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 的大作中提到】
: 【 以下文字转载自 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 的大作中提到】
: 【 以下文字转载自 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 的大作中提到】
: 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 的大作中提到】
: 这个东东我写过个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 的大作中提到】
: 用 C 吧
: 这俩文件不大呀，整数的话，A 80MB B 160MB 空间就都读内存里了
: 然后就傻 loop 也不会花太久
|
y***d
发帖数: 2330 | 12 你是说 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一次大概几分钟时间，但我要做很多次这个事情，想请教大家更好的方法。
|
G*****h
发帖数: 33134 | 13 有序序列，二分搜索好了
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
发帖数: 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 的大作中提到】
: 你是说 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 的大作中提到】
: 有序序列，二分搜索好了
:
: monte
|
y***d
发帖数: 2330 | 16 俺完全没看懂... 当我没说...
【在 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
发帖数: 4671 | 17 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
作业一律零分。
【在 i***r 的大作中提到】
: 谢谢大家，我确实是要搜SNP，chromosome都排好序了，sort过了，稍微优化了一下，
: 一会儿晚上看看跑了多少。
: 包子已发，谢谢
: PS 这里有没有new jersey附近做生物信息的？pm我大家认识下？
|
o**n
发帖数: 1249 | 18 请教我前面那个例子怎么避免二重循环或在matlab里怎么能做的更有效些？谢谢！
【在 l***y 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
|
G*****h
发帖数: 33134 | 19 matlab 也就在学校用用
不会也没事
【在 l***y 的大作中提到】
: 竟然用 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 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(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 的大作中提到】
: 我去查查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 的大作中提到】
: 有序序列，二分搜索好了
:
: 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 的大作中提到】
: 借机问一个问题。我有一个较大的二维矩阵，是(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 的大作中提到】
: 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 的大作中提到】
: 谢谢大家，我确实是要搜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 的大作中提到】
: 这里的人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 的大作中提到】
: 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 的大作中提到】
: 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 的大作中提到】
: 多谢思考我的问题和读我的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 的大作中提到】
: 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 的大作中提到】
: 【 以下文字转载自 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 的大作中提到】
: 竟然用 matlab 做大规模的多重循环？这都是在哪里学的 matlab？要是在我的班上，
: 作业一律零分。
|
