text:
function [pairwiseInteractions, pairwiseSolutions] = simulatePairwiseInteractions(pairwiseModelFolder, varargin)
% This function predicts the outcome of pairwise simulations in every
% combination from a given list of pairwise models. The pairwise models
% need to be created first with the function joinModelsPairwiseFromList.
% This script requires the COBRA Toolbox function solveCobraLP. Due to the
% coupling constraints on the pairwise models, the simulations cannot
% currently be run with optimizeCbModel.
%
% Below is a description of all possible consequences for the two joined
% organisms. Please note that the outcomes depend on the two genome-scale
% reconstructions joined and are highly dependent on the applied
% constraints.
%
% * Competition: both organisms grow slower in co-growth than separately
%   (same outcome for both).
%
% * Parasitism: one organism grows faster in co-growth than separately, while
%   the other grows slower in co-growth than separately. Outcome for
%   faster-growing organism: Parasitism_Taker, for slower-growing organism:
%   Parasitism_Giver
%
% * Amensalism: one organism's growth is unaffected by co-growth,  while
%   the other grows slower in co-growth than separately. Outcome for
%   unaffected organism: Amensalism_Unaffected, for slower-growing organism:
%   Amensalism_Affected
%
% * Neutralism: both organisms' growths are unaffected by co-growth (same
%   outcome for both)
%
% * Commensalism: one organism's growth is unaffected by co-growth,  while
%   the other grows fatser in co-growth than separately. Outcome for
%   unaffected organism: Commensalism_Giver, for slower-growing organism:
%   Commensalism_Taker
%
% * Mutualism: both organisms growth faster in co-growth than separately
%   (same outcome for both)
%
% Please note: the function takes the name of the folder containing
% previously created pairwise models as the input, e.g.:
% pairwiseInteractions, pairwiseSolutions] =
% simulatePairwiseInteractions('/Users/almut.heinken/Documents/PairwiseModels');
%
% USAGE:
%     [pairwiseInteractions, pairwiseSolutions] = simulatePairwiseInteractions(pairwiseModelFolder, varargin)
%
% INPUTS:
%     pairwiseModelFolder:     Folder where pairwise models and pairwise model info file are located
%
% OPTIONAL INPUTS:
%     inputDiet:               Cell array of strings with three columns containing
%                              exchange reaction model abbreviations, lower bounds,
%                              and upper bounds.
%                              If no diet is input then the input pairwise models
%                              will be used with unchanged constraints.
%     saveSolutionsFlag:       If true, flux solutions are stored (may result in
%                              large output files)
%     numWorkers:              Number of workers in parallel pool if desired
%     sigD:                    Difference in growth rate that counts as significant
%                              change (by default: 10%)
%
% OUTPUTS:
%     pairwiseInteractions:    Table with computed pairwise and single growth
%                              rates for all entered microbe-microbe models
%
% OPTIONAL OUTPUT:
%     pairwiseSolutions:       Table with all computed flux solutions saved
%                              (may result in large output files)
%
% .. Authors:
%       - Almut Heinken, 02/2018
%       - Almut Heinken, 02/2020: Inputs changed for more efficient computation

parser = inputParser();  % Parse input parameters
parser.addRequired('pairwiseModelFolder', @ischar);
parser.addParameter('inputDiet', {}, @iscell)
parser.addParameter('saveSolutionsFlag', false, @(x) isnumeric(x) || islogical(x))
parser.addParameter('numWorkers', 0, @(x) isnumeric(x))
parser.addParameter('sigD', 0.1, @(x) isnumeric(x))

parser.parse(pairwiseModelFolder, varargin{:});

pairwiseModelFolder = parser.Results.pairwiseModelFolder;
inputDiet = parser.Results.inputDiet;
saveSolutionsFlag = parser.Results.saveSolutionsFlag;
numWorkers = parser.Results.numWorkers;
sigD = parser.Results.sigD;

%% initialize COBRA Toolbox and parallel pool
global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end
solver = CBT_LP_SOLVER;

if numWorkers > 0
    % with parallelization
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        parpool(numWorkers)
    end
end
environment = getEnvironment();

pairwiseSolutions = {};
if saveSolutionsFlag == true
    pairwiseSolutions{1, 1} = 'pairedModelID';
    pairwiseSolutions{1, 2} = 'PairwiseSolution';
    pairwiseSolutions{1, 3} = 'SingleModel1Solution';
    pairwiseSolutions{1, 4} = 'SingleModel2Solution';
end

pairwiseInteractions{1, 1} = 'pairedModelID';
pairwiseInteractions{1, 2} = 'ModelID1';
pairwiseInteractions{1, 3} = 'ModelID2';
pairwiseInteractions{1, 4} = 'pairedGrowth_Model1';
pairwiseInteractions{1, 5} = 'pairedGrowth_Model2';
pairwiseInteractions{1, 6} = 'singleGrowth_Model1';
pairwiseInteractions{1, 7} = 'singleGrowth_Model2';
pairwiseInteractions{1, 8} = 'Outcome_Model1';
pairwiseInteractions{1, 9} = 'Outcome_Model2';
pairwiseInteractions{1, 10} = 'Total_Outcome';

growthRates={};
solutionsTmp={};

if isfile([pairwiseModelFolder filesep 'pairedModelInfo.mat'])
load([pairwiseModelFolder filesep 'pairedModelInfo.mat'],'pairedModelInfo')
else
    error('Pairwise models have not been created. Please run joinModelsPairwiseFromList first.')
end

info=pairedModelInfo;
parfor i = 1:size(info, 1)
    restoreEnvironment(environment);
    changeCobraSolver(solver, 'LP', 0, -1);
    changeCobraSolverParams('LP', 'logFile', 0);

    % load the model
    filename=[pairwiseModelFolder filesep info{i,1}];
    pairedModel=readCbModel(filename);
    pairedModelOrg=pairedModel;
    % if a diet was input
    if ~isempty(inputDiet)
        pairedModel = useDiet(pairedModel, inputDiet);
    end
    % for each paired model, set both biomass objective functions as
    % objectives
    biomass1 = strcat(info{i, 2}, '_', info{i, 3});
    biomass2 = strcat(info{i, 4}, '_', info{i, 5});
    model1biomass = find(ismember(pairedModel.rxns, biomass1));
    pairedModel.c(model1biomass, 1) = 1;
    model2biomass = find(ismember(pairedModel.rxns, biomass2));
    pairedModel.c(model2biomass, 1) = 1;
    % enforce minimal growth
    pairedModel = changeRxnBounds(pairedModel, biomass1, 0.001, 'l');
    pairedModel = changeRxnBounds(pairedModel, biomass2, 0.001, 'l');
    % calculate joint biomass
    solutionPaired = solveCobraLP(buildLPproblemFromModel(pairedModel,false));
    % separate growth
    % silence model 2 and optimize model 1
    % load the model again to avoid errors
    pairedModel=pairedModelOrg;
    % if a diet was input
    if ~isempty(inputDiet)
        pairedModel = useDiet(pairedModel, inputDiet);
    end
    pairedModel = changeObjective(pairedModel, biomass1);
    % disable flux through the second model
    pairedModel = changeRxnBounds(pairedModel, pairedModel.rxns(strmatch(strcat(info{i, 4}, '_'), pairedModel.rxns)), 0, 'b');
    % calculate single biomass
    solutionSingle1 = solveCobraLP(buildLPproblemFromModel(pairedModel,false));
    % silence model 1 and optimize model 2
    % load the model again to avoid errors
    pairedModel=pairedModelOrg;
    % if a diet was input
    if ~isempty(inputDiet)
        pairedModel = useDiet(pairedModel, inputDiet);
    end
    pairedModel = changeObjective(pairedModel, biomass2);
    % disable flux through the first model
    pairedModel = changeRxnBounds(pairedModel, pairedModel.rxns(strmatch(strcat(info{i, 2}, '_'), pairedModel.rxns)), 0, 'b');
    % calculate single biomass
    solutionSingle2 = solveCobraLP(buildLPproblemFromModel(pairedModel,false));

    % save results temporarily
    if solutionPaired.stat==1
        growthRates{i}{1,1} = solutionPaired.full(model1biomass);
        growthRates{i}{1,2} = solutionPaired.full(model2biomass);
    else
        growthRates{i}{1,1} = 0;
        growthRates{i}{1,2} = 0;
    end
    if solutionSingle1.stat==1
        growthRates{i}{1,3} = solutionSingle1.full(model1biomass);
    else
        growthRates{i}{1,3} = 0;
    end
    if solutionSingle2.stat==1
        growthRates{i}{1,4} = solutionSingle2.full(model2biomass);
    else
        growthRates{i}{1,4} = 0;
    end

    solutionsTmp{i}{1,1} = solutionPaired;
    solutionsTmp{i}{1,2} = solutionSingle1;
    solutionsTmp{i}{1,3} = solutionSingle2;
end

% backup results for large number of simulations
if size(pairedModelInfo, 1) > 1000
save('growthRates','growthRates','-v7.3');
save('solutionsTmp','solutionsTmp','-v7.3');
end

for i = 1:size(pairedModelInfo, 1)
    solPaired = solutionsTmp{i}{1,1};
    solSingle1 = solutionsTmp{i}{1,2};
    solSingle2 = solutionsTmp{i}{1,3};

    grPaired1 = growthRates{i}{1,1};
    grPaired2 = growthRates{i}{1,2};
    grSingle1 = growthRates{i}{1,3};
    grSingle2 = growthRates{i}{1,4};

    % interpret computed growth rates-only if solutions are feasible
    if solPaired.stat ~= 0 && solSingle1.stat ~= 0 && solSingle2.stat ~= 0
        [iAFirstModel, iASecondModel, iATotal] = analyzePairwiseInteractions(grPaired1, grPaired2, grSingle1, grSingle2, sigD);
    end
    %% list all results
    % fill out the results table with the model information
    pairwiseInteractions{i + 1, 1} = pairedModelInfo{i, 1};
    pairwiseInteractions{i + 1, 2} = pairedModelInfo{i, 2};
    pairwiseInteractions{i + 1, 3} = pairedModelInfo{i, 4};
    % combined growth
    pairwiseInteractions{i + 1, 4} = grPaired1;
    pairwiseInteractions{i + 1, 5} = grPaired2;
    pairwiseInteractions{i + 1, 6} = grSingle1;
    pairwiseInteractions{i + 1, 7} = grSingle2;
    % make sure that infeasible solutions are flagged
    if solPaired.stat == 0 || solSingle1.stat == 0 || solSingle2.stat == 0
        pairwiseInteractions{i + 1, 8} = 'Infeasible';
        pairwiseInteractions{i + 1, 9} = 'Infeasible';
        pairwiseInteractions{i + 1, 10} = 'Infeasible';
    else
        % save the results
        pairwiseInteractions{i + 1, 8} = iAFirstModel;
        pairwiseInteractions{i + 1, 9} = iASecondModel;
        pairwiseInteractions{i + 1, 10} = iATotal;
    end
    %% store the solutions if storeSolutionFlag==true
    if saveSolutionsFlag == true
        pairwiseSolutions{i + 1, 1} = pairedModelInfo{i, 1};
        pairwiseSolutions{i + 1, 2} = solPaired.full;
        pairwiseSolutions{i + 1, 3} = solSingle1.full;
        pairwiseSolutions{i + 1, 4} = solSingle2.full;
    end
end

end

function [iAFirstOrg, iASecondOrg, iATotal] = analyzePairwiseInteractions(pairedGrowthOrg1, pairedGrowthOrg2, singleGrowthOrg1, singleGrowthOrg2, sigD)
% This function evaluates the outcome of pairwise growth of two organisms
% compared with the two organisms separately. There are six possible
% outcomes of the interaction, and nine possible outcomes from the
% perspective of each organism.

