text:
function vw = motionCompSetArgs(vw)
%
%   gb 14/04/05
%
%   vw = motionCompSetArgs(vw)
%
% Runs an interface to set the arguments of motionCompMutualInfMeanInit
% Also computes the MSE and the MI between consecutive frames of the new
% dataType.
%
global dataTYPES HOMEDIR

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creates the figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nScans = numberScans(vw);
nTrans = 4;

xs = 1.8;
ys = 1.8;

botMargin = 0.2;
butWidth = 8;
editWidthY = 1.4;
fontSize = 9;

height11 = nTrans + 2 + botMargin*2 + .5;
height21 = nScans + 2 + botMargin*2 + .5;
heightEdge = 1;
heightSep = 2;
height = height11 + height21 + heightSep + 2*heightEdge;

width1 = 22;
width2 = 22;
width3 = 35;
widthAlgo = width1;
widthSep = 5;
widthEdge = 1;
width = width1 + width2 + width3 + 2*widthSep + 2*widthEdge;

h = figure('MenuBar','none',...
    'Name','Choose parameters of the transformation',...
    'Units','char',...
    'Resize','off',...
    'NumberTitle','off',...
    'Position',[60,(60 - height*ys), width*xs,height*ys]);

hControl = cell(1,9);

scanNames = {};
for i = 1:nScans
    scanNames = [scanNames,{['Scan' num2str(i)]}];
end
hControl{1} = fillFigure('Scans',scanNames,h,widthEdge,height11 + heightSep + heightEdge,width1);

hControl{2} = fillFigure('Algorithms',{'Rigid (MI)','Non linear (MSE)','Apply GLM','Consecutive frame error'},...
    h,widthEdge,heightEdge,widthAlgo,'',0);

scanNames = [scanNames,{'Anatomy'}];
x = widthEdge + width1 + widthSep;
y = height11 + heightSep + heightEdge + nScans;
fillFigure('BaseScan',{},h,x,y,width2,'',0);
hControl{3} = uicontrol('Style','popupmenu',...
                'Units','char',...
                'String',scanNames,...
                'Position',[(x + widthSep)*xs,(y)*xs,(width1 - 2*widthSep)*xs,editWidthY],...
                'HorizontalAlignment','center',...
                'FontSize',fontSize);

y = y - nScans;
fillFigure('Original data type name',{},h,x,y,width2,'',0);
dataNames = {};
for i = 1:length(dataTYPES)
    dataNames = [dataNames,{dataTYPES(i).name}];
end

hControl{4} = uicontrol('Style','popupmenu',...
                'Units','char',...
                'String',dataNames,...
                'Position',[(x + widthSep)*xs,(y)*xs,(width1 - 2*widthSep)*xs,editWidthY],...
                'HorizontalAlignment','center',...
                'FontSize',fontSize);

y = heightEdge + nTrans;
fillFigure('New data type name',{},h,x,y,width2,'',0);

hControl{5} = uicontrol('Style','edit',...
                'Units','char',...
                'String','',...
                'Position',[(x + widthSep)*xs,(y)*xs,(width1 - 2*widthSep)*xs,editWidthY],...
                'HorizontalAlignment','center',...
                'FontSize',fontSize,...
                'Enable','Off');

x = widthEdge + width1 + width2 + 2*widthSep;
y = height11 + heightSep + heightEdge + nScans;
fillFigure('ROI',{},h,x,y,width3,'',0);

ROInames = {'None','ROIocc','ROIfront'};
for i = 1:length(vw.ROIs)
    ROInames = [ROInames,{vw.ROIs(i).name}];
end
ROIfiles = dir(roiDir(vw));
for i = 3:length(ROIfiles)
    if isempty(find(strcmp(ROInames,ROIfiles(i).name(1:end - 4)), 1))
        ROInames = [ROInames, {ROIfiles(i).name(1:end - 4)}];
    end
end

hControl{6} = uicontrol('Style','popupmenu',...
                'Units','char',...
                'String',ROInames,...
                'Position',[(x + widthSep)*xs,(y)*xs,(width1 - 2*widthSep)*xs,editWidthY],...
                'HorizontalAlignment','center',...
                'FontSize',fontSize);

x = width - widthEdge - width3/2;
nSlices = numberSlices(vw);
hControl{7} = createSlider('Slice',[(x + widthSep/2)*xs y*ys (width3/2 - 3/2*widthSep)*xs editWidthY],[1 nSlices]);

x = widthEdge + width1 + width2 + width3/4;
y = heightEdge + heightSep;
hControl{8} = gca;
set(hControl{8},'Units','char');
set(hControl{8},'Position',[(x + widthSep)*xs,(y)*xs,(width3 - 2*widthSep)*xs,(nScans + nTrans + heightSep - 2*botMargin)*ys]);

str = [vw.ui.displayMode,'Mode'];
modeStr=['vw.ui.',str];
mode = eval(modeStr);
cmap = mode.cmap;
numGrays = mode.numGrays;
numColors = mode.numColors;
clipMode = mode.clipMode;

% image(vw.ui.image);
% colormap(cmap);
% axis image;
% axis off;

x = widthEdge + width1 + width2 + width3 - butWidth - botMargin;
y = heightEdge + botMargin;
uicontrol('Style','pushbutton',...
    'String','Cancel',...
    'Units','char',...
    'Position',[x*xs,y*ys,butWidth*xs,ys],...
    'CallBack','uiresume',...
    'FontSize',fontSize,...
    'UserData','Cancel');

x = x - butWidth - widthEdge;
hControl{9} = uicontrol('Style','pushbutton',...
                'String','OK',...
                'Enable','Off',...
                'Units','char',...
                'Position',[x*xs,y*ys,butWidth*xs,ys],...
                'CallBack','uiresume',...
                'FontSize',fontSize,...
                'UserData','OK');

hControl{10} = 0;

for i = 1:(length(hControl) - 3)
    for j = 1:length(hControl{i})
        set(hControl{i}(j),'CallBack','uiresume');
        set(hControl{i}(j),'UserData','Update');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main Loop
%
% Update the arguments whenever the user adds
% something to the figure
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

response = 'Update';
while strcmp(response,'Update')
    uiwait

    % determine which button was hit.
    response = get(gco,'UserData');

    if isempty(response)
        close(get(0,'currentFigure'));
        fprintf('Motion Compensation Canceled by the user');
        return;
    end

    try
        if strcmp(get(gco,'Style'),'edit')
            hControl{end} = 1;
        end
    end

    % Update the arguments
    [scans,ROI,baseScan,newDataType,currentDataType,rigid,nonLinear,glmIndex,consError] = updateHControl(vw,hControl);

    % Update the new slice showed. If an ROI has been loaded, loads a slice
    % where the ROI is not zero.
    if ~ieNotDefined('ROI')
        if sum(sum(abs(ROI(:,:,get(hControl{7}(1),'Value'))))) == 0
            mn = mean(mean(ROI,1),2);
            displaySlices = find(mn > 0);
            [minimum,index] = min(abs(displaySlices - get(hControl{7}(1),'Value')));
            if ~isempty(displaySlices)
                set(hControl{7}(1),'Value',displaySlices(index));
                set(hControl{7}(2),'String',num2str(displaySlices(index)));
            else
                msgBox('This ROI is empty!','Warning!');
            end
        end
    end

    vw = viewSet(vw, 'Current Slice', get(hControl{7}(1),'Value'));

    % Recomputes the shown image
    vw = recomputeImage(vw,numGrays,numColors,clipMode);

    size1 = sliceDims(vw,1);
    size2 = size(vw.anat);
    T = maketform('affine',[size2(2)/size1(2) 0 0; 0 size2(1)/size1(1) 0; 0 0 1]);

    im = vw.ui.image;

    if ~ieNotDefined('ROI')

        im = cmap(double(im)+1);
        roi = ROI(:,:,str2num(get(hControl{7}(2),'String')));
        roi = imtransform(roi,T,'nearest');
        roi = [zeros(1,size(roi,2) + 1);zeros(size(roi,1),1),roi];

        im = repmat(im,[1 1 3]);
        roi = repmat(roi,[1 1 3]);
        color = [1 0 0];
        colorIndex = repmat(reshape(color,[1 1 length(color)]),[size(im,1) size(im,2) 1]);

        transparency = 2/3;
        im = (im + transparency)/(1 + transparency).*roi.*colorIndex + im.*(1 - roi);

    end
	image(im);
    colormap(cmap);
    axis image;
    axis off;

end

if ~ieNotDefined('ROI')
    ROIname = ['ROI_' newDataType];
end
close(h)
vw = refreshScreen(vw);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Once the ok button is pressed, the motion compensation algorithm is
% called. It can be followed by the glm and the computation of the
% consecutive frame errors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(response,'OK')
    vw = motionCompMutualInfMeanInit(vw,scans,ROI,baseScan,newDataType,currentDataType,rigid,nonLinear);

    % Runs the GLM
    if glmIndex
        er_runSelxavgBlock(vw,intersect(union(baseScan,scans),1:6),1);
    end

	% Creates the movie of the mean maps
%     try,	motionCompMeanMovie(vw);  end

    % Computes the consecutive-frame error
    if consError,
        try
            for i = 1:length(dataTYPES)
                if isequal(dataTYPES(i).name,currentDataType)
                    vw = viewSet(vw,'currentDataType',i);
                    break
                end
            end

            motionCompPlotMSE(vw,'',0);
            motionCompPlotMI(vw,'',0);

            if ~ieNotDefined('ROIname')
                motionCompPlotMSE(vw,ROIname,0);
                motionCompPlotMI(vw,ROIname,0);
            end

            if exist(fullfile(roiDir(vw),'ROIdef.mat'),'file')
                pathMSE = fullfile(dataDir(vw),['MSE_' vw.sessionCode '_' currentDataType '_ROIdef.mat']);
                if exist(pathMSE,'file')
                    delete(pathMSE);
                end
                pathMI = fullfile(dataDir(vw),['MI_' vw.sessionCode '_' currentDataType '_ROIdef.mat']);
                if exist(pathMI,'file')
                    delete(pathMI);
                end

                motionCompPlotMSE(vw,'ROIdef',0);
                motionCompPlotMI(vw,'ROIdef',0);
            end

            for i = 1:length(dataTYPES)
                if isequal(dataTYPES(i).name,newDataType)
                    vw = viewSet(vw,'currentDataType',i);
                    break
                end
            end

            motionCompPlotMSE(vw,'',0);
            motionCompPlotMI(vw,'',0);

            if ~ieNotDefined('ROIname')
                motionCompPlotMSE(vw,ROIname,0);
                motionCompPlotMI(vw,ROIname,0);
            end

            if exist(fullfile(roiDir(vw),'ROIdef.mat'),'file')
                pathMSE = fullfile(dataDir(vw),['MSE_' vw.sessionCode '_' newDataType '_ROIdef.mat']);
                if exist(pathMSE,'file')
                    delete(pathMSE);
                end
                pathMI = fullfile(dataDir(vw),['MI_' vw.sessionCode '_' newDataType '_ROIdef.mat']);
                if exist(pathMI,'file')
                    delete(pathMI);
                end

                motionCompPlotMSE(vw,'ROIdef',0);
                motionCompPlotMI(vw,'ROIdef',0);

            end
        end
    end
end

return

% /-----------------------------------------------------------------------/ %
function h_button = fillFigure(headerStr,optionStr,h,xoff,yoff,width,height,displayButtons)
% function reply = fillFigure(headerStr,optionStr)
%
% Button Dialog Box
%
%  Allows the user to toggle and select options in
%  the cell array string 'optionStr'.  Reply is a
%  boolean vector with length(optionStr) that indicates
%  selected options. reply is empty if the user
%  chooses to cancel.
%
%  Example:
%  reply = buttondlg('pick it',{'this','that','the other'})
%
%4/11/98 gmb   Wrote it
% 15/07/02  fwc     it will now make multiple columns of options if the number is large
% 2002.07.24  rfd & fwc: fixed minor bug where short checkbox labels were invisible.
% 2004.12.16  arw: Added option to select / deselect all
%
% gb 04/18/05
%
% Several options are added in order to make it more flexible and to
% integrate this window to an existing figure

OptionsPerColumn=11; % max number of options in one column

if nargin<2
    disp('Error: "bottondlg" requires two inputs');
    return
end

if ieNotDefined('xoff')
    xoff = 0;
end

if ieNotDefined('yoff')
    yoff = 0;
end

if ieNotDefined('displayButtons')
    displayButtons = 1;
end

if isunix
    fontSize = 10;
else
    fontSize = 9;
end

if strcmp(class(optionStr),'cell')
    optionStr = char(optionStr);
end

