i***r 发帖数: 1035 | 1 【 以下文字转载自 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也可以考虑。。。
多谢了,解决问题的发包子 |
G*****h 发帖数: 33134 | 2 用 C 吧
这俩文件不大呀,整数的话,A 80MB B 160MB 空间就都读内存里了
然后就傻 loop 也不会花太久
【在 i***r 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 【 以下文字转载自 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)
|
w****w 发帖数: 521 | 3 用perl,方便且较快。
把第一个文件的数分成区间,做个hash,key是区间,value是个list,即所有在区间里
的数。然后把hash sort一下,读第二个文件,再到hash里找区间,处理一下重叠的情
况就行了。 |
i***r 发帖数: 1035 | 4 perl不会。。。可能最后要用2楼的方法。
我的文件实际上要复杂些,除了数字还有gene序列在后面,找到交叉后我要去找对应的
gene信息。
我现在是一行一行的scan A(然后去和B的每一行找交叉),可能这就是程序慢的主要
因素?
包子随后奉上,先多谢了 |
l******n 发帖数: 1683 | 5 这个活一定不能用matlab, 那会慢死. 另外A的任意行前两个数字构成的区间都不重叠
么?
【在 i***r 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 【 以下文字转载自 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)
|
w*r 发帖数: 2421 | |
w****w 发帖数: 521 | 7 给SNP annotate?你这个是单个chromosome还是whole genome?
【在 i***r 的大作中提到】![](/moin_static193/solenoid/img/up.png) : perl不会。。。可能最后要用2楼的方法。 : 我的文件实际上要复杂些,除了数字还有gene序列在后面,找到交叉后我要去找对应的 : gene信息。 : 我现在是一行一行的scan A(然后去和B的每一行找交叉),可能这就是程序慢的主要 : 因素? : 包子随后奉上,先多谢了
|
w****w 发帖数: 521 | 8 这个东东我写过个C# library,把gene按chromosome分好后排序,然后再搜SNP,不应该
超过一个小时。 |
n******7 发帖数: 12463 | 9 恩,感觉lz是要做这个
我也做过类似的,sort是关键,还有就是chr by chr来做
【在 w****w 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 这个东东我写过个C# library,把gene按chromosome分好后排序,然后再搜SNP,不应该 : 超过一个小时。
|
i***r 发帖数: 1035 | 10 谢谢大家,我确实是要搜SNP,chromosome都排好序了,sort过了,稍微优化了一下,
一会儿晚上看看跑了多少。
包子已发,谢谢
PS 这里有没有new jersey附近做生物信息的?pm我大家认识下? |
|
|
o**n 发帖数: 1249 | 11 借机问一个问题。我有一个较大的二维矩阵,是(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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 用 C 吧 : 这俩文件不大呀,整数的话,A 80MB B 160MB 空间就都读内存里了 : 然后就傻 loop 也不会花太久
|
y***d 发帖数: 2330 | 12 你是说 loop grid 花时间?弄个 dirty set 就行了吧 ...都不用,算 xy 的 grid 时
候记录 index 就成了?
monte
【在 o**n 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 借机问一个问题。我有一个较大的二维矩阵,是(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 发帖数: 33134 | 13 有序序列,二分搜索好了
monte
【在 o**n 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 借机问一个问题。我有一个较大的二维矩阵,是(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 发帖数: 1249 | 14 我去查查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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 你是说 loop grid 花时间?弄个 dirty set 就行了吧 ...都不用,算 xy 的 grid 时 : 候记录 index 就成了? : : monte
|
o**n 发帖数: 1249 | 15 我猜如果我用matlab里矩阵的比较时,比如A(i,j)>B(i,j)返回所有满足条件的数据行
列index,应该用的是很优化的算发,我的理解是matlab建议用户尽量用矩阵形式避免
用loop就是他内部优化的结果。但我要对每个grid点和原数据都做一次这个比较,所以
费时,请参见我贴出的code。
【在 G*****h 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 有序序列,二分搜索好了 : : monte
|
y***d 发帖数: 2330 | 16 俺完全没看懂... 当我没说...
【在 o**n 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 我去查查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 发帖数: 4671 | 17 竟然用 matlab 做大规模的多重循环?这都是在哪里学的 matlab?要是在我的班上,
作业一律零分。
【在 i***r 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 谢谢大家,我确实是要搜SNP,chromosome都排好序了,sort过了,稍微优化了一下, : 一会儿晚上看看跑了多少。 : 包子已发,谢谢 : PS 这里有没有new jersey附近做生物信息的?pm我大家认识下?
|
o**n 发帖数: 1249 | 18 请教我前面那个例子怎么避免二重循环或在matlab里怎么能做的更有效些?谢谢!
【在 l***y 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 竟然用 matlab 做大规模的多重循环?这都是在哪里学的 matlab?要是在我的班上, : 作业一律零分。
|
G*****h 发帖数: 33134 | 19 matlab 也就在学校用用
不会也没事
【在 l***y 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 竟然用 matlab 做大规模的多重循环?这都是在哪里学的 matlab?要是在我的班上, : 作业一律零分。
|
j**u 发帖数: 6059 | 20 这里的人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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 借机问一个问题。我有一个较大的二维矩阵,是(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 发帖数: 6059 | 21 有点糊涂了,你最终的目的是什么?
【在 o**n 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 我去查查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 发帖数: 1759 | 22 no need to bother with quadtree
float / mesh_grid_size = index of containing cell
【在 G*****h 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 有序序列,二分搜索好了 : : monte
|
G*****7 发帖数: 1759 | 23 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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 借机问一个问题。我有一个较大的二维矩阵,是(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 发帖数: 1759 | 24 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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 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 发帖数: 1759 | 25 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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 谢谢大家,我确实是要搜SNP,chromosome都排好序了,sort过了,稍微优化了一下, : 一会儿晚上看看跑了多少。 : 包子已发,谢谢 : PS 这里有没有new jersey附近做生物信息的?pm我大家认识下?
|
o**n 发帖数: 1249 | 26 多谢思考我的问题和读我的code,先包子一个。我觉得你的方法好像比我的code那个好
一些,我去试试看。
我的问题表达得不太好,让很多人糊涂了。我再试着说一遍:
有一个二维面上离散化的物理量A=A(i,j),0
过一些过程后对应到另一个二维面上,每个点在新二维面上坐标变了且与在初始面上的
坐标没有单调关系,i.e.每个A(i,j)在新面的坐标为(x(i,j),y(i,j))。我现在需要知
道这个物理量在新的二维面上的分布。所以在新面上重新做grid,这样就不可避免有多
个点落在同一个grid格子里,如果是这样,就把这个格子里的所有点的物理量值相加,
最后再做一个二维内插平滑化,就能得出在新面上的物理量分布了。其实是个monte
carlo的东西,当然在初始面离散的越细(点越多),最后在新面上的分布就越准确连
续平滑一些。
希望我说明白了。
【在 j**u 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 这里的人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 发帖数: 1249 | 27 thanks for your code! baozi first. I'll read the code and try it.
【在 G*****7 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 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 发帖数: 1249 | 28 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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 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 发帖数: 1759 | 29 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 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 多谢思考我的问题和读我的code,先包子一个。我觉得你的方法好像比我的code那个好 : 一些,我去试试看。 : 我的问题表达得不太好,让很多人糊涂了。我再试着说一遍: : 有一个二维面上离散化的物理量A=A(i,j),0: 过一些过程后对应到另一个二维面上,每个点在新二维面上坐标变了且与在初始面上的 : 坐标没有单调关系,i.e.每个A(i,j)在新面的坐标为(x(i,j),y(i,j))。我现在需要知 : 道这个物理量在新的二维面上的分布。所以在新面上重新做grid,这样就不可避免有多 : 个点落在同一个grid格子里,如果是这样,就把这个格子里的所有点的物理量值相加, : 最后再做一个二维内插平滑化,就能得出在新面上的物理量分布了。其实是个monte : carlo的东西,当然在初始面离散的越细(点越多),最后在新面上的分布就越准确连
|
o**n 发帖数: 1249 | 30 thanks for your advise!
the
布,
【在 G*****7 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 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 发帖数: 131 | 31 bedtools
【在 i***r 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 【 以下文字转载自 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 发帖数: 110 | 32 河马新马甲?
【在 l***y 的大作中提到】![](/moin_static193/solenoid/img/up.png) : 竟然用 matlab 做大规模的多重循环?这都是在哪里学的 matlab?要是在我的班上, : 作业一律零分。
|