text:
function [X, fX, i] = minimize_lbfgsb(X, f, length, varargin)

% Minimize a differentiable multivariate function using quasi Newton.
%
% Usage: [X, fX, i] = minimize_lbfgsb(X, f, length, P1, P2, P3, ... )
%
% X       initial guess; may be of any type, including struct and cell array
% f       the name or pointer to the function to be minimized. The function
%         f must return two arguments, the value of the function, and it's
%         partial derivatives wrt the elements of X. The partial derivative
%         must have the same type as X.
% length  length of the run; if it is positive, it gives the maximum number of
%         line searches, if negative its absolute gives the maximum allowed
%         number of function evaluations. Optionally, length can have a second
%         component, which will indicate the reduction in function value to be
%         expected in the first line-search (defaults to 1.0).
% P1, P2  ... parameters are passed to the function f.
%
% X       the returned solution
% fX      vector of function values indicating progress made
% i       number of iterations (line searches or function evaluations,
%         depending on the sign of "length") used at termination.
%
% The function returns when either its length is up, or if no further progress
% can be made (ie, we are at a (local) minimum, or so close that due to
% numerical problems, we cannot get any closer). NOTE: If the function
% terminates within a few iterations, it could be an indication that the
% function values and derivatives are not consistent (ie, there may be a bug in
% the implementation of your "f" function).
%
% Copyright (C) 2010 by Hannes Nickisch, 2010-02-05

% global variables serve as communication interface between calls
global minimize_lbfgsb_iteration_number
global minimize_lbfgsb_objective
global minimize_lbfgsb_gradient
global minimize_lbfgsb_X

% init global variables
minimize_lbfgsb_iteration_number = 0;
minimize_lbfgsb_objective = Inf;
minimize_lbfgsb_gradient = 0*unwrap2vec(X);
minimize_lbfgsb_X = X;

X0 = X;
lb = -Inf*ones(size(unwrap2vec(X0)));
ub =  Inf*ones(size(unwrap2vec(X0)));
maxiter = abs(length); % max number of iterations

% no callback routine used so far
% m is the number of saved vectors used to estimate the Hessian
% factr is the precision 1e-12
X = lbfgsb( unwrap2vec(X0), lb, ub, 'minimize_lbfgsb_objfun', ...
                                'minimize_lbfgsb_gradfun', ...
         {f,X0,varargin{:}}, [], ...
         'maxiter',maxiter, 'm',4, 'factr',1e-12, 'pgtol',1e-5);
i = minimize_lbfgsb_iteration_number;
fX = minimize_lbfgsb_objective;
X = rewrap(X0,X);

% clear global variables
clear minimize_lbfgsb_iteration_number
clear minimize_lbfgsb_objective
clear minimize_lbfgsb_gradient
clear minimize_lbfgsb_X

% Extract the numerical values from "s" into the column vector "v". The
% variable "s" can be of any type, including struct and cell array.
% Non-numerical elements are ignored. See also the reverse rewrap.m.
function v = unwrap2vec(s)
  v = [];
  if isnumeric(s)
    v = s(:);                       % numeric values are recast to column vector
  elseif isstruct(s)
    v = unwrap2vec(struct2cell(orderfields(s)));% alphabetize, conv to cell, recurse
  elseif iscell(s)
    for i = 1:numel(s)            % cell array elements are handled sequentially
      v = [v; unwrap2vec(s{i})];
    end
  end                                                  % other types are ignored

% Map the numerical elements in the vector "v" onto the variables "s" which can
% be of any type. The number of numerical elements must match; on exit "v"
% should be empty. Non-numerical entries are just copied. See also unwrap2vec.m.
function [s v] = rewrap(s, v)
  if isnumeric(s)
    if numel(v) < numel(s)
      error('The vector for conversion contains too few elements')
    end
    s = reshape(v(1:numel(s)), size(s));           % numeric values are reshaped
    v = v(numel(s)+1:end);                       % remaining arguments passed on
  elseif isstruct(s)
    [s p] = orderfields(s); p(p) = 1:numel(p);     % alphabetize, store ordering
    [t v] = rewrap(struct2cell(s), v);                % convert to cell, recurse
    s = orderfields(cell2struct(t,fieldnames(s),1),p); % conv to struct, reorder
  elseif iscell(s)
    for i = 1:numel(s)            % cell array elements are handled sequentially
      [s{i} v] = rewrap(s{i}, v);
    end
  end                                            % other types are not processed

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [l, L_rf, L_sf, L_k, L_seg, L_n] = ...
    retroProjHmgLinFromPinHoleOnRob(Rf, Sf, k, seg, n)

% RETROPROJHMGLINFROMPINHOLEONROB retroprj Hmg Line from pinhole on robot.

%   Copyright 2009 Teresa Vidal.

if nargout == 1

    ls = invPinHoleHmgLin(k, seg, n) ;
    lr = fromFrameHmgLin(Sf, ls);
    l  = fromFrameHmgLin(Rf, lr);

else % Jacobians requested

    [ls, LS_seg, LS_n, LS_k] = invPinHoleHmgLin(seg, n, k) ;
    [lr, LR_sf, LR_ls]       = fromFrameHmgLin(Sf, ls);
    [l, L_rf, L_lr]          = fromFrameHmgLin(Rf, lr);

    L_sf  = L_lr*LR_sf;
    L_ls  = L_lr*LR_ls;
    L_k   = L_ls*LS_k;
    L_seg = L_ls*LS_seg;
    L_n   = L_ls*LS_n;

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [mapStruct] = modifyReactionsMetabolites(map, rxnList, metList, newColor, newAreaWidth)
% Modifies the color and areaWidth of reactions from a given list as input
% and the color of the corresponding metabolites from a given list as
% input. The colors and areaWidth are given as inputs and only metabolites
% present in the given reactions list will be colored.
%
% USAGE:
%
%   [mapStruct] = modifyReactionsMetabolites(map, rxnList, metList, newColor, newAreaWidth)
%
% INPUTS:
%   map:            Matlab structure of the map obtained from the
%                   function "transformXML2Map".
%   rxnList:        List of reaction names as a string array
%   metList:        List of metabolite names as a string array
%   newColor:       Color chosen for reaction lines and metabolites
%                   given as a string with the corresponding real name.
%                   Possible names can be found in the function
%                   "createColorsMap.m".
%   newAreaWidth:   Width size for the reaction lines. Can be given as
%                   a string or a double.
%
% OUTPUT:
%   mapStruct:      Updated map structure with the changed areaWidth and
%                   color of the reactions and their corresponding
%                   metabolites.
%
% .. Author: - N.Sompairac - Institut Curie, Paris, 25/07/2017

    colors = createColorsMap;  % Create a Color map with corresponding colors names and their HTML code

    mapStruct = map;

    % Get the indexes of the needed reactions to color
    rxnIndexList = find(ismember(mapStruct.rxnName, rxnList));

    % Initialise a list that will contain aliases of molecules implicated in
    % all the needed reactions
    rxnAliasList = {};

    % Loop over reactions and change the corresponding values
    for rxn = rxnIndexList'

        % Change de color of the reaction
        mapStruct.rxnColor{rxn} = colors(newColor);

        % Change de areaWidth of the reaction
        mapStruct.rxnWidth{rxn} = newAreaWidth;

        % Get the list of aliases involved in this reaction
        % Loop over base reactants
        for x = 1:length(mapStruct.rxnBaseReactantID{rxn})
            rxnAliasList = [rxnAliasList, mapStruct.rxnBaseReactantAlias{rxn}{x}];
        end

        % Loop over reactants
        for x = 1:length(mapStruct.rxnReactantID{rxn})
            rxnAliasList = [rxnAliasList, mapStruct.rxnReactantAlias{rxn}{x}];
        end

        % Loop over base products
        for x = 1:length(mapStruct.rxnBaseProductID{rxn})
            rxnAliasList = [rxnAliasList, mapStruct.rxnBaseProductAlias{rxn}{x}];
        end

        % Loop over products
        for x = 1:length(mapStruct.rxnProductID{rxn})
            rxnAliasList = [rxnAliasList, mapStruct.rxnProductAlias{rxn}{x}];
        end
    end

    % Get the indexes of the needed molecules to color
    % Get the corresponding IDs of the species based on their Names
    specIdList = mapStruct.specID(ismember(mapStruct.specName, metList));
    % Get the correspoding Aliases of the molecules based on the species IDs
    mapAliasList = mapStruct.molAlias(ismember(mapStruct.molID, specIdList));
    % Get only the Aliases of the molecules implicated in the needed reactions
    neededAliasList = mapAliasList(ismember(mapAliasList, rxnAliasList));
    % Get the corresponding Indexes of the molecules in the needed reactions
    molIndexList = find(ismember(mapStruct.molAlias, neededAliasList));

    % Loop over molecules and change the color
    for x = molIndexList'
        mapStruct.molColor{x} = colors(newColor);
    end

end

author: opencobra

repo: cobratoolbox


text:
function varargout = McCabeThieleGUI(varargin)
% MCCABETHIELEGUI M-file for McCabeThieleGUI.fig
%      MCCABETHIELEGUI, by itself, creates a new MCCABETHIELEGUI or raises the existing
%      singleton*.
%
%      H = MCCABETHIELEGUI returns the handle to a new MCCABETHIELEGUI or the handle to
%      the existing singleton*.
%
%      MCCABETHIELEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MCCABETHIELEGUI.M with the given input arguments.
%
%      MCCABETHIELEGUI('Property','Value',...) creates a new MCCABETHIELEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before McCabeThieleGUI_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to McCabeThieleGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Last Modified by GUIDE v2.5 06-Jul-2010 12:55:42

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @McCabeThieleGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @McCabeThieleGUI_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

% --- Executes just before McCabeThieleGUI is made visible.
function McCabeThieleGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% 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 McCabeThieleGUI (see VARARGIN)

set(handles.mixturepanel_buttongroup,'SelectionChangeFcn',@mixturepanel_buttongroup_SelectionChangeFcn);

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = McCabeThieleGUI_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;
%% Mixture Panel

function mixturepanel_buttongroup_SelectionChangeFcn(hObject, eventdata)

handles = guidata(hObject);

