text:
function i4vec_run_count_test ( )

%*****************************************************************************80
%
%% I4VEC_RUN_COUNT_TEST tests I4VEC_RUN_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 20;
  test_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_RUN_COUNT_TEST\n' );
  fprintf ( 1, '  I4VEC_RUN_COUNT counts runs in an I4VEC\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ' Run Count        Sequence\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for test = 1 : test_num

    [ a, seed ] = i4vec_uniform_ab ( n, 0, 1, seed );

    run_count = i4vec_run_count ( n, a );

    fprintf ( 1, '  %8d        ', run_count )
    for i = 1 : n
      fprintf ( 1, '%2d', a(i) );
    end
    fprintf ( 1, '\n' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [AngData] = fft_angle(Xcube,fft_Ang,Is_Windowed)

Nr=size(Xcube,1);   %%%length of Chirp
Ne=size(Xcube,2);   %%%length of receiver
Nd=size(Xcube,3);   %%%length of chirp loop

% win = taylorwin(Ne,5,-60);
% win = win/norm(win);
for i = 1:Nd
    for j = 1:Nr
        if Is_Windowed
            win_xcube = reshape(Xcube(j,:,i),Ne,1).*taylorwin(Ne);
        else
            win_xcube = reshape(Xcube(j,:,i),Ne,1).*1;
        end
        AngData(j,:,i) = fftshift(fft(win_xcube,fft_Ang));
    end
end
end

author: Xiangyu-Gao


text:
% LORENZ system
% Visualization for dependency on noise level

clear all, close all, clc

SystemModel = 'LORENZ';

figpath = '../FIGURES/'; mkdir(figpath)
datapath = ['../DATA/EX_',SystemModel,'_Dependencies/'];
addpath('../utils');

%% Paramaters
NARXtraining = 'trainbr'; %'trainlm',trainbr
InputSignalType = 'sphs'; %sphs,sine2
Ntrain_vec = 3000; %1000;
N_LENGTHS = length(Ntrain_vec);

eta_vec = [0.01 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5];
N_ETA = length(eta_vec);

Nmodels = N_ETA*N_LENGTHS;

Nr = 50;
Nt = 1000; % length(tv);
Nvar = 3;
%% Load all models
ResultsALL(1:Nmodels,3) = struct('err', zeros(Nr,1), 'errM', zeros(Nr,1), 'xA', zeros(Nt,Nvar), 'xB', zeros(Nt,Nvar,Nr),'Ttraining',zeros(Nr,1));

count = 0;
for iN = 1:N_LENGTHS
    for iNoise = 1:N_ETA
        count = count + 1;

        %         ModelName = 'DMDc';
        %         datapath1 = [datapath,'DMDc/'];
        %         filename = fullfile(datapath1,['EX_LOTKA_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS.mat']);
        %         load(filename);
        %         ResultsALL(count,1) = Results;
        ResultsALL(count,1) = struct('err', zeros(Nr,1), 'errM', zeros(Nr,1), 'xA', zeros(Nt,Nvar), 'xB', zeros(Nt,Nvar,Nr),'Ttraining',zeros(Nr,1));

        ModelName = 'SINDYc';
        datapath1 = [datapath,'SINDYc/'];
        filename = fullfile(datapath1,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS.mat']);
        load(filename);
        ResultsALL(count,2) = Results;

        ModelName = 'NARX';
        datapath1 = [datapath,'NARX/',NARXtraining,'/'];
        filename = fullfile(datapath1,['EX_',SystemModel,'_SI_',ModelName,'_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_Eta',sprintf('%03g',100*eta_vec(iNoise)),'_STATS.mat']);
        load(filename);
        ResultsALL(count,3) = Results;
    end
end

%%
clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
symbols = {'o', 'd', 's'};
xx = 1:0.1:11;%0:0.01:0.6;
yy = zeros(3,length(xx));
err = zeros(Nr,N_ETA,3);
for iM = 1:3
    for i = 1:N_ETA, err(1:Nr,i,iM) = ResultsALL(i,iM).err(1:Nr); end
    if iM == 2
        xm = zeros(1,N_ETA);
        for k = 1:N_ETA
            TF = isnan(err(1:Nr,k,iM));
            xm(k) = median(err(TF==0,k,iM),1);
        end
    else
        xm = median(err(1:Nr,:,iM),1);
    end
    yy(iM,:) = spline(1:11,xm,xx); %eta_vec(:)
end
figure; hold on
data = err(1:Nr,eta_vec==0.4,2);% sine 0.4, sphs 0.45: 278.4062
% ylimregion = median(data(isnan(data)==0))
ylimregion = 300;
fillh1 = fill([0.1 11.9 11.9 0.1], [10 10 ylimregion.*ones(1,2)],0.9*ones(1,3));
fillh1.EdgeColor = 0.9*ones(1,3); fillh1.FaceAlpha = 0.5;
%ph(1)=plot(xx,yy(1,:),'-r','LineWidth',1); hold on,
ph(2)=plot(xx(1:end),yy(2,1:end),'-g','LineWidth',1);
ph(3)=plot(xx(1:end),yy(3,1:end),'-','Color',[0.7,0.7,1],'LineWidth',1);
ph(1) = [];
% plot([0,12],max(err(1:Nr,eta_vec==0.3,2)).*ones(1,2),'-k')
for iM = 2:3
    boxplot(err(:,:,iM),'PlotStyle','compact','Colors',ccolors(iM,:), 'Symbol', symbols{iM}, ... %, 'Labels', num2str(eta_vec'), eta_vec
        'Widths',0.1); hold on
    delete(findobj(gca,'Type','text'))
end
set(gca,'xtick',[1:2:11],'xticklabel',num2str(eta_vec(1:2:end)'),'Position',[0.2 0.22 0.75 0.73])
xlim([0 12]), ylim([0 1300])
xt = xlabel('Eta'); xt.Position = [115 -20 -0.1];
ylabel('MSE')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noise_',NARXtraining,'_noleg.eps']);

% l1 = legend(ph,'DMDc','SINDYc','NARX');
l1 = legend(ph,'SINDYc','NARX');
set(l1,'Location','NorthWest')
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noise_',NARXtraining,'.eps']);

%%
clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
symbols = {'o', 'd', 's'};
xx = 1:0.1:11;%0:0.01:0.6;
yy = zeros(3,length(xx));
err = zeros(Nr,N_ETA,3);
for iM = 1:3
    for i = 1:N_ETA, err(1:Nr,i,iM) = ResultsALL(i,iM).errM(1:Nr); end
    if iM == 2
        xm = zeros(1,N_ETA);
        for k = 1:N_ETA
            TF = isnan(err(1:Nr,k,iM));
            xm(k) = median(err(TF==0,k,iM),1);
        end
    else
        xm = median(err(1:Nr,:,iM),1);
    end
    yy(iM,:) = spline(1:11,xm,xx); %eta_vec(:)
end
figure; hold on
fillh1 = fill([0.1 11.9 11.9 0.1], [10 10 ylimregion.*ones(1,2)],0.9*ones(1,3));
fillh1.EdgeColor = 0.9*ones(1,3); fillh1.FaceAlpha = 0.5;
% ph(1)=plot(xx,yy(1,:),'-r','LineWidth',1); hold on,
ph(2)=plot(xx(1:end),yy(2,1:end),'-g','LineWidth',1);
ph(3)=plot(xx(1:end),yy(3,1:end),'-','Color',[0.7,0.7,1],'LineWidth',1);
ph(1) = [];
% plot([0,12],max(err(1:Nr,eta_vec==0.3,2)).*ones(1,2),'-k')
for iM = 2:3
    boxplot(err(:,:,iM),'PlotStyle','compact','Colors',ccolors(iM,:), 'Symbol', symbols{iM}, ... %, 'Labels', num2str(eta_vec'), eta_vec
        'Widths',0.1); hold on
    delete(findobj(gca,'Type','text'))
end
set(gca,'xtick',[1:2:11],'xticklabel',num2str(eta_vec(1:2:end)'),'Position',[0.2 0.22 0.75 0.73])
xlim([0 12]), ylim([0 1300])
xt = xlabel('Eta'); xt.Position = [115 -20 -0.1];
ylabel('MSE')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noiseM_',NARXtraining,'_noleg.eps']);

% l1 = legend(ph,'DMDc','SINDYc','NARX');
l1 = legend(ph,'SINDYc','NARX');
set(l1,'Location','NorthWest')
% l1.Position = [l1.Position(1)-0.06 l1.Position(2)-0.07 l1.Position(3)-0.01 l1.Position(4)];
% l1.Position = [l1.Position(1)+0.08 l1.Position(2)-0.08 l1.Position(3)-0.01 l1.Position(4)];
l1.Position = [l1.Position(1)-0.01 l1.Position(2)+0.02 l1.Position(3)-0.01 l1.Position(4)];
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noiseM_',NARXtraining,'.eps']);

%% Normalized error
clear LegLoc
LOG_SCALE = 0;
data2plot = zeros(Nr,N_ETA,3);
% for iN = 1:N_LENGTHS
for iNoise = 1:N_ETA
    for iR = 1:Nr
        for iM = 1:3
            data2plot(iR,iNoise,iM) = mean( sum( abs(ResultsALL(iNoise,iM).xA - ResultsALL(iNoise,iM).xB(:,:,iR))./abs(ResultsALL(iNoise,iM).xA) ,2)  );
            %data2plot(iR,iNoise,iM) = mean( abs(sum( (ResultsALL(iNoise,iM).xA - ResultsALL(iNoise,iM).xB(:,:,iR))./ResultsALL(iNoise,iM).xA ,2)) );
        end
    end
end
% end

PostName = 'RelErr';
ytext = 'Avg. Rel. Error';
yaxlim = [0 14];
shaded_region = 0;
VIZ_ERROR_STATS

%%
iNoise = 1; iM = 3;
figure, hold on
plot(ResultsALL(iNoise,iM).xA,'-k')
for i = 1:Nr
    plot(ResultsALL(iNoise,iM).xB(:,:,iR),'--r');%,'Color',[0.7,0.7,0.7])
end
%%
ModelSelection = {'DMDc','SINDYc', 'NARX'};
dt = 0.01;
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
tspan = [1:Nt].*dt+30;
iM = 3;
data2plot = zeros(Nt,N_ETA,3);

for iM = 2:3
    PostName = ['_TS_',ModelSelection{iM}];
    for iNoise = 1:N_ETA
        f1 = figure('visible','off');hold on, box on
        xshift = [-35,0,+20];
        for iVar = 1:Nvar
            %data2plot(iT,iNoise,iM) = (ResultsALL(iNoise,iM).xA - ResultsALL(iNoise,iM).xB(:,:,iR));
            %boxplot(squeeze(ResultsALL(iNoise,iM).xB(1:10:end,1,:))')
            Xstats = squeeze(ResultsALL(iNoise,iM).xB(1:1:end,iVar,:))';
            Ystats = prctile(Xstats,[25 50 75],1);

            X=[tspan(1:end),fliplr(tspan(1:end))];                %#create continuous x value array for plotting
            Y=[Ystats(1,:)+xshift(iVar),fliplr(Ystats(3,:))+xshift(iVar)];              %#create y values for out and then back
            fillh1 = fill(X,Y,ccolors(iM,:));                  %#plot filled area
            fillh1.EdgeColor = ccolors(iM,:);
            fillh1.FaceAlpha = 0.3;
            fillh1.EdgeAlpha = 0.3;

            plot(tspan,Ystats(2,:)+xshift(iVar),'-','Color',ccolors(iM,:),'LineWidth',1)
            %                 plot(tspan,Ystats(1,:),'--','Color',ccolors(iM,:))
            %                 plot(tspan,Ystats(3,:),'--','Color',ccolors(iM,:))
            plot(tspan,ResultsALL(iNoise,iM).xA(1:1:end,iVar)+xshift(iVar),'-k','LineWidth',1)%,'Color',[0.7,0.7,0.7])
        end
        ylim([-60 65])
        ylabel('xi'), xlabel('Time')
        set(gca,'LineWidth',1, 'FontSize',14)
        set(gcf,'Position',[100 100 300 200])
        set(gcf,'PaperPositionMode','auto'),
        drawnow
        print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_',SystemModel,'_PREDPERF_',InputSignalType,'_N',sprintf('%04g',Ntrain_vec(iN)),'_noise',sprintf('%03g',100*eta_vec(iNoise)),'_',NARXtraining,'_',PostName,'_noleg.eps']);
        close(f1);

    end
end

%% Training time
LOG_SCALE = 1;
data2plot = zeros(Nr,N_ETA,3);
% for iN = 1:N_LENGTHS
for iNoise = 1:N_ETA
    for iR = 1:Nr
        for iM = 1:3
            data2plot(iR,iNoise,iM) = ResultsALL(iNoise,iM).Ttraining(iR);
        end
    end
end
% end

clear YSCALE_SET
PostName = 'TrainTime';
ytext = 'Training Time [s]';
yaxlim = [0 max(data2plot(:))];
shaded_region = 0;
LegLoc = 'NorthEast';
VIZ_ERROR_STATS

% Zoom
% PostName = 'TrainTimeZOOM';
% yaxlim = [0 0.04];
% VIZ_ERROR_STATS
%% Prediction horizon
dt = 0.01;
clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
ERR_RADIUS = 0.1;
PredHor_Ball = zeros(Nr,N_ETA,Nmodels);
PredHor_Max = zeros(Nr,N_ETA,Nmodels);
%Err_SMAPE = zeros(Nr,N_ETA,Nmodels);
Err_Ratio = zeros(Nr,N_ETA,Nmodels);
for iM = 2:3
    for iNoise = 1:N_ETA
        for iR = 1:Nr
            tmp =  (ResultsALL(iNoise,iM).xA - ResultsALL(iNoise,iM).xB(:,:,iR))./ResultsALL(iNoise,iM).xA ;
            TF = abs(tmp)>0.1;
            idx = [];
            for iVar = 1:Nvar
                idx = [idx ; find(TF(:,iVar)==1,1,'first')];
            end
            PredHor_Max(iR,iNoise,iM) = max(idx);

%             tmp2 = sqrt(sum(tmp.^2,2));
            tmp2 = sqrt(sum(abs((ResultsALL(iNoise,iM).xA - ResultsALL(iNoise,iM).xB(:,:,iR))).^2,2));
            idx = find(tmp2>3,1,'first');
            PredHor_Ball(iR,iNoise,iM) = idx-1;

            %Err_SMAPE(iR,iNoise,iM) =  sum(sum(abs(sum(ResultsALL(iNoise,iM).xB(:,:,iR)-ResultsALL(iNoise,iM).xA)./sum(ResultsALL(iNoise,iM).xA + ResultsALL(iNoise,iM).xB(:,:,iR)))./Nt))./Nvar;
            Err_Ratio(iR,iNoise,iM) =  sum(sum( abs((ResultsALL(iNoise,iM).xB(:,:,iR)./ResultsALL(iNoise,iM).xA))  )./Nt)./Nvar;

        end
    end
end

%%
data2plot = (PredHor_Ball).*dt;
PostName = 'PredHoriz_Ball';
ytext = 'Pred. Horizon';
yaxlim = [10^-3 10^1];
%yaxlim = [0 10];
LOG_SCALE = 1;
YSCALE_SET = 1;
shaded_region = 0;
LegLoc = 'NorthEast';
VIZ_ERROR_STATS

return
data2plot = PredHor_Max;
PostName = 'PredHoriz_Max';
ytext = 'Pred. Horizon';
yaxlim = [10^-1 10^3];
LOG_SCALE = 1;
YSCALE_SET = 1;
shaded_region = 0;
LegLoc = 'NorthEast';
VIZ_ERROR_STATS

data2plot = Err_Ratio;
PostName = 'Err_Ratio';
ytext = 'Error';
yaxlim = [10^-1 10^1];
LOG_SCALE = 1;
YSCALE_SET = 1;
shaded_region = 0;
LegLoc = 'NorthEast';
VIZ_ERROR_STATS

author: eurika-kaiser

repo: SINDY-MPC


text:
function n = numel(varargin)
% returns the number of tensors

n = 1;

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function loadScenarioDD(filename_DD)
%
% Author:
%   Leonhard Hermansdorfer      Date: 03.12.2020
%
% Description:
%   function used to load a scenario into the raceline data dictionary which is then used
%   for trajectory planning emulation.
%
% Input:
%   filename_DD:        name of data dictionary where to load scenario from;
%                       this file must be created via loadRaceline()

%% read csv-file which contains the raceline

    % add .sldd to data dictionary name if not already included
    if ~contains(filename_DD, ".sldd")
        filename_DD = strcat(filename_DD,'.sldd');
    end

    % raise error when no raceline file with specified name is found
    if isempty(which(filename_DD))
        error(strcat('No raceline data dictionary found with name: "', filename_DD, '"!'))
    end

%% assign struct to data dictionary

    DDname = "raceline.sldd";

    try

        datadict_raceline = Simulink.data.dictionary.open(DDname);
        datadict_scenario = Simulink.data.dictionary.open(filename_DD);

        % delete old data dictionary content
        dDataSectObj_racelineDD = getSection(datadict_raceline,'Design Data');
        list_struct_entries = find(dDataSectObj_racelineDD, 'DataSource', DDname);

        % delete all data dict entries of type struct
        for i=1:length(list_struct_entries)
            deleteEntry(dDataSectObj_racelineDD, list_struct_entries(i).Name)
        end

        % write raceline, start position and raceline file info to data dictionary
        disp(strcat('Write scenario of "', filename_DD, '" to raceline data dictionary...'));

        % get data of scenario data dictionary
      	dDataSectObj_scenarioDD = getSection(datadict_scenario,'Design Data');
        list_struct_entries = find(dDataSectObj_scenarioDD, 'DataSource', filename_DD);

     	% copy all data dict entries to target data dict
        for i=1:length(list_struct_entries)
            temp = getEntry(dDataSectObj_scenarioDD, list_struct_entries(i).Name);
            addEntry(dDataSectObj_racelineDD, list_struct_entries(i).Name, getValue(temp));
        end

        addEntry(dDataSectObj_racelineDD, 'current_scenario', filename_DD)

        % save and close data dictionary
        datadict_raceline.saveChanges()
        datadict_raceline.close();
        datadict_scenario.close();

    catch e
        warning(['Something went wrong during configuration of ' ...
                'raceline data dictionary']);
        disp('******************************************');
        disp('Exception: ');
        disp(getReport(e))
        disp('******************************************');

    end

end

author: TUMFTM


text:
function [  ] = pipeline_release(  )
    caffe.reset_all();
end

repo: fashion-landmarks


text:
function [cmap] = buildcmap(colors)
% [cmap]=buildcmap(colors)
%
% This function can be used to build your own custom colormaps. Imagine if
% you want to display rainfall distribution map. You want a colormap which
% ideally brings rainfall in mind, which is not achiveved by colormaps such
% as winter, cool or jet and such. A gradient of white to blue will do the
% task, but you might also use a more complex gradient (such as
% white+blue+red or colors='wbr'). This function can be use to build any
% colormap using main colors rgbcmyk. In image processing, w (white) can be
% used as the first color so that in the output, the background (usually
% with 0 values) appears white. In the example of rainfall map, 'wb' will
% produce a rainfall density map where the background (if its DN values are
% 0) will appear as white.
%
% Inputs:
%  colors: string (char) of color codes, any sequence of rgbcmywk
%  representing different colors (such as 'b' for blue) is acceptable. If a
%  gradient of white to blue is needed, colors would be 'wb'; a rainbow of
%  white+blue+red+green would be 'wbrg'.
%
% Example:
%  [cmap]=buildcmap('wygbr');
% %try the output cmap:
% im=imread('cameraman.tif');
% imshow(im), colorbar
% colormap(cmap) %will use the output colormap
%
% First version: 14 Feb. 2013
% sohrabinia.m@gmail.com
%--------------------------------------------------------------------------

if nargin<1
    colors='wrgbcmyk';
end

if ~ischar(colors)
    error(['Error! colors must be a variable of type char with '...
        'color-names, such as ''r'', ''g'', etc., '...
        'type ''help buildcmap'' for more info']);
end

ncolors=length(colors)-1;

bins=round(255/ncolors);
% diff1=255-bins*ncolors;

vec=zeros(300,3);

switch colors(1)
    case 'w'
        vec(1,:)=1;
    case 'r'
        vec(1,:)=[1 0 0];
    case 'g'
        vec(1,:)=[0 1 0];
    case 'b'
        vec(1,:)=[0 0 1];
    case 'c'
        vec(1,:)=[0 1 1];
    case 'm'
        vec(1,:)=[1 0 1];
    case 'y'
        vec(1,:)=[1 1 0];
    case 'k'
        vec(1,:)=[0 0 0];
end

for i=1:ncolors
 beG=(i-1)*bins+1;
 enD=i*bins+1; %beG,enD
 switch colors(i+1)
     case 'w'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';%colors(i+1),beG,enD,
     case 'r'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';%colors(i+1),beG,enD
     case 'g'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';%colors(i+1),beG,enD
     case 'b'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';%colors(i+1),beG,enD
     case 'c'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';%colors(i+1),beG,enD
     case 'm'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';
     case 'y'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';
     case 'k'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';
 end
end
cmap=vec(1:bins*ncolors,:);
end %end of buildcmap

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function p = sphere_stereograph2_inverse ( q, focus, center )

%*****************************************************************************80
%
%% SPHERE_STEREOGRAPH2_INVERSE computes stereographic preimages of points.
%
%  Discussion:
%
%    We start with a sphere of center C.
%
%    F is a point on the sphere which is the focus of the mapping,
%    and the antipodal point 2*C-F is the point of tangency
%    to the sphere of a plane.
%
%    For any point Q on the plane, the stereographic inverse projection
%    P of the point is defined by drawing the line from F through Q, and
%    computing P as the intersection of this line with the sphere.
%
%    The spatial dimension M is arbitrary, but should be at least 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    C F Marcus,
%    The stereographic projection in vector notation,
%    Mathematics Magazine,
%    Volume 39, Number 2, March 1966, pages 100-102.
%
%  Parameters:
%
%    Input, real Q(M,N), the points, which are presumed to lie
%    on the plane.
%
%    Input, real FOCUS(M,1), the coordinates of the focus point.
%
%    Input, real CENTER(M,1), the coordinates of the center of the sphere.
%
%    Output, real P(M,N), the stereographic inverse projections of the points.
%
  [ m, n ] = size ( q );

  ff = repmat ( focus, 1, n );

  s(1,1:n) = 2.0 * ( center - focus )' * ( q - ff ) ./ sum ( ( q - ff ).^2, 1 );

  ss = repmat ( s, m, 1 );

  p = ss .* q + ( 1.0 - ss ) .* ff;

  return
end


text:
% Copyright 2017 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     https://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Demonstration of the image restoration experiments conducted in
% Y. Romano, M. Elad, and P. Milanfar, "The Little Engine that Could:
% Regularization by Denoising (RED)", submitted to SIAM Journal on Imaging
% Sciences, 2016. https://arxiv.org/abs/1611.02862
%
% This example reads a ground-truth image, degrades the image by
% first blurring or downscaling it, followed by an addition of random white
% Gaussian noise. Then it calls to RED in order to restore the image.
% This example compares the input and output PSNR, shows and saves the
% results. The suggested image-adaptive Laplacian-regularization functional
% is minimized using the Fixed-Point, ADMM, and Steepest Descent methods.
% Please refer to the paper for more details.
%
% The following are the degradation models that this example handles:
% 'UniformBlur'  - 9X9 uniform psf with noise-level equal to sqrt(2)
% 'GaussianBlur' - 25X25 Gaussian psf with std 1.6 and noise-level
%                  equal to sqrt(2)
% 'Downscale'    - 7X7 Gaussian psf with std 1.6 and noise-level equal to 5
%
% The denoising engine is TNRD: Yunjin Chen, and Thomas Pock, "Trainable
% Nonlinear Reaction Diffusion: A Flexible Framework for Fast and Effective
% Image Restoration", IEEE TPAMI 2016. The code is available in
% http://www.icg.tugraz.at/Members/Chenyunjin/about-yunjin-chen
% Note: Enable parallel pool to reduce runtime.
%
% The degradation process is similar to the one suggested in NCSR paper:
% Weisheng Dong, Lei Zhang, Guangming Shi, and Xin Li "Nonlocally
% Centralized Sparse Representation for Image Restoration", IEEE-TIP, 2013.
% The code is available in http://www4.comp.polyu.edu.hk/~cslzhang/NCSR.htm
%

clc;
clear;
close all;

% configure the path
% denoising functions
addpath(genpath('./tnrd_denoising/'));
% SD, FP, and ADMM methods
addpath(genpath('./minimizers/'));
% contains the default params
addpath(genpath('./parameters/'));
% contains basic functions
addpath(genpath('./helper_functions/'));
% test images for the debluring and super resolution problems,
% taken from NCSR software package
addpath(genpath('./test_images/'));

% set light_mode = true to run the code in a sub optimal but faster mode
% set light_mode = false to obtain the results reported in the RED paper
light_mode = false;

if light_mode
    fprintf('Running in light mode. ');
    fprintf('Turn off to obatain the results reported in RED paper.\n');
else
    fprintf('Light mode option is off. ');
    fprintf('Reproducing the result in RED paper.\n');
end

%% read the original image

file_name = 'starfish.tif';

fprintf('Reading %s image...', file_name);
orig_im = imread(['./test_images/' file_name]);
orig_im = double(orig_im);

fprintf(' Done.\n');

%% define the degradation model

% choose the secenrio: 'UniformBlur', 'GaussianBlur', or 'Downscale'
degradation_model = 'UniformBlur';

fprintf('Test case: %s degradation model.\n', degradation_model);

switch degradation_model
    case 'UniformBlur'
        % noise level
        input_sigma = sqrt(2);
        % filter size
        psf_sz = 9;
        % create uniform filter
        psf = fspecial('average', psf_sz);
        % use fft to solve a system of linear equations in closed form
        use_fft = true;
        % create a function-handle to blur the image
        ForwardFunc = ...
            @(in_im) imfilter(in_im,psf,'conv','same','circular');
        % the psf is symmetric, i.e., the ForwardFunc and BackwardFunc
        % are the same
        BackwardFunc = ForwardFunc;
        % special initialization (e.g. the output of other method)
        % set to identity mapping
        InitEstFunc = @(in_im) in_im;

    case 'GaussianBlur'
        % noise level
        input_sigma = sqrt(2);
        % filter size
        psf_sz = 25;
        % std of the Gaussian filter
        gaussian_std = 1.6;
        % create gaussian filter
        psf = fspecial('gaussian', psf_sz, gaussian_std);
        % use fft to solve a system of linear equations in closed form
        use_fft = true;
        % create a function handle to blur the image
        ForwardFunc = ...
            @(in_im) imfilter(in_im,psf,'conv','same','circular');
        % the psf is symmetric, i.e., the ForwardFunc and BackwardFunc
        % are the same
        BackwardFunc = ForwardFunc;
        % special initialization (e.g. the output of other method)
        % set to identity mapping
        InitEstFunc = @(in_im) in_im;

    case 'Downscale'
        % noise level
        input_sigma = 5;
        % filter size
        psf_sz = 7;
        % std of the Gaussian filter
        gaussian_std = 1.6;
        % create gaussian filter
        psf = fspecial('gaussian', psf_sz, gaussian_std);
        % scaling factor
        scale = 3;

        % compute the size of the low-res image
        lr_im_sz = [ceil(size(orig_im,1)/scale),...
                    ceil(size(orig_im,2)/scale)];
        % create the degradation operator
        H = CreateBlurAndDecimationOperator(scale,lr_im_sz,psf);
        % downscale
        ForwardFunc = @(in_im) reshape(H*in_im(:),lr_im_sz);
        % upscale
        BackwardFunc = @(in_im) reshape(H'*in_im(:),scale*lr_im_sz);
        % special initialization (e.g. the output of other method)
        % use bicubic upscaler
        InitEstFunc = @(in_im) imresize(in_im,scale,'bicubic');

    otherwise
        error('Degradation model is not defined');
end

%% degrade the original image

switch degradation_model
    case {'UniformBlur', 'GaussianBlur'}
        fprintf('Blurring...');
        % blur each channel using the ForwardFunc
        input_im = zeros( size(orig_im) );
        for ch_id = 1:size(orig_im,3)
            input_im(:,:,ch_id) = ForwardFunc(orig_im(:,:,ch_id));
        end
        % use 'seed' = 0 to be consistent with the experiments in NCSR
        randn('seed', 0);

    case 'Downscale'
        fprintf('Downscaling...');
        % blur the image, similar to the degradation process of NCSR
        input_im = Blur(orig_im, psf);
        % decimate
        input_im = input_im(1:scale:end,1:scale:end,:);
        % use 'state' = 0 to be consistent with the experiments in NCSR
        randn('state', 0);

    otherwise
        error('Degradation model is not defined');
end

% add noise
fprintf(' Adding noise...');
input_im = input_im + input_sigma*randn(size(input_im));

% convert to YCbCr color space if needed
input_luma_im = PrepareImage(input_im);
orig_luma_im = PrepareImage(orig_im);

if strcmp(degradation_model,'Downscale')
    % upscale using bicubic
    input_im = imresize(input_im,scale,'bicubic');
    input_im = input_im(1:size(orig_im,1), 1:size(orig_im,2), :);
end
fprintf(' Done.\n');
psnr_input = ComputePSNR(orig_im, input_im);

%% minimize the Laplacian regularization functional via Fixed Point

fprintf('Restoring using RED: Fixed-Point method\n');

switch degradation_model
    case 'UniformBlur'
        params_fp = GetUniformDeblurFPParams(light_mode, psf, use_fft);
    case 'GaussianBlur'
        params_fp = GetGaussianDeblurFPParams(light_mode, psf, use_fft);
    case 'Downscale'
        assert(exist('use_fft','var') == 0);
        params_fp = GetSuperResFPParams(light_mode);
    otherwise
        error('Degradation model is not defined');
end

[est_fp_im, psnr_fp] = RunFP(input_luma_im,...
                             ForwardFunc,...
                             BackwardFunc,...
                             InitEstFunc,...
                             input_sigma,...
                             params_fp,...
                             orig_luma_im);
out_fp_im = MergeChannels(input_im,est_fp_im);

fprintf('Done.\n');

%% minimize the Laplacian regularization functional via ADMM

fprintf('Restoring using RED: ADMM method\n');

switch degradation_model
    case 'UniformBlur'
        params_admm = GetUniformDeblurADMMParams(light_mode, psf, use_fft);
    case 'GaussianBlur'
        params_admm = GetGaussianDeblurADMMParams(light_mode, psf, use_fft);
    case 'Downscale'
        assert(exist('use_fft','var') == 0);
        params_admm = GetSuperResADMMParams(light_mode);
    otherwise
        error('Degradation model is not defined');
end

[est_admm_im, psnr_admm] = RunADMM(input_luma_im,...
                                   ForwardFunc,...
                                   BackwardFunc,...
                                   InitEstFunc,...
                                   input_sigma,...
                                   params_admm,...
                                   orig_luma_im);
out_admm_im = MergeChannels(input_im,est_admm_im);

fprintf('Done.\n');

%% minimize the laplacian regularization functional via Steepest Descent

fprintf('Restoring using RED: Steepest-Descent method\n');

switch degradation_model
    case 'UniformBlur'
        params_sd = GetUniformDeblurSDParams(light_mode);
    case 'GaussianBlur'
        params_sd = GetGaussianDeblurSDParams(light_mode);
    case 'Downscale'
        params_sd = GetSuperResSDParams(light_mode);
    otherwise
        error('Degradation model is not defined');
end

[est_sd_im, psnr_sd] = RunSD(input_luma_im,...
                             ForwardFunc,...
                             BackwardFunc,...
                             InitEstFunc,...
                             input_sigma,...
                             params_sd,...
                             orig_luma_im);
% convert back to rgb if needed
out_sd_im = MergeChannels(input_im,est_sd_im);

fprintf('Done.\n');

%% display final results

fprintf('Image name %s \n', file_name);
fprintf('Input PSNR = %f \n', psnr_input);
fprintf('RED: Fixed-Point PSNR = %f \n', psnr_fp);
fprintf('RED: ADMM PSNR = %f \n', psnr_admm);
fprintf('RED: Steepest-Decent PSNR = %f \n', psnr_sd);

%% write images

if ~exist('./results/','dir')
    mkdir('./results/');
end

fprintf('Writing the images to ./results...');

imwrite(uint8(input_im),['./results/input_' file_name]);
imwrite(uint8(out_fp_im),['./results/est_fp_' file_name]);
imwrite(uint8(out_admm_im),['./results/est_admm_' file_name]);
imwrite(uint8(out_sd_im),['./results/est_sd_' file_name]);

fprintf(' Done.\n');

author: google

repo: RED

save_path: github-repos/MATLAB/google-RED


text:
function [ pts ] = genRansacTestPoints( ptNum,outlrRatio,inlrStd,inlrCoef )
%GENRANSACTESTPOINTS Generate the points used by RANSAC function
%   PTS = GENRANSACTESTPOINTS(PTNUM,OUTLRRATIO,INLRSTD,INLRCOEF) PTS is
%   2*PTNUM, including PTNUM points, among which ROUND(OUTLRRATIO*PTNUM)
%   are outliers, others are inliers.
%		The inliers are around the line: y = INLRCOEF(1)*x + INLRCOEF(2),
%   INLRSTD is the standard deviation of, the dist between inliers and the
%	line. The outliers

outlrNum = round(outlrRatio*ptNum);
inlrNum = ptNum-outlrNum;

k = inlrCoef(1);
b = inlrCoef(2);
X = (rand(1,inlrNum)-.5)*ptNum; % X is in [-ptNum/2,ptNum/2]
Y = k*X+b;

% add noise for inliers
dist = randn(1,inlrNum)*inlrStd;
theta = atan(k);
X = X+dist*(-sin(theta));
Y = Y+dist*cos(theta);
inlrs = [X;Y];

outlrs = (rand(2,outlrNum)-.5)*ptNum;
% outlrs = (rand(2,outlrNum)-[ones(1,outlrNum)*.5;ones(1,outlrNum)*.1])*ptNum;
pts = [inlrs,outlrs];

end

author: Sable

repo: mcbench-benchmarks


text:
%IPROFILE Extract pixels along a line
%
% V = IPROFILE(IM, P1, P2) is a vector of pixel values extracted from the
% image IM (HxWxP) between the points P1 (2x1) and P2 (2x1).  V (NxP) has
% one row for each point along the line and the row is the pixel value
% which will be a vector for a multi-plane image.
%
% [P,UV] = IPROFILE(IM, P1, P2) as above but also returns the coordinates of
% the pixels for each point along the line.  Each row of UV is the pixel
% coordinate (u,v) for the corresponding row of P.
%
% Notes::
% - The Bresenham algorithm is used to find points along the line.
%
% See also BRESENHAM, ILINE.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function [p,uv] = iprofile(c, p1, p2)

    if nargin == 2
        % p1 is a set of points, in columns
        p = [];
        for i=2:numcols(p1)
            p = [p; iprofile(c, p1(:,i-1), p1(:,i))];
        end
        return
    end
    % coordinates must be integers
    p1 = round(p1); p2 = round(p2);

    points = bresenham(p1, p2);

    p = [];
    for point = points'
        p = [p; c(point(2), point(1), :)];
    end

    if nargout > 1
        uv = squeeze(points');
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [ scores, valid ] = objectClassifier(room_obj, objects, model, ClassNames )
%OBJECT_CLASSIFIER Summary of this function goes here
%   Detailed explanation goes here
% ClassNames = model.ClassNames;
% assert(length(ClassNames)==28, 'Classifier involve class number less than 28.');
% fprintf('object classifier: %d categories\n', length(ClassNames));

room_points = room_obj.out_points_w;
room_xyz = [min(room_points, [], 1) max(room_points, [],  1)];
total_fea_dim = comp_object_3d_feature();
num_object = length(objects);
features = zeros(num_object, total_fea_dim);
valid = true(num_object,1);

for oid = 1:num_object
    obj_points = objects(oid).out_points_w;
    obj_xyz = [min(obj_points, [], 1) max(obj_points, [], 1)];
    test_min = (room_xyz(1:3) - obj_xyz(1:3));
    norm_min = (room_xyz(1:3) - obj_xyz(1:3))./(obj_xyz(4:6)-obj_xyz(1:3)+0.00001);
    test_max = (obj_xyz(4:6) - room_xyz(4:6));
    norm_max = (obj_xyz(4:6) - room_xyz(4:6))./(obj_xyz(4:6)-obj_xyz(1:3)+0.00001);
    if any(norm_min>0.2) || any(norm_max>0.2) || any(test_min>20) || any(test_max>20)
        valid(oid) = false;
        continue;
    end

    obj_x = objects(oid).x_w;
    obj_ori_size = obj_x(4:6)';
    obj_angle = obj_x(7);
    features(oid,:) = comp_object_3d_feature(obj_xyz, room_xyz, obj_ori_size, obj_angle);
end

[~,scores,~] = predict(model, features);
if ~exist('ClassNames','var') || isempty(ClassNames)
    ClassNames = model.ClassNames;
end
I = zeros(1,length(ClassNames));
for i = 1:length(ClassNames)
    I(i) = str2double(ClassNames{i});
end
scores(:,I) = scores;

end

author: zouchuhang

repo: LayoutNet


text:
function [hfun] = limhfn2(vert,tria,hfun,dhdx)
%LIMHFN2 impose gradient limits on a discrete mesh-size fun-
%ction defined over a 2-simplex triangulation.
%   [HFUN] = LIMHFN2(VERT,TRIA,HFUN,DHDX) returns a "gradie-
%   nt-limited" function HFUN, defined over a triangulation
%   {VERT,TRIA}. HFUN is a T-by-1 vector of function values,
%   VERT is a V-by-2 array of XY coordinates and TRIA is a
%   T-by-3 array of triangles. Each row of TRIA
%   defines a triangle, such that VERT(TRIA(II,1),:), VERT(
%   TRIA(II,2),:) and VERT(TRIA(II,3),:) are the coordinates
%   of the II-TH triangle. DHDX is a scalar gradient-limit.
%   HFUN is "limited" to control variation over the elements
%   in the triangulation, such that (HFUN(V2)-HFUN(V1))/LL<=
%   DHDX, where {V1,V2} are the vertices of a given triangle
%   edge and LL is the edge-length. Limits are enforced exh-
%   austively over all edges.
%
%   See also TRIHFN2, LFSHFN2

%   This function is based on a very simplified version of:
%   Persson, P.O. "Mesh size functions for implicit geometr-
%   ies and PDE-based gradient limiting." Engineering with
%   Computers 22 (2006): 95-109.

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 18/04/2017

%---------------------------------------------- basic checks
    if ( ~isnumeric(vert) || ...
         ~isnumeric(tria) || ...
         ~isnumeric(hfun) || ...
         ~isnumeric(dhdx) )
        error('limhfn2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(vert) ~= +2 || ...
        ndims(tria) ~= +2 || ...
        ndims(hfun) ~= +2 || ...
        numel(dhdx) ~= +1 )
        error('limhfn2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(vert,2)~= +2 || ...
        size(tria,2) < +3 || ...
        size(hfun,2)~= +1 || ...
        size(vert,1)~= size(hfun,1) )
        error('limhfn2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nvrt = size(vert,1) ;

%---------------------------------------------- basic checks
    if (min(min(tria(:,1:3))) < +1 || ...
            max(max(tria(:,1:3))) > nvrt )
        error('limhfn2:invalidInputArgument', ...
            'Invalid TRIA input array.') ;
    end

%-------------------- impose gradient limits over mesh edges
   [edge,tria] = tricon2(tria);

    evec = vert(edge(:,2),:) - ...
           vert(edge(:,1),:) ;
    elen = sqrt(sum(evec.^2,2)) ;

%-------------------- impose gradient limits over edge-graph
   [hfun] = limgrad( ...
         edge,elen,hfun,dhdx,sqrt(nvrt)) ;

end

author: dengwirda


text:
function strides = compute_strides(bnet)
% COMPUTE_STRIDES For each CPT and each variable in that CPT,
% returns the stride of that variable.  So in future, we can
% quickly extract a slice of the CPT.
%
% The return value is a 2d array, where strides(i,j) contains the
% stride of the jth variable in the ith CPT.  Cell arrays would
% have saved space but they are slower.
%

num_cpts = size(bnet.CPD, 2);
max_cpt_dim = 1 + max(sum(bnet.dag));
strides = zeros(num_cpts, max_cpt_dim);

for i = 1:num_cpts
  c = CPT(bnet, i);
  siz = size(CPT(bnet, i));

  % Deal with the special case of a 1-d array separately
  if siz(2) == 1
    dim = 1;
  else
    dim = size(siz, 2);
  end

  strides(i, 1:dim ) = [1 cumprod(siz(1:dim-1))];
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function dlx = mat2dl_complex(x)
% Convert x into a particular data structure to store complex numbers
%
% function dlx = mat2dl_complex(x)
%
% The iput x can be defined recursively by arrays, structs and cells. Each
% part of x should contain complex numbers. The function converts each
% part of x into a struct containing dlarrays with fields real and imag
% which indicate the real and imaginary part of the stored complex numbers.
%
% See also: dl2mat_complex, manoptADhelp

% This file is part of Manopt: www.manopt.org.
% Original author: Xiaowen Jiang, July. 31, 2021.
% Contributors: Nicolas Boumal
% Change log:

    if ~isstruct(x) && ~iscell(x) && ~isnumeric(x)
        up = MException('manopt:autodiff:mat2dl_complex', ...
                    'mat2dl_complex should only accept a struct, a cell or a numeric array');
        throw(up);
    end

    % recursively convert each part of x into a particular struct
    if isstruct(x)
        dlx = mat2dl_struct(x);
    elseif iscell(x)
        dlx = mat2dl_cell(x);
    else
        xreal = real(x);
        ximag = imag(x);
        dlx.real = dlarray(xreal);
        dlx.imag = dlarray(ximag);
    end

    % convert x into a particular dlarray struct if x is a struct
    function dlx = mat2dl_struct(x)
        elems = fieldnames(x);
        nelems = numel(elems);
        for ii = 1:nelems
            if isstruct(x.(elems{ii}))
                dlx.(elems{ii}) = mat2dl_struct(x.(elems{ii}));
            elseif iscell(x.(elems{ii}))
                dlx.(elems{ii}) = mat2dl_cell(x.(elems{ii}));
            else
                dlx.(elems{ii}) = struct();
                xreal = real(x.(elems{ii}));
                ximag = imag(x.(elems{ii}));
                dlx.(elems{ii}).real = dlarray(xreal);
                dlx.(elems{ii}).imag = dlarray(ximag);
            end
        end
    end

    % convert x into a particular dlarray struct if x is a cell
    function dlx = mat2dl_cell(x)
        ncell = length(x);
        for ii = 1:ncell
            if isstruct(x{ii})
                dlx{ii} = mat2dl_struct(x{ii});
            elseif iscell(x{ii})
                dlx{ii} = mat2dl_cell(x{ii});
            else
                xreal = real(x{ii});
                ximag = imag(x{ii});
                dlx{ii} = struct();
                dlx{ii}.real = dlarray(xreal);
                dlx{ii}.imag = dlarray(ximag);
            end
        end
    end
end

author: NicolasBoumal

repo: manopt


text:
function str = cvx_replicate_structure( str, sz )
%CVX_REPLICATE_STRUCTURE Apply matrix structure to N-D arrays.

szs  = size( str );
nmat = prod( sz );
omat = ones( 1, nmat );
[ r, c, v ] = find( str );
nelm = length( r );
oelm = ones( 1, nelm );
r = r( : );
nvec = ( 0 : nmat - 1 ) * szs( 1 );
r = r( :, omat ) + nvec( oelm, : );
c = c( : );
c = c( :, omat );
nvec = ( 0 : nmat - 1 ) * szs( 2 );
c = c( :, omat ) + nvec( oelm, : );
v = v( : );
v = v( :, omat );
str = sparse( r, c, v, nmat * szs( 1 ), nmat * szs( 2 ) );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function inspect_qsubcellfun5

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

% this test script pertains to
% http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1255

%%

y = qsubcellfun(@inspect_qsubcellfun5_fn, {1, 2, 3}, 'memreq', 0, 'timreq', 500)
y = qsubcellfun(@inspect_qsubcellfun5_fn, {1, 2, 3}, 'memreq', 0, 'timreq', 500, 'stack', 3)

%%

fc = qsubcompile(@inspect_qsubcellfun5_fn)
y = qsubcellfun(fc, {1, 2, 3}, 'memreq', 0, 'timreq', 500)

%%

y = qsubcellfun(@inspect_qsubcellfun5_fn, {1, 2, 3}, 'memreq', 0, 'timreq', 500, 'compile', true, 'stack', 1)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear all, close all, clc

% J = @(u,t)(25-(5-(u-t))^2);

J = @(u,t)(25-(5-(u)).^2);
y0 = J(0,0);  % u = 0

% Extremum Seeking Control Parameters
freq = 10*2*pi; % sample frequency
dt = 1/freq;
T = 10; % total period of simulation (in seconds)
A = .2;  % amplitude
omega = 10*2*pi; % 10 Hz
phase = 0;
K = 5;   % integration gain

% High pass filter (Butterworth filter)
butterorder=1;
butterfreq=2;  % in Hz for 'high'
[b,a] = butter(butterorder,butterfreq*dt*2,'high')
ys = zeros(1,butterorder+1)+y0;
HPF=zeros(1,butterorder+1);

uhat=u;
for i=1:T/dt
    t = (i-1)*dt;
    yvals(i)=J(u,t);

    for k=1:butterorder
        ys(k) = ys(k+1);
        HPF(k) = HPF(k+1);
    end
    ys(butterorder+1) = yvals(i);
    HPFnew = 0;
    for k=1:butterorder+1
        HPFnew = HPFnew + b(k)*ys(butterorder+2-k);
    end
    for k=2:butterorder+1
        HPFnew = HPFnew - a(k)*HPF(butterorder+2-k);
    end
    HPF(butterorder+1) = HPFnew;

    xi = HPFnew*sin(omega*t + phase);
    uhat = uhat + xi*K*dt;
    u = uhat + A*sin(omega*t + phase);
    uhats(i) = uhat;
    uvals(i) = u;
end

%%
figure
subplot(2,1,1)
plot(t,uvals,t,uhats,'LineWidth',1.2)
l1=legend('$u$','$\hat{u}$')
set(l1,'interpreter','latex','Location','NorthWest')
grid on
subplot(2,1,2)
plot(t,yvals,'LineWidth',1.2)
ylim([-1 26])
grid on

set(gcf,'Position',[100 100 500 350])
set(gcf,'PaperPositionMode','auto')
% print('-depsc2', '-loose', '../../../figures/ESC_Response');

author: dynamicslab


text:
function obj = Infill_Maximal_Distance(x, sample_x)
obj = - min(pdist2(x, sample_x),[],2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [] = fseriesdemo(f,b,NN)
% FSERIESDEMO Plots a function f and it's Fourier series approximation.
% FSERIESDEMO(f,b,N) plots the function f and it's N term Fourier series
% approximation from zero to b.  FSERIESDEMO uses the FFT to approximate
% the Trigonometric Fourier Series coefficients.  FSERIESDEMO only works
% for the function from 0 to b.  Input N must be an integer >=2.
% Example:
%
%           f = @(x) sin(x)./x+(cos(x)).^2-exp(x/8)+x.^(x/20);
%           fseriesdemo(f,20,100)  % Approx. f on (0,20) with 100 terms.
%
%
% Author:  Matt Fig
% Contact:  popkenai@yahoo.com

% First some argument checking.
if nargin<3
   error ('Not enough arguments. See help.')
end

if ~isreal(NN) || ~isreal(b) || b<=0 || NN < 2 || floor(NN)~=NN
    error ('2nd and 3rd args must be positive real scalars. See help.')
end

n = 2*NN;
T = b-eps+(b-eps)/(n-2); % Get the appropriate endpoint for fft.
x = linspace(eps,T,n+1); % Use this x in fft.
x(end) = [];   % But not the endpoint.
fun = f(x);  % Get discrete points.
FUNC = fft(fun);  % Call fft
FUNC = [conj(FUNC(NN+1)) FUNC(NN+2:end) FUNC(1:NN+1)]/n; % Shift/scale.
A0 = FUNC(NN+1);   % First cosine term.
AN = 2*real(FUNC(NN+2:end));  % General cosine term.
BN = -2*imag(FUNC(NN+2:end));  % General sine term.

p = app1(x,A0,AN,BN,T);
% Now on to the plotting.
xx = linspace(eps,x(end),b*100);
xx(end) = [];
fun = f(xx);
hand = plot(xx,fun,'r',x,p,'b');
set(hand,'linewidth',2)
grid on;
axis tight;
legend('User Function',[num2str(NN),' Term Approximation'],0)

function answ = app1(x,A0,AN,BN,T)
ii = 1:length(AN);   % The indeces.
lgh = length(x);
[ii x] = meshgrid(ii,x); % Arrays for angles.
AN = repmat(AN,lgh,1);   % Create arrays for approximation.
BN = repmat(BN,lgh,1);
theta = ii.*x.*(2*pi)/T; % The angles.
answ = (sum(AN.*cos(theta)+BN.*sin(theta),2)+A0)'; % the F approximation.


text:
function [Y, G] = rbfOut(model, X);

% RBFOUT Output of an RBF model.
% FORMAT
% DESC gives the output of a radial basis function model, the function is
% a wrapper for rbffwd.
% ARG model : the model for which the output is required.
% ARG X : the input data for which the output is required.
% RETURN Y : the output.
%
% FORMAT
% DESC gives the output of a radial basis function model.
% ARG model : the model for which the output is required.
% ARG X : the input data for which the output is required.
% RETURN Y : the output.
% RETURN G : the hidden layer activations.
%
% SEEALSO : rbffwd, rbf, modelOut
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2007, 2008

% MLTOOLS

  if nargout > 1
    [Y, G] = rbffwd(model, X);
  else
    Y = rbffwd(model, X);
  end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function estimation_results = REKF_SLAM(data, NumberOfSteps)
% R-EKF SLAM

% load pre-given data: odometry and observations
if nargin < 1
    load('./data.mat');
end

data_matrix = data.state;

odom_sigma = data.odom_sigma;
obsv_sigma = data.obsv_sigma;

% odoCov = data.odom_cov;   % constant variable
% obsCov = data.obse_cov;   % constant variable

%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% In every step, all elements of Estimation_X will be changed %%%%%%%%%%%%
Estimation_X.orientation = data.poses.orientation(1:3,1:3);
Estimation_X.position    = data.poses.position(:,1);
Estimation_X.cov         = sparse(6,6);
Estimation_X.landmarks   = [];       % the landmarks observed until this step (included), 4*N format, the 4-th row is the index
Estimation_X0.IndexObservedNew=[];
Estimation_X0.IndexObservedAlreadyThis=[];
%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%

% Initialize
if nargin < 2
    NumberOfSteps = max(data_matrix(:,4));  % step instead of pose,  hence, it does not include pose 0
elseif NumberOfSteps > max(data_matrix(:,4))
    NumberOfSteps = max(data_matrix(:,4));
end
estimation_results = cell(1, NumberOfSteps+1);
estimation_results{1} = Estimation_X;
row_idx = (data_matrix(:, end) <= NumberOfSteps+1);
data_matrix = data_matrix(row_idx, :);

for i = 0:NumberOfSteps
    IndexOfCurrentStepInDataMatrix = find(data_matrix(:,4) == i);
    m = size(IndexOfCurrentStepInDataMatrix, 1);
    if ( mod(i, 50) == 0 )
        disp(['Processing pose ', int2str(i)]);
    end
    % det(Estimation_X.cov)
    if i==NumberOfSteps-1
    a=1;
    end

    if i ~= NumberOfSteps
        OdometryFromThis2Next = data_matrix(IndexOfCurrentStepInDataMatrix(m-5):IndexOfCurrentStepInDataMatrix(m),1);
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 1 ),...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 3 ),...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 5 )];
           [Estimation_X] = REKF_update(Estimation_X, CameraMeasurementThis, obsv_sigma );
        end

        estimation_results{i+1} = Estimation_X;

%        propagation using odometry info
        [Estimation_X] = REKF_propagate(Estimation_X, OdometryFromThis2Next, odom_sigma );

    else
        a=2;
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 1 ) , ...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 3 ),...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 5 )];
            CameraMeasurementThis=CameraMeasurementThis(1:end-6,:);

            [Estimation_X] = REKF_update(Estimation_X, CameraMeasurementThis, obsv_sigma );
        end
        estimation_results{i+1} = Estimation_X;
    end
end
clearvars -except estimation_results

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function test_values_test ( )

%*****************************************************************************80
%
%% TEST_VALUES_TEST tests the TEST_VALUES library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_VALUES_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the TEST_VALUES library.\n' );

  abram0_values_test ( );
  abram1_values_test ( );
  abram2_values_test ( );
  agm_values_test ( );
  airy_ai_values_test ( );
  airy_ai_int_values_test ( );
  airy_ai_prime_values_test ( );
  airy_bi_values_test ( );
  airy_bi_int_values_test ( );
  airy_bi_prime_values_test ( );
  airy_cai_values_test ( );
  airy_cbi_values_test ( );
  airy_gi_values_test ( );
  airy_hi_values_test ( );
  arccos_values_test ( );
  arccosh_values_test ( );
  arcsin_values_test ( );
  arcsinh_values_test ( );
  arctan_values_test ( );
  arctan_int_values_test ( );
  arctan2_values_test ( );
  arctanh_values_test ( );

  bei0_values_test ( );
  bei1_values_test ( );
  bell_values_test ( );
  ber0_values_test ( );
  ber1_values_test ( );
  bernoulli_number_values_test ( );
  bernoulli_poly_values_test ( );
  bernstein_poly_01_values_test ( );
  bessel_i0_values_test ( );
  bessel_i0_int_values_test ( );
  bessel_i0_spherical_values_test ( );
  bessel_i1_values_test ( );
  bessel_i1_spherical_values_test ( );
  bessel_in_values_test ( );
  bessel_ix_values_test ( );
  bessel_j0_values_test ( );
  bessel_j0_int_values_test ( );
  bessel_j0_spherical_values_test ( );
  bessel_j1_values_test ( );
  bessel_j1_spherical_values_test ( );
  bessel_jn_values_test ( );
  bessel_jx_values_test ( );
  bessel_k0_values_test ( );
  bessel_k0_int_values_test ( );
  bessel_k1_values_test ( );
  bessel_kn_values_test ( );
  bessel_kx_values_test ( );
  bessel_y0_values_test ( );
  bessel_y0_int_values_test ( );
  bessel_y0_spherical_values_test ( );
  bessel_y1_values_test ( );
  bessel_y1_spherical_values_test ( );
  bessel_yn_values_test ( );
  bessel_yx_values_test ( );
  beta_values_test ( );
  beta_cdf_values_test ( );
  beta_inc_values_test ( );
  beta_log_values_test ( );
  beta_noncentral_cdf_values_test ( );
  binomial_values_test ( );
  binomial_cdf_values_test ( );
  bivariate_normal_cdf_values_test ( );

  catalan_values_test ( );
  cauchy_cdf_values_test ( );
  cbrt_values_test ( );
  cheby_t_poly_values_test ( );
  cheby_u_poly_values_test ( );
  cheby_v_poly_values_test ( );
  cheby_w_poly_values_test ( );
  chi_values_test ( );
  chi_square_cdf_values_test ( );
  chi_square_noncentral_cdf_values_test ( );
  ci_values_test ( )
  cin_values_test ( );
  cinh_values_test ( );
  clausen_values_test ( );
  clebsch_gordan_values_test ( );
  collatz_count_values_test ( );
  cos_values_test ( );
  cos_degree_values_test ( );
  cos_power_int_values_test ( );
  cosh_values_test ( );
  cot_values_test ( );
  cp_values_test ( );

  dawson_values_test ( );
  debye1_values_test ( );
  debye2_values_test ( );
  debye3_values_test ( );
  debye4_values_test ( );
  dedekind_sum_values_test ( );
  dielectric_values_test ( );
  dilogarithm_values_test ( );

  e1_values_test ( );
  ei_values_test ( );
  elliptic_ea_values_test ( );
  elliptic_em_values_test ( );
  elliptic_ka_values_test ( );
  elliptic_km_values_test ( );
  erf_values_test ( );
  erfc_values_test ( );
  euler_number_values_test ( );
  euler_poly_values_test ( );
  exp_values_test ( );
  exp3_int_values_test ( );
  exponential_cdf_values_test ( );
  extreme_values_cdf_values_test ( );

  f_cdf_values_test ( );
  f_noncentral_cdf_values_test ( );
  fresnel_cos_values_test ( );
  fresnel_sin_values_test ( );
  frobenius_number_data_values_test ( );
  frobenius_number_order_values_test ( );
  frobenius_number_order2_values_test ( );

  gamma_values_test ( );
  gamma_cdf_values_test ( );
  gamma_inc_values_test ( );
  gamma_inc_p_values_test ( );
  gamma_inc_q_values_test ( );
  gamma_inc_tricomi_values_test ( );
  gamma_log_values_test ( );
  gegenbauer_poly_values_test ( );
  geometric_cdf_values_test ( );
  goodwin_values_test ( );
  gud_values_test ( );

  hermite_function_values_test ( );
  hermite_poly_phys_values_test ( );
  hermite_poly_prob_values_test ( );
  hyper_1f1_values_test ( );
  hyper_2f1_values_test ( );
  hypergeometric_cdf_values_test ( );
  hypergeometric_pdf_values_test ( );
  hypergeometric_u_values_test ( );

  i0ml0_values_test ( );
  i1ml1_values_test ( );
  i4_factorial_values_test ( );
  i4_factorial2_values_test ( );
  i4_fall_values_test ( );
  i4_rise_values_test ( );
  int_values_test ( );

  jacobi_cn_values_test ( );
  jacobi_dn_values_test ( );
  jacobi_poly_values_test ( );
  jacobi_sn_values_test ( );
  jed_ce_values_test ( );
  jed_mjd_values_test ( );
  jed_rd_values_test ( );
  jed_weekday_values_test ( );

  kei0_values_test ( );
  kei1_values_test ( );
  ker0_values_test ( );
  ker1_values_test ( );

  laguerre_associated_values_test ( );
  laguerre_general_values_test ( );
  laguerre_polynomial_values_test ( );
  lambert_w_values_test ( );
  laplace_cdf_values_test ( );
  legendre_associated_values_test ( );
  legendre_associated_normalized_values_test ( );
  legendre_associated_normalized_sphere_values_test ( );
  legendre_poly_values_test ( );
  legendre_function_q_values_test ( );
  lerch_values_test ( );
  lobachevsky_values_test ( );
  lobatto_polynomial_values_test ( );
  lobatto_polynomial_derivatives_test ( );
  log_values_test ( );
  log_normal_cdf_values_test ( );
  log_series_cdf_values_test ( );
  log10_values_test ( );
  logarithmic_integral_values_test ( );
  logistic_cdf_values_test ( );

  mertens_values_test ( );
  moebius_values_test ( );

  negative_binomial_cdf_values_test ( );
  nine_j_values_test ( );
  normal_cdf_values_test ( );
  normal_01_cdf_values_test ( );

  omega_values_test ( );
  owen_values_test ( );

  partition_count_values_test ( );
  partition_distinct_count_values_test ( );
  phi_values_test ( );
  pi_values_test ( );
  poisson_cdf_values_test ( );
  polylogarithm_values_test ( );
  prandtl_values_test ( );
  prime_values_test ( );
  psat_values_test ( );
  psi_values_test ( );

  r8_factorial_values_test ( );
  r8_factorial_log_values_test ( );
  r8_factorial2_values_test ( );
  r8_fall_values_test ( );
  r8_rise_values_test ( );
  rayleigh_cdf_values_test ( );

  secvir_values_test ( );
  shi_values_test ( );
  si_values_test ( );
  sigma_values_test ( );
  sin_values_test ( );
  sin_degree_values_test ( );
  sin_power_int_values_test ( );
  sinh_values_test ( );
  six_j_values_test ( );
  sound_values_test ( );
  sphere_unit_area_values_test ( );
  sphere_unit_volume_values_test ( );
  spherical_harmonic_values_test ( );
  sqrt_values_test ( );
  stirling1_values_test ( );
  stirling2_values_test ( );
  stromgen_values_test ( );
  struve_h0_values_test ( );
  struve_h1_values_test ( );
  struve_l0_values_test ( );
  struve_l1_values_test ( );
  student_cdf_values_test ( );
  student_noncentral_cdf_values_test ( );
  subfactorial_values_test ( );
  surten_values_test ( );
  synch1_values_test ( );
  synch2_values_test ( );

  tan_values_test ( );
  tanh_values_test ( );
  tau_values_test ( );
  thercon_values_test ( );
  three_j_values_test ( );
  tran02_values_test ( );
  tran03_values_test ( );
  tran04_values_test ( );
  tran05_values_test ( );
  tran06_values_test ( );
  tran07_values_test ( );
  tran08_values_test ( );
  tran09_values_test ( );
  trigamma_values_test ( );
  truncated_normal_ab_cdf_values_test ( );
  truncated_normal_ab_pdf_values_test ( );
  truncated_normal_a_cdf_values_test ( );
  truncated_normal_a_pdf_values_test ( );
  truncated_normal_b_cdf_values_test ( );
  truncated_normal_b_pdf_values_test ( );
  tsat_values_test ( );

  van_der_corput_values_test ( );
  viscosity_values_test ( );
  von_mises_cdf_values_test ( );

  weekday_values_test ( );
  weibull_cdf_values_test ( );

  zeta_values_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_VALUES_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer


text:
%% Example of using KST class for interfacing with KUKA iiwa robots
% This example script is used to show how to utilise the networking
% functions integrated into the class KST of the KUKA Sunrise Toolbox

% Before using the script:
% 1- You need to establish a network between the PC and the robot
%       controller
% 2- Synchronise the application "MatlabToolboxServer" into the robot
%       controller using the SunriseWorkbench, from more info consult the
%       document "Import KST to SunriseWorkbench.pdf"
% 3- Run the application "MatlabToolboxServer" from the teach pendant of
%       the robot.
% 4- Run the following script from MATLAB on a PC that is paired to the
%       controller of the robot using an ethernet cable connected to the
%       X66 connector.

% Note: When you run the "MatlabToolboxServer" app, you have 60 seconds to
% connect to the robot from MATLAB, if the limit time is exceeded without
% establishing a connection from MATLAB, the "MatlabToolboxServer" app
% turns off automatically. You have to manyally restart the
% "MatlabToolboxServer" app before connecting to the robot again.

% Exceptions, to use the instruction "net_pingIIWA", or the general porpuse
% functions it is not required to run the "MatlabToolboxServer" app.
% Yet a pairing the PC and the controller on ethernet is still important to
% run the function "net_pingIIWA".

% Copyright Mohammad SAFEEA, 30th-June-2018

close all;clear;clc;
warning('off')

%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object
%% Check if the robot is reachable
rep=iiwa.net_pingIIWA();
disp(rep);

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
  return;
end
pause(1);
disp('This script demonstrates the networking functions of KST')

%% Read EEF position from the controller
Pos=iiwa.getEEFPos();
disp('EEF position and orientation is:')
disp(Pos)

%% turn off the server
iiwa.net_turnOffServer();
warning('on')

repo: KST-Kuka-Sunrise-Toolbox


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logistic regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell you
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

k = num_labels;
for i = 1:k

    % Set Initial theta
    initial_theta = zeros(n + 1, 1);

    % Set options for fminunc
    options = optimset('GradObj', 'on', 'MaxIter', 50);

    % Run fmincg to obtain the optimal theta
    % This function will return theta and the cost
    [initial_theta] = fmincg (@(t)(lrCostFunction(t, X, (y == i), lambda)),initial_theta, options);

    all_theta(i,:) = initial_theta';
end

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

end

author: zzlyw

repo: machine-learning-exercises


text:
function [r,lon,lat] = SHMapToGrid(vec,res,lmax,location)

% [r,lon,lat] = SHMapToGrid(vec,res[,lmax,location])
%
% Maps a vector of spherical harmonic coefficients to grid.
% To use with multi-layer vectors, specify the array lmax.
% Grid is generated with the default resolution 10 degrees.
% By default, the spherical harmonics are computed at the
% centers of the cells. This behaviour can be modified by
% setting location='corner' and any value of res in degrees

if nargin < 2
    res = 10;
end

if nargin < 3
    lmax = SHn2lm(length(vec));
end

if nargin < 4
    location = 'center';
end

if strcmp(location,'center')
    shift = res/2;
    lon=linspace(shift,360-shift,360/res);
    lat=linspace(90-shift,-90+shift,180/res)';
    colat = 90 - lat;
elseif strcmp(location,'corner')
    lon=linspace(0,360,360/res+1);
    lat=linspace(90,-90,180/res+1)';
    colat = 90 - lat;
else
    error('location argument must be "center" or "corner"');
end

nlon=length(lon);
nlat=length(lat);
nrad=length(lmax);

r=zeros(nlat,nlon,nrad);

for i=1:nlon
    for j=1:nlat
        [SHvec,i1,i2] = SHCreateYVec(lmax,lon(i),colat(j),'deg');
        for k=1:nrad
            r(j,i,k)= vec(i1(k):i2(k))'*SHvec(i1(k):i2(k));
        end
    end
end

% if nargout > 1
%     [lon,lat] = meshgrid(lon,lat);
% end


text:
function [track_info,pos_info_val] = bz_olypherInfo(data,round_to,smoothing)
% USAGE
%        [track_info,pos_info_val] = bz_olypherInfo(data,round_to,smoothing)
%
% INPUTS
%         data - matrix (M x N x D) where M is the number of cells to analyze,
%                N is the number of trials for each cell, and D is the number
%                of time bins
%         round_to - integer value that data is discritized to. A value of
%                    2 means all data will be rounded to nearest 2 (i.e.
%                    2,4,8...)
%         smoothing - 0 no smoothing, else smooth with N bins
% OUTPUTS
%         track_info - matrix (M x D-2) of information scores across all
%                      trials or behavior windows (N)
%         pos_info_val - matrix (M x N x D) of all information values
%                        that are calculated
%
% this function calculates the information carried in the firing rate of single
% neurons per spatial/temporal bin
%
%  Written by David Tingley
%  UCSD Cognitive Neuroscience
%  1/15/12

%% TODO
% - convert to varargin with input parser
% - add 'exclude' input to remove 0's from info calculation
% -

M = size(data,1);
N = size(data,2);
D = size(data,3);

if M == 0
    M = 1
end
if N == 0
    N = 1
end
if D == 0
    D = 1
end
pos_info_val = zeros(M,N,D);
a = N*D;

%% Rounding
if smoothing ~= 0
    for i = 1 : M
        for k = 1:N
            data(i,k,:) = smooth(squeeze(data(i,k,:)),smoothing).*smoothing;
        end
    end
end
data = round(data./round_to)*round_to;

%% Info Analysis

for i = 1 : M
      for x = 1 : D
         for k = 1:N

                q = data(i,k,x);

                pKx = ((length(find(data(i,:,x) == q)))/N);

                pK = ((length(find(data(i,:,:) == q)))/a);

                if pK == 0 || pKx == 0 || pKx < pK
                    pos_info_val(i,k,x) = pos_info_val(i,k,x);
                else
                    pos_info_val(i,k,x) = pos_info_val(i,k,x) + (pKx*log2(pKx/pK));
                end

         end
     end

end

for i = 1:M
track_info(i,:) = sum(pos_info_val(i,:,2:end-1),2);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%This Matlab script can be used to reproduce Figure 7.3 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.01 (Last edited: 2019-03-27)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Number of BS antennas
M = 100;

%Define the pilot reuse factor
f = 2;

%Select the number of setups with random UE locations
nbrOfSetups = 100;

%Select the number of channel realizations per setup
nbrOfRealizations = 200;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Define noise figure at BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Set the range of Delta values in the power control policy of (7.10)
DeltadB = 0:10:20;

%Generate uncorrelated covariance matrices
Runcorr = repmat(eye(M),[1 1 K L L]);

%Prepare to save simulation results
sumSE_MR = zeros(L*K*nbrOfSetups,length(DeltadB));
sumSE_ZF = zeros(L*K*nbrOfSetups,length(DeltadB));
sumSE_SMMSE = zeros(L*K*nbrOfSetups,length(DeltadB));
sumSE_RZF = zeros(L*K*nbrOfSetups,length(DeltadB));
sumSE_MMMSE = zeros(L*K*nbrOfSetups,length(DeltadB));

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,K,M,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoiseOriginal = channelGaindB - noiseVariancedBm;

    %Go through all values of Delta
    for s = 1:length(DeltadB)

        %Extract the average channel gains before power control
        channelGainOverNoise = channelGainOverNoiseOriginal;

        %Go through all cells
        for j = 1:L

            %Compute beta_j,min in the power control policy of (7.10)
            betajMin = min(channelGainOverNoiseOriginal(:,j,j));

            %Scale the average channel gains by applying the power control
            %policy of (7.10). Note that we are including this power
            %control here so that we can then view it as if all UEs of
            %transmitting at maximum power
            differenceSNR = channelGainOverNoiseOriginal(:,j,j)-betajMin;
            backoff = differenceSNR-DeltadB(s);
            backoff(backoff<0) = 0;

            channelGainOverNoise(:,j,:) = channelGainOverNoiseOriginal(:,j,:)-repmat(backoff,[1 1 L]);

        end

        %Generate channel realizations with estimates and estimation
        %error correlation matrices
        [Hhat,C,tau_p,Rscaled] = functionChannelEstimates(R,channelGainOverNoise,nbrOfRealizations,M,K,L,p,f);

        %Compute SEs using Theorem 4.1
        [SE_MR,SE_ZF,SE_SMMSE,SE_RZF,SE_MMMSE] = functionComputeSE_UL(Hhat,C,Rscaled,tau_c,tau_p,nbrOfRealizations,M,K,L,p);

        %Save results
        sumSE_MR(1+(n-1)*K*L:n*K*L,s) = SE_MR(:);
        sumSE_ZF(1+(n-1)*K*L:n*K*L,s) = SE_ZF(:);
        sumSE_SMMSE(1+(n-1)*K*L:n*K*L,s) = SE_SMMSE(:);
        sumSE_RZF(1+(n-1)*K*L:n*K*L,s) = SE_RZF(:);
        sumSE_MMMSE(1+(n-1)*K*L:n*K*L,s) = SE_MMMSE(:);

        %Delete large matrices
        clear Hhat C Rscaled;

    end

    %Delete large matrices
    clear R;

end

%% Plot the simulation results

figure;
hold on; box on;

plot(sort(sumSE_MR(:,1)),linspace(0,1,K*L*nbrOfSetups),'k-','LineWidth',1);
plot(sort(sumSE_MR(:,2)),linspace(0,1,K*L*nbrOfSetups),'b-.','LineWidth',1);
plot(sort(sumSE_MR(:,3)),linspace(0,1,K*L*nbrOfSetups),'r--','LineWidth',1);
xlabel('SE per UE [bit/s/Hz]');
ylabel('CDF');
legend('\Delta=0 dB','\Delta=10 dB','\Delta=20 dB','Location','SouthEast');
xlim([0 8]);

figure;
hold on; box on;

plot(sort(sumSE_RZF(:,1)),linspace(0,1,K*L*nbrOfSetups),'k-','LineWidth',1);
plot(sort(sumSE_RZF(:,2)),linspace(0,1,K*L*nbrOfSetups),'b-.','LineWidth',1);
plot(sort(sumSE_RZF(:,3)),linspace(0,1,K*L*nbrOfSetups),'r--','LineWidth',1);
xlabel('SE per UE [bit/s/Hz]');
ylabel('CDF');
legend('\Delta=0 dB','\Delta=10 dB','\Delta=20 dB','Location','SouthEast');
xlim([0 8]);

figure;
hold on; box on;

plot(sort(sumSE_MMMSE(:,1)),linspace(0,1,K*L*nbrOfSetups),'k-','LineWidth',1);
plot(sort(sumSE_MMMSE(:,2)),linspace(0,1,K*L*nbrOfSetups),'b-.','LineWidth',1);
plot(sort(sumSE_MMMSE(:,3)),linspace(0,1,K*L*nbrOfSetups),'r--','LineWidth',1);
xlabel('SE per UE [bit/s/Hz]');
ylabel('CDF');
legend('\Delta=0 dB','\Delta=10 dB','\Delta=20 dB','Location','SouthEast');
xlim([0 8]);

author: emilbjornson

repo: massivemimobook


text:
% VL_HOMKERMAP Homogeneous kernel map
%   V = VL_HOMKERMAP(X, N) computes a 2*N+1 dimensional approximated
%   kernel map for the Chi2 kernel. X is an array of data points. Each
%   point is expanded into a vector of dimension 2*N+1 and saved to
%   the output V. The expanded feature vectors are stacked along the
%   first dimension, so that the output array V has the same
%   dimensions of the input array X except for the first one, which is
%   2*N+1 times larger.
%
%   The function accepts the following options:
%
%   KChi2::
%     Compute the map for the Chi2 kernel.
%
%   KINTERS::
%     Compute the map for the intersection kernel.
%
%   KL1::
%     Same as KINTERS, but deprecated as the name is not fully
%     accurate.
%
%   KJS::
%     Compute the map for the JS (Jensen-Shannon) kernel.
%
%   Period:: [automatically tuned]
%     Set the period of the kernel specturm. The approximation is
%     based on periodicizing the kernel specturm. If not specified,
%     the period is automatically set based on the heuristic described
%     in [2].
%
%   Window:: [RECTANGULAR]
%     Set the window used to truncate the spectrum before The window
%     can be either RECTANGULAR or UNIFORM window. See [2] and the API
%     documentation for details.
%
%   Gamma:: [1]
%     Set the homogeneity degree of the kernel. The standard kernels
%     are 1-homogeneous, but sometimes smaller values perform better
%     in applications. See [2] for details.
%
%   Example::
%     The following code results in approximatively the same
%     similarities matrices between points X and Y:
%
%       x = rand(10,1) ;
%       y = rand(10,100) ;
%       psix = vl_homkermap(x, 3) ;
%       psiy = vl_homkermap(y, 3) ;
%       figure(1) ; clf ;
%       ker = vl_alldist(x, y, 'kchi2') ;
%       ker_ = psix' * psiy ;
%       plot([ker ; ker_]') ;
%
%   Note::
%     The homogeneous kernels K(X,Y) are normally defined for
%     non-negative data only. VL_HOMKERMAP defines them for both
%     positive and negative data by using the definition
%     SIGN(X)SIGN(Y)K(ABS(X),ABS(Y)) -- note that other extensions are
%     possible as well (see [2]).
%
%   REFERENCES::
%     [1] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     Proc. CVPR, 2010.
%
%     [2] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     PAMI, 2011 (submitted).
%
%   See also: VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: Cloud-CV

repo: object-proposals


text:
function [stat, res, stat0] = get_validregions(stat0, res0, clustrules)

ix = 0;

iclust = res0.iclust;
% ipixbad = true(size(res0.M));
for k = 1:length(stat0)

    if ~isempty(stat0(k).region)

%         vs = [stat0(k).region.V];
%         criterion = sum(vs > mean(vs(:)) * clustrules.parent.PixelFractionThreshold);
        if 1 %criterion<clustrules.parent.MaxRegions

            igood = [stat0(k).region.mrs]./[stat0(k).region.mrs0] <clustrules.Compact;
            igood = igood & [stat0(k).region.npix]>clustrules.MinNpix & ...
                [stat0(k).region.npix]<clustrules.MaxNpix;
%             igood = [stat0(k).region.npix]>clustrules.MinNpix;

            igood = find(igood);
            for j = 1:length(igood)
                ix = ix+1;
                stat(ix)             = stat0(k).region(igood(j));

                iclust(stat(ix).ipix) = ix + numel(stat0);

                ismmb = ismember(stat0(k).ipix, stat(ix).ipix);
                stat0(k).ipix(ismmb) = [];
            end
        end
    end
    stat0(k).mrs   = Inf;

end

for k = 1:length(stat0)
    stat0(k).igood = 0;
end
for k = 1:length(stat)
    stat(k).igood = 1;
end

stat(1).region = [];
stat0(1).V = [];

stat0(1).parent = [];

res = res0;
res.iclust = iclust;

% stat(1).ipix = find(res.iclust(:)==1);
% stat(1).mrs = Inf;
% stat(1).mrs0 = Inf;
% stat(1).npix = numel(stat(1).ipix);
% stat(1).lambda = res.lambda(stat(1).ipix);
% stat(1).V = Inf;

stat = [stat0 stat];

author: cortex-lab


text:
validationMatrix = [1     1     0     0     0
                    1     0     1     0     0
                    0     0     0     0     0
                    0     0     0     1     0
                    0     0     0     0     1
                    0     0     0     0     0
                    0     0     0     0     1
                    0     0     0     0     0];

clusterer = NaiveClustererX();

clusters_test = clusterer.cluster(validationMatrix);

load('eval.mat');
assert(isequaln(clusters_test,clusters_eval));

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function  [frames, names, height,width,nframe ] = readAllFrames(options)

    frameFiles = imdir(fullfile(options.infolder));
    nframe = length( frameFiles);
    frames = cell( nframe, 1 );
    names = cell( nframe, 1 );
    for index = 1: nframe
        [~, frameName] = fileparts(frameFiles(index).name);
        if exist(fullfile(options.infolder, [frameName '.png']),'file')
            frame = imread(fullfile(options.infolder, [frameName '.png']));
        elseif exist(fullfile(options.infolder, [frameName '.jpg']),'file')
            frame = imread(fullfile(options.infolder, [frameName '.jpg']));
        elseif exist(fullfile(options.infolder, [frameName '.bmp']),'file')
            frame = imread(fullfile(options.infolder, [frameName '.bmp']));
        end
        frames{ index } = double(frame);
        names{ index } =  frameName;
    end
    [ height,width ] = size(frame(:,:,1));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

steps = [ 0.01 0.03 0.1 0.3 1 3 10 30 ];
minError = Inf;
minC = Inf;
minSigma = Inf;

for i = 1:length(steps)
    for j = 1:length(steps)
        currentC = steps(i);
        currentSigma = steps(j);
        model = svmTrain(X, y, currentC, @(x1, x2) gaussianKernel(x1, x2, currentSigma));
        predictions = svmPredict(model, Xval);
        error = mean(double(predictions ~= yval));

        if error < minError
            minError = error;
            minC = currentC;
            minSigma = currentSigma;
        end
    end
end

C = minC;
sigma = minSigma;

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

end

repo: MachineLearning


text:
%  Figure 10.45      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% fig10_45.m is a script to generate Fig. 10.45 the step response for the SRL
% design for the aircraft with state feedback augmenting the
% inner loop stabilization
clf;
%system matrices
ftq =[-0.0064    0.0263         0  -32.2000      0;
   -0.0941   -0.6240  761.1400 -196.2000         0;
   -0.0002   -0.0015   -4.4120  -12.4800         0;
         0         0    1.0000         0         0;
         0   -1.0000         0  830.0000         0];
g=[0; -32.7; -2.08; 0; 0];
h=[0 0 0 0 1];
j=0;
k=[-.0011 .0016 -.0833 -1.613 -.0010];
fc=ftq-g*k;
% dc gain
dcgain=-h*inv(fc)*g;
r=100/dcgain;
t=0:.3:30;
% step response
y=step(fc,g*r,h,j,1,t);
plot(t,y);
v=[0, 30, 0, 100];
axis(v);
xlabel('Time (sec)');
ylabel('Altitude, h (ft)');
grid;

title( 'Fig. 10.45 Step response of the altitude autopilot')

author: Sable

repo: mcbench-benchmarks


text:
%% Author: epokh
%% Website: www.epokh.org/drupy
%% This software is under GPL

clf
clear
%%Cartesian manipulator with an RPY wrist
px=2.3;
py=1.7;
pz=7.4;

fia=32;
fio=178;
fin=4;
Tcartesian=Tras(px,py,pz);
Torient=RPY(fia,fio,fin);

%%The trasformation matrix is:
Tend=Tcartesian*Torient;
figure(1);
plot3(0,0,0,'r');
plotT(Tend);
title('Cartesian manipulator example');
%%Cylindrical manipulator with an Euler wrist
Tcyl=Tcyl(62,8.2,5.2)*Euler(32,15,17);
figure(2);
plot3(0,0,0,'r');
plotT(Tcyl);
title('Cylindrical manipulator example');
%%Spherical maninupaltor with an Euler an RPY wrist
Tsfer=Tsfer(40,50,10)*RPY(10,0,10);
figure(3);
plot3(0,0,0,'r');
plotT(Tsfer);
title('Spherical manipulator example');

%%Invert the Euler transformation
euwrist=Euler(32,15,17);
[fi1,fio,fi2]=invEuler(euwrist);
fprintf('Invers euler angles %d %d %d \n',fi1,fio,fi2);
rpywrist=RPY(5,10,15);
[fia,fio,fin]=invRPY(rpywrist);
fprintf('Invers rpy angles %d %d %d \n',fia,fio,fin);


text:
function [varargout] = nanmax(varargin)
%NANMAX Maximum value, ignoring NaNs.
%   M = NANMAX(A) returns the maximum of A with NaNs treated as missing.
%   For vectors, M is the largest non-NaN element in A.  For matrices, M is
%   a row vector containing the maximum non-NaN element from each column.
%   For N-D arrays, NANMAX operates along the first non-singleton
%   dimension.
%
%   [M,NDX] = NANMAX(A) returns the indices of the maximum values in A.  If
%   the values along the first non-singleton dimension contain more than
%   one maximal element, the index of the first one is returned.
%
%   M = NANMAX(A,B) returns an array the same size as A and B with the
%   largest elements taken from A or B.  Either one can be a scalar.
%
%   [M,NDX] = NANMAX(A,[],DIM) operates along the dimension DIM.
%
%   See also MAX, NANMIN, NANMEAN, NANMEDIAN, NANMIN, NANVAR, NANSTD.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.8.1 $  $Date: 2010/03/16 00:15:49 $

% Call [m,ndx] = max(a,b) with as many inputs and outputs as needed
[varargout{1:nargout}] = max(varargin{:});

author: marianux

repo: ecg-kit

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


text:
function [metabolite_structure,IDsAdded] = parseKeggWebpage(metabolite_structure,startSearch,endSearch)

annotationSource = 'Kegg website';
annotationType = 'automatic';

Mets = fieldnames(metabolite_structure);
a = 1;
IDsAdded = '';

if ~exist('endSearch','var')
    endSearch = length(Mets);
end

if ~exist('startSearch','var')
    startSearch = 1;
end

% I am not going to collect pubchemId's as kegg has deposited its own id's
% to pubchem in the past and I prefer the "real" pubchem id's
mapping ={
   'casRegistry'   '<nobr>CAS:'
   'cheBIId' 'https://www.ebi.ac.uk/chebi/'
  %  'pubchemId' 'https://pubchem.ncbi.nlm.nih.gov/'
    'chembl'    'https://www.ebi.ac.uk/chembldb'
   'knapsack'  'http://kanaya.naist.jp/knapsack_jsp'
    };

for i = startSearch : endSearch%size(Mets,1)

    if ~isempty(metabolite_structure.(Mets{i}).VMHId) && isempty(find(isnan(metabolite_structure.(Mets{i}).VMHId),1))

        if(~isempty(metabolite_structure.(Mets{i}).keggId) && isempty(find(isnan(metabolite_structure.(Mets{i}).keggId),1)))

            try

                url = strcat('https://www.genome.jp/dbget-bin/www_bget?',metabolite_structure.(Mets{i}).keggId);
                syst = urlread(url);

                for k = 1 : size(mapping,1)
                    [metabolite_structure,idNew] = getData(metabolite_structure,syst,Mets{i}, mapping(k,:),IDsAdded);
                    if ~isempty(idNew) && isempty(find(isnan(idNew),1))
                        if isempty(metabolite_structure.(Mets{i}).(mapping{k,1})) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).(mapping{k,1})),1))
                            metabolite_structure.(Mets{i}).(mapping{k,1}) = (idNew);
                            metabolite_structure.(Mets{i}).(strcat(mapping{k,1},'_source')) = [annotationSource,':',annotationType,':',datestr(now)];
                            IDsAdded{a,1} = Mets{i};
                            IDsAdded{a,2} = mapping{k,1};
                            IDsAdded{a,3} = char(idNew);
                            a = a+1;
                        end
                    end
                end
            catch
                continue;
            end
        end
    end
end

function [metabolite_structure,idNew] = getData(metabolite_structure,syst,met, map,IDsAdded)
a = size(IDsAdded,1)+1;
idNew = '';
try
    startSearchvalue=strfind(syst,map{1,2});
    string = syst(startSearchvalue(1):startSearchvalue(1)+200);

    if contains(map{1,1},'cas')
        string = regexprep(string,'<nobr>CAS: </nobr></td><td>','');
        [tok,rem] = strtok(string,'<');
        idNew = tok;
    elseif contains(map{1,1},'pubchemId')
        [tok,rem] = strtok(string,'=');
        [tok2,rem2] = strtok(rem,'>');
        idNew = regexprep(tok2,'=','');
        idNew = regexprep(idNew,'\W$','');
    elseif contains(map{1,1},'cheBIId')
        [tok,rem] = strtok(string,'=');
        [tok2,rem2] = strtok(rem,'>');
        idNew = regexprep(tok2,'=','');
        idNew = regexprep(idNew,'CHEBI:','');
        idNew = regexprep(idNew,'\W$','');
    elseif contains(map{1,1},'chembl') ||  contains(map{1,1},'knapsack')
        [tok,rem] = strtok(string,'>');
        [tok2,rem2] = strtok(rem,'<');
        idNew = regexprep(tok2,'>','');

    end

end

author: opencobra

repo: cobratoolbox


text:
function x = touch( p, x, iseq )
global cvx___
if nargin < 3, iseq = false; end

if isa( x, 'cvx' ),
    p  = index( p );
    b  = cvx_basis( x );
    y  = any( b, 2 );
    if iseq,
    	cvx___.canslack( y ) = false;
    end
    v  = cvx___.problems( p ).t_variable;
    nv = size( v, 1 );
    ny = length( y );
    if ny < nv,
        y( nv, : ) = 0;
    elseif nv < ny,
        y = y( 1 : nv, : );
    end
    cvx___.problems( p ).t_variable = v | y;
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function term  = sdlfmjMeanCompute(sdlfmKern, t , option)

% SDLFMJMEANCOMPUTE Jolt mean for the switching dynamical LFM model.
% Computes the terms $r_d$ and $q_d$ that appear in the mean function
% associated with the switching dynamical LFM model. If the mean function
% is mu(t), then
%
%     mu(t) = r_d(t)y_d(t_0) + q_d(t)\dot{y}_d(t_0),
%
% where $y_d(t_0)$ is the initial condition associated to the position and
% $\dot{y}_d(t_0)$ is the initial condition associated to the velocity.
%
% FORMAT
% DESC
% ARG sdlfmKern : switching dynamical LFM kernel structure with the
% parameters.
% ARG t : input times for which the mean is to be computed.
% RETURN term : the value of $r_d$.
%
% FORMAT
% DESC
% Computes the terms that appear in the mean function associated with the
% switching dynamical LFM model.
% ARG sdlfmKern : switching dynamical LFM kernel structure with the
% parameters.
% ARG t : input times for which the mean is to be computed.
% ARG option : indicates which term of the mean should be computed. Option
% 'Pos' computes the term $r_d$ and option 'Vel' computes $q_d$ that
% accompanies the initial condition of the velocity.
% RETURN term : the value of $r_d$ or $q_d$ depending of option.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 3
    option = 'Pos';
end

alpha = sdlfmKern.damper/(2*sdlfmKern.mass);
omega = sqrt(sdlfmKern.spring/sdlfmKern.mass-alpha^2);
freq = omega*t;

switch option
    case 'Pos'
        term = (alpha^2/omega + omega)*exp(-alpha*t).*...
            ((omega^2 - alpha^2)*sin(freq) + 2*alpha*omega*cos(freq));
    case 'Vel'
        term = exp(-alpha*t).*((3*alpha*omega - alpha^3/omega)*sin(freq)...
            +(3*alpha^2- omega^2)*cos(freq));
    otherwise
     error('No recognized option')
end


text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)

