text:
function varargout = optionpricegui2(varargin)
% OPTIONPRICEGUI2 M-file for optionpricegui2.fig
%      OPTIONPRICEGUI2, by itself, creates a new OPTIONPRICEGUI2 or raises the existing
%      singleton*.
%
%      H = OPTIONPRICEGUI2 returns the handle to a new OPTIONPRICEGUI2 or the handle to
%      the existing singleton*.
%
%      OPTIONPRICEGUI2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in OPTIONPRICEGUI2.M with the given input arguments.
%
%      OPTIONPRICEGUI2('Property','Value',...) creates a new OPTIONPRICEGUI2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before optionpricegui2_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to optionpricegui2_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

% Edit the above text to modify the response to help optionpricegui2

% Last Modified by GUIDE v2.5 14-Aug-2008 14:18:46

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

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

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

set(handles.editStart, 'String', datestr(today, 2))
set(handles.editEnd, 'String', datestr(today + 91, 2))

handles = runSimulation(handles);

% This is an awkward place to put the rotate3d command.  If you place it in
% the runSimulation function, however, it tends to toggle rotate3d on and
% off after each simulation.  Here, it gets set to "on" once and stays that
% way.
rotate3d(handles.surfAxes)

% Update handles structure
guidata(hObject, handles);
end

% UIWAIT makes optionpricegui2 wait for user response (see UIRESUME)
% uiwait(handles.mainFigure);

% --- Outputs from this function are returned to the command line.
function varargout = optionpricegui2_OutputFcn(hObject, eventdata, handles) %#ok
% 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;
end

