text:
function h = qqPlot(varargin)
%QQPLOT Quantile-quantile plot with patch option
%
%   IOSR.STATISTICS.QQPLOT(Y) displays a quantile-quantile plot of the
%   sample quantiles of Y versus theoretical quantiles from a normal
%   distribution. If the distribution of Y is normal, the plot will be
%   close to linear.
%
%   IOSR.STATISTICS.QQPLOT(X,Y) displays a quantile-quantile plot of two
%   samples. If the samples come from the same distribution, the plot will
%   be linear.
%
%   The inputs X and Y should be numeric and have an equal number of
%   elements; every element is treated as a member of the sample.
%
%   The plot displays the sample data with the plot symbol 'x'.
%   Superimposed on the plot is a dashed straight line connecting the first
%   and third quartiles.
%
%   IOSR.STATISTICS.QQPLOT(...,MODE) allows the appearance of the plot to
%   be configured. With MODE='line' (default), the plot appears as
%   described above. With MODE='patch', the data are plotted as a patch
%   object, with the area bound by the x-distribution and the linear fit
%   shaded grey. With mode='both' the two appearances are combined.
%
%   IOSR.STATISTICS.QQPLOT(...,MODE,METHOD) and
%   IOSR.STATISTICS.QQPLOT(...,[],METHOD) allows the method for calculating
%   the quartiles, used for the fit line, to be specified. The default is
%   'R-8'. Type 'help iosr.statistics.quantile' for more information. The
%   latter form of the function call uses the default mode.
%
%   H = IOSR.STATISTICS.QQPLOT(...) returns a two- or three-element vector
%   of handles to the plotted object. The nature of the handles depends
%   upon the mode. In all cases, the first handle is to the sample data,
%   the second handle is to the fit line. With MODE='patch' or MODE='both',
%   there is third handle to the patch object.
%
%   Example
%
%      % Display Q-Q plots for the rand and randn functions
%      figure
%      subplot(2,1,1)
%      iosr.statistics.qqPlot(rand(20),'patch')
%      subplot(2,1,2)
%      h = iosr.statistics.qqPlot(randn(20),'patch');
%      set(h(3),'FaceColor','r') % change fill color
%
%   See also IOSR.STATISTICS.QUANTILE, IOSR.STATISTICS.BOXPLOT.

%   Copyright 2016 University of Surrey.

    %% determine X and Y

    IXn = cellfun(@(x) isnumeric(x) & ~isempty(x),varargin);

    switch sum(IXn)
        case 0
            error('iosr:qqPlot:noData','No input data specified')
        case 1
            % compare to normal distrbution
            Y = get_input_sample(varargin,IXn);
            p = (.5:length(Y))/length(Y);
            X = sqrt(2)*erfinv(2*p - 1);
            x_label = 'Standard normal quantiles';
            y_label = 'Sample quantiles';
        case 2
            % compare to input data distribution
            Y = get_input_sample(varargin,find(IXn,1,'last'));
            X = get_input_sample(varargin,find(IXn,1,'first'));
            assert(isequal(size(X),size(Y)), 'iosr:quantile:invalidInput', 'Input data must be the same size')
            x_label = 'X quantiles';
            y_label = 'Y quantiles';
        otherwise
            error('iosr:qqPlot:unkonwnInput','Unknown input specified')
    end

    %% determine mode and method

    % find inputs
    IXc = cellfun(@(x) ischar(x) | isempty(x),varargin);
    switch sum(IXc)
        case 0
            mode = [];
            method = [];
        case 1
            mode = varargin{IXc};
            method = [];
        case 2
            mode = varargin{find(IXc,1,'first')};
            method = varargin{find(IXc,1,'last')};
        otherwise
            error('iosr:qqPlot:unknownString','Unknown string specified')
    end

    % defaults
    if isempty(mode)
        mode = 'line';
    end
    if isempty(method)
        method = 'R-8';
    end

    %% calculate fit to first and third quartile

    % quartiles
    q1x = iosr.statistics.quantile(X,.25,[],method);
    q3x = iosr.statistics.quantile(X,.75,[],method);
    q1y = iosr.statistics.quantile(Y,.25,[],method);
    q3y = iosr.statistics.quantile(Y,.75,[],method);

    % slope
    slope = (q3y-q1y)./(q3x-q1x);
    centerx = (q1x+q3x)/2;
    centery = (q1y+q3y)/2;

    % fit
    maxx = max(X);
    minx = min(X);
    maxy = centery + slope.*(maxx - centerx);
    miny = centery - slope.*(centerx - minx);

    % lines
    X_fit = linspace(minx,maxx,length(X));
    Y_fit = linspace(miny,maxy,length(X));

    %% plot data

    hold on

    if strcmpi(mode,'patch') || strcmpi(mode,'both')
        Hp = patch([X fliplr(X_fit)],[Y fliplr(Y_fit)],[0.5 0.5 0.5]);
        set(Hp,'edgecolor','none')
    else
        Hp = NaN;
    end

    switch lower(mode)
        case 'line'
            linestyle = 'xk';
        case 'patch'
            linestyle = '-k';
        case 'both'
            linestyle = 'xk';
        otherwise
            error('iosr:qqPlot:unknownMode','Unknown mode specified')
    end

    H = plot(X,Y,linestyle,X_fit,Y_fit,'--k');

    hold off

    % axis labels
    xlabel(x_label)
    ylabel(y_label)

    box on; axis tight
    set(gca,'layer','top')

    % return handle
    if nargout>0
        h = H;
        if isobject(Hp) || ishandle(Hp)
            h = [h; Hp];
        end
    end

end

function Z = get_input_sample(z,IX)
%GET_INPUT_SAMPLE get sample, order, and convert to vector
    Z = z{IX};
    Z = sort(Z(:))';
end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function square_grid_test03 ( )

%*****************************************************************************80
%
%% SQUARE_GRID_TEST03 uses a square with different sizes in each dimension.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 August 2014
%
%  Author:
%
%    John Burkardt
%
  a = [   0.0, -2.0 ];
  b = [ +10.0, +2.0 ];
  c = [ 3, 4 ];
  ns = [ 3, 3 ];

  n = ns(1) * ns(2);

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_GRID_TEST03\n' );
  fprintf ( 1, '  Create a grid using SQUARE_GRID.\n' );
  fprintf ( 1, '  Use a different physical size in every dimension.\n' );
  fprintf ( 1, '  Number of grid points N = %d\n', n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    NS     C      A         B\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 2
    fprintf ( 1, '  %4d  %4d  %4d  %8.4f  %8.4f\n', i, ns(i), c(i), a(i), b(i) );
  end

  x = square_grid ( n, ns, a, b, c );
  r8mat_transpose_print ( 2, n, x, '  Grid points:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,fval,exitflag,output]=fminsearchbnd3(fun,x0,LB,UB,options,varargin)
% FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation
% usage: x=FMINSEARCHBND(fun,x0)
% usage: x=FMINSEARCHBND(fun,x0,LB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options)
% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...)
% usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...)
%
% arguments:
%  fun, x0, options - see the help for FMINSEARCH
%
%  LB - lower bound vector or array, must be the same size as x0
%
%       If no lower bounds exist for one of the variables, then
%       supply -inf for that variable.
%
%       If no lower bounds at all, then LB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
%  UB - upper bound vector or array, must be the same size as x0
%
%       If no upper bounds exist for one of the variables, then
%       supply +inf for that variable.
%
%       If no upper bounds at all, then UB may be left empty.
%
%       Variables may be fixed in value by setting the corresponding
%       lower and upper bounds to exactly the same value.
%
% Notes:
%
%  If options is supplied, then TolX will apply to the transformed
%  variables. All other FMINSEARCH parameters should be unaffected.
%
%  Variables which are constrained by both a lower and an upper
%  bound will use a sin transformation. Those constrained by
%  only a lower or an upper bound will use a quadratic
%  transformation, and unconstrained variables will be left alone.
%
%  Variables may be fixed by setting their respective bounds equal.
%  In this case, the problem will be reduced in size for FMINSEARCH.
%
%  The bounds are inclusive inequalities, which admit the
%  boundary values themselves, but will not permit ANY function
%  evaluations outside the bounds. These constraints are strictly
%  followed.
%
%  If your problem has an EXCLUSIVE (strict) constraint which will
%  not admit evaluation at the bound itself, then you must provide
%  a slightly offset bound. An example of this is a function which
%  contains the log of one of its parameters. If you constrain the
%  variable to have a lower bound of zero, then FMINSEARCHBND may
%  try to evaluate the function exactly at zero.
%
%
% Example usage:
% rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%
% fminsearch(rosen,[3 3])     % unconstrained
% ans =
%    1.0000    1.0000
%
% fminsearchbnd(rosen,[3 3],[2 2],[])     % constrained
% ans =
%    2.0000    4.0000
%
% See test_main.m for other examples of use.
%
%
% See also: fminsearch, fminspleas
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 4
% Release date: 7/23/06

% size checks
xsize = size(x0);
x0 = x0(:);
n=length(x0);

if (nargin<3) || isempty(LB)
  LB = repmat(-inf,n,1);
else
  LB = LB(:);
end
if (nargin<4) || isempty(UB)
  UB = repmat(inf,n,1);
else
  UB = UB(:);
end

if (n~=length(LB)) || (n~=length(UB))
  error 'x0 is incompatible in size with either LB or UB.'
end

% set default options if necessary
if (nargin<5) || isempty(options)
  options = optimset('fminsearch');
end

% stuff into a struct to pass around
params.args = varargin;
params.LB = LB;
params.UB = UB;
params.fun = fun;
params.n = n;
params.OutputFcn = [];

% 0 --> unconstrained variable
% 1 --> lower bound only
% 2 --> upper bound only
% 3 --> dual finite bounds
% 4 --> fixed variable
params.BoundClass = zeros(n,1);
for i=1:n
  k = isfinite(LB(i)) + 2*isfinite(UB(i));
  params.BoundClass(i) = k;
  if (k==3) && (LB(i)==UB(i))
    params.BoundClass(i) = 4;
  end
end

% transform starting values into their unconstrained
% surrogates. Check for infeasible starting guesses.
x0u = x0;
k=1;
for i = 1:n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      if x0(i)<=LB(i)
        % infeasible starting value. Use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(x0(i) - LB(i));
      end

      % increment k
      k=k+1;
    case 2
      % upper bound only
      if x0(i)>=UB(i)
        % infeasible starting value. use bound.
        x0u(k) = 0;
      else
        x0u(k) = sqrt(UB(i) - x0(i));
      end

      % increment k
      k=k+1;
    case 3
      % lower and upper bounds
      if x0(i)<=LB(i)
        % infeasible starting value
        x0u(k) = -pi/2;
      elseif x0(i)>=UB(i)
        % infeasible starting value
        x0u(k) = pi/2;
      else
        x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1;
        % shift by 2*pi to avoid problems at zero in fminsearch
        % otherwise, the initial simplex is vanishingly small
        x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k))));
      end

      % increment k
      k=k+1;
    case 0
      % unconstrained variable. x0u(i) is set.
      x0u(k) = x0(i);

      % increment k
      k=k+1;
    case 4
      % fixed variable. drop it before fminsearch sees it.
      % k is not incremented for this variable.
  end