% this is a simpler demo file for testing on your own images.

function demo_test_simple_voc()

rng('shuffle');
addpath('./my_utils');
dir_matConvNet='../libs/matconvnet/matlab';
run(fullfile(dir_matConvNet, 'vl_setupnn.m'));

run_config=[];

run_config.use_gpu=true;
% run_config.use_gpu=false;
run_config.gpu_idx=1;

% result dir:
result_name=['runner_result_dir' datestr(now, 'YYYYmmDDHHMMSS')];
result_dir=fullfile('../cache_data', 'test_examples_voc', result_name);
mkdir_notexist(result_dir);

% the folder that contains testing images:
img_data_dir='../datasets/example_imgs_voc';

% using a trained model which is trained on VOC 2012
run_config.trained_model_path='../model_trained/refinenet_res101_voc2012.mat';
run_config.class_info=gen_class_info_voc();

% for trained model, control the size of input images
run_config.input_img_short_edge_min=450;
run_config.input_img_short_edge_max=600;

runner_info=prepare_runner_test_simple(run_config);

img_filenames=my_list_file(img_data_dir);
img_num=length(img_filenames);
for img_idx=1:img_num
    task_info=[];
    task_info.img_dir=img_data_dir;
    task_info.img_filename=img_filenames{img_idx};
    task_result=runner_info.run_task_fn(runner_info, task_info);

    [~, img_name]=fileparts(task_info.img_filename);
    one_cache_file=fullfile(result_dir, [img_name '.png']);
    fprintf('save predict mask:%s\n', one_cache_file);
    imwrite(task_result.mask_data, run_config.class_info.mask_cmap, one_cache_file);
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function payload_removed_zeros = removeZeros(payload)
    for j = 1:size(payload, 2)
        if ~any(all(payload(:, j), 2))
            payload = payload(:, 1:j-1);
            break;
        end
    end
    payload_removed_zeros = payload;
