text:
function unannotatedRxns = getUnannotatedReactionsFromPubSeedSpreadsheets(infoFilePath,inputDataFolder,spreadsheetFolder)
% Part of the DEMETER pipeline. Prepares input file for the refinement
% based on comparative genomics analyses.
% Gets all the reactions that were not found in the respective organism
% through comparative genomics to remove them from the draft
% reconstructions.
%
% USAGE:
%   unannotatedRxns = getUnannotatedReactionsFromPubSeedSpreadsheets(infoFilePath,inputDataFolder,spreadsheetFolder)
%
% INPUTS
% infoFilePath          File with information on reconstructions to refine
% inputDataFolder       Folder to save propagated data to (default: folder
%                       in current path called "InputData")
% spreadsheetFolder     Folder with comparative genomics data retrieved
%                       from PubSEED in spreadsheet format if available.
%                       For an example of the required format, see
%                       cobratoolbox/papers/2021_demeter/exampleSpreadsheets.
%
% .. Authors:
%       - Almut Heinken, 06/2020

% get PubSEED IDs of new organisms to reconstruct
infoFile = readInputTableForPipeline(infoFilePath);

% load reactions
reactions = readInputTableForPipeline('InReactions.txt');

% get all spreadsheets
dInfo = dir(spreadsheetFolder);
fileList={dInfo.name};
fileList=fileList';
fileList(find(strcmp(fileList(:,1),'.')),:)=[];
fileList(find(strcmp(fileList(:,1),'..')),:)=[];

unannotatedRxns={};
cnt=1;

for i=1:length(fileList)
    i
    % read spreadsheets
    spreadsheet = readInputTableForPipeline([spreadsheetFolder filesep fileList{i}]);

    for j=2:size(spreadsheet,1)
        % replace PubSeed with AGORA Model IDs
        % some entries in the comparative genomics spreadsheet have no
        % reconstruction -> skip
        if ~isempty(find(strcmp(infoFile(:,2),spreadsheet{j,1})))
            spreadsheet{j,1}=infoFile{find(strcmp(infoFile(:,2),spreadsheet{j,1})),1};
            for k=2:size(spreadsheet,2)
                % get the ones that are empty so genoem annotation was not
                % found
                % no drug genes with many associated reactions
                if isempty(spreadsheet{j,k}) && ~any(strcmp(spreadsheet{1,k},{'eUidA','cUidA','UidB','UidP','UidABC'}))
                    getRxns=reactions(find(strcmp(reactions(:,1),spreadsheet{1,k})),2);
                    for r=1:length(getRxns)
                        unannotatedRxns{cnt,1}=spreadsheet{j,1};
                        unannotatedRxns{cnt,2}=getRxns{r};
                        cnt=cnt+1;
                    end
                end
            end
        end
    end
end

% read comparative genomics data
genomeAnnotation = readInputTableForPipeline([inputDataFolder filesep 'gapfilledGenomeAnnotation.txt']);

% read translated SEED reactions
translateRxns = readInputTableForPipeline('ReactionTranslationTable.txt');

% remove duplicate reactions for organisms, gap-filled reactions present in
% comparative genomics spreadsheets
delArray=[];
cnt=1;
orgs=unique(unannotatedRxns(:,1));
for i=1:length(orgs)
    i
    getRxns=unannotatedRxns(find(strcmp(unannotatedRxns(:,1),orgs{i})),2);
    getRxnInds=find(strcmp(unannotatedRxns(:,1),orgs{i}));
    [uniqueRxns, ~, J]=unique(getRxns);
    cntRxns = histc(J, 1:numel(uniqueRxns));
    duplRxns=uniqueRxns(cntRxns>1);
    if ~isempty(duplRxns)
        for j=1:length(duplRxns)
            findInds=getRxnInds(find(strcmp(getRxns(:,1),duplRxns{j})),1);
            for k=2:length(findInds)
                delArray(cnt,1)=findInds(k);
                cnt=cnt+1;
            end
        end
    end
    annRxns=genomeAnnotation(find(strcmp(orgs{i},genomeAnnotation(:,1))),2);
    [~,keptRxnInd]=intersect(getRxns,annRxns);
    for j=1:length(keptRxnInd)
        delArray(cnt,1)=getRxnInds(keptRxnInd(j));
        cnt=cnt+1;
    end

    % remove reactions not present in KBase translation table
    [~,notTranslated]=setdiff(getRxns,translateRxns(:,2));
    for j=1:length(notTranslated)
        delArray(cnt,1)=getRxnInds(notTranslated(j));
        cnt=cnt+1;
    end
end
unannotatedRxns(delArray,:)=[];

end

author: opencobra

repo: cobratoolbox


text:
%% ADP for sensorimotor control

%% Simulate the DF
MDF = DFSimulator();
simNF(MDF);
simDF(MDF);
simAL(MDF);
simAE(MDF);
showStiffness(MDF);
%% Simulate the VF
MVF = VFSimulator();
simNF(MVF);
simVF(MVF);
simAL(MVF);
simAE(MVF);
showStiffness(MVF);
%% Valide Fitts Law
%NF
disp('Simulating the NF...')
wdd = 2.2:0.1:6;
r=0.5./(2.^wdd);
T=r;
for i=1:length(r) % generating each single trial
    disp(['Simulating the NF...', 'Trial: #', num2str(i), ' (', num2str(length(r)-i), 'trials left)'])
    % T(i) = fitts_movement_duration_NF(r(i));
    MDF.reset();
    T(i) = MDF.getMovementDuration(r(i));
end
figure
subplot(321)
fit1l = polyfit(wdd,T,1);
plot(wdd,T,'+',wdd,fit1l(2)+fit1l(1).*wdd,'linewidth',1.5);
axis([1,7,0,0.8]);
xlabel('\fontsize{12}log_2(2d/s)')
ylabel('\fontsize{12}Movement time t_f (sec)')
title('\fontsize{12}A');
legend('\fontsize{12}Movement duration of each trial', '\fontsize{12}Least-squares fit using the log law');
subplot(322)
fit1p = polyfit(log(0.25./r),log(T),1);
plot(log(0.25./r),log(T),'+',log(0.25./r),fit1p(2)+fit1p(1).*log(0.25./r),'linewidth',1.5);
legend('\fontsize{12}Movement duration  of each trial', '\fontsize{12}Least-squares fit using the power law');
axis([0,4,-2,0]);
xlabel('\fontsize{12}ln(d/s)');
ylabel('\fontsize{12}ln(t_f)');
title('\fontsize{12}B');
% VF
Tvf = zeros(size(r));
for i=1:length(r)
    disp(['Simulating the VF...', 'Trial: #', num2str(i), ' (', num2str(length(r)-i), 'trials left)']);
    % Tvf(i) = fitts_movement_duration_VF(r(i));
    Tvf(i) = MVF.getPostLearningMovementDuration(r(i));
end
%figure
%subplot(121)
subplot(323)
fit2l = polyfit(wdd,Tvf,1);
plot(wdd,Tvf,'+',wdd,fit2l(2)+fit2l(1).*wdd,'linewidth',1.5);
axis([1,7,0,1.5]);
xlabel('\fontsize{12}log_2(2d/s)');
ylabel('\fontsize{12}Movement time t_f (sec)');
title('\fontsize{12}C');
legend('\fontsize{12}Movement duration of each trial', '\fontsize{12}Least-squares fit using the log law');
%subplot(122)
subplot(324)
fit2p = polyfit(log(0.25./r),log(Tvf),1);
plot(log(0.25./r),log(Tvf),'+',log(0.25./r),fit2p(2)+fit2p(1).*log(0.25./r),'linewidth',1.5);
legend('\fontsize{12}Movement duration  of each trial', '\fontsize{12}Least-squares fit using the power law');
axis([0,4,-2,0.7]);
xlabel('\fontsize{12}ln(d/s)');
ylabel('\fontsize{12}ln(t_f)');
title('\fontsize{12}D');
% DF
Tdf = zeros(size(r));
for i=1:length(r)
    disp(['Simulating the DF...', 'Trial: #', num2str(i), ' (', num2str(length(r)-i), 'trials left)']);
    MDF.K = MDF.Ko;
    Tdf(i) = MDF.getMovementDuration(r(i));
end
%
%figure
%subplot(121);
subplot(325)
fit3l = polyfit(wdd,Tdf,1);
plot(wdd,Tdf,'+',wdd,fit3l(2)+fit3l(1).*wdd,'linewidth',1.5);
axis([1,7,0,0.8]);
xlabel('\fontsize{12}log_2(2d/s)');
ylabel('\fontsize{12}Movement time t_f (sec)');
title('\fontsize{12}E');
legend('\fontsize{12}Movement duration of each trial', '\fontsize{12}Least-squares fit using the log law');
%subplot(122);
subplot(326)
fit3p = polyfit(log(0.25./r),log(Tdf),1);
plot(log(0.25./r),log(Tdf),'+',log(0.25./r),fit3p(2)+fit3p(1).*log(0.25./r),'linewidth',1.5);
legend('\fontsize{12}Movement duration  of each trial', '\fontsize{12}Least-squares fit using the power law');
axis([0,4,-2,0]);
xlabel('\fontsize{12}ln(d/s)');
ylabel('\fontsize{12}ln(t_f)');
title('\fontsize{12}F');
% Display the fitted parameters
disp('Fitting results:')
disp(['NF Log   Law:', 'a=',num2str(fit1l(1)), ' b=', num2str(fit1l(2))]);
disp(['NF Power Law:', 'a=',num2str(fit1p(1)), ' b=', num2str(fit1p(2))]);
disp(['VF Log   Law:', 'a=',num2str(fit2l(1)), ' b=', num2str(fit2l(2))]);
disp(['VF Power Law:', 'a=',num2str(fit2p(1)), ' b=', num2str(fit2p(2))]);
disp(['DF Log   Law:', 'a=',num2str(fit3l(1)), ' b=', num2str(fit3l(2))]);
disp(['DF Power Law:', 'a=',num2str(fit3p(1)), ' b=', num2str(fit3p(2))]);

author: yu-jiang

repo: radpbook

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


text:
% Graph of Buckley-Leverett flux function

% Theory in Section 9.4 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% 	See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program generates Fig. 9.10 in the book

		% ...............Input.........................................
a = 0.5;	% Parameter in flux function on page 371
		% ............End of input.....................................
x = 0.0:0.05:1.0; den = x.*x + 0.5*(1-x).^2; f = (x.*x)./den;
figure(1), clf, plot(x,f,'-')

repo: Aero-matlab


text:
function printSolutionSummary(algname, problem, inliers, info)

[true_positive_rate, false_positive_rate, precision, recall] = ...
    detectionStats(problem, inliers);

outliers = setdiff(1:problem.N, inliers);

fprintf('%s\n', algname);
fprintf('   Iterations: %d\n', info.Iterations);
fprintf('   Correctly estimated outliers: %d / %d\n', ...
    length(intersect(problem.outliers, outliers)), ...
    length(problem.outliers));
if(length(outliers) ~= length(problem.outliers))
    fprintf('   Num. est. outliers: %d\n', length(outliers));
end
fprintf('   True positive rate: %.3f%%\n', ...
    true_positive_rate * 100);
fprintf('   False positive rate: %.3f%%\n', ...
    false_positive_rate * 100);
fprintf('   Precision: %.3f %%\n', precision * 100);
fprintf('   Recall: %.3f %%\n', recall * 100);
if isfield(info, 'time')
    fprintf('   Exec. time: %.2f sec.\n', info.time);
end

end

author: MIT-SPARK

repo: GNC-and-ADAPT


text:
function check = bernoulli_check ( a )

