i*r
2 楼
【 以下文字转载自 Programming 讨论区 】
发信人: iiiir (哎呀我最牛), 信区: Programming
标 题: 【包子求助】20M*20M的loop怎么搞？
发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
我有2个文件A和B，各自约20 million lines
现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
之间。
比较直观的是，我写2个loop
for i=1:length(A)
for k=1:length(B)
do what I want to check...
end
end
问题是20M*20M似乎太长，我的电脑化了2个小时，才算到i=600左右，也就是跑了600*
20M loops
请问大虾们，怎么整。我用的是matlab，C也可以用，python也可以考虑。。。
多谢了，解决问题的发包子
发信人: iiiir (哎呀我最牛), 信区: Programming
标 题: 【包子求助】20M*20M的loop怎么搞？
发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
我有2个文件A和B，各自约20 million lines
现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
之间。
比较直观的是，我写2个loop
for i=1:length(A)
for k=1:length(B)
do what I want to check...
end
end
问题是20M*20M似乎太长，我的电脑化了2个小时，才算到i=600左右，也就是跑了600*
20M loops
请问大虾们，怎么整。我用的是matlab，C也可以用，python也可以考虑。。。
多谢了，解决问题的发包子
t*r
3 楼
没认出来，：（
G*h
4 楼
用 C 吧
这俩文件不大呀，整数的话，A 80MB B 160MB 空间就都读内存里了
然后就傻 loop 也不会花太久
【在 i***r 的大作中提到】
: 【 以下文字转载自 Programming 讨论区 】
: 发信人: iiiir (哎呀我最牛), 信区: Programming
: 标 题: 【包子求助】20M*20M的loop怎么搞？
: 发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
: 我有2个文件A和B，各自约20 million lines
: 现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
: 之间。
: 比较直观的是，我写2个loop
: for i=1:length(A)
: for k=1:length(B)
这俩文件不大呀，整数的话，A 80MB B 160MB 空间就都读内存里了
然后就傻 loop 也不会花太久
【在 i***r 的大作中提到】
: 【 以下文字转载自 Programming 讨论区 】
: 发信人: iiiir (哎呀我最牛), 信区: Programming
: 标 题: 【包子求助】20M*20M的loop怎么搞？
: 发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
: 我有2个文件A和B，各自约20 million lines
: 现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
: 之间。
: 比较直观的是，我写2个loop
: for i=1:length(A)
: for k=1:length(B)
d*y
5 楼
一开始出来也没认出来。
下面一集看着有点像。
搜了一把，还真是她。
老得像大妈了。
【在 t******r 的大作中提到】
: 没认出来，：（
下面一集看着有点像。
搜了一把，还真是她。
老得像大妈了。
【在 t******r 的大作中提到】
: 没认出来，：（
w*w
6 楼
用perl,方便且较快。
把第一个文件的数分成区间，做个hash，key是区间，value是个list,即所有在区间里
的数。然后把hash sort一下，读第二个文件，再到hash里找区间，处理一下重叠的情
况就行了。
把第一个文件的数分成区间，做个hash，key是区间，value是个list,即所有在区间里
的数。然后把hash sort一下，读第二个文件，再到hash里找区间，处理一下重叠的情
况就行了。
l*r
7 楼
还好吧，不算太大妈
话说大妈就是女子的归宿啊
【在 d********y 的大作中提到】
: 一开始出来也没认出来。
: 下面一集看着有点像。
: 搜了一把，还真是她。
: 老得像大妈了。
话说大妈就是女子的归宿啊
【在 d********y 的大作中提到】
: 一开始出来也没认出来。
: 下面一集看着有点像。
: 搜了一把，还真是她。
: 老得像大妈了。
i*r
8 楼
perl不会。。。可能最后要用2楼的方法。
我的文件实际上要复杂些，除了数字还有gene序列在后面，找到交叉后我要去找对应的
gene信息。
我现在是一行一行的scan A（然后去和B的每一行找交叉），可能这就是程序慢的主要
因素？
包子随后奉上，先多谢了
我的文件实际上要复杂些，除了数字还有gene序列在后面，找到交叉后我要去找对应的
gene信息。
我现在是一行一行的scan A（然后去和B的每一行找交叉），可能这就是程序慢的主要
因素？
包子随后奉上，先多谢了
L*7
9 楼
我还没看到她出现
不过她的长相是容易见老那种
不过她的长相是容易见老那种
l*n
10 楼
这个活一定不能用matlab, 那会慢死. 另外A的任意行前两个数字构成的区间都不重叠
么?
【在 i***r 的大作中提到】
: 【 以下文字转载自 Programming 讨论区 】
: 发信人: iiiir (哎呀我最牛), 信区: Programming
: 标 题: 【包子求助】20M*20M的loop怎么搞？
: 发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
: 我有2个文件A和B，各自约20 million lines
: 现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
: 之间。
: 比较直观的是，我写2个loop
: for i=1:length(A)
: for k=1:length(B)
么?
【在 i***r 的大作中提到】
: 【 以下文字转载自 Programming 讨论区 】
: 发信人: iiiir (哎呀我最牛), 信区: Programming
: 标 题: 【包子求助】20M*20M的loop怎么搞？
: 发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
: 我有2个文件A和B，各自约20 million lines
: 现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
: 之间。
: 比较直观的是，我写2个loop
: for i=1:length(A)
: for k=1:length(B)
r*y
11 楼
当年忘了看她那个片子了 很pp！
【在 d********y 的大作中提到】
: 那个独立电影投资人珊莉。
: 真老了。
【在 d********y 的大作中提到】
: 那个独立电影投资人珊莉。
: 真老了。
s*z
13 楼
还好吧，毕竟她已经41了~
当年的惊鸿仙子真是pp啊
【在 d********y 的大作中提到】
: 一开始出来也没认出来。
: 下面一集看着有点像。
: 搜了一把，还真是她。
: 老得像大妈了。
当年的惊鸿仙子真是pp啊
【在 d********y 的大作中提到】
: 一开始出来也没认出来。
: 下面一集看着有点像。
: 搜了一把，还真是她。
: 老得像大妈了。
w*w
14 楼
给SNP annotate？你这个是单个chromosome还是whole genome?
【在 i***r 的大作中提到】
: perl不会。。。可能最后要用2楼的方法。
: 我的文件实际上要复杂些，除了数字还有gene序列在后面，找到交叉后我要去找对应的
: gene信息。
: 我现在是一行一行的scan A（然后去和B的每一行找交叉），可能这就是程序慢的主要
: 因素？
: 包子随后奉上，先多谢了
【在 i***r 的大作中提到】
: perl不会。。。可能最后要用2楼的方法。
: 我的文件实际上要复杂些，除了数字还有gene序列在后面，找到交叉后我要去找对应的
: gene信息。
: 我现在是一行一行的scan A（然后去和B的每一行找交叉），可能这就是程序慢的主要
: 因素？
: 包子随后奉上，先多谢了
w*w
15 楼
这个东东我写过个C# library,把gene按chromosome分好后排序，然后再搜SNP,不应该
超过一个小时。
超过一个小时。
n*7
16 楼
恩，感觉lz是要做这个
我也做过类似的，sort是关键，还有就是chr by chr来做
【在 w****w 的大作中提到】
: 这个东东我写过个C# library,把gene按chromosome分好后排序，然后再搜SNP,不应该
: 超过一个小时。
我也做过类似的，sort是关键，还有就是chr by chr来做
【在 w****w 的大作中提到】
: 这个东东我写过个C# library,把gene按chromosome分好后排序，然后再搜SNP,不应该
: 超过一个小时。
i*r
17 楼
谢谢大家，我确实是要搜SNP，chromosome都排好序了，sort过了，稍微优化了一下，
一会儿晚上看看跑了多少。
包子已发，谢谢
PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
一会儿晚上看看跑了多少。
包子已发，谢谢
PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
o*n
18 楼
借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
【在 G*****h 的大作中提到】
: 用 C 吧
: 这俩文件不大呀，整数的话，A 80MB B 160MB 空间就都读内存里了
: 然后就傻 loop 也不会花太久
有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
【在 G*****h 的大作中提到】
: 用 C 吧
: 这俩文件不大呀，整数的话，A 80MB B 160MB 空间就都读内存里了
: 然后就傻 loop 也不会花太久
y*d
19 楼
你是说 loop grid 花时间？弄个 dirty set 就行了吧 ...都不用，算 xy 的 grid 时
候记录 index 就成了？
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
候记录 index 就成了？
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
G*h
20 楼
有序序列，二分搜索好了
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
o*n
21 楼
我去查查dirty set。
我没太搞懂你说的第二个
我现在是两层loop grid，用matlab的矩阵比较把index搞出来，类似
for jx=1:nx
ix = xt>x(jx) & xt<=x(jx)+dx;
for jy=1:ny
iy = yt>y(jy) & yt<=y(jy)+dy;
ixy = ix&iy;
#use ixy to process data ...
end
end
xt,yt是(x,y)坐标的二维数据nx*ny，其实是我有一个矩阵数据A，每个数据点对应一个
(x,y)坐标，A(i,j) -> (x,y)。上面code中的x,y是新建的均匀grid，每个间距dx,dy。
我不太理解你说的意思。
【在 y***d 的大作中提到】
: 你是说 loop grid 花时间？弄个 dirty set 就行了吧 ...都不用，算 xy 的 grid 时
: 候记录 index 就成了？
:
: monte
我没太搞懂你说的第二个
我现在是两层loop grid，用matlab的矩阵比较把index搞出来，类似
for jx=1:nx
ix = xt>x(jx) & xt<=x(jx)+dx;
for jy=1:ny
iy = yt>y(jy) & yt<=y(jy)+dy;
ixy = ix&iy;
#use ixy to process data ...
end
end
xt,yt是(x,y)坐标的二维数据nx*ny，其实是我有一个矩阵数据A，每个数据点对应一个
(x,y)坐标，A(i,j) -> (x,y)。上面code中的x,y是新建的均匀grid，每个间距dx,dy。
我不太理解你说的意思。
【在 y***d 的大作中提到】
: 你是说 loop grid 花时间？弄个 dirty set 就行了吧 ...都不用，算 xy 的 grid 时
: 候记录 index 就成了？
:
: monte
o*n
22 楼
我猜如果我用matlab里矩阵的比较时，比如A(i,j)>B(i,j)返回所有满足条件的数据行
列index，应该用的是很优化的算发，我的理解是matlab建议用户尽量用矩阵形式避免
用loop就是他内部优化的结果。但我要对每个grid点和原数据都做一次这个比较，所以
费时，请参见我贴出的code。
【在 G*****h 的大作中提到】
: 有序序列，二分搜索好了
:
: monte
列index，应该用的是很优化的算发，我的理解是matlab建议用户尽量用矩阵形式避免
用loop就是他内部优化的结果。但我要对每个grid点和原数据都做一次这个比较，所以
费时，请参见我贴出的code。
【在 G*****h 的大作中提到】
: 有序序列，二分搜索好了
:
: monte
y*d
23 楼
俺完全没看懂... 当我没说...
【在 o**n 的大作中提到】
: 我去查查dirty set。
: 我没太搞懂你说的第二个
: 我现在是两层loop grid，用matlab的矩阵比较把index搞出来，类似
: for jx=1:nx
: ix = xt>x(jx) & xt<=x(jx)+dx;
: for jy=1:ny
: iy = yt>y(jy) & yt<=y(jy)+dy;
: ixy = ix&iy;
: #use ixy to process data ...
: end
【在 o**n 的大作中提到】
: 我去查查dirty set。
: 我没太搞懂你说的第二个
: 我现在是两层loop grid，用matlab的矩阵比较把index搞出来，类似
: for jx=1:nx
: ix = xt>x(jx) & xt<=x(jx)+dx;
: for jy=1:ny
: iy = yt>y(jy) & yt<=y(jy)+dy;
: ixy = ix&iy;
: #use ixy to process data ...
: end
l*y
24 楼
竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
作业一律零分。
【在 i***r 的大作中提到】
: 谢谢大家，我确实是要搜SNP，chromosome都排好序了，sort过了，稍微优化了一下，
: 一会儿晚上看看跑了多少。
: 包子已发，谢谢
: PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
作业一律零分。
【在 i***r 的大作中提到】
: 谢谢大家，我确实是要搜SNP，chromosome都排好序了，sort过了，稍微优化了一下，
: 一会儿晚上看看跑了多少。
: 包子已发，谢谢
: PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
o*n
25 楼
请教我前面那个例子怎么避免二重循环或在matlab里怎么能做的更有效些？谢谢！
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
G*h
26 楼
matlab 也就在学校用用
不会也没事
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
不会也没事
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
j*u
27 楼
这里的人C专家众多，但是matlab的编程方式很不一样，很多时候是靠点小聪明，而不
是什么算法的。
你的问题，我的理解是平面上很多点，你知道每个点的坐标，现在在平面上画一个棋盘
格子，然后需要知道每个格子里面是那些点？如果是这样的话，可以比较简单的用
matlab的矩阵来进行计算：
1. 每个点的x和y是可以分开考虑的（因为x轴垂直于y轴）
2. 根据你的棋盘格子离散化x和y矩阵，例如x=x0+nx*dx，这个nx就是x这个点在棋盘格
子里的位置，nx是正整数
3. 新得倒的nx和ny矩阵就是你(x,y)坐标矩阵在棋盘格子的新坐标矩阵。利用逻辑矩阵
，很容易找出棋盘格子里任意一格包含的所有点的坐标
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
是什么算法的。
你的问题，我的理解是平面上很多点，你知道每个点的坐标，现在在平面上画一个棋盘
格子，然后需要知道每个格子里面是那些点？如果是这样的话，可以比较简单的用
matlab的矩阵来进行计算：
1. 每个点的x和y是可以分开考虑的（因为x轴垂直于y轴）
2. 根据你的棋盘格子离散化x和y矩阵，例如x=x0+nx*dx，这个nx就是x这个点在棋盘格
子里的位置，nx是正整数
3. 新得倒的nx和ny矩阵就是你(x,y)坐标矩阵在棋盘格子的新坐标矩阵。利用逻辑矩阵
，很容易找出棋盘格子里任意一格包含的所有点的坐标
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
j*u
28 楼
有点糊涂了，你最终的目的是什么？
【在 o**n 的大作中提到】
: 我去查查dirty set。
: 我没太搞懂你说的第二个
: 我现在是两层loop grid，用matlab的矩阵比较把index搞出来，类似
: for jx=1:nx
: ix = xt>x(jx) & xt<=x(jx)+dx;
: for jy=1:ny
: iy = yt>y(jy) & yt<=y(jy)+dy;
: ixy = ix&iy;
: #use ixy to process data ...
: end
【在 o**n 的大作中提到】
: 我去查查dirty set。
: 我没太搞懂你说的第二个
: 我现在是两层loop grid，用matlab的矩阵比较把index搞出来，类似
: for jx=1:nx
: ix = xt>x(jx) & xt<=x(jx)+dx;
: for jy=1:ny
: iy = yt>y(jy) & yt<=y(jy)+dy;
: ixy = ix&iy;
: #use ixy to process data ...
: end
G*7
29 楼
no need to bother with quadtree
float / mesh_grid_size = index of containing cell
【在 G*****h 的大作中提到】
: 有序序列，二分搜索好了
:
: monte
float / mesh_grid_size = index of containing cell
【在 G*****h 的大作中提到】
: 有序序列，二分搜索好了
:
: monte
G*7
30 楼
function fast_grid_lookup
%% prepare the data
% randomly generate the query points
num_points = 20e6;
points = rand(num_points,2);
% specify the evenly-spaced mesh grid
num_divs = 1028; % per side
num_cells = num_divs^2;
grid_spacing = 1/(num_divs-1);
% you do not have to instantiate the grid points by
% grid = meshgrid(linspace(0, 1, num_div), linspace(0, 1, num_div));
%% find the enclosing cell of each point
tic;
points_in_cell = ceil(points/grid_spacing)+1; % damn you, 1-based matlab
toc;
%% one-shot query: which points are within a given cell,
% for instance, the cell at the center
i=num_divs/2;
j=num_divs/2;
cell_idx = [i, j];
tic;
points_idx_1 = find( points_in_cell(:,1)==cell_idx(1) & points_in_cell(:,2)=
=cell_idx(2) );
toc;
%% amalgamate the cell-point membership into a sparse matrix
tic;
cell_indices = sub2ind([num_divs,num_divs], points_in_cell(:,1), points_in_
cell(:,2));
point_indices = 1:num_points;
indicator = ones(1, num_points);
cell_point = sparse(cell_indices, point_indices, indicator, num_cells, num_
points);
toc;
% query like this
tic;
points_idx_2 = cell_point(sub2ind([num_divs, num_divs], i, j), :);
toc;
assert( all(transpose(points_idx_1) == find(points_idx_2)) );
end
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
%% prepare the data
% randomly generate the query points
num_points = 20e6;
points = rand(num_points,2);
% specify the evenly-spaced mesh grid
num_divs = 1028; % per side
num_cells = num_divs^2;
grid_spacing = 1/(num_divs-1);
% you do not have to instantiate the grid points by
% grid = meshgrid(linspace(0, 1, num_div), linspace(0, 1, num_div));
%% find the enclosing cell of each point
tic;
points_in_cell = ceil(points/grid_spacing)+1; % damn you, 1-based matlab
toc;
%% one-shot query: which points are within a given cell,
% for instance, the cell at the center
i=num_divs/2;
j=num_divs/2;
cell_idx = [i, j];
tic;
points_idx_1 = find( points_in_cell(:,1)==cell_idx(1) & points_in_cell(:,2)=
=cell_idx(2) );
toc;
%% amalgamate the cell-point membership into a sparse matrix
tic;
cell_indices = sub2ind([num_divs,num_divs], points_in_cell(:,1), points_in_
cell(:,2));
point_indices = 1:num_points;
indicator = ones(1, num_points);
cell_point = sparse(cell_indices, point_indices, indicator, num_cells, num_
points);
toc;
% query like this
tic;
points_idx_2 = cell_point(sub2ind([num_divs, num_divs], i, j), :);
toc;
assert( all(transpose(points_idx_1) == find(points_idx_2)) );
end
monte
【在 o**n 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(x,y)坐标矩阵，float type，x,y都没
: 有规则顺序。我再建一个x,y mesh grid的矩阵，间距是均匀的或有一定规律比如线性
: 增加等，要把前面那个乱的(x,y)矩阵放入grid里，会有多个(x,y)落在同一个grid的情
: 况，把每个grid中有的(x,y)在原矩阵中的行列index记录下来。其实就是个二维monte
: carlo的东西。我现在也是用loop硬来，目前用matlab，用C肯定快一些。数组也不是特
: 别大，run一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
G*7
31 楼
1) for-free. the time taken is minimal. for 20 million points in 1024^2
grids, it took me < 5 seconds.
2) trivial to extend to cases with dx!=dy.
【在 G*****7 的大作中提到】
: function fast_grid_lookup
: %% prepare the data
: % randomly generate the query points
: num_points = 20e6;
: points = rand(num_points,2);
: % specify the evenly-spaced mesh grid
: num_divs = 1028; % per side
: num_cells = num_divs^2;
: grid_spacing = 1/(num_divs-1);
: % you do not have to instantiate the grid points by
grids, it took me < 5 seconds.
2) trivial to extend to cases with dx!=dy.
【在 G*****7 的大作中提到】
: function fast_grid_lookup
: %% prepare the data
: % randomly generate the query points
: num_points = 20e6;
: points = rand(num_points,2);
: % specify the evenly-spaced mesh grid
: num_divs = 1028; % per side
: num_cells = num_divs^2;
: grid_spacing = 1/(num_divs-1);
: % you do not have to instantiate the grid points by
G*7
32 楼
i m gonna repeat what ppl said in programming:
1) if you have to explicitly do for-loop and linearly search for the
interval membership, use mex to embed c/cpp code in matlab.
2) if you got time, check out the search acceleration data structures
provided by, say cgal, such as segment tree, and write mex+cpp code that
takes advantage of that.
【在 i***r 的大作中提到】
: 谢谢大家，我确实是要搜SNP，chromosome都排好序了，sort过了，稍微优化了一下，
: 一会儿晚上看看跑了多少。
: 包子已发，谢谢
: PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
1) if you have to explicitly do for-loop and linearly search for the
interval membership, use mex to embed c/cpp code in matlab.
2) if you got time, check out the search acceleration data structures
provided by, say cgal, such as segment tree, and write mex+cpp code that
takes advantage of that.
【在 i***r 的大作中提到】
: 谢谢大家，我确实是要搜SNP，chromosome都排好序了，sort过了，稍微优化了一下，
: 一会儿晚上看看跑了多少。
: 包子已发，谢谢
: PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
o*n
33 楼
多谢思考我的问题和读我的code，先包子一个。我觉得你的方法好像比我的code那个好
一些，我去试试看。
我的问题表达得不太好，让很多人糊涂了。我再试着说一遍：
有一个二维面上离散化的物理量A=A(i,j)，0过一些过程后对应到另一个二维面上，每个点在新二维面上坐标变了且与在初始面上的
坐标没有单调关系，i.e.每个A(i,j)在新面的坐标为(x(i,j),y(i,j))。我现在需要知
道这个物理量在新的二维面上的分布。所以在新面上重新做grid，这样就不可避免有多
个点落在同一个grid格子里，如果是这样，就把这个格子里的所有点的物理量值相加，
最后再做一个二维内插平滑化，就能得出在新面上的物理量分布了。其实是个monte
carlo的东西，当然在初始面离散的越细（点越多），最后在新面上的分布就越准确连
续平滑一些。
希望我说明白了。
【在 j**u 的大作中提到】
: 这里的人C专家众多，但是matlab的编程方式很不一样，很多时候是靠点小聪明，而不
: 是什么算法的。
: 你的问题，我的理解是平面上很多点，你知道每个点的坐标，现在在平面上画一个棋盘
: 格子，然后需要知道每个格子里面是那些点？如果是这样的话，可以比较简单的用
: matlab的矩阵来进行计算：
: 1. 每个点的x和y是可以分开考虑的（因为x轴垂直于y轴）
: 2. 根据你的棋盘格子离散化x和y矩阵，例如x=x0+nx*dx，这个nx就是x这个点在棋盘格
: 子里的位置，nx是正整数
: 3. 新得倒的nx和ny矩阵就是你(x,y)坐标矩阵在棋盘格子的新坐标矩阵。利用逻辑矩阵
: ，很容易找出棋盘格子里任意一格包含的所有点的坐标
一些，我去试试看。
我的问题表达得不太好，让很多人糊涂了。我再试着说一遍：
有一个二维面上离散化的物理量A=A(i,j)，0过一些过程后对应到另一个二维面上，每个点在新二维面上坐标变了且与在初始面上的
坐标没有单调关系，i.e.每个A(i,j)在新面的坐标为(x(i,j),y(i,j))。我现在需要知
道这个物理量在新的二维面上的分布。所以在新面上重新做grid，这样就不可避免有多
个点落在同一个grid格子里，如果是这样，就把这个格子里的所有点的物理量值相加，
最后再做一个二维内插平滑化，就能得出在新面上的物理量分布了。其实是个monte
carlo的东西，当然在初始面离散的越细（点越多），最后在新面上的分布就越准确连
续平滑一些。
希望我说明白了。
【在 j**u 的大作中提到】
: 这里的人C专家众多，但是matlab的编程方式很不一样，很多时候是靠点小聪明，而不
: 是什么算法的。
: 你的问题，我的理解是平面上很多点，你知道每个点的坐标，现在在平面上画一个棋盘
: 格子，然后需要知道每个格子里面是那些点？如果是这样的话，可以比较简单的用
: matlab的矩阵来进行计算：
: 1. 每个点的x和y是可以分开考虑的（因为x轴垂直于y轴）
: 2. 根据你的棋盘格子离散化x和y矩阵，例如x=x0+nx*dx，这个nx就是x这个点在棋盘格
: 子里的位置，nx是正整数
: 3. 新得倒的nx和ny矩阵就是你(x,y)坐标矩阵在棋盘格子的新坐标矩阵。利用逻辑矩阵
: ，很容易找出棋盘格子里任意一格包含的所有点的坐标
o*n
34 楼
thanks for your code! baozi first. I'll read the code and try it.
【在 G*****7 的大作中提到】
: function fast_grid_lookup
: %% prepare the data
: % randomly generate the query points
: num_points = 20e6;
: points = rand(num_points,2);
: % specify the evenly-spaced mesh grid
: num_divs = 1028; % per side
: num_cells = num_divs^2;
: grid_spacing = 1/(num_divs-1);
: % you do not have to instantiate the grid points by
【在 G*****7 的大作中提到】
: function fast_grid_lookup
: %% prepare the data
: % randomly generate the query points
: num_points = 20e6;
: points = rand(num_points,2);
: % specify the evenly-spaced mesh grid
: num_divs = 1028; % per side
: num_cells = num_divs^2;
: grid_spacing = 1/(num_divs-1);
: % you do not have to instantiate the grid points by
o*n
35 楼
just looked through your code. I think you and jzxu used the same concept
which is a much better way than what I did. And your code can be simply
revised to deal with non-uniform grid. Thanks for your help!!
【在 G*****7 的大作中提到】
: function fast_grid_lookup
: %% prepare the data
: % randomly generate the query points
: num_points = 20e6;
: points = rand(num_points,2);
: % specify the evenly-spaced mesh grid
: num_divs = 1028; % per side
: num_cells = num_divs^2;
: grid_spacing = 1/(num_divs-1);
: % you do not have to instantiate the grid points by
which is a much better way than what I did. And your code can be simply
revised to deal with non-uniform grid. Thanks for your help!!
【在 G*****7 的大作中提到】
: function fast_grid_lookup
: %% prepare the data
: % randomly generate the query points
: num_points = 20e6;
: points = rand(num_points,2);
: % specify the evenly-spaced mesh grid
: num_divs = 1028; % per side
: num_cells = num_divs^2;
: grid_spacing = 1/(num_divs-1);
: % you do not have to instantiate the grid points by
G*7
36 楼
a better way to explain your problem to a programmer is to abstract out the
non-essentials (离散点经过一些过程, 物理量值相加,二维内插平滑化, 物理量分布,
monte carlo...), and tell them what is your input (points and grid) and
what is the desired output (cell-point membership).
【在 o**n 的大作中提到】
: 多谢思考我的问题和读我的code，先包子一个。我觉得你的方法好像比我的code那个好
: 一些，我去试试看。
: 我的问题表达得不太好，让很多人糊涂了。我再试着说一遍：
: 有一个二维面上离散化的物理量A=A(i,j)，0: 过一些过程后对应到另一个二维面上，每个点在新二维面上坐标变了且与在初始面上的
: 坐标没有单调关系，i.e.每个A(i,j)在新面的坐标为(x(i,j),y(i,j))。我现在需要知
: 道这个物理量在新的二维面上的分布。所以在新面上重新做grid，这样就不可避免有多
: 个点落在同一个grid格子里，如果是这样，就把这个格子里的所有点的物理量值相加，
: 最后再做一个二维内插平滑化，就能得出在新面上的物理量分布了。其实是个monte
: carlo的东西，当然在初始面离散的越细（点越多），最后在新面上的分布就越准确连
non-essentials (离散点经过一些过程, 物理量值相加,二维内插平滑化, 物理量分布,
monte carlo...), and tell them what is your input (points and grid) and
what is the desired output (cell-point membership).
【在 o**n 的大作中提到】
: 多谢思考我的问题和读我的code，先包子一个。我觉得你的方法好像比我的code那个好
: 一些，我去试试看。
: 我的问题表达得不太好，让很多人糊涂了。我再试着说一遍：
: 有一个二维面上离散化的物理量A=A(i,j)，0: 过一些过程后对应到另一个二维面上，每个点在新二维面上坐标变了且与在初始面上的
: 坐标没有单调关系，i.e.每个A(i,j)在新面的坐标为(x(i,j),y(i,j))。我现在需要知
: 道这个物理量在新的二维面上的分布。所以在新面上重新做grid，这样就不可避免有多
: 个点落在同一个grid格子里，如果是这样，就把这个格子里的所有点的物理量值相加，
: 最后再做一个二维内插平滑化，就能得出在新面上的物理量分布了。其实是个monte
: carlo的东西，当然在初始面离散的越细（点越多），最后在新面上的分布就越准确连
o*n
37 楼
thanks for your advise!
the
布,
【在 G*****7 的大作中提到】
: a better way to explain your problem to a programmer is to abstract out the
: non-essentials (离散点经过一些过程, 物理量值相加,二维内插平滑化, 物理量分布,
: monte carlo...), and tell them what is your input (points and grid) and
: what is the desired output (cell-point membership).
the
布,
【在 G*****7 的大作中提到】
: a better way to explain your problem to a programmer is to abstract out the
: non-essentials (离散点经过一些过程, 物理量值相加,二维内插平滑化, 物理量分布,
: monte carlo...), and tell them what is your input (points and grid) and
: what is the desired output (cell-point membership).
z*r
38 楼
bedtools
【在 i***r 的大作中提到】
: 【 以下文字转载自 Programming 讨论区 】
: 发信人: iiiir (哎呀我最牛), 信区: Programming
: 标 题: 【包子求助】20M*20M的loop怎么搞？
: 发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
: 我有2个文件A和B，各自约20 million lines
: 现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
: 之间。
: 比较直观的是，我写2个loop
: for i=1:length(A)
: for k=1:length(B)
【在 i***r 的大作中提到】
: 【 以下文字转载自 Programming 讨论区 】
: 发信人: iiiir (哎呀我最牛), 信区: Programming
: 标 题: 【包子求助】20M*20M的loop怎么搞？
: 发信站: BBS 未名空间站 (Fri Sep 14 14:51:39 2012, 美东)
: 我有2个文件A和B，各自约20 million lines
: 现在我要check的是，B的每一行中的第一个数字，是否落入A的某一行的第一第二数字
: 之间。
: 比较直观的是，我写2个loop
: for i=1:length(A)
: for k=1:length(B)
N8
39 楼
河马新马甲?
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
相关阅读
GNU screen 下面 vim 没有颜色？请高手解答grep和ls的两个问题linux下能处理plain tex文档的软件都有什么？ubuntu 8.10对video memory要求很高吗?土问平衡disk i/o的方法Gnome又开始搞global menu bar了啊Gnome-Do 打开一个软件在docky模式下比classic模式下慢多了free linux book死机reset后到底需不需要force fsch check?我觉得dreamlinux 挺不错的哪个 desktop 能自定义这几个热键？Ubuntu 或者 Debian based 的哪个 distro 最省资源？subversion 初学者问题dreamLinux 3.5 终于出了AWN里面能设windows list么请问我应该安装哪个linux的dist？perl能不能一次把一个str中的a替换成x，b替换成y？adobe flash player 10PCMan X的字体/编码设置问题Ubuntu 9.04's New Login Screen