end

author: UMich-BipedLab


text:
function pass = test_jacobian( )
% Test jacobian

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% Check jacobian of an empty diskfunv is an empty diskfun.
u = diskfunv;
f = jacobian(u);
pass(1) = isempty(f) & isa(f,'diskfun');

% Check definition:
F = diskfunv(@(x,y) cos(pi*x.*y), @(x,y) sin(pi*y).*cos(pi*x));
Fx = diffx(F);
Fy = diffy(F);
pass(2) = ( norm(jacobian(F) - (Fx(1).*Fy(2)-Fy(1).*Fx(2)) ) < tol );

%check against true value
true = diskfun(@(x,y)-pi^2*y.*sin(pi*x.*y).*cos(pi*y).*cos(pi.*x)-...
    pi^2*x.*sin(pi*x.*y).*sin(pi*y).*sin(pi*x));
pass(3) = ( norm(jacobian(F)-true) < 4*tol ) ;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function cub = CubatureVolumeMesh2D(CubatureOrder)

% function cub = CubatureVolumeMesh2D(CubatureOrder)
% purpose: build cubature nodes, weights and geometric factors for all elements

Globals2D;

% set up cubature nodes
[cub.R,cub.S,cub.W, cub.Ncub] = Cubature2D(CubatureOrder);

