Split non-bijective function into two bijective regions

8 views (last 30 days)
I have a model which generates a large number of samples (i.e. a set of x values and a set of y values). The output roughly follows the shape of the curve pictured, so that there are two y values for each x value. The data is not sorted by x or y.
I would like to split this dataset into two bijective datasets, representing points above and below the marked X. (This will then allow me to interpolate between points on the upper portion of the curve, where otherwise interpolation would bounce up and down between the two parts of the curve.)
My initial approach was to find the co-ordinates of the marked X, xX and yX, and assign the two datasets according to y > yX and y < yX. However, there are a handful of points on the bottom portion of the curve that can sometimes have y > yX, so this doesn't work.
Can anyone suggest a potential solution to this problem?
  4 Comments
Dan
Dan on 3 Oct 2023
Thanks for the comments.
Dyuman: The marked point X has an x co-ordinate which is the maximum value of x (i.e. max(x))
Bruno: This is not my actual data but is a simple version that has the essence of the problem. My approach would put the points (1,2) and (2.5,1.8) into the top half of the curve
x = [4;2;3;6;6.5;4.5;1;1.5;3.5;2.5;5;6.2;2.8;6.3];
y = [5.5;6.5;1;4;1.5;1;2;7;6;1.8;5;1;6.2;3];
scatter(x,y)

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 3 Oct 2023
Edited: Bruno Luong on 3 Oct 2023
x = [4;2;3;6;6.5;4.5;1;1.5;3.5;2.5;5;6.2;2.8;6.3];
y = [5.5;6.5;1;4;1.5;1;2;7;6;1.8;5;1;6.2;3];
% TSP to order the data, See here
% https://www.mathworks.com/help/optim/ug/travelling-salesman-problem.html
nStops = length(x);
idxs = nchoosek(1:nStops,2);
% Bruno comment: you might normalize x and y before comuting dist if they have different units
dist = hypot(x(idxs(:,1)) - x(idxs(:,2)), ...
y(idxs(:,1)) - y(idxs(:,2)));
lendist = length(dist);
G = graph(idxs(:,1),idxs(:,2));
Aeq = spalloc(nStops,length(idxs),nStops*(nStops-1)); % Allocate a sparse matrix
for ii = 1:nStops
whichIdxs = (idxs == ii); % Find the trips that include stop ii
whichIdxs = sparse(sum(whichIdxs,2)); % Include trips where ii is at either end
Aeq(ii,:) = whichIdxs'; % Include in the constraint matrix
end
beq = 2*ones(nStops,1);
intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);
opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);
x_tsp = logical(round(x_tsp));
Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
%%
[~,minloc]=min(x);
[~,maxloc]=max(x);
p1 = Gsol.shortestpath(minloc,maxloc);
A=Gsol.adjacency;
A(p1(1),p1(2))=0;
A(p1(2),p1(1))=0;
Gtmp = graph(A);
p2 = Gtmp.shortestpath(minloc,maxloc);
p2(1) = [];
x1 = x(p1); y1 = y(p1);
x2 = x(p2); y2 = y(p2);
figure()
plot(x1,y1,'Linewidth',2)
hold on
plot(x2,y2,'Linewidth',2)
  4 Comments
Bruno Luong
Bruno Luong on 3 Oct 2023
The challenge is reorder the data, selecting the turning point is easy.
Torsten
Torsten on 3 Oct 2023
You are right (as always) - I didn't read the question carefully.

Sign in to comment.

More Answers (0)

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!