text:
function [stats, TABLE] = SummStats(TEMP,p)
% =======================================================================
% Computes descriptive stats of a time series (or a matrix of time series)
% =======================================================================
% [stats, TABLE] = SummStats(TEMP,p)
% -----------------------------------------------------------------------
% INPUT
%   - TEMP  : T obs (rows) x N series (columns)
%   - p = order of time polynomial in the ADF null-hypothesis.
%         p ~exist, does not perform ADF test
%         p = -1, no deterministic part
%         p =  0, for constant term
%         p =  1, for constant plus time-trend
%         p >  1, for higher order polynomial
% -----------------------------------------------------------------------
% OUTPUT
%   - stats: [Obs, mean, median, max, min, StDev, AutoCorr, Skew, Kurt,...
%       ADF(4), ADF(4) Conf, ADF(8), ADF(8) Conf]
%   - TABLE: cell, matrix with the names of the stats & the stats
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015

[~, c] = size(TEMP);

% Compute summary statistics NaN are not considered (like in Excel)
mean    = nanmean(TEMP);
median  = nanmedian(TEMP);
max     = nanmax(TEMP);
min     = nanmin(TEMP);
StDev   = nanstd(TEMP);
CV      = StDev./mean;

% Compute autocorrelation coefficient AutoCorr on the maximum number of
% observation available per series
for i=1:c
    X = TEMP(:,i);
    if isnan(X)==1
        Obs(1,i)  = NaN;
        AutoCorr(1,i)  = NaN;
        Skew(1,i) = NaN;
        Kurt(1,i) = NaN;
        ADF4(1,i) = NaN;
        ADF4conf(1,i) = NaN;
        ADF8(1,i) = NaN;
        ADF8conf(1,i) = NaN;
    else
        X(any(isnan(X),2),:) = [];
        size_X    = size(X);
        Obs(1,i)  = size_X(1);
        AutoCorr(1,i)  = corr(X(2:end),X(1:end-1));
        Skew(1,i) = skewness(X);
        Kurt(1,i) = kurtosis(X);
        % If less than 20 observations don't do ADF4
        if length(X)<20 || ~exist('p','var')
            ADF4(1,i) = NaN;
            ADF4conf(1,i) = NaN;
        else
            test = adf(X,p,4);
            ADF4(1,i) = test.adf;
            if test.adf<test.crit(1)    % cannot reject the null of I(0)
                ADF4conf(1,i) = 0.01;
            elseif test.adf<test.crit(2)
                ADF4conf(1,i) = 0.05;
            elseif test.adf<test.crit(3)
                ADF4conf(1,i) = 0.10;
            else                        % reject the null of I(0) => I(1)
                ADF4conf(1,i) = 0;
            end
        end
        % If less than 40 observations don't do ADF8
        if length(X)<40 || ~exist('p','var')
            ADF8(1,i) = NaN;
            ADF8conf(1,i) = NaN;
        else
            test = adf(X,p,8);
            ADF8(1,i) = test.adf;
            if test.adf<test.crit(1)
                ADF8conf(1,i) = 0.01;
            elseif test.adf<test.crit(2)
                ADF8conf(1,i) = 0.05;
            elseif test.adf<test.crit(3)
                ADF8conf(1,i) = 0.10;
            else
                ADF8conf(1,i) = 0;
            end
        end
    end
end

% Create matrix with all summary stats
stats = [Obs; mean; median; max; min; StDev; AutoCorr; Skew; Kurt; ADF4; ADF4conf; ADF8; ADF8conf; CV];

% Create titles
titles ={'Obs'; 'Mean'; 'Median'; 'Max'; 'Min'; 'St. Dev.'; 'Auto Corr.'; 'Skew.'; 'Kurt.'; 'ADF(4)'; 'ADF(4) Conf.'; 'ADF(8)'; 'ADF(8) Conf.'; 'Coeff. Variation'};

% Create cell
TABLE = TabPrint(stats,[],titles);

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [data,units] = compute_min_wing_area(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = min(trx(fly).wing_areal_mm,trx(fly).wing_arear_mm);

end
units = parseunits('mm^2');

author: kristinbranson

repo: JAABA


text:
% homodyne_recon_demo.m
%
% - Demonstrates use of homodyne_recon.m, which implements the homodyne
% reconstruction method outlined in Doug Noll's 1991 IEEE T-MI paper
% "Homodyne detection in magnetic resonance imaging"
% - Figure shows the decrease in RMSE that results from sampling more.
% - Reconstructed image can be rectangular, with odd or even dimensions.
%
% 2012-06-15, Mai Le
% 2012-09-26 Jeff Fessler tweaked

% set up "true" brain image
%f.dir = [path_find_dir('mri') '/../data/mri/'];
%f.xtrue = [f.dir 'brainweb_t1.jpg'];
%mag_true = double(imread(f.xtrue))'; % true image magnitude
mag_true = 25*ellipse_im(256);
mag_true = mag_true(16:end-16,3:end-2); % make it rectangle to stress test
dims = size(mag_true);

% simulate linear phase over object
[xx yy] = ndgrid(-dims(1)/2:dims(1)/2-1,-dims(2)/2:dims(2)/2-1);
ph_max = pi/2;
ph = ph_max*xx/(dims(1)/2) + ph_max*yy/(dims(2)/2);
image = mag_true .* exp(1i*ph); % modulate to make it complex
im_fft = fftshift(fft2(image));

im plc 2 3
im(1, mag_true)

for homodyne_direction = [1 2] % partial k-space can be in direction 1 or 2
	% # of samples in "half" kspace:
	nhalf = ceil((dims(homodyne_direction)+1)/2);
	overlaps = dims(homodyne_direction) - nhalf; % maximum # of extra rows
	overlaps = [0:5:(overlaps-1) overlaps];
	nover = numel(overlaps);

	nrms_im = zeros(nover,1);
	nrms_ph_demod = zeros(nover,1);

	for ii = 1:nover
		overlap = overlaps(ii);
		if (homodyne_direction == 1)
			k1 = 1:min(dims(1), nhalf+overlap);
			partial_kspace = im_fft(k1,:);
		elseif (homodyne_direction == 2)
			k2 = 1:min(dims(2), nhalf+overlap);
			partial_kspace = im_fft(:,k2);
		else
			fail 'bug'
		end

		f.sigma = 800; % complex AWGN std dev
		noisy_partial_kspace = partial_kspace + f.sigma * ...
			(randn(size(partial_kspace)) + 1i*randn(size(partial_kspace)));

		[recon_ph_demod, recon_im, lp_im, full_kspace] = ...
			homodyne_recon(noisy_partial_kspace, dims(1), dims(2), ...
				overlap, 'direction', homodyne_direction);

		nrms_fun = @(x,y) 100 * nrms(x(:), y(:));
		nrms_im(ii) = nrms_fun(abs(recon_im), mag_true);
		nrms_ph_demod(ii) = nrms_fun(abs(recon_ph_demod), mag_true);

		if overlap <= 15
			clim = minmax(mag_true)';
			im(2, abs(recon_im), clim)
			titlef('Conventional, overlap=%d', overlap)
			im(3, abs(recon_ph_demod), clim)
			titlef('Demodulated, overlap=%d', overlap)
			im(4, log(abs(partial_kspace)), clim)
			im(5, abs(recon_im) - mag_true, [-9 9], 'Error')
			im(6, abs(recon_ph_demod) - mag_true, [-9 9], 'Error')
			drawnow
	%		prompt
		end
	end

	f.snr = sqrt(mean(abs(partial_kspace(:)).^2)) / ...
		sqrt(mean(abs(noisy_partial_kspace(:) - partial_kspace(:)).^2));
	f.snr = 20*log10(f.snr);
	pr f.snr

	im subplot 4
	plot(overlaps, nrms_im, 'g-+', overlaps, nrms_ph_demod, 'b-o')
	axis tight
	xtick([0 15 max(overlaps)])
	title('NRMSE (%)')
	legend('without demodulation', 'phase demodulated');
	xlabel('rows acquired beyond 1/2 of kspace');
	ylabel('NRMS error of reconstructed image (%)');
prompt
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [] = mc_test_pdfs(MM,VV,desired_steps)
% Function tests and plots the Gaussian mixture in the desired simulation
% steps.
%
%% Syntax
% mc_test_pdfs(MM,VV,desired_steps)
%
%% Description
% Function tests and plots the Gaussian mixture described by the vector of Gaussian
% means and variances in the desired simulation steps.
%
% Input:
% * MM            ... the matrix of mean values of Gaussian components, the matrix of dimensions
%                       ksteps x Nsamples
% * VV            ... the matrix of associated variances of Gaussian components
% * desired_steps ... the vector with indices of steps, where we want to test
%                       the distributions
%
% See Also:
% demo_example_gp_simulation
%
% Examples:
% demo_example_gp_simulation
%

%%
% * Written by K. Azman, 2007.
%

% variables
MM = MM';
VV = VV';

Nsamples = size(MM,1);
mu = mean(MM);
s2 = mean(VV) + mean((MM-repmat(mu,Nsamples,1)).^2);
mu = mu';
s2 = s2';

% "calculation"
STD = sqrt(VV);
Xwide = 200; % resolution = no. of samples between xmin and xmax
Xmin = min(MM-4*STD);
Xmax = max(MM+4*STD);
dX = (Xmax-Xmin)/(Xwide-1);

% XX,YY ... matrices for storing the GP model's pdfs
% XX(:,k) ... yvalues at time step k
% YY(:,k) ... density values at time step k and coresponding y values in XX(:,k)
% sum(YY(:,k)) = 1 for all k-s
XX = zeros(Xwide,size(VV,2));
YY = zeros(size(XX));

% for desired steps
for jj = 1:length(desired_steps)
    kk = desired_steps(jj);
    XX(:,kk) = Xmin(kk):dX(kk):Xmax(kk);
    % ... and for every realisation
    ytemp = zeros(size(YY,1),1);

    for ii=1:Nsamples
        ytemp = normpdf(XX(:,kk), MM(ii,kk), STD(ii,kk));
        YY(:,kk) = YY(:,kk) + ytemp;
    end
    YY(:,kk) = YY(:,kk)/Nsamples;

    yapprox = normpdf(XX(:,kk),mu(kk),sqrt(s2(kk)));

    % plot pdfs
    figure(10000+kk);

    plot(XX(:,kk),YY(:,kk),'Color',[0.6 0.6 0.6],'LineWidth',4);
    hold on;
    plot(XX(:,kk),yapprox,'LineStyle','--','Color',[0 0 1],'LineWidth',4);
     for ii=1:Nsamples
        ytemp = normpdf(XX(:,kk), MM(ii,kk), STD(ii,kk));
        plot(XX(:,kk),ytemp,'Color',[0.9 0.9 0.9]);
    end
    plot(XX(:,kk),YY(:,kk),'Color',[0.6  0.6 0.6],'LineWidth',4);
    plot(XX(:,kk),yapprox,'LineStyle','--','Color',[0 0 1],'LineWidth',4);
    hold off;
    grid;
    legend('true','GP approx','samples')

    title(strcat(['test MC pdfs, simul step=', num2str(kk)]));

    disp('press sth to continue');
    pause;

end

author: Dynamic-Systems-and-GP

repo: GPdyn


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

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

% Last Modified by GUIDE v2.5 23-Dec-2020 21:34:33
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GINavPlot_OpeningFcn, ...
                   'gui_OutputFcn',  @GINavPlot_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

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

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

H=get(0,'ScreenSize'); h=get(gcf,'Position');
x=H(3)/2-h(3)/2; y=H(4)/2-h(4)/2; h(1)=x; h(2)=y;
set(gcf,'Position',h);

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

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes on button press in solfile.
function solfile_Callback(hObject, eventdata, handles)
% hObject    handle to solfile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global glc
fullpath=which('GINavExe.m');
[path,~,~]=fileparts(fullpath);
result_path=[path,glc.sep,'result',glc.sep];
[filename,pathname] = uigetfile('*.pos','Select solution file',result_path);
if any(filename)
    fullname = [pathname filename];
    if ~exist(fullname,'file')
        error('The solution file does not exist!!!');
    end
    idx=find(filename=='.');
    if any(idx)
        suffix=filename(idx+1:end);
    else
        error('Incorrect file suffix name!!!');
    end
    switch suffix
        case 'mat'
            load(fullname);
        case 'pos'
            solution=readGINavsol(fullname);
        case 'pva'
            solution=readGINavsol(fullname);
    end
    handles.solution=solution;
    guidata(hObject,handles)
else
    fprintf('Please reselect solution file!!!\n');
end

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

% Hints: contents = cellstr(get(hObject,'String')) returns plot_any contents as cell array
%        contents{get(hObject,'Value')} returns selected item from plot_any
solution=handles.solution;
str = get(handles.plot_any,'String');
val = get(handles.plot_any,'Value');
switch str{val}
    case 'Trajectory'
        plot_trajectory(solution);
    case 'Position'
        plot_position(solution);
    case 'Velocity'
        plot_velocity(solution);
    case 'Number of satellite'
        plot_nsat(solution);
    case 'AR ratio factor'
        plot_ARratio(solution);
end

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

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

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

% Hints: get(hObject,'String') returns contents of staname as text
%        str2double(get(hObject,'String')) returns contents of staname as a double
station_name=get(handles.staname,'String');
if any(station_name)
    handles.station_name=station_name;
    guidata(hObject,handles);
else
    fprintf('Please re-enter the station name!!!\n');
end

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

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

% --- Executes on button press in SNX.
function SNX_Callback(hObject, eventdata, handles)
% hObject    handle to SNX (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename,pathname] = uigetfile({'*.snx';'*.SNX'},'Select SNX file');
if any(filename)
    snxfilename = [pathname filename];
    if ~exist(snxfilename,'file')
        error('The SNX file does not exist!!!');
    end
    handles.snxfilename=snxfilename;
    guidata(hObject,handles);
else
   fprintf('Please reselect SNX file!!!\n');
end

% --- Executes on button press in plotPPP.
function plotPPP_Callback(hObject, eventdata, handles)
% hObject    handle to plotPPP (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
station_name=handles.station_name;
snxfilename=handles.snxfilename;
solution=handles.solution;
plot_ppp_err(solution,station_name,snxfilename);

% --- Executes on button press in ref.
function ref_Callback(hObject, eventdata, handles)
% hObject    handle to ref (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename,pathname] = uigetfile('*.mat','Select pos reference data');
if any(filename)
    fullname = [pathname filename];
    if ~exist(fullname,'file')
        error('The pos reference file does not exist!!!');
    end
    load(fullname);
    handles.reference=reference;
    guidata(hObject,handles)
else
    fprintf('Please reselect pos ref file!!!\n');
end

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

% Hints: contents = cellstr(get(hObject,'String')) returns plot_err contents as cell array
%        contents{get(hObject,'Value')} returns selected item from plot_err
solution  = handles.solution;
reference = handles.reference;
str = get(handles.plot_err,'String');
val = get(handles.plot_err,'Value');
switch str{val}
    case 'Position error'
        plot_pva_err(solution,reference,1)
    case 'Velocity error'
        plot_pva_err(solution,reference,2)
    case 'Attitude error'
        plot_pva_err(solution,reference,3)
end

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

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

repo: GINav


text:
%% -- mat2PDB.m --
%
% this function creates a PDB from coordinate data. Represent all inputs as
% a structure field for it to be read. The output format is as given in
% online documentation (as of July 2012 when writing this program)
% http://www.wwpdb.org/documentation/format33/sect9.html#ATOM
%
% Make sure all data input is one-dimensional with the same length. If
% they are not the same length, the program ignores user input, states
% an error, and uses defaults. All string inputs should be in cell-format.
% Keep in mind the "element" and "charge" inputs are strings in
% cell-format, not numbers.
%
%
% -- required inputs (3) --
%
% input value        meaning
%
% input.X            orthagonal X coordinate data (angstroms)
% input.Y            orthagonal Y coordinate data (angstroms)
% input.Z            orthagonal Z coordinate data (angstroms)
%
% -- optional inputs (12): generates defaults when not user-specified --
%
% input value        meaning                           default value
%
% input.outfile      output file name                 "mat2PDB.pdb"
% input.recordName   output record name of atoms      "ATOM"
% input.atomNum      atom serial number                sequential number
% input.atomName     name of atoms                    "OW" (water oxygen)
% input.altLoc       alt. location indicator          " "
% input.resName      name of residue                  "SOL" (water)
%
% input.chainID      protein chain identifier         "A"
% input.resNum       residue sequence number           sequential number
% input.occupancy    occupancy factor                 "1.00"
% input.betaFactor   beta factor, temperature         "0.00"
% input.element      element symbol                   "O" (oxygen)
% input.charge       atomic charge                    " "
%
%
% -- example uses --
%
% % translates both X and Y coordinates of 3IJU.pdb by 5 angstroms
% PDBdata = pdb2mat('3IJU.pdb');
% PDBdata.X = PDBdata.X + 5;
% PDBdata.Y = PDBdata.Y + 5;
% mat2pdb(PDBdata)
%
% % make a PDB with 30 atoms in random places within a 10 angstrom box
% data.X = rand(1,20)*10;
% data.Y = rand(1,20)*10;
% data.Z = rand(1,20)*10;
% mat2pdb(data)
%
%

function mat2pdb(input)
%% review XYZ coordinate data

% coordinate data is required! Checking XYZ input
if ~isfield(input, 'X') || ~isfield(input, 'Y') || ~isfield(input, 'Z')
    fprintf('we need xyz coordinate data to make a PDB!!\n\texiting...\n');
    return;
end
X = input.X;
Y = input.Y;
Z = input.Z;
if length(X) ~= length(Y) || length(X) ~= length(Z)
    fprintf('xyz coordinate data is not of equal lengths!\n\texiting...\n');
    return;
end

%% review optional data inputs

% in case optional data data not given, fill in blanks
if ~isfield(input, 'outfile')
    input.outfile = 'mat2PDB.pdb';
end
if ~isfield(input, 'recordName')
    input.recordName = cell(1,length(X));
    input.recordName(1:end) = {'ATOM'};
end
if ~isfield(input, 'atomNum')
    input.atomNum = 1:length(X);
end
if ~isfield(input, 'atomName')
    input.atomName = cell(1,length(X));
    input.atomName(1:end) = {'OW'};
end
if ~isfield(input, 'altLoc')
    input.altLoc = cell(1,length(X));
    input.altLoc(1:end) = {' '};
end
if ~isfield(input, 'resName')
    input.resName = cell(1,length(X));
    input.resName(1:end) = {'SOL'};
end
if ~isfield(input, 'chainID')
    input.chainID = cell(1,length(X));
    input.chainID(1:end) = {'A'};
end
if ~isfield(input, 'resNum')
    input.resNum = 1:length(X);
end
if ~isfield(input, 'occupancy')
    input.occupancy = ones(1,length(X));
end
if ~isfield(input, 'betaFactor')
    input.betaFactor = zeros(1, length(X));
end
if ~isfield(input, 'element')
    input.element = cell(1,length(X));
    input.element(1:end) = {'O'};
end
if ~isfield(input, 'charge')
    input.charge = cell(1,length(X));
    input.charge(1:end) = {' '};
end

outfile    = input.outfile;
recordName = input.recordName;
atomNum    = input.atomNum;
atomName   = input.atomName;
altLoc     = input.altLoc;
resName    = input.resName;
chainID    = input.chainID;
resNum     = input.resNum;
occupancy  = input.occupancy;
betaFactor = input.betaFactor;
element    = input.element;
charge     = input.charge;

%% remove faulty inputs

if length(recordName) ~= length(X)
    fprintf('recordName input is not the correct length!\n\tignoring user input\n');
    recordName = cell(1,length(X));
    recordName(1:end) = {'ATOM'};
end
if length(atomNum) ~= length(X)
    fprintf('atom serial number input is not the correct length!\n\tignoring user input\n');
    atomNum = 1:length(X);
end
if length(atomName) ~= length(X)
    fprintf('atom name input is not the correct length!\n\tignoring user input\n');
    atomName = cell(1,length(X));
    atomName(1:end) = {'OW'};
end
if length(altLoc) ~= length(X)
    fprintf('alternate location input is not the correct length!\n\tignoring user input\n');
    altLoc = cell(1,length(X));
    altLoc(1:end) = {' '};
end
if length(resName) ~= length(X)
    fprintf('residue name input is not the correct length!\n\tignoring user input\n');
    resName = cell(1,length(X));
    resName(1:end) = {'SOL'};
end
if length(chainID) ~= length(X)
    fprintf('chain ID input is not the correct length!\n\tignoring user input\n');
    chainID = cell(1,length(X));
    chainID(1:end) = {'A'};
end
if length(resNum) ~= length(X)
    fprintf('residue number input is not the correct length!\n\tignoring user input\n');
    resNum = 1:length(X);
end
if length(occupancy) ~= length(X)
    fprintf('occupancy input is not the correct length!\n\tignoring user input\n');
    occupancy = ones(1,length(X));
end
if length(betaFactor) ~= length(X)
    fprintf('beta factor input is not the correct length!\n\tignoring user input\n');
    betaFactor = zeros(1, length(X));
end
if length(element) ~= length(X)
    fprintf('element symbol input is not the correct length!\n\tignoring user input\n');
    element = cell(1,length(X));
    element(1:end) = {'O'};
end
if length(charge) ~= length(X)
    fprintf('charge input is not the correct length!\n\tignoring user input\n');
    charge = cell(1,length(X));
    charge(1:end) = {' '};
end

% fix atomName spacing
for n = 1:length(atomName)
    atomName(n) = {sprintf('%-3s',cell2mat(atomName(n)))};
end

%% create PDB

% open file
fprintf('outputting PDB in file %s\n', outfile);
FILE = fopen(outfile, 'w');

% output data
for n = 1:length(atomNum)

    % standard PDB output line
    fprintf( FILE, '%-6s%5u%5s%1.1s%3s %1.1s%4u%12.3f%8.3f%8.3f%6.2f%6.2f%12s%2s\n', ...
        cell2mat(recordName(n)), atomNum(n), cell2mat(atomName(n)), ...
        cell2mat(altLoc(n)), cell2mat(resName(n)), cell2mat(chainID(n)), ...
        resNum(n), X(n), Y(n), Z(n), occupancy(n), betaFactor(n), ...
        cell2mat(element(n)), cell2mat(charge(n)));

    % output progress in terminal
    if ~mod(n,400)
        fprintf('   %6.2f%%', 100*n / length(atomNum));
        if ~mod(n, 4000)
            fprintf('\n');
        end
    end

end
fprintf( FILE, 'END\n');

% close file
fprintf('   %6.2f%%\n    done! closing file...\n', 100);

fclose(FILE);

end

author: opencobra

repo: cobratoolbox


text:
function res = badchannels(this, varargin)
% Method for getting/setting bad channels
% FORMAT res = badchannels(this)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: badchannels.m 5933 2014-03-28 13:22:28Z vladimir $

if length(varargin) == 2 && isnumeric(varargin{1})
    % make sure that the two inputs for set are the same length
    if ~(length(varargin{2}) == 1 || (length(varargin{1}) == length(varargin{2})))
        error('Use either same vector length or scalar for value');
    end
end

if numel(varargin) >= 1  && (isnumeric(varargin{1})  && ~isempty(varargin{1}))
    if ~(all(varargin{1} >= 1) && all(varargin{1} <= nchannels(this)))
        error('Channel number out of range.');
    end
end

if numel(varargin) >= 2 && (isnumeric(varargin{1})  && ~isempty(varargin{1}))
    ubad = unique(varargin{2});
    if isempty(ubad) || ~all(ismember(ubad, [0 1]))
        error('Illegal bad flags (should be 0 or 1)');
    end
end

if this.montage.Mind == 0
    res = getset(this, 'channels', 'bad', varargin{:});
elseif numel(varargin) >= 2
    this.montage.M(this.montage.Mind) = getset(this.montage.M(this.montage.Mind), 'channels', 'bad', varargin{:});
    res = this;
else
    res = getset(this.montage.M(this.montage.Mind), 'channels', 'bad', varargin{:});
end

% Return channel indices if called without arguments and [0, 1] if called
if numel(varargin) <= 1 % get
    if iscell(res)
        res = [res{:}];
    end
    res = logical(res);
    if isempty(varargin)
        res = find(res);
    end
end

author: spm


text:
function [imuData, frameRange] = loadImuData(imuDataFolder, imageTimestamps)
% loadImageData Read captured image data into memory. Tries to load
% individual images or a saved mat file if one exists.

    %Read IMU data
    oxtsData = loadOxtsliteData(imuDataFolder);
    % load IMU Data
    v_index = 9:11; % 12:14 body frame, 15:17 FLU frame
    a_index = 15:17; % 15:17 FLU frame
    omega_index = 21:23; % 18:20 body frame, 21:23 FLU frame

    dateStrings = loadTimestamps([imuDataFolder '/oxts']);
    timestamps = zeros(1, length(dateStrings));
    for i = 1:length(dateStrings)
        timestamps(i) =  datenum_to_unixtime(datenum(dateStrings(i))) + 0.00001;
    end

    frameRange = 1:length(timestamps);
    frameNum = length(frameRange);
    imuData.timestamps = zeros(1, frameNum);
    imuData.measVel = zeros(3, frameNum);
    imuData.measOrient = zeros(4, frameNum);
    imuData.measOmega = zeros(3,frameNum);
    imuData.measAccel = zeros(3, frameNum);
    imuData.initialVelocity = zeros(3,1);

    % for all oxts packets do
    for meas_i=1:length(frameRange)

      i = frameRange(meas_i);
      % if there is no data
      if isempty(oxtsData{i})
        continue;
      end

     imuData.timestamps(1,meas_i) = timestamps(i);
     imuData.measOmega(:,meas_i) = oxtsData{i}(omega_index);

      rx = oxtsData{i}(4); % roll
      ry = oxtsData{i}(5); % pitch
      rz = oxtsData{i}(6); % heading
      Rx = [1 0 0; 0 cos(rx) -sin(rx); 0 sin(rx) cos(rx)]; % base => nav  (level oxts => rotated oxts)
      Ry = [cos(ry) 0 sin(ry); 0 1 0; -sin(ry) 0 cos(ry)]; % base => nav  (level oxts => rotated oxts)
      Rz = [cos(rz) -sin(rz) 0; sin(rz) cos(rz) 0; 0 0 1]; % base => nav  (level oxts => rotated oxts)
      R  = Rz*Ry*Rx;

      imuData.measOrient(:,meas_i) = rotMatToQuat(R);
      imuData.measVel(:,meas_i) =  oxtsData{i}(v_index);
      imuData.measAccel(:,meas_i) =  oxtsData{i}(a_index)' - R'*[0;0;9.81];

      if meas_i == 1
          imuData.initialVelocity = getRnb(oxtsData{i})'*[ oxtsData{i}(8); oxtsData{i}(7); 0; ];
      end
    end

    function dn = datenum_to_unixtime( date_num )
      dn =  (date_num - 719529)*86400;         %# 719529 == datenum(1970,1,1)
    end

end

repo: MSCKF


text:
function f = prolong(f, nOut)
%PROLONG   Manually adjust the number of points used in a TRIGTECH.
%   G = PROLONG(F, N) returns a TRIGTECH G where LENGTH(G) = N and G represents
%   the same function as F but using more or less coefficients than F.
%
%   If N < LENGTH(F) the representation is compressed by chopping
%   coefficients, which may result in a loss of accuracy.
%
%   If N > LENGTH(F) the coefficients are padded with zeros.
%
% See also ALIAS.

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

% Get the number of coefficients.
n = length(f);

% Return if nOut == n
if ( nOut == n )
    return
end

% Get coefficients
coeffs = f.coeffs;

% If n is even extend to coeffs to n+1
if ( mod(n,2) == 0 )
    coeffs = [.5*coeffs(1,:);coeffs(2:end,:);.5*coeffs(1,:)];
    n = n + 1;
end

% Return if nOut == n
if ( nOut == n )
    f.coeffs = coeffs;
    f.values = f.coeffs2vals(f.coeffs);
    f.values(:,f.isReal) = real(f.values(:,f.isReal));
    return
end

% Pad with zeros when nOut > n:
if ( nOut > n )
    kup = ceil((nOut-n)/2);
    kdown = floor((nOut-n)/2);
    coeffs = [zeros(kup, size(coeffs, 2)); coeffs; zeros(kdown, size(coeffs,2))];
    f.coeffs = coeffs;
    f.values = f.coeffs2vals(f.coeffs);
    f.values(:,f.isReal) = real(f.values(:,f.isReal));
    return
end

% Chop coefficients when nOut < n:
if ( nOut < n )
    kup = floor((n-nOut)/2);
    kdown = ceil((n-nOut)/2);
    coeffs(end-kdown+1:end,:) = [];
    coeffs(1:kup,:) = [];
    if ( kup < kdown )
        coeffs(1,:) = 2*coeffs(1,:); % scale coefficients in even case
    end
    f.coeffs = coeffs;
    f.values = f.coeffs2vals(f.coeffs);
    f.values(:,f.isReal) = real(f.values(:,f.isReal));
    return
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [f1,f2,DialedNum] = DecodedFrequencies(E);

%Computes the energy of the Signal
%Input Variables
%    E    = Energy Matrix of the Filtered Signals
%
% Output Variables
%    f1   = first decoded frequency
%    f2   = second decoded frequency

Fmat = [697 770 852 941 1209 1336 1477 1633];

LowEnergyMatrix = E(1:4);
HighEnergyMatrix = E(5:8);

[e1,e1Index]=max(LowEnergyMatrix);
[e2,e2Index]=max(HighEnergyMatrix);

f1 = Fmat(e1Index);
f2 = Fmat(4+e2Index);

if        f1 == Fmat(1) && f2 == Fmat(5)
   DialedNum = '1';
   elseif f1 == Fmat(1) && f2 == Fmat(6)
   DialedNum = '2';
   elseif f1 == Fmat(1) && f2 == Fmat(7)
   DialedNum = '3';
   elseif f1 == Fmat(1) && f2 == Fmat(8)
   DialedNum = 'A';
   elseif f1 == Fmat(2) && f2 == Fmat(5)
   DialedNum = '4';
   elseif f1 == Fmat(2) && f2 == Fmat(6)
   DialedNum = '5';
   elseif f1 == Fmat(2) && f2 == Fmat(7)
   DialedNum = '6';
   elseif f1 == Fmat(2) && f2 == Fmat(8)
   DialedNum = 'B';
   elseif f1 == Fmat(3) && f2 == Fmat(5)
   DialedNum = '7';
   elseif f1 == Fmat(3) && f2 == Fmat(6)
   DialedNum = '8';
   elseif f1 == Fmat(3) && f2 == Fmat(7)
   DialedNum = '9';
   elseif f1 == Fmat(3) && f2 == Fmat(8)
   DialedNum = 'C';
   elseif f1 == Fmat(4) && f2 == Fmat(5)
   DialedNum = '*';
   elseif f1 == Fmat(4) && f2 == Fmat(6)
   DialedNum = '0';
   elseif f1 == Fmat(4) && f2 == Fmat(7)
   DialedNum = '#';
   elseif f1 == Fmat(4) && f2 == Fmat(8)
   DialedNum = 'D';
else
    DialedNum = 'p';
end

author: Sable

repo: mcbench-benchmarks


text:
function [ ft ] = gsp_time_translate(G, f,tau,param)
%GSP_TIME_TRANSLATE Generalized time-vertex translation of the signal f to the node i and time t
%   Usage: ft = gsp_time_translate(G, f,tau,param);
%
%   Input parameters
%       G   : Time-Vertex Graph structure
%       f   : Time-Vertex signal
%       tau : Time location
%   Output parameters
%       ft  : Translated signal
%
%   This function translate the time-vertex signal *f* at time tau.
%
%   Additional parameters
%   ---------------------
%   * *param.boundary*  : Time boundary condition for the translation:
%                                                'periodic' (fft),  (default)
%                                                'reflecting' (dct),
%                                                'absorbing' (zero-padding).
%
%
% Author :  Francesco Grassi

if nargin<5
    param=struct;
end

if ~isfield(param,'boundary'), param.boundary='periodic'; end

[N,T] = size(f);

error('Change the boundary!')

switch param.boundary
    case 'periodic'
        delta = gsp_delta(G.jtv.T,tau).';
        fhat = fft(f,[],2);
        operator = repmat(fft(delta,[],2),N,1);
        ft = sqrt(G.N)*ifft(fhat .* operator,[],2);

    case 'symmetric'
        delta = gsp_delta(G.jtv.T,tau).';
        fhat = dct(f.').';
        operator = repmat(dct(delta.').',N,1);
        ft = sqrt(G.N)*idct( (fhat .* operator).' ).';

    case 'absorbing'

        ind = round(-tau*G.jtv.fs);
        ft = circshift(f,[0 ind]);

        if tau >0

            ft = [ft(:,1:end+ind) zeros(N,-ind)];

        else

            ft = [zeros(N,ind-1) ft(:,ind+1:end)];

        end

    otherwise
        error('Unknown boundary condition');

end

end

repo: gspbox


text:
% DEMOILVARGPLVM4 Run variational GPLVM on oil data.

% VARGPLVM

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

dataSetName = 'oil';
experimentNo = 4;
printDiagram = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = vargplvmOptions('dtcvar');
%options.kern = {'rbfard2', 'bias', 'white'};
options.kern = 'rbfardjit';
options.numActive = 50;
%options.tieParam = 'tied';

options.optimiser = 'scg';
latentDim = 10;
d = size(Y, 2);

% demo using the variational inference method for the gplvm model
model = vargplvmCreate(latentDim, d, Y, options);
%
model = vargplvmParamInit(model, model.m, model.X);
model.vardist.covars = 0.5*ones(size(model.vardist.covars)) + 0.001*randn(size(model.vardist.covars));
model.learnBeta=1;

% Optimise the model.
iters = 1800;
display = 1;

model.beta = 1/((1/100 * var(model.m(:))));
model.learnBeta = false; model.learnSigmaf = false; model.initVardist = true;
model = vargplvmOptimise(model, display, 500);
model.learnBeta = true; model.learnSigmaf = true; model.initVardist = false;

model = vargplvmOptimise(model, display, iters);

capName = dataSetName;
capName(1) = upper(capName(1));
modelType = model.type;
modelType(1) = upper(modelType(1));
save(['dem' capName modelType num2str(experimentNo) '.mat'], 'model');

% order wrt to the inputScales
mm = vargplvmReduceModel(model,2);
%% plot the two largest twe latent dimensions
if exist('printDiagram') & printDiagram
  lvmPrintPlot(mm, lbls, capName, experimentNo);
  bar(model.kern.inputScales);
end
errors = fgplvmNearestNeighbour(mm, lbls);
fprintf('# Vargplvm errors in the 2-D projection: %d\n', errors)

author: SheffieldML

repo: vargplvm


text:
function r = gamrnd(a,b,m,n);
%GAMRND Random matrices from gamma distribution.
%   R = GAMRND(A,B) returns a matrix of random numbers chosen
%   from the gamma distribution with parameters A and B.
%   The size of R is the common size of A and B if both are matrices.
%   If either parameter is a scalar, the size of R is the size of the other
%   parameter. Alternatively, R = GAMRND(A,B,M,N) returns an M by N matrix.
%
%   Some references refer to the gamma distribution
%   with a single parameter. This corresponds to GAMRND
%   with B = 1. (See Devroye, pages 401-402.)

%   GAMRND uses a rejection or an inversion method depending on the
%   value of A.

%   References:
%      [1]  L. Devroye, "Non-Uniform Random Variate Generation",
%      Springer-Verlag, 1986

%   B.A. Jones 2-1-93
%   Copyright (c) 1993-98 by The MathWorks, Inc.
%   $Revision: 1.1.1.1 $  $Date: 2005/04/26 02:29:18 $

if nargin < 2,
   error('Requires at least two input arguments.');
end

if nargin == 2
   [errorcode rows columns] = rndcheck(2,2,a,b);
end

if nargin == 3
   [errorcode rows columns] = rndcheck(3,2,a,b,m);
end

if nargin == 4
   [errorcode rows columns] = rndcheck(4,2,a,b,m,n);
end

if errorcode > 0
   error('Size information is inconsistent.');
end

% Initialize r to zero.
lth = rows*columns;
r = zeros(lth,1);
a = a(:); b = b(:);

scalara = (length(a) == 1);
if scalara
   a = a*ones(lth,1);
end

scalarb = (length(b) == 1);
if scalarb
   b = b*ones(lth,1);
end

% If a == 1, then gamma is exponential. (Devroye, page 405).
k = find(a == 1);
if any(k)
   r(k) = -b(k) .* log(rand(size(k)));
end

k = find(a < 1 & a > 0);
% (Devroye, page 418 Johnk's generator)
if any(k)
  c = zeros(lth,1);
  d = zeros(lth,1);
  c(k) = 1 ./ a(k);
  d(k) = 1 ./ (1 - a(k));
  accept = k;
  while ~isempty(accept)
    u = rand(size(accept));
    v = rand(size(accept));
    x = u .^ c(accept);
    y = v .^ d(accept);
    k1 = find((x + y) <= 1);
    if ~isempty(k1)
      e = -log(rand(size(k1)));
      r(accept(k1)) = e .* x(k1) ./ (x(k1) + y(k1));
      accept(k1) = [];
    end
  end
  r(k) = r(k) .* b(k);
end

% Use a rejection method for a > 1.
k = find(a > 1);
% (Devroye, page 410 Best's algorithm)
bb = zeros(size(a));
c  = bb;
if any(k)
  bb(k) = a(k) - 1;
  c(k) = 3 * a(k) - 3/4;
  accept = k;
  count = 1;
  while ~isempty(accept)
    m = length(accept);
    u = rand(m,1);
    v = rand(m,1);
    w = u .* (1 - u);
    y = sqrt(c(accept) ./ w) .* (u - 0.5);
    x = bb(accept) + y;
    k1 = find(x >= 0);
    if ~isempty(k1)
      z = 64 * (w .^ 3) .* (v .^ 2);
      k2 = (z(k1) <= (1 - 2 * (y(k1) .^2) ./ x(k1)));
      k3 = k1(find(k2));
      r(accept(k3)) = x(k3);
      k4 = k1(find(~k2));
      k5 = k4(find(log(z(k4)) <= (2*(bb(accept(k4)).*log(x(k4)./bb(accept(k4)))-y(k4)))));
      r(accept(k5)) = x(k5);
      omit = [k3; k5];
      accept(omit) = [];
    end
  end
  r(k) = r(k) .* b(k);
end

% Return NaN if a or b is not positive.
r(b <= 0 | a <= 0) = NaN;

r = reshape(r,rows,columns);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [fx] = f_metaToM(x,P,u,inF)
% evolution function for meta-learning (k-ToM vs sequence)
% function [fx] = f_metaToM(x,P,u,inF)
% This function update the belief of an agent who has to identify whether
% he faces an intentional agent (k-ToM) or an inanimate patter (sequence).
% Learning here derives from VB-Laplace update rules of the agent's
% posterior belief. practically speaking, this function (i) wraps modified
% evolution functions of both k-ToM learning and sequence learning, and
% (ii) updates the probability Pi that the agent is facing an intentional
% agent given a new observation u (the probability that the agent faces an
% inanimate pattern is 1-Pi).
% IN:
%   - x: hidden states (see indexing in inF.indlev)
%   - P: evolution params:
%   P(1)= agent's prior opponent's (log-) volatility on opponent's params
%   P(2)= agent's (invsigmoid-) dilution coefficient
%   - u: sequence of past actions:
%       u(1)= opponent's last move
%       u(2)= learner's last move
%       u(3:K+2) = sequence of K past opponent's moves
%   - inF: input structure (see prepare_metaToM.m)
% OUT:
%   - fx: updated hidden states
% [see RecToMfunction.m and f_BSL.m]

ot = u(1); % opponent's last move
fx = NaN(size(x)); % initialize updated states

% 1- update P(agent=kToM)
Pi0 = VBA_sigmoid(x(inF.meta.indx)); % prior P(agent=1)
% partial forgetting of prior belief on opponent's type?
if inF.meta.diluteP
    dc = VBA_sigmoid(P(inF.meta.indP)); % dilution coefficient
    Pi0 = (1-dc).*Pi0 + dc./2;
end

% Get k-ToM's likelihood, ie P(o|k-ToM).
xktom = x(inF.ktom.indx);
inFktom = inF.ktom.inF;
ntotPar = inFktom.indParev+inFktom.indParobs; % total number of params
level = inFktom.lev; % depth of k-ToM's recursive beliefs (k=level)
indlev = defIndlev(level,ntotPar); % states indexing
if level==0 % 0-ToM [should be useless]
    mx = xktom(1); % E[log-odds of P(o=1)]
    Vx = exp(xktom(2)); % V[log-odds of P(o=1)]
    Els1 = VBA_Elogsig(mx,Vx);
    Els0 = VBA_Elogsig(-mx,Vx);
    ELL = ot.*Els1 + (1-ot).*Els0;
    h_kToM = exp(ELL); % P(o|k-ToM)
else
    Pk = VBA_sigmoid(x(1:(level-1))); % P(k'), with k'=0,...,k-1
    Pk = [Pk;max(0,1-sum(Pk))]; % insert last P(k'=k-1)
    f = zeros(level,1); % E[x(theta)]
    Vx = zeros(level,1); % V[x(theta)]
    for j=1:level % loop over possible opponent's levels  (k'=j-1)
        f(j) = xktom(indlev(j).f); % E[x(theta)|k'=j-1]
        df = xktom(indlev(j).df); % d[x(theta)]/dtheta for k'=j-1
        Sig = exp(xktom(indlev(j).Par(2:2:2*ntotPar))); % V[theta|k'=j-1]
        Vx(j) = sum(Sig.*df.^2); % V[x(theta)|k'=j-1]
    end
    Els1 = VBA_Elogsig(f,Vx);
    Els0 = VBA_Elogsig(-f,Vx);
    ELL = ot.*Els1 + (1-ot).*Els0;
    h_kToM = exp(Pk'*ELL); % P(o|k-ToM)
end

% Get sequence's likelihood, ie P(o|seq).
xseq = x(inF.seq.indx);
inFseq = inF.seq.inF;
useq = u; % remove agent's previous move
useq(2) = [];
K = inFseq.K; % sequence depth
yb = useq(2:K+1); % previous outcomes
if VBA_isWeird (yb)
    h_seq = 1/2;
else
    if K >0
        indSeq = bin2dec(num2str(yb'))+1; % index of sequence of previous outcomes
    else
        indSeq = 1;
    end
    m = xseq(indSeq);
    v = exp(xseq((2^K)+indSeq));
    Els1 = VBA_Elogsig(m,v);
    Els0 = VBA_Elogsig(-m,v);
    ELL = ot.*Els1 + (1-ot).*Els0;
    h_seq = exp(ELL); % P(o|seq)
end

% VB update of P(agent=kToM)
Pi = Pi0.*h_kToM./(Pi0.*h_kToM+(1-Pi0).*h_seq);
fx(inF.meta.indx) = VBA_sigmoid(Pi, 'inverse', true);

% 2- update k-ToM belief
inFktom.metaweight = Pi;
Par_ktom = P(inF.ktom.indP);
uktom = u(1:2);
fx(inF.ktom.indx) = RecToMfunction(xktom,Par_ktom,uktom,inFktom);

% 3- update seq learner belief
inFseq.metaweight = 1-Pi;
Par_seq = P(inF.seq.indP);
fx(inF.seq.indx) = f_BSL(xseq,Par_seq,useq,inFseq);

author: MBB-team

repo: VBA-toolbox


text:
function b = ismonom(a)
% function B = ismonom(A)
%
% DESCRIPTION
%   Returns true for a monomial or a vector or matrix of monomials.
%
% INPUTS
%   A: polynomial
%
% OUTPUTS
%   B: 1 if A is a monomial or a vector or matrix of monomials.  Returns
%      0 otherwise.
%
% SYNTAX
%   B = ismonom(A);

% 1/29/2008: PJS  Initial Coding

if isa(a,'double') && a==1
    b = true;
    return;
elseif ~isa(a,'polynomial')
    b = false;
    return;
end

% acoef will be Nterms-by-(Nrows*Ncols)
acoef = a.coefficient;

if all(nonzeros(acoef)==1) && all(sum(acoef,1)==1)
    b = true;
else
    b = false;
end

author: yu-jiang

repo: radpbook

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


text:
function varargout = goldsteinprice(X)
% Goldstein-Price function
%
%   GOLDSTEINPRICE([x1, x2]) returns the value of the value of the
%   Goldstein-Price function at the specified points. [x1] and [x2]
%   may be vectors. The search domain is
%
%               -2 < x_i < 2
%
%   The golbal minimum is
%
%               f(x1, x2) = f(0, -1) = 3.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;  % # dims
        varargout{2} = [-2, -2]; % LB
        varargout{3} = [+2, +2]; % UB
        varargout{4} = [0, -1]; % solution
        varargout{5} = 3; % function value at solution

    % otherwise, output function value
    else

        % keep values in the search domain
        X(X < -2) = inf;  X(X > 2) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = (1  + (x1 + x2 + 1).^2.*(19 - 14*x1 +  3*x1.^2 - 14*x2 +  6*x1.*x2 +  3*x2.^2)).*...
            (30 + (2*x1 - 3*x2).^2.*(18 - 32*x1 + 12*x1.^2 + 48*x2 - 36*x1.*x2 + 27*x2.^2));

    end

end


text:
function ft_sourcewrite(cfg, source)

% FT_SOURCEWRITE exports source-reconstructed results to gifti or nifti format file.
% The appropriate output file depends on whether the source locations are described by
% on a cortically constrained sheet (gifti) or by a regular 3D lattice (nifti).
%
% Use as
%  ft_sourcewrite(cfg, source)
% where source is a source structure obtained from FT_SOURCEANALYSIS and
% cfg is a structure that should contain
%
%   cfg.filename  = string, filename without the extension
%   cfg.filetype  = string, can be 'nifti', 'gifti' or 'cifti' (default is automatic)
%   cfg.parameter = string, functional parameter to be written to file
%   cfg.precision = string, can be 'single', 'double', etc.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this the input data will be read from a *.mat
% file on disk. This mat file should contain only a single variable,
% corresponding with the input data structure.
%
% See also FT_SOURCEANALYSIS, FT_SOURCEDESCRIPTIVES, FT_VOLUMEWRITE

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar source
ft_preamble provenance source

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% ensure that the required options are present
cfg = ft_checkconfig(cfg, 'required', {'parameter', 'filename'});

% get the options
cfg.parameter      = ft_getopt(cfg, 'parameter');
cfg.filename       = ft_getopt(cfg, 'filename');
cfg.filetype       = ft_getopt(cfg, 'filetype');        % the default is determined further down
cfg.brainstructure = ft_getopt(cfg, 'brainstructure');  % is used for cifti
cfg.parcellation   = ft_getopt(cfg, 'parcellation');    % is used for cifti
cfg.precision      = ft_getopt(cfg, 'precision');       % is used for cifti

% check if the input data is valid for this function
if strcmp(cfg.filetype, 'cifti')

  % keep the transformation matrix
  if isfield(source, 'transform')
    transform = source.transform;
  elseif isfield(source, 'brainordinate') && isfield(source.brainordinate, 'transform')
    transform = source.brainordinate.transform;
  else
    transform = [];
  end

  if isfield(source, 'brainordinate')
    % it is a parcellated source representation, i.e. the main structure one channel for each parcel
    brainordinate = source.brainordinate;
    source = rmfield(source, 'brainordinate');

    % split them and check individually
    source        = ft_checkdata(source, 'datatype', {'timelock', 'freq', 'chan'}, 'feedback', 'yes');
    brainordinate = ft_checkdata(brainordinate, 'datatype', 'parcellation', 'parcellationstyle', 'indexed', 'hasunit', 'yes');

    % merge them again
    source = copyfields(brainordinate, source, setdiff(fieldnames(brainordinate), {'cfg'}));
  else
    source = ft_checkdata(source, 'datatype', 'source', 'hasunit', true, 'feedback', 'yes');
  end

  % keep the transformation matrix
  if ~isempty(transform)
    source.transform = transform;
  end

end % if cifti

if isempty(cfg.filetype)
  if isfield(source, 'dim')
    % source positions are on a regular 3D lattice, save as nifti
    cfg.filetype = 'nifti';
  elseif isfield(source, 'tri')
    % there is a specification of a 2D cortical sheet, save as gifti
    cfg.filetype = 'gifti';
  else
    ft_error('the input data does not look like a 2D sheet, nor as a 3D regular volume');
  end
end

switch (cfg.filetype)
  case 'nifti'
    if numel(cfg.filename)<=4 || ~strcmp(cfg.filename(end-3:end), '.nii');
      cfg.filename = cat(2, cfg.filename, '.nii');
    end

    % convert functional data into 4D
    dat = getsubfield(source, cfg.parameter);
    dat = reshape(dat, source.dim(1), source.dim(2), source.dim(3), []);

    if ~isfield(source, 'transform')
      source.transform = pos2transform(source.pos);
    end

    % this is a bit of a kludge, but ft_write_mri can write 3D and 4D nifti
    ft_write_mri(cfg.filename, dat, 'dataformat', 'nifti', 'transform', source.transform);

  case 'gifti'
    if numel(cfg.filename)<=4 || ~strcmp(cfg.filename(end-3:end), '.gii');
      cfg.filename = cat(2, cfg.filename, '.gii');
    end

    % this is a bit of a kludge, but ft_write_headshape can write gifti including data
    ft_write_headshape(cfg.filename, source, 'data', getsubfield(source, cfg.parameter), 'format', 'gifti');

  case 'cifti'
    % brainstructure should represent the global anatomical structure, such as CortexLeft, Thalamus, etc.
    % parcellation should represent the detailled parcellation, such as BA1, BA2, BA3, etc.
    ft_write_cifti(cfg.filename, source, 'parameter', cfg.parameter, 'brainstructure', cfg.brainstructure, 'parcellation', cfg.parcellation, 'precision', cfg.precision);

  otherwise
    ft_error('unsupported output format "%s"', cfg.filetype);
end % switch filetype

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous source
ft_postamble provenance

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function u = ldivide(a,b)
%LDIVIDE      Slope elementwise left division  a .\ b
%

% written  11/02/05     S.M. Rump
%

  u = b ./ a;

repo: OpenMAS


text:
function [ value, en ] = wew_a ( x )

%*****************************************************************************80
%
%% WEW_A estimates Lambert's W function.
%
%  Discussion:
%
%    For a given X, this routine estimates the solution W of Lambert's
%    equation:
%
%      X = W * EXP ( W )
%
%    This routine has higher accuracy than WEW_B.
%
%  Modified:
%
%    11 June 2014
%
%  Reference:
%
%    Fred Fritsch, R Shafer, W Crowley,
%    Algorithm 443: Solution of the transcendental equation w e^w = x,
%    Communications of the ACM,
%    October 1973, Volume 16, Number 2, pages 123-124.
%
%  Parameters:
%
%    Input, real X, the argument of W(X)
%
%    Output, real VALUE, the estimated value of W(X).
%
%    Output, real EN, the last relative correction to W(X).
%
  c1 = 4.0 / 3.0;
  c2 = 7.0 / 3.0;
  c3 = 5.0 / 6.0;
  c4 = 2.0 / 3.0;
%
%  Initial guess.
%
  f = log ( x );

  if ( x <= 6.46 )

    wn = x .* ( 1.0 + c1 * x ) ./ ( 1.0 + x .* ( c2 + c3 * x ) );
    zn = f - wn - log ( wn );

  else

    wn = f;
    zn = - log ( wn );

  end
%
%  Iteration 1.
%
  temp = 1.0 + wn;
  y = 2.0 * temp .* ( temp + c4 * zn ) - zn;
  wn = wn .* ( 1.0 + zn .* y ./ ( temp .* ( y - zn ) ) );
%
%  Iteration 2.
%
  zn = f - wn - log ( wn );
  temp = 1.0 + wn;
  temp2 = temp + c4 * zn;
  en = zn .* temp2 ./ ( temp .* temp2 - 0.5 * zn );
  wn = wn .* ( 1.0 + en );

  value = wn;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ value, x ] = randlc_jump ( x, k )

%*****************************************************************************80
%
%% RANDLC_JUMP returns the K-th element of a uniform pseudorandom sequence.
%
%  Discussion:
%
%    The sequence uses the linear congruential generator:
%
%      X(K+1) = A * X(K)  mod 2^46
%
%    The K-th element, which can be represented as
%
%      X(K) = A^K * X(0)  mod 2^46
%
%    is computed directly using the binary algorithm for exponentiation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Bailey, Eric Barszcz, John Barton, D Browning, Robert Carter,
%    Leonardo Dagum, Rod Fatoohi,
%    Samuel Fineberg, Paul Frederickson, Thomas Lasinski, Robert Schreiber,
%    Horst Simon, V Venkatakrishnan, Sisira Weeratunga,
%    The NAS Parallel Benchmarks,
%    RNR Technical Report RNR-94-007,
%    March 1994.
%
%    Donald Knuth,
%    The Art of Computer Programming,
%    Volume 2, Seminumerical Algorithms,
%    Third Edition,
%    Addison Wesley, 1997,
%    ISBN: 0201896842,
%    LC: QA76.6.K64.
%
%  Parameters:
%
%    Input, uint64 X, the initial seed (with index 0).
%
%    Input, integer K, the index of the desired value.
%
%    Output, real VALUE, the K-th value in the sequence.
%
%    Output, uint64 X, the K-th seed.
%
  persistent a;
  persistent a1;
  persistent a2;
  persistent ks;
  persistent r23;
  persistent r46;
  persistent t23;
  persistent t46;
%
%  If this is the first call, compute
%
%    R23 = 2 ^ -23,
%    R46 = 2 ^ -46,
%    T23 = 2 ^ 23,
%    T46 = 2 ^ 46.
%
%  These are computed in loops, rather than by merely using the power operator,
%  in order to insure that the results are exact on all systems.
%
  if ( isempty ( ks ) )

    r23 = 1.0;
    r46 = 1.0;
    t23 = 1.0;
    t46 = 1.0;

    for i = 1 : 23
      r23 = 0.5 * r23;
      t23 = 2.0 * t23;
    end

    for i = 1 : 46
      r46 = 0.5 * r46;
      t46 = 2.0 * t46;
    end

    a = 1220703125.0;
%
%  Break A into two parts such that A = 2^23 * A1 + A2.
%
    t1 = r23 * a;
    a1 = floor ( t1 );
    a2 = a - t23 * a1;

    ks = 1;

  end

  if ( k < 0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'RANDLC_JUMP - Fatal error!\n' );
    fprintf ( 1, '  K < 0.\n' );
    error ( 'RANDLC_JUMP - Fatal error!' )

  elseif ( k == 0 )

%
%  Find M so that K < 2^M.
%
  else

    m = 1;
    twom = 2;
    while ( twom <= k )
      twom = twom * 2;
      m = m + 1;
    end

    b = a;
    b1 = a1;
    b2 = a2;

    for i = 1 : m

      j = floor ( k / 2 );
%
%  Replace X by A * X, if appropriate.
%
      if ( 2 * j ~= k )

        t1 = r23 * x;
        x1 = floor ( t1 );
        x2 = x - t23 * x1;

        t1 = b1 * x2 + b2 * x1;
        t2 = floor ( r23 * t1 );
        z = t1 - t23 * t2;

        t3 = t23 * z + b2 * x2;
        t4 = floor ( r46 * t3 );
        x = t3 - t46 * t4;

      end
%
%  Replace A by A * A mod 2^46.
%
      t1 = r23 * b;
      x1 = floor ( t1 );
      x2 = b - t23 * x1;

      t1 = b1 * x2 + b2 * x1;
      t2 = floor ( r23 * t1 );
      z = t1 - t23 * t2;

      t3 = t23 * z + b2 * x2;
      t4 = floor ( r46 * t3 );
      b = t3 - t46 * t4;
%
%  Update A1, A2.
%
      t1 = r23 * b;
      b1 = floor ( t1 );
      b2 = b - t23 * b1;

      k = j;

    end

  end

  value = r46 * x;

  return
end


text:
function exampleFilter
% function exampleFilter
%
% Example of how the Kalman filter performs for real data

load raw

%Apply filter
k=Kalman_Stack_Filter(raw);
k75=Kalman_Stack_Filter(raw,0.75);

%Set up image window
minMax=[min(raw(:)),max(raw(:))];
clf, colormap gray

subplot(1,3,1)
imagesc(raw(:,:,1))
title('original')
set(gca,'clim',minMax), axis off equal

subplot(1,3,2)
imagesc(k(:,:,1))
title('filtered, gain=0.5')
set(gca,'clim',minMax), axis off equal

subplot(1,3,3)
imagesc(k75(:,:,1))
title('filtered, gain=0.75')
set(gca,'clim',minMax), axis off equal

%Loop movie
disp('crtl-c to stop movie')
while 1
    for i=1:size(k,3)

        subplot(1,3,1)
        set(get(gca,'children'),'CData',raw(:,:,i))

        subplot(1,3,2)
        set(get(gca,'children'),'CData',k(:,:,i))

        subplot(1,3,3)
        set(get(gca,'children'),'CData',k75(:,:,i))

        pause(0.05)
        drawnow
    end
end


text:
function toggletoolbar(varargin)
%TOGGLETOOLBAR Toggle distribution fit plot toolbar on or off

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

if (nargin>0 & ishandle(varargin{1}) & ...
               isequal(get(varargin{1},'Type'),'figure'))
   dffig = varargin{1};
else
   dffig = gcbf;
end

tbstate = get(dffig,'toolbar');
h = findall(dffig,'Type','uitoolbar');
if isequal(tbstate,'none') | isempty(h)
   % Create toolbar for the first time
   set(dffig,'toolbar','figure');
   dfadjusttoolbar(dffig);
elseif nargin>1 & isequal(varargin{2},'on')
   % Hide toolbar
   set(h,'Visible','on');
else
   % Show toolbar
   set(h,'Visible','off');
end

author: zouchuhang

repo: LayoutNet


text:
function [zgMax, bgMax, bgOriented] = getLocalGradients(imName)
	paths = getPaths();
	typ = {'bgOriented', 'bgMax', 'zgMax'};
	for t = 1:3,
		fileName = fullfile(paths.gradientsForRecognition, typ{t}, strcat(imName, '.mat'));
		dt = load(fileName);
		eval(sprintf('%s = dt.signal;', typ{t}));
	end
	ng1 = zgMax(:,:,[1:4 4]);
	ng2 = zgMax(:,:,[5:8 8]);
	dg = zgMax(:,:,[9:12 12]);
	zgMax = cat(3, ng1, ng2, dg);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function p20_default ( dim_num )

%*****************************************************************************80
%
%% P20_DEFAULT sets default values for problem 20.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p20_i4 ( 'D', '*', i4 );

  r8 = 0.0;
  p20_r8 ( 'D', '*', r8 );

  return
end


text:
function varargout = setxor(varargin)
% Adds support for 'stable' input flag with one output
% Also, z = vertcat(z, y); is used instead of z = [z; y]; because of an
% Octave bug when vertically concatenating a 1x0 char array with another
% char array; see https://savannah.gnu.org/bugs/index.php?52542. It seems
% to be avoided using `vertcat` instead of `;`
if iscell(varargin{1}) && ~iscell(varargin{2}), varargin{2} = {varargin{2}}; end
if iscell(varargin{2}) && ~iscell(varargin{1}), varargin{1} = {varargin{1}}; end
if nargin>=3 && strcmp(varargin{end},'stable')
    % This is almost like setdiff:
    if strcmp(varargin{3},'rows')
        y = [varargin{1}; varargin{2}];
        [a, ~, x] = unique(y,'rows');
        x1 = x(1:size(varargin{1},1));
        x2 = x(size(varargin{1},1)+1:end);
    else
        y = [varargin{1}(:); varargin{2}(:)];
        [a, ~, x] = unique(y);
        x1 = x(1:numel(varargin{1}));
        x2 = x(numel(varargin{1})+1:end);
    end
    ind = ~ismember(x1, x2);
    x = x1(ind);
    if ~isempty(x)
        x = x(~any(triu(bsxfun(@eq, x, x.'),1)));
    end
    if strcmp(varargin{3},'rows')
        y = a(x,:);
    else
        y = a(x);
        y = y(:);
    end
    % Keep partial output:
    z = y;
    % Same with inputs reversed:
    if strcmp(varargin{3},'rows')
        y = [varargin{2}; varargin{1}];
        [a, ~, x] = unique(y,'rows');
        x2 = x(1:size(varargin{2},1));
        x1 = x(size(varargin{2},1)+1:end);
    else
        y = [varargin{2}(:); varargin{1}(:)];
        [a, ~, x] = unique(y);
        x2 = x(1:numel(varargin{2}));
        x1 = x(numel(varargin{2})+1:end);
    end
    ind = ~ismember(x2, x1);
    x = x2(ind);
    if ~isempty(x)
        x = x(~any(triu(bsxfun(@eq, x, x.'),1)));
    end
    if strcmp(varargin{3},'rows')
        y = a(x,:);
    else
        y = a(x);
        y = y(:);
    end
    % Join partial outputs:
    z = vertcat(z, y);
    if ~strcmp(varargin{3},'rows') && isrow(varargin{1}) && isrow(varargin{2})
        z = z.';
    end
    varargout{1} = z;
else
    if strcmp(varargin{end},'sorted'), varargin(end) = []; end
    varargout = cell(1,max(nargout,1)); % if called without outputs, nargout will be zero. In that case we return one output
    [varargout{:}] = builtin('setxor', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


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

    properties
        stage LaunchVehicleStage

        stageStateToSet(1,1) logical = true;

        inheritStageState(1,1) logical = true;
        inheritStageStateFrom(1,1) InheritStateEnum = InheritStateEnum.InheritFromLastState;
        inheritStageStateFromEvent LaunchVehicleEvent = LaunchVehicleEvent.empty(1,0);
    end

    methods
        function obj = SetKinematicStateStageState(stage)
            obj.stage = stage;
        end

        function updateStateLogEntry(obj, newStateLogEntry, stateLog)
            stageState = obj.getStageState(newStateLogEntry);

            if(obj.inheritStageState)
                if(obj.inheritStageStateFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
                   not(isempty(obj.inheritStageStateFromEvent)) && ...
                   not(isempty(obj.inheritStageStateFromEvent.getEventNum())) && ...
                   obj.inheritStageStateFromEvent.getEventNum() > 0)

                    evtStateLog = stateLog.getLastStateLogForEvent(obj.inheritStageStateFromEvent);
                    if(not(isempty(evtStateLog)))
                        evtStageState = obj.getStageState(evtStateLog);
                        stageState.active = evtStageState.active;
                    end
                end
            else
                stageState.active = obj.stageStateToSet;
            end
        end

        function listboxStr = getListboxStr(obj)
            listboxStr = {};
            for(i=1:length(obj))
                listboxStr{i} = obj(i).stage.name; %#ok<AGROW>
            end
        end
    end

    methods(Access=private)
        function stageState = getStageState(obj, stateLogEntry)
            stageStates = stateLogEntry.stageStates;
            stageState = stageStates([stageStates.stage] == obj.stage);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function h = ldivide( f, g )
%.\   Pointwise SEPARABLEAPPROX left array divide.
%   F.\G if G is a SEPARABLEAPPROX and F is a double this returns (1/F)*G
%
%   F.\G if G is a double and F is a SEPARABLEAPPROX this returns G\F, but this does
%   not work if F becomes numerically close to zero.
%
%   F.\G we do not allow F and G to both be SEPARABLEAPPROX objects.
%
%   F.\G is the same as the command ldivide(F,G)

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

h = rdivide( g, f );

end


text:
function op = prox_l1pos( q )
%PROX_L1POS    L1 norm, restricted to x >= 0
%    OP = PROX_L1( q ) implements the nonsmooth function
%        OP(X) = norm(q.*X,1) + indicator_{ X >= 0 }
%    Q is optional; if omitted, Q=1 is assumed. But if Q is supplied,
%    then it must be a positive real scalar (or must be same size as X).

% New in v1.0d

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) ||  any( q(:) < 0 ) || all(q(:)==0) %|| numel( q ) ~= 1
	error( 'Argument must be positive.' );
end

op = tfocs_prox( @(x)f(x,q), @(x,t)prox_f(x,t,q), 'vector');
end

function v = f(x,q)
    if any( x(:) < 0 )
        v = Inf;
    elseif isscalar(q)
        v = q*sum( x(:) );
    else
        v = sum( q(:).*x(:) );
    end
end

% The proximity operator is a simplified version of shrinkage:
function x = prox_f(x,t,q)
    x   = max( 0, x - t*q );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
% Jason Joseph Rebello
% Carnegie Mellon University (Jan 2012 - May 2013)
% MS in Electrical & Computer Engineering
% K Means Algorithm with the application to image compression

% This program uses the K means clustering algorithm to group the pixels
% in an image in order to provide image compression.

clear all;
clc;
close all;

fprintf('K means clustering algorithm used for image compression\n\n');
t = cputime;

%% Read the image
fprintf('Reading image');
I = imread('bird_small.png');
%imshow(I);
I = (double(I))/255;
fprintf('...done\n\n');

%% Declare and Initialize Variabels
fprintf('Initializing variables');
K = 16; % number of clusters
imgSize = size(I); % get size of image
iterCentroids = 10; % number of times K means runs to find the best centroid
iterKMeans = 10; % number of times K means runs with different initial centroids
fprintf('...done\n\n');

%% Get input
fprintf('Formatting input');
X = reshape(I, imgSize(1) * imgSize(2), 3); % resize into (total pixel x features)
fprintf('...done\n\n');

%% Run K Means
for i=1:iterKMeans

    fprintf(' ********* Running K means iteration %d ***********\n\n',i);
    [centroids cost idx] = runKMeans(X, K, iterCentroids);
    fprintf('Cost after %d iteration : %f\n\n',i,cost);

    if i==1
        bestCentroids = centroids;
        bestCost = cost;
        bestidx = idx;
    elseif (i>1 && cost<bestCost) % stores the best clustering
        bestCentroids = centroids;
        bestCost = cost;
        bestidx = idx;
    end
    fprintf('Best cost : %f\n\n',bestCost);

XCompressed = centroids(idx,:);

% Reshape the recovered image into proper dimensions
XCompressed = reshape(XCompressed, imgSize(1), imgSize(2), 3);
imshow(XCompressed); % display final compressed image for each iteration
pause(1);

end

%% Display original and best compressed image
XCompressed = bestCentroids(bestidx,:);

% Reshape the recovered image into proper dimensions
XCompressed = reshape(XCompressed, imgSize(1), imgSize(2), 3);
displayImage(I, XCompressed, K);

fprintf('Program executed in %f seconds or %f minutes\n\n', cputime-t, (cputime-t)/60);


text:
function [f,g] = funObj_mfvi_exact(v, y, X, lambda)
% Primal objective function (Omega = inv(Sigma)) for mean-field
% f = 0.5*(logdet(V) - logdet(Sigma) - tr(V*SigmaInv) - (m-mu)'*SigmaInv*(m-mu)
%     + L) - sum_d fb(mbar_d, vbar_d)
% where mbar = X*m, vbar = diag(X*V*X')
%
% Written by Emtiyaz,
% Modified by Wu Lin
  [D L] = size(X);

  Omega = diag(lambda);
  %Extract mean, Cholesky and bias
  m = v(1:L);
  sigma = v(L+1:end);
  U = diag(sigma);

  % compute V
  V = diag(sigma.^2);

  % compute kl and its gradient
  kl = 0.5*(2*sum(log(diag(U))) + sum(log(lambda(2:end))) - trace(V*Omega) -m'*(Omega*m) + L);

  % contribution from the bound
  mbar = X*m;
  vbar = sum(X.*(V*X')',2); % diag(X*V*X') efficient
  [fb, gmb, gvb] = E_log_p('bernoulli_logit', y, mbar, vbar, []);
  fb = -fb;
  gm_lvb = X'*(-gmb);
  gV_lvb = X'*bsxfun(@times, -gvb, X); %efficient X'*diag(gllp/2)*X;

  % final
  f = kl - sum(fb);
  gm = - Omega*m - gm_lvb;
  gU = diag(1./diag(U)) - diag(U*Omega) - diag(2*U*gV_lvb);

  g=[gm(:); diag(gU)];

  % return
  g=-g;
  f=-f;

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function fit = lf_censor(x,y,cens,varargin)
%
% Censored local regression using normal assumption.
% Must provide x, y and cens.
% All other arguments to locfit() can be provided, with the
% exception of weights.
%
% NEED: Kaplan Meier Estimate. Iterations are fixed.
%

lfc_y = y;
unc = find(~cens);

for i = 0:3
  fit = locfit(x,lfc_y,varargin{:});
  fh = fitted(fit);

  rs = rsum(fit);
  df0 = rs(1);
  df1 = rs(2);

  rdf = sum(1-cens) - 2*df0 + df1;
  sigma = sqrt(sum( (y-fh).*(lfc_y-fh) / rdf));
  sr = (y-fh)/sigma;
  lfc_y = fh + sigma*normpdf(sr)./normcdf(-sr);
  lfc_y(unc) = y(unc);
end;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function a = maxij ( m, n )

%*****************************************************************************80
%
%% MAXIJ returns the MAXIJ matrix.
%
%  Discussion:
%
%    This matrix is occasionally known as the "Boothroyd MAX" matrix.
%
%  Formula:
%
%    A(I,J) = max(I,J)
%
%  Example:
%
%    N = 5
%
%    1 2 3 4 5
%    2 2 3 4 5
%    3 3 3 4 5
%    4 4 4 4 5
%    5 5 5 5 5
%
%  Properties:
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    The inverse of A is tridiagonal.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.13,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 42,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns
%    of the matrix.
%
%    Output, real A(M,N), the matrix.
%
  for i = 1 : m
    for j = 1 : n
      a(i,j) = max ( i, j );
    end
  end

  return
end


text:
% Compiles all mex routines that are part of toolbox.
%
% USAGE
%  toolboxCompile
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also
%
% Piotr's Image&Video Toolbox      Version 3.01
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

disp('Compiling.......................................');
rd=fileparts(mfilename('fullpath')); rd=rd(1:end-9);

% general compile options (can make architecture specific)
opts = {'-output'};
if(exist('OCTAVE_VERSION','builtin')), opts={'-o'}; end

% if you get warnings on linux, you can set the gcc version using:
% opts = {'CXX=g++-4.1' 'CC=g++-4.1' 'LD=g++-4.1' '-l' ...
%   'mwlapack' '-l' 'mwblas' '-output' };

% compile c functions
fs={'assignToBins1','histc2c','ktHistcRgb_c','ktComputeW_c',...
  'nlfiltersep_max','nlfiltersep_sum','imtransform2_c',...
  'meanShift1','fernsInds1'};
ds=[repmat({'images'},1,7),repmat({'classify'},1,2)];
for i=1:length(fs), fs{i}=[rd '/' ds{i} '/private/' fs{i}]; end
for i=1:length(fs), mex([fs{i} '.c'],opts{:},[fs{i} '.' mexext]); end

% compile c++ functions
try
  fs={'convConst','gradientMex','rgbConvertMex','imResampleMex',...
    'imPadMex','opticalFlowHsMex','forestFindThr','forestInds'};
  ds=[repmat({'channels'},1,5) 'images' repmat({'classify'},1,2)];
  for i=1:length(fs), fs{i}=[rd '/' ds{i} '/private/' fs{i}]; end
  for i=1:length(fs), mex([fs{i} '.cpp'],opts{:},[fs{i} '.' mexext]); end
catch ME
  fprintf(['C++ mex failed, likely due to lack of a C++ compiler.\n' ...
    'Run ''mex -setup'' to specify a C++ compiler if available.\n'...
    'Or, one can specify a specific C++ explicitly (see mex help).\n']);
end

disp('..................................Done Compiling');

author: kristinbranson

repo: JAABA


text:
function tests = test_yale_faces
  tests = functiontests(localfunctions);
end

function setupOnce(testCase)
    close all;
end

function teardownOnce(testCase)
    close all;
end

function test_basic(testCase)
    yf = spx.data.image.YaleFaces();
    yf.load();
    verifyEqual(testCase, yf.ImageHeight, 192);
    verifyEqual(testCase, yf.ImageWidth, 168);
    verifyEqual(testCase, yf.num_subjects, 38);
    verifyEqual(testCase, yf.num_total_images, 64*38);
    verifyEqual(testCase, yf.image_size, 192*168);
    yf.describe;
end

function test1(testCase)
    yf = spx.data.image.YaleFaces();
    yf.load();
    yf.describe();
    mf = spx.graphics.Figures();
    mf.new_figure('Example faces');
    canvas = yf.create_random_canvas();
    imshow(canvas);
    colormap(gray);
    axis image;
    axis off;

    yf.resize_all(42, 48);
    mf.new_figure('Subject 1');
    canvas = yf.create_subject_canvas(1);
    imshow(canvas);
    colormap(gray);
    axis image;
    axis off;
end

author: indigits

repo: sparse-plex


text:
% IMU information (todo: move this to the yaml file)
update_rate = 400;

% xsens mti-G-700 (page 32 and 33)
% https://www.xsens.com/download/usermanual/MTi_usermanual.pdf
gyroscope_noise_density = 0.01; %datasheet => deg/s/sqrt(Hz)
accelerometer_noise_density = 60; %datasheet => µg/sqrt(Hz)

% Convert to kalibr format
gyroscope_noise_density = pi/180*gyroscope_noise_density; %convert to => rad/s/sqrt(Hz)
accelerometer_noise_density = 9.81*accelerometer_noise_density/1000000; %convert to => m/s^2/sqrt(Hz)

%% Output results for yaml file
fprintf('\n\n=================================\n');
fprintf('rostopic: /imu0\n');
fprintf('update_rate: %.2f\n',update_rate);
fprintf('#Accelerometers\n');
fprintf('accelerometer_noise_density: %.5f.\n',accelerometer_noise_density);
fprintf('accelerometer_random_walk: %.5f.\n',0.0);
fprintf('#Gyroscopes\n');
fprintf('gyroscope_noise_density: %.6f.\n',gyroscope_noise_density);
fprintf('gyroscope_random_walk: %.6f.\n',0.0);
fprintf('=================================\n');

author: rpng


text:
function str = encodevarname(str,varargin)
%
%    newname = encodevarname(name)
%
%    Encode an invalid variable name using a hex-format for bi-directional
%    conversions.

%    This function is sensitive to the default charset
%    settings in MATLAB, please call feature('DefaultCharacterSet','utf8')
%    to set the encoding to UTF-8 before calling this function.
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        name: a string, can be either a valid or invalid variable name
%
%    output:
%        newname: a valid variable name by converting the leading non-ascii
%              letter into "x0xHH_" and non-ascii letters into "_0xHH_"
%              format, where HH is the ascii (or Unicode) value of the
%              character.
%
%              if the encoded variable name CAN NOT be longer than 63, i.e.
%              the maximum variable name specified by namelengthmax, and
%              one uses the output of this function as a struct or variable
%              name, the name will be trucated at 63. Please consider using
%              the name as a containers.Map key, which does not have such
%              limit.
%
%    example:
%        encodevarname('_a')   % returns x0x5F_a
%        encodevarname('a_')   % returns a_ as it is a valid variable name
%        encodevarname('变量')  % returns 'x0xE58F98__0xE9878F_'
%
%    this file is part of EasyH5 Toolbox: https://github.com/fangq/easyh5
%
%    License: GPLv3 or 3-clause BSD license, see https://github.com/fangq/easyh5 for details
%

    if(~isvarname(str(1)))
        if(exist('unicode2native','builtin'))
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',char(str(1))+0,str(2:end));
        end
    end
    if(isvarname(str))
        return;
    end
    if(exist('unicode2native','builtin'))
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        cpos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(cpos))
            return;
        end
        str0=str;
        pos0=[0 cpos(:)' length(str)];
        str='';
        for i=1:length(cpos)
            str=[str str0(pos0(i)+1:cpos(i)-1) sprintf('_0x%X_',str0(cpos(i))+0)];
        end
        if(cpos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end
end

author: brainstorm-tools


text:
function depth =denoise(depth,camera_D)

XYZcamera(:,:,1)=camera_D.X .* depth;
XYZcamera(:,:,2)=camera_D.Y .* depth;
XYZcamera(:,:,3)=depth .* (~isnan(camera_D.X) & ~isnan(camera_D.Y));

%figure;imagesc(XYZcamera(:,:,3)); axis equal; axis tight;

X = reshape(XYZcamera,[],3);

validID = find(X(:,3)~=0);

X = X(validID,:);

%points2ply('original.ply', X);

% denoise

kdOBJ = KDTreeSearcher(X);

[~,mindist] = knnsearch(kdOBJ,X,'K',4);

threshold1 = 0.1;
threshold2 = 0.2;
threshold3 = 0.25;
threaholdFurthest = 8; %15;
threaholdClosest = 0.5;

goodones = mindist(:,2)<threshold1 & mindist(:,3)<threshold2 & mindist(:,4)<threshold3 & X(:,3)<threaholdFurthest & X(:,3)>threaholdClosest ;

XYZcamera( 2*camera_D.width*camera_D.height + validID(~goodones)) = 0;
X = X(goodones,:);

%figure;imagesc(XYZcamera(:,:,3)); axis equal; axis tight;

%points2ply('filter1.ply', X);

depth3Max = ordfilt2(XYZcamera(:,:,3),23,true(5));

outliers = XYZcamera(:,:,3) > (depth3Max + 0.5) & (depth3Max>0);
XYZcamera( 2*camera_D.width*camera_D.height + find(outliers)) = 0;

%figure;imagesc(XYZcamera(:,:,3)); axis equal; axis tight;

%depthMed = medfilt2(XYZcamera(:,:,3), [5 5]);
%XYZcamera(:,:,3)

%X = reshape(XYZcamera,[],3);
%validID = find(X(:,3)~=0);
%X = X(validID,:);
%points2ply('filter2.ply', X);

depth = XYZcamera(:,:,3);

author: jianxiongxiao

repo: ProfXkit


text:
function cx = zscal ( n, ca, cx, incx )

%*****************************************************************************80
%
%% ZSCAL scales a complex vector by a constant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, complex CA, the multiplier.
%
%    Input, complex CX(*), the vector to be scaled.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Output, complex CX(*), the scaled vector.
%
  cx(1:incx:1+(n-1)*incx) = ca * cx(1:incx:1+(n-1)*incx);

  return
end