text:
function [] = knn(train_data, test_data, train_target, test_target, k)
    classification_accuracy = 0;
    for i = 1:size(test_data, 1)
        %============================================
        %         Calculating Euclidean Distance
        %============================================
        D = test_data(i, :) - train_data(: , :);
        D = D.^2;
        dist_mat = sum(D, 2);
        dist_mat = sqrt(dist_mat);
        dist = [dist_mat train_target];
        %============================================
        %  Sorting Row according to minimum distance
        %============================================
        dist = sortrows(dist, 1);
        %============================================
        %       If K value is 1 Print Results
        %============================================
        if k == 1
            k_neighbours = dist(k, :);
            predicted = k_neighbours(1, 2);
            true = test_target(i, 1);
            if true == predicted
                accuracy = 1;
                classification_accuracy = classification_accuracy + accuracy;
            else
                accuracy = 0;
            end
            fprintf('ID=%5d, predicted=%3d, true=%3d, accuracy=%4.2f \n', i, predicted, true, accuracy)
        %===============================================
        % Else K value is greater then 1 Print Results
        %===============================================
        else
            k_neighbours = dist(1:k, :);
            if size(unique(k_neighbours(:, 2))) == 1
                predicted = unique(k_neighbours(:, 2));
            elseif size(unique(k_neighbours(:, 2))) == k
                predicted = k_neighbours(1, 2);
            else
                predicted = mode(k_neighbours(:, 2))
            end
            true = test_target(i, 1);
            if true == predicted
                accuracy = 1;
                classification_accuracy = classification_accuracy + accuracy;
            else
                accuracy = 0;
            end
            fprintf('ID=%5d, predicted=%3d, true=%3d, accuracy=%4.2f \n', i, predicted, true, accuracy)
        end
    end
    fprintf('classification_accuracy=%6.4f \n', classification_accuracy/size(test_target, 1))
end

repo: Machine-Learning-Models


text:
function spp_params = spp_load_pooling_params(spp_params_def)
% spp_params = spp_load_pooling_params(spp_params_def)
%   spp_params -- parameters for spatial pyramid pooling for spatial model
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

[~, ~, ext] = fileparts(spp_params_def);
if isempty(ext)
    spp_params_def = [spp_params_def, '.m'];
end
assert(exist(spp_params_def, 'file') ~= 0);

% change folder to avoid too long path for eval()
cur_dir = pwd;
[spp_def_dir, spp_def_file] = fileparts(spp_params_def);

cd(spp_def_dir);
spp_params = eval(spp_def_file);
cd(cur_dir);

spp_params.expected_scale = @spm_expected_scale;
spp_params.response_boxes = @spm_response_boxes;

end

author: ShaoqingRen


text:
function [clustCent,data2cluster,cluster2dataCell] = mean_shift_algorithm(dataPts,bandwidth,threshold);

[dim_dataset,no_data_dataset] = size(dataPts);
no_cluster = 0;         bandSq = (bandwidth).^2;
dataset_index = 1:no_data_dataset;
threshold_convergence =  threshold;% bandwidth(1,1)*exp(-3);  % this condition is taken from a reference paper on the web
tracking_array= false(1,no_data_dataset); % tracking: if a points been seen already
no_point_4_initial = no_data_dataset; % number of points to posibaly use as initilization points
clusterV = zeros(1,no_data_dataset,'uint16'); % cluster allotment method
clustMembsCell = [];
clustCent = []; % center of clust

while no_point_4_initial
    % choosing random data set for and then converging it
    random_index = ceil( (no_point_4_initial-1e-6)*rand); % pick a random seed point
    random_data_point = dataset_index(random_index); % use this point as start of mean
    mean_cur = dataPts(:,random_data_point);  % intilize mean to this points location
    cluster_members = []; % points that will get added to this cluster
    cluster_mem = zeros(1,no_data_dataset,'uint16');
    % convergence loop
    while true
        squ_Euclidean_distance(1,:) = sum(bsxfun(@minus,mean_cur(1:2,:),dataPts(1:2,:)).^2); % dist squared from mean to all points still active
        squ_Euclidean_distance(2,:) = sum(bsxfun(@minus,mean_cur(3:5,:),dataPts(3:5,:)).^2);
        kernel_range_datapoint_index = find(squ_Euclidean_distance(1,:) < bandSq(1,1)  );% points within bandWidth
        kernel_range_datapoint_index= find(squ_Euclidean_distance(2,:) < bandSq(1,2) );
        cluster_mem(kernel_range_datapoint_index) = cluster_mem(kernel_range_datapoint_index)+1;
        mean_previous = mean_cur; % save the old mean
        mean_cur = gaussian_kernel(dataPts(:,kernel_range_datapoint_index),sqrt(squ_Euclidean_distance(1,kernel_range_datapoint_index)),sqrt(squ_Euclidean_distance(2,kernel_range_datapoint_index)),bandwidth); % compute the new mean
        cluster_members = [cluster_members kernel_range_datapoint_index]; % add any point within bandWidth to the cluster
        tracking_array(cluster_members) = true; % mark that these points have been visited
        % converging condition
        if norm(mean_cur-mean_previous) < threshold_convergence
            join_cluster = 0;
            for cno = 1:no_cluster
                dist1 = norm(mean_cur(1:2)-clustCent(1:2,cno)); % spatial
                dist2 = norm(mean_cur(3:5)-clustCent(3:5,cno)); %range
                if( dist1 < bandwidth(1,1) &&dist2<bandwidth(1,2)) % condition to join the kernel
                    join_cluster = cno;
                    break;
                end
            end

            if join_cluster > 0
                nc = numel(cluster_members);
                no = numel(clustMembsCell{join_cluster});
                nw = [nc;no]/(nc+no);
                clustMembsCell{join_cluster} = unique([clustMembsCell{join_cluster},cluster_members]);
                clustCent(:,join_cluster) = mean_cur*nw(1) + mean_previous*nw(2);
                clusterV(join_cluster,:) = clusterV(join_cluster,:) + cluster_mem;
            else
                no_cluster = no_cluster+1;
                clustCent(:,no_cluster) = mean_cur;
                clustMembsCell{no_cluster} = cluster_members;
                clusterV(no_cluster,:) = cluster_mem;
            end
            break;
        end
    end
    dataset_index = find(~tracking_array);
    no_point_4_initial = length(dataset_index);
end
[~,data2cluster] = max(clusterV,[],1);
if nargout > 2
    cluster2dataCell = cell(no_cluster,1);
    for cno = 1:no_cluster
        cluster_members = find(data2cluster == cno);
        cluster2dataCell{cno} = cluster_members;
    end
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function T = auxstructure3(elem)
%% AUXSTRUCTURE3 auxiliary structure for a 3-D triangulation.
%
%  T = AUXSTRUCTURE3(elem) constucts the indices map between elements,
%  faces, edges and nodes, and the boundary edge information. T is a struct
%  data.
%
%  T.neighor(1:NT,1:4): the indices map of neighbor information of elements,
%  where neighbor(t,i) is the global index of the element oppoiste to the
%  i-th vertex of the t-th element.
%
%  T.elem2face(1:NT,1:4): the indices map from elements to faces, where
%  elem2face(t,i) is the global index of the  edge opposite to the i-th
%  vertex of t-th elemment.
%
%  T.face(1:NF,1:3): faces, where face(k,i) is the global index of the i-th
%  vetex of the k-th face, and face(k,1)<face(k,2)<face(k,3).
%
%  T.bdFace(1:Nbd,1:3): boundary faces with positive oritentation, where
%  bdFace(k,i) is the global index of the i-th vetex of the k-th boundary
%  face. The positive oritentation means that the order of the three
%  vertices of bdFace(k,:) satisfy the right-handed system and normal
%  points the outside of the domain. Note that this requires elem is
%  positive ordered, i.e., the signed volume of each tetrahedron is
%  positive. If not, use elem = fixorder3(node,elem) to fix the order.
%
%  T.face2elem(1:NF,1:4): the indices map from faces to elements, where
%  face2elem(k,1:2) are two global indices of the elements sharing the k-th
%  face, and face2elem(k,3:4) are local indices of e to edge2elem(k,1:2).
%
%  To save space all the data type in T is uint32. When use them as a input
%  of sparse(i,j,s,m,n), please change them into double type.
%
%  See also auxstructure.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1);
totalFace = uint32([elem(:,[2 3 4]); elem(:,[1 4 3]); ...
                    elem(:,[1 2 4]); elem(:,[1 3 2])]);
matlabversion = version;
if str2double(matlabversion(end-5:end-2)) > 2012
    [face, i2, j] = unique(sort(totalFace,2),'rows','lagacy');
else
    [face, i2, j] = unique(sort(totalFace,2),'rows');
end
i1(j(4*NT:-1:1)) = 4*NT:-1:1; i1 = i1';
k1 = ceil(i1/NT); t1 = i1 - NT*(k1-1);
k2 = ceil(i2/NT); t2 = i2 - NT*(k2-1);
idx = (i1 ~= i2);
neighbor = uint32(accumarray([[t1(idx),k1(idx)];[t2,k2]],[t2(idx);t1],[NT 4]));
elem2face = uint32(reshape(j,NT,4));
face2elem = uint32([t1,t2,k1,k2]);
bdElem = t1(t1 == t2);
bdk1 = k1(t1 == t2);
% consistent ordering of faces is used.
bdFace = [elem(bdElem(bdk1==1),[2 3 4]); elem(bdElem(bdk1==2),[1 4 3]);...
          elem(bdElem(bdk1==3),[1 2 4]); elem(bdElem(bdk1==4),[1 3 2])];
T = struct('neighbor',neighbor,'elem2face',elem2face,'face',uint32(face),...
           'face2elem',face2elem,'bdElem',bdElem,'bdFace',uint32(bdFace));

repo: Aero-matlab


text:
classdef KSPTOT_AppOptions < handle
    %KSPTOT_AppOptions Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %KSPTOT Options
        bodiesinifile
        rtshostname
        timesystem char
        gravparamtype char

        plotmaxdeltav double
        porkchopptsaxes double
        porkchopnumsynperiods double
        departplotnumoptiters double
        quant2opt char

        %LVD Options
        lvdrecentmsncasefiles(1,:) char = '';
    end

    events
        GravParamTypeUpdated
        TimeSystemUpdated
    end

    methods
        function obj = KSPTOT_AppOptions()

        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function mono_print_test ( )

%*****************************************************************************80
%
%% MONO_PRINT_TEST tests MONO_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_PRINT_TEST\n' );
  fprintf ( 1, '  MONO_PRINT can print out a monomial.\n' );
  fprintf ( 1, '\n' );

  m = 1;
  f = [ 5 ];
  mono_print ( m, f, '  Monomial [5]:' );

  m = 1;
  f = [ -5 ];
  mono_print ( m, f, '  Monomial [5]:' );

  m = 4;
  f = [ 2, 1, 0, 3 ];
  mono_print ( m, f, '  Monomial [2,1,0,3]:' );

  m = 3;
  f = [ 17, -3, 199 ];
  mono_print ( m, f, '  Monomial [17,-3,199]:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%klas_bayes 'Performs Bayesian maximum likelihood classification '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros las_bayes.pane file
%
% Parameters:
% InputFile: i 'Data object ', required: 'Input data object'
% InputFile: p 'Pre-Classification', required: 'File containing input predetermined classification'
% InputFile: ppf 'Prior Probabilities File (ASCII)', optional: 'Prior probabilities files (ASCII)'
% OutputFile: o 'Classified data object', required: 'Resulting output classified data object'
% OutputFile: mlo 'Maximum Likelihood object', optional: 'output maximum likelihood object'
% OutputFile: cso 'Chi-Square object', optional: 'output Chi-Square tail-rejection areas'
%
% Example: [o, mlo, cso] = klas_bayes({i, p, ppf}, {'i','';'p','';'ppf','';'o','';'mlo','';'cso',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% las_bayes - Performs Bayesian maximum likelihood classification
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% llas_bayes
%
%  RESTRICTIONS
%
%  REFERENCES
% This program is modeled after the Land Analysis System (LAS)
% program "bayes".
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc.  All rights reserved.
%

function varargout = klas_bayes(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = klas_bayes(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'p', '__input';'ppf', '__input';'o', '__output';'mlo', '__output';'cso', '__output'};
maxval={0,0,1,0,1,1};
minval={0,0,1,0,1,1};
istoggle=[0,0,1,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','OutputFile','OutputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'las_bayes"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [rhsHx, rhsHy, rhsEz] = MaxwellCurvedRHS2D(cinfo, Hx,Hy,Ez)

% function [rhsHx, rhsHy, rhsEz] = MaxwellCurvedRHS2D(cinfo, Hx,Hy,Ez)
% Purpose  : Evaluate RHS flux in 2D Maxwell TM form

Globals2D;

[rhsHx,rhsHy,rhsEz] = MaxwellRHS2D(Hx, Hy, Ez);

% correct residuals at each curved element
Ncinfo = length(cinfo);
for n=1:Ncinfo

  % for each curved element computed L2 derivatives via cubature
  cur = cinfo(n); k1 = cur.elmt; cDx = cur.Dx; cDy = cur.Dy;

  rhsHx(:,k1) = -cDy*Ez(:,k1);
  rhsHy(:,k1) =  cDx*Ez(:,k1);
  rhsEz(:,k1) =  cDx*Hy(:,k1) - cDy*Hx(:,k1);

  % for each face of each curved element use Gauss quadrature based lifts
  for f1=1:Nfaces
    k2 = EToE(k1,f1);
    gnx = cur.gnx(:,f1); gny = cur.gny(:,f1);
    gVM = cur.gVM(:,:,f1); gVP = cur.gVP(:,:,f1);
    glift = cur.glift(:,:,f1);

    % compute difference of solution traces at Gauss nodes
    gdHx = gVM*Hx(:,k1) - gVP*Hx(:,k2);
    gdHy = gVM*Hy(:,k1) - gVP*Hy(:,k2);
    gdEz = gVM*Ez(:,k1) - gVP*Ez(:,k2);

    % correct jump at Gauss nodes on domain boundary faces
    if(k1==k2)
      gdHx = 0*gdHx; gdHy = 0*gdHy; gdEz = 2*gVM*Ez(:,k1);
    end

    % perform upwinding
    gndotdH =  gnx.*gdHx+gny.*gdHy;
    fluxHx =  gny.*gdEz + gndotdH.*gnx-gdHx;
    fluxHy = -gnx.*gdEz + gndotdH.*gny-gdHy;
    fluxEz = -gnx.*gdHy + gny.*gdHx   -gdEz;

    % lift flux terms using Gauss based lift operator
    rhsHx(:,k1) = rhsHx(:,k1) + glift*fluxHx/2;
    rhsHy(:,k1) = rhsHy(:,k1) + glift*fluxHy/2;
    rhsEz(:,k1) = rhsEz(:,k1) + glift*fluxEz/2;
  end
end
return;


text:
function scale = p25_scale ( neqn )

%*****************************************************************************80
%
%% P25_SCALE returns scale factors for problem p25.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = [ 1.0; 1.0 ];

  return
end


text:
function [eventLog] = ma_executeCoast_nBody_goto_soi_trans(initialState, eventNum, forceModel, soiSkipIds, massLoss, maxPropTime, events, celBodyData)
%ma_executeCoast_nBody_goto_soi_trans Summary of this function goes here
%   Detailed explanation goes here

    eventLog = ma_executeCoast_nBody_goto_dt(maxPropTime, initialState, eventNum, forceModel, true, soiSkipIds, massLoss, maxPropTime, events, false, celBodyData);
end


text:
function vw = roiSetVertInds(vw, prefs)
%  Add a field to each current ROI with the mesh vertices needed to show
%  the ROI. By storing these values in the ROI field, they do not have to
%  be computed each time the mesh is updated. This significantly speeds the
%  update process.
%
% vw = roiSetVertInds([vw], [prefs])
%
% See roiRemoveVertInds.m
%
% Example: vw = roiSetVertInds(vw);
%
% 8/2009: JW

%% Variable check
if ~exist('vw', 'var'),
    % this takes a LONG time!
    vw = getCurView;
end

if ~exist('prefs', 'var'),
    prefs = mrmPreferences;
end

msh  = viewGet(vw, 'currentmeshnumber');
if isempty(msh), warning('need a mesh to get ROI vertex indices'); return; end

%% Set the ROI vert indices
nROIs = length(viewGet(vw, 'ROIs'));
verbose = prefsVerboseCheck;

for roi = 1:nROIs

    % Because this code can be slow, display some progress...
    if verbose
        fprintf('[%s]: Adding vertex indices to ROI %d\n', mfilename, roi);
    end

    roiVertInds = viewGet(vw, 'roiVertInds', roi, prefs);

    vw = viewSet(vw, 'roivertinds', roiVertInds, roi, prefs);

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function X = update_particles(F_update, Xstd_pos, Xstd_vec, X)

N = size(X, 2);

X = F_update * X;

X(1:2,:) = X(1:2,:) + Xstd_pos * randn(2, N);
X(3:4,:) = X(3:4,:) + Xstd_vec * randn(2, N);

author: Sable

repo: mcbench-benchmarks


text:
function r = randpick(x)
% RANDPICK - Pick element from x randomly
%            If x is matrix, pick row from x randomly.

%   Author: Aki Vehtari <Aki.Vehtari@hut.fi>
%   Last modified: 2004-09-07 11:25:23 EEST

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if any(size(x)==1)
  r=x(floor(rand.*length(x)+1));
else
  r=x(floor(rand.*size(x,1)+1),:);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = vec2pdfb(c, s)
% VEC2PDFB   Convert the vector form to the output structure of the PDFB
%
%       y = vec2pdfb(c, s)
%
% Input:
%   c:  1-D vector that contains all PDFB coefficients
%   s:  structure of PDFB output
%
% Output:
%   y:  PDFB coefficients in cell vector format that can be used in pdfbrec
%
% See also:	PDFB2VEC, PDFBREC

% Copy the coefficients from c to y according to the structure s
n = s(end, 1);      % number of pyramidal layers
y = cell(1, n);

% Variable that keep the current position
pos = prod(s(1, 3:4));

y{1} = reshape(c(1:pos), s(1, 3:4));

% Used for row index of s
ind = 1;

for l = 2:n
    % Number of directional subbands in this layer
    nd = length(find(s(:, 1) == l));

    y{l} = cell(1, nd);

    for d = 1:nd
        % Size of this subband
        p = s(ind + d, 3);
        q = s(ind + d, 4);
        ss = p * q;

        y{l}{d} = reshape(c(pos+[1:ss]), [p, q]);
        pos = pos + ss;
    end

    ind = ind + nd;
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [A,idxList]=quickSort(A,byCol,gtCompareFunc)
%%QUICKSORT Sort an array, a matrix, or a cell array in ascending order
%           using a non-recursive implementation of the quicksort
%           algorithm. The quicksort algorithm tends to have better average
%           performance than other techniques, being O(n*log(n)) averge and
%           best case, with a small multiplicative constant, but has O(n^2)
%           worst case performance, which is worse than heapsort.  A custom
%           comparison function can be used so that one can, for example,
%           sort strings or other things. Changing the direction of the
%           comparison function (from greater than to less than) changes
%           the sorting order from ascending to descending. Note that
%           quickSort is not a stable sorting algorithm, meaning that the
%           order of items having the same value might change.
%
%INPUTS: A  An array, 2D matrix, or linear cell array that is to be sorted.
%           When given a matrix, the comparison function gtCompareFunc must
%           be provided so that it is clear how the columns are compared
%           (or the rows if byCol is false).
%     byCol A boolean value indicating whether sorting should be performed
%           by row or by column. This must be specified when a matrix is
%           passed as as a linear array could be confused with a matrix
%           with just one row/column. The default if this parameter is
%           omitted and the input is 1-dimensional is whatever would sort
%           over the 1D array. If the input is 2D, then the default is
%           true. This parameter is not used when cell arrays are passed
%           (an empty matrix can be passed for this) as only linear cell
%           arrays are supported, so no ambiguity between 1D and 2D inputs
%           would exist.
%gtCompareFunc  A function handle that performs a greather-than comparison
%               of two entries in A. This lets one define custom comparison
%               operations. Providing a less-than comparison for this will
%               cause the list to be sorted in descending order rather than
%               ascending order. When sorting an array or matrix, the
%               the function handle takes inputs of the form
%               gtCompareFunc(A(:,i),A(:,j)) if byCol=true and with
%               reversed row and column indices if byCol=false. When
%               comparing cell arrays, it must handle inputs of the form
%               gtCompareFunc(A{i},A{j}). The default if this parameter is
%               omitted is @(x1,x2)(x1>x2);
%
%OUTPUTS: A The sorted array/ matrix/ cell array. If gtCompareFunc
%           performs a greater-than comparison, then it is in increasing
%           order. Otherwise it is in decreasing order.
%   idxList The indices of the original elements with respect to the sorted
%           order. For example, if the input A is an array, then
%           A(idxList) on the input A will give the sorted output A.
%
%The quicksort algorithm is based on the description given in Chapter 7 of
%[1]. However, it has been modified to eliminate the recursion.
%
%Note that sorting large matrices can be slow as each element (so an entire
%row/ column) is copied during the search. Thus, if one is just sorting
%according to a particular row, it makes sense ot get the idxList for that
%one row and then use it to sort everything else.
%
%REFERENCES:
%[1] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein,
%    Introduction to Algorithms, 2nd ed. Cambridge, MA: The MIT Press,
%    2001.
%
%December 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %If an empty matrix is passed, return an empty matrix.
    if(isempty(A))
        idxList=[];
        return;
    end

    if(nargin<2||isempty(byCol))
        %If the input is 1-dimensional, then make the byCol parameter go
        %over whichever dimension would sort it. Otherwise, make it go over
        %columns.
        if(size(A,1)==1)
            byCol=true;
        elseif(size(A,2)==1)
            byCol=false;
        else
            byCol=true;
        end
    end

    if(nargin<3||isempty(gtCompareFunc))
        gtCompareFunc=@(x1,x2)(x1>x2);
    end

    isACellArray=isa(A,'cell');
    if(~isACellArray)
        %If A is supposed to be sorted by row.
        if(byCol==false)
            A=A';
        end

        numPoints=size(A,2);
    else
        numPoints=length(A);
    end

    %The index list is not used in the algorithm for sorting; it is just
    %computed in case the user wants it as a return variable.
    idxList=1:numPoints;

    %stack(:,i) holds the lower and upper bounds of the partiion being
    %sorted at the ith level. The upper bound on the stack space in
    %the worst-case scenario is where every single point visited
    %creates two partitions. However, since each time it splits, it
    %deletes itself, the limit on the memory is ceil(numPoints/2) and
    %not numPoints.
    stack=zeros(2,ceil(numPoints/2));

    p=1;
    r=numPoints;
    stackIdx=1;
    stack(:,stackIdx)=[p;r];
    while(stackIdx>0)
        %Pop the p and r values off of the stack.
        p=stack(1,stackIdx);
        r=stack(2,stackIdx);
        stackIdx=stackIdx-1;

        %Partition the array, getting the correct location of the pivor
        %element.
        if(~isACellArray)
            [q,A,idxList]=partition(A,idxList,p,r,gtCompareFunc);
        else
            [q,A,idxList]=partitionCell(A,idxList,p,r,gtCompareFunc);
        end

        %If there are elements on the left side, then push them onto
        %the stack.
        if(p<q-1)
            stackIdx=stackIdx+1;
            stack(:,stackIdx)=[p;q-1];
        end

        %If there are elements on the right side, then push them onto
        %the stack.
        if(q+1<r)
            stackIdx=stackIdx+1;
            stack(:,stackIdx)=[q+1;r];
        end
    end

    %Restore the orientation if it was supposed to be sorted by row.
    if(~isACellArray&&byCol==false)
        A=A';
    end
end

function [q,A,idxList]=partition(A,idxList,p,r,gtCompareFunc)
%Rearrange the subarray A(p:r) in place.
    x=A(:,r);
    i=p-1;
    for j=p:(r-1)
        if(~gtCompareFunc(A(:,j),x))
            i=i+1;
            %Swap A(:,i) and A(:,j)
            temp=A(:,i);
            A(:,i)=A(:,j);
            A(:,j)=temp;

            %Swap the indices in idxList to keep track of changes in the
            %ordering in A.
            temp=idxList(i);
            idxList(i)=idxList(j);
            idxList(j)=temp;
        end
    end

    %Swap A(:,i+1) and A(:,r)
    temp=A(:,i+1);
    A(:,i+1)=A(:,r);
    A(:,r)=temp;

    %Record the swap in the index list.
    temp=idxList(i+1);
    idxList(i+1)=idxList(r);
    idxList(r)=temp;

    %The new pivot point to return.
    q=i+1;
end

function [q,A,idxList]=partitionCell(A,idxList,p,r,gtCompareFunc)
%This is the same as the partition function, but has been modified to index
%cell arrays.
    x=A{r};
    i=p-1;
    for j=p:(r-1)
        if(~gtCompareFunc(A{j},x))
            i=i+1;
            %Swap A{i} and A{j}
            temp=A{i};
            A{i}=A{j};
            A{j}=temp;
            %Record the swap in the index list.
            temp=idxList(i);
            idxList(i)=idxList(j);
            idxList(j)=temp;
        end
    end

    %Swap A{i+1} and A{r}
    temp=A{i+1};
    A{i+1}=A{r};
    A{r}=temp;

    %Record the swap in the index list.
    temp=idxList(i+1);
    idxList(i+1)=idxList(r);
    idxList(r)=temp;

    %The new pivot point to return.
    q=i+1;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [g1, g2, g3] = sdlfmaXsdlfmKernGradientBlockILJ(lfmKern1, lfmKern2, ...
    t1, t2, i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
    g2Mean, gsp1Mean, gsp2Mean, typeMeanParam, typeMeanSwitching, ...
    typeKernParam1, typeKernParam2, typeKernSwitching1, typeKernSwitching2)

% SDLFMAXSDLFMKERNGRADIENTBLOCKILJ
% FORMAT
% DESC computes the gradients of the parameters for system 1 and system 2
% when i is lower than j. System 1 represents an acceleration while system
% 2 represents a position.
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG covGrad : partial derivatives of the objective function wrt portion
% of the corresponding kernel matrix
% ARG g1Mean : gradients of the parameter of the system 1 obtained from the
% part of the function that uses the funcitons accommpanying the initial
% conditions.
% ARG g2Mean : gradients of the parameter of the system 2 obtained from the
% part of the function that uses the funcitons accommpanying the initial
% conditions.
% ARG gsp1Mean : gradient of the switching point of the system 1 obtained
% from the part of the function that uses the funcitons accommpanying the
% initial conditions.
% ARG gsp2Mean : gradient of the switching point of the system 2 obtained
% from the part of the function that uses the funcitons accommpanying the
% initial conditions.
% ARG typeMeanParam : specify the mean functions used to compute this part
% of the kernel
% ARG typeMeanSwitching : specify the functions used to compute the
% gradients of the swicthing points in part thst uses mean functions
% ARG typeKernParam1 : specify the first kernel function used to compute
% this part of the kernel
% ARG typeKernParam2 : specify the second kernel function used to compute
% this part of the kernel
% ARG typeKernSwitching1 : specify the functions used to compute the
% gradients of the swicthing points in both sides of the kernel function 1
% ARG typeKernSwitching2 : specify the functions used to compute the
% gradients of the swicthing points in both sides of the kernel function 2
% RETURN g1 : gradients of parameters for the system 1
% RETURN g2 : gradients of parameters for the system 2
% RETURN g3 : gradients of switching points
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin < 14
    typeMeanParam = 'sdlfm';
    typeMeanSwitching = 'sdlfmv';
    typeKernParam1 = 'lfmXlfm';
    typeKernParam2 = 'lfmvXlfm';
    typeKernSwitching1= {'lfmvXlfm', 'lfmvXlfm'};
    typeKernSwitching2 ={'lfmvXlfmv', 'lfmaXlfm'};
end

g3 = [];

fhandleMeanParam = str2func([typeMeanParam 'MeanCompute']);
fhandleMeanGradParam = str2func([typeMeanParam 'MeanGradient']);
fhandleMeanGradSwitching = str2func([typeMeanSwitching 'MeanCompute']);
fhandleKernPosPos = str2func([typeKernParam1 'KernCompute']);
fhandleKernGradPosPos = str2func([typeKernParam1 'KernGradient']);
fhandleKernGradSwitchingPosPos1 = str2func([typeKernSwitching1{1} 'KernCompute']);
fhandleKernGradSwitchingPosPos2 = str2func([typeKernSwitching1{2} 'KernCompute']);
fhandleKernVelPos = str2func([typeKernParam2 'KernCompute']);
fhandleKernGradVelPos = str2func([typeKernParam2 'KernGradient']);
fhandleKernGradSwitchingVelPos1 = str2func([typeKernSwitching2{1} 'KernCompute']);
fhandleKernGradSwitchingVelPos2 = str2func([typeKernSwitching2{2} 'KernCompute']);

c2 = fhandleMeanParam(lfmKern2(1), t2, 'Pos');
e2 = fhandleMeanParam(lfmKern2(1), t2, 'Vel');

if isempty(generalConst{i,j})
    coeffPosPos = c2;
    coeffPosVel = e2;
else
    coeffPosPos = generalConst{i,j}(1,1)*c2 + generalConst{i,j}(2,1)*e2;
    coeffPosVel = generalConst{i,j}(1,2)*c2 + generalConst{i,j}(2,2)*e2;
end
g1Kern = zeros(length(lfmKern1), 5);
g2Kern = zeros(length(lfmKern1), 5);
PosPos = zeros(length(t1),1);
PosVel = zeros(length(t1),1);
for k=1:length(lfmKern1)
    PosPos = PosPos + fhandleKernPosPos(lfmKern1(k), lfmKern2(k), t1,...
        lfmKern1(k).limit);
    [g1KLocal, g2KLocal] = fhandleKernGradPosPos(lfmKern1(k), lfmKern2(k), ...
        t1, lfmKern1(k).limit, covGrad, coeffPosPos.');
    g1Kern(k,:) = g1Kern(k, :) + g1KLocal;
    g2Kern(k,:) = g2Kern(k, :) + g2KLocal;
    PosVel = PosVel + fhandleKernVelPos(lfmKern1(k), lfmKern2(k), ...
        t1, lfmKern1(k).limit);
    [g1KLocal, g2KLocal] = fhandleKernGradVelPos(lfmKern1(k), lfmKern2(k), ...
        t1, lfmKern1(k).limit, covGrad, coeffPosVel.');
    g1Kern(k,:) = g1Kern(k, :) + g1KLocal;
    g2Kern(k,:) = g2Kern(k, :) + g2KLocal;
end
[gcAlpha, gcOmega] = fhandleMeanGradParam(lfmKern2(1), t2, 'Pos');
[geAlpha, geOmega] = fhandleMeanGradParam(lfmKern2(1), t2, 'Vel');
[gradAlpha, gradOmega] = getLocalGradAlphaOmega(lfmKern2);
g2Pos = fhandleMeanGradSwitching(lfmKern2(1), t2, 'Pos');
h2Vel = fhandleMeanGradSwitching(lfmKern2(1), t2, 'Vel');
gsp1Kern = zeros(1, length(lfmKern1));
gsp2Kern = zeros(1, length(lfmKern1));
% This means that are only two switching points involved, t1
% and t0 appear in k, like k_{ff}(t1-t0,t'-t0) and
% k_{mf}(t1-t0,t'-t0). The other points appear in c1(t - t1)
% and e1(t - t1)
for k=1:length(lfmKern1)
    temp = fhandleKernGradSwitchingPosPos1(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosPos.').*covGrad));
    %temp = fhandleKernGradSwitchingPosPos2(lfmKern2(k), lfmKern1(k), lfmKern1(k).limit, t1)';
    temp = fhandleKernGradSwitchingPosPos2(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosPos.').*covGrad));
    gsp2Kern(k) = gsp2Kern(k) + sum(sum((temp*coeffPosPos.').*covGrad));
    temp = fhandleKernGradSwitchingVelPos1(lfmKern1(k), lfmKern2(k), t1, lfmKern1(k).limit);
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosVel.').*covGrad));
    temp = fhandleKernGradSwitchingVelPos2(lfmKern2(k), lfmKern1(k), lfmKern1(k).limit, t1).';
    gsp1Kern(k) = gsp1Kern(k) - sum(sum((temp*coeffPosVel.').*covGrad));
    gsp2Kern(k) = gsp2Kern(k) + sum(sum((temp*coeffPosVel.').*covGrad));
end
if isempty(generalConst{i,j})
    matGradAlpha = PosPos*gcAlpha.' + PosVel*geAlpha.';
    matGradOmega = PosPos*gcOmega.' + PosVel*geOmega.';
    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    matGradSp2 = PosPos*g2Pos.' + PosVel*h2Vel.';
    gsp2 = - sum(sum(matGradSp2.*covGrad));
    g3(i) = gsp1Mean + sum(gsp1Kern);                % switching points 1
    g3(j) = gsp2Mean + sum(gsp2Kern) + gsp2;         % switching points 2
else
    constGradAlpha = generalConstGrad{1}{i,j};
    constGradOmega = generalConstGrad{2}{i,j};
    constVal = generalConst{i,j};
    % Derivative wrt parameters
    matGradAlpha = PosPos*((constVal(1,1)*gcAlpha.' + constGradAlpha(1,1)*c2.') ...
        + constVal(2,1)*geAlpha.' + constGradAlpha(2,1)*e2.') ...
        + PosVel*((constVal(1,2)*gcAlpha.' + constGradAlpha(1,2)*c2.') ...
        + constVal(2,2)*geAlpha.' + constGradAlpha(2,2)*e2.');

    matGradOmega = PosPos*((constVal(1,1)*gcOmega.' + constGradOmega(1,1)*c2.') ...
        + constVal(2,1)*geOmega.' + constGradOmega(2,1)*e2.') ...
        + PosVel*((constVal(1,2)*gcOmega.' + constGradOmega(1,2)*c2.') ...
        + constVal(2,2)*geOmega.' + constGradOmega(2,2)*e2.');

    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    % Derivative wrt switching points
    % Firts, notice that gsp2Kern doesn't correspond to the switching point
    % of the interval j (or say 2). It's just the derivative of the inner
    % switching point that lead to the actual switching point for the
    % interval j. So we rename it first.
    gspInt = sum(gsp2Kern);
    % Compute the derivative of the swicthing point j, not appearing in the
    % constant
    matGradSp2 = PosPos*(constVal(1,1)*g2Pos.' + constVal(2,1)*h2Vel.') + ...
        PosVel*(constVal(1,2)*g2Pos.' + constVal(2,2)*h2Vel.');
    gsp2 = - sum(sum(matGradSp2.*covGrad));
    % Compute the derivatives for all the switching points apearing in the
    % constant
    constGradSPoint = generalConstGrad{3}{i,j};
    numberSP = size(constGradSPoint,2);
    gspInBetween = zeros(1, numberSP);
    for k=1:numberSP
        temp = PosPos*(constGradSPoint(1,k)*c2.' + constGradSPoint(3,k)*e2.') + ...
            PosVel*(constGradSPoint(2,k)*c2.' + constGradSPoint(4,k)*e2.');
        gspInBetween(k) = sum(sum(temp.*covGrad));
    end
    % Assign derivatives wrt all other switching points, with a correction
    % for the derivative of the innermost switching point in the constant
    gspInBetween(end) =  gspInBetween(end) + gspInt;
    g3(i) = gsp1Mean + sum(gsp1Kern);           % switching point 1
    g3(j) = gsp2Mean + gsp2 + gspInBetween(1);  % switching point 2
    g3(i+1:j-1) = fliplr(gspInBetween(2:end));
end
% Assign derivatives wrt first system
g1{1} = g1Mean + sum(g1Kern(:,1:3), 1);          % mass 1, spring 1, damper 1
g1{2} = g1Kern(:,4).';                            % inverse widths
g1{3} = g1Kern(:,5).';                            % sensitivities 1
% Assign derivatives wrt second system
g2{1} = g2Mean + sum(g2Kern(:,1:3), 1) + gCoeff; % mass 2, spring 2, damper 2
g2{2} = g2Kern(:,4).';                            % inverse widths
g2{3} = g2Kern(:,5).';                            % sensitivities 2

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% compute the eccentricity of a 2D shape

path(path, 'toolbox/');
path(path, 'data/');

rep = 'images/';
name = 'bird';
name = 'apple';
name = 'device';
name = 'camel';
name = 'cavern';
name = 'chicken';
name = 'giraffe';

method = 'ecc';
method = 'mean';
metric = 'euclidean';
metric = 'geodesic';
options.metric = metric;
options.method = method;

% exponent for the computation
s = Inf;
if strcmp(method, 'mean')
    s = 1;
end

rep = 'results/eccentricity/';
if ~exist(rep)
    mkdir(rep);
end
testname = [method '-' metric];

n = 256;
M = rescale( load_image(name,n), 0,1 );
M = double(M>0.5);

% make sure pixels on the boundary are black
if M(1)==1
    M = 1-M;
end

% compare the geodesic distance to the euclidean distance
clf;
imagesc(M); axis image; axis off;
title('click on a point inside the shape');
% [y,x] = ginput(1);
x = 50; y = 50;
start_points = round([x y]');
W = ones(n);
L = zeros(n)-Inf; L(M==1) = +Inf;
options.constraint_map = [];
[D0,S,Q] = perform_fast_marching(W, start_points, options);
options.constraint_map = L;
[D1,S,Q] = perform_fast_marching(W, start_points, options);
D0(M==0) = 0; D1(M==0) = 0;

clf;
subplot(1,2,1);
hold on;
display_eccentricity(D0, 1);
h = plot(y,x, 'r.');
set(h, 'MarkerSize', 20);
axis ij;
title('Euclidean distance');
subplot(1,2,2);
hold on;
display_eccentricity(D1, 1);
h = plot(y,x, 'r.');
set(h, 'MarkerSize', 20);
axis ij;
title('Geodesic distance');

% saveas(gcf, [rep name '-geodesic-dist.eps'], 'epsc');

options.s = s;
options.nb_samples = 300; % subamples
[E,Ind] = compute_eccentricity_transform(M, options);

clf;
subplot(1,2,1);
display_eccentricity(E, 1);
title('Eccentricity Transform');

if strcmp(method, 'ecc')
    [ep,I,J] = unique( Ind(:) ); ep(ep==Inf) = [];
    [eccx,eccy] = ind2sub([n n], ep);
    J = reshape(J,n,n);

    subplot(1,2,2);
    hold on;
    imagesc(J);
    h = plot(eccy,eccx,'r.');
    set(h, 'MarkerSize', 20);
    hold off;
    title('Eccentric points');
    axis image; axis off; axis ij;

    % saveas(gcf, [rep name '-' testname '.png'], 'png');
    % saveas(gcf, [rep name '-' testname '.eps'], 'epsc');

    % subplot(2,2,3);
    % display_eccentricity(E, 2);
end

% display the histograms
clf;
I = find(E>0);
v = rescale(E(I));
hist(v, 60);
title('Histogram');
axis tight;
% saveas(gcf, [rep name '-histogram.eps'], 'epsc');
saveas(gcf, [rep name '-histogram-' method '.png'], 'png');

%% save as pretty images
[Y,X] = meshgrid(1:n,1:n);
r = 0.02*n;
% save shape
warning off;
imwrite( 1-M, [rep name '-shape.png'], 'png' );
% save eccentricity
E0 = E; E0(M==0) = Inf;
E0 = convert_distance_color(E0);
% place a dot at the center
[tmp,i] = min(E(:) + (E(:)==0)*1e9); [x,y] = ind2sub([n n], i);
I = find( (X-x).^2 + (Y-y).^2 <= r^2 );
for i=1:3
    A = E0(:,:,i); A(I) = i==1;
    E0(:,:,i) = A;
    A = E0(:,:,i); A(I) = i==1;
    E0(:,:,i) = A;
end
imwrite( E0, [rep name '-shape-' method '.png'], 'png' );
if strcmp(method, 'ecc')
    % save eccentric points
    E0 = Ind;
    [tmp,a,E0] = unique(E0); E0 = reshape(E0,n,n);
    E0(M==0) = Inf;
    E0 = convert_distance_color(E0);
    for k=1:length(eccx)
        I = find( (X-eccx(k)).^2 + (Y-eccy(k)).^2 <= r^2 );
        for i=1:3
            A = E0(:,:,i); A(I) = i==1;
            E0(:,:,i) = A;
            A = E0(:,:,i); A(I) = i==1;
            E0(:,:,i) = A;
        end
    end
    imwrite( E0, [rep name '-shape-eccpoints.png'], 'png' );
    warning on;
end

author: alecjacobson

repo: gptoolbox


text:
filename='Gripping_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.75;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function hFigBlur = blurFigure(hFig, state)
% blurFigure blurs and prevents interaction on a figure window
%
% Syntax:
%    hFigBlur = blurFigure(hFig, state)
%
% Description:
%    blurFigure(hFig) blurs figure hFig and prevents interaction with it.
%    The only interaction possible is with user-created controls on the
%    blurring panel (see below).
%
%    hFigBlur = blurFigure(hFig) returns the overlaid blurred figure pane.
%    This is useful to present a progress bar or other GUI controls, for
%    user interaction during the blur phase.
%
%    blurFigure(hFig,STATE) sets the blur status of figure hFig to STATE,
%    where state is 'on','off',true or false (default='on'/true).
%
%    blurFigure(hFig,'on') or blurFigure(hFig,true) is the same as:
%    blurFigure(hFig).
%
%    blurFigure(hFig,'off') or blurFigure(hFig,false) is the same as:
%    close(hFigBlur).
%
%    blurFigure('demo') displays a simple demo of the blurring
%
% Input parameters:  (all parameters are optional)
%
%    hFig -  (default=gcf) Handle(s) of the modified figure(s).
%            If component handle(s) is/are specified, then the containing
%            figure(s) will be inferred and used.
%
%    state - (default='on'/true) blurring flag: 'on','off',true or false
%
% Examples:
%    hFigBlur = blurFigure(hFig);       % blur hFig (alternative #1)
%    hFigBlur = blurFigure(hFig,true);  % blur hFig (alternative #2)
%
%    blurFigure(hFig,false);       % un-blur hFig (alternative #1)
%    blurFigure(hFig,'off');       % un-blur hFig (alternative #2)
%    close(hFigBlur);              % un-blur hFig (alternative #3)
%    delete(hFigBlur);             % un-blur hFig (alternative #4)
%
%    blurFigure('demo');           % blur demo with progress bar etc.
%
%    hFigBlur = blurFigure(hFig);  % add a blur with a cancel button
%    uicontrol('parent',hFigBlur, 'string','Cancel', 'Callback','close(gcbf)');
%
% Technical Description:
%    http://UndocumentedMatlab.com/blog/blurred-matlab-figure-window
%
% Bugs and suggestions:
%    Please send to Yair Altman (altmany at gmail dot com)
%
% Warning:
%    This code heavily relies on undocumented and unsupported Matlab functionality.
%    It works on Matlab 7.9 (R2009b) and higher, but use at your own risk!
%
% Change log:
%    2011-03-07: First version posted on Matlab's File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/?term=authorid%3A27420
%    2011-03-08: Minor fix to the demo
%    2011-10-14: Fix for R2011b
%
% See also:
%    enableDisableFig, setFigTransparency, getJFrame (all of them on the File Exchange)

% Programmed by Yair M. Altman: altmany(at)gmail.com
% $Revision: 1.2 $  $Date: 2011/11/14 03:13:27 $

  % Require Java 1.6.0_10 to run
  if ~usejava('awt')
      error([mfilename ' requires Java to run']);
  elseif ~ismethod('com.sun.awt.AWTUtilities','setWindowOpacity')
      error([mfilename ' requires a newer Java engine (1.6.0_10 minimum) to run']);
  end

  % Set default input parameters values
  if nargin < 1,  hFig = gcf;    end
  if nargin < 2,  state = true;  end

  % Check for demo request
  if ischar(hFig)
      if strcmpi(hFig,'demo')
          f1 = figure;
          try oldWarn = warning('off','MATLAB:uitree:MigratingFunction'); catch, end
          h1 = uitree('root','c:\');  %#ok unused
          drawnow;
          try h1.getTree.expandRow(0); catch, end
          try warning(oldWarn); catch, end
          h2 = uicontrol('string','click me!', 'units','pixel', 'pos',[300,50,100,20]);  %#ok unused
          hax = axes('parent',gcf,'units','pixel','pos',[230,100,300,300]);  %#ok unused
          surf(peaks);
          set(gcf,'ToolBar','figure');  % restore the toolbar that was removed by the uicontrol() call above
          hFigBlurTemp = blurFigure(f1);
          h4 = uicontrol('parent',hFigBlurTemp, 'style','text', 'units','pixel', 'pos',[130,85,390,80], 'string','Processing - please wait...', 'FontSize',12, 'FontWeight','bold','Fore','red','Back','yellow');  %#ok unused
          [h5a,h5b]=javacomponent('javax.swing.JProgressBar',[180,115,310,20],hFigBlurTemp); h5a.setValue(67);  %#ok unused
          h6 = uicontrol('parent',hFigBlurTemp, 'string','Cancel', 'pos',[280,90,100,20], 'Callback','close(gcbf)');  %#ok unused
          return;
      else
          error([mfilename ' requires a figure handle or the string ''demo''']);
      end
  end

  % Check valid state values
  if ischar(state)
      if ~any(strcmpi(state,{'on','off'}))
          error([mfilename ' requires a STATE value of ''on'', ''off'', true or false']);
      else
          state = strcmpi(state,'on');
      end
  elseif ~isscalar(state)
      error([mfilename ' requires a scalar STATE value of ''on'', ''off'', true or false']);
  else
      state = logical(state);
  end

  % Get unique figure handles
  if ~all(ishghandle(hFig))
      error('hFig must be a valid GUI handle or array of handles');
  end
  hFigCell = num2cell(zeros(1,length(hFig)));
  for hIdx = 1 : length(hFig)
      if hFig(hIdx) ~= 0
          hFigCell{hIdx} = ancestor(hFig(hIdx),'figure');
      end
  end
  [unigueHandles,uniqueIdx,foldedIdx] = unique(hFig,'first');  %#ok unused
  hFigBlur_internal = zeros(size(hFig));

  % Loop over all requested figures
  for figIdx = 1 : length(hFigCell)

      % Get the root Java frame
      try
          % The following is only used for the internal sanity checks on the target figure
          jff = getJFrame(hFigCell{figIdx});
      catch
          errMsg = regexprep(lasterr,'Error using [^\n]+\n','');  %#ok
          error(errMsg);
      end

      % If this figure has already been processed - bail out
      if ~ismember(figIdx,uniqueIdx)
          hFigBlur_internal(figIdx) = hFigBlur_internal(foldedIdx(figIdx));
          continue;
      end

      hFigBlurTemp = getappdata(hFigCell{figIdx}, 'blur_hFig');

      % If blurring is requested
      if state

          % If this figure is already blurred, bail out
          if ~isempty(hFigBlurTemp) && ishghandle(hFigBlurTemp),  continue;  end

          % Overlap the current frame with an identical empty figure frame
          hProps = get(hFigCell{figIdx});
          hFigBlur_internal(figIdx) = figure('Units',hProps.Units, 'Pos',hProps.Position, 'Color',hProps.Color, ...
                                             'Menu',hProps.MenuBar, 'Toolbar',hProps.ToolBar, ...
                                             'IntegerHandle','off','HandleVisibility','off', ...
                                             'Name','', 'NumberTitle','off', ...
                                             'DockControls','off','Resize','off');

          % Remove all toolbar sub-components from the new figure
          hToolbar = findall(hFigCell{figIdx},'tag','FigureToolBar');
          if isempty(hToolbar)
              set(hFigBlur_internal(figIdx),'ToolBar','none');  % might be 'auto' without a visible toolbar, so we need this...
          else
              hToolbar = findall(hFigBlur_internal(figIdx),'tag','FigureToolBar');
              hComponents = findall(hToolbar);
              delete(hComponents(2:end));  % keep the toolbar without any sub-components, to preserve the bar area
          end

          % Remove all menubar sub-components from the new figure
          hMenuBar = findall(hFigCell{figIdx},'type','uimenu','visible','on');
          if isempty(hMenuBar)
              set(hFigBlur_internal(figIdx),'Menu','none');  % might be 'figure' without any visible menu item, so we need this...
          else
              hComponents = findall(hFigBlur_internal(figIdx),'type','uimenu');
              set(hComponents,'Label','');  % keep the menubar without any visible sub-components, to preserve the bar area
          end

          % Blur by setting the new frame's opacity to translucent
          jNewFig = getJFrame(hFigBlur_internal(figIdx));
          com.sun.awt.AWTUtilities.setWindowOpacity(jNewFig,0.8);

          % Add a hidden handle to the corresponding figure in each of the figures
          setappdata(hFigBlur_internal(figIdx), 'blur_hFig',hFigCell{figIdx});
          setappdata(hFigCell{figIdx}, 'blur_hFig',hFigBlur_internal(figIdx));

          % Prevent resizing during blur
          setappdata(hFigBlur_internal(figIdx), 'blur_originalResize', hProps.Resize);
          set(hFigCell{figIdx}, 'Resize','off');

          % Prevent closing the figures via the 'X' button
          set(hFigBlur_internal(figIdx), 'CloseRequestFcn',@closeRequestFcn)

          % Attach listeners to the original and new frames: Position, Visibility, Docking, Close
          hFigures = [hFigCell{figIdx}, hFigBlur_internal(figIdx)];
          hLink(1) = linkprop(hFigures,'Units');
          hLink(2) = linkprop(hFigures,'Position');
          hLink(3) = linkprop(hFigures,'Visible');
          hLink(4) = linkprop(hFigures,'WindowStyle');
          hLink(5) = handle.listener(hFigCell{figIdx},'ObjectBeingDestroyed',{@deleteFcn,hFigBlur_internal(figIdx)});
          hLink(6) = handle.listener(hFigBlur_internal(figIdx),'ObjectBeingDestroyed',@unblurFcn);
          hLink(7) = handle.listener(hFigCell{figIdx},findprop(handle(hFigCell{figIdx}),'Visible'),'PropertyPostSet',{@visibilityFcn,hFigBlur_internal(figIdx)});

          % Persist the listener handles in the blur figure's hidden ApplicationData
          setappdata(hFigBlur_internal(figIdx), 'blur_hListeners',hLink);

          % Prevent Mininization
          try
              % This prevents minimization better, but also prevents interaction with
              % any user-generated uicontrols in the blur window - not good!
              %jNewFig.setEnabled(false);

              % This simply restores the window back to place after a minimization
              % - ugly but better than nothing...
              hjNewFig = handle(jNewFig,'CallbackProperties');
              set(hjNewFig, 'WindowIconifiedCallback',@minimizedFcn);
          catch
              % might fail in some future Matlab version - never mind...
          end

          % Prevent direct activation/invocation of the original figure
          try
              hjFig = handle(jff,'CallbackProperties');
              set(hjFig, 'WindowActivatedCallback',{@visibilityFcnHelper,hFigBlur_internal(figIdx)});
          catch
              % might fail in some future Matlab version - never mind...
          end

      else  % Unbluring was requested
          try delete(hFigBlurTemp); catch, end
      end
  end

  % Return the blurring figure pane handle(s), if requested
  if nargout && state
      hFigBlur = hFigBlur_internal;
  end

end  % blurFigure

%% Get the root Java frame (up to 10 tries, to wait for figure to become responsive)
function jframe = getJFrame(hFig)

  drawnow;

  % Ensure that hFig is a non-empty handle...
  if isempty(hFig)
      error(['Cannot retrieve the figure handle for handle ' num2str(hFig)]);
  end

  % Check for the desktop handle
  if hFig == 0
      error('Blurring the Desktop is not supported');
  end

  % Check whether the figure is invisible
  if strcmpi(get(hFig,'Visible'),'off')
      error(['Cannot blur non-visible figure ' num2str(hFig)]);
  end

  % Check whether the figure is docked
  if strcmpi(get(hFig,'WindowStyle'),'docked')
      error(['Cannot blur docked figure ' num2str(hFig)]);
  end

  % Retrieve the figure window (JFrame) handle
  jframe = [];
  maxTries = 10;
  while maxTries > 0
      try
          % Get the figure's underlying Java frame
          jf = get(handle(hFig),'JavaFrame');

          % Get the Java frame's root frame handle
          %jframe = jf.getFigurePanelContainer.getComponent(0).getRootPane.getParent;
          try
              jframe = jf.fFigureClient.getWindow;  % equivalent to above...
          catch
              jframe = jf.fHG1Client.getWindow;  % equivalent to above...
          end
          if ~isempty(jframe)
              break;
          else
              maxTries = maxTries - 1;
              drawnow; pause(0.1);
          end
      catch
          maxTries = maxTries - 1;
          drawnow; pause(0.1);
      end
  end
  if isempty(jframe)
      error(['Cannot retrieve the Java Frame reference for figure ' num2str(hFig)]);
  end

end  % getJFrame

%% Close request callback function
function closeRequestFcn(hFig, varargin)
  stk = dbstack(1);
  if ~isempty(stk) && strcmpi(stk(1).file,'close.m')
      delete(hFig);
  else
      % came here by trying to close the figure by clicking the 'X' - ignore!
  end
end  % closeRequestFcn

%% Delete blur figure callback function
function deleteFcn(hObj, eventData, hFig, varargin)  %#ok unused
  delete(hFig);
end  % deleteFcn

%% Restore original resize behavior upon un-bluring
function unblurFcn(hFigBlur, varargin)
  try
      originalResizeMode = getappdata(hFigBlur, 'blur_originalResize');
      hFig = getappdata(hFigBlur, 'blur_hFig');
      set(hFig, 'Resize',originalResizeMode);
  catch
      % never mind...
  end
end  % unblurFcn

%% Minimization callback function
function minimizedFcn(hObj, varargin)
  % Restore to previous position (prevent minimization)
  try hObj.setMinimized(false); catch, end
end  % minimizedFcn

%% Gained-focus callback function
function visibilityFcn(hFig, eventData, hFigBlur, varargin)  %#ok unused
  % Set the blurring window on top
  try
      if strcmpi(eventData.NewValue,'on')
          start(timer('ExecutionMode','singleShot', 'StartDelay',0.05, 'TasksToExecute',1, 'TimerFcn', {@visibilityFcnHelper,hFigBlur}));
      end
  catch
      % never mind...
  end
end  % minimizedFcn

%% visibility helper function
function visibilityFcnHelper(hTimer,hEventData,hFigBlur)  %#ok unused
  try
      jNewFig = getJFrame(hFigBlur);
      %figure(hFigBlur);
      jNewFig.requestFocus;
      com.sun.awt.AWTUtilities.setWindowOpacity(jNewFig,0.8);
      hjNewFig = handle(jNewFig,'CallbackProperties');
      set(hjNewFig, 'WindowIconifiedCallback',@minimizedFcn);
      stop(hTimer)
      delete(hTimer);
  catch
      % never mind...
  end
end  % visibilityFcnHelper


text:
function [status, fsHome] = mrvFreeSurferExists(varargin)
% Check whether FREESURFER_HOME is in the matlab environment
%
% Syntax
%   status = mrvFreeSurferExist;
%
% Returns true if FREESURFER_HOME is in the environment.
%
% Inputs
% Optional key/value
%   'verbose' - Prints out the FREESURFER_HOME variable
%
% Returns
% Wandell, Vistasoft Team, 2018
%
% See also
%  mrvFreeSurferConfg (NYI)

%%
p = inputParser;
p.addParameter('verbose',false,@islogical);
p.parse(varargin{:});

%%
fsHome = getenv('FREESURFER_HOME');

if isempty(fsHome),     status = 0;
else,                   status = 1;
end

if p.Results.verbose
    fprintf('\nFREESURFER_HOME is "%s"\n',fsHome);
end

end


text:
function printPostSolver(fid, expInfo)
%PRINTPOSTSOLVER   Print commands after solution has been found
%
% Calling sequence:
%   PRINTPOSTSOLVER(FID, EXPINFO)
% where
%   FID:        ID of a file-writing stream.
%   EXPINFO:    Struct containing information for printing the problem.

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

% Extract info from the EXPINFO struct:
sol = expInfo.sol;
indVarName = expInfo.indVarName;
deInput = expInfo.deInput;
s = expInfo.s;

% Plotting. Deal with scalar and system problems separately.
if ( numel(deInput) == 1 )
    % Scalar problems.
    fprintf(fid, '\n%%%% Plot the solution.\n');
    fprintf(fid, 'waterfall(%s, %s)\n', sol, ...
        indVarName{2});
    fprintf(fid, 'xlabel(''%s''), ylabel(''%s'')', ...
        indVarName{1},indVarName{2});

else
    % Coupled systems.
    fprintf(fid, '\n%%%% Plot the solution components.');
    for k = 1:numel(deInput)
        fprintf(fid, '\nfigure\n');
        fprintf(fid, 'waterfall(%s, %s)\n', s{k}, indVarName{2});
        fprintf(fid, 'xlabel(''%s''), ylabel(''%s''), title(''%s'')', ...
            indVarName{1},indVarName{2},s{k});
    end

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
mex pochhammer.c ../lightspeed/util.obj -I../lightspeed
mex di_pochhammer.c ../lightspeed/util.obj -I../lightspeed
%mex tri_pochhammer.c ../lightspeed/util.obj -I../lightspeed
%mex s_derivatives.c ../lightspeed/util.obj -I../lightspeed

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function P = optimizer(self,x,u)
%OPTIMIZER  Container for optimization problem
%
%   OPT = OPTIMIZER(Problem,x,u) exports an object that contains
%   precompiled numerical data to be solved for varying arguments
%   x, returning the optimal value of the expression u.
%
%   SEE OPTIMIZER for more info
%
%   Example
%
%    The following problem creates an LP with varying upper and lower
%    bounds on the decision variable.
%
%    The optimizing argument is obtained by indexing (with {}) the optimizer
%    object with the point of interest. The argument should be a column
%    vector (if the argument has a width larger than 1, YALMIP assumes that
%    the optimal solution should be computed in several points)
%
%     A = randn(10,3);
%     b = rand(10,1)*19;
%     c = randn(3,1);
%
%     z = sdpvar(3,1);
%     sdpvar UB LB
%
%     Constraints = [A*z < b, LB < z < UB];
%     Objective = c'*z;
%     P = optproblem(Constraints,Objective);
%     % We want the optimal z as a function of [LB;UB]
%     optZ = optimizer(P,[LB; UB],z);
%
%     % Compute the optimal z when LB=1, UB = 3;
%     zopt = optZ{[1; 3]}
%
%     % Compute two solutions, one for (LB,UB) [1;3] and one for (LB,UB) [2;6]
%     zopt = optZ{[[1; 3], [2;6]]}
%
%     A second output argument can be used to catch infeasibility
%     [zopt,infeasible] = optZ{[1; 3]}

P = optimizer(self.Constraints,self.Objective,self.Options,x,u);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% Test file for ADCHEBFUN times

function pass = test_innerProduct

% Initialise pass vector
pass = zeros(3, 5);

% Tolerance for Taylor testing
tolOrder = 1e-2;
tolDiff = 1e-14;

% The function handle we're working with:
func = @innerProduct;

% Compare values
[err, lin] = adchebfun.valueTestingBinary(func);

% Confirm that the error returned is close to zero. Note that the error won't be
% exact zero, due to slightly different mechanisms used in
% CHEBFUN/INNERPRODUCT() and ADCHEBFUN/INNERPRODUCT() for computing the value of
% the inner products.
pass(1, :) = ( err < tolDiff );

% Taylor testing
[order1, order2, nDiff2] = adchebfun.taylorTestingBinary(func);

% We expect all elements of ORDER1 to be close to 1. Since depending on the
% combinations of variables (ADCHEBFUN, CHEBFUN and SCALAR), we should expect
% some of ORDER2 to be noise (for linear operations), in which case NDIFF2
% should have values close to machine epsilon. For the nonlinear operations,
% ORDER2 will be close having all elements value close to 2.

% Check whether we get the expected results for linear and nonlinear
% operations.
linearOpResults = ( ((max(abs(order1(:, 2:end))) - 1) < tolOrder) & ...
    (max(abs(nDiff2(:, 2:end))) < tolDiff) );

nonlinearOpResults = ( (max(abs(order1(:, 1) - 1)) < tolOrder) && ...
    (max(abs(order2(:, 1) - 2)) < tolOrder) );

% Concatenate results
pass(2, :) = [nonlinearOpResults, linearOpResults];

% Check whether we received correct linearity information.
pass(3, :) = ( lin == [0, 1, 1, 1, 1]);
end


text:
function test_bug1607

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_channelrepair ft_topoplotER

% During preprocessing I lost some channels which I got back through
% ft_channelrepair after which the channel order changed per subject.
% I realigned the channel order and the associated data structure of each subject
% to one reference order after which the topoplots look completely different.
% They look even worse than before realignment. Anybody got a clue where the
% confusion arises?
% Attached a script and two data files, one with reference data and one with the
% to be aligned data.

% solution
% was not a bug, instead jonas was indexing the wrong way around:
%
%     nERPdata_nD_left{isubject}.avg(:,:) = ERPdata_nD_left{isubject}.avg(loc,:);     % realign data to ref channel order
% instead of
%     nERPdata_nD_left{isubject}.avg(loc,:) = ERPdata_nD_left{isubject}.avg(:,:);     % realign data to ref channel order

ERPdata_nD_left = [];
reference_labels = [];

%load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1607/06_control_ICA_clean.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1607/ERPdata.mat'));
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1607/reference.mat'));

% reference_labels = ICA_clean.label;

nERPdata_nD_left = ERPdata_nD_left;
% nERPdata_D_left = ERPdata_D_left;
% nERPdata_nD_right = ERPdata_nD_right;
% nERPdata_D_right = ERPdata_D_right;

for isubject = 1:1%length(ERPdata_nD_left)
%     disp(num2str(isubject))
    [a loc] = ismember(ERPdata_nD_left{isubject}.label, reference_labels);  % compare the channel labeling of each subject with a reference
%     for isens = 1 : 64
%         if loc(isens) ~= isens
%             disp(['changing subj ' num2str(isubject) ' sensnr ' num2str(isens)]);
%             break;
%         end
%     end
%     disp(num2str(any(find((ERPdata_nD_left{isubject}.avg == ERPdata_nD_left{isubject}.avg(loc,:))==0))));
    nERPdata_nD_left{isubject}.label = reference_labels;    % relabel the channels according to the reference
%     nERPdata_nD_left{isubject}.avg(:,:) = ERPdata_nD_left{isubject}.avg(loc,:);     % realign data to ref channel order
    nERPdata_nD_left{isubject}.avg(loc,:) = ERPdata_nD_left{isubject}.avg(:,:);     % realign data to ref channel order
    nERPdata_nD_left{isubject}.dof(:,:) = ERPdata_nD_left{isubject}.dof(loc,:);     % realign data to ref channel order

end
% after realigning the data the plots look completely different, compare
% old and new data

cfg = [];
cfg.layout      = 'biosemi64.lay';
cfg.xlim        = [.04 .06];
cfg.channel     = 'all';
cfg.interactive = 'no';
cfg.highlight   = 'yes';
cfg.highlightchannel = reference_labels(end-6:end);

figure; ft_topoplotER(cfg,ERPdata_nD_left{1});      % plot original data
figure; ft_topoplotER(cfg,nERPdata_nD_left{1});     % plot corrected data


text:
function ecef = llh2ecef(llh)
% llh2ecef: converts from navigation coordinates (latitude, longitude and
% altitude) to ECEF coordinates.
%
% INPUTS
%   llh: Nx3 LLH coordinates [lat, lon, h] (rad, rad, m).
%
% OUTPUTS
%   ned: Nx3 NED coordinates [X Y Z] (m, m, m).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%   Guowei Cai et al. Unmanned Rotorcraft Systems. Springer. 2011.
% Eq. 2.22, p. 31.
%
% Version: 001
% Date:    2019/01/16
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

% Preallocate

ecef = zeros(size(llh));

lat = llh(:,1);
lon = llh(:,2);
h   = llh(:,3);

[~,RN] = radius(lat);

e = 0.0818191908426;    % Eccentricity

slat = sin(lat);
clat = cos(lat);
clon = cos(lon);
slon = sin(lon);

ecef(:,1) = (RN + h) .* clat .* clon;
ecef(:,2) = (RN + h) .* clat .* slon;
ecef(:,3) = (RN *(1-e^2) + h) .* slat;

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function X = hnormalise(X)
% 	hnormalise   Normalise homogeneous coordinates to have 1 as the last component
% 		Xn = hnormalise(X)
%
%		Inputs:
%				X		ndims x npts array, where each column is a vector of length ndims,
%							the homogeneous coordinates of a point in ndims-1 space
%
%		Outputs:
%				Xn		ndims x npts array, where each column represents the same point as in X
%							but normalised to have the last component equal to 1

%	Author:		Oliver Whyte <oliver.whyte@ens.fr>
%	Date:		November 2011
%	Copyright:	2011, Oliver Whyte
%	Reference:  O. Whyte, J. Sivic and A. Zisserman. "Deblurring Shaken and Partially Saturated Images". In Proc. CPCV Workshop at ICCV, 2011.
%	URL:		http://www.di.ens.fr/willow/research/saturation/

for i=1:size(X,1)-1
	X(i,:) = X(i,:)./X(end,:);
end
X(end,:) = 1;

author: IVRL

repo: Kernel-Modeling-Super-Resolution


text:
% LOGDET_KRON_LDLCHOL - Computes the log-determinant of a matrix using its
% Cholesky factor from SuiteSparse.
%
% Y = LOGDET_LDLCHOL(LD1,LD2)
%
% Computes the log-determinant of the matrix X whose Cholesky factor is
% the triangular matrix U. U can be lower or upper triangular.
%
% More generally, X could be a product of two instances of U and/or U'. For
% instance, X=U'*U' or X=U'*U.

% Last modified 2010-11-09
% Copyright (c) Jaakko Luttinen

function y = logdet_ldlchol(LD1,LD2)

N1 = size(LD1,1);
N2 = size(LD2,1);

y = N2*logdet_tri(LD1) + N1*logdet_tri(LD2);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [f] = spm_fx_poly(x,v,P)
% Normal (bilinear) form equation of motion
% FORMAT [f] = spm_fx_poly(x,v,P)
% x      - state vector
% v      - exogenous cause
% P      - free parameters
%
% f      - dx/dt
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_poly.m 3878 2010-05-07 19:53:54Z karl $

% compute Jacobian from blinear terms
%--------------------------------------------------------------------------
x     = spm_vec(x);
J     = P.A;
for i = 1:length(P.B)
    J = J + P.B{i}*x(i);
end
for i = 1:length(P.C)
    J = J + P.C{i}*v(i);
end
f     = J*x;

author: spm


text:
function [v,beta]=HouseholderVec(x,forceSign)
%%HOUSEHOLDERVEC Given an mX1 column vector x, compute mX1 column
%          vector v and scalar beta such that P=eye(m,m)-beta*v*v' is
%          orthogonal and P*x=(+/-)norm(x)*e1 for real x, where e1 is a
%          vector with a 1 in the first entry and zeros elsewhere. If
%          forceSign is true, then it will be such that the + sign is
%          always chosen. For complex x,
%          P*x=(+/-)exp(1j*angle(x(1)))*norm(x)*e1, whereby the sign in
%          front can vary. The v vector is chosen such that v(1)=1. If a
%          scalar value is provided, then v=1, beta=0 are returned.
%
%INPUTS: x An mX1 vector, real or complex.
% forceSign If true and x is real, then P*x always equals +norm(x)*e1.
%          Otherwise, it can be + or -. Often, assumptions on recovering
%          matrices from factored form implicitly assume that
%          forceSign=false. The default if omitted or an empty matrix is
%          passed is false.
%
%OUTPUTS: v An mX1 vector such that P=eye(m)-beta*v*v' is an orthogonal
%           matrix and P*x having the only nonzero element be the first
%           one. If a scalar input is given, then v=0 is returned.
%      beta The real scalar beta as mentioned with v. beta= 2/(v'*v). If a
%           scalar input is given, then beta=0 is returned.
%
%For real vectors x, the Householder vector algorithm of Section 5.1.1 of
%[1] is used. For complex vectors, the Householder algorithm of Section
%5.1.13 is used.
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(forceSign))
    forceSign=false;
end

m=length(x);

if(isempty(x))
   v=[];
   beta=[];
   return;
end

if(m==1)
   v=1;
   beta=0;
   return;
end

if(isreal(x))%The real case from Chapter 5.1.3 in [1].
    sigma=x(2:m)'*x(2:m);
    v=[1;x(2:m)];
    if(sigma==0&&x(1)>=0)
        if(forceSign)
            beta=0;
        else
            beta=2;
        end

        return;
    elseif(sigma==0&&x(1)<0)
        beta=2;%Sign corrected from that in the book.
        return;
    else
        mu=sqrt(x(1).*x(1)+sigma);
        %The book has <=; we are using < as it does not change things and
        %for some data types, such as Intervals, it is easier to overload <
        %than <=.
        if(x(1)<0)
            v(1)=x(1)-mu;
        else
            v(1)=-sigma/(x(1)+mu);
        end
        beta=2*v(1).*v(1)./(sigma+v(1).*v(1));
        v=v./v(1);
    end
else%The complex case from Chapter 5.1.13
    angVal=x(1)/abs(x(1));
    %Deal with NaNs due to x(1)=0.
    if(~isfinite(angVal))
        angVal=1;
    end

    x2Norm=norm(x,2);

    v=zeros(m,1);

    v(2:end)=x(2:end);
    v1=x(1)+angVal*x2Norm;
    v2=x(1)-angVal*x2Norm;

    %Of v1 and v2, choose the one with the largest norm.
    if(abs(v1)>abs(v2))
        v(1)=v1;
    else
        v(1)=v2;
    end

    v=v/v(1);
    beta=2/(v'*v);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Save and close the Simulink model through Matlab command line.
% It also closes the associate static GUI

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('[closeModel]: closing the Simulink model...')

% save and close the Simulink model
save_system('jointsControl.mdl');
close_system('jointsControl.mdl');

% close all figures
close all

% remove paths (optional)
rmpath('../../library/matlab-gui');
rmpath('./src-static-gui');

disp('[closeModel]: done.')

author: robotology

repo: whole-body-controllers


text:
function CLim= visutil_unifyCLim(h, CLim)
%VISUTIL_UNIFYCLIM - Unify the CLim in all axes pointed to by handles h or
% in the current figure if no handles are given
%
%Synopsis:
% CLIM= visutil_unifyCLim
% CLIM= visutil_unifyCLim(H)
% CLIM= visutil_unifyCLim(H, CLIM)
%
%Input:
% H: List of axis handles
% CLIM: CLim value to set, otherwise the color limit range among all axes
%       taken

% blanker@cs.tu-berlin.de, 09/2009

if ~exist('h','var') || isempty(h),
  h= get(gcf, 'children');
end
hax= findobj(h, 'Type','axes', 'Tag','');
if ~exist('CLim','var') || isempty(CLim)
    for hi= 1:length(hax),
        cl(hi,:)= get(hax(hi), 'CLim');
    end
    CLim= [min(cl(:,1)) max(cl(:,2))];
end
set(hax, 'CLim', CLim);

author: bbci


text:
function job = calcParent2Child(job, varargin)
% compute optimal parent to child orientation relationship
%
% The function |calcParent2Child| uses the parent to child orientation
% relationship stored in |job.p2c| as a starting point for an iterative
% process to find a parent to child orientation
% relationship that best possible fits to the child to child
% misorientations in the measured grain data |job.grains|.
%
% Syntax
%
%   % find optimal parent to child orientation relationship
%   job.calcParent2Child
%
%   % display distribtion of the misfit
%   histogram(job.calcGBFit ./ degree)
%
% Input
%  job - @parentGrainReconstructor
%
% Output
%  job.p2c - fitted parent to child orientation relationship
%
% Options
%  c2c - consider only child to child misorientations
%  p2c - consider only parent to child misorientations
%  peakFitting - use peak fitting algorithm
%  quantile  - consider only misorientation within this quantile to the current p2c guess (default 0.9)
%  threshold - only consider misorientations that are within this threshold of the current parent to child OR guess
%
% References
%
% * Tuomo Nyyssönen, <https://www.researchgate.net/deref/http%3A%2F%2Fdx.doi.org%2F10.1007%2Fs11661-018-4904-9?_sg%5B0%5D=gRJGzFvY4PyFk-FFoOIj2jDqqumCsy3e8TU6qDnJoVtZaeUoXjzpsGmpe3TDKsNukQYQX9AtKGniFzbdpymYvzYwhg.5jfOl5Ohgg7pW_6yACRXN3QiR-oTn8UsxZjTbJoS_XqwSaaB7r8NgifJyjSES2iXP6iOVx57sy8HC4q2XyZZaA
% Crystallography, Morphology, and Martensite Transformation of Prior
% Austenite in Intercritically Annealed High-Aluminum Steel>
%
% See also
% calcParent2Child
%

noOpt = ~check_option(varargin,{'p2c','c2c'});
threshold = get_option(varargin,'threshold',5*degree);
quant = get_option(varargin,'quantile',0.5);

if ~isempty(job.p2c), p2c0 = job.p2c; end
p2c0 = getClass(varargin,'orientation',p2c0);

% get p2c from parent 2 child OR
if nnz(job.ebsdPrior.phaseId==job.parentPhaseId) > 0.01 * length(job.ebsdPrior) ...
    && (noOpt || check_option(varargin,'p2c'))

  p2cPairs = neighbors(job.grains(job.csParent),job.grains(job.csChild));

  p2c = inv(job.grains(p2cPairs(:,2)).meanOrientation) .* ...
    job.grains(p2cPairs(:,1)).meanOrientation;

  for k = 1:3
    omega = angle(p2c,p2c0);
    ind = omega < min(threshold, quantile(omega, quant));
    p2c0 = mean(p2c(ind));
  end

  if check_option(varargin,'peakFitting')
    mdf = calcDensity(p2c(ind),'halfwidth',1*degree);
    p2c0 = steepestDescent(mdf,p2c0);
  end

  p2cData = p2c0;

else % some default parent2child orientation relationship

  p2cData = [];

end

% consider also child to child
if (noOpt && angle(p2c0)>5*degree) || check_option(varargin,'c2c')

  % get neighbouring grain pairs
  [c2cPairs, oriChild] = getC2CPairs(job, varargin{:});

  % compute c2c misorientation
  mori = inv(oriChild(:,1)) .* oriChild(:,2);

  % ignore pairs with misorientation angle smaller then 5 degree
  mori(mori.angle < 5 * degree) = [];

  if ~isempty(job.p2c), p2c0 = job.p2c; end
  p2c0 = getClass(varargin,'orientation',p2c0);

  % compute an optimal parent to child orientation relationship
  if check_option(varargin,'v3')

   p2c = fitP2C(mori,p2c0);

  elseif check_option(varargin,'v2')

    p2c = calcParent2Child2(mori,p2c0);

  else

    p2c = calcParent2Child(mori,p2c0,varargin{:});

  end

  % combine p2c and c2c orientation relationships
  if ~isempty(p2cData) && ~check_option(varargin,'noP2C')

    p2c = mean([p2c,p2cData],'weights',[length(c2cPairs),length(p2cPairs)]);

  end

else

  p2c = p2c0;

end

% compute new variantMap
if ~isempty(job.p2c) && length(p2c.variants) == length(job.p2c.variants)
  [~,new2old] = min(angle_outer(p2c.variants,job.p2c.variants,'noSym1'),[],2);
  p2c.opt.variantMap = new2old;
end

% update p2c
job.p2c = p2c;

end

author: mtex-toolbox

repo: mtex

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


text:
function [year,month,day,hour,minute,second]=TT2Cal(Jul1,Jul2,givenDayFrac)
%%TT2CAL Convert terrestrial times (TT) given as a two-part pseudo-Julian
%        dates to a date in terms of the Gregorian calendar in years,
%        months, days, hours, minutes and seconds or in terms of the
%        Gregorian calendar in years, months, days and the fraction of a
%        day.
%
%INPUTS: Jul1, Jul2 Vectors or matrices of the two parts of Julian dates
%                   given in TT. The units of the date are days. The full
%                   date is the sum of both terms. The date is broken into
%                   two parts to provide more bits of precision. It does
%                   not matter how the date is split.
%       giveDayFrac An optional boolean variable specifying whether the
%                   output should be as years months, days and a fraction
%                   of a day. If this parameter is omitted, the default is
%                   false. That means that the output will be given as
%                   years, months, days, hours, minutes, and seconds.
%
%OUTPUTS: Regardless of the value of giveDayFrac, the first three outputs
%         are the same. They are:
%           year A vector or matrix of years in the Gregorian calendar
%                under UTC time, one for for each Julian date.
%          month A vector or matrix of months in the Gregorian calendar
%                under UTC time, one for each Julian date. 1<=month<=12
%            day A vector or matrix of days in the Gregorian calendar
%                under UTC time, one for each Julian date. Days count from
%                1.
%          If giveDayFrac is omitted or is false, then three additional
%          outputs are
%           hour  A vector or matrix of hours in the Gregorian calendar
%                 under UTC time, one for each Julian date. 0<=hour<=23.
%          minute A vector or matrix of minutes in the Gregorian calendar
%                 under UTC time, one for each Julian date.
%                 0<=minute<=59.
%          second A vector or matrix of seconds in the Gregorian calendar
%                 under UTC time, one for each Julian date. This includes
%                 the possibility of a leap second on the day in question.
%          On the other hand, if giveDayFrac is true, then the one
%          additional output is
%          dayFrac (Given as the hour output) A vector or matrix of values
%                  >=0 and <1 indicating the fraction of the day elapsed,
%                  one for each Julian date.
%          and requesting more than four outputs will cause an error.
%
%This function just calls the function TT2UTC and then UTC2Cal.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(givenDayFrac))
   givenDayFrac=false;
end

[Jul1,Jul2]=TT2UTC(Jul1,Jul2);

if(givenDayFrac==false)
	[year,month,day,hour,minute,second]=UTC2Cal(Jul1,Jul2,givenDayFrac);
else
    if(nargout>4)
        error('Wrong number of outputs')
    end
    [year,month,day,hour]=UTC2Cal(Jul1,Jul2,true);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% sinusoidal sequence

 n=0:20;
 x=2*cos(1/2*n+pi/4);
 stem(n,x)
 legend('x[n]')
 grid

 figure
 y=2*cos(pi/6*n+pi/4);
 stem(n,y)
 legend('y[n]')
 grid

%sampling
figure
t=0:0.01:10;
x=cos(7*t);

Ts1=pi/7;
ts1=0:Ts1:10;
xs1=cos(7*ts1);

Ts2=pi/4;
ts2=0:Ts2:10;
xs2=cos(7*ts2);

plot(t,x,ts1,xs1,':o',ts2, xs2,':+')
legend('x(t)','x[n],T_s=\pi/7','x[n],T_s=\pi/4')

author: Sable


text:
function y = objfun_grad(x)

y = [(2/3)*x(1);x(2)];

repo: Optimization-Theory-and-Algorithm


text:
function status = SetupRawDir(dirs,anatSeriesNames,makeIFileFlag)
% Set up the directories containing Pfiles, Anat files
%
%  status = SetupRawDir([dirs,anatSeriesNames,makeIFileFlag]);
%
% This routine is useful if you've dumped the Pfiles and DICOM files (using
% sftp and get11, respectively) in a single directory.
%
% dirs: specify the directory, or a cell specifying many directories.
% Default is current working directory.
%
% anatSeriesNames: a cell specifying the names of each
% anatomical series, e.g. {'3plane' 'Inplane' 'SS'}.
% Anatomical series are identified as directories with
% numeric names like '001' or '099'. The script
% will rename the series in numeric order according to the
% name. Leave empty to not rename that series. (E.g.,
% {'Inplane' '' '' 'SS'} ignores the 2 middle series.)
% enter the number 0 to get a default:
% {'Localizer' 'Inplane' 'SS'}
% leave blank and it will prompt for each series #.
%
% makeIFileFlag: if 1, use dcm2gen to convert DICOM images
% to genesis I-files for each anatomical subdir.
%
% Analogous to the old script simptrans3.sh [now defunct],
% set up the Raw directory for a session/sessions, with
% the appropriate subdirectories.
%
% ras 12/01/04.
status = 0;

if ~exist('dirs','var') || isempty(dirs)
    dirs = pwd;
end

if ~exist('makeIFileFlag','var') || isempty(makeIFileFlag)
    makeIFileFlag = 0;
end

if ~exist('anatSeriesNames','var') || isempty(anatSeriesNames)
    anatSeriesNames = []; % flag to check for each dir
elseif isequal(anatSeriesNames,0)
    anatSeriesNames = {'Localizer' 'Inplane' 'SS'};
end

if iscell(dirs)
    for i = 1:length(dirs)
        SetupRawDir(dirs{i},anatSeriesNames,makeIFileFlag);
    end
    return
end

callingDir = pwd;

cd(dirs);

% make the Raw dir and subdirs
if ~exist('Raw','dir') || ~exist('Raw/Anatomy','dir') ...
    || ~exist('Raw/Pfiles','dir')
    fprintf('Making Raw directory and subdirs...\n');
%     unix('mkdir Raw Raw/Anatomy Raw/Pfiles');
    mkdir('Raw');
    mkdir('Raw','Anatomy');
    mkdir('Raw','Pfiles');
end

% get list of anatomical series dirs
anatSeries = dir('0*');
if ~isempty(anatSeries)
    isdir = [anatSeries.isdir];
    names = {anatSeries.name};
    for i = 1:length(names)
        nums(i) = str2num(names{i});
    end
    isnum = ismember(nums,(1:99));
    nums = nums(isnum==1);
    names = names(isnum==1);
    anatSeries = anatSeries(isdir==1 & isnum==1);
end

% move anat series to Raw/Anatomy/
fprintf('Moving anatomical series dirs to Raw/Anatomy...\n');
for i = 1:length(anatSeries);
%     cmd = sprintf('mv %s Raw/Anatomy/',anatSeries(i).name);
%     unix(cmd);
    movefile(anatSeries(i).name, 'Raw/Anatomy');
end

% make sure we have a list of names for anat series
if isempty(anatSeriesNames)
    for i = 1:length(anatSeries);
        msg = sprintf('Enter name for series %i [enter to leave as is]: ',nums(i));
        anatSeriesNames{i} = input(msg,'s');
    end
elseif isequal(anatSeriesNames,0)
    anatSeriesNames = {'Localizer' 'Inplane' 'SS'};
    anatSeriesNames(end+1:length(anatSeries)) = '';
end

% rename anat series
for i = 1:length(anatSeries)
    if ~isempty(anatSeriesNames{i})
%         cmd = sprintf('mv Raw/Anatomy/%s Raw/Anatomy/%s',...
%                        anatSeries(i).name,anatSeriesNames{i});
%         unix(cmd);
        cd Raw/Anatomy
        movefile(anatSeries(i).name, anatSeriesNames{i});
        cd ../..
    end
end

% get list of Pfiles, move to Raw/Pfiles
pfiles = dir('P*.7*');
pnames = {pfiles.name};
fprintf('Moving Pfiles and Pmags to Raw/Pfiles...\n');
for i = 1:length(pnames)
%     cmd = sprintf('mv %s Raw/Pfiles/',pnames{i});
%     unix(cmd);
    movefile(pnames{i}, 'Raw/Pfiles/');
end

% move .img files (from rtfmri) to Raw/Pfiles/imgFiles
imgfiles = dir('P*.img');
if ~isempty(imgfiles)
    fprintf('Making Raw/Pfiles/imgFiles and moving .img files there...\n');
%     unix('mkdir Raw/Pfiles/imgFiles');
    mkdir('Raw/Pfiles', 'imgFiles');
    for i = 1:length(imgfiles)
%         cmd = sprintf('mv %s Raw/Pfiles/imgFiles/',imgfiles(i).name);
%         unix(cmd);
        movefile(imgfiles(i).name, 'Raw/Pfiles/imgFiles/');
    end
end

% get list of E-file headers, compare to pfile names
% (we want to store 'junk' e-files that lack corresponding
% p-files in Raw/Pfiles/extraEfiles
efiles = dir('E*.7');
if ~isempty(efiles)
    enames = {efiles.name};
    for i = 1:length(enames)
        enums(i) = str2num(enames{i}(12:16));
    end
    for i = 1:length(pnames)
        pnums(i) = str2num(pnames{i}(2:6));
    end
    pnums = unique(pnums);
    ind = ismember(enums,pnums);
    goodEfiles = find(ind==1);
    extraEfiles  = find(ind==0);

    % move 'good' efiles to Raw/Pfiles
    for i = goodEfiles
%         cmd = sprintf('mv %s Raw/Pfiles/',enames{i});
%         unix(cmd);
        movefile(enames{i}, 'Raw/Pfiles');
    end

    % move 'extra' efiles to Raw/Pfiles/extraEfiles
    for i = extraEfiles
          if ~exist('Raw/Pfiles/extraEfiles', 'dir')
              mkdir('Raw/Pfiles/', 'extraEfiles');
          end
%         cmd = sprintf('mv %s Raw/Pfiles/extraEfiles/',enames{i});
%         unix(cmd);
          movefile(enames{i}, 'Raw/Pfiles/extraEfiles');
    end
end

% read the first I-file found to get exam #
fprintf('Trying to get exam #...');
try
    tmp = dir('Raw/Anatomy');
    pth = fullfile('Raw','Anatomy',tmp(3).name);
    tmp = dir(fullfile(pth,'I*'));
    [img hdr] = ReadMRImage(fullfile(pth,tmp(1).name));
%     cmd = sprintf('echo %i >Raw/examNum.txt',...
%                    hdr.exam(1).ex_no);
%     unix(cmd);
    fid = fopen('Raw/examNum.txt', 'w');
    fprintf(fid, '%s', num2str(hdr.exam(1).ex_no));
    fclose(fid);
    fprintf('succeeded.\n');
catch
    fprintf('Couldn''t get exam #\n');
end

% convert DICOM to Genesis I-file, if selected
if makeIFileFlag==1 && isunix
    dcm2genPath = fullfile(RAID,'dataTools','dcm2gen');
    cd('Raw/Anatomy');
    subdirs = dir;
    subdirs = subdirs([subdirs.isdir]==1);
    for i = 1:length(subdirs)
        cmd = sprintf('%s %s/*',dcm2genPath,subdirs(i).name);
        unix(cmd);
    end
end

cd(callingDir);

status = 1;

return

author: vistalab

repo: vistasoft


text:
% Linear Diophantine Equation
%Given three integers a, b, c representing a linear equation of the form : ax + by = c.
%To find if the equation has a solution such that x and y are both integral values.
%Used in programming to find the exact solution exists or not.

function retval = linear_diophantine_eqn (a,b,c)
   if mod(c,gcd(a,b))==0
     retval=true;            % 1 represent yes it exist
   else
     retval= false;          %  0 represent no it doesn't exist

end

author: TheAlgorithms

repo: MATLAB-Octave