if abs(1 - (pairedGrowthOrg1 / singleGrowthOrg1)) < sigD
    % first microbe unaffected - all possible cases resulting froms
    % second microbe's growth
    if abs(1 - (pairedGrowthOrg2 / singleGrowthOrg2)) < sigD
        % second microbe unaffected
        iAFirstOrg = 'Neutralism';
        iASecondOrg = 'Neutralism';
        iATotal = 'Neutralism';
    elseif abs((pairedGrowthOrg2 / singleGrowthOrg2)) > 1 + sigD
        % second microbe grows better
        iAFirstOrg = 'Commensalism_Giver';
        iASecondOrg = 'Commensalism_Taker';
        iATotal = 'Commensalism';
    elseif abs((singleGrowthOrg2 / pairedGrowthOrg2)) > 1 + sigD
        % second microbe grows slower
        iAFirstOrg = 'Amensalism_Unaffected';
        iASecondOrg = 'Amensalism_Affected';
        iATotal = 'Amensalism';
    else
        % if no case fits - needs inspection!
        iAFirstOrg = 'No_Result';
        iASecondOrg = 'No_Result';
        iATotal = 'No_Result';
    end
elseif abs((pairedGrowthOrg1 / singleGrowthOrg1)) > 1 + sigD
    % first microbe grows better - all possible cases resulting froms
    % second microbe's growth
    if abs(1 - (pairedGrowthOrg2 / singleGrowthOrg2)) < sigD
        % second microbe unaffected
        iAFirstOrg = 'Commensalism_Taker';
        iASecondOrg = 'Commensalism_Giver';
        iATotal = 'Commensalism';
    elseif abs((pairedGrowthOrg2 / singleGrowthOrg2)) > 1 + sigD
        % second microbe grows better
        iAFirstOrg = 'Mutualism';
        iASecondOrg = 'Mutualism';
        iATotal = 'Mutualism';
    elseif abs((singleGrowthOrg2 / pairedGrowthOrg2)) > 1 + sigD
        % second microbe grows slower
        iAFirstOrg = 'Parasitism_Taker';
        iASecondOrg = 'Parasitism_Giver';
        iATotal = 'Parasitism';
    else
        % if no case fits - needs inspection!
        iAFirstOrg = 'No_Result';
        iASecondOrg = 'No_Result';
        iATotal = 'No_Result';
    end
elseif abs((singleGrowthOrg1 / pairedGrowthOrg1)) > 1 + sigD
    % first microbe grows slower - all possible cases resulting froms
    % second microbe's growth
    if abs(1 - (pairedGrowthOrg2 / singleGrowthOrg2)) < sigD
        % second microbe unaffected
        iAFirstOrg = 'Amensalism_Affected';
        iASecondOrg = 'Amensalism_Unaffected';
        iATotal = 'Amensalism';
    elseif abs((pairedGrowthOrg2 / singleGrowthOrg2)) > 1 + sigD
        % second microbe grows better
        iAFirstOrg = 'Parasitism_Giver';
        iASecondOrg = 'Parasitism_Taker';
        iATotal = 'Parasitism';
    elseif abs((singleGrowthOrg2 / pairedGrowthOrg2)) > 1 + sigD
        % second microbe grows slower
        iAFirstOrg = 'Competition';
        iASecondOrg = 'Competition';
        iATotal = 'Competition';
    else
        % if no case fits - needs inspection!
        iAFirstOrg = 'No_Result';
        iASecondOrg = 'No_Result';
        iATotal = 'No_Result';
    end
else
    % if no case fits - needs inspection!
    iAFirstOrg = 'No_Result';
    iASecondOrg = 'No_Result';
    iATotal = 'No_Result';
end
end

author: opencobra

repo: cobratoolbox


text:
function vis_combined_scoremap(expidx, img_idx, ends)

p = exp_params(expidx);
load(p.testGT)

im_fn = annolist(img_idx).image.name;
[~,im_name,~] = fileparts(im_fn);
im = imread(im_fn);

scmap_name = fullfile(p.unary_scoremap_dir, [im_name '.mat']);
load(scmap_name, 'scoremaps');

colors = [1 0 1; 1 1 0; 0 1 1; 1 0 0; 0 1 0; 0 0 1; 1 1 1; 1 1 0; 1 0 1; 0 1 1; 1 0 0; 0 1 0; 0 0 1; 1 1 1];

pw_map_img = im;
un_map_img = im;

for end_ = ends

    if true
        scmap = [];

        for start = 1:14
            if start ~= end_
                current = visualise_pairwise_probabilities(expidx, img_idx, 1, start, end_, false);
                if isempty(current)
                    continue;
                end
                if isempty(scmap)
                    scmap = current;
                else
                    scmap = scmap .* current;
                end
            end
        end
        %save('~/pw_scmap.mat', 'scmap');
    else
        %load('~/pw_scmap.mat', 'scmap');
    end

    %figure(1);
    %imagesc(scmap);
    %colorbar;

    scmap_rescaled = scmap/max(scmap(:));
    pw_map_img = vis_scoremap_on_img(pw_map_img, scmap_rescaled, colors(end_,:));

    unary_map = scoremaps(:, :, end_);
    un_map_img = vis_scoremap_on_img(un_map_img, unary_map, colors(end_,:));
end

out_dir = '/tmp/';

figure(1);
fname = fullfile(out_dir, [num2str(img_idx) '_cidxs_all.png']);
showsaveimg(pw_map_img, fname);

figure(2);
fname = fullfile(out_dir, [num2str(img_idx) '_unary.png']);
showsaveimg(un_map_img, fname);

end

function showsaveimg(im, fname)
imshow(im);
axis off;
set(gca, 'LooseInset', get(gca, 'TightInset'));
print(gcf,'-dpng', fname);
end

function final_img = vis_scoremap_on_img(im, scoremap, color1)
scoremap_img = imresize(scoremap, [size(im, 1), size(im, 2)], 'bicubic');

c = reshape(color1, 1, 1, 3);
color_rep = repmat(c, size(im, 1), size(im, 2), 1);
color_rep = color_rep*255;

scoremap_img = repmat(scoremap_img, 1, 1, 3);
%figure(4);
%imagesc(pw_prob_img);

final_img = (scoremap_img) .* double(color_rep) + (1-scoremap_img) .* double(im);
final_img = uint8(final_img);

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [Flag,jAngle]=dialog_jointAnlgesInput()
% Dialog box used by the user to enter the values of the joint angles of the robot.

% Copyright: Mohammad SAFEEA, 2018-July-20

global daHandleof_dialog_jointAnlgesInput;
global dialog_jointAnlgesInput_button_hit;
global dialog_jointAnlgesInput_flag;

dialog_jointAnlgesInput_button_hit=0;
dialog_jointAnlgesInput_flag=false;
dialog_jointAnlgesInput_jangle={};
%% Interface container
dialog_width=540;
dialog_height=500;
dialog_x0=10;
dialog_y0=100;
daCaption='Joint angles (degree)';
jointAnglesDialogHandle=figure('Name',daCaption,'rend','painters','pos',[dialog_x0 dialog_y0 dialog_width dialog_height]);
daHandleof_dialog_jointAnlgesInput=jointAnglesDialogHandle;
%% Button dimension
enter_buttons_separator=40;
button_width=100;
button_height=50;
%% Add ok button
fcn=@callback_dialog_jointAnlgesInput_okBtn;
btn_ok_h= uicontrol(jointAnglesDialogHandle,'Style', 'pushbutton', 'String', 'OK',...
                'Units','pixels',...
                'Callback', fcn);
x0=round((dialog_width-enter_buttons_separator-button_width*2)/2);
xtemp=x0;
y0=enter_buttons_separator;
ytemp=y0;
btn_ok_h.HorizontalAlignment='center';
btn_ok_h.Position=[x0 y0 button_width button_height];
%% Add cancel button
fcn=@callback_dialog_jointAnlgesInput_cancelBtn;
btn_cancel_h= uicontrol(jointAnglesDialogHandle,'Style', 'pushbutton', 'String', 'Cancel',...
                'Units','pixels',...
                'Callback', fcn);
x0=xtemp+enter_buttons_separator+button_width;
y0=ytemp;
btn_cancel_h.HorizontalAlignment='center';
btn_cancel_h.Position=[x0 y0 button_width button_height];
%% Dimensions of lables and sliders,
labelsNum=7;
enter_controls_horizontal_margin=10;
x0=enter_controls_horizontal_margin;
y0=ytemp+button_height;

labels_width=100;
enter_labels_vertical_margin=30;
labels_height=(dialog_height-enter_labels_vertical_margin*(labelsNum+1)-y0)/labelsNum;

grandSliderWidth=280;
bigSliderWidth=250;
smallSliderWidth=220;
sliderHeight=labels_height;
slider_mid_coordinate=275;
%% Sliders handles
sliders_handles=[];
angle_value_texts_handles=[];
%% Add lables and sliders,
h=daHandleof_dialog_jointAnlgesInput;
labels_text_array={'J1','J2','J3','J4','J5','J6','J7'};
for i=1:7
    labels_text=labels_text_array{7-i+1};
    lbl_handle= uicontrol(h,'Style', 'text', 'String', labels_text,...
            'Units','pixels');
    x_labels_0=x0;
    y_labels_0=y0+(enter_labels_vertical_margin)*i+(labels_height)*(i-1);
    lbl_handle.HorizontalAlignment='center';
    lbl_handle.Position=[x_labels_0 y_labels_0 labels_width labels_height];
    lbl_handle.Enable='inactive';

    %% Add slider and its labels
    if i==1
        slider_width=grandSliderWidth;
        min=-175;
        max=175;
    elseif mod(i,2)==1
        slider_width=bigSliderWidth;
        min=-170;
        max=170;
    else
        slider_width=smallSliderWidth;
        min=-120;
        max=120;
    end
    x_Slider=getXCoordiantesFromMidWidth(slider_mid_coordinate,slider_width);
    y_Slider=y_labels_0;
    silder_handle = uicontrol('Parent',h,'Style','slider',...
         'Units','pixels', ...
        'Position',[x_Slider,y_Slider,slider_width,sliderHeight],...
              'min',min, 'max',max);
    sliders_handles{i}=silder_handle;
          %% Min labels
    minmaxLablesWidth=50;
    minmaxLabelHeight=sliderHeight;
    min_labels_x=x_Slider - minmaxLablesWidth;
    tempPos=[min_labels_x,y_Slider,minmaxLablesWidth,minmaxLabelHeight];
    bl1 = uicontrol(h,'Style','text',...
        'Units','pixels', ...
        'Position',tempPos,...
        'String',[num2str(min),char(176)]);

          %% Max labels
    max_labels_x=x_Slider +slider_width ;
    tempPos=[max_labels_x,y_Slider,minmaxLablesWidth,minmaxLabelHeight];
    bl1 = uicontrol(h,'Style','text',...
        'Units','pixels', ...
        'Position',tempPos,...
        'String',[num2str(max),char(176)]);

    %% Angle value text
    labelsWidth=80;
    angle_val_label_x=dialog_width-labelsWidth ;
    tempPos=[angle_val_label_x,...
        y_Slider,...
        labelsWidth,minmaxLabelHeight];
    angle_value_text_handle = uicontrol(h,'Style','text',...
        'Units','pixels', ...
        'Position',tempPos,...
        'String',num2str(0.0));
           angle_value_texts_handles{i} =angle_value_text_handle;
end

tempJointAngles={1,2,3,4,5,6,7};
dialog_jointAnlgesInput_jangle={0,0,0,0,0,0,0};
while(dialog_jointAnlgesInput_button_hit==0)
    figure(daHandleof_dialog_jointAnlgesInput)
    pause(0.1);
    try
    for i=1:7
        % get value of the slider
        temp_handle=sliders_handles{i};
        val=get(temp_handle,'Value');
        tempJointAngles{7-i+1}=val*pi/180;
        % through the value into the textbox
        temp_handle=angle_value_texts_handles{i} ;
%         numString=num2str(val)
        numString=sprintf('%0.2f',val);
        set(temp_handle,'String',numString);
        dialog_jointAnlgesInput_jangle=tempJointAngles;
    end
    catch
        break;
    end
