1

I have a 3D density function q(x,y,z) that I am trying to plot in Matlab 8.3.0.532 (R2014a).

The domain of my function starts at a and ends at b, with uniform spacing ds. I want to plot the density on a ternary surface plot, where each dimension in the plot represents the proportion of x,y,z at a given point. For example, if I have a unit of density on the domain at q(1,1,1) and another unit of density on the domain at q(17,17,17), in both cases there is equal proportions of x,y,z and I will therefore have two units of density on my ternary surface plot at coordinates (1/3,1/3,1/3). I have code that works using ternsurf. The problem is that the number of proportion points grows exponentially fast with the size of the domain. At the moment I can only plot a domain of size 10 (in each dimension) with unit spacing (ds = 1). However, I need a much larger domain than this (size 100 in each dimension) and much smaller than unit spacing (ideally as small as 0.1) - this would lead to 100^3 * (1/0.1)^3 points on the grid, which Matlab just cannot handle. Does anyone have any ideas about how to somehow bin the density function by the 3D proportions to reduce the number of points?

My working code with example:

a = 0; % start of domain
b = 10; % end of domain
ds = 1; % spacing

[x, y, z] = ndgrid((a:ds:b)); % generate 3D independent variables
n = size(x);

q = zeros(n); % generate 3D dependent variable with some norm distributed density
for i = 1:n(1)
    for j = 1:n(2)
        for k = 1:n(2)
            q(i,j,k) = exp(-(((x(i,j,k) - 10)^2 + (y(i,j,k) - 10)^2 + (z(i,j,k) - 10)^2) / 20));
        end
    end
end

Total = x + y + z; % calculate the total of x,y,z at every point in the domain
x = x ./ Total; % find the proportion of x at every point in the domain
y = y ./ Total; % find the proportion of y at every point in the domain
z = z ./ Total; % find the proportion of z at every point in the domain
x(isnan(x)) = 0; % set coordinate (0,0,0) to 0
y(isnan(y)) = 0; % set coordinate (0,0,0) to 0
z(isnan(z)) = 0; % set coordinate (0,0,0) to 0

xP = reshape(x,[1, numel(x)]); % create a vector of the proportions of x
yP = reshape(y,[1, numel(y)]); % create a vector of the proportions of y
zP = reshape(z,[1, numel(z)]); % create a vector of the proportions of z

q = reshape(q,[1, numel(q)]);  % create a vector of the dependent variable q

ternsurf(xP, yP, q);  % plot the ternary surface of q against proportions
shading(gca, 'interp');
colorbar
view(2)
Improve this question

2 Answers 2

Reset to default
0

I believe you meant n(3) in your innermost loop. Here are a few tips:

1) Loose the loops:

q = exp(- ((x - 10).^2 + (y - 10).^2 + (z - 10).^2) / 20);

2) Loose the reshapes:

xP = x(:); yP = y(:); zP = z(:);

3) Check Total once, instead of doing three checks on x,y,z:

Total = x + y + z; % calculate the total of x,y,z at every point in the domain
Total( abs(Total) < eps ) = 1;
x = x ./ Total; % find the proportion of x at every point in the domain
y = y ./ Total; % find the proportion of y at every point in the domain
z = z ./ Total; % find the proportion of z at every point in the domain

PS: I just recognized your name.. it's Jonathan ;)

Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

Hey Jonathan, thanks, all good tips! Any ideas about subsequently binning my 'proportion space' to reduce the number of points? I think it's an interesting mathematical problem actually! I can think of how to do it with for loops and a lot of if else statements, but that will be very unwieldy and slow. Cheers!
I don't see anything about binning in the OP. If there is a subquestion you should include a clear description, possibly with a dummy code that doesn't run but looks like what you're trying to achieve :)
0

Discretization method probably depends on use of your plot, maybe it make sense to clarify your question from this point of view. Overall, you probably struggling with an "Out of memory" error, a couple of relevant tricks are described here http://www.mathworks.nl/help/matlab/matlab_prog/resolving-out-of-memory-errors.html?s_tid=doc_12b?refresh=true#brh72ex-52 . Of course, they work only up to certain size of arrays.

A more generic solution is too save parts of arrays on hard drive, it makes processing slower but it'll work. E.g., you can define several q functions with the scale-specific ngrids (e.g. ngridOrder0=[0:10:100], ngridOrder10=[1:1:9], ngridOrder11=[11:1:19], etc... ), and write an accessor function which will load/save the relevant grid and q function depending on the part of the plot you're looking.

Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.