end
% if any of the unknowns were fixed, then we need to shorten
% x0u now.
if k<=n
  x0u(k:n) = [];
end

% were all the variables fixed?
if isempty(x0u)
  % All variables were fixed. quit immediately, setting the
  % appropriate parameters, then return.

  % undo the variable transformations into the original space
  x = xtransform(x0u,params);

  % final reshape
  x = reshape(x,xsize);

  % stuff fval with the final value
  fval = feval(params.fun,x,params.args{:});

  % fminsearchbnd was not called
  exitflag = 0;

  output.iterations = 0;
  output.funcount = 1;
  output.algorithm = 'fminsearch';
  output.message = 'All variables were held fixed by the applied bounds';

  % return with no call at all to fminsearch
  return
end

% Check for an outputfcn. If there is any, then substitute my
% own wrapper function.
if ~isempty(options.OutputFcn)
  params.OutputFcn = options.OutputFcn;
  options.OutputFcn = @outfun_wrapper;
end

% now we can call fminsearch, but with our own
% intra-objective function.
[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params);

% undo the variable transformations into the original space
x = xtransform(xu,params);

% final reshape
x = reshape(x,xsize);

% Use a nested function as the OutputFcn wrapper
  function stop = outfun_wrapper(x,varargin);
    % we need to transform x first
    xtrans = xtransform(x,params);

    % then call the user supplied OutputFcn
    stop = params.OutputFcn(xtrans,varargin{1:(end-1)});

  end

end % mainline end

% ======================================
% ========= begin subfunctions =========
% ======================================
function fval = intrafun(x,params)
% transform variables, then call original function

% transform
xtrans = xtransform(x,params);

% and call fun
fval = feval(params.fun,xtrans,params.args{:});

end % sub function intrafun end

% ======================================
function xtrans = xtransform(x,params)
% converts unconstrained variables into their original domains

xtrans = zeros(1,params.n);
% k allows some variables to be fixed, thus dropped from the
% optimization.
k=1;
for i = 1:params.n
  switch params.BoundClass(i)
    case 1
      % lower bound only
      xtrans(i) = params.LB(i) + x(k).^2;

      k=k+1;
    case 2
      % upper bound only
      xtrans(i) = params.UB(i) - x(k).^2;

      k=k+1;
    case 3
      % lower and upper bounds
      xtrans(i) = (sin(x(k))+1)/2;
      xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i);
      % just in case of any floating point problems
      xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i)));

      k=k+1;
    case 4
      % fixed variable, bounds are equal, set it at either bound
      xtrans(i) = params.LB(i);
    case 0
      % unconstrained variable.
      xtrans(i) = x(k);

      k=k+1;
  end
end

end % sub function xtransform end

author: goGPS-Project


text:
function [str, tag, cind] = gencode(item, tag, tagctx)

% function [str, tag, cind] = gencode(item, tag, tagctx)
% Generate code to recreate a cfg_item object. This function calls the
% cfg_item specific gencode_item method. Instead of creating one large,
% deeply nested struct/cell array, it creates separate variables for each
% cfg_item.
%
% Input arguments:
% item - MATLAB variable to generate code for (the variable itself, not its
%        name)
% tag     - optional: name of the variable, i.e. what will be displayed left
%           of the '=' sign. This can also be a valid struct/cell array
%           reference, like 'x(2).y'. If not provided, inputname(1) will be
%           used.
% tagctx  - optional: variable names not to be used (e.g. keywords,
%           reserved variables). A cell array of strings.
%
% Output arguments:
% str  - cellstr containing code lines to reproduce
% tag  - name of the generated variable
% cind - index into str to the line where the variable assignment is coded
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: gencode.m 3355 2009-09-04 09:37:35Z volkmar $

rev = '$Rev: 3355 $'; %#ok

if nargin < 2
    tag = inputname(1);
end;
if nargin < 3
    tagctx = {};
end
stoptag = tag;
tropts = cfg_tropts({{}},1,inf,1,inf,true);
[str, tag, cind] = gencode_item(item, tag, tagctx, stoptag, tropts);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%EASY4     The ephemeris file is reformated by rinexe. The observation type
%               is selected and the observation files for master and rover are opened.
%               The preliminary (X,Y,Z) position of the master receiver is read from
%               the header.
%               The reading position in the file is moved to the first epoch and
%               sow is computed. sow in master and rover files are synchronized.
%               We test if all SVs have an ephemris and SVs and their
%               observations are sorted to match. Next we call the function
%               baseline which computes the components of the baseline.
%              Some results are printed in workspace and some figures are plotted.

% The code could be shortened by introducing functions. However, for a
% teaching situation I found it more illustrative to compile a script.

%RINEX version 3.03

%Kai Borre 27-07-2002
%Copyright (c) by Kai Borre
%$Revision: 1.0 $  $Date: 2002/07/27  $
% Total revison 2.0, January 13, 2016

% Read RINEX ephemerides file and convert to internal Matlab format
rinexe('log_24h.15n','eph.dat');
Eph = get_eph('eph.dat');

% Selection of observation type
ss = 'C1W'%;

% Open the master observation file
ofile1 = 'log_24h.15o';
fid1 = fopen(ofile1,'rt');

% Open the rover observation file
ofile2 = 'log_r.15o'; % log_24h.15o
fid2 = fopen(ofile2,'rt');

coo = [];
ant_delta = [];
linjer = 0;

% Gobbling the master header, and collecting useful information
while 1
    linjer = linjer +1;
    line = fgetl(fid1);
    answer = strfind(line,'END OF HEADER');
    if  ~isempty(answer), break; end;
    if (line == -1), eof = 1; break; end;

    answer = strfind(line,'REC #');
    if ~isempty(answer)
        var = textscan(fid1,'%s %d14','Delimiter','\n');
        var1 = var{1}{1};
        cox = str2double(var1(1:14));
        coy = str2double(var1(15:28));
        coz = str2double(var1(29:42));
    end

    answer = strfind(line,'ANT #');
    if ~isempty(answer)
        var = textscan(fid1,'%s %d14','Delimiter','\n');
        var1 = var{1}{1};
        ax = str2double(var1(1:14));
        ay = str2double(var1(15:28));
        az = str2double(var1(29:42));
    end

    answer = strfind(line,'SYS / # / OBS TYPES');
    if ~isempty(answer)
        tline1 = strsplit(line);
        line = fgetl(fid1);
        tline2 = strsplit(line);
        tt1 = horzcat(tline1,tline2);
        i = strcmp(tt1,ss); % if tt equals ss, i =1, else 0
        ii = find(i == 1) ;
        ii = ii-2;
        if ii > 15,  ii = ii-7, end;
        % the cell array of strings tline1 originally contains two strings
        % which describe the system and number of observation types.
        % Both tline1 and tline2 terminates with six additional strings.
        % An extra string appears at  the start of tline2; it originates
        % from concatenation of the two lines. The indexing does not
        % change even if you empty the cells. They remain as empty
        % cells and keep a place
        obs_col = ii%;
    end;
    answer = strfind(line,'INTERVAL');
    if ~isempty(answer)
        interval = strtok(line);
        int = str2double(interval);
    end;
end % end reading header

% the string arrays for the tline1 and tline2 contain an integer after the
% carrier phase observation. We must account for this by the following
% correctional table
if         strcmp(ss(2:3), '1C'), obs_col = obs_col +1;
elseif  strcmp(ss(2:3), '1W'), obs_col = obs_col+2;
elseif strcmp(ss(2:3), '2X'), obs_col = obs_col +3;
else  strcmp(ss(2:3), '2W'), obs_col = obs_col +4;
end

Pos = [];
dt1 = [];
dt2 = [];
Gdop = [];
Omc = [];
bases = [];
epoch = 0;

while  1
    epoch = epoch +1;
    time1 = 0;
    sats1 = [];
    time2 = 0;
    sats2 = [];

    % Reading first observation record in master file.

    %  We read the first line in every epoch in the master file.
    % Output is sow nd number of SVs: NoSvs1.
    [time,~] = textscan(fid1,'%s %d8','Delimiter','\n');
    tid = time{1}{1};
    year = str2double(tid(3:6));
    month = str2double(tid(8:9));
    day = str2double(tid(11:12));
    hour = str2double(tid(14:15));
    minute = str2double(tid(17:18));
    second = str2double(tid(20:29));
    % static = str2double(tid(31:32));
    NoSvs1 = str2double(tid(34:36));
    dt1  = str2double(tid(38:56));
    h = hour+minute/60+second/3600;
    jd = julday(year, month, day, h);
    [~, sec_of_week] = gps_time(jd);
    time1 = sec_of_week; % sow1

    Obs1 = zeros(NoSvs1,length(obs_col));
    % Reading NoSvs1 lines of observations in the master file
    for i = 1:NoSvs1
        [obs1,~] = textscan(fid1,'%s %d8','Delimiter','\n');
        obs1x = obs1{1}{1} ;
        obs1y = strsplit(obs1x);
        sat1 = obs1y{1} ;
        sats1(i) = str2double(sat1(2:3));
        Obs1(i,:) = str2double(obs1y(obs_col));
    end

    % Next we test if all observed SVs have an ephemeris.
    % sats1 contains the SVs as read in the observation lines.
    % The intersect command delivers Sats in sorted order!
    % Therefore we must be careful in the follwing manipulations
    Sats = intersect(sats1,Eph(1,:));

    %The command ismember does not change the sequence of entries in sats1
    lia = ismember(sats1,Sats);
    % A 0 (zero) in lia indicates that a SV has been deleted. We delete the
    % corresponding row in the observations
    sats1(lia==0) = [];
    Obs1(lia==0) = [];
    NoSV1 = length(sats1);

    % dimensions
    % time1  % 1 x 1
    % NoSV1  % 1 x 1
    % sats1 % NoSV1 x 1
    % Obs1 % NoSV1 x 1

    % the whole story repated for the rover with index 2
    % we start by gobbeling the header
    if epoch == 1
        for qq = 1:linjer+ 3
            fgetl(fid2);
        end
    end

    % We read the header line for the first epoch in the rover file.
    % It gives us sow and number of SVs.
    [time,~] = textscan(fid2,'%s %d8','Delimiter','\n');
    tid = time{1}{1};
    year = str2double(tid(3:6));
    month = str2double(tid(8:9));
    day = str2double(tid(11:12));
    hour = str2double(tid(14:15));
    minute = str2double(tid(17:18));
    second = str2double(tid(20:29));
    % static = str2double(tid(31:32));
    NoSvs2 = str2double(tid(34:36));
    dt2  = str2double(tid(38:56));
    h = hour+minute/60+second/3600;
    jd = julday(year, month, day, h);
    [~, sec_of_week] = gps_time(jd);
    time2 = sec_of_week; % sow2

    % Reading the observations in the first epch in the rover file
    Obs2 = zeros(NoSvs2,length(obs_col));
    for i = 1:NoSvs2
        [obs2,~] = textscan(fid2,'%s %d8','Delimiter','\n');
        obs2y = obs2{1}{1};
        obs2 = strsplit(obs2y);
        sat2 = obs2{1};
        sats2(i,:) = str2double(sat2(2:3));
        Obs2(i,:) = str2double(obs2(obs_col));
    end

    % Next we test if all observed sats have an ephemeris.
    % sats contains the SVs as read in the observation lines
    % The intersect command delivers Sats in sorted order!
    % Therefore we must be careful in the follwing manipulations
    Sats = intersect(sats2,Eph(1,:));

    % The command ismember does not change the sequence of entries in sats
    lia = ismember(sats2,Sats);
    % A 0 (zero) in lia indicates that a SV has been deleted. We delete the
    % corresponding row in the observations
    sats2(lia==0) = [];
    Obs2(lia==0) = [];
    NoSV2 = length(sats2);

    % dimensions
    % time2  % 1 x 1
    % NoSV2  % 1 x 1
    % sats2   % NoSV2 x 1
    % Obs2  % NoSV2 x  1

    if epoch == 1
        % Establish synchronous reading of fid1 and fid2
        % The epoch interval is int seconds
        if  time1 > time2
            for i = 1: round(time1-time2)
                 [time,~] = textscan(fid2,'%s %d8','Delimiter','\n');
                tid = time{1}{1};
                year = str2double(tid(3:6));
                month = str2double(tid(8:9));
                day = str2double(tid(11:12));
                hour = str2double(tid(14:15));
                minute = str2double(tid(17:18));
                second = str2double(tid(20:29));
                %static = str2double(tid(31:32));
                NoSvs2 = str2double(tid(34:36));
                dt2  = str2double(tid(38:56));
                h = hour+minute/60+second/3600;
                jd = julday(year, month, day, h);
                [~, sec_of_week] = gps_time(jd);
                time2 = sec_of_week; % sow2

                for qq = 1:NoSvs2
                    fgetl(fid2);
                end
            end % i
        end % if time1 ...

        if  time1 < time2
            for i = 1: round(time2-time1)
                [time, ~] = textscan(fid1,'%s %d8','Delimiter','\n');
                tid = time{1}{1};
                year = str2double(tid(3:6));
                month = str2double(tid(8:9));
                day = str2double(tid(11:12));
                hour = str2double(tid(14:15));
                minute = str2double(tid(17:18));
                second = str2double(tid(20:29));
                static = str2double(tid(31:32));
                NoSvs1 = str2double(tid(34:36));
                dt1  = str2double(tid(38:56));
                h = hour+minute/60+second/3600;
                jd = julday(year, month, day, h);
                [~, sec_of_week] = gps_time(jd);
                time1 = sec_of_week; % sow1

                for qq = 1:NoSvs1
                    fgetl(fid1);
                end
            end % i
        end % if time1
        fprintf('First Common Epoch Time1 :  %8.0f\n', time1)
        fprintf('First Common Epoch Time2 :  %8.0f\n', time2)
    end % synchronization

    % By chance sats1 and sats 2 are the same. Otherwise we need to establish
    % a matching between sats1 and sats2 and likewise for Obs1 and Obs2

    % master observations: obs1, and rover observations: obs2
    [omc,base] = baseline([cox;coy;coz],Obs1,Obs2,sats1,time1,Eph);
    Omc = [Omc, omc];
    bases = [bases base];

    [pos, el, gdop] = recpo_ls(Obs1,sats1,time1,Eph);
    Gdop = [Gdop gdop];
    Pos = [Pos pos];
    if feof(fid1) == 1, break, end;
    if feof(fid2) == 1, break; end;