function editSpot_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editSpot (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 editSpot as text
%        str2double(get(hObject,'String')) returns contents of editSpot as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editStrike_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editStrike (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 editStrike as text
%        str2double(get(hObject,'String')) returns contents of editStrike as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editRFR_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editRFR (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 editRFR as text
%        str2double(get(hObject,'String')) returns contents of editRFR as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editStart_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editStart (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 editStart as text
%        str2double(get(hObject,'String')) returns contents of editStart as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editEnd_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editEnd (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 editEnd as text
%        str2double(get(hObject,'String')) returns contents of editEnd as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editVol_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editVol (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 editVol as text
%        str2double(get(hObject,'String')) returns contents of editVol as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editYield_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editYield (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 editYield as text
%        str2double(get(hObject,'String')) returns contents of editYield as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

end

function editIter_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editIter (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 editIter as text
%        str2double(get(hObject,'String')) returns contents of editIter as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

end

function editValue_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editValue (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 editValue as text
%        str2double(get(hObject,'String')) returns contents of editValue as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

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

% Hints: contents = get(hObject,'String') returns popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1

% Some options only require a single strike price; others, two.  The
% butterfly option requires one strike price and a spread.
contents = get(hObject,'String');
switch contents{get(hObject,'Value')}
    case {'Call' 'Put' 'Straddle'}
        handles = oneStrike(handles);
    case {'Strangle' 'Bull Spread' 'Bear Spread'}
        handles = twoStrikes(handles);
    case {'Butterfly'}
        handles = butterfly(handles);
end
guidata(hObject, handles);
end
% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)%#ok
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editStrike1_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editStrike1 (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 editStrike1 as text
%        str2double(get(hObject,'String')) returns contents of editStrike1 as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editStrike2_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editStrike2 (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 editStrike2 as text
%        str2double(get(hObject,'String')) returns contents of editStrike2 as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

function editSpread_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to editSpread (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 editSpread as text
%        str2double(get(hObject,'String')) returns contents of editSpread as a double
end

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

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
end

% --- Executes on button press in buttonRun.
function buttonRun_Callback(hObject, eventdata, handles)%#ok
% hObject    handle to buttonRun (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles = runSimulation(handles);
guidata(hObject, handles);
end

% --- Helper functions follow ---
% --- oneStrike, twoStrikes: toggle between displaying one and two strike
% prices in the UI panel ---
function handles = oneStrike(handles)
set(handles.text3, 'String', 'Strike Price')
set(handles.editStrike, 'Visible', 'on')
set(handles.editStrike1, 'Visible', 'off')
set(handles.editStrike2, 'Visible', 'off')
set(handles.editSpread, 'Visible', 'off')
set(handles.text11, 'Visible', 'off')
end

function handles = twoStrikes(handles)
set(handles.text3, 'String', 'Strike Prices')
set(handles.editStrike, 'Visible', 'off')
set(handles.editStrike1, 'Visible', 'on')
set(handles.editStrike2, 'Visible', 'on')
set(handles.editSpread, 'Visible', 'off')
set(handles.text11, 'Visible', 'off')
end

function handles = butterfly(handles)
set(handles.text3, 'String', 'Strike Price')
set(handles.editStrike, 'Visible', 'on')
set(handles.editStrike1, 'Visible', 'off')
set(handles.editStrike2, 'Visible', 'off')
set(handles.editSpread, 'Visible', 'on')
set(handles.text11, 'Visible', 'on')
end

%--- runSimulation: the main control program for the code ---
function handles = runSimulation(handles)
% This function has four steps:
% Step 1: Collect inputs and find option value
% Step 2: Plot the option value surface
% Step 3: Run a Monte Carlo simulation
% Step 4: Plot the Monte Carlo simulation
%
% In addition, there are two nested functions that allow you to incorporate
% a WindowButtonMotionFcn and a WindowButtonDownFcn on the Monte Carlo
% Paths axes.

% Step 1: Collect inputs and find option value
try
    spot = str2double(get(handles.editSpot, 'String'));
    rate = str2double(get(handles.editRFR, 'String'));
    vol  = str2double(get(handles.editVol, 'String'));
    yield= str2double(get(handles.editYield, 'String'));
    type = get(handles.popupmenu1,'String');
    type = type{get(handles.popupmenu1,'Value')};

    % No harm in getting Strike2 and Spread, even if they aren't used.
    strike2 = str2double(get(handles.editStrike2, 'String'));
    spread  = str2double(get(handles.editSpread, 'String'));

    startDate = get(handles.editStart, 'String');
    endDate = get(handles.editEnd, 'String');
    % I only include this extra TRY/CATCH loop in order to provide a better
    % error message if the invalid dates are applied: YEARFRAC itself will
    % do all of the appropriate error checking.
    try
        time = yearfrac(startDate, endDate);
    catch %#ok<CTCH>
        error('optionPriceGui:InvalidDates', ...
            'Dates must be in a valid format')
    end

    iter = str2double(get(handles.editIter, 'String'));
    if int32(iter) ~= iter || iter <= 0
        error('optionPriceGui:InvalidIterations',...
            'Iteration count must be a positive integer.')
    end

    switch type
        case {'Call' 'Put' 'Straddle' 'Butterfly'}
            strike = str2double(get(handles.editStrike, 'String'));
        case {'Strangle' 'Bull Spread' 'Bear Spread'}
            strike = str2double(get(handles.editStrike1, 'String'));
    end

    value = optPriceVal(type, spot, strike, rate, time, vol, yield, strike2, spread);
    set(handles.editValue, 'String', num2str(value,'%4.3f'))

catch ME
    errordlg(ME.message, 'Invalid Input', 'modal')
    set(handles.editValue, 'String', 'error')
    return
end

% Step 2: Now that we have all the data we need, let's plot the option
% value surface.

% A 51x51 surface plot should look nice and not be too taxing on the
% graphics processor, but you can adjust that here.
xPoints = 51;
yPoints = 51;

T = linspace(0, time, xPoints);

% The spot price axis span is a little tricky to determine in advance.  We
% will go with +/- 10% from the strike price in the case of a put, call, or
% straddle option, +/- 50% of the strike prices' span in the case of a
% strangle, bull spread, or bear spread option, and +/- 40% of the spread
% from a butterfly option.  Of course, we will truncate if these go below a
% price of 0.
switch type
    case {'Call' 'Put' 'Straddle'}
        lPrice = 0.9 * strike;
        uPrice = 1.1 * strike;
    case {'Strangle' 'Bull Spread' 'Bear Spread'}
        span = abs(strike2 - strike);
        if strike2 > strike
            lPrice = strike - 0.5 * span;
            uPrice = strike2+ 0.5 * span;
        else
            lPrice = strike2- 0.5 * span;
            uPrice = strike + 0.5 * span;
        end
    case 'Butterfly'
        lPrice = (1 - 1.4 * spread) * strike;
        uPrice = (1 + 1.4 * spread) * strike;
end
lPrice(lPrice < 0) = 0;
S = linspace(lPrice, uPrice, yPoints);
[S,T] = meshgrid(S,T);

% Get the option value for the surface
V = optPriceVal(type, S, strike, rate, T, vol, yield, strike2, spread);

% Plot the 3D surface of changing option value with spot price and time
surf(handles.surfAxes, S, T, V);
box(handles.surfAxes, 'on')
shading(handles.surfAxes, 'interp')
axis(handles.surfAxes, 'tight')
xlabel(handles.surfAxes, 'Underlying ($)','fontweight','bold')
ylabel(handles.surfAxes, 'Time to Expiry (Years)','fontweight','bold')
zlabel(handles.surfAxes, 'Option Value ($)','fontweight','bold')

% Step 3: Let's now run a Monte Carlo simulation
% We will run a single step of the simulation per day-- this is fairly
% coarse-grained, but it gets the message across.

% Assume geometric Brownian motion and apply Ito's lemma to simplify the
% Monte Carlo equation. (see: eqn T6.9 of Kevin Dowd's _An Introduction to
% Market Risk Measurement_, p. 224)
numDays = daysact(startDate, endDate);
TT = sqrt(time/numDays) * ones(numDays, iter);
rNum = randn(numDays, iter);

BT = TT.*rNum;
% And we pad this by an initial row of zeros to handle the initial point at
% day 0.
TT = [zeros(1, iter); TT];
BT = [zeros(1, iter); BT];

ST = spot * exp(cumsum((rate - 0.5 * vol^2) * TT + vol * BT));

% Step 4: Plot the Monte Carlo simulation
% Create line/bar object for histogram of monte-carlo price distributions
set(handles.axesHist,'yaxislocation','right','xaxislocation','top',...
    'xtick',[],'xlimmode','auto')

% Only pick at most 1000 lines from the Monte-Carlo simulation for plotting
% purposes; all data is still used for the underlying math, however.
if size(ST,2)>1000
    ind = randperm(size(ST,2))';
    STplot = ST(:,ind(1:1000));
else
    STplot = ST;
end

% Plot the monte-carlo underlying price paths
handles.pricePaths = plot(handles.axesMonte, (0:numDays), STplot,...
    'color', [.7,.7,.7]);

% Create the tracking line.  Rather, create an object to hold the tracking
% line as it is generated by the windowbuttonmotion function.
handles.trackLine = line('XData', NaN, 'YData', NaN,...
    'Parent', handles.axesMonte, 'erasemode', 'xor', 'Color', 'b', 'LineWidth', 2);
handles.trackText = text('parent', handles.axesMonte, 'erasemode','xor','edgecolor','k',...
    'background',[1,1,.87],'verticalalign','bot','fontsize',14, 'visible', 'off');
% handles.pdfLine = line('XData', NaN, 'YData', NaN,...
%     'Parent', handles.axesHist, 'erasemode', 'xor', 'Color', 'r');

% Plot the mean of the underlying monte-carlo prices paths.
STmean = mean(ST,2);
handles.meanLine = line('XData', (0:numDays), 'YData', STmean, ...
    'Parent', handles.axesMonte, 'Marker', '.' ,'MarkerSize', 10, 'Color', 'r');

% Annotate the monte-carlo plots
xlabel(handles.axesMonte, 'Time (Days)','fontweight','bold')
ylabel(handles.axesMonte, 'Underlying ($)','fontweight','bold')
axis(handles.axesMonte, 'tight')
grid(handles.axesMonte, 'on')
legend(handles.axesMonte, ...
    {'$$S_{t+dt} = S_te^{(r-\frac{{\sigma^2}}{2})dt+\sigma\epsilon\sqrt{dt}}$$'}, ...
    'fontsize',14,'interpreter','latex', 'Location', 'NorthWest');

set(handles.mainFigure,'windowbuttonmotionfcn',@wbm)

% Nested function: WBM for the window button motion function
    function wbm(src, event)
        drawnow
        ylimit = get(handles.axesMonte, 'ylim');
        xlimit = get(handles.axesMonte, 'xlim');
        pos =  get(handles.axesMonte, 'currentpoint');
        x = pos(1,1);
        y = pos(1,2);
        % check if the mouse is out of bounds of the axes
        if x < xlimit(1) || x > xlimit(2) || y < ylimit(1) || y>ylimit(2)
            return
        end

        % create a histogram of the underlying price paths at the current point
        binsC = linspace(ylimit(1), ylimit(2)+eps(ylimit(2)), 21);
        n = histc(ST(round(x)+1, :), binsC);
        n = n./max(n);
        barh(handles.axesHist, binsC, n, 'histc')
        set(handles.axesHist, 'ylim', ylimit)
        set(handles.axesHist,'YTick',[])
        set(handles.axesHist,'XTick',[])
        set(handles.trackLine,'xdata',[x,x],'ydata',ylim)

        % Create a maximum likelihood estimate of the histogram, assuming a
        % lognormal distribution.
        params = lognfit(ST(round(x)+1, :));
        pdfPts = ylimit(1) : ylimit(2);
        pdfVals = lognpdf(pdfPts, params(1), params(2));
        line('Parent', handles.axesHist, ...
            'XData', pdfVals./max(pdfVals), 'YData', pdfPts,...
            'erasemode', 'xor', 'color', 'b')

        % get the current datestr and mean value of the mouse position
        currentDate = datestr(round(x + datenum(startDate)), 2);
        currentMean = STmean(round(x)+1);

        str = sprintf(['Date: ' currentDate '\nAvg Price: $%.2f'], currentMean);
        set(handles.trackText,...
            'position',[x + diff(xlimit)*.02, y + diff(ylimit)*.02, 0],...
            'string', str, 'visible', 'on')

    end
end


text:
classdef ZXH_CF5 < PROBLEM
% <multi/many> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Zhou, Xiang, and He

%------------------------------- Reference --------------------------------
% Y. Zhou, Y. Xiang, and X. He, Constrained multiobjective optimization:
% Test problem construction and performance evaluations, IEEE Transactions
% on Evolutionary Computation, 2021, 25(1): 172-186.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        k;  % Number of constrained variables
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M+10;  end
            obj.lower    = zeros(1,obj.D) + 1e-10;
            obj.upper    = ones(1,obj.D)  - 1e-10;
            obj.encoding = ones(1,obj.D);
            if obj.M <= 3
                obj.k = obj.M - 1;
            elseif obj.M > 3 && obj.M <= 8
                obj.k = floor(obj.M/2);
            else
                obj.k = 3;
            end
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            PopDec = varargin{1};
            OptX   = 0.2;
            [N,D]  = size(PopDec);
            M      = obj.M;
            % Step 1: Compute cumsum
            Sx = cumsum(PopDec(:,1:M).^2,2,'reverse');
            % Step 2: Compute theta
            THETA = 2/pi*atan(sqrt(Sx(:,2:end))./PopDec(:,1:M-1));
            % Step 3: Calculate Rosenbrock function
            h = sum(100*((PopDec(:,M+1:end-1)-OptX).^2-(PopDec(:,M+2:end)-OptX)).^2+(PopDec(:,M+1:end-1)-OptX).^2,2);
            % Step 4: Compute T_
            T = (1 - Sx(:,1)).^2 + h;
            % Step 5: Objectives (convex)
            G      = 1-[ones(N,1) cumprod(sin(pi/2*THETA),2)] .* [cos(pi/2*THETA) ones(N,1)];
            PopObj = G .* repmat((1+T),1,M);
            % Step 6: Constraints
            PopCon(:,1) = Sx(:,1) + h - 1;
            for i = 1 : obj.k
                PopCon(:,i+1) = max(1/4-THETA(:,i),THETA(:,i)-3/4);
            end
            Population = SOLUTION(varargin{1},PopObj,PopCon,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
            R = 1 - R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
            T = zeros(size(R));
            for i = obj.M-1 : -1 : 1
                T(:,i) = atan((1-R(:,i+1))./(1-R(:,i))./cos(T(:,i+1)));
            end
            THETA = T(:,1:obj.k)*2/pi;
            Valid = all(THETA>=1/4&THETA<=3/4,2);
            R     = R(Valid,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a  = linspace(0,pi/2,30)';
                R  = {1-sin(a)*cos(a'),1-sin(a)*sin(a'),1-cos(a)*ones(size(a'))};
                T2 = atan((1-R{3})./(1-R{2}));
                T1 = asin((1-R{2})./cos(T2));
                THETA = cat(3,T1,T2)*2/pi;
                Valid = all(THETA>=1/4&THETA<=3/4,3);
                R{1}(~Valid) = nan;
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function metricsout = assess_sign_flips(netmats,flips)
% function for comparing different sign flip solutions
if min(unique(flips))==0
    flips = 1-2*flips;
end
[N,nch] = size(flips);
nlags = size(netmats,1)/nch;
for iSj=1:N
    flips_to_do = repelem(flips(iSj,:),1,nlags);
    flipped_netmats(iSj,:,:) = netmats(:,:,iSj).*(flips_to_do'*flips_to_do);
end
offdiagblocks = triu(ones(nlags*nch)-repelem(eye(nch),nlags,nlags));
values = flipped_netmats(:,logical(offdiagblocks));

metricsout = [];
metricsout.abscorr = mean(abs(sum(values)));
metricsout.corr = corr(values');
metricsout.corrmean = sum(sum(triu(metricsout.corr,1)))./((N.^2-N)/2);
end

author: OHBA-analysis

repo: HMM-MAR


text:
function w = clean(w)
% CLEAN Clean up waveform object(s)
%   w = clean(w) will detrend waveforms (in a smart way, aware of NaN
%   values marking missing values), then use fillgaps to mark bounded NaNs
%   with linear-interpolated values, and also mark bounded NaNs with zeroes
%   (now okay, since trend has been removed, so no weird startup effects).
%   It then removes non-linear trends using a 20-s (0.05 Hz) highpass
%   filter.

% Glenn Thompson March 28, 2017

    for c=1:numel(w)

        if ~isempty(w(c))

            % remove spikes of length 1
            w(c) = medfilt1(w(c), 3);

            % smart detrend
            w(c) = detrend(w(c));

            % % fill gaps to get rid of NaNs marking missing values, so we can filter
            w(c) = fillgaps(w(c), 'interp');

            % highpass data at 20s - to remove non-linear trends
            f = filterobject('h',0.05,2);
            w(c) = filtfilt(f,w(c));

        end
    end

end

author: geoscience-community-codes

repo: GISMO


text:
function Wrot = whiteningFromCovariance(CC)
% takes as input the matrix CC of channel pairwise correlations
% outputs a symmetric rotation matrix (also Nchan by Nchan) that rotates
% the data onto uncorrelated, unit-norm axes

[E, D] 	= svd(CC); % covariance eigendecomposition (same as svd for positive-definite matrix)
D       = diag(D); % take the non-zero values from the diagonal
eps 	= 1e-6;
Wrot 	= E * diag(1./(D + eps).^.5) * E'; % this is the symmetric whitening matrix (ZCA transform)

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function xs = searchGauss(x,gpstruct,LB,UB,optimState,options)
%SEARCHGAUSS Multivariate normal random search.

if nargin == 0
    xs = 'gauss';
    return;
end

MeshSize = optimState.meshsize;
SearchFactor = optimState.searchfactor;

nvars = length(x);

% Diagonal covariance matrix for random search
sigma = 0.5*MeshSize*SearchFactor*eye(nvars);

% Random draw from multivariate normal
xs = bsxfun(@plus, x, randn(options.Nsearch,nvars)*sigma);

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function f = tan(f)
%TAN   Tangent of a CHEBFUN3T object.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty(f) )
    return
end

op = @(x,y,z) tan(feval(f, x, y, z));     % Resample
f = chebfun3t(op, f.domain);              % Call constructor.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
x = uiuc_sample;
x = x(1:256, 1:256);
filt_opt.J = 6;
Wop = wavelet_factory_2d(size(x), filt_opt);
%%
Sx = scat(x, Wop);
%%
scatter_meta(Sx{3}.meta, 'j', 1, 'j', 2);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function env_testslave(varargin)
% Run as a test slave: run tests for the toolbox whenever it has changed.
%
% In:
%   ScanInterval : the interval at which the toolbox is scanned for changes, in seconds
%                  (default: 30)
%
%   WaitPeriod : the period that must have passed between the last change and a re-test
%                (default: 120)
%
%   ScanDirectories : the directories and files to scan for changes until re-tests are triggered
%                     (default: {'bcilab:/code', 'bcilab:/*.m'})
%
%   TestDirectory : the directory where the test cases are located
%                   (default: 'bcilab:/resources/testcases')
%
%   IssueTest : function to be called to issue a test run (default: @env_test_bcilab)
%               (assumed to return data related to the test summary)
%
%   PostTest : function to be called after a test run (default: [])
%              (called with the results of IssueTest)
%
%   LogFile : name/location of the logfile, if any
%             (default: 'bcilab:/build/buildslave.log')
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-03-17

% read options
o = arg_define(varargin, ...
    arg({'scan_interval','ScanInterval'},15,[0.1 15 3600 1000000],'Interval between scans for changes. In seconds.'), ...
    arg({'wait_period','WaitPeriod'},120,[0 60 300 3600],'Wait this long for further changes until rebuild.'), ...
    arg({'proj_dir','ProjectDirectory'},'bcilab:/build',[],'Name/location of the test project.'), ...
    arg({'scan_dirs','ScanDirectories'},{'bcilab:/code','bcilab:/*.m'},[],'Directories/files to scan for changes.'), ...
    arg({'issue_test','IssueTest'},'@env_test_bcilab',[],'Function called to run new tests. Assumed to return data related to the test summary.','type','expression'), ...
    arg({'post_test','PostTest'},[],[],'Function called after successful tests.','type','expression'), ...
    arg({'log','LogFile'},sprintf('bcilab:/resources/testcases/log_%s.log',hlp_hostname),[],'Name/location of log file.'), ...
    arg({'report','ReportFile'},'bcilab:/resources/testcases/report_%s-%s.log',[],'Name/location of build report file.'));

% sanitize inputs
if ischar(o.scan_dirs)
    o.scan_dirs = {o.scan_dirs}; end
for d=1:length(o.scan_dirs)
    o.scan_dirs{d} = env_translatepath(o.scan_dirs{d}); end
o.log = env_translatepath(o.log);
o.report = env_translatepath(o.report);
if ischar(o.post_test)
    o.post_test = str2func(o.post_test); end

% names of marker files
mrk_untested = [o.proj_dir 'untested'];
mrk_testing = [o.proj_dir 'testing'];
mrk_tested = [o.proj_dir 'tested'];

% turn off a few warnings
if ispc
    warning off MATLAB:FILEATTRIB:SyntaxWarning; end
warning off MATLAB:DELETE:FileNotFound;

% run...
quicklog(o.log,'===== Now running as test slave =====');
while 1
    % get the date of the last test (if any)
    testdate = timestamp(mrk_tested);

    % get the most recent date of scanned directories
    codedate = max(cellfun(@timestamp,o.scan_dirs));

    % check if we are still up-to-date
    if codedate > testdate
        quicklog(o.log,'out of sync');
        rmfile(mrk_tested);
        makefile(mrk_untested);
        rmfile(mrk_testing);
    end

    % check if wait period expired ...
    if (codedate - testdate)*24*60*60 > o.wait_period
        % issue a re-test...
        quicklog(o.log,'re-testing...');
        makefile(mrk_testing);
        rmfile(mrk_untested);
        rmfile(mrk_tested);

        % issue a re-test
        try
            results = o.issue_test();
        catch e
            quicklog(o.log,evalc('env_handleerror(e)'));
            rethrow(e);
        end

        quicklog(o.log,'finished.');

        % check whether the code has been edited during the build
        newdate = max(cellfun(@timestamp,o.scan_dirs));
        if (newdate > codedate) && codedate
            quicklog(o.log,'code has been edited concurrently.');
            % if so, revert to untested
            makefile(mrk_untested);
            rmfile(mrk_tested);
            rmfile(mrk_testing);
        else
            quicklog(o.log,'testing completed.');
            if ~isempty(o.post_test)
                % run the post-testing step
                try
                    quicklog(o.log,'running post-testing step...');
                    o.post_test(results);
                    quicklog(o.log,'post-testing step completed successfully...');
                catch e
                    quicklog(o.log,evalc('env_handleerror(e)'));
                    rethrow(e);
                end
            end
            % test completed successfully
            makefile(mrk_tested);
            rmfile(mrk_untested);
            rmfile(mrk_testing);
            % write/update report, if applicable
            ds = datestr(clock); ds(ds == ':' | ds == ' ') = '-';
            if ischar(results)
                quicklog(sprintf(o.report,hlp_hostname,ds),results); end
        end
    else
        if codedate > testdate
            quicklog(o.log,'waiting for edits to finish...'); end
        pause(o.scan_interval);
    end
end

% calc the most recent time stamp of a given file system location
function ts = timestamp(loc)
if any(loc=='*') || ~exist(loc,'dir')
    % assume that it's a file mask / reference
    infos = dir(loc);
else
    % assume that it's a directory reference
    infos = cellfun(@dir,hlp_split(genpath(loc),pathsep),'UniformOutput',false);
    infos = vertcat(infos{:});
end
if ~isempty(infos)
    % take the maximum time stamp
    if ~isfield(infos,'datenum')
        [infos.datenum] = celldeal(cellfun(@datenum,{infos.date},'UniformOutput',false)); end
    ts = max([infos.datenum]);
else
    % no file: set to the beginning of time...
    ts = 0;
end

% create a status marker file...
function makefile(name)
try
    fid = fopen(name,'w+');
    fclose(fid);
    fileattrib(name,'+w','a');
catch
    if fid ~= -1
        try fclose(fid); catch,end
    end
end

% remove a status marker file...
function rmfile(name)
try
    delete(name);
catch
end

repo: REST


text:
function r = cache_file_is_there_and_newer_than(cache_file,main_file)

r = false;
main_info = dir(main_file);
main_date = datenum(main_info.date);
cache_info = dir(cache_file);
if(~isempty(cache_info)) %if there is a cache file (and it is newer than the main file)
    cache_date = datenum(cache_info.date);
    if  cache_date>main_date
        r = true;
    end
end

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function i4mat_max_test ( )

%*****************************************************************************80
%
%% I4MAT_MAX_TEST tests I4MAT_MAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 7;
  b = 0;
  c = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4MAT_MAX_TEST\n' );
  fprintf ( 1, '  I4MAT_MAX returns the maximum;\n' );

  seed = 123456789;

  [ a, seed ] = i4mat_uniform_ab ( m, n, b, c, seed );

  i4mat_print ( m, n, a, '  Random array:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum entry = %d', i4mat_max ( m, n, a ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% To plot Basin of Attraction for a given set of equations using Newton
% Raphson Method
% Governing Equations are:
% x^3-y = 0 ; y^3 - x = 0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Warning : On running this the workspace memory will be deleted. Save if
% any data present before running the code !!
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%--------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula, PhD                           |
%                   Structural Mechanics Laboratory                       |
%                   Indira Gandhi Center for Atomic Research              |
%                   India                                                 |
% E-mail : allwayzitzme@gmail.com                                         |
% web-link: https://sites.google.com/site/kolukulasivasrinivas/           |
%--------------------------------------------------------------------------
% Version 1 : 19 September 2013

clc ; clear all
warning('off') % To off the warning which shows "Matrix is close to singular
               % badly scaled" when algorithm passes through a point where the Jacobian
               % matrix is singular

% The roots of the given governing equations
r1 = [-1 ;-1] ;
r2 = [0 ;0] ;
r3 = [1 ;1] ;
% Initial conditions
x = linspace(-2,2,200) ;
y = linspace(-2,2,200) ;
% Initialize the required matrices
Xr1 = [] ; Xr2 = [] ; Xr3 = [] ; Xr4 = [] ;
tic
for i = 1:length(x)
     for j = 1:length(y)
          X0 = [x(i);y(j)] ;
          % Solve the system of Equations using Newton's Method
          X = NewtonRaphson(X0) ;
          % Locating the initial conditions according to error
          if norm(X-r1)<1e-8
             Xr1 = [X0 Xr1]  ;
          elseif norm(X-r2)<1e-8
             Xr2 = [X0 Xr2] ;
          elseif norm(X-r3)<1e-8
             Xr3 = [X0 Xr3] ;
          else               % if not close to any of the roots
             Xr4 = [X0 Xr4] ;
          end

     end
end
toc
warning('on') % Remove the warning off constraint

% Initialize figure
figure
set(gcf,'color','w')
hold on
plot(Xr1(1,:),Xr1(2,:),'.','color','r') ;
plot(Xr2(1,:),Xr2(2,:),'.','color','b') ;
plot(Xr3(1,:),Xr3(2,:),'.','color','g') ;
plot(Xr4(1,:),Xr4(2,:),'.','color','k') ;
title('Basin of attraction for f(x,y) = x^3-y = 0 and y^3-x=0')

author: Sable

repo: mcbench-benchmarks


text:
function F = sign(F, pref)
%SIGN   Sign function of a CHEBFUN.
%   G = SIGN(F) returns a piecewise constant CHEBFUN G such that G(x) = 1 in the
%   interval where F(x) > 0, G(x) = -1 in the interval where F(x) < 0 and G(x) =
%   0 in the interval where F(x) = 0. Breakpoints in G are introduced at zeros
%   of F.
%
%   For the nonzero values of complex F, SIGN(F) = F./ABS(F)
%
% See also ABS, HEAVISIDE, ROOTS.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Deal with the empty case:
if ( isempty(F) )
    return
end

if ( nargin < 2 )
    pref = chebfunpref();
end

for k = 1:numel(F)
    F(k) = signColumn(F(k), pref);
end

end

function g = signColumn(f, pref)

% Add breaks at the appropriate roots of f:
g = addBreaksAtRoots(f, pref);

% Call SIGN on each of the FUNs: (result will be smooth)
for k = 1:numel(g.funs)
    g.funs{k} = sign(g.funs{k}, pref);
end

% Take the sign of the pointValues in the first row:
g.pointValues = sign(g.pointValues);

% Remove unnecessary breakpoints:
g = merge(g, pref);

end


text:
function [ objects ] = rectangle2hypothesis( xyzBox, rule, room_points )
%RECTANGLE2HYPOTHESIS Popup rectangles to cuboids
%   Detailed explanation goes here
COUNTER_CLOCKWISE = true;
FITTINGRULES;
uni_room_points = zeros(8,3);
uni_room_points(POINTMAPPING{3},1) = room_points(:,1);
uni_room_points(POINTMAPPING{3},2) = room_points(:,2);
uni_room_points(POINTMAPPING{3},3) = room_points(:,3);

I = find(uni_room_points(1:4,1)<0 & uni_room_points(1:4,2)<0);
ID = rem([I I+1 I+2 I+3]-1,4) + 1;
uni_room_points = uni_room_points([ID ID+4],:);
wall_points = zeros(4,3,6);
wall_points(:,:,1) = uni_room_points([1 4 8 5],:);
wall_points(:,:,2) = uni_room_points([3 2 6 7],:);
wall_points(:,:,3) = uni_room_points([2 1 5 6],:);
wall_points(:,:,4) = uni_room_points([4 3 7 8],:);
wall_points(:,:,5) = uni_room_points([1 2 3 4],:);
wall_points(:,:,6) = uni_room_points([8 7 6 5],:);

LRDU_SURF = zeros(6,4);
LRDU_SURF(1,:) = [3 4 5 6];
LRDU_SURF(2,:) = [4 3 5 6];
LRDU_SURF(3,:) = [2 1 5 6];
LRDU_SURF(4,:) = [1 2 5 6];
LRDU_SURF(5,:) = [1 2 4 3];
LRDU_SURF(6,:) = [1 2 3 4];

normals = [1 0 0; -1 0 0; 0 1 0; 0 -1 0; 0 0 1; 0 0 -1];
normal_ids = [1 1 2 2 3 3];
main_angle = [pi/2 -pi/2 -pi 0 0];

objects = repmat(struct('out_points_w',[],'name',[],'type',[],'points',[],'x_w',[]),0,1);
object = struct('out_points_w',[],'name','other','type',[],'points',[],'x_w',zeros(7,1));
xyzBox = reshape(xyzBox', 3, 4)';

%% if left two points in lhs wall, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,1)), xyzBox([1 4],:), 0 );
if all(inside)
    LID = LRDU_SURF(rule,1);
    wall_point = wall_points(:,:,LID);
    point4 = LineFaceIntersection(wall_point(1,:), normals(LID,:), [0 0 0], xyzBox(1,:));
    point3 = LineFaceIntersection(wall_point(1,:), normals(LID,:), [0 0 0], xyzBox(4,:));
    point1 = LineFaceIntersection(point4, normals(rule,:), [0 0 0], xyzBox(2,:));
    point2 = LineFaceIntersection(point3, normals(rule,:), [0 0 0], xyzBox(3,:));
    point5 = point1;    point6 = point2;
    point5(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point6(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 9;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point6-point2);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

%% if right two points in rhs wall, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,2)), xyzBox([2 3],:), 0 );
if all(inside)
    RID = LRDU_SURF(rule,2);
    wall_point = wall_points(:,:,RID);
    point5 = LineFaceIntersection(wall_point(1,:), normals(RID,:), [0 0 0], xyzBox(2,:));
    point6 = LineFaceIntersection(wall_point(1,:), normals(RID,:), [0 0 0], xyzBox(3,:));
    point1 = LineFaceIntersection(point5, normals(rule,:), [0 0 0], xyzBox(1,:));
    point2 = LineFaceIntersection(point6, normals(rule,:), [0 0 0], xyzBox(4,:));
    point4 = point1;    point3 = point2;
    point4(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point3(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 9;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point6-point2);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

%% if bottom two points in floor, add a hypothesis
[inside,~,~] = insideCone( wall_points(:,:,LRDU_SURF(rule,3)), xyzBox([4 3],:), 0 );
if all(inside)
    DID = LRDU_SURF(rule,3);
    wall_point = wall_points(:,:,DID);
    point3 = LineFaceIntersection(wall_point(1,:), normals(DID,:), [0 0 0], xyzBox(4,:));
    point2 = LineFaceIntersection(wall_point(1,:), normals(DID,:), [0 0 0], xyzBox(3,:));
    point4 = LineFaceIntersection(point3, normals(rule,:), [0 0 0], xyzBox(1,:));
    point1 = LineFaceIntersection(point2, normals(rule,:), [0 0 0], xyzBox(2,:));
    point5 = point4;    point6 = point1;
    point5(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);
    point6(normal_ids(rule)) = wall_points(1,normal_ids(rule),rule);

    object.out_points_w = [point1; point2; point3; point4; point5; point6];
    object.type = 13;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point5-point4);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

%% if all four points inside the wall, add a rectangle
[inside,~,~] = insideCone( wall_points(:,:,rule), xyzBox, 0 );
if all(inside)
    wall_point = wall_points(:,:,rule);
    point4 = LineFaceIntersection(wall_point(1,:), normals(rule,:), [0 0 0], xyzBox(1,:));
    point1 = LineFaceIntersection(wall_point(1,:), normals(rule,:), [0 0 0], xyzBox(2,:));
    point2 = LineFaceIntersection(wall_point(1,:), normals(rule,:), [0 0 0], xyzBox(3,:));
    point3 = LineFaceIntersection(wall_point(1,:), normals(rule,:), [0 0 0], xyzBox(4,:));
    point5 = point4;
    point8 = point1;
    point7 = point2;
    point6 = point3;

    object.out_points_w = [point1; point2; point3; point4; point5; point6; point7; point8];
    object.type = 1;
    object.points = object.out_points_w ./ repmat(sqrt(sum(object.out_points_w.^2,2)),1,3);
    object.x_w(1:3) = point3;
    object.x_w(4) = norm(point2-point3);
    object.x_w(5) = norm(point7-point2);
    object.x_w(6) = norm(point4-point3);
    object.x_w(7) = main_angle(rule);
    objects(end+1,1) = object;
end

end

author: zouchuhang

repo: LayoutNet


text:
function [Dtest,modelF,allF]=spm_eeg_invertiter(Dtest,Npatchiter,funcname,patchind)

%  Function to perform several MSP type inversions with different
%  pseudo-randomly selected priors- in this case single cortical patches
%
% Npatchiter: number of iterations
% funcname is name of MSP alogorithm: current (spm_eeg_invert) or classic (spm_eeg_invert_classic)
% patchind is an optional list of indices of vertices which will be patch
% centres. patchind will have size Npatchiter*Np (where Np is number of patches set in
% inverse.Np )
% if  Dtest{1}.inv{val}.inverse.mergeflag==1 then merges Npatchiter posterior current
% distributions, else replaces posterior with best of the iterations.
% __________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging
%
% Gareth Barnes
% $Id: spm_eeg_invertiter.m 6367 2015-03-09 15:02:25Z gareth $

if nargin<2,
    Npatchiter=[];
end;

if nargin<3,
    error('No function call specified');
end;

if nargin<4,
    patchind=[];
end;
if isempty(Npatchiter),
    Npatchiter=16;
end;

if numel(Dtest)>1,
    error('only works with single datasets at the moment');
end;

val=Dtest{1}.val;

Nvert=size(Dtest{1}.inv{val}.mesh.tess_mni.vert,1);
Np=Dtest{1}.inv{val}.inverse.Np;

allF=zeros(Npatchiter,1);

fprintf('Checking leadfields');

[L Dtest{1}] = spm_eeg_lgainmat(Dtest{1});  % Generate/load lead field- this stops it being done at each iteration

if isempty(patchind),
    disp('Reseting random number seed ! and then generating random patch centres ');
    rand('state',0);

    %% make random patch centers, except on the first iteration when we keep to a fixed stepsize
    for f=1:Npatchiter,
        tmp=randperm(Nvert);
        allIp(f).Ip=tmp(1:Np);
    end;
    allIp(1).Ip=[]; %% fix the step size of the first patch set (uses default MSP indices)
else
    for f=1:Npatchiter,
        allIp(f).Ip=patchind(f,:);
    end;
end;

modelF=[];
bestF=-Inf;
for patchiter=1:Npatchiter,
    %Din=Dtest{1};
    %Dout=Din;
    Ip=allIp(patchiter).Ip;
    disp(sprintf('patchiter %d/%d',patchiter,Npatchiter));
    switch funcname,
        case 'Classic',
            Dtest{1}.inv{val}.inverse.Ip=Ip;
            Dtest{1}    = spm_eeg_invert_classic(Dtest{1});
        case 'Current',
            warning('Patch centres are currently fixed for this algorithm (iteration will have no effect!)');
            Dtest{1}   = spm_eeg_invert(Dtest{1}); %
            %Dout.inv{val}.inverse.Ip=Ip;
        otherwise
            error('unknown function');
    end;
    if (Dtest{1}.inv{val}.inverse.mergeflag==1), %% will need all posteriors to merge them (could consider getting rid of worst ones here to save memory
        modelF(patchiter).inverse=Dtest{1}.inv{val}.inverse;
    else %% just keep track of best inversion
        if Dtest{1}.inv{val}.inverse.F>bestF,
            modelF(1).inverse=Dtest{1}.inv{val}.inverse;
            bestF=Dtest{1}.inv{val}.inverse.F;
        end;
    end;
    allF(patchiter)=Dtest{1}.inv{val}.inverse.F;
end; % for patchiter

%% will use these to merge posteriors or take best one

[bestF,bestind]=max(allF);
disp('model evidences relative to maximum:')

sort(allF-bestF)

%% for 1 iteration or for best solution
Dtest{1}.inv{val}.inverse=modelF(1).inverse; %% return best model (only 1 model saved in this option)
Dtest{1}.inv{val}.inverse.allF=allF;

if Npatchiter>1, %% keep iterations if more than 1

    if (Dtest{1}.inv{val}.inverse.mergeflag==1)
        Dtest{1}.inv{val}.inverse=modelF(bestind).inverse; % start with best model so far
        Qpriors=sparse(zeros(Npatchiter,size(modelF(1).inverse.qC,1)));
        for patchiter=1:Npatchiter,
            Qpriors(patchiter,:)=modelF(patchiter).inverse.qC;
        end;
        disp('Merging posterior distributions..');
        ugainfiles=Dtest{1}.inv{val}.gainmat;
        surfind=ones(Npatchiter,1);
        %% now mix the posteriors
        [Dtest{1}] = spm_eeg_invert_classic_mix(Dtest{1},val,Qpriors,surfind,ugainfiles);

        Dtest{1}.inv{val}.comment{1}=sprintf('Merged posterior of %d solutions',Npatchiter);
        mixF=Dtest{1}.inv{val}.inverse.F;
        if mixF-bestF<0
            warning('Merged solution is worse than the best');
        end;
        disp(sprintf('Improvement (when +ve) over best iteration %3.2f',mixF-bestF));

    else % NO merge - just take the best

        disp('Using best patch set to current estimate');
        Dtest{1}.inv{val}.comment{1}=sprintf('Best F of %d solutions',Npatchiter);
    end; % if BMA

end;

spm_eeg_invert_display(Dtest{1});

author: spm


text:
function x0 = randomObjFBASol(model, initArgs)
% Solves an FBA problem with a random objective function
%
% USAGE:
%
%    x0 = randomObjSol(model, initArgs)
%
% INPUTS:
%    model:        COBRA model structure
%    initArgs:     Cell array containing the following data:
%
%                    1.  osenseStr - Maximize ('max') / minimize ('min')
%                    2.  minObjFrac - Minimum initial objective fraction
%                    3.  minObjValue - Minimum initial objective value (opt)
%                        (Default = `minObjFrac*sol.f`)
%
% OUTPUT:
%    x0:           solution of FBA problem
%
% .. Author: - Markus Herrgard

osenseStr = initArgs{1};
minObjFrac = initArgs{2};

if (length(initArgs) < 3)
    solOpt = optimizeCbModel(model,osenseStr);
    model.lb(model.c==1) = minObjFrac*solOpt.f;
else
    model.lb(model.c==1) = initArgs{3};
end

nRxns = length(model.rxns);

model.c = rand(nRxns,1)-0.5;
sol = optimizeCbModel(model,osenseStr);
x0 = sol.x;

author: opencobra

repo: cobratoolbox


text:
function test_suite = test_edgeLength
% One-line description here, please.
%   output = testEdgeLength(input)
%
%   Example
%   testEdgeLength
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSimple(testCase) %#ok<*DEFNU>

p1 = [2 3];
p2 = [5 7];
edge = createEdge(p1, p2);
len = edgeLength(edge);
testCase.assertEqual(5, len, 'AbsTol', .01);

function testArray(testCase)

p1 = [2 3;2 3;3 2;5 7];
p2 = [5 7;5 7;7 5;2 3];
edge = createEdge(p1, p2);
len = edgeLength(edge);

testCase.assertEqual(4, size(len, 1));
testCase.assertEqual([5;5;5;5], len, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [res] = fiff_pick_channels_evoked(orig,include,exclude)
%
% [res] = fiff_pick_channels_evoked(orig,include,exclude)
%
% Pick desired channels from evoked-response data
%
% orig      - The original data
% include   - Channels to include (if empty, include all available)
% exclude   - Channels to exclude (if empty, do not exclude any)
%
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.9  2006/10/04 20:12:37  msh
%   Added fiff_read_evoked_all
%   Modified fiff_pick_channels_evoked to handle multiple data sets
%   Added bad channel writing to fiff_write_evoked
%
%   Revision 1.8  2006/06/29 22:12:28  msh
%   Fixed errors in channel picking
%
%   Revision 1.7  2006/04/26 00:14:24  msh
%   Return empty structure from fiff_pick_channels evoked and fiff_pick_types_evoked if there is no match.
%   Accuracy checking was incorrect in mne_add_coil_defs
%
%   Revision 1.6  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.5  2006/04/22 10:59:30  msh
%   Added fiff_pick_info
%
%   Revision 1.4  2006/04/20 23:33:20  msh
%   Separated fiff_pick_channels and fiff_pick_types from the _evoked versions
%
%   Revision 1.3  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.2  2006/04/14 15:49:49  msh
%   Improved the channel selection code and added ch_names to measurement info.
%
%   Revision 1.1  2006/04/14 00:45:42  msh
%   Added channel picking and fiff_invert_transform
%
%

me='MNE:fiff_pick_channels_evoked';

if nargin == 1
    res = orig;
    return;
elseif nargin == 2
    exclude = [];
elseif nargin ~= 3
    error(me,'Incorrect number of arguments');
end

if isempty(include) && isempty(exclude)
    res = orig;
    return;
end

sel = fiff_pick_channels(orig.info.ch_names,include,exclude);
if isempty(sel)
    res = [];
    fprintf(1,'Warning : No channels match the selection.\n');
    return;
end

res = orig;
%
%   Modify the measurement info
%
res.info = fiff_pick_info(res.info,sel);
%
%   Create the reduced data set
%
for k = 1:length(res.evoked)
    res.evoked(k).epochs = res.evoked(k).epochs(sel,:);
end

return;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function complainif_toomanyargs(fnargin,limit,callfun)

if nargin<3
    callfun = mfilename;
end

if fnargin>limit
   error('%s: Too many input arguments.',upper(callfun));
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function opts = mergeOptions(opts1, opts2)
% Merges two options structures with one having precedence over the other.
%
% function opts = mergeOptions(opts1, opts2)
%
% input: opts1 and opts2 are two structures.
% output: opts is a structure containing all fields of opts1 and opts2.
% Whenever a field is present in both opts1 and opts2, it is the value in
% opts2 that is kept.
%
% The typical usage is to have opts1 contain default options and opts2
% contain user-specified options that overwrite the defaults.
%
% See also: getGlobalDefaults

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    if isempty(opts1)
        opts1 = struct();
    end
    if isempty(opts2)
        opts2 = struct();
    end

    opts = opts1;
    fields = fieldnames(opts2);
%     for i = 1 : length(fields)
%         opts.(fields{i}) = opts2.(fields{i});
%     end

    for i = 1 : length(fields)
        if strcmp(fields{i}, 'max_epoch') || strcmp(fields{i}, 'verbose')
            if ischar(opts2.(fields{i}))
                opts.(fields{i}) = str2double(opts2.(fields{i}));
            else
                opts.(fields{i}) = opts2.(fields{i});
            end
        else
            opts.(fields{i}) = opts2.(fields{i});
        end

    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function plot_updatePlot(hFig, idxPC)

% PC in UserData of figure -> PC -----------------------------------------------

PC = hFig.UserData.PC{idxPC};

% Delete points and normals of previous plot -----------------------------------

% Delete points
if ~iscell(PC.U.hPlot)
    delete(PC.U.hPlot);
else
    for j = 1:numel(PC.U.hPlot)
        delete(PC.U.hPlot{j});
    end
end

% Delete normals
delete(PC.U.hPlotNormals);

% Consider parameter 'MaxPoints' -----------------------------------------------

if sum(PC.act) > PC.U.p.MaxPoints
    PC.U.percentOfVisiblePoints = PC.U.p.MaxPoints/sum(PC.act)*100;
    idxAct = find(PC.act); % indices of active points
    idxRandom = randperm(numel(idxAct), PC.U.p.MaxPoints); % indices of randomly selected points
    act4plot = false(PC.noPoints,1);
    act4plot(idxAct(idxRandom)) = true; % reactivate randomly selected points
else
    PC.U.percentOfVisiblePoints = 100;
    act4plot = PC.act;
end

% Plot -------------------------------------------------------------------------

if PC.U.Visible && sum(act4plot) > 0

    % Axes handle
    hAxes = findobj('Parent', hFig, 'Type', 'Axes');

    % Colorize points according to an attribute
    if numel(PC.U.p.Color) > 1 && strcmp(PC.U.p.Color(1:2), 'A.') && ~strcmpi(PC.U.p.Color, 'A.RGB') % if p.Color begins with 'A.'

        % Values of attribute to plot
        attribute2plot = PC.U.p.Color(3:end); % attribute name
        if strcmpi(attribute2plot, 'x') % attribute = x coordinate
            A = PC.X(act4plot,1);
        elseif strcmpi(attribute2plot, 'y') % attribute = y coordinate
            A = PC.X(act4plot,2);
        elseif strcmpi(attribute2plot, 'z') % attribute = z coordinate
            A = PC.X(act4plot,3);
        else
            A = PC.A.(attribute2plot)(act4plot);
        end

        % Plot!
        PC.U.hPlot = scatter3(hAxes, PC.X(act4plot,1), PC.X(act4plot,2), PC.X(act4plot,3), PC.U.p.MarkerSize^2, A, '.');

        % Colormap
        colormap(PC.U.p.ColormapName);

        % CAxisLim
        if strcmpi(PC.U.p.ColormapName, 'classpal') % for visualization of classes with maximum distiguishable colors

            caxis([-0.5 50.5]); % for no. of classes = 50, see classpal.m

        else % usual case

            if isempty(PC.U.p.CAxisLim)
                CAxisLim = [min(A) max(A)];
            else
                CAxisLim = PC.U.p.CAxisLim;
            end
            caxis(CAxisLim);

        end

        % % Colorbar
        % hColorbar = findobj('Type', 'Colorbar', 'Parent', gcf);
        % if ~isempty(hColorbar) % applies if min and max of CAxisLim is equal (e.g. if attribute contains only a single value)
        %     PC.U.p.CAxisLim = [hColorbar.Limits(1) hColorbar.Limits(2)];
        %     hColorbar.Label.String = PC.label;
        %     hColorbar.Label.Interpreter = 'none';
        % end

    % RGB plot
    elseif strcmpi(PC.U.p.Color, 'A.RGB')

        colorbar off;

        % Check how colors are saved
        if isfield(PC.A, 'r') && isfield(PC.A, 'g') && isfield(PC.A, 'b')
            r = PC.A.r;
            g = PC.A.g;
            b = PC.A.b;
        elseif isfield(PC.A, 'red') && isfield(PC.A, 'green') && isfield(PC.A, 'blue')
            r = PC.A.red;
            g = PC.A.green;
            b = PC.A.blue;
        elseif isfield(PC.A, 'Red') && isfield(PC.A, 'Green') && isfield(PC.A, 'Blue')
            r = PC.A.Red;
            g = PC.A.Green;
            b = PC.A.Blue;
        elseif isfield(PC.A, 'diffuse_red') && isfield(PC.A, 'diffuse_green') && isfield(PC.A, 'diffuse_blue')
            r = PC.A.diffuse_red;
            g = PC.A.diffuse_green;
            b = PC.A.diffuse_blue;
        end

        % Check if colors were found
        if ~exist('r', 'var') | ~exist('g', 'var') | ~exist('b', 'var')
            errordlg('No RGB attributes found!');
            return;
        end

        % Conversion of colors
        if max(r) <= 1
            % do nothing
        elseif max(r) > 1 && max(r) <= 2^8-1
            r = r/(2^8-1); g = g/(2^8-1); b = b/(2^8-1);
        elseif max(r) > 255 && max(r) <= 2^16-1
            r = r/(2^16-1); g = g/(2^16-1); b = b/(2^16-1);
        end

        % Plot!
        PC.U.hPlot = scatter3(hAxes, PC.X(act4plot,1), PC.X(act4plot,2), PC.X(act4plot,3), PC.U.p.MarkerSize^2, [r(act4plot) g(act4plot) b(act4plot)], '.');

    % Unicolor plot
    else

        colorbar off;

        % Plot!
        if strcmpi(PC.U.p.Color, 'by PC')
            nColor = size(hFig.UserData.colors,1);
            idxColor = mod(idxPC,nColor); if idxColor == 0, idxColor = nColor; end % get color index
            PC.U.p.Color = hFig.UserData.colors{idxColor,1};
        end

        PC.U.hPlot = plot3(hAxes, PC.X(act4plot,1), PC.X(act4plot,2), PC.X(act4plot,3), '.', 'MarkerSize', PC.U.p.MarkerSize, 'Color', PC.U.p.Color);

    end

    % Plot normals
    if hFig.UserData.PC{idxPC}.U.plotNormals
        PC.U.hPlotNormals = PC.plotNormals('Scale', hFig.UserData.options.NormalsScale, 'Color', hFig.UserData.options.NormalsColor, 'Axes', hAxes); % plot normals!
    end

end

% PC -> UserData of figure -----------------------------------------------------

hFig.UserData.PC{idxPC} = PC;

% Update menus -----------------------------------------------------------------

plot_updateMenus(hFig)

% Update histogram -------------------------------------------------------------

plot_updateHisto(hFig)

end

author: pglira


text:
function K=centKer(K)

%function K=centKer(K)
%
% Centeres the matrix K
%
%INPUTS
% K = the kernel matrix to be centered
%
%OUTPUTS
% Kc = the centered kernel matrix
%
%
%For more info, see www.kernel-methods.net

% original kernel matrix stored in variable K
% output uses the same variable K
% K is of dimension ell x ell
% D is a row vector storing the column averages of K
% E is the average of all the entries of K
ell = size(K,1);
D = sum(K) / ell;
E = sum(D) / ell;
J = ones(ell,1) * D;
K = K - J - J' + E;


text:
%% Code to grab data, parse, adjust timing
% Adam Werries 2016, see Apache 2.0 license.
addpath('GrovesCode');
addpath('Utilities');
addpath('Tuning');

if ~exist('last_applanix_dir', 'var') || ~ischar(last_applanix_dir)
   last_applanix_dir = pwd;
end
text_files = {'*.txt;*.csv;*.log','Data files (*.txt,*.csv,*.log)'; '*.*', 'All Files (*.*)'};
%% Import Applanix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ground truth import
disp('Please select Applanix log file')
[applanix_file, applanix_path] = uigetfile(text_files, 'Select Applanix Log File', last_applanix_dir);
last_applanix_dir = applanix_path;

gt = csvread([applanix_path applanix_file]);
gt_t = gt(:,1);
gt_vel = gt(:,8:10);
gt_speed = gt(:,21);
gt_rot = gt(:,11:13);
gt_accel = gt(:,14:16);
gt_xyz = gt(:,27:29);

%% Import gps data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Format:  1           2         3         4         5           6     7       8
%         start_time, end_time, gps_time, fix_mode, num_sat,    lat,  long,   elev,
%         9       10       11       12       13       14       15    16    17    18    19
%         ecef_x, ecef_y,  ecef_z,  ecef_vx, ecef_vy, ecef_vz, gdop, pdop, hdop, vdop, tdop
if ~exist('last_lowcost_dir', 'var') || ~ischar(last_lowcost_dir)
   last_lowcost_dir = last_applanix_dir;
end
disp('Please select Skytraq log file')
[skytraq_file, skytraq_path] = uigetfile(text_files, 'Select Skytraq Log File', last_lowcost_dir);
last_lowcost_dir = skytraq_path;

gps = csvread([skytraq_path skytraq_file],1,0);
gps_time = gps(:,1);
% grab lat-long
lla = gps(:,6:8);
% Convert latlon to UTM
[gps_x,gps_y,utmzone] = deg2utm(lla(:,1),lla(:,2));
gps_h = -lla(:,3);

%% Import IMU data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Format:  1                2              3        4        5        6       7       8
%         meas_start_time, meas_end_time, accel_x, accel_y, accel_z, gyro_x, gyro_y, gyro_z
disp('Please select IMU log file')
[imu_file, imu_path] = uigetfile(text_files, 'Select IMU Log File', last_lowcost_dir);
imu = csvread([imu_path imu_file]);
% only using end_time here
imu = imu(:,2:end);
% convert g's to m/s^2 using local gravity estimate
imu(:,2:4) = imu(:,2:4)*9.80097;
% rotate sensor readings 180 deg about the x-axis
rotation = [1   0         0;
            0   cosd(180) -sind(180)
            0   sind(180) cosd(180)];
imu(:,2:4) = (rotation*imu(:,2:4)')';
imu(:,5:7) = (rotation*imu(:,5:7)')';
%% Correct time racnges %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
timelag = -0.0963798;
[start_time,Is] = max([min(gps_time),min(gt_t),min(imu(:,1))]);
[end_time,Ie] = min([max(gps_time),max(gt_t),max(imu(:,1))]);
end_time = end_time - start_time;
gt_t = gt_t - start_time;
gps_time = gps_time - start_time + timelag;
imu(:,1) = imu(:,1) - start_time;
% cut off extra gps data at end, fold time back in
gps_mask = gps_time < end_time & gps_time > 0;
gps_time = gps_time(gps_mask);
gps = gps(gps_mask,:);
lla = lla(gps_mask,:);
gps(:,1) = gps_time;
% cut off extra IMU data at end
imu = imu(imu(:,1) < end_time, :);
imu = imu(imu(:,1) > 0, :);
imu_time = imu(:,1);

%% Generate filter time, configurable epochs
epoch = 0.04;
filter_time = 0:epoch:imu_time(end);

%% Generate ground truth

min_x = min([min(gps_x),min(gt_xyz(:,2))]);
min_y = min([min(gps_y),min(gt_xyz(:,1))]);
min_z = min([min(gps_h),min(gt_xyz(:,3))]);
gps_x = gps_x(gps_mask) - min_x;
gps_y = gps_y(gps_mask) - min_y;
gps_h = gps_h(gps_mask);
ground_truth = zeros(length(gps_time), 6);
ground_truth_full = zeros(length(filter_time), 9);

ground_truth(:,1) =  (interp1(gt_t,gt_xyz(:,2),gps_time)-min_x)';
ground_truth(:,2) =  (interp1(gt_t,gt_xyz(:,1),gps_time)-min_y)';
ground_truth(:,3) =  interp1(gt_t,gt_xyz(:,3),gps_time)';
ground_truth(:,4) =  interp1(gt_t,gt_rot(:,1),gps_time)';
ground_truth(:,5) =  interp1(gt_t,gt_rot(:,2),gps_time)';
ground_truth(:,6) =  interp1(gt_t,gt_rot(:,3),gps_time)';

ground_truth_full(:,1) = (interp1(gt_t,gt_xyz(:,2),filter_time)-min_x)';
ground_truth_full(:,2) = (interp1(gt_t,gt_xyz(:,1),filter_time)-min_y)';
ground_truth_full(:,3) = interp1(gt_t,gt_xyz(:,3),filter_time)';
ground_truth_full(:,4) = interp1(gt_t,gt_rot(:,1),filter_time)';
ground_truth_full(:,5) = interp1(gt_t,gt_rot(:,2),filter_time)';
ground_truth_full(:,6) = interp1(gt_t,gt_rot(:,3),filter_time)';
ground_truth_full(:,7) = interp1(gt_t,gt_vel(:,1),filter_time)';
ground_truth_full(:,8) = interp1(gt_t,gt_vel(:,2),filter_time)';
ground_truth_full(:,9) = interp1(gt_t,gt_vel(:,3),filter_time)';

author: awerries

repo: kalman-localization


text:
function y = nanstd(x, flag, dim)
% nanstd - Standard deviation ignoring NaNs
%
% Synopsis:
%   y = nanstd(x)
%   y = nanstd(x,flag)
%   y = nanstd(x,flag,dim)
%
% Arguments:
%   x: Matrix or vector
%   flag: Logical. Default: false (normalize by N-1). If true: normalize by N
%       (produce the second moment about the mean)
%   dim: Dimension along which to operate
%
% Returns:
%   y: Standard deviation along the chosen dimension
%
% Description:
%   nanstd(X) returns the standard deviation treating NaNs
%   as missing values.
%   For vectors, nanstd(X) is the standard deviation of the
%   non-NaN elements in X.  For matrices, nanstd(X) is a row
%   vector containing the standard deviation of each column,
%   ignoring NaNs.
%
% Examples:
%   nanstd([1 2 NaN]) returns 0.7071 (the same result as std([1 2])
%   nanstd([1 2 NaN; 3 4 5], 1, 1) returns [1 1 0]
%
% See also: nanmean,nansum
%

% Author(s): Anton Schwaighofer, Feb 2005

if nargin<3,
  % Operate along the first non-singleton dimension
  dim = min(find(size(x)~=1));
  if isempty(dim),
    dim = 1;
  end
end
if nargin<2,
  flag = 0;
end
if isempty(x),
  y = 0/0;
  return;
end
tile = ones(1, max(ndims(x), dim));
tile(dim) = size(x,dim);
% Effective number of non-NaNs
effSize = size(x,dim)-sum(isnan(x),dim);
% Remove mean, also ignoring NaNs
xc = x - repmat(nanmean(x, dim), tile);
if flag,
  y = sqrt(nansum(conj(xc).*xc, dim)./effSize);
else
  y = sqrt(nansum(conj(xc).*xc, dim)./(effSize-1));
  % Simulate the behaviour of std: if only one sample, output 0
  y(effSize==1) = 0;
end

author: bbci


text:
function p = exp_part(x,idx)
% Get a specific part of an expression; can be more efficient than exp_parts.
% Parts = exp_part(Expression, Index)
%
% In:
%   Expression  : any expression
%   Index       : index of the part;
%                  * 0 is the expression's head
%                  * 1..n is the expression's parts, n being exp_length(Expression)
%                  * negative indices count parts from the end
%                  * index vectors denote paths into the expression tree
%
% Out:
%   Part        : the selected part of the expression, as in exp_parts
%
% See also:
%   exp_head, exp_parts, exp_length
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-26

while isfield(x,'tracking') && isfield(x.tracking,'expression')
    x = x.tracking.expression; end

try
    if length(idx)==1
        if idx > 0
            % positive indices behave as exp_parts(x){idx} (if that could be written in MATLAB...)
            if isfield(x,{'head','parts'})
                p = x.parts{idx};
            elseif has_function_symbol(x)
                p = get_function_symbol(x);
            elseif idx == 1
                p = x;
            else
                error; %#ok<LTARG>
            end
        elseif idx == 0
            % part 0 is the head
            p = exp_head(x);
        elseif idx < 0
            % negative indices count from the end (-1 is the last part)
            p = exp_part(x,exp_length(x)+idx+1);
        end
    else
        % index vectors specify paths into the expression
        p = exp_part(exp_part(x,idx(1)),idx(2:end));
    end
catch
    error('out-of-range part accessed.');
end


text:
function stats = FieldPSP(stims,lfp,varargin)

%FieldPSP - Measure field EPSPs (waveforms, slope, amplitude) across time.
%
%  Slope is determined as the maximum slope during the time window around stimulation.
%  Amplitude is determined as the maximum amplitude during the same time window.
%
%  USAGE
%
%    stats = FieldPSP(stims,lfp,<options>)
%
%    stims          stimulation timestamps
%    lfp            local field potential samples
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'direction'   direction of the field EPSP, either 'up' (recording
%                   from str. pyr., default) or 'down' (recording from str.
%                   radiatum)
%     'durations'   durations before and after synchronizing events for each
%                   trial (in s) (default = [-0.01 0.02])
%     'slope'       where max slope should be sought (default = [0.001 0.005])
%     'amplitude'   where max amplitude should be sought
%                   (default = [0.005 0.010])
%     'mode'        compute slope and amplitude running averages using a
%                   fixed time window ('time'), or a window with a fixed
%                   number of stimulations ('count', default)
%     'window'      window length (in s or in counts) for running averages
%                   (default = 100)
%     'show'        either 'on' or 'off' (default)
%    =========================================================================
%
%  OUTPUT
%
%    stats.amplitude.t         times of max amplitude (relative to stimulation)
%    stats.amplitude.v         values of max amplitudes
%    stats.slope.t             times of max slopes (relative to stimulation)
%    stats.slope.v             values of max slopes
%    stats.psp.t               times of fPSP samples (relative to stimulation)
%    stats.psp.v               fPSP waveforms, one per stimulation
%

% Copyright (C) 2010-2012 by Michaël Zugaro
%
% This program 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.

% Default values
parent = [];
durations = [-0.01 0.02];
mode = 'count';
window = 100;
direction = 'up';
stats = [];
show = 'off';
slopeWindow = [0.001 0.005];
amplitudeWindow = [0.005 0.010];

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help FieldPSP'' for details).');
end

% Check parameter sizes
if ~isdvector(stims),
	error('Parameter ''stims'' is not a vector (type ''help FieldPSP'' for details).');
end
if ~isdmatrix(lfp),
	error('Parameter ''lfp'' is not a matrix (type ''help FieldPSP'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help FieldPSP'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'direction',
			direction = lower(varargin{i+1});
			if ~isstring_FMAT(direction,'up','down'),
				error('Incorrect value for property ''direction'' (type ''help FieldPSP'' for details).');
			end
		case 'durations',
			durations = varargin{i+1};
			if ~isdvector(durations,'#2','<'),
				error('Incorrect value for property ''durations'' (type ''help FieldPSP'' for details).');
			end
		case 'slope',
			slopeWindow = varargin{i+1};
			if ~isdvector(slopeWindow,'#2','<'),
				error('Incorrect value for property ''slope'' (type ''help FieldPSP'' for details).');
			end
		case 'amplitude',
			amplitudeWindow = varargin{i+1};
			if ~isdvector(amplitudeWindow,'#2','<'),
				error('Incorrect value for property ''amplitude'' (type ''help FieldPSP'' for details).');
			end
		case 'mode',
			mode = lower(varargin{i+1});
			if ~isstring_FMAT(mode,'time','count'),
				error('Incorrect value for property ''mode'' (type ''help FieldPSP'' for details).');
			end
		case 'parent',
			parent = varargin{i+1};
			if ~ishandle(parent),
				error('Incorrect value for property ''parent'' (type ''help FieldPSP'' for details).');
			end
		case 'window',
			window = varargin{i+1};
			if ~isdscalar(window,'>0'),
				error('Incorrect value for property ''window'' (type ''help FieldPSP'' for details).');
			end
		case 'show',
			show = lower(varargin{i+1});
			if ~isstring_FMAT(show,'on','off'),
				error('Incorrect value for property ''show'' (type ''help FieldPSP'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help FieldPSP'' for details).']);
	end
end

up = strcmp(direction,'up');
show = strcmp(show,'on');

% Some analyses depend on the method (fixed time vs fixed count), but others do not
% To simplify the code, we will use two variables:
%  - stims0 uses time, whatever the method
%  - stims  uses either time or count, depending on the method
stims0 = stims;
if strcmp(mode,'count'),
	stims = (1:length(stims))';
end

% Constants
nStims = length(stims);
height = max(lfp(:,2))-min(lfp(:,2));
nSlices = 12;
start = stims(1);
stop = stims(end);
window = (stop-start)/nSlices;
slices = (0:nSlices-1)*window+start;
nStims = length(stims);
frequency = 1/median(diff(lfp(:,1)));
smooth = frequency/1250;
nStimsPlottedPerSlice = 5;
nStimsPerSlice = nStims/nSlices;
nSkip = floor(nStims/(nSlices*(nStimsPlottedPerSlice-1))); % plot one stim every n
colors = Bright(nStimsPlottedPerSlice);

% Compute slope and amplitude for each stimulation

% Resynchronize all fPSPs (use only 15ms after stims)
[sync,index] = Sync(lfp,stims0,'durations',durations);

% Preallocate large matrix to store all fPSPs
n = sum(index==1);
stats.psp.v = nan(nStims,n);
stats.psp.t = linspace(durations(1),durations(2),n);
% Compute slope and amplitude for each stimulation
if show, figure; end
for i = 1:nStims,
	% Select the appropriate data, store in fPSP matrix, smooth and differentiate
	s = sync(index==i,:);
	stats.psp.v(i,:) = interp1(s(:,1),s(:,2),stats.psp.t);
	s(:,2) = Smooth(s(:,2),smooth);
	% Compute amplitude as the first local maximum (or minimum) in the appropriate time interval
	in = InIntervals(s(:,1),amplitudeWindow);
	if up,
		extrema = find(IsExtremum(s)&in);
	else
		extrema = find(IsExtremum(s,'mode','minima')&in);
	end
	if isempty(extrema),
		amplitude(i) = nan;
		tAmplitude(i) = nan;
		slope(i) = nan;
		tSlope(i) = nan;
		continue;
	end
	extremum = extrema(1);
	tAmplitude(i) = s(extremum,1);
	amplitude(i) = s(extremum,2);
	% Compute slope as the maximum slope in the appropriate time interval
	% (set all values outside this time interval to nan to make sure the max will be in the interval)
	ds = Diff(s);
	ds(~InIntervals(ds(:,1),slopeWindow),2) = nan;
	if up,
		[slope(i),j] = nanmax(ds(:,2));
	else
		[slope(i),j] = nanmin(ds(:,2));
	end
	tSlope(i) = ds(j,1);
	% Plot slope and amplitude
	slice = floor(i/nStimsPerSlice);
	stimInSlice = floor(i-1-slice*nStimsPerSlice);
	if show && rem(stimInSlice,nSkip) == 0 && stimInSlice >= 0,
		k = ceil(nSlices*i/nStims);
		l = floor(stimInSlice/nSkip)+1;
		if l <= nStimsPlottedPerSlice,
			SquareSubplot(nSlices,k);
			xlabel(num2str(slices(k)));
			hold on;
			PlotXY(s,'color',colors(l,:),'linewidth',3);
			PlotSlope(tAmplitude(i),amplitude(i),0,0.005,'r');
			PlotSlope(tSlope(i),s(j,2),slope(i),0.005,'k');
		end
	end
end

% Make sure all subplots use the same scale
if show, AdjustAxes('y','uniform'); end

stats.amplitude.t = tAmplitude;
stats.amplitude.v = amplitude;
stats.slope.t = tSlope;
stats.slope.v = slope;

% Additional figures: mean fPSPs, slopes and amplitudes (scatterplots + running averages), CVs, ISIs and short time CCGs

if show,
	% Mean fPSPs
	figure;
	for i = 1:nSlices,
		if strcmp(mode,'count'),
			start = floor((i-1)*nStimsPerSlice)+1;
			stop = floor(i*nStimsPerSlice);
			ok = start:stop;
		else
			ok = find(InIntervals(stims,[0 window]+start+(i-1)*window));
		end
		m = mean(stats.psp.v(ok,:));
		e = sem(stats.psp.v(ok,:));
		SquareSubplot(nSlices,i);
		xlabel(num2str(slices(i)));
		hold on;
		PlotMean(stats.psp.t,m,m-e/2,m+e/2,':');
	end
	% Make sure all subplots use the same scale
	AdjustAxes('y','uniform');

	figure;
	% 'Rectify' slope and amplitude
	if ~up,
		amplitude = -amplitude;
		slope = -slope;
	end
	if strcmp(mode,'time'),
		xLabel = 'Time (s)';
	else
		xLabel = 'Count';
	end
	% Slope (scatterplot)
	subplot(2,3,1);
	hold on;
	plot(stims,slope,'.','markersize',1);
	% Slope (running average)
	[ts,ms,es] = RunningAverage(stims,slope,'window',window,'overlap',0.5*window);
	PlotMean(ts,ms,es(:,1),es(:,2),':','k');
	PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
	xlim([stims(1,1) stims(end,1)]);
	ylabel('Slope (a.u.)');
	set(gca,'xtick',[]);
	% Slope time (scatterplot)
	subplot(2,3,2);
	hold on;
	plot(stims,tSlope*1000,'.','markersize',1);
	% Slope time (running average)
	[ts,ms,es] = RunningAverage(stims,tSlope*1000,'window',window,'overlap',0.5*window);
	PlotMean(ts,ms,es(:,1),es(:,2),':','k');
	PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
	xlim([stims(1,1) stims(end,1)]);
	ylabel('Slope time (ms)');
	set(gca,'xtick',[]);

	% Amplitude (scatter)
	subplot(2,3,4);
	hold on;
	plot(stims,amplitude,'.','markersize',1);
	% Amplitude (running average)
	[ta,ma,ea] = RunningAverage(stims,amplitude,'window',window,'overlap',0.5*window);
	PlotMean(ta,ma,ea(:,1),ea(:,2),':','k');
	PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
	xlim([stims(1,1) stims(end,1)]);
	xlabel(xLabel);
	ylabel('Amplitude (a.u.)');
	% Amplitude time (scatter)
	subplot(2,3,5);
	hold on;
	plot(stims,tAmplitude*1000,'.','markersize',1);
	% Amplitude time (running average)
	[ta,ma,ea] = RunningAverage(stims,tAmplitude*1000,'window',window,'overlap',0.5*window);
	xlim([stims(1,1) stims(end,1)]);
	PlotMean(ta,ma,ea(:,1),ea(:,2),':','k');
	PlotHVLines(slices,'v','color',[0.8 0.8 0.8]);
	ylabel('Amplitude time (ms)');
	xlabel(xLabel);

	% ISIs
	subplot(2,3,3);
	ds = diff(stims0(:,1));
	x = 0:0.1:10;
	h = hist(ds,x);
	h = h/sum(h);
	bar(x,h);
	xlabel('ISI');
	xlim([0 10]);
	ylabel('Count');
	% CV and CV2
	binSize = 0.1;
	smooth = 10;
	for i = 1:20,
		cv(i) = CV(stims0,'order',i,'binSize',binSize,'smooth',smooth,'method','fixed');
	end
	cv2 = CV(stims0,'measure','cv2');
	subplot(2,3,6);
	plot(cv,'+-');
	hold on;
	plot(cv2,'r+');
	xlabel('nth ISI');
	ylabel('CV & CV2');

	% Stims autocorrelograms
	figure;
	[ccg,x,y] = ShortTimeCCG(stims0,'min',1,'mode','norm','smooth',[2 0]);
	PlotShortTimeCCG(ccg,'x',x,'y',y);
	clim([0 0.01]);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% EMGSignal.m
% This code is about the EMG signal amplitude in each session for all subjects.
% Please set the download directory of your computer for visualization.
% Please add the bbci toolbox in the 'Reference_toolbox' folder

edit EMGSignal_Session1
edit EMGSignal_Session2
edit EMGSignal_Session3

author: PatternRecognition

repo: OpenBMI


text:
classdef MyODEMR < ALGORITHM
% <multi/many> <real/integer>
% Many-objective differential evolution with mutation restriction
% nP --- 500 --- Number of reference points for IGD calculation

%------------------------------- Reference --------------------------------
% R. Denysiuk, L. Costa, and I. E. Santo, Many-objective optimization using
% differential evolution with variable-wise mutation restriction,
% Proceedings of the Annual Conference on Genetic and Evolutionary
% Computation, 2013, 591-598.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Roman Denysiuk

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            nP = Algorithm.ParameterSet(500);

            %% Generate hyperplane
            P = UniformPoint(nP,Problem.M);

            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                Offspring  = Operator(Problem,Population(1:Problem.N),Population(randi(Problem.N,1,Problem.N)),Population(randi(Problem.N,1,Problem.N)));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N,P);
            end
        end
    end
end

author: BIMK

repo: PlatEMO


text:
function [norm_im, mask] = segmentation(or_im, blksze, thr)

    fun = inline('std(x(:))*ones(size(x))');

    std_devim = blkproc(or_im, [blksze blksze], fun);

    mask = std_devim > thr;
    mask_i = find(mask);

    or_im = or_im - mean(or_im(mask_i));
    norm_im = or_im / std(or_im(mask_i));

    %imshow(normim);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [dx, dxGrad] = dynamics(x,u,p)
% [dx, dxGrad] = dynamics(x,u,p)
%
% Computes the dynamics (and gradients) for the simple pendulum
%

q = x(1,:);
dq = x(2,:);

    k = p.k;    c = p.c;
    ddq = -c*dq - k*sin(q) + u;
    dx = [dq;ddq];

if nargout == 2   % Analytic gradients
    nTime = length(u);

    dqGrad = zeros(1,4,nTime); %4 = [time + angle + rate + torque];
    dqGrad(1,3,:) = 1; %gradient dq wrt dq

    ddqGrad = zeros(1,4,nTime);  %4 = [time + angle + rate + torque];
    ddqGrad(1,2,:) = -k*cos(q);   %gradient ddq wrt q
    ddqGrad(1,3,:) = -c;  %gradient ddq wrt dq
    ddqGrad(1,4,:) = 1;  %gradient ddq wrt u

    dxGrad = cat(1, dqGrad, ddqGrad);

end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function ok = dfcanplotdata(ds,dffig)
%DFCANPLOTDATA Determine if we can plot data in the current plot

%   $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:35:19 $
%   Copyright 2003-2004 The MathWorks, Inc.

% No problem unless we have non-positive data
ok = true;
xlimits = xlim(ds);
xmin = xlimits(1);
if xmin>0
   return
end

% Get handle to control containing the distribution list
if nargin<2
   dffig = dfgetset('dffig');
end
hsel = getappdata(dffig,'selectioncontrols');

% No problem unless we have a probability plot
h = hsel(3);   % handle of display type control
choice = get(hsel(3),'Value');
ftypes = getappdata(hsel(3),'codenames');
ftype = ftypes{choice};
if ~isequal(ftype, 'probplot')
   return
end

% No problem unless distribution has support that excludes some data
ax = get(dffig,'CurrentAxes');
distspec = getappdata(ax,'DistSpec');
lobnd = distspec.support(1);
strict = ~distspec.closedbound(1);
if strict && lobnd>=xmin
   ok = false;
elseif ~strict && lobnd>xmin
   ok = false;
end


text:
function r = mpower(p,q)
% MEAS/MPOWER  Implement p^q for meas.

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

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


text:
function res = Profile_Horn_res(p,  doseProfile1, xIndexStart, xIndexEnd, yIndex, s, ICsigma, hornOffAxis, fractionEF, varargin)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% p(1) := calibration constant = 25. (to scale DPM dose to 1cGy/MU at
% depth = 1.5cm for field size 10x10cm2, at 100SSD.

% p(2) := horn at atan( 2.5/100);
% p(3) := horn at atan( 5.0/100);
% p(4) := horn at atan( 7.5/100);
% p(5) := horn at atan(10.0/100);
% p(6) := horn at atan(15.0/100);
% p(7) := horn at atan(20.0/100);
% p(8) := horn at atan(28.5/100);

if length(varargin) == 2
    IM = varargin{1};
    FS = varargin{2};
else
    % When called by 'fimsearch(FUN, x0)', varargin's elements turns into cell from struct.
    IM = varargin{1}{1};
    FS = varargin{1}{2};
end

% Get the horn based on
hornCoef = [0, p(2:end)];
offAxis = [0:0.1:28.5];
test = interp1(hornOffAxis, hornCoef, offAxis);

% get off-axis-distance at 100 cm SSD.
PBOffAxis = sqrt(IM.beams.xPBPosV.^2 + IM.beams.yPBPosV.^2);
w_field = 1 + interp1(offAxis, test, PBOffAxis);
% Correction for isotropic derive
% IM2isoCorrection = (cos(atan(PBOffAxis/IM.beams.isodistance)).^2);
% w_field = w_field .* IM2isoCorrection;

% Now this is only Primary photon contribution.
IM.beamlets = IM.beamletsPrimary;
tic; [dose3D] = getIMDose(IM, w_field, 1); toc
% Use Issam's code to do the denoising.
% Denoising part.
%! dose3D = max(dose3D(:))*anisodiff3d_miao_exp(dose3D/max(dose3D(:)), 0.03, 4);

 % Get the Extra-focal contribution.
IM.beamlets = IM.beamletsFF;
tic; [dose3D_FF] = getIMDose(IM,[], 1); toc
% Denoising.
%! dose3D_FF = max(dose3D_FF(:))*anisodiff3d_miao_exp(dose3D_FF/max(dose3D_FF(:)), 0.03, 4);

% Add dose3D and dose3D_FF together.
dose3D = dose3D + fractionEF * dose3D_FF;

% Scale dose up by p(1) to match with the measurements.
doseProfile1_DPM = p(1)*dose3D(yIndex, xIndexStart:xIndexEnd, int16(s));

%Start the function.
%Calculate the residue of the PDD and dose profile difference between DPM and the measurement.
diffsqr = sum((doseProfile1 - doseProfile1_DPM).^2);
res = sqrt(diffsqr)
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [theta,HH,C,P,hrf,fit,e,param] = Anneal_Logit(theta0,t,tc,Run)
%
% [theta,HH,C,P] = Anneal_Logit(theta0,t,tc,Run)
%
% Estimate inverse logit (IL) HRF model using Simulated Annealing
% Creates fitted curve - 3 logistic functions to be summed together - from parameter estimates
%
% INPUT: theta0, t, tc, Run
% Run = stick function
% tc = time course
% t = vector of time points
% theta0 = initial value for the parameter vector
%
% By Martin Lindquist and Tor Wager
% Edited 12/12/06
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Initial values

iter = 15000;                               % Number of iterations
theta = theta0;                             % Set initial value for the parameter vector
h0 = cost(theta0,t,tc,Run);                 % Calculate cost of initial estimate
LB = [0.05, 1, 0, 0.05, 5, 0, 10];      % Lower bounds for parameters
UB = [10, 15, 5, 10, 15, 5, 30];           % Upper bounds for parameters

%
% These values may need tweaking depending on the individual situation.
%

r1= 0.001;                                 % A parameters
r1b= 0.001;                                % A parameters
r2 = 0.05;                                 % T parameters
r3 = 0.001;                                % delta parameters

t1 = [1 4];
t1b = [6];
t2 = [2 5 7];
t3 = [3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

u = zeros(1,7);

HH = zeros(1+iter,7);               % Keep track of theta_i
HH(1,:) = theta0;
P = zeros(1+iter,1);
C = zeros(1+iter,1);                % Keep track of the cost function
C(1) = h0;

cnt = 0;
for i=1:iter,

    T = 100/log(1+i);    %Temperature function (may require tweaking)
    th = zeros(1,7);
    ind = 0;

    % Choose a new candidate solution theta_{i+1}, based on a random perturbation of the current solution of theta_{i}.
    % Check new parameters are within accepted bounds
    while ( (sum((LB-th)>0) + sum((th-UB)>0)) > 0),

        % Perturb solution

        u(t1) = normrnd(0,r1,1,2);
        u(t1b) = normrnd(0,r1b,1,1);
        u(t2) = normrnd(0,r2,1,3);
        u(t3) = normrnd(0,r3,1,1);

        % Update solution
        th = theta + u;
        ind = ind + 1;

        if(ind > 500),
            warning('stuck!');
            return;
        end;
    end;

    h = cost(th,t,tc,Run);
    C(i+1) = h;
    delta = h - h0;

    % Determine whether to update the parameter vector.
    if (unifrnd(0,1) < min(exp(-delta/T),1)),
        theta = th;
        h0=h;
        cnt = cnt+1;
    end;

    HH(i+1,:) = theta;
    P(i+1) = min(exp(-delta/T),1);

end;

%cnt/iter

[a,b] = min(C);
theta = HH(b,:);
%h

% Additional outputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get HRF for final model
if nargout > 4
    hrf = Get_Logit(theta(1:7),t);                   % Calculate HRF estimate (fit, given theta)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convolve HRF and stick function
if nargout > 5
    len = length(Run);
    fit = conv(Run, hrf);
    fit = fit(1:len);
    e = tc - fit;
end

if nargout > 7
    [param] = get_parameters_logit(hrf,t,theta);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [c1,c2]=lox(p1,p2,bounds,genInfo,ops)
% Linearorder crossover takes two parents P1,P2 and performs linear order
% crossover for permutation strings.
%
% function [c1,c2] = linearOrderXover(p1,p2,bounds,Ops)
% p1      - the first parent ( [solution string function value] )
% p2      - the second parent ( [solution string function value] )
% bounds  - the bounds matrix for the solution space
% Ops     - Options matrix for simple crossover [gen #SimpXovers].

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% This program 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 1, or (at your option)
% any later version.
%
% This program 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
sz=size(p1,2)-1;
c1=p1(1:sz);%zeros(1,sz);
c2=p2(1:sz);%zeros(1,sz);
cut1=round(rand*(sz-1)+1.5);
cut2=round(rand*(sz-1)+1.5);
while cut2 == cut1
   cut2 = round(rand*sz + 0.5);
end
if cut1 > cut2
   t = cut1; cut1 = cut2; cut2 = t;
end
for i=cut1:cut2
  c1=strrep(c1,p2(i),-1);
  c2=strrep(c2,p1(i),-1);
end
g1=find(c1>-1);
g2=find(c2>-1);
c1=[c1(g1(1:(cut1-1))) p2(cut1:cut2) c1(g1(cut1:end)) p1(end)];
c2=[c2(g2(1:(cut1-1))) p1(cut1:cut2) c2(g2(cut1:end)) p2(end)];

author: HuangCongQing