switch get(eventdata.NewValue,'Tag')   % Get Tag of selected object
    case 'methanolwater_radiobutton'
        %equilibrium data
        xMW = '[0    0.0200    0.0400    0.0600    0.0800    0.1000    0.1500    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    0.9500    1.0000]';
        yMW = '[0    0.1340    0.2300    0.3040    0.3650    0.4180    0.5170    0.5790    0.6650    0.7290    0.7790    0.8250    0.8700    0.9150    0.9580    0.9790    1.0000]';
        set(handles.userdefined_input,'Enable','Off')
        set(handles.userdefined_input2,'Enable','Off')
        set(handles.userdefined_input,'String',xMW)
        set(handles.userdefined_input2,'String',yMW)
        %hold off
        %methanol-water mixture
        %polynomial regression and plot of equilibrium data
        %p = polyfit(xMW,yMW,6);
        %x = 0:.1:100;
        %y = polyval(p,x);
        %plot(xMW,yMW)
        %hold on
        %axis([0 1 0 1])
        %plot(x,y)
        %plot(x,x)
        %xlabel('Liquid Mole Fraction of Methanol')
        %ylabel('Vapor Mole Fraction of Methanol')
        %title('McCabe-Thiele Diagram for Methanol-Water')

    case 'ethanolwater_radiobutton'
        %equilibrium data
        xEW = '[0    0.0190    0.0721    0.0966    0.1238    0.1661    0.2337    0.2608    0.3273    0.3965    0.5198    0.5732    0.6763    0.7472    0.8943    1.0000]';
        yEW = '[0    0.1700    0.3891    0.4375    0.4704    0.5089    0.5445    0.5580    0.5826    0.6122    0.6599    0.6841    0.7385    0.7815    0.8943    1.0000]';
        xeq = num2str(xEW);
        yeq = num2str(yEW);
        set(handles.userdefined_input,'Enable','Off')
        set(handles.userdefined_input2,'Enable','Off')
        set(handles.userdefined_input,'String',xeq)
        set(handles.userdefined_input2,'String',yeq)
        %hold off
        %ethanol-water mixture
        %polynomial regression and plot of equilibrium data
        %p = polyfit(xEW,yEW,8);
        %x = 0:.1:100;
        %y = polyval(p,x);
        %plot(x,y)
        %plot(xEW,yEW)
        %hold on
        %axis([0 1 0 1])
        %plot(x,x)
        %xlabel('Liquid Mole Fraction of Ethanol')
        %ylabel('Vapor Mole Fraction of Ethanol')
        %title('McCabe-Thiele Diagram for Ethanol-Water')

    case 'acetonethanol_radiobutton'
        %equilibrium data
        xAE = '[0    0.1000    0.1500    0.2000    0.2500    0.3000    0.3500    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    1.0000]';
        yAE = '[0    0.2620    0.3480    0.4170    0.4780    0.5240    0.5660    0.6050    0.6740    0.7390    0.8020    0.8650    0.9290    1.0000]';
        set(handles.userdefined_input,'Enable','Off')
        set(handles.userdefined_input2,'Enable','Off')
        set(handles.userdefined_input,'String',xAE)
        set(handles.userdefined_input2,'String',yAE)
        %hold off
        %acetone-ethanol mixture
        %polynomial regression and plot of equilibrium data
        %p = polyfit(xAE,yAE,6);
        %x = 0:.1:100;
        %y = polyval(p,x);
        %plot(xAE,yAE)
        %hold on
        %plot(x,y)
        %plot(x,x)
        %axis([0 1 0 1])
        %xlabel('Liquid Mole Fraction of Acetone')
        %ylabel('Vapor Mole Fraction of Acetone')
        %title('McCabe-Thiele Diagram for Acetone-Ethanol')

    case 'userdefined_radiobutton'
        hold off
       %User Defined Input
        choice = menu('Input Equilibrium Data as a Row Vector or as a Relative Volatility (alpha)','Vector','Alpha');
        switch choice
            case 1
                set(handles.userdefined_input,'Enable','On')
                set(handles.userdefined_input2,'Enable','On')
                set(handles.userdefined_input,'String','[Liquid Mole Fractions]')
                set(handles.userdefined_input2,'String','[Vapor Mole Fractions]')
                warndlg('Enlose your vectors with brackets and single quotes', 'Warning');
                %disp('You must enter your own equilibrium data in vector form ')
                %MVC = input('Enter the name of the more volatile component ','s');
                %LVC = input('Enter the name of the less volatile component ','s');
                %xU = input('Enter liquid mole fraction data as a row vector (use brackets) \n');
                %a = length(xU);
                %fprintf('You have entered %1.0f "x" data points, you must enter an equal number of y data points \n',a)
                %yU = input('Enter vapor phase mole fraction data as a row vector (use brackets) \n');
                %b = length(yU);
                %if a == b
                    %plot(xU,yU,'o')
                    %hold on
                    %p = polyfit(xU,yU,6);
                    %x = 0:.1:100;
                    %y = polyval(p,x);
                    %plot(x,y)
                    %plot(x,x)
                    %x_label = sprintf('Liquid Phase Mole Fraction of %s',MVC);
                    %y_label = sprintf('Vapor Phase Mole Fraction of %s',MVC);
                    %t_itle = sprintf('McCabe-Thiele Diagram for %s-%s',MVC,LVC);
                    %xlabel(x_label)
                    %ylabel(y_label)
                    %title(t_itle)
                %else
                    %disp('Equilibrium data does not match, vectors must be equal lengths')
                %end
            case 2
                set(handles.userdefined_input,'Enable','On')
                set(handles.userdefined_input2,'String','[]')
                set(handles.userdefined_input,'String','<Enter Relative Volatility>')
        end
end
%updates the handles structure
guidata(hObject, handles);
function userdefined_input_Callback(hObject, eventdata, handles)
% --- Executes during object creation, after setting all properties.
function userdefined_input_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function userdefined_input2_Callback(hObject, eventdata, handles)
% --- Executes during object creation, after setting all properties.
function userdefined_input2_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%% Column Pressure Checkbox

%Column Pressure Input
function pressureinput_Callback(hObject, eventdata, handles)
% Hints: get(hObject,'String') returns contents of pressureinput as text
%        str2double(get(hObject,'String')) returns contents of pressureinput as a double

% --- Executes during object creation, after setting all properties.
function pressureinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function columnpressure_checkbox_Callback(hObject, eventdata, handles)

checkbox1_Status = get(handles.columnpressure_checkbox,'Value');
if checkbox1_Status == 0
    set(handles.pressureinput,'String',' ')
else
    set(handles.pressureinput,'Enable','On')
    set(handles.pressureinput,'String','<Enter Column Pressure>')
end

%% Feed
%% Feed Flow Rate
%Feed Flow Rate Input
function feedflowinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function feedflowinput_CreateFcn(hObject, eventdata, handles)

    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
    end

%Feed Flow Rate Checkbox
function feedflowrate_checkbox_Callback(hObject, eventdata, handles)

checkbox2_Status = get(handles.feedflowrate_checkbox,'Value');

if checkbox2_Status == 0
    set(handles.feedflowinput,'Enable','Off')
    set(handles.feedflowinput,'String',' ');
else
    set(handles.feedflowinput,'Enable','On')
    set(handles.feedflowinput,'String','<Enter Feed Flow Rate>')
end

%% Feed Composition of the Volatile Component Input
function feedcompvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function feedcompvolatile_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Feed Composition of the Less Volatile Component Input
function feedcomplessvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function feedcomplessvolatile_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Feed Composition Checkbox
function feedcomposition_checkbox_Callback(hObject, eventdata, handles)

checkbox3_Status = get(handles.feedcomposition_checkbox,'Value');

if checkbox3_Status == 0
    set(handles.feedcompvolatile,'Enable','Off')
    set(handles.feedcomplessvolatile,'Enable','Off')
    set(handles.feedcompvolatile,'String',' ')
    set(handles.feedcomplessvolatile,'String',' ')
else
    set(handles.feedcompvolatile,'Enable','On')
    set(handles.feedcomplessvolatile,'Enable','On')
    set(handles.feedcompvolatile,'String','<Volatile>')
    set(handles.feedcomplessvolatile,'String','<Non-Volatile>')
end

%% Feed Enthalpy
%Feed Enthalpy Input
function feedenthalpyinput_Callback(hObject, eventdata, handles)
function feedenthalpyinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Feed Enthalpy Checkbox
function feedenthalpy_Callback(hObject, eventdata, handles)

checkbox4_Status = get(handles.feedenthalpy,'Value');

if checkbox4_Status == 0
    set(handles.feedenthalpyinput,'Enable','Off')
    set(handles.feedenthalpyinput,'String',' ')
else
    set(handles.feedenthalpyinput,'Enable','On')
    set(handles.feedenthalpyinput,'String','<Feed Enthalpy>')
end

%Reboiler Duty Input
function reboilerdutyinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function reboilerdutyinput_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% Reboiler Duty Checkbox
function reboilerdutycheckbox_Callback(hObject, eventdata, handles)

checkbox7_Status = get(handles.reboilerdutycheckbox,'Value');
if checkbox7_Status == 0
    set(handles.reboilerdutyinput,'Enable','Off')
    set(handles.reboilerdutyinput,'String',' ')
else
    set(handles.reboilerdutyinput,'Enable','On')
    set(handles.reboilerdutyinput,'String','<Reboiler Duty>')
end

%% Distillate composition

%Distillate Composition, Volatile Component Input
function distcompvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function distcompvolatile_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Distillate Composition of the less Volatile Component Input
function distcomplessvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function distcomplessvolatile_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Distillate composition checkbox
function distillatecomposition_checkbox_Callback(hObject, eventdata, handles)

checkbox9_Status = get(handles.distillatecomposition_checkbox,'Value');

if checkbox9_Status == 0
    set(handles.distcompvolatile,'Enable','Off')
    set(handles.distcomplessvolatile,'Enable','Off')
    set(handles.distcompvolatile,'String',' ')
    set(handles.distcomplessvolatile,'String',' ')
else
    set(handles.distcompvolatile,'Enable','On')
    set(handles.distcomplessvolatile,'Enable','On')
    set(handles.distcompvolatile,'String','<Volatile>')
    set(handles.distcomplessvolatile,'String','<Non-Volatile>')
end

%% Reflux Ratio
%Reflux Ratio Input
function refluxratio_input_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function refluxratio_input_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function minimumrefluxinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function minimumrefluxinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Reflux Ratio Checkbox
function refluxratio_checkbox_Callback(hObject, eventdata, handles)

checkbox11_Status = get(handles.refluxratio_checkbox,'Value');

if checkbox11_Status == 0
    set(handles.refluxratio_input,'Enable','Off')
    set(handles.refluxratio_input,'String',' ')
    set(handles.minimumrefluxinput,'Enable','Off')
    set(handles.minimumrefluxinput,'String',' ')
elseif checkbox11_Status == 1
    choice = menu('Are you specifying minimum reflux or actual reflux','Actual','Minumum');
    switch choice
        case 1
            set(handles.refluxratio_input,'Enable','On')
            set(handles.refluxratio_input,'String','<Reflux Ratio>')
            set(handles.minimumrefluxinput,'Enable','Off')
            set(handles.minimumrefluxinput,'String',' ')
        case 2
            Rmin = str2double(get(handles.minimumrefluxinput,'String'));
            if isnan(Rmin) == 1
                set(handles.refluxratio_input,'Enable','On')
                set(handles.refluxratio_input,'String','<Operating Factor>')
                set(handles.minimumrefluxinput,'Enable','On')
                set(handles.minimumrefluxinput,'String','<Minumum Reflux>')
            elseif isnan(Rmin) == 0
                set(handles.refluxratio_input,'Enable','On')
                set(handles.refluxratio_input,'String','<Operating Factor>')
                set(handles.minimumrefluxinput,'Enable','Off')
            end
    end

end

%% Bottoms Composition
%Bottoms Composition of the Less Volatile Component Input
function botcomplessvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function botcomplessvolatile_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Bottoms Composition of the Volatile Component Input
function botcompvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function botcompvolatile_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Bottoms Composition Checkbox
function bottomscomposition_checkbox_Callback(hObject, eventdata, handles)

checkbox10_Status = get(handles.bottomscomposition_checkbox,'Value');

if checkbox10_Status == 0
    set(handles.botcomplessvolatile,'Enable','Off')
    set(handles.botcompvolatile,'Enable','Off')
    set(handles.botcomplessvolatile,'String',' ')
    set(handles.botcompvolatile,'String',' ')
elseif checkbox10_Status == 1
    set(handles.botcomplessvolatile,'Enable','On')
    set(handles.botcompvolatile,'Enable','On')
    set(handles.botcompvolatile,'String','<Volatile>')
    set(handles.botcomplessvolatile,'String','<Non-Volatile>')
end

%% Optimum Feed Plate
%Optimum Feed Plate Checkbox
function optimumfeedplate_Callback(hObject, eventdata, handles)

%% Distillate Flow Rate
%Distillate Flow Rate Input
function distillateflowrateinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function distillateflowrateinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Distillate Flow Rate Checkbox
function distillateflowrate_checkbox_Callback(hObject, eventdata, handles)
checkbox14_Status = get(handles.distillateflowrate_checkbox,'Value');

if checkbox14_Status == 0
    set(handles.distillateflowrateinput,'Enable','Off')
    set(handles.distillateflowrateinput,'String',' ')
else
    set(handles.distillateflowrateinput,'Enable','On')
    set(handles.distillateflowrateinput,'String','<Enter Distillate Flow Rate>')
end

%% Fractional Recoveries

%Fractional Recovery of the Volatile Component in the Distillate Input
function fracrecoveryvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function fracrecoveryvolatile_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Fractional Recovery of the Less Volatile Component in the Bottoms Input
function fracrecoverylessvolatile_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function fracrecoverylessvolatile_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Fractional Recoveries Checkbox
function fractionalrecoveries_checkbox_Callback(hObject, eventdata, handles)

checkbox12_Status = get(handles.fractionalrecoveries_checkbox,'Value');

if checkbox12_Status == 0
    set(handles.fracrecoveryvolatile,'Enable','Off')
    set(handles.fracrecoverylessvolatile,'Enable','Off')
    set(handles.fracrecoveryvolatile,'String','0')
    set(handles.fracrecoverylessvolatile,'String','0')
else
    set(handles.fracrecoveryvolatile,'Enable','On')
    set(handles.fracrecoverylessvolatile,'Enable','On')
    set(handles.fracrecoveryvolatile,'String','<Volatile>')
    set(handles.fracrecoverylessvolatile,'String','<Non-Volatile>')
end

%% Bottoms Flow Rate
function bottomsflowrateinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function bottomsflowrateinput_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Bottoms Flow Rate Checkbox
function bottomsflowrate_checkbox_Callback(hObject, eventdata, handles)
checkbox15_Status = get(handles.bottomsflowrate_checkbox,'Value');
if checkbox15_Status == 0
    set(handles.bottomsflowrateinput,'Enable','Off')
    set(handles.bottomsflowrateinput,'String',' ')
else
    set(handles.bottomsflowrateinput,'Enable','On')
    set(handles.bottomsflowrateinput,'String','<Enter Bottoms Flow Rate>')
end

%% Boilup Ratio Checkbox

%Boilup Ratio Input
function boilupratioinput_Callback(hObject, eventdata, handles)

function boilupratioinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function boilupratio_checkbox_Callback(hObject, eventdata, handles)

checkbox15_Status = get(handles.boilupratio_checkbox,'Value');

if checkbox15_Status == 0
    set(handles.boilupratioinput,'Enable','Off')
    set(handles.boilupratioinput,'String',' ')
elseif checkbox15_Status == 1
    set(handles.boilupratioinput,'Enable','On')
    set(handles.boilupratioinput,'String','<Enter Boilup Ratio>')
elseif checkbox15_Status == 0 && checkbox11_Status == 1 && checkbox6_Status == 1 && checkbox9_Status == 1
    set(handles.boilupratioinput,'Enable','Off')
    set(handles.boilupratioinput,'String',' ')
    xi = (-(Q-1)*(1-R/(R+1))*xd-zf)/((Q-1)*R/(R+1)-Q);
    yi = (zf+xd*Q/R)/(1+Q/R);
    xs = linspace(xb,xi);
    ys = linspace(xb,yi);
    plot(xs,ys)
    hold on
end

%% Feed Quality

%Feed Quality Input
function feedqualityinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function feedqualityinput_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Feed Quality Checkbox
function feedquality_checkbox_Callback(hObject, eventdata, handles)

checkbox6_Status = get(handles.feedquality_checkbox,'Value');

if checkbox6_Status == 0
    set(handles.feedqualityinput,'Enable','Off')
    set(handles.feedqualityinput,'String',' ')
else
    set(handles.feedqualityinput,'Enable','On')
    set(handles.feedqualityinput,'String','<Enter Feed Quality>')
end

%Condenser Duty Input
function condenserdutyinput_Callback(hObject, eventdata, handles)
function condenserdutyinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Condenser Duty Checkbox
function condenserdutycheckbox_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function condenserdutycheckbox_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

checkbox5_Status = get(handles.condenserdutycheckbox,'Value');
if checkbox5_Status == 0
    set(handles.condenserdutyinput,'Enable','Off')
    set(handles.condenserdutyinput,'String',' ')
else
    set(handles.condenserdutyinput,'Enable','On')
    set(handles.condenserdutyinput,'String','<Enter Condenser Duty>')
end

%Reflux Enthalpy Input
function refluxenthalpyinput_Callback(hObject, eventdata, handles)
% --- Executes during object creation, after setting all properties.
function refluxenthalpyinput_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
%Reflux Enthalpy Checkbox
function refluxenthalpy_Callback(hObject, eventdata, handles)

checkbox8_Status = get(handles.refluxenthalpy,'Value');
if checkbox8_Status == 0
    set(handles.refluxenthalpyinput,'Enable','Off')
    set(handles.refluxenthalpyinput,'String',' ')
else
    set(handles.refluxenthalpyinput,'Enable','On')
    set(handles.refluxenthalpyinput,'String','<Reflux Enthalpy>')
end

%Boilup Enthalpy Input
function boilupenthalpyinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function boilupenthalpyinput_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Boilup Enthalpy Checkbox
function boilupenthalpycheckbox_Callback(hObject, eventdata, handles)
checkbox20_Status = get(handles.boilupenthalpycheckbox,'Value');
if checkbox20_Status == 0
    set(handles.boilupenthalpyinput,'Enable','Off')
    set(handles.boilupenthalpyinput,'String',' ')
else
    set(handles.boilupenthalpyinput,'Enable','On')
    set(handles.boilupenthalpyinput,'String','<Boilup Enthalpy>')
end

%Distillate Enthalpy
function distillateenthalpyinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function distillateenthalpyinput_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in distillateenthalpycheckbox.
function distillateenthalpycheckbox_Callback(hObject, eventdata, handles)
checkbox21_Status = get(handles.distillateenthalpycheckbox,'Value');
if checkbox21_Status == 0
    set(handles.distillateenthalpyinput,'Enable','Off')
    set(handles.distillateenthalpyinput,'String',' ')
else
    set(handles.distillateenthalpyinput,'Enable','On')
    set(handles.distillateenthalpyinput,'String','<Distillate Enthalpy>')
end

%Bottoms enthalpy
function bottomsenthalpyinput_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function bottomsenthalpyinput_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in bottomsenthalpycheckbox.
function bottomsenthalpycheckbox_Callback(hObject, eventdata, handles)
checkbox22_Status = get(handles.bottomsenthalpycheckbox,'Value');
if checkbox22_Status == 0
    set(handles.bottomsenthalpyinput,'Enable','Off')
    set(handles.bottomsenthalpyinput,'String',' ')
else
    set(handles.bottomsenthalpyinput,'Enable','On')
    set(handles.bottomsenthalpyinput,'String','<Bottoms Enthalpy>')
end

% --- Executes on button press in startbutton.
function startbutton_Callback(hObject, eventdata, handles)

% Get Values

Alpha = str2double(get(handles.userdefined_input,'String'));

P = str2double(get(handles.pressureinput,'String'));
F = str2double(get(handles.feedflowinput,'String'));
zf = str2double(get(handles.feedcompvolatile,'String'));
zfnv = str2double(get(handles.feedcomplessvolatile,'String'));
xd = str2double(get(handles.distcompvolatile,'String'));
xdnv = str2double(get(handles.distcomplessvolatile,'String'));
xb = str2double(get(handles.botcompvolatile,'String'));
xbnv = str2double(get(handles.botcomplessvolatile,'String'));
D = str2double(get(handles.distillateflowrateinput,'String'));
fv = str2double(get(handles.fracrecoveryvolatile,'String'));
fnv = str2double(get(handles.fracrecoverylessvolatile,'String'));
B = str2double(get(handles.bottomsflowrateinput,'String'));
S = str2double(get(handles.boilupratioinput,'String'));
Q = str2double(get(handles.feedqualityinput,'String'));
R = str2double(get(handles.refluxratio_input,'String'));
Rmin = str2double(get(handles.minimumrefluxinput,'String'));
hf = str2double(get(handles.feedenthalpyinput,'String'));
hD = str2double(get(handles.refluxenthalpyinput,'String'));
hB = str2double(get(handles.boilupenthalpyinput,'String'));
Qc = str2double(get(handles.condenserdutyinput,'String'));
Qr = str2double(get(handles.reboilerdutyinput,'String'));
HD = str2double(get(handles.distillateenthalpyinput,'String'));
HB = str2double(get(handles.bottomsenthalpyinput,'String'));