end % while
fclose all;

me1 = mean(bases,2);
spread1 = std(bases,1,2);
fprintf('\nBaseline Components as Computed From %2.0f Epochs:', epoch)
fprintf('\n\nX: %12.3f  Y: %12.3f  Z: %12.3f\n\n', me1(1,1),me1(2,1),me1(3,1))

figure(1);
plot((bases(:,2:end)-bases(:,2)*ones(1,epoch-1))','linewidth',.25)
title(['Variation of Baseline Components Over ',int2str(epoch),' Epochs'],'fontsize',14)
legend('X','Y','Z')
xlabel(['Epochs ', int2str(int), ' s interval'], 'fontsize',14)
ylabel('[m]','fontsize',14)
set(gca,'fontsize',14)
legend
print -dpdf easy41

figure(2);
plot(Omc(:,2:end)')
xlabel('Epochs')
ylabel('[m]','fontsize',14)
title({'Observed minus Computed Values'; 'for All SVs in Each Epoch'},'fontsize',14)
set(gca,'fontsize',14)
print -dpdf easy42

figure(3);
plot(Gdop,'linewidth',1)
axis([1 length(Gdop) 0 5])
title('GDOP')
print -dpdf easy43

%%%%%%%%%%%%%% end easy4.m %%%%%%%%%%%%%%%

author: yandld


text:
function [GAmp,GTime]=GyTrapezoid(p)

tStart=p.tStart; % start time ms
tEnd=p.tEnd;     % end time ms
tRamp=p.tRamp;   % ramp duration time ms
sRamp=p.sRamp;   % ramp steps
GyAmp=p.GyAmp;   % Gy amplitude
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

[GAmp,GTime]=StdTrap(tStart-tRamp, ...
                     tEnd+tRamp,   ...
                     tStart,               ...
                     tEnd,                 ...
                     GyAmp,max(2,sRamp),2,max(2,sRamp));

[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% Assumes depth == 0
function point3D = proj2Dto3D(point2D, cameraParams, rotationMatrix, translationVector)
    N = size(point2D, 1);
    if N ~= 2
        point2D = point2D';
    end
    if size(point2D, 1) == 2
        point2D(3, :) = 1;
    end
    numPoints = size(point2D, 2);

    % if length(point2D) == 2
    %     point2D(3) = 1;
    % end
    % if size(point2D, 2) > 1
    %     point2D = point2D';
    % end
    K = cameraParams.IntrinsicMatrix';
    z_small = inv(rotationMatrix) * inv(K) * point2D;
    lamda = repmat(-translationVector(3), 1, numPoints) ./ z_small(3, :);

    translationVector = translationVector';
    % multiply lamda*z_small, i.e. lamda(i)*z_small(:, i)
    B = cellfun(@(x, y)(translationVector + x * y), num2cell(lamda), num2cell(z_small, 1), 'uniformoutput', false);
    point3D = cell2mat(B);
end

author: kavli-ntnu


text:
function [noisy, enhanced, clean, enhanced_wav, noisySTFT, mask, variance] = RunEnhanceNN(Data, layer, para)

output = FeatureTree2(Data, para, layer);

noisySTFT = gather(output{1}{1});
noisy = gather(output{1}{2});
enhanced = gather(output{1}{3});

enhanced_wav = abs2wav(exp(enhanced(1:257,:)/2)', angle(noisySTFT)', 400, 240);

if isempty(para.test.clean_idx)
    clean  = [];
else
    clean = gather(output{1}{para.test.clean_idx});
end

if isempty(para.test.mask_idx)
    mask = [];
else
    mask = output{1}{para.test.mask_idx};
end

if isempty(para.test.var_idx)
    variance = [];
else
    variance = output{1}{para.test.var_idx};
    if 0
        enhanced = enhanced - variance.^0.5;
        enhanced_wav = abs2wav(exp(enhanced(1:257,:)/2)', angle(noisySTFT)', 400, 240);
    end
    % [enhancedML,LL] = delta2static_ML(enhanced', variance', 2, 257);
end

end

author: singaxiong

repo: SignalGraph


text:
%--- help for mcmc/autocorrplot ---
%
%  Plots autocorrelations of a given parameter
%
%  ::
%
%     hdl=autocorrplot(obj,pname)
%     hdl=autocorrplot(obj,pname,chain_id)
%     hdl=autocorrplot(obj,pname,chain_id,order)
%
%  Args:
%
%     obj (mcmc object): mcmc object
%
%     pname (string): parameter name
%
%     chain_id (integer | {[]}): choice of chain for which to plot the
%       autocorrelation. If empty, all chains are used.
%
%     order (integer | {40}): maximum order of the autocorrelation
%
%  Returns:
%     :
%
%     - **hdl** [integer]: handle to the plot
%
%

author: jmaih


text:
function[varargout]=gulf4plot(varargin)
%GULF4PLOT  A four-panel circulation plot for the Gulf of Mexico.
%
%   GULF4PLOT(LAT,LON,CV) makes a four-panel quiver plot in the Gulf of
%   Mexico with color shading being the speed of the mean flow.  CV is the
%   complex-valued velocity U+iV of size LENGTH(LAT) x LENGTH(LON) x 4.
%
%   The plot has 2 rows and 2 columns, with a single colorbar centered
%   beneath the four plots.
%
%   GULF4PLOT(LAT,LON,CV,X) instead uses X, of the same size as CV, for the
%   color shading.
%
%   GULF4PLOT(...,CAX) uses CAX, a length 2 array, for the color axis.
%
%   GULF4PLOT(...,CAX,CLABEL,LABELS), where LABELS is a cell array of four
%   strings, prints LABELS{i} after the letter label of the ith subplot.
%
%   LAT, LON, CV, and X can also be cell arrays of four elements.  This is
%   useful if the grid size is not the same in all four plots.
%
%   Usage: h=gulf4plot(lat,lon,cv);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2020 J.M. Lilly --- type 'help jlab_license' for details

% if strcmp(varargin{1}, '--t')
%     gulf4plot_test,return
% end

tweakmap_thin=['plot([-84.5+1i*21.5 -84.5+1i*22],''k--'',''linewidth'',0.5);hold on;'...
    'plot([-86.9+1i*21.5 -84.5+1i*21.5],''k--'',''linewidth'',0.5);' setstr(10) ...
    'jfig axis|[-98.35 -81.5 18 30.75] latratio|25 topoplot|gray '...
    'eval|topoplot([],[],-5:1:-1,''0.5E'') eval|topoplot([],[],-1/2,''0.5k'') eval|topoplot([],[],-0.005,''0.5E'') ' ...
    'ticksout portrait fontsize|[11 10 10 10] '];

cax=[];
clabel=[];
labels={'','','',''};
alpha=1/80;

if length(varargin{1})==1
    alpha=varargin{1};
    varargin=varargin(2:end);
end

lat=varargin{1};
lon=varargin{2};
cv=varargin{3};
x=varargin{4};
if length(varargin)>4
    cax=varargin{5};
end
if length(varargin)>5
    clabel=varargin{6};
end
if length(varargin)>6
    labels=varargin{7};
end

if isempty(x)
    if ~iscell(cv)
        x=abs(cv);
    else
        for i=1:length(cv)
            x{i}=abs(cv{i});
        end
    end
end

if ~iscell(x)
    xcell{1}=x(:,:,1);
    xcell{2}=x(:,:,2);
    xcell{3}=x(:,:,3);
    xcell{4}=x(:,:,4);
else
    xcell=x;
end

if ~iscell(cv)&&~isempty(cv)
    cvcell{1}=cv(:,:,1);
    cvcell{2}=cv(:,:,2);
    cvcell{3}=cv(:,:,3);
    cvcell{4}=cv(:,:,4);
else
    cvcell=cv;
end

if ~iscell(lat)
    latcell{1}=lat;
    latcell{2}=lat;
    latcell{3}=lat;
    latcell{4}=lat;
else
    latcell=lat;
end

if ~iscell(lon)
    loncell{1}=lon;
    loncell{2}=lon;
    loncell{3}=lon;
    loncell{4}=lon;
else
    loncell=lon;
end

figure
subplot(3,2,1)
jpcolor(loncell{1},latcell{1},xcell{1}),
hold on,eval(tweakmap_thin),
ar=get(gca,'dataaspectratio');
[long,latg]=meshgrid(loncell{1},latcell{1});
if ~isempty(cvcell)
    quiver([long(:);-90.2],[latg(:);20.5],alpha*[vcolon(real(cvcell{1}));30]*ar(1),alpha*[vcolon(imag(cvcell{1}));0],0,'k');hold on
end
text(-98,30.155,['(a) ' labels{1}])
text(-90.2,19.95,'30 cm/s')
%--------------------------------------------------------------------------
for i=2:4
    subplot(3,2,i)
    jpcolor(loncell{i},latcell{i},xcell{i}),
    hold on,eval(tweakmap_thin),
    ar=get(gca,'dataaspectratio');
    [long,latg]=meshgrid(loncell{i},latcell{i});
    if ~isempty(cvcell)
        h=quiver(long,latg,alpha*real(cvcell{i})*ar(1),alpha*imag(cvcell{i}),0,'k');hold on
    end
    text(-98,30.155,['(' setstr(97+i-1) ') ' labels{i}])
end
%--------------------------------------------------------------------------
if ~isempty(cax)
    for i=1:4,subplot(3,2,i),caxis(cax),end
end
h1=packfig(3,2,'both');
delete(h1(5:6))
orient tall
fontsize 8 8 8 8
set(gcf,'paperposition',[0.25 0.25 8 9.5])
axes(h1(3)),set(gca,'xticklabelmode','auto')
axes(h1(4)),set(gca,'xticklabelmode','auto')

axes(h1(3))
hc=colorbar('South');
if ~isempty(clabel)
    hc.Label.String=clabel;
end
pos=hc.Position;
hc.Position=[0.32 0.34 0.4 pos(4)/2];
set(hc,'AxisLocation','out')

varargout{1}=h1;

function[]=gulf4plot_test

%reporttest('GULF4PLOT',aresame())

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function varargout = mars_struct(action, varargin)
% Multifunction function for manipulating structures
%
% To help the exposition a bit:
% 'fill' in a name, means that values empty or missing
% in one structure are fetched from another
%
% 'merge' means simply that missing fields are added, with
% values, from a second structure (but not filled if empty)
%
% Each function needs to deal with the case of empty arguments
%
% FORMAT c = mars_struct('fillafromb', a, b, fieldns, flags)
% fills structure fields empty or missing in a from those present in b
% a, b are structures
% fieldns (optional) is cell array of field names to fill from in b
% c is returned structure
% Is recursive, will fill struct fields from struct fields
% flags may contain 'f', which Force fills a from b (all non empty
% fields in b overwrite those in a)
% flags may also contain 'r', which Restricts fields to write from b, to
% those that are already present in a
%
% FORMAT [c, d] = mars_struct('split', a, b)
% split structure a into two, according to fields in b
% so that c becomes a structure which contains the fields
% in a, that are also present in b, and d contains the fields
% in a that are not present in b.  b can be a structure
% or a cell array of fieldnames
%
% FORMAT [d] = mars_struct('strip', a, b)
% strips all fields present in b from those in a,
% returning denuded structure as d. b can be a structure
% or a cell array of fieldnames.  'strip' is just 'split'
% but returning only the second argument
%
% FORMAT c = mars_struct('merge', a, b)
% merges structure a and b (fields present in b added to a)
%
% FORMAT [c,d] = mars_struct('ffillsplit', a, b)
% force fill, followed by split
% All fields from a, that are also present in b, and not empty in b,
% are replaced with the values in b; the result is returned as c
% Any fields present in a, but not present in b, are returned in d
%
% FORMAT c = mars_struct('ffillmerge', a, b)
% force fill followed by merge
% performs 'ffillsplit' on a and b, then merges a and b
% All fields present in a or b are returned in c, but
% any fields present in both, now have the value from b
%
% FORMAT [c d] = mars_struct('splitmerge', a, b)
% performs 'split' on a and b, creating c and d
% then merges c with b.
% d contains fields in a that were not present in b
% c contains fields present in both, or just in b
%
% FORMAT z = mars_struct('isthere', a, b [, c [, d ...])
% returns 1 if field named in b is present in a
% and field value is not empty.
% The call is recursive if more than two arguments are passed
% Thus with structure s = struct('one', struct('two', 3))
% mars_struct('isthere', s, 'one', 'two') returns 1
%
% FORMAT z = mars_struct('getifthere', a, b [, c [, d ...])
% returns value of field named in b from a or [] if absent
% Call is recursive, like 'isthere' above.
%
% FORMAT strs = mars_struct('celldisp', a)
% returns output like disp(a) as a cell array
% Useful for printing text description of structure
%__________________________________________________________________________

% Matthew Brett
% $Id: mars_struct.m 6623 2015-12-03 18:38:08Z guillaume $

if nargin < 1
    error('Action needed');
end
if nargin < 2
    error('Must specify structure')
end
if nargin < 3
    varargin = {varargin{:} []};
end
[a,b] = deal(varargin{1:2});

switch lower(action)
    case 'fillafromb'
        % Return for empty passed args
        if isempty(a), varargout = {b}; return, end
        if isempty(b), varargout = {a}; return, end
        if nargin < 4, fieldns = []; else fieldns = varargin{3}; end
        if isempty(fieldns)
            if ~isstruct(b), error('Need struct as 2nd argument'); end
            fieldns = fieldnames(b);
        end
        if nargin < 5, flags = ''; else flags = varargin{4}; end
        if isempty(flags), flags = ' ';end

        if ischar(fieldns), fieldns=cellstr(fieldns);end

        af = fieldnames(a)';
        bf = fieldns';

        % classify fields 0 = a~b, 1 = a&b, 2=b~a
        cf = af;
        ftype = ismember(af, bf);
        if ~any(flags == 'r')
            b_not_a = find(~ismember(bf, af));
            cf =  {cf{:} bf{b_not_a}};
            ftype = [ftype ones(1, length(b_not_a))*2];
        end

        % cope with arrays of structures
        alen = prod(size(a));
        blen = prod(size(b));
        maxlen = max(alen, blen);

        for si=1:maxlen
            ctmp = [];
            for i=1:length(cf)
                fn = cf{i};
                switch ftype(i)
                    case 0 % a~b
                        fval = getfield(a(si), fn);
                    case 1 % shared field
                        bfc = getfield(b(si), fn);
                        if isempty(getfield(a(si), fn)) || ... % a field is empty
                                (any(flags == 'f' & ~isempty(bfc)))% or force fill
                            fval = bfc;
                        else % field not empty, could be struct -> recurse
                            fval = getfield(a(si),fn);
                            if isstruct(fval) && isstruct(bfc)
                                fval = mars_struct('fillafromb',fval,bfc);
                            end
                        end
                    case 2 % b~a
                        fval = getfield(b(si), fn);
                    case 3 % no field information, see below
                        fval = [];
                end
                if isempty(ctmp)
                    ctmp = struct(fn, fval);
                else
                    ctmp = setfield(ctmp, fn, fval);
                end
            end
            c(si) = ctmp;

            if si == blen % reached end of bs, rest of b~a fields are empty
                ftype = (ftype == 2) * 3;
            elseif si == alen % end of a's rest of a~b fields are empty
                ftype = (ftype == 0) * 2 + 1;
            end

        end
        varargout = {c};

    case 'split'
        if isempty(a), varargout = {a,a}; return, end
        if isempty(b), varargout = {b,a}; return, end
        d = a;
        c = [];

        if ischar(b), b = {b};end
        if isstruct(b), b = fieldnames(b);end

        for bf = b(:)'
            if isfield(a, bf{1})
                c = setfield(c, bf{1}, getfield(a, bf{1}));
                d = rmfield(d, bf{1});
            end
        end
        varargout = {c, d};

    case 'strip'
        [c,d] = mars_struct('split', a, b);
        varargout = {d};

    case 'merge'
        if isempty(a), varargout = {b}; return, end
        if isempty(b), varargout = {a}; return, end
        c = a;

        for bf = fieldnames(b)';
            if ~isfield(a, bf{1})
                c = setfield(c, bf{1}, getfield(b, bf{1}));
            end
        end
        varargout = {c};

    case 'ffillsplit'
        if isempty(a) || isempty(b)
            % Nothing in common, return unchanged
            varargout = {a, b}; return
        end
        c = a; d = b;

        cf = fieldnames(c);
        for i=1:length(cf)
            if isfield(d, cf{i})
                dfc = getfield(d,cf{i});
                if ~isempty(dfc)
                    c = setfield(c, cf{i}, dfc);
                end
                d = rmfield(d, cf{i});
            end
        end
        varargout = {c,d};

    case 'ffillmerge'
        [a,b] = mars_struct('ffillsplit', a, b);
        varargout = {mars_struct('merge', a, b)};

    case 'splitmerge'
        [a,c] = mars_struct('split', a, b);
        varargout = {mars_struct('merge', a, b) c};

    case 'isthere'
        if isempty(a), varargout = {0}; return, end
        c = mars_struct('getifthere', varargin{:});
        varargout = {~isempty(c)};

    case 'getifthere'
        if isempty(a), varargout = {[]}; return, end
        if isempty(b), varargout = {[]}; return, end
        for v = 2:nargin-1
            b = varargin{v};
            if ~isfield(a, b)
                varargout = {[]};
                return
            end
            a = getfield(a, b);
        end
        varargout = {a};

    case 'celldisp'
        if isempty(a), varargout = {{}}; return, end
        af = fieldnames(a);
        c  = {};
        pad_len = size(char(af), 2) + 4;
        pad_str = ['%' num2str(pad_len) 's: %s'];
        for f = 1:length(af)
            d     = getfield(a, af{f});
            cls   = class(d);
            sz    = size(d);
            szstr = sprintf('%dx', size(d));
            szstr(end) = [];
            switch cls
                case 'char'
                case {'double', 'float'}
                    d = ['['  num2str(d) ']'];
                otherwise
                    d = sprintf('[%s %s]', szstr, cls);
            end
            c{f} = sprintf(pad_str, af{f}, d);
        end
        varargout = {c};

    otherwise
        error(['Suspicious action was ' action]);
end % switch

author: spm


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% Find indices of positive and negative examples
pos = find(y==1); neg = find(y==0);

% Plot examples
plot(X(pos,1), X(pos,2), 'r+', 'LineWidth', 2, 'MarkerSize', 7);
plot(X(neg,1), X(neg,2), 'ko', 'MarkerFaceColor', 'yellow', 'MarkerSize', 7);

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

hold off;

end

author: yhyap

repo: machine-learning-coursera


text:
function linpack_c_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests CCHUD and CTRSL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  p = 20;
  ldr = p;
  ldz = p;
  nz = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  Hermitian matrix\n' );
  fprintf ( 1, '  CCHUD updates a Cholesky decomposition.\n' );
  fprintf ( 1, '  CTRSL solves a triangular linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this example, we use CCHUD to solve a\n' );
  fprintf ( 1, '  least squares problem R * b = z.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of equations is P = %d\n', p );
%
%  Initialize.
%
  r(1:p,1:p) = 0.0;
  z(1:p,1:nz) = 0.0;

  for i = 1 : p
    x(i) = complex ( i, mod ( i, 2 ) );
  end
%
%  Use CCHUD to form R, Z and RHO by adding X and Y a row at a time.
%  X is a row of the least squares matrix and Y the right hand side.
%
  seed = 123456789;

  for i = 1 : p
    [ row, seed ] = c4vec_uniform_01 ( p, seed );
    y(1) = row(1:p) * transpose ( x(1:p) );
    rho(1) = 0.0;
    [ r, z, rho, c, s ] = cchud ( r, ldr, p, row, z, ldz, nz, y, rho );

  end
%
%  Generate the least squares solution, b = inverse ( R ) * Z.
%
  for j = 1 : nz

    b(1:p) = z(1:p,j);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  RHS #%d\n', j );
    fprintf ( 1, '\n' );

    for i = 1 : p
      if ( i <= 5 | p-5 < i )
        fprintf ( 1, '  %8d  (%8f  %8f)\n', i, real ( b(i) ), imag ( b(i) ) );
      end
      if ( i == 5 )
        fprintf ( 1, '  ......  ..............\n' );
      end
    end

    job = 01;

    [ b, info ] = ctrsl ( r, ldr, p, b, job );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solution vector #%d\n', j );
    fprintf ( 1, '  (Should be (1,1) (2,0), (3,1) (4,0) ...)\n' );
    fprintf ( 1, '\n' );

    for i = 1 : p
      if ( i <= 5 | p-5 < i )
        fprintf ( 1, '  %8d  (%8f  %8f)\n', i, real ( b(i) ), imag ( b(i) ) );
      end
      if ( i == 5 )
        fprintf ( 1, '  ......  ..............\n' );
      end
    end

  end

  return
end


text:
function [BstMriFile, sMri, Messages] = import_mri(iSubject, MriFile, FileFormat, isInteractive, isAutoAdjust, Comment, Labels)
% IMPORT_MRI: Import a MRI file in a subject of the Brainstorm database
%
% USAGE: [BstMriFile, sMri, Messages] = import_mri(iSubject, MriFile, FileFormat='ALL', isInteractive=0, isAutoAdjust=1, Comment=[], Labels=[])
%               BstMriFiles = import_mri(iSubject, MriFiles, ...)   % Import multiple volumes at once
%
% INPUT:
%    - iSubject  : Indice of the subject where to import the MRI
%                  If iSubject=0 : import MRI in default subject
%    - MriFile   : Full filename of the MRI to import (format is autodetected)
%                  => if not specified : file to import is asked to the user
%    - FileFormat : String, one on the file formats in in_mri
%    - isInteractive : If 1, importation will be interactive (MRI is displayed after loading)
%    - isAutoAdjust  : If isInteractive=0 and isAutoAdjust=1, relice/resample automatically without user confirmation
%    - Comment       : Comment of the output file
%    - Labels        : Labels attached to this file (cell array {Nlabels x 3}: {index, text, RGB})
% OUTPUT:
%    - BstMriFile : Full path to the new file if success, [] if error
%    - sMri       : Brainstorm MRI structure
%    - Messages   : String, messages reported by this function

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2023

%% ===== PARSE INPUTS =====
if (nargin < 3) || isempty(FileFormat)
    FileFormat = 'ALL';
end
if (nargin < 4) || isempty(isInteractive)
    isInteractive = 0;
end
if (nargin < 5) || isempty(isAutoAdjust)
    isAutoAdjust = 1;
end
if (nargin < 6) || isempty(Comment)
    Comment = [];
end
if (nargin < 7) || isempty(Labels)
    Labels = [];
end
% Initialize returned variables
BstMriFile = [];
sMri = [];
Messages = [];
% Get Protocol information
ProtocolInfo     = bst_get('ProtocolInfo');
ProtocolSubjects = bst_get('ProtocolSubjects');
% Default subject
if (iSubject == 0)
	sSubject = ProtocolSubjects.DefaultSubject;
% Normal subject
else
    sSubject = ProtocolSubjects.Subject(iSubject);
end

%% ===== SELECT MRI FILE =====
% If MRI file to load was not defined : open a dialog box to select it
if isempty(MriFile)
    % Get last used directories
    LastUsedDirs = bst_get('LastUsedDirs');
    % Get last used format
    DefaultFormats = bst_get('DefaultFormats');
    if isempty(DefaultFormats.MriIn)
        DefaultFormats.MriIn = 'ALL';
    end
    % Get MRI file
    [MriFile, FileFormat, FileFilter] = java_getfile( 'open', ...
        'Import MRI...', ...              % Window title
        LastUsedDirs.ImportAnat, ...      % Default directory
        'multiple', 'files_and_dirs', ... % Selection mode
        bst_get('FileFilters', 'mri'), ...
        DefaultFormats.MriIn);
    % If no file was selected: exit
    if isempty(MriFile)
        return
    end
    % Expand file selection (if inputs are folders)
    MriFile = file_expand_selection(FileFilter, MriFile);
    if isempty(MriFile)
        error(['No ' FileFormat ' file in the selected directories.']);
    end
    % Save default import directory
    LastUsedDirs.ImportAnat = bst_fileparts(MriFile{1});
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default import format
    DefaultFormats.MriIn = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
end

%% ===== DICOM CONVERTER =====
TmpDir = [];
if strcmpi(FileFormat, 'DICOM-SPM')
    % Convert DICOM to NII
    DicomFiles = MriFile;
    TmpDir = bst_get('BrainstormTmpDir', 0, 'dicom');
    MriFile = in_mri_dicom_spm(DicomFiles, TmpDir, isInteractive);
    if isempty(MriFile)
        return;
    end
    FileFormat = 'Nifti1';
end

%% ===== LOOP ON MULTIPLE MRI =====
if iscell(MriFile) && (length(MriFile) == 1)
    MriFile = MriFile{1};
elseif iscell(MriFile) && ~strcmpi(FileFormat, 'SPM-TPM')
    % Only allow multiple import if there is already a MRI
    if isempty(sSubject.Anatomy)
        error(['You must import the first MRI in the subject folder separately.' 10 'Please select only one volume at a time.']);
    end
    % Initialize returned values
    nFiles = length(MriFile);
    BstMriFile = cell(1, nFiles);
    sMri = cell(1, nFiles);
    % Import all volumes without supervision
    for i = 1:nFiles
        [BstMriFile{i}, sMri{i}] = import_mri(iSubject, MriFile{i}, FileFormat, isInteractive, isAutoAdjust);
    end
    % All the files are imported: exit
    return;
end

%% ===== LOAD MRI FILE =====
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'Import MRI', 'Loading MRI file...');
end
% MNI / Atlas?
isMni = ismember(FileFormat, {'ALL-MNI', 'ALL-MNI-ATLAS'});
isAtlas = ismember(FileFormat, {'ALL-ATLAS', 'ALL-MNI-ATLAS', 'SPM-TPM'});
% Load MRI
isNormalize = 0;
sMri = in_mri(MriFile, FileFormat, isInteractive && ~isMni, isNormalize);
if isempty(sMri)
    bst_progress('stop');
    return
end
% History: File name
if iscell(MriFile)
    sMri = bst_history('add', sMri, 'import', ['Import from: ' MriFile{1}]);
else
    sMri = bst_history('add', sMri, 'import', ['Import from: ' MriFile]);
end

%% ===== DELETE TEMPORARY FILES =====
if ~isempty(TmpDir)
    file_delete(TmpDir, 1, 1);
end

%% ===== GET ATLAS LABELS =====
% Try to get associated labels
if isempty(Labels) && ~iscell(MriFile)
    Labels = mri_getlabels(MriFile, sMri, isAtlas);
end
% Save labels in the file structure
if ~isempty(Labels)   % Labels were found in the input folder
    sMri.Labels = Labels;
    tagAtlas = '_volatlas';
    isAtlas = 1;
elseif isAtlas    % Volume was explicitly imported as an atlas
    tagAtlas = '_volatlas';
else
    tagAtlas = '';
end
% Get atlas comment
if isAtlas && isempty(Comment) && ~iscell(MriFile)
    [fPath, fBase, fExt] = bst_fileparts(MriFile);
    switch (fBase)
        case 'aseg'
            Comment = 'ASEG';
        case 'aparc+aseg'
            Comment = 'Desikan-Killiany';
        case 'aparc.a2009s+aseg'
            Comment = 'Destrieux';
        case {'aparc.DKTatlas+aseg', 'aparc.mapped+aseg'}  % FreeSurfer, FastSurfer
            Comment = 'DKT';
    end
end

%% ===== MANAGE MULTIPLE MRI =====
fileTag = '';
% Add new anatomy
iAnatomy = length(sSubject.Anatomy) + 1;
% If add an extra MRI: read the first one to check that they are compatible
if (iAnatomy > 1) && (isInteractive || isAutoAdjust)
    % Load the reference MRI (the first one)
    refMriFile = sSubject.Anatomy(1).FileName;
    sMriRef = in_mri_bst(refMriFile);
    % Adding an MNI volume to an existing subject
    if isMni
        sMri = mri_reslice_mni(sMri, sMriRef, isAtlas);
        isSameSize = 1;
        errMsg = '';
    % Regular coregistration options between volumes
    else
        % If some transformation where made to the intial volume: apply them to the new one ?
        if isfield(sMriRef, 'InitTransf') && ~isempty(sMriRef.InitTransf) && any(ismember(sMriRef.InitTransf(:,1), {'permute', 'flipdim'}))
            if ~isInteractive || java_dialog('confirm', ['A transformation was applied to the reference MRI.' 10 10 'Do you want to apply the same transformation to this new volume?' 10 10], 'Import MRI')
                % Apply step by step all the transformations that have been applied to the original MRI
                for it = 1:size(sMriRef.InitTransf,1)
                    ttype = sMriRef.InitTransf{it,1};
                    val   = sMriRef.InitTransf{it,2};
                    switch (ttype)
                        case 'permute'
                            sMri.Cube = permute(sMri.Cube, [val, 4]);
                            sMri.Voxsize = sMri.Voxsize(val);
                        case 'flipdim'
                            sMri.Cube = bst_flip(sMri.Cube, val(1));
                    end
                end
                % Modifying the volume disables the option "Reslice"
                isResliceDisabled = 1;
            else
                isResliceDisabled = 0;
            end
        else
            isResliceDisabled = 0;
        end

        % === ASK REGISTRATION METHOD ===
        % Get volumes dimensions
        refSize = size(sMriRef.Cube(:,:,:,1));
        newSize = size(sMri.Cube(:,:,:,1));
        isSameSize = all(refSize == newSize) && all(round(sMriRef.Voxsize(1:3) .* 1000) == round(sMri.Voxsize(1:3) .* 1000));
        % Ask what operation to perform with this MRI
        if isInteractive
            % Initialize list of options to register this new MRI with the existing one
            strOptions = 'How to register the new volume with the reference image?
';
            cellOptions = {};
            % Register with the SPM
            strOptions = [strOptions, '
- SPM:   Coregister the two volumes with SPM (requires SPM toolbox).'];
            cellOptions{end+1} = 'SPM';
            % Register with the MNI transformation
            strOptions = [strOptions, '
- MNI:   Compute the MNI transformation for both volumes (inaccurate).'];
            cellOptions{end+1} = 'MNI';
            % Skip registration
            strOptions = [strOptions, '
- Ignore:   The two volumes are already registered.'];
            cellOptions{end+1} = 'Ignore';
            % Ask user to make a choice
            RegMethod = java_dialog('question', [strOptions '

'], 'Import MRI', [], cellOptions, 'Reg+reslice');
        % In non-interactive mode: ignore if possible, or use the first option available
        else
            RegMethod = 'Ignore';
        end
        % User aborted the import
        if isempty(RegMethod)
            sMri = [];
            bst_progress('stop');
            return;
        end

        % === ASK RESLICE ===
        if isInteractive && (~strcmpi(RegMethod, 'Ignore') || ...
            (isfield(sMriRef, 'InitTransf') && ~isempty(sMriRef.InitTransf) && ismember('vox2ras', sMriRef.InitTransf(:,1)) && ...
             isfield(sMri,    'InitTransf') && ~isempty(sMri.InitTransf)    && ismember('vox2ras', sMri.InitTransf(:,1)) && ...
             ~isResliceDisabled))
            % If the volumes don't have the same size, add a warning
            if ~isSameSize
                strSizeWarn = '
The two volumes have different sizes: if you answer no here,
you will not be able to overlay them in the same figure.';
            else
                strSizeWarn = [];
            end
            % Ask to reslice
            isReslice = java_dialog('confirm', [...
                'Reslice the volume?

' ...
                'This operation rewrites the new MRI to match the alignment,
size and resolution of the original volume.' ...
                strSizeWarn ...
                '

'], 'Import MRI');
        % In non-interactive mode: never reslice
        else
            isReslice = 0;
        end

        % === REGISTRATION ===
        switch (RegMethod)
            case 'MNI'
                % Register the new MRI on the existing one using the MNI transformation (+ RESLICE)
                [sMri, errMsg, fileTag] = mri_coregister(sMri, sMriRef, 'mni', isReslice, isAtlas);
            case 'SPM'
                % Register the new MRI on the existing one using SPM + RESLICE
                [sMri, errMsg, fileTag] = mri_coregister(sMri, sMriRef, 'spm', isReslice, isAtlas);
            case 'Ignore'
                if isReslice
                    % Register the new MRI on the existing one using the transformation in the input files (files already registered)
                    [sMri, errMsg, fileTag] = mri_reslice(sMri, sMriRef, 'vox2ras', 'vox2ras', isAtlas);
                else
                    % Just copy the fiducials from the reference MRI
                    [sMri, errMsg, fileTag] = mri_coregister(sMri, sMriRef, 'vox2ras', isReslice, isAtlas);
                    % Transform error in warning
                    if ~isempty(errMsg) && ~isempty(sMri) && isSameSize && ~isReslice
                        disp(['BST> Warning: ' errMsg]);
                        errMsg = [];
                    end
                end
                % Copy the old SCS and NCS fields to the new file (only if registered)
                if isSameSize || isReslice
                    sMri.SCS = sMriRef.SCS;
                    %sMri.NCS = sMriRef.NCS;
                end
        end
    end
    % Stop in case of error
    if ~isempty(errMsg)
        if isInteractive
            bst_error(errMsg, [RegMethod ' MRI'], 0);
            sMri = [];
            bst_progress('stop');
            return;
        else
            error(errMsg);
        end
    end
end

%% ===== SAVE MRI IN BRAINSTORM FORMAT =====
% Add a Comment field in MRI structure, if it does not exist yet
if ~isempty(Comment)
    sMri.Comment = Comment;
    importedBaseName = file_standardize(Comment);
else
    if ~isfield(sMri, 'Comment') || isempty(sMri.Comment)
        sMri.Comment = 'MRI';
    end
    % Use filename as comment
    if (iAnatomy > 1) || isInteractive || ~isAutoAdjust
        [fPath, fBase, fExt] = bst_fileparts(MriFile);
        fBase = strrep(fBase, '.nii', '');
        if isMni
            sMri.Comment = file_unique(fBase, {sSubject.Anatomy.Comment});
        else
            sMri.Comment = file_unique([fBase, fileTag], {sSubject.Anatomy.Comment});
        end
    end
    % Add MNI tag
    if isMni
        if isfield(sMri, 'NCS') && isfield(sMri.NCS, 'y_method') && ~isempty(sMri.NCS.y_method)
            sMri.Comment = [sMri.Comment ' (MNI-' sMri.NCS.y_method ')'];
        elseif isfield(sMri, 'NCS') && isfield(sMri.NCS, 'y') && isfield(sMri.NCS, 'iy') && ~isempty(sMri.NCS.y) && ~isempty(sMri.NCS.iy)
            sMri.Comment = [sMri.Comment ' (MNI-nonlin)'];
        elseif isfield(sMri, 'NCS') && isfield(sMri.NCS, 'R') && isfield(sMri.NCS, 'T') && ~isempty(sMri.NCS.R) && ~isempty(sMri.NCS.T)
            sMri.Comment = [sMri.Comment ' (MNI-linear)'];
        else
            sMri.Comment = [sMri.Comment ' (MNI)'];
        end
    end
    % Get imported base name
    [tmp__, importedBaseName] = bst_fileparts(MriFile);
    importedBaseName = strrep(importedBaseName, 'subjectimage_', '');
    importedBaseName = strrep(importedBaseName, '_subjectimage', '');
    importedBaseName = strrep(importedBaseName, '.nii', '');
end

%% ===== SAVE FILE =====
% Get subject subdirectory
subjectSubDir = bst_fileparts(sSubject.FileName);
% Produce a default anatomy filename
BstMriFile = bst_fullfile(ProtocolInfo.SUBJECTS, subjectSubDir, ['subjectimage_' importedBaseName fileTag tagAtlas '.mat']);
% Make this filename unique
BstMriFile = file_unique(BstMriFile);
% Save new MRI in Brainstorm format
sMri = out_mri_bst(sMri, BstMriFile);

%% ===== REFERENCE NEW MRI IN DATABASE ======
% New anatomy structure
sSubject.Anatomy(iAnatomy) = db_template('Anatomy');
sSubject.Anatomy(iAnatomy).FileName = file_short(BstMriFile);
sSubject.Anatomy(iAnatomy).Comment  = sMri.Comment;
% Default anatomy: do not change
if isempty(sSubject.iAnatomy)
    sSubject.iAnatomy = iAnatomy;
end
% Default subject
if (iSubject == 0)
	ProtocolSubjects.DefaultSubject = sSubject;
% Normal subject
else
    ProtocolSubjects.Subject(iSubject) = sSubject;
end
bst_set('ProtocolSubjects', ProtocolSubjects);
% Save first MRI as permanent default
if (iAnatomy == 1)
    db_surface_default(iSubject, 'Anatomy', iAnatomy, 0);
end

%% ===== UPDATE GUI =====
% Refresh tree
panel_protocols('UpdateNode', 'Subject', iSubject);
panel_protocols('SelectNode', [], 'subject', iSubject, -1 );
% Save database
db_save();
% Unload MRI (if a MRI with the same name was previously loaded)
bst_memory('UnloadMri', BstMriFile);

%% ===== MRI VIEWER =====
if isInteractive
    % First MRI: Edit fiducials
    if (iAnatomy == 1)
        % MRI Visualization and selection of fiducials (in order to align surfaces/MRI)
        hFig = view_mri(BstMriFile, 'EditMri');
        drawnow;
        bst_progress('stop');
        % Wait for the MRI Viewer to be closed
        if ishandle(hFig)
            waitfor(hFig);
        end
    % Other volumes: Display registration
    else
        % If volumes are registered
        if isSameSize || isReslice
            % Open the second volume as an overlay of the first one
            hFig = view_mri(refMriFile, BstMriFile);
            % Set the amplitude threshold to 30%
            if ~isAtlas
                panel_surface('SetDataThreshold', hFig, 1, 0.3);
            end
        else
            hFig = view_mri(BstMriFile);
        end
    end
else
    if ~isProgress
        bst_progress('stop');
    end
end

author: brainstorm-tools


text:
function varargout = tucker(f)
%TUCKER   SLICE-TUCKER expansion of a CHEBFUN3 object.
%   [CORE, COLS, ROWS, TUBES] = TUCKER(F) returns the core tensor CORE and
%   the three factor quasimatrices COLS, ROWS, and TUBES in the low-rank
%   representation of a CHEBFUN3 object F. The factor quasimatrices are of
%   size Inf-by-length(F, 1), Inf-by-length(F, 2) and Inf-by-length(F, 3),
%   respectively and we have
%
%   F(x,y,z) = CORE x_1 COLS(x,:,:) x_2 ROWS(:,y,:) x_3 TUBES(:,:,z).
%
%   CORE = TUCKER(F) returns the core tensor used in the construction of F.
%
% See also CHEBFUN2/CDR.

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

if ( isempty(f) )
    varargout = cell(1, nargout);
    return
end

% Get the low rank representation for f.
fCore = f.core;
fCols = f.cols;
fRows = f.rows;
fTubes = f.tubes;

% Output:
if ( nargout <= 1 )
    varargout = {fCore};
else
    % ST decomposition
    varargout = {fCore, fCols, fRows, fTubes};
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function coords = gsp_isomap(G, dim, param)
%GSP_ISOMAP isomap
%   Usage: coords = gsp_isomap(G, dim, param);
%          coords = gsp_isomap(G, dim);
%
%   Input parameters
%         G         : Graph
%         dim       : Dimensionality of the embedding
%         param     : Structure of optional parameters
%
%   Output parameters
%         coords    : Coordinates of the embedding
%
%   This function uses the weight matrix of a graph G, in order to compute
%   a *dim* -dimensional embedding (output coordinates). The algorithm used
%   is Isomap. Warning, this function might not work if the
%   graph is not connected.
%
%   *param* is a structure with optional parameters:
%
%   * *param.tol*    : Tolerance for the spectral gap (default 1e-6).
%   * *param.kernel* : The kernel used to create the graph weight matrix:
%     + 'exp'        : exponential kernel ($e^(frac{-d_{ij}}{sigma^2})$)
%     + '1/x'        : inverse of x kernel ($frac{1}{sigma+d_{ij}}$)
%     + '1/x^2'      : inverse of x^2 kernel ($frac{1}{(sigma+d_{ij})^2}$)
%     + 'resistance' : Resistance distance.
%   * *param.k*      : Max number of nearest neighbors. If not defined, the
%
%   References: tenenbaum2000global
%
%   See also: gsp_update_coordinates gsp_laplacian_eigenmaps gsp_lle
%
%   Demo: gsp_demo_graph_embedding
%
%   References: tenenbaum2000global
%

% Authors : Dion O. E. Tzamarias
% Date    : 20/11/2015

if nargin<3
    param = struct;
end

if ~isfield(param,'kernel'), param.kernel = 'exp'; end

D = gsp_weight2distance(G,param.kernel);

% Shortest paths
d = zeros(G.N);
for ii = 1:G.N
    d(ii, :) = dijkstra(D, ii);
end
D = (d.^2);

if any(isinf(D))
    warning('Disconnected graph!');
    D(isinf(D)) = max(D(not(isinf(D)))) * 10;
end

% Perform MDS on Dijkstra Distance matrix D
B = 0.5*(repmat(sum(D,2),1,length(D))/length(D)...
    + repmat(sum(D,2).',length(D),1)/length(D) ...
    - D - repmat(sum(sum(D,1)),...
    length(D),length(D))/length(D).^2);

[coords, e] = eigs(B);
e = diag(e);
[e , ind] = sort(e, 'descend');
coords = coords(:,ind(1:dim));
coords = coords * diag(sqrt(e(1:dim)));

end

repo: gspbox


text:
classdef nntopology < nntest
  properties (TestParameter)
    topology = {'sequential', 'diamond'}
    conserveMemory = {false, true}
  end

  methods (Test)
    function testEval(test, topology, conserveMemory)
      if strcmp(test.currentDataType, 'double'), return ; end

      x = Input() ;

      switch topology
      case 'sequential'
        % a sequence of 3 layers that can be worked out by hand.

        % network input
        x_value = [1, -1] ;
        x_value = reshape(single(x_value), 1, 1, 2) ;

        % expected output
        y_value = [1, 2] ;
        y_value = reshape(single(y_value), 1, 1, 2) ;

        % expected derivative
        x_der = [-1, 0] ;
        x_der = reshape(single(x_der), 1, 1, 2) ;

        % derivative for backprop
        y_der = ones(1, 1, 2, 'single') ;

        % convs will just be the identity, or its negative
        w_value = reshape(eye(2, 'single'), 1, 1, 2, 2) ;

        y = vl_nnconv(x, Param('value', w_value), []) ;
        y = vl_nnrelu(y) ;
        y = vl_nnconv(y, Param('value', -w_value), Param('value', single([2, 2]))) ;

      case 'diamond'
        % input branches out into 2 middle layers, which are then joined
        % again. tests correct accumulation of derivatives, and math
        % operators.

        % network input
        x_value = [3; 10] ;

        % expected output
        y_value = [9; 30] ;

        % expected derivative
        x_der = [3; 3] ;

        % derivative for backprop
        y_der = ones(2, 1) ;

        % ensure they're matrix multiplies, instead of scalar
        x1 = 2 * eye(2) * x ;
        x2 = 5 * eye(2) * x ;
        y = x2 - x1 ;

      otherwise
        error('Unknown topology.')
      end

      % name layers, create net and set input
      Layer.workspaceNames() ;
      net = Net(y, 'conserveMemory', conserveMemory) ;

      % handle GPU
      if strcmp(test.currentDevice, 'gpu')
        gpuDevice(1) ;
        net.move('gpu') ;
        y_der = gpuArray(y_der) ;
      end

      % run forward and backward
      net.eval({'x', x_value}, 'normal', y_der) ;

      % check output and input derivatives
      disp('Output:') ;
      disp(squeeze(net.getValue(y))) ;
      disp('Input derivative:') ;
      disp(squeeze(net.getDer(x))) ;

      test.eq(net.getValue(y), y_value) ;
      test.eq(net.getDer(x), x_der) ;
    end
  end
end

repo: MFT


text:
function [ f ] = getEEF_Force( t )
%% This function is used to get the force acting at the endeffector of the KUKA iiwa 7 R 800.

%% Syntax:
% [ f ] = getEEF_Force( t )

%% About:
% This function is used to get the endeffector force of the KUKA iiwa 7 R 800.
% The force is measured at the media flange of the robot. The components of
% the force vector are described in the base refernce frame of the robot.

%% Arreguments:
% t: is the TCP/IP connection
% f: is 1x3 cell array. Representing the  X,Y and Z components of the force acting on the
% endeffector the forces are measured in (Newton).

% Copyright, Mohammad SAFEEA, 3rd of May 2017

theCommand='Eef_force';
fprintf(t, theCommand);
message=fgets(t);
%disp(message)
[P1,N]=getDoubleFromString(message);

for i=1:3
    f{i}=P1{i};
end

end

function [jPos,j]=getDoubleFromString(message)
n=max(max(size(message)));
j=0;
numString=[];
for i=1:n
    if message(i)=='_'
        j=j+1;
        jPos{j}=str2num(numString);
        numString=[];
    else
        numString=[numString,message(i)];
    end
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% SOM Toolbox
% Version 2.0beta, May 30 2002
%
% Copyright 1997-2000 by
% Esa Alhoniemi, Johan Himberg, Juha Parhankangas and Juha Vesanto
% Contributed files may contain copyrights of their own.
%
% SOM Toolbox comes with ABSOLUTELY NO WARRANTY; for details
% see License.txt in the program package. This is free software,
% and you are welcome to redistribute it under certain conditions;
% see License.txt for details.
%
%
% Demos
%
%            som_demo1   SOM Toolbox demo 1: basic properties
%            som_demo2   SOM Toolbox demo 2: basic usage
%            som_demo3   SOM Toolbox demo 3: visualization
%            som_demo4   SOM Toolbox demo 4: data analysis
%
% Creation of structs
%
%              som_set   create & set (& check) values to structs
%             som_info   print out information on a given struct
%      som_data_struct   create & initialize a data struct
%       som_map_struct   create & initialize a map struct
%     som_topol_struct   create & initialize a topology struct
%     som_train_struct   create & initialize a train struct
%         som_clstruct   create a cluster struct
%            som_clset   set properties in a cluster struct
%            som_clget   get stuff from a cluster struct
%
% Struct conversion and file I/O
%
%           som_vs1to2   converts a version 1.0 struct to version 2.0 struct
%           som_vs2to1   converts a version 2.0 struct to version 1.0 struct
%        som_read_data   reads a (SOM_PAK format) ASCII data file
%       som_write_data   writes a SOM_PAK format codebook file
%        som_write_cod   writes a SOM_PAK format data file
%         som_read_cod   reads a SOM_PAK format codebook file
%
% Data preprocessing
%
%        som_normalize   normalize data set
%      som_denormalize   denormalize data set
%    som_norm_variable   (de)normalize one variable
%           preprocess   preprocessing GUI
%
% Initialization and training functions
%
%             som_make   create, initialize and train a SOM
%         som_randinit   random initialization algorithm
%          som_lininit   linear initialization algorithm
%         som_seqtrain   sequential training algorithm
%       som_batchtrain   batch training algorithm
%              som_gui   SOM initialization and training GUI
%       som_prototrain   a simple version of sequential training: easy to modify
%
% Clustering algorithms
%
%           som_kmeans   k-means algorithm (was earlier kmeans)
%      kmeans_clusters   try and evaluate several k-means clusterings
%           neural_gas   neural gas vector quantization algorithm
%          som_linkage   hierarchical clustering algorithms
%        som_cllinkage   hierarchical clustering of SOM
%       som_dmatminima   local minima from distance (or U-) matrix
%     som_dmatclusters   distance (or U-) matrix based clustering
%         som_clspread   spreads clusters to unassinged map units
%           som_cldist   calculate distances between clusters
%         som_gapindex   gap validity index of clustering
%             db_index   Davies-Bouldin validity index of clustering
%
% Supervised/classification algorithms
%
%       som_supervised   supervised SOM algorithm
%                 lvq1   LVQ1 algorithm
%                 lvq3   LVQ3 algorithm
%                  knn   k-NN classification algorithm
%              knn_old   k-NN classification algorithm (old version)
%
% SOM error measures
%
%          som_quality   quantization and topographic error of SOM
%       som_distortion   SOM distortion measure
%      som_distortion3   elements of the SOM distortion measure
%
% Auxiliary functions
%
%             som_bmus   calculates BMUs for given data vectors
%         som_eucdist2   pairwise squared euclidian distances between vectors
%            som_mdist   calculates pairwise distances between vectors
%           som_divide   extract subsets of data based on map
%            som_label   give labels to map units
%        som_label2num   rcodes string data labels to interger class labels
%        som_autolabel   automatically labels the SOM based on given data
%      som_unit_coords   calculates coordinates in output space for map units
%       som_unit_dists   distances in output space between map units
%      som_unit_neighs   units in 1-neighborhood for each map unit
%     som_neighborhood   calculates neighborhood matrix for the given map
%        som_neighbors   calculates different kinds of neighborhoods
%           som_neighf   calculates neighborhood function values
%           som_select   GUI for manual selection of map units
%     som_estimate_gmm   create Gaussian mixture model on top of SOM
%  som_probability_gmm   evaluate Gaussian mixture model
%          som_ind2sub   from linear index to subscript index
%          som_sub2ind   from subscript index to linear index
%          som_ind2cod   from linear index to SOM_PAK linear index
%          som_cod2ind   from SOM_linear index to SOM_PAK linear index
%             nanstats   mean, std and median which ignore NaNs
%   som_modify_dataset   add, remove, or extract samples and components
%         som_fillnans   fill NaNs in a data set based on given SOM
%            som_stats   statistics of a data set
%           som_drmake   calculate descriptive rules for a cluster
%           som_dreval   evaluate descriptive rules for a cluster
%         som_drsignif   rule significance measures
%
% Using SOM_PAK from Matlab
%
%      som_sompaktrain   uses SOM_PAK to train a map
%           sompak_gui   GUI for using SOM_PAK from Matlab
%          sompak_init   call SOM_PAK's initialization programs from Matlab
%      sompak_init_gui   GUI for using SOM_PAK's initialization from Matlab
%    sompak_rb_control   an auxiliary function for sompak_*_gui functions.
%        sompak_sammon   call SOM_PAK's Sammon program from Matlab
%    sompak_sammon_gui   GUI for using SOM_PAK's Sammon program from Matlab
%         sompak_train   call SOM_PAK's training program from Matlab
%     sompak_train_gui   GUI for using SOM_PAK's training program from Matlab
%
% Visualization
%
%             som_show   basic visualization
%         som_show_add   add labels, hits and trajectories
%       som_show_clear   remove extra markers
%       som_recolorbar   refresh/reconfigure colorbars
%         som_show_gui   GUI for using som_show and associated functions
%             som_grid   visualization of SOM grid
%           som_cplane   component planes and U-matrices
%         som_barplane   bar chart visualization of map
%         som_pieplane   pie chart visualization of map
%        som_plotplane   plot chart visualization of map
%       som_trajectory   launches a GUI for presenting comet-trajectories
%       som_dendrogram   visualization of clustering tree
%       som_plotmatrix   pairwise scatter plots and histograms
%    som_order_cplanes   order and visualize the component planes
%           som_clplot   plots of clusters (based on cluster struct)
% som_projections_plot   projections plots (see som_projections)
%       som_stats_plot   plots of statistics (see som_stats)
%
% Auxiliary functions for visualization
%
%                 hits   calculates hits, or sum of values for each map unit
%             som_hits   calculates the response of data on the map
%             som_umat   calculates the U-matrix
%                  cca   curvilinear component analysis projection algorithm
%              pcaproj   principal component projection algorithm
%               sammon   Sammon's mapping projection algorithm
%       som_connection   connection matrix for map
%       som_vis_coords   map unit coordinates used in visualizations
%        som_colorcode   create color coding for map/2D data
%         som_bmucolor   colors of the BMUs from a given map color code
%        som_normcolor   simulate indexed colormap
%     som_clustercolor   color coding which depends on clustering structure
%      som_kmeanscolor   color coding according to k-means clustering
%     som_kmeanscolor2   a newer version of the som_kmeanscolor function
%       som_fuzzycolor   a fuzzy color coding
%         som_coloring   a SOM-based color coding
%      som_projections   calculates a default set of projections
%
% Report generation stuff
%
%     som_table_struct   creates a table struct
%     som_table_modify   modifies a table struct
%      som_table_print   print a table in various formats
%            rep_utils   various utilities for printing report elements
%      som_stats_table   a table of data set statistics
%     som_stats_report   report on data set statistics
%
% Low level routines used by visualization functions
%
%            vis_patch   defines hexagonal and rectangular patches
%    vis_som_show_data   returns UserData and subplot handles stored by som_show.m
%        vis_valuetype   used for type checks
%         vis_footnote   adds a movable text to the current figure
%          vis_trajgui   the actual GUI started by som_trajectory.m
% vis_PlaneAxisProperties   set axis properties in visualization functions
% vis_footnoteButtonDownFcn   callback function for vis_footnote.m
%     vis_planeGetArgs   converts topol struct to lattice, msize argument pair
%    vis_show_gui_comp   internal function used by som_show_gui.m
%    vis_show_gui_tool   internal function used by som_show_gui.m
%
% Other
%
%           somtoolbox   this file
%            iris.data   IRIS data set (used in demos)
%          License.txt   GNU General Public License
%        Copyright.txt   Copyright notice

author: ilarinieminen

repo: SOM-Toolbox


text:
function GUI = sessionGUI_addSession(session);
%
% Add a new session to the current mrVista study.
%
%  GUI = sessionGUI_addSession(<session=dialog>);
%
% ras, 07/2006.
mrGlobals2;

if notDefined('session')            % dialog
    p = 'Select a first session to add to this study ';
    session = sessionSelectDialog(1, p);
end

% evaluate relative to the current study's parent dir, if specified
n = GUI.settings.study;
if ~isempty(GUI.studies(n).studyDir)
    session = fullfile(eval(GUI.studies(n).studyDir), session);
else
    session = fullpath(session);
end

%%%%% add to sepcified study list:
otherSessions = GUI.studies(n).sessions;
if ~isempty(cellfind(otherSessions, session))
	ok = setdiff(1:length(otherSessions), cellfind(otherSessions,session));
	otherSessions = otherSessions(ok);
end
GUI.studies(n).sessions = [session otherSessions];

if n ~= 1
    %%%%% also add to recent sessions list:
    % first, if this session was already there, remove it:
    GUI.studies(1).sessions = setdiff(GUI.studies(1).sessions, session);

    % now, add to top of the list:
	otherSessions = GUI.studies(1).sessions;
	if ~isempty(cellfind(otherSessions, session))
		ok = setdiff(1:length(otherSessions), cellfind(otherSessions,session));
		otherSessions = otherSessions(ok);
	end

    GUI.studies(1).sessions = [session otherSessions];
end

% restrict to a max # of recent sessions
maxRecentSessions = 50;
if length(GUI.studies(1).sessions) > maxRecentSessions
    GUI.studies(1).sessions = GUI.studies(1).sessions(1:maxRecentSessions);
end

studySave(GUI.studies);

% select the new session
GUI = sessionGUI_selectSession(session);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [G]=gsp_logo()
%GSP_LOGO Initialize a graph with the GSP logo
%   Usage:  G = gsp_logo();
%
%   Output parameters:
%         G     : Graph structure.
%
%   'gsp_logo()' initializes a graph structure containing the GSP logo
%
%   Example:::
%
%          G = gsp_logo();
%          gsp_plot_graph(G);
%

% Author : Johan Paratte
% Test:

load logogsp

G.W = W;
G.coords = coords;
G.info.idx_g = idx_g;
G.info.idx_s = idx_s;
G.info.idx_p = idx_p;
%G.plotting.vertex_color = [200 136/255.0 204/255.0];
%G.plotting.edge_color = [0 136/255.0 204/255.0];
G.plotting.vertex_size = 20;
G.limits = [0 640 -400 0];
G = gsp_graph_default_parameters(G);

end


text:
function [xyz_out, indx, min_distance] = select_coordinates_near_regions(region_obj, xyz, cutoff_in_mm)
% xyz_out = select_coordinates_near_regions(region_obj, xyz, cutoff_in_mm)
%
% xyz = r x 2 matrix of XYZ mm coordinates, e.g., DB.xyz from meta-analysis
% cutoff_in_mm: xyz coordinates within this distance from any coordinate in the region object will be saved

xyz1 = cat(2, region_obj.XYZmm)'; % region object

% xyz1 : q x 3 matrix of XYZ mm coordinates
% xyz : 3 x r matrix of XYZ mm coordinates

D = dist(xyz1, xyz');

min_distance = min(D, [], 1)';

indx = min_distance < cutoff_in_mm;

xyz_out = xyz(indx, :);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function loglike = gppca_loglikelihood(Y, Qgp, Ytrain)
% loglike = gppca_loglikelihood(Y, Qgp)
%
% Returns posterior predictive loglikelihood.

% Get covariance matrices (needs the training data)
D = cols(Qgp.W);
if isempty(Qgp.CovWp) || isempty(Qgp.CovXp)
  disp('Must run GP PCA to get covariance matrices.');
  Q = vbgppcamv(Ytrain,D,Qgp.inW,Qgp.inX,Qgp.covfuncW,Qgp.logthetaW,Qgp.covfuncX, ...
                Qgp.logthetaX,'init',Qgp, 'maxiter',1, 'updatepseudox', ...
                false, 'updatepseudow', false, 'initpseudow', {Qgp.pseudoW}, ...
                'initpseudox', {Qgp.pseudoX}, 'loglikelihood', false, ...
                'updatehyper', false, 'reconstruct', false);
  Qgp.CovXp = Q.CovXp;
  Qgp.CovWp = Q.CovWp;
end

% Remove empty rows/columns
Obs = ~isnan(Y);
rowrm = (colsum(Obs) == 0);
colrm = (rowsum(Obs) == 0);
if any(rowrm)
  Qgp.inW(:,rowrm) = [];
  Y(rowrm,:) = [];
  Obs(rowrm,:) = [];
end
if any(colrm)
  Qgp.inX(:,colrm) = [];
  Y(:,colrm) = [];
  Obs(:,colrm) = [];
end

% Sample X and W
N = 10; % number of samples
W = zeros([rows(Y), D, N]);
X = zeros([D, cols(Y), N]);
for d=1:cols(Qgp.W)
  W(:,d,:) = reshape(gppredrnd(Qgp.inW, Qgp.pseudoW{d}, Qgp.Wp{d}, ...
                               Qgp.CovWp{d}, Qgp.logthetaW{d}, Qgp.covfuncW{d}, ...
                               N), [rows(W), 1, N]);
  X(d,:,:) = reshape(gppredrnd(Qgp.inX, Qgp.pseudoX{d}, Qgp.Xp{d}, ...
                               Qgp.CovXp{d}, Qgp.logthetaX{d}, Qgp.covfuncX{d}, ...
                               N), [1, cols(X), N]);

  fprintf('Sampled dimension %d/%d\n', d, D);
end

% Sample tau
tau = gamrnd(Qgp.a_tau, 1/Qgp.b_tau, [N,1]);

% Estimate predictive loglikelihood by averaging over samples
loglikes = zeros(N,1);
loglike = 0;
[I,J] = find(Obs);
for k=1:length(I)
  for n=1:N
    mu = W(I(k),:,n)*X(:,J(k),n); % reconstruct
    loglikes(n) = sum(norm_lpdf(Y(I(k),J(k)), mu, sqrt(1/tau(n))));
  end
  loglike = loglike + log( mean(exp(loglikes)) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = gppredrnd(inputpred, input, postmu, postCov, logtheta, covfunc, N)

I = eye(length(postmu));

if ~iscell(covfunc), covfunc = {covfunc}; end

% Calculate the covariance matrices
Kpp = feval(covfunc{:}, logtheta, input, input);
Kpp = Kpp + 1e-6*I;
Lp = chol(Kpp, 'lower');

Kxp = feval(covfunc{:}, logtheta, inputpred, input);

Kxx = feval(covfunc{:}, logtheta, inputpred, inputpred);

R = solve_tril(Lp, Kxp');
S = solve_triu(Lp', R);

mu = R' * solve_tril(Lp, postmu);
% $$$ szKxx = size(Kxx)
% $$$ szS = size(S)
% $$$ szKpp = size(Kpp)
% $$$ szPostCov = size(postCov)
Cov = Kxx - S' * (Kpp - postCov) * S + 1e-6*eye(length(mu));

X = zeros(length(Kxx), N);
for i=1:N
  X(:,i) = mymvnrnd(mu, Cov);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%DISNORM Trainable mapping for dissimilarity matrix normalization
%
% 	V = DISNORM(D,OPT)
% 	V = D*DISNORM([],OPT)
% 	V = D*DISNORM(OPT)
%   F = E*V
%
% INPUT
%   D 	 NxN dissimilarity matrix or dataset, which defines the norm
%   E    Matrix to be normalized, e.g. D itself
% 	OPT  'max' : maximum dissimilarity is set to 1 by global rescaling
% 		   'mean': average dissimilarity is set to 1 by global rescaling (default)
%
% OUTPUT
%   V 	Trained mapping
%   F  	Normalized dissimilarity data
%
% DEFAULT
%   OPT = 'mean'
%
% DESCRIPTION
% Operation on dissimilarity matrices, like the computation of classifiers
% in dissimilarity space, may depend on the scaling of the dissimilarities
% (a single scalar for the entire matrix). This routine computes a scaling
% for a giving matrix, e.g. a training set and applies it to other
% matrices, e.g. the same training set or based on a test set.
%
% Note that DNORM = (DISNORM*MAPEX); F = D*DNORM is equivalent to
% F = D*DISNORM(D). So the user defined DNORM can be treated as a fixed
% mapping for normalizing dissimilarity matrices by themselves.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, MAPEX

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function out = disnorm(varargin)

	argin = shiftargin(varargin,'char');
  argin = setdefaults(argin,[],'mean');
  if mapping_task(argin,'definition')
    % call like U=disnorm or U=disnorm('mean') or U=disnorm([],'mean')
    out = define_mapping(argin,'untrained');
    out = setname(out,'Disnorm');
  elseif mapping_task(argin,'training')
    % call like W=disnorm(D,'mean') or W=D*U
    [D,opt] = deal(argin{:});
    if ~isdataset(D)
      D = prdataset(D,1);
      D = setfeatlab(D,getlabels(D));
    end
    if strcmpi(opt,'mean')
      n = size(D,1);
      m = sum(sum(+D))/(n*(n-1));
    elseif strcmpi(opt,'max')
      m = max(D(:));
    else
      error('Wrong OPT')
    end
    out = prmapping(mfilename,'trained',{m},[],size(D,2),size(D,2));
  elseif mapping_task(argin,'trained execution')
    % call like E=D*W
    [D,W] = deal(argin{:});
    m = getdata(W,1);
    out = D./m;
  else
    error('Illegal input')
  end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit