text:
function YLim= visutil_selectYLim(h, varargin)

props = {'Policy',          'auto'          '!CHAR(auto tightest tight)';
         'TightenBorder'     0.03        	'DOUBLE';
         'Symmetrize'       0               '!BOOL';
         'SetLim'           1               '!BOOL';
         };

if nargin==0,
  YLim= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);
misc_checkType(h,'!GRAPHICS');

switch(opt.Policy),
 case 'auto',
  YLim= get(h, 'YLim');
 case 'tightest',
  visutil_backaxes(h);
  axis('tight');
  YLim= get(h, 'YLim');
 case 'tight',
  visutil_backaxes(h);
  axis('tight');
  yl= get(h, 'YLim');
  % add border not to make it too tight:
  yl= yl + [-1 1]*opt.TightenBorder*diff(yl);
  % determine nicer limits
  dig= floor(log10(diff(yl)));
  if diff(yl)>1,
    dig= max(1, dig);
  end
  YLim= [util_trunc(yl(1),-dig+1,'floor') util_trunc(yl(2),-dig+1,'ceil')];
end

if opt.Symmetrize,
  ma= max(abs(YLim));
  YLim= [-ma ma];
end

if opt.SetLim,
  set(h, 'YLim',YLim);
end

if nargout==0,
  clear YLim;
end

author: bbci


text:
%%*******************************************************
%% fapproblem:
%%
%% (primal prob.) min  Tr C*X
%%                s.t. diag(X) = e
%%                    -Eij * X = 2/(k-1)  if (i,j) in  U
%%                    -Eij * X <= 2/(k-1) if (i,j) in GU
%%  where Eij = ei*ej' + ej*ei'
%%  b = [e; 2/(k-1)*ones(mU+mGU,1)]
%%  C = (1/2)*Diag(We) - (k-1)/(2k) *L(G)
%%-------------------------------------------------------
%% [blk,Avec,C,b] = fapproblme(fname);
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*******************************************************

 function [blk,Avec,C,b,kparm] = fapread(fname)

%%
%% read fap data
%%
   if exist(fname)
      fid = fopen(fname,'r');
   elseif exist([fname,'.dat']);
      fid = fopen([fname,'.dat'],'r');
   else
      fprintf('** Problem not found. \n');
      blk = []; Avec = []; C = []; b = []; kparm = [];
      return;
   end
   [tmpr,count] = fscanf(fid,'%c');
   datavec = sscanf(tmpr,'%f'); clear tmpr;
   n = datavec(1);
   numedges = datavec(2);
   kparm = datavec(3);
   datavec = datavec(4:length(datavec));
   len = length(datavec);
   if (len ~= 3*numedges)
      error(' fapread: numedges and data do not match.');
   end
   I = datavec(1:3:len);
   J = datavec(2:3:len);
   w = datavec(3:3:len);

   idxU = find(w==1000);
   IU = I(idxU); JU = J(idxU); wU = w(idxU);
   U = spconvert([IU JU wU; n n 0]);
   U = U + U';

   idx2 = find(w~=1000);
   I2 = I(idx2); J2 = J(idx2); w2 = w(idx2);
   GU = spconvert([I2 J2 w2; n n 0]);
   GU = GU + GU';

   fclose(fid);
%%
%% blk, Avec, C, b
%%
    n = length(U);
    mU = nnz(triu(U,1));
    mGU = nnz(triu(GU,1));
    m = mGU + mU + n;

    b = [ones(n,1); (2/(kparm-1))*ones(mU+mGU,1)];
    LG = diag(GU*ones(n,1))-GU;
    C{1,1} = 0.5*(diag(GU*ones(n,1))) - ((kparm-1)/(2*kparm))*LG;
    C{2,1} = zeros(mGU,1);

    blk{1,1} = 's';  blk{1,2} = n;
    blk{2,1} = 'l';  blk{2,2} = mGU;
%%
%%
%%
    r2 = sqrt(2);
    I = zeros(m,1); J = zeros(m,1); w = zeros(m,1);
    cnt = 0;
    e = [1:n]';
    I(1:n) = e.*(e+1)/2;
    J(1:n) = e;
    w = ones(n,1);
    cnt = cnt+n;
    for i = 1:n
        idx = find(U(i,i+1:n));
        idx = idx+i;      %% adjust index.
        len = length(idx);
        I(cnt+[1:len]) = i + idx.*(idx-1)/2;
        J(cnt+[1:len]) = cnt+[1:len]';
        w(cnt+[1:len]) = -r2*ones(len,1);
        cnt = cnt + len;
    end
    for i = 1:n
        idx = find(GU(i,i+1:n));
        idx = idx+i;      %% adjust index.
        len = length(idx);
        I(cnt+[1:len]) = i + idx.*(idx-1)/2;
        J(cnt+[1:len]) = cnt+[1:len]';
        w(cnt+[1:len]) = -r2*ones(len,1);
        cnt = cnt + len;
    end
    Avec{1,1} = spconvert([I,J,w; n*(n+1)/2, m, 0]);
%%
    I = [1:mGU]'; J = [m-mGU+1:m]';  w = ones(mGU,1);
    Avec{2,1} = spconvert([I, J, w; mGU, m, 0]);
%%
%%*******************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function siftroot = hlp_getSiftRoot
% get the root path for SIFT

siftroot = fileparts(utl_whichfile('StartSIFT'));

repo: REST


text:
% PENDING
%
%
%
%
%
%

function [ERP serror] = importerpss(filename, format, transpose)

serror = 0;
if nargin<2
      error('ERPLAB says: error, filename and pathname are needed as inputs.')
end

%%%taken from asc2erp.m%%%
fullname = filename;
fid = fopen(fullname);

%used to verify whether it is an ERPSS file
C = textscan(fid, '%s', 5);
fclose(fid);

G = [C{1:end}];

disp('hola')

%if strcmpi(G{5}, 'nchans') % then it is ERPSS

      % %       #
      % %       #	binno          1
      % %       #
      % %       bindesc="Cue LSF "
      % %       condesc="Cue LSF "
      % %       npoints      613
      % %       sampleperiod      1953
      % %       presampling     199219
      % %       sums        1181
      % %       procfuncs   "avg"
      % %       arejects    819
      % %

      %
      % Getting bin descriptors
      %
      fid = fopen(fullname);
      BD  = textscan(fid, '%s', 'delimiter','\n');
      fclose(fid);

      BDesc = regexpi(BD{:}, 'bindesc="(.*)"','tokens');
      BDesc = BDesc(~cellfun(@isempty, BDesc));
      BDesc = [BDesc{:}];
      BDesc = [BDesc{:}];
      % %       C   = regexprep(C{:}, '"(.*)"','"${strrep($1,'' '',''_'')}"');
      % %
      % %       [pathproc fnameproc ext] = fileparts(fullname);
      % %
      % %       fullfile(pathproc, [fnameproc '_procerplab' ext])
      % %
      % %       fid = fopen(fullfile(pathproc, [fnameproc '_procerplab' ext]), 'w');
      % %       fprintf(fid, '%s\n', C{:});
      % %       fclose(fid);
      % %
      % %       fid = fopen([ fnameproc '_procerplab' ext]);
      % %       C   = textscan(fid, '%s');
      % %       fclose(fid);

      fid = fopen(fullname);
      C   = textscan(fid, '%s');
      fclose(fid);

      D   = [C{1:end}];
      m   = length(D);
      count=1; x=0; a=0; b=0; c=0; d=0; e=0; f=0; j=1; k=0;

      %Explicit Format (header for bins and channels)
      if format==0  % explicit
            fprintf('\nPlease wait, loading ERPset from ERPSS ascii file (explicit format)...\n');
            i=1;

            while i<=m-5
                  %x = D{i};
                  a = D{i};
                  b = D{i+1};
                  d = str2num(D{i+3});
                  f = D{i+5};
                  err = 0;

                  if (a=='#' & b=='#' & d>=0 & d<=256 & f=='#' )
                        k = i+6;
                        j = 1;

                        while (D{k} ~= '#' & err==0)
                              E(j,1,count) = str2num(D{k});
                              j = j+1;
                              k = k+1;
                              if k==m
                                    err=1;
                              end
                        end
                        count = count+1; % bin counter
                        i = k-1;
                  else
                        i = i+1;
                  end
            end

            %to get very last data point
            E(j,1, count-1) = str2num(D{m});
            %fclose(fid);

            %finding channel numbers and bin numbers
            ERP       = buildERPstruct;
            numofchan = str2num(D{6});
            numofbins = ((count-1)/numofchan)-1;

            for binnumber=1:numofbins
                  k = (binnumber*numofchan) + 1;
                  for n=1:numofchan
                        ERP.bindata(n,:,binnumber) = E(:,1,k)';
                        k = k+1;
                  end
            end
            %Need to fill ERP.bindata, numofchans, numofbins
      elseif format ==1     % implicit
            fprintf('\nPlease wait, loading ERPset from ERPSS ascii file (implicit format, ');

            %rows are channels, columns are samples
            if transpose ==0
                  fprintf('untransposed)...\n')
                  i=1;

                  while i<=m-5
                        %x = D{i};
                        a = D{i};
                        b = D{i+2};
                        d = str2num(D{i+3});
                        err = 0;

                        if (strcmp(a,'procfuncs') & strcmp(b,'arejects') & d>=0) % JLC: "d" might not necessarily be zero
                              k=i+4;
                              j=1;

                              while (D{k} ~= '#' & err==0)
                                    E(j,1,count) = str2num(D{k});
                                    j = j+1;
                                    k = k+1;
                                    if k==m
                                          err=1;
                                    end
                              end
                              count = count+1;
                              i = k-1;
                        else
                              i = i+1;
                        end
                  end
                  %to get very last data point
                  E(j,1, count-1) = str2num(D{m});
                  %fclose(fid);

                  %finding channel numbers and bin numbers
                  ERP = buildERPstruct;
                  numofchan  = str2num(D{6});
                  numofbins  = count-1;
                  sigperchan = j/numofchan; %number of signals per channel
                  k=1;
                  m=1;

                  for binnumber = 1:numofbins
                        for n = 1:numofchan
                              for j=1:sigperchan
                                    x = E(m,1,k);
                                    ERP.bindata(n,j,binnumber) = x;
                                    m = m+1;
                              end
                        end
                        k = k+1;
                        m=1;
                  end
                  %columns are channels, rows are samples
            else
                  fprintf('transposed)...\n')
                  i=1;

                  while i<=m-5
                        %x = D{i};
                        a = D{i};
                        b = D{i+2};
                        d = str2num(D{i+3});
                        err =0;

                        if (strcmp(a,'procfuncs') & strcmp(b,'arejects') & d>=0) %JLC: "d" might not necessarily be zero
                              k = i+4;
                              j=1;

                              while (D{k} ~= '#' & err==0)
                                    E(j,1,count) = str2num(D{k});
                                    j = j+1;
                                    k = k+1;
                                    if k==m
                                          err=1;
                                    end
                              end
                              count = count+1;
                              i = k-1;
                        else
                              i = i+1;
                        end
                  end
                  %to get very last data point
                  E(j,1, count-1) = str2num(D{m});
                  %fclose(fid);

                  %finding channel numbers and bin numbers
                  ERP = buildERPstruct;
                  numofchan  = str2num(D{6});
                  numofbins  = (count-1);
                  sigperchan = j/numofchan; %number of signals per channel
                  k=1;

                  for binnumber=1:numofbins
                        for n=1:numofchan
                              m=n;
                              for j=1:sigperchan
                                    x = E(m,1,k);
                                    ERP.bindata(n,j,binnumber) = x;
                                    m = m+numofchan;
                              end
                        end
                        k = k+1;
                        m=1;
                  end
            end
      end

      ERP.nbin     = numofbins; %number of bins
      ERP.nchan    = numofchan; %number of channels
      ERP.binerror = [];
      ERP.isfilt   = 0;
      [a numberofpoints] = size(ERP.bindata(1,:,1));
      ERP.pnts  = numberofpoints;

      % ################# NO contrast against OPTIONS #################### JLC
      %ERP.xmin  = xlim(1)/1000;
      %ERP.xmax  = xlim(2)/1000;

      %error if fs and estimated fs are not within 10% of each other
      %sdistance = (xlim(2)-xlim(1))/numberofpoints;  %distance between points
      %est_fs=1/sdistance*1000;
      %srdiff = abs((fs-est_fs)/est_fs);
      %if srdiff>0 && srdiff>=0.1
      %        serror=3;
      %end
      %###########################  JLC

      %now if the srate from the erpss is offf another error
      %now take the value of ERP.srate= the one from ERPSS

      for i=1:500 %because the sample period will be located somewhere in between here
            if strcmp(D{i},'sampleperiod')
                  sam_per  = str2num(D{i+1});
                  fs_ERPSS = round((1/sam_per)*10^6);
                  break % JLC
            end
      end

      % #######  JLC  ##############
      for i=1:500 %because the sample period will be located somewhere in between here
            if strcmp(D{i},'presampling')
                  pre_sam = str2num(D{i+1});
                  pre_sam = -pre_sam/1000; % usec to msec
                  break % JLC
            end
      end

      %Create ERP structure

      % JLC Remainder:  presampling     200000
      %
      % @Eric : ERPSS text file has ALL info you need to load it. It is NOT necessary to contrast fs, xmin, xmax, pnts, etc
      % against what the user entered in OPTIONS at the IMPORT GUI. OPTIONS must work only for universal text.
      %
      % ####################

      %####### no fs test for ERPSS #######  JLC
      %if abs(((fs-fs_ERPSS)/fs_ERPSS)*100)>10
      %        serror=3;
      %end
      %###################################

      %ERP.srate     = fs;
      %ERP.srate=(1/sdistance)*1000;
      ERP.srate = fs_ERPSS;  % JLC

      %times= round(linspace(xlim(1), xlim(2), numberofpoints));

      % JLC  ######
      times = pre_sam:(1000/fs_ERPSS):numberofpoints*(1000/fs_ERPSS)+(pre_sam)-1;
      ERP.times = times;
      ERP.xmin  = min(times)/1000;
      ERP.xmax  = max(times)/1000;
      % JLC ###########

      %channel locations
      fid = fopen(fullname);
      C   = textscan(fid, '%d %s', 256, 'headerLines', 6);
      fclose(fid);

      channelnames = cellstr(C{2});
      [a b] = size(channelnames);

      for i=1:a
            k = channelnames{i};
            j = regexp(k,'"(\w*)"', 'tokens');
            channelnames{i} = char(j{:});
      end

      b = int2str(C{1});
      channelnumbers = cellstr(b);

      for i=1:a
            x = channelnumbers{i};
            tempchannum = str2num(x);
            tempchannum = tempchannum+1;
            channelnumbers{i} = int2str(tempchannum);
      end
      if isfield(ERP, 'chanlocs')
            ERP=rmfield(ERP, 'chanlocs');
      end
      %channelnames
      for i=1:a
            ERP.chanlocs(1,i) = struct('labels', channelnames{i});
      end
      %ERP.erpname
      g = D{1};
      g = strtrim(strrep(g,' ','_'));
      g = strrep(g,'-','_');
      %g(ismember_bc2(double(g),[45])) = ''; %removes '-'   %%If name starts with numbers, will not work

      a = regexp(g,'"(.*)"', 'tokens');
      a = strrep(a{:},'"','');
      %     if double(b(1))==48:57
      %         strcat('A_', b);
      %     end

      ERP.erpname = char(a{:});    %char(a{:}) converts cell to character

