%% Problems:  10.12.2017, spmd_minmaxproblems2.m
%    
%% (1) Max and min-points
%
% Find max and min-values and locations by decomposing domain
%  into labs and working in parallel with all labs available.
%
%% (2) Bisection ?, maybe not, n-section, yes!
%
%%  First example: 
clear
close all
format compact
%% Define  function f 
%
f=@(x) x.*sin(x);

%% Find local min and max on I=[-pi N*pi], take N=nlabs+1
% 
% In this case we know zeros of f: n*pi .
%% Open pool:
nlabs=2; % Nr. of workers on laptop
nlabs=20; % Triton max
I=[0, (nlabs+1)*pi];
ezplot(f,I);grid on;shg
hold on
%
% parpool(nlabs)   % Remove comment if no pool is open
%%  Some rules of pools: parpool(nlabs), spmd ... end, delete(gcp)
%
% parpool(4)   % This is just to show the "labintervals"
spmd 
    a=-pi+(labindex-1)*2*pi;
    b=pi+(labindex-1)*2*pi;
     fprintf('Subinterval: [%-4g, %-4g]\n',a,b)  % Remove comment for
    % nlabs=4 (eg.)
end
% Note1: If you started pool for nlabs=24 (eg.), you have to close it before
% changing to parpool(4) (eg.)  
% To close pool: 
%    delete(gcp)
%
% Note 2: If you have a pool open, you can do several spmd-blocks one after
% another. So the end stetement od yourcspmd-block doesn't close the pool.
%
%% Composite to ordinary variables (workers -> client)
A=[a{:}];B=[b{:}];
plot(A,f(A),'*k',B,f(B),'ok')
% Lab-boundaries a,b are zeros of f this time (zeros of sin).
%%  Find min and max on each lab.
n=1000;
disp([num2str(n*nlabs) ' calculations, each lab doing ' ... 
    num2str(n) ' simultaneously' ])
spmd
    x=linspace(a,b,n);
    y=f(x);
    % [ymax,maxI]=max(y); % Max on each lab and index of the 1'st
    % [ymin,minI]=min(y); % Min on each lab and its index of the 1'st 
    ymax=max(y); % Max values on each lab
    maxI=ymax==y; % Logical indexvector of (all) max-values in each lab
    minI=ymin==y; % Logical indexvector of (all) min-values in each lab
    xmax=x(maxI); % x-values of max's
    xmin=x(minI); % x-values of min's
end  
%% Error estimate
L=2*pi     % "Labinterval's" length
abserr=L/n % upper bound for abs error

%% Global max and min points
% gmaxx=[xmax{1} xmax{2}]
% gmaxy=[ymax{1} ymax{2}]
% gminx=[xmin{1} xmin{2}]
% gminy=[ymin{1} ymin{2}]
%% Composite to ordinary variables (results from workers to client):
Xmax=[xmax{:}];
Ymax=[ymax{:}]; 
Xmin=[xmin{:}]; 
Ymin=[ymin{:}]; 
%% Check lab-borders
sum(Xmax==A)
sum(Xmin==A)
% Both return 0, as they should, because the labintervals are [k-1,k+1]*pi.
plot(Xmin,Ymin,'*r',Xmax,Ymax,'*k')
%% Increase accuracy!  
%
%  Repeat the same procedure using labs with small radii around min/max-points.
nlabs=24; % Triton max  %as before
n=1000; % as before
% parpool(nlabs)
spmd
    am=xmin-abserr;bm=xmin+abserr;
    aM=xmax-abserr;bM=xmax+abserr;
    xm=linspace(am,bm,n);
    xM=linspace(aM,bM,n);
    ym=f(xm);yM=f(xM);
    [ymax,maxI]=max(yM); % Max on each lab and its index
    [ymin,minI]=min(ym); % Min on each lab and its index
    % In this case the shorter form using max/min is ok, as we know
    % uniqueness on the small intertevals. (But this is hardly more
    % efficient though.)
    xmax2=xM(maxI);       % x-values of max's
    ymax2=f(xmax2);
    xmin2=xm(minI);       % x-values of min's
    ymin2=f(xmin2);
end  
Xmax2=[xmax2{:}];
Ymax2=[ymax2{:}]; 
Xmin2=[xmin2{:}]; 
Ymin2=[ymin2{:}]; 
%% Error:
abserr2=abserr/n
max(abs(Xmax2-Xmax))
max(abs(Ymax2-Ymax))
%
%% In some cases need to remove extra points
% Here's one attempt:
%{
aa=[a{:}];bb=[b{:}];
Axmax=[aa' gmaxx']
pick=~(Axmax(:,1)==Axmax(:,2));
xmaxtrue=gmaxx(pick);
ymaxtrue=gmaxy(pick);
%}
%% Removal of a single value is easy:
t=-2:2; T=[t t t t]
% Remove 0's
T(T~=0)
% or
T(T==0)=[]
%% Zeros of f using brackets we just found
Br=[Ymax' Ymin']
spmd
    ab=[Xmax(labindex) Xmin(labindex)];
    z=fzero(f,ab);
end        
Z=[z{:}]
% Here of course we know the result.
Z-(1:nlabs)*pi
% YES !!!
%% 








