text:
function results = run_STRCF(seq, res_path, bSaveImage)

% Feature specific parameters
hog_params.cell_size = 4;
hog_params.compressed_dim = 10;
hog_params.nDim = 31;

grayscale_params.colorspace='gray';
grayscale_params.cell_size = 4;

cn_params.tablename = 'CNnorm';
cn_params.useForGray = false;
cn_params.cell_size = 4;
cn_params.nDim = 10;

% Which features to include
params.t_features = {
    struct('getFeature',@get_colorspace, 'fparams',grayscale_params),...
    struct('getFeature',@get_fhog,'fparams',hog_params),...
    struct('getFeature',@get_table_feature, 'fparams',cn_params),...
};

% Global feature parameters1s
params.t_global.cell_size = 4;                  % Feature cell size

% Image sample parameters
params.search_area_shape = 'square';    % The shape of the samples
params.search_area_scale = 5;         % The scaling of the target size to get the search area
params.min_image_sample_size = 150^2;   % Minimum area of image samples
params.max_image_sample_size = 200^2;   % Maximum area of image samples

% Spatial regularization window_parameters
params.feature_downsample_ratio = [4]; %  Feature downsample ratio
params.reg_window_max = 1e5;           % The maximum value of the regularization window
params.reg_window_min = 1e-3;           % the minimum value of the regularization window

% Detection parameters
params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection score
params.clamp_position = false;          % Clamp the target position to be inside the image

% Learning parameters
params.output_sigma_factor = 1/16;		% Label function sigma
params.temporal_regularization_factor = [15 15]; % The temporal regularization parameters

% ADMM parameters
params.max_iterations = [2 2];
params.init_penalty_factor = [1 1];
params.max_penalty_factor = [0.1, 0.1];
params.penalty_scale_step = [10, 10];

% Scale parameters for the translation model
% Only used if: params.use_scale_filter = false
params.number_of_scales = 5;            % Number of scales to run the detector
params.scale_step = 1.01;               % The scale factor

% Visualization
params.visualization = 1;               % Visualiza tracking and detection scores

% GPU
params.use_gpu = false;                 % Enable GPU or not
params.gpu_id = [];                     % Set the GPU id, or leave empty to use default

% Initialize
params.seq = seq;

% Run tracker
results = tracker(params);

repo: STRCF


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

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

% Last Modified by GUIDE v2.5 30-Nov-2017 13:55:03

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

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

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

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

% Update handles structure
guidata(hObject, handles);
global FilterType;
FilterType=0;
global std_rls;
global std_lms;
global std_fxlms;
global std_nlms;
global std_ap;
std_rls=0;
std_lms=0;
std_fxlms=0;
std_nlms=0;
std_ap=0;
global corr_rls;
global corr_lms;
global corr_fxlms;
global corr_nlms;
global corr_ap;
corr_rls=0;
corr_lms=0;
corr_fxlms=0;
corr_nlms=0;
corr_ap=0;

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

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

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

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global x;
global y;
global res;
global Fs;
global std_rls;
global std_lms;
global std_fxlms;
global std_nlms;
global std_ap;
global corr_rls;
global corr_lms;
global corr_fxlms;
global corr_nlms;
global corr_ap;
[x,Fs] = audioread('handel.wav');
[y,Fs]=audioread('handel_echo.wav');
length=get(handles.edit1,'String');
length=str2num(length);
ff=get(handles.edit2,'String');
ff=str2double(ff);
mu=get(handles.edit3,'String');
mu=str2double(mu);
global FilterType;
switch FilterType
    case 0
        ha=dsp.RLSFilter('Length',length,'ForgettingFactor',ff);
        [res,e] = ha(x,y);
        std_rls=std(e,1);
        corr_rls=corr(x,res);
    case 1
        ha=dsp.LMSFilter('Length',length,'StepSize',mu);
        [res,e] = ha(x,y);
        std_lms=std(e,1);
        corr_lms=corr(x,res);
    case 2
        ha=dsp.FilteredXLMSFilter('Length',length,'StepSize',mu);
        [res,e] = ha(x,y);
        std_fxlms=std(e,1);
        corr_fxlms=corr(x,res);
    case 3
        ha=dsp.LMSFilter('Length',length,'Method','Normalized LMS','StepSize',mu);
        [res,e] = ha(x,y);
        std_nlms=std(e,1);
        corr_nlms=corr(x,res);
    case 4
        ha=dsp.AffineProjectionFilter('Length',length,'StepSize',mu);
        [res,e] = ha(x,y);
        std_ap=std(e,1);
        corr_ap=corr(x,res);
end
audiowrite('handle Regenerated.wav',res,Fs);
set(handles.pushbutton4,'Enable','on');
set(handles.radiobutton1,'Enable','on');
set(handles.radiobutton1,'Value',0);
%axes(handles.axes1);
%spectrogram(x);
plot(handles.axes1,x);
title(handles.axes1,'Original Signal');
xlabel(handles.axes1,'Time Index'); ylabel(handles.axes1,'Signal Value');
plot(handles.axes2,y);
title(handles.axes2,'Echoed Signal');
xlabel(handles.axes2,'Time Index'); ylabel(handles.axes2,'Signal Value');
plot(handles.axes3,res);
title(handles.axes3,'Regenerated Signal');
xlabel(handles.axes3,'Time Index'); ylabel(handles.axes3,'Signal Value');
std_all=[std_rls,std_lms,std_nlms,std_fxlms,std_ap];
corr_all=[corr_rls,corr_lms,corr_nlms,corr_fxlms,corr_ap];
%bar(handles.axes4,mse_all);
bar(handles.axes4,std_all,0.4);
gap1=get(handles.axes4,'YTick');
ylength1=get(handles.axes4,'YLim');
ylim(handles.axes4,[0 ylength1(1,2)*1.1]);
if std_rls>0
    text(handles.axes4,1-0.3,std_rls+gap1(1,2)*0.3,'RLS');
end
if std_lms>0
    text(handles.axes4,2-0.3,std_lms+gap1(1,2)*0.3,'LMS');
end
if std_fxlms>0
    text(handles.axes4,4-0.3,std_fxlms+gap1(1,2)*0.3,'FX');
end
if std_nlms>0
    text(handles.axes4,3-0.4,std_nlms+gap1(1,2)*0.3,'NLMS');
end
if std_ap>0
    text(handles.axes4,5-0.3,std_ap+gap1(1,2)*0.3,'Affine');
end
title(handles.axes4,{'Standard deviation';'the smaller the better'});
xlabel(handles.axes4,'Filter type'); ylabel(handles.axes4,'Value');

bar(handles.axes5,corr_all,0.4);
gap2=get(handles.axes5,'YTick');
ylength2=get(handles.axes5,'YLim');
ylim(handles.axes5,[0 ylength2(1,2)*1.1]);
if corr_rls>0
    text(handles.axes5,1-0.3,corr_rls+gap2(1,2)*0.3,'RLS');
end
if corr_lms>0
    text(handles.axes5,2-0.3,corr_lms+gap2(1,2)*0.3,'LMS');
end
if corr_fxlms>0
    text(handles.axes5,4-0.3,corr_fxlms+gap2(1,2)*0.3,'FX');
end
if corr_nlms>0
    text(handles.axes5,3-0.5,corr_nlms+gap2(1,2)*0.3,'NLMS');
end
if corr_ap>0
    text(handles.axes5,5-0.3,corr_ap+gap2(1,2)*0.3,'Affine');
end
title(handles.axes5,{'Correlation coefficient';'the bigger the better'});
xlabel(handles.axes5,'Filter type'); ylabel(handles.axes5,'Value');

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global x;
global Fs;
sound(x,Fs);
pause(10);

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global y;
global Fs;
sound(y,Fs);
pause(10);

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global res;
global Fs;
sound(res,Fs);
pause(10);

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

% Hint: get(hObject,'Value') returns toggle state of radiobutton1
global x;
global y;
global res;
if get(handles.radiobutton1,'Value')
    axes(handles.axes1);
    spectrogram(x);
    title('Original Signal');
    axes(handles.axes2);
    spectrogram(y);
    title('Echoed Signal');
    axes(handles.axes3);
    spectrogram(res);
    title('Regenerated Signal');
else
    plot(handles.axes1,x);
    title(handles.axes1,'Original Signal');
    xlabel(handles.axes1,'Time Index'); ylabel(handles.axes1,'Signal Value');
    plot(handles.axes2,y);
    title(handles.axes2,'Echoed Signal');
    xlabel(handles.axes2,'Time Index'); ylabel(handles.axes2,'Signal Value');
    plot(handles.axes3,res);
    title(handles.axes3,'Regenerated Signal');
    xlabel(handles.axes3,'Time Index'); ylabel(handles.axes3,'Signal Value');
end
guidata(hObject,handles)