% %       %bindescriptor name
% %       if format ==0 % explicit
% %             i=1; j=1; k=1;
% %
% %             while i<=(m-numberofpoints)
% %                   if strcmp(D{i}, 'Bin') & strcmp(D{i+1}, 'desc:')
% %                         i=i+2;
% %
% %                         while strcmp(D{i},'#') ~= 1
% %                               binnames{j,k} = [blanks(1) D{i} ];
% %                               k = k+1;
% %                               i = i+1;
% %                         end
% %                         j = j+1;
% %                         k=1;
% %                   else
% %                         i = i+1;
% %                   end
% %             end
% %             for i=1:numofbins
% %                   ERP.bindescr{1,i} = strtrim(strcat(binnames{i+1,:}));
% %             end
% %       else % implicit
% %             i=1;j=1;
% %             m = length(D);
% %
% %             while i<= m %JLC
% %                   if D{i}=='#' & strcmp(D{i+1}, 'binno') & D{i+3}=='#'
% %                         i = i+4;
% %                         h = D{i}
% %                         %h = strtrim(strrep(h,' ','_'));
% %                         h = strrep(h,'-','_')
% %                         %h(ismember_bc2(double(h),[45])) = '';   %removes '-'
% %                         %h
% %                         b = regexp(h,'bindesc=(.*)', 'tokens')  %if there is a space between binname and last ", will not work
% %                         b = strrep(b{:},'"','')
% %
% %                         ERP.bindescr{1,j} = char(b{:});
% %                         j = j+1;
% %                   else
% %                         i = i+1;
% %                   end
% %             end
% %       end
      ERP.bindescr = BDesc;
      ERP = sorterpstruct(ERP);

      % Implicit format (bin header but no channel header)

% else %if not ERPSS file
%       msgboxText =  ['Something went wrong. '...
%             'Please, verify the file format.'...
%             ' For ERPSS text file, please check the organization of the data values,'...
%             ' as well as channel headers between channels.'];
%       %title = 'ERPLAB: pop_importerpss() error:';
%       error('Error:importerpss', msgboxText);
%       return
% end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%----------------------------------------------------------------------
%Find the largest eigenvalue
%----------------------------------------------------------------------

function [ opts ] = find_max_eigenv ( matrix )

[v,s]     =       eig(matrix);
sort_s    =       sort(diag(s),'descend');
opts      =       sort_s(1);

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function out = get(f, prop)
%GET   GET method for the DELTAFUN class.
%   P = GET(F, PROP) returns the property P specified in the string PROP from
%   the DELTAFUN object F.  The string PROP may be the name of a DELTAFUN
%   property (see the DELTAFUN classdef file for a list) or any of the
%   following strings:
%       'DELTAS'       - Return the locations and magnitude of delta functions
%                        as a single matrix with the first row corresponding to
%                        the locations of the delta functions, and the
%                        magnitudes appended below the first row.
%   If PROP is a string other than those specified above, GET(F, PROP) returns
%   the result of GET(F.FUNPART, PROP).
%
% See also DELTAFUN.

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

switch prop
    case fieldnames(f)
        out = f.(prop);
    case 'deltas'
        out = [f.deltaLoc ; f.deltaMag];
    otherwise
        out = get(f.funPart, prop);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [bestEpsilon bestF1] = selectThreshold(yval, pval)
%SELECTTHRESHOLD Find the best threshold (epsilon) to use for selecting
%outliers
%   [bestEpsilon bestF1] = SELECTTHRESHOLD(yval, pval) finds the best
%   threshold to use for selecting outliers based on the results from a
%   validation set (pval) and the ground truth (yval).
%

bestEpsilon = 0;
bestF1 = 0;
F1 = 0;

stepsize = (max(pval) - min(pval)) / 1000;
for epsilon = min(pval):stepsize:max(pval)

    % ====================== YOUR CODE HERE ======================
    % Instructions: Compute the F1 score of choosing epsilon as the
    %               threshold and place the value in F1. The code at the
    %               end of the loop will compare the F1 score for this
    %               choice of epsilon and set it to be the best epsilon if
    %               it is better than the current choice of epsilon.
    %
    % Note: You can use predictions = (pval < epsilon) to get a binary vector
    %       of 0's and 1's of the outlier predictions

    tp = sum( ((pval<epsilon)==1) & (yval==1) );
    fp = sum( ((pval<epsilon)==1) & (yval==0) );
    fn = sum( ((pval<epsilon)==0) & (yval==1) );
    prec = tp / (tp+fp);
    rec  = tp / (tp+fn);
    F1 = 2 * prec * rec / (prec+rec);

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

    if F1 > bestF1
       bestF1 = F1;
       bestEpsilon = epsilon;
    end
end

end

author: zlotus


text:
% PURPOSE: demo of chowlin_co()
%          Temporal disaggregation with indicators.

% 		  Chow-Lin method, rho derived via Cochrane-Orcutt
%-------------------------------------------------------------
% USAGE: chowlin_co_d
%-------------------------------------------------------------

close all; clear all; clc;

% Low-frequency data: Bournay-Laroque (1979)
% Y: Output. Textile industries at 1956 prices. Unit: FF.
% Sample: 1949-1959

Y = [18664
    21049
    22472
    21831
    23011
    24045
    24621
    26037
    28195
    27702
    27273 ];

% High-frequency data: Bournay-Laroque (1979)
% x: Index of Industrial Production. Textile industries. Unit: base 1952=100
% Sample: 1949.1 - 1960.4
 x = [94
    93
    92
    95
    101
    104
    108
    108
    107
    111
    108
    106
    105
    99
    97
    99
    97
    101
    105
    108
    108
    110
    114
    112
    109
    107
    105
    106
    108
    113
    113
    120
    124
    127
    128
    128
    130
    124
    119
    112
    106
    117
    118
    123
    124
    125
    125
    128 ];

% ---------------------------------------------

% Inputs for td library

% Type of aggregation
ta=1;

% Frequency conversion
sc=4;

% Intercept
opC = -1;

% Name of ASCII file for output
file_sal='td.sal';

% Calling the function: output is loaded in a structure called res
res=chowlin_co(Y,x,ta,sc,opC);

% Calling printing function
tdprint(res,file_sal);

edit td.sal;

% Calling graph function
tdplot(res);

author: Sable

repo: mcbench-benchmarks


text:
function t_pf_ac(quiet)
%T_PF_AC  Tests for AC power flow solvers.

%   MATPOWER
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

%%  alg         name                            check       opts
cfg = {
    {'NR',      'Newton (default, power-polar)',[]          []  },
    {'NR-SP',   'Newton (power-polar)',         [],         []  },
    {'NR-SC',   'Newton (power-cartesian)',     [],         []  },
    {'NR-SH',   'Newton (power-hybrid)',        [],         []  },
    {'NR-IP',   'Newton (current-polar)',       [],         []  },
    {'NR-IC',   'Newton (current-cartesian)',   [],         []  },
    {'NR-IH',   'Newton (current-hybrid)',      [],         []  },
    {'FDXB',    'Fast Decoupled (XB)',          [],         []  },
    {'FDBX',    'Fast Decoupled (BX)',          [],         []  },
    {'GS',      'Gauss-Seidel',                 [],         []  },
    {'ZG',      'Implicit Z-bus Gauss',         [],         []  },
};
if have_feature('mp_core')
    cfg{end+1} = {'FSOLVE',  'fsolve (power-polar)',         'fsolve',   []  };
    cfg{end+1} = {'FSOLVE',  'fsolve (power-cartesian)',     'fsolve',   {'pf.v_cartesian', 1}  };
    cfg{end+1} = {'FSOLVE',  'fsolve (current-polar)',       'fsolve',   {'pf.current_balance', 1, 'pf.tol', 1e-10}  };
    cfg{end+1} = {'FSOLVE',  'fsolve (current-cartesian)',   'fsolve',   {'pf.v_cartesian', 1, 'pf.current_balance', 1, 'pf.tol', 1e-10}  };
%     cfg{end+1} = {'FSOLVE',  'fsolve (power-polar)',         'fsolve',   struct('Algorithm', 'trust-region-dogleg')  },
%     cfg{end+1} = {'FSOLVE',  'fsolve (power-polar)',         'fsolve',   struct('Algorithm', 'trust-region-reflective')  },
%     cfg{end+1} = {'FSOLVE',  'fsolve (power-polar)',         'fsolve',   struct('Algorithm', 'levenberg-marquardt', 'TolX', 1e-11) },
end
% %%  alg         name                            check       opts
% cfg = {
%     {'NR',      'Newton (default, power-polar)',[]          []  },
% };

ntests = 48;
t_begin(length(cfg)*ntests + 6, quiet);

casefile = 't_case9_pf';
if quiet
    verbose = 0;
else
    verbose = 1;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end
mpopt0 = mpoption('out.all', 0, 'pf.tol', 1e-9, 'verbose', 0);
mpopt0 = mpoption(mpopt0, 'verbose', verbose);

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%% network with islands
mpc0 = loadcase(casefile);
mpc0.gen(1, PG) = 60;
mpc0.gen(1, QG) = 10;
mpc0.gen(1, [PMIN PMAX QMIN QMAX PG QG]) = mpc0.gen(1, [PMIN PMAX QMIN QMAX PG QG]) / 2;
mpc0.gen = [mpc0.gen(1, :); mpc0.gen];
mpc1 = mpc0;
mpc  = mpc0;
nb = size(mpc.bus, 1);
mpc1.bus(:, BUS_I)      = mpc1.bus(:, BUS_I) + nb;
mpc1.branch(:, F_BUS)   = mpc1.branch(:, F_BUS) + nb;
mpc1.branch(:, T_BUS)   = mpc1.branch(:, T_BUS) + nb;
mpc1.gen(:, GEN_BUS)    = mpc1.gen(:, GEN_BUS) + nb;
mpc.bus         = [mpc.bus; mpc1.bus];
mpc.branch      = [mpc.branch; mpc1.branch];
mpc.gen         = [mpc.gen; mpc1.gen];
mpc1 = mpc;

%%-----  AC power flow  -----
%% get solved AC power flow case from MAT-file
load soln9_pf;      %% defines bus_soln, gen_soln, branch_soln
soln_vg = load('soln9_pf_vg');

%% run AC PF
for k = 1:length(cfg)
  if ~isempty(cfg{k}{3}) && ~have_fcn(cfg{k}{3})
    t_skip(ntests, sprintf('%s not available', cfg{k}{3}));
  else
    t = sprintf('AC PF - %s : ', cfg{k}{2});
    mpopt = mpoption(mpopt0, 'pf.alg', cfg{k}{1});
    if ~isempty(cfg{k}{4}) && iscell(cfg{k}{4})
        mpopt = mpoption(mpopt, cfg{k}{4}{:});
    end
    [baseMVA, bus, gen, branch, success, et] = runpf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(bus, bus_soln, 6, [t 'bus']);
    t_is(gen, gen_soln, 6, [t 'gen']);
    t_is(branch, branch_soln, 6, [t 'branch']);

    r = runpf(casefile, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.bus, bus_soln, 6, [t 'bus']);
    t_is(r.gen, gen_soln, 6, [t 'gen']);
    t_is(r.branch, branch_soln, 6, [t 'branch']);

    %% check when Vg ~= 1
    t = sprintf('%s - Vg ~= 1 : ', cfg{k}{1});
    mpopt = mpoption(mpopt, 'verbose', 0);
    mpc = loadcase(casefile);
    mpc.gen(:, VG) = [1.04; 1.025; 1.025];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.bus, soln_vg.bus_soln, 6, [t 'bus']);
    t_is(r.gen, soln_vg.gen_soln, 6, [t 'gen']);
    t_is(r.branch, soln_vg.branch_soln, 6, [t 'branch']);

    %% check Qg distribution, when Qmin = Qmax
    t = sprintf('%s - check Qg : ', cfg{k}{1});
    mpc = loadcase(casefile);
    mpc.gen(1, [QG QMIN QMAX]) = [10 20 20];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1, QG), 24.07, 2, [t 'single gen, Qmin = Qmax']);

    mpc.gen = [mpc.gen(1, :); mpc.gen];
    mpc.gen(1, [QMIN QMAX]) = [10 10];
    mpc.gen(2, [QMIN QMAX]) = [0 50];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [10; 14.07], 2, [t '2 gens, Qmin = Qmax for one']);

    mpc.gen(1, [QMIN QMAX]) = [10 10];
    mpc.gen(2, [QMIN QMAX]) = [-50 -50];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [42.03; -17.97], 2, [t '2 gens, Qmin = Qmax for both']);

    mpc.gen(1, [QMIN QMAX]) = [0 50];
    mpc.gen(2, [QMIN QMAX]) = [0 100];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [8.02; 16.05], 2, [t '2 gens, proportional']);

    mpc.gen(1, [QMIN QMAX]) = [-50 0];
    mpc.gen(2, [QMIN QMAX]) = [50 150];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [-50+8.02; 50+16.05], 2, [t '2 gens, proportional']);

    mpc.gen(1, [QMIN QMAX]) = [-50 Inf];
    mpc.gen(2, [QMIN QMAX]) = [50 150];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [-31.61; 55.68], 2, [t '2 gens, one infinite range']);

    mpc.gen(1, [QMIN QMAX]) = [-50 Inf];
    mpc.gen(2, [QMIN QMAX]) = [50 Inf];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [-33.12; 57.18], 2, [t '2 gens, both infinite range']);

    mpc.gen(1, [QMIN QMAX]) = [-50 Inf];
    mpc.gen(2, [QMIN QMAX]) = [-Inf 150];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [76.07; -52], 2, [t '2 gens, both infinite range']);

    mpc.gen(1, [QMIN QMAX]) = [-Inf Inf];
    mpc.gen(2, [QMIN QMAX]) = [-Inf Inf];
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(1:2, QG), [12.03; 12.03], 2, [t '2 gens, both infinite range']);

    t = sprintf('%s - reactive generation allocation : ', cfg{k}{1});
    mpc = loadcase(casefile);
    %% generator data
    %	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
    mpc.gen = [
		1	0	0	300	-300	1	100	1	250	10	0	0	0	0	0	0	0	0	0	0	0;
		2	54	0	0	-5	1	100	1	300	10	0	0	0	0	0	0	0	0	0	0	0;
		2	54	0	5	-5	1	100	1	300	10	0	0	0	0	0	0	0	0	0	0	0;
		2	55	0	25	 10	1	100	1	300	10	0	0	0	0	0	0	0	0	0	0	0;
		30	25	1	300	-300	1	100	1	270	10	0	0	0	0	0	0	0	0	0	0	0;
		30	30	2	300	-300	1	100	1	270	10	0	0	0	0	0	0	0	0	0	0	0;
		30	30	-3	300	-300	1	100	1	270	10	0	0	0	0	0	0	0	0	0	0	0;
    ];
    mpc.bus(3, BUS_TYPE) = PQ;
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.gen(2:4, QG), [-5; -5; 10] + [1; 2; 3]*1.989129794, 7, [t 'PV bus']);
    t_is(r.gen(5:7, QG), [1; 2; -3], 8, [t 'PQ bus']);

    %% network with islands
    t = sprintf('%s - network w/islands : AC PF : ', cfg{k}{1});
    mpc = mpc1;
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.bus( 1:9,  VA), bus_soln(:, VA), 7, [t 'voltage angles 1']);
    t_is(r.bus(10:18, VA), bus_soln(:, VA), 7, [t 'voltage angles 2']);
    Pg = [gen_soln(1, PG)-30; 30; gen_soln(2:3, PG)];
    t_is(r.gen(1:4, PG), Pg, 6, [t 'active power generation 1']);
    t_is(r.gen(5:8, PG), Pg, 6, [t 'active power generation 2']);

    t = sprintf('%s - all buses isolated : ', cfg{k}{1});
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = runpf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end

    %% case 14 with Q limits
    t = sprintf('%s - pf.enforce_q_lims == 0 : ', cfg{k}{1});
    mpc = loadcase('case14');
    mpc.gen(1, QMIN) = -10;
    mpc.gen(:, QMAX) = [10; 30; 29; 15; 15];
    bt0 = mpc.bus(:, BUS_TYPE);
    bt = bt0;
    mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 0);
    r = runpf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.bus(:, BUS_TYPE), bt, 12, [t 'bus type']);
    t_is(r.gen(:, QG), [-16.549300542; 43.557100134; 25.075348495; 12.730944405; 17.623451366], 6, [t 'Qg']);

    t = sprintf('%s - pf.enforce_q_lims == 1 : ', cfg{k}{1});
    mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 1);
    r = runpf(mpc, mpopt);
    bt = bt0;
    bt([1 2 3 8]) = [PQ PQ REF PQ];
    t_is(r.success, 0, 12, [t 'success = 0']);
    t_is(r.bus(:, BUS_TYPE), bt, 12, [t 'bus type']);
    t_is(r.gen(:, QG), [-10; 30; 31.608422873; 16.420423190; 15], 4, [t 'Qg']);

    t = sprintf('%s - pf.enforce_q_lims == 2 : ', cfg{k}{1});
    mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 2);
    r = runpf(mpc, mpopt);
    bt = bt0;
    bt([1 2 3 6 8]) = [REF PQ PQ PQ PQ];
    t_ok(r.success, [t 'success']);
    t_is(r.bus(:, BUS_TYPE), bt, 12, [t 'bus type']);
    t_is(r.gen(:, QG), [-6.30936644; 30; 29; 15; 15], 5, [t 'Qg']);
  end
end

t = sprintf('%s - pf.enforce_q_lims == 1 : ', 'NR-SP');
mpopt = mpoption(mpopt, 'pf.alg', 'NR-SP', 'pf.tol', 1e-9);
mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 1);
r = runpf('case14', mpopt);
t_ok(r.success, [t 'success']);
t_is(r.iterations, 5, 12, [t 'iterations']);
t_is(r.bus(1, VA), 0, 12, [t 'ref bus Va = 0']);

t = sprintf('%s - pf.enforce_q_lims == 1 : ', 'NR-IC');
mpopt = mpoption(mpopt, 'pf.alg', 'NR-SP', 'pf.tol', 1e-9);
mpopt = mpoption(mpopt, 'pf.enforce_q_lims', 1);
r = runpf('case14', mpopt);
t_ok(r.success, [t 'success']);
t_is(r.iterations, 5, 12, [t 'iterations']);
t_is(r.bus(1, VA), 0, 12, [t 'ref bus Va = 0']);

t_end;

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [svolpulse, spulse, spulse_per_vol, verbose] = tapas_physio_get_onsets_from_locs(t, VOLLOCS, LOCS, sqpar, verbose)
% creates timing vectors of found
%
%   [svolpulse, spulse, spulse_per_vol] = tapas_physio_get_onsets_from_locs(t, VOLLOCS,
%                                           LOCS, sqpar, verbose);
%
% IN
%   t           time vector of logfile (in seconds)
%   VOLLOCS     volume scan events (in samples of logfile)
%   LOCS        all slice (including volume) scan events (in samples of logfile)
%   sqpar                   - sequence timing parameters
%           .Nslices        - number of slices per volume in fMRI scan
%           .NslicesPerBeat - usually equals Nslices, unless you trigger with the heart beat
%           .TR             - repetition time in seconds
%           .Ndummies       - number of dummy volumes
%           .Nscans         - number of full volumes saved (volumes in nifti file,
%                             usually rows in your design matrix)
%           .Nprep          - number of non-dummy, volume like preparation pulses
%                             before 1st dummy scan. If set, logfile is read from beginning,
%                             otherwise volumes are counted from last detected volume in the logfile
%           .time_slice_to_slice - time between the acquisition of 2 subsequent
%                             slices; typically TR/Nslices or
%                             minTR/Nslices, if minimal temporal slice
%                             spacing was chosen
%            onset_slice    - slice whose scan onset determines the adjustment of the
%                             regressor timing to a particular slice for the whole volume
%
%   verbose             physio.verbose, See also tapas_physio_new
% OUT
%   svolpulse   vector of volume scan pulse events (in seconds from logfile start)
%   spulse      vector of scan pulse events (in seconds from logfile start)
%   spulse_per_vol  cell of slice scan events; one cell element per volume
% EXAMPLE
%   [ons_secs.svolpulse, ons_secs.spulse] = tapas_physio_get_onsets_from_locs(ons_secs.t, VOLLOCS, LOCS);
%
%   See also tapas_physio_main_create_regressors tapas_physio_read_physlogfiles_philips

% Author: Lars Kasper
% Created: 2013-02-16
%
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

Nscans          = sqpar.Nscans;
Ndummies        = sqpar.Ndummies;
NslicesPerBeat  = sqpar.NslicesPerBeat;
Nslices         = sqpar.Nslices;
do_count_from_start = isfield(sqpar, 'Nprep') && ~isempty(sqpar.Nprep);
if do_count_from_start
    Nprep = sqpar.Nprep;
end

%% remove slice scan events for last, incompletely acquired volume
if VOLLOCS(end) ~= LOCS(end-NslicesPerBeat+1)
    LOCS(LOCS>=VOLLOCS(end))=[];
    VOLLOCS(end) = [];
end

Nallvols = length(VOLLOCS);

hasWarned = false;

for v = 1:Nallvols-1
    SLICELOCS{v} = LOCS(intersect(find(LOCS>=VOLLOCS(v), Nslices, 'first'), ...
        find(LOCS<VOLLOCS(v+1))));
    if length(SLICELOCS{v})~=Nslices && ~hasWarned
        verbose = tapas_physio_log(sprintf(...
            'First volume with missing slice events (%d: %d instead of %d found)\n', ...
            Nallvols, length(SLICELOCS{v}), Nslices), ...
            verbose, 1);
        hasWarned = true;
    end

end

SLICELOCS{Nallvols} = LOCS(find(LOCS>=VOLLOCS(Nallvols), Nslices, 'first'));

if length(SLICELOCS{Nallvols})~=Nslices && ~hasWarned
    verbose = tapas_physio_log(sprintf(...
        'First volume with missing slice events (%d: %d instead of %d found)\n', ...
        Nallvols, length(SLICELOCS{Nallvols}), Nslices), ...
        verbose, 1);
end

if verbose.level >= 3
    stringTitle =  'Sync: Slice bundles belonging to 1 volume';
    verbose.fig_handles(end+1) = tapas_physio_get_default_fig_params();
    set(gcf, 'Name', stringTitle);
    for v=1:Nallvols-1, stem(t(SLICELOCS{v}),ones(size(SLICELOCS{v})));hold all;end
    title(stringTitle);
    xlabel('t (seconds since SCANPHYSLOG-start)');
end