%xie = 1/2*(-Q+zf*Alpha+Q*Alpha-zf-Alpha+(2*Q*Alpha^2*zf+Q^2-2*Q*zf+2*zf*Alpha-2*zf^2*Alpha+zf^2*Alpha^2-2*zf*Alpha^2+zf^2+2*Q*Alpha-2*Q^2*Alpha+Alpha^2+Q^2*Alpha^2-2*Q*Alpha^2)^(1/2))/Q/(-1+Alpha);
%yie = xie*Alpha/(1-xie+xie*Alpha);
xi = (-(Q-1)*(1-R/(R+1))*xd-zf)/((Q-1)*R/(R+1)-Q);
yi = (zf+xd*Q/R)/(1+Q/R);

%Error Tests

%Feed Quality & Feed Line
checkbox6_Status = get(handles.feedquality_checkbox,'Value');
if checkbox6_Status == 0
    errordlg('Feed quality must be specified')
    return
elseif checkbox6_Status == 1 && isnan(Q) == 1
    errordlg('Feed quality must be specified')
    return
elseif checkbox6_Status == 1 && isnan(Q) == 0
    ye = 0:.001:1;
    if Q > -inf && Q <1
            xf = 0:.001:1;
    elseif Q > 1
            xf = zf:.001:1;
    end
    if Q ~= 1
        yf = yfeed(xf,Q,zf);
    elseif Q == 1
        yf = 0:.001:1;
        xf = zf*ones(length(yf));
    end
end

%Equilibrium Data
if any(isnan(Alpha)) == 1
    xeq = str2num(get(handles.userdefined_input,'String'));
    yeq = str2num(get(handles.userdefined_input2,'String'));
    %Error Test on Equilibrium Vector Length
    if length(yeq) ~= length(xeq) && isempty(yeq) == 1 && length(xeq) > 1
        errordlg('Equilibrium data does not match')
        return
    elseif length(yeq) == length(xeq) && length(xeq) < 10
        errordlg('Too few data points will give a poor equilibrium curve')
    else
        k = 1; %k is a constant, metric for stepping off equilibrium stages
    end
elseif isnan(Alpha) == 0 && length(Alpha) == 1
    yeq = [];
    if Alpha <= 1
        errordlg('Relative volatility must be greater than 1')
        return
    elseif Alpha > 1
        k = 2;
    end
end
if k == 1
    xx = 0:.001:1;
    yy = spline(xeq,yeq,xx);
    if Q == 0
        xx = .001:.001:1.001;
        yy = spline(xeq,yeq,xx);
        diffe = abs(yy-zf);
        [a b] = min(diffe) ;
        xie = xx(b);
        yie = zf;
    elseif Q == 1
        yie = spline(xeq,yeq,zf);
        xie = zf;
    elseif Q~=0 && Q~=1
        diffe = abs(yy-yf);
        [a b] = min(diffe);
        xie = xx(b);
        yie = yy(b);
    end
elseif k == 2
    xie = 1/2*(-Q+zf*Alpha+Q*Alpha-zf-Alpha+(2*Q*Alpha^2*zf+Q^2-2*Q*zf+2*zf*Alpha-2*zf^2*Alpha+zf^2*Alpha^2-2*zf*Alpha^2+zf^2+2*Q*Alpha-2*Q^2*Alpha+Alpha^2+Q^2*Alpha^2-2*Q*Alpha^2)^(1/2))/Q/(-1+Alpha);
    yie = xie*Alpha/(1-xie+xie*Alpha);
end

%Overall Mass Balance
checkbox2_Status = get(handles.feedflowrate_checkbox,'Value');
checkbox17_Status = get(handles.distillateflowrate_checkbox,'Value');
checkbox18_Status = get(handles.bottomsflowrate_checkbox,'Value');
l = [checkbox2_Status checkbox17_Status checkbox18_Status]; %l is a status vector for mass flow streams
if any([F D B]< 0)
    errordlg('Mass flow rates must be greater than zero')
    return
elseif sum(l) == 3 && D + B ~= F
    errordlg('The mass balance is incorrect')
    return
elseif sum(l) == 2
    if isnan(F) == 0
        if isnan(B) == 0
            D = F - B;
            Ds = num2str(D);
            set(handles.distillateflowrateinput,'String',Ds)
        elseif isnan(D) == 0
            B = F - D;
            Bs = num2str(B);
            set(handles.bottomsflowrateinput,'String',Bs)
        end
    elseif isnan(F) == 1
        F = B + D;
        Fs = num2str(F);
        set(handles.feedflowinput,'String',Fs)
    end
end

%Feed Composition
checkbox3_Status = get(handles.feedcomposition_checkbox,'Value');
if checkbox3_Status == 1
    if isnan(zf) == 1 && isnan(zfnv) == 1
        errordlg('Feed composition must be a real number')
        return
    elseif isnan(zf) == 0 && isnan(zfnv) == 1
        zfnv = 1 - zf;
        zfnv = num2str(zfnv);
        set(handles.feedcomplessvolatile,'String',zfnv)
    elseif isnan(zf) == 1 && isnan(zfnv) == 0
        zf = 1 - zfnv;
        zf = num2str(zf);
        set(handles.feedcompvolatile,'String',zf)
    elseif isnan(zf) == 0 && isnan(zfnv) == 0
        if zf + zfnv ~= 1
            errordlg('Mole fractions must add up to 1')
            return
        elseif zf < 0 || zf > 1 || zfnv < 0 || zfnv > 1
            errordlg('Mole fractions must be between 0 and 1')
            return
        end
    end
end

%Distillate Composition
checkbox9_Status = get(handles.distillatecomposition_checkbox,'Value');
if checkbox9_Status == 1
    if isnan(xd) && isnan(xdnv)
        errordlg('Distillate composition must be a real number')
        return
    elseif isnan(xd) == 0 && isnan(xdnv) == 1
        xdnv = 1 - xd;
        xdnvs = num2str(xdnv);
        set(handles.distcomplessvolatile,'String',xdnvs)
    elseif isnan(xd) == 1 && isnan(xdnv) == 0
        xd = 1 - xdnv;
        xds = num2str(xd);
        set(handles.distcompvolatile,'String',xds)
    elseif isnan(xd) == 0 && isnan(xdnv) == 0
        if xd + xdnv ~= 1
            errordlg('Mole fractions must add up to 1')
            return
        elseif xd < 0 || xd > 1 || xdnv < 0 || xdnv > 1
            errordlg('Mole fractions must be between 0 and 1')
            return
        end
    end
end

%Bottoms Composition
checkbox10_Status = get(handles.bottomscomposition_checkbox,'Value');
if checkbox10_Status == 1
    if isnan(xb) == 1 && isnan(xbnv) == 1
        errordlg('Bottoms composition must be a real number')
        return
    elseif isnan(xb) == 0 && isnan(xbnv) == 1
        xbnv = 1 - xb;
        xbnvs = num2str(xbnv);
        set(handles.botcomplessvolatile,'String',xbnvs)
    elseif isnan(xb) == 1 && isnan(xbnv) == 0
        xb = 1 - xbnv;
        xbs = num2str(xb);
        set(handles.botcompvolatile,'String',xbs)
    elseif isnan(xb) == 0 && isnan(xbnv) == 0
        if xb + xbnv ~= 1
            errordlg('Mole fractions must add up to 1')
            return
        elseif xb < 0 || xb > 1 || xbnv < 0 || xbnv > 1
            errordlg('Mole fractions must be between 0 and 1')
            return
        end
    end
end

%Test to see if the mass balance has converged
%lk is a status vector for feed stream compositions
lk = [checkbox3_Status checkbox9_Status checkbox10_Status];
if sum(lk) == 3
    if sum(l) == 1
        if isnan(F) == 0
            if isnan(D) == 1 && isnan(B) == 1
                D = F*(zf-xb)/(xd-xb);
                Ds = num2str(D);
                set(handles.distillateflowrateinput,'String',Ds)
                B = F - D;
                Bs = num2str(B);
                set(handles.bottomsflowrateinput,'String',Bs)
                fv = D*xd/F/zf;
                fvs = num2str(fv);
                set(handles.fracrecoveryvolatile,'String',fvs)
                fnv = B*(1-xb)/F/(1-zf);
                fnvs = num2str(fnv);
                set(handles.fracrecoverylessvolatile,'String',fnvs)
            end
        end
    elseif sum(l) == 2
        if isnan(F) == 0
            if isnan(D) == 0
                Dk = F*(zf-xb)/(xd-xb);
                Bk = F - Dk;
                if D ~= Dk
                    Ds = num2str(Dk);
                    set(handles.distillateflowrateinput,'String',Ds)
                    Bk = F-Dk;
                    Bs = num2str(Bk);
                    set(handles.bottomsflowrateinput,'String',Bs)
                    fv = D*xd/F/zf;
                    fvs = num2str(fv);
                    set(handles.fracrecoveryvolatile,'String',fvs)
                    fnv = B*(1-xb)/F/(1-zf);
                    fnvs = num2str(fnv);
                    set(handles.fracrecoverylessvolatile,'String',fnvs)
                end
                if B ~= Bk
                    Bs = num2str(Bk);
                    set(handles.bottomsflowrateinput,'String',Bs)
                    Dk = F-Bk;
                    Ds = num2str(Dk);
                    set(handles.distillateflowrateinput,'String',Ds)
                    fv = D*xd/F/zf;
                    fvs = num2str(fv);
                    set(handles.fracrecoveryvolatile,'String',fvs)
                    fnv = B*(1-xb)/F/(1-zf);
                    fnvs = num2str(fnv);
                    set(handles.fracrecoverylessvolatile,'String',fnvs)
                end
            end
        end
    end
elseif sum(lk) == 1 && isnan(zf) == 0
checkbox12_Status = get(handles.fractionalrecoveries_checkbox,'Value');
    if sum(l) == 2 || sum(l) == 3
        %Fractional Recoveries
        if  checkbox12_Status == 1
            if any(isnan([fv fnv])) == 1
                errordlg('Fractional recovery must be a real number')
                return
            elseif any([fv fnv]>=1) || any([fv fnv]<0)
                errordlg('Fractional recovery must be between 0 and 1')
                return
            elseif any([fv fnv]==1)
                errordlg('100 percent recovery is physically impossible, try a different value')
                return
            else
                xd = fv/D*F*zf;
                xds = num2str(xd);
                set(handles.distcompvolatile,'String',xds)
                xb = 1-fnv/(F-D)*F*(1-zf);
                xbs = num2str(xb);
                set(handles.botcompvolatile,'String',xbs)
            end
        end
    elseif checkbox12_Status == 0
        errordlg('The system is underspecified')
        return
    end
