function varargout = RootFinding(varargin)
% ROOTFINDING -- Demo of several root-finding algorithms.
%       Simply type 'RootFinding' to run.
%

% Last Modified by GUIDE v2.5 21-Aug-2015 14:38:38

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @RootFinding_OpeningFcn, ...
                   'gui_OutputFcn',  @RootFinding_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

%#ok<*INUSL>
%#ok<*INUSD>
%#ok<*DEFNU>

% --- Executes just before RootFinding is made visible.
function RootFinding_OpeningFcn(hObject, eventdata, handles, varargin)  %#ok<*INUSL>
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to RootFinding (see VARARGIN)

global EvalError

EvalError = [];

% Create x- and y-axes for figure.
axes(handles.FunctionAxes);
hold all;
xAxis = line([-10000 10000],[0 0]);
yAxis = line([0 0], [-10000 10000]);
set([xAxis yAxis], 'Color', [0.6 0.6 0.6], 'LineWidth', 1.5);

handles.xAxis = xAxis;
handles.yAxis = yAxis;

% This variable will store the handle to the plotted function.
handles.fPlot = [];

% Style structure for plotted function
handles.fPlotStyle.LineWidth = 2;
handles.fPlotStyle.Color = [0 0.39 0.66];

% Handle to error text (if any)
handles.fErrorText = [];

% Style structure for error text.
handles.fErrorTextStyle.HorizontalAlignment = 'center';
handles.fErrorTextStyle.FontSize = 20.0;
handles.fErrorTextStyle.Color = [1 0 0];

% Create root region shading rectangle.
handles.fRootRegionRect = TransparentRect([0 1], [0 1], 'y', 0.5);
set(handles.fRootRegionRect, 'ZData', [1 1 1 1]);

% Create root region
handles.fRootRegion = plot3(0,0,2);
%set(handles.fRootRegion, 'XData', [], 'YData', []);
set(handles.fRootRegion, 'ButtonDownFcn', @RRMouseDown);

% Styles for root region.
handles.fRootRegionBaseStyle.MarkerSize = 12;
handles.fRootRegionBaseStyle.MarkerFaceColor = [0.54 0.76 0.97];
handles.fRootRegionBaseStyle.MarkerEdgeColor = [0.17 0.25 0.79];
handles.fRootRegionBaseStyle.Marker = 'o';
handles.fRootRegionBaseStyle.LineWidth = 2;
handles.fRootRegionShowLineStyle.LineStyle = '-';
handles.fRootRegionHideLineStyle.LineStyle = 'none';

set(handles.fRootRegion, handles.fRootRegionBaseStyle);

% Handle to array of the plotted error lineseries
handles.ePlot = [];

% Style structure for plotted error
handles.ePlotStyle(1).LineWidth = 2;
handles.ePlotStyle(2).LineWidth = 2;
handles.ePlotStyle(3).LineWidth = 2;
handles.ePlotStyle(4).LineWidth = 2;
handles.ePlotStyle(1).Color = [0.23 0.44 0.34];
handles.ePlotStyle(2).Color = [0.51 0.38 0.48];
handles.ePlotStyle(3).Color = [0.39 0.47 0.64];
handles.ePlotStyle(4).Color = [0.87 0.49 0.00];
handles.ePlotStyle(1).Marker = 'x';
handles.ePlotStyle(2).Marker = 's';
handles.ePlotStyle(3).Marker = 'd';
handles.ePlotStyle(4).Marker = 'o';


% Variable which stores what the user is currently dragging ('', 'a', or
% 'b')
handles.Dragging = '';

% Flag for new user input (used to halt 'Run' mode)
handles.GotUI = false;

% Collection of checkboxes for errors.
handles.AllErrorCBs = [handles.ErrBisectionCB
                       handles.ErrRegulaFalsiCB
                       handles.ErrSecantCB
                       handles.ErrNewtonCB];

% Number of iterations to use
handles.Iterations = 30;

% Choose default command line output for RootFinding
handles.output = hObject;

handles.MethodIndex = 1;
handles.Method = '';
handles.Function = '';
handles.History = [];
handles.Roots = [];

set(handles.ErrorLabelAxes, 'Color', 'none');
axes(handles.ErrorLabelAxes);
t = text(0,0,'Error (in x)');
ErrorLabelStyle.FontSize = 14;
ErrorLabelStyle.Rotation = 90;
ErrorLabelStyle.Color = [1 1 1];
ErrorLabelStyle.HorizontalAlignment = 'center';
set(t, ErrorLabelStyle);

%set(get(handles.ErrorAxes, 'XLabel'), 'String', 'Iteration');
%set(get(handles.ErrorAxes, 'YLabel'), 'String', 'Error in x');

% Update handles structure
guidata(hObject, handles);

UpdateFunction;
RedrawFunction;
UpdateMethod;
ResetStep;




% UIWAIT makes RootFinding wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = RootFinding_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes during object deletion, before destroying properties.
function RootFinding_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

GotUI;

% --- Executes on mouse motion over figure - except title and menu.
function RootFinding_WindowButtonMotionFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

try
    if ~isempty(handles.Dragging)
        RRDragUpdate;
    end
catch me %#ok<NASGU>
    %disp(me.message);
    % Apparently this function can be called even without user input
    % when window is created, in which case we don't have a 'Dragging'
    % field yet.
end

% --- Executes on mouse press over figure background, over a disabled or
% --- inactive control, or over an axes background.
function RootFinding_WindowButtonUpFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isempty(handles.Dragging)
    RRDragDone;
end

% --- Executes on selection change in MethodPopup.
function MethodPopup_Callback(hObject, eventdata, handles)
% hObject    handle to MethodPopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

GotUI;
UpdateMethod;

% --- Executes during object creation, after setting all properties.
function MethodPopup_CreateFcn(hObject, eventdata, handles)
% hObject    handle to MethodPopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called


% --- Executes on selection change in FunctionPopup.
function FunctionPopup_Callback(hObject, eventdata, handles) 
% hObject    handle to FunctionPopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

GotUI;
UpdateFunction;
RedrawFunction;

% --- Executes during object creation, after setting all properties.
function FunctionPopup_CreateFcn(hObject, eventdata, handles)
% hObject    handle to FunctionPopup (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called


function CustomFunction_Callback(hObject, eventdata, handles)
% hObject    handle to CustomFunction (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of CustomFunction as text
%        str2double(get(hObject,'String')) returns contents of CustomFunction as a double
GotUI;
UpdateFunction;
RedrawFunction;


% --- Executes during object creation, after setting all properties.
function CustomFunction_CreateFcn(hObject, eventdata, handles)
% hObject    handle to CustomFunction (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --- Executes on key press with focus on CustomFunction and none of its controls.
function CustomFunction_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to CustomFunction (see GCBO)
% eventdata  structure with the following fields (see UICONTROL)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)


% --- Executes on button press in RunButton.
function RunButton_Callback(hObject, eventdata, handles)
% hObject    handle to RunButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

GotUI;
if strcmp(get(hObject, 'String'), 'Run');
    set(hObject, 'String', 'Stop');
    Run(1.5);
end
set(hObject, 'String', 'Run');

% --- Executes on button press in StepButton.
function StepButton_Callback(hObject, eventdata, handles)
% hObject    handle to StepButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
GotUI;
Step;



% --- Executes on button press in ErrBisectionCB.
function ErrBisectionCB_Callback(hObject, eventdata, handles)
% hObject    handle to ErrBisectionCB (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
UpdateError;


% --- Executes on button press in ErrSecantCB.
function ErrSecantCB_Callback(hObject, eventdata, handles)
% hObject    handle to ErrSecantCB (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
UpdateError;



% --- Executes on button press in ErrSecantCB.
function ErrRegulaFalsiCB_Callback(hObject, eventdata, handles)
% hObject    handle to ErrSecantCB (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
UpdateError;



% --- Executes on button press in ErrNewtonCB.
function ErrNewtonCB_Callback(hObject, eventdata, handles)
% hObject    handle to ErrNewtonCB (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
UpdateError;




function UpdateError

handles = guidata(gcf);

axes(handles.ErrorAxes);
hold all;

delete(handles.ePlot); handles.ePlot = [];

show = get(handles.AllErrorCBs, 'Value');
show = logical([show{:}]);
ePlot = nan(length(show), 1);

for i = 1:length(show)
    if show(i)
        ePlot(i) = semilogy(log(handles.Errors(i,:)));
        set(ePlot(i), handles.ePlotStyle(i));
    end
end
ePlot(isnan(ePlot)) = [];

handles.ePlot = ePlot;

guidata(gcf, handles);



function GotUI

handles = guidata(gcf);
handles.GotUI = true;
guidata(gcf, handles);

function UpdateMethod

handles = guidata(gcf);

mpop = handles.MethodPopup;
contents = cellstr(get(mpop, 'String'));
midx = get(mpop, 'Value');
mname = contents{midx};

handles.MethodIndex = midx;
handles.Method = mname;

% In case the 'root not bracketed' error is shown
handles = ClearErrorText(handles);

guidata(gcf, handles);
ResetStep;

%DefaultRR;


function UpdateFunction

global EvalError

EvalError = [];

handles = guidata(gcf);

fpop = handles.FunctionPopup;
contents = cellstr(get(fpop, 'String'));
fname = contents{get(fpop, 'Value')};
fname(isspace(fname)) = [];                     % Delete spaces from name

custom = strcmp(fname, 'Custom');
if custom,
    vis = 'on';
else
    vis = 'off';
end
set(handles.CustomFunction, 'Visible', vis);

handles.Function = fname;

if custom,
    uicontrol(handles.CustomFunction);
    try
        handles.Function = eval(['@(x) ' get(handles.CustomFunction, 'String')]);
    catch me %#ok<NASGU>
        handles.Function = @Invalid;
    end
end

guidata(gcf, handles);

DefaultRR;


function DefaultRR

handles = guidata(gcf);

a = -4; b = 4;
if ischar(handles.Function),
    switch handles.Function
        case 'BesselJ0'
            a = 0; b = 4;
        case 'SmoothedStep'
            a = -1.3; b = 2;
        case 'Kepler'
            a = -1.5; b = 5;
        case 'DoubleRoot'
            a = -2; b = 3;
    end
end

handles.nextA = a;
handles.nextB = b;
handles.nextO = (a+b)/2;
fa = Eval(a);
fb = Eval(b);

guidata(gcf, handles);

RRControl(true,false,[a b],[fa fb]);
RRUpdated;

function RRUpdated

RunMethods;
UpdateError;
ResetStep;

function RRMouseDown(rr, eventData)

GotUI;

handles = guidata(gcbo);

% In case the 'root not bracketed' error is shown
handles = ClearErrorText(handles);

cp = get(handles.FunctionAxes, 'CurrentPoint');

if strcmp(handles.Method, 'Newton')
    handles.Dragging = 'o';
else
    x = cp(1,1);
    if abs(x - handles.nextA) < abs(x - handles.nextB),
        handles.Dragging = 'a';
    else
        handles.Dragging = 'b';
    end
end

guidata(gcbo, handles);

function RRDragUpdate

handles = guidata(gcbo);
cp = get(handles.FunctionAxes, 'CurrentPoint');
x = cp(1,1);

switch handles.Dragging
    case 'a'
        handles.nextA = x;
    case 'b'
        handles.nextB = x;
    case 'o'            % Newton's method: single point.
        handles.nextO = x;
end

if strcmp(handles.Dragging, 'o')
    h = 1e-5;
    xd = [x-h x+h];
    doLine = true;
    doRect = false;
else
    xd = [handles.nextA handles.nextB];
    doLine = false;
    doRect = true;
end

yd = Eval(xd);

RRControl(doRect, doLine, xd, yd);

guidata(gcbo, handles);

function RRDragDone

% Update final location
RRDragUpdate;

% Stop dragging
handles = guidata(gcbo);
handles.Dragging = '';
guidata(gcbo, handles);

% Update the RR.
RRUpdated;

function RRRedraw

handles = guidata(gcf);

if isempty(handles.History), return; end

it = handles.StepIteration;
isMidStep = it > floor(it);

midx = handles.MethodIndex;
method = handles.Method;
hist = handles.History(midx,floor(it),:);

if isnan(hist(1)),
    handles = MakeErrorText(handles, 'Root is not bracketed. ');
    a = handles.nextA;
    b = handles.nextB;
    fa = Eval(a);
    fb = Eval(b);
    RRControl(true, false, [a b], [fa fb]);
    guidata(gcf, handles);
    return;
end

% In case the 'root not bracketed' error is shown
handles = ClearErrorText(handles);
guidata(gcf,handles);

if strcmp(method, 'Newton')
    RRControl(false, true, hist([1 3]), [hist(4) 0]);
elseif isMidStep,
    if strcmp(method, 'Bisection'),
        RRControl(true, false, hist([1 3 2]), hist([4 6 5]));
    else
        RRControl(true, true, hist([1 3 2]), [hist(4) 0 hist(5)]);            
    end
else
    RRControl(true, false, hist(1:2), hist(4:5));
end



function RedrawFunction

global EvalError

handles = guidata(gcf);

% Break if nothing to draw
if isempty(handles.Function), return; end

% Delete old function
delete(handles.fPlot); handles.fPlot = [];
handles = ClearErrorText(handles);

axes(handles.FunctionAxes); hold on;
%handles.fPlot = ezplot(handles.Function, xlim());
xl = xlim();
xs = linspace(xl(1), xl(2), 500);

% Evaluate function.
ys = Eval(xs);

if ~isempty(EvalError)
    handles = MakeErrorText(handles,{'Error evaluating function:'; EvalError});
end

% Plot
d = max(max(ys) - min(ys), 1e-8);
ylim([min(ys) - 0.3*d, max(ys) + 0.3*d]);

h = plot(xs,ys);
set(h, handles.fPlotStyle);
handles.fPlot = h;

guidata(gcf, handles);


function handles = MakeErrorText(handles, errmsg)

axes(handles.FunctionAxes);
t = text(mean(xlim),mean(ylim),100,errmsg);
set(t, handles.fErrorTextStyle);
handles.fErrorText = t;


function handles = ClearErrorText(handles)

delete(handles.fErrorText); handles.fErrorText = [];


% Create a semi-transparent rectangle.
%   x = [xmin xmax]
%   y = [ymin ymax]
%   c: Color specification ([r g b] or single-letter spec: 'y', 'k', ...)
%   a: Alpha value (0 <= a <= 1)
function [xc,yc] = RectXYToVertices(x, y)
xc = [x(1) x(2) x(2) x(1)];
yc = [y(1) y(1) y(2) y(2)];

function h = TransparentRect(x, y, c, a)
[xc,yc] = RectXYToVertices(x,y);
h = fill(xc,yc,c);
alpha(h, a);
set(h, 'EdgeColor', 'none');

function MoveRect(h, x, y)
[xc,yc] = RectXYToVertices(x,y);
set(h, 'XData', xc, 'YData', yc);

% Root region display
function RRControl(show, doLine, x, y)
handles = guidata(gcf);
rr = handles.fRootRegion;
rrRect = handles.fRootRegionRect;

    
if show, vis = 'on'; else vis = 'off'; end

set(rrRect, 'Visible', vis);

if doLine,
    set(rr, handles.fRootRegionShowLineStyle);
    
    % Extend line.
    xl = get(handles.FunctionAxes, 'XLim');
    x0 = xl(1) - 1;
    x1 = xl(2) + 1;
    m = (y(end)-y(1))/(x(end)-x(1));
    y0 = y(1) + (x0 - x(1)) * m;
    y1 = y(1) + (x1 - x(1)) * m;
    ex = [x0 x(:).' x1];
    ey = [y0 y(:).' y1];
else
    set(rr, handles.fRootRegionHideLineStyle);
    ex = x;
    ey = y;
end

% Put root region at z = 2.
set(rr, 'XData', ex, 'YData', ey, 'ZData', ones(size(ex))*2);

axes(handles.FunctionAxes);
MoveRect(rrRect, [min(x) max(x)], [-1000000 1000000]);




function y = Eval(x)

global EvalError

handles = guidata(gcf);

try
    y = feval(handles.Function, x);
catch me %#ok<NASGU>
    % OK, not vectorized, do it one at a time
    try
        y = zeros(size(x));
        for i = 1:length(x)
            y(i) = feval(handles.Function, x(i));
        end
    catch me
        % Error evaluating function.
        EvalError = [' ' me.message ' '];
        y = nan(size(x));
    end
end


%
% Root-finding algorithms
%

function history = SecantMethod(a,b,it)

fa = Eval(a);
fb = Eval(b);
history = nan(it,6);

for i = 1:it
    c = (a*fb - b*fa) / (fb - fa);
    fc = Eval(c);

    if a < b,
        history(i,:) = [a b c fa fb fc];
    else
        history(i,:) = [b a c fb fa fc];        
    end

    % Drop a and move c->b->a
    a = b; fa = fb;
    b = c; fb = fc;    
end

function history = RegulaFalsi(a,b,it)

fa = Eval(a);
fb = Eval(b);
history = nan(it,6);

for i = 1:it
    c = (a*fb - b*fa) / (fb - fa);
    fc = Eval(c);

    history(i,:) = [a b c fa fb fc];

    if fb*fc >= 0,
        b = c; fb = fc;
    else
        a = c; fa = fc;
    end
end

function history = BisectionMethod(a,b,it)

fa = Eval(a);
fb = Eval(b);
history = nan(it,6);

for i = 1:it
    c = (a+b)/2;
    fc = Eval(c);

    history(i,:) = [a b c fa fb fc];

    if fb*fc >= 0,
        b = c; fb = fc;
    else
        a = c; fa = fc;
    end
end

function [d,fxmh,fxph] = ApproxDerivative(x)

h = 1e-5;
fxph = Eval(x+h);
fxmh = Eval(x-h);
d = (fxph - fxmh) / (2*h);


function history = NewtonMethod(o,it)

history = nan(it,6);
x = o;
fx = Eval(x);
h = 1e-5;

for i = 1:it
    [fpx,fxmh,fxph] = ApproxDerivative(x);
    xn = x - fx/fpx;
    fxn = Eval(xn);
    history(i,:) = [x-h x+h xn fxmh fxph fxn];
    x = xn;
    fx = fxn;
end

function RunMethods

handles = guidata(gcf);

a = handles.nextA;
b = handles.nextB;
o = handles.nextO;
it = handles.Iterations;
fa = Eval(a);
fb = Eval(b);

history = nan(4,it,6);

if (fa*fb) <= 0,
    history(1,:,:) = BisectionMethod(a,b,it);
    history(2,:,:) = RegulaFalsi(a,b,it);
end
history(3,:,:) = SecantMethod(a,b,it);
history(4,:,:) = NewtonMethod(o,it);

handles.History = history;
guidata(gcf, handles);

r = GetRoots;

errs = bsxfun(@minus, history(:,:,3), reshape(r, 1, 1, []));
errs = min(abs(errs),[],3);

handles.Errors = errs;
guidata(gcf, handles);


% Find "all" relevant roots using MATLAB's fzero.
function r = GetRoots

handles = guidata(gcf);

% Get a bunch of initial guesses for roots.
a = handles.nextA;
b = handles.nextB;
x0s = linspace(a,b,11);                         % Evenly spaced points in the region.
x0s = [x0s handles.History(:,end,3).'];           % Also include the final results from each root-finding algorithm.
x0s(~isfinite(x0s)) = [];                       % Remove any nonfinites

r = zeros(size(x0s));
for i = 1:numel(x0s)
    r(i) = fzero(@Eval, x0s(i));
end
r(~isfinite(r)) = [];                           % Remove any failed roots.
r = unique(r);







function it = Step

handles = guidata(gcf);

it = handles.StepIteration;

% Abort if we are at the end.
if it >= handles.Iterations + 1, return; end

if strcmp(handles.Method, 'Newton')
    inc = 1;
else
    inc = 0.5;
end

handles.StepIteration = it + inc;
guidata(gcf,handles);

RRRedraw;

function ResetStep

handles = guidata(gcf);
handles.StepIteration = 1;
guidata(gcf,handles);

RRRedraw;



function Run(delay)

handles = guidata(gcf);
isStarted = handles.StepIteration > 0;
ResetStep;
RRRedraw;
handles = guidata(gcf);
handles.GotUI = false;
guidata(gcf, handles);

i = 0;

if isStarted, pause(delay); end

% Test for isstruct(handles) because on figure deletion, handles becomes
%  a non-structure.
while isstruct(handles) && i <= handles.Iterations && ~handles.GotUI,
    i = Step();
    pause(delay);
    handles = guidata(gcf);
end


%
% Built-in functions for UI
%

function y = BesselJ0(x)
y = besselj(0,x);

function y = SmoothedStep(x)
y = (0.5*(sign(x)+1) .* atan(x.^8)/(pi/2)) - 0.5;

function y = Kepler(x)
M = pi/2;
e = 0.5;
y = M-(x - e*sin(x));

function y = DoubleRoot(x)
y = 0.9 - exp(-(x-1.43).^2);

function y = OffcenteredCubic(x)
y = (x-0.3).^3;

function Invalid(~)
error('Could not parse input');