How to form new meshgrid based on output condition of the vectorization of original meshgrid?

I apply a function to a meshgrid. If the output of the function f_nml2 is greater than say 20e3 I need to get rid of the value. Easy enough.

[N1,M1,L1]=meshgrid(0:N,0:M,0:L_1);

f_nml2=0.5*c_o*sqrt((N1./L).^2+(M1./W).^2+((2.*L1+1)./(2*H)).^2);

f_nml2=f_nml2(f_nml2<20e3)

However I required a new meshgrid N2,M2,L2 which has the corresponding indexes removed from N1,M1,L1 if the condition f_nml2>20e3 is met.

How to I remove combinations from the original meshgrid N1,M1,L1 which cause f_nml2 to be greater than 20e3?

Many thanks for you help

[Matt J]

You realize, don't you, that if you remove data from N1,M1,L1, leaving nothing behind to fill the hole, then there is no way for them to keep their original 3D shape?

[Henry Way]

Yes this is precisely my problem. The reason I am wanting to do this is that N1,M1,L1 have size 163x3726x95 which takes a lot of memory and computational time when manipulating.

I have come across your ndSparse class, if I made all values greater than 20e3 = 0 and then make it sparse, this may be a work around. Would this work, and if so please could you show me how to use your ndSparse to do this?

many thanks

[Henry Way]

I have tried the ndSparse(f_nml2) and it increase memory size from 461576880Bytes to 479652600Bytes.

[Matt J]

I have tried the ndSparse(f_nml2) and it increase memory size from 461576880Bytes to 479652600Bytes.

Probably because f_nml2 doesn't contain very many zeros.

What exactly is the motivation for keeping the 3D shape of N2,M2,L2? If you are replacing the unwanted elements with a non-index, it's hard to see what use it would be to you.

[Henry Way]

When making everything greater than 20e3 == 0, there are a total of 27718872 zeros in total.

My motivation is that the permutations of N,M and L are very important. This is acoustical modal theory, where I am calculating natural frequencies in a room (a water tank in this case). Which is dependent on a given permutation of N,M and L. 0,0,0;0,1,0;0,0,1 etc up to a maximum value of N=3725 M=162 L=94. Each permutation is known as a mode number.

[Matt J]

When making everything greater than 20e3 == 0, there are a total of 27718872 zeros in total.

That's not enough to benefit from sparse representation. You need well over 2/3 the total number of elements to be zero to benefit.

I think we need to see a typical computation that you would perform with N,M,L. I still don't see why their 3D shape is important. If the permutations are what matters, why not just keep a 3-column list of them?

[Henry Way]

%Calculate normalisation factors

epsilon1=N1;
epsilon2=M1;
epsilon3=L1;

epsilon1(epsilon1~=0)=2;
epsilon1(epsilon1==0)=1;

epsilon2(epsilon2~=0)=2;
epsilon2(epsilon2==0)=1;

epsilon3(epsilon3~=0)=2;
epsilon3(epsilon3==0)=1;

%Calculate pressure at source and microphone location

a=pi*x/L;
b=pi*y/W;
c=pi*z/(2*H);
a1=pi*xs/L;
b1=pi*ys/W;
c1=pi*zs/(2*H);
Psi_xyz=sqrt(epsilon1.*epsilon2.*epsilon3).*cos(N1.*(a)).*cos(M1.*(b).*cos((2.*L1+1).*(c)));
Psi_source=sqrt(epsilon1.*epsilon2.*epsilon3).*cos(N1.*(a1)).*cos(M1.*(b1)).*cos((2.*L1+1).*(c1));
zeta=c_o/2.*((epsilon1.*(wall_refl_coef(1,1)+wall_refl_coef(2,1))./L)+(epsilon2.*(wall_refl_coef(1,2)+wall_refl_coef(2,2))./W)+(epsilon3.*(wall_refl_coef(1,3)+wall_refl_coef(2,3))./H));

Psi_xyz1=reshape(Psi_xyz,[],1);
Psi_source1=reshape(Psi_source,[],1);
zeta=reshape(zeta,[],1);

%%Compute a_n for each mode. 

V=L*W*H;
omega_driving=driving_freq*2*pi;
res_freq=(2*pi).*res_freq;
q_source=0.00003;

a_n_numer=Psi_source1.*(omega_driving*rho*(c_o^2));%
a_n_denom=(V).*((2.*zeta.*res_freq.*omega_driving)+(1i*((omega_driving^2)-(res_freq.^2))));
a_n=q_source*a_n_numer./a_n_denom;

%%Calculate pressure at hydrophone location due to every mode contribution

Pressure=sum(a_n.*Psi_xyz1);

[Matt J]

Like I said, I see nothing in the calculations that utilizes the 3D nature of N1,M1,L1. In fact, you ultimately reshape() everything computed from them to column vectors anyway.

[Guillaume]

And in fact, keeping N1, M1, L1 together as a Nx3 matrix would also simplify all these equations. e.g.:

NML = [N(:), M(:), L(:)];
epsilon = (NML ~= 0) + 1;  %just one line instead of 6!
psi_xyz = sqrt(prod(epsilon, 2)) .* prod(cos((NML .* [1 1 2] + [0 0 1]) .* [a b c]), 2);

[Henry Way]

Your right. Many thanks for your help Matt and Guillaume.

Answer 1

[Matt J]

 keep=(f_nml2<20e3);

 N2=N1(keep);
 M2=M1(keep);
 L2=L1(keep);

[Henry Way]

Thank you. However I need N2,M2,L2 to remain a 3D matrix, so I would need to do a reshape and fill in the unwanted indices with something that will reduce memory and computational time.

Answer 2

[Matt J]

You can save about a factor of 6 in memory consumption by converting N1,M1,L1 to integer types

>>  M1=cast(M1,'uint8'); N1=cast(N1,'uint16'); L1=cast(L1,'uint8');

    Name   Size            Kilobytes     Class    Attributes

    L1     3727x164x96         57303     uint8              
    M1     3727x164x96         57303     uint8              
    N1     3727x164x96        114606     uint16

[Henry Way]

Thank you. However i do not think you can sqrt or cos an uint type. I would have to convert back to double before sqrt and cos, which brings me back to the same memory and performance problem.

[Matt J]

Well, you can achieve some savings by converting to 'single', if you can tolerate the reduced precision. That would be a factor of 2 on top of the factor of 2 you save by discarding elements, assuming you now agree that you don't need 3D arrays for N,M,L.

Answer 3

[Matt J]

None of the computations you've shown in fact require that you use meshgrid at all, assuming you have R2016b or higher. Just define N1, M1, L1 as vectors:

   N1=0:N;
   M1=(0:M).';
   L1=reshape(0:L,1,1,[]);

and all your computations are still valid. Even if you don't have R2016b, you can still achieve the same with bsxfun() albeit more tediously.