end

%Fractional Recoveries
checkbox12_Status = get(handles.fractionalrecoveries_checkbox,'Value');
if checkbox12_Status == 0
    if all(isnan([F D B zf xd xb])) == 0
        fv = D*xd/F/zf;
        fvs = num2str(fv);
        set(handles.fracrecoveryvolatile,'String',fvs)
        fnv = B*(1-xb)/F/(1-zf);
        fnvs = num2str(fnv);
        set(handles.fracrecoverylessvolatile,'String',fnvs)
    end
end

%Feed Enthalpy
checkbox4_Status = get(handles.feedenthalpy,'Value');
if checkbox4_Status == 1 && isnan(hf) == 1
    errordlg('Feed enthalpy must be a real number')
    return
end

%Condenser Duty
checkbox5_Status = get(handles.condenserdutycheckbox,'Value');
if checkbox5_Status == 1 && isnan(Qc) == 1
    errordlg('Condenser duty must be a real number')
    return
end

%Reboiler Duty
checkbox7_Status = get(handles.reboilerdutycheckbox,'Value');
if checkbox7_Status == 1 && isnan(Qr) == 1
    errordlg('Reboiler duty must be a real number')
    return
end

%Reflux Enthalpy
checkbox8_Status = get(handles.refluxenthalpy,'Value');
if checkbox8_Status == 1 && isnan(hD) == 1
    errordlg('Reflux enthalpy must be a real number')
    return
end

%Reflux Ratio
checkbox11_Status = get(handles.refluxratio_checkbox,'Value');
if checkbox11_Status == 1
    if isnan(R) == 1 && isnan(Rmin) == 1
        errordlg('Reflux ratio must be a real number')
        return
    elseif R < 0
        errordlg('Reflux ratio must be greater than 0')
        return
    elseif Rmin < 0
        errordlg('Minumum reflux must be greater than 0')
        return
    elseif isnan(Rmin) == 0 && isnan(R) == 1
        errordlg('You must specify an operating factor')
        return
    end
elseif checkbox11_Status == 0
    if isnan(xd) == 1 && isnan(Q) == 1
        errordlg('You must specify feed quality and distillate composition')
        return
    elseif isnan(xd) == 0 && isnan(Q) == 1
        errordlg('You must specify feed quality')
        return
    elseif isnan(xd) == 1 && isnan(Q) == 0
        errordlg('You must specify distillate composition')
        return
    elseif isnan(Rmin) == 1 && isnan(Q) == 0 && isnan(xd) == 0
        if isnan(Alpha) == 0
            if isnan(R) == 1
                Emin = (xd-yie)/(xd-xie);
                Rmin = Emin/(1-Emin);
                Rmins = num2str(Rmin);
                set(handles.minimumrefluxinput,'String',Rmins)
                errordlg('Minimum reflux has been calculated and returned but you must specify an operating factor')
                return
            elseif isnan(R) == 0
                Emin = (xd-yie)/(xd-xie);
                Rmin = Emin/(1-Emin);
                Rmins = num2str(Rmin);
                R = R*Rmin;
                Rs = num2str(R);
                set(handles.minimumrefluxinput,'String',Rmins)
                set(handles.refluxratio_input,'String',Rs)
            elseif isnan(Rmin) == 0 && isnan(R) == 1
                errordlg('You must specify an operating factor')
                return
            end
        elseif isnan(Alpha) == 1 && k ==1
            if isnan(R) == 1
                Emin = (xd-yie)/(xd-xie);
                Rmin = Emin/(1-Emin);
                Rmins = num2str(Rmin);
                set(handles.minimumrefluxinput,'String',Rmins)
                errordlg('Minimum reflux has been calculated and returned but you must specify an operating factor')
                return
            elseif isnan(R) == 0
                Emin = (xd-yie)/(xd-xie);
                Rmin = Emin/(1-Emin);
                Rmins = num2str(Rmin);
                R = R*Rmin;
                Rs = num2str(R);
                set(handles.minimumrefluxinput,'String',Rmins)
                set(handles.refluxratio_input,'String',Rs)
            elseif isnan(Rmin) == 0 && isnan(R) == 1
                errordlg('You must specify an operating factor')
                return
            end
        end
    end
end

%Optimum Feed Plate
checkbox12_Status = get(handles.optimumfeedplate,'Value');
if checkbox12_Status == 0
    errordlg('This version does not support non-optimum feed plate')
    return
end

%Boilup Ratio
if S < 0
    errordlg('Boilup ratio must be greater than 0')
    return
end

%Boilup Enthalpy
checkbox20_Status = get(handles.boilupenthalpycheckbox,'Value');
if checkbox20_Status == 1 && isnan(hB) == 1
    errordlg('Boilup enthalpy must be a real number')
    return
end

%Pressure
if P < 0
    errordlg('Pressure must be greater than 0')
end

%Energy Balance
%checkbox4_Status = get(handles.feedenthalpy,'Value');
%checkbox5_Status = get(handles.condenserdutycheckbox,'Value');
%checkbox7_Status = get(handles.reboilerdutycheckbox,'Value');
%checkbox8_Status = get(handles.refluxenthalpy,'Value');
%checkbox20_Status = get(handles.boilupenthalpycheckbox,'Value');
%checkbox21_Status = get(handles.distillateenthalpycheckbox,'Value');
%checkbox22_Status = get(handles.bottomsenthalpycheckbox,'Value');
%jkl = [checkbox4_Status checkbox5_Status checkbox7_Status checkbox8_Status checkbox20_Status checkbox21_Status checkbox22_Status]; %Status vector
%Condenser Energy Balance
if isnan(Qc) == 0
    if isnan(hD) == 1
        hD = Qc/D/(R+1)*3600;
        hDs = num2str(hD);
        set(handles.refluxenthalpyinput,'String',hDs)
    end
elseif isnan(hD) == 0
    if isnan(Qc) == 1
        Qc = hB*D*(R+1)/3600;
        Qcs = num2str(Qc);
        set(handles.condenserdutyinput,'String',Qcs)
    end
end
%Reboiler Energy Balance
if isnan(Qr) == 0
    if isnan(hB) == 1
        hB = Qr/B/S*3600;
        hBs = num2str(hB);
        set(handles.boilupenthalpyinput,'String',hBs)
    end
elseif isnan(hB) == 0
    if isnan(Qr) == 1
        Qr = B*S*hB/3600;
        Qrs = num2str(Qr);
        set(handles.reboilerdutyinput,'String',Qrs)
    end
end
%jklm = [isnan(F) isnan(hf) isnan(Qc) isnan(Qr) isnan(B) isnan(HB) isnan(D) isnan(HD)];
%if sum(jklm) == 0
%    if F*hf + Qc + Qr - B*HB - D*HD >= 100
%        errordlg('The energy balance is not correct')
%        return
%    end
%end

%Equilibrium Curves
if k ==1 %Vector
    xx = 0:.001:1;
    yy = spline(xeq,yeq,xx);
    hold off
    plot(xx,yy)
    axis([0 1 0 1])
    hold on
    plot(xx,xx)
elseif k == 2 %Alpha
    ye = 0:.001:1;
    xe = equilib(ye,Alpha);
    hold off
    plot(xe,ye)
    axis([0 1 0 1])
    hold on
    plot(ye,ye)
end

%Plotting the Operating Lines
ytopop = ytop(ye,R,xd);
ybotop = ybot(ye,S,xb);
plot(xf,yf);
plot(ye,ytopop)
plot(ye,ybotop)
plot(ye,ye)
plot(xb,xb,'o')
text(xb+.025,xb-.025,'x_B')
plot(xd,xd,'o')
text(xd+.025,xd-.025,'x_D')
plot(zf,zf,'o')

% Computing the intersection of feed line and operating lines
plot(xi,yi,'o')
text(xi-.1,yi+.025,'Feed Line')
text(xb-.025,xb+.05,'Stripping Line')
text(xd-.25,xd-.05,'Enriching Line')
xlabel('Liquid Phase Mole Fraction, x')
ylabel('Vapor Phase Mole Fraction, y')
title('McCabe-Thiele Diagram for Binary Column Distillation')
if k == 2
    yi2 = interp1 (xe,ye,xi);
    if yi > yi2
        errordlg ('The distillation is not possible. Try a different operating condition.')
        return
    end
elseif k == 1
    yi2 = interp1(xx,yy,xi);
    if yi > yi2
        errordlg('The distillation is not possible. Try a different operating condition.')
        return
    end
end

% Rectifying section
if k == 2 %user defined systems
    i = 1;
    xp(1) = xd;
    yp(1) = xd;
    y = xd;
    while xp (i) > xi
        xp(i+1)= equilib(y,Alpha);
        yp(i+1)= R/(R+1)*xp (i+1)+xd/(R+1);
        y = yp (i+1);
        a = linspace(xp(i),xp(i+1));
        b = yp(i)*ones(length(a));
        plot(a,b)
        text ((xp(i+1)-.025),(yp(i)+.025),num2str (i))
        if xp (i+1) > xi
            c = linspace(yp(i),yp(i+1));
            d = xp(i+1)*ones(length(c));
            plot(d,c)
        end
        i = i+1;
    end
    % Stripping section
    xs = linspace(xb,xi);
    ys = linspace(xb,yi);
    plot(xs,ys)
    ss = (yi-xb)/(xi-xb);
    S = 1/(ss-1);
    S = num2str(S);
    set(handles.boilupratioinput,'String',S)
    yp(i) = ss*(xp(i)-xb)+xb;
    y = yp(i);
    a = linspace(yp(i-1),yp(i));
    b = xp(i)*ones(length(a));
    plot(b,a)
    while xp (i) > xb
        xp(i+1) = equilib(y,Alpha);
        yp(i+1) = ss*(xp(i+1)-xb)+xb;
        y = yp(i+1);
        a = linspace(xp(i),xp(i+1));
        b = yp(i)*ones(length(a));
        plot(a,b)
        text ((xp(i+1)-.025),(yp(i)+.025),num2str (i))
        if xp (i+1) > xb
            a = linspace(yp(i),yp(i+1));
            b = xp(i+1)*ones(length(a));
            plot(b,a)
        end
        i = i+1;
    end