nOptions = size(optionStr,1);

ncols=ceil(nOptions/OptionsPerColumn);
nOptionsPerColumn=max(ceil(nOptions/ncols),0);

%scale factors for x and y axis coordinates
xs = 1.8;
ys = 1.8;

%default sizes
butWidth=10;
botMargin = 0.2;

if ieNotDefined('height')
    height = nOptionsPerColumn+2+botMargin*2+.5;
end

% If we don't have a minimum colwidth (5 in this case), short strings don't show up.
colwidth=max(size(optionStr,2),5);%

if ieNotDefined('width')
    width = ncols*colwidth+2;
    width = max(width,butWidth+2);
end

%open the figure
if ieNotDefined('h')
    h = figure('MenuBar','none',...
    'Units','char',...
    'Resize','off',...
    'NumberTitle','off',...
    'Position',[20, 10, width*xs,height*ys]);
end
bkColor = get(h,'Color');

%Display title text inside a frame
x = 1;
% y = nOptions+1+botMargin;
y = nOptionsPerColumn+1+botMargin;

x = x + xoff;
y = y + yoff;

uicontrol('Style','frame',...
    'Units','char',...
    'String',headerStr,...
    'Position',[x*xs,(y+.2)*ys,(width-2)*xs,ys*1.3],...
    'HorizontalAlignment','center',...
    'FontSize',fontSize);

uicontrol('Style','text',...
    'Units','char',...
    'String',headerStr,...
    'Position',[(x+.25)*xs,(y+.3)*ys,(width-2.5)*xs,ys*.9],...
    'HorizontalAlignment','center',...
    'FontSize',fontSize);

%Display the radio buttons
y = y+botMargin/2;
y0=y;
c=0;

for optionNum=1:nOptions
    y = y-1;
    h_button(optionNum) = ...
        uicontrol('Style','checkbox',...
        'Units','char',...
        'String',optionStr(optionNum,:),...
        'BackgroundColor',bkColor,...
        'Position',[x*xs+c*colwidth*xs,y*ys,colwidth*xs+(c+1)*colwidth*xs,ys],...
        'HorizontalAlignment','left',...
        'FontSize',fontSize);
    if optionNum>=(c+1)*nOptionsPerColumn
        c=c+1;
        y=y0;
    end
end

if displayButtons
	%Display the OK/Cancel buttons
	x = width - butWidth - 3;
	x = x + xoff;
	y = botMargin + yoff;

	sabutton=uicontrol('Style','pushbutton',...
        'String','SelAll',...
        'Units','char',...
        'Position',[x*xs+2,y*ys,butWidth*xs/2,ys],...
        'FontSize',fontSize-4,...
        'Callback',['set([',num2str(h_button,30),'],''Value'',1);uiresume'],...
        'UserData','Update');
	cabutton=uicontrol('Style','pushbutton',...
        'String','ClrAll',...
        'Units','char',...
        'Position',[x*xs+2+butWidth,y*ys,butWidth*xs/2,ys],...
        'FontSize',fontSize-4,...
        'Callback',['set([',num2str(h_button,30),'],''Value'',0);uiresume'],...
        'UserData','Update');
end

return;

% /------------------------------------------------------------------------------/ %
function [scans,ROI,baseScan,newDataType,currentDataType,rigid,nonLinear,glmIndex,consError] = updateHControl(vw,hControl)
%
% Updates the arguments from the current figure
%

nScans = length(hControl{1});
scans = [];
for i = 1:nScans
    if get(hControl{1}(i),'Value')
        scans = [scans i];
    end
end

scanNames = get(hControl{3},'String');
baseScan = get(hControl{3},'Value');

dataNames = get(hControl{4},'String');
currentDataType = dataNames{get(hControl{4},'Value')};
newDataType = get(hControl{5},'String');

rigid = get(hControl{2}(1),'Value');
nonLinear = get(hControl{2}(2),'Value');
glmIndex = get(hControl{2}(3),'Value');
consError = get(hControl{2}(4),'Value');

if gco == hControl{7}(2)
    val = str2num(get(hControl{7}(2),'String'));
    val = min(get(hControl{7}(1),'Max'),val);
    val = max(get(hControl{7}(1),'Min'),val);
    val = round(val);

    set(hControl{7}(1),'Value',val);
    set(hControl{7}(2),'String',num2str(val));
else
    val = get(hControl{7}(1),'Value');
    val = min(get(hControl{7}(1),'Max'),val);
    val = max(get(hControl{7}(1),'Min'),val);
    val = round(val);

    set(hControl{7}(1),'Value',val);
    set(hControl{7}(2),'String',num2str(val));
end

if (rigid | nonLinear) & ~isempty(setdiff(scans,baseScan))
    set(hControl{5},'Enable','On');
    set(hControl{end - 1},'Enable','On');

    if hControl{end} == 0;
        newDataName = '';
        if rigid
            newDataName = [newDataName 'Rigid'];
            if nonLinear
                newDataName = [newDataName '+'];
            end
        end

        if nonLinear
            newDataName = [newDataName 'NL'];
        end

        newDataName = [newDataName '_' scanNames{baseScan}];

        set(hControl{5},'String',newDataName);
    end

else
    set(hControl{5},'Enable','Off');
    set(hControl{5},'String','');
    set(hControl{end - 1},'Enable','Off');
    hControl{end} = 0;
end

ROINames = get(hControl{6},'String');
switch ROINames{get(hControl{6},'Value')}

    case 'None'
        ROI = '';

    case 'ROIocc'
        ROI = zeros([sliceDims(vw,1) numberSlices(vw)]);
        ROI((end - 50):(end - 20),30:(end - 30),15:19) = ones(size(ROI((end - 50):(end - 20),30:(end - 30),15:19)));

    case 'ROIfront'
        ROI = zeros([sliceDims(vw,1) numberSlices(vw)]);
        ROI((20:end - 20),30:(end - 30),2:8) = ones(size(ROI((20:end - 20),30:(end - 30),2:8)));

    otherwise
        ROI = motionCompGetROI(vw,ROINames{get(hControl{6},'Value')});

end

if baseScan == (nScans + 1)
    baseScan = 0;
end

return

% /-----------------------------------------------------------------------------/ %
function handle = createSlider(name,position,range)
%
% Creates a slider at a position and with a range specified by the user
%

color = get(gcf,'Color');

% Make slider
sliderHandle = ...
    uicontrol('Style','slider',...
    'Units','char',...
    'Position',position,...
    'min',range(1),...
    'max',range(2),...
    'SliderStep',[1/(range(2) - range(1)) 1/(range(2) - range(1))],...
    'val',floor(range(2)/2));

% Make label
labelOffset = [2,1.8,0,0];
labelPos = position - labelOffset;
labelPos(4) = 1.3;     % don't encroach on the slider

uicontrol('Style','text',...
          'Units','char',...
          'Position',labelPos,...
          'BackgroundColor',color,...
          'String',[name ':']);

editOffset = [7 - position(3),1.8,0,0];
editPos = position - editOffset;
editPos(3) = 5;
editPos(4) = 1.5;    % don't encroach on the slider
editHandle = ...
    uicontrol('Style','edit',...
    'Units','char',...
    'Position',editPos,...
    'String',num2str(get(sliderHandle,'Value')),...
    'BackgroundColor',color);

handle = [sliderHandle editHandle];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Copyright 2013 The MathWorks, Inc
%% Setup Visuals
if m == 1
    sScope = dsp.TimeScope('SampleRate',fs,...
        'TimeSpan',10/prf,'Grid',true,...
        'LayoutDimensions',[1 2],'MaximizeAxes','off', ...
        'Position',[371 646 1000 409],'NumInputPorts',2);
    set(sScope,'ActiveDisplay',2,'YLabel','Magnitude','Title','Collected Signal')
    set(sScope,'ActiveDisplay',1,'YLabel','Magnitude','Title','Transmitted Signal')
    show(sScope);
    %% Range Doppler Map
    figure('WindowStyle','docked')
    hrdmap = imagesc(sgrid,rgrid,abs(rdmap));
    xlabel('Speed (m/s)'); ylabel('Range (m)'); title('Range Doppler Map');

    %% Detection and Range Estimation
    figure('WindowStyle','docked')
    tgtrange = [NaN NaN NaN]; pmax = tgtrange;
    hold on
    for n=1:3,
        plot(sTgtMotion{n}.InitialPosition(1)*ones(2,1),[0 7e-5],'r:')
        htext(n) = text(tgtrange(n),1.05*pmax(n),int2str(tgtrange(n))); %#ok<*SAGROW>
    end
    legend('Initial Range')
    xlabel('Range (m)'); ylabel('Magnitude'); title ('Estimated Range')
    hbar = sqrt(threshold)*ones(numel(fast_time),1);
    hline = plot(range_gates,[hbar hbar]); % Threshold
    offset = numel(sMFilt.Coefficients)-1;
end

%% Stream Signals
step(sScope,abs(s),abs(rsig)); % Ctrl + A to scale axis limits
set(hrdmap,'CData',abs(fliplr(rdmap)))
drawnow
set(hline(2),'YData',[intpulses(offset:numel(fast_time)); NaN*ones(offset-1,1)])
for n=1:3,
    set(htext(n),'String',int2str(tgtrange(n)),'Position',[tgtrange(n) 1.05*pmax(n)]);
end

author: Sable

repo: mcbench-benchmarks


text:
%% Mean-Shift Video Tracking
% by Sylvain Bernhardt
% July 2008
%% Description
% Insert in the image I a rectangle
% which size is H,W and thickness is
% expressed in pixels.
%
% I = Draw_target(x,y,W,H,I,thick)

function I = Draw_target(x,y,W,H,I,thick)
R=1; G=2; B=3;
I(y:y+thick-1,x:x+W,R) = 255;
I(y:y+thick-1,x:x+W,G) = 0;
I(y:y+thick-1,x:x+W,B) = 0;

I(y+H-thick+1:y+H,x:x+W,R) = 255;
I(y+H-thick+1:y+H,x:x+W,G) = 0;
I(y+H-thick+1:y+H,x:x+W,B) = 0;

I(y:y+H,x:x+thick-1,R) = 255;
I(y:y+H,x:x+thick-1,G) = 0;
I(y:y+H,x:x+thick-1,B) = 0;

I(y:y+H,x+W-thick+1:x+W,R) = 255;
I(y:y+H,x+W-thick+1:x+W,G) = 0;
I(y:y+H,x+W-thick+1:x+W,B) = 0;


text:
function [x0, lb, ub] = ma_formOptim_X0_LB_UB(variables)
%ma_formOptim_X0_LB_UB Summary of this function goes here
%   Detailed explanation goes here

    lb = [];
    ub = [];
    x0 = [];
    for(varC = variables) %#ok<*NO4LP>
        var = varC{1};
        lb = [lb,var.lb]; %#ok<*AGROW>
        ub = [ub,var.ub];
        x0 = [x0,var.x0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [varargout] = minus(varargin)

% MINUS compile the missing mex file on the fly

% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);
  mex(mexsrc);
  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was just created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function V = erasure_symbol_merge(W)
cnt = 0;
N = size(W, 2);
for i = N/2 : -1  : 1
    if W(1, i) == W(2, i)
        cnt = cnt + 1;
    else
        break;
    end