end
jAngle=dialog_jointAnlgesInput_jangle;
Flag=dialog_jointAnlgesInput_flag;

clear daHandleof_dialog_jointAnlgesInput;
clear dialog_jointAnlgesInput_button_hit;
clear dialog_jointAnlgesInput_flag;
clear dialog_jointAnlgesInput_jangle;
end

function callback_dialog_jointAnlgesInput_okBtn(PushButton, EventData)
    global daHandleof_dialog_jointAnlgesInput;
    global dialog_jointAnlgesInput_button_hit;
    global dialog_jointAnlgesInput_flag;

    close(daHandleof_dialog_jointAnlgesInput);
    dialog_jointAnlgesInput_flag=true;
    dialog_jointAnlgesInput_button_hit=1;
end

function callback_dialog_jointAnlgesInput_cancelBtn(PushButton, EventData)
    global daHandleof_dialog_jointAnlgesInput;
    global dialog_jointAnlgesInput_button_hit;
    global dialog_jointAnlgesInput_flag;

    close(daHandleof_dialog_jointAnlgesInput);
    dialog_jointAnlgesInput_flag=false;
    dialog_jointAnlgesInput_button_hit=1;
end

function x=getXCoordiantesFromMidWidth(mid,width)
x=mid-width/2;
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function geometry_test085 ()