% evaluate generalized Vandermonde of Lagrange interpolant functions at cubature nodes
cub.V  = InterpMatrix2D(cub.R, cub.S);

% evaluate local derivatives of Lagrange interpolants at cubature nodes
[cub.Dr,cub.Ds] = Dmatrices2D(N,cub.R,cub.S,V);

% evaluate the geometric factors at the cubature nodes
[cub.rx,cub.sx,cub.ry,cub.sy,cub.J] = GeometricFactors2D(x,y, cub.Dr,cub.Ds);

% custom mass matrix per element
cub.mmCHOL = zeros(Np, Np, K); cub.mm = zeros(Np, Np, K);
for k=1:K
  cub.mm(:,:,k)     = cub.V'*diag(cub.J(:,k).*cub.W(:))*cub.V;
  cub.mmCHOL(:,:,k) = chol(cub.mm(:,:,k));
end

% incorporate weights and Jacobian
cub.w = cub.W; cub.W = cub.W*ones(1,K); cub.W = cub.W.*cub.J;

% compute coordinates of cubature nodes
cub.x = cub.V*x; cub.y = cub.V*y;
return

repo: Aero-matlab


text:
function [snrdB,ptotdB,psigdB,pnoisedB] = calcSNR(vout,f,fB,w,N)
% SNR calculation in the time domain (P. Malcovati, S. Brigati)
% function [snrdB,ptotdB,psigdB,pnoisedB] = calcSNR(vout,f,fB,w,N)
% vout: Sigma-Delta bit-stream taken at the modulator output
% f:    Normalized signal frequency (fs -> 1)
% fB:	Base-band frequency bins
% w:	windowing vector
% N:    samples number
%
% snrdB:     SNR in dB
% ptotdB:    Bit-stream power spectral density (vector)
% psigdB:    Extracted signal power spectral density (vector)
% pnoisedB:  Noise power spectral density (vector)
%
fB=ceil(fB);
signal=(N/sum(w))*sinusx(vout(1:N).*w,f,N);	% Extracts sinusoidal signal
noise=vout(1:N)-signal;			            % Extracts noise components
stot=((abs(fft((vout(1:N).*w)'))).^2);		% Bit-stream PSD
ssignal=(abs(fft((signal(1:N).*w)'))).^2;	% Signal PSD
snoise=(abs(fft((noise(1:N).*w)'))).^2;		% Noise PSD
pwsignal=sum(ssignal(1:fB));                % Signal power
pwnoise=sum(snoise(1:fB));		            % Noise power
snr=pwsignal/pwnoise;
snrdB=dbp(snr);
norm=sum(stot)/sum(vout(1:N).^2)*N;								% PSD normalization
if nargout > 1
	ptot=stot/norm;
	ptotdB=dbp(ptot);
end

if nargout > 2
	psig=ssignal/norm;
	psigdB=dbp(psig);
end

if nargout > 3
	pnoise=snoise/norm;
	pnoisedB=dbp(pnoise);
end


text:
% Exercise 4.57: Capacity of a communication channel
% Boyd & Vandenberghe "Convex Optimization"
% Joëlle Skaf - 04/24/08
%
% We consider a discrete memoryless communication channel, with input
% X(t) \in {1,...,n}, and output Y(t) \in {1,...,m}, for t = 1,2,...
% The relation between the input and output is given statistically:
%           p_ij = Prob(Y(t)=i|X(t)=j), i=1,...,m,  j=1,...,n
% The matrix P is called the channel transition matrix.
% The channel capacity C is given by
%           C = sup{ I(X;Y) | x >= 0, sum(x) = 1},
% I(X;Y) is the mutual information between X and Y, and it can be shown
% that:     I(X;Y) = c'*x - sum_{i=1}^m y_i*log_2(y_i)
% where     c_j = sum_{i=1}^m p_ij*log_2(p_ij), j=1,...,m

% Input data
rand('state', 0);
n = 15;
m = 10;
P = rand(m,n);
P = P./repmat(sum(P),m,1);
c = sum(P.*log2(P))';

% Channel capacity
cvx_begin
    variable x(n)
    y = P*x;
    maximize (c'*x + sum(entr(y))/log(2))
    x >= 0;
    sum(x) == 1;
cvx_end
C = cvx_optval;

% Results
display(['The channel capacity is: ' num2str(C) ' bits.'])

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function db_close

% DB_CLOSE closes the connection to the database
%
% Use as
%   db_close

% Copyright (C) 2007, Robert Oostenveld
%
% $Log: db_close.m,v $
% Revision 1.1  2007/11/07 10:50:45  roboos
% created helper functions for easy access to a MySQL database table using a structure for representing the data
%

mysql('close');

author: fieldtrip

repo: fieldtrip


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%

for i = 1:length(lambda_vec),
    lambda = lambda_vec(i);
    theta = trainLinearReg(X, y, lambda);
    error_train(i) = linearRegCostFunction(X, y, theta, 0);
    error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end;

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

end