- #1
hoffmann
- 70
- 0
Hi everyone. I'm working on a piece of code and I'm having a bit of trouble with array pre-allocation. Here it is:
function [react_in1 react_in2 react_out1 react_out2]=oneVtwoXone(masses,prec);
prec = 1e-4;
N=length(masses);
masses=[0;masses]; % adding dummy metabolite with zero mass
masses=masses(:);
%mass2=masses*ones(size(masses))'+ones(size(masses))*masses';
%mass2=mass2(:);
mass2=zeros(N*(N+1)/2,1); % list of mass pairs
met1=mass2; % list of the first metabolite in a pair
met2=mass2; % same for the second metabolite
k=0;
for i=1:N
for j=i:N
k=k+1;
met1(k)=i;
met2(k)=j;
mass2(k)=masses(i)+masses(j);
end
end
disp(['Number of distinct pairs is' length(mass2) ])
[mass2, ord]=sort(mass2); % sorted mass pairs
met1=met1(ord); % metabolites in the same order as sorted mass pairs
met2=met2(ord);
dm=diff(mass2);
while(length(find(dm <= mass2(1:length(dm))*prec)) > 0) % while have metabolite pairs that have about the same mass
ind=find(dm<mass2(length(mass2)-length(dm)+1:end)*prec);
react_in1=[react_in1 met1(ind)];
react_in2=[react_in2 met2(ind)];
react_out1=[react_out1 met1(ind+1)];
react_out2=[react_out2 met2(ind+1)];
dm = dm(1:end-1) + dm(2:end);
end
disp(['Number of conforming reactions is' length(react_in1) ])
The source of my trouble is in the second while loop. The react_in1 etc. takes up too much memory and I need to figure out how to preallocate a matrix to hold values approximately 10000, and how to delete the excess zeros if there is any unused space.
Can anyone help?
function [react_in1 react_in2 react_out1 react_out2]=oneVtwoXone(masses,prec);
prec = 1e-4;
N=length(masses);
masses=[0;masses]; % adding dummy metabolite with zero mass
masses=masses(:);
%mass2=masses*ones(size(masses))'+ones(size(masses))*masses';
%mass2=mass2(:);
mass2=zeros(N*(N+1)/2,1); % list of mass pairs
met1=mass2; % list of the first metabolite in a pair
met2=mass2; % same for the second metabolite
k=0;
for i=1:N
for j=i:N
k=k+1;
met1(k)=i;
met2(k)=j;
mass2(k)=masses(i)+masses(j);
end
end
disp(['Number of distinct pairs is' length(mass2) ])
[mass2, ord]=sort(mass2); % sorted mass pairs
met1=met1(ord); % metabolites in the same order as sorted mass pairs
met2=met2(ord);
dm=diff(mass2);
while(length(find(dm <= mass2(1:length(dm))*prec)) > 0) % while have metabolite pairs that have about the same mass
ind=find(dm<mass2(length(mass2)-length(dm)+1:end)*prec);
react_in1=[react_in1 met1(ind)];
react_in2=[react_in2 met2(ind)];
react_out1=[react_out1 met1(ind+1)];
react_out2=[react_out2 met2(ind+1)];
dm = dm(1:end-1) + dm(2:end);
end
disp(['Number of conforming reactions is' length(react_in1) ])
The source of my trouble is in the second while loop. The react_in1 etc. takes up too much memory and I need to figure out how to preallocate a matrix to hold values approximately 10000, and how to delete the excess zeros if there is any unused space.
Can anyone help?