elseif k ==1 %predefined equilibrium systems
    i = 1;
    xp(1) = xd;
    yp(1) = xd;
    y = xd;
    while xp(i) > xi
        xp(i+1)= predefeq(xeq,yeq,y);
        yp(i+1)= R/(R+1)*xp(i+1)+xd/(R+1);
        y = yp(i+1);
        a = linspace(xp(i),xp(i+1));
        b = yp(i)*ones(length(a));
        plot(a,b)
        text ((xp(i+1)-.025),(yp(i)+.025),num2str (i))
        if xp (i+1) > xi
            c = linspace(yp(i),yp(i+1));
            d = xp(i+1)*ones(length(c));
            plot(d,c)
        end
        i = i+1;
    end
    % Stripping section
    xs = linspace(xb,xi);
    ys = linspace(xb,yi);
    plot(xs,ys)
    ss = (yi-xb)/(xi-xb);
    S = 1/(ss-1);
    S = num2str(S);
    set(handles.boilupratioinput,'String',S)
    yp(i) = ss*(xp(i)-xb)+xb;
    y = yp(i);
    a = linspace(yp(i-1),yp(i));
    b = xp(i)*ones(length(a));
    plot(b,a)
    while xp (i) > xb
        xp(i+1) = predefeq(xeq,yeq,y);
        yp(i+1) = ss*(xp(i+1)-xb)+xb;
        y = yp(i+1);
        a = linspace(xp(i),xp(i+1));
        b = yp(i)*ones(length(a));
        plot(a,b)
        text ((xp(i+1)-.025),(yp(i)+.025),num2str (i))
        if xp (i+1) > xb
            a = linspace(yp(i),yp(i+1));
            b = xp(i+1)*ones(length(a));
            plot(b,a)
        end
        i = i+1;
    end
end

% Functions called with start button: operating & equilibrium lines
function yf = yfeed(xf,Q,zf)
yf = Q/(Q-1)*xf+zf/(1-Q);

function xa = equilib(ya,Alpha)
xa = ya./(Alpha-ya*(Alpha-1));

function ybotop = ybot(x,S,xb)
ybotop = (S+1)/S*x-xb/S;

function ytopop = ytop(x,R,xd)
ytopop = R/(R+1)*x+xd/(R+1);

function backl = predefeq(xeq,yeq,val)
xx = 0:.001:1;
yy = spline(xeq,yeq,xx);
ym = val*ones(1,length(xx));
diffe = abs(yy-ym);
[a b] = min(diffe);
xk = xx(b);
backl = xk;

%System of Nonlinear Equations -- Will be functional in a later version

%User Input Vector
%L = R*D;
%V = D*(R+1);
%Vbar = S*B;
%Lbar = B*(S+1);
%Input_Vector = [F,zf,zfnv,xd,xdnv,xb,xbnv,D,fv,fnv,B,S,Q,R,Rmin,hf,hD,hB,Qc,Qr,L,V,Lbar,Vbar,xie,yie,xi,yi];
%for k = 1:20
%    if isnan(Input_Vector(k)) == 0
%        v(k) = Input_Vector(k);
%    end
%end

%v(1) = F;
%v(2) = zf;
%v(3) = zfnv;
%v(4) = xd;
%v(5) = xdnv;
%v(6) = xb;
%v(7) = xbnv;
%v(8) = D;
%v(9) = fv;
%v(10) = fnv;
%v(11) = B;
%v(12) = S;
%v(13) = Q;
%v(14) = R;
%v(15) = Rmin;
%v(16) = hf;
%v(17) = hD;
%v(18) = hB;
%v(19) = Qc;
%v(20) = Qr;
%v(21) = L;
%v(22) = V;
%v(23) = Lbar;
%v(24) = Vbar;
%v(25) = xie;
%v(26) = yie;
%v(27) = xi;
%v(28) = yi;

%F = @(v) [v(10)-v(11).*v(7)./v(1)./v(3); v(9)-v(8).*v(4)./v(1)./v(2); v(11)+v(24)-v(23);
%    v(21)+v(8)-v(22); v(23)./v(24)-1./v(12)+1; v(22)./v(21)-1./v(14)-1; v(11)+v(8)-v(1);
%    v(11).*v(6)+v(8).*v(4)-v(1).*v(2); v(2)+v(3)-1; v(6)+v(7)-1; v(4)+v(5)-1;
%    v(1).*v(16)+v(19)+v(20)-v(8).*v(17)-v(11).*v(18); v(15)./(1+v(15))-(v(4)-v(26))./(v(4)-v(25));
%    v(12)-1./(((v(28)-v(6))./(v(27)-v(6)))-1); v(14)./(1+v(14))-(v(4)-v(28))./(v(4)-v(27))];

%InitialGuess = [50;.5;.5;.5;.5;.5;.5;25;.75;.75;25;2;.5;3;1.5;1;1;1;1;1;1;1;1;1;.5;.5;.5;.5];
%Options = optimset('Display','iter');
%XY = fsolve(F, InitialGuess, Options);

%fnv - B*xbnv/F/zfnv
%fv - D*xd/F/zf;
%B + Vbar - Lbar;
%L + D - V;
%Lbar/Vbar - 1/S + 1;
%V/L - 1/R - 1;
%B + D - F;
%B*xb + D*xd - F*zf;
%zf + zfnv - 1;
%xb + xbnv - 1;
%xd + xdnv - 1;
%F*hf + Qc + Qr - D*hD - B*hB;
%Rmin/(1+Rmin) - (xd-yie)/(xd-xie);
%S - 1/((yi-xb)/(xi-xb)-1);
%R/(1+R)-(xd-yi)/(xd-xi);

author: Sable

repo: mcbench-benchmarks


text:
function D = genDeltaTransform(nFr, delta_order, type, precision)
if nargin<4
    precision = 'single';
end
if nargin<3
    type = 1;
end

if 0
    D = sparse(nFr, nFr+delta_order*2);
    for t=1 : nFr
        for j=1:delta_order
            D(t,t-j+delta_order) = -j;
            D(t,t+j+delta_order) = j;
        end
    end
else
    D = zeros(nFr,nFr+delta_order*2, precision);
    for j=1:delta_order
        t = 1:nFr;
        idx = (t+j+delta_order-1)*nFr + t;
        D(idx) = j;
        idx = (t-j+delta_order-1)*nFr + t;
        D(idx) = -j;
    end
end

