% Finding the feasible region and optimal solution of the (u, u^3) example.
thmin = -1.1;
thmax = 1.1;
% Need to sample with higher density for this boundary.
theta = linspace(thmin, thmax, 1001);
% Points on unit circle.
x = cos(pi()*theta);
y = sin(pi()*theta);
% Weird curvy branch for numerator = 0.
numcurves = NaN(3, length(theta));
for i = 1:length(theta)
r = roots([3, -3*x(i), -3*x(i)^2, -x(i)^3 + 4*y(i)]);
r(abs(imag(r)) > 1e-6) = NaN();
numcurves(:,i) = real(r);
end
% Nice branch where denominator = 0.
dencurve = -x;
% Need to do some magic to plot the curve properly because it's a multivalued
% function. Loop through the points, finding the next closest one at each step.
% This works because the curve is not self-intersecting and we know one of the
% endpoints. We bring the nice denominator curve along for the ride so that
% the region is easier to plot later on.
[inum, jnum] = find(~isnan(numcurves));
z = [numcurves(sub2ind(size(numcurves), inum, jnum))'; dencurve(jnum); ...
theta(jnum)];
zscale = max(z, [], 2) - min(z, [], 2);
z = bsxfun(@rdivide, z, zscale);
for i = 2:size(z, 2)
dist = [1, 0, 1]*bsxfun(@minus, z(:,i:end), z(:,i - 1)).^2;
[~, inext] = min(dist);
inext = inext + i - 1; % Correct for offset.
if inext ~= i
z(:,[i, inext]) = z(:,[inext, i]);
end
end
z = bsxfun(@times, z, zscale); % Rescale.
ufeas1 = max(z(1:2,:)); % Top piece.
ufeas2 = min(z(1:2,:)); % Bottom piece.
thfeas = z(3,:);
% Now find the minimizer. We use a coarser gridding here.
theta = linspace(thmin, thmax, 241);
urange = linspace(-2, 2, 1001);
[th, u] = meshgrid(theta, urange);
Vplus = cost(u, th, 1);
Vminus = cost(u, th, -1);
[V, firstsign] = min(cat(3, Vplus, Vminus), [], 3);
firstsign = 3 - 2*firstsign;
% Find local minima.
Vmid = V(2:end-1,:); % Slice out edges.
localmin = (Vmid < V(1:end-2,:)) & (Vmid < V(3:end,:));
localmin = [false(size(theta)); localmin; false(size(theta))]; % Add edges back.
ilocal = find(localmin(:));
Vlocal = V(ilocal);
thlocal = th(ilocal);
ulocal = u(ilocal);
% Find the absolute minima. Run fminsearch to refine a bit.
[Vmin, imin] = min(V);
umin = urange(imin);
for i = 1:length(umin)
[umin(i), Vmin(i)] = fminsearch(@(u) cost(u, theta(i), ...
firstsign(imin(i),i)), umin(i));
end
% Find lines of discontinuity.
[upieces, ujumps] = getdiscont(theta, umin, 2);
[Vpieces, Vjumps] = getdiscont(theta, Vmin, 2);
% Make a plot of u.
figure();
hold('on');
ufill = [ufeas1, ufeas2(end:-1:1)];
thfill = [thfeas, thfeas(end:-1:1)];
fill(thfill, ufill, 'r', 'edgecolor', 'r');
plot(thlocal, ulocal, 'sb', upieces(:,1), upieces(:,2), '-g', ...
ujumps(:,1), ujumps(:,2), ':g');
legend('Infeasible', 'Local Optima', 'Global Optima');
xlabel('\theta/\pi');
ylabel('Input');
% Make a plot of V.
figure();
hold('on');
semilogy(thlocal, Vlocal, 'sr', Vpieces(:,1), Vpieces(:,2), '-g', ...
Vjumps(:,1), Vjumps(:,2), ':g');
xlabel('\theta/\pi');
ylabel('Cost');
% Now save data.
feas = [thfeas', ufeas1', ufeas2'];
save('circle.dat', 'upieces', 'ujumps', 'Vpieces', 'Vjumps', 'feas');