end
W_erasure = W(:, N/2 - cnt + 1 : N/2 + cnt);
erasure_probability = sum(W_erasure(1, :));
middle = erasure_probability/2 * ones(2, 2);
W_left = W(:, 1 : N/2 - cnt);
W_right = W(:, N/2 + cnt + 1 : end);
V = [W_left middle W_right];
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
%% THE OPTIMAL RECIPROCAL COLLISION AVOIDANCE (agent_2D_ORCA.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This agent class provides support for the RVO2 library for ORCA-based
% Reciprocal Velocity Obstacle collision avoidance.

% Author: James A. Douthwaite

classdef agent_2D_ORCA < agent_2D_RVO
%% INITIALISE THE AGENT SPECIFIC PARAMETERS
    properties
        % AGENT PARAMETERS
%         radius = 2;          % Declare (independant of the VO varients)
%         neighbourDist = 15;  %(s)
%         maxNeighbours = 10;
        timeHorizon     = 10;
        timeHorizonObst = 10;
        RVO_EPSILON     = 1E-5;

        % PARAMETERS UNIQUE TO RVO2
%         timeHorizon     = 10;
%         timeHorizonObst = 10;

        % ORCA-LINE CONTAINERS
        orcaLines = struct('point',[],'direction',[]);
    end
%%  CLASS METHODS
    methods
        % CONSTRUCTION METHOD
        function obj = agent_2D_ORCA(varargin)
            % This function is to construct the agent object using the
            % object defintions held in the 'objectDefinition' base class.
            % INPUTS:
            % namestr - Assigned name string
            % OUTPUTS:
            % obj     - The constructed object

            % CALL THE SUPERCLASS CONSTRUCTOR
            obj@agent_2D_RVO(varargin);

            % SENSORS
            % Sensor properties will be inherited from 'agent_2D_VO' and
            % 'agent_VO' unless explictly overridden here.

            % //////////////// Check for user overrides ///////////////////
            % - It is assumed that overrides to the properties are provided
            %   via the varargin structure.
            [obj] = obj.ApplyUserOverrides(varargin);
            % /////////////////////////////////////////////////////////////
        end
    end
    % //////////////////////// OPENMAS FUNCTIONS ///////////////////////////
    methods
        % MAPPING OF OPENMAS TO THE RVO LIBRARY
        function [headingVector,speed] = GetAvoidanceCorrection(obj,dt,desiredVelocity,visualiseProblem)
            % This function computes the mapping of the simulation inputs
            % for the RVO library functions and returns the calculated
            % optimal velocity.

            % AGENT KNOWLEDGE (2D)
            [p_i,v_i,r_i] = obj.GetAgentMeasurements();

            % PLOT INPUT VARIABLES
            inputDebugPlot = 0;
            if inputDebugPlot
                figureHandle = figure(1);
                hold on; grid on; box on;
                axis equal;
                [figureHandle] = obj.getObjectScene(knownObstacles,figureHandle);           % Plot the objects
            end

            % GET AGENT DATA
            agentIDs = [obj.MEMORY([obj.MEMORY.type] == OMAS_objectType.agent).objectID];

            % MOVE THROUGH THE PRIORITISED OBSTACLE SET
            agentNeighbours = cell(numel(agentIDs),1);
            agentNeighbourLength = 0;
            for item = 1:numel(agentIDs)
                % Get object data from memory structure
                p_j = obj.GetLastMeasurementByID(agentIDs(item),'position');
                v_j = obj.GetLastMeasurementByID(agentIDs(item),'velocity');
                r_j = obj.GetLastMeasurementByID(agentIDs(item),'radius');

                % NEIGHBOUR CONDITIONS
                neighbourConditionA = item < obj.maxNeighbours;            % Maximum number of neighbours
                neighbourConditionB = norm(p_j) < obj.neighbourDist;       % [CONFIRMED]
                neighbourConditionC = ~any(isnan(v_j));                    % Wait for a valid velocity reading
                if ~neighbourConditionB || ~neighbourConditionC
                    continue
                end

                p_j = p_j + p_i;
                v_j = v_j + v_i;

                % BUILD THE RVO AGENT DESCRIPTION
                agentRef = struct('position',p_j,...
                                  'velocity',v_j,...
                                  'unitDir',obj.nullVelocityCheck(v_j),...
                                  'radius',r_j,...
                                  'id',agentIDs(item));

                % BUILD LIST OF AGENTS
%                 agentNeighbours = horzcat(agentNeighbours,agentRef);
                agentNeighbours{item} = agentRef;
                agentNeighbourLength  = agentNeighbourLength + 1;
            end
            % REMOVE EMPTY CELLS
            agentNeighbours = agentNeighbours(~cellfun('isempty',agentNeighbours));

            % OBSTACLE DATA
            obstacleIDs  = [obj.MEMORY([obj.MEMORY.type] == OMAS_objectType.obstacle).objectID];

            % BUILD THE REQUIRED OBSTACLE STRUCTURES
            obstacleNeighbours = cell(numel(obstacleIDs),1);
            for item = 1:numel(obstacleIDs)
                % Get object data from memory structure
                p_j = obj.GetLastMeasurementByID(obstacleIDs(item),'position');
                v_j = obj.GetLastMeasurementByID(obstacleIDs(item),'velocity');

                % PROCESS THE OBSTACLE'S VERTICES
                obstacleRef.position = p_j + p_i;
                obstacleRef.velocity = v_j + v_i;
%               obstacleRef.geometry.vertices = obstacleRef.geometry.vertices + [agentPosition;0]';
                g_j = obj.GetLastMeasurementByID(obstacleIDs(item),'geometry');

                % [TO-DO] PASS THE OBSTACLE GEOMETRY
                obstacleNeighbours{item} = obj.createRVOobstacleVertexList(g_j);
            end
            % REMOVE EMPTY CELLS
            obstacleNeighbours = obstacleNeighbours(~cellfun('isempty',obstacleNeighbours));

            % COMPUTE THE AGENT'S NEW VELOCITY
            % This calls the libraries 'computeNewVelocity' function
            [obj,avoidanceVelocity] = obj.computeNewVelocity(dt,desiredVelocity,agentNeighbours,obstacleNeighbours);

            % SPECIAL CASE- VELOCITY MAGNITUDE IS ZERO
            speed = norm(avoidanceVelocity);
            headingVector = avoidanceVelocity/speed;
            if isnan(headingVector)
                headingVector = [1;0];  % Retain previous heading
            end

        end
    end
    %% FUNCTIONS IN THE "agent.cpp" FILE
    methods
        % COMPUTE THE AGENT VELOCITY
        function [obj,avoidanceVelocity] = computeNewVelocity(obj,dt,prefVelocity,agentNeighbours,obstacleNeighbours)
            % This function computes the new velocity of this agent from the
            % RVO2 library.

            % DEFINE THE AGENT PARAMETERS
            [p_a,v_a,r_a] = obj.GetAgentMeasurements();

            % MAP THE PREFERRED VELOCITY TO NEW FRAME
            prefVelocity = -prefVelocity;

            % RESET THE OCRA-LINES VECTOR FROM THE LAST TIMESTEP
            obj.orcaLines = [];

            % //////////// GET THE OBSTACLE ORCA LINES ////////////////////
            [obst_orcaLines] = obj.getObstacleORCAlines(p_a,v_a,r_a,obj.timeHorizonObst,obstacleNeighbours);
            obj.orcaLines = [obj.orcaLines,obst_orcaLines];
            numObstLines = numel(obj.orcaLines);

            % ////////////// GET THE AGENT ORCA LINES /////////////////////
            [agnt_orcaLines] = obj.getAgentORCAlines(dt,p_a,v_a,r_a,obj.timeHorizon,agentNeighbours);
            obj.orcaLines = [obj.orcaLines;agnt_orcaLines];

            % //////////////// DEBUG PLOTS ////////////////////////////////
%             if ~isempty(obj.orcaLines)
%                 % PLOT THE ORCA-LINES
%                 h = figure(1);
%                 ax = axes(h);
%                 hold on; grid on; box on;
%                 axis equal;
%
%                 g_a = obj.GEOMETRY.vertices + [p_a;0]';
%                 patch(ax,...
%                     'Faces',obj.GEOMETRY.faces,...
%                     'Vertices',g_a);
%
%                 % The objects vertices
%                 objectVertices = obstacleNeighbours{1};
%                 for i = 1:numel(objectVertices)
%                     scatter(ax,objectVertices(i).point(1),objectVertices(i).point(2),'r*');
%                 end
%
%                 for i = 1:numel(obj.orcaLines)
%                     lineVar = obj.orcaLines(i);
%                     len = 5;
%                     X = [lineVar.point(1);len*lineVar.direction(1)];
%                     Y = [lineVar.point(2);len*lineVar.direction(2)];
%                     plot(ax,X,Y,'b','marker','o','lineWidth',2);
%                 end
%                 close(h);
%             end

            % CALL THE SECOND LINEAR PROGRAM
            [obj,lineObjFail,optimalVelocity] = obj.linearProgram2(obj.orcaLines,obj.v_max,prefVelocity,logical(false));

            % IF SOME ORCA-LINES FAIL
            % This function will generate direction based ORCA lines where
            % the previous lines failed.
            if lineObjFail < numel(obj.orcaLines)
                [obj,optimalVelocity] = obj.linearProgram3(obj.orcaLines,numObstLines,lineObjFail,obj.v_max,optimalVelocity);
            end
            % REASSIGN OUTPUT VELOCITY
%             avoidanceVelocity = optimalVelocity;
            avoidanceVelocity = -optimalVelocity;
        end
        % PROCESS THE VERTEX SET
        function [vertexSet] = createRVOobstacleVertexList(obj,observedObject)
            % The obstacle consideration expects the obstacles to be of
            % polygonal form. Each obstacle is defined by a subset of
            % points that resemble its vertices.

            % We assume that each obstacle defines a set of obstacle
            % vertices. In the RVO2 library, the assumption is made
            % that the vertices are given in counter-clockwise order.

            %/*
            %* Add (polygonal) obstacles, specifying their vertices in counterclockwise
            %* order.
            %*/

            % THE OBSTACLES VERTICES
            observedVertices = observedObject.geometry.vertices;

            % GET THE VERTICES PROJECTED ONTO THE X/Y PLANE
            [planarVertices] = OMAS_graphics.geometryPlanarProjection([0;0;1],...
                                                                      observedObject.position,...
                                                                      observedVertices);
            % SORT THE VERTICES CCW ABOUT THE +ve Z-AXIS
            [verticesCCW] = OMAS_graphics.sortVerticesCCW([0;0;1],...
                                                          observedObject.position,...
                                                          planarVertices);
            % GETTING THE 2D VERTEX SET
            if size(verticesCCW,2) > 2
               verticesCCW = verticesCCW(:,1:2);
            end
            % BUILD THE RVO POLY-GONAL OSTACLE DESCRIPTION
            [vertexSet] = obj.defineRVOVertexObstacleSet(verticesCCW);
        end
        % GET THE OBSTACLE ORCA-LINE VECTOR % <<<<<<<<< THE FINAL PROBLEM
        function [obst_orcalines] = getObstacleORCAlines(obj,pa,va,ra,timeHorizonObst,obstacleNeighbours)
            % This function calculates the the orca-lines for static
            % obstacles using the prinicples of the VO.
            % INPUTS:
            % obj - Agent self-reference
            % dt  - The simulation timestep

            % The obstacle neighbour list is composed of the discrete
            % vertex sets. Each neighbour is a set of vertices with their
            % links to each other specified.

            invTimeHorizonObst = 1/timeHorizonObst;

            % ORCA-LINE CONTAINER
            obst_orcalines = [];

            if isempty(obstacleNeighbours)
                return
            end

            % FOR EACH NEIGHBOURING OBSTACLE
            for n = 1:numel(obstacleNeighbours)
                % THE SUBSET OF VERTICES
                obstacleVertexSet = obstacleNeighbours{n};
                % [JD] - Small notation change. The obstacleVertexSet is
                % now a list of structures that defines the relationship
                % between each vertex.

                % FOR EACH VERTEX IN THE NEIGHBOUR SET
                for i = 1:numel(obstacleVertexSet)
                    % GET THE OBSTACLE
                    obstacle1 = obstacleVertexSet(i);
                    obstacle2 = obstacleVertexSet(obstacle1.nextObstacle);
                    % RELATIVE POSITION OF THE POINTS
                    relativePosition1 = obstacle1.point - pa;
                    relativePosition2 = obstacle2.point - pa;

                    % Check if velocity obstacle of obstacle is already taken care of by
                    % previously constructed obstacle ORCA lines.
                    alreadyCovered = logical(false);

                    if ~isempty(obst_orcalines)
                        for j = 1:numel(obj.orcaLines)
                            conditionA = obj.det(invTimeHorizonObst*relativePosition1 - obst_orcalines(j).point, obst_orcalines(j).direction) - invTimeHorizonObst*ra >= -obj.RVO_EPSILON;
                            conditionB = obj.det(invTimeHorizonObst*relativePosition2 - obst_orcalines(j).point, obst_orcalines(j).direction) - invTimeHorizonObst*ra >= -obj.RVO_EPSILON;
                            if conditionA && conditionB
                                alreadyCovered = logical(true);
                                break;
                            end
                        end
                    end
                    % If covered, move to the next obstacle
                    if alreadyCovered
                        continue
                    end

                    % Obstacle is not yet covered, check for collisions
                    distSq1 = obj.absSq(relativePosition1);
                    distSq2 = obj.absSq(relativePosition2);

                    radiusSq = ra^2;

                    obstacleVector = obstacle2.point - obstacle1.point;
                    s = obj.dot(-relativePosition1,obstacleVector)/obj.absSq(obstacleVector);
                    distSqlineObj = obj.absSq(-relativePosition1 - s*obstacleVector);

                    % ////////////// CREATE ORCA (lineObj) ////////////////
                    lineObj = struct('point',[],'direction',[]);

                    if s < 0 && distSq1 <= radiusSq
                        % Collision with left vertex. Ignore if non-convex
                        if obstacle1.isConvex
                            lineObj.point = [0;0];
                            direction = [-relativePosition1(2);relativePosition1(1)];  % [y;x]
                            lineObj.direction = obj.normalise(direction);
                            obst_orcalines = [obst_orcalines;lineObj];
                        end
                        continue
                    elseif s > 1 && distSq2 <= radiusSq
                        % Collision with right vertex. Ignore if non-convex or
                        % if it will be taken care of by neighbouring obstacle.
                        if obstacle2.isConvex && obj.det(relativePosition2,obstacle2.unitDir) >= 0
                            lineObj.point = [0;0];
                            direction = [-relativePosition2(2);relativePosition2(1)];
                            lineObj.direction = obj.normalise(direction);
                            obst_orcalines = [obst_orcalines;lineObj];
                        end
                        continue
                    elseif s >= 0 && s < 1 && distSqlineObj <= radiusSq
                        % Collision with obstacle segment
                        lineObj.point = [0;0];
                        lineObj.direction = -obstacle1.unitDir;
                        obst_orcalines = [obst_orcalines;lineObj];
                        continue
                    end

%                     /*
%                      * No collision.
%                      * Compute legs. When obliquely viewed, both legs can come from a single
%                      * vertex. Legs extend cut-off line when nonconvex vertex.
%                      */

                    if s < 0 && distSqlineObj <= radiusSq
                        % Obstacle viewed obliquely so that left vertex defines
                        % velocity obstacle
                        if ~obstacle1.isConvex
                            continue;   % Ignore obstacle
                        end
                        obstacle2 = obstacle1;

                        leg1 = sqrt(distSq1 - radiusSq);
                        leftLegDirection  = [relativePosition1(1)*leg1 - relativePosition1(2)*ra;  relativePosition1(1)*ra + relativePosition1(2)*leg1] / distSq1;
                        rightLegDirection = [relativePosition1(1)*leg1 + relativePosition1(2)*ra; -relativePosition1(1)*ra + relativePosition1(2)*leg1] / distSq1; % Assuming [x;y]
                    elseif s > 1 && distSqlineObj <= radiusSq
                        % Obstacle viewed obliquely so that right vertex defines
                        % velocity obstacle.
                        if ~obstacle2.isConvex
                            % Ignore obstacle
                            continue
                        end
                        obstacle1 = obstacle2;

                        leg2 = sqrt(distSq2 - radiusSq);
                        leftLegDirection  = [relativePosition2(1)*leg2 - relativePosition2(2)*ra;  relativePosition2(1)*ra + relativePosition2(2)*leg2] / distSq2;
                        rightLegDirection = [relativePosition2(1)*leg2 + relativePosition2(2)*ra; -relativePosition2(1)*ra + relativePosition2(2)*leg2] / distSq2; % Assuming [x;y]
                    else
                        % Unusual situation
                        if obstacle1.isConvex
                            leg1 = sqrt(distSq1 - radiusSq);
                            leftLegDirection = [relativePosition1(1)*leg1 - relativePosition1(2)*ra; relativePosition1(1)*ra + relativePosition1(2)*leg1] / distSq1;
                        else
                            % Left vertex non-convex; left leg extends cut-off
                            % lineObj
                            leftLegDirection = -obstacle1.unitDir;
                        end

                        if obstacle2.isConvex
                            leg2 = sqrt(distSq2 - radiusSq);
                            rightLegDirection = [relativePosition2(1)*leg2 + relativePosition2(2)*ra; -relativePosition2(1)*ra + relativePosition2(2)*leg2] / distSq2;
                        else
                            % Right vertex non-convex; right leg extends cut-off
                            % lineObj
                            rightLegDirection = obstacle1.unitDir;
                        end
                    end

                    % Legs can never point into neighbouring edge when convex
                    % vertex, take cutoff-lineObj of neighbouring edge instead. If
                    % velocity projected on "foreign" leg, no constraint is
                    % added.

                    leftNeighbour = obstacleVertexSet(obstacle1.prevObstacle);
                    isLeftLegForeign  = logical(false);
                    isRightLegForeign = logical(false);

                    if obstacle1.isConvex && obj.det(leftLegDirection,-leftNeighbour.unitDir) >= 0
                        % Left leg points into obstacle.
                        leftLegDirection = -leftNeighbour.unitDir;
                        isLeftLegForeign = logical(true);
                    end
                    if obstacle2.isConvex && obj.det(rightLegDirection, obstacle2.unitDir) <= 0
                        % Right leg points into obstacle.
                        rightLegDirection = obstacle2.unitDir;
                        isRightLegForeign = logical(true);
                    end
                    % COMPUTE CUT-OFF CENTERS
                    leftCutoff  = invTimeHorizonObst*(obstacle1.point - pa);
                    rightCutoff = invTimeHorizonObst*(obstacle2.point - pa);
                    cutoffVec   = rightCutoff - leftCutoff;

                    % Project current velocity on velocity obstacle.

                    % Check if current velocity is projected on cutoff circles
                    if obstacle1.id == obstacle2.id
                        t = 0.5;
                    else
                        t = obj.dot((va - leftCutoff),cutoffVec) / obj.absSq(cutoffVec);
                    end

                    tLeft  = obj.dot((va - leftCutoff),leftLegDirection);
                    tRight = obj.dot((va - rightCutoff),rightLegDirection);

                    if ((t < 0 && tLeft < 0) || (obstacle1.id == obstacle2.id && tLeft < 0 && tRight < 0))
                        % Project on left cut-off circle.
                        unitW = obj.normalise(va - leftCutoff);

                        lineObj.direction = [unitW(2); -unitW(1)];
                        lineObj.point = leftCutoff + ra*invTimeHorizonObst*unitW;
                        obst_orcalines = [obst_orcalines;lineObj];
                        continue
                    elseif (t > 1 && tRight < 0)
                        % Project on right cut-off circle.
                        unitW = obj.normalise(va - rightCutoff);

                        lineObj.direction = [unitW(2); -unitW(1)];
                        lineObj.point = rightCutoff + ra*invTimeHorizonObst*unitW;
                        obst_orcalines = [obst_orcalines;lineObj];
                        continue;
                    end

                    % Project on left leg, right leg, or cut-off lineObj, whichever is closest
                    % to velocity.
                    % CALCULATE distSqCutoff
                    if t < 0 || t > 1 || obstacle1.id == obstacle2.id
                        distSqCutOff = inf;
                    else
                        distSqCutOff = obj.absSq(va - (leftCutoff + t*cutoffVec));
                    end

                    % CALCULATE distSqLeft
                    if (tLeft < 0)
                        distSqLeft = inf;
                    else
                        distSqLeft = obj.absSq(va - (leftCutoff + tLeft*leftLegDirection));
                    end
                    % CALCULATE distSqRight
                    if (tRight < 0)
                        distSqRight = inf;
                    else
                        distSqRight = obj.absSq(va - (rightCutoff + tRight*rightLegDirection));
                    end

                    if distSqCutOff <= distSqLeft && distSqCutOff <= distSqRight
                        % Project on cut-off lineObj.
                        lineObj.direction = -obstacle1.unitDir;
                        lineObj.point = leftCutoff + ra*invTimeHorizonObst*[-lineObj.direction(2); lineObj.direction(1)];
                        obst_orcalines = [obst_orcalines;lineObj];
                        continue;
                    elseif distSqLeft <= distSqRight
                        % Project on left leg.
                        if isLeftLegForeign
                            continue
                        end
                        lineObj.direction = leftLegDirection;
                        lineObj.point = leftCutoff + ra*invTimeHorizonObst*[-lineObj.direction(2); lineObj.direction(1)];
                        obst_orcalines = [obst_orcalines;lineObj];
                        continue;
                    else
                        % Project on right leg. */
                        if isRightLegForeign
                            continue;
                        end
                        lineObj.direction = -rightLegDirection;
                        lineObj.point = rightCutoff + ra*invTimeHorizonObst*[-lineObj.direction(2); lineObj.direction(1)];
                        obst_orcalines = [obst_orcalines;lineObj];
                        continue;
                    end

                end
            end
        end
        % GET THE AGENT ORCA-LINE VECTOR
        function [agnt_orcalines] = getAgentORCAlines(obj,dt,pa,va,ra,timehorizon,agentNeighbours)
            % This function computes the orca lines using the principles of
            % the RVO method.
            % INPUTS:
            % obj - Agent reference, for use of local functions.
            % dt  - The simulation timestep.
            % pa  - The agents global position
            % va  - The agents global velocity
            % ra  - The agents characteristic radius
            % timeHorizon - The time for which avoidance should be assured.
            % agentNeighbours - The list of considered agent-neighbours.

            % ORCA-LINE CONTAINER
            agnt_orcalines = [];

            if isempty(agentNeighbours)
                return
            end

            % BEGIN CREATING ORCA LINES FOR THE AGENTS
            invTimeHorizon  = 1/timehorizon;

            for i = 1:numel(agentNeighbours)
                % CALCULATE THE OTHER AGENT PROPERTIES
%                 other = agentNeighbours(i);
                %                     relativePosition = other.position - pa;
                %                     relativeVelocity = va - other.velocity;
%                 agentNeighbours{i}
                relativePosition = -(agentNeighbours{i}.position - pa);
                relativeVelocity = -(va - agentNeighbours{i}.velocity);

                distSq = obj.dot(relativePosition,relativePosition);
                combinedRadius = ra + agentNeighbours{i}.radius;
                combinedRadiusSq = combinedRadius^2;

                % CREATE ORCA lineObj
                lineObj = struct('point',[],'direction',[]);

                % NO COLLISION
                if distSq > combinedRadiusSq
                    w = relativeVelocity - invTimeHorizon*relativePosition; % Vector from cutoff center to relative velocity
                    wLengthSq   = obj.dot(w,w); % Square length
                    dotProduct1 = obj.dot(w,relativePosition);

                    if dotProduct1 < 0 && dotProduct1^2 > combinedRadiusSq*wLengthSq
                        % Project on cut-off circle
                        wLength = sqrt(wLengthSq);
                        unitW = w/wLength;

                        lineObj.direction = [unitW(2); -unitW(1)];
                        u = (combinedRadius*invTimeHorizon - wLength)*unitW;

                    else
                        % Project on legs
                        leg = sqrt(distSq - combinedRadiusSq);

                        if obj.det(relativePosition, w) > 0
                            % Project on left leg
                            lineObj.direction = ([relativePosition(1)*leg - relativePosition(2)*combinedRadius;
                                relativePosition(1)*combinedRadius + relativePosition(2)*leg]) / distSq;
                        else
                            % Project on right leg
                            lineObj.direction = -([relativePosition(1)*leg + relativePosition(2)*combinedRadius;
                                -relativePosition(1)*combinedRadius + relativePosition(2)*leg]) / distSq;
                        end

                        dotProduct2 = obj.dot(relativeVelocity,lineObj.direction);
                        u = dotProduct2*lineObj.direction - relativeVelocity;
                    end
                    % COLLISION
                else
                    % Collision. Project on cut-off circle of time timeStep.
                    invTimeStep = 1/dt;
                    % Vector from cutoff center to relative velocity. */
                    w = relativeVelocity - invTimeStep*relativePosition;
                    wLength = sqrt(obj.dot(w,w));
                    unitW = w/wLength;

                    lineObj.direction = [unitW(2);-unitW(1)];
                    u = (combinedRadius*invTimeStep - wLength)*unitW;
                end
                % lineObj.point = va + 0.5*u;  % Equation 5 (definition of ORCA_AB)
                lineObj.point  = -va + 0.5*u;  % Equation 5 (definition of ORCA_AB)
                agnt_orcalines = [agnt_orcalines;lineObj];
            end
        end

        % LINEAR PROGRAM ONE (CORRECT)
        function [obj,flag,result_LP1] = linearProgram1(obj,lineVector,lineNo,constraintRadius,optVelocity,directionOpt,currentResult)
            % Solves a one-dimensional lineObjar program on a specified lineObj
            % subject to linear constraints defined by lineObjs and a circular
            % constraint.
            % INPUTS:
            % obj          - The agent object
            % lineObjs        - A vector of lineObjs defining the lineObjar constraints.
            % lineObjNo       - The specified lineObj constraint
            % radius       - The radius of the circulat constraint.
            % optVelocity  - The optimization velocity.
            % directionOpt - True if the direction should be optimized.
            % result       - A reference to the result of the lineObjar program.
            % OUTPUT:
            % flag - True if successful.

            % THE PROBLEM IS THAT THE INDEX NUMBER 'lineNo' is different
            % between linearProgram1 and linearProgram2

            dotProduct   = dot(lineVector(lineNo).point,lineVector(lineNo).direction);
            discriminant = dotProduct^2 + constraintRadius^2 - dot(lineVector(lineNo).point,lineVector(lineNo).point);

            % Check if max speed circle fully invalidates lineObj constraint "lineObjNo"
            if discriminant < 0
                flag = logical(false);
                result_LP1 = currentResult;
                return
            end

            sqrtDiscriminant = sqrt(discriminant);
            tLeft  = -dotProduct - sqrtDiscriminant;
            tRight = -dotProduct + sqrtDiscriminant;

            % CHANGED TO PREVENT LOOP FROM ENTERING WHEN THE FAILURE OCCURS
            % ON THE FIRST ORCA-LINE
            i = 1;
            while (i < lineNo)
                denominator = obj.det(lineVector(lineNo).direction,lineVector(i).direction);
                numerator   = obj.det(lineVector(i).direction,(lineVector(lineNo).point - lineVector(i).point));

                if sqrt(dot(denominator,denominator)) <= obj.RVO_EPSILON    % The norm is less than the RVO constant
                    % lineObjs "lineObjNo" and "i" are (almost) parallel.
                    if numerator < 0
                        flag = logical(false);
                        result_LP1 = currentResult;
                        return
                    else
                        i = i + 1;                                          % Increment the counter before moving to next loop
                        continue
                    end
                end

                t = numerator/denominator;
                if denominator >= 0
                    tRight = min([tRight;t]);                               % lineObj "i" bounds lineObj "lineObjNo" on the right
                else
                    tLeft = max([tLeft,t]);                                 % lineObj "i" bounds lineObj "lineObjNo" on the left
                end

                if tLeft > tRight
                    flag = logical(false);
                    result_LP1 = currentResult;
                    return
                end
                i = i + 1;                                                  % Increment the counter before moving to next loop
            end

            % OPTIMISE DIRECTION & VELOCITY
            if directionOpt
                % Optimize direction
                if dot(optVelocity,lineVector(lineNo).direction) > 0
                    % Take right extreme
                    currentResult = lineVector(lineNo).point + tRight*lineVector(lineNo).direction;
                else
                    % Take left extreme
                    currentResult = lineVector(lineNo).point + tLeft*lineVector(lineNo).direction;
                end
            else
                % Optimise closest point
                t = dot(lineVector(lineNo).direction,(optVelocity - lineVector(lineNo).point));
                if t < tLeft
                    currentResult = lineVector(lineNo).point + tLeft*lineVector(lineNo).direction;
                elseif t > tRight
                    currentResult = lineVector(lineNo).point + tRight*lineVector(lineNo).direction;
                else
                    currentResult = lineVector(lineNo).point + t*lineVector(lineNo).direction;
                end
            end
            % DEFINE OUTPUT PARAMETERS
            result_LP1 = currentResult;
            flag = logical(true);

%             obj.LP1 = currentResult; % <<< REMOVE ME LATER
        end
        % LINEAR PROGRAM TWO
        function [obj,lineObjFlag,result_LP2]  = linearProgram2(obj,lineVector,constraintRadius,optVelocity,directionOpt)
            % Solves a two-dimensional lineObjar program subject to lineObjar
            % constraints defined by lineObjs and a circular constraint.
            % INPUTS:
            % obj          - The agent object
            % lineObjs        - lineObjs defining the lineObjar constraints.
            % radius       - The radius of the circular constraint.
            % optVelocity  - The optimization velocity.
            % directionOpt - True if the direction should be optimized.
            % result       - A reference to the result of the lineObjar program.
            % OUTPUT:
            % size_t       - The number of the lineObj it fails on, and the number of lineObjs if successful.

            if directionOpt
                % Optimise direction. Noe that the optimzation velocity is
                % of unit length in this case.
                currentResult = optVelocity*constraintRadius;
            elseif dot(optVelocity,optVelocity) > constraintRadius^2
                % Optimize closest point and ouside circle.
                currentResult = obj.normalise(optVelocity)*constraintRadius;
            else
                % Optimise closest point and inside circle
                currentResult = optVelocity;
            end

            % MOVE THROUGH THE lineObjS AND COMPARE THE lineObjAR CONSTRAINTS
            for i = 1:numel(lineVector)

                if obj.det(lineVector(i).direction,(lineVector(i).point - currentResult)) > 0
                    % "result" does not satisfy constraint "i". Compute new optimal "result".
                    tempResult = currentResult;

                    % EVALUATE lineObj USING
                    [obj,flag,currentResult] = obj.linearProgram1(lineVector,i,constraintRadius,optVelocity,directionOpt,currentResult);

                    if ~flag
                        % RETURN THE TEMP RESULT
                        result_LP2 = tempResult;
                        lineObjFlag = i;
                        return
                    end
                end

            end

            % RETURN THE FINAL RESULT OF THE LINEAR PROGRAM
            result_LP2 = currentResult;
            % RETURN THE TOTAL NUMBER OF LINES
            lineObjFlag = numel(lineVector); % The number of the line it fails on, and the number of lines if successful.
        end
        % LINEAR PROGRAM THREE
        function [obj,result_LP3] = linearProgram3(obj,lineVector,numObstLines,beginLine,constraintRadius,currentResult)
           % Solves a two-dimensional lineObjar program subject to lineObjar
           % constraints defined by lineObjs and a circular constraint.
           % INPUTS:
           % obj             - The agent object
           % lineObjs        - Linears defining the lineObjar constraints.
           % numObstlineObjs - Count of obstacle lineObjs.
           % beginlineObj    - The lineObj on which the 2-d linear program failed.
           % radius          - The radius of the circular constraint.
           % result          - A reference to the result of the linear program.

           % INITIAL DISTANCE TO "BEST CASE" POINT
           distance = 0;

           % IF THERE ARE NO OBSTACLES, ONLY AGENTS
           if numObstLines == 0
               numObstLines = 1; % Correct for c++ indexing
           end

           for i = beginLine:numel(lineVector)

               if obj.det(lineVector(i).direction,(lineVector(i).point - currentResult)) > distance
                   % Results does not satisfy constraint of lineObj 1
                   projlineObjs = [];
                   for j = numObstLines:i
                       % CREATE LINE OBJECT
                       newLine = struct('point',[],'direction',[]);

                       determinant = obj.det(lineVector(i).direction,lineVector(j).direction);
                       if sqrt(dot(determinant,determinant)) <= obj.RVO_EPSILON
                           % lineObj "i" and lineObj "j" are parallel.
                           if dot(lineVector(i).direction,lineVector(j).direction) > 0
                               % lineObj "i" and lineObj "j" point in the same
                               % direction
                               continue;
                           else
                               % lineObj "i" and lineObj "j" point in opposite
                               % directions
                               newLine.point = 0.5*(lineVector(i).point + lineVector(j).point);
                           end
                       else
                           newLine.point = lineVector(i).point + (obj.det(lineVector(j).direction,...
                                          (lineVector(i).point - lineVector(j).point)) / determinant) * lineVector(i).direction;
                       end

%                        newLine.direction = obj.normalise(lineVector(j).direction - lineVector(i).direction);
                       direction = lineVector(j).direction - lineVector(i).direction;
                       newLine.direction = direction/norm(direction);
                       projlineObjs = [projlineObjs;newLine];              % Append the new line to the projection lines
                   end
                   % HOLD ONTO THE CURRENT RESULT AS THE TEMP RESULT

                   tempResult = currentResult;

                   flowVector = [-lineVector(i).direction(2); lineVector(i).direction(1)]; % Vector in the best direction "going with the flow"

                   % RE-EVALUATE THE NEW LINE OBJECTS
                   [obj,lineObjFlag,currentResult] = obj.linearProgram2(projlineObjs,constraintRadius,flowVector,logical(true));
%                    obj.LP2 = currentResult;
                   if lineObjFlag < numel(projlineObjs)
                       % This should in principle not happen.  The result is by definition
                       % already in the feasible region of this linear program. If it fails,
                       % it is due to small floating point error, and the current result is
                       % kept.
                       currentResult = tempResult;
                   end
                   % The distance to the new optimal flow line
                   distance = obj.det(lineVector(i).direction,(lineVector(i).point - currentResult));
               end
           end
           % REDEFINE THE RESULT
           result_LP3 = currentResult;
%            obj.LP3 = result_LP3;
       end
    end
    %/////////////////////// "agent.cpp" FUNCTIONS ////////////////////////
    methods (Static)
        % AGENT UPDATE FUNCTION (CORRECT)
        function [position,velocity] = update(position,newVelocity)
            % INPUT HANDLING
            assert(numel(position) == numel(newVelocity),'Position and velocity updates must be of the same dimensions.');
            % UPDATE THE AGENT'S VELOCITY
            velocity = newVelocity;
            position = position + dt*velocity; % Update in accordance to the RVO2 library
        end
    end

    % ////////////////////// "vector_2.h" FUNCTIONS ///////////////////////
    methods (Static)
        % COMPUTES THE LENGTH OF A 2D VECTOR
        function [absLength] = abs(v)
           absLength = sqrt(dot(v,v));
        end
        % RETURN THE SQUARE LENGTH OF A 2D VECTOR
        function [absSq] = absSq(v)
           absSq = dot(v,v);
        end
        % COMPUTE THE DETERMININANT OF A 2D SQUARE MATRIX
        function [det_uv] = det(u,v)
            det_uv = u(1)*v(2) - u(2)*v(1);
        end
        % THE DOT PRODUCT
        function [dot_uv] = dot(u,v)
           dot_uv = u(1)*v(1) + u(2)*v(2);
        end
        % NORMALIZE A VECTOR
        function [unitVector] = normalise(v)
           unitVector = v/sqrt(dot(v,v));
        end
    end

    %/////////////////////// HELPER FUNCTIONS [JD] ////////////////////////
    methods (Static)
        % BUILD THE VERTEX-OBSTACLE LIST
        function [vertexSet] = defineRVOVertexObstacleSet(vertexData)
            % This function takes a list of points and contstucts a linked
            % obstacle list description. The vertices are assumed ordered
            % in terms of their connectivity (i.e. clockwise).

            % VERTEX STRUCTURE
            obstacleTemp = struct(...
                'point',zeros(2,1),...
                'nextObstacle',0,...
                'prevObstacle',0,...
                'unitDir',zeros(2,1),...
                'isConvex',boolean(true),...
                'id',0);

            numberOfVertices = size(vertexData,1);

            % LIST OF OBSTACLES
            vertexSet = repmat(obstacleTemp,numberOfVertices, 1 );
            % MOVE THROUGH THE VERTEX SET
            for i = 1:numberOfVertices
                % INDEX OF i+1
                ind_i_plus = mod(i,numberOfVertices) + 1;
                ind_i_plus2 = mod(i+1,numberOfVertices) + 1;
                % ALLOCATE POINT
                vertexSet(i).point = vertexData(i,:)';                     % Assign obstacle
                % SET NEXT AND PREVIOUS OBSTACLES/VERTEX POINTS
                vertexSet(i).nextObstacle = ind_i_plus;
                vertexSet(ind_i_plus).prevObstacle = i;
                % UNIT DIRECTION
                lineDirection = vertexData(ind_i_plus,:)' - vertexData(i,:)';
                vertexSet(i).unitDir = OMAS_geometry.unit(lineDirection);
                % IS CONVEX
                if numberOfVertices == 2
                    vertexSet(ind_i_plus).isConvex = 1;
                else
                    % CHECK IF NEXT POINT IS ON LEFT (IS CONVEX)
                    isLeft = OMAS_geometry.leftOf(vertexData(i,:),vertexData(ind_i_plus,:),vertexData(ind_i_plus2,:));
                    % IS POINT LEFT OF
                    vertexSet(ind_i_plus).isConvex = isLeft <= 0;          % If list is clockwise "<=0" , else list is not ">=0"
                end
                % SET THE INDEX OF POINT
                vertexSet(i).id = i;
            end
        end
    end
end

repo: OpenMAS


text:
function [obs,nobs,obss]=searchobsr(obss)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%search rover observations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global gls

if obss.idx>=obss.n
    obs=NaN;nobs=-1;return
end

obss.idx=obss.idx+1;
obs_tmp=obss.data(obss.idx,:);

time0=obss.data(:,1)+obss.data(:,2);
time1=obs_tmp(1)+obs_tmp(2);
idx=time0==time1;  pos=find(time0==time1);
if any(idx)
    obs0=obss.data(idx,:);  nobs=size(obs0,1); obss.idx=pos(end);
else
    nobs=0;obs=NaN;return;
end

obs=repmat(gls.obs_tmp,nobs,1);
for i=1:nobs
    obs(i).time.time=obs0(i,1);
    obs(i).time.sec=obs0(i,2);
    obs(i).sat=obs0(i,3);
    obs(i).P=obs0(i,4:6);
    obs(i).L=obs0(i,7:9);
    obs(i).D=obs0(i,10:12);
    obs(i).S=obs0(i,13:15);
    obs(i).LLI=obs0(i,16:18);
    obs(i).code=obs0(i,19:21);
end

return

repo: GINav


text:
function [KB] = TB2KB(TB)
% Convert computery things from terabytes to kilobytes.
% Chad A. Greene 2012
KB = TB*1073741824 ;


text:
function blkStruct = slblocks

blkStruct.Name = 'Real-Time Pacer';
blkStruct.OpenFcn = 'realtime_pacer_lib';
blkStruct.MaskInitialization = '';

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
%% AKAZE local features matching
%
% In this demo, we will learn how to use AKAZE local features to detect and
% match keypoints on two images. We will find keypoints on a pair of images
% with given homography matrix, match them and count the number of inliers
% (i.e. matches that fit in the given homography).
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/db/d70/tutorial_akaze_matching.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/features2D/AKAZE_match.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/xfeatures2D/LATCH_match.cpp>
%

%% Description
%
% You can find expanded version of this example
% <https://github.com/pablofdezalc/test_kaze_akaze_opencv here>.
%
% We are going to use images 1 and 3 from _Graffity_ sequence of Oxford
% dataset.
%
% <<https://docs.opencv.org/3.2.0/graf.png>>
%
% Homography is given by a 3-by-3 matrix shown below.
%

%% Code
% Options
inlier_threshold = 2.5;  % Distance threshold to identify inliers
nn_match_ratio = 0.8;    % Nearest neighbor

%%
% Load grayscale images
fnames = {
    fullfile(mexopencv.root(), 'test', 'graf1.png')
    fullfile(mexopencv.root(), 'test', 'graf3.png')
};
imgs = cell(size(fnames));
for i=1:numel(fnames)
    % downloading if necessary
    if exist(fnames{i}, 'file') ~= 2
        disp('Downloading image...')
        [~,name,ext] = fileparts(fnames{i});
        baseURL = 'https://cdn.rawgit.com/opencv/opencv/3.2.0/samples/data/';
        urlwrite([baseURL, name, ext], fnames{i});
    end
    % read image
    imgs{i} = cv.imread(fnames{i}, 'Grayscale',true);
end
[img1, img2] = deal(imgs{:});
whos img1 img2

%%
% ground-truth homography
H = [
    7.6285898e-01  -2.9922929e-01   2.2567123e+02
    3.3443473e-01   1.0143901e+00  -7.6999973e+01
    3.4663091e-04  -1.4364524e-05   1.0000000e+00
];
display(H)

%%
% Detect keypoints and compute descriptors using AKAZE
if true
    name = 'AKAZE';
    akaze = cv.AKAZE();
    [kpts1, desc1] = akaze.detectAndCompute(img1);
    [kpts2, desc2] = akaze.detectAndCompute(img2);
else
    % requires xfeatures2d from opencv_contrib
    name = 'LATCH';
    orb = cv.ORB('MaxFeatures',10000);
    kpts1 = orb.detect(img1);
    kpts2 = orb.detect(img2);
    latch = cv.LATCH();
    desc1 = latch.compute(img1, kpts1);
    desc2 = latch.compute(img2, kpts2);
end
whos kpts1 kpts2 desc1 desc2

%%
% Use brute-force matcher to find 2-nn matches.
% We use Hamming distance, because AKAZE uses binary descriptor by default.
matcher = cv.DescriptorMatcher('BruteForce-Hamming');
matches = matcher.knnMatch(desc1, desc2, 2);
whos matches

%%
% Use 2-nn matches to find correct keypoint matches.
% If the closest match is |ratio| closer than the second closest one, then the
% match is correct.
idx = cellfun(@(m) m(1).distance < nn_match_ratio * m(2).distance, matches);
matches = cellfun(@(m) m(1), matches(idx));
kp1 = kpts1([matches.queryIdx] + 1);
kp2 = kpts2([matches.trainIdx] + 1);
whos matches kp1 kp2

%%
% Check if our matches fit in the homography model.
% If the distance from first keypoint's projection to the second keypoint
% is less than threshold, then it it fits in the homography.
pts1 = cat(1, kp1.pt);
pts1(:,3) = 1;
pts1 = (H * pts1.').';
pts1 = bsxfun(@rdivide, pts1(:,1:2), pts1(:,3));
pts2 = cat(1, kp2.pt);
d = sqrt(sum((pts1 - pts2).^2, 2));
idx = d < inlier_threshold;
whos pts1 pts2
fprintf('%d inliers\n', nnz(idx));

%%
% We create a new set of matches for the inliers, because it is required
% by the drawing function.
kp1_good = kp1(idx);
kp2_good = kp2(idx);
matches_good = matches(idx);
for i=1:numel(matches_good)
    matches_good(i).queryIdx = i-1;
    matches_good(i).trainIdx = i-1;
end
whos kp1_good kp2_good matches_good

%%
% Show the result
res = cv.drawMatches(img1, kp1_good, img2, kp2_good, matches_good);
imshow(res), title(name)

%%
% Print some statistics
fprintf('%s Matching Results:\n', name);
fprintf('Keypoints 1: %d\n', numel(kpts1));
fprintf('Keypoints 2: %d\n', numel(kpts2));
fprintf('Matches: %d\n', numel(matches));
fprintf('Inliers: %d\n', numel(matches_good));
fprintf('Inliers Ratio: %f\n', numel(matches_good) / numel(matches));

%% References
%
% * *(AKAZE)* Pablo F Alcantarilla, Jesus Nuevo, and Adrien Bartoli. "Fast
%   explicit diffusion for accelerated features in nonlinear scale spaces".
%   Trans. Pattern Anal. Machine Intell, 34(7):1281-1298, 2011.
% * *(LATCH)* Gil Levi and Tal Hassner, "LATCH: Learned Arrangements of
%   Three Patch Codes", arXiv preprint arXiv:1501.03719, 15 Jan. 2015
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef CEC2010_F14 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- 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
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{14};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 1000;
            obj.upper    = zeros(1,obj.D) + 1000;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = 1 + PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(100*(Z(:,1:end-1).^2-Z(:,2:end)).^2+(Z(:,1:end-1)-1).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = sum((-Z.*cos(sqrt(abs(Z)))),2) - size(Z,2);
            PopCon(:,2) = sum((Z.*cos(sqrt(abs(Z)))),2) - size(Z,2);
            PopCon(:,3) = sum((Z.*sin(sqrt(abs(Z)))),2) - 10*size(Z,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function t = tsurf(F,V,varargin)
  % TSURF trisurf wrapper, easily plot triangle meshes with(out) face or vertex
  % indices. Attaches callbacks so that click-and-holding on the mesh and then
  % pressing 'm' launches meshplot (if available)
  %
  % t = tsurf(F,V);
  % t = tsurf(F,V,'ParameterName',ParameterValue, ...)
  %
  % Inputs:
  %   F  list of faces #F x 3
  %   V  vertex positiosn #V x 3 or #V x 2
  %   Optional:
  %     'VertexIndices' followed by
  %     'FaceIndices' followed by
  %                   0 -> off
  %                   1 -> text and grey background
  %                  -1 -> text and grey background (decremented by 1)
  %                  >1 -> text
  %     ... additional options passed on to trisurf
  % Outputs:
  %  t  handle to trisurf object
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % Example:
  %   % Compose with set function to set trisurf parameters to achieve
  %   % sharp isointervals
  %   tsurf(F,V, ...
  %     'FaceColor','interp', ...
  %     'FaceLighting','phong', ...
  %     'EdgeColor',[0.2 0.2 0.2]);
  %   colormap(flag(8))
  %
  % See also: trisurf
  %

  if issparse(V)
      V = full(V);
  end

  vertex_indices = 0;
  face_indices = 0;
  tets = [];
  buttondownfcn = 'default';

  v = 1;
  while v<=numel(varargin) && ischar(varargin{v})
    switch varargin{v}
    case 'VertexIndices'
      v = v+1;
      assert(v<=numel(varargin));
      vertex_indices = varargin{v};
    case 'FaceIndices'
      v = v+1;
      assert(v<=numel(varargin));
      face_indices = varargin{v};
    case 'Tets'
      v = v+1;
      assert(v<=numel(varargin));
      tets = varargin{v};
    case 'ButtonDownFcn'
      v = v+1;
      assert(v<=numel(varargin));
      buttondownfcn = varargin{v};
    otherwise
      break;
    end
    v = v+1;
  end

  % number of vertices
  n = size(V,1);

  % nuymber of dimensions
  dim = size(V,2);

  if(dim==2 || (dim ==3 && sum(abs(V(:,3))) == 0))
    V = [V(:,1) V(:,2) 0*V(:,1)];
    dim = 2;
  elseif ( (dim>3 || dim<2 ))
      if dim==4 && all(V(:,4)==1.0)
          V = V(:,1:3);
      else
    error('V must be #V x 3 or #V x 2');
    return;
      end
  end

  %tets = size(F,2) ==4 && (size(F,1)*4 > 1.01*size(boundary_faces(F),1));
  if isempty(tets)
    tets = false;
  if size(F,2) == 4
    VV = bsxfun(@minus,V,min(V))*max(max(V)-min(V));
    tets = sum(volume(VV,F))>1e-10;
    if ~tets
      Ftri = [F(:,[1 2 3]);F(:,[1 3 4])];
      Itri = repmat(1:size(F,1),1,2);
    end
  else
    tets = false;
    Ftri = F;
    Itri = 1:size(F,1);
  end
  end

  FI = 1:size(F,1);
  if face_indices<0
    FI = FI-1;
  end
  if tets
    t_copy = tetramesh(F,V,'FaceAlpha',0.5);
    FC = barycenter(V,F);
    if(abs(face_indices)==1)
      text(FC(:,1),FC(:,2),FC(:,3),num2str((FI)'),'BackgroundColor',[.7 .7 .7]);
    elseif(face_indices)
      text(FC(:,1),FC(:,2),FC(:,3),num2str((FI)'));
    end
    set(gcf,'Renderer','OpenGL');
  else
    t_copy = trisurf(F,V(:,1),V(:,2),V(:,3));
    FC = barycenter(V,F);
    if(abs(face_indices)==1)
      text(FC(:,1),FC(:,2),FC(:,3),num2str((FI)'),'BackgroundColor',[.7 .7 .7]);
    elseif(face_indices)
      text(FC(:,1),FC(:,2),FC(:,3),num2str((FI)'));
    end
  end

  % if 2d then set to view (x,y) plane
  if( dim == 2)
    view(2);
  end

  if vertex_indices
    visible = reshape(unique(F),[],1);
    VI = visible - (vertex_indices<0);
    if(abs(vertex_indices)==1)
      text(V(visible,1),V(visible,2),V(visible,3),num2str(VI),'BackgroundColor',[.8 .8 .8]);
    elseif(vertex_indices)
      text(V(visible,1),V(visible,2),V(visible,3),num2str(VI));
    end
  end
  % uncomment these to switch to a better 3d surface viewing mode
  %axis equal; axis vis3d;
  %axis(reshape([min(V(:,1:dim));max(V(:,1:dim))],1,dim*2));
  if ...
    strcmp(get(gca,'XLimMode'),'auto') && ...
    strcmp(get(gca,'YLimMode'),'auto') && ...
    strcmp(get(gca,'ZLimMode'),'auto')
    axis tight;
  end

  % surgery on varargin
  cdata_v = find(cellfun(@(c) ischar(c) && strcmpi(c,'cdata'),varargin),1,'last');
  if ~isempty(cdata_v) && (cdata_v+1)<=numel(varargin)
    varargin{cdata_v+1} = double(varargin{cdata_v+1});
  end

  if v<=numel(varargin)
    set(t_copy,varargin{v:end});
  end

  % Only output if called with output variable, prevents ans = ... output on
  % terminal
  if nargout>=1
    t = t_copy;
  end

  % subversively set up the callbacks so that if the user clicks and holds on
  % the mesh then hits m, meshplot will open with this mesh
  switch buttondownfcn
  case 'default'
    set(t_copy,'buttondownfcn',@ondown);
  case 'none'
  otherwise
  end

  function ondown(src,ev)
    if exist('point_mesh_squared_distance','file')==3
      [~,ci,C] = point_mesh_squared_distance(ev.IntersectionPoint,V,Ftri);
      warning off;
        B = barycentric_coordinates(C,V(Ftri(ci,1),:),V(Ftri(ci,2),:),V(Ftri(ci,3),:));
      warning on;
      on_vertex = sum(B<0.15)==2;
      color = [.7 .5 .5];
      if on_vertex
        ci = Ftri(ci,max(B)==B);
        color = [.5 .5 .7];
        C = V(ci,:);
      else
        ci = Itri(ci);
        C = FC(ci,:);
      end
      text_h = text(C(:,1),C(:,2),C(:,3),num2str(ci),'BackgroundColor',color);
    end
    set(gcf,'windowbuttonupfcn',@onup);
    set(gcf,'keypressfcn',@onkeypress);

    function onkeypress(src,ev)
      switch ev.Character
      case 'm'
        V3 = t_copy.Vertices;
        if size(V3,2) == 2
          V3(:,3) = 0*V(:,1);
        end
        fprintf('Opening mesh in viewmesh...\n');
        writeOBJ('.vm.obj',V3,F);
        system('/usr/local/bin/viewmesh .vm.obj');
      otherwise
        warning(['Unknown key: ' ev.Character]);
      end
    end

    function onup(src,ev)
      if exist('text_h','var') && ishandle(text_h)
        delete(text_h);
      end
      set(gcf,'windowbuttonupfcn','');
      set(gcf,'keypressfcn',      '');
      set(gcf,'windowbuttonmotionfcn','');
    end
  end

end


text:
%====================================================%
%===========ECEN 5793: DIP===========================%
%===========PROJECT 1================================%
%====================================================%
function A3 = imageresize(B3,p,m)
% Zoom or shrink a image
% Input:image B3: gray-image or RGB image
%       p: ratio, zooming (p>1), shrink(p<1)
%       m: mode; =0:(default) nearest neighbor interpolation
%                =1: Bilinear interpolation
%                =2: Bicubic interpolation
% Output: image A3
% @Trung Duong, trungd@okstate.edu
% Jan 25, 2009
%====================================================%
% Input checking and default values
error(nargchk(2, 3, nargin, 'struct'));
if nargin < 3, m = 0; end

dimB = length(size(B3));
if dimB== 2 % GRAY-IMAGE INPUT:
    A3 = gray_resize(B3,p,m);
elseif dimB== 3 % COLOR RGB-IMAGE INPUT:
    AR = gray_resize(B3(:,:,1),p,m);
    AG = gray_resize(B3(:,:,2),p,m);
    AB = gray_resize(B3(:,:,3),p,m);
    A3 = zeros([size(AR) 3]);
    A3(:,:,1) = AR;
    A3(:,:,2) = AG;
    A3(:,:,3) = AB;
else
    error('Improper input image');
end

%**************************************************%
%==================================================%
% Sub-Function: Resize a gray-image
% Same input argument with imageresize()
function A = gray_resize(B,p,m)
% Initialize new-grid and output
[N,M] = size(B);
xp = 1:1/p:N+1/p; yp = 1:1/p:M+1/p;
A = zeros(length(xp),length(yp));
% Symmetric Padding
npad = 3;
B = sym_pad(B,npad);
switch m
    case 0 % Nearest neighbor interpolation
        U = round(xp); V = round(yp);
        U(find(U<1)) = 1; V(find(V<1)) = 1;
        U(find(U>N)) = N; V(find(V>M)) = M;
        A = B(U+npad,V+npad);
    case 1 % Bilinear interpolation
        % Floor of (xp,yp)
        xf = floor(xp);     yf = floor(yp);
        % Distance to top-left neighbors
        [XF,YF] = ndgrid(xf,yf);
        [XP,YP] = ndgrid(xp,yp);
        u = XP - XF;         v = YP  - YF;
        % Change xf, yf for new padding image
        xf = xf + npad; yf = yf + npad;
        % Interpolation
        A = (1-u).*(1-v).*B(xf,yf)   + ...
                 (1-u).*v    .*B(xf,yf+1) + ...
                 u   .*(1-v) .*B(xf+1,yf) + ...
                 u   .*v     .*B(xf+1,yf+1);
   case 2
        % Floor of (xp,yp)
        xf = floor(xp);yf = floor(yp);
        % Distance to top-left neighbors
        [XF,YF] = ndgrid(xf,yf);
        [XP,YP] = ndgrid(xp,yp);
        u = XP - XF;        v = YP - YF;
        % Change xf, yf for new padding image
        xf = xf + npad; yf = yf + npad;
        % Interpolation: 16 neighbors
        for i = -1:2
            for j = -1:2
                if i==-1, sgi = -1; else sgi = 1; end
                if j==-1, sgj = -1; else sgj = 1; end
                A = A + mex_hat(sgi*(i-u)).*mex_hat(sgj*(j-v)).*B(xf+i,yf+j);
            end
        end
   otherwise
        error('Undefined Interpolation method');
end
%**************************************************%
%==================================================%
% Sub-Function: Mexican-hat kernel
function hx = mex_hat(x)
    hx = zeros(size(x));
    x = abs(x);
    ind1 = find(x<=1);  ind2 = find(x>1 & x<=2);

    hx(ind1) = 1 - 2*x(ind1).^2 + x(ind1).^3;
    hx(ind2) = 4 - 8*x(ind2) + 5*x(ind2).^2 - x(ind2).^3;
% END of sub-function
%==================================================%
% Sub-Function: symmetric padding, by default 2 pixels
% Input: gray image
function Bp = sym_pad(B,n)
    Bp = zeros(size(B)+2*n);
    Bp(n+1:end-n,n+1:end-n) = B;
    % Padding symmetrically 4 boundaries
    Bp(n:-1:1,n+1:end-n) = B(1:n,:);
    Bp(n+1:end-n,n:-1:1) = B(:,1:n);
    Bp(end-n+1:end,n+1:end-n) = B(end:-1:end-n+1,:);
    Bp(n+1:end-n,end-n+1:end) = B(:,end:-1:end-n+1);
% END of sub-function
%==================================================%


text:
function t = calcInitial( Deriv, X, xi)
%CALCINITIAL Calculate the initial guess for the line search.
%
%   Wrapper function for calcInitial_mex.c
%
%   See also calcGradient, calcFunction
%

%   GeomCG Tensor Completion. Copyright 2013 by
%   Michael Steinlechner
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

    t = -calcInitial_mex( Deriv.subs', Deriv.vals, ...
                          X.core.data, X.U{1}', X.U{2}', X.U{3}',...
                          xi.Y_tilde.data, xi.U1_tilde', xi.U2_tilde', xi.U3_tilde');

end

author: andrewssobral


text:
function x=fhtnat(data)
% The function implement the 1D natural(Hadamard) ordered fast Hadamard transform,
% wchich can be used in signal processing, pattern recongnition and Genetic alogorithms.
% This algorithm uses a Cooley-Tukey type signal flow graph and is implemented in N log2 N
% additions and subtractions. Data sequence length must be an integer power of 2.
%
% The inverse transform is the same as the forward transform except for the multiplication factor N.
%
% Example:
% x=[1 2 1 1];
% W=fhtnat(x);
%
% Author: Gylson Thomas
% e-mail: gylson_thomas@yahoo.com
% Asst. Professor, Electrical and Electronics Engineering Dept.
% MES College of Engineering Kuttippuram,
% Kerala, India, December 2006.
% copyright 2006.
% Reference: N.Ahmed, K.R. Rao, "Orthogonal Transformations for
% Digital Signal Processing" Spring Verlag, New York 1975. page-111.
N = length(data);
x=data;
L=log2(N);
k1=N; k2=1; k3=N/2;
for i1=1:L  %In-place iterations begins here
    L1=1;
    for i2=1:k2
        for i3=1:k3
            i=i3+L1-1; j=i+k3;
            temp1= x(i); temp2 = x(j);
            x(i) = temp1 + temp2;
            x(j) = temp1 - temp2;
        end
            L1=L1+k1;
    end
        k1 = k1/2;  k2 = k2*2;  k3 = k3/2;
end
x=inv(N)*x; %Delete this line for inverse transform

function x=fhtdya(data)
% The function implement the 1D dyadic (Paley) ordered fast Hadamard transform,
% wchich can be used in signal processing, pattern recongnition and Genetic alogorithms.
% This algorithm uses a Cooley-Tukey type signal flow graph and is implemented in N log2 N
% additions and subtractions. Data sequence length must be an integer power of 2.
%
% The inverse transform is the same as the forward transform except for the multiplication factor N.
%
% Example:
% x=[1 2 1 1];
% W=fhtdya(x);
%
% Author: Gylson Thomas
% e-mail: gylson_thomas@yahoo.com
% Asst. Professor, Electrical and Electronics Engineering Dept.
% MES College of Engineering Kuttippuram,
% Kerala, India, December 2006.
% copyright 2006.
% Reference: N.Ahmed, K.R. Rao, "Orthogonal Transformations for
% Digital Signal Processing" Spring Verlag, New York 1975. page-111.

N = length(data);
x=bitrevorder(data);
L=log2(N);
k1=N; k2=1; k3=N/2;
for i1=1:L  %In-place iteration begins here
    L1=1;
    for i2=1:k2
        for i3=1:k3
            i=i3+L1-1; j=i+k3;
            temp1= x(i); temp2 = x(j);
            x(i) = temp1 + temp2;
            x(j) = temp1 - temp2;
        end
            L1=L1+k1;
    end
        k1 = k1/2;  k2 = k2*2;  k3 = k3/2;
end
x=inv(N)*x; %Delete this line for inverse transform


text:
function hBodySurf = rts_initOrbPlot(hFig, orbitDispAxes, bodyInfo)
%rts_initOrbPlot Summary of this function goes here
%   Detailed explanation goes here

    hold(orbitDispAxes,'on');
    set(orbitDispAxes,'XTickLabel',[]);
    set(orbitDispAxes,'YTickLabel',[]);
    set(orbitDispAxes,'ZTickLabel',[]);
    grid(orbitDispAxes,'on');
    view(orbitDispAxes,3);

    if(~isempty(bodyInfo))
        dRad = bodyInfo.radius;
        [X,Y,Z] = sphere(30);
        hold(orbitDispAxes,'on');
        hBodySurf = surf(orbitDispAxes, dRad*X,dRad*Y,dRad*Z);
        hold(orbitDispAxes,'on');
        colormap(orbitDispAxes,bodyInfo.bodycolor);
    else
        hBodySurf = [];
    end
    axis(orbitDispAxes, 'equal');
    hold(orbitDispAxes,'off');
end


text:
%% Generate Test Problems + Results for opti_Install_Test
% Assumes you are in base opti directory!
% All simple, quick problems

%% Linear Programming
clc
clear
no = 5;
lp_tprob = cell(no,1);
lp_sval = zeros(no,1);
i = 1;

%LP 1
f = -[-1, 2]';
A = [2, 1;-4, 4];
b = [5, 5]';
lp_tprob{i} = optiprob('f',f,'ineq',A,b);
lp_sval(i) = -3.75;
i = i + 1;

%LP 2
f = -[50, 100];
A = [10, 5;4, 10; 1, 1.5];
b = [2500, 2000, 450]';
lp_tprob{i} = optiprob('f',f,'ineq',A,b);
lp_sval(i) = -21875;
i = i + 1;

%LP 3
f = [2, 3, 7, 7];
A = [1, 1, -2, -5;-1, 2, 1, 4];
b = [2, -3]';
e = [1, 1];
lb = zeros(4,1);
ub = [30 100 20 1]';
lp_tprob{i} = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub);
lp_sval(i) = 4;
i = i + 1;

%LP 4
f = [1, 2, 3, 7, 8, 8];
A = [5, -3, 2, -3, -1, 2; -1, 0, 2, 1, 3, -3;1, 2, -1, 0, 5, -1];
b = [-5, -1, 3]';
e = [1, 1, 1];
lb = zeros(6,1);
ub = 10*ones(6,1);
lp_tprob{i} = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub);
lp_sval(i) = 3;
i = i + 1;

%LP 5
f = -[1 2 3]';
A = [-1,1,1; 1,-3,1];
b = [20,30]';
Aeq = [1 1 1];
beq = 40;
lb = [0 0 0]';
ub = [40 inf inf]';
lp_tprob{i} = optiprob('grad',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub);
lp_sval(i) = -97.5;
i = i + 1;

%LP Clean Up + Save
clear f A b e Aeq beq lb ub i no
save 'Utilities/Install/Test Results/lp_test_results.mat'

%% Mixed Integer Linear Programming
clc
clear
no = 5;
milp_tprob = cell(no,1);
milp_sval = zeros(no,1);
i = 1;

%MILP 1
f = -[-1, 2]';
A = [2, 1;-4, 4];
b = [5, 5]';
xint = 'II';
milp_tprob{i} = optiprob('f',f,'ineq',A,b,'int',xint);
milp_sval(i) = -3;
i = i + 1;

%MILP 2
f = [2, 3, 7, 7];
A = [1, 1, -2, -5;-1, 2, 1, 4];
b = [2, 3]';
e = [1, 1];
lb = zeros(4,1);
ub = [30 100 20 1]';
xint = 'CICI';
milp_tprob{i} = optiprob('f',f,'mix',A,b,e,'bounds',lb,ub,'int',xint);
milp_sval(i) = 6;
i = i + 1;

%MILP 3
f = -[592, 381, 273, 55, 48, 37, 23];
A = [3534, 2356, 1767, 589, 528, 451, 304];
b = 119567;
lb = zeros(7,1);
ub = [100 50 33 20 77 44 20]';
xint = 'IIIIIII';
milp_tprob{i} = optiprob('f',f,'ineq',A,b,'bounds',lb,ub,'int',xint);
milp_sval(i) = -19979;
i = i + 1;

%MILP 4
f = -[6 5]';
A = [-3,5; 6,4; 3, -5; -6, -4];
b = [6;9;1;3];
xint = 'BC';
milp_tprob{i} = optiprob('grad',f,'ineq',A,b,'int',xint);
milp_sval(i) = -9.75;
i = i + 1;

%MILP 5
f = -[1 2 3 1]';
A = [-1 1 1 10; 1 -3 1 0];
b = [20;30];
Aeq = [0 1 0 -3.5];
beq = 0;
lb = [0 0 0 2]';
ub = [40 inf inf 3]';
xint = 'CCCI';
%Build & Options
milp_tprob{i} = optiprob('grad',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'int',xint);
milp_sval(i) = -122.5;

clear f A b e Aeq beq lb ub xint i no
save 'Utilities/Install/Test Results/milp_test_results.mat'

%% Quadratic Programming
clc
clear
no = 3;
qp_tprob = cell(no,1);
qp_sval = zeros(no,1);
i = 1;

%QP 1
H = eye(3);
f = -[2 3 1]';
A = [1 1 1;3 -2 -3; 1 -3 2];
b = [1;1;1];
qp_tprob{i} = optiprob('H',H,'f',f,'ineq',A,b);
qp_sval(i) = -2.83333333333333;
i = i + 1;

%QP 2
H = [1 -1; -1 2];
f = -[2 6]';
A = [1 1; -1 2; 2 1];
b = [2; 2; 3];
lb = [0;0];
qp_tprob{i} = optiprob('qp',H,f,'ineq',A,b,'lb',lb);
qp_sval(i) = -8.22222220552525;
i = i + 1;

%QP 3
H = [1 -1; -1 2];
f = -[2 6]';
A = [1 1; -1 2; 2 1];
b = [2; 2; 3];
Aeq = [1 1.5];
beq = 2;
lb = [0;0];
ub = [10;10];
qp_tprob{i} = optiprob('qp',H,f,'eq',Aeq,beq,'ineq',A,b,'bounds',lb,ub);
qp_sval(i) = -6.41379310344827;
i = i + 1;

clear H f A b e Aeq beq lb ub i no
save 'Utilities/Install/Test Results/qp_test_results.mat'

%% MI Quadratic Programming
clc
clear
no = 2;
miqp_tprob = cell(no,1);
miqp_sval = zeros(no,1);
i = 1;

%MIQP 1
H = [1 -1; -1 2];
f = -[2 6]';
A = [1 1; -1 2];
b = [3; 3.5];
miqp_tprob{i} = optiprob('hess',H,'grad',f,'ineq',A,b,'bounds',[0;0],[10;10],'int',1);
miqp_sval(i) = -11.5;
i = i + 1;

%MIQP 2
H = eye(3);
f = -[2 3 1]';
A = [1 1 1;3 -2 -3; 1 -3 2];
b = [1;1;1];
miqp_tprob{i} = optiprob('hess',H,'grad',f,'ineq',A,b,'int','CIC');
miqp_sval(i) = -2.75;
i = i + 1;

clear H f A b e Aeq beq lb ub i no
save 'Utilities/Install/Test Results/miqp_test_results.mat'

%% SDP Programming
clc
clear
no = 5;
sdp_tprob = cell(no,1);
sdp_sval = zeros(no,1);
i = 1;

%SDP1
f = 1;
A = eye(2);
C = -[0 sqrt(2); sqrt(2) 0];
sdp = sparse([C(:) A(:)]);
sdp_tprob{i} = optiprob('f',f,'sdcone',sdp);
sdp_sval(i) = sqrt(2);
i = i + 1;
sdp = [];

%SDP2
f = [1 1]';
C = -[0 2; 2 0];
A0 = [1 0; 0 0];
A1 = [0 0; 0 1];
sdp = sparse([C(:) A0(:) A1(:)]);
lb = [0;0];
ub = [10;10];
sdp_tprob{i} = optiprob('f',f,'sdcone',sdp,'bounds',lb,ub);
sdp_sval(i) = 4;
i = i + 1;
sdp = [];

%SDP3
f = [1 0 0 0]';
%SDP Constraint1 [x2 x3;x3 x4] <= x1*eye(2)
C = zeros(2);
A0 = eye(2);
A1 = -[1 0; 0 0];
A2 = -[0 1; 1 0];
A3 = -[0 0; 0 1];
sdp{1} = sparse([C(:) A0(:) A1(:) A2(:) A3(:)]);
%SDP Constraint2 [x2 x3; x3 x4] >= [1 0.2; 0.2 1]
C = [1 0.2; 0.2 1];
A0 = zeros(2);
A1 = [1 0; 0 0];
A2 = [0 1; 1 0];
A3 = [0 0; 0 1];
sdp{2} = sparse([C(:) A0(:) A1(:) A2(:) A3(:)]);
sdp_tprob{i} = optiprob('f',f,'sdcone',sdp);
sdp_sval(i) = 1.2;
i = i + 1;
sdp = [];

%SDP4
f = [1 1 1]';
C = -[0 1 2; 1 0 3; 2 3 100];
A0 = [1 0 0; 0 0 0; 0 0 0];
A1 = [0 0 0; 0 1 0; 0 0 0];
A2 = zeros(3);
sdp = sparse([C(:) A0(:) A1(:) A2(:)]);
lb = [10;0;0];
ub = [1000;1000;1000];
sdp_tprob{i} = optiprob('f',f,'sdcone',sdp,'bounds',lb,ub);
sdp_sval(i) = 10.1787148490962;
i = i + 1;
sdp = [];

%SDP5
f = -[1 2];
%SDP Constraint 1
C = [2 1; 1 2];
A0 = -[3 1; 1 3];
A1 = -zeros(2);
sdp{1} = sparse([C(:) A0(:) A1(:)]);
%SDP Constraint 2
C = [3 0 1; 0 2 0; 1 0 3];
A0 = -zeros(3);
A1 = -[3 0 1; 0 4 0; 1 0 5];
sdp{2} = sparse([C(:) A0(:) A1(:)]);
%SDP Constraint 3
C = zeros(2);
A0 = -[1 0; 0 0];
A1 = -[0 0; 0 1];
sdp{3} = sparse([C(:) A0(:) A1(:)]);
sdp_tprob{i} = optiprob('f',f,'sdcone',sdp);
sdp_sval(i) = 2.74999999779937;

clear f A0 A1 A2 A3 sdp lb ub i no
save 'Utilities/Install/Test Results/sdp_test_results.mat'

%% Nonlinear Least Squares
clc
clear
no = 5;
nls_tprob = cell(no,1);
nls_sval = zeros(no,1);
i = 1;

%NLS 1
fun = @(x) [10*(x(2)-x(1)^2); 1 - x(1)];
ydata = zeros(2,1);
x0 = [-1.2;1];
nls_tprob{i} = optiprob('fun',fun,'ydata',ydata,'x0',x0);
nls_sval(i) = 0;
i = i + 1;

%NLS 2
fun = @(x) [1.5 - x(1)*(1 - x(2)); 2.25 - x(1)*(1 - x(2)^2); 2.625 - x(1)*(1 - x(2)^3) ];
ydata = zeros(3,1);
x0 = [1;1];
nls_tprob{i} = optiprob('fun',fun,'ydata',ydata,'x0',x0);
nls_sval(i) = 0;
i = i + 1;

%NLS 3
fun = @(x) [1e4*x(1)*x(2) - 1;
            exp(-x(1)) + exp(-x(2)) - 1.0001];
ydata = zeros(2,1);
x0 = [0;1];
nls_tprob{i} = optiprob('fun',fun,'ydata',ydata,'x0',x0);
nls_sval(i) = 0;
i = i + 1;

%NLS 4
ii = (1:10)';
fun = @(x) 2 + 2*ii - (exp(0.2578*ii) + exp(0.2578*ii));
ydata = zeros(10,1);
x0 = [0.3;0.4];
nls_tprob{i} = optiprob('fun',fun,'ydata',ydata,'x0',x0);
nls_sval(i) = 124.36226865;
i = i + 1;

%NLS 5
fun = @(x) [10*(x(2)-x(1)^2);
            1 - x(1);
            sqrt(90)*(x(4)-x(3)^2);
            1 - x(3);
            sqrt(10)*(x(2) + x(4) - 2);
            10^(-0.5)*(x(2)-x(4))];
x0 = [3;-1;-3;-1];
ydata = zeros(6,1);
nls_tprob{i} = optiprob('fun',fun,'ydata',ydata,'x0',x0);
nls_sval(i) = 0;

clear fun ydata x0 ii i no
save 'Utilities/Install/Test Results/nls_test_results.mat'

%% Nonlinear Programming
clc
clear
no = 3;
nlp_tprob = cell(no,1);
nlp_sval = zeros(no,1);
i = 1;

%NLP 1
fun = @(x) 100*(x(2) - x(1)^2)^2 + (1 - x(1))^2;
grad = @(x)[[2*x(1)-400*x(1)*(-x(1)^2+x(2))-2];[-200*x(1)^2+200*x(2)]];
lb = [-inf; -1.5];
ub = [100; 100];
x0 = [-2; 1];
nlp_tprob{i} = optiprob('obj',fun,'grad',grad,'bounds',lb,ub,'x0',x0);
nlp_sval(i) = 0;
i = i + 1;

%NLP 2
fun = @(x) 1/3*(x(1) + 1)^3 + x(2);
grad = @(x)[[(x(1)+1)^2];[1]];
lb = [1;0];
ub = [100;100];
x0 = [1.125;0.125];
nlp_tprob{i} = optiprob('obj',fun,'grad',grad,'bounds',lb,ub,'x0',x0);
nlp_sval(i) = 8/3;
i = i + 1;

%NLP 3
fun = @(x) sin(x(1) + x(2)) + (x(1) - x(2))^2 - 1.5*x(1) + 2.5*x(2) + 1;
grad = @(x)[[2*x(1)-2*x(2)+cos(x(1)+x(2))-1.5];[2*x(2)-2*x(1)+cos(x(1)+x(2))+2.5]];
lb = [-1.5;-3];
ub = [4;3];
x0 = [0;0];
nlp_tprob{i} = optiprob('obj',fun,'grad',grad,'bounds',lb,ub,'x0',x0);
nlp_sval(i) = -0.5*sqrt(3)-pi/3;
i = i + 1;

clear fun grad lb ub x0 i no
save 'Utilities/Install/Test Results/nlp_test_results.mat'

%% MI Nonlinear Programming
clc
clear
no = 1;
minlp_tprob = cell(no,1);
minlp_sval = zeros(no,1);
i = 1;

%MINLP 1
obj = @(x) (x(1) - 5)^2 + x(2)^2 - 25;
nlcon = @(x) -x(1)^2 + x(2)-0.5;
nlrhs = 0;
nle = 1;
x0 = [4;0];
minlp_tprob{i} = optiprob('obj',obj,'ndec',2,'ivars',[1 2],'nlmix',nlcon,nlrhs,nle,'x0',x0);
minlp_sval(i) = -5;
i = i + 1;

clear obj A b nlcon nlrhs nle lb ub x0 i no
save 'Utilities/Install/Test Results/minlp_test_results.mat'

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function a = permute(a,order)
%PERMUTE      Permute array dimensions
%
%Implements the Matlab function "permute", similar syntax and semantics
%

% written  09/17/08     S.M. Rump
%

  if isa(a,'intval')
    if a.complex
      a.mid = permute(a.mid,order);
      if ~isequal(a.rad,0)
        a.rad = permute(a.rad,order);
      end
    else
      a.inf = permute(a.inf,order);
      a.sup = permute(a.sup,order);
    end
  else
    error('invalid call of permute for intervals')
  end