D2 = D(:, delta_order+1:end-delta_order);
if type==1
    D2(:,1) = sum(D(:,1:delta_order+1)')';
    D2(:,end) = sum(D(:,end-delta_order:end)')';
end

D = D2 / sum((1:delta_order).^2) * 0.5;

author: singaxiong

repo: SignalGraph


text:
classdef ShapeTransformer < handle
    %SHAPETRANSFORMER  Base class for shape transformation algorithms
    %
    % See also: cv.ShapeTransformer.ShapeTransformer,
    %  cv.ShapeContextDistanceExtractor, tpaps
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent, Hidden)
        % The regularization parameter for relaxing the exact interpolation
        % requirements of the TPS algorithm.
        % Beta value of the regularization parameter, default 0
        % (only applicable for ThinPlateSplineShapeTransformer)
        RegularizationParameter
        % default true
        % (only applicable for AffineTransformer)
        FullAffine
    end

    %% Constructor/destructor
    methods
        function this = ShapeTransformer(transformerType, varargin)
            %SHAPETRANSFORMER  Constructor
            %
            %     obj = cv.ShapeTransformer(transformerType)
            %     obj = cv.ShapeTransformer(transformerType, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __transformerType__ an algorithm that defines the aligning
            %   transformation. One of:
            %   * __ThinPlateSplineShapeTransformer__ Definition of the
            %     transformation occupied in the paper [Bookstein89].
            %   * __AffineTransformer__ Wrapper class for the OpenCV Affine
            %     Transformation algorithm. See cv.estimateRigidTransform,
            %     cv.warpAffine, and cv.transform
            %
            % ## Options
            % The following are options for the various algorithms:
            %
            % ### `ThinPlateSplineShapeTransformer`
            % * __RegularizationParameter__ The regularization parameter for
            %   relaxing the exact interpolation requirements of the TPS
            %   algorithm. default 0
            %
            % ### `AffineTransformer`
            % * __FullAffine__ see cv.estimateRigidTransform, default true
            %
            % ## References
            % [Bookstein89]:
            % > "Principal Warps: Thin-Plate Splines and Decomposition of
            % > Deformations", by F.L. Bookstein (PAMI 1989)
            %
            % See also: cv.ShapeTransformer.estimateTransformation
            %
            this.id = ShapeTransformer_(0, 'new', transformerType, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.ShapeTransformer
            %
            if isempty(this.id), return; end
            ShapeTransformer_(this.id, 'delete');
        end
    end

    %% ShapeTransformer
    methods
        function estimateTransformation(this, transformingShape, targetShape, matches)
            %ESTIMATETRANSFORMATION  Estimate the transformation parameters of the current transformer algorithm, based on point matches
            %
            %     obj.estimateTransformation(transformingShape, targetShape, matches)
            %
            % ## Input
            % * __transformingShape__ Contour defining first shape. A numeric
            %   Nx2/Nx1x2/1xNx2 array or a cell-array of 2D points
            %   `{[x,y], ...}`
            % * __targetShape__ Contour defining second shape (target). Same
            %   format as `transformingShape`.
            % * __matches__ Standard vector of matches between points.
            %
            % See also: cv.ShapeTransformer.applyTransformation,
            %  cv.ShapeTransformer.warpImage
            %
            ShapeTransformer_(this.id, 'estimateTransformation', transformingShape, targetShape, matches);
        end

        function [cost, output] = applyTransformation(this, input)
            %APPLYTRANSFORMATION  Apply a transformation, given a pre-estimated transformation parameters
            %
            %     [cost, output] = obj.applyTransformation(input)
            %
            % ## Input
            % * __input__ Contour (set of points) to apply the transformation.
            %
            % ## Output
            % * __cost__ Transformation cost.
            % * __output__ Output contour.
            %
            % See also: cv.ShapeTransformer.estimateTransformation
            %
            [cost, output] = ShapeTransformer_(this.id, 'applyTransformation', input);
        end

        function output = warpImage(this, transformingImage, varargin)
            %WARPIMAGE  Apply a transformation, given a pre-estimated transformation parameters, to an Image
            %
            %     output = obj.warpImage(transformingImage)
            %     output = obj.warpImage(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __transformingImage__ Input image.
            %
            % ## Output
            % * __output__ Output image.
            %
            % ## Options
            % * __Interpolation__ Image interpolation method. default 'Linear'
            % * __BorderType__ border style. default 'Constant'
            % * __BorderValue__ border value. default 0
            %
            % See cv.remap or cv.warpAffine for options description.
            %
            % See also: cv.ShapeTransformer.estimateTransformation
            %
            output = ShapeTransformer_(this.id, 'warpImage', transformingImage, varargin{:});
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.ShapeTransformer.empty, cv.ShapeTransformer.load
            %
            ShapeTransformer_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.ShapeTransformer.clear, cv.ShapeTransformer.load
            %
            b = ShapeTransformer_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.ShapeTransformer.load
            %
            ShapeTransformer_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.ShapeTransformer.save
            %
            ShapeTransformer_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.ShapeTransformer.save, cv.ShapeTransformer.load
            %
            name = ShapeTransformer_(this.id, 'getDefaultName');
        end
    end

    %% Getters/Setters
    methods
        function value = get.RegularizationParameter(this)
            value = ShapeTransformer_(this.id, 'get', 'RegularizationParameter');
        end
        function set.RegularizationParameter(this, value)
            ShapeTransformer_(this.id, 'set', 'RegularizationParameter', value);
        end

        function value = get.FullAffine(this)
            value = ShapeTransformer_(this.id, 'get', 'FullAffine');
        end
        function set.FullAffine(this, value)
            ShapeTransformer_(this.id, 'set', 'FullAffine', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sde_install(opt)
%SDE_INSTALL  Add or remove SDETools from Matlab search path and save path.
%   SDE_INSTALL adds the SDETools directory (where this function is located) to
%   the Matlab search path, saves the path, and prints the help for the toolbox.
%
%   SDE_INSTALL('remove') uninstalls SDETools by removing the SDETools directory
%   from the Matlab search path and saving the path. If SDETools is not
%   installed (on the path), a warning is issued.
%
%   See also: PATH, ADDPATH, RMPATH, SAVEPATH

%   Andrew D. Horchler, horchler @ gmail . com, Created 8-12-13
%   Revision: 1.2, 11-16-13

if nargin == 0 || any(strcmp(opt,{'add','install','addpath'}))
    addpath(fileparts(mfilename('fullpath')));
    status = true;
elseif any(strcmp(opt,{'remove','uninstall','rmpath'}))
    rmpath(fileparts(mfilename('fullpath')));
    status = false;
else
    error('SDETools:sde_install:UnknownOption',...
         ['Input argument must be the string ''install'' to install or '...
          '''remove'' to uninstall.']);
end

if savepath
    error('SDETools:sde_install:SavePathError',...
          'Unable to save pathdef.m file.');
end
rehash('toolbox');
clear('sde_install');

if status
    fprintf(1,'\n SDETools installed.\n\n');
    help('SDETools');
else
    fprintf(1,'\n SDETools uninstalled.\n');
end

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools


text:
function [username,password]=GetAuthentication(defaultuser)
%GetAuthentication prompts a username and password from a user and hides the
% password input by *****
%
%   [user,password] = GetAuthentication;
%   [user,password] = GetAuthentication(defaultuser);
%
% arguments:
%   defaultuser - string for default name
%
% results:
%   username - string for the username
%   password - password as a string
%
% Created by Felix Ruhnow, MPI-CBG Dresden
% Version 1.00 - 20th February 2009
%

if nargin==0
    defaultuser='';
end

hAuth.fig = figure('Menubar','none','Units','normalized','Resize','off','NumberTitle','off', ...
                   'Name','Authentication','Position',[0.4 0.4 0.2 0.2],'WindowStyle','normal');

uicontrol('Parent',hAuth.fig,'Style','text','Enable','inactive','Units','normalized','Position',[0 0 1 1], ...
          'FontSize',12);

uicontrol('Parent',hAuth.fig,'Style','text','Enable','inactive','Units','normalized','Position',[0.1 0.8 0.8 0.1], ...
                      'FontSize',12,'String','Username:','HorizontalAlignment','left');

hAuth.eUsername = uicontrol('Parent',hAuth.fig,'Style','edit','Tag','username','Units','normalized','Position',[0.1 0.675 0.8 0.125], ...
                       'FontSize',12,'String',defaultuser,'BackGroundColor','white','HorizontalAlignment','left');

uicontrol('Parent',hAuth.fig,'Style','text','Enable','inactive','Units','normalized','Position',[0.1 0.5 0.8 0.1], ...
          'FontSize',12,'String','Password:','HorizontalAlignment','left');

hAuth.ePassword = uicontrol('Parent',hAuth.fig,'Style','edit','Tag','password','Units','normalized','Position',[0.1 0.375 0.8 0.125], ...
                            'FontSize',12,'String','','BackGroundColor','white','HorizontalAlignment','left');

uicontrol('Parent',hAuth.fig,'Style','pushbutton','Tag','OK','Units','normalized','Position',[0.1 0.05 0.35 0.2], ...
                            'FontSize',12,'String','OK','Callback','uiresume;');

uicontrol('Parent',hAuth.fig,'Style','pushbutton','Tag','Cancel','Units','normalized','Position',[0.55 0.05 0.35 0.2], ...
                            'FontSize',12,'String','Cancel','Callback',@AbortAuthentication);

set(hAuth.fig,'CloseRequestFcn',@AbortAuthentication)
set(hAuth.ePassword,'KeypressFcn',@PasswordKeyPress)

setappdata(0,'hAuth',hAuth);
uicontrol(hAuth.eUsername);
uiwait;

username = get(hAuth.eUsername,'String');
password = get(hAuth.ePassword,'UserData');
delete(hAuth.fig);

function PasswordKeyPress(hObject,event)
hAuth = getappdata(0,'hAuth');
password = get(hAuth.ePassword,'UserData');
switch event.Key
   case 'backspace'
      password = password(1:end-1);
   case 'return'
      uiresume;
      return;
   otherwise
      password = [password event.Character];
end
set(hAuth.ePassword,'UserData',password)
set(hAuth.ePassword,'String',char('*'*sign(password)))

function AbortAuthentication(hObject,event)
hAuth = getappdata(0,'hAuth');
set(hAuth.eUsername,'String','');
set(hAuth.ePassword,'UserData','');
uiresume;


text:
function [ bgU, status, OPTS ] = bgtraining( videopath, OPTIONS, max_cycles, subsampling, training_size, FPS_ONLY , TRAIN_FRAME)
%BGTRAINING Summary of this function goes here
%   Detailed explanation goes here
thresh = 0.1;

FILE_EXT        = '.bmp';
OPTS            = struct(); % initial a empty struct for OPTS
status.init     = 0;   % status of grasta at each iteration
U_hat           = zeros(1);

%
Dir_lists = dir(videopath);
Video_Length = min(training_size,length(Dir_lists));

if TRAIN_FRAME == 1, % Random selection
    p = randperm(length(Dir_lists)); % length(Dir_lists)-1500 for hall
    frames_idx = p(1:Video_Length);
else                % Select all frames
    frames_idx = 1:Video_Length;
end

if ~FPS_ONLY,
    h_fg = subplot(2,2,1);set(gca,'nextplot','replacechildren');title('GRASTA-Foreground');
    h_training_bg = subplot(2,2,2);set(gca,'nextplot','replacechildren');title('Training Background');
    h_img = subplot(2,2,3);set(gca,'nextplot','replacechildren');title('Training video');
end

for outiter = 1:max_cycles,
    video_order = randperm(Video_Length);  %  1:Video_Length; %
    t_start = tic;
    frame_count  = 0;
    for i=1:Video_Length,
        if Dir_lists(frames_idx(video_order(i))).isdir,
            continue;
        end
        if isempty(strfind(Dir_lists(frames_idx(video_order(i))).name,FILE_EXT)) ,
            continue;
        end

        frame_count = frame_count+1;
        fname = [videopath Dir_lists(frames_idx(video_order(i))).name];

        % prepare the image
        I = imread(fname);
        I = double(rgb2gray(I));

        if frame_count==1,
            [rows,cols]     = size(I);
            VW_ROWS         = rows;
            VW_COLS         = cols; %ceil(cols * VW_RATIO);
            DIM             = VW_ROWS * VW_COLS;
            OPTIONS.DIM_M   = DIM; % video ambient dimension
        end

        I = I/max(max(I));

        % random subsampling the frame I
        M = round(subsampling * DIM);
        p = randperm(DIM);
        idx = p(1:M)';

        I_Omega = I(idx);

        [U_hat, status, OPTS] = grasta_stream(I_Omega, idx, U_hat, status, OPTIONS, OPTS);

        if mod(i,1) == 0 && ~FPS_ONLY,
            bg_img = reshape(U_hat * status.w * status.SCALE, VW_ROWS,VW_COLS);
            axes(h_training_bg); imagesc(bg_img);colormap gray;axis off;axis ij ;

            o_img = reshape( I ,VW_ROWS,VW_COLS );
            axes(h_img); imagesc(o_img);colormap gray;axis off;axis ij ;

            s_hat = I(:) - U_hat * status.w * status.SCALE;
            s_hat = fg_thresholding(s_hat,thresh);

            s_img = reshape(s_hat,VW_ROWS,VW_COLS);
            axes(h_fg); imagesc(s_img);colormap gray;axis off;axis ij ;

        end
    end
    t_end = toc(t_start);
    fprintf('Training %d/%d: %.2f seconds, %.2f fps, grasta_t %.2e \n',...
        outiter, max_cycles,t_end, frame_count/t_end,status.grasta_t);
end

bgU = U_hat;
end

author: andrewssobral

repo: lrslibrary


text:
function r = mtimes(p,q)
% MEAS/MTIMES  Implement p * q for meas.

% make a meas called r
r = meas();

% give it the right entries
r.value = p.value * q.value;
r.error = r.value * sqrt((p.error/p.value)^2 + (q.error/q.value)^2);


text:
% trial using t-test

range_fish = 8:11;
numFish = length(range_fish);
TF_fishrange = cell(1,numFish);
% also save clusID of matching pairs???
Pairs_AllClusAllFish = []; % each row: fish#1, clusID1, fish#2, clusID2

for i_fishnum = 1:numFish,
    % choose reference fish
    i_fish_ref = range_fish(i_fishnum);
    U_ref = unique(FC{i_fish_ref}.gIX);
    numClus_ref = length(U_ref);
    TF = zeros(numClus_ref,numFish-1);

    % cycle through clusters in reference fish
    for i_ref_clus = 1:numClus_ref,
        refclusID = U_ref(i_ref_clus);
        IX_ref = find(FC{i_fish_ref}.gIX == refclusID);
        Ref_xyz = FC{i_fish_ref}.xyz_norm(IX_ref,:);

        % for current cluster, compare with closest cluster in each other fish
        for i_testnum = 1:numFish,
            i_fish_test = range_fish(i_testnum);
            if i_fish_test ~= i_fish_ref,
                U = unique(FC{i_fish_test}.gIX);
                numClus = length(U);

                % find testclusID that is closest to current ref-clus in anat-space
                refcoord = mean(Ref.xyz_n,1);
                testcoords = zeros(numClus,3);
                Dist = zeros(numClus,1);
                for i = 1:numClus,
                    IX = find(FC{i_fish_test}.gIX==U(i));
                    testcoords(i,:) = mean(FC{i_fish_test}.xyz_norm(IX,:),1);
                    Dist(i) = pdist([refcoord;testcoords(i,:)]);
                end
                [~,ix] = min(Dist);
                testclusID = U(ix);

                % perform t-test between this pair of clusters
                IX = find(FC{i_fish_test}.gIX==testclusID);
                Test_xyz = FC{i_fish_test}.xyz_norm(IX,:);
                [h,p] = ttest2(Ref_xyz,Test_xyz,'Vartype','unequal');
                if max(h) == 0,
                    % clusters are the same, collect
                    TF(i_ref_clus,i_testnum) = 1;
                    if i_fish_ref<i_fish_test,
                        Pairs_AllClusAllFish = [Pairs_AllClusAllFish; ...
                            [i_fish_ref, refclusID, i_fish_test, testclusID]];
                    end
                else
                    TF(i_ref_clus,i_testnum) = 1-sum(h)/3;
                end
            else % don't compare within same fish
                TF(i_ref_clus,i_testnum) = NaN;
            end
        end
    end
    TF_fishrange{i_fishnum} = TF;
end

author: xiuyechen

repo: FishExplorer


text:
function plasma_matrix_test ( n )

%*****************************************************************************80
%
%% PLASMA_MATRIX_TEST tests PLASMA_MATRIX.
%
%  Discussion:
%
%    This program shows how a MATLAB sparse matrix (and possibly a right
%    hand side vector) can be stored into a Harwell-Boeing file,
%    and later retrieved.
%
%    Harwell-Boeing files are useful as a means of storing sparse matrices,
%    especially when data is created with one program and needed by another.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 January 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PLASMA_MATRIX_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Create a large sparse plasma matrix and right hand side.\n' );
  fprintf ( 1, '  Store the matrix and right hand side in a Harwell-Boeing file.\n' );
  fprintf ( 1, '  Then retrieve the information.\n' );

  plasma_matrix_test01 ( 5 );
  plasma_matrix_test01 ( 100 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PLASMA_MATRIX_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function msh = msh_2d (breaks, qn, qw, geo, varargin)

warning ('geopdes:obsolete', 'The class MSH_2D is obsolete. Using MSH_CARTESIAN instead')

msh = msh_cartesian (breaks, qn, qw, geo, varargin{:});

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
classdef EventPlottingMethodEnum < matlab.mixin.SetGet
    %EventPlottingMethodEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        PlotContinuous('Plot Continuous');
        SkipFirstState('Skip First State');
        DoNotPlot('Do Not Plot');
    end

    properties
        name
    end

    methods
        function obj = EventPlottingMethodEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('EventPlottingMethodEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('EventPlottingMethodEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('EventPlottingMethodEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [BIscores]=mvg_superpixelBBintegralScore(superpixels, windows, config)

%% Default settings
if nargin<3
    config.probDist='Gaussian'; % 'Dist' 'Gaussian' 'DistGaussian'
    config.GaussSigma=3;
    config.exponent=2;
    config.dilate=3;
    config.verbose=4;
end

%% Find bounding boxes and make boundary mask image
lbs=unique(superpixels);
numSpix=length(lbs);
spixBox=zeros(numSpix,4);
for i=1:numSpix
    [rw,cl]=find(superpixels==lbs(i));
    spixBox(i,:)=[min(cl),min(rw),max(cl),max(rw)];
end

%% Make boundary mask
boundaryMask=zeros(size(superpixels));
for i=1:numSpix
    boundaryMask(spixBox(i,2):spixBox(i,4),spixBox(i,1))=1;
    boundaryMask(spixBox(i,2):spixBox(i,4),spixBox(i,3))=1;
    boundaryMask(spixBox(i,2),spixBox(i,1):spixBox(i,3))=1;
    boundaryMask(spixBox(i,4),spixBox(i,1):spixBox(i,3))=1;
end

%% Dilate
if config.dilate>0
    structEl=strel('disk',config.dilate);
    boundaryMask=imdilate(boundaryMask,structEl);
end

%% Make distance transform based probability
switch config.probDist
    case 'Dist'
        probFun=1-bwdist(boundaryMask);
        probFun=(probFun-min(probFun(:)))/max((probFun(:)-min(probFun(:))));
        probFun=probFun.^config.exponent;
    case 'Gaussian'
        GaussianFun=makeGaussian_(config.GaussSigma*[1 1]);
        probFun=conv2(double(boundaryMask),GaussianFun,'same');
    case 'DistGaussian'
        probFun=1-bwdist(boundaryMask);
        probFun=(probFun-min(probFun(:)))/max((probFun(:)-min(probFun(:))));
        probFun=probFun.^config.exponent;
        GaussianFun=makeGaussian(config.GaussSigma*[1 1]);
        probFun=conv2(double(probFun),GaussianFun,'same');
    otherwise
        error('Unknown method');
end

%% Integrate over window borders to get BI score
BIscores=integrateOverWindowBorder_(probFun,windows);

%%%%%%%%%%%%%%%%%%%%%%%%
% Additional functions %
%%%%%%%%%%%%%%%%%%%%%%%%

%%% Integrate over the windows border %%%
function [windowIntegral,winMeasure]=integrateOverWindowBorder_(probFun,windows)

%% Initialize
if size(windows,2)>4
    windows=windows(:,1:4);
end
numWindow=size(windows,1);
intLineImg=zeros(size(probFun,1),size(probFun,2),2);
windowIntegral=zeros(numWindow,1);

%% Form line integral image
intLineImg(:,:,1)=cumsum(probFun,1); % Row sums
intLineImg(:,:,2)=cumsum(probFun,2); % Column sums
intLineImg=[zeros(size(intLineImg,1),1,2),intLineImg]; % Add zero column
intLineImg=[zeros(1,size(intLineImg,2),2);intLineImg]; % Add zero row

%% If given windows are in normalized coordinates, denormalize them
if max(windows(:))<1.000001
    imgCol=size(probFun,2);
    imgRow=size(probFun,1);
    windows(:,[1,3])=windows(:,[1,3])*(imgCol-1)+1;
    windows(:,[2,4])=windows(:,[2,4])*(imgRow-1)+1;
end

%% Round windows to integer coordinates
windows=max(round(windows(:,1:4)),1);

%% Loop over windows and compute bounding box scores
for i=1:numWindow
    % Get critical points
    Xmin=windows(i,1);
    Ymin=windows(i,2);
    Xmax=windows(i,3);
    Ymax=windows(i,4);

    % Compute sums (start each one pixel further to avoid taking corners twice)
    topSum=intLineImg(Ymin+1,Xmax+1,2)-intLineImg(Ymin+1,Xmin+1,2);
    bottomSum=intLineImg(Ymax+1,Xmax+1,2)-intLineImg(Ymax+1,Xmin+1,2);
    leftSum=intLineImg(Ymax+1,Xmin+1,1)-intLineImg(Ymin+1,Xmin+1,1);
    rightSum=intLineImg(Ymax+1,Xmax+1,1)-intLineImg(Ymin+1,Xmax+1,1);
    %% For reference the full sums, which count corner points twice
    %topSum=intLineImg(Ymin+1,Xmax+1,2)-intLineImg(Ymin+1,Xmin,2);
    %bottomSum=intLineImg(Ymax+1,Xmax+1,2)-intLineImg(Ymax+1,Xmin,2);
    %leftSum=intLineImg(Ymax+1,Xmin+1,1)-intLineImg(Ymin,Xmin+1,1);
    %rightSum=intLineImg(Ymax+1,Xmax+1,1)-intLineImg(Ymin,Xmax+1,1);

    % Assing sum over bounding box to window score
    windowIntegral(i)=topSum+bottomSum+leftSum+rightSum;

end

%% If two output arguments are required, return also box size and edge length
if nargin>1
    boxWidth=windows(:,3)-windows(:,1)+1;
    boxHeight=windows(:,4)-windows(:,2)+1;

    winMeasure.Perimeter=2*boxWidth+2*boxHeight-4; % Need to subtract extra corners (that's why -4)
    winMeasure.Area=boxWidth.*boxHeight;
end

%%% Make Gaussian funtion %%%
function [GaussianFun]=makeGaussian_(GaussianSigma,WindowRadius)

%% Return trivial case
if min(GaussianSigma)<eps
    GaussianFun=1;
    return;
end

%% Default settings
if nargin<2
    %WindowRadius=max(round(2*GaussianSigma),1);
    WindowRadius=round(2*GaussianSigma);
end

%% Initialize
dim=length(GaussianSigma); % How many dimensions in output Gaussian (3D is max and ordering is row, column, and third dimension).
GaussianSigma=2*GaussianSigma.^2; %Turn sigma from standard deviation to variance (2* is to include 1/(2*sigma^2) already here).
% If only one value for size is given and more dimensions are required, use same size for all dimensions.
if length(WindowRadius)==1 && dim>1
    WindowRadius=WindowRadius(1)*ones(1,dim);
end

%% Generate Gaussian function to the required dimensions
if dim==1
    d1=-WindowRadius(1):WindowRadius(1); % spatial coordinates
    GaussianFun=exp(-((d1.^2)/GaussianSigma(1))); % Gaussian values

elseif dim==2
    [d2,d1]=meshgrid(-WindowRadius(2):WindowRadius(2),-WindowRadius(1):WindowRadius(1)); % spatial coordinates
    GaussianFun=exp(-((d1.^2)/GaussianSigma(1)+(d2.^2)/GaussianSigma(2))); % Gaussian values

elseif dim==3
    [d2,d1,d3]=meshgrid(-WindowRadius(2):WindowRadius(2),-WindowRadius(1):WindowRadius(1),-WindowRadius(3):WindowRadius(3)); % spatial coordinates
    GaussianFun=exp(-((d1.^2)/GaussianSigma(1)+(d2.^2)/GaussianSigma(2)+(d3.^2)/GaussianSigma(3))); % Gaussian values

else
    error('Not implemented');
end

%% Normalize data to have sum equal to fun.
GaussianFun=GaussianFun/sum(GaussianFun(:));

author: Cloud-CV

repo: object-proposals