%*****************************************************************************80
%
%% BERNOULLI_CHECK checks the parameter of the Bernoulli CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the PDF.
%    0.0D+00 <= A <= 1.0.
%
%    Output, logical CHECK, is TRUE if the parameters are OK.
%
  if ( a < 0.0 | 1.0 < a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BERNOULLI_CHECK - Fatal error!\n' );
    fprintf ( 1, '  A < 0 or 1 < A.\n' );
    check = 0;
  else
    check = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [status,datafiles] = writelevel(filename,datafiles)

% print count of renamings.
%datafiles=cell(10,10);

s=signal(filename);
str=history(s);
if ~isempty(datafiles)
  m=length(datafiles(:,1))+1;
else
  m=1;
end
n=0;
status=-1;
lasttoken='';
for i=1:length(str(:,1))
  [token,rem]=strtok(str(i,:),'>');
  rem=rem(2:length(rem));
  [token,rem]=strtok(rem,'<');
  if token & ~strcmp(lasttoken,token)
    lasttoken=token;
    if n
      m=m+1;
      for i1=1:n
	datafiles(m,i1)=datafiles(m-1,i1);
      end
    end
    n=n+1;
    datafiles(m,n)=cellstr(token);
    status=0;
  end
end

[path,fname,ext] = fileparts(filename);

if ~strcmp(lasttoken,fname)
    if n
      m=m+1;
      for i1=1:n
	datafiles(m,i1)=datafiles(m-1,i1);
      end
    end
    n=n+1;
    datafiles(m,n)=cellstr(fname);
end

datafiles(cellfun('isempty',datafiles))=cellstr('');

return;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function params = get_params_value()
% constant parameters
params.c = physconst('LightSpeed');% Speed of light in air (m/s)
params.fc = 77e9; % Center frequency (Hz)
params.lambda = params.c/params.fc;
params.Rx = 4;
params.Tx = 2;

% configuration parameters
params.Fs = 4*10^6;
params.sweepSlope = 21.0017e12;
params.samples = 128;
params.loop = 255;

params.Tc = 120e-6; % us
params.fft_Rang = 134; % 134=>128
params.fft_Vel = 256;
params.fft_Ang = 128;
params.num_crop = 3;
params.max_value = 1e+04; % data WITH 1843

% Creat grid table
freq_res = params.Fs/params.fft_Rang; % range_grid
freq_grid = (0:params.fft_Rang-1).'*freq_res;
params.rng_grid = freq_grid*params.c/params.sweepSlope/2; % d=frediff_grid*c/sweepSlope/2;

w = linspace(-1,1,params.fft_Ang); % angle_grid
params.agl_grid = asin(w)*180/pi; % [-1,1]->[-pi/2,pi/2]

% velocity_grid
dop_grid = fftshiftfreqgrid(params.fft_Vel,1/params.Tc); % now fs is equal to 1/Tc
params.vel_grid = dop_grid*params.lambda/2;   % unit: m/s, v = lamda/4*[-fs,fs], dopgrid = [-fs/2,fs/2]

end

author: Xiangyu-Gao


text:
function [metabolite_structure,IDsAdded,InchiKeyList,InchiStringList] = generateInchiFromMol(metabolite_structure,folder, inchiKey, smiles,formula,startSearch,endSearch)
% This function adds inchistrings to the metabolite structure entries if
% mol files are available (in folder) and no inchistring is present in
% metabolite structure for a given entry. It also calculated inchikeys and
% smiles if desired.
%
% INPUT
% metabolite_structure  metabolite structure
% folder                folder containing the mol files
% inchiKey              if 1 (default), then inchiKey will be generated
%                       from mol file
% smiles                if 1 (default), then smiles will be generated
%                       from mol file
% formula               if 1 (default), then inchiKey will be generated
%                       from inchistring file
% startSearch           specify where the search should start in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% endSearch             specify where the search should end in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
%
% OUTPUT
% metabolite_structure  Updated metabolite structure
% IDsAddeed             List of IDs added to the metabolite structure
% InchiKeyList          List of inchikeys added
% InchiStringList       List of inchistrings added
%
%
% Ines Thiele, 2020-2021

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

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

% folder    folder containing all mol files
annotationSource = ['Mol files from ' folder];
annotationType = 'automatic';

Mets = fieldnames(metabolite_structure);

if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(Mets);
end
cnt = 1;
for i = startSearch : endSearch

    VMHId{cnt,1} = Mets{i};
    if isfield(metabolite_structure.(Mets{i}),'VMHId')
        VMHId{cnt,2} = metabolite_structure.(Mets{i}).VMHId;
    else
        VMHId{cnt,2} = Mets{i};
        VMHId{cnt,2} =regexprep(  VMHId{i,2},'VMH_','');
        VMHId{cnt,2} =regexprep(  VMHId{i,2},'__','_');
    end
    cnt = cnt + 1;
end
a = 1;
IDsAdded = ''; a =  1;
InchiStringList = ''; b =1;
InchiKeyList = ''; c= 1;
%folder = 'FileDumps/molFiles/molFiles_AGORA2/';
files = dir(folder);
space =' ';

% this could be accelerated
for i = 1 :length(files)
    met = regexprep(files(i).name,'.mol','');
    match = VMHId(find(ismember(VMHId(:,2),met),1),1);

    if contains(files(i).name,'.mol') && ~isempty(match)

        metabolite_structure.(match{1}).hasmolfile = 1;
        metabolite_structure.(match{1}).hasmolfile_source = [annotationSource,':',annotationType,':',datestr(now)];

        if isempty(metabolite_structure.(match{1}).inchiString) || length(find(isnan(metabolite_structure.(match{1}).inchiString)))>0
            if ismac
                [status, result]=system(strcat('/usr/local/bin/obabel',[space folder filesep files(i).name],' -o inchi -h' ));
            elseif ispc % pc
                [status, result]=system(strcat('obabel',[space folder filesep files(i).name],' -o inchi -h' ));
            elseif isunix
                [status, result]=system(strcat('obabel',[space folder filesep files(i).name],' -o inchi -h' ));
            end

            if contains(result,'1 molecule converted')
                %result = regexprep(result,'1 molecule converted','');
                [tok,rem] = strtok(result,'=');
                result = regexprep(rem,'1 molecule converted','');
                result = regexprep(result,'\n','');
                if strfind(result,'==')
                    tokr = split(result,'==');
                    result = tokr{1};
                end
                if ~isempty(result)
                     inchiString = strcat('InChI',result);
                else
                    inchiString = {};
                end
                % find metabolite in metabolite_structure
                if ~isempty(inchiString)
                    inchiString = strcat('InChI',result);
                    metabolite_structure.(match{1}).inchiString = inchiString;
                    metabolite_structure.(match{1}).inchiString_source = [annotationSource,':',annotationType,':',datestr(now)];
                    IDsAdded{a,1} = match{1};
                    IDsAdded{a,2} = 'inchiString';
                    IDsAdded{a,3} = metabolite_structure.(match{1}).inchiString;
                    a = a + 1;

                    InchiStringList{b,1} = match{1};
                    InchiStringList{b,2} = strcat(folder,files(i).name);
                    InchiStringList{b,3} = metabolite_structure.(match{1}).inchiString;
                    InchiStringList{b,4} = metabolite_structure.(match{1}).inchiString_source;
                    b = b + 1;
                end

                if inchiKey
                    if(isempty(metabolite_structure.(match{1}).inchiKey) || length(find(isnan(metabolite_structure.(match{1}).inchiKey)))>0)
                        if    ~isempty(metabolite_structure.(match{1}).inchiString) &&  length(find(isnan(metabolite_structure.(match{1}).inchiString)))==0
                            inchiString =  metabolite_structure.(match{1}).inchiString;
                            [result] = convertInchiString2format(inchiString,'inchiKey');
                            metabolite_structure.(match{1}).inchiKey = result;
                            metabolite_structure.(match{1}).inchiKey_source = [annotationSource,':',annotationType,':',datestr(now)];
                            IDsAdded{a,1} = match{1};
                            IDsAdded{a,2} = 'inchiKey';
                            IDsAdded{a,3} = metabolite_structure.(match{1}).inchiKey;
                            a = a + 1;

                            InchiKeyList{c,1} = match{1};
                            InchiKeyList{c,2} = strcat(folder,files(i).name);
                            InchiKeyList{c,3} = metabolite_structure.(match{1}).inchiKey;
                            InchiKeyList{c,4} = metabolite_structure.(match{1}).inchiKey_source;
                            c = c + 1;
                        end
                    end
                end
                if smiles
                    if (isempty(metabolite_structure.(match{1}).smile) || length(find(isnan(metabolite_structure.(match{1}).smile)))>0)
                        if ~isempty(metabolite_structure.(match{1}).inchiString) &&  length(find(isnan(metabolite_structure.(match{1}).inchiString)))==0
                            inchiString =  metabolite_structure.(match{1}).inchiString;

                            [result] = convertInchiString2format(inchiString,'smiles');
                            result = regexprep(result,'\W$','');
                            metabolite_structure.(match{1}).smile = result;
                            metabolite_structure.(match{1}).smile_source = [annotationSource,':',annotationType,':',datestr(now)];
                            IDsAdded{a,1} = match{1};
                            IDsAdded{a,2} = 'smiles';
                            IDsAdded{a,3} = metabolite_structure.(match{1}).smile;
                            a = a + 1;
                        end
                    end
                end
            end
        end
    end
end
if formula
    [metabolite_structure] = addMetFormulaCharge(metabolite_structure,startSearch,endSearch);
end


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren and Li Xu, "On Vectorization of Deep Convolutional Neural Networks for Vision Tasks",
% The 29th AAAI Conference on Artificial Intelligence (AAAI-15). Austin, Texas, USA, January 25-30, 2015
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
function mnist_configure()
    global config;
    % general configurations for both training and testing
    config.input_size = [32 32];
    config.chs = 1;
    config.forward_pass_scheme = {'conv_v', 'pool', 'conv_v', 'pool', 'conv_v', 'full', 'full', 'full', 'out'};
    config.nonlinearity = 'relu';   % 'relu', 'tanh', 'sigmoid'
    config.output_activation = 'softmax';   % 'softmax', 'inherit'
    config.cost_function = 'cross entropy'; % 'cross entropy', 'L2 norm'
    config.kernel_size = [5 5; 5 5; 5 5];
    config.conv_hidden_size = [6 20 150];
    config.full_hidden_size = [150 150];
    config.output_size = [1 1 10];
    config.batch_size = 100;
    config.compute_device = 'GPU';

    % the following items are only for training
    config.learning_rate = 0.01;
    config.weight_range = 2;
    config.decay = 5e-7 / 10;
    config.normalize_init_weights = 1;
    config.dropout_full_layer = 1;
    config.optimization = 'adagrad';
end

author: jimmy-ren


text:
function [fg,filename] = mtrImportFibers(filename,xform,bEndOnly, coordspace)
%
% [fg,filename] = mtrImportFibers(filename, [xform], [bEndOnly], [coordspace])
%
% Imports non-mrDiffusion fiber/path file formats. Currently supports:
%
% 1. PDB (pathway data base)
% 2. Bfloat (ConTrack/Camino)
% 3. tck (MRtrix)
%
% Parameters
% ----------
% filename: Full path to the file containing the fiber group. Currently
% supports:
%   1. PDB (pathway data base)
%   2. Bfloat (ConTrack/Camino)
%   3. MRtrix (MRtrix fibers)
% xform: A transformation matrix to apply to the coordinates of all the
% fibers in the fiber group, optional.
% bEndOnly: For Bfloat files, whether to read in only the end points,
% optional, default to read the full fibers.
% coordspace: string {'acpc' | 'img' | ...}, the coordinate space in which
% the fiber coordinates are defined. optional. Defaults to an empty matrix
%
% 2011.07.21 - LMP changed the call for Bfloat file types to
% dtiLoadBfloatPaths
%
% 2009 Stanford, Vista

if notDefined('bEndOnly'), bEndOnly = 0; end
if notDefined('xform'),    xform = eye(4); end

%persistent defaultPath;
defaultPath = pwd;

if (~exist('filename','var') || isempty(filename))
    if(isempty(defaultPath)), filename = pwd;
    else                      filename = defaultPath;
    end
end

if (~exist('xform','var') || isempty(xform)), xform = eye(4); end

% Read the file name from the user
if(isdir(filename)) || ~exist(filename,'file')
    if(~isempty(defaultPath)), filename = defaultPath; end
    [f,p] = uigetfile({'*.pdb','ConTrack PDB *.pdb';'*.Bfloat','ConTrack/Camino *.Bfloat';'*.*','all files'},'Select fiber file...',filename);
    if(isnumeric(f)), disp('Cancel.'); return; end
    % defaultPath = p;
    filename = fullfile(p,f);
end

%
[p,f,e] = fileparts(filename); %#ok<ASGLU>
if(isempty(e)), e = '.pdb'; end

fg = [];
switch(e)
    case '.pdb',
        %disp('Trying to import PDB paths...');

        % Read the whole file into a string, str.  The str is decoded in
        % the function below.
        fid = fopen(filename, 'r');
        str = fread(fid, inf,'uint8');
        fclose(fid);

        % MetroTrack is the old name.  It reads the pdb file.
        mt = dtiLoadMetrotracPathsFromStr(str,xform);

        % Copy the loaded data in the mrDiffusion format
        if ~isempty(mt)
            fg = dtiNewFiberGroup;
            fg.fibers = mt.pathways;
            fg.name = f;
            fg.colorRgb = [200 200 100];
            fg.pathwayInfo = mt.pathwayInfo;
            for ss = 1:length(mt.statHeader)
                statstruct.name = mt.statHeader(ss).agg_name;
                statstruct.uid=mt.statHeader(ss).uid;
                statstruct.ile=mt.statHeader(ss).is_luminance_encoding;
                statstruct.icpp=mt.statHeader(ss).is_computed_per_point;
                statstruct.ivs=mt.statHeader(ss).is_viewable_stat;
                statstruct.agg=mt.statHeader(ss).agg_name;
                statstruct.lname=mt.statHeader(ss).local_name;
                for pp = 1:length(mt.pathwayInfo)
                    statstruct.stat(pp) = mt.pathwayInfo(pp).pathStat(ss);
                end
                fg.params{ss} =  statstruct;
            end
        end

    case '.Bfloat',
        mt = dtiLoadBfloatPaths(filename,xform,bEndOnly);
        if ~isempty(mt)
            fg = fgCreate;
            fg.fibers = mt.pathways;
            fg.name = f;
            fg.colorRgb = [200 200 100];
            for ss = 1:length(mt.statHeader)
                statstruct.name = mt.statHeader(ss).agg_name;
                for pp = 1:length(mt.pathwayInfo)
                    statstruct.stat(pp) = mt.pathwayInfo(pp).pathStat(ss);
                end
                fg.params{ss} =  statstruct;
            end
        end

    case '.tck'
        fg = dtiImportFibersMrtrix(filename);

    otherwise,
        error('unknown format: %s\n',e);
end

% These are always going to be in acpc space, so there we go:
if ~notDefined('coordspace')
    fgSet(fg, 'coordspace', coordspace);
else
    fgSet(fg, 'coordspace', []);
end

fprintf('Fibers loaded from file %s\n',f);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% first casadi test for mpc fpr mobile robots
clear all
close all
clc

% CasADi v3.1.1
addpath('C:\Users\mehre\OneDrive\Desktop\CasADi\casadi-windows-matlabR2016a-v3.4.5')
import casadi.*

T = 0.2; %[s]
N = 10; % prediction horizon
rob_diam = 0.3;

v_max = 0.6; v_min = -v_max;
omega_max = pi/4; omega_min = -omega_max;

x = SX.sym('x'); y = SX.sym('y'); theta = SX.sym('theta');
states = [x;y;theta]; n_states = length(states);

v = SX.sym('v'); omega = SX.sym('omega');
controls = [v;omega]; n_controls = length(controls);
rhs = [v*cos(theta);v*sin(theta);omega]; % system r.h.s

f = Function('f',{states,controls},{rhs}); % nonlinear mapping function f(x,u)
U = SX.sym('U',n_controls,N); % Decision variables (controls)
P = SX.sym('P',n_states + n_states);
% parameters (which include at the initial state of the robot and the reference state)

X = SX.sym('X',n_states,(N+1));
% A vector that represents the states over the optimization problem.

obj = 0; % Objective function
g = [];  % constraints vector

Q = zeros(3,3); Q(1,1) = 1;Q(2,2) = 5;Q(3,3) = 0.1; % weighing matrices (states)
R = zeros(2,2); R(1,1) = 0.5; R(2,2) = 0.05; % weighing matrices (controls)

st  = X(:,1); % initial state
g = [g;st-P(1:3)]; % initial condition constraints
for k = 1:N
    st = X(:,k);  con = U(:,k);
    obj = obj+(st-P(4:6))'*Q*(st-P(4:6)) + con'*R*con; % calculate obj
    st_next = X(:,k+1);
    f_value = f(st,con);
    st_next_euler = st+ (T*f_value);
    g = [g;st_next-st_next_euler]; % compute constraints
end
% make the decision variable one column  vector
OPT_variables = [reshape(X,3*(N+1),1);reshape(U,2*N,1)];

nlp_prob = struct('f', obj, 'x', OPT_variables, 'g', g, 'p', P);

opts = struct;
opts.ipopt.max_iter = 2000;
opts.ipopt.print_level =0;%0,3
opts.print_time = 0;
opts.ipopt.acceptable_tol =1e-8;
opts.ipopt.acceptable_obj_change_tol = 1e-6;

solver = nlpsol('solver', 'ipopt', nlp_prob,opts);

args = struct;

args.lbg(1:3*(N+1)) = 0;
args.ubg(1:3*(N+1)) = 0;

args.lbx(1:3:3*(N+1),1) = -2; %state x lower bound
args.ubx(1:3:3*(N+1),1) = 2; %state x upper bound
args.lbx(2:3:3*(N+1),1) = -2; %state y lower bound
args.ubx(2:3:3*(N+1),1) = 2; %state y upper bound
args.lbx(3:3:3*(N+1),1) = -inf; %state theta lower bound
args.ubx(3:3:3*(N+1),1) = inf; %state theta upper bound

args.lbx(3*(N+1)+1:2:3*(N+1)+2*N,1) = v_min; %v lower bound
args.ubx(3*(N+1)+1:2:3*(N+1)+2*N,1) = v_max; %v upper bound
args.lbx(3*(N+1)+2:2:3*(N+1)+2*N,1) = omega_min; %omega lower bound
args.ubx(3*(N+1)+2:2:3*(N+1)+2*N,1) = omega_max; %omega upper bound
%----------------------------------------------
% ALL OF THE ABOVE IS JUST A PROBLEM SET UP

% THE SIMULATION LOOP SHOULD START FROM HERE
%-------------------------------------------
t0 = 0;
x0 = [0.1 ; 0.1 ; 0.0];    % initial condition.
xs = [1.5 ; 1.5 ; 0.0]; % Reference posture.

xx(:,1) = x0; % xx contains the history of states
t(1) = t0;

u0 = zeros(N,2);        % two control inputs for each robot
X0 = repmat(x0,1,N+1)'; % initialization of the states decision variables

sim_tim = 20; % total sampling times

% Start MPC
mpciter = 0;
xx1 = [];
u_cl=[];

% the main simulaton loop... it works as long as the error is greater
% than 10^-6 and the number of mpc steps is less than its maximum
% value.
tic
while(norm((x0-xs),2) > 0.05 && mpciter < sim_tim / T)
    args.p   = [x0;xs]; % set the values of the parameters vector
    % initial value of the optimization variables
    args.x0  = [reshape(X0',3*(N+1),1);reshape(u0',2*N,1)];
    sol = solver('x0', args.x0, 'lbx', args.lbx, 'ubx', args.ubx,...
        'lbg', args.lbg, 'ubg', args.ubg,'p',args.p);
    u = reshape(full(sol.x(3*(N+1)+1:end))',2,N)'; % get controls only from the solution
    xx1(:,1:3,mpciter+1)= reshape(full(sol.x(1:3*(N+1)))',3,N+1)'; % get solution TRAJECTORY
    u_cl= [u_cl ; u(1,:)];
    t(mpciter+1) = t0;
    % Apply the control and shift the solution
    [t0, x0, u0] = shift(T, t0, x0, u,f);
    xx(:,mpciter+2) = x0;
    X0 = reshape(full(sol.x(1:3*(N+1)))',3,N+1)'; % get solution TRAJECTORY
    % Shift trajectory to initialize the next step
    X0 = [X0(2:end,:);X0(end,:)];
    mpciter
    mpciter = mpciter + 1;
end;
toc

ss_error = norm((x0-xs),2)
%Draw_MPC_point_stabilization_v1 (t,xx,xx1,u_cl,xs,N,rob_diam)
%-----------------------------------------
%-----------------------------------------
%-----------------------------------------
%    Start MHE implementation from here
%-----------------------------------------
%-----------------------------------------
%-----------------------------------------
% plot the ground truth
figure(1)
subplot(311)
plot(t,xx(1,1:end-1),'b','linewidth',1.5); axis([0 t(end) 0 1.8]);hold on
ylabel('x (m)')
grid on
subplot(312)
plot(t,xx(2,1:end-1),'b','linewidth',1.5); axis([0 t(end) 0 1.8]);hold on
ylabel('y (m)')
grid on
subplot(313)
plot(t,xx(3,1:end-1),'b','linewidth',1.5); axis([0 t(end) -pi/4 pi/2]);hold on
xlabel('time (seconds)')
ylabel('\theta (rad)')
grid on

% Synthesize the measurments
con_cov = diag([0.005 deg2rad(2)]).^2;
meas_cov = diag([0.1 deg2rad(2)]).^2;

r = [];
alpha = [];
for k = 1: length(xx(1,:))-1
    r = [r; sqrt(xx(1,k)^2+xx(2,k)^2)  + sqrt(meas_cov(1,1))*randn(1)];
    alpha = [alpha; atan(xx(2,k)/xx(1,k))      + sqrt(meas_cov(2,2))*randn(1)];
end
y_measurements = [ r , alpha ];

% Plot the cartesian coordinates from the measurements used
figure(1)
subplot(311)
plot(t,r.*cos(alpha),'r','linewidth',1.5); hold on
grid on
legend('Ground Truth','Measurement')
subplot(312)
plot(t,r.*sin(alpha),'r','linewidth',1.5); hold on
grid on

% plot the ground truth mesurements VS the noisy measurements
figure(2)
subplot(211)
plot(t,sqrt(xx(1,1:end-1).^2+xx(2,1:end-1).^2),'b','linewidth',1.5); hold on
plot(t,r,'r','linewidth',1.5); axis([0 t(end) -0.2 3]); hold on
ylabel('Range: [ r (m) ]')
grid on
legend('Ground Truth','Measurement')
subplot(212)
plot(t,atan(xx(2,1:end-1)./xx(1,1:end-1)),'b','linewidth',1.5); hold on
plot(t,alpha,'r','linewidth',1.5); axis([0 t(end) 0.2 1]); hold on
ylabel('Bearing: [ \alpha (rad) ]')
grid on

% The following two matrices contain what we know about the system, i.e.
% the nominal control actions applied (measured) and the range and bearing
% measurements.
%-------------------------------------------------------------------------
u_cl;
y_measurements;

% Let's now formulate our MHE problem
%------------------------------------
T = 0.2; %[s]
N_MHE = 6; % prediction horizon

v_max = 0.6; v_min = -v_max;
omega_max = pi/4; omega_min = -omega_max;

x = SX.sym('x'); y = SX.sym('y'); theta = SX.sym('theta');
states = [x;y;theta]; n_states = length(states);
v = SX.sym('v'); omega = SX.sym('omega');
controls = [v;omega]; n_controls = length(controls);
rhs = [v*cos(theta);v*sin(theta);omega]; % system r.h.s
f = Function('f',{states,controls},{rhs}); % MOTION MODEL

r = SX.sym('r'); alpha = SX.sym('alpha'); % range and bearing
measurement_rhs = [sqrt(x^2+y^2); atan(y/x)];
h = Function('h',{states},{measurement_rhs}); % MEASUREMENT MODEL
%y_tilde = [r;alpha];

% Decision variables
U = SX.sym('U',n_controls,N_MHE);    %(controls)
X = SX.sym('X',n_states,(N_MHE+1));  %(states) [remember multiple shooting]

P = SX.sym('P', 2 , N_MHE + (N_MHE+1));
% parameters (include r and alpha measurements as well as controls measurements)

V = inv(sqrt(meas_cov)); % weighing matrices (output)  y_tilde - y
W = inv(sqrt(con_cov)); % weighing matrices (input)   u_tilde - u

obj = 0; % Objective function
g = [];  % constraints vector
for k = 1:N_MHE+1
    st = X(:,k);
    h_x = h(st);
    y_tilde = P(:,k);
    obj = obj+ (y_tilde-h_x)' * V * (y_tilde-h_x); % calculate obj
end

for k = 1:N_MHE
    con = U(:,k);
    u_tilde = P(:,N_MHE+ k);
    obj = obj+ (u_tilde-con)' * W * (u_tilde-con); % calculate obj
end

% multiple shooting constraints
for k = 1:N_MHE
    st = X(:,k);  con = U(:,k);
    st_next = X(:,k+1);
    f_value = f(st,con);
    st_next_euler = st+ (T*f_value);
    g = [g;st_next-st_next_euler]; % compute constraints
end

% make the decision variable one column  vector
OPT_variables = [reshape(X,3*(N_MHE+1),1);reshape(U,2*N_MHE,1)];

nlp_mhe = struct('f', obj, 'x', OPT_variables, 'g', g, 'p', P);

opts = struct;
opts.ipopt.max_iter = 2000;
opts.ipopt.print_level =0;%0,3
opts.print_time = 0;
opts.ipopt.acceptable_tol =1e-8;
opts.ipopt.acceptable_obj_change_tol = 1e-6;

solver = nlpsol('solver', 'ipopt', nlp_mhe,opts);

args = struct;

args.lbg(1:3*(N_MHE)) = 0;
args.ubg(1:3*(N_MHE)) = 0;

args.lbx(1:3:3*(N_MHE+1),1) = -2; %state x lower bound
args.ubx(1:3:3*(N_MHE+1),1) = 2; %state x upper bound
args.lbx(2:3:3*(N_MHE+1),1) = -2; %state y lower bound
args.ubx(2:3:3*(N_MHE+1),1) = 2; %state y upper bound
args.lbx(3:3:3*(N_MHE+1),1) = -pi/2; %state theta lower bound
args.ubx(3:3:3*(N_MHE+1),1) = pi/2; %state theta upper bound

args.lbx(3*(N_MHE+1)+1:2:3*(N_MHE+1)+2*N_MHE,1) = v_min; %v lower bound
args.ubx(3*(N_MHE+1)+1:2:3*(N_MHE+1)+2*N_MHE,1) = v_max; %v upper bound
args.lbx(3*(N_MHE+1)+2:2:3*(N_MHE+1)+2*N_MHE,1) = omega_min; %omega lower bound
args.ubx(3*(N_MHE+1)+2:2:3*(N_MHE+1)+2*N_MHE,1) = omega_max; %omega upper bound
%----------------------------------------------
% ALL OF THE ABOVE IS JUST A PROBLEM SET UP

% MHE Simulation loop starts here
%------------------------------------------
X_estimate = []; % X_estimate contains the MHE estimate of the states
U_estimate = []; % U_estimate contains the MHE estimate of the controls

U0 = zeros(N_MHE,2);   % two control inputs for each robot
X0 = zeros(N_MHE+1,3); % initialization of the states decision variables

% Start MHE
mheiter = 0;

U0 = u_cl(1:N_MHE,:); % initialize the control actions by the measured
% initialize the states from the measured range and bearing
X0(:,1:2) = [y_measurements(1:N_MHE+1,1).*cos(y_measurements(1:N_MHE+1,2)),...
    y_measurements(1:N_MHE+1,1).*sin(y_measurements(1:N_MHE+1,2))];

tic
for k = 1: size(y_measurements,1) - (N_MHE)
    mheiter = k
    % Get the measurements window and put it as parameters in p
    args.p   = [y_measurements(k:k+N_MHE,:)',u_cl(k:k+N_MHE-1,:)'];
    % initial value of the optimization variables
    args.x0  = [reshape(X0',3*(N_MHE+1),1);reshape(U0',2*N_MHE,1)];
    sol = solver('x0', args.x0, 'lbx', args.lbx, 'ubx', args.ubx,...
        'lbg', args.lbg, 'ubg', args.ubg,'p',args.p);
    U_sol = reshape(full(sol.x(3*(N_MHE+1)+1:end))',2,N_MHE)'; % get controls only from the solution
    X_sol = reshape(full(sol.x(1:3*(N_MHE+1)))',3,N_MHE+1)'; % get solution TRAJECTORY
    X_estimate = [X_estimate;X_sol(N_MHE+1,:)];
    U_estimate = [U_estimate;U_sol(N_MHE,:)];

    % Shift trajectory to initialize the next step
    X0 = [X_sol(2:end,:);X_sol(end,:)];
    U0 = [U_sol(2:end,:);U_sol(end,:)];
end;
toc

figure(1)
subplot(311)
plot(t(N_MHE+1:end),X_estimate(:,1),'g','linewidth',1.5); hold on
legend('Ground Truth','Measurement','MHE')

subplot(312)
plot(t(N_MHE+1:end),X_estimate(:,2),'g','linewidth',1.5); axis([0 t(end) 0 1.8]);hold on

subplot(313)
plot(t(N_MHE+1:end),X_estimate(:,3),'g','linewidth',1.5); axis([0 t(end) -pi/4 pi/2]);hold on

author: MMehrez


text:
function mM = minandmax2est(f, N)
%MINANDMAX2EST   Estimates the minimum and maximum of a SEPARABLEAPPROX.
%   mM = MINANDMAX2EST(F) returns estimates for the minimum and maximum of the
%   SEPARABLEAPPROX F over its domain.  mM is a vector of length 2 such that
%   mM(1) is the estimated minimum and mM(2) is the estimated maximum.
%
%   mM = MINANDMAX2EST(F, N) returns estimates for the minimum and maximum of
%   the SEPARABLEAPPROX F over its domain, based on samples on an N by N grid
%   (N = 33 by default).
%
% See also MINANDMAX2.

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

if ( isempty(f) )
    mM = [0, 0];
    return
end

if ( ( nargin < 2 ) || isempty(N) )
    % Default to N = 33:
    N = 33;
end

% Sample f on an appropriate grid:
vals = sample(f, N, N);

% Make result a column vector:
vals = vals(:);

% Get min and max:
mM = [ min(vals), max(vals) ];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Obs = projectLmk(Rob,Sen,Lmk,Obs,Opt)

% PROJECTLMK  Project landmark estimate into sensor's measurement space.
%   Obs = PROJECTLMK(Rob,Sen,Lmk,Obs) projects the landmark Lmk into sensor
%   Sen mounted on robot Rob, and updates the information of the
%   observation structure Obs. The observation model is determined from
%   Sen.type and Lmk.type. It is an error if no model exists for the chosen
%   Sen-Lmk pair.
%
%   The updated fields in Obs are:
%       .sid        % sensor ID
%       .lid        % landmark ID
%       .ltype      % landmark type
%       .vis        % flag: true if landmark is visible
%       .meas.R     % measurement noise cov. matrix
%       .exp.e      % expectation's mean
%       .exp.E      % expectation's covariances matrix
%       .exp.um     % expectation's uncertainty measure
%       .Jac.E_r    % Jacobian wrt robot frame
%       .Jac.E_s    % Jacobian wrt sensor frame
%       .Jac.E_l    % Jacobian wrt landmark state
%
%   See also OBSERVEKNOWNLMKS.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

global Map

% PREVIOUS TASKS
% get landmark range and mean
lr = Lmk.state.r ;        % lmk range in Map
switch Map.type
    case 'ekf'
        l  = Map.x(lr);   % lmk mean
    case 'graph'
        l = Lmk.state.x;  % lmk mean
    otherwise
        error('??? Unknown Map type ''%s''.',Map.type)
end

% PROJECTION FUNCTION
% explore all sensor and landmark types
switch Sen.type

    case {'pinHole'} % camera pinHole

        switch Lmk.type

            case {'eucPnt'} % euclidean point

                % Point3D --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projEucPntIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l) ;

                vis = isVisible(e,depth,Sen.par.imSize);

            case {'idpPnt'} % inverse depth point

                % IDP --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projIdpPntIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l) ;

                vis = isVisible(e,depth,Sen.par.imSize);

            case {'hmgPnt'} % euclidean point

                % Point3D --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projHmgPntIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l) ;

                vis = isVisible(e,depth,Sen.par.imSize);

            case {'ahmPnt'}

                % Point3D --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projAhmPntIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l);

                vis = isVisible(e,depth,Sen.par.imSize);

            case {'fhmPnt'}

                % Point3D --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projFhmPntIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l);

                vis = isVisible(e,depth,Sen.par.imSize);

            case {'plkLin'}

                % Plucker line --> homogeneous line (value and Jacs)
                [e, v, E_rf, E_sf, E_k, E_l] = ...
                    projPlkLinIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    l); % expectation e is a homogeneous line

                % normalize wrt director vector e(1:2):
                ine12 = 1/norm(e(1:2));
                e     = e    * ine12;
                E_rf  = E_rf * ine12;
                E_sf  = E_sf * ine12;
                % E_k   = E_k  * ine12;
                E_l   = E_l  * ine12;

                % 3d Segment from Plucker line and abscissas
                [si,SI_l] = pluckerSegment(l,[Lmk.par.endp.t]);

                % projected segment
                [s, d, S_rf, S_sf, S_k, S_si] = projSegLinIntoPinHoleOnRob(...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    si);

                % segment visibility
                [s,vis] = visibleSegment(s,d,Sen.par.imSize);

            case 'aplLin'

                % Anchored Plucker line --> homogeneous line (value and Jacs)
                [e, v, E_rf, E_sf, E_k, E_l] = ...
                    projAplLinIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    l); % expectation e is a homogeneous line

                % normalize wrt director vector e(1:2):
                ine12 = 1/norm(e(1:2));
                e     = e    * ine12;
                E_rf  = E_rf * ine12;
                E_sf  = E_sf * ine12;
                % E_k   = E_k  * ine12;
                E_l   = E_l  * ine12;

                % 3d Segment from Plucker line and abscissas
                [si,SI_l] = aPluckerSegment(l,[Lmk.par.endp.t]);

                % projected segment
                [s, d, S_rf, S_sf, S_k, S_si] = projSegLinIntoPinHoleOnRob(...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    si);

                % segment visibility
                [s,vis] = visibleSegment(s,d,Sen.par.imSize);

            case 'idpLin'

                % IDP line --> homogeneous line (value and Jacs)
                [s, v, S_rf, S_sf, S_k, S_l] = ...
                    projIdpLinIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    l); % expectation s is a 2d segment
                [e, E_s] = seg2hmgLin(s); % expectation e is a homogeneous line
                E_rf = E_s*S_rf;
                E_sf = E_s*S_sf;
                E_l  = E_s*S_l;

                % normalize wrt director vector e(1:2):
                ine12 = 1/norm(e(1:2));
                e     = e    * ine12;
                E_rf  = E_rf * ine12;
                E_sf  = E_sf * ine12;
                % E_k   = E_k  * ine12;
                E_l   = E_l  * ine12;

                % 3d Segment from Plucker line and abscissas
                [si,SI_l] = idpLinSegment(l,[Lmk.par.endp.t]);

                % projected segment
                [s, d, S_rf, S_sf, S_k, S_si] = projSegLinIntoPinHoleOnRob(...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    si);

                % segment visibility
                [s,vis] = visibleSegment(s,d,Sen.par.imSize,0,Opt.obs.lines.minLength);

            case 'hmgLin'

                % HMG line --> homogeneous line (value and Jacs)
                [s, v, S_rf, S_sf, S_k, S_l] = ...
                    projHmgLinIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    l); % expectation s is a 2d segment
                [e, E_s] = seg2hmgLin(s); % expectation e is a homogeneous line
                E_rf = E_s*S_rf;
                E_sf = E_s*S_sf;
                E_l  = E_s*S_l;

                % normalize wrt director vector e(1:2):
                ine12 = 1/norm(e(1:2));
                e     = e    * ine12;
                E_rf  = E_rf * ine12;
                E_sf  = E_sf * ine12;
                % E_k   = E_k  * ine12;
                E_l   = E_l  * ine12;

                % 3d Segment from Plucker line and abscissas
                [si,SI_l] = hmgLinSegment(l,[Lmk.par.endp.t]);

                % projected segment
                [s, d, S_rf, S_sf, S_k, S_si] = projSegLinIntoPinHoleOnRob(...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    si);

                % segment visibility
                [s,vis] = visibleSegment(s,d,Sen.par.imSize,0,Opt.obs.lines.minLength);

            case 'ahmLin'

                % AHM line --> homogeneous line (value and Jacs)
                [s, v, S_rf, S_sf, S_k, S_l] = ...
                    projAhmLinIntoPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    l); % expectation s is a 2d segment
                [e, E_s] = seg2hmgLin(s); % expectation e is a homogeneous line
                E_rf = E_s*S_rf;
                E_sf = E_s*S_sf;
                E_l  = E_s*S_l;

                % normalize wrt director vector e(1:2):
                ine12 = 1/norm(e(1:2));
                e     = e    * ine12;
                E_rf  = E_rf * ine12;
                E_sf  = E_sf * ine12;
                % E_k   = E_k  * ine12;
                E_l   = E_l  * ine12;

                % 3d Segment from Plucker line and abscissas
                [si,SI_l] = ahmLinSegment(l,[Lmk.par.endp.t]);

                % projected segment
                [s, d, S_rf, S_sf, S_k, S_si] = projSegLinIntoPinHoleOnRob(...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    si);

                % segment visibility
                [s,vis] = visibleSegment(s,d,Sen.par.imSize,0,Opt.obs.lines.minLength);

            otherwise % unknown landmark type for pin hole sensor
                error('??? Unknown landmark type ''%s'' for sensor ''%s''.', Lmk.type, Sen.type);

        end

    case 'pinHoleDepth'  % Pin hole sensor with depth information
        switch Lmk.type
            case 'eucPnt'

                % Point3D --> pixel+depth -(value and Jacobians)-
                [e, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projEucPntIntoPhdOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l) ;

                vis = isVisible(e(1:2,:),e(3,:),Sen.par.imSize);

            otherwise
                error('??? Unknown landmark type ''%s'' for sensor ''%s''.',Lmk.type,Sen.type);
        end

    case {'omniCam'} % Omnidirectional camera

        switch Lmk.type

            case {'eucPnt'} % euclidean point

                % Point3D --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projEucPntIntoOmniCamOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l) ;

                vis = isVisible(e,depth,Sen.par.imSize);

            case {'ahmPnt'}

                % Point3D --> pixel -(value and Jacobians)-
                [e, depth, E_rf, E_sf, E_k, E_d, E_l] = ...
                    projAhmPntIntoOmniCamOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.d, ...
                    l);

                vis = isVisible(e,depth,Sen.par.imSize);

            otherwise % unknown landmark type for pin hole sensor
                error('??? Unknown landmark type ''%s'' for sensor ''%s''.',Lmk.type,Sen.type);

        end

    otherwise % unknown Sensor type
        error('??? Unknown sensor type ''%s''.',Sen.type);

