|
| 1 | +n = 1000; |
| 2 | + |
| 3 | +%initialize variables needed |
| 4 | +h = 0.1; |
| 5 | +e = 0.1; |
| 6 | + |
| 7 | +%Take the input value and create input amount of points |
| 8 | +points = rand(n,2); |
| 9 | +%create all the points within the unit circle |
| 10 | +theta = 0:0.01:2*pi; |
| 11 | +p = [cos(theta'),sin(theta')]; |
| 12 | +%Take only half the points in the heap and keep the points that are less |
| 13 | +%than 1 to create a simplex |
| 14 | +O = points(:,2) - 2*points(:,1) > 0; |
| 15 | +points(O,:) = []; |
| 16 | +O = points(:,2) + 2*points(:,1) > 2; |
| 17 | +points(O,:) = []; |
| 18 | +%draw the points |
| 19 | +hull = scatter(points(:,1),points(:,2),'yellow'); |
| 20 | +hold on; |
| 21 | + |
| 22 | +%Having the graph instantiated, we will dot p and all the points in the |
| 23 | +%simplex. Function takes 2 matrixes. |
| 24 | +maxU = u(p,points); |
| 25 | + |
| 26 | +%Next we want to determine the approximate half spaces of the simplex. This |
| 27 | +%function will take in 3 arguements: p, points, and maxU-h. |
| 28 | +arg3 = maxU-h; |
| 29 | +[result, index] = Nh(p,points,arg3) |
| 30 | + |
| 31 | +%Then we will sort the indexes in descending order |
| 32 | +[g,I] = sort(index, 'descend') |
| 33 | + |
| 34 | +N = zeros(length(points),1); |
| 35 | +ph = zeros(length(points),1); |
| 36 | +for k = 1:length(points) |
| 37 | + pcount = zeros(length(p),1); |
| 38 | + a = Inf; |
| 39 | + for i = 1:length(p) |
| 40 | + pxDot = p(i,1)*points(k,1)+p(i,2)*points(k,2); |
| 41 | + for j = 1:length(points) |
| 42 | + pyDot = p(i,1)*points(j,1)+p(i,2)*points(j,2); |
| 43 | + if pyDot >= pxDot - h |
| 44 | + pcount(i) = pcount(i) + 1; |
| 45 | + end |
| 46 | + end |
| 47 | + end |
| 48 | + [a, i] = min(pcount); |
| 49 | + N(k) = a; |
| 50 | + ph(k) = i; |
| 51 | +end |
| 52 | + |
| 53 | +[Nhs, J] = sort(N) |
| 54 | +J(1:20) |
| 55 | + |
| 56 | +dtheta = 0.5; |
| 57 | +thetaList = zeros(0); |
| 58 | +x = -0.1:0.01:1.1; |
| 59 | +j = 1; |
| 60 | +while length(thetaList) <= 2 |
| 61 | + y = (maxU(I(j)) - p(I(j),1)*x) / p(I(j),2); |
| 62 | + ans = theta(I(j)); |
| 63 | + a = 1; |
| 64 | + for k = 1:length(thetaList) |
| 65 | + if abs(ans-thetaList(k)) <= dtheta |
| 66 | + a = 0; |
| 67 | + end |
| 68 | + end |
| 69 | + if a == 1 |
| 70 | + pause; |
| 71 | + plot(x,y,'blue') |
| 72 | + thetaList(end+1) = ans; |
| 73 | + end |
| 74 | + j = j + 1; |
| 75 | +end |
| 76 | +pointsList = zeros(0); |
| 77 | + |
| 78 | +i = 1; |
| 79 | +interior = zeros(3,1); |
| 80 | +while length(pointsList) <= 2 |
| 81 | + ans1 = points(I(i)); |
| 82 | + a = 1; |
| 83 | + for j = 1:length(pointsList) |
| 84 | + if abs(ans1-pointsList(j)) <= dtheta |
| 85 | + a = 0; |
| 86 | + end |
| 87 | + end |
| 88 | + if a == 1 |
| 89 | + scatter(points(J(i),1),points(J(i),2),'red') |
| 90 | + pointsList(end+1) = ans1; |
| 91 | + interior(i,1) = points(J(i),1); |
| 92 | + interior(i,2) = points(J(i),2); |
| 93 | + end |
| 94 | + for z = 1:length(interior) %Trying to isolate 3 unique points. |
| 95 | + if points(z,1) == 0 |
| 96 | + points(z,1) = []; |
| 97 | + end |
| 98 | + if points(z,2) == 0 |
| 99 | + points(z,2) = []; |
| 100 | + end |
| 101 | + end |
| 102 | + i = i+1; |
| 103 | +end |
0 commit comments