%*****************************************************************************80
%
%% TEST085 tests PLANE_EXP_PRO3.
%
%  Discussion:
%
%    Projection is ( -1, 1, 1 ).
%    Projection is ( 4, 5, -8 ).
%    Projection is ( 0.33, 0.33, 0.33).
%    Projection is ( 5.33, -1.66, -2.66 ).
%    Projection is ( -1, 1, 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  ntest = 5;

  p1 = [ 1.0, 0.0, 0.0 ];
  p2 = [ 0.0, 1.0, 0.0 ];
  p3 = [ 0.0, 0.0, 1.0 ];
  po = [ ...
    0.0,  2.0,  2.0;  ...
    4.0,  5.0, -8.0;  ...
    0.25, 0.25, 0.25; ...
    5.0, -2.0, -3.0;  ...
   -2.0,  0.0,  0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST085\n' );
  fprintf ( 1, '  PLANE_EXP_PRO3 projects an object point\n' );
  fprintf ( 1, '    orthographically into a plane.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '            PO                PP\n' );
  fprintf ( 1, '\n' );

  pp = plane_exp_pro3 ( p1, p2, p3, ntest, po );

  for i = 1 : ntest
    fprintf ( 1, '  %10f  %10f  %10f  %10f  %10f  %10f\n', ...
      po(1:dim_num,i), pp(1:dim_num,i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [map] = data2map(data)
%
%  Convert a data item (a cell or a vec) into a map by mapping values into their indices
%  Thang Luong @ 2013, 2015 <lmthang@stanford.edu>
%
  if isempty(data)
    map = containers.Map();
  else
    if iscell(data)
      map = containers.Map(data, 1:length(data));
    else
      map = containers.Map(num2cell(data), 1:length(data));
    end
  end

author: lmthang


text:
function rotMat = axisAngle2RotMat(axis,angle)
% Function to get the rotation matrix for a given axis and angle where x
% is a column vector, such that:
%    xNew = rotMat*x
%
% Alternate use where x is a row vector:
%    xNew = x*rotMat'
%
% Parameters
% ------------
%   axis : 1 x 3 float vector
%       Axis about which to rotate the point x
%   angle : float
%       Rotation angle (radian)
%
% Returns
% ------------
%   rotMat : 3 x 3 float vector
%       Rotation matrix from the input axis and angle
%

rotMat = zeros(3);
rotMat(1,1) = axis(1)*axis(1)*(1-cos(angle)) + cos(angle);
rotMat(1,2) = axis(2)*axis(1)*(1-cos(angle)) - axis(3)*sin(angle);
rotMat(1,3) = axis(3)*axis(1)*(1-cos(angle)) + axis(2)*sin(angle);
rotMat(2,1) = axis(1)*axis(2)*(1-cos(angle)) + axis(3)*sin(angle);
rotMat(2,2) = axis(2)*axis(2)*(1-cos(angle)) + cos(angle);
rotMat(2,3) = axis(3)*axis(2)*(1-cos(angle)) - axis(1)*sin(angle);
rotMat(3,1) = axis(1)*axis(3)*(1-cos(angle)) - axis(2)*sin(angle);
rotMat(3,2) = axis(2)*axis(3)*(1-cos(angle)) + axis(1)*sin(angle);
rotMat(3,3) = axis(3)*axis(3)*(1-cos(angle)) + cos(angle);

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
classdef TestMorphologyEx
    %TestMorphologyEx

    properties (Constant)
        img = uint8([...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 1 1 1 0 0 0 0;...
            0 0 0 1 1 1 0 0 0 0;...
            0 0 0 1 1 1 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0;...
        ]);
    end

    methods (Static)
        function test_1
            result = cv.morphologyEx(TestMorphologyEx.img, 'Tophat');
            validateattributes(result, {class(TestMorphologyEx.img)}, ...
                {'size',size(TestMorphologyEx.img)});
        end

        function test_2
            result = cv.morphologyEx(TestMorphologyEx.img, 'Tophat', ...
                'Element',eye(3), 'Anchor',[1,1], 'Iterations',1, ...
                'BorderType','Constant', 'BorderValue',255);
            validateattributes(result, {class(TestMorphologyEx.img)}, ...
                {'size',size(TestMorphologyEx.img)});
        end

        function test_error_argnum
            try
                cv.morphologyEx();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [selected] = ft_select_point3d(bnd, varargin)

% FT_SELECT_POINT3D helper function for selecting one or multiple points on a 3D mesh
% using the mouse. It returns a list of the [x y z] coordinates of the selected
% points.
%
% Use as
%   [selected] = ft_select_point3d(bnd, ...)
%
% Optional input arguments should come in key-value pairs and can include
%   'multiple'    = true/false, make multiple selections, pressing "q" on the keyboard finalizes the selection (default = false)
%   'nearest'     = true/false (default = true)
%   'marker'      = character or empty, for example '.', 'o' or 'x' (default = [])
%   'markersize'  = scalar, the size of the marker (default = 10)
%   'markercolor' = character, for example 'r', 'b' or 'g' (default = 'k')
%
% Example
%   [pos, tri] = mesh_sphere(162);
%   bnd.pos = pos;
%   bnd.tri = tri;
%   ft_plot_mesh(bnd)
%   camlight
%   ... do something here
%
% See also FT_SELECT_BOX, FT_SELECT_CHANNEL, FT_SELECT_POINT, FT_SELECT_RANGE, FT_SELECT_VOXEL

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

% get optional input arguments
nearest     = ft_getopt(varargin, 'nearest',  true);
multiple    = ft_getopt(varargin, 'multiple', false);
marker      = ft_getopt(varargin, 'marker', []);
markersize  = ft_getopt(varargin, 'markersize', 10);
markercolor  = ft_getopt(varargin, 'markercolor', 'k');

% ensure that it is boolean
nearest  = istrue(nearest);
multiple = istrue(multiple);

% get the object handles
h = get(gca, 'children');

% select the correct objects
iscorrect = false(size(h));
for i=1:length(h)
  try
    pos = get(h(i),'vertices');
    tri = get(h(i),'faces');
    if ~isempty(bnd) && isequal(bnd.pos, pos) && isequal(bnd.tri, tri)
      % it is the same object that the user has plotted before
      iscorrect(i) = true;
    elseif isempty(bnd)
      % assume that it is the same object that the user has plotted before
      iscorrect(i) = true;
    end
  end
end
h = h(iscorrect);

if isempty(h) && ~isempty(bnd)
  figure
  ft_plot_mesh(bnd);
  camlight
  selected = ft_select_point3d(bnd, varargin{:});
  return
end

if length(h)>1
  ft_warning('using the first patch object in the figure');
  h = h(1);
end

selected = zeros(0,3);

% everything is added to the current figure
holdflag = ishold;
if ~holdflag
  hold on
end

done = false;
az = 0;
el = 0;
view(az,el);

while ~done
  k = waitforbuttonpress;
  [p, v, vi, facev, facei] = select3d(h);
  if k == 1 %checks if waitforbuttonpress was a key
    key = get(gcf,'CurrentCharacter'); % which key was pressed (if any)?
    if strcmp(key, 'q')
      % finished selecting points
      done = true;
    elseif strcmp(key, 'r')
      % remove last point
      if ~isempty(selected)
        if ~isempty(marker)
          delete(findobj('marker', '*'));
          hs = plot3(selected(1:end-1,1), selected(1:end-1,2), selected(1:end-1,3), [markercolor marker]);
          set(hs, 'MarkerSize', markersize);
        end
        fprintf('removed latest point at [%f %f %f]\n', selected(end,1), selected(end,2), selected(end,3));
        selected = selected(1:end-1,:);
      end
    elseif strcmp(key,'+')
      zoom(1.1)
    elseif strcmp(key,'-')
      zoom(0.9)
    elseif strcmp(key,'w')
      az = az+6;
      view(az,el)
    elseif strcmp(key,'a')
      el = el+6;
      view(az,el)
    elseif strcmp(key,'s')
      az = az-6;
      view(az,el)
    elseif strcmp(key,'d')
      el = el-6;
      view(az,el)
    end

  elseif ~isempty(p)
    % a new point was selected
    if nearest
      selected(end+1,:) = v;
    else
      selected(end+1,:) = p;
    end % if nearest

    if multiple
      fprintf('selected points at\n');
      disp(selected);
    else
      fprintf('selected point at [%f %f %f]\n', selected(end,1), selected(end,2), selected(end,3));
    end

    if ~isempty(marker)
      hs = plot3(selected(end,1), selected(end,2), selected(end,3), [markercolor marker]);
      set(hs, 'MarkerSize', markersize);
    end
  end

  if ~multiple
    done = true;
  end
end

if ~holdflag
  hold off
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function infrasoundEvent = xcorr3C(wevent, infrasoundEvent, make_figures, figureOutDirectory, pretrigger)
%XCORR3C Cross-correlation an event recorded on 3 infrasound components
% infrasoundEvent = xcorr3C(infrasoundEvent)
%   Input:
%       wevent - a cell array where each component is a vector of 3
%           waveform objects, 1 per infrasound channel
%       infrasoundEvent is a structure containing two elements:
%           FirstArrivalTime
%           LastArrivalTime
%       make_figures - if true, a figure is generated for each xcorr pair
%
%   Output:
%       infrasoundEvent with some additional elements added
%           maxCorr - a 3x3 array of the maximum cross correlation values
%           secsDiff - a 3x3 array of the time lags corresponding to
%                      maxCorr
%           meanSecsDiff - the mean of secsDiff for non-diagonal components
%
%       each component is cross correlated with each component, hence 3x3

    %% correlate
    % loop through infrasound channels
    % take a 0.3-second snippet starting 0.1s before FirstArrivalTime, till 0.2s
    % after it
    % correlate this against the whole wevent for each infrasound
    % trace
    % this should result in a correlation matrix
    % from this record the time lag matrix for each infrasound channel against each other
    % infrasound channel
    disp('CORRELATION ...')
    disp('_______________')
    numEvents = numel(infrasoundEvent);
    for eventNumber=1:numEvents
        fprintf('- processing event %d of %d\n', eventNumber, numEvents);
        haystacks = wevent{eventNumber};
        infrasoundEvent(eventNumber).maxCorr = eye(3);
        infrasoundEvent(eventNumber).secsDiff = eye(3);
        precorrtime = 0.1; % NEEDLE seconds of data to add before first arrival
        postcorrtime = 0.2; % NEEDLE seconds of data to add after first arrival
        %postcorrtime = 0.4;
        for chanNumber=1:3
            needle = extract(haystacks(chanNumber), 'time', infrasoundEvent(eventNumber).FirstArrivalTime-precorrtime/86400, infrasoundEvent(eventNumber).FirstArrivalTime+postcorrtime/86400);
            needle_data = detrend(get(needle, 'data'));
%             % upsample by a factor of 8
%             nx = get(needle,'timevector');
%             nxx = nx(1):(nx(2) - nx(1))/8:nx(end);
%             nyy = spline(nx,needle_data,nxx);
            for haystackNum = 1:3
                fprintf('  - looking for needle %d in haystack %d\n', chanNumber, haystackNum);
                haystack = haystacks(haystackNum);
                haystack_data = get(haystack,'data');
%                 % upsample by a factor of 8
%                 hx = get(haystack,'timevector');
%                 hxx = hx(1):(hx(2) - hx(1))/8:hx(end);
%                 hyy = spline(hx,haystack_data,hxx);
                [acor,lag] = xcorr(needle_data, haystack_data);
%                 [acor,lag] = xcorr(nyy, hyy);
                cxx0 = sum(abs(needle_data).^2);
                cyy0 = sum(abs(haystack_data).^2);
%                  cxx0 = sum(abs(nyy).^2);
%                  cyy0 = sum(abs(hyy).^2);
                scale = sqrt(cxx0*cyy0);
                acor = acor./scale;
                [m,I] = max(abs(acor));
                infrasoundEvent(eventNumber).maxCorr(chanNumber,haystackNum) = m;
                infrasoundEvent(eventNumber).secsDiff(chanNumber,haystackNum) = lag(I)/get(haystack,'freq') + pretrigger - precorrtime;
                if make_figures
                    figure;
                    subplot(3,1,1),plot(haystack,'axeshandle',gca);
                    subplot(3,1,2),plot(needle,'axeshandle',gca);
                    subplot(3,1,3),plot(lag,acor);
                    outfile = sprintf('%s/xcorr_infrasoundEvent%03d_%d_%d.png',figureOutDirectory,eventNumber,chanNumber,haystackNum);
                    feval('print', '-dpng', outfile);
                    close
                end
            end
        end
        infrasoundEvent(eventNumber).meanCorr = mean(infrasoundEvent(eventNumber).maxCorr([2 3 4 6 7 8]));
        infrasoundEvent(eventNumber).stdCorr = std(infrasoundEvent(eventNumber).maxCorr([2 3 4 6 7 8]));
        infrasoundEvent(eventNumber).meanSecsDiff = mean(infrasoundEvent(eventNumber).secsDiff([2 3 4 6 7 8]));
        infrasoundEvent(eventNumber).stdSecsDiff = std(infrasoundEvent(eventNumber).secsDiff([2 3 4 6 7 8]));
    end
end

author: geoscience-community-codes

repo: GISMO


text:
%%********************************************************
%% Arrow:
%%
%% Fx = Arrow(pblk,f,x,options);
%%
%% if options == 0;
%%    Fx = Arr(F)*x
%% if options == 1;
%%    Fx = Arr(F)^{-1}*x
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%********************************************************

   function   Fx = Arrow(pblk,f,x,options);

   if nargin == 3; options = 0; end;

   s = 1 + [0, cumsum(pblk{2})];
   idx1 = s(1:length(pblk{2}));
   if options == 0
      inprod = mexqops(pblk{2},f,x,1);
      Fx  = mexqops(pblk{2},f(idx1),x,3) + mexqops(pblk{2},x(idx1),f,3);
      Fx(idx1) = inprod;
   else
      gamf2 = mexqops(pblk{2},f,f,2);
      gamprod = mexqops(pblk{2},f,x,2);
      alpha = gamprod./gamf2;
      Fx = mexqops(pblk{2},1./f(idx1),x,3) - mexqops(pblk{2},alpha./f(idx1),f,3);
      Fx(idx1) = alpha;
   end
%%
%%********************************************************

repo: REST


function [posterior,prior] = p2p(hh, shrinkage,prior0,olsreg,F,G,Fo,positions_nylags,position_constant)
%
%********************************************************
% Conjugate Prior: MN-IW for Direct Methods
%********************************************************

% settings
if isempty(position_constant) == 0
    nx = 1;
else
    nx = 0;
end

ny       = size(G,2);
nylags   = size(F,1);
iIminusF = inv(eye(nylags) - F);

% constructing the prior mean
Fhh        = F^hh;
Fohh       = (iIminusF * (eye(nylags)-Fhh)) * Fo;
if hh < 40
    priorSigmaScale  = zeros(ny);
    for j = 0 : hh
        priorSigmaScale = priorSigmaScale + G' * F^(hh-j)* G * prior0.Sigma.scale * G' * F^(hh-j)' * G;
    end
else
    % solves x-a*x*a'=b for b (and then x) symmetrical
    % function [x,info]=lyapunov_symm(a,b)
    [priorSigmaScale0,~] = lyapunov_symm(F,G*prior0.Sigma.scale*G' - F^(hh+1)* G * prior0.Sigma.scale * G' * F^(hh+1)');
    % max(max(abs(priorSigmaScale-G'*priorSigmaScale0*G)))
    priorSigmaScale      = G'*priorSigmaScale0*G;
end

prior.BetaMean   = [Fhh(1:ny,:)'; Fohh(1 : ny,1)'];
prior.BetaVar    = prior0.Phi.cov * 1/ shrinkage;

prior.df  = prior0.Sigma.df;            % usually number of regressors minus the 2
prior.XXi = inv( prior.BetaVar );       % V^{-1}
prior.S   = priorSigmaScale;            % Sigma0

% retrieve the OLS
B_   = olsreg.beta([positions_nylags position_constant], :);
XX_  = olsreg.X(:,[positions_nylags position_constant])' * olsreg.X(:,[positions_nylags position_constant]);
E_   = olsreg.error;
XXp_ = XX_ + prior.XXi;

% construct the posterior
posterior.df     = olsreg.N - nylags - nx + prior0.Sigma.df;
posterior.XXi    = inv( XXp_ );
posterior.PhiHat = posterior.XXi * (XX_ * B_ + prior.XXi * prior.BetaMean);
posterior.S     =  ...
    E_'* E_ + priorSigmaScale + prior.BetaMean' * prior.XXi * prior.BetaMean + ...
    B_' * XX_ * B_ - posterior.PhiHat' * XXp_ * posterior.PhiHat;
% FF    = prior.Sigma.scale + (posterior1.PhiHat - prior1.BetaMean)'* prior1.XXi * (posterior1.PhiHat - prior1.BetaMean) ...
%     + posterior1.E_'* posterior1.E_;

posterior.B_   = B_ ;
posterior.XX_  = XX_;       %olsreg.X(:,[positions_nylags position_constant])' * olsreg.X(:,[positions_nylags position_constant]);
posterior.E_   = E_;        %olsreg.error';
posterior.XXp_ = XXp_;      %XX_ + prior.XXi;
posterior.U_   = olsreg.Y - olsreg.X(:,[positions_nylags position_constant]) * posterior.PhiHat;% posterior errors;


text:
% std_setcomps2cell - convert .sets and .comps to cell array. The .sets and
%                     .comps format is useful for GUI but the cell array
%                     format is used for plotting and statistics.
%
% Usage:
%   [ struct setinds allinds ] = std_setcomps2cell(STUDY, clustind);
%   [ struct setinds allinds ] = std_setcomps2cell(STUDY, sets, comps);
%
% Author: Arnaud Delorme, CERCO/CNRS, UCSD, 2009-

% Copyright (C) Arnaud Delorme, arno@salk.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

function [ tmpstruct setinds allinds ] = std_setcomps2cell(STUDY, sets, comps, generateerror)

if nargin < 4
    generateerror = 0;
end;
if nargin < 3
    tmpstruct = STUDY.cluster(sets);
    sets  = tmpstruct.sets;
    comps = tmpstruct.comps; % old format
else
    tmpstruct     = [];
end;
comps   = repmat(comps, [size(sets,1) 1]);
oldsets = sets;
sets    = reshape(sets , 1, size(sets ,1)*size(sets ,2));
comps   = reshape(comps, 1, size(comps,1)*size(comps,2));

% get indices for all groups and conditions
% -----------------------------------------
setinfo       = STUDY.design(STUDY.currentdesign).cell;
allconditions = STUDY.design(STUDY.currentdesign).variable(1).value;
allgroups     = STUDY.design(STUDY.currentdesign).variable(2).value;
nc = max(length(allconditions),1);
ng = max(length(allgroups),    1);
allinds = cell( nc, ng );
setinds = cell( nc, ng );

for index = 1:length(setinfo)
    % get index of independent variables
    % ----------------------------------
    condind = std_indvarmatch( setinfo(index).value{1}, allconditions);
    grpind  = std_indvarmatch( setinfo(index).value{2}, allgroups    );
    if isempty(allconditions), condind = 1; end;
    if isempty(allgroups),     grpind  = 1; end;

    % get the position in sets where the dataset is
    % if several datasets check that they all have the same
    % ICA and component index
    % -----------------------
    datind  = setinfo(index).dataset;
    ind     = find(datind(1) == sets);
    if ~isempty(ind) && length(datind) > 1
        [ind1 ind2] = find(datind(1) == oldsets);
        columnica   = oldsets(:,ind2(1));
        if ~all(ismember(datind, columnica));
            disp('Warning: ***** change STUDY design as it combines datasets with different ICA decompositions');
        end;
    end;

    allinds{ condind, grpind } = [ allinds{ condind, grpind } comps(ind) ];
    setinds{ condind, grpind } = [ setinds{ condind, grpind } repmat(index, [1 length(ind)]) ];
end;
tmpstruct.allinds = allinds;
tmpstruct.setinds = setinds;

if generateerror && isempty(setinds{1})
    error( [ 'Some datasets not included in preclustering' 10 ...
             'because of partial STUDY design. You need to' 10 ...
             'use a STUDY design that includes all datasets.' ]);
end;

author: PatternRecognition

repo: OpenBMI


text:
%FtsZRing_Test
% FtsZ ring test class.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 1/4/2011
classdef FtsZRing_Test < edu.stanford.covert.cell.sim.CellStateTestCase
    methods
        function this = FtsZRing_Test(name)
            this = this@edu.stanford.covert.cell.sim.CellStateTestCase(name);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function particles = correction_step(particles, z)

% Weight the particles according to the current map of the particle
% and the landmark observations z.
% z: struct array containing the landmark observations.
% Each observation z(j) has an id z(j).id, a range z(j).range, and a bearing z(j).bearing
% The vector observedLandmarks indicates which landmarks have been observed
% at some point by the robot.

% Number of particles
numParticles = length(particles);

% Number of measurements in this time step
m = size(z, 2);

% TODO: Construct the sensor noise matrix Q_t (2 x 2)
Q_t = 0.01*eye(2);
% process each particle
for i = 1:numParticles
  robot = particles(i).pose;
  % process each measurement
  for j = 1:m
    % Get the id of the landmark corresponding to the j-th observation
    % particles(i).landmarks(l) is the EKF for this landmark
    l = z(j).id;

    % The (2x2) EKF of the landmark is given by
    % its mean particles(i).landmarks(l).mu
    % and by its covariance particles(i).landmarks(l).sigma

    % If the landmark is observed for the first time:
    if (particles(i).landmarks(l).observed == false)

      % TODO: Initialize its position based on the measurement and the current robot pose:

      particles(i).landmarks(l).mu = [robot(1)+z(j).range*cos(z(j).bearing + robot(3)); robot(2) + z(j).range*sin(z(j).bearing + robot(3))];

      % get the Jacobian with respect to the landmark position
      [h, H] = measurement_model(particles(i), z(j));

      % TODO: initialize the EKF for this landmark

      particles(i).landmarks(l).sigma = H\Q_t*inv(H)';

      % Indicate that this landmark has been observed
      particles(i).landmarks(l).observed = true;

    else

      % get the expected measurement
      [expectedZ, H] = measurement_model(particles(i), z(j));

      % TODO: compute the measurement covariance

	Q = H*particles(i).landmarks(l).sigma*H' + Q_t;

      % TODO: calculate the Kalman gain

	K = particles(i).landmarks(l).sigma*H'/Q;

      % TODO: compute the error between the z and expectedZ (remember to normalize the angle)

	z_diff = [z(j).range;z(j).bearing] - expectedZ;
	z_diff(2) = normalize_angle(z_diff(2));

      % TODO: update the mean and covariance of the EKF for this landmark

	particles(i).landmarks(l).mu = particles(i).landmarks(l).mu + K*z_diff;
	particles(i).landmarks(l).sigma = (eye(2) - K*H)*particles(i).landmarks(l).sigma;

      % TODO: compute the likelihood of this observation, multiply with the former weight
      %       to account for observing several features in one time step
	particles(i).weight = particles(i).weight*1/sqrt(det(2*pi*Q))*exp(-1/2*z_diff'/Q*z_diff);
    end

  end % measurement loop
end % particle loop

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


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

    enumeration
        SoIChunk('View By SoI Grouping')
        All('View All')
    end

    properties
        name(1,:) char = '';
    end

    methods
        function obj = ViewEventsTypeEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('ViewEventsTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('ViewEventsTypeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('ViewEventsTypeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function []=sb_identify_good_pixels()
%SB_IDENTIFY_GOOD_PIXELS
%   Do the closure tests to identify 'good pixels', i.e. ones that are
%   likely unwrapped properly.
%   Ekbal Hussain @ Uni Leeds
%   19/11/2013
%
%   ======================================================================
%   06/2014 AH: Modify for StaMPS and fix
%   04/2015 EH: Speed optimisations
%   09/2015 EH/DB: Account for Nan values
%   ======================================================================
logit;
fprintf('Identifying good pixels...\n')

load psver
psname=['ps',num2str(psver)];
phuwname=['phuw_sb',num2str(psver)];
loopname=['phuw_loops',num2str(psver)];
goodname=['phuw_good',num2str(psver)];

sb_make_closure_loops;

ps=load(psname);
load(loopname);
load(phuwname,'ph_uw')

good_pixels = false(ps.n_ps,ps.n_ifg);

for m=1:size(intfg_loops,1)

    %Sum positive interferograms in each loop
    positive = find(intfg_loops(m,:)==1);
    positive_ints = zeros(size(ph_uw,1),1);
    for p=1:size(positive,2)
       positive_ints = positive_ints + ph_uw(:,positive(1,p));
    end

    %Sum negative interferograms in each loop
    negative = find(intfg_loops(m,:)==-1);
    negative_ints = zeros(size(ph_uw,1),1);
    for p=1:size(negative,2)
       negative_ints = negative_ints + ph_uw(:,negative(1,p));
    end

    ints_used = [positive, negative];

    %Remove the average residual (closest 2*pi value) to ensure each int is
    %in the same reference
    average_resid = nanmean(positive_ints-negative_ints);

    average_resid = 2*pi*(round(average_resid/2/pi));

    phase_resid = positive_ints-negative_ints-average_resid;
    %hist(positive_ints-negative_ints,100)
    %a = phase_resid(find(phase_resid(:,1)<=1 & phase_resid(:,1)>=-1),:);
    %size(a,1)

    %Index each interferogram with the pixels that are good=1 and bad=0
    %A good pixel is good if: -1<closure residual<1
    good_pixels(phase_resid(:,1)<=1 & phase_resid(:,1)>=-1,ints_used)=true;
end

save(goodname,'good_pixels')
logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% Mathematics Q2876283
% https://math.stackexchange.com/questions/2876283
% Least Square with Optimization of Triangular Matrix
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     12/08/2018
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

numRowsA    = 6;
numColsA    = 20;

numRowsB    = 10;
numColsB    = numColsA;

stepSize        = 0.0025;
numIterations   = 500;

%% Generate Data

mA = randn([numRowsA, numColsA]);
mB = randn([numRowsB, numColsB]);

%% Projection onto Lower Triangular Matrices Set

mY = randn([numRowsA, numRowsA]);

cvx_begin('quiet')
    cvx_precision('best');
    variable mX(numRowsA, numRowsA) lower triangular
    minimize( norm(mX - mY, 'fro') )
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mX(:).'), ' ]']);
disp([' ']);

mX = tril(mY);

disp([' ']);
disp(['Projection Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(norm(mX - mY, 'fro'))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mX(:).'), ' ]']);
disp([' ']);

%% Solution by CVX

cvx_begin('quiet')
    cvx_precision('best');
    variable mX(numRowsA, numRowsB) lower triangular;
    minimize( 0.5 * pow_pos(norm( mX * mB - mA, 'fro' ), 2) );
cvx_end

disp([' ']);
disp(['CVX Solution Summary']);
disp(['The CVX Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(cvx_optval)]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mX(:).'), ' ]']);
disp([' ']);

%% Solution by Projected Gradient Descent

hObjFun = @(mX) (0.5 * (norm( (mX * mB) - mA, 'fro' ) ^ 2));
vObjVal = zeros([numIterations, 1]);

mBB = mB * mB.';
mAB = mA * mB.';

mX         = mAB * pinv(mBB); %<! Initialization by the Least Squares Solution
mX         = tril(mX);
vObjVal(1) = hObjFun(mX);

for ii = 2:numIterations

    mG = (mX * mBB) - mAB;
    mX = mX - (stepSize * mG);

    % Projection Step
    mX = tril(mX);

    vObjVal(ii) = hObjFun(mX);
end

disp([' ']);
disp(['Projected Gradient Descent Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(vObjVal(numIterations))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(mX(:).'), ' ]']);
disp([' ']);

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineSeries = plot(1:numIterations, [vObjVal, cvx_optval * ones([numIterations, 1])]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['Objective Function Value vs. Iteration']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Objective Function Value', ...
    'FontSize', fontSizeAxis);
set(hAxes, 'XLim', [1, numIterations]);
hLegend = ClickableLegend({['Projected Gradient Descent'], ['Optimal Value (CVX)']});
set(hAxes, 'LooseInset', [0.07, 0.07, 0.07, 0.07]);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
%function []=mcmc()

clear all; clc;

%% MCMC options
nloop=110; % number of iterations for mcmcm
burnin=10;  % number of initial iterations to discard
thin=10;    % save every nth iteration
M=4;        %number of group components desired

%% set starting values for Z, S_BAR, SIGMA_S B_BAR SIGMA_B MU
load real_data;
[Z_st S_BAR_st SIGMA_S_st B_BAR_st SIGMA_B_st N_k_st N_k1k2_st]=get_initial_values_new(N,M_i,S,B,Q,M);
Z=Z_st; % Mi x M matrices of group indicators
S_BAR=S_BAR_st; % cluster centroid locations
SIGMA_S=SIGMA_S_st; % cluster centroid variances
B_BAR=B_BAR_st; % group level connectivities
SIGMA_B=SIGMA_B_st; % variances of connectivities
N_k=N_k_st; % number of components that belong to each cluster
N_k1k2=N_k1k2_st; % pairwise counts of group level clusters
MU=ones(M,1)/M; %probabilities of membership in group-level clusters.

%% initialize arrays of posterior draws
keep=round((nloop-burnin)/thin);
Z_array=cell(keep,1);
S_BAR_array=cell(keep,1);
SIGMA_S_array=cell(keep,1);
B_BAR_array=cell(keep,1);
SIGMA_B_array=cell(keep,1);
MU_array=cell(keep,1);
iter_array=0;

%%Hyperparameters
c=10000;                    %'c' in the paper
eta=median(N_k_st);         % 'eta' in paper
SS=eta*mean(SIGMA_S,3);     % 'S' in the paper
p1=1;                       % 'a' in the paper
p2=1;                       % 'b' in the paper

%% run MCMC

for iter=1:nloop

    % Draw S_BAR (group centroid locations)
    for k=1:M
        mu_s_bar_k=zeros(3,1);
        Sigma_s_bar_k=inv(N_k(k)*inv(SIGMA_S(:,:,k))+inv(SS));
        for i=1:N
            if max(Z{i}(:,k))==1
                j_k=find(Z{i}(:,k)==1);
                for j=1:length(j_k)
                    mu_s_bar_k=mu_s_bar_k+S{i}(j_k(j),:)';
                end
            end
        end
        mu_s_bar_k=Sigma_s_bar_k*inv(SIGMA_S(:,:,k))*mu_s_bar_k;
        S_BAR(k,:)=mvnrnd(mu_s_bar_k,Sigma_s_bar_k);
    end

    % Draw SIGMA_S (group centoid covariance matrices)
    for k=1:M
        eta_k=eta+.5*N_k(k);
        SS_k=SS;
        for i=1:N
            if max(Z{i}(:,k))==1
                j_k=find(Z{i}(:,k)==1);
                for j=1:length(j_k)
                    SS_k=SS_k+.5*(S{i}(j_k(j),:)-S_BAR(k,:))'*(S{i}(j_k(j),:)-S_BAR(k,:));
                end
            end
        end
        SS_k=.5*(SS_k+SS_k');
        inv_SS_k=inv(SS_k);
        inv_SS_k=.5*(inv_SS_k+inv_SS_k');
        SIGMA_S(:,:,k)=inv(wishrnd(inv_SS_k,eta_k));
    end

    % Draw B_BAR (group mean connectivites)
    for k1=1:M
        for k2=1:M
            mu_b_bar_k1k2=zeros(Q,1);
            Sigma_sq_b_k1k2=1/(1/c+N_k1k2(k1,k2)/SIGMA_B(k1,k2))*eye(Q);
            for i=1:N
                n_k1=sum(Z{i}(:,k1));
                n_k2=sum(Z{i}(:,k2));
                if(and(n_k1>=1,n_k2>=1))
                    j_k1=find(Z{i}(:,k1)==1);
                    j_k2=find(Z{i}(:,k2)==1);
                    for j1=1:n_k1
                        for j2=1:n_k2
                            mu_b_bar_k1k2=mu_b_bar_k1k2+B{i}{j_k1(j1),j_k2(j2)}';
                        end
                    end
                end
            end
            mu_b_bar_k1k2=Sigma_sq_b_k1k2*mu_b_bar_k1k2/SIGMA_B(k1,k2);
            B_BAR{k1,k2}=mvnrnd(mu_b_bar_k1k2,Sigma_sq_b_k1k2);
        end
    end

    % Draw SIGMA_B (group connectivity variances)
    for k1=1:M
        for k2=1:M
            p1_k1k2=p1+.5*Q*N_k1k2(k1,k2);
            p2_k1k2=p2;
            for i=1:N
                n_k1=sum(Z{i}(:,k1));
                n_k2=sum(Z{i}(:,k2));
                if(and(n_k1>=1,n_k2>=1))
                    j_k1=find(Z{i}(:,k1)==1);
                    j_k2=find(Z{i}(:,k2)==1);
                    for j1=1:n_k1
                        for j2=1:n_k2
                            p2_k1k2=p2_k1k2+.5*sum((B{i}{j_k1(j1),j_k2(j2)}-B_BAR{k1,k2}).^2);
                        end
                    end
                end
            end
            SIGMA_B(k1,k2)=1/gamrnd(p1_k1k2,1/p2_k1k2);
        end
    end

    % Draw Z (indicators of group membership)
    ind=randperm(N);
    for i=ind
        ind_i=randperm(M_i(i));
        for j=ind_i
            Log_p_ij=zeros(M,1);
            for k=1:M
                Sigma_s_ijk=SIGMA_S(:,:,k);
                s_ijk=S{i}(j,:)-S_BAR(k,:);
                log_p_s_ijk=-.5*log(det(Sigma_s_ijk))-.5*s_ijk*inv(Sigma_s_ijk)*s_ijk';
                log_p_b_ijk=0;
                for j1=1:M_i(i)
                    if j==j1
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(k,k)))...
                        -.5*sum((B{i}{j,j}-B_BAR{k,k}).^2)/SIGMA_B(k,k);
                    end
                    if and(not(j==j1),find(Z{i}(j1,:)==1)==k)
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(k,k)))...
                        -.5*sum((B{i}{j1,j1}-B_BAR{find(Z{i}(j1,:)==1),find(Z{i}(j1,:)==1)}).^2)/...
                            SIGMA_B(find(Z{i}(j,:)==1),find(Z{i}(j,:)==1));
                    end
                    if not(find(Z{i}(j,:)==1)==find(Z{i}(j1,:)==1))
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(k,find(Z{i}(j1,:)==1))))...
                        -.5*sum((B{i}{j,j1}-B_BAR{k,find(Z{i}(j1,:)==1)}).^2)/...
                              SIGMA_B(k,find(Z{i}(j1,:)==1));
                    end
                end
                for j2=1:M_i(i)
                    if not(find(Z{i}(j2,:)==1)==find(Z{i}(j,:)==1))
                        log_p_b_ijk=log_p_b_ijk-.5*Q*log(det(SIGMA_B(find(Z{i}(j2,:)==1),k)))...
                        -.5*sum((B{i}{j2,j}-B_BAR{find(Z{i}(j2,:)==1),k}).^2)/...
                                SIGMA_B(find(Z{i}(j2,:)==1),k);
                    end
                end
                Log_p_ij(k)=MU(k)+log_p_s_ijk+log_p_b_ijk;
            end
            Log_p_ij=Log_p_ij-max(Log_p_ij);
            P_ij=exp(Log_p_ij)/sum(exp(Log_p_ij));
            v=rand;
            Z{i}(j,:)=0*Z{i}(j,:);
            if v<=P_ij(1)
                   Z{i}(j,1)=1;
                   log_p_ij=Log_p_ij(1);
            end
            for k=1:(M-1)
               if(and(v>sum(P_ij(1:k)),v<=sum(P_ij(1:(k+1)))))
                   Z{i}(j,k+1)=1;
                   log_p_ij=Log_p_ij(k);
               end
            end
        end
    end

    % Draw MU (probabilities of clusters)
    N_k1k2=zeros(M);
    for k1=1:M
        for k2=1:M
            for i=1:N
                n_k1=sum(Z{i}(:,k1));
                n_k2=sum(Z{i}(:,k2));
                if(and(n_k1>=1,n_k2>=1))
                    if not(k1==k2)
                        N_k1k2(k1,k2)=N_k1k2(k1,k2)+n_k1*n_k2;
                    end
                    if k1==k2
                        N_k1k2(k1,k2)=N_k1k2(k1,k2)+n_k1;
                    end
                end
            end
        end
    end
    N_k=diag(N_k1k2);
    MU=drchrnd(1/M*ones(1,M)+N_k',1);

    % Posterior Draw arrays
    if and(iter>=burnin,thin*round(iter/thin)==iter)
        iter_array=iter_array+1;
        Z_array{iter_array}=Z;
        S_BAR_array{iter_array}=S_BAR;
        SIGMA_S_array{iter_array}=SIGMA_S;
        B_BAR_array{iter_array}=B_BAR;
        SIGMA_B_array{iter_array}=SIGMA_B;
        MU_array{iter_array}=MU;

        %compute and print out intermediate estimates of source locations
        S_BAR_mean=S_BAR_array{1};
        for j=2:iter_array
            S_BAR_mean=S_BAR_mean+S_BAR_array{j};
        end
        S_BAR_mean=S_BAR_mean/iter_array;
        iter
        [S_BAR_st zeros(M,1) S_BAR_mean]

        %[sqrt(sum(min(pdist2(S_BAR_true,S_BAR_st)').^2))/M sqrt(sum(min(pdist2(S_BAR_true,S_BAR_mean)').^2))/M]

        %compute and print intermediate indicator probs for random subject
        i=discretesample(repmat(1,N,1)/N,1);
        Z_i_mean=Z_array{1}{i};
        for j=2:iter_array
           Z_i_mean=Z_i_mean+Z_array{j}{i};
        end
        Z_st{i}
        Z_i_mean=Z_i_mean/iter_array
    end

end

save mcmc_outputs nloop M Z_array S_BAR_array SIGMA_S_array B_BAR_array SIGMA_B_array MU_array iter_array


text:
clear;
clc;
close all;

A = [1 1 1 1; 2 5  7 8]';
B = [1 2 3 3]';

X = inv((A'*A))*A'*B

X = pinv(A)*B

[U,S,V] = svd(A);
PINV_A = V*pinv(S)*U';
X = PINV_A*B

X = U*U'*B

author: yandld


text:
%% Transient diffusion equation
%% PDE and boundary conditions
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc

%% Define the domain and create a mesh structure
L = 50;  % domain length
Nx = 20; % number of cells
m = createMesh3D(Nx,Nx,Nx, L,L,L);
%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure
BC.left.a(:) = 0; BC.left.b(:)=1; BC.left.c(:)=0; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
BC.top.a(:) = 0; BC.top.b(:)=1; BC.top.c(:)=0; % top boundary
BC.bottom.a(:) = 0; BC.bottom.b(:)=1; BC.bottom.c(:)=0; % bottom boundary
%% define the transfer coeffs
D_val = 1;
D = createCellVariable(m, D_val);
alfa = createCellVariable(m, 1);
%% define initial values
c_init = 1;
c_old = createCellVariable(m, c_init,BC); % initial values
c = c_old; % assign the old value of the cells to the current values
%% loop
dt = 1; % time step
final_t = 50;
for t=dt:dt:final_t
    [M_trans, RHS_trans] = transientTerm(c_old, dt, alfa);
    Dave = harmonicMean(D);
    Mdiff = diffusionTerm(Dave);
    [Mbc, RHSbc] = boundaryCondition(BC);
    M = M_trans-Mdiff+Mbc;
    RHS = RHS_trans+RHSbc;
    c = solvePDE(m,M, RHS);
    c_old = c;
end
%% visualization
visualizeCells(c)

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [params, result] = ortho_default(params, W, w)
%  DSS No orthogonalization function
%  Does mere normalisation
%   W = orthof(W)     for symmetric dss
%   w = orthof(W, w)  for deflation dss
%     W Matrix with projection vectors as rows. For deflation
%       algorithm only previously calculated projections are given.
%     w Currently iterated projection. Only for deflation algorithm.

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Licensed under the Creative Commons Attribution-NonCommercial-ShareAlike
% License. http://creativecommons.org/licenses/by-nc-sa/2.0/.
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<2
    params.name = 'No orthogonalization';
    params.description = 'Description of this function.';
    return;
end

disp('no orthogonalisation!');
if nargin>2
  % per component orthogonalization
  result = w / norm(w);
else
  % symmetric orthogonalization
  result = diag(1./ sqrt(diag(W * W'))) * W;
end


text:
function refactor_fcn_name(fcnname, newfcnname, topdir, domove)
% refactor the name of a function, changing all references to the function
% name in the path and moving the function file to the new named file
%
% Syntax
%
% refactor_fcn_name(fcnname, newfcnname)
% refactor_fcn_name(fcnname, newfcnname, topdir)
% refactor_fcn_name(fcnname, newfcnname, topdir, domove)
%
% Description
%
% refactor_fcn_name(fcnname, newfcnname) finds all uses of the function
% name in fcnname and replaces it with the string in newfcnname in the
% entire matlab path. The function must be function in an m-file on the
% Matlab path. By default the mfile containing the fuction is also then
% renamed to the new name.
%
% refactor_fcn_name(fcnname, newfcnname, topdir) performs the same action
% but searching the folder provided in 'topdir' and all it's subdirectories
% rather than the entire Matlab path.
%
% refactor_fcn_name(fcnname, newfcnname, topdir, domove) performs the same
% action but the 'domove' flag determines whether the function file is
% actually moved or not. If this evaluates to true the file is moved to the
% mfile with the new name, if false nothing is moved.
%
%
% See also: regexprepfile
%

% Created by Richard Crozier 2013

    if ~ischar(fcnname)
        error('fcnname must be a string');
    end

    if ~ischar(newfcnname)
        error('newfcnname must be a string');
    end

    % refactor_fcn_name
    loc = which(fcnname);

    if isempty(loc)
        error('fcnname is not a function or class on the path')
    end

    [pathstr, ~, ext] = fileparts(which(fcnname));

    if isempty(strcmpi(ext, '.m'))
       error('function is not an m-file on the path.')
    end

    if nargin < 3 || isempty(topdir)
        thepath = path2cell(path);
    else
        if ~ischar(topdir)
            error('topdir must be a string');
        end
        if exist(topdir, 'file') ~= 7
            error('supplied directory name does not exist.')
        end
        thepath = path2cell(genpath(topdir));
    end

    if nargin < 4
        domove = true;
    end

    for indi = 1:numel(thepath)

        mfiles = dir([ thepath{indi}, '\*.m' ]);

        for indii = 1:numel(mfiles)
            % replace string in file using regexprepfile, with a regex
            % looking for word boundaries to avoid functions with similar
            % substrings
            regexprepfile(fullfile(thepath{indi}, mfiles(indii).name), ['\<', fcnname, '\>'], newfcnname);
        end

    end

    if domove
        % finally move the function file to the new m-file with the correct
        % name,
        if ispc
            % On windows we can't use the built-in movefile function as the
            % file system is not case sensitive, and will refuse to move
            % the file to what it beleives is the same location, so we take
            % the roundabout route

            % move to a temporary location
            tmpdest = [tempname, '_', newfcnname, ext];
            movefile(fullfile(pathstr, [fcnname, ext]), tmpdest);

            % move it again to where we really want it
            movefile(tmpdest, fullfile(pathstr, [newfcnname, ext]));

        else
            % on other systems, just move it!
            movefile(fullfile(pathstr, [fcnname, ext]), fullfile(pathstr, [newfcnname, ext]));
        end

    end

end

function pathlist = path2cell(pathstr)
%PATH2CELL Convert search path to cell array.
%
%   PATH2CELL returns a cell array where each element is a directory
%   in the search path.
%
%   PATH2CELL(MYPATH) converts MYPATH rather than the search path.
%
%   Empty directories are removed, so under UNIX, PATH2CELL('foo::bar')
%   returns {'foo', 'bar'} rather than {'foo', '', 'bar'}.

%   Author:      Peter J. Acklam
%   Time-stamp:  2001-10-18 21:23:19 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
%    error(nargchk(0, 1, nargin));

   % use MATLAB's search path if no input path is given
   if ~nargin
      pathstr = path;
   end

   k = find(pathstr == pathsep);            % find path separators
   k = [0 k length(pathstr)+1];             % find directory boundaries
   ndirs = length(k)-1;                     % number of directories
   pathlist = cell(0);                      % initialize output argument
   for i = 1 : ndirs
      dir = pathstr(k(i)+1 : k(i+1)-1);     % get i'th directory
      if ~isempty(dir)                      % if non-empty...
         pathlist{end+1,1} = dir;           % ...append to list
      end
   end
end

author: Sable

repo: mcbench-benchmarks


text:
function dat_comb = searchlight_applymask_collate(dat2, file_list)
% Estimate local pattern weight on train data using SVR and searchlight and
% apply weights to the test dataset
%
% :Usage:
% ::
%
%     dat_comb = searchlight_applymask(dat2, file_list)
%
% :Inputs:
%
%   **dat2:**
%        fmri_data test object
%
%   **file_list:**
%        Cellarray of list of file distributed in parallel.  Make sure
%        it is sorted correctly (e.g., sort_nat())
%
% :Output:
%
%   **dat_comb:**
%        An fmri_data object with searchlight weights applied to
%        test dataset.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2015  Luke Chang & Wani Woo
%
%     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/>.
% ..
%
% :Example:
% ::
%
%    dat_comb = searchlight_applymask_collate(dat2, file_list);

if ~isa(dat2,'fmri_data') % Check inputs
    dat1 = fmri_data(dat2); dat2 = remove_empty(dat2);
end

if ~iscell(file_list)
    error('Make sure searchlight file list is a cellarray and that files are on the path')
end

% Collate searchlight .mat files
vox_num = size(dat2.dat,1);
dist_n = length(file_list);
dat_comb = dat2;
for i = 1:dist_n
    load(file_list{i});
    dist_indx = select_voxels(i, dist_n, vox_num);
    dat_comb.dat(dist_indx,:) = dat.dat;
end
end

% ========== SUBFUNCTION ===========

function dist_indx = select_voxels(dist_id, dist_n, vox_num)
% preparation of distribution indx

dist_indx = false(vox_num,1);

unit_num = ceil(vox_num/dist_n);
unit = true(unit_num,1);

start_point = (unit_num.*(dist_id-1)+1);
end_point = min(start_point+unit_num-1, vox_num);
dist_indx(start_point:end_point) = unit(1:end_point-start_point+1);

dist_indx = logical(dist_indx);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = subsref(v,s)
%overloads subsref

switch s(1).type
  case '()'

    v.x = subsref(v.x,s(1));
    v.y = subsref(v.y,s(1));
    v.z = subsref(v.z,s(1));

    v.opt = struct;

    if numel(s)>1
      [varargout{1:nargout}] = builtin('subsref',v,s(2:end));
    else
      varargout{1} = v;
    end
  otherwise

    [varargout{1:nargout}] = builtin('subsref',v,s);

end
end

author: mtex-toolbox

repo: mtex

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


text:
function S = SOAFinit(w0,mu,leak)

% SOAFinit       Initialize Parameter Structure for the SOAF-LMS Algorithm
%
% Arguments:
% w0             Coefficients of FIR filter at start (@n=1)
% mu             Step size for LMS algorithm
% leak           Leaky factor:
%                  leak = 0 for conventional LMS
%                  0 < leak <1/step for leaky LMS algorithm
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

if nargin < 3                 % Set default to conventional LMS algorithm
	leak = 0;
end

% Assign structure fields

S.coeffs        = w0(:);      % Coefficients of FIR filter
S.step          = mu;         % Step size
S.leakage       = leak;       % Leakage factor
S.iter          = 0;          % Iteration count
S.AdaptStart    = length(w0); % Running effects
M               = length(w0); % Length of filter
S.beta          = 1/M;        % Forgetting factor

for j=1:M,                    % DCT-transform matrix
   S.T(1,j)=1/sqrt(M);
   for i=2:M,
     S.T(i,j)=sqrt(2/M)*cos(pi*(i-1)*(2*j-1)/2/M);
   end
end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
% gammaScanScript
global planC
indexS = planC{end};

[newXgrid, newYgrid, newZgrid] = getScanXYZVals(planC{indexS.scan}(1));
deltaX = abs(newXgrid(2) - newXgrid(1));
deltaY = abs(newYgrid(2) - newYgrid(1));
deltaZ = abs(newZgrid(2) - newZgrid(1));
doseAgreement = 15;
distAgreement = 0.3;
thresholdAbsolute = 4000;
doseArray1 = single(planC{indexS.scan}(1).scanArray);
doseArray2 = single(planC{indexS.scan}(2).scanArray);
gammaM = gammaDose3d(doseArray1, doseArray2, [deltaX deltaY deltaZ], doseAgreement, distAgreement, [], thresholdAbsolute);

newDoseNum = length(planC{indexS.dose}) + 1;
%Remove old caching info.
planC{indexS.dose}(newDoseNum).cachedMask = [];
planC{indexS.dose}(newDoseNum).cachedColor = [];
planC{indexS.dose}(newDoseNum).cachedTime = [];
%Set coordinates.
planC{indexS.dose}(newDoseNum).sizeOfDimension1 = length(newXgrid);
planC{indexS.dose}(newDoseNum).sizeOfDimension2 = length(newYgrid);
planC{indexS.dose}(newDoseNum).sizeOfDimension3 = length(newZgrid);
planC{indexS.dose}(newDoseNum).horizontalGridInterval = newXgrid(2)-newXgrid(1);
planC{indexS.dose}(newDoseNum).verticalGridInterval = newYgrid(2)-newYgrid(1);
planC{indexS.dose}(newDoseNum).depthGridInterval = newZgrid(2)-newZgrid(1);
planC{indexS.dose}(newDoseNum).coord1OFFirstPoint = newXgrid(1);
planC{indexS.dose}(newDoseNum).coord2OFFirstPoint = newYgrid(1);
planC{indexS.dose}(newDoseNum).coord3OfFirstPoint = newZgrid(1);
planC{indexS.dose}(newDoseNum).zValues = newZgrid;
planC{indexS.dose}(newDoseNum).doseUnits = 'Gy';
planC{indexS.dose}(newDoseNum).doseArray = gammaM;
planC{indexS.dose}(newDoseNum).doseUID = createUID('dose');
planC{indexS.dose}(newDoseNum).fractionGroupID = 'Gamma Scan 4mm, 50HU';
%Switch to new dose
sliceCallBack('selectDose', num2str(newDoseNum));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function CNR =Compute_VM(IM,sigma1,sigma2)

%  This new visibility metric is to calculate the Contrast-to-Noise Ratio
%  of noise image estimated by Gaussian kernel with sigma1 and sigma2 respectively.
%  Written by zhengguo dai in Beijing JX Digital Wave Co.ltd
%  Email: zhgdai@126.com
%  Note: This new metric should be tested by more pictures.

epsilon=1e-1;
half_size1 = ceil( -norminv( epsilon/2, 0, sigma1 ) );
size1 = 2 * half_size1 + 1 ;
half_size2 = ceil( -norminv( epsilon/2, 0, sigma2 ) );
size2 = 2 * half_size2 + 1 ;
gaussian1 = fspecial( 'gaussian', size1, sigma1 ) ;
gaussian2 = fspecial( 'gaussian', size2, sigma2 ) ;
IM = abs(IM - imfilter( IM, gaussian1, 'replicate' )) ;
noise = IM - imfilter( IM, gaussian1, 'replicate' ) ;
lsd = sqrt(imfilter( noise.^2, gaussian2, 'replicate' ) ) ;
lsd =  lsd./max( lsd( : ) ) .*255 ; % normalize
% figure; imshow(lsd,[])
CNR = sum(lsd(:))/(size(lsd,1)*size(lsd,2));


text:
function [ qy, qty, b, rsd, xb, info ] = zqrsl ( x, ldx, n, k, qraux, y, job )

%*****************************************************************************80
%
%% ZQRSL solves, transforms or projects systems factored by ZQRDC.
%
%  Discussion:
%
%    The routine applies the output of ZQRDC to compute coordinate
%    transformations, projections, and least squares solutions.
%
%    For K <= min ( N, P ), let XK be the matrix
%
%      XK = ( X(IPVT(1)), X(IPVT(2)), ... ,X(IPVT(k)) )
%
%    formed from columnns IPVT(1), ... ,IPVT(K) of the original
%    N by P matrix X that was input to ZQRDC (if no pivoting was
%    done, XK consists of the first K columns of X in their
%    original order).  ZQRDC produces a factored unitary matrix Q
%    and an upper triangular matrix R such that
%
%      XK = Q * ( R )
%               ( 0 )
%
%    This information is contained in coded form in the arrays
%    X and QRAUX.
%
%    The parameters QY, QTY, B, RSD, and XB are not referenced
%    if their computation is not requested and in this case
%    can be replaced by dummy variables in the calling program.
%
%    To save storage, the user may in some cases use the same
%    array for different parameters in the calling sequence.  A
%    frequently occuring example is when one wishes to compute
%    any of B, RSD, or XB and does not need Y or QTY.  In this
%    case one may identify Y, QTY, and one of B, RSD, or XB, while
%    providing separate arrays for anything else that is to be
%    computed.  Thus the calling sequence
%
%      call zqrsl ( x, ldx, n, k, qraux, y, dum, y, b, y, dum, 110, info )
%
%    will result in the computation of B and RSD, with RSD
%    overwriting Y.  More generally, each item in the following
%    list contains groups of permissible identifications for
%    a single callinng sequence.
%
%    1. ( Y, QTY, B )   ( RSD )      ( XB )  ( QY )
%    2. ( Y, QTY, RSD ) ( B )        ( XB )  ( QY )
%    3. ( Y, QTY, XB )  ( B )        ( RSD ) ( QY )
%    4. ( Y, QY )       ( QTY, B )   ( RSD ) ( XB )
%    5. ( Y, QY )       ( QTY, RSD ) ( B )   ( XB )
%    6. ( Y, QY )       ( QTY, XB )  ( B )   ( RSD )
%
%    In any group the value returned in the array allocated to
%    the group corresponds to the last member of the group.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex X(LDX,P), the output of ZQRDC.
%
%    Input, integer LDX, the leading dimension of X.
%
%    Input, integer N, the number of rows of the matrix XK, which
%    must have the same value as N in ZQRDC.
%
%    Input, integer K, the number of columns of the matrix XK.  K must not
%    be greater than min ( N, P), where P is the same as in the calling
%    sequence to ZQRDC.
%
%    Input, complex QRAUX(P), the auxiliary output from ZQRDC.
%
%    Input, complex Y(N), a vector that is to be manipulated by ZQRSL.
%
%    Input, integer JOB, specifies what is to be computed.  JOB has
%    the decimal expansion ABCDE, meaning:
%    if A /= 0, compute QY.
%    if B, D, D, or E /= 0, compute QTY.
%    if C /= 0, compute B.
%    if D /= 0, compute RSD.
%    if E /= 0, compute XB.
%    A request to compute B, RSD, or XB automatically triggers the
%    computation of QTY, for which an array must be provided in the
%    calling sequence.
%
%    Output, complex QY(N), contains Q*Y, if it has been requested.
%
%    Output, complex QTY(N), contains hermitian(Q)*Y, if it has
%    been requested.  Here hermitian(Q) is the conjugate transpose
%    of the matrix Q.
%
%    Output, complex B(K), the solution of the least squares problem
%      minimize norm2 ( Y - XK * B ),
%    if it has been requested.  If pivoting was requested in ZQRDC,
%    the J-th component of B will be associated with column IPVT(J)
%    of the original matrix X that was input into ZQRDC.
%
%    Output, complex RSD(N), the least squares residual Y - XK*B,
%    if it has been requested.  RSD is also the orthogonal projection
%    of Y onto the orthogonal complement of the column space of XK.
%
%    Output, complex XB(N), the least squares approximation XK*N,
%    if its computation has been requested.  XB is also the orthogonal
%    projection of Y onto the column space of X.
%
%    Output, integer INFO, is zero unless the computation of B has
%    been requested and R is exactly singular.  In this case, INFO is the
%    index of the first zero diagonal element of R and B is left unaltered.
%
  qy = [];
  qty = [];
  b = [];
  rsd = [];
  xb = [];
  info = 0;
%
%  Determine what is to be computed.
%
  cqy =  ( floor (       job / 10000        ) ~= 0 );
  cqty = (         mod ( job,  10000 )        ~= 0 );
  cb =   ( floor ( mod ( job,   1000 ) /100 ) ~= 0 );
  cr =   ( floor ( mod ( job,    100 ) / 10 ) ~= 0 );
  cxb =  (         mod ( job,     10 )        ~= 0 );

  ju = min ( k, n - 1 );
%
%  Special action when N=1.
%
  if ( ju == 0 )

    if ( cqy )
      qy(1) = y(1);
    end

    if ( cqty )
      qty(1) = y(1);
    end

    if ( cxb )
      xb(1) = y(1);
    end

    if ( cb )
      if ( zabs1 ( x(1,1) ) == 0.0 )
        info = 1;
      else
        b(1) = y(1) / x(1,1);
      end
    end

    if ( cr )
      rsd(1) = 0.0;
    end

    return

  end
%
%  Set up to compute QY or QTY.
%
  if ( cqy )
    qy(1:n) = y(1:n);
  end

  if ( cqty )
    qty(1:n) = y(1:n);
  end
%
%  Compute QY.
%
  if ( cqy )

    for jj = 1 : ju

      j = ju - jj + 1;

      if ( zabs1 ( qraux(j) ) ~= 0.0 )
        temp = x(j,j);
        x(j,j) = qraux(j);
        t = - ( qy(j:n) * conj ( x(j:n,j) ) ) / x(j,j);
        qy(j:n) = qy(j:n) + t * transpose ( x(j:n,j) );
        x(j,j) = temp;
      end

    end

  end
%
%  Compute hermitian ( A ) * Y.
%
  if ( cqty )
    for j = 1 : ju
      if ( zabs1 ( qraux(j) ) ~= 0.0 )
        temp = x(j,j);
        x(j,j) = qraux(j);
        t = - ( qty(j:n) * conj ( x(j:n,j) ) ) / x(j,j);
        qty(j:n) = qty(j:n) + t * x(j:n,j);
        x(j,j) = temp;
      end
    end
  end
%
%  Set up to compute B, RSD, or XB.
%
  if ( cb )
    b(1:k) = qty(1:k);
  end

  kp1 = k + 1;

  if ( cxb )
    xb(1:k) = qty(1:k);
  end

  if ( cr & k < n )
    rsd(k+1:n) = qty(k+1:n);
  end

  if ( cxb )
    xb(k+1:n) = 0.0;
  end

  if ( cr )
    rsd(1:k) = 0.0;
  end
%
%  Compute B.
%
  if ( cb )

    for jj = 1 : k

      j = k - jj + 1;

      if ( zabs1 ( x(j,j) ) == 0.0 )
        info = j;
        break
      end

      b(j) = b(j) / x(j,j);

      if ( j ~= 1 )
        t = -b(j);
        b(1:j-1) = b(1:j-1) + t * x(1:j-1,j);
      end

    end

  end

  if ( cr | cxb )
%
%  Compute RSD or XB as required.
%
    for jj = 1 : ju

      j = ju - jj + 1;

      if ( zabs1 ( qraux(j) ) ~= 0.0 )

        temp = x(j,j);
        x(j,j) = qraux(j);

        if ( cr )
          t = - ( rsd(j:n) * conj ( x(j:n,j) ) ) / x(j,j);
          rsd(j:n) = rsd(j:n) + t * x(j:n,j);
        end

        if ( cxb )
          t = - ( xb(j:n) * conj ( x(j:n,j) ) ) / x(j,j);
          xb(j:n) = xb(j:n) + t * x(j:n,j);
        end

        x(j,j) = temp;

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%XMLDEMO3 Demonstrate how to convert an XMLtree in a simple structure
%
%   Description
%   This script demonstrates the use of the xmltree class to
%   convert an XMLtree (when possible) in a simple structure.
%   This can only be performed when the XML file is simple enough
%   (one element cannot have more than one chardata child).

%   Copyright (C) 2003  Guillaume Flandin

clc;
disp('This demonstration illustrates the use of the xmltree class')
disp('to convert an xmltree in a simple Matlab structure.')
disp(' ')
disp('Press any key to open the XML file ''example.xml''.')
pause;

disp('>> t = xmltree(''example.xml'')')
try,
	t = xmltree('example.xml');
catch,
	disp('Please first launch xmldemo1')
end

disp(' ')
disp('Press any key to display the XML file.')
pause;

save(t);

disp(' ')
disp('Convert the xmltree in a Matlab structure.')
disp(' ')
disp('Press any key to continue.')
disp(' ')
pause; clc;

disp('>> s = convert(tree);')
s = convert(t);

disp(' ')
disp('Display the content of the structure.')
disp(' ')
disp('Press any key to continue.')
disp(' ')
pause;

disp('>> s.entry');
s.entry
disp('>> s.entry.address');
s.entry.address
disp('>> s.entry.address.institute');
s.entry.address.institute

disp(' ');
disp('Press any key to end.')
pause; clc; close all;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% function [A,truth]=create_SBM(N,q,c,epsi,com_size)
%
% This creates an SBM graph with
% Inputs:
% - N nodes,
% - q communites of sizes listed in com_size,
% - an average degree of c,
% - and a difficulty \epsilon=epsi.
% Outputs:
% - G the graph structure
%
% Copyright (C) 2016 Nicolas Tremblay, Gilles Puy.
% This file is part of the CSCbox (Compressive Spectral Clustering toolbox)
%
% The CSCbox 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.
%
% The CSCbox 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/>.
%
% If you use this toolbox please kindly cite
%     N. Tremblay, G. Puy, R. Gribonval and P. Vandergheynst.
%     Compressive Spectral Clustering.
%     ArXiv e-prints:1602.02018 Feb. 2016.

function [G]=create_SBM(N,q,c,epsi,com_size)

if sum(com_size)~=N
    error('error in create_SBM.m : sum(com_size)~=N!')
end

pin=(q*c)/(N-q+(q-1)*epsi*N);
pout=(q*c*epsi)/(N-q+(q-1)*epsi*N);

rng('shuffle');

% total # of intra community edges is taken from Bernouilli distrib
I=[]; J=[]; truth=zeros(N,1);
for k=1:q
    truth(sum(com_size(1:k-1))+1:sum(com_size(1:k)))=k;
    Numedgesk=my_binornd(com_size(k).*(com_size(k)-1)/2,pin); %normal approx
    while Numedgesk>com_size(k).*(com_size(k)-1)/2
        Numedgesk=my_binornd(com_size(k).*(com_size(k)-1)/2,pin); %normal approx
    end
    edges = randperm(com_size(k).*(com_size(k)-1)/2,Numedgesk);
    %edges=datasample([1:com_size(k).*(com_size(k)-1)./2],Numedgesk,'Replace',false);
    [Inew, Jnew] = ind2sub4up(edges);
    I=[I;Inew+sum(com_size(1:k-1))];
    J=[J;Jnew+sum(com_size(1:k-1))];
end

% for k=1:q-1
%     for kk=k+1:q
%         Numedgesk=my_binornd(com_size(k)*com_size(kk),pout); %normal approx
%         edges = randperm(com_size(k)*com_size(kk),Numedgesk);
%         %Numedgesk=binornd(com_size(k)*com_size(kk),pout);
%         %edges=datasample([1:com_size(k)*com_size(kk)],Numedgesk,'Replace',false);
%         [Inew, Jnew] = ind2sub([com_size(k),com_size(kk)],edges);
%         I=[I;Inew'+sum(com_size(1:k-1))];
%         J=[J;Jnew'+sum(com_size(1:kk-1))];
%     end
% end

for k=1:q-1
    Numedgesk=my_binornd(com_size(k)*sum(com_size(k+1:end)),pout); %normal approx
    edges = randperm(com_size(k)*sum(com_size(k+1:end)),Numedgesk);
    %Numedgesk=binornd(com_size(k)*sum(com_size(k+1:end)),pout);
    %edges=datasample([1:com_size(k)*sum(com_size(k+1:end))],Numedgesk,'Replace',false);
    [Inew, Jnew] = ind2sub([com_size(k),sum(com_size(k+1:end))],edges);
    I=[I;Inew'+sum(com_size(1:k-1))];
    J=[J;Jnew'+sum(com_size(1:k))];
end

A=sparse(I,J,1,N,N);
A=A+A';
G=gsp_graph(A);
G.type='SBM_CSC';
G.truth = truth;
G.k=q;

repo: gspbox


text:
function [J,Msp,Nsp]=mtfftpt(data,tapers,nfft,t,f,findx)
% Multi-taper fourier transform for point process given as times
%
% Usage:
% [J,Msp,Nsp]=mtfftpt (data,tapers,nfft,t,f,findx) - all arguments required
% Input:
%       data        (struct array of times with dimension channels/trials;
%                   also takes in 1d array of spike times as a column vector)
%       tapers      (precalculated tapers from dpss)
%       nfft        (length of padded data)
%       t           (time points at which tapers are calculated)
%       f           (frequencies of evaluation)
%       findx       (index corresponding to frequencies f)
% Output:
%       J (fft in form frequency index x taper index x channels/trials)
%       Msp (number of spikes per sample in each channel)
%       Nsp (number of spikes in each channel)
if nargin < 6; error('Need all input arguments'); end;
if isstruct(data); C=length(data); else C=1; end% number of channels
K=size(tapers,2); % number of tapers
nfreq=length(f); % number of frequencies
if nfreq~=length(findx); error('frequency information (last two arguments) inconsistent'); end;
H=fft(tapers,nfft,1);  % fft of tapers
H=H(findx,:); % restrict fft of tapers to required frequencies
w=2*pi*f; % angular frequencies at which ft is to be evaluated
Nsp=zeros(1,C); Msp=zeros(1,C);
for ch=1:C;
  if isstruct(data);
     fnames=fieldnames(data);
     eval(['dtmp=data(ch).' fnames{1} ';'])
     indx=find(dtmp>=min(t)&dtmp<=max(t));
     if ~isempty(indx); dtmp=dtmp(indx);
     end;
  else
     dtmp=data;
     indx=find(dtmp>=min(t)&dtmp<=max(t));
     if ~isempty(indx); dtmp=dtmp(indx);
     end;
  end;
  Nsp(ch)=length(dtmp);
  Msp(ch)=Nsp(ch)/length(t);
  if Msp(ch)~=0;
      data_proj=interp1(t',tapers,dtmp);
      exponential=exp(-i*w'*(dtmp-t(1))');
      J(:,:,ch)=exponential*data_proj-H*Msp(ch);
  else
      J(1:nfreq,1:K,ch)=0;
  end;
end;


text:
function qwv_2d_test02 ( )

%*****************************************************************************80
%
%% QWV_2D_TEST02 tests QWV_2D for Padua points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 May 2014
%
%  Author:
%
%    John Burkardt
%
  a = -1.0;
  b = +1.0;
  c = -1.0;
  d = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWV_2D_TEST02:\n' );
  fprintf ( 1, '  Compute the weights associated with an interpolatory\n' );
  fprintf ( 1, '  quadrature rule defined by N=(T+1)*(T+2)/2 points,\n' );
  fprintf ( 1, '  exact for polynomials of total degree T or less.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X Interval = [%g,%g]\n', a, b );
  fprintf ( 1, '  Y Interval = [%g,%g]\n', c, d );

  for t = 0 : 10

    n = ( ( t + 1 ) * ( t + 2 ) ) / 2;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Degree T = %d\n', t );
    fprintf ( 1, '  Number of points N = %d\n', n );

    [ x, y ] = padua_point_set ( t );
%
%  Compute the weights.
%
    w = qwv_2d ( t, n, a, b, c, d, x, y );

    r8vec_print_16 ( n, w, '  Weights:' );

  end

  return
end


text:
function [f, df, ddf] = rosenbrock(x);

% rosenbrock.m This function returns the function value, partial derivatives
% and Hessian of the (general dimension) rosenbrock function, given by:
%
%       f(x) = sum_{i=1:D-1} 100*(x(i+1) - x(i)^2)^2 + (1-x(i))^2
%
% where D is the dimension of x. The true minimum is 0 at x = (1 1 ... 1).
%
% Carl Edward Rasmussen, 2001-07-21.

D = length(x);
f = sum(100*(x(2:D)-x(1:D-1).^2).^2 + (1-x(1:D-1)).^2);

if nargout > 1
  df = zeros(D, 1);
  df(1:D-1) = - 400*x(1:D-1).*(x(2:D)-x(1:D-1).^2) - 2*(1-x(1:D-1));
  df(2:D) = df(2:D) + 200*(x(2:D)-x(1:D-1).^2);
end

if nargout > 2
  ddf = zeros(D,D);
  ddf(1:D-1,1:D-1) = diag(-400*x(2:D) + 1200*x(1:D-1).^2 + 2);
  ddf(2:D,2:D) = ddf(2:D,2:D) + 200*eye(D-1);
  ddf = ddf - diag(400*x(1:D-1),1) - diag(400*x(1:D-1),-1);
end

repo: pTVreg


text:
function printFuncName = getStandardOutputFuncName()
global logger

if isempty(logger)
    printFuncName = 'fprintf';
elseif isa(logger, 'Logger')
    if ~logger.IsOpen()
        logger.Open();
    end
    if ~logger.IsOpen()
        printFuncName = 'fprintf';
    else
        printFuncName = 'logger.Write';
    end
else
    printFuncName = 'logger.Write';
end

author: BUNPC