try
    if do_count_from_start
        REALLOCS = cell2mat(SLICELOCS(Nprep+Ndummies+(1:Nscans))');
        DUMMYLOCS = cell2mat(SLICELOCS(Nprep+(1:Ndummies))');
    else
        REALLOCS = cell2mat(SLICELOCS(end-Nscans+1:end)');
        DUMMYLOCS = cell2mat(SLICELOCS(end-Nscans-Ndummies+1:end-Nscans)');
    end
catch err
    if exist('REALLOCS', 'var')
        verbose = tapas_physio_log(...
            sprintf('Number of dummy slice event bundles does not match Ndummies. Check parameter sqpar.Ndummies or set to 0.\n\t (Matlab error: %s)', ...
            err.message), verbose, 2);
    else
        verbose = tapas_physio_log(...
            sprintf('Number of slice event bundles does not match Nscans. Check parameters sqpar.Nscans, Ndummies and Nprep.\n\t (Matlab error: %s)', ...
            err.message), verbose, 2);
    end
end

DUMMYLOCS = reshape(DUMMYLOCS,length(DUMMYLOCS),1);
REALLOCS = reshape(REALLOCS,length(REALLOCS),1);

ons.acq_slice_all = [DUMMYLOCS; REALLOCS];

if do_count_from_start
    ons.acq_vol_all   = VOLLOCS(Nprep+(1:Nscans+Ndummies));
    ons.acq_slice_per_vol_all = SLICELOCS(Nprep+(1:Nscans+Ndummies));
else
    ons.acq_vol_all   = VOLLOCS(end-Ndummies-Nscans+1:end);
    ons.acq_slice_per_vol_all = SLICELOCS(end-Nscans-Ndummies+1:end);
end

%% set scan time (in phys time scale) corresponding to index, where valid index exists
% (i.e. for scan events occuring during phys logfile

spulse = NaN(size(ons.acq_slice_all));
svolpulse = NaN(size(ons.acq_vol_all));

iValidSlice = find(~isnan(ons.acq_slice_all));
spulse(iValidSlice) = t(ons.acq_slice_all(iValidSlice));

iValidVol = find(~isnan(ons.acq_vol_all));
svolpulse(iValidVol) = t(ons.acq_vol_all(iValidVol));

spulse_per_vol = cellfun(@(x) t(x), ons.acq_slice_per_vol_all, 'UniformOutput', false);

author: translationalneuromodeling

repo: tapas


text:
function [sbc, fpe, logdp, np] = arord(R, m, mcor, ne, pmin, pmax)
%ARORD	Evaluates criteria for selecting the order of an AR model.
%
%  [SBC,FPE]=ARORD(R,m,mcor,ne,pmin,pmax) returns approximate values
%  of the order selection criteria SBC and FPE for models of order
%  pmin:pmax. The input matrix R is the upper triangular factor in the
%  QR factorization of the AR model; m is the dimension of the state
%  vectors; the flag mcor indicates whether or not an intercept vector
%  is being fitted; and ne is the number of block equations of size m
%  used in the estimation. The returned values of the order selection
%  criteria are approximate in that in evaluating a selection
%  criterion for an AR model of order p < pmax, pmax-p initial values
%  of the given time series are ignored.
%
%  ARORD is called by ARFIT.
%
%  See also ARFIT, ARQR.

%  For testing purposes, ARORD also returns the vectors logdp and np,
%  containing the logarithms of the determinants of the (scaled)
%  covariance matrix estimates and the number of parameter vectors at
%  each order pmin:pmax.

%  Modified 17-Dec-99
%  Author: Tapio Schneider
%          tapio@gps.caltech.edu

  imax 	  = pmax-pmin+1;        % maximum index of output vectors

  % initialize output vectors
  sbc     = zeros(1, imax);     % Schwarz's Bayesian Criterion
  fpe     = zeros(1, imax);     % log of Akaike's Final Prediction Error
  logdp   = zeros(1, imax);     % determinant of (scaled) covariance matrix
  np      = zeros(1, imax);     % number of parameter vectors of length m
  np(imax)= m*pmax+mcor;

  % Get lower right triangle R22 of R:
  %
  %   | R11  R12 |
  % R=|          |
  %   | 0    R22 |
  %
  R22     = R(np(imax)+1 : np(imax)+m, np(imax)+1 : np(imax)+m);

  % From R22, get inverse of residual cross-product matrix for model
  % of order pmax
  invR22  = inv(R22);
  Mp      = invR22*invR22';

  % For order selection, get determinant of residual cross-product matrix
  %       logdp = log det(residual cross-product matrix)
  logdp(imax) = 2.*log(abs(prod(diag(R22))));

  % Compute approximate order selection criteria for models of
  % order pmin:pmax
  i = imax;
  for p = pmax:-1:pmin
    np(i)      = m*p + mcor;	% number of parameter vectors of length m
   if p < pmax
      % Downdate determinant of residual cross-product matrix
      % Rp: Part of R to be added to Cholesky factor of covariance matrix
      Rp       = R(np(i)+1:np(i)+m, np(imax)+1:np(imax)+m);

      % Get Mp, the downdated inverse of the residual cross-product
      % matrix, using the Woodbury formula
      L        = chol(eye(m) + Rp*Mp*Rp')';
      N        = L \ Rp*Mp;
      Mp       = Mp - N'*N;

      % Get downdated logarithm of determinant
      logdp(i) = logdp(i+1) + 2.* log(abs(prod(diag(L))));
   end

   % Schwarz's Bayesian Criterion
   sbc(i) = logdp(i)/m - log(ne) * (ne-np(i))/ne;

   % logarithm of Akaike's Final Prediction Error
   fpe(i) = logdp(i)/m - log(ne*(ne-np(i))/(ne+np(i)));

   % Modified Schwarz criterion (MSC):
   % msc(i) = logdp(i)/m - (log(ne) - 2.5) * (1 - 2.5*np(i)/(ne-np(i)));

   i      = i-1;                % go to next lower order
end


text:
function D = spm_eeg_correct_sensor_data(S)
% Remove artefacts from the data based on their topography
% FORMAT D = spm_eeg_correct_sensor_data(S)
%
% S      - input structure (optional)
% (optional) fields of S:
%   S.D    - MEEG object or filename of M/EEG mat-file
%   S.mode - 'SSP': simple projection
%          - 'Berg': the method of Berg (see the reference below)
% Output:
% D      - MEEG object (also written on disk)
%
% Implements:
%   Berg P, Scherg M.
%   A multiple source approach to the correction of eye artifacts.
%   Electroencephalogr Clin Neurophysiol. 1994 Mar;90(3):229-41.
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_correct_sensor_data.m 7701 2019-11-21 21:50:17Z vladimir $

SVNrev = '$Rev: 7701 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','Correct sensor data');

if ~isfield(S, 'mode') && isfield(S, 'correction'),  S.mode  = S.correction;  end
if ~isfield(S, 'prefix'),                            S.prefix   = 'T';        end

%-Get MEEG object
%--------------------------------------------------------------------------
D = spm_eeg_load(S.D);

inputfile = fullfile(D);

if ~any(D.sconfounds)
    D = spm_eeg_spatial_confounds(S);
    if ~any(D.sconfounds)
        return;
    end
end

[mod, list] = modality(D, 1, 1);

A = {};
if isequal(mod, 'Multimodal')
    sconf = getfield(D, 'sconfounds');

    for i = 1:numel(list)
        chanind = indchantype(D, list{i}, 'GOOD');
        [sel1, sel2] = spm_match_str(chanlabels(D, chanind), sconf.label);

        if any(sconf.bad(sel2))
            error(['Channels ' sprintf('%s ', sconf.label{sel2}) ' should be set to bad.']);
        end

        A{i} = sconf.coeff(sel2, :);
    end
else
    A = {D.sconfounds};
    list = {mod};
end

Dorig = D;

for i = 1:numel(A)
    label = D.chanlabels(indchantype(D, list{i}, 'GOOD'))';

    montage = [];
    montage.labelorg = label;
    montage.labelnew = label;

    montage.chantypeold = lower(D.chantype(D.indchannel(label)))';
    montage.chantypenew  = lower(montage.chantypeold);

    montage.chanunitold = D.units(D.indchannel(label))';
    montage.chanunitnew  = montage.chanunitold;

    if size(A{i}, 1)~=numel(label)
        error('Spatial confound vector does not match the channels.');
    end

    if isequal(lower(S.mode), 'berg')
        % These are the locations taken from the file BR_Brain Regions_LR.bsa
        % in BESA distribution, transformed to Tailarach coordinates using
        % BESA simulator and then to MNI template space using tal2icbm_spm
        % from http://brainmap.org/icbm2tal/
        sources = [
            -47.3124    6.4922  -10.3381
            -49.1870  -38.7590   -4.3574
            -35.2804   38.7888   21.0944
            -41.0574  -16.4018   46.0553
            -34.9993  -70.8348   20.8450
            0.8786   60.5733   -5.2797
            1.5716   38.5344   44.5853
            1.9792  -16.5380   65.3503
            1.8525  -71.0130   44.3363
            1.2689  -91.6233   -5.6259
            49.1987    6.8326  -12.0156
            51.4597  -38.4040   -6.1068
            37.8063   39.0466   19.8240
            44.5032  -16.1000   44.5681
            38.0874  -70.5770   19.5747
            ];

        [D, ok] = check(D, 'sensfid');

        if ~ok
            if check(D, 'basic')
                error(['The requested file is not ready for source reconstruction.'...
                    'Use prep to specify sensors and fiducials.']);
            else
                error('The meeg file is corrupt or incomplete');
            end
        end

        %-Find or prepare head model
        %==================================================================

        if ~isfield(D, 'val')
            D.val = 1;
        end

        if ~isfield(D, 'inv') || ~iscell(D.inv) ||...
                ~(isfield(D.inv{D.val}, 'forward') && isfield(D.inv{D.val}, 'datareg')) ||...
                ~isa(D.inv{D.val}.mesh.tess_ctx, 'char') % detects old version of the struct
            D = spm_eeg_inv_mesh_ui(D, D.val);
            D = spm_eeg_inv_datareg_ui(D, D.val);
            D = spm_eeg_inv_forward_ui(D, D.val);

            save(D);
        end

        fwd = spm_eeg_inv_get_vol_sens(D, D.val, 'MNI-aligned', 'inv', list{i});

        [vol, sens] = ft_prepare_vol_sens(fwd.(list{i}).vol, fwd.(list{i}).sens, 'channel', label);

        L = ft_compute_leadfield(spm_eeg_inv_transform_points(inv(fwd.transforms.toMNI), sources), sens, vol);
        %[L, D] = spm_eeg_lgainmat(D, [], label);

        B = spm_svd(L*L', 0.01);

        lim = min(0.5*size(L, 1), 45); % 45 is the number of dipoles BESA would use.

        if size(B, 2) > lim;
            B = B(:, 1:lim);
        end

        SX = full([A{i} B]);

        SXi = pinv(SX);
        SXi = SXi(1:size(A{i}, 2), :);

        montage.tra = eye(size(A{i}, 1)) - A{i}*SXi;
    else
        montage.tra = eye(size(A{i}, 1)) - A{i}*pinv(A{i});
    end

    S1   = [];
    S1.D = D;
    S1.montage = montage;
    S1.keepothers = true;
    S1.updatehistory  = false;

    Dnew = spm_eeg_montage(S1);

    if isfield(D,'inv')
        Dnew.inv = D.inv;
    end

    if i>1
        delete(D);
    end

    D = Dnew;
end

%-Change the channel order to the original order
%==========================================================================
tra = eye(D.nchannels);
montage = [];
montage.labelorg = D.chanlabels';
montage.labelnew = Dorig.chanlabels';

montage.chantypeold  = lower(D.chantype)';
montage.chantypenew  = lower(Dorig.chantype)';

montage.chanunitold  = D.units';
montage.chanunitnew  = Dorig.units';

[sel1, sel2]  = spm_match_str(montage.labelnew, montage.labelorg);

montage.tra = tra(sel2, :);

S1   = [];
S1.D = D;
S1.montage = montage;
S1.keepothers = true;
S1.updatehistory  = 0;

Dnew = spm_eeg_montage(S1);

delete(D);
D = Dnew;

if ~isempty(badchannels(Dorig))
    D = badchannels(D, badchannels(Dorig), 1);
end

D = D.history(mfilename, S);
save(D);

D = move(D, spm_file(inputfile, 'prefix', S.prefix));

%-Cleanup
%--------------------------------------------------------------------------
spm('FigName', 'Correct sensor data: done');

author: spm


text:
function FDM=aggregateFDM(FDM,z,k,n)

for i=1:n
        tmp1=[];     tmp2=[];     tmp3=[];     tmp4=[];
    for j=1:k
        tmp1=[tmp1 FDM{j}{z,i}(1)];
        tmp2=[tmp2 FDM{j}{z,i}(2)];
        tmp3=[tmp3 FDM{j}{z,i}(3)];
        tmp4=[tmp4 FDM{j}{z,i}(4)];
    end
   Wj1a(i)=min(tmp1);
   Wj2a(i)=1/k*sum(tmp2);
   Wj3a(i)=1/k*sum(tmp3);
   Wj4a(i)=max(tmp4);

end
FDMtmp=[Wj1a; Wj2a; Wj3a; Wj4a];
for i=1:n
    FDM2{:,i}=FDMtmp(:,i)';
end
FDM=FDM2;


text:
classdef DisjointSetM < DisjointSet
%%DISJOINTSETM A disjoint set class that can be used to partition a set of
%              targets into clusters while keeping track of which
%              measurements are associated with which cluster.
%
%This subclass generalizes the DisjointSet class to keep track of which
%measurements caused the partitioning of the targets. The set size is the
%number of ordered targets and a target is identifiable by its index number
%in the set. The class is created for a fixed number of targets and
%measurements. The main point of this class is to facilitate the creation
%of ClusterSet classes for the targets that cluster together and for the
%measurements that are associated with each cluster.
%
%November 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

properties
    numMeasInTree
    measToNode
end
methods
    function newSet=DisjointSetM(numTar,numMeas)
        %The constructor method.
        newSet=newSet@DisjointSet(numTar);
        newSet.numMeasInTree=zeros(numTar,1);
        newSet.measToNode=zeros(numMeas,1);
    end

    function unionFromList(DSObj,uList,measNum)
    %%UNIONFROMLIST Merge clusters of targets from a given common
    %               measurement. NOTE: Calling unionFromList twice with
    %               the same measurement number will produce incorrect
    %               results.
    %
    %INPUTS: DSObj  The implicitly passed calling object.
    %        uList  A one-dimensional list of indices of targets to
    %               merge into one cluster
    %       measNum The index of the measurement.
    %
    %The uList passed for a particular measNum should a list of the
    %indices of ALL of the targets that gate with the measurement.
    %
    %November 2013 David F. Crouse, Naval Research Laboratory, Washington
    %D.C.

        numInList=length(uList);
        if(numInList==0)
            return;
        end

        rootIdx1=DSObj.find(uList(1));
        DSObj.numMeasInTree(rootIdx1)=DSObj.numMeasInTree(rootIdx1)+1;
        DSObj.measToNode(measNum)=rootIdx1;
        for curTar=2:numInList
            rootIdx2=DSObj.find(uList(curTar));
            rootIdx1=DSObj.unionRoots(rootIdx1,rootIdx2);
        end
    end

    function unionFromBoolList(DSObj,boolList,measNum)
    %%UNIONFROMBOOLLIST Merge targets specified by a boolean array for
    %               a given measurement index. NOTE: Calling
    %               unionFromBoolList twice with the same measurement
    %               number will produce incorrect results.
    %
    %INPUTS: DSObj  The implicitly passed calling object.
    %     boolList  A length numTarX1 or 1XnumTar boolean array
    %               indicating whether which targets are in the same
    %               cluster due to the specified measurement.
    %       measNum The index of the measurement.
    %
    %November 2013 David F. Crouse, Naval Research Laboratory, Washington
    %D.C.

        numTar=length(DSObj.setArray);
        %Check for input validity.
        assert(numTar==length(DSObj.setArray))

        tarIdxCur=find(boolList,1);

        if(isempty(tarIdxCur))
            %If no target gates with the specified measurement.
            return;
        end

        rootIdx1=DSObj.find(tarIdxCur);
        DSObj.numMeasInTree(rootIdx1)=DSObj.numMeasInTree(rootIdx1)+1;
        DSObj.measToNode(measNum)=rootIdx1;
        for curTarIdx=(tarIdxCur+1):numTar
            if(boolList(curTarIdx))
                rootIdx2=DSObj.find(curTarIdx);
                rootIdx1=DSObj.unionRoots(rootIdx1,rootIdx2);
            end
        end
    end

    function unionFromBinMat(DSObj,binMat)
    %UNIONFROMBINMAT Given a binary matrix specifying which
    %                measurements gate with which targets, cluster
    %                together all targets that gate with each other.
    %
    %INPUTS: DSObj The implicitly passed calling object.
    %       binMat A numTarXnumMeas binary matrix where numTar is the
    %              same as the length of setArray in DSObj and numMeas
    %              is the number of measurements with which the
    %              DisjointSet object was created.
    %
    %The measurement number is assumed from the column index when
    %clustering. This function should not be called multiple times with
    %different binary matrices as it will produce incorrect results.
    %Rather, multiple binary matrices can be OR-ed together prior to
    %calling this method.
    %
    %November 2013 David F. Crouse, Naval Research Laboratory, Washington
    %D.C.

        numTar=size(binMat,1);
        numMeas=size(binMat,2);

        %Check for input validity
        assert(numTar==length(DSObj.setArray))
        if(~isempty(DSObj.measToNode))
            assert(numMeas==length(DSObj.measToNode))
        end

        for curMeas=1:numMeas
            rootIdx1=0;
            for curTar=1:numTar
                if(binMat(curTar,curMeas)~=0)
                    if(rootIdx1==0)
                        rootIdx1=DSObj.find(curTar);

                        %Do measurement bookeeping.
                        DSObj.measToNode(curMeas)=curTar;
                        DSObj.numMeasInTree(rootIdx1)=DSObj.numMeasInTree(rootIdx1)+1;
                    else
                        rootIdx2=DSObj.find(curTar);
                        rootIdx1=DSObj.unionRoots(rootIdx1,rootIdx2);
                    end

                    %Do measurement bookeeping.
                    DSObj.measToNode(curMeas)=rootIdx1;
                end
            end
        end
    end

    function [newCSet,newCSMeasSet,ungatedMeas]=createClusterSet(DSObj)
    %%CREATECLUSTERSET Turn the DisjointSet object into a ClusterSet so
    %                  that targets in each cluster can be easily
    %                  addressed. Also produce a corresponding
    %                  ClusterSet for the clusters of measurements that
    %                  are associated with each cluster of targets and
    %                  produce a list of which measurements do not gate
    %                  with anything.
    %
    %OUTPUTS: newCSet A ClusterSet object representing the clustering
    %                 given by this DisjointSetM object.
    %    newCSMeasSet A ClusterSet object of the measurements that go
    %                 with each cluster of targets.
    %     ungatedMeas An array of indices of measurements that do not
    %                 gate with anything.
    %
    %newCSet(c,:) is the set of targets in cluster c and
    %newCSMeasSet(c,:) is the corresponding set of measurements.
    %
    %November 2013 David F. Crouse, Naval Research Laboratory, Washington
    %D.C.

        numTar=length(DSObj.setArray);
        numMeas=length(DSObj.measToNode);

        rootIdxList=zeros(DSObj.numberOfSets,1);
        curNumAdded=zeros(DSObj.numberOfSets,1);

        %This holds the numbers of targets per cluster.
        rootRankList=zeros(DSObj.numberOfSets,1);
        %This holds the numbers of measurements per cluster.
        rootRankListMeas=zeros(DSObj.numberOfSets,1);
        clusterEls=zeros(numTar,1);
        clusterElsMeas=zeros(numMeas,1);

        %First, record the indices and ranks of all of the root nodes.
        rootNodesFound=0;
        for curItem=1:numTar
            if(DSObj.setArray(curItem)<0)
                rootNodesFound=rootNodesFound+1;
                rootIdxList(rootNodesFound)=curItem;
                rootRankList(rootNodesFound)=DSObj.numInTree(curItem);
                rootRankListMeas(rootNodesFound)=DSObj.numMeasInTree(curItem);

                if(rootNodesFound==DSObj.numberOfSets)
                    break;
                end
            end
        end

        %Next, sort the root nodes by their index, so that children can
        %be efficiently assigned to them.
        [rootIdxList,idx] = sort(rootIdxList,1,'ascend');
        rootRankList=rootRankList(idx);
        %The measurement clusters must correspond to the target
        %clusters.
        rootRankListMeas=rootRankListMeas(idx);

        offsetArray=zeros(DSObj.numberOfSets,1);
        offsetArray(2:end)=cumsum(rootRankList(1:(end-1)));
        for curItem=1:numTar
            %Find the root node for the given node.
            rootNode=DSObj.find(curItem);
            %Find the index of the cluster corresponding to that root
            %node.
            [~, clusterIdx]=binSearch(rootIdxList,rootNode);
            %Add the node to the cluster.
            clusterEls(1+offsetArray(clusterIdx)+curNumAdded(clusterIdx))=curItem;
            curNumAdded(clusterIdx)=curNumAdded(clusterIdx)+1;
        end
        %Create the cluster set for the measurements
        newCSet=ClusterSet(clusterEls,rootRankList,offsetArray);
        %Create an array of the maximum possible size to hold the
        %ungated measurements.
        ungatedMeas=zeros(numMeas,1);
        numUngatedMeas=0;

        %Fill the measurement clusters in the same manner as done
        %for the regular clusters.
        curNumAdded=zeros(DSObj.numberOfSets,1);
        offsetArrayMeas=zeros(DSObj.numberOfSets,1);
        offsetArrayMeas(2:end)=cumsum(rootRankListMeas(1:(end-1)));
        for curItem=1:numMeas
            if(DSObj.measToNode(curItem)>0)
                rootNode=DSObj.find(DSObj.measToNode(curItem));
                %Find the index of the cluster corresponding to that root
                %node.
                [~, clusterIdx]=binSearch(rootIdxList,rootNode);
                %Add the node to the cluster.
                clusterElsMeas(1+offsetArrayMeas(clusterIdx)+curNumAdded(clusterIdx))=curItem;
                curNumAdded(clusterIdx)=curNumAdded(clusterIdx)+1;
            else
                numUngatedMeas=numUngatedMeas+1;
                ungatedMeas(numUngatedMeas)=curItem;
            end
        end

        %Shorten the array to the actual number of ungated
        %measurements.
        ungatedMeas=ungatedMeas(1:numUngatedMeas);

        newCSMeasSet=ClusterSet(clusterElsMeas,rootRankListMeas,offsetArrayMeas);
    end

    function mergedRootIdx=unionRoots(DSObj,rootIdx1,rootIdx2)
    %%unionRoots Merge two sets given their root nodes, keeping
    %            track of the number of measurements associated
    %            with both sets.
    %
    %INPUTS: DSObj      The implicitly passed calling object.
    %        rootIdx1   The root index of the first set to merge.
    %        rootIdx2   The root index of the second set to merge.
    %
    %OUTPUTS: mergedRootIdx The index of the root of the merged
    %                       cluster. It will either be rootIdx1 or
    %                       rootIdx2.
    %
    %This is a union-by-rank algorithm. The merged cluster has a number
    %of associated measurements equal to the sum of the number of
    %measurements in both clusters. This method should not be called
    %with indices that do not correspond to root nodes.
    %
    %November 2013 David F. Crouse, Naval Research Laboratory, Washington
    %D.C.

        %If the roots are the same, then there is nothing to merge.
        if(rootIdx1==rootIdx2)
            %Measurement bookeeping
            DSObj.numMeasInTree(rootIdx1)=DSObj.numMeasInTree(rootIdx1);
            mergedRootIdx=rootIdx1;
            return;
        end

        %If root 2 is deeper, than make root2 the new root.
        if(DSObj.setArray(rootIdx2)< DSObj.setArray(rootIdx1))
            DSObj.setArray(rootIdx1)=rootIdx2;
            DSObj.numInTree(rootIdx2)=DSObj.numInTree(rootIdx1)+DSObj.numInTree(rootIdx2);
            mergedRootIdx=rootIdx2;

            DSObj.numMeasInTree(rootIdx2)=DSObj.numMeasInTree(rootIdx1)+DSObj.numMeasInTree(rootIdx2);
        else
            %If they are the same height, then update the height of
            %rootIdx1.
            if(DSObj.setArray(rootIdx1)== DSObj.setArray(rootIdx2))
                DSObj.setArray(rootIdx1)=DSObj.setArray(rootIdx1)-1;
            end
            DSObj.setArray(rootIdx2)=rootIdx1;
            DSObj.numInTree(rootIdx1)=DSObj.numInTree(rootIdx1)+DSObj.numInTree(rootIdx2);
            mergedRootIdx=rootIdx1;

            DSObj.numMeasInTree(rootIdx1)=DSObj.numMeasInTree(rootIdx2)+DSObj.numMeasInTree(rootIdx1);
        end

        %Merging reduces the number of sets by one.
        DSObj.numberOfSets=DSObj.numberOfSets-1;
    end
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.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y=evalSinCosSeries(c,cosX,sinX)
%%EVALSINCOSSERIES Evaluate the cosine series
%                  y=sum_{i=1}^Nc(i)*cos((i-1)*x)
%                  or the sine series
%                  y=sum_{i=1}^(N+1)c(i)*sin(i*x)
%                  efficiently using Clenshaw summation.
%
%INPUTS: c An NX1 or a 1XN vector of the coefficients in the sum.
%     cosX The value cos(x). This value is always required. To evaluate
%          multiple sums at once, this can be a vector or a matrix.
%     sinX If this value is provided, it is sin(x) and the sum evaluated by
%          this function will be a sine series. Otherwise, if this
%          parameter is omitted or an empty matrix is passed, this function
%          will evaluate a cosine series. If this parameter is provided, it
%          must be the same size as cosX.
%
%OUTPUTS: y The value of the sum. If cosX was a matrix, then this will be a
%           matrix.
%
%This function chooses the correct parameters and calls the function
%evalClenshawRecurSeries.
%
%EXAMPLE:
% N=200;
% c=rand(N,1);
% x=2*pi*rand(4,4);
% cosX=cos(x);
% yF=evalSinCosSeries(c,cosX);
% ySum=zeros(4,4);
% for k=1:N
%    ySum=ySum+c(k)*cos((k-1)*x);
% end
% relativeError=max(max(abs(ySum-yF)./yF))
% %This will typically be on the order of 1e-12 or less, due simply to
% %finite precision differences.
% %Similarly, if we wanted the sine series.
% sinX=sin(x);
% yF=evalSinCosSeries(c,cosX,sinX);
% ySum=zeros(4,4);
% for k=1:N
%     ySum=ySum+c(k)*sin(k*x);
% end
% relativeError=max(max(abs(ySum-yF)./yF))
% %This will again be on the order of 1e-12 or less, due simply to finite
% %precision differences.
%
%July 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

alphaVal=2*cosX;
betaVal=-1;

if(nargin<3||isempty(sinX))
    %A cosine series.
    F0=ones(size(cosX));%=cos(0)
    F1=cosX;
    y=evalClenshawRecurSeries(c,alphaVal,betaVal,F0,F1);
else%A sine series
    F0=sinX;
    F1=2*sinX.*cosX;
    y=evalClenshawRecurSeries(c,alphaVal,betaVal,F0,F1);
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:
function result = getResult(s,Y)
nCorrect = 0;
for i = 1:length(s)
    if s(i) > 0.5 & Y(i) ~= -1
        nCorrect = nCorrect + 1;
    end
    if s(i) < 0.5 & Y(i) ~= 1
         nCorrect = nCorrect + 1;
    end
end
result = nCorrect/length(s);

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function OverlapRatio= overlap(ROI1,ROI2,ImageSize)

Image=zeros(ImageSize);

Xpoint=mod(double(ROI1(1,:))',ImageSize(1));
Ypoint=mod(double(ROI1(2,:))',ImageSize(2));
Xpoint(Xpoint==0)=ImageSize(1);
Ypoint(Ypoint==0)=ImageSize(2);
for m=1:1:length(Xpoint)
      Image(Ypoint(m),Xpoint(m))=1;
end

Xpoint=mod(double(ROI2(1,:))',ImageSize(1));
Ypoint=mod(double(ROI2(2,:))',ImageSize(2));
Xpoint(Xpoint==0)=ImageSize(1);
Ypoint(Ypoint==0)=ImageSize(2);
for m=1:1:length(Xpoint)
      Image(Ypoint(m),Xpoint(m))=Image(Ypoint(m),Xpoint(m))+1;
end

OverlapRatio(1)=length(find(Image(:)==2))./size(ROI1,2);
OverlapRatio(2)=length(find(Image(:)==2))./size(ROI2,2);
OverlapRatio(3)=length(find(Image(:)==2))./(length(find(Image(:)==2))+length(find(Image(:)==1)));

end

author: kavli-ntnu


text:
function edges=meshedge(elem,varargin)
%
% edges=meshedge(elem,opt)
%
% return all edges in a surface or volumetric mesh
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2011/02/26
%
% input:
%    elem:  element table of a mesh (support N-d space element)
%    opt: optional input, giving the additional options. If opt
%         is a struct, it can have the following field:
%       opt.nodeorder: if 1, assuming the elem node indices is in CCW
%                      orientation; 0 use nchoosek() output to order edges
%         you can replace opt by a series of ('param', value) pairs.
%
% output:
%    edge:  edge list; each row is an edge, specified by the starting and
%           ending node indices, the total edge number is
%           size(elem,1) x nchoosek(size(elem,2),2). All edges are ordered
%           by looping through each element first.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

dim=size(elem);
edgeid=nchoosek(1:dim(2),2);
len=size(edgeid,1);
edges=zeros(dim(1)*len,2);
for i=0:len-1
    edges((i*dim(1)+1):((i+1)*dim(1)),:)=[elem(:,edgeid(i+1,1)) elem(:,edgeid(i+1,2))];
end


text:
function [mSynCat] = create_syncat(mCatalog,nCdf,fMu,fSigma, fMinMagIn,bPlot)
% function [mSynCat] = create_syncat(mCatalog,nCdf,fMu,fSigma,fMinMagIn,bPlot);
% --------------------------------------------------
% Creates a synthetic catalog for the part below Mc with a Normal CDF.
% Usage:
% Before using this function, create a defined synthetic catalog according
% to a power-law behavior. It's assumed that the smallest magnitude bin is
% the bin of the magnitude of completeness. If fMu, fSigma and fMinMagIn are
% missing, enter interactively.
%
% Incoming variable:
% mCatalog : EQ catalog
% nCdf : Choice of cumulative distribution function
%        1 : Normal cumulative distribution function
%        2 : Lognormal cumulative distribution function
%        3 : Weibull CDF
%
% fMu    : Mu of Normal / Lognormal CDF
% fSigma : Sigma of Normal / Lognormal CDF
% fMinMagIn : Minimum magnitude for synthetic catalog
% bPlot : Plot historgram, 0 =no plot, 1 = plot
% Output:
% mSynCat : Synthetic catalog
% mCat : Catalog of events below Mc
%
% Authour: J. Woessner, j.woessner@sed.ethz.ch
% last update: 10.07.04

% Initialize
mCat = [];

if ~exist('bPlot','var')
    bPlot = 0;
end

% Check for input values
if nargin < 3
    % Calculate probablities for Normal CDF
    prompt = {'Enter mu:','Enter sigma:','Minimum magnitude:'};
    dlg_title = 'Parameters for normal CDF';
    num_lines= 1;
    def     = {'0.8','0.4','0'};
    answer  = inputdlg(prompt,dlg_title,num_lines,def);
    fMu = str2double(answer(1));
    fSigma = str2double(answer(2));
    fMinMagIn = str2double(answer(3));
end

% Calculate FMD
[mFMDC, mFMD] = calc_FMD(mCatalog);

% Find first bin (Mc bin) with data
vSel = (max(mFMD(2,:)) == mFMD(2,:));
fN_Mc = mFMD(2,vSel);

fMinBin = mFMD(1,vSel);
vMagstep = fMinMagIn:0.1:fMinBin-0.1;

% Choose CDF
switch nCdf
    case 1
        vProb = normcdf(vMagstep,fMu, fSigma);
    case 2
        vProb = logncdf(vMagstep,fMu,fSigma);
    case 3
        vProb = wblcdf(vMagstep,fMu,fSigma);
    otherwise
        disp('Check nCdf parameter!');
        return;
end

vProb = vProb';
vMagstep = vMagstep';

% Calculate number of EQs in bins
vN = round(vProb(:,1)*fN_Mc);
mData = [vMagstep vN];
vMag = [];
nCount=1;
for nMag=fMinMagIn:0.1:fMinBin-0.1
    %fM = repmat(nMag,mData(floor(abs(nMag)*10+1),2),1);
    fM = repmat(nMag,mData(nCount,2),1);
    vMag = [vMag; fM];
    nCount=nCount+1;
end

% Choose appropriate number of events from the catalog, change magnitudes
% and add to original catalog
vInd = round(rand(length(vMag),1)*length(mCatalog(:,1)));
% Avoid zeros
vIndice = find(vInd == 0);
vInd(vIndice) = vInd(vIndice)+1;

for n=1:length(vInd)
    mCat = [mCat; mCatalog(vInd(n),:)];
end
%mCat = mCatalog(vSel,:);
%mCat = mCatalog(1:length(vMag),:);
mCat(:,6) = vMag;
mSynCat = [mCat;mCatalog];

if (exist('bPlot','var') & bPlot == 1)
    %Plot result
    fMaxMag = max(mCatalog(:,6));
    figure_w_normalized_uicontrolunits('tag','maghist');
    histogram(mSynCat(:,6),fMinMagIn:0.1:fMaxMag);
    xlabel('Magnitude');
    ylabel('Frequency');
end

author: CelsoReyes


text:
function result = set_concat(varargin)
% Concatenate continuous signals across time.
% Result = set_joinepos(Set1, Set2, ...)
%
% In:
%   SetK   : The k'th data set to concatenate.
%
% Out:
%   Result : A new data set that is the concatenation of all input sets. The following changes are made:
%            * .data and all other time-series fields are concatenated across time (2nd dimension)
%            * .event is joined and .latency fields are updated appropriately
%            * .xmax/.pnts are updated
%
% Notes:
%   This function returns a new data set with meta-data set to that of the first input set, and the
%   time series fields joined across all sets. No checks for meta-data consistency are done. There
%   is a heavy-duty function for merging inconsistent sets called set_merge, which can merge cats
%   and dogs. This function does not attempt to keep miscellaneous EEGLAB meta-data consistent,
%   including: setname,filename,filepath,subject,group,condition,session,comments,urevent,reject,stats,history,etc
%
% Examples:
%   % concatenate data sets eegA, eegB and eegC across time
%   eeg = set_concat(eegA,eegB,eegC)
%
% See also:
%   set_joinepos, set_merge
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-03-31
dp;

% set_concat_version<1.0> -- for the cache

if ~exp_beginfun('editing') return; end

declare_properties('name','Concatenate','independent_channels',true,'independent_trials',false);

% input validation
for k=1:length(varargin)
    utl_check_fields(varargin{k},{'data','event','pnts','xmin','xmax','srate'},'input','signal');
    if ~isempty(varargin{k}.event)
        if ~isfield(varargin{k}.event,'latency')
            error('A dataset passed to set_concat is lacking the .event.latency field.'); end
        latency_numels = cellfun('prodofsize',{varargin{k}.event.latency});
        if any(latency_numels == 0)
            error('One or more of the events in the given data set have an empty .latency field, which is not permitted.'); end
        if any(latency_numels ~= 1)
            error('One or more of the events in the given data set have a .latency value that is not a scalar, which is not permitted.'); end
    end
    if isfield(varargin{k},'epoch') && ~isempty(varargin{k}.epoch)
        error('Only continuous data can be concatenated with set_concat -- use set_joinepos for epoched data.'); end
end

if ~isempty(varargin)
    result = varargin{1};
    if length(varargin) > 1
        % concatenate time series fields
        for field = utl_timeseries_fields(result)
            data = cellfun(@(x)x.(field{1}),varargin,'UniformOutput',false);
            try
                result.(field{1}) = cat(2,data{:});
            catch e
                % concatenation failed: produce a reasonable error message
                if ~isempty(data)
                    sizes = cellfun('size',data,1);
                    sizes = sizes(sizes ~= 0);
                    if length(sizes) > 1 && ~all(sizes==sizes(1))
                        error('The time-series field .%s must have the same number of channels in each data set.',field{1}); end
                    if any(cellfun('size',data,3) > 1)
                        error('One or more of the datasets passed to set_concat were epoched -- use set_joinepos to concatenate epoched data.'); end
                end
                size_info = hlp_tostring(cellfun(@size,data,'UniformOutput',false));
                error('Concatenation of time-series fields failed with error: %s (the data sizes were %s -- make sure that they are mutually compatible).',e.message,size_info);
            end
            if isempty(result.(field{1}))
                result.(field{1}) = []; end
        end
        % count events, epochs and samples in each set
        event_count = cellfun(@(x)length(x.event),varargin);
        sample_count = cellfun(@(x)x.pnts,varargin);
        % concatenate .event and .epoch fields
        event = cellfun(@(x)x.event,varargin,'UniformOutput',false); result.event = [event{:}];
        % shift event latencies based on cumulative sample counts
        if ~isempty(result.event)
            [result.event.latency] = arraydeal([result.event.latency]+replicate(cumsum(sample_count)-sample_count,event_count)); end
        % update misc fields
        [result.nbchan,result.pnts,result.trials,extra_dims] = size(result.data); %#ok<ASGLU,NASGU>
        result.xmax = result.xmin + (result.pnts-1)/result.srate;
    end
else
    result = struct('setname','','filename','','filepath','','subject','','group','','condition','','session',[],'comments','','nbchan',0,...
        'trials',0,'pnts',0,'srate',1,'xmin',0,'xmax',0,'times',[],'data',[],'icaact',[],'icawinv',[],'icasphere',[],'icaweights',[], ...
        'icachansind',[],'chanlocs',[],'urchanlocs',[],'chaninfo',[],'ref',[],'event',[],'urevent',[],'eventdescription',{{}}, ...
        'epoch',[],'epochdescription',{{}},'reject',[],'stats',[],'specdata',[],'specicaact',[],'splinefile','','icasplinefile','', ...
        'dipfit',[],'history','','saved','no','etc',[]);
end

exp_endfun;

function result = replicate(values,counts)
% Replicate each element Values(k) by Count(k) times.
result = zeros(1,sum(counts));
k = 0;
for p=find(counts)
    result(k+(1:counts(p))) = values(p);
    k = k+counts(p);
end


text:
function [err,time,solver,eqn,node,elem] = afemPoisson(mesh,pde,option,varargin)

%% Check input arguments
if isfield(mesh,'node') && isfield(mesh,'elem')
    node = mesh.node;
    elem = double(mesh.elem);
end
if ~exist('node','var') || ~exist('elem','var')
    % default mesh: Lshape
    [node,elem] = squaremesh([-1,1,-1,1],1);
    [node,elem] = delmesh(node,elem,'x>0 & y<0');
end
if isfield(mesh,'bdFlag')
    bdFlag = mesh.bdFlag;
else % default boundary condition
    bdFlag = setboundary(node,elem,'Dirichlet');
end
if ~exist('option','var'), option = []; end
if ~exist('pde','var')
    pde = Lshapedata;                          % default data
end

%% Parameters
option = afemoption(option,2);
maxIt = option.maxIt;
refType = option.refType;
elemType = option.elemType;
theta = option.theta;

%% Initialize err
errL2 = zeros(maxIt,1);   errH1 = zeros(maxIt,1); erreta = zeros(maxIt,1);
erruIuh = zeros(maxIt,1); errMax = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
estimateTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1);

%% Generate an initial mesh
for k = 1:option.L0
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    end
end

%%  Adaptive Finite Element Method
% *SOLVE* -> *ESTIMATE* -> *MARK* -> *REFINE*
for k = 1:maxIt
    % Step 1: SOLVE
    switch elemType
        case 'P1'     % piecewise linear function P1 element
            [soln,eqn,info] = Poisson(node,elem,bdFlag,pde,option);
        case 'CR'     % piecewise linear function CR element
            [soln,eqn,info] = PoissonCR(node,elem,bdFlag,pde,option);
        case 'P2'     % piecewise quadratic function
            [soln,eqn,info] = PoissonP2(node,elem,bdFlag,pde,option);
        case 'WG'     % weak Galerkin element
            [soln,eqn,info] = PoissonWG(node,elem,bdFlag,pde,option);
    end
    uh = soln.u;
    % compute error
    t = cputime;
    if isfield(pde,'Du')
        if ~isfield(pde,'d')
            pde.d = [];
        end
        if isfield(soln,'Du') && ~isempty(soln.Du) % Du is in the output
            errH1(k) = getH1error(node,elem,pde.Du,soln.Du,pde.d);
        else
            errH1(k) = getH1error(node,elem,pde.Du,soln.Du,pde.d);
        end
    end
    if isfield(pde,'exactu')
        errL2(k) = getL2error(node,elem,pde.exactu,uh);
        % interpolation
        switch elemType
            case 'P1'
                uI = Lagrangeinterpolate(pde.exactu,node,elem);
            case 'Q1'
                uI = Lagrangeinterpolate(pde.exactu,node,elem);
            case 'CR'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'CR',eqn.edge);
            case 'P2'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'P2',eqn.edge);
            case 'P3'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'P3',eqn.edge);
            case 'WG'
                uI = Lagrangeinterpolate(pde.exactu,node,elem,'WG',eqn.edge);
        end
        erruIuh(k) = sqrt((uh-uI)'*eqn.A*(uh-uI));
        errMax(k) = max(abs(uh-uI));
    end
    errTime(k) = cputime - t;
    % record time
    solverTime(k) = info.solverTime;
    assembleTime(k) = info.assembleTime;
    if option.printlevel>1
        fprintf('Time to compute the error %4.2g s \n H1 err %4.2g    L2err %4.2g \n',...
            errTime(k),errH1(k), errL2(k));
    end
    % record solver information
    itStep(k) = info.itStep;
    stopErr(k) = info.stopErr;
    flag(k) = info.flag;
    % plot
    N(k) = length(uh);
    if  strcmp(elemType,'WG') % modify size for WG
        if ~isfield(option,'reducesystem') || (option.reducesystem == 1)
            N(k) = N(k) - size(elem,1); % reduced system
        end
    end
    if option.plotflag % show mesh and solution
       figure(1);
       showresult(node,elem,uh,option.viewangle);
    end
    % Step 2: ESTIMATE
    t = cputime;
    switch option.estType
        case 'recovery' % recovery type
            eta = estimaterecovery(node,elem,uh);
        case 'residual' % residual type
            switch elemType
                case 'P1'
                    eta = estimateresidual(node,elem,uh,pde,bdFlag);
                case 'WG'
                    eta = estimateresidualWG(node,elem,uh,soln.Du,pde);
            end
    end
    erreta(k) = sqrt(sum(eta.^2));
    estimateTime(k) = cputime - t;
    % Step 3: MARK
    switch option.markType
        case 'L2'
            markedElem = mark(elem,eta,theta);
        case 'MAX'
            markedElem = mark(elem,eta,theta,'MAX');
    end
    % Step 4: REFINE
    t = cputime;
    if (N(k) > option.maxN) || (k == maxIt) % no refinement
        break;
    end
    [node,elem,bdFlag,HB,tree] = bisect(node,elem,markedElem,bdFlag); %#ok<ASGLU>
    % Step 4.2: COARSEN
    if option.coarsenflag
        eta = eleminterpolate(eta,tree);
        markedElem = mark(elem,eta,0.25*theta,'COARSEN');
        [node,elem,bdFlag] = coarsen(node,elem,markedElem,bdFlag);
    end
    meshTime(k) = cputime - t;
end

%% Plot convergence rates
if option.rateflag
    if ~isfield(option,'rateshift')
        option.rateshift = 10;
    end
    figure;
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.55,0.4]);
    subplot(1,2,1)
    showrate2(N(1:k),errH1(1:k),option.rateshift,'-*','|| Du-Du_h||',...
              N(1:k),errL2(1:k),option.rateshift,'k-+','|| u-u_h||');
    subplot(1,2,2)
    showrate2(N(1:k),erruIuh(1:k),option.rateshift,'m-+','||Du_I-Du_h||',...
              N(1:k),errMax(1:k),option.rateshift,'r-*','||u_I-u_h||_{\infty}');
end

%% Output
err = struct('N',N(1:k),'H1',errH1(1:k),'L2',errL2(1:k),...
             'uIuhH1',erruIuh(1:k),'uIuhMax',errMax(1:k),'eta',erreta(1:k));
time = struct('N',N(1:k),'err',errTime(1:k),'solver',solverTime(1:k), ...
              'assemble',assembleTime(1:k),'mesh',meshTime(1:k),...
              'estimate',estimateTime(1:k));
solver = struct('N',N(1:k),'itStep',itStep(1:k),'time',solverTime(1:k),...
                'stopErr',stopErr(1:k),'flag',flag(1:k));

%% Display error and CPU time
if option.dispflag
    if ~isfield(option,'dispspace')
       option.dispspace = 5;    % display error for every 5 iterations
    end
    idx = 1:option.dispspace:k;

    fprintf('\n');
    disp('Table: Error')
    colname = {'#Dof','||u-u_h||','||Du-Du_h||','||DuI-Du_h||','||uI-u_h||_{max}','eta'};
    disptable(colname,err.N(idx),[],err.L2(idx),'%0.5e',err.H1(idx),'%0.5e',...
          err.uIuhH1(idx),'%0.5e',err.uIuhMax(idx),'%0.5e',err.eta(idx),'%0.5e');
%
%     disp('Table: CPU time')
%     colname = {'#Dof','Assemble','Solve','Error','Mesh'};
%     disptable(colname,time.N,[],time.assemble,'%0.2e',time.solver,'%0.2e',...
%                       time.err,'%0.2e',time.mesh,'%0.2e');
    fprintf('\n');
end

repo: ifem


text:
function varargout=imagehomog(varargin)
%
% For calling details please see v_imagehomog.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_imagehomog.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: imagehomog.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_imagehomog(varargin{:});
else
    v_imagehomog(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [germs, germPaths] = centroidalVoronoi2d(germs, poly, varargin)
%CENTROIDALVORONOI2D Centroidal Voronoi tesselation within a polygon.
%
%   PTS = centroidalVoronoi2d(NPTS, POLY)
%   Generate points in a polygon based on centroidal voronoi tesselation.
%   Centroidal germs can be computed by using the Llyod's algorithm:
%   1) initial germs are chosen at random within polygon
%   2) voronoi polygon of the germs is computed
%   3) the centroids of each domain are computed, and used as germs of the
%   next iteration
%
%   [PTS, PATHLIST] = centroidalVoronoi2d(NPTS, POLY)
%   Also returns the path of each germs at each iteration. The result
%   PATHLIST is a cell array with as many cells as the number of germs,
%   containing in each cell the successive positions of the germ.
%
%   PTS = centroidalVoronoi2d(.., PARAM, VALUE)
%   Specify one or several optional arguments. PARAM can be one of:
%   * 'nIter'   specifies the number of iterations of the algorithm
%       (default is 50)
%   * 'verbose' display iteration number. Default is false.
%
%   Example
%     poly = ellipseToPolygon([50 50 40 30 20], 200);
%     nGerms = 100;
%     germs = centroidalVoronoi2d(nGerms, poly);
%     figure; hold on;
%     drawPolygon(poly, 'k');
%     drawPoint(germs, 'bo');
%     axis equal; axis([0 100 10 90]);
%     % extract regions of the CVD
%     box = polygonBounds(poly);
%     [n, e] = boundedVoronoi2d(box, germs);
%     [n2, e2] = clipGraphPolygon(n, e, poly);
%     drawGraphEdges(n2, e2, 'b');
%
%   See also
%   graphs, boundedVoronoi2d, centroidalVoronoi2d_MC
%
%   Rewritten from programs found in
%   http://people.scs.fsu.edu/~burkardt/m_src/cvt/cvt.html
%
%   Reference:
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2012-02-23, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

%% Parse input arguments

% Number of germs
if isscalar(germs)
    nGerms = germs;
    germs = [];
else
    nGerms = size(germs, 1);
end

% Number of iterations
nIter = 50;

verbose = false;

keepPaths = nargout > 1;

while length(varargin) > 1
    paramName = varargin{1};
    switch lower(paramName)
        case 'verbose'
            verbose = varargin{2};
        case 'niter'
            nIter = varargin{2};

        otherwise
            error(['Unknown parameter name: ' paramName]);
    end

    varargin(1:2) = [];
end

%% Initialisations

% bounding box of polygon
bbox = polygonBounds(poly);

% init germs if needed
if isempty(germs)
    germs = generatePointsInPoly(nGerms);
end
germIters = cell(nIter, 1);

%% Iteration of the Lloyd algorithm

for i = 1:nIter
     if verbose
        disp(sprintf('Iteration: %d/%d', i, nIter)); %#ok<DSPS>
    end

    if keepPaths
        germIters{i} = germs;
    end

    % Compute Clipped Voronoi diagram of germs
    if verbose
        disp('  compute Voronoi Diagram');
    end
    [n, e, f] = boundedVoronoi2d(bbox, germs);
    [n2, e2, f2] = clipMesh2dPolygon(n, e, f, poly); %#ok<ASGLU>

    % update the position of each germ
    if verbose
        disp('  compute centroids');
    end
    for iGerm = 1:nGerms
        polygon = n2(f2{iGerm}, :);
        germs(iGerm,:) = polygonCentroid(polygon);
    end

end

%% Evenutally compute germs trajectories

if nargout > 1
    % init
    germPaths = cell(nGerms, 1);
    path = zeros(nIter+1, 2);

    % Iteration on germs
    for i = 1:nGerms

        % create path corresponding to germ
        for j = 1:nIter
            pts = germIters{j};
            path(j,:) = pts(i,:);
        end
        path(nIter+1, :) = germs(i,:);

        germPaths{i} = path;
    end
end

function pts = generatePointsInPoly(nPts)
    % extreme coordinates
    xmin = bbox(1);  xmax = bbox(2);
    ymin = bbox(3);  ymax = bbox(4);

    % compute size of box
    dx = xmax - xmin;
    dy = ymax - ymin;

    % allocate memory for result
    pts = zeros(nPts, 2);

    % iterate until all points have been sampled within the polygon
    ind = (1:nPts)';
    while ~isempty(ind)
        NI = length(ind);
        x = rand(NI, 1) * dx + xmin;
        y = rand(NI, 1) * dy + ymin;
        pts(ind, :) = [x y];

        ind = ind(~polygonContains(poly, pts(ind, :)));
    end
end

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function view = fixPhase(view, scanNum, sliceNum, shift, scale)
%
% function view = fixPhase(view, [scanNum], [sliceNum], [shift], [scale])
%
%
% Adds specified offset and scale to the phases of the specified coranal.
% Will prompt for offset & scale if not given of empty.
% The current slice will be used if sliceNum not given or empty.
%
% HISTORY:
%   2002.03.08 RFD (bob@white.stanford.edu) wrote it.

global dataTYPES;

if(~exist('scanNum','var') | isempty(scanNum))
    scanNum = getCurScan(view);
end

if(~exist('sliceNum','var') | isempty(sliceNum))
    sliceNum = viewGet(view, 'Current Slice');
end

oldShift = dataTYPES(view.curDataType).atlasParams(scanNum).phaseShift(sliceNum);
oldScale = dataTYPES(view.curDataType).atlasParams(scanNum).phaseScale(sliceNum);

oldPh = view.ph;
% undo the old shift & scale so that we can start anew
view.ph{scanNum}(:,:,sliceNum) = mod((view.ph{scanNum}(:,:,sliceNum) - oldShift)/oldScale, 2*pi);
%refreshView(view);

if(~exist('shift','var')) shift = []; end
if(~exist('scale','var')) scale = []; end
if(isempty(shift) | isempty(scale))
    if(isempty(shift))
        def{1} = num2str(oldShift);
    else
        def{1} = num2str(shift);
    end
    if(isempty(scale))
        def{2} = num2str(oldScale);
    else
        def{2} = num2str(scale);
    end

    ok = 0;
    while(~ok)
        answer = inputdlg({'Phase shift:','Scale Factor:'}, 'Fix Phase', 1, def);
        if(isempty(answer))
            % Cancel all changes
            view.ph = oldPh;
            refreshView(view);
            ok = 1;
            return;
        end
        shift = eval(answer{1});
        scale = eval(answer{2});
        view.ph{scanNum}(:,:,sliceNum) = view.ph{scanNum}(:,:,sliceNum) * scale + shift;
        view.ph{scanNum}(:,:,sliceNum) = phaseStandardize(view.ph{scanNum}(:,:,sliceNum));
        view.ph{scanNum}(:,:,sliceNum) = mod(view.ph{scanNum}(:,:,sliceNum), 2*pi);
        refreshView(view);
        ans = questdlg('Is this OK?', ...
            'Confirm Phase Adjustment', ...
            'Yes','No- try again','Cancel all changes','Yes');

        switch ans,
        case 'Yes',
            saveCorAnal(view);
            ok = 1;
        case 'No- try again',
            refreshView(view);
            def = answer;
            ok = 0;
        case 'Cancel all changes',
            view.ph = oldPh;
            refreshView(view);
            ok = -1;
            return;
        end % switch
    end
else
    view.ph{scanNum}(:,:,sliceNum) = view.ph{scanNum}(:,:,sliceNum) * scale + shift;
    view.ph{scanNum}(:,:,sliceNum) = phaseStandardize(view.ph{scanNum}(:,:,sliceNum));
    view.ph{scanNum}(:,:,sliceNum) = mod(view.ph{scanNum}(:,:,sliceNum), 2*pi);
    refreshView(view);
    saveCorAnal(view);
end
dataTYPES(view.curDataType).atlasParams(scanNum).phaseShift(sliceNum) = shift;
dataTYPES(view.curDataType).atlasParams(scanNum).phaseScale(sliceNum) = scale;
saveSession;
return;

function ph = phaseStandardize(ph)
% wrap neagative phases back up to the 2pi end

l = ph<0;
if(~isempty(l))
    ph(l) = ph(l)+2*pi;
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function box_behnken_test02 ( )

%*****************************************************************************80
%
%% BOX_BEHNKEN_TEST02 tests R8MAT_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2012
%
%  Author:
%
%    John Burkardt
%
  dim_num = 4;

  range = [ ...
    0.0, 0.0, 0.0, 0.0; ...
    1.0, 1.0, 1.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_BEHNKEN_TEST02\n' );
  fprintf ( 1, '  R8MAT_WRITE writes a Box-Behnken dataset\n' );
  fprintf ( 1, '  to a file.\n' );

  r8mat_transpose_print ( dim_num, 2, range, '  The ranges:' );

  x_num = box_behnken_size ( dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  the Box-Behnken design is of size %d\n', x_num );

  x = box_behnken ( dim_num, x_num, range );

  file_out_name = 'box_behnken_04_33.txt';

  r8mat_write ( file_out_name, dim_num, x_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data was written to the file "%s".\n', file_out_name );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function v = null(N, prefs, nullity)
%NULL   Null space of a linear CHEBOP.
%   Z = NULL(N) returns a CHEBMATRIX with orthonormal columns which span the
%   null space of the linear CHEBOP N. That is, N(Z) has negligible elements,
%   SIZE(Z, 2) is the nullity of N, and Z'*Z = I. N may contain linear boundary
%   conditions, but they will be treated as homogeneous. The nullity is
%   determined by comparing the differential order of the system and the number
%   of supplied boundary conditions.
%
%   NULL(N, PREFS) allows additional preferences to be passed via the
%   CHEBOPPREF, PREFS.
%
%   NULL(n, PREFS, K) or NULL(N, K) attempts to find K null vectors. If the
%   nullity of N is determined to be less than K then a warning is thrown. This
%   is useful in situations where the nullity is known in advance and the
%   algorithm struggles to determine it automatically.
%
%   Systems of equations are not yet supported.
%
% Example 1:
%   N = chebop(@(u) diff(u), [0, pi]);
%   V = null(N);
%   norm(N(V))
%
% Example 2:
%   N = chebop(@(x, u) 0.2*diff(u, 3) - sin(3*x).*diff(u));
%   N.rbc = 1;
%   V = null(N)
%   norm(N(V))
%   plot(V)
%
% See also CHEBOP/SVDS, CHEBOP/EIGS, NULL.

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

if ( nargin < 3 )
    nullity = [];
end
if ( nargin < 2 )
    prefs = cheboppref();
elseif ( ~isa(prefs, 'cheboppref') )
    nullity = prefs;
    prefs = cheboppref();
end

% Linearize and check whether the CHEBOP is linear:
[L, ignored, fail] = linop(N); %#ok<ASGLU>
if ( fail )
    error('CHEBFUN:CHEBOP:null:nonlinear', ...
        ['The operator appears to be nonlinear.\n', ...
         'NULL() supports only linear CHEBOP instances.']);
end

% Determine the discretization:
prefs = determineDiscretization(N, L, prefs);

% Call LINOP/NULL:
v = null(L, prefs, nullity);

% Return a CHEBFUN rather than a CHEBMATRIX for scalar problems:
if ( isa(v, 'chebmatrix') && all(size(v, 1) == 1) )
    v = chebfun(v);
end

end


text:
% Multiple plots using plotPub

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;
Vm = 10;
phi = pi/4;

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v1 = Vm*sin(th);
v2 = Vm*sin(th - phi);
v3 = Vm*sin(th - phi*2);

%% plot and settings
plt = Plot(t*1E3, v1, t*1E3, v2, t*1E3, v3);

plt.XLabel = 'Time, t (ms)'; % xlabel
plt.YLabel = 'Voltage, V (V)'; %ylabel
plt.YTick = [-10, 0, 10];
plt.YLim = [-11, 11];

% Save? comment the following line if you do not want to save
plt.export('plotMultiple.png');

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
%  Figure 3.26      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.26

clf;
zeta=.5;
k=1/zeta;
den=[1 1 1];

a=10;
num=[k/a 1];
t=0:.1:10;
y1=step(num,den,t);

a=4;
num=[k/a 1];
t=0:.1:10;
y2=step(num,den,t);

a=2;
num=[k/a 1];
t=0:.1:10;
y3=step(num,den,t);

a=1;
num=[k/a 1];
t=0:.1:10;
y4=step(num,den,t);

plot(t,y1,'-',t,y2,'-',t,y3,'-',t,y4,'-'),grid
title('Fig. 3.26 Step response with \xi = 0.5')
xlabel('\omega_n t')
ylabel('Step response of H(s)')
text(.1,.9,'\alpha=')
text(.6,.9,'1')
text(1.1,.9,'2')
text(1.5,.85,'4')
text(1.8,.8,'100')
% grid
nicegrid