% --- Executes on selection change in popupmenu1.
function popupmenu1_Callback(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Determine the selected data set.
global FilterType;
set(handles.pushbutton4,'Enable','off');
set(handles.radiobutton1,'Enable','off');
str = get(hObject, 'String');
val = get(hObject,'Value');
% Set current data to the selected data set.
switch str{val}
    case 'RLS' % User selects peaks.
        FilterType=0;
        set(handles.edit1,'String','4');
        set(handles.edit2,'String','1');
        set(handles.edit2,'Enable','on');
        set(handles.edit3,'Enable','off');
    case 'LMS' % User selects membrane.
        FilterType=1;
        set(handles.edit1,'String','2');
        set(handles.edit3,'String','0.022');
        set(handles.edit2,'Enable','off');
        set(handles.edit3,'Enable','on');
    case 'FilteredXLMS' % User selects sinc.
        FilterType=2;
        set(handles.edit1,'String','2');
        set(handles.edit3,'String','0.002');
        set(handles.edit2,'Enable','off');
        set(handles.edit3,'Enable','on');
    case 'NLMS' % User selects membrane.
        FilterType=3;
        set(handles.edit1,'String','4');
        set(handles.edit3,'String','0.022');
        set(handles.edit2,'Enable','off');
        set(handles.edit3,'Enable','on');
    case 'AffineProjectionFilter' % User selects membrane.
        FilterType=4;
        set(handles.edit1,'String','2');
        set(handles.edit3,'String','0.022');
        set(handles.edit2,'Enable','off');
        set(handles.edit3,'Enable','on');
end
% Save the handles structure.
guidata(hObject,handles)

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

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

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

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

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double
value = str2double(get(handles.edit1,'String'));
flag = value > 0
if ~flag
   set(handles.edit1,'String','');
   warndlg('Value should be positive !! Please Enter Again');
% else
%    msgbox('Got the value, On it!');
end
set(handles.pushbutton4,'Enable','off');
set(handles.radiobutton1,'Enable','off');

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

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

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

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double
value = str2double(get(handles.edit2,'String'));
flag = value >= 0 && value <= 1
if ~flag
   set(handles.edit2,'String','');
   warndlg('Value should be between 0 and 1 !! Please Enter Again');
% else
%    msgbox('Got the value, On it!');
end
set(handles.pushbutton4,'Enable','off');
set(handles.radiobutton1,'Enable','off');

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

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

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

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double
value = str2double(get(handles.edit3,'String'));
flag = value > 0
if ~flag
   set(handles.edit3,'String','');
   warndlg('Value should be positive !! Please Enter Again');
% else
%    msgbox('Got the value, On it!');
end
set(handles.pushbutton4,'Enable','off');
set(handles.radiobutton1,'Enable','off');

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

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

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

% Hint: place code in OpeningFcn to populate axes1

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
plot(handles.axes1,5);% To override colorbar of spectrogram
cla(handles.axes1);
reset(handles.axes1);
plot(handles.axes2,5);
cla(handles.axes2);
reset(handles.axes2);
plot(handles.axes3,5);
cla(handles.axes3);
reset(handles.axes3);
cla(handles.axes4);
reset(handles.axes4);
cla(handles.axes5);
reset(handles.axes5);
set(handles.popupmenu1,'Value',1);
set(handles.edit2,'Enable','on');
set(handles.edit3,'Enable','off');
set(handles.pushbutton4,'Enable','off');
set(handles.radiobutton1,'Enable','off');
global FilterType;
FilterType=0;
global std_rls;
global std_lms;
global std_fxlms;
global std_nlms;
global std_ap;
std_rls=0;
std_lms=0;
std_fxlms=0;
std_nlms=0;
std_ap=0;
global corr_rls;
global corr_lms;
global corr_fxlms;
global corr_nlms;
global corr_ap;
corr_rls=0;
corr_lms=0;
corr_fxlms=0;
corr_nlms=0;
corr_ap=0;
set(handles.edit1,'String','4');
set(handles.edit2,'String','1');
set(handles.edit3,'String','0.022');

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

author: Spritea

repo: AEC

save_path: github-repos/MATLAB/Spritea-AEC


text:
function [D]=d2a2(A,B)
%D2A2 Pairwise squared L2 distance along 2nd axis.
% D = D2A2(A,B) Computes squared L2 (Euclidean) distance
% between all pairs of d-dimensional points in A and B.
%
% Inputs:
% A - m-by-d matrix of m d-dimensional points;
% B - n-by-d matrix of n d-dimensional points.
%
% Outputs:
% D - m-by-n matrix of pairwise distances.
%
% See also D2A1.
if nargin < 2
D = full(A*A');
d = diag(D);
D = bsxfun(@minus, d, 2*D);
D = bsxfun(@plus, d', D);
D = max(D, 0);
else
D = full(A*B');
D = bsxfun(@minus, sum(B.^2,2)', 2*D);
D = bsxfun(@plus, sum(A.^2,2), D);
D = max(D, 0);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% system process model in ecef frame with dcm formulation, Cse is the
% transformation between the sensor and e-frame, with the assumption that
% the sensor is also the body frame, we have Cse=Cbe.
%implement this function as sys_llh_phipsi_v000
function [STM Qd]=sys_ecef_dcm_v001(xyz_imu, Cse, acc, gyro, dt, imutype, modelNo)
% for modelNo=5 and imutype=5
% position   (1-3)
% velocity   (4-6)
% attitude   (7-9)
% acc bias drift
% gyro bias drift
% acc scale factor
% gyro scale factor
% acc turn on bias
% gyro turn on bias
% for modelNo=1 and imutype=5 the acc and gyro turn on bias are removed
WIE_E=7292115e-11;  %Earth rotation rate
g = 9.7803267714; % nominal g on earth surface
R=6317000; % earth average radius
%system disturbance coefs
N=zeros(9,6);
N(7:9,4:6)=-Cse; %attitude
N(4:6,1:3)=Cse; %velocity

%system matrix
A=zeros(9);
acc_e=Cse*acc;

%Position
A(1,4)=1;
A(2,5)=1;
A(3,6)=1;

%Velocity
unit=xyz_imu/norm(xyz_imu,2);
A(4:6,1:3)=-g*R^2/norm(xyz_imu,2)^3*(eye(3)-3*(unit*unit'))-skew([0;0;WIE_E])*skew([0;0;WIE_E]);
A(4:6,4:6)=-2*skew([0;0;WIE_E]);
A(4:6,7:9)=skew(acc_e);

%Attitude
A(7:9,7:9)=-skew([0;0;WIE_E]);
Anav=A;
Nnav=N;
% X(k+1) = ffun[X(k),U(k),V(k)]
% X(k+1) = Ak*X(k)+Gk*Vk

%%%%Imu error model parameters
[Aimu_d, Qimu_d, Cimu, Rimu]=imu_err_model_v001(acc, gyro, dt, imutype, modelNo);

%%%%Combine and discretize nav and imu models
% this discretization can also be accomplished by Loan's matrix exponential
% method, see sys_metric_phipsi_v000.m
Anav_d=eye(9)+dt*Anav;  %Use 1st order taylor series to discretize Anav
Qnav=Nnav*Rimu*Nnav';
Qnav_d=dt/2*(Anav_d*Qnav+Qnav*Anav_d');      %Use trapezoidal rule to discretize Rimu

STM=zeros(9+size(Aimu_d,1));
STM(1:9,1:9)=Anav_d;
STM(1:9,10:end)=Nnav*Cimu*dt;
STM(10:end,10:end)=Aimu_d;

Qd=zeros(9+size(Aimu_d,1));
Qd(1:9,1:9)=Qnav_d;
Qd(10:end,10:end)=Qimu_d;
Qd(1:9,10:end)=Nnav*Cimu*Qimu_d*dt/2;
Qd(10:end,1:9)=Qd(1:9,10:end)';
end

repo: ekfmonoslam


text:
function [upsilonaa, upsilonap] = lfmaaComputeUpsilonMatrix(gamma, sigma2, t1, t2, mode)

% LFMAACOMPUTEUPSILONMATRIX Upsilon matrix acce. accel. with t1, t2 limits
% FORMAT
% DESC computes a portion of the LFM kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

sigma = sqrt(sigma2);
gridt1 = repmat(t1, 1, length(t2));
gridt2 = repmat(t2', length(t1), 1);
timeGrid = gridt1 - gridt2;

if mode==0
    if nargout > 1
        upsilonap = lfmapComputeUpsilonMatrix(gamma, sigma2, t1, t2,1);
        upsilonaa = gamma^2*upsilonap + (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2).* ...
            ((gamma + 2*timeGrid/sigma2).*(2/sigma2 - (2*timeGrid/sigma2).^2) + 8*timeGrid/sigma^4);
    else
        upsilonaa = gamma^2*lfmapComputeUpsilonMatrix(gamma, sigma2, t1, t2,1) ...
            + (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2).* ...
            ((gamma + 2*timeGrid/sigma2).*(2/sigma2 - (2*timeGrid/sigma2).^2) + 8*timeGrid/sigma^4);
    end
else
    if nargout > 1
        upsilonap = lfmapComputeUpsilonMatrix(gamma, sigma2, t1, t2, 0);
        upsilonaa = gamma^2*upsilonap + (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2).* ...
            ((gamma + 2*timeGrid/sigma2).*(2/sigma2 - (2*timeGrid/sigma2).^2) + 8*timeGrid/sigma^4) ...
            + (2*gamma^2/(sqrt(pi)*sigma))*exp(-gamma*t1)*((gamma - 2*t2/sigma2).*exp(-(t2.^2)/sigma2)).';

    else
        upsilonaa = gamma^2*lfmapComputeUpsilonMatrix(gamma, sigma2, t1, t2, 0) ...
            + (2/(sqrt(pi)*sigma))*exp(-(timeGrid.^2)./sigma2).* ...
            ((gamma + 2*timeGrid/sigma2).*(2/sigma2 - (2*timeGrid/sigma2).^2) + 8*timeGrid/sigma^4) ...
            + (2*gamma^2/(sqrt(pi)*sigma))*exp(-gamma*t1)*((gamma - 2*t2/sigma2).*exp(-(t2.^2)/sigma2)).';
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = removeMeshEars(varargin)
%REMOVEMESHEARS Remove vertices that are connected to only one face.
%
%   [V, F] = removeMeshEars(V, F)
%   [V, F] = removeMeshEars(MESH)
%   Remove vertices that are connected to only one face. This removes also
%   "pending" faces.
%   Note that if the mesh has boundary, this may remove some regular faces
%   located on the boundary.
%
%   Example
%   removeMeshEars
%
%   See also
%     meshes3d, ensureManifoldMesh
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-01-08, using Matlab 8.6.0.267246 (R2015b)
% Copyright 2019-2022 INRA - Cepia Software Platform

[vertices, faces] = parseMeshData(varargin{:});

nVertices = size(vertices, 1);

% for each vertex, determine the number of faces it belongs to
vertexDegree = zeros(nVertices, 1);
for iv = 1:nVertices
    vertexDegree(iv) = sum(sum(faces == iv, 2) > 0);
end

% remove vertices with degree 1
inds = find(vertexDegree == 1);
[vertices, faces] = removeMeshVertices(vertices, faces, inds);

%% Format output

varargout = formatMeshOutput(nargout, vertices, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

    delta = (1 / m * (X * theta - y)' * X)';

    theta = theta - alpha * delta;

    % ============================================================

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

author: benoitvallon

repo: coursera-machine-learning


text:
function dep = dep_add(cdep, dep, ntgt_input, njtsubs)

% augment cdep tsubs references, and add them to dependency list
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: dep_add.m 2512 2008-12-01 13:21:29Z volkmar $

rev = '$Rev: 2512 $'; %#ok

for k = 1:numel(cdep)
    cdep(k).tgt_input = [ntgt_input cdep(k).tgt_input];
    cdep(k).jtsubs = [njtsubs cdep(k).jtsubs];
end;
if isempty(dep)
    dep = cdep(:);
else
    dep = [dep(:); cdep(:)];
end;

author: spm


text:
function varargout = equationPoints(disc)
%EQUATIONPOINTS   Points at which collocation is enforced.
%   In CHEBCOLLOC2, functions are discretized at 2nd kind points but equations
%   are enforced at 1st kind points, to avoid duplication at boundaries.

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

pointsFun = @(n) chebtech1.chebpts(n);
[varargout{1:nargout}] = valsDiscretization.points(disc, pointsFun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ama] = loadama(filename);

% LOADAMA read an inverted A-matrix and associated geometry information
% from an ama file that was written by Tom Oostendorp's DIPOLI
%
% Use as
%   [ama] = loadama(filename)
%
% See also LOADTRI, LOADMAT

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: loadama.m 2885 2011-02-16 09:41:58Z roboos $

fid = fopen(filename, 'rb', 'ieee-le');

version = fread(fid, 1, 'int');
if version~=10
  error(sprintf('%s is either not an inverted A matrix, or one of an old version', filename));
end

mode = fread(fid, 1, 'int');
ngeo = fread(fid, 1, 'int');

totpnt = 0;
totdhk = 0;
nrow   = 0;

% read the boundaries
geo  = [];
for i=1:ngeo
  geo(i).name    = char(fread(fid, [1 80], 'uchar'));
  geo(i).npnt    = fread(fid, 1, 'int');
  geo(i).pnt     = fread(fid, [3 geo(i).npnt], 'float')';
  geo(i).ndhk    = fread(fid, 1, 'int');
  geo(i).dhk     = fread(fid, [3 geo(i).ndhk], 'int')' + 1;  % Matlab indexing starts at 1
  geo(i).sigmam  = fread(fid, 1, 'float');
  geo(i).sigmap  = fread(fid, 1, 'float');
  geo(i).geocon  = fread(fid, ngeo, 'int');
  geo(i).deflat  = fread(fid, ngeo, 'float');
  totpnt = totpnt + geo(i).npnt;
  totdhk = totdhk + geo(i).ndhk;
end

% read the electrodes
if mode~=1
  elec.name    = char(fread(fid, [1 80], 'uchar'));
  elec.npnt    = fread(fid, 1, 'int');
  for i=1:(elec.npnt+1)
    elec.el(i).dhk  = fread(fid, 1, 'int') + 1; % Matlab indexing starts at 1
    elec.el(i).la   = fread(fid, 1, 'float');
    elec.el(i).mu   = fread(fid, 1, 'float');
    elec.el(i).name = char(fread(fid, [1 10], 'char'));
    % the ELECTRODE c-structure is padded to word boundaries, i.e. to 4 bytes
    dum = fread(fid, 2, 'char');
  end
  elec.vertex  = fread(fid, 1, 'int');
  elec.surface = fread(fid, 1, 'int');
  nrow = nrow + elec.npnt;
else
  elec = [];
end

% read the gradiometers
if mode~=0
  error('gradiometers not yet implemented');
else
  grad = [];
end

% read the inverted A-matrix
bi = fread(fid, [totpnt nrow], 'float')';

% read the isolated source compartment information, if present
iso_sur    = fread(fid, 1, 'int') + 1;  % Matlab indexing starts at 1
inner_only = fread(fid, 1, 'int');
if iso_sur~=0
  iso_totpnt = geo(iso_sur).npnt;
  iso_b      = fread(fid, [iso_totpnt iso_totpnt], 'float')';
else
  iso_b = [];
end

fclose(fid);

% put all local variables into a structure, this is a bit unusual programming style
% the output structure is messy, but contains all relevant information
tmp = whos;
ama = [];
for i=1:length(tmp)
  if isempty(strmatch(tmp(i).name, {'tmp', 'fid', 'ans', 'handles'}))
    ama = setfield(ama, tmp(i).name, eval(tmp(i).name));
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function LMstats(D, HOMEIMAGES)

D = LMvalidobjects(D);
%objectnames = lower(LMobjectnames(D));
[D, objectnames] = LMcreateObjectIndexField(D);
objectnames = lower(objectnames);
Nimages = length(D);
Nobjects = length(objectnames);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Collect data:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
instancecounts = sparse(Nobjects, Nimages);
areacounts = sparse(Nobjects, Nimages);
labelingdate = [];
imagendx = [];
objndx = [];

p = 0;
for n = 1:Nimages
    Nimages - n
    if isfield(D(n).annotation, 'object')
        objectsinimage = {D(n).annotation.object.name};

        %a = imfinfo(strrep(fullfile(HOMEIMAGES, D(n).annotation.folder, D(n).annotation.filename), '\', '/'));
        %nrows = a.Width;
        %ncols = a.Height;

        ncols = 0; nrows = 0;
        [TF, ndx] = ismember(strtrim(lower(objectsinimage)), objectnames);
        for m = 1:length(ndx)
            [X,Y] = getLMpolygon(D(n).annotation.object(m).polygon);
            area = polyarea(double(X),double(Y)); % ignores intersections

            instancecounts(ndx(m),n) = instancecounts(ndx(m),n)+1;
            areacounts(ndx(m),n) = areacounts(ndx(m),n) + area;
            nrows = max(nrows, max(Y(:)));
            ncols = max(ncols, max(X(:)));
        end
        areacounts(:,n) = areacounts(:,n) / double(nrows*ncols);

        % Collect date
        if isfield(D(n).annotation.object(1), 'date')
            for m = 1:length(D(n).annotation.object)
                d = D(n).annotation.object(m).date;
                if ~isempty(d)
                    p = p + 1;
                    labelingdate{p} = d(1:11);
                    imagendx(p) = n;
                    objndx(p) = D(n).annotation.object(m).namendx;
                end
            end
        end
    end
end
Labelingdate = datenum(labelingdate);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OBJECTS
% 1) distribution object labels
frequency = full(sum(instancecounts,2));
[sortedFreq, ndx] = sort(frequency, 'descend');

% Show object counts sorted by frequency
nn = min(40, length(ndx));
figure
subplot(121)
loglog(sortedFreq)
xlabel('rank')
ylabel('counts')
title('Distribution of objects')
axis('square'); axis('tight')
subplot(122)
barh(frequency(ndx(1:nn)))
set(gca, 'YTick', 1:nn)
set(gca, 'YtickLabel', objectnames(ndx(1:nn)))
axis([0 max(frequency)+5 0 nn+1])
grid on
title(sprintf('Most frequent obj. Number polygons: %d, number classes: %d', sum(frequency), length(frequency)))
axis('tight')

% histogram of number of objects per image
valid = find(full(sum(areacounts,1)>.97));
imagesperobject = full(sum(instancecounts,1));
figure
hist(imagesperobject(valid), [1:200])
xlabel('number of objects/image')
ylabel('counts')
title('Histogram of number of objects per image')

% 2) evolution dataset/time
lastDate = max(Labelingdate);
year = [2005 2006 2007 2008 2009 2010]
months = [1:12]
t = 0; dt = [];
counts = [];
countUniqueImages = [];
countObjectClasses = [];
for y = year
    for m = months
        t = t+1;
        d = (datenum(y,m,31));
        if d<lastDate

            dt{t} = datestr(d);
            j = find(Labelingdate<d);
            counts(t) = length(j);
            countUniqueImages(t) = length(unique(imagendx(j)));
            countObjectClasses(t) = length(unique(objndx(j)));
        else
            break
        end
    end
end

if length(dt)>48
init = 5;
figure
subplot(311)
plot(counts(init:end))
set(gca,'Xtick',[1 12 24 36 48]-init+1)
set(gca,'XTickLabel', dt([1 12 24 36 48]));
subplot(312)
plot(countUniqueImages(init:end),'r')
set(gca,'Xtick',[1 12 24 36 48]-init+1)
set(gca,'XTickLabel', dt([1 12 24 36 48]));
subplot(313)
plot(countObjectClasses(init:end),'g')
set(gca,'Xtick',[1 12 24 36 48]-init+1)
set(gca,'XTickLabel', dt([1 12 24 36 48]));
end

% 3) perplexity labels / proportion new labels added every day

% SCENES
% 4) distribution of scenes: if a scene is defined as the 5 largest
% objects: draw distribution of scene types.
Nobj = 4;
valid = find(full(sum(areacounts>0.05,1)>0) & full(sum(areacounts,1)>.97));
valid = find(full(sum(areacounts,1)>.97));
% for each scene collect only the 5 largest objects
scenes = sparse(Nobjects, Nimages);
for i = valid
    s = full(areacounts(:,i));
    jj = find(s>0);
    [ss,j] = sort(s(jj),'descend');

    L = min(length(j),Nobj);
    j = jj(j(1:L));

    scenes(j,i) = 100-(0:L-1);
end
scenes = scenes(:,valid)';

colors = 'rgcb'; k = 0;
figure
for ngram = [1 2 4 8]
    k = k+1;
    S = (scenes>100-ngram);
    % remove scenes with less than ngram objects
    c = sum(S,2);
    svalid = find(c>=ngram);

    [scenesunique, ns, j] = unique(S(svalid,:), 'rows');
    ns = valid(svalid(ns)); % recover initial indexing
    scenesunique = full(scenesunique');
    ju = unique(j);
    counts = hist(j, ju);
    [counts, ndx] = sort(counts, 'descend');

    loglog(counts, colors(k))
    hold on
    title('Distribution of scenes: a scene is a unique N-gram')
    axis('square'); axis('tight')
    drawnow
end
axis([1 3000 .9 3000])

% Create figure with most common scene within each scene type: look at the
% most frequent distribution of the other objects. Take the image closer to
% the mean. Make a figure with 10 different scene types.
ngram = 4;
S = (scenes>100-ngram);
c = sum(S,2);
svalid = find(c>=ngram);
[scenesunique, ns, j] = unique(S(svalid,:), 'rows');
ns = valid(svalid(ns)); % recover initial indexing
scenesunique = full(scenesunique');
ju = unique(j);
counts = hist(j, ju);
[counts, ndx] = sort(counts, 'descend');

figure
k=0;
for n = 1:1:49
    k = k+1;
    % get list of large objects in image
    o = find(scenesunique(:,ndx(n)));
    list = sprintf('%s,',objectnames{o});
    list = sprintf('%s (%d)', list(1:end-1), counts(n));

    % get typical scene from set
    img = LMimread(D, ns(ndx(n)), HOMEIMAGES);
    thumb = LMsceneThumbnail(D(ns(ndx(n))).annotation, HOMEIMAGES, img);

    subplottight(7,7,k,.1)
    imshow(thumb)
    title(list)
    drawnow
end
% 5) perplexity scenes

% TYPICALITY
% Show the most typical scene. Or run Kmeans using the distances and chose the 10 central images. The one with the lowest perplexity.

author: CSAILVision

repo: LabelMeToolbox


text:
% DEMHIGHFIVE1 Demonstration of hierarchical GP-LVM on walking and running data.
%
%	Description:
%

%	Copyright (c) 2007 Neil D. Lawrence
% 	demHighFive1.m version 1.1

if exist('diaryFile')
    diary(diaryFile)
end

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

hsvargplvm_init;

dataSetName = 'highFive';
capName = dataSetName;
capName(1) = upper(capName(1));
experimentNo = 1;
dirSep = filesep;
baseDir = datasetsDirectory;

%--- Load data
try
    load([baseDir 'dem' dataSetName]);
catch
    [void, errid] = lasterr;
    if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        skelA = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '20' dirSep '20.asf']);
        [YA, skelA] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '20' dirSep '20_11.amc'], skelA);
        seqInd = [50:4:113 114:155 156:4:size(YA, 1)];
        YA = YA(seqInd, :);
        %    YA(:, [4:end]) = asind(sind(YA(:, [4:end])));
        skelB = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '21' dirSep '21.asf']);
        [YB, skelB] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '21' dirSep '21_11.amc'], skelB);
        YB = YB(seqInd, :);
        %    YB(:, [4:end]) = asind(sind(YB(:, [4:end])));
        save([baseDir 'dem' dataSetName], 'YA', 'YB', 'skelA', 'skelB', ...
            'seqInd');
    else
        error(lasterr);
    end
end

Yall{1} = YA;
Yall{2} = YB;
clear('YA','YB');

%-- Set up model
numberOfDatasets = length(Yall);
globalOpt.indPoints = min(globalOpt.indPoints, size(Yall{1},1));

%-- Load datasets
for i=1:numberOfDatasets
    Y = Yall{i};
    dims{i} = size(Y,2);
    N{i} = size(Y,1);
    indTr = globalOpt.indTr;
    if indTr == -1
        indTr = 1:N{i};
    end
    if ~exist('Yts')
        indTs = setdiff(1:size(Y,1), indTr);
        Yts{i} = Y(indTs,:);
    end
    Ytr{i} = Y(indTr,:);

    t{i} = linspace(0, 2*pi, size(Y, 1)+1)'; t{i} = t{i}(1:end-1, 1);
    timeStampsTraining{i} = t{i}(indTr,1); %timeStampsTest = t(indTs,1);
end

for i=2:numberOfDatasets
    if N{i} ~= N{i-1}
        error('The number of observations in each dataset must be the same!');
    end
end

clear('Y')

%%--- Optimise %%--------- TO FIX (from this point and below) (this is now taken from svargplvm)
options = svargplvmOptions(Ytr, globalOpt, labelsTrain);

if ~isempty(globalOpt.dynamicsConstrainType)
    for i=1:numberOfDatasets
        % Set up dynamcis (i.e. back-constraints) model
        optionsDyn{i}.type = 'vargpTime';
        optionsDyn{i}.inverseWidth=30;
        %   optionsDyn.vardistCovars = vardistCovarsMult;
        optionsDyn{i}.initX = globalOpt.initX;
        optionsDyn{i}.constrainType = globalOpt.dynamicsConstrainType;

        if exist('timeStampsTraining')
            optionsDyn{i}.t = timeStampsTraining;
        end
        if exist('labelsTrain') && ~isempty(labelsTrain)
            optionsDyn{i}.labels = labelsTrain;
        end
    end
else
    optionsDyn= [];
end

model = svargplvmModelCreate(Ytr, globalOpt, options, optionsDyn);
if exist('diaryFile')
    model.diaryFile = diaryFile;
end

model.globalOpt = globalOpt;
model.options = options;

%%%% TEMP
if exist('whiteVar')
    model.dynamics.kern.comp{2}.variance = whiteVar;
end
%%%%

% Force kernel computations
params = svargplvmExtractParam(model);
model = svargplvmExpandParam(model, params);

%%
%fprintf('# Median of vardist. covars: %d \n',median(median(model.vardist.covars)));
%fprintf('# Min of vardist. covars: %d \n',min(min(model.vardist.covars)));
%fprintf('# Max of vardist. covars: %d \n',max(max(model.vardist.covars)));

if displayIters
    model = svargplvmOptimiseModel(model);
else
    model = svargplvmOptimiseModelNoDisplay(model);
end

%--------

%%
colordef white
ax = hgplvmHierarchicalVisualise(model, visualiseNodes, [], [0.03 ...
    0.5 0.03 0.03])
tar = get(ax, 'cameratarget');
pos = get(ax, 'cameraposition');
newPos = tar + (rotationMatrix(0, -pi/8, 3*pi/2)*(pos - tar)')';
set(ax, 'cameraposition', newPos)
set(ax, 'xlim', [-20 25]);
set(ax, 'ylim', [-15 8])
set(ax, 'visible', 'off')

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


% dubinsParameters
%   - Find Dubin's parameters between two configurations
%
% input is:
%   start_node  - [wn_s, we_s, wd_s, chi_s, 0, 0]
%   end_node    - [wn_e, wn_e, wd_e, chi_e, 0, 0]
%   R           - minimum turn radius
%
% output is:
%   dubinspath  - a matlab structure with the following fields
%       dubinspath.ps   - the start position in re^3
%       dubinspath.chis - the start course angle
%       dubinspath.pe   - the end position in re^3
%       dubinspath.chie - the end course angle
%       dubinspath.R    - turn radius
%       dubinspath.L    - length of the Dubins path
%       dubinspath.cs   - center of the start circle
%       dubinspath.lams - direction of the start circle
%       dubinspath.ce   - center of the end circle
%       dubinspath.lame - direction of the end circle
%       dubinspath.w1   - vector in re^3 defining half plane H1
%       dubinspath.q1   - unit vector in re^3 along straight line path
%       dubinspath.w2   - vector in re^3 defining position of half plane H2
%       dubinspath.w3   - vector in re^3 defining position of half plane H3
%       dubinspath.q3   - unit vector defining direction of half plane H3
%

function dubinspath = dubinsParameters(start_node, end_node, R)

  ell = norm(start_node(1:2)-end_node(1:2));
  if ell<2*R,
      disp('The distance between nodes must be larger than 2R.');
      dubinspath = [];
  else

    ps   = start_node(1:3);
    chis = start_node(4);
    pe   = end_node(1:3);
    chie = end_node(4);

    crs = ps' + R*rotz(pi/2)*[cos(chis); sin(chis); 0];
    cls = ps' + R*rotz(-pi/2)*[cos(chis); sin(chis); 0];
    cre = pe' + R*rotz(pi/2)*[cos(chie); sin(chie); 0];
    cle = pe' + R*rotz(-pi/2)*[cos(chie); sin(chie); 0];

    % compute L1
    connecting_line = cre-crs;
    north_unit = [1; 0; 0];
%     theta = acos(dot(connecting_line,north_unit)/norm(connecting_line));
    theta = atan2(connecting_line(2), connecting_line(1));
    L1 = norm(crs-cre) + R*mod(2*pi + mod(theta-pi/2, 2*pi) - mod(chis-pi/2, 2*pi), 2*pi) ...
         + R*mod(2*pi + mod(chie-pi/2, 2*pi) - mod(theta-pi/2, 2*pi), 2*pi);
    % compute L2
    connecting_line = cle-crs;
    ell = norm(connecting_line);
%     theta = acos(dot(connecting_line,north_unit)/norm(connecting_line));
    theta = atan2(connecting_line(2), connecting_line(1));
    theta2 = theta - pi/2 + asin(2*R/ell);
    if isreal(theta2)==0,
      L2 = 9999;
    else
      L2 = sqrt(ell^2-4*R^2) + R*mod(2*pi + mod(theta2, 2*pi) - mod(chis-pi/2, 2*pi), 2*pi) ...
         + R*mod(2*pi + mod(theta2+pi, 2*pi) - mod(chie+pi/2, 2*pi), 2*pi);
    end
    % compute L3
    connecting_line = cre-cls;
    ell = norm(connecting_line);
%     theta = acos(dot(connecting_line,north_unit)/norm(connecting_line));
    theta = atan2(connecting_line(2), connecting_line(1));
    theta2 = acos(2*R/ell);
    if isreal(theta2)==0,
      L3 = 9999;
    else
      L3 = sqrt(ell^2-4*R^2) + R*mod(2*pi + mod(chis+pi/2, 2*pi) - mod(theta+theta2, 2*pi), 2*pi) ...
         + R*mod(2*pi + mod(chie-pi/2, 2*pi) - mod(theta+theta2-pi, 2*pi), 2*pi);
    end
    % compute L4
    connecting_line = cle-cls;
%     theta = acos(dot(connecting_line,north_unit)/norm(connecting_line));
    theta = atan2(connecting_line(2), connecting_line(1));
    L4 = norm(cls-cle) + R*mod(2*pi + mod(chis+pi/2, 2*pi) - mod(theta+pi/2, 2*pi), 2*pi) ...
         + R*mod(2*pi + mod(theta+pi/2, 2*pi) - mod(chie+pi/2, 2*pi), 2*pi);
    % L is the minimum distance
    [L,idx] = min([L1,L2,L3,L4]);
    e1 = [1; 0; 0];     % north unit vector
    switch(idx),
        case 1,
            cs = crs;
            lams = 1;
            ce = cre;
            lame = 1;
            q1 = (ce-cs)/norm(ce-cs);
            w1 = cs + R*rotz(-pi/2)*q1;
            w2 = ce + R*rotz(-pi/2)*q1;
        case 2,
            cs = crs;
            lams = 1;
            ce = cle;
            lame = -1;
            ell = norm(ce-cs);
            connecting_line = ce - cs;
            theta = atan2(connecting_line(2), connecting_line(1));
%             theta = acos(dot(ce-cs,e1)/ell);
            theta2 = theta - pi/2 + asin(2*R/ell);
            q1 = rotz(theta2+pi/2)*e1;
            w1 = cs + R*rotz(theta2)*e1;
            w2 = ce + R*rotz(theta2+pi)*e1;
        case 3,
            cs = cls;
            lams = -1;
            ce = cre;
            lame = 1;
            ell = norm(ce-cs);
            connecting_line = ce - cs;
%             theta = acos(dot(ce-cs,e1)/ell);
            theta = atan2(connecting_line(2), connecting_line(1));
            theta2 = acos(2*R/ell);
            q1 = rotz(theta+theta2-pi/2)*e1;
            w1 = cs + R*rotz(theta+theta2)*e1;
            w2 = ce + R*rotz(theta+theta2-pi)*e1;
         case 4,
            cs = cls;
            lams = -1;
            ce = cle;
            lame = -1;
            q1 = (ce-cs)/norm(ce-cs);
            w1 = cs + R*rotz(pi/2)*q1;
            w2 = ce + R*rotz(pi/2)*q1;
    end
    w3 = pe';
    q3 = rotz(chie)*e1;

    % assign path variables
    dubinspath.ps   = ps;
    dubinspath.chis = chis;
    dubinspath.pe   = pe;
    dubinspath.chie = chie;
    dubinspath.R    = R;
    dubinspath.L    = L;
    dubinspath.cs   = cs;
    dubinspath.lams = lams;
    dubinspath.ce   = ce;
    dubinspath.lame = lame;
    dubinspath.w1   = w1;
    dubinspath.q1   = q1;
    dubinspath.w2   = w2;
    dubinspath.w3   = w3;
    dubinspath.q3   = q3;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% rotz(theta)
%%   rotation matrix about the z axis.
function R = rotz(theta)
    R = [...
        cos(theta), -sin(theta), 0;...
        sin(theta), cos(theta), 0;...
        0, 0, 1;...
        ];
end


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = (1/m) * (X'*X);
[U, S, V] = svd(Sigma);

% =========================================================================

end

author: zlotus


text:
function hammersley_dataset ( )

%*****************************************************************************80
%
%% HAMMERSLEY_DATASET generates a Hammersley dataset and writes it to a file.
%
%  Discussion:
%
%    This program is meant to be used interactively.  It's also
%    possible to prepare a simple input file beforehand and use it
%    in batch mode.
%
%    The program requests input values from the user:
%
%    * NDIM, the spatial dimension,
%    * N, the number of points to generate,
%    * STEP, the index of the first subsequence element to be computed.
%    * SEED(1:NDIM), the Hammersley sequence index corresponding to STEP = 0.
%    * LEAP(1:NDIM), the successive jumps in the Hammersley sequence.
%    * BASE(1:NDIM), the bases (usually distinct primes or -N).
%
%    The program generates the data, writes it to the file
%
%      hammersley_NDIM_N.txt
%
%    where "NDIM" and "N" are the numeric values specified by the user,
%    and then asks the user for more input.   To indicate that no further
%    computations are desired, it is enough to input a nonsensical
%    value, such as -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2004
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Generate a Hammersley dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program is meant to be used interactively.\n' );
  fprintf ( 1, '  It is also possible to prepare a simple input \n' );
  fprintf ( 1, '  file beforehand and use it in batch mode.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * NDIM, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * STEP, the index of the first subsequence element.\n' );
  fprintf ( 1, '  * SEED(1:NDIM), the sequence element\n' );
  fprintf ( 1, '    corresponding to STEP = 0\n' );
  fprintf ( 1, '  * LEAP(1:NDIM), the successive jumps in the sequence.\n' );
  fprintf ( 1, '  * BASE(1:M), the bases, usually distinct primes\n' );
  fprintf ( 1, '    or -N (to generate values like J/N).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program generates the data, writes it to the file\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    hammersley_NDIM_N.txt\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where "NDIM" and "N" are the numeric values specified\n' );
  fprintf ( 1, '  by the user, and then asks the user for more input.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To indicate that no further computations are \n' );
  fprintf ( 1, '  desired, it is enough to input a nonsensical value, \n' );
  fprintf ( 1, '  such as -1.\n' );

  while ( 1 )

    fprintf ( 1, '  *\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '*  Ready to generate a new dataset:\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '  *\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NDIM is the spatial dimension.\n' );
    fprintf ( 1, '  (Try "2" if you have no preference.\n' );
    fprintf ( 1, '  Any value less than 1 terminates execution.\n' );
    ndim = [];
    ndim = input ( '  Enter NDIM:  ' );

    fprintf ( 1, '  User input NDIM = %d\n', ndim );

    if ( ~halham_ndim_check ( ndim ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
      fprintf ( 1, '  The input value of NDIM = %d\n', ndim );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  N is the number of points.\n' );
    fprintf ( 1, '  (Try "25" if you have no preference.\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    n = [];
    n = input ( '  Enter N:  ' );

    fprintf ( 1, '  User input N = %d\n', n );

    if ( ~halham_n_check ( n ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
      fprintf ( 1, '  The input value of N = %d\n', n );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  STEP is the index of the first subsequence element.\n' );
    fprintf ( 1, '  (Try "0" or "1" if you have no preference.\n' );
    fprintf ( 1, '  (Any value less than 0 terminates execution.)\n' );
    step = [];
    step = input ( '  Enter STEP:  ' );

    fprintf ( 1, '  User input STEP = %d\n', step );

    if ( ~halham_step_check ( step ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
      fprintf ( 1, '  The input value of STEP = %d\n', step );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SEED(1:NDIM) is the starting element index\n' );
    fprintf ( 1, '  for each coordinate.\n' );
    fprintf ( 1, '  (Try "[0,0,...,0]" if you have no preference.\n' );
    fprintf ( 1, '  (Any value less than 0 terminates execution.)\n' );
    seed = [];
    seed(1:ndim) = input ( '  Enter SEED(1:NDIM):  ' );

    fprintf ( 1, '\n' );
    i4vec_transpose_print ( ndim, seed, 'SEED:' );

    if ( ~halham_seed_check ( ndim, seed ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
      fprintf ( 1, '  The negative input value of at least one entry of\n' );
      fprintf ( 1, '  SEED is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  LEAP(1:NDIM) is the leaping multiplier\n' );
    fprintf ( 1, '  for each coordinate.\n' );
    fprintf ( 1, '  (Try "[1,1,...,1]" if you have no preference.\n' );
    fprintf ( 1, '  (Another choice is a prime bigger than all the bases.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    leap = [];
    leap(1:ndim) = input ( '  Enter LEAP(1:NDIM):  ' );

    fprintf ( 1, '\n' );
    i4vec_transpose_print ( ndim, leap, 'LEAP:' );

    if ( ~halham_leap_check ( ndim, leap ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
      fprintf ( 1, '  The input value of at least one entry of\n' );
      fprintf ( 1, '  LEAP is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  BASE(1:NDIM) is the base for each coordinate,\n' );
    fprintf ( 1, '  usually distinct primes or -N for values like J/N.\n' );
    fprintf ( 1, '  (Try "[-N,2,3,5,7,11,13,...]" if you have no preference.\n' );
    fprintf ( 1, '  (Any value of 0 or 1 terminates execution.)\n' );
    base = [];
    base(1:ndim) = input ( '  Enter BASE:  ' );

    fprintf ( 1, '\n' );
    i4vec_transpose_print ( ndim, base, 'BASE:' );

    if ( ~hammersley_base_check ( ndim, base ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HAMMERSLEY_DATASET\n' );
      fprintf ( 1, '  The input value of at least one entry of\n' );
      fprintf ( 1, '  BASE is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
    end

    r = i4_to_hammersley_sequence ( ndim, n, step, seed, leap, base );

    file_out_name = ...
      strcat ( 'hammersley_', num2str ( ndim ), '_', num2str ( n ), '.txt' );

    halham_write ( ndim, n, step, seed, leap, base, r, file_out_name );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  The data was written to the file "%s".\n', ...
      file_out_name );

  end

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% evaluate simulated data - ICA and NMF
function separ(sep0, offset0, path_data, path_res, prename, niter, savethis, hint, sep_how)
% sep_how: string i - ica, n - nmf -> method for separation
if ~exist('hint', 'var')
    hint = 0;
end

if ~exist('sep_how', 'var')
    sep_how = 'in';
end

fprintf('Separation of components... \n')
for rr = 1: length(offset0)
    fprintf('\n%g: ',rr)
    for ll=1 : length(sep0)
        fprintf('.')
        p.namedir = [prename num2str(100*sep0(ll)) 'offset_' num2str(offset0(rr))];
        cd ([path_data p.namedir])
        for mm=1:niter

            namefile = [p.namedir '-iter_' num2str(mm)];
            load ([namefile '.mat'])
            %         ims(psf);
            %         SaveImageFULL('psf', 'pf');

            if sum(sep_how == 'i')>0 %ICA
                [icasig{mm}, A{mm}, W{mm}] = fastica (dveccr, 'numOfIC', 2, 'g', 'tanh');
                icapixICA{mm} = reshape(A{mm},32, 32, 2);
            end
            if sum(sep_how == 'n')>0 %NMF
                ncomp = 2; %number of components to be separated
                if hint
% % %                     dvec_ind = double(squeeze(reshape(array2im(dpixc_ind), p.nx*p.ny, 1, 2))); % vectors of resized images
%                     [out, bg(mm), bg_im]=backgroundoffset(dpixc);
                    [out, bg(mm), bg_im]=backgroundoffset(dpixc, 'no', 5, 20, 8); %empirical values...
                    dvec_bg = bg(mm)*ones(1, p.nx*p.ny);
%                     dvec_bg = p.offset*ones(1, p.nx*p.ny); %changed for offset 10...

% % %                     blinkmatrand = rand(p.Nt, ncomp);
                    blinkmatrand = blinkmat';
                    winit = [blinkmatrand,ones(p.Nt,1)];             %random weights will be assigned to firts two and bg fixed

%                     hinit = [dvec_ind; dvec_bg];       %original 'true' points + background
                    hinit = [rand(ncomp, p.nx*p.ny); dvec_bg];
                    ncomp = ncomp+1; %background added
                    [w{mm},h{mm}]=nmf(double(dveccr'),ncomp+1,1,winit,hinit, [3], [3]);

                else
                    [w{mm},h{mm}]=nmf(double(dveccr'),ncomp,1);
                end
                icapixNMF{mm} = shiftdim(reshape(h{mm},ncomp,32,32), 1);
            end
            %             imstiled(icapixICA{mm});
            %             SaveImageFULL([p.namedir 'ICA_' num2str(mm)], 'p');
            %             imstiled(icapixNMF{mm});
            %             SaveImageFULL([p.namedir 'NMF_' num2str(mm)], 'p');
            close all

        end

        p.path_data = path_data;
        p.path = path_res;

        if savethis == 1
            fprintf('saving data \n');
            if ~(strcmp('p.path', 'p.path_data')) %not identical
                mkdir ([p.path p.namedir]);
                cd ([p.path p.namedir]);
            end
            save ([p.namedir '_separ'])
            writedata([],[],p,[p.namedir '_param'])
        end
    end
end

fprintf('\n')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%PREX_PLOTC  PRTools example on the dataset scatter and classifier plot
help prex_plotc
echo on
              % Generate Higleyman data
A = gendath([100 100]);
              % Split the data into the training and test sets
[C,D] = gendat(A,[20 20]);
              % Compute classifiers
w1 = ldc(C);        % linear
w2 = qdc(C);        % quadratic
w3 = parzenc(C);    % Parzen
w4 = dtc(C);        % decision tree
              % Compute and display errors
              % Store classifiers in a cell
W = {w1,w2,w3,w4};
              % Plot errors
disp(D*W*testc);
              % Plot the data and classifiers
figure
              % Make a scatter-plot
scatterd(A);
              % Plot classifiers
plotc({w1,w2,w3,w4});
echo off

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function res = matrix_max_percent(mat)
    res = floor(mat/max(mat(:))*100);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function D = compute_histogram_distance(H, options)

% compute_histogram_distance - compute distance between histograms
%
%   D = compute_histogram_distance(H, options);
%
%   H(:,i) is the ith histogram.
%	D(i,g) is the distance between histogram i and j.
%
%	options.histmetric is the metric used to compute the distance d(g,h) between two histograms.
%	We denote by G and H the cumulative distribution, i.e. G=cumsum(g)
%		'l2'  -> d(g,h)^2 = sum_i (g(i)-h(i))^2
%		'l1'  -> d(g,h)   = sum_i abs(g(i)-h(i))
%		'cl2' -> d(g,h)^2 = sum_i (G(i)-H(i))^2
%		'cl1' -> d(g,h)   = sum_i abs(G(i)-H(i))
%		'chi2' -> d(g,h) = sum_i (G(i)-H(i))^2 / (G(i)+H(i))
%		'bhatta' -> d(g,h) = 1 - sum_i sqrt(G(i)*H(i))
%   options.sigma is a pre-smoothing factor, counted in pixels.
%
%   Copyright (c) 2007 Gabriel Peyre

if isfield(options, 'sigma')
	sigma = options.sigma;
else
	sigma = 1.2;
end
if isfield(options, 'histmetric')
	histmetric = options.histmetric;
else
	histmetric = 'histmetric';
end

%% smooth the histograms
n = size(H,1);
m = size(H,2);
h = compute_gaussian_filter( 21,sigma/(2*n),n);
for i=1:m
    H(:,i) = perform_convolution(H(:,i),h);
end

% make them sum to 1
H = max(H,0);
H = H ./ repmat( sum(H,1), [n 1] );

%% compute distance
switch lower(histmetric)
    case {'cl2' 'cl1'}
        D = distance_matrix(cumsum(H), histmetric(2:end) );
    otherwise
        D = distance_matrix(H, histmetric );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = distance_matrix(X,metric)

n = size(X,1); % dimension
p = size(X,2);

A = repmat( reshape(X', [p 1 n] ), [1 p]);
B = permute(A, [2 1 3]);

switch metric
    case 'l1'
        D = sum( abs(A-B), 3 );
    case 'l2'
        D = sqrt( sum( (A-B).^2, 3 ) );
    case 'chi2'
        a = A+B; a(a<eps) = 1;
        D = sum( ((A-B).^2) ./ a, 3 );
    case 'bhatta'
        D = 1 - sum( sqrt(A.*B), 3 );
    otherwise
        error('Unknown method');
end

author: gpeyre

repo: matlab-toolboxes


text:
clear, clc;

% This is an example for running the function fusedLeastR
%
%% Problem
%
%  min  1/2 || A x - y||^2 + rho * ||x||_1 +
%       opts.fusedPenalty * sum_i |x_i-x_{i+1}|
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% [1]  Jun Liu, Lei Yuan, and Jieping Ye, An Efficient Algorithm for
%      a Class of Fused Lasso Problems, KDD, 2010.
%
%% ------------   History --------------------
%
% First version on 14 November 2009.
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                   % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/fusedLasso;

m=100;  n=1000;    % The data matrix is of size m x n

% for reproducibility
randNum=1;

% ---------------------- Generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);       % the data matrix

randn('state',(randNum-1)*3+2);
xOrin=randn(n,1);

randn('state',(randNum-1)*3+3);
noise=randn(m,1);
y=A*xOrin +...
    noise*0.01;     % the response

rho=0.001;          % the regularization parameter
                    % it is a ratio between (0,1), if .rFlag=1

%----------------------- Set optional items ------------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% termination criterion
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=500;   % maximum number of iterations

% normalization
opts.nFlag=0;       % without normalization

% regularization
opts.rFlag=1;       % the input parameter 'rho' is a ratio in (0, 1)
%opts.rsL2=0.01;    % the squared two norm term

% fused penalty
opts.fusedPenalty=0.01;

% line search
opts.lFlag=0;

%----------------------- Run the code LeastR -----------------------
tic;
[x1, funVal1, ValueL1]= fusedLeastR(A, y, rho, opts);
toc;

figure;
plot(funVal1,'-.b');
xlabel('Iteration (i)');
ylabel('The objective function value');

author: zzstefan

repo: BrainNetClass


text:
function [x,fval,exitflag,info] = opti_cbc(H,f,A,rl,ru,lb,ub,xtype,sos,x0,opts)
%OPTI_CBC Solve a MILP or MIQP using CBC
%
%   min f'*x      subject to:     rl <= A*x <= ru
%    x                            lb <= x <= ub
%                                 for i = 1..n: xi in Z
%                                 for j = 1..m: xj in {0,1}
%
%   x = opti_cbc([],f,A,rl,ru,lb,ub,xtype) solves a MILP where f is the
%   objective vector, A,rl,ru are the linear constraints, lb,ub are the
%   bounds and xtype is a string of integer variables ('C', 'I', 'B')
%
%   x = opti_cbc(H,f,A,rl,ru,lb,ub,xtype) solves a MIQP where H and f are
%   the objective matrix and vector respectively, A,rl,ru are the linear
%   constraints, lb,ub are the bounds and xtype is a string of integer
%   variables ('C', 'I', 'B') [NOT CURRENTLY SUPPORTED]
%
%   x = opti_cbc([H,...,xtype,sos) sos is a structure with fields type,
%   index, weight for SOS.
%
%   x = opti_cbc(H,...,sos,x0) uses x0 to warm start the solver.
%
%   x = opti_cbc(H,...,x0,opts) uses opts to pass cbcset options to the
%   solver.
%
%   [x,fval,exitflag,info] = opti_cbc(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR CBC USING THE MEX INTERFACE
%   See supplied Eclipse Public License

%   Copyright (C) 2012 Jonathan Currie (IPL)

t = tic;

% Handle missing arguments
if nargin < 11, opts = optiset; end
if nargin < 10, x0 = []; end
if nargin < 9, sos = []; end
if nargin < 8, xtype = repmat('C',size(f)); end
if nargin < 7, ub = []; end
if nargin < 6, lb = []; end
if nargin < 5, error('You must supply at least 5 arguments to opti_cbc'); end

warn = strcmpi(opts.warnings,'all');

%Add in cbcset settings
if(isfield(opts,'solverOpts') && ~isempty(opts.solverOpts))
    popts = cbcset(opts.solverOpts);
else
    popts = cbcset;
end
%Add in options from optiset
popts.maxnodes = opts.maxnodes;
popts.maxtime = opts.maxtime;
popts.primalTol = opts.tolrfun;
popts.dualTol = opts.tolrfun;
popts.intTol = opts.tolint;
%Add objective bias
if(isfield(opts,'objbias')), popts.objbias = opts.objbias; end
%Setup display level
popts.display = dispLevel(opts.display,1,1);
popts.optiver = optiver;

%Check sparsity
if(~issparse(A))
    if(warn), optiwarn('OPTI:NotSparseA','The A matrix should be sparse, correcting: [sparse(A)]'); end
    A = sparse(A);
end
if(~issparse(H))
    if(~isempty(H) && warn), optiwarn('OPTI:NotSparseH','The H matrix should be sparse, correcting: [sparse(H)]'); end
    H = sparse(H);
end
%Check Sym Tril
if(any(any(triu(H,1) ~= 0)))
    if(warn), optiwarn('OPTI:NotTrilH','The H matrix should be Symmetric TRIL, correcting: [tril(H)]'); end
    H = tril(H);
end
%Check SOS
if(~isempty(sos))
    if(~isstruct(sos) || ~isfield(sos,'type') || ~isfield(sos,'index') || ~isfield(sos,'weight'))
        error('SOS constraints must be a structure with fields ''type'', ''index'' and ''weight''');
    end
    if(length(sos.type) == 1)
        sos.index = {sos.index};
        sos.weight = {sos.weight};
    end
end

%MEX contains error checking
[x,fval,status,iter,cobj] = cbc(H, f, A, rl, ru, lb, ub, xtype, sos, x0, popts);

%Assign Outputs
info.Nodes = iter;
info.AbsGap = abs(cobj-fval);
info.RelGap = abs(cobj-fval)/(1e-1 + abs(fval));
info.Time = toc(t);
info.Algorithm = 'CBC: Branch and Cut using CLP';

switch(status)
    case 0
        info.Status = 'Integer Optimal';
        exitflag = 1;
    case 1
        info.Status = 'Linear Relaxation Infeasible';
        exitflag = -1;
    case 2
        info.Status = 'Gap Reached';
        exitflag = 1;
    case 3
        info.Status = 'Maximum Nodes Reached';
        exitflag = 0;
    case 4
        info.Status = 'Maximum Time Reached';
        exitflag = 0;
    case 5
        info.Status = 'User Exited';
        exitflag = -5;
    case 6
        info.Status = 'Number of Solutions Reached';
        exitflag = 0;
    case 7
        info.Status = 'Linear Relaxation Unbounded';
        exitflag = -2;
	case 8
        info.Status = 'Proven Infeasible';
        exitflag = -1;
    otherwise
        info.Status = 'Unknown Termination';
        exitflag = -3;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function saveHDF5(fileName, model, params)
  %hdf5write(fileName, 'params', params)
  hdf5write(fileName, 'dec_wordvec', model.W_emb_tgt);
  hdf5write(fileName, 'dec_vocab', params.tgtVocab, 'WriteMode', 'append');

  if params.isBi
    hdf5write(fileName, 'enc_wordvec', model.W_emb_src, 'WriteMode', 'append');
    hdf5write(fileName, 'enc_vocab', params.tgtVocab, 'WriteMode', 'append');
  end

  % ip weight
  hdf5write(fileName, 'ip_weight', model.W_soft, 'WriteMode', 'append');

  % recurrent params
  input_col_ranges = 1:params.lstmSize;
  hidden_col_ranges = params.lstmSize+1:2*params.lstmSize;
  for dd=1:params.numLayers
    prefix_tgt = ['dec_d', num2str(dd), '_lstm:'];

    % input
    row_ranges = 1:params.lstmSize;
    hdf5write(fileName, [prefix_tgt 'input_gate'], [model.W_tgt{dd}(row_ranges,hidden_col_ranges) model.W_tgt{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

    % forget
    row_ranges = params.lstmSize+1:2*params.lstmSize;
    hdf5write(fileName, [prefix_tgt 'forget_gate'], [model.W_tgt{dd}(row_ranges,hidden_col_ranges) model.W_tgt{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

    % output
    row_ranges = 2*params.lstmSize+1:3*params.lstmSize;
    hdf5write(fileName, [prefix_tgt 'output_gate'], [model.W_tgt{dd}(row_ranges,hidden_col_ranges) model.W_tgt{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

    % signal
    row_ranges = 3*params.lstmSize+1:4*params.lstmSize;
    hdf5write(fileName, [prefix_tgt 'input_value'], [model.W_tgt{dd}(row_ranges,hidden_col_ranges) model.W_tgt{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

    if params.isBi
      prefix_src = ['enc_d', num2str(dd), '_lstm:'];

      % input
      row_ranges = 1:params.lstmSize;
      hdf5write(fileName, [prefix_src 'input_gate'], [model.W_src{dd}(row_ranges,hidden_col_ranges) model.W_src{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

      % forget
      row_ranges = params.lstmSize+1:2*params.lstmSize;
      hdf5write(fileName, [prefix_src 'forget_gate'], [model.W_src{dd}(row_ranges,hidden_col_ranges) model.W_src{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

      % output
      row_ranges = 2*params.lstmSize+1:3*params.lstmSize;
      hdf5write(fileName, [prefix_src 'output_gate'], [model.W_src{dd}(row_ranges,hidden_col_ranges) model.W_src{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');

      % signal
      row_ranges = 3*params.lstmSize+1:4*params.lstmSize;
      hdf5write(fileName, [prefix_src 'input_value'], [model.W_src{dd}(row_ranges,hidden_col_ranges) model.W_src{dd}(row_ranges, input_col_ranges)], 'WriteMode', 'append');
    end
  end

end

author: lmthang


text:
function imResult = blendMode_Multiply(A, B, offsetW, offsetH)
%% Multiply blending mode: multiplies the numbers for each pixel of the top
%   layer with the corresponding pixel for the bottom layer. The result is a
%   darker picture.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_Multiply));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

C = A .* B;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end

author: Sable

repo: mcbench-benchmarks


text:
function output = MSE(predicted, actual)
 output = mean(sum(((predictions-actual).^2)'));
end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [dLon_deg,dLat_deg,dz_km,maxZ_km,minZ_km] = request_3dgrid_params(tit)
    % prompt for grid spacing.
    % uses catalog "primeCatalog" for default depth limits

    dx = 0.1;
    dy = 0.1 ;
    dz = 5.00 ;
    ZG=ZmapGlobal.Data;

    def= {dx, dy, dz, max(ZG.primeCatalog.Depth), min(ZG.primeCatalog.Depth)}; % as numbers
    defstr = cellfun(@num2str,def,'UniformOutput',false); % converted to strings

    if ~exist('tit','var') || isempty(tit)
        tit ='Three dimesional analysis params';
    end
    prompt={ 'Spacing in Longitude (dx in [deg])',...
        'Spacing in Latitude  (dy in [deg])',...
        'Spacing in Depth    (dz in [km ])',...
        'Depth Range: deep limit [km] ',...
        'Depth Range: shallow limit [km] ',...
        };

    ni2 = inputdlg(prompt,tit,1,defstr); %as strings
    ni2 = cellfun(@str2double,ni2); %converted to numbers

    dLon_deg = ni2(1);
    dLat_deg = ni2(2);
    dz_km = ni2(3);
    maxZ_km = ni2(4);
    minZ_km = ni2(5);
end

author: CelsoReyes


text:
function f = sqrt( f )
%SQRT   Square root.
%   SQRT(F) returns the square root of a positive SEPARABLEAPPROX F.

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

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

if ( isreal(f) )
    % Positive/negative test.
    bool = singleSignTest(f);  % Returns TRUE if there is no sign change.
    if ( ~bool )
        error('CHEBFUN:SEPARABLEAPPROX:sqrt:notSmooth', ...
            'Sign change detected. Unable to represent the result.');
    end
end

f = compose(f, @sqrt);

end


text:
function [au] = ly2au(ly)
% Convert length from light years to astronomical units.
% Chad A. Greene 2012
au = ly*63241.07708807;


text:
function errorbar_tick(h,w,xtype)
%ERRORBAR_TICK Adjust the width of errorbars
%   ERRORBAR_TICK(H) adjust the width of error bars with handle H.
%      Error bars width is given as a ratio of X axis length (1/80).
%   ERRORBAR_TICK(H,W) adjust the width of error bars with handle H.
%      The input W is given as a ratio of X axis length (1/W). The result
%      is independent of the x-axis units. A ratio between 20 and 80 is usually fine.
%   ERRORBAR_TICK(H,W,'UNITS') adjust the width of error bars with handle H.
%      The input W is given in the units of the current x-axis.
%
%   See also ERRORBAR
%

% Author: Arnaud Laurent
% Creation : Jan 29th 2009
% MATLAB version: R2007a
%
% Notes: This function was created from a post on the french forum :
% http://www.developpez.net/forums/f148/environnements-developpement/matlab/
% Author : Jerome Briot (Dut)
%   http://www.mathworks.com/matlabcentral/newsreader/author/94805
%   http://www.developpez.net/forums/u125006/dut/
% It was further modified by Arnaud Laurent and Jerome Briot.

% Check numbers of arguments
error(nargchk(1,3,nargin))

% Check for the use of V6 flag ( even if it is depreciated ;) )
flagtype = get(h,'type');

% Check number of arguments and provide missing values
if nargin==1
	w = 80;
end

if nargin<3
   xtype = 'ratio';
end

% Calculate width of error bars
if ~strcmpi(xtype,'units')
    dx = diff(get(gca,'XLim'));	% Retrieve x limits from current axis
    w = dx/w;                   % Errorbar width
end

% Plot error bars
if strcmpi(flagtype,'hggroup') % ERRORBAR(...)

    hh=get(h,'children');		% Retrieve info from errorbar plot
    x = get(hh(2),'xdata');		% Get xdata from errorbar plot

    x(4:9:end) = x(1:9:end)-w/2;	% Change xdata with respect to ratio
    x(7:9:end) = x(1:9:end)-w/2;
    x(5:9:end) = x(1:9:end)+w/2;
    x(8:9:end) = x(1:9:end)+w/2;

    set(hh(2),'xdata',x(:))	% Change error bars on the figure

else  % ERRORBAR('V6',...)

    x = get(h(1),'xdata');		% Get xdata from errorbar plot

    x(4:9:end) = x(1:9:end)-w/2;	% Change xdata with respect to the chosen ratio
    x(7:9:end) = x(1:9:end)-w/2;
    x(5:9:end) = x(1:9:end)+w/2;
    x(8:9:end) = x(1:9:end)+w/2;

    set(h(1),'xdata',x(:))	% Change error bars on the figure

end


text:
classdef HomogenizedVarComputer < handle

    properties (Access = public)
        Cref
        dCref
        C
        dC
        rho
        drho
        Pp
        dPp
    end

    methods (Access = public, Static)

        function obj = create(cParams)
            f = HomogenizedVarComputerFactory();
            obj = f.create(cParams);
        end

    end

    methods (Access = public, Abstract)
        computeCtensor(obj)
        computePtensor(obj)
        computeDensity(obj)
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function A = probBasedRandomGraph(N,P)
%PROBBASEDRANDOMGRAPH Generates a graph with random connections whose
%                     probability is based on pair-wise probabilities.
%
%INPUTS:
% N: The number of nodes in the graph.
% P: A N-by-N matrix containing probabilities of pair-wise connections. If
%    an undirected graph is desired, this should be an upper triangular
%    matrix. Otherwise, a directed graph is assumed.
%
%OUTPUT:
% A: The N-by-N adjacency matrix for the realized graph. This has 2 on the
%    diagonal for self-loops to count both ends of the edges.
%
%EXAMPLE 1: Generate a random digraph using randomly generated pairwise
%           probabilities.
% rng(42)
% P = rand(10);
% A = probBasedRandomGraph(10,P);
% g = digraph(A);
% plot(g,"MarkerSize",indegree(g))
%
%EXAMPLE 2: Generate a random undirected graph using randomly generated
%           pairwise probabilities.
% rng(42)
% P = triu(rand(10));
% A = probBasedRandomGraph(10,P);
% g = graph(A);
% plot(g,"MarkerSize",degree(g))
%
%August 2022 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if any(P>1 | P<0)
    error('Edge probabilities must be in [0,1]')
end

if istriu(P)
    directed = false;
elseif istril(P)
    P = P';
    directed = false;
else
    directed = true;
end

r = rand(N);
A = zeros(N);
if directed
    for i = 1:N
        for j = 1:N
            if r(i,j)<P(i,j)
                if i~=j
                    A(i,j) = 1;
                else
                    A(i,j) = 2;
                end
            end
        end
    end
else %undirected
    for i = 1:N
        for j = i:N
            if r(i,j)<P(i,j)
                A(i,j) = 1;
            end
        end
    end
    A = A+A';
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Zica A T mu] = myICA(Z,r,dispFlag)
%--------------------------------------------------------------------------
% Syntax:       Zica = myICA(Z,r);
%               Zica = myICA(Z,r,dispFlag);
%               [Zica A T mu] = myICA(Z,r);
%               [Zica A T mu] = myICA(Z,r,dispFlag);
%
% Inputs:       Z is an (d x n) matrix containing n samples of an
%               d-dimensional random vector
%
%               r is the desired number of independent components
%
%               [OPTIONAL] dispFlag = {true false} sets the stdout print
%               state. The default value is dispFlag = true
%
% Outputs:      Zica is an (r x n) matrix containing the r independent
%               components - scaled to variance 1 - of the input samples
%
%               A and T are the ICA transformation matrices such that
%               Zr = T \ pinv(A) * Zica + repmat(mu,1,n);
%               is the r-dimensional ICA approximation of Z
%
%               mu is the (d x 1) sample mean of Z
%
% Description:  This function performs independent component analysis (ICA)
%               on the input samples using the FastICA algorithm with
%               Gaussian negentropy
%
% Author:       Brian Moore
%               brimoor@umich.edu
%
% Date:         April 26, 2015
%--------------------------------------------------------------------------

% Knobs
eps = 1e-6;         % Convergence criteria
maxSamples = 1000;  % Max # data points in sample mean calculation
maxIters = 100;     % Maximum # iterations

% Parse display flag
dispFlag = (nargin < 3) || dispFlag;

% Center and whiten input data
[Zcw T mu] = myCenterAndWhiten(Z);

% Parse whitened data
[d n] = size(Zcw);
if (n > maxSamples)
    % Truncate data for sample mean calculations
    Zcwt = Zcw(:,randperm(n,maxSamples));
else
    % Full data
    Zcwt = Zcw;
end

% Random initial weights
normRows = @(X) bsxfun(@times,X,1 ./ sqrt(sum(X.^2,2)));
W = normRows(rand(r,d));

% FastICA w/ Gaussian negentropy
k = 0;
err = inf;
while (err > eps) && (k < maxIters)
    % Increment counter
    k = k + 1;

    % Update weights
    Wlast = W; % Save last weights
    Sk = permute(Wlast * Zcwt,[1 3 2]);
    G = Sk .* exp(-0.5 * Sk.^2);
    Gp = Sk .* G;
    W = mean(bsxfun(@times,G,permute(Zcwt,[3 1 2])),3) + bsxfun(@times,mean(Gp,3),Wlast);
    W = normRows(W);

    % Decorrelate weights
    [U,S,~] = svd(W,'econ');
    W = U * diag(1 ./ diag(S)) * U' * W;

    % Update error
    err = max(1 - dot(W,Wlast,2));

    % Display progress
    if dispFlag == true
        sprintf('Iteration %i: max(1 - <w%i,w%i>) = %.4g\n',k,k,k - 1,err);
    end
end

% Transformation matrix
A = W;

% Independent components
Zica = A * Zcw;


text:
function zernike_poly_coef_test ( )

%*****************************************************************************80
%
%% ZERNIKE_POLY_COEF_TEST tests ZERNIKE_POLY_COEF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ZERNIKE_POLY_COEF_TEST\n' );
  fprintf ( 1, '  ZERNIKE_POLY_COEF determines the Zernike \n' );
  fprintf ( 1, '  polynomial coefficients.\n' );

  for m = 0 : n

    c = zernike_poly_coef ( m, n );

    r8poly_print ( n, c, '  Zernike polynomial' );

  end

  return
end


text:
%% Housekeeping
% In this version, one can choose to run model 1 or model 2 by selecting
% the appropriate flag "ncp_shocks". If ncp_shocks==false, we have model 1
% otherwise we have model 2.
% Model 1: nonstationary and cointegrated
% - Neutral technology shocks
% - Investment-specific technology shocks
% Model 2: nonstationary and cointegrated
% - Neutral technology shocks
% - Investment-specific technology shocks
% - Preference shocks
% cointegrated
clc

%% RISE the model
ncp_shocks=true;
m=rise('sgusea12',...
    'rise_flags',struct('CountryNames',{{'H','F'}},'ncp_shocks',ncp_shocks),...
    'saveas',true);

%% Assign steady state file

m=set(m,'steady_state_file','sstate_model12b');

%% create parameters

[p,priors]=create_parameters12(ncp_shocks);

m=set(m,'parameters',p);

%% Bring in the data

data=create_data();

%% Estibrate the model
clc
ms=estimate(m,'data',data,'estim_priors',priors);
%%
[mfilt,loglik]=filter(m,'data',data);

author: jmaih


text:
function [G, G_ind, related, n_genes_KO, G_time] = buildGmatrix(model_name, model, separate_isoform, numWorkers, printLevel)
% Build the G matrix required for the calculation of genetic Minimal Cut
% Sets (gMCSs).
%
% USAGE:
%
%    [G, G_ind, related, n_genes_KO, G_time] = buildGmatrix(model_name, model_struct, separate_isoform, numWorkers, printLevel)
%
% INPUTS:
%    model_name:          Name of the metabolic model under study.
%    model_struct:        Metabolic model structure (COBRA Toolbox format).
%    separate_isoform:    Character used to discriminate different isoforms of a gene.
%
% OPTIONAL INPUTS:
%    numWorkers:        Maximum number of workers
%                       * 0 - maximum provided by the system (automatic)
%                       * 1 - sequential
%                       * 2+ - parallel
%    printLevel:        show the reactions created in models.
%                       * 0 - shows nothing
%                       * 1 - shows progress by reactions (default)
%                       * 2+ - shows everything (reaction and network generation)
%
% OUTPUTS:
%    G:             G matrix.
%    G_ind:         Gene knockouts associated with each row in the G matrix.
%    related:       Relationships among rows of the G matrix.
%    n_genes_KO:    Number of genes whose knockout is added by each row of the G matrix.
%    G_time:        Calculation times for each of the steps.
%
% EXAMPLE:
%
%    [G, G_ind, related, n_genes_KO, G_time] = buildGmatrix('Recon2.v04', modelR204, '.')
%
% .. Authors:
%       - Inigo Apaolaza, 16/11/2017, University of Navarra, TECNUN School of Engineering.
%       - Luis V. Valcarcel, 19/11/2017, University of Navarra, TECNUN School of Engineering.
%       - Francisco J. Planes, 20/11/2017, University of Navarra, TECNUN School of Engineering.
%       - Inigo Apaolaza, 10/04/2018, University of Navarra, TECNUN School of Engineering.

time_a = tic;
% Generate name for temporary folder
global CBTDIR
tmpFolderName = [CBTDIR filesep '.tmp'];
if ~exist(tmpFolderName,'dir')  % Create directories if needed
    mkdir(tmpFolderName)
end
if ~exist([tmpFolderName filesep 'rxn_level_gMCSs'],'dir')
    mkdir([tmpFolderName filesep 'rxn_level_gMCSs'])
end
if ~exist([tmpFolderName filesep 'rxn_level_gMCSs_by_rxn'],'dir')
    mkdir([tmpFolderName filesep 'rxn_level_gMCSs_by_rxn'])
end
if ~exist([tmpFolderName filesep 'rxn_level_models'],'dir')
    mkdir([tmpFolderName filesep 'rxn_level_models'])
end

% Analyze the GPR rules in order to set the strategy for the calculation of
% the knockouts for the G matrix
grRules = model.grRules;
rxnGeneMat = model.rxnGeneMat;
n_rxns = size(rxnGeneMat, 1);

rxns_0_genes = sum(rxnGeneMat, 2) == 0;
rxns_1_gene = sum(rxnGeneMat, 2) == 1;
rxns_or = cellfun(@strfind, grRules, repmat({'or'}, n_rxns, 1), 'UniformOutput', false);
rxns_or = ~cellfun(@isempty, rxns_or) & sum(rxnGeneMat, 2) > 1;
rxns_and = cellfun(@strfind, grRules, repmat({'and'}, n_rxns, 1), 'UniformOutput', false);
rxns_and = ~cellfun(@isempty, rxns_and) & sum(rxnGeneMat, 2) > 1;
rxns_only_or = rxns_or & ~rxns_and;
rxns_only_and = ~rxns_or & rxns_and;
rxns_or_and = rxns_or & rxns_and;

n_rxns_0_genes = sum(rxns_0_genes);
n_rxns_1_gene = sum(rxns_1_gene);
n_rxns_only_or = sum(rxns_only_or);
n_rxns_only_and = sum(rxns_only_and);
n_rxns_or_and = sum(rxns_or_and);
n_rxns_total = n_rxns_0_genes+n_rxns_1_gene+n_rxns_only_or+n_rxns_only_and+n_rxns_or_and;

summary_1 = {'n_rxns_0_genes', n_rxns_0_genes; 'n_rxns_1_gene', n_rxns_1_gene;
    'n_rxns_only_or', n_rxns_only_or; 'n_rxns_only_and', n_rxns_only_and;
    'n_rxns_or_and', n_rxns_or_and; 'n_rxns_total', n_rxns_total};
if printLevel >=1
    fprintf('\nG MATRIX - Summary\n')
    disp(summary_1)
end
ini_time = toc(time_a);

% Step 1 - Reactions with 1 gene
% clc
if printLevel >=1
    disp('G MATRIX - STEP 1');
end
time_b = tic;
act_rxnGeneMat = rxnGeneMat(rxns_1_gene, :);
not_delete_cols = sum(act_rxnGeneMat, 1) ~= 0;
act_rxnGeneMat = act_rxnGeneMat(:, not_delete_cols);
G_ind_1 = model.genes(not_delete_cols);
n_KO_1 = length(G_ind_1);
tmp_G_1 = act_rxnGeneMat';
G_1 = zeros(n_KO_1, n_rxns);
G_1(:, rxns_1_gene) = tmp_G_1;
G_time(1, 1) = toc(time_b);

% Step2 - Reactions with more than one gene and only OR rules
% clc
if printLevel >=1
    disp('G MATRIX - STEP 2');
end
time_c = tic;
act_rxnGeneMat = rxnGeneMat(rxns_only_or, :);
pos_rxns_only_or = find(rxns_only_or);
for i = 1:n_rxns_only_or
    G_ind_2{i, 1} = model.genes(logical(act_rxnGeneMat(i, :)))';
end
tmp_G_2 = eye(n_rxns_only_or);
G_2 = zeros(n_rxns_only_or, n_rxns);
G_2(:, rxns_only_or) = tmp_G_2;
G_time(2, 1) = toc(time_c);

% Reactions with more than one gene and only AND rules
% clc
if printLevel >=1
    disp('G MATRIX - STEP 3');
end
time_d = tic;
act_rxnGeneMat = rxnGeneMat(rxns_only_and, :);
not_delete_cols = sum(act_rxnGeneMat, 1) ~= 0;
act_rxnGeneMat = act_rxnGeneMat(:, not_delete_cols);
G_ind_3 = model.genes(not_delete_cols);
n_KO_3 = length(G_ind_3);
tmp_G_3 = act_rxnGeneMat';
G_3 = zeros(n_KO_3, n_rxns);
G_3(:, rxns_only_and) = tmp_G_3;
G_time(3, 1) = toc(time_d);

% Reactions with more than one gene and both OR and AND rules
time_e = tic;
search_filename = [tmpFolderName filesep 'rxn_level_models' filesep 'rxn_level_' model_name '_and_or.mat'];
if exist(search_filename,'file')
    load(search_filename);
else
    pos_rxns_or_and = find(rxns_or_and);
    [models_or_and, rxnNumGenes_or_and] = GPR2models(model, pos_rxns_or_and, separate_isoform, numWorkers, printLevel);
    save(search_filename, 'models_or_and', 'rxnNumGenes_or_and');
end

search_filename_2 = [tmpFolderName filesep 'rxn_level_gMCSs' filesep 'rxn_level_gMCSs_' model_name '.mat'];
if exist(search_filename_2,'file')
    load(search_filename_2);
else
    target_b = 1e-3;
    n_mcs = 100000000000;
    timelimit = 5*60;
    pos_rxns_or_and = find(rxns_or_and);
    if printLevel >=1
        disp('G MATRIX - STEP 4');
        showprogress(0);
    end
    for i = 1:n_rxns_or_and
        if printLevel >=1
%             disp([num2str(i),' of ', num2str(n_rxns_or_and)]);
            showprogress(i/n_rxns_or_and);
        end
        search_filename_3 = [tmpFolderName filesep 'rxn_level_gMCSs_by_rxn' filesep 'rxn_level_gMCSs_' model_name '_rxn' num2str(pos_rxns_or_and(i)) '.mat'];
        if exist(search_filename_3,'file')
            load(search_filename_3);
            mcs{i, 1} = act_mcs;
        else
            act_model = models_or_and{i};
            nbio = find(act_model.c);
            rxns = act_model.rxns;
            tmp = repmat({'DM_'}, length(rxns), 1);
            DM = cellfun(@strfind, rxns, tmp, 'UniformOutput', false);
            DM = ~cellfun(@isempty, DM);
            n_DM = sum(DM);
            DM = rxns(find(DM));
%             options.rxn_set = DM;
% %             options.timelimit = timelimit;
%             options.target_b = target_b;
%             options.printLevel = 0;
            max_len_mcs = length(DM);
            [act_mcs, act_mcs_time] = calculateMCS(act_model, n_mcs, max_len_mcs,...
                'rxn_set', DM,...
                'timelimit', timelimit,...
                'target_b', target_b,...
                'printLevel', 0);
            mcs{i, 1} = act_mcs;
            mcs_time{i, 1} = act_mcs_time;
            save(search_filename_3, 'act_mcs', 'act_mcs_time');
        end
    end
    save(search_filename_2);
end

k = 0;
for i = 1:n_rxns_or_and
    load([tmpFolderName filesep 'rxn_level_gMCSs_by_rxn' filesep 'rxn_level_gMCSs_' model_name '_rxn' num2str(pos_rxns_or_and(i)) '.mat']);
    n_act_mcs = length(act_mcs);

    for j = 1:n_act_mcs
        act_G_ind = act_mcs{j};
        if ~iscell(act_G_ind) && isnan(act_G_ind)
        else
            act_G_ind = cellfun(@strrep, act_G_ind, repmat({'DM_'}, length(act_G_ind), 1), repmat({''}, length(act_G_ind), 1), 'UniformOutput', false);
            act_G_ind = cellfun(@strtok, act_G_ind, repmat({separate_isoform}, length(act_G_ind), 1), 'UniformOutput', false);
            k = k+1;
            G_4(k, :) = zeros(1, n_rxns);
            G_4(k, pos_rxns_or_and(i)) = 1;
            G_ind_4{k, 1} = act_G_ind';
        end
    end
end
G_time(4, 1) = toc(time_e);

% Delete isoforms in order to work at the gene level
time_f = tic;
G_ind_1 = cellfun(@strtok, G_ind_1, repmat({separate_isoform}, length(G_ind_1), 1), 'UniformOutput', false);

n_KO_2 = length(G_ind_2);
for i = 1:n_KO_2
    act_G_ind_2 = G_ind_2{i};
    act_G_ind_2 = cellfun(@strtok, act_G_ind_2, repmat({separate_isoform}, 1, length(act_G_ind_2)), 'UniformOutput', false);
    G_ind_2{i} = unique(act_G_ind_2);
end
G_ind_3 = cellfun(@strtok, G_ind_3, repmat({separate_isoform}, length(G_ind_3), 1), 'UniformOutput', false);

% Delete repeats
if printLevel >=1
    disp('G MATRIX - Delete Repeats');
end
tmp_G = [];
try tmp_G = [tmp_G; G_1]; end
try tmp_G = [tmp_G; G_2]; end
try tmp_G = [tmp_G; G_3]; end
try tmp_G = [tmp_G; G_4]; end
tmp_G_ind = [];
try tmp_G_ind = [tmp_G_ind; G_ind_1]; end
try tmp_G_ind = [tmp_G_ind; G_ind_2]; end
try tmp_G_ind = [tmp_G_ind; G_ind_3]; end
try tmp_G_ind = [tmp_G_ind; G_ind_4]; end
n_tmp_G_ind = length(tmp_G_ind);

k = 0;
clear G G_ind
for i = 1:n_tmp_G_ind
    if i == 1
        k = k+1;
        G(k, :) = tmp_G(i, :);
        if ~iscell(tmp_G_ind{i})
            G_ind{k, 1} = tmp_G_ind(i);
        else
            G_ind{k, 1} = sort(tmp_G_ind{i});
        end
    else
        if ~iscell(tmp_G_ind{i})
            act_G_ind = tmp_G_ind(i);
        else
            act_G_ind = sort(tmp_G_ind{i});
        end

        pos_equal = cellfun(@isequal, G_ind, repmat({act_G_ind}, length(G_ind), 1));
        if sum(pos_equal) > 0
            G(pos_equal, :) = G(pos_equal, :) + tmp_G(i, :);
        else
            k = k+1;
            G(k, :) = tmp_G(i, :);
            G_ind{k, 1} = act_G_ind;
        end
    end
end

% Fill the G matrix with reactions which are knocked out by a given KO
% without being a gMCS
n_genes_KO = cellfun(@length, G_ind);
[n_genes_KO, ind] = sort(n_genes_KO, 'ascend');
G_ind = G_ind(ind);
G = G(ind, :);
n_G_ind = length(G_ind);
for i = 1:n_G_ind
    act_G_ind = G_ind{i};
    n_act_G_ind = length(act_G_ind);
    pos = find(n_genes_KO > n_act_G_ind);
    greater_G_ind = G_ind(pos);
    n_greater_G_ind = length(pos);
    for j = 1:n_greater_G_ind
        if sum(ismember(G_ind{pos(j)}, act_G_ind)) == n_act_G_ind
            G(pos(j), :) = G(pos(j), :) + G(i, :);
        end
    end
end
G = double(G>0);

% Check the interconnections between KOs
if printLevel >=1
    disp('G MATRIX - Check Relations');
end
k = 0;
n_genes_KO = cellfun(@length, G_ind);
for i = 2:n_G_ind
    act_G_ind = G_ind(i);
    n_act_G_ind = length(act_G_ind{:});
    pos = find(n_genes_KO < n_act_G_ind);
    if ~isempty(pos)
        pos = pos(end);
        for j = 1:pos
            act_G_ind_2 = G_ind(j);
            tmp = ismember(act_G_ind{:}, act_G_ind_2{:});
            n_act_G_ind_2 = length(act_G_ind_2{:});
            if sum(tmp) == n_act_G_ind_2
                k = k+1;
                related(k, 1) = i;
                related(k, 2) = j;
            end
        end
    end
end

if exist('related')
    un_related = unique(related(:, 1));
    n_un_related = length(un_related);
    for i = 1:n_un_related
        act_KO = un_related(i);
        ind = find(related(:, 1) == act_KO);
        act_related = related(ind, 2);
        all_genes = [G_ind{act_related}];
        un_all_genes = unique(all_genes);
        n_un_all_genes = length(un_all_genes);
        n_genes_KO(act_KO) = n_genes_KO(act_KO)-n_un_all_genes;
    end
else
    related = NaN;
end

final_filename = [pwd filesep 'G_' model_name '.mat'];
G_time(5, 1) = toc(time_f)+ini_time;
save(final_filename, 'G', 'G_ind', 'related', 'n_genes_KO', 'G_time');
if printLevel >=1
    disp('The G Matrix has been successfully calculated');
end
rmdir(tmpFolderName, 's');
end

author: opencobra

repo: cobratoolbox


text:
function [newVOI]=dicomrt_closevoi(VOI)
% dicomrt_closevoi(VOI)
%
% Close VOI contours.
%
% VOI contains VOI contours.
%
% See also: dicomrt_close3dvoi
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check case and set-up some parameters and variables
[voi_temp,type,label]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(voi_temp);

newVOI=VOI;
for k=1:size(VOI,1) % loop through the number of VOIs
    for j=1:size(VOI{k,2},1) % loop through the number of slices/VOI
        temp=VOI{k,2}{j};
        isequalx=isequal(temp(1,1,1),temp(end,1,1));
        isequaly=isequal(temp(1,2,1),temp(end,2,1));
        if or(~isequalx,~isequaly)
            temp(end+1,:)=temp(1,:);
            newVOI{k,2}{j}=temp;
        end
    end
end

[newVOI]=dicomrt_restorevarformat(voi_temp,newVOI);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR