text:
function lbl = removeBorderRegions(lbl, varargin)
%REMOVEBORDERREGIONS Remove regions on the border of an image
%
%   LBL2 = removeBorderRegions(LBL);
%   LBL is a labeled image, or a binary image. Image can be 2D or 3D.
%   In the case of binary image, the image is labeled using 4-adjacency for
%   2D images, or 6-adjacency for 3D images.
%   The regions which touch the border of the image are removed, and the
%   remaining regions are re-labeled. Result has the same class than the
%   original image.
%
%   LBL2 = removeBorderRegions(LBL, BORDERS);
%   Specifies which borders need to be removed. BORDERS is a cell array of
%   strings, with each string corresponding to a border:
%   'left':     all pixels with x=1
%   'right':    all pixels with x=last
%   'top':      all pixels with y=1
%   'bottom':   all pixels with y=last
%
%   For 3D images, 2 additonal borders can be removed:
%   'front':    all pixels with z=1
%   'back':     all pixels with z=last
%
%   LBL2 = removeBorderRegions(LBL, BORDERS, RELABEL);
%   where RELABEL is a logical, specifies wheter remaining regions should
%   be relabeled or not (default is TRUE).
%
%   See also:
%   bwmorph
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 23/04/2007
%

% HISTORY :
%   19/07/2007: add possibility to specify borders
%   07/08/2007: relabeling is now optional
%   01/02/2008: return binImg image if input is binImg

warning('matImage:deprecated', ...
    'function "removeBorderRegion" has been replaced by "imKillBorders"');

% default values
binImg  = false;
relabel = true;

% If image is binImg, convert to labeled image
if islogical(lbl)
    binImg = true;
    relabel = false;
    if ndims(lbl)==2 %#ok<ISMAT>
        lbl = bwlabel(lbl, 4);
    elseif ndims(lbl)==3
        lbl = bwlabeln(lbl, 6);
    end
end

% check if relabel flag is present
if ~isempty(varargin)
    if islogical(varargin{end})
        relabel = varargin{end};
        varargin(end) = [];
    end
end

% find borders to remove
borders = {'top', 'bottom', 'left', 'right'};
if ~isempty(varargin)
    if iscell(varargin{1})
        borders = varargin{1};
    else
        borders = varargin;
    end
end

% create mask for borders
mask = false(size(lbl));
if ndims(lbl)==2 %#ok<ISMAT>
    for i=1:length(borders)
        border = borders{i};
        if strcmp(border, 'left')
            mask(:,1) = true;
        elseif strcmp(border, 'right')
            mask(:,end) = true;
        elseif strcmp(border, 'top')
            mask(1,:) = true;
        elseif strcmp(border, 'bottom')
            mask(end,:) = true;
        end
    end
else
    for i=1:length(borders)
        border = borders{i};
        if strcmp(border, 'left')
            mask(:,1,:) = true;
        elseif strcmp(border, 'right')
            mask(:,end,:) = true;
        elseif strcmp(border, 'top')
            mask(1,:,:) = true;
        elseif strcmp(border, 'bottom')
            mask(end,:,:) = true;
        elseif strcmp(border, 'front')
            mask(:,:,1) = true;
        elseif strcmp(border, 'back')
            mask(:,:,end) = true;
        end
    end
end

% find labels of regions touching borders
lbl2 = unique(lbl(mask));

% remove label '0', which corresponds to background
if lbl2(1)==0
    lbl2 = lbl2(2:length(lbl2));
end

% set labels of border regions to 0
indices = ismember(lbl, lbl2);
lbl(indices)=0;

% Eventually relabel remaining regions by shifting labels
if relabel && ~isempty(lbl2)
    % For remaining regions, shift the labels
    for i=1:length(lbl2)-1
        indices = lbl>lbl2(i) & lbl<=lbl2(i+1);
        lbl(indices) = imsubtract(lbl(indices), i);
    end

    % Shift labels for regions whose index is greater than last index of LBL2.
    indices = lbl>lbl2(length(lbl2));
    lbl(indices) = imsubtract(lbl(indices), length(lbl2));
end

% convert to binImg image if needed
if binImg
    lbl = lbl>0;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

m = size(X  , 1);

for i = 1 : m,
	T = [];
	for j = 1 : K,
		T = [T ; X(i,:)];
	end
	[Max , idx(i)] = min(sum((T - centroids).^2 , 2));
end

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function linplus_test60 ( )

%*****************************************************************************80
%
%% TEST60 tests R8UT_DET, R8UT_INVERSE, R8UT_MXM, R8UT_RANDOM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST60\n' );
  fprintf ( 1, '  For an upper triangular matrix,\n' );
  fprintf ( 1, '  R8UT_DET computes the determinant.\n' );
  fprintf ( 1, '  R8UT_INVERSE computes the inverse.\n' );
  fprintf ( 1, '  R8UT_MXM computes matrix products.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  [ a, seed ] = r8ut_random ( n, n, seed );

  r8ut_print ( n, n, a, '  The matrix A:' );
%
%  Compute the determinant.
%
  det = r8ut_det ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Determinant is %f\n', det );
%
%  Compute the inverse matrix B.
%
  b = r8ut_inverse ( n, a );

  r8ut_print ( n, n, b, '  The inverse matrix B:' );
%
%  Check
%
  c = r8ut_mxm ( n, a, b );

  r8ut_print ( n, n, c, '  The product A * B:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lpcb
global lowpass hlp

if get(hlp,'value')
    lowpass=true;
else
    lowpass=false;
end

author: Sable

repo: mcbench-benchmarks


text:
function weights = get_w(a)

%
%               w = get_w()
%
%   Returns the weight vector from the svm
%   AE: i have changed this procedure to include non linear feature
%   selection with RFE. When non-linear, it gives  a score vector for each
%   feature (used in RFE) which is not the margin but is sorted as the margin

weights=a.w;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef Jab < handle
  % .jab files

  methods (Static)

    function J = load(jabfile)
      % J = load(jabfile)
      %
      % J: Modernized Macguffin object

      assert(ischar(jabfile) && exist(jabfile,'file')==2,...
        'Cannot find jabfile ''%s''.',jabfile);

      J = loadAnonymous(jabfile);
      if isstruct(J)
        J = Macguffin(J);
      end
      J.modernize();
    end

    function save(J,jabfile)
      saveAnonymous(jabfile,J);
    end

    function savePrompt(J,jabfile)
      if exist('jabfile','var')==0
        jabfile = [];
      end
      [fname,pname] = uiputfile('*.jab','Save',jabfile);
      fname = fullfile(pname,fname);
      Jab.save(J,fname);
    end

    function [tfsuccess,jabfiles] = uiGetJabFiles(varargin)
      % [tfsuccess,jabfiles] = uiGetJabFiles(p1,v1,...)
      %
      % Nice(r) GUI for multiple jab selection.
      %
      % Optional PVs:
      %  * promptstr, char.
      %
      % tfsuccess: if false, user canceled GUI
      % jabfiles: cellstr

      promptstr = myparse(varargin,...
        'promptstr','Select jab files');

      %jabpath = ExpPP.loadConfigVal('jabpath');
      jabpath = [];
      if isempty(jabpath)
        jabpath = pwd;
      end
      jabfiles = uipickfiles('Prompt',promptstr,...
        'FilterSpec',jabpath,'Type',{'*.jab','JAB-files'},'SmartAdd',true);
      if ~iscell(jabfiles) || isempty(jabfiles), % canceled
        tfsuccess = false;
        jabfiles = [];
      else
        tfsuccess = true;
        %jabpath = fileparts(jabfiles{1});
        %ExpPP.saveConfigVal('jabpath',jabpath);
      end
    end

    function clearLabels(jabfile,realbeh)
      % Clear all labels for realbeh and No-realbeh

      J = Jab.load(jabfile);

      jabrealbehs = J.behaviors.names(1:J.behaviors.nbeh);
      tf = strcmpi(realbeh,jabrealbehs);
      assert(any(tf),'Specified behavior ''%s'' not present in jabfile.',realbeh);
      assert(nnz(tf)==1);
      realbeh = jabrealbehs{tf}; % case could be different

      tfMultiCls = J.behaviors.nbeh>1;
      if tfMultiCls
        nobeh = Labels.noBehaviorName(realbeh);
      else
        nobeh = 'None';
      end

      J.labels = Labels.clearLabels(J.labels,realbeh,realbeh);
      J.labels = Labels.clearLabels(J.labels,nobeh,realbeh);

      J.savePrompt(J,jabfile);
    end

    function tf = isMultiClassifier(jabfile)
      % tf = isMultiClassifier(jabfile)
      J = Jab.load(jabfile);
      tf = J.isMultiClassifier();
    end

    function rmClassifier(jabfile,clsnames)
      % rmClassifier(jabfile,clsnames)
      %
      % jab: char, jab filename
      % clsname: char or cellstr

      J = Jab.load(jabfile);
      Jbeh = J.behaviors;
      Jclsnames = Jbeh.names(1:Jbeh.nbeh);

      if exist('clsnames','var')>0
        if ischar(clsnames)
          clsnames = cellstr(clsnames);
        end
        assert(iscellstr(clsnames),'Expected ''clsname'' argument to be a char or cellstr.');
      else
        [sel,ok] = listdlg('PromptString','Select classifiers to remove:',...
          'ListString',Jclsnames,...
          'ListSize',[240 240]);
        if ~ok
          return;
        end
        clsnames = Jclsnames(sel);
      end

      [tf,iClsRm] = ismember(clsnames,Jclsnames);
      if ~all(tf)
        clsmissing = clsnames(~tf);
        error('Jab:rmClassifier','Classifier(s) %s not present in jabfile.',...
          civilizedStringFromCellArrayOfStrings(clsmissing));
      end

      nClsRm = numel(iClsRm);
      nClsNew = Jbeh.nbeh - nClsRm;
      if nClsNew==0
        error('Jab:rmClassifier','Cannot remove all classifiers from a .jab file.');
      end

      iLblsRm = [iClsRm iClsRm+Jbeh.nbeh];
      origBehs = J.behaviors.names(iClsRm);
      origNoBehs = J.behaviors.names(iClsRm+Jbeh.nbeh);
      J.behaviors.names(iLblsRm) = [];
      J.behaviors.labelcolors(iLblsRm,:) = [];
      J.behaviors.nbeh = nClsNew;
      J.file.scorefilename(iClsRm) = [];
      J.windowFeaturesParams(iClsRm) = [];
      for iCls = 1:nClsRm
        J.labels = Labels.removeClassifier(J.labels,origBehs{iCls},origNoBehs{iCls});
        J.gtLabels = Labels.removeClassifier(J.gtLabels,origBehs{iCls},origNoBehs{iCls});
      end
      J.classifierStuff(iClsRm) = [];

      % Quirk, no-beh name for single-classifier projs
      if nClsNew==1
        oldNoBehavior = J.behaviors.names{2};
        newNoBehavior = 'None';
        J.behaviors.names{2} = newNoBehavior;

        J.labels = Labels.renameBehaviorRaw(J.labels,oldNoBehavior,newNoBehavior);
        J.gtLabels = Labels.renameBehaviorRaw(J.gtLabels,oldNoBehavior,newNoBehavior);

        % ALXXX to be fixed
        warning('Jab:rmClassifier',...
          'Going from multi-classifier to single-classifier project. Label importance may not be set correctly for GT mode.');
      end

      Jab.savePrompt(J,jabfile);
    end

    function tf = formatChangedBetweenVersions(v1,v2)
      % v2 should be later-than-or-equal-to-v1.
      %
      % This method introduced in version 0.6.0. Jabs saved in prior
      % versions have version prop enclosed in double-cell array. We avoid
      % the complications for now and just return true if v2 is 0.6.0 and
      % v1 is not. This will need to be updated with subsequent versions.

      tf = strcmp(v2,'0.6.0') && ~strcmp(v1,'0.6.0');
    end

    function merge(jabfiles,jabout)
      % jabMerge(jabfiles,jabout)
      % Merge two or more jab files.
      %
      % Behaviors with the same name (precisely, case-sensitive, etc) are
      % treated as the same classifier. Experiment names/directories should
      % be in a consistent format.
      %
      % jabfiles: optional. cellstr of jab filenames
      % jabout: optional. output jab filename

      if ~exist('jabfiles','var') || isempty(jabfiles)
        [tfsuccess,jabfiles] = Jab.uiGetJabFiles('promptstr','Select jab files to combine');
        if ~tfsuccess
          return;
        end
      end

      if ischar(jabfiles)
        jabfiles = cellstr(jabfiles);
      end
      assert(iscellstr(jabfiles),'Expected ''jabfiles'' to be a cellstr of jab filenames.');

      Js = cellfun(@Jab.load,jabfiles,'uni',0);
      Js = cat(1,Js{:});
      Jmerge = Macguffin(Js);

      % come up with a proposed name for combined jab
      if ~exist('jabout','var') || isempty(jabout)
        MAXFILENAMELENGTH = 45;
        behnames = Labels.verifyBehaviorNames(Jmerge.behaviors.names);
        combjabname = '';
        for i = 1:numel(behnames)
          combjabname = [combjabname behnames{i} '.']; %#ok<AGROW>
          if numel(combjabname)>MAXFILENAMELENGTH
            combjabname = [combjabname 'etc.']; %#ok<AGROW>
            break;
          end
        end
        combjabname = [combjabname 'jab'];

        %jabpath = ExpPP.loadConfigVal('jabpath');
        jabpath = [];
        if isempty(jabpath)
          jabpath = pwd;
        end
        [filename,pathname] = ...
          uiputfile({'*.jab','JAABA files (*.jab)'},'Save combined jabfile',fullfile(jabpath,combjabname));
        if ~ischar(filename),
          % user hit cancel
          return;
        end

        jabout = fullfile(pathname,filename);
      else
        if exist(jabout,'file')
          uiwait(warndlg(sprintf('Output file ''%s'' exists and will be overwritten.',jabout)));
        end
      end

      Jab.save(Jmerge,jabout);
      %ExpPP.saveConfigVal('jabpath',fileparts(jabout));
    end

    function [scorefiles,behaviors] = jabfileScoresBehs(jabname)
      % [scorefiles,behaviors] = jabfileScoresBehs(jabname)
      % scorefile: cellstr of score files
      % behaviors: cellstr of "real" behaviors (doesn't include No_<beh>)

      x = loadAnonymous(jabname);
      scorefiles = x.file.scorefilename;
      if ischar(scorefiles)
        scorefiles = cellstr(scorefiles);
      end
      behaviors = Labels.verifyBehaviorNames(x.behaviors.names);
      assert(numel(scorefiles)==numel(behaviors));
    end

  end

end

author: kristinbranson

repo: JAABA


text:
% This script assumes that the variables under interest of analysis are:
%
% 1) SAVED FROM SIMULINK as STRUCTURE WITH TIME format;
%
% 2) contains the suffix "_DATA" in the name.
%
% Then, the script reads the workspace (or loads a .mat file), looks for
% these variables and allows the user to choose which variable to plot.
%
%

%% DATA TO LOAD (comment if data are in the workspace)
% load('EXAMPLE.mat')

%% SETUP VISUALIZER
close all
clc

addpath('./src_plot')

%%utility variables

% labels and font size
xAxisLabel  = 'Time [s]';
lineSize    = 2;

% print figures to .png
printFigure = false;

% dock all figures together
dockFigures = true;

% desired init and end time in the figures [s]
timeInit    = 0;
timeLimit   = 120;

if dockFigures

    set(0,'DefaultFigureWindowStyle','docked')
else
    set(0,'DefaultFigureWindowStyle','normal')
end

%% Get dataset from workspace
dataVariables = whos('*DATA');
plotList      = {};

for k = 1:length(dataVariables)

    plotList{k} = dataVariables(k).name; %#ok<SAGROW>
end

[plotNumber, ~] = listdlg('PromptString',  'Choose the variables to plot:', ...
                          'ListString',     plotList, ...
                          'SelectionMode', 'multiple', ...
                          'ListSize',       [250 150]);

% regenerate the plot list with the data selected by user
plotList = plotList(plotNumber);

%% Select which time to use for the plot (if Yarp time is available)
timeList = {'simulink time'};

if exist('time_Yarp','var')

    timeList{2} = 'yarp time';
end

if length(timeList) > 1

   [timeNumber, ~] = listdlg('PromptString',  'Choose the time vector to use:', ...
                             'ListString',     timeList, ...
                             'SelectionMode', 'single', ...
                             'ListSize',       [250 150]);
else
    timeNumber = 1;
end

%% Process data to plot
figNumber = 0;

for k = 1:length(plotList)

    currentData = eval(plotList{k});

    if timeNumber == 1

        time = currentData.time;

    elseif timeNumber == 2

        time = time_Yarp.signals.values;
    end

    % get the user defined time length
    timeLimitIndex = sum(time <= timeLimit);
    timeInitIndex  = sum(time <= timeInit);
    tVector        = time(timeInitIndex:timeLimitIndex);

    dataStructure  = currentData.signals;
    figNumber      = figNumber + 1;

    for kk = 1:length(currentData.signals)

        valuesToPlot = squeeze(currentData.signals(kk).values);

        % twist data dimensions if necessary
        if size(valuesToPlot,1) < length(time)

            valuesToPlot = valuesToPlot';
        end

        valuesToPlot = valuesToPlot(timeInitIndex:timeLimitIndex,:);
        yAxisLabel   = currentData.signals(kk).label;

        if kk == 1

            titleLabel = plotList{k}(1:end-5);
        else
            titleLabel = '';
        end

        figure(figNumber)
        subplot(length(currentData.signals),1,kk)
        custom_plot(tVector, valuesToPlot, xAxisLabel, yAxisLabel, titleLabel, ...
                    lineSize, '', figNumber, printFigure, 0);
    end
end

rmpath('./src_plot')

author: robotology

repo: whole-body-controllers


text:
clear all;
close all;
clc
%%
load('ExperimentInformation.mat');
load('GridCellAnalysis.mat');
load('NAT.mat');
%% paremater settings
AngleSmooth=2;
AngleBinsize=3; % cm
SpeedThreadhold=2.5; %cm/mm
MinEventCount=100;
MinSNR=3;
Shuffling=1000;
Shuffling_mininterval=30; %second
%% shuffling
SelectedFrame_raw=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell); %% which frame in each session should be used to do analyses for each cell
SelectedFrame_filtered=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell); %% which frame in each session should be used to do analyses for each cell, filtered by some cratirial
TurningCurve_whole=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell);
TurningCurve_firsthalf=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell);
TurningCurve_secondthalf=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell);
TurningCurveStat_whole=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell);
TurningCurveStat_firsthalf=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell);
TurningCurveStat_secondthalf=cell(ExperimentInformation.Session,ExperimentInformation.TotalCell);
CellForAnalysis=cell(ExperimentInformation.Session,1);
MVL_shuffled=zeros(ExperimentInformation.TotalCell,Shuffling+3,ExperimentInformation.Session);
Correlation_shuffled=zeros(ExperimentInformation.TotalCell,Shuffling+3,ExperimentInformation.Session);