end % sensor type

% COVARIANCES
if strcmp(Map.type,'ekf')

    % Rob-Sen-Lmk range and Jacobian
    if Sen.frameInMap
        rslr  = [Rob.frame.r ; Sen.frame.r ; lr]; % range of robot, sensor, and landmark
        E_rsl = [E_rf E_sf E_l];
    else
        rslr  = [Rob.frame.r ; lr];               % range of robot and landmark
        E_rsl = [E_rf E_l];
    end

    % Expectation covariances matrix
    E = E_rsl*Map.P(rslr,rslr)*E_rsl' ;

else
    E = [];
end

% Other parameters
switch Lmk.type(4:6)
    case 'Lin'
        % for lines, project endpoints with covariances:

        % compute endpoints
        Obs.par.endp(1).e = s(1:2);
        Obs.par.endp(2).e = s(3:4);

        if Map.type == 'ekf'
            % Rob-Sen-Lmk Jacobian of projected segment
            if Sen.frameInMap
                S_rsl = [S_rf S_sf S_si*SI_l];
            else
                S_rsl = [S_rf S_si*SI_l];
            end
            % compute covariances
            S = S_rsl*Map.P(rslr,rslr)*S_rsl'; % segment covariance
            Obs.par.endp(1).E = S(1:2,1:2);
            Obs.par.endp(2).E = S(3:4,3:4);
        end

end

% UPDATE OBS STRUCTURE
Obs.sid     = Sen.id ;
Obs.lid     = Lmk.id ;
Obs.ltype   = Lmk.type ;
Obs.vis     = vis ;
Obs.exp.e   = e ;
Obs.exp.E   = E ;
Obs.exp.um  = det(E);  % uncertainty measure proportional to ellipsoid area
Obs.Jac.E_r = E_rf;
Obs.Jac.E_s = E_sf;
Obs.Jac.E_l = E_l;

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% clc - clear command window.
%
% USAGE
%  c
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also CLC, CC, CCC
%
% Piotr's Image&Video Toolbox      Version 1.5
% 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]

clc

author: zouchuhang

repo: LayoutNet


text:
function [params, names] = gammaPriorExtractParam(prior)

% GAMMAPRIOREXTRACTPARAM Extract params from gamma prior structure.

% PRIOR

params = [prior.a prior.b];
if nargout > 1
  names = {'gamma a', 'gamma b'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Test file for @classicfun/plus.m

function pass = test_plus(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

singPref = pref;
singPref.blowup = true;

% Set a domain for BNDFUN.
data.domain = [-2 7];

% Generate a few random points to use as test values.
seedRNG(6178);
x = diff(data.domain) * rand(100, 1) + data.domain(1);

% A random number to use as an arbitrary additive constant.
alpha = -0.194758928283640 + 0.075474485412665i;

%%
% Check operation in the face of empty arguments.
f = bndfun();
g = bndfun(@(x) x, data, pref);
pass(1) = (isempty(f + f) && isempty(f + g) && isempty(g + f));

%%
% Check addition with scalars.
f_op = @(x) sin(x);
f = bndfun(f_op, data, pref);
pass(2:3) = test_add_function_to_scalar(f, f_op, alpha, x);

%%
% Check addition of two BNDFUN objects.
f_op = @(x) zeros(size(x));
f = bndfun(f_op, data, pref);
pass(4:5) = test_add_function_to_function(f, f_op, f, f_op, x);

f_op = @(x) exp(x) - 1;
f = bndfun(f_op, data, pref);

g_op = @(x) 1./(1 + x.^2);
g = bndfun(g_op, data, pref);
pass(6:7) = test_add_function_to_function(f, f_op, g, g_op, x);

g_op = @(x) cos(1e4*x);
g = bndfun(g_op, data, pref);
pass(8:9) = test_add_function_to_function(f, f_op, g, g_op, x);

g_op = @(t) sinh(t*exp(2*pi*1i/6));
g = bndfun(g_op, data, pref);
pass(10:11) = test_add_function_to_function(f, f_op, g, g_op, x);

%%
% Check operation for array-valued BNDFUN objects.
f_op = @(x) [zeros(size(x)) zeros(size(x)) zeros(size(x))];
f = bndfun(f_op, data, pref);
pass(12:13) = test_add_function_to_function(f, f_op, f, f_op, x);

f_op = @(x) [sin(x) cos(x) exp(x)];
f = bndfun(f_op, data, pref);
pass(14:15) = test_add_function_to_scalar(f, f_op, alpha, x);

g_op = @(x) [cosh(x) airy(1i*x) sinh(x)];
g = bndfun(g_op, data, pref);
pass(16:17) = test_add_function_to_function(f, f_op, g, g_op, x);

%%
% This should fail with a dimension mismatch error.
g_op = @(x) sin(x);
g = bndfun(g_op, data, pref);
try
    h = f + g; %#ok<NASGU>
    if ( verLessThan('matlab', '9.1') )
        pass(18) = false;
    else
        pass(18) = true;
    end
catch ME
    if ( verLessThan('matlab', '9.1') )
        pass(18) = strcmp(ME.message, 'Matrix dimensions must agree.');
    else
        pass(18) = false;
    end
end

%%
% Check that direct construction and PLUS give comparable results.
tol = 10*eps;
f = bndfun(@(x) x, data, pref);
g = bndfun(@(x) cos(x) - 1, data, pref);
h1 = f + g;
h2 = bndfun(@(x) x + cos(x) - 1, data, pref);
pass(19) = norm(feval(h1, x) - feval(h2, x), inf) < 2*tol;

%%
% Check that adding a BNDFUN to an unhappy BNDFUN gives an unhappy
% result.
f = bndfun(@(x) cos(x + 1), data);    % Happy
g = bndfun(@(x) sqrt(x + 1), data);   % Unhappy
h = f + g;  % Add unhappy to happy.
pass(20) = (~get(g, 'ishappy')) && (~get(h, 'ishappy'));
h = g + f;  % Add happy to unhappy.
pass(21) = (~get(g, 'ishappy')) && (~get(h, 'ishappy'));

%%
% Test on singular BNDFUN.
pow = -1;
op1 = @(x) (x - data.domain(2)).^pow.*sin(x);
op2 = @(x) (x - data.domain(2)).^pow.*cos(3*x);
singData = data;
singData.exponents = [0 pow];
f = bndfun(op1, singData, singPref);
g = bndfun(op2, singData, singPref);
h = f + g;
vals_h = feval(h, x);
op = @(x)  (x - data.domain(2)).^pow.*(sin(x)+cos(3*x));
h_exact = op(x);
pass(22) = ( norm(vals_h-h_exact, inf) < 1e3*max(eps, ...
    eps)*norm(h_exact, inf) );

%% Test for UNBNDFUN:

% Functions on [-inf inf]:

% Set the domain:
data.domain = [-Inf Inf];
domCheck = [-1e2 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

opf = @(x) exp(-x.^2);
opg = @(x) x.^2.*exp(-x.^2);
oph = @(x) exp(-x.^2) + x.^2.*exp(-x.^2);
f = unbndfun(opf, data);
g = unbndfun(opg, data);
h = f + g;
hVals = feval(h, x);
hExact = oph(x);
err = hVals - hExact;
pass(23) = norm(err, inf) < 1e1*eps*get(h,'vscale');

end

%%
% Test the addition of a BNDFUN F, specified by F_OP, to a scalar ALPHA using
% a grid of points X in [a  b] for testing samples.
function result = test_add_function_to_scalar(f, f_op, alpha, x)
    g1 = f + alpha;
    g2 = alpha + f;
    result(1) = isequal(g1, g2);
    g_exact = @(x) f_op(x) + alpha;
    result(2) = norm(feval(g1, x) - g_exact(x), inf) < ...
        10*max(get(g1, 'vscale')*eps);
end

%%
% Test the addition of two BNDFUN objects F and G, specified by F_OP and
% G_OP, using a grid of points X in [-1  1] for testing samples.
function result = test_add_function_to_function(f, f_op, g, g_op, x)
    h1 = f + g;
    h2 = g + f;
    result(1) = isequal(h1, h2);
    h_exact = @(x) f_op(x) + g_op(x);
    result(2) = norm(feval(h1, x) - h_exact(x), inf) <= ...
        100*max(get(h1, 'vscale')*eps);

end


text:
function res = GRAPPA(kData,kCalib,kSize,lambda, dispp)
% res = GRAPPA(kData,kCalib,kSize,lambda [, disp)
%
% This is a GRAPPA reconstruction algorithm that supports
% arbitrary Cartesian sampling. However, the implementation
% is highly inefficient in Matlab because it uses for loops.
% This implementation is very similar to the GE ARC implementation.
%
% The reconstruction looks at a neighborhood of a point and
% does a calibration according to the neighborhood to synthesize
% the missing point. This is a k-space varying interpolation.
% A sampling configuration is stored in a list, and retrieved
% when needed to accelerate the reconstruction (a bit)
%
% Inputs:
%       kData     -   [Size x, Size y, num coils] 2D multi-coil k-space data to reconstruct from.
%                   Make sure that the missing entries have exact zeros in them.
%       kCalib    -   calibration data (fully sampled k-space)
%       kSize     - size of the 2D GRAPPA kernel [kx, ky]
%       lambda    - Tykhonov regularization for the kernel calibration.
%       dispp      - Figure number to display images as they are
%                   reconstructed
% Outputs:
%       res       - k-space data where missing entries have been filled in.
%
% Example:
%   [x,y] = meshgrid(linspace(0,1,128));
%   % Generate fake Sensitivity maps
%   sMaps = cat(3,x.^2,1-x.^2,y.^2,1-y.^2);
%   % generate 4 coil phantom
%   imgs = repmat(phantom(128),[1,1,4]).*sMaps;
%   DATA = fft2c(imgs);
%   % crop 20x20 window from the center of k-space for calibration
%   kCalib = crop(DATA,[20,20,4]);
%
%   %calibrate a kernel
%   kSize = [5,5];
%   coils = 4;
%
%   % undersample by a factor of 2
%   DATA(1:2:end,2:2:end,:) = 0;
%   DATA(2:2:end,1:2:end,:) = 0;
%
%   %reconstruct:
%   [res] = GRAPPA(DATA,kCalib, kSize, 0.01);
%   figure, imshow(cat(2,sos(imgs), 2*sos(ifft2c(DATA)), sos(ifft2c(res))),[]);
%   title('full,  zero-fill,   result')
%
%
%
% (c) Michael Lustig 2008

if nargin < 5
    dispp = 0;
end

pe = size(kData,2); fe = size(kData,1); coils = size(kData,3); % get sizes

res = kData*0;
%[AtA] = corrMatrix(kCalib,kSize); % build coil correlation matrix
AtA = dat2AtA(kCalib, kSize); % build coil calibrating matrix

for n=1:coils
	disp(sprintf('reconstructing coil %d',n));
	res(:,:,n) = ARC(kData, AtA,kSize, n,lambda); % reconstruct single coil image
    if dispp ~=0
        figure(dispp), imshow3(abs(ifft2c(res)),[]); drawnow
    end

end

function [res] = ARC(kData, AtA, kSize, c,lambda);
[sx,sy,nCoil] = size(kData);

kData = zpad(kData,[sx+kSize(1)-1, sy+kSize(2)-1,nCoil]);

dummyK = zeros(kSize(1),kSize(2),nCoil); dummyK((end+1)/2,(end+1)/2,c) = 1;
idxy = find(dummyK);

res = zeros(sx,sy);

MaxListLen = 100;
LIST = zeros(kSize(1)*kSize(2)*nCoil,MaxListLen);
KEY =  zeros(kSize(1)*kSize(2)*nCoil,MaxListLen);
count = 0;

%H = waitbar(0);
for y = 1:sy
	for x=1:sx
	%	waitbar((x + (y-1)*sx)/sx/sy,H);
		tmp = kData(x:x+kSize(1)-1,y:y+kSize(2)-1,:);
		pat = abs(tmp)>0;
		if pat(idxy) | sum(pat)==0
			res(x,y) = tmp(idxy);
		else
			key = pat(:);
            idx = 0;
			for nn=1:size(KEY,2);
				if sum(key==KEY(:,nn))==length(key)
				   idx = nn;
				   break;
			   	end
			end
			if idx == 0
				count = count + 1;
				kernel = calibrate(AtA,kSize,nCoil,c,lambda,pat);
				KEY(:,mod(count,MaxListLen)+1) = key(:);
				LIST(:,mod(count,MaxListLen)+1) = kernel(:);
				%disp('add another key');size(KEY,2)
			else
				kernel = LIST(:,idx);
			end
			res(x,y) = sum(kernel(:).*tmp(:));
		end

	end
end

%close(H);

author: thomaskuestner


text:
function gdrawnow
% function gdrawnow
% ------------------------------------------------------------------------
% The |gdrawnow| function is similar to the |drawnow| command but also
% activates the |vcw| widget if present for the current figure window.
%
%
% Kevin Mattheus Moerman
%
% Change log:
% 2020/05/11 Created
% 2020/05/11 Added triple call to drawnow which appears to fix MATLAB bug
% when multiple figures are used and plotting occurs in the wrong figure
% window.
% 2020/05/11 Added new axis colorbar handling
% ------------------------------------------------------------------------
%%
% Force full drawnow
drawnow; drawnow;  drawnow; %Triple to avoid MATLAB bug

%% Activate vcw

hf=gcf;

if isfield(hf.UserData,'vcw')

    switch hf.UserData.vcw.buttonHandle.State
        case 'on'
            H=findobj(hf,'Type','colorbar'); %Handle set for all colorbars in figure;
            if ~isempty(H)
                if ~isempty(hf.UserData.vcw.colorbarHandles)
                    Hp=hf.UserData.vcw.colorbarHandles;
                    logicSetManual=~ismember(H,Hp);
                    if nnz(logicSetManual)>0
                        H_new=H(logicSetManual);
                        hf.UserData.vcw.colorbarHandles=[hf.UserData.vcw.colorbarHandles H_new];
                        hf.UserData.vcw.colorbarLocSet(end:end+numel(H_new))={H_new.Location};
                        for q=1:1:numel(H_new) %Loop over colorbar handles and set location property
                            set(H_new(q),'Location','manual');
                        end
                    end
                end
            end
        case 'off'
            hf.UserData.vcw.buttonHandle.State='on';
    end

    %Set default views for new axes and attempt to hide interactive toolbars
    h = findobj(hf, 'Type', 'axes', '-depth', 1)'; %All axis handles
    if ~isempty(h)
        for hNow = h
            %Set default views if not present
            if ~isfield(hNow.UserData,'defaultView')
                hNow.UserData.defaultView=camview(hNow);
            end

            %Try to hide toolbar
            try
                hNow.Toolbar.Visible = 'off';
            catch
            end
        end
    end

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function r = plus(p,q)
% MEAS/PLUS  Implement p + q for meas.

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

% give it the right entries
r.value = p.value + q.value;
r.error = sqrt(p.error^2 + q.error^2);

author: Sable

repo: mcbench-benchmarks


text:
function[total_cost,grad]=softmax(result,docbatch,parameter)
    total_cost=0;
    grad.soft_W=zeroMatrix(size(parameter.soft_W));
    step_size=1;
    num_word=0;
    N=size(docbatch.target_sen_matrix,1);
    zeroState=zeroMatrix([parameter.hidden,N]);
    for ll=1:parameter.layer_num
        grad.source_h{ll,1}=zeroState;
        grad.source_c{ll,1}=zeroState;
        for sen_tt=1:length(result.Target_sen)-1
            grad.target_sen_h{ll,sen_tt}=zeroState;
            grad.target_sen_c{ll,sen_tt}=zeroState;
        end
    end
    for sen_tt=1:length(result.Target_sen)
        Word_List=docbatch.target_word{sen_tt}.Word;
        Word_Mask=docbatch.target_word{sen_tt}.Mask;
        num_word=num_word+length(find(Word_Mask==1));
        N=size(Word_List,1);
        T=size(Word_List,2);
        target_sen=result.Target_sen{sen_tt};
        N_examples=size(Word_List,1)*size(Word_List,2);
        predict_Words=reshape(Word_List,1,N_examples);
        mask=reshape(Word_Mask,1,N_examples);
        h_t=[];
        if sen_tt==1
            h_t=[h_t,result.source_sen{parameter.layer_num,1}];
        else
            dim=size(result.h_t_target_sen);
            h_t=[h_t,result.h_t_target_sen{dim(1),sen_tt-1}];
        end
        dim=size(target_sen.h_t_target_word);
        h_t=[h_t,[target_sen.h_t_target_word{parameter.layer_num,1:dim(2)-1}]];
        [cost,grad_softmax_h]=batchSoftmax(h_t,mask,predict_Words,parameter);
        total_cost=total_cost+cost;
        grad.soft_W=grad.soft_W+grad_softmax_h.soft_W;
        if sen_tt==1
            grad.source_h{parameter.layer_num,1}=grad_softmax_h.h(:,1:N);
        else
            grad.target_sen_h{parameter.layer_num,sen_tt-1}=grad_softmax_h.h(:,1:N);
        end
        for i=1:T-1
            grad.ht{sen_tt}{1,i}=grad_softmax_h.h(:,N*i+1:N*(i+1));
        end
    end
    total_cost=total_cost/N;
    grad.soft_W=grad.soft_W/N;
    clear predict_Words; clear mask;
    clear grad_softmax_h;
end

function[cost,softmax_grad]=batchSoftmax(h_t,mask,predict_Words,parameter)
    unmaskedIds=find(mask==1);
    scores=parameter.soft_W*h_t;
    mx = max(scores,[],1);
    scores=bsxfun(@minus,scores,mx);
    scores=exp(scores);
    norms = sum(scores, 1);
    if length(find(mask==0))==0
        scores=bsxfun(@rdivide, scores, norms);
    else
        scores=bsxfun(@times,scores, mask./norms);
    end
    scoreIndices = sub2ind(size(scores),predict_Words(unmaskedIds),unmaskedIds);
    cost=sum(-log(scores(scoreIndices)));
    scores(scoreIndices) =scores(scoreIndices) - 1;
    softmax_grad.soft_W=scores*h_t';  %(N_word*examples)*(examples*diemsnion)=N_word*diemsnion;
    softmax_grad.h=(scores'*parameter.soft_W)';%(diemsnion*N_word)*(N_word*examples)=dimension*examples
    clear scores;
    clear norms;
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
function Yhat = LRSDL_buildYhat(Y, D, X, D_range, Y_range)
    % Yhat = [Yhat1 ... Yhatj ... YhatC]
    % Yhat_c = Y_c - Dc Xcc;
    C = numel(Y_range) - 1;
    Yhat = zeros(size(Y));
    for c = 1: C
        Yc = get_block_col(Y, c, Y_range);
        Dc = get_block_col(D, c, D_range);
        Xcc = get_block(X, c, c, D_range, Y_range);
        Yhat(:, Y_range(c) + 1: Y_range(c+1)) = Yc-Dc*Xcc;
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [RF] = trainRF_table(table,Y,categories,nBoot,cost)

[trainSets,~] = buildBootSet(Y,nBoot,'NoAdjust');

RF = applyEnsembleRF_table(table(trainSets(:,1),:),Y(trainSets(:,1)),categories,cost);
for n = 2:nBoot
   RFtemp = applyEnsembleRF_table(table(trainSets(:,n),:),Y(trainSets(:,n)),categories,cost);
   RF = append(RF,RFtemp);
end

end

author: mvallieres

repo: radiomics


text:
function UNew = fluidDirichlet2D(varargin);
% fluidDirichlet2D: solve fluid registraion in 2D with Dirichlet
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[DU,F,mu,lambda,PixSize,NumPix,HX,HY] = parse_inputs(varargin{:});

% construct filters that implement discretized Navier-Lame equations
d1 = [1;-2;1]/(PixSize(1)^2);
d2 = [1 -2 1]/(PixSize(2)^2);
d12 = [1 0 -1;0 0 0;-1 0 1]/(4*PixSize(1)*PixSize(2));

[A11,A22] = deal(zeros(3,3));
A11(:,2) = A11(:,2) + (lambda+2*mu)*d1;
A11(2,:) = A11(2,:) + mu*d2;
A22(:,2) = A22(:,2) + mu*d1;
A22(2,:) = A22(2,:) + (lambda+2*mu)*d2;

A12 = d12*(lambda+mu)/4;
A21 = A12;

% multiply force field by adjoint of Navier-Lame equations
Fnew = zeros(NumPix(1),NumPix(2),2);
Fnew(:,:,1) = imfilter(F(:,:,1),A22,'replicate') - imfilter(F(:,:,2),A12,'replicate');
Fnew(:,:,2) = imfilter(F(:,:,2),A11,'replicate') - imfilter(F(:,:,1),A21,'replicate');

% compute sine transform of new force field
FnewF1 = imag(fft(imag(fft(Fnew(:,:,1),2*NumPix(1)-2,1)),2*NumPix(2)-2,2));
FnewF2 = imag(fft(imag(fft(Fnew(:,:,2),2*NumPix(1)-2,1)),2*NumPix(2)-2,2));
FnewF1 = FnewF1(1:NumPix(1),1:NumPix(2));
FnewF2 = FnewF2(1:NumPix(1),1:NumPix(2));

% construct images of coordinates scaled by pi/(N or M)
[alpha,beta] = ndgrid(pi*(0:(NumPix(1)-1))/(NumPix(1)-1),pi*(0:(NumPix(2)-1))/(NumPix(2)-1));

% construct LHS factor
LHSfactor = mu.*(lambda+2*mu).*(2*cos(alpha) + 2*cos(beta) - 4).^2;

% set origin term to 1, as DC term does not matter
LHSfactor(1,1) = 1;

% solve for FFT of V
VF1 = FnewF1./LHSfactor;
VF2 = FnewF2./LHSfactor;

% perform inverse DST
V1 = imag(ifft(imag(ifft(VF1,2*NumPix(1)-2,1)),2*NumPix(2)-2,2));
V2 = imag(ifft(imag(ifft(VF2,2*NumPix(1)-2,1)),2*NumPix(2)-2,2));

% crop and concatenate
V = cat(3,V1(1:NumPix(1),1:NumPix(2)),V2(1:NumPix(1),1:NumPix(2)));

% construct estimate of transformation Jacobian
J = zeros(NumPix(1),NumPix(2),2,2);
J(:,:,1,1) = 1 - imfilter(V(:,:,1),HX,'replicate','same');
J(:,:,2,1) = -imfilter(V(:,:,1),HY,'replicate','same');
J(:,:,1,2) = -imfilter(V(:,:,2),HX,'replicate','same');
J(:,:,2,2) = 1 - imfilter(V(:,:,2),HY,'replicate','same');

% now perform Euler integration to construct new displacements
UNew = zeros(NumPix(1),NumPix(2),2);
UNew(:,:,1) = J(:,:,1,1).*V(:,:,1) + J(:,:,1,2).*V(:,:,2);
UNew(:,:,2) = J(:,:,2,1).*V(:,:,1) + J(:,:,2,2).*V(:,:,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DU,F,mu,lambda,PixSize,NumPix,HX,HY] = parse_inputs(varargin);

% get arguments
F = varargin{2};
PixSize = varargin{4}(1:2);
NumPix = [varargin{5} varargin{6}];
mu = varargin{8};
lambda = varargin{9};
DU = varargin{11};
HX = varargin{12};
HY = varargin{13};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
 function y = mtimes_block(ob, x, iblock, nblock)
%function y = mtimes_block(ob, x, iblock, nblock)
% y = G(i'th block) * x	or y = G'(i'th block) * x
% in either case the project data will be "small"
% iblock is 1,...,nblock

% support 'exists' option for seeing if this routine is available
if nargin == 2 & ischar(x) & streq(x, 'exists')
	y = 1;
	return
end

if nargin ~= 4, ir_usage, end

nb = ob.nb;
na = ob.na;

%
% forward projection
%
if ~ob.is_transpose
	if ob.apower ~= 1, error notdone, end

	% if needed, expand concise column
	if ob.is_masked
		idim = size(x);
		np = sum(ob.mask(:));
		nxy = numel(ob.mask);
		if idim(1) ~= np
			error 'size mismatch'
		end
		x = embed(x, ob.mask);
		x = reshape(x, nxy, idim(2));
	end

	y = double(wtfmex('chat', ob.chat, 'proj,block', single(x), ...
		int32(iblock-1), int32(nblock)));

	% fix: extract the relevant columns - should do in wtfmex?
	ia = iblock:nblock:na;
	nv = length(ia);

	if ndims(y) == 3			% [nb,na,nz]
		y = y(:,ia,:);			% [nb,nv,nz]

	elseif size(y,1) == nb*na		% [nba,nz]
		nz = size(y,2);
		y = reshape(y, nb, na, nz);	% [nb,na,nz]
		y = y(:,ia,:);			% [nb,nv,nz]
		y = reshape(y, nb*nv, nz);	% [nb*nv,nz]

	elseif size(y,1) == nb			% [nb,na]
		y = y(:,ia);			% [nb,nv]
	else
		error size
	end

%
% backprojection
%
else
	if ob.apower ~= 1, error notdone, end

	ia = iblock:nblock:na;
	nv = length(ia);

	if ndims(x) == 3		% [nb,nv,nz]
		error todo

	elseif size(x,1) == nb*nv	% [nb*nv,nz]
		nz = size(x,2);
		tmp = zeros(nb, na, nz);
		tmp(:,ia,:) = reshape(x, [nb nv nz]);
		x = reshape(tmp, [nb*na nz]);

	elseif size(x,1) == nb		% [nb,nv]
		error todo

	else
		error bug
	end

	y = double(wtfmex('chat', ob.chat, 'back,block', single(x), ...
		int32(iblock-1), int32(nblock)));

	if ob.is_masked
		if size(y,1) == numel(ob.mask)
			y = y(ob.mask,:);	% [nxy,nz] -> [np,nz]
		else
			y = y(ob.mask);		% [nx,ny] -> [np]
		end
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function cheby_u_poly_test ( )

%*****************************************************************************80
%
%% CHEBY_U_POLY_TEST tests CHEBY_U_POLY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2015
%
%  Author:
%
%    John Burkardt
%
  n_max = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_U_POLY_TEST:\n' );
  fprintf ( 1, '  CHEBY_U_POLY evaluates the Chebyshev T polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N      X        Exact F       U(N)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = cheby_u_poly_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = cheby_u_poly ( 1, n, x );

    fprintf ( 1, '  %6d  %8f  %12f  %12f\n', n, x, fx, fx2(n+1) );

  end

  return
end

author: johannesgerer


text:
function points = robotlaser_as_cartesian(rl, maxRange = 15, subsample = false)

numBeams = length(rl.ranges);
maxRange=min(maxRange, rl.maximum_range);
% apply the max range
idx = rl.ranges<maxRange & rl.ranges>0;

if (subsample)
	idx(2:2:end) = 0;
endif

angles = linspace(rl.start_angle, rl.start_angle + numBeams*rl.angular_resolution, numBeams)(idx);
points = [rl.ranges(idx) .* cos(angles); rl.ranges(idx) .* sin(angles); ones(1, length(angles))];
transf = v2t(rl.laser_offset);

% apply the laser offset
points = transf * points;

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


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
    %   logisitc 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
    n = size(X, 1);
    d = size(X, 2);

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

    % Add ones to the X data matrix
    X = [ones(n, 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 use
    %       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);
    %

    % for optimization of theta
    options = optimset('GradObj', 'on', 'MaxIter', 50);

    % for each class label
    for k = 1 : num_labels

        % initialize theta, and minimize value for cost function
        init_theta = zeros(d + 1, 1);
        theta = fmincg(@(t)(lrCostFunction(t, X, (y == k), lambda)), ...
            init_theta, options);
        all_theta(k, :) = theta';
    end
end

author: worldveil

repo: coursera-ml


text:
function value = epn_glg_monomial_integral ( n, expon, alpha )

%*****************************************************************************80
%
%% EPN_GLG_MONOMIAL_INTEGRAL: integral of monomial with GLG weight on EPN.
%
%  Discussion:
%
%    EPN_GLG is the N-dimensional positive space [0,+oo)^N with generalized
%    Laguerre weight function:
%
%      w(alpha;x) = product ( 1 <= i <= n ) x(i)^alpha exp ( - x(i) )
%
%    value = integral ( EPN )
%      product ( 1 <= i <= n ) x(I)^expon(i) x(i)^alpha exp ( - x(i) ) dx(i)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, integer EXPON(N), the exponents.
%
%    Input, real ALPHA, the exponent of X in the weight function.
%    -1.0 < ALPHA.
%
%    Output, real VALUE, the value of the integral.
%
  value = 1.0;
  for i = 1 : n
    value2 = ep1_glg_monomial_integral ( expon(i), alpha );
    value = value * value2;
  end

  return
end


text:
function [ bx, ndim ] = cvx_sparsify( bx, tt, mode )

global cvx___
ox = isa( bx, 'cvx' );
if ox,
    sx = size( bx );
    bx = cvx_basis( bx );
end
if isempty( tt ),
    by = bx;
elseif ~any( tt ),
    return
else
    by = bx( :, tt );
end
if isempty( mode ),
    by = bx;
else
    [ xR, by ] = cvx_bcompress( by, mode, 0 );
end
nx = size( by, 2 );
ndim = cvx_newvar( nx, cvx_classify_mex( by, cvx___.classes ) );
bz = sparse( ndim, 1 : nx, 1 );
nz = size(bz,1);
by( nz, end ) = 0;
cvx_newcnstr( by - bz, true );
if ~isempty( mode ),
    bz = bz * xR;
end
if isempty( tt ),
    bx = bz;
else
    bx( 1:nz, tt ) = bz;
end
if ox,
    bx = cvx( sx, bx );
end

% 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


text:
% This m-file takes care of synthetic division.
% By giving one polynomial and one root this function returns
% the polynomial formed with the other roots of the given polynomial excluding the given root.
% Keerthi Venkateswara Rao
function coeff_second=syn_division(coeff_function,fun_root_new)
order_fun=size((coeff_function),2);
coeff_second=0;
for index=1:size((coeff_function),2)-1
    if index==1
        coeff_second(index)=coeff_function(index);
    else
        coeff_second(index)=coeff_function(index)+fun_root_new*coeff_second(index-1);
    end
end


text:
function [ x, y ] = i4_swap ( x, y )

%*****************************************************************************80
%
%% I4_SWAP swaps two I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, Y, two values to interchange.
%
%    Output, integer X, Y, the interchanged values.
%
  z = x;
  x = y;
  y = z;

  return
end

repo: jburkardt-m


text:
function im_opp = rgb_to_opp(im)
% Transforms image into opponent color space
% Code by Ross Girshick

if ~isa(im, 'double')
  im = double(im);
end

im_opp = zeros(size(im));
% Color saliency boosting? c1*0.850, c2*0.524, c3*0.065
im_opp(:,:,1) = (im(:,:,1) - im(:,:,2)) / sqrt(2);
im_opp(:,:,2) = (im(:,:,1) + im(:,:,2) - 2*im(:,:,3)) / sqrt(6);
im_opp(:,:,3) = sum(im,3) / sqrt(3);
%min(im_opp(:)) % about -50
%max(im_opp(:)) % about 450
im_opp = max(0, min(255, im_opp)); % is this a good idea?

author: Cloud-CV

repo: object-proposals


text:
function [stiff, diinsy, cols, sysmat] = sb_calc_stiff(vol)

% SB_CALC_STIFF
%
% $Id$

if(~(size(vol.pos,2)==3))
    if(size(vol.pos,1)==3)
        node = vol.pos';
        warning('Dimension of vol.pos should be #nodes x 3!')
    else
        error('vol.pos has wrong dimension!')
    end
else
    node = vol.pos;
end
npnt = size(node,1);
npnt = int32(npnt);

if isfield(vol,'tet')
    if size(vol.tet,1) == 4
        mele = size(vol.tet,1);
        elem = vol.tet;
    elseif size(vol.tet,2) == 4
        mele = size(vol.tet,2);
        elem = vol.tet';
    else
        error('vol.tet has wrong dimensions!')
    end
    elem = [elem; zeros(4,size(elem,2))];
elseif isfield(vol,'hex')
    if size(vol.hex,1) == 8
        mele = size(vol.hex,1);
        elem = vol.hex;
    elseif size(vol.hex,2) == 8
        mele = size(vol.hex,2);
        elem = vol.hex';
    else
        error('vol.hex has wrong dimensions!')
    end
else
    error('Could not find connectivity information!')
end

if min(min(elem(1:mele,:))) == 0
    elem = elem + 1;
    warning('Numbering of nodes in vol.tet/vol.hex must start at 1 (Fortran numbering)!')
elseif min(min(elem(1:mele,:))) < 0
    error('No negative indices for conectivity information allowed!')
end

if isfield(vol,'cond') && isfield(vol,'tissue') && isfield(vol,'tissuelabel')
    if length(vol.tissuelabel) == length(vol.cond)
         if length(vol.tissue) == size(elem,2)
             cond = zeros(size(elem,2),1);
             numlabels = length(vol.tissuelabel);
             for i=1:numlabels
                 cond(vol.tissue == i) = vol.cond(i);
             end
        else
            error('Dimensions of vol.tet or vol.hex and vol.tissue do not fit!');
        end
    else
        error('Dimensions of vol.cond and entries of vol.tissuelabel do not fit!');
    end
end

mele = int32(mele);
elem = int32(elem);

% check whether the nodes have right orientation

if isfield(vol,'tet')
  if ~sb_test_ori(node,elem(1:4,:)')
    error('Elements have wrong orientation, consider exchanging node 3 and 4');
  end
elseif isfield(vol,'hex')
  if ~sb_test_ori(node,elem')
    error('Elements have wrong orientation or are degenerated');
  end
end

try
  [diinsy,cols,sysmat] = calc_stiff_matrix_val(node,elem,cond,mele);
catch err
  if ispc && strcmp(err.identifier,'MATLAB:invalidMEXFile')
    error('Error executing mex-file. Microsoft Visual C++ 2008 Redistributables and Intel Visual Fortran Redistributables are required.')
  else
    rethrow(err)
  end
end
npnt   = double(npnt);
diinsy = double(diinsy);
cols   = double(cols);
rows   = sb_sparse_to_mat(diinsy);
stiff  = sparse(rows,cols,sysmat,npnt,npnt,length(sysmat));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function args = hlp_struct2varargin(struc,varargin)
% Convert a struct into a sequence of name-value pairs; inverse to hlp_varargin2struct.
% Args = hlp_struct2varargin(Struct,Options...)
%
% In:
%   Struct: a 1x1 structure array
%   Options: optional name-value pairs;
%            'suppress': suppress the output names listed in the following cell array
%            'rewrite' : rewrite some names in the struct, if present (executed after suppress and before restrict)
%                        specified as a cell array of {oldname,newname,oldname,newname, ....}
%            'restrict': restrict the output names to those listed in the following cell array
%
% Out:
%  Args: a list of name-value pairs
%
% Examples:
%   % in the given cell array of name-value pairs, remove occurrences of myarg1 and myarg2 parameters
%   hlp_struct2varargin(options,'suppress',{'myarg1','myarg2'});
%
%   % rewrite any occurrences of a 'channel' parameter into 'chns', and likewise for samplerate
%   hlp_struct2varargin(options,'rewrite',{'channels','chns','samplerate','srate'});
%
%   % restrict the given name-value pairs to those whose name is either test, arg1, or arg2.
%   hlp_struct2varargin(options,'restrict',{'test','arg1','arg2'});
%
% See also:
%   hlp_varargin2struct
%
%                               Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                               2010-03-28

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

if isempty(struc)
    args = {};
    return;
end

opts = hlp_varargin2struct(varargin,'restrict',[],'suppress',[],'rewrite',[]);

fields = fieldnames(struc)';
data = struct2cell(struc)';
% suppress some of the original field names
if ~isempty(opts.suppress)
    [fields,I] = setdiff(fields,opts.suppress);
    data = data(I);
end
% rewrite some of the original field names into new field names
for c=1:2:length(opts.rewrite)
    fields(strcmp(fields,opts.rewrite{c})) = opts.rewrite(c+1); end
% restrict to a subset of old/new field names
if ~isempty(opts.restrict)
    [fields,I,J] = intersect(fields,opts.restrict);  %#ok<NASGU>
    data = data(I);
end

args = vertcat(fields,data);
args = args(:)';

repo: REST


text:
function y = sinc(x, c)
%SINC   The function SIN(PI*X)/(PI*X).
%
%   SINC(X) returns the SIN(PI*X)/(PI*X) function which is defined as
%
%      1                 if x = 0
%      sin(pi*x)/(pi*x)  if 0 < |x| < infinity
%      0                 if |x| = infinity
%
%   SINC(X, C) returns the SIN(C*X)/(C*X) function, where C > 0, which is a
%   generalization of the above and is defined as
%
%      1                 if x = 0
%      sin(c*x)/(c*x)    if 0 < |x| < infinity
%      0                 if |x| = infinity
%
%   See also SIN.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-20 08:45:14 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   nargsin = nargin;

   % check number of input arguments
   error(nargchk(1, 2, narsgin));

   if nargsin < 2
      c = pi;
   else
      if any(size(c) ~= 0) | (c <= 0)
         error('Second argument must be a positive scalar.');
      end
   end

   y = ones(size(x));           % initialize output
   i = x ~= 0;                  % find non-zero elements
   t = c * x(i);                % precompute C*X
   y(i) = sin(t) ./ t;          % compute SIN(C*X)/(C*X)
   y(isinf(x)) = 0;             % 0 not NaN when X = +/-Inf

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function targets=FSEOF(model,biomassRxn,targetRxn,iterations,coefficient,outputFile)
% FSEOF: implements the algorithm of Flux Scanning based on Enforced Objective Flux.
%
%   model           a model structure
%   biomassRxn      string with reaction ID of the biomass formation or
%                   growth reaction
%   targetRxn       string with reaction ID of target reaction
%   iterations      double indicating number of iterations (opt, default 10)
%   coefficient     double indicating ratio of optimal target reaction
%                   flux, must be less than 1 (opt, default 0.9)
%   outputFile      string with output filename (opt, default prints to
%                   command window)
%
%   targets         structure with target identifying information for each reaction
%       logical     logical array indicating FSEOF identified target reaction
%       slope       double array with FSEOF calculated slopes for each reaction

%OUTPUTS
%   This function writes an tab-delimited file or prints to command window. If an output
%   has been specified (targets), it will also generate a structure indicating for
%   each reaction whether it is identified by FSEOF as a target and the slope of the
%   reaction when switching from biomass formation to product formation.
%
%   Usage: targets=FSEOF(model,biomassRxn,targetRxn,iterations,coefficient,outputFile)

biomassRxn=char(biomassRxn);
targetRxn=char(targetRxn);

if nargin<4
    iterations=10;
    coefficient=0.9;
end

if nargin <5
    coefficient=0.9;
end

if nargin == 6
    output=1;
else
    output=0;
end

%Find out the maximum theoretical yield of target reaction
model=setParam(model,'obj',targetRxn,1);
sol=solveLP(model,1);
targetMax=abs(sol.f*coefficient);   % 90 percent of the theoretical yield

model=setParam(model,'obj',biomassRxn,1);

fseof.results=zeros(length(model.rxns),iterations);
fseof.target=zeros(length(model.rxns),1);
rxnDirection=zeros(length(model.rxns),1);

%Enforce objective flux iteratively
for i=1:iterations
    n=i*targetMax/iterations;
    model=setParam(model,'lb',targetRxn,n);

    sol=solveLP(model,1);

    fseof.results(:,i)=sol.x;

    %Loop through all fluxes and identify the ones that increased upon the
    %enforced objective flux
    for j=1:length(fseof.results)
        if fseof.results(j,1) > 0   %Check the positive fluxes

            if i == 1   %The initial round
                rxnDirection(j,1)=1;
                fseof.target(j,1)=1;
            else

                if (fseof.results(j,i) > fseof.results(j,i-1)) & fseof.target(j,1)
                    fseof.target(j,1)=1;
                else
                    fseof.target(j,1)=0;
                end
            end

        elseif fseof.results(j,1) < 0 %Check the negative fluxes

            if i == 1   %The initial round
                rxnDirection(j,1)=-1;
                fseof.target(j,1)=1;
            else
                if (fseof.results(j,i) < fseof.results(j,i-1)) & fseof.target(j,1)
                    fseof.target(j,1)=1;
                else
                    fseof.target(j,1)=0;
                end
            end

        end

    end
end

%Generating output
formatSpec='%s\t%s\t%s\t%s\t%s\t%s\t%s\n';
if output == 1    %Output to a file
    outputFile=char(outputFile);
    fid=fopen(outputFile,'w');
    fprintf(fid,formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystems','Direction','Gr Rule');
else              %Output to screen
    fprintf(formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystems','Direction','Gr Rule');
end

for num=1:length(fseof.target)
    if fseof.target(num,1) == 1
        A0=num2str(abs(fseof.results(num,iterations)-fseof.results(num,1))/abs(targetMax-targetMax/iterations)); %Slope calculation
        A1=num2str(num);                                  %row ID
        A2=char(model.rxns(num));                         %enzyme ID
        A3=char(model.rxnNames(num));                     %enzyme Name
        if isfield(model,'subSystems') && ~isempty(model.subSystems{num});
            A4=char(strjoin(model.subSystems{num,1},';'));                   %Subsystems
        else
            A4='';
        end
        A5=num2str(model.rev(num)*rxnDirection(num,1));   %reaction Dirction
        A6=char(model.grRules(num));                      %Gr Rule
        if output == 1    %Output to a file
            fprintf(fid,formatSpec,A0,A1,A2,A3,A4,A5,A6);
        else              %Output screen
            fprintf(formatSpec,A0,A1,A2,A3,A4,A5,A6);
        end
    end
end

if output == 1    %Output to a file
    fclose(fid);
end

if nargout == 1
    targets.logical=logical(fseof.target);
    targets.slope=abs(fseof.results(:,iterations)-fseof.results(:,1))/abs(targetMax-targetMax/iterations); %Slope calculation
end
end

author: SysBioChalmers

repo: RAVEN


text:
function jd = julday(y,m,d,h)
% JULDAY  Conversion of date as given by
%         y ... year (four digits)
%         m ... month
%         d ... day
%         h ... hour and fraction hereof
%         The conversion is only valid in the time span
%         from March 1, 1900 to February 28, 2100

%  For further information see
%  Meeus, Jean (1991) Astronomical Algorithms,
%         Willmann-Bell, Richmond, Virginia, p. 59--62

%  Written by Kai Borre
%  February 14, 2001

      if m <= 2, y = y-1; m = m+12; end
      jd = floor(365.25*(y+4716))+floor(30.6001*(m+1))+d+h/24-1537.5;
%      mjd = jd-2400000.5;
%%%%%%% end julday.m  %%%%%%%%%%%%

author: yandld


text:
% Equality constrained norm minimization.
%
% This script constructs a random equality-constrained norm minimization
% problem and solves it using CVX. You can also change p to +2 or +Inf
% to produce different results. Alternatively, you an replace
%     norm( A * x - b, p )
% with
%     norm_largest( A * x - b, 'largest', p )
% for 1 <= p <= 2 * n.

% Generate data
p = 1;
n = 10; m = 2*n; q=0.5*n;
A = randn(m,n);
b = randn(m,1);
C = randn(q,n);
d = randn(q,1);

% Create and solve problem
cvx_begin
   variable x(n)
   dual variable y
   minimize( norm( A * x - b, p ) )
   subject to
        y : C * x == d;
cvx_end

% Display results
disp( sprintf( 'norm(A*x-b,%g):', p ) );
disp( [ '   ans   =   ', sprintf( '%7.4f', norm(A*x-b,p) ) ] );
disp( 'Optimal vector:' );
disp( [ '   x     = [ ', sprintf( '%7.4f ', x ), ']' ] );
disp( 'Residual vector:' );
disp( [ '   A*x-b = [ ', sprintf( '%7.4f ', A*x-b ), ']' ] );
disp( 'Equality constraints:' );
disp( [ '   C*x   = [ ', sprintf( '%7.4f ', C*x ), ']' ] );
disp( [ '   d     = [ ', sprintf( '%7.4f ', d   ), ']' ] );
disp( 'Lagrange multiplier for C*x==d:' );
disp( [ '   y     = [ ', sprintf( '%7.4f ', y ), ']' ] );

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% WSUM - Weighted sum of a matrix along one dimension.
%
% Y = WSUM(X,W,DIM)
%
% X   : matrix
% W   : vector of weights
% DIM : dimension to sum along (default: first non-singleton dimension)
%
% When X is a 2D-matrix, one can simply use W'*X or X*W for more
% efficient performace.

% Last modified 2010-06-09
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function y = wsum(X,w,dim)

% warning('This function seems to be inefficient..')

if ~isvector(w)
  error('Weights must be given as a vector.');
end

if nargin < 3
  % Default dimension: first non-singleton dimension
  dim = find(size(X)>1, 1);
  if isempty(dim)
    dim = 1;
  end
end

if size(X,dim) ~= length(w)
  error(['The length of the weight vector does not match with the size ' ...
         'of the matrix.']);
end

s = ones(1,max(dim,2));
s(dim) = length(w);
y = sum(bsxfun(@times, X, reshape(w,s)), dim);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function status=odetpbar(t,y,flag)
    persistent tf tstart;

    if isempty(flag)
        % Integration steps
        ts=mean(t);
        progress=100*ts/tf;
        textprogressbar(progress);
        status=0;
    else
        switch flag
            case 'init'     % Initializing progress bar
                tstart=tic;
                tf=max(t);
                textprogressbar('ODE integration: ');
            case 'done'     % Finishing status function
                tf=[];
                textprogressbar('');
                display([ '   Integration time: ' num2str(toc(tstart))]);
                tstart=[];
            otherwise
                error('odetpbar:UnknownError',...
                    'Unknown error has occured');
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function res = matricize( U, mode )
    %MATRICIZE Matricize 3D Matlab array.
    %   A = MATRICIZE(U, MODE) matricizes the 3D Matlab array U along the
    %   specified mode MODE. Higher dimensions than 3 are not supported.
    %
    %   See also TENSORIZE, TENSORPROD_TTEMPS, UNFOLD.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    d = size(U);
    % pad with 1 for the last dim (annoying)
    if length(d) == 2
        d = [d, 1];
    end

    switch mode
        case 1
            res = reshape( U, [d(1), d(2)*d(3)] );
        case 2
            res = reshape( permute( U, [2, 1, 3]), [d(2), d(1)*d(3)] );
        case 3
            res = transpose( reshape( U, [d(1)*d(2), d(3)] ) );
        otherwise
            error('Invalid mode input in function matricize')
    end
end

author: NicolasBoumal

repo: manopt