for j=1:1:ExperimentInformation.Session
    k=1;
    for i=1:1:ExperimentInformation.TotalCell
       SelectedFrame_raw=find(~isnan(NAT{1,j}(:,4*i+10)));
        SelectedFrame_filtered=intersect(find(~isnan(NAT{1,j}(:,4*i+10))),find(NAT{1,j}(:,6)==1));% filter out the frames with speed valid
        SelectedFrame_filtered=intersect(SelectedFrame_filtered,find(NAT{1,j}(:,5)>SpeedThreadhold));% filter out the frames with speed threadhold
        SelectedFrame_filtered_firstHalf=SelectedFrame_filtered(find(SelectedFrame_filtered<=ExperimentInformation.FrameinEachSession(j)/2));
        SelectedFrame_filtered_secondHalf=SelectedFrame_filtered(find(SelectedFrame_filtered>ExperimentInformation.FrameinEachSession(j)/2));
        HeadDirectionTrain=NAT{1,j}(SelectedFrame_filtered,4);
        HeadDirectionTrain_firstHalf=NAT{1,j}(SelectedFrame_filtered_firstHalf,4);
        HeadDirectionTrain_secondHalf=NAT{1,j}(SelectedFrame_filtered_secondHalf,4);
        Event_filtered=intersect(find(NAT{1,j}(:,4*i+12)>0),find(NAT{1,j}(:,6)==1));% filter out the frames with speed valid
        Event_filtered=intersect(Event_filtered,find(NAT{1,j}(:,5)>SpeedThreadhold));% filter out the frames with speed threadholds
        Event_firstHalf=NAT{1,j}(SelectedFrame_filtered_firstHalf,[1 4*i+12]);
        Event_secondHalf=NAT{1,j}(SelectedFrame_filtered_secondHalf,[1 4*i+12]);
        if length(Event_filtered)>MinEventCount && ExperimentInformation.CellSNR(i)>3 && ~ismember(i,ExperimentInformation.RepeatCell)
            EventTrain=NAT{1,j}(SelectedFrame_filtered,[1 4*i+12]);
            Event_firstHalf=NAT{1,j}(SelectedFrame_filtered_firstHalf,[1 4*i+12]);
            Event_secondHalf=NAT{1,j}(SelectedFrame_filtered_secondHalf,[1 4*i+12]);
            CellForAnalysis{j,1}(k)=i;
            for n=1:1:Shuffling
                xmin=Shuffling_mininterval*ExperimentInformation.Trackingframerate;
                xmax=size(length(SelectedFrame_filtered),1)-xmin;
                ShiftFrame=round(xmin+rand(1,1)*(xmax-xmin));
                EventTrain_shuffled=circshift(EventTrain,ShiftFrame);
                EventTrain_shuffled(:,1)=EventTrain(:,1);
                EventTrain_shuffled_firstHalf=EventTrain_shuffled(1:size(HeadDirectionTrain_firstHalf,1),:);
                EventTrain_shuffled_secondHalf=EventTrain_shuffled(size(HeadDirectionTrain_firstHalf,1)+1:end,:);
                TurningCurve_shuffled=SpatialTuning_BNT.TurnningCurver_power(HeadDirectionTrain,...
                    EventTrain_shuffled(:,2),...
                    AngleBinsize,...
                    1/(ExperimentInformation.Trackingframerate/ExperimentInformation.ImagingPlane),...
                    AngleSmooth);
                TurnningCurveStastics_shuffled=SpatialTuning_BNT.tcStatistics(TurningCurve_shuffled,AngleBinsize, 50);
                MVL_shuffled(i,n,j)=TurnningCurveStastics_shuffled.r;

                TurningCurve_firstHalf_shuffled=SpatialTuning_BNT.TurnningCurver_power(HeadDirectionTrain_firstHalf,...
                    EventTrain_shuffled_firstHalf(:,2),...
                    AngleBinsize,...
                    1/(ExperimentInformation.Trackingframerate/ExperimentInformation.ImagingPlane),...
                    AngleSmooth);
                %                 TurnningCurveStastics_firstHalf=SpatialTuning_BNT.tcStatistics(TurningCurve_firstHalf,AngleBinsize, 50);

                TurningCurve_secondHalf_shuffled=SpatialTuning_BNT.TurnningCurver_power(HeadDirectionTrain_secondHalf,...
                    EventTrain_shuffled_secondHalf(:,2),...
                    AngleBinsize,...
                    1/(ExperimentInformation.Trackingframerate/ExperimentInformation.ImagingPlane),...
                    AngleSmooth);
                %                 TurnningCurveStastics_secondHalf=SpatialTuning_BNT.tcStatistics(TurningCurve_secondHalf,AngleBinsize, 50);

                Correlation_shuffled(i,n,j)=corr(TurningCurve_firstHalf_shuffled(:,2),TurningCurve_secondHalf_shuffled(:,2));
                MVL_shuffled(i,Shuffling+1,j)=prctile(MVL_shuffled(i,1:Shuffling,j),95);
                MVL_shuffled(i,Shuffling+2,j)=prctile(MVL_shuffled(i,1:Shuffling,j),99);
                Correlation_shuffled(i,Shuffling+1,j)=prctile(Correlation_shuffled(i,1:Shuffling,j),95);
                Correlation_shuffled(i,Shuffling+2,j)=prctile(Correlation_shuffled(i,1:Shuffling,j),99);
                TurningCurve_whole{j,i}=SpatialTuning_BNT.TurnningCurver_power(HeadDirectionTrain,...
                    EventTrain(:,2),...
                    AngleBinsize,...
                    1/(ExperimentInformation.Trackingframerate/ExperimentInformation.ImagingPlane),...
                    AngleSmooth);
                TurningCurveStat_whole{j,i}=SpatialTuning_BNT.tcStatistics(TurningCurve_whole{j,i}, AngleBinsize, 50);
                MVL_shuffled(i,Shuffling+3,j)=TurningCurveStat_whole{j,i}.r;
                TurningCurve_firsthalf{j,i}=SpatialTuning_BNT.TurnningCurver_power(HeadDirectionTrain_firstHalf,...
                    Event_firstHalf(:,2),...
                    AngleBinsize,...
                    1/(ExperimentInformation.Trackingframerate/ExperimentInformation.ImagingPlane),...
                    AngleSmooth);
                TurningCurveStat_firsthalf{j,i}=SpatialTuning_BNT.tcStatistics(TurningCurve_firsthalf{j,i}, AngleBinsize, 50);
                TurningCurve_secondthalf{j,i}=SpatialTuning_BNT.TurnningCurver_power(HeadDirectionTrain_secondHalf,...
                    Event_secondHalf(:,2),...
                    AngleBinsize,...
                    1/(ExperimentInformation.Trackingframerate/ExperimentInformation.ImagingPlane),...
                    AngleSmooth);
                TurningCurveStat_secondthalf{j,i}=SpatialTuning_BNT.tcStatistics(TurningCurve_secondthalf{j,i},AngleBinsize, 50);
                Correlation_shuffled(i,Shuffling+3,j)=corr(TurningCurve_firsthalf{j,i}(:,2),TurningCurve_secondthalf{j,i}(:,2));
                disp([num2str(i),'-',num2str(n)]);
            end
            k=k+1;
        else
        end
    end
end
%% identify HD cells
IsHDCell=cell(ExperimentInformation.Session,1);
for j=1:1:ExperimentInformation.Session
    g=1;
    for i=1:1:ExperimentInformation.TotalCell
            if MVL_shuffled(i,Shuffling+3,j)>MVL_shuffled(i,Shuffling+1,j) && Correlation_shuffled(i,Shuffling+3,j)>Correlation_shuffled(i,Shuffling+1,j)
                IsHDCell{j,1}(g)=i;
                g=g+1;
            else
            end
    end
end
%% save HD cell information
HDCellAnalysis=struct;
HDCellAnalysis.IsHDCell=IsHDCell;
HDCellAnalysis.TurningCurve_whole=TurningCurve_whole;
HDCellAnalysis.TurningCurve_firsthalf=TurningCurve_firsthalf;
HDCellAnalysis.TurningCurve_secondthalf=TurningCurve_secondthalf;
HDCellAnalysis.TurningCurveStat_whole=TurningCurveStat_whole;
HDCellAnalysis.TurningCurveStat_firsthalf=TurningCurveStat_firsthalf;
HDCellAnalysis.TurningCurveStat_secondthalf=TurningCurveStat_secondthalf;
HDCellAnalysis.CellForAnalysis=CellForAnalysis;
HDCellAnalysis.MVL_shuffled=MVL_shuffled;
HDCellAnalysis.Correlation_shuffled=Correlation_shuffled;
HDCellAnalysis.Shuffling=Shuffling;
HDCellAnalysis.AngleSmooth=AngleSmooth;
HDCellAnalysis.AngleBinsize=AngleBinsize; % cm
HDCellAnalysis.SpeedThreadhold=SpeedThreadhold; %cm/mm
HDCellAnalysis.MinEventCount=MinEventCount;
HDCellAnalysis.MinSNR=MinSNR;
HDCellAnalysis.Shuffling_mininterval=Shuffling_mininterval; %second
%% plot  whole-session calcium activity maps, colored by head-direction, for all HD cells
close all
figure
x0=-100;
y0=-100;
width=3440;
height=1200;
set(gcf,'position',[x0,y0,width,height])
j=1;
coluum=20;
S=0;
for k=1:1:size(HDCellAnalysis.IsHDCell{j,1},2)
    i=HDCellAnalysis.IsHDCell{j,1}(k);
    P=mod(k,coluum);
    if P==0
        P=coluum;

    elseif P==1
        S=S+1;
    else
    end
    SelectedFrame_filtered=intersect(find(~isnan(NAT{1,j}(:,4*i+10))),find(NAT{1,j}(:,6)==1));% filter out the frames with speed valid
    SelectedFrame_filtered=intersect(SelectedFrame_filtered,find(NAT{1,j}(:,5)>2.5));% filter out the frames with speed threadhold
    Event_filtered=intersect(find(NAT{1,j}(:,4*i+12)>0),find(NAT{1,j}(:,6)==1));% filter out the frames with speed valid
    Event_filtered=intersect(Event_filtered,find(NAT{1,j}(:,5)>2.5));% filter out the frames with speed threadhold
    Bestshift=GridCellAnalysis.GridBestShift(i);
    Position=NAT{1,j}(Event_filtered,2:3);
    Position(:,1)=Position(:,1)+ (Bestshift*cos(NAT{1,j}(Event_filtered,4) * pi/180));
    Position(:,2)=Position(:,2)+ (Bestshift*sin(NAT{1,j}(Event_filtered,4) * pi/180));

    Event=NAT{1,j}(Event_filtered,4*i+12);
    Direction=NAT{1,j}(Event_filtered,4);
    Max=max(Event(:));
    scatter(Position(:,1)+100*(P-1),Position(:,2)-100*(S-1),10*(sqrt(Event./Max)),Direction,'filled','MarkerFaceAlpha',0.8)
    colormap(gca,hsv)
    caxis([0 360]);
%     colorbar;
    ylim([-100*(round(size(HDCellAnalysis.IsHDCell{j,1},2)/coluum)+0.5),50])
    xlim([-50,100*coluum])
    daspect([1 1 1]);
    hold on
    axis off
end
set(gca,'color',[1 1 1]);
set(gcf,'color',[1 1 1]);

%% plot  whole-session HD tuning for all HD cells
close all
figure
x0=10;
y0=10;
width=4000;
height=1400;
set(gcf,'position',[x0,y0,width,height])
j=1;
for k=1:1:size(HDCellAnalysis.IsHDCell{j,1},2)
    i=HDCellAnalysis.IsHDCell{j,1}(k);
    subplot(9,20,k,'align')
    Occupancy=preprocessing.circularSmooth(TurningCurve_whole{j,i}(:,3),AngleSmooth);
    Tuning=TurningCurve_whole{j,i}(:,2);
    MaxTime=max(Occupancy);
    MaxTuning=max(Tuning);
    Occupancy=0.5*Occupancy./MaxTime;
    Tuning=Tuning./MaxTuning;
    H=polarplot([TurningCurve_whole{j,i}(:,1)/180*pi;TurningCurve_whole{j,i}(1,1)/180*pi],[Occupancy;Occupancy(1)],'color',[0.6 0.6 0.6],'LineWidth',1);
    hold on
    H=polarplot([TurningCurve_whole{j,i}(:,1)/180*pi;TurningCurve_whole{j,i}(1,1)/180*pi],[Tuning;Tuning(1)],'color',[0 0 0],'LineWidth',1.5);
    Ax = gca;
    Ax.RTick = [];
    Ax.RTickLabel = [];
    Ax.ThetaTickLabel = [];
    Ax.ThetaTick  = [0 90 180 270];
    Ax.RLim=[0 1];
    AX.FontSize=0;
    Ax.ThetaAxis.Color = [1 1 1];
    AX.RColor = [1 1 1];
    set(gca,'GridColor',[0 0 0]);
    set(gca,'GridAlpha',1);
    set(gca,'RGrid','off');
    set(H,'LineWidth',1);
    AX.RAxis.Visible='off';
%     title(['P-time: ',num2str(MaxTime,'%.2f'),' s',', P-tuning: ',num2str(MaxTuning,'%.2f')])
end
set(gca,'color',[1 1 1]);
set(gcf,'color',[1 1 1]);

%%
save ([ExperimentInformation.RawDataAddress,'\HDCellAnalysis.mat'],'HDCellAnalysis','-v7.3');
disp('HD cell analysis was done!');

author: kavli-ntnu


text:
%% A 10x10 tridiagonal matrix, with 2 on the diagonal, -1 on the off diagonal.

A = blktridiag(2,-1,-1,10);

% The sparsity pattern is correct
spy(A)

% and the elements are as designated
full(A)

%% A lower block bidiagonal matrix with replicated blocks

% with 2x2 blocks of ones on the main diagonal, and
% 2x2 blocks of twos on the sub-diagonal

A = blktridiag(ones(2),2*ones(2),zeros(2),5);

spy(A)
full(A)

%% A block tridiagonal matrix with replicated blocks

Amd = reshape(1:9,3,3);
Asub = reshape(11:19,3,3);
Asup = reshape(21:29,3,3);
A = blktridiag(Amd,Asub,Asup,4);

spy(A)
full(A)

%% A tridiagonal matrix with random elements

Amd = rand(1,1,7);
Asub = rand(1,1,6);
Asup = rand(1,1,6);
A = blktridiag(Amd,Asub,Asup);

spy(A)
full(A)

%% A block tridiagonal matrix with distinct elements

Amd = reshape(1:27,[3 3 3]);
Asub = reshape(101:118,[3 3 2]);
Asup = reshape(201:218,[3 3 2]);
A = blktridiag(Amd,Asub,Asup);

spy(A)
full(A)

%% A block tridiagonal matrix with 2x3 fixed non-square blocks

Amd = rand(2,3);
Asub = 2*ones(2,3);
Asup = ones(2,3);
A = blktridiag(Amd,Asub,Asup,3);

spy(A)
full(A)

%% A block tridiagonal matrix with varying 3x2 non-square blocks
Amd = reshape(1:18,[3 2 3]);
Asub = reshape(101:112,[3 2 2]);
Asup = reshape(201:212,[3 2 2]);
A = blktridiag(Amd,Asub,Asup);

spy(A)
full(A)


text:
function mesh_write_fs_curv(p)

% mesh_write_fs_curv - Save mesh CData to FreeSurfer curvature file
%
% USEAGE: mesh_write_fs_curv(p)
%
% Write a binary curvature file for each mesh in p.mesh.data.
% If any cells of p.mesh.data.meshtype are empty, they are skipped.
%
% This function outputs an *.curv curvature file from the
% vertex scalar data contained in p.mesh.data.Cdata - if it
% contains a timeseries, the output values are taken from
% the column specified in p.mesh.samplePoint.
%
% All output is handled by fs_write_curv.
%
% See the FreeSurfer website at http://surfer.nmr.mgh.harvard.edu/
% for more information.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:58 $

% Licence:  GNU GPL, no implied or express warranties
% History:  10/2002, Darren.Weber_at_radiology.ucsf.edu
%           01/2003, Darren.Weber_at_radiology.ucsf.edu
%                    now using fs_write_surf and fs_write_wfile
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\nMESH_WRITE_FS_CURV...\n');

if ~exist('p','var'),
    error('...no input p struct.\n');
elseif isempty(p),
    error('...input p struct is empty.\n');
elseif isempty(p.mesh.data),
    error('...input p struct has no mesh data.\n');
end

[path,name,ext] = fileparts(strcat(p.mesh.path,filesep,p.mesh.file));
file = fullfile(path,[name ext]);

fprintf('...writing FreeSurfer curvatures to:\n\t%s\n',path);

tic;

Meshes = p.mesh.data.meshtype;

[path,name,ext] = fileparts(file);

for i=1:size(Meshes,2),

    if Meshes{i},

        % check for mesh Cdata
        if size(p.mesh.data.Cdata,2) >= i,
            if ~isempty(p.mesh.data.Cdata{1}),
                if size(p.mesh.data.Cdata{i},2) > 1,
                    % Obtain the Cdata at the selected time point
                    curv = p.mesh.data.Cdata{i}(:,p.mesh.samplePoint);
                else
                    curv = p.mesh.data.Cdata{i};
                end

                Nfaces = size(p.mesh.data.faces{i},1);

                outputExt = [ext,'.curv'];
                outputName = strcat(name,'.',Meshes{i});
                outputFile = fullfile(path,[outputName outputExt]);

                fprintf('...writing curvature: %s\n',[outputName outputExt]);
                fs_write_curv(outputFile,curv,Nfaces);
            else
                msg = sprintf('\np.mesh.data.Cdata{%d} is empty\n',i);
                warning(msg);
            end
        else
            msg = sprintf('\nno p.mesh.data.Cdata{%d}\n',i);
            warning(msg);
        end
    end
end

t=toc; fprintf('...done (%5.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
function s = structmerge(varargin)
% s = structmerge(s1,s2,s3,...)
% Merge structures into one combined structure.

s = struct();
for c = 1:numel(varargin)
    stmp = varargin{c};
    assert(isscalar(stmp) && isstruct(stmp),...
        'All input arguments must be scalar structures.');
    tmpflds = fieldnames(stmp);
    for fld = tmpflds(:)'
        fld = fld{1}; %#ok<FXSET>
        if isfield(s,fld)
            warning('structmerge:overlappingFields',...
                'Overwriting field ''%s''.',fld);
        end
        s.(fld) = stmp.(fld);
    end
end


text:
function fx = p12_fun ( option, nvar, x )

%*****************************************************************************80
%
%% P12_FUN evaluates the function for problem 12.
%
%  Title:
%
%    Materially nonlinear problem.
%
%  Description:
%
%    The problem is the two point boundary value problem
%
%      U'' + LAMBDA * SIN ( U + U**2 + U**3 ) = 0
%
%    with boundary conditions
%
%      U(0) = 0.0
%      U(1) = 0.0
%
%    U is approximated by piecewise polynomials whose coefficients are
%    the unknowns U(1), ..., U(NVAR-1), and the value of LAMBDA is
%    stored as U(NVAR).
%
%  Options:
%
%    OPTION  Polynomials   Continuity
%      1     linear         1
%      2     cubic          1
%      3     cubic          2
%      4     quintic        1
%      5     quintic        2
%      6     quintic        3
%
%    All options use 8 intervals.
%
%  Comments:
%
%    The current program has zero as solution for all X(nvar).
%    Must find bifurcation branch and jump on to it.
%    Perhaps add X(nvar+1) a perturbation to right hand side.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ivo Babuska, Werner Rheinboldt,
%    Reliable Error Estimations and Mesh Adaptation for the Finite
%    Element Method,
%    in International Conference on Computational Methods
%    in Nonlinear Mechanics,
%    edited by John Oden,
%    Elsevier, 1980,
%    ISBN: 0444853820,
%    LC: QA808.I57.
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the function.
%
%    Output, real FX(NVAR-1), the value of the function at X.
%
  nbco = 1;
  nbcz = 1;
  nint = 8;
  maxpolys = 6;

  bcone(1) = 0.0;
  bczero(1) = 0.0;

  fx(1:nvar-1) = 0.0;

  if ( option == 1 )
    npolys = 2;
    nderiv = 1;
  elseif ( option == 2 )
    npolys = 4;
    nderiv = 1;
  elseif ( option == 3 )
    npolys = 4;
    nderiv = 2;
  elseif ( option == 4 )
    npolys = 6;
    nderiv = 1;
  elseif ( option == 5 )
    npolys = 6;
    nderiv = 2;
  elseif ( option == 6 )
    npolys = 6;;
    nderiv = 3;
  end

  nvary = nint * npolys;
%
%  Get the Gauss quadrature rule.
%
  [ gcoef, gpoint ] = p12_gauss8 ( );
%
%  Set up the terms A * Y involving the bivariate form
%
%  For each interval I:
%
  for i = 1 : nint

    iskip = ( i - 1 ) * npolys;
    xl = ( i - 1 ) / nint;
    xr = i / nint;
    dtdx = 2.0 / ( xr - xl );
%
%  For each Gauss point, J, evaluate the integrand.
%
    for j = 1 : 8

      t = gpoint(j);
      coef = gcoef(j) * ( xr - xl ) / 2.0;
      [ pl, pld ] = p12_legendre_val ( t, dtdx, npolys );

      u = 0.0;
      uprym = 0.0;
      for k = 1 : npolys
        u = u + x(iskip+k) * pl(k);
        uprym = uprym + x(iskip+k) * pld(k);
      end

      phi = - uprym;
      psi = x(nvar) * sin ( u * ( 1.0 + u * ( 1.0 + u ) ) );
      lskip = iskip;
%
%  Project onto each test function L.
%
      for l = 1 : npolys
        ieqn = lskip + l;
        fx(ieqn) = fx(ieqn) + coef * ( psi * pl(l) + phi * pld(l) );
      end

      lskip = lskip + npolys;

    end

  end
%
%  2. Add the terms B * Z for the continuity of the test functions.
%
%  For each interval I:
%
  for i = 1 : nint

    if ( i == 1 )
      ncl = nvary;
    else
      ncl = nvary + nbcz + ( i - 2 ) * nderiv;
    end

    ncr = nvary + nbcz + ( i - 1 ) * nderiv;
    xl = ( i - 1 ) / nint;
    xr = i / nint;
    dtdx = 2.0 / ( xr - xl );
%
%  Count conditions at the left endpoint, LHIL, and at right, LHIR.
%  If we are in the first or last interval, one of
%  these will be boundary conditions.
%
    if ( i == 1 )
      lhil = nbcz;
    else
      lhil = nderiv;
    end

    if ( i == nint )
      lhir = nbco;
    else
      lhir = nderiv;
    end
%
%  For each test function PL(K):
%
    for k = 1 : npolys

      s = r8_mop ( k + 1 );
      ieqn = ( i - 1 ) * npolys + k;
%
%  Apply the boundary conditions.
%
      h2i = 1.0;
      for l = 1 : lhil
        s = - s;
        ivar = ncl + l;
        fx(ieqn) = fx(ieqn) + s * x(ivar) * h2i * p12_theta ( l, k );
        h2i = h2i * dtdx;
      end

      h2i = 1.0;
      for l = 1 : lhir
        ivar = ncr + l;
        fx(ieqn) = fx(ieqn) + x(ivar) * h2i * p12_theta ( l, k );
        h2i = h2i * dtdx;
      end

    end

  end
%
%  3. Create the C * Y terms for U and its derivatives.
%  One equation is generated for component and condition.
%
  npsum = 0;
  dtdxr = 0.0;
  dtdxl = 0.0;
%
%  For each node:
%
  ndsum = nvary;

  for i = 1 : nint + 1

    if ( 1 < i )
      xl = ( i - 2 ) / nint;
    end

    xc = ( i - 1 ) / nint;

    if ( i < nint + 1 )
      xr = i / ( nint );
    end

    if ( xc ~= xl )
      dtdxl = 2.0 / ( xc - xl );
    end

    if ( xr ~= xc )
      dtdxr = 2.0 / ( xr - xc );
    end

    h2il = 1.0;
    h2ir = 1.0;
%
%  Count the conditions:
%
    if ( i == 1 )
      khi = nbcz;
    elseif ( i < nint + 1 )
      khi = nderiv;
    elseif ( i == nint + 1 )
      khi = nbco;
    end

    for k = 1 : khi

      s = r8_mop ( k + 1 );
%
%  Set up the term from the left hand interval.
%
      ieqn = ndsum + k;

      if ( i == 1 )

        fx(ieqn) = fx(ieqn) + bczero(k);

      else

        for l = 1 : npolys
          ivar = npsum + l - npolys;
          fx(ieqn) = fx(ieqn) + x(ivar) * h2il * p12_theta ( k, l );
        end

      end
%
%  Set up the term from the right hand interval.
%
      if ( i == nint + 1 )

        fx(ieqn) = fx(ieqn) - bcone(k);

      else

        for l = 1 : npolys
          ivar = npsum + l;
          s = - s;
          fx(ieqn) = fx(ieqn) + s * x(ivar) * h2ir * p12_theta(k,l);
        end
      end

      h2il = h2il * dtdxl;
      h2ir = h2ir * dtdxr;

    end

    ndsum = ndsum + khi;
    npsum = npsum + npolys;

  end

  return
end


text:
function cleanmacpollution(loc)
%CLEANMACPOLLUTION deletes the annoying search trace files generated by MAC.
%   CLEANMACPOLLUTION(WHERE) deletes the annoying search trace generated
%   by MAC OS for all subdirectorys under WHERE.
%
%   MAC OS has the bad habit of creating garbage finder data file in every
%   directory it touches. Usually this means your flashdrive will be
%   polluted by numerous unwanted files once you've plugged it into a mac.
%   On OS 9 these garbages are FINDER.DAT; on OS X they are .DS_Store;
%   Use this little tool to remove those unwanted files;

% Written by Siyi Deng; 12-19-2007;

if nargin < 1 || isempty(loc), loc = uigetdir; end
if ~loc, return; end
s = genpath(loc);
p = strfind(s,pathsep);
q = [[1,p(1:end-1)+1];p-1];
mc = {'FINDER.DAT','.DS_Store'}; % OS 9, OS X;
for v = 1:numel(mc)
    for k = 1:size(q,2)
        thePath = [s(q(1,k):q(2,k)),filesep];
        if exist([thePath,mc{v}],'file')
            disp(['Deleting ',thePath,mc{v},'...']);
            delete([thePath,mc{v}]);
        end
        if exist([thePath,'private'],'dir')
            if exist([thePath,'private',filesep,mc{v}],'file')
                disp(['Deleting ',thePath,'private',filesep,mc{v},'...']);
                delete([thePath,'private',filesep,mc{v}]);
            end
        end
    end
end

end % CLEANMACPOLLUTION;


text:
%  Function gets activations of network (net) when it infers on an image
%  (img), for only specified layers (layers_to_use)
%
%  @authors: ufukefe, kutalmisince
%  Created on March 23, 2021
%  @Middle East Technical University, Center for Image Analysis
%  Last Edited on July 1, 2021

function activations = GetActivations(img,net,layers_to_use)

    % Check the type of the input image
    if ~isa(img, 'single')
        img = single(img);
    end

    % Forward pass
    net.eval({'x0',gpuArray(img)});

    % Insert activations inside cell array
    activations = cell(size(layers_to_use,1),2);

    for i=1:size(layers_to_use,1)
        activations{i,1} = layers_to_use{i};

        activations{i,2} = ...
        gather(net.vars(getVarIndex(net,net.layers(getLayerIndex(net, layers_to_use{i})).outputs)).value);
    end

    % Reset network for memory issues
    reset(net);

end

author: ufukefe

repo: DFM

save_path: github-repos/MATLAB/ufukefe-DFM


text:
%Copyright 2008-2009  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
%Copyright 2008-2009  David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
%
%THE BSD LICENSE
%
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions
%are met:
%
%1. Redistributions of source code must retain the above copyright
%   notice, this list of conditions and the following disclaimer.
%2. Redistributions in binary form must reproduce the above copyright
%   notice, this list of conditions and the following disclaimer in the
%   documentation and/or other materials provided with the distribution.
%
%THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
%IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
%OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
%IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
%INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
%NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
%DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
%THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
%(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
%THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function index = flann_load_index(filename, dataset)
%FLANN_LOAD_INDEX  Loads an index from disk
%
% Marius Muja, March 2009

    index = nearest_neighbors('load_index', filename, dataset);
end

repo: gspbox


text:
% Mathematics Q43918
% https://dsp.stackexchange.com/questions/questions/43918
% Determine the Unit Impulse Response of ARMA Filter for 8 First Indices
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     25/09/2017
%   *   First release.

%% General Parameters

run('InitScript.m');

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

generateFigures = OFF;

%% Simulation Parameters

numSamples = 8;

% Filter Coefficients
% See MATLAB's Filter Documentation for the sign.
vA = [1; -0.5; -(1 / 8)];
vB = [0; 0; 0.5];

%% Generate Data

% Unit Impulse
vX = zeros([numSamples, 1]);
vX(1) = 1;

%% Simulate Result

vY = filter(vB, vA, vX);

%% Display Results

hFigure         = figure('Position', figPosDefault);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries  = line([0:(numSamples - 1)], [vX, vY]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['ARMA Filter Unit Impulse Response']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Sample Indices']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Value']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Unite Impulse'], ['ARMA Filter Response']});

%% Restore Defaults

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

author: RoyiAvital

repo: StackExchangeCodes


text:
function julian = jed_to_yearcount_julian ( jed )

%*****************************************************************************80
%
%% JED_TO_YEARCOUNT_JULIAN converts a JED to a Julian year count.
%
%  Discussion:
%
%    An average year in the Julian calendar is exactly 365.25 days long.
%    This calculation counts the number of average Julian years from
%    the beginning of the year 2000.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, real JULIAN, the Julian year.
%
  year_length = 365.25;

  jed_epoch = epoch_to_jed_y2k ( );

  julian = 2000.0 + ( jed - jed_epoch ) / year_length;

  return
end


text:
function [ grid_num, point_num ] = gl_grids_minmax_size ( dim_num, q_min, ...
  q_max )

%*****************************************************************************80
%
%% GL_GRIDS_MINMAX_SIZE returns sizes for the GL_GRIDS_MINMAX.
%
%  Discussion:
%
%    This routine can be used to determine the necessary size to be
%    allocated to arrays GRID_ORDER and GRID_POINT in a call to
%    GL_GRIDS_MINMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer Q_MIN, Q_MAX, the minimum and maximum values of
%    Q, the sum of the orders in each spatial coordinate.
%
%    Output, integer GRID_NUM, the number of Gauss-Legendre
%    grids whose Q value is between Q_MIN and Q_MAX.
%
%    Output, integer POINT_NUM, the total number of points in the grids.
%

%
%  Determine the total number of points that will be generated
%  by "going through the motions".
%
  point_num = 0;
  grid_num = 0;

  for q = q_min : q_max

    more = 0;
    order_1d = [];

    while ( 1 )

      [ order_1d, more ] = compnz_next ( q, dim_num, order_1d, more );

      order_nd = prod ( order_1d(1:dim_num) );

      point_num = point_num + order_nd;

      grid_num = grid_num + 1;

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function out = feval(N, varargin)
%FEVAL   Evaluate the operator of the CHEBOP at a CHEBFUN or CHEBMATRIX.
%   OUT = FEVAL(N, U) for a CHEBFUN or CHEBMATRIX U applies the CHEBOP N to U,
%   i.e., it returns N(U). Here, N.OP should be of the form @(u) diff(u,2) + ...
%   If N.op is of the form @(x, u) diff(u,2) + ... then an x variable is
%   instantiated internally and included automatically, however, this is not
%   the preferred syntax and may not be supported in future releases.
%
%   OUT = FEVAL(N, X, U) for the CHEBFUN X and CHEBFUN or CHEBMATRIX U applies
%   the CHEBOP N to X and  U, i.e., it returns N(X, U) where N.OP has the form
%   @(x, u) diff(u,2) + .... Here, X should be the dependent variable on
%   N.DOMAIN.
%
%   OUT = FEVAL(N, X, U1, U2, ..., UM) for a CHEBFUN X and CHEBFUN or CHEMBATRIX
%   objects U1, ..., UM applies the CHEBOP N to the functions Uk; i.e., it
%   returns N(X, U1, U2, ..., UM) where N.OP has the form @(x, u1, u2, ..., um).
%   Note that for systems of equations, X _must_ be included in N.OP.
%
%   OUT = FEVAL(N, X, U) where U is a CHEBMATRIX of M entries and N.OP has the
%   form @(X, U1, U2, ..., UM) is equivalent FEVAL(N, X, U{1}, ..., U{M}).
%   Again, OUT = FEVAL(N, U) will also work in this situation, but this is not
%   the preferred syntax.
%
%   OUT = FEVAL(N, DIM) returns an DIM-point discretization of the linear
%   operator N. If N is not linear an error is thrown. OUT = FEVAL(N, DIM,
%   'oldschool') uses boundary bordering to deal wiuth the boundary conditions,
%   rather than the rectangular projection approach of hale and Driscoll. Note
%   that this syntax exists only the support ATAP and doesn't necessarily give a
%   clear picture of the discretizations now being used in the Chebfun release.
%   CHEBOP/MATRIX, which is the preferred syntax for this functionality,
%   provides further details
%
% See also CHEBOP/SUBSREF, LINOP/MTIMES, CHEBOP/MATRIX.

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

% Support for calling a linear CHEBOP with a numerical input to get its
% discretization. This is deprecated
if ( (nargin == 2) && isnumeric(varargin{1}) )
    warning('CHEBFUN:CHEBOP:feval:deprecated', ...
        ['FEVAL(N, DIM) or N(DIM) exists only to provide backwards\n', ...
        'compatibility with ATAP. The preferred method for visualizing a\n', ...
        'discretization of a linear CHEBOP is MATRIX(N, DIM). Note, however,\n', ...
        'that these may not give the same result due to changes in how\n', ...
        'CHEBOP discretizes differential operators.'])
    warning('off', 'CHEBFUN:CHEBOP:feval:deprecated');

    % We cannot support boundary conditions in this way.
    throwBCwarning = false;
    if ( ~isempty(N.lbc) )
        N.lbc = [];
        throwBCwarning = true;
    end
    if ( ~isempty(N.rbc) )
        N.rbc = [];
        throwBCwarning = true;
    end
    if ( ~isempty(N.bc) )
        N.bc = [];
        throwBCwarning = true;
    end
    if ( throwBCwarning )
        warning('CHEBFUN:CHEBOP:feval:BCs', ...
            'Boundary conditions are not supported in FEVAL(N, DIM).')
    end

    % Because the native V5 matrix is rectangular, we have to do some resizing
    % to recapture V4 behavior. Note: this may break for systems, piecewise
    % cases.
    n = varargin{1};
    L = linop(N);
    pref = cheboppref();
    disc = pref.discretization;
    if ( isa(disc, 'char') && strcmpi(disc, 'values') )
        pref.discretization = @chebcolloc2;
    elseif ( isa(disc, 'function_handle') && ~isequal(disc, @chebcolloc2) )
        % We only support CHEBCOLLOC2 discretizations!
        error('CHEBFUN:CHEBOP:feval:notColloc2', ...
            ['FEVAL(N, DIM) only supports CHEBCOLLOC2 discretizations.\n', ...
             'Use MATRIX(N, DIM) or change the discretization in CHEBOPPREF.']);
    end
    A = matrix(L, n, pref); % has n rows but maybe more columns

    % We want an n by n result. We have that A is (n-d) x n where d =
    % L.diffOrder. Also, the output of A is at 1st kind points. We need to do:
    %  (map from n 1st kind to n 2nd kind) * A * (map from n 2nd kind to
    %  n+d 2nd kind).

    % Note that we don't have to translate/scale to the domain for barymat
    % matrices that go between grids.
    d = L.diffOrder;
    [x1,~,w1] = chebtech1.chebpts( n );
    x2 = chebtech2.chebpts( n );
    x3 = chebtech2.chebpts( n + d );

    out = barymat(x2, x1, w1) * A * barymat(x3, x2);
    return

elseif ((nargin == 3) && ischar(varargin{2}) ...
        && strcmpi(varargin{2}, 'oldschool'))

    % TODO: Is this still correct?
    out = matrix(N, varargin{:});
    return
end

% We must expand CHEBMATRIX entries out to a cell for {:} to work below.
isChebMat = cellfun(@(u) isa(u, 'chebmatrix'), varargin);
if ( any(isChebMat) )
    args = {};
    for k = 1:numel(varargin)
        % Append variables from the kth input:
        if ( isChebMat(k) )
            args = [args , varargin{k}.blocks.']; %#ok<AGROW>
        else
            args = [args , varargin(k)];          %#ok<AGROW>
        end
    end
else
    args = varargin;
end

% How many input arguments are there to N.op?
numberOfInputs = nargin(N);

% If no arguments, return empty:
if ( numberOfInputs == 0 )
    out = [];
    return
end

if ( numberOfInputs == 1)
    % If N has one input arguments, either we have a scalar problem, or the
    % problem is specified on chebmatrix format, e.g.,
    %   N.op = @(u) [ diff(u{1},2) + u{2}; u{1} + diff(u{2}];
    % Here, importantly, x does not appear in the argument list for N.op.
    u = varargin{1};

    % If we have a scalar problem, but U is still a CHEBMATRIX, we need to
    % extract the BLOCK of U in order to be evaluate N. If on the other hand, U
    % has more than one block, but NUMBEROFINPUTS is still equal to 1 (which got
    % us here in the first place), we must be dealing with a CHEBOP N, whose OP
    % is specified on CHEBMATRIX format, e.g.
    %   N.op = @(u) [diff(u{1}) + u{2}; u{1} + diff(u{2})];
    if ( isa(u, 'chebmatrix') && max(size(u)) == 1 )
        u = u.blocks{1};
    end

    out = N.op(u);

elseif ( numberOfInputs == 2 )
    % If N has two input arguments, either we have a scalar problem, or the
    % problem is specified on chebmatrix format, e.g.,
    %   N.op = @(x, u) [ diff(u{1},2) + u{2}; u{1} + diff(u{2}];
    % Here, importantly, x must appear in the argument list for N.op.

    % Did we not get the x variable passed in as argument?
    if ( numel(varargin) == 1 )
        u = varargin{1};

        % Construct the independent variable X.
        x = chebfun(@(x) x, N.domain);
        % If the CHEBFUN U passed in is a quasimatrix, we need to tile the
        % independent variable X to have the same dimensions as U:
        x = repmat(x, 1, size(u,2));

    elseif ( numel(varargin) == numberOfInputs )
        % Got passed both X and U.
        x = varargin{1};
        u = varargin{2};

    else
        error('CHEBFUN:CHEBOP:feval:numInputs', ...
            'Unexpected number of input arguments.')
    end

    % If we have a scalar problem, but U is still a CHEBMATRIX, we need to
    % extract the BLOCK of U in order to be evaluate N. If on the other
    % hand, U has more than one block, but NUMBEROFINPUTS is still less than
    % or equal to 2 (which got us here in the first place), we must be
    % dealing with a CHEBOP N, whose OP is specified on CHEBMATRIX format,
    % e.g.,
    %   N.op = @(x,u) [diff(u{1}) + u{2}; u{1} + diff(u{2})];
    if ( isa(u, 'chebmatrix') && max(size(u)) == 1 )
        u = u.blocks{1};
    end

    % Evaluate the operator!
    out = N.op(x, u);

else
    % The operator is specified on the form
    %   N.op = @(x, u, v) = [diff(u,2) + v; u + diff(v)]

    % Count the number of RHSs:
    numCols = max(max(cellfun(@(v) size(v, 2), varargin)));

    % We must expand CHEBMATRIX entries out to a cell for {:} to work below.
    isChebMat = cellfun(@(u) isa(u, 'chebmatrix'), varargin);
    if ( any(isChebMat) )
        args = {};
        for k = 1:numel(varargin)
            % Append variables from the kth input:
            if ( isChebMat(k) )
                args = [args , varargin{k}.blocks.']; %#ok<AGROW>
            else
                args = [args , varargin(k)];          %#ok<AGROW>
            end
        end
        % ARGS need to be vertically concatenated for operator to be evaluated
        % correctly.
        args = args.';
    else
        % ARGS need to be vertically concatenated for operator to be evaluated
        % correctly.
        args = varargin.';
    end

    if ( numCols > 1 )
        % Deal with multiple RHS
        out = cell(1, numCols);
        for k = 1:numCols
            out{k} = doEval(N, args(:,k));
        end
        out = horzcat(out{:});
    else
        out = doEval(N, args);
    end

end

% If the operator is written in a function file or nested function, then it
% is natural for it to be returned in cell form. Convert it to a chebmatrix:
if ( iscell(out) )
    out = vertcat(out{:}); % TODO: is this correct for multiple RHSs?
end

end

function out = doEval(N, args)

    numberOfInputs = nargin(N);
    if ( numel(args) == numberOfInputs - 1 )
        % Check if we need to include an x (independent variable):
        x = chebfun(@(x) x, N.domain);
        args = [ {x} ; args ];
    end
    % Evaluate the operator:
    out = N.op(args{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function found = existDataType(dataTypeName,dataTypes,fullMatch)
%
%      found = existDataType([dataTypeName],[dataTypes],[fullMatch])
%
%   Check whether a dataType with dataTypeName exists in the mrSESSION structure.
%
%   If a full match to the dataType name is not found, then we return
%   found=0.
%
%   If a full match is found, mrSESSION.dataTypes(found) is the dataType.
%
%   If a partial match (fullMatch = 0) is requested, then existDataType() can
%   find several matches (i.e., whenever the first part of the data type name matches
%   the dataTypeName string). In this case, found returns a list of indices
%   whose names match the dataTypeName.  If no matches are found, then
%   found = 0 is returned.
%
%   Examples:
%       existDataType('Atlases',dataTypes,0) returns a list of partial matches
%        mrSESSION.dataTypes(found) all begin with the phrase 'Atlases', such as
%        'Atlases-1'.
%
% dataTypes: optional dataTypes structure. Default is to use the global dataTYPES.
%
% djh, 3/2001

global dataTYPES

if ~exist('dataTypes','var') || isempty(dataTypes),  dataTypes = dataTYPES; end
if ~exist('fullMatch','var'), fullMatch = 1; end

ndataTypes = length(dataTypes);

if fullMatch
    found = 0;
    for ii=1:ndataTypes
        if strcmp(dtGet(dataTypes(ii), 'Name'),dataTypeName);
            found = ii;
            return;
        end
    end
else
    found = zeros(1,ndataTypes);
    for ii=1:ndataTypes
        found(ii) =  strncmp(dtGet(dataTypes(ii), 'Name'),dataTypeName,length(dataTypeName));
    end
    found = find(found);
    if isempty(found), found = 0; end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function military = ch_to_military ( ch )

%*****************************************************************************80
%
%% CH_TO_MILITARY converts an ASCII character to a Military code word.
%
%  Example:
%
%    'A'  'Alpha'
%    'B'  'Bravo'
%    'Z'  'Zulu'
%    'a'  'alpha'
%    '7'  '7'
%    '%'  '%'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CH, the ASCII character.
%
%    Output, string MILITARY, the military code word.
%    If CH is not an alphabetic letter, then MILITARY is simply set equal to CH.
%
    code = [ ...
    'alpha   '; ...
    'bravo   '; ...
    'charlie '; ...
    'delta   '; ...
    'echo    '; ...
    'foxtrot '; ...
    'golf    '; ...
    'hotel   '; ...
    'india   '; ...
    'juliet  '; ...
    'kilo    '; ...
    'lima    '; ...
    'mike    '; ...
    'november'; ...
    'oscar   '; ...
    'papa    '; ...
    'quebec  '; ...
    'romeo   '; ...
    'sierra  '; ...
    'tango   '; ...
    'uniform '; ...
    'victor  '; ...
    'whiskey '; ...
    'x-ray   '; ...
    'yankee  '; ...
    'zulu    ' ];

  if ( 'A' <= ch && ch <= 'Z' )
    military = code( ch - 'A' + 1, : );
    military(1) = ch_cap ( military(1) );
  elseif ( 'a' <= ch && ch <= 'z' )
    military = code( ch - 'a' + 1, : );
  else
    military = ch;
  end

  return
end


text:
function [ a, info ] = cpofa ( a, lda, n )

%*****************************************************************************80
%
%% CPOFA factors a complex hermitian positive definite matrix.
%
%  Discussion:
%
%    CPOFA is usually called by CPOCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%    (time for CPOCO) = (1 + 18/N) * (time for CPOFA).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 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 A(LDA,N), the hermitian matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); an upper triangular matrix R so that
%      A = hermitian(R)*R
%    where hermitian(R) is the conjugate transpose.  The strict lower
%    triangle is unaltered.  If INFO /= 0, the factorization is not
%    complete.  Only the diagonal and upper triangle are used.
%
%    Output, integer INFO.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is
%    not positive definite.
%
  info = 0;

  for j = 1 : n

    s = 0.0;

    for k = 1 : j-1
      t = a(k,j) - a(1:k-1,k)' * a(1:k-1,j);
      t = t / a(k,k);
      a(k,j) = t;
      s = s + real ( t * conj ( t ) );
    end

    s = real ( a(j,j) ) - s;

    if ( s <= 0.0 | imag ( a(j,j) ) ~= 0.0 )
      info = j;
      break
    end

    a(j,j) = sqrt ( s );

  end

  return
end


text:
%% Copyright (C) 2014-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym log (@var{x})
%% Symbolic log function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = log (x)
%%   @result{} y = (sym) log(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = log(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('log', x);
end

%!error log (sym(1), 2)
%!assert (isequaln (log (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = log(x);
%! f2 = log(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = log(A);
%! f2 = log(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = log (d);
%! f = log (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function [rRelmTest] = relm_NTest(vRatesH, vRatesN, nNumberSimulation, fMagThreshold, bOptimized, bDrawFigure)
% function [rRelmTest] = relm_NTest(vRatesH, vRatesN, nNumberSimulation, fMagThreshold, bOptimized, bDrawFigure)
% --------------------------------------------------------------------------------------------------------------
% Computation of the N-test for the RELM framework
%
% Input parameters:
%   vRatesH                       Matrix with rates of the test hypothesis
%   vRatesN                       Matrix with rates of the null hypothesis
%   nNumberSimulation             Number of random simulations
%   fMagThreshold                 Magnitude threshold (Use only bins with magnitude >= threshold
%   bOptimized                    0 (default): use a for loop, 1: matrix-wise calculation (needs a lot of memory)
%   bDrawFigure                   Draw the cumulative density plot after testing (default: off)
%
% Output paramters:
%   rRelmTest.fAlpha              Alpha-value of the cumulative density
%   rRelmTest.fBeta               Beta-value of the cumulative density
%   rRelmTest.vSimValues_H        Vector containing the sorted simulated numbers of events for the test hypothesis
%   rRelmTest.vSimValues_N        Vector containing the sorted simulated numbers of events for the null hypothesis
%   rRelmTest.nNumberSimulation   Number of random simulations
%   rRelmTest.fObservedData       Observed total number of events
%
% Danijel Schorlemmer
% October 4, 2002

% Exit on empty rate matrices
if isempty(vRatesH)  ||  isempty(vRatesN)
  rRelmTest.fAlpha = nan;
  rRelmTest.fBeta = nan;
  rRelmTest.vSimValues_H = nan;
  rRelmTest.vSimValues_N = nan;
  rRelmTest.nNumberSimulation = nan;
  rRelmTest.fObservedData = nan;
  return;
end

if ~exist('bDrawFigure')
  bDrawFigure = 0;
end

if ~exist('bOptimized')
  bOptimized = 0;
end

% Randomize
rand('state',sum(100*clock));

% Get the necessary data from the rate matrices and weight them properly
[vLambdaH, vLambdaN, vNumberQuake] = relm_PrepareData(vRatesH, vRatesN, fMagThreshold);
nNumberQuake = sum(vNumberQuake);

% Get the number of bins (rows)
[nRow, nColumn] = size(vLambdaH);

if bOptimized
  % Create the random numbers for the simulation
  vRandom = rand(nRow, nNumberSimulation);

  % Replicate the rate vectors
  vLambdaH = repmat(vLambdaH, 1, nNumberSimulation);
  vLambdaN = repmat(vLambdaN, 1, nNumberSimulation);

  % Compute the simulated number of events and sum them up
  vNum_H = poissinv(vRandom, vLambdaH);
  vNum_N = poissinv(vRandom, vLambdaN);
  vSimNum_H = sum(vNum_H);
  vSimNum_N = sum(vNum_N);
else
  % Create empty vectors for the total number of events
  vSimNum_H = [];
  vSimNum_N = [];

  % Loop over the simulations
  for nCnt = 1:nNumberSimulation
    % Create the random numbers for the simulation
    vRandom = rand(nRow, 1);

    % Compute the simulated number of events and sum them up
    vNum_H = poissinv(vRandom, vLambdaH);
    vNum_N = poissinv(vRandom, vLambdaN);
    vSimNum_H = [vSimNum_H; nansum(vNum_H)];
    vSimNum_N = [vSimNum_N; nansum(vNum_N)];
  end
end

% Sort them for the cumulative density plot
rRelmTest.vSimValues_H = sort(vSimNum_H);
rRelmTest.vSimValues_N = sort(vSimNum_N);

% Compute Alpha and Beta and store the important parameters
rRelmTest.fAlpha = sum(rRelmTest.vSimValues_N > nNumberQuake)/nNumberSimulation;
rRelmTest.fBeta = sum(rRelmTest.vSimValues_H < nNumberQuake)/nNumberSimulation;
rRelmTest.nNumberSimulation = nNumberSimulation;
rRelmTest.fObservedData = nNumberQuake;

if bDrawFigure
  relm_PaintCumPlot(rRelmTest, 'Number of earthquakes');
end

author: CelsoReyes


text:
function jmp = estimateJump(series,jmp_idx,same_slope,mode)
    %Description: estimate jmp in a series doing a dtrend using medians
    %INPUT:
    % series : the series
    % jmp_idx : where the jmp occurs
    % mode : 1 joint detrend , 2: detrend separately
    if nargin < 3
        same_slope = true;
        mode = 'median';
    end
    if nargin < 4
        mode = 'median';
    end
    if strcmp(mode, 'median')
        diff_pre = diff(series(1:jmp_idx-1));
        diff_post = diff(series(jmp_idx:end));
        if same_slope
            trnd = median([diff_pre ;diff_post],'omitnan');
            series = series - trnd * [1:length(series)]';
        else
            trend1 = median([diff_pre],'omitnan');
            trend2 = median([diff_post],'omitnan');
            %series(1:jmp_idx-1) = series(1:jmp_idx-1) - trend1 * [1 : jmp_idx-1]';
            %series(jmp_idx:end) = series(jmp_idx:end) - trend2 * [1 : (length(series) - jmp_idx + 1)]';
        end
        jmp = median(series(jmp_idx:end),'omitnan') - median(series(1:jmp_idx-1),'omitnan')  ;
    elseif strcmp(mode, 'ls')
        valid_idx_l = ~isnan(series);
        valid_idx = find(valid_idx_l);
        idx_bf = valid_idx < jmp_idx;
        idx_aft = valid_idx >= jmp_idx;
        n_valid = sum(valid_idx_l);
        if same_slope
            A= zeros(n_valid,3);
            A(:,1) = find(valid_idx_l);
            A(idx_bf,2) = 1;
            A(idx_aft,3) = 1;
            x = A\series(valid_idx_l);
            jmp = x(3) -x(2);
        else
            A= zeros(n_valid,4);
            A(idx_bf, 1) = valid_idx(idx_bf);
            A(idx_aft, 2) = valid_idx(idx_aft);
            A(idx_bf, 3) = 1;
            A(idx_aft, 4) = 1;
            x = A\series(valid_idx_l);
            jmp = x(4) -x(3);
        end
    end
end

author: goGPS-Project


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        Sample Data and Metch Registration Sessions        %
%                                                           %
%           by Qianqian Fang <q.fang at neu.edu>            %
%                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first of all, make sure you've already add the metch root
% folder to your matlab/octave search path list

% load the sample data, where
%   no: the node coordinates of a surface mesh
%   el: the surface triangles
%   pt: the point cloud to be registered

load sampledata

% start the metch-gui to perform the registration

%== Description of the workflow ==
%
% 1. when the GUI pops up, it will display the mesh and the points,
%    you can rotate both plots so that you can identify the matching
%    features
% 2. switch on "Select" mode, then, click on a land-mark point on the point
%    plot, when a data-tip shows up, click "Add Selected" button
% 3. click on the corresponding position on the mesh, and click
%    "Add Selected"
% 4. repeat the above for at least 4 point pairs (you can select more);
%    if you want to change views, switch off "Select" box and rotate;
%    after rotation, switch on "Select" box again
% 5. click "Initialize": this will create the initial mapping using the
%    selected point pairs
% 6. click "Optimize": this will fit the surface with the whole point cloud
% 7. click "Proj2Mesh": this will project the fitted point clouds onto the
%    mesh
% 8. you can quit the GUI by hit "Close", your results will be saved to reg
% 9. close the window

if(exist('OCTAVE_VERSION')~=0)
        reg=metchgui_one(no,el,pt);
else
        reg=metchgui(no,el,pt);
end

author: fangq


text:
classdef KosDragCoeffientModel < AbstractDragCoefficientModel
    %kOSDragCoefficientModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        dataFile(1,:) char

        machNum(:,1) double
        aoa(:,1) double
        sideslip(:,1) double
        data double

        giDragCube
        giOtherDrag
    end

    properties(Access=private)
        giRnCorr
    end

    properties(Constant)
        enum = DragCoefficientModelEnum.kOSModel
    end

    methods
        function obj = KosDragCoeffientModel(dataFile)
            obj.dataFile = dataFile;

            obj.dataFile = dataFile;

            if(isfile(dataFile))
                obj.createGriddedInterpFromFile();
            else
                obj.machNum = [0; 1];
                obj.aoa = [0; deg2rad(1)];
                obj.sideslip = [0; deg2rad(1)];

                rows = combvec(obj.machNum', obj.aoa', obj.sideslip')';
                rows(:,4:5) = 0;
                obj.data = rows;

                obj.createGriddedInterpFromData();
            end

            %see Profile.ks in the Project-Atmospheric-Drag for details of
            %this correction factor
            %https://github.com/Ren0k/Project-Atmospheric-Drag/blob/main/CCAT/DragProfile/LIB/Profile.ks
            rnCorrData = [0.00, 4.0000;
                            0.0001 ,3.0000;
                            0.01   ,2.0000;
                            0.10   ,1.2000;
                            1.00   ,1.0000;
                            100.0  ,1.0000;
                            200.0  ,0.8200;
                            500.0  ,0.8600;
                            1000.0 ,0.9000;
                            10000  ,0.9500];
            obj.giRnCorr = griddedInterpolant(rnCorrData(:,1), rnCorrData(:,2), 'pchip', 'nearest');
        end

        function createGriddedInterpFromFile(obj)
            obj.data = readmatrix(obj.dataFile);
            obj.data(:,2:3) = deg2rad(obj.data(:,2:3)); %convert aoa and sideslip to radian from deg

            obj.data = sortrows(obj.data, [3 2 1]);
            obj.machNum = unique(obj.data(:,1));
            obj.aoa = unique(obj.data(:,2));
            obj.sideslip = unique(obj.data(:,3));

            obj.createGriddedInterpFromData();
        end

        function createGriddedInterpFromData(obj)
            %drag cube Model
            c = obj.data(:,4);
            V = reshape(c, [length(obj.machNum), length(obj.aoa), length(obj.sideslip)]);

            gridVecs = {obj.machNum,obj.aoa,obj.sideslip};

            obj.giDragCube = griddedInterpolant(gridVecs,V, "linear", "nearest");

            %other drag Model
            c = obj.data(:,5);
            V = reshape(c, [length(obj.machNum), length(obj.aoa), length(obj.sideslip)]);

            gridVecs = {obj.machNum,obj.aoa,obj.sideslip};

            obj.giOtherDrag = griddedInterpolant(gridVecs,V, "linear", "nearest");
        end

        function clearData(obj)
            obj.machNum = [];
            obj.aoa = [];
            obj.sideslip = [];
            obj.data = [];
            obj.giDragCube = [];
        end

        function CdA = getDragCoeff(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEFMag, totalAoA, aoa, sideslip)
            arguments
                obj(1,1) KosDragCoeffientModel
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                mass(1,1) double
                altitude(1,1) double
                pressure(1,1) double
                density(1,1) double
                vVectECEFMag(1,1) double
                totalAoA(1,1) double = 0;
                aoa(1,1) double = 0;
                sideslip(1,1) double = 0;
            end

            if(density > 0)
                sideslip = angleNegPiToPi(sideslip);

                %see useProfile.ks in the Project-Atmospheric-Drag for details
                %of this algorithm
                %https://github.com/Ren0k/Project-Atmospheric-Drag/blob/main/CCAT/useProfile.ks
                reynoldsNumber = density*vVectECEFMag;
                reynoldsCorrection = obj.giRnCorr(reynoldsNumber);

                pressurePa = pressure*1000;
                speedSound = sqrt(1.4 * pressurePa / density); %m/s
                speedMS = vVectECEFMag*1000;                   %m/s
                thisMachNum = speedMS / speedSound;

                dragCubeCdA = obj.giDragCube(thisMachNum, aoa, sideslip);
                otherDragCdA = obj.giOtherDrag(thisMachNum, aoa, sideslip);

                CdA = ((dragCubeCdA*reynoldsCorrection)+otherDragCdA);
            else
                CdA = 0;
            end
        end

        function tf = usesTotalAoA(obj)
            tf = false;
        end

        function tf = usesAoaAndSideslip(obj)
            tf = true;
        end

        function useTf = openEditDialog(obj, lvdData)
            out = AppDesignerGUIOutput({false});
            lvd_EditKosDragPropertiesGUI_App(obj, lvdData, out);
            useTf = out.output{1};
        end

        function plotDragEnvelope(obj, hAx, machNum)
            arguments
                obj(1,1) KosDragCoeffientModel
                hAx = [];
                machNum(1,1) double = 0;
            end

            if(isempty(obj.machNum) || isempty(obj.aoa) || isempty(obj.sideslip))
                warning('No kOS drag coefficient data loaded: cannot plot.');
                return;
            end

            if(isempty(hAx))
                hAx = axes(figure());
            end

            dragCubeData = obj.data(:,4);

            maxCdA = max(dragCubeData);
            minCdA = min(dragCubeData);
            levels = linspace(minCdA, maxCdA, 50);

            warning("off",'MATLAB:contour:ConstantData');

            allVect = combvec(obj.aoa(:)', obj.sideslip(:)')';
            allVect = [machNum*ones(height(allVect), 1), allVect];
            cd = obj.giDragCube(allVect(:,1), allVect(:,2), allVect(:,3));

            [~,hContour] = contour(hAx, rad2deg(obj.aoa), rad2deg(obj.sideslip), reshape(cd, [length(obj.aoa), length(obj.sideslip)]), levels, 'Fill','on');
            xlabel(hAx, 'Angle of Attack [deg]');
            ylabel(hAx, 'Sideslip Angle [deg]');
            title(hAx, sprintf('Cd*A for Mach Number = %0.3f', machNum));
            grid(hAx, 'on');

            hContour.DataTipTemplate.DataTipRows(1).Label = 'AoA (deg)';
            hContour.DataTipTemplate.DataTipRows(2).Label = 'Sideslip (deg)';
            hContour.DataTipTemplate.DataTipRows(3).Label = 'Cd*A (m^2)';

            hContour.DataTipTemplate.DataTipRows(1).Format = '%0.3f';
            hContour.DataTipTemplate.DataTipRows(2).Format = '%0.3f';
            hContour.DataTipTemplate.DataTipRows(3).Format = '%0.3f';

            warning("on",'MATLAB:contour:ConstantData');

            hC = colorbar(hAx);
            hC.Label.String = 'Cd*A [m^2]';
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ labels, valid_ids, vid_ids  ] = extract_SEMAINE_labels( SEMAINE_dir, recs, aus )
%EXTRACT_SEMAINE_LABELS Summary of this function goes here
%   Detailed explanation goes here

    % Get the right eaf file

    aus_SEMAINE = [2 12 17 25 28 45];

    inds_to_use = [];

    for i=1:numel(aus)

        inds_to_use = cat(1, inds_to_use, find(aus_SEMAINE == aus(i)));

    end

    labels = cell(numel(recs), 1);
    valid_ids = cell(numel(recs), 1);
    vid_ids = zeros(numel(recs), 2);

    for i=1:numel(recs)

        file = dir([SEMAINE_dir, '/', recs{i}, '/*.eaf']);

        vid_ids(i,:) = dlmread([SEMAINE_dir, '/', recs{i}, '.txt'], ' ');

        xml_file = [SEMAINE_dir, recs{i}, '\' file.name];
        [root_xml, name_xml, ~] = fileparts(xml_file);
        m_file = [root_xml, name_xml, '.mat'];

        if(~exist(m_file, 'file'))
            activations = ParseSEMAINEAnnotations([SEMAINE_dir, recs{i}, '\' file.name]);
            save(m_file, 'activations');
        else
            load(m_file);
        end
        if(size(activations,1) < vid_ids(i,2))
            vid_ids(i,2) = size(activations,1);
            if(vid_ids(i,2) > 2999)
                vid_ids(i,1) = vid_ids(i,2) - 2999;
            end
        end

        labels{i} = activations(vid_ids(i,1)+1:vid_ids(i,2), 1 + inds_to_use);

        % all indices in SEMAINE are valid
        valid_ids{i} = ones(size(labels{i},1),1);

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function stla_write ( output_file_name, node_num, face_num, node_xyz, ...
  face_node, face_normal )

%*****************************************************************************80
%
%% STLA_WRITE writes graphics information to an ASCII StereoLithography file.
%
%  Example:
%
%    solid MYSOLID
%      facet normal 0.4 0.4 0.2
%        outerloop
%          vertex  1.0 2.1 3.2
%          vertex  2.1 3.7 4.5
%          vertex  3.1 4.5 6.7
%        end loop
%      end facet
%      ...
%      facet normal 0.2 0.2 0.4
%        outerloop
%          vertex  2.0 2.3 3.4
%          vertex  3.1 3.2 6.5
%          vertex  4.1 5.5 9.0
%        end loop
%      end facet
%    end solid MYSOLID
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string OUTPUT_FILE_NAME, the name of the output file.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the node coordinates.
%
%    Input, integer FACE_NODE(3,FACE_NUM), the nodes making faces.
%
%    Input, integer FACE_MAX, the maximum number of faces.
%
%    Input, real FACE_NORMAL(3,FACE_NUM), the normal vector at each face.
%
  text_num = 0;
%
%  Open the file.
%
  iunit = fopen ( output_file_name, 'wt' );

  if ( iunit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', output_file_name );
    error ( 'STLA_WRITE - Fatal error!' );
  end

  fprintf ( iunit, 'solid MYSOLID\n' );
  text_num = text_num + 1;

  for face = 1 : face_num

    fprintf ( iunit, '  facet normal  %14f  %14f  %14f\n', ...
      face_normal(1:3,face) );
    text_num = text_num + 1;

    fprintf ( iunit, '    outer loop\n' );
    text_num = text_num + 1;

    for vertex = 1 : 3

      node = face_node(vertex,face);

      fprintf ( iunit, '      vertex    %14f  %14f  %14f\n', ...
        node_xyz(1:3,node) );
      text_num = text_num + 1;

    end

    fprintf ( iunit, '    end loop\n' );
    text_num = text_num + 1;
    fprintf ( iunit, '  end facet\n' );
    text_num = text_num + 1;

  end

  fprintf ( iunit, 'end solid MYSOLID\n' );
  text_num = text_num + 1;

  fclose ( iunit );

  return
end


[ H, ok, score ] = HM_ransac(X1, X2, 500, 0.1);

figure; clf;
marker_size = 15;
imshow([im1,im2], 'border', 'tight'); hold on;
plot(X1(1,ok), X1(2,ok), 'b.', 'MarkerSize', marker_size);
plot(X2(1,ok)+imsize1(2), X2(2,ok), 'b.', 'MarkerSize', marker_size);
plot(X1(1,~ok), X1(2,~ok), 'r.', 'MarkerSize', marker_size);
plot(X2(1,~ok)+imsize1(2), X2(2,~ok), 'r.', 'MarkerSize', marker_size);

% h = line([X1(1,ok) ; X2(1,ok)+imsize1(2)], [X1(2,ok) ; X2(2,ok)], 'LineWidth', 2) ;

X1_ok = X1(:,ok);
X2_ok = X2(:,ok);

% bundle adjustment on the normallized matching data
options = optimoptions('lsqnonlin', 'Algorithm','levenberg-marquardt', 'Display','final',...
    'MaxFunEvals',2000, 'MaxIter',1e3, 'TolFun',1e-6, 'TolX',1e-6, 'Jacobian','off');
[paras_init,resnorm_init,residual_init,exitflag_init,output_init] = lsqnonlin(...
    @(p)residual_KR_robust(double(X1_ok), double(X2_ok), imsize1, imsize2, p, 1000), [1000 1000 0 0 0],...
    [],[],options);
[paras,resnorm,residual,exitflag,output] = lsqnonlin(...
    @(p)residual_KR_robust(double(X1_ok), double(X2_ok), imsize1, imsize2, p, 10), paras_init,...
    [],[],options);
k1 = paras(1);
k2 = paras(2);
theta = paras(3:5);

K1 = [k1, 0, imsize1(2)/2;
     0, k1, imsize1(1)/2;
     0,  0, 1];
K2 = [k2, 0, imsize2(2)/2;
      0, k2, imsize2(1)/2;
      0,  0, 1];
theta_m = [0         -theta(3) theta(2)
           theta(3)  0         -theta(1)
           -theta(2) theta(1)  0];
R = expm(theta_m);

H = K2*R/K1;
M1 = K1; M2 = K2; D1 = 0; D2 = 0;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Main problem is minimize \| A x \|_1 subject to A x = b
% We introduce z = A x - b
% x \in R^n A \in R^{m x n} b \in R^m
% m < n
% ADMM terms
% minimize f(x) + g(z) subject to A x + B z = c
% f(x) : {x \in R^n | Ax = b}
% g(z) : \| z \|_1
% minimize f(x) + \| z \|_1 subject to x - z = 0
% A : 1
% B : -1
% c : 0
%
% x update x  = proj (z - u) over the set {Ax = b}
% z update: z = shrinkage(x + u)
% u update : u = u + x - z
%
% Primal residual r = A x  + B z - c :  x - z
% Dual residual s = rho A^T B (z_new - z_old) = rho (z_old - z_new)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, history] = bp(A, b, options)
    if nargin < 2 || nargin > 3
        error('Invalid arguments');
    end
    if nargin < 3
        options = struct;
    end
    % weight for the quadratic penalty term
    rho = 1;
    if isfield(options, 'rho')
        rho = options.rho;
    end
    % weight for the relaxation
    alpha = 1;
    if isfield(options, 'alpha')
        alpha = options.alpha;
    end
    % size of the matrix
    [M, N] = size(A);
    % Main optimization variable
    if isfield(options, 'x')
        x = options.x;
    else
        x = zeros(N, 1);
    end
    % Auxiliary optimization variable z = Ax - b
    % Goal of ADMM iterations is to bring z closer to Ax - b
    if isfield(options, 'z')
        z = options.z;
    else
        z = zeros(N, 1);
    end
    % scaled Lagrangian variable
    if isfield(options, 'u')
        u = options.u;
    else
        u = zeros(N, 1);
    end
    % verbosity
    verbose = 0;
    if isfield(options, 'verbose')
        verbose = options.verbose;
    end
    % Maximum number of ADMM iterations
    max_iterations = 1000;
    if isfield(options, 'max_iterations')
        max_iterations = options.max_iterations;
    end
    % absolute tolerance
    eps_abs = 1e-4;
    if isfield(options, 'absolute_tolerance')
        eps_abs = options.absolute_tolerance;
    end
    % relative tolerance
    eps_rel = 1e-2;
    if isfield(options, 'relative_tolerance')
        eps_rel = options.relative_tolerance;
    end
    %  Compute and cache factorizations
    AAt = A*A';
    P = eye(N) - A' * (AAt\A);
    % pseudo-inverse based solution of equation Ax=b
    q = A' * (AAt \b);
    % import relevant functions
    import spx.opt.shrinkage;
    if verbose
        fprintf('%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
          'r norm', 'eps pri', 's norm', 'eps dual', 'objective');
    end
    % perform ADMM iterations
    for k=1:max_iterations
        % preserve old value of z
        z_old = z;
        % update x as a projection on to the set {Ax = b}
        % we add the projection of (z-u) to the null-space of A and add
        % it to the pseudo-inverse solution of Ax = b
        x = P*(z-u) + q;
        % relaxation
        if alpha ~= 1
            x_hat = alpha*x + (1-alpha)*z;
        else
            x_hat = x;
        end
        % update z
        z = shrinkage(x_hat + u, 1/rho);
        % update u
        u = u + (x_hat - z);
        % primal residual
        r = x - z;
        r_norm = norm(r);
        % dual residual
        s = rho * (z_old - z);
        s_norm = norm(s);
        % upper bound on primal residual
        eps_primal = sqrt(N) * eps_abs + eps_rel * max([norm(x), norm(z)]);
        eps_dual = sqrt(N) * eps_abs + eps_rel * norm(rho*u);
        if verbose
            % value of objective function
            objective_value = objective(x);
            fprintf('%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
            r_norm, eps_primal, s_norm, eps_dual, objective_value);
        end
        if r_norm < eps_primal && s_norm < eps_dual
            % We have achieved convergence
            break;
        end
    end
end

% objective value
function obj = objective(x)
    obj = norm(x,1);
end

author: indigits

repo: sparse-plex


text:
% ----------------------------------------------------------------------------
% function hfssAssignSlave(fid, Name, ObjName, iUStart, iUEnd, Units,
%                              Master, [Phi], [Theta], [ReverseV])
%
% Description :
% -------------
% Creates the necessary VB Script to assign a Slave Boundary to a given Object.
%
% Parameters :
% ------------
% fid        - file identifier of the HFSS script file.
% Name       - name of the slave boundary (appears under 'Boundaries' in
%              HFSS).
% ObjName    - name of the (sheet-like) object to which the slave boundary
%              is to be assigned.
% iUStart    - (vector) starting point of the U vector. Specify as
%              [x, y, z].
% iUEnd      - (vector) ending point of the U vector. Specify as
%              [x, y, z].
% Units      - specify as 'meter', 'in', 'cm' (defined in HFSS).
% Master     - name of the master boundary to which it is assigned.
% [Phi]      - (degrees, optional) scan angle phi (defauls to 0deg).
% [Theta]    - (degrees, optional) scan angle theta (defaults to 0deg).
% [ReverseV] - (boolean, optional) reverses vector V (defaults to
%              true if the U vector points +y, and to false elsewise)
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
% hfssAssignSlave(fid, 'Slave', 'Sheet', [-width/2, 0, 0], ...
%	                 [width/2, 0, 0], 'meter', 'Master', 0, 0, false);
%
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% 21-May-2013: *Initial release (PAG).
% 01-Sep-2020: *Fixed example (DRP).
% 08-Sep-2020: *Added phi and theta as input variables (DRP).
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Written by Pablo Alcon Garcia
% pabloalcongarcia@gmail.com / palcon@tsc.uniovi.es
% 21 May 2013
% ----------------------------------------------------------------------------
function hfssAssignSlave(fid, Name, ObjName, iUStart, iUEnd, Units, ...
    Master, Phi, Theta, ReverseV)

% arguments processor.
if (nargin < 10)
    if iUEnd(2)~=iUStart(2)
        ReverseV = true;
    else
        ReverseV = false;
    end
end
if (nargin < 9)
    Theta = 0;
end
if (nargin < 8)
    Phi = 0;
end
if (nargin < 7)
	error('Insufficient # of arguments !');
end

if ReverseV
    ReverseV = 'true';
else
    ReverseV = 'false';
end

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'Set oModule = oDesign.GetModule("BoundarySetup")\n');

% Parameters
fprintf(fid, 'oModule.AssignSlave _\n');
fprintf(fid, 'Array("NAME:%s", _\n', Name);
fprintf(fid, '\t"Objects:=", Array("%s"), _\n', ObjName);
fprintf(fid, '\tArray("NAME:CoordSysVector", "Origin:=", _\n');
fprintf(fid, '\t\tArray("%f%s", "%f%s", "%f%s"), _\n', ...
        iUStart(1), Units, iUStart(2), Units, iUStart(3), Units);
fprintf(fid, '\t\t"UPos:=", Array("%f%s", "%f%s", "%f%s") _\n', ...
        iUEnd(1), Units, iUEnd(2), Units, iUEnd(3), Units);
fprintf(fid, '\t\t), _\n');
fprintf(fid, '\t"ReverseV:=", %s, _\n', ReverseV);
fprintf(fid, '\t"Master:=", "%s", _\n', Master);
fprintf(fid, '\t"UseScanAngles:=", true, _\n');
fprintf(fid, '\t\t"Phi:=", "%fdeg", "Theta:=", "%fdeg" _\n', Phi, Theta);
fprintf(fid, '\t)\n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function geometry_test080 ()

%*****************************************************************************80
%
%% TEST080 tests POLYGON_DIAMETER_2D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 4;
  dim_num = 2;

  diameter_exact = 2.0;
  v = [ ...
    1.0, 0.0; ...
    2.0, 1.0; ...
    1.0, 2.0; ...
    0.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST080\n' );
  fprintf ( 1, '  For a polygon in 2D:\n' );
  fprintf ( 1, '  POLYGON_DIAMETER_2D computes the diameter;\n' );

  r8mat_transpose_print ( dim_num, n, v, '  The polygon vertices:' );

  diameter = polygon_diameter_2d ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Diameter ( computed ) %f\n', diameter );
  fprintf ( 1, '  Diameter ( exact )    %f\n', diameter_exact );

  return
end

author: johannesgerer