text:
function [output] = ft_volumelookup(cfg, volume)

% FT_VOLUMELOOKUP can be used in to combine an anatomical or functional
% atlas with the source reconstruction results. You can use it for forward
% and reverse lookup.
%
% Given the region of interest (ROI) as anatomical or functional label, it
% looks up the locations and creates a mask (as a binary volume) based on
% the label. Given the ROI as point in the brain, it creates a sphere or
% box around that point. In these two case the function is to be used as:
%   mask = ft_volumelookup(cfg, volume)
%
% Given a binary volume that indicates a ROI or a point of interest (POI),
% it looks up the corresponding anatomical or functional labels from the
% atlas. In this case the function is to be used as:
%    labels = ft_volumelookup(cfg, volume)
%
% In both cases the input volume can be:
%   mri    is the output of FT_READ_MRI source is the output of FT_SOURCEANALYSIS
%   stat   is the output of FT_SOURCESTATISTICS
%
% The configuration options for a mask according to an atlas:
%   cfg.atlas               = string, filename of atlas to use, see FT_READ_ATLAS
%   cfg.roi                 = string or cell-array of strings, ROI from anatomical atlas
%
% The configuration options for a spherical/box mask around a POI:
%   cfg.roi                 = Nx3 vector, coordinates of the POI
%   cfg.sphere              = radius of each sphere in cm/mm dep on unit of input
%   cfg.box                 = Nx3 vector, size of each box in cm/mm dep on unit of input
%   cfg.round2nearestvoxel  = 'yes' or 'no' (default = 'no'), voxel closest to point of interest is calculated
%                             and box/sphere is centered around coordinates of that voxel
%
% The configuration options for labels from a mask:
%   cfg.atlas               = string, filename of atlas to use, see FT_READ_ATLAS
%   cfg.maskparameter       = string, field in volume to be looked up, data in field should be logical
%   cfg.minqueryrange       = number, should be odd and <= to maxqueryrange (default = 1)
%   cfg.maxqueryrange       = number, should be odd and >= to minqueryrange (default = 1)
%
% The configuration options for labels around POI:
%   cfg.output              = 'single' always outputs one label; if several POI are provided, they are considered together as describing a ROI (default)
%                             'multiple' outputs one label per POI (e.g., choose to get labels for different electrodes)
%   cfg.roi                 = Nx3 vector, coordinates of the POI
%   cfg.atlas               = string, filename of atlas to use, see FT_READ_ATLAS
%   cfg.minqueryrange       = number, should be odd and <= to maxqueryrange (default = 1)
%   cfg.maxqueryrange       = number, should be odd and >= to minqueryrange (default = 1)
%   cfg.querymethod         = 'sphere' searches voxels around the ROI in a sphere (default)
%                           = 'cube' searches voxels around the ROI in a cube
%   cfg.round2nearestvoxel  = 'yes' or 'no', voxel closest to POI is calculated (default = 'yes')
%
% The label output has a field "names", a field "count" and a field "usedqueryrange".
% To get a list of areas of the given mask you can do for instance:
%      [tmp ind] = sort(labels.count,1,'descend');
%      sel = find(tmp);
%      for j = 1:length(sel)
%        found_areas{j,1} = [num2str(labels.count(ind(j))) ': ' labels.name{ind(j)}];
%      end
% In the "found_areas" variable you can then see how many times which labels are
% found. Note that in the AFNI brick one location can have 2 labels.
%
% Dependent on the input coordinates and the coordinates of the atlas, the
% input MRI is transformed betweem MNI and Talairach-Tournoux coordinates
% See http://www.mrc-cbu.cam.ac.uk/Imaging/Common/mnispace.shtml for more details.
%
% See http://www.fieldtriptoolbox.org/template/atlas for a list of templates and
% atlasses that are included in the FieldTrip release.
%
% See also FT_READ_ATLAS, FT_SOURCEPLOT

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar volume
ft_preamble provenance volume

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
% additional checking of the input data is done further down
volume = ft_checkdata(volume, 'datatype', {'volume', 'source'}, 'feedback', 'yes', 'hasunit', 'yes', 'hascoordsys', 'yes');

% the handling of the default cfg options is done further down
cfg.minqueryrange      = ft_getopt(cfg,'minqueryrange', 1);
cfg.maxqueryrange      = ft_getopt(cfg,'maxqueryrange', 1);
cfg.output             = ft_getopt(cfg,'output', []); % in future, cfg.output could be extended to support both 'label' and 'mask'

% ensure that old and unsupported options are not being relied on by the end-user's script
% instead of specifying cfg.coordsys, the user should specify the coordsys in the data
cfg = ft_checkconfig(cfg, 'forbidden', {'units', 'coordsys', 'inputcoord', 'inputcoordsys', 'coordinates'});
cfg = ft_checkconfig(cfg, 'renamedval', {'output', 'label', 'single'});

roi2mask   = 0;
mask2label = 0;
roi2label  = 0;
if isfield(cfg, 'roi') && ~isempty(cfg.output)
  roi2label = 1;
elseif isfield(cfg, 'roi')
  roi2mask = 1;
elseif isfield(cfg, 'maskparameter')
  mask2label = 1;
else
  ft_error('you should either specify cfg.roi, or cfg.maskparameter')
end

if roi2mask
  % only for volume data
  volume = ft_checkdata(volume, 'datatype', 'volume');

  cfg.round2nearestvoxel = ft_getopt(cfg, 'round2nearestvoxel', 'no');

  isatlas = iscell(cfg.roi) || ischar(cfg.roi);
  ispoi   = isnumeric(cfg.roi);
  if isatlas+ispoi ~= 1
    ft_error('do not understand cfg.roi')
  end

  if isatlas
    ft_checkconfig(cfg, 'forbidden', {'sphere', 'box'}, 'required',  {'atlas'});
  elseif ispoi
    ft_checkconfig(cfg, 'forbidden', {'atlas'});
    if isempty(ft_getopt(cfg, 'sphere')) && isempty(ft_getopt(cfg, 'box'))
      ft_error('you should either specify cfg.sphere or cfg.box')
    end
  end

elseif mask2label || roi2label
  % convert to source representation, this is easier to work with
  volume = ft_checkdata(volume, 'datatype', 'source');
  ft_checkconfig(cfg, 'required', {'atlas'});

  if cfg.minqueryrange > cfg.maxqueryrange
    ft_error('maxqueryrange should be superior or equal to minqueryrange');
  end

  if rem(cfg.minqueryrange, 2) == 0 || rem(cfg.maxqueryrange, 2) == 0
    ft_error('incorrect query range, should be odd numbers');
  end

  cfg.round2nearestvoxel = ft_getopt(cfg, 'round2nearestvoxel', 'yes');
  cfg.querymethod        = ft_getopt(cfg, 'querymethod', 'sphere');
end

if roi2mask
  % determine location of each anatomical voxel in its own voxel coordinates
  dim = volume.dim;
  i = 1:dim(1);
  j = 1:dim(2);
  k = 1:dim(3);
  [I, J, K] = ndgrid(i, j, k);
  ijk = [I(:) J(:) K(:) ones(prod(dim),1)]';
  % determine location of each anatomical voxel in head coordinates
  xyz = volume.transform * ijk; % note that this is 4xN

  if isatlas
    if ischar(cfg.atlas)
      % assume it to represent a filename
      atlas = ft_read_atlas(cfg.atlas);
    else
      % assume cfg.atlas to be a struct
      atlas = cfg.atlas;
    end

    % determine which field(s) to use to look up the labels,
    % and whether these are boolean or indexed
    fn = fieldnames(atlas);
    isboolean = false(numel(fn),1);
    isindexed = false(numel(fn),1);
    for i=1:length(fn)
      if islogical(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
        isboolean(i) = true;
      elseif isnumeric(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
        isindexed(i) = true;
      end
    end
    if any(isindexed)
      % let the indexed prevail
      fn = fn(isindexed);
      isindexed = 1;
    elseif any(isboolean)
      % use the boolean
      fn = fn(isboolean);
      isindexed = 0;
    end

    if ischar(cfg.roi)
      cfg.roi = {cfg.roi};
    end

    if isindexed
      sel = zeros(0,2);
      for m = 1:length(fn)
        for i = 1:length(cfg.roi)
          tmp = find(strcmp(cfg.roi{i}, atlas.([fn{m},'label'])));
          sel = [sel; tmp m*ones(numel(tmp),1)];
        end
      end
      fprintf('found %d matching anatomical labels\n', size(sel,1));

      % this is to accommodate for multiple parcellations:
      % the brick refers to the parcellationname
      % the value refers to the value within the given parcellation
      brick = sel(:,2);
      value = sel(:,1);

      % convert between MNI head coordinates and TAL head coordinates
      % coordinates should be expressed compatible with the atlas
      if     strcmp(volume.coordsys, 'mni') && strcmp(atlas.coordsys, 'tal')
        xyz(1:3,:) = mni2tal(xyz(1:3,:));
      elseif strcmp(volume.coordsys, 'mni') && strcmp(atlas.coordsys, 'mni')
        % nothing to do
      elseif strcmp(volume.coordsys, 'tal') && strcmp(atlas.coordsys, 'tal')
        % nothing to do
      elseif strcmp(volume.coordsys, 'tal') && strcmp(atlas.coordsys, 'mni')
        xyz(1:3,:) = tal2mni(xyz(1:3,:));
      elseif ~strcmp(volume.coordsys, atlas.coordsys)
        ft_error('the mismatch between the coordinate system in the atlas and the coordinate system in the data cannot be resolved');
      end

      % determine location of each anatomical voxel in atlas voxel coordinates
      ijk = atlas.transform \ xyz;
      ijk = round(ijk(1:3,:))';

      inside_vol = ijk(:,1)>=1 & ijk(:,1)<=atlas.dim(1) & ...
        ijk(:,2)>=1 & ijk(:,2)<=atlas.dim(2) & ...
        ijk(:,3)>=1 & ijk(:,3)<=atlas.dim(3);
      inside_vol = find(inside_vol);

      % convert the selection inside the atlas volume into linear indices
      ind = sub2ind(atlas.dim, ijk(inside_vol,1), ijk(inside_vol,2), ijk(inside_vol,3));

      brick_val = cell(1,numel(brick));
      % search the bricks for the value of each voxel
      for i=1:numel(brick_val)
        brick_val{i} = zeros(prod(dim),1);
        brick_val{i}(inside_vol) = atlas.(fn{brick(i)})(ind);
      end

      mask = zeros(prod(dim),1);
      for i=1:numel(brick_val)
        %fprintf('constructing mask for %s\n', atlas.descr.name{sel(i)});
        mask = mask | (brick_val{i}==value(i));
      end
    else
      ft_error('support for atlases that have a probabilistic segmentationstyle is not supported yet');
      % NOTE: this may be very straightforward indeed: the mask is just the logical "or" of the specified ROIs.
    end

  elseif ispoi

    if istrue(cfg.round2nearestvoxel)
      for i=1:size(cfg.roi,1)
        cfg.roi(i,:) = poi2voi(cfg.roi(i,:), xyz);
      end
    end

    % sphere(s)
    if isfield(cfg, 'sphere')
      mask = zeros(1,prod(dim));
      for i=1:size(cfg.roi,1)
        dist = sqrt( (xyz(1,:) - cfg.roi(i,1)).^2 + (xyz(2,:) - cfg.roi(i,2)).^2 + (xyz(3,:) - cfg.roi(i,3)).^2 );
        mask = mask | (dist <= cfg.sphere(i));
      end
      % box(es)
    elseif isfield(cfg, 'box')
      mask = zeros(1, prod(dim));
      for i=1:size(cfg.roi,1)
        mask = mask | ...
          (xyz(1,:) <= (cfg.roi(i,1) + cfg.box(i,1)./2) & xyz(1,:) >= (cfg.roi(i,1) - cfg.box(i,1)./2)) & ...
          (xyz(2,:) <= (cfg.roi(i,2) + cfg.box(i,2)./2) & xyz(2,:) >= (cfg.roi(i,2) - cfg.box(i,2)./2)) & ...
          (xyz(3,:) <= (cfg.roi(i,3) + cfg.box(i,3)./2) & xyz(3,:) >= (cfg.roi(i,3) - cfg.box(i,3)./2));
      end
    end
  end

  mask = reshape(mask, dim);
  fprintf('%i voxels in mask, which is %.3f %% of total volume\n', sum(mask(:)), 100*mean(mask(:)));
  output = mask;

elseif mask2label || roi2label
  if ischar(cfg.atlas)
    % assume it to represent a filename
    atlas = ft_read_atlas(cfg.atlas);
  else
    % assume cfg.atlas to be a struct
    atlas = cfg.atlas;
  end

  % determine which field(s) to use to look up the labels,
  % and whether these are boolean or indexed
  fn = fieldnames(atlas);
  isboolean = false(numel(fn),1);
  isindexed = false(numel(fn),1);
  for i=1:length(fn)
    if islogical(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
      isboolean(i) = true;
    elseif isnumeric(atlas.(fn{i})) && isequal(size(atlas.(fn{i})), atlas.dim)
      isindexed(i) = true;
    end
  end
  if any(isindexed)
    % let the indexed prevail
    fn = fn(isindexed);
    isindexed = 1;
  elseif any(isboolean)
    % use the boolean
    fn = fn(isboolean);
    isindexed = 0;
  end

  labels.name = cell(0,1);
  for k = 1:numel(fn)
    % ensure that they are concatenated as column
    tmp = atlas.([fn{k},'label']);
    labels.name = cat(1, labels.name(:), tmp(:));
  end
  labels.name{end+1} = 'no_label_found';
  labels.count = zeros(length(labels.name),1);
  for iLab = 1:length(labels.name)
    labels.usedqueryrange{iLab} = [];
  end

  if mask2label
    sel = find(volume.(cfg.maskparameter)(:));
  elseif roi2label
    if istrue(cfg.round2nearestvoxel)
      % determine location of each anatomical voxel in head coordinates
      xyz = [volume.pos ones(size(volume.pos,1),1)]'; % note that this is 4xN
      nSel = size(cfg.roi, 1);
      for i=1:nSel
        cfg.roi(i,:) = poi2voi(cfg.roi(i,:), xyz);
      end
    end % round2nearestvoxel
    sel = zeros(size(cfg.roi,1), 1);
    for i = 1:size(cfg.roi,1)
        sel(i) = find(volume.pos(:, 1) == cfg.roi(i, 1) & volume.pos(:, 2) == cfg.roi(i, 2) & volume.pos(:, 3) == cfg.roi(i, 3));
    end
  end
  if strcmp(cfg.output, 'multiple')
      labels = repmat(labels, length(sel), 1);
  end
  for iVox = 1:length(sel)
      label = {};
      for qr = cfg.minqueryrange:2:cfg.maxqueryrange
          if isempty(label)
              label = atlas_lookup(atlas, volume.pos(sel(iVox), :), 'coordsys', volume.coordsys, 'queryrange', qr, 'method', cfg.querymethod);
              usedQR = qr;
          end
      end

      if isempty(label)
          label = {'no_label_found'};
      elseif length(label) == 1
          label = {label};
      end

      if strcmp(cfg.output, 'multiple')
          iLabOut = iVox;
      else
          iLabOut = 1;
      end
      ind_lab = [];
      for iLab = 1:length(label)
          ind_lab = find(strcmp(label{iLab}, labels(iLabOut).name));
          labels(iLabOut).count(ind_lab) = labels(iLabOut).count(ind_lab)+1; % labels.count should give the number of times a label was found within a query range
      end

      %     labels.count(ind_lab) = labels.count(ind_lab) + (1/length(ind_lab));
      %     ^this gives each label a weight depending on the number of
      %     labels found within the query range. Using this method, all labels
      %     that were found will have the same number listed for labels.count,
      %     which defeats the point of the labels.count field as I understand it
      for iFoundLab = 1:length(ind_lab)
          if isempty(labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)})
              labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)} = usedQR;
          else
              labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)} = [labels(iLabOut).usedqueryrange{ind_lab(iFoundLab)} usedQR];
          end
      end
  end %iVox

  output = labels;

end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble provenance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION point of interest to voxel of interest
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function voi = poi2voi(poi, xyz)
xmin = min(abs(xyz(1,:) - poi(1))); xcl = round(abs(xyz(1,:) - poi(1))) == round(xmin);
ymin = min(abs(xyz(2,:) - poi(2))); ycl = round(abs(xyz(2,:) - poi(2))) == round(ymin);
zmin = min(abs(xyz(3,:) - poi(3))); zcl = round(abs(xyz(3,:) - poi(3))) == round(zmin);
xyzcls = xcl + ycl + zcl; ind_voi = xyzcls == 3;
if sum(ind_voi) > 1
  fprintf('%i voxels at same distance of poi, taking first voxel\n', sum(ind_voi))
  ind_voi_temp = find(ind_voi); ind_voi_temp = ind_voi_temp(1);
  ind_voi = zeros(size(ind_voi));
  ind_voi(ind_voi_temp) = 1;
  ind_voi = logical(ind_voi);
end
voi = xyz(1:3,ind_voi);
fprintf('coordinates of voi: %.1f  %.1f %.1f\n', voi(1), voi(2), voi(3));


text:
function kern = ouKernExpandParam(kern, params)

% OUKERNEXPANDPARAM Create kernel structure from OU kernel's parameters
% (see ouKernCompute or ouKernParamInit for a more detailed description of
% the OU kernel).
% FORMAT
% DESC returns a Ornstein-Uhlenbeck kernel kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : ouKernParamInit, ouKernExtractParam, kernExpandParam
%
% COPYRIGHT : David Luengo, 2009

% KERN

kern.decay = params(1);
kern.variance = params(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%kkpad 'Pad Data with a Constant Value'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kpad.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data file'
% Integer: wsize 'Padded Width ', default: 512: 'Size of width dimension after padding'
% Integer: hsize 'Padded Height ', default: 512: 'Size of height dimension after padding'
% Integer: dsize 'Padded Depth ', default: 1: 'Size of depth dimension after padding'
% Integer: tsize 'Padded Time Size ', default: 1: 'Size of time dimension after padding'
% Integer: esize 'Padded Element Size', default: 1: 'Size of element dimension after padding'
% Integer: woff 'Width Offset ', default: 0: 'Offset of original object in width dimension'
% Integer: hoff 'Height Offset ', default: 0: 'Offset of original object in height dimension'
% Integer: doff 'Depth Offset ', default: 0: 'Offset of original object in depth dimension'
% Integer: toff 'Time Offset ', default: 0: 'Offset of original object in time dimension'
% Integer: eoff 'Element Offset ', default: 0: 'Offset of original object in elements dimension'
% Toggle: wrap 'Fill Outside Data Regions by WRAPPING', default: 0: 'fill outside data regions by WRAPPING'
% OutputFile: o 'Output', required: 'Resulting output file'
%
% Example: o = kkpad(i, {'i','';'wsize',512;'hsize',512;'dsize',1;'tsize',1;'esize',1;'woff',0;'hoff',0;'doff',0;'toff',0;'eoff',0;'wrap',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kpad - Pad Data with a Constant Value
%
%  DESCRIPTION
% "Pad Object" (kpad) performs the padding of a data object using the pad
% values defined by the Real and Imaginary padding value options.
%
%  "Size"
% The output data object will be the size specified by the final
% padded size parameters, Padded Width, Padded Height, etc (wsize,
% hsize, dsize, tsize, esize).
%
%  "Offset Values"
% The offset parameters (woff, hoff, doff, toff, eoff) define where the
% origin of the source object will lie with respect to the output object.
% Specifying an offset greater than zero can be used to generate a border
% around the original data.  Negative offset values will result in
% truncation of data lying in the region (0-woff, 0-hoff, 0-doff, 0-toff,
% 0-eoff).  For example, if an offset of (-1,-1,0,0,0) is specified, the
% object will effectively be shifted by -1 in the width and height dimensions,
% and that data truncated in the destination object.  If the offset plus
% the original data size is greater than the specified destination size,
% truncation of the source data will occur.
% If the input object is larger than the final padded size specified,
% the output data object will be truncated to the specified size.
%
%  "Pad Values"
% The Real Pad Value (real) specifies the pad value that will be assigned
% to the real component of a complex pair, and the Imaginary Pad Value
% (imag) specifies the imaginary component of the pair.
% If the input data object type is not complex, the imaginary pad value
% will be ignored.  For input data objects that are not of type double,
% the specified pad values are cast to the same type as the input data object.
%
%  "Wrapping"
% The regions outside the original data bounds can be filled also by wrapping
% data from the opposite side of the dataset by specifying the Wrap flag.
% This will override the Pad Value functionality.
%
%  "Mask Data & Validity Option"
% The validity mask corresponding to padded data will be zero if
% the user specifies that padded data added by this program be marked
% as invalid (valid FALSE).  Likewise, if the user specifies that padded
% data is valid (valid TRUE), the mask corresponding to padded data will
% be non-zero.  In the first case (invalid), if no mask exists in the
% source object, on will be created.  In the valid case, no mask will be
% created.  Validity mask data corresponding to the original source object
% data is transferred to the destination unaltered.
%
%  "Location & Time Data"
% If curvilinear or rectilinear time or location data
% exist, it is unclear how to alter the dimensionality,
% so the pad operation will fail.  If uniform location data exist
% padding occurs, but the location information is not altered in
% any way.
%
%  "Map Data"
% The pad operator does not modify map data.
%
%  "Failure Modes"
% The input object must contain value or mask data.
%
%
%
%  EXAMPLES
%
% kpad -i input_file -o output_file- Pad Data with a Constant Value
% Will create a zero padded output object, using the default
% size and zero offset.
%
% kpad -i input_file -o output_file -wsize 512 -hsize 256 -real 200 -imag 100- Pad Data with a Constant Value
% This padding operation will produce an output object with a
% width of 512 and height of 256.  If the input data is complex,
% padding for the real component will be 200, and 100 for the
% imaginary component.   If the data is not complex, the imaginary
% value is ignored.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% kpad will fail if the source object contains curvilinear or rectilinear
% location or time data.  If uniform time or location data exist, it is
% transferred, but not modified.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Eldar reconstruction

N = 3;                  % Nth order nonuniform sampling
TQ = 1;                 % Nyquist Period
L=20;                   % Resolution factor for fraction delay

T = [2*TQ 3*TQ 6*TQ];      % Decimation Periods
K = 0.5*lcm(2*T(1), 2*T(2)); K = 0.5*lcm(2*K, 2*T(3))/TQ;
capT = K*TQ; M = capT./T;
capM = lcm(M(1), M(2)); capM = lcm(capM, M(3));
excess = ceil((K-1)/capM);
maxf = K/(excess*capM+1);
TQ1 = maxf*TQ;
K1 = K/maxf;

ML = 400; % number of slices
w_c = 0.85;
NS = 100;  % Number of Sinusoids

LF = capM*2*K1+1;  %359,159,239          % min length of LF should be capM*2*K1
n = -(LF-1)/2:1:(LF-1)/2;
Hd = firpm(LF-1,[0 w_c],[0 w_c*pi],'differentiator');
delayV = (LF-1)/2;

std = [1e-6 1e-5 1e-4 1e-3 1e-2 1e-1];%5*1e-1];

serE = zeros(size(std));

MCruns = 25;
MCruns1 = 25;

for tt = 1:length(std)
    aa = 0;
    display(tt);
    for rrr = 1:MCruns1
        taus = [0 1.1+std(tt)*randn 2.2+std(tt)*randn]*TQ;
        if or(taus(2)==TQ,or(taus(2)==2*TQ,or(taus(3)==2*TQ,or(taus(2)==taus(3),taus(2)==taus(1)))))
            aa = aa+1;
            display(taus);
            continue;
        end

        tausI = sort([taus(1) taus(2) taus(3) T(1)+taus(1) T(2)+taus(2) 2*T(1)+taus(1)]);
        a = zeros(1,K);
        for p = 1:K
            a(p) = 1;
            for q = 1:K
                if q ~= p
                        a(p) = a(p)/sin(pi*(tausI(p)-tausI(q))/capT);
                end;
            end;
        end;

        for pp = 1:MCruns;
            Frq = rand(1,NS)*w_c/2;
            Amp = rand(1,NS)/(sqrt(NS)*2);
            Phi = rand(1,NS)*2*pi;
            inputN = zeros(1,ML*K1);
            input = zeros(1,ML*K);
            for k = 1:NS
              inputN = inputN + Amp(k)*sin(2*pi*Frq(k)*(0:ML*K1-1)*TQ1+Phi(k));
              input = input + Amp(k)*sin(2*pi*Frq(k)*(0:ML*K-1)*TQ+Phi(k));
            end;
        end

        y = zeros(N,ML*K1);
        for p = 1:N
            tau = taus(p)+(0:ML*M(p)-1)*T(p);
            x1 = zeros(1,ML*M(p));
            for k = 1:NS
                x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
            end;

            y1 = upsample(x1,K1);

            LFE = M(p)*capM*2*K1+1;  %359,159,239          % min length of LF should be capM*2*K1
            nE = -(LFE-1)/2:1:(LFE-1)/2;
            h = sinc((nE/K1)-(taus(p)/T(p))).*kaiser_mine1(LFE,18,-K1*(taus(p)/T(p)));
            aaa = ones(1,M(p));
            bb = ones(M(p),LFE);
            for l = 1:M(p)
                for q = 1:N
                    if q ~= p
                        aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
                        bb(l,:) = bb(l,:).*sin(pi*((nE*TQ1/M(p))-taus(q)+(l-1)*T(p))/T(q));
                    end;
                end;
                bb(l,:) = bb(l,:)/aaa(l);
            end;

            bbb = zeros(M(p),LFE);
            bbn = zeros(M(p),LFE);
            for m = 1:M(p)
                for l = 1:M(p)
                    bbb(l,:) = bb(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
                end
                bbb = sum(bbb,1);
                bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
            end;
            bbn = sum(bbn,1);
            bbn = bbn.*h/M(p);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            y1 = conv(y1,bbn);
            delay = (length(h)-1)/2;
            y1 = y1(1+delay:M(p):end-delay);
            y(p,:) = y1;
        end;
        y = real(sum(y,1));
        x = inputN;
        y = y(160:end-60);
        x = x(160:end-60);
        serE(tt) = serE(tt)+20*log10(norm(x,2)/norm(y-x,2));
    end
end

serE = serE/(MCruns*(MCruns1-aa));

plot(std,serE);
xlabel('Standard Deviation (\sigma)','fontsize',14,'fontweight','b');
ylabel('SNR in dB','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function fusion = qfuser_backoff(w,scores,scrQ,ndx,backoff_scores)
% A wrapper around qfuser_linear that substitutes backoff scores
% for any trials that are missing scores from the quality fuser.
% The backoff scores will typically be the scores from the
% non-quality fusion.
% Inputs:
%   w: The trained quality fusion weights.
%   scores: The system scores.
%   scrQ: The quality measures.
%   ndx: A Key or Ndx object indicating trials.
%   backoff_scores: A vector of scores to use if missing from
%     quality fusion output.
% Outputs:
%   fusion: A vector of scores; one for each trial.  Scores will be
%     from quality fusion unless they are missing in which case
%     they will be copied over from backoff_scores.

assert(isa(scrQ,'Quality'))
assert(isa(ndx,'Ndx')||isa(ndx,'Key'))
assert(isvector(backoff_scores))

if isa(ndx,'Ndx')
    trials = ndx.trialmask;
else
    trials = ndx.tar | ndx.non;
end

[m,n] = size(scores);
assert(n==sum(trials(:)));

% Get fused scores from quality fuser.
fusion = qfuser_linear(w,scores,scrQ,ndx);

% Check for missing scores.
assert(all(size(trials)==size(scrQ.scoremask)));
missing = trials & ~scrQ.scoremask;
missing = missing(trials(:));
nmiss = sum(missing(:));

% Use backoff scores for missing scores.
if nmiss > 0
    log_warning('Quality missing for %i of %i trials, substituting backoff scores.\n',nmiss,n);
    fusion(missing) = backoff_scores(missing);
end

end

author: nesl


text:
function [A, vA, vB, bb_rel] = crop_borders(A, bcol, padding, crop_amounts)
%CROP_BORDERS Crop the borders of an image or stack of images
%
%   [B, vA, vB, bb_rel] = crop_borders(A, bcol, [padding])
%
%IN:
%   A - HxWxCxN stack of images.
%   bcol - Cx1 background colour vector.
%   padding - scalar indicating how much padding to have in relation to
%             the cropped-image-size (0<=padding<=1). Default: 0
%   crop_amounts - 4-element vector of crop amounts: [top,right,bottom,left]
%             where NaN/Inf indicate auto-cropping, 0 means no cropping,
%             and any other value mean cropping in pixel amounts.
%
%OUT:
%   B - JxKxCxN cropped stack of images.
%   vA     - coordinates in A that contain the cropped image
%   vB     - coordinates in B where the cropped version of A is placed
%   bb_rel - relative bounding box (used for eps-cropping)

%{
% 06/03/15: Improved image cropping thanks to Oscar Hartogensis
% 08/06/15: Fixed issue #76: case of transparent figure bgcolor
% 21/02/16: Enabled specifying non-automated crop amounts
% 04/04/16: Fix per Luiz Carvalho for old Matlab releases
% 23/10/16: Fixed issue #175: there used to be a 1px minimal padding in case of crop, now removed
% 15/05/22: Fixed EPS bounding box (issue #356)
%}

    if nargin < 3
        padding = 0;
    end
    if nargin < 4
        crop_amounts = nan(1,4);  % =auto-cropping
    end
    crop_amounts(end+1:4) = NaN;  % fill missing values with NaN

    [h, w, c, n] = size(A);
    if isempty(bcol)  % case of transparent bgcolor
        bcol = A(ceil(end/2),1,:,1);
    end
    if isscalar(bcol)
        bcol = bcol(ones(c, 1));
    end

    % Crop margin from left
    if ~isfinite(crop_amounts(4))
        bail = false;
        for l = 1:w
            for a = 1:c
                if ~all(col(A(:,l,a,:)) == bcol(a))
                    bail = true;
                    break;
                end
            end
            if bail
                break;
            end
        end
    else
        l = 1 + abs(crop_amounts(4));
    end

    % Crop margin from right
    if ~isfinite(crop_amounts(2))
        bcol = A(ceil(end/2),w,:,1);
        bail = false;
        for r = w:-1:l
            for a = 1:c
                if ~all(col(A(:,r,a,:)) == bcol(a))
                    bail = true;
                    break;
                end
            end
            if bail
                break;
            end
        end
    else
        r = w - abs(crop_amounts(2));
    end

    % Crop margin from top
    if ~isfinite(crop_amounts(1))
        bcol = A(1,ceil(end/2),:,1);
        bail = false;
        for t = 1:h
            for a = 1:c
                if ~all(col(A(t,:,a,:)) == bcol(a))
                    bail = true;
                    break;
                end
            end
            if bail
                break;
            end
        end
    else
        t = 1 + abs(crop_amounts(1));
    end

    % Crop margin from bottom
    bcol = A(h,ceil(end/2),:,1);
    if ~isfinite(crop_amounts(3))
        bail = false;
        for b = h:-1:t
            for a = 1:c
                if ~all(col(A(b,:,a,:)) == bcol(a))
                    bail = true;
                    break;
                end
            end
            if bail
                break;
            end
        end
    else
        b = h - abs(crop_amounts(3));
    end

    if padding == 0  % no padding
        % Issue #175: there used to be a 1px minimal padding in case of crop, now removed
        %{
        if ~isequal([t b l r], [1 h 1 w]) % Check if we're actually croppping
            padding = 1; % Leave one boundary pixel to avoid bleeding on resize
            bcol(:) = nan;  % make the 1px padding transparent
        end
        %}
    elseif abs(padding) < 1  % pad value is a relative fraction of image size
        padding = sign(padding)*round(mean([b-t r-l])*abs(padding)); % ADJUST PADDING
    else  % pad value is in units of 1/72" points
        padding = round(padding);  % fix cases of non-integer pad value
    end

    if padding > 0  % extra padding
        % Create an empty image, containing the background color, that has the
        % cropped image size plus the padded border
        B = repmat(bcol,[(b-t)+1+padding*2,(r-l)+1+padding*2,1,n]);  % Fix per Luiz Carvalho
        % vA - coordinates in A that contain the cropped image
        vA = [t b l r];
        % vB - coordinates in B where the cropped version of A will be placed
        vB = [padding+1, (b-t)+1+padding, padding+1, (r-l)+1+padding];
        % Place the original image in the empty image
        B(vB(1):vB(2), vB(3):vB(4), :, :) = A(vA(1):vA(2), vA(3):vA(4), :, :);
        A = B;
    else  % extra cropping
        vA = [t-padding b+padding l-padding r+padding];
        A = A(vA(1):vA(2), vA(3):vA(4), :, :);
        vB = [NaN NaN NaN NaN];
    end

    % For EPS cropping, determine the relative BoundingBox - bb_rel
    bb_pixels = [l-1 h-b-1 r+1 h-t+1];  %[LowerLeftXY, UpperRightXY]
    bb_pixels(bb_pixels<0) = 0;
    bb_pixels = min(bb_pixels, [w h w h]);
    bb_rel = bb_pixels ./ [w h w h];
end

function A = col(A)
    A = A(:);
end

author: altmany


text:
function [K, sK] = heatXrbfhKernCompute(heatKern, rbfhKern, x1, x2)

% HEATXRBFHKERNCOMPUTE Cross kernel between a HEAT and a RBF kernels.
% FORMAT
% DESC computes cross kernel terms between a HEAT kernel and a RBF kernel
% for the multiple output kernel.
% ARG heatKern : the kernel structure associated with the HEAT kernel.
% ARG rbfKern : the kernel structure associated with the RBFH kernel.
% ARG x1 : inputs for which kernel is to be computed. First column represent
% the time points, while the second column represents the spatial points.
% Entries with Inf indicate missing values.
% RETURN K : block of values from kernel matrix.
% RETURN sK : unscaled kernel matrix
%
% FORMAT
% DESC computes cross kernel terms between a HEAT kernel and a RBF kernel
% for the multiple output kernel.
% ARG heatKern : the kernel structure associated with the HEAT kernel.
% ARG rbfKern : the kernel structure associated with the RBFH kernel.
% ARG x1 : row inputs for which kernel is to be computed. First column
% corresponds to time points and the second column corresponds to spatial
% points. Entries with Inf indicate missing values.
% ARG x2 : column inputs for which kernel is to be computed. First column
% corresponds to time points and the second column corresponds to spatial
% points. Entries with Inf indicate missing values.
% RETURN k : block of values from kernel matrix.
% RETURN sK : unscaled kernel matrix
%
% SEEALSO : multiKernParamInit, multiKernCompute, heatKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    x2 = x1;
end
if size(x1, 2) ~= 2 || size(x2, 2) ~= 2
    error('Input can only have two columns');
end
if (heatKern.inverseWidthTime ~= rbfhKern.inverseWidthTime) || ...
        (heatKern.inverseWidthSpace ~= rbfhKern.inverseWidthSpace)
    error('Kernels cannot be cross combined if they have different inverse widths.')
end

% Split the domain into time domain and spatial domain and account for
% missing values. If there are no missing values the computation of the
% kernel is a pointwise prodruct, otherwise it is a kronecker product.
t1 = x1(x1(:,1)~=Inf,1);
t2 = x2(x2(:,1)~=Inf,1);
s1 = x1(x1(:,2)~=Inf,2);
s2 = x2(x2(:,2)~=Inf,2);
if (length(t1) == length(s1)) && (length(t2) == length(s2))
    ut1 = unique(t1);
    ut2 = unique(t2);
    us1 = unique(s1);
    us2 = unique(s2);
    if (length(ut1)*length(us1) == length(t1)) && ...
            (length(ut2)*length(us2) == length(t2))
        t1 = ut1; s1 = us1; t2 = ut2; s2 = us2;
        isPointwise = false;
        K = zeros(length(t1)*length(s1), length(t2)*length(s2));
    else
        isPointwise = true;
        K = zeros(length(t1), length(t2));
    end
else
    isPointwise = false;
    K = zeros(length(t1)*length(s1), length(t2)*length(s2));
end

% Although this is done in heatKernExpandParam.m, we do it here again as a
% precaution.

heatKern.sim.inverseWidth = heatKern.inverseWidthTime;
rbfhKern.rbf.inverseWidth = rbfhKern.inverseWidthTime;

sigmax = sqrt(2/heatKern.inverseWidthSpace);
lengthX = heatKern.lengthX;
nterms = heatKern.nTerms;
decay = heatKern.decay;
diff = heatKern.diffusion;

% Precompute some terms
w = ((1:nterms)*(pi/lengthX))';
gamma = sqrt(-1)*w;
beta = decay + diff*(w.^2);
cK = 2/lengthX;

if heatKern.includeIC
   error('Not implemented yet')
else
    if isPointwise
        for i=1:nterms
            heatKern.sim.decay = beta(i);
            Kt = simXrbfKernCompute(heatKern.sim, rbfhKern.rbf, t1, t2);
            Ks = srbfhKernCompute(sigmax, lengthX, s1, s2, w, gamma, i);
            K = K + Kt.*Ks;
        end
    else
        for i=1:nterms
            heatKern.sim.decay = beta(i);
            Kt = simXrbfKernCompute(heatKern.sim, rbfhKern.rbf, t1, t2);
            Ks = srbfhKernCompute(sigmax, lengthX, s1, s2, w, gamma, i);
            K = K + kron(Kt,Ks);
        end
    end
    sK = cK*K;
    K = heatKern.sensitivity*sK;
end


text:
function [p,s]=v_minspane(x)
%V_MINSPANE calculate minimum spanning tree using euclidean distance [p,s]=X
%
% Inputs:  x(n,d)    d-dimensional data points, one per row
%
% Outputs: p(n-1,1)  indices of the parent of each node within the tree
%                      (data point n is the root)
%          s(n-1,1)  list of the edges in ascending order of euclidean
%                      distance. Thus the shortest edge goes from node
%                      s(1) to node p(s(1))
%
% The minimum spanning tree (or shortest spanning tree ) defines a set of
% n-1 links that interconnect n points (or nodes) with the minimum total
% length. We represent these links in the form of a tree with node n as
% the root. Each node (except node n) has a unique parent node that is
% given in the output vector p; it is possible for several nodes to share
% the same parent.
% It can be useful to know which of the links are the longest, so the output
% argument lists them in ascending order. s could be calculated directly by
% calculaing the length, l, of each link as follows:
%       l=sqrt(sum((x(1:n-1,:) - x(p,:)).^2),2);
%       [v,s]=sort(l);

%      Copyright (C) Mike Brookes 2000-2009
%      Version: $Id: v_minspane.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[np,nd]=size(x);

% first do delauny tessalation to find feasible edges

t=delaunayn(x); % nd+1 vertices per polytope
nt=size(t,1);
cn=v_choosenk(nd+1,2);
nf=size(cn,1);
ee=zeros(nt,nf,2);     % space for all the edges
ee(:,:,1)=t(:,cn(:,1));
ee(:,:,2)=t(:,cn(:,2));
ee=reshape(ee,[],2);
mk=ee(:,1)>ee(:,2);
ee(mk,:)=ee(mk,2:-1:1); % make all edges in ascending order
% ees=sparse(ee(:,1),ee(:,2),1);      % remove duplicates
[er,ec]=find(sparse(ee(:,1),ee(:,2),1));      % remove duplicates
ee=[er,ec];
ne=size(ee,1);

% now apply Kruskal shortest spanning tree algorithm

sz=sum((x(ee(:,1),:)-x(ee(:,2),:)).^2,2);     % length^2 of each edge
[vz,mz]=sort(sz);
ee=ee(mz,:);        % sort edges into ascenging length order
ts=ones(ne,1);      % size of each component
tp=zeros(np,1);     % root node links
ei=zeros(np-1,1);     % index of shortest spanning tree edges
k=0;
for i=1:ne
    i1=ee(i,1);
    j=tp(i1);
    while j         % find root node for ee(i,1)
        i1=j;
        j=tp(i1);
    end
    i2=ee(i,2);
    j=tp(i2);
    while j         % find root node for ee(i,2)
        i2=j;
        j=tp(i2);
    end
    if i1~=i2       % if they are different, merge them
        k=k+1;
        ei(k)=i;    % add to the shortest spanning tree
        if ts(i1)>ts(i2)
            tp(i2)=i1;          % make i2 a sub tree of i1
            ts(i1)=ts(i1)+ts(i2);
        else
            tp(i1)=i2;          % make i1 a sub tree of i2
            ts(i2)=ts(i1)+ts(i2);
        end
    end
end
ee=ee(ei,:);        % refine the edges to include only those in the minimum spanning tree

% now arrange as a tree with point np as the head

eet=sparse(ee(:,1),ee(:,2),1:np-1,np,np);
p=zeros(np-1,1);              % points to parent node
s=zeros(np-1,1);              % sorted index
chn=np;                                 % start with the root node of the tree
[rf,cf]=find(eet(:,chn));               % find any nodes that connect to it
[rg,cg]=find(eet(chn,:));
while ~isempty(rf) || ~isempty(rg)
    p(rf)=chn(cf);                      % set the parents
    rcf=rf(:)+np*(chn(cf(:))-1);
    s(eet(rcf))=rf(:);
    eet(rcf)=0;                         % delete the edges linking them to their parents
    p(cg)=chn(rg);
    rcg=chn(rg(:))+np*(cg(:)-1);
    s(eet(rcg))=cg(:);
    eet(rcg)=0;
    chn=[rf(:); cg(:)];                 % now search for their children
    [rf,cf]=find(eet(:,chn));
    [rg,cg]=find(eet(chn,:));
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function area = polygonArea(poly, varargin)
%POLYGONAREA Compute the signed area of a polygon.
%
%   A = polygonArea(POINTS);
%   Compute area of a polygon defined by POINTS. POINTS is a N-by-2 array
%   of double containing coordinates of vertices.
%
%   Vertices of the polygon are supposed to be oriented Counter-Clockwise
%   (CCW). In this case, the signed area is positive.
%   If vertices are oriented Clockwise (CW), the signed area is negative.
%
%   If polygon is self-crossing, the result is undefined.
%
%   Examples
%     % compute area of a simple shape
%     poly = [10 10;30 10;30 20;10 20];
%     area = polygonArea(poly)
%     area =
%         200
%
%     % compute area of CW polygon
%     area2 = polygonArea(poly(end:-1:1, :))
%     area2 =
%         -200
%
%     % Computes area of a paper hen
%     x = [0 10 20  0 -10 -20 -10 -10  0];
%     y = [0  0 10 10  20  10  10  0 -10];
%     poly = [x' y'];
%     area = polygonArea(poly)
%     area =
%        400
%
%     % Area of unit square with 25% hole
%     pccw = [0 0; 1 0; 1 1; 0 1];
%     pcw = pccw([1 4 3 2], :) * .5 + .25;
%     polygonArea ([pccw; nan(1,2); pcw])
%     ans =
%        0.75
%
%   References
%   algo adapted from P. Bourke web page
%   http://paulbourke.net/geometry/polygonmesh/
%
%   See also:
%   polygons2d, polygonCentroid, polygonSecondAreaMoments, triangleArea
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 05/05/2004.
%

%   HISTORY
%   25/04/2005: add support for multiple polygons
%   12/10/2007: update doc

%% Process special cases

% in case of polygon sets, computes the sum of polygon areas
if iscell(poly)
    area = 0;
    for i = 1:length(poly)
        area = area + polygonArea(poly{i});
    end
    return;
end

% check there are enough points
if size(poly, 1) < 2
    area = 0;
    return;
end

% case of polygons with holes -> computes the sum of areas
if any(isnan(poly))
    area = sum(polygonArea(splitPolygons(poly)));
    return;
end

%% Process single polygons or single rings

% extract coordinates
if nargin == 1
    % polygon given as N-by-2 array
    px = poly(:, 1);
    py = poly(:, 2);

elseif nargin == 2
    % poylgon given as two N-by-1 arrays
    px = poly;
    py = varargin{1};
end

% indices of next vertices
N = length(px);
iNext = [2:N 1];

% compute area (vectorized version)
area = sum(px .* py(iNext) - px(iNext) .* py) / 2;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function V = load_images(params)
  sequence_name = params.sequences_name(params.current_sequence).name;
  path = fullfile(params.sequences_path,sequence_name,params.sequences_format);
  list = dir(fullfile(path,params.sequences_ext));
  names = char({list.name}');
  frames = size(names,1);
  V = [];
  for i = 1:frames
    if(mod(i,10) == 0) fprintf('%d ',i); end
    I = imread(fullfile(path,names(i,:)));
    I = im2double(I);
    V(:,:,:,i) = I;
  end
  disp(i);
end

author: andrewssobral


text:
classdef PerformanceTestsSuite < handle

    methods

        function obj = PerformanceTestsSuite()
            path = './tests/Source/PerformanceTests/PerformanceTests.m';
            suite = matlab.unittest.TestSuite.fromFile(path, 'Tag','Performance');
            results = suite.run;
            table(results)
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function w = ymdf_to_weekday_common ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_WEEKDAY_COMMON returns the weekday of a Common YMDF date.
%
%  Discussion:
%
%    The "common" calendar is meant to be the calendar which is Julian up to
%    day JED = 2299160, and Gregorian from day JED = 2299161 and after.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, integer W, is the week day number of the date, with
%    1 for Sunday, through 7 for Saturday.
%
  jed = ymdf_to_jed_common ( y, m, d, f );

  [ w, f2 ] = jed_to_weekday ( jed );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function calcFeatureImportanceRF_Old(pathRF,nPerms,seed)

startpath = pwd;

cd(pathRF), load('training'), load('testing')
nameOutcomes = fieldnames(training.outcomes); nOutcomes = numel(nameOutcomes);
fSetNames = fieldnames(training.textures.(nameOutcomes{1})); nFset = numel(fSetNames);

nInst = numel(testing.outcomes.(nameOutcomes{1})); % Number of testing patients
perms = zeros(nInst,nPerms);
rng(seed)
for p = 1:nPerms
    perms(:,p) = datasample(1:nInst,nInst,'Replace',false)';
end

for o = 1:nOutcomes
    outcome = testing.outcomes.(nameOutcomes{o});
    for f = 1:nFset
        results = load(['testResults_',fSetNames{f},'_',nameOutcomes{o}]); results = struct2cell(results); results = results{1};
        AUC_baseline = results.AUC;
        RF = load(['RF_',fSetNames{f},'_',nameOutcomes{o}]); RF = struct2cell(RF); RF = RF{1};
        indClinic = training.clinical.bestAdd.(nameOutcomes{o}).(fSetNames{f});
        text = testing.textures.(nameOutcomes{o}).(fSetNames{f}); nText = size(text,2);
        tableTest = [text,testing.clinical.table(:,indClinic)];
        varNames = tableTest.Properties.VariableNames'; nVar = numel(varNames);
        percentAUCdecrease = zeros(nVar,1);
        for v = 1:nVar
            tableTemp = tableTest;
            percentAUCdecrease_temp  = 0;
            for p = 1:nPerms
                tableTemp(:,v) = tableTest(perms(:,p),v);
                [prob] = predictRF(tableTemp,RF);
                [AUC,~,~,~] = calcPerformMetrics(prob,outcome,0.5);
                percentAUCdecrease_temp = percentAUCdecrease_temp + (AUC - AUC_baseline)/AUC_baseline;
            end
            percentAUCdecrease_temp = percentAUCdecrease_temp/nPerms;
            percentAUCdecrease(v) = percentAUCdecrease_temp;
        end
        [percentAUCdecrease,ind] = sort(percentAUCdecrease,'descend');
        variableImportance.(nameOutcomes{o}).(fSetNames{f}).percentAUCdecrease = percentAUCdecrease;
        varNames = varNames(ind);
        variableImportance.(nameOutcomes{o}).(fSetNames{f}).varNames = varNames;
    end
end
save('testingVariableImportance','variableImportance')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function X = ctransposefactor(X)
if ~isempty(X.midfactors)
    for i = 1:length(X.midfactors)
        temp = X.leftfactors{i};
        X.leftfactors{i} = X.rightfactors{i}';
        X.rightfactors{i} = temp';
        X.midfactors{i} = X.midfactors{i}';
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [y] = mci_approach_gen (P,M,U)
% Approach to limit model
% FORMAT [y] = mci_approach_gen (P,M,U)
%
% P         parameters
% M,U       as usual
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_approach_gen.m 6548 2015-09-11 12:39:47Z will $

V=exp(P(1));
tau=exp(P(2));
t=U.X;

y=-60+V*(1-exp(-t/tau));

author: spm


text:
function pass = test_get( pref )
% Test GET.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e2 * pref.cheb2Prefs.chebfun2eps;
j = 1;

f = diskfun(@(x,y) 1 + sin(pi*x.*y) + sin(pi*x));
g = diskfun(@(x,y) cos(1-x.*y));
F = diskfunv(f,g);
Fc = F.components;
G = F.';
pass(j) = norm( Fc{1} - f ) < tol; j = j + 1;
pass(j) = norm( Fc{2} - g ) < tol; j = j+1;
pass(j) = norm( 2 - F.nComponents ) < tol; j = j + 1;
pass(j) = F.isTransposed==0; j = j + 1;
pass(j) = G.isTransposed==1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function title = p15_title ( )

%*****************************************************************************80
%
%% P15_TITLE returns the title for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'log(x) / ( 1 + 100 x^2 )';

  return
end


text:
function out = run_sqrs(ecg,HRVparams,rs)
%out = run_sqrs(ecg,s,rs)
%	OVERVIEW:
%       QRS onset detector
%   INPUT:
%       ecg = single channel of ECG in digital values
%       s   = settings struct
%       rs  = resampling option; 1 = resample original signal, 0 = don't
%   OUTPUT:
%       out = sample points of onset of each QRS complex
%   DEPENDENCIES & LIBRARIES:
%   REFERENCE:
%	REPO:
%       https://github.com/cliffordlab/hrv_toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Converted from sqrs.c at:
%       http://www.physionet.org/physiotools/wfdb/app/sqrs.c
%       to  sqrs.m (Matlab)     by J. Perry,    July        2006
%           sqrs.c	(C)         by	G.B. Moody  27 October  1990
%
%	Last revised by John:  25 February 2006
%	Last revised by Gari:  07 February 2007
%
%   03-06-2017
%   Edited by Adriana Vest
%   Now requires settings struct HRVparams, which defines debug mode and
%   sampling frequency, and resampling option rs.
%   LICENSE AND COPYRIGHT:
%       See Below
%
% if debug>0, then plot final data (default; debug=0)
%
% This algorithm is sensitive to high frequency noise, so you
% might want to use an FIR band-pass filter first.
%
% Avaialble under the GPL (see below)
%-------------------------------------------------------------------------------
% sqrs: Single-channel QRS detector
% Copyright (C) 1990-2006 George B. Moody
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
% 02111-1307, USA.
%
% You may contact the author by e-mail (george@mit.edu) or postal
% mail (MIT Room E25-505A, Cambridge, MA 02139 USA).  For updates to
% this software, please visit PhysioNet (http://www.physionet.org/).
% _________________________________________________________________
%
% The detector algorithm is based on example 10 in the WFDB
% Programmer's Guide, which in turn is based on a Pascal program
% written by W.A.H. Engelse and C. Zeelenberg, "A single scan
% algorithm for QRS-detection and feature extraction", Computers in
% Cardiology 6:37-42 (1979).  `sqrs' does not include the feature
% extraction capability of the Pascal program.  The output of `sqrs'
% is an annotation file (with annotator name `qrs') in which all
% detected beats are labelled normal;  the annotation file may also
% contain `artifact' annotations at locations which `sqrs' believes
% are noise-corrupted.
%
% 'sqrs' has been optimized for adult human ECGs. For other ECGs,
% it may be necessary to experiment with the input sampling
% frequency and the time constants indicated below.
%
% This program is provided as an example only, and is not intended
% for any clinical application.  At the time the algorithm was
% originally published, its performance was typical of
% state-of-the-art QRS detectors.  Recent designs, particularly
% those that can analyze two or more input signals, may exhibit
% significantly better performance.
%%

if nargin<2
    error('Incorrect number of input arguments: must provide ECG signal and HRV paramters struct')
end
if nargin < 3
    rs = 1;
end

%debug = s.debug;
debug = 0;
fs = HRVparams.Fs;

if ~rs
    %% 1. Use Fs of Input
    time = 0;
    now = 10;

    %freq = 256;
    freq = fs;
    ms160 = ceil(0.16*freq);
    ms200 = ceil(0.2*freq);
    s2 = ceil(2*freq);
    scmin = 500;
    % number of ADC units corresponding to scmin microvolts
    % scmin = muvadu(signal, scmin);
    scmax = 10 * scmin;
    slopecrit = 10 * scmin;
    maxslope = 0;
    nslope = 0;
    out = [];

    while (now < length(ecg)) %   && (to == 0)
        filter = [1 4 6 4 1 -1 -4 -6 -4 -1] * ecg((now-9):now);
        if (mod(time, s2) == 0)
	    % Adjust slope
            if (nslope == 0)
	      slopecrit = max(slopecrit - slopecrit/16, scmin);
            elseif (nslope >= 5)
	      slopecrit = min(slopecrit + slopecrit/16, scmax);
            end
        end
        if (nslope == 0  && abs(filter) > slopecrit)
            nslope = nslope + 1;
	    maxtime = ms160;
	    if (filter > 0)
	      sign = 1;
	    else
	      sign = -1;
	    end
            qtime = time;
        end
        if (nslope ~= 0)
            if (filter * sign < -slopecrit)
                sign = -sign;
                nslope = nslope + 1;
                if (nslope > 4)
                  maxtime = ms200;
                else
                  maxtime = ms160;
                end
            elseif (filter * sign > slopecrit &&  abs(filter) > maxslope)
                maxslope = abs(filter);
            end

            if (maxtime < 0)
                if (2 <= nslope && nslope <= 4)
                    slopecrit = slopecrit + ((maxslope/4) - slopecrit)/8;
                    if (slopecrit < scmin)
                        slopecrit = scmin;
                    elseif (slopecrit > scmax)
                        slopecrit = scmax;
                    end
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = NORMAL;
                    time = 0;
                elseif (nslope >= 5)
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = ARFCT;
                end
                nslope = 0;
            end
            maxtime = maxtime - 1;
        end
	time = time + 1;
	now = now + 1;
    end

    out=out-1; % adjust for 1 sample offset problem.

    if debug > 0
       plot(ecg,'b');
       hold on;
       plot(out,ecg(out),'m*');
    end
elseif rs
    %% 2. Resample at 256 Hz
    % These time constants may need adjustment for pediatric or
    % small mammal ECGs.
    ecg_rs = resample(ecg,256,fs);

    time = 0;
    now = 10;

    freq = 256;
    ms160 = ceil(0.16*freq);
    ms200 = ceil(0.2*freq);
    s2 = ceil(2*freq);
    %scmin = 500;
    % ANV changed this to work with real units
    scmin = 500;
    % number of ADC units corresponding to scmin microvolts
    % scmin = muvadu(ecg_rs, scmin);
    scmax = 10 * scmin;
    slopecrit = 10 * scmin;
    maxslope = 0;
    nslope = 0;
    out = [];

    while (now < length(ecg_rs)) %   && (to == 0)
        filter = [1 4 6 4 1 -1 -4 -6 -4 -1] * ecg_rs((now-9):now);
        if (mod(time, s2) == 0)
            % Adjust slope
            if (nslope == 0)
                slopecrit = max(slopecrit - slopecrit/16, scmin);
            elseif (nslope >= 5)
                slopecrit = min(slopecrit + slopecrit/16, scmax);
            end
        end
        if (nslope == 0  && abs(filter) > slopecrit)
            nslope = nslope + 1;
            maxtime = ms160;
            if (filter > 0)
                sign = 1;
            else
                sign = -1;
            end
            qtime = time;
        end
        if (nslope ~= 0)
            if (filter * sign < -slopecrit)
                sign = -sign;
                nslope = nslope + 1;
                if (nslope > 4)
                    maxtime = ms200;
                else
                  maxtime = ms160;
                end
            elseif (filter * sign > slopecrit &&  abs(filter) > maxslope)
                maxslope = abs(filter);
            end
            if (maxtime < 0)
                if (2 <= nslope && nslope <= 4)
                    slopecrit = slopecrit + ((maxslope/4) - slopecrit)/8;
                    if (slopecrit < scmin)
                        slopecrit = scmin;
                    elseif (slopecrit > scmax)
                        slopecrit = scmax;
                    end
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = NORMAL;
                    time = 0;
                elseif (nslope >= 5)
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = ARFCT;
                end
                nslope = 0;
            end
            maxtime = maxtime - 1;
        end
	time = time + 1;
	now = now + 1;
    end

    out=out-1; % adjust for 1 sample offset problem.

    if debug > 0
       plot(ecg_rs,'b');
       hold on;
       plot(out,ecg_rs(out),'m*');
    end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function ft_defaults

% FT_DEFAULTS (ending with "s") sets some general settings in the global variable
% ft_default (without the "s") and takes care of the required path settings. You can
% call this function in your startup.m script. This function is also called at the
% begin of all FieldTrip functions.
%
% The global configuration defaults are stored in the global "ft_default" structure.
% The ft_checkconfig function that is called by many FieldTrip functions will merge
% these global configuration defaults with the cfg ctructure that you pass to
% the FieldTrip function that you are calling.
%
% The global options and their default values are
%   ft_default.checkconfig       = string, can be 'pedantic', 'loose', 'silent' (default = 'loose')
%   ft_default.checkpath         = string, can be 'pedantic', 'once', 'no' (default = 'pedantic')
%   ft_default.checksize         = number in bytes, can be inf (default = 1e5)
%   ft_default.showlogo          = string, can be 'yes' or 'no' (default = 'yes')
%   ft_default.showcallinfo      = string, can be 'yes' or 'no' (default = 'yes')
%   ft_default.trackcallinfo     = string, can be 'yes' or 'no' (default = 'yes')
%   ft_default.trackusage        = false, or string with salt for one-way encryption of identifying information (by default this is enabled and an automatic salt is created)
%   ft_default.trackdatainfo     = string, can be 'yes' or 'no' (default = 'no')
%   ft_default.outputfilepresent = string, can be 'keep', 'overwrite', 'error' (default = 'overwrite')
%   ft_default.debug             = string, can be 'display', 'displayonerror', 'displayonsuccess', 'save', 'saveonerror', saveonsuccess' or 'no' (default = 'no')
%   ft_default.toolbox.signal    = string, can be 'compat' or 'matlab' (default is automatic, see below)
%   ft_default.toolbox.stats     = string, can be 'compat' or 'matlab' (default is automatic, see below)
%   ft_default.toolbox.images    = string, can be 'compat' or 'matlab' (default is automatic, see below)
%   ft_default.reproducescript   = string, directory to which the script and intermediate data are written (default = [])
%
% If you want to overrule these default settings, you can add something like this in your startup.m script
%   ft_defaults
%   global ft_default
%   ft_default.option1 = value1
%   ft_default.option2 = value2
%
% The toolbox option for signal, stats and images allows you to specify whether you
% want to use the original version from MathWorks or a compatible drop-in to be used.
% When you use the Radboud University license server, i.e. at the Donders, the
% default is 'compat'. This has the advantage that you do not need a license for
% these toolboxes; we do not have that many licenses and parallel computations on our
% Donders compute cluster would otherwise use all licenses. In all other cases, the
% default is 'matlab' when the toolbox is available, and 'compat' when it is not
% available.
%
% See also FT_HASTOOLBOX, FT_CHECKCONFIG, FT_TRACKUSAGE, LICENSE

% undocumented options
%   ft_default.siunits        = 'yes' or 'no'

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

global ft_default

persistent initialized
persistent checkpath

if isempty(initialized)
  initialized = false;
end

if isempty(checkpath)
  checkpath = false;
end

% ft_warning is located in fieldtrip/utilities, which may not be on the path yet
if ~exist('ft_warning', 'file')
  ft_warning = @warning;
end

% locate the file with the persistent FieldTrip preferences
fieldtripprefs = fullfile(prefdir, 'fieldtripprefs.mat');
if exist(fieldtripprefs, 'file')
  prefs       = load(fieldtripprefs); % the file contains multiple fields
  ft_default  = mergestruct(ft_default, prefs);
end

% Set the defaults in a global variable, ft_checkconfig will copy these over into the local configuration.
% NOTE ft_getopt might not be available on the path at this moment and can therefore not yet be used.
% NOTE all options here should be explicitly listed as allowed in ft_checkconfig

if ~isfield(ft_default, 'checkconfig'),       ft_default.checkconfig    = 'loose';    end % pedantic, loose, silent
if ~isfield(ft_default, 'checkpath'),         ft_default.checkpath      = 'pedantic'; end % pedantic, once, no
if ~isfield(ft_default, 'checksize'),         ft_default.checksize      = 1e5;        end % number in bytes, can be inf
if ~isfield(ft_default, 'showlogo'),          ft_default.showlogo       = 'yes';      end % yes or no, this is relevant for SPM and EEGLAB
if ~isfield(ft_default, 'showcallinfo'),      ft_default.showcallinfo   = 'yes';      end % yes or no, this is used in ft_pre/postamble_provenance
if ~isfield(ft_default, 'debug'),             ft_default.debug          = 'no';       end % no, save, saveonerror, display, displayonerror, this is used in ft_pre/postamble_debug
if ~isfield(ft_default, 'outputfilepresent'), ft_default.outputfilepresent = 'overwrite'; end % can be keep, overwrite, error

% These options allow to disable parts of the provenance
if ~isfield(ft_default, 'trackcallinfo'),  ft_default.trackcallinfo  = 'yes';    end % yes or no
if ~isfield(ft_default, 'trackdatainfo'),  ft_default.trackdatainfo  = 'no';     end % yes or no
if ~isfield(ft_default, 'tracktimeinfo'),  ft_default.tracktimeinfo  = 'yes';    end % yes or no
if ~isfield(ft_default, 'trackmeminfo')
  if ispc()
    % don't track memory usage info under Windows; this does not work (yet)
    ft_default.trackmeminfo   = 'no';
  else
    ft_default.trackmeminfo   = 'yes';
  end
end

if strcmp(license(), '38957')
  % within the Donders and Radboud University the default is compat, since we do not have that many licenses
  default_images = 'compat';
  default_stats  = 'compat';
  default_signal = 'compat';
else
  % in this case it will still check whether the Mathworks toolbox is available and if not, the compatible version will be added
  default_images = 'matlab';
  default_stats  = 'matlab';
  default_signal = 'matlab';
end

% These options allow to prefer the MATLAB toolbox implementations ('matlab') over the drop-in replacements ('compat').
if ~isfield(ft_default, 'toolbox'), ft_default.toolbox  = []; end
if ~isfield(ft_default.toolbox, 'images'), ft_default.toolbox.images  = default_images; end % matlab or compat
if ~isfield(ft_default.toolbox, 'stats') , ft_default.toolbox.stats   = default_stats;  end % matlab or compat
if ~isfield(ft_default.toolbox, 'signal'), ft_default.toolbox.signal  = default_signal; end % matlab or compat

% Some people mess up their path settings and then have stuff on the path that should not be there.
% The following will issue a warning
switch ft_default.checkpath
  case 'pedantic'
    % check every time
    checkIncorrectPath();
  case 'once'
    % check only once
    if ~checkpath
      checkIncorrectPath();
      checkpath = true;
    end
  case 'no'
    % do not check
end % case

% Check whether this ft_defaults function has already been executed. Note that we
% should not use ft_default itself directly, because the user might have set stuff
% in that struct already prior to ft_defaults being called for the first time.
if initialized && exist('ft_hastoolbox', 'file')
  return;
end

if isfield(ft_default, 'toolbox') && isfield(ft_default.toolbox, 'cleanup')
  prevcleanup = ft_default.toolbox.cleanup;
else
  prevcleanup = {};
end

if strcmp(ft_default.showlogo, 'yes')
  % show a welcome message
  fprintf([
    '-------------------------------------------------------------------------------------------' ...
    '\nFieldTrip is developed by members and collaborators of the Donders Institute for Brain,' ...
    '\nCognition and Behaviour at Radboud University, Nijmegen, the Netherlands.' ...
    '\n'...
    '\n                          --------------------------'...
    '\n                        /                            \\'...
    '\n                     ------------------------------------'...
    '\n                    /                                    \\'...
    '\n          -------------------------------------------------'...
    '\n         /                            /\\/\\/\\/\\/\\ '...
    '\n         ---------------------------------------------------'...
    '\n                  |        F  i  e  l  d  T  r  i  p       |'...
    '\n                  ------------------------------------------'...
    '\n                   \\                                      /'...
    '\n                     ------------------------------------'...
    '\n                          \\            /'...
    '\n                            ----------'...
    '\n' ...
    '\nPlease cite the FieldTrip reference paper when you have used FieldTrip in your study.' ...
    '\nRobert Oostenveld, Pascal Fries, Eric Maris, and Jan-Mathijs Schoffelen. FieldTrip: Open' ...
    '\nSource Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data.' ...
    '\nComputational Intelligence and Neuroscience, vol. 2011, Article ID 156869, 9 pages, 2011.' ...
    '\ndoi:10.1155/2011/156869.' ...
    '\n-------------------------------------------------------------------------------------------\n'
    ]);
end % if showlogo

% Ensure that the path containing ft_defaults is on the path.
% This allows people to do "cd path_to_fieldtrip; ft_defaults"
ftPath = fileparts(mfilename('fullpath')); % get the full path to this function, strip away 'ft_defaults'
ftPath = strrep(ftPath, '\', '\\');
if isempty(regexp(path, [ftPath pathsep '|' ftPath '$'], 'once'))
  ft_warning('FieldTrip is not yet on your MATLAB path, adding %s', strrep(ftPath, '\\', '\'));
  addpath(ftPath);
end

if ~isdeployed

  if isempty(which('ft_hastoolbox')) || isempty(which('ft_warning'))
    % the fieldtrip/utilities directory contains the ft_hastoolbox and ft_warning
    % functions, which are required for the remainder of this script
    addpath(fullfile(fileparts(which('ft_defaults')), 'utilities'));
  end

  % Some people mess up their path settings and then have different versions of certain toolboxes on the path.
  % The following will issue a warning
  checkMultipleToolbox('FieldTrip',           'ft_defaults.m');
  checkMultipleToolbox('spm',                 'spm.m');
  checkMultipleToolbox('mne',                 'fiff_copy_tree.m');
  checkMultipleToolbox('eeglab',              'eeglab2fieldtrip.m');
  checkMultipleToolbox('dipoli',              'write_tri.m');
  checkMultipleToolbox('eeprobe',             'read_eep_avr.mexa64');
  checkMultipleToolbox('yokogawa',            'GetMeg160ChannelInfoM.p');
  checkMultipleToolbox('simbio',              'sb_compile_vista.m');
  checkMultipleToolbox('fns',                 'fns_region_read.m');
  checkMultipleToolbox('bemcp',               'bem_Cii_cst.mexa64');
  checkMultipleToolbox('bci2000',             'load_bcidat.m');
  checkMultipleToolbox('openmeeg',            'openmeeg_helper.m');
  checkMultipleToolbox('freesurfer',          'vox2ras_ksolve.m');
  checkMultipleToolbox('fastica',             'fastica.m');
  checkMultipleToolbox('besa',                'readBESAmul.m');
  checkMultipleToolbox('neuroshare',          'ns_GetAnalogData.m');
  checkMultipleToolbox('ctf',                 'setCTFDataBalance.m');
  checkMultipleToolbox('afni',                'WriteBrikHEAD.m');
  checkMultipleToolbox('gifti',               '@gifti/display.m');
  checkMultipleToolbox('sqdproject',          'sqdread.m');
  checkMultipleToolbox('xml4mat',             'xml2mat.m');
  checkMultipleToolbox('cca',                 'ccabss.m');
  checkMultipleToolbox('bsmart',              'armorf.m');
  checkMultipleToolbox('iso2mesh',            'iso2meshver.m');
  checkMultipleToolbox('bct',                 'degrees_und.m');
  checkMultipleToolbox('yokogawa_meg_reader', 'getYkgwHdrEvent.p');
  checkMultipleToolbox('biosig',              'sopen.m');
  checkMultipleToolbox('icasso',              'icassoEst.m');

  try
    % external/signal contains alternative implementations of some signal processing functions
    external_signal = fullfile(fileparts(which('ft_defaults')), 'external', 'signal');
    if ~ft_platform_supports('signal') || ~strcmp(ft_default.toolbox.signal, 'matlab') || ~ft_hastoolbox('signal')
      addpath(external_signal);
    elseif contains(path, external_signal)
      rmpath(external_signal);
    end
  end

  try
    % external/stats contains alternative implementations of some statistics functions
    external_stats = fullfile(fileparts(which('ft_defaults')), 'external', 'stats');
    if ~ft_platform_supports('stats') || ~strcmp(ft_default.toolbox.stats, 'matlab') || ~ft_hastoolbox('stats')
      addpath(external_stats);
    elseif contains(path, external_stats)
      rmpath(external_stats);
    end
  end

  try
    % external/images contains alternative implementations of some image processing functions
    external_images = fullfile(fileparts(which('ft_defaults')), 'external', 'images');
    if ~ft_platform_supports('images') || ~strcmp(ft_default.toolbox.images, 'matlab') || ~ft_hastoolbox('images')
      addpath(external_images);
    elseif contains(path, external_images)
      rmpath(external_images);
    end
  end

  try
    % this directory contains various functions that were obtained from elsewere, e.g. MATLAB file exchange
    ft_hastoolbox('fileexchange', 3, 1); % not required
  end

  try
    % these directories deal with compatibility with older MATLAB versions
    if ft_platform_supports('matlabversion', -inf, '2008a'), ft_hastoolbox('compat/matlablt2008b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2008b'), ft_hastoolbox('compat/matlablt2009a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2009a'), ft_hastoolbox('compat/matlablt2009b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2009b'), ft_hastoolbox('compat/matlablt2010a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2010a'), ft_hastoolbox('compat/matlablt2010b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2010b'), ft_hastoolbox('compat/matlablt2011a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2011a'), ft_hastoolbox('compat/matlablt2011b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2011b'), ft_hastoolbox('compat/matlablt2012a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2012a'), ft_hastoolbox('compat/matlablt2012b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2012b'), ft_hastoolbox('compat/matlablt2013a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2013a'), ft_hastoolbox('compat/matlablt2013b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2013b'), ft_hastoolbox('compat/matlablt2014a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2014a'), ft_hastoolbox('compat/matlablt2014b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2014d'), ft_hastoolbox('compat/matlablt2015a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2015a'), ft_hastoolbox('compat/matlablt2015b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2015b'), ft_hastoolbox('compat/matlablt2016a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2016a'), ft_hastoolbox('compat/matlablt2016b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2016b'), ft_hastoolbox('compat/matlablt2017a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2017a'), ft_hastoolbox('compat/matlablt2017b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2017b'), ft_hastoolbox('compat/matlablt2018a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2018a'), ft_hastoolbox('compat/matlablt2018b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2018b'), ft_hastoolbox('compat/matlablt2019a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2019a'), ft_hastoolbox('compat/matlablt2019b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2019b'), ft_hastoolbox('compat/matlablt2020a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2020a'), ft_hastoolbox('compat/matlablt2020b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2020b'), ft_hastoolbox('compat/matlablt2021a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2021a'), ft_hastoolbox('compat/matlablt2021b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2021b'), ft_hastoolbox('compat/matlablt2022a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2022a'), ft_hastoolbox('compat/matlablt2022b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2022b'), ft_hastoolbox('compat/matlablt2023a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2023a'), ft_hastoolbox('compat/matlablt2023b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2023b'), ft_hastoolbox('compat/matlablt2024a', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2024a'), ft_hastoolbox('compat/matlablt2024b', 3, 1); end
    if ft_platform_supports('matlabversion', -inf, '2024b'), ft_hastoolbox('compat/matlablt2025a', 3, 1); end
    % this deals with compatibility with all OCTAVE versions
    if ft_platform_supports('octaveversion', -inf, +inf),    ft_hastoolbox('compat/octave', 3, 1); end
  end

  try
    % these contains template layouts, neighbour structures, MRIs and cortical meshes
    ft_hastoolbox('template/layout',      1, 1);
    ft_hastoolbox('template/anatomy',     1, 1);
    ft_hastoolbox('template/headmodel',   1, 1);
    ft_hastoolbox('template/electrode',   1, 1);
    ft_hastoolbox('template/neighbours',  1, 1);
    ft_hastoolbox('template/sourcemodel', 1, 1);
  end

  try
    % this is used in ft_statistics
    ft_hastoolbox('statfun', 1, 1);
  end

  try
    % this is used in ft_definetrial
    ft_hastoolbox('trialfun', 1, 1);
  end

  try
    % this contains the low-level reading functions
    ft_hastoolbox('fileio', 1, 1);
  end

  try
    % this is for filtering etc. on time-series data
    ft_hastoolbox('preproc', 1, 1);
  end

  try
    % this contains forward models for the EEG and MEG volume conductor
    ft_hastoolbox('forward', 1, 1);
  end

  try
    % this contains inverse source estimation methods
    ft_hastoolbox('inverse', 1, 1);
  end

  try
    % this contains intermediate-level plotting functions, e.g. multiplots and 3-d objects
    ft_hastoolbox('plotting', 1, 1);
  end

  try
    % this contains intermediate-level functions for spectral analysis
    ft_hastoolbox('specest', 1, 1);
  end

  try
    % this contains the functions to compute connectivity metrics
    ft_hastoolbox('connectivity', 1, 1);
  end

  try
    % this contains test scripts
    ft_hastoolbox('test', 1, 1);
  end

  try
    % this contains the functions for spike and spike-field analysis
    ft_hastoolbox('contrib/spike', 1, 1);
  end

  try
    % this contains user contributed functions
    ft_hastoolbox('contrib/misc', 1, 1);
  end

  try
    % this contains specific code and examples for realtime processing
    ft_hastoolbox('realtime/example', 3, 1);    % not required
    ft_hastoolbox('realtime/online_mri', 3, 1); % not required
    ft_hastoolbox('realtime/online_meg', 3, 1); % not required
    ft_hastoolbox('realtime/online_eeg', 3, 1); % not required
  end

end

% the toolboxes added by this function should not be removed by FT_POSTAMBLE_HASTOOLBOX
ft_default.toolbox.cleanup = prevcleanup;

% track the usage of this function, this only happens once at startup
ft_trackusage('startup');

% remember that the function has executed in a persistent variable
initialized = true;

end % function ft_default

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function checkMultipleToolbox(toolbox, keyfile)

persistent warned
if isempty(warned)
  warned = false;
end

if ~ft_platform_supports('which-all')
  return;
end

list = which(keyfile, '-all');
if length(list)>1
  ft_warning('Multiple versions of %s on your path will confuse FieldTrip', toolbox);
  if ~warned % only throw the following warnings once
    warned = true;
    for i=1:length(list)
      ft_warning('one version of %s is found here: %s', toolbox, list{i});
    end
  end
  ft_warning('You probably used addpath(genpath(''path_to_fieldtrip'')), this can lead to unexpected behavior. See http://www.fieldtriptoolbox.org/faq/should_i_add_fieldtrip_with_all_subdirectories_to_my_matlab_path');
end
end % function checkMultipleToolbox

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function checkIncorrectPath
p = fileparts(mfilename('fullpath'));
incorrect = fullfile(p, 'compat', 'incorrect');
if ~isempty(strfind(path, incorrect))
  ft_warning('Your path is set up incorrectly. You probably used addpath(genpath(''path_to_fieldtrip'')), this can lead to unexpected behavior. See http://www.fieldtriptoolbox.org/faq/should_i_add_fieldtrip_with_all_subdirectories_to_my_matlab_path');
end
end % function checkIncorrectPath

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Fit a piece-wise linear regression model.
% Here is the model
%
%  X \
%  | |
%  Q |
%  | /
%  Y
%
% where all arcs point down.
% We condition everything on X, so X is a root node. Q is a softmax, and Y is a linear Gaussian.
% Q is hidden, X and Y are observed.

X = 1;
Q = 2;
Y = 3;
dag = zeros(3,3);
dag(X,[Q Y]) = 1;
dag(Q,Y) = 1;
ns = [1 2 1]; % make X and Y scalars, and have 2 experts
dnodes = [2];
onodes = [1 3];
bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'observed', onodes);

IRLS_iter = 10;
clamped = 0;

bnet.CPD{1} = root_CPD(bnet, 1);

% start with good initial params
w = [-5 5];  % w(:,i) is the normal vector to the i'th decisions boundary
b = [0 0];  % b(i) is the offset (bias) to the i'th decisions boundary

mu = [0 0];
sigma = 1;
Sigma = repmat(sigma*eye(ns(Y)), [ns(Y) ns(Y) ns(Q)]);
W = [-1 1];
W2 = reshape(W, [ns(Y) ns(X) ns(Q)]);

bnet.CPD{2} = softmax_CPD(bnet, 2, w, b,  clamped, IRLS_iter);
bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', mu, 'cov', Sigma, 'weights', W2);

engine = jtree_inf_engine(bnet);

evidence = cell(1,3);
evidence{X} = 0.68;

engine = enter_evidence(engine, evidence);

m = marginal_nodes(engine, Y);
m.mu

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function train_id_net_vgg16(varargin)
% -------------------------------------------------------------------------
% Part 4.1: prepare the data
% -------------------------------------------------------------------------

imdb = load('./dataset/CUHK-PEDES-prepare/url_data.mat');
imdb = imdb.imdb;
imdb.images.set(imdb.images.set==2)=3;
load('./dataset/CUHK-PEDES-prepare/cuhk_word2.mat');
imdb.charcnn = wordcnn;
%imdb.charmean = mean(imdb.charcnn(:,:,:,imdb.images.set==1),4);
% -------------------------------------------------------------------------
% Part 4.2: initialize a CNN architecture
% -------------------------------------------------------------------------
net = cuhk_word2_Rankloss_vgg16();
net.conserveMemory = true;
im_mean = imdb.rgbMean;
net.meta.normalization.averageImage = im_mean;
%net.meta.normalization.charmean = imdb.charmean;
% -------------------------------------------------------------------------
% Part 4.3: train and evaluate the CNN

% -------------------------------------------------------------------------
opts.train.averageImage = net.meta.normalization.averageImage;
opts.train.batchSize = 32;
opts.train.continue = true;
opts.train.gpus = 1;
opts.train.prefetch = false ;
opts.train.nesterovUpdate = true ;
opts.train.expDir = './data/vgg16_cuhk_batch32_Rankloss_2:1:0.5_margin1_re';
opts.train.derOutputs = {'objective_f',2,'objective_img',1,'objective_txt',0.5} ;
%opts.train.gamma = 0.9;
opts.train.momentum = 0.9;
%opts.train.constraint = 100;
opts.train.learningRate = [0.1*ones(1,40),0.01*ones(1,20)] ;
opts.train.weightDecay = 0.0001;
opts.train.numEpochs = numel(opts.train.learningRate) ;
[opts, ~] = vl_argparse(opts.train, varargin) ;
% Call training function in MatConvNet
[net,info] = cnn_train_dag_batchsize(net, imdb, @getBatch,opts) ;

% --------------------------------------------------------------------
function inputs = getBatch(imdb,batch,opts)
% --------------------------------------------------------------------
batchsize = numel(batch);
half = batchsize/2;
label_img =  imdb.images.label(batch(1:half));
txt_batch = zeros(1,batchsize,'single');
for i= 1:half  % Yp ~ Xp
    txt_batch(i) = rand_same_class(imdb, label_img(i));
end
for i= half+1:batchsize   % Yn ~ Xp
    txt_batch(i) = rand_diff_class2(imdb, label_img(i-half));   %negative pair
end
for i= half+1:batchsize   % Xn ~ Xp
    if(imdb.images.label(batch(i))== imdb.images.label(batch(i-half)))
        batch(i) = rand_diff_class3(imdb, label_img(i-half));
    end
end
%-- img data
im_url = imdb.images.data(batch) ;
im = vl_imreadjpeg(im_url,'Pack','Resize',[224,224],'Flip',...
    'CropLocation','random','CropSize',[0.8,1],...
    'Interpolation', 'bicubic','NumThreads',16,... %'Brightness', double(0.1*imdb.rgbCovariance),...
    'SubtractAverage',imdb.rgbMean,...
    'CropAnisotropy',[3/4,4/3]);
oim = im{1};
label_img =  imdb.images.label(batch);
%-- txt data
label_txt =  imdb.images.label2(txt_batch);
txt = single(imdb.charcnn(:,txt_batch));
txtinput = zeros(1,56,7263,batchsize,'single');
for i=1:batchsize
    len = sum(txt(:,i)>0);
    location = randi(57-len);
    for j=1:len
        v = txt(j,i);
       txtinput(1,location,v,i)=1;
       location = location+1;
    end
end
txtinput = gpuArray(txtinput);
%--
inputs = {'x0',gpuArray(oim),'data2',txtinput,'label_img',label_img,'label_txt',label_txt};

author: layumi

repo: Image-Text-Embedding


text:
function roomhyp = sample_roomtype3(n, lines, vp, imgwidth, imgheight)

% 4 ways to create hypotheses
% roomtype1 + 3,top,left
% roomtype1 + 2,top
% roomtype6 + 3,bot,left
% roomtype6 + 2,bot

lc = [lines.lineclass]; % 0,1,2,3
lr = [lines.leftorright]; % 1 if right, -1 if left, 0 if neither
tb = [lines.above_horizon]; % 1 if above, -1 if below, 0 if neither

twotop = find((lc==2) & (tb==1));
twobot = find((lc==2) & (tb==-1));
threetopleft = find((lc==3) & (tb==1) & (lr==-1));
threebotleft = find((lc==3) & (tb==-1) & (lr==-1));

linecode{1} = twotop;
linecode{2} = twobot;
linecode{3} = threetopleft;
linecode{4} = threebotleft;

% recipe: [roomtypepart linecode oldcornerid newcornerid]
recipe(1,:) = [1 3 3 1];
recipe(2,:) = [1 1 3 1];
recipe(3,:) = [6 4 1 3];
recipe(4,:) = [6 2 1 3];

% roomhyp = [];
count = 0;
num_continue_without_progress = 0;
while count < n && num_continue_without_progress < 30
% pick recipe
    cur_recipe = recipe(randsample(size(recipe,1), 1), :);

    % sample partial room
    if cur_recipe(1) == 1
        roomhyppart = sample_roomtype1(1, lines, vp, imgwidth, imgheight);
    elseif cur_recipe(1) == 6
        roomhyppart = sample_roomtype6(1, lines, vp, imgwidth, imgheight);
    end

    if isempty(roomhyppart)
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    % sample lines
    l1 = linecode{cur_recipe(2)};
    if length(l1)<1
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end
    ls1 = l1(randsample(length(l1), 1));

    % new corner
    oldpoint = roomhyppart.corner(cur_recipe(3)).pt;
    [newpoint degen] = line_intersect(vp{1},oldpoint,...
        lines(ls1).point1, lines(ls1).point2);
    if degen==1
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end
    % check if inside img
    MARGIN = 5;
    if ~is_in_image(newpoint, imgwidth, imgheight, MARGIN)
        num_continue_without_progress = num_continue_without_progress + 1;
        continue;
    end

    % TODO: check location of new line sample
    %

    % add to hypothesis
    count = count + 1;
    num_continue_without_progress = 0;
    roomhyp(count) = roomhyppart;
    roomhyp(count).corner(cur_recipe(4)).pt = newpoint;
    roomhyp(count).type = 3;

%     global img;
%     disp_vanish(img, lines(ls1), vp);
%     plot(oldpoint(1), oldpoint(2), 'rx', 'MarkerSize',10, 'LineWidth',2);
%     plot(newpoint(1), newpoint(2), 'bx', 'MarkerSize',10, 'LineWidth',2);
%     pause;
%     close;
end

if ~exist('roomhyp','var')
    roomhyp = [];
end

author: zouchuhang

repo: LayoutNet


text:
function [ r, info ] = r8po_fa ( n, a )

%*****************************************************************************80
%
%% R8PO_FA factors a R8PO matrix.
%
%  Discussion:
%
%    The R8PO storage format is appropriate for a symmetric positive definite
%    matrix and its inverse.  (The Cholesky factor of a R8PO matrix is an
%    upper triangular matrix, so it will be in R8GE storage format.)
%
%    Only the diagonal and upper triangle of the square array are used.
%    This same storage scheme is used when the matrix is factored by
%    R8PO_FA, or inverted by R8PO_INVERSE.  For clarity, the lower triangle
%    is set to zero.
%
%    The positive definite symmetric matrix A has a Cholesky factorization
%    of the form:
%
%      A = R' * R
%
%    where R is an upper triangular matrix with positive elements on
%    its diagonal.  This routine overwrites the matrix A with its
%    factor R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the matrix in R8PO storage.
%
%    Output, real R(N,N), the Cholesky factor R in R8GE storage.
%
%    Output, integer INFO, error flag.
%    0, normal return.
%    K, error condition.  The principal minor of order K is not
%    positive definite, and the factorization was not completed.
%
  r(1:n,1:n) = a(1:n,1:n);

  for j = 1 : n

    for k = 1 : j - 1
      t = 0.0;
      for i = 1 : k-1
        t = t + r(i,k) * r(i,j);
      end
      r(k,j) = ( r(k,j) - t ) / r(k,k);
    end

    t = 0.0;
    for i = 1 : j - 1
      t = t + r(i,j)^2;
    end

    s = r(j,j) - t;

    if ( s <= 0.0 )
      info = j;
      return;
    end

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

  end

  info = 0;
%
%  Since the Cholesky factor is stored in R8GE format, be sure to
%  zero out the lower triangle.
%
  for i = 1 : n
    for j = 1 : i-1
      r(i,j) = 0.0;
    end
  end

  return
end


text:
%PlotBrowser_Test
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 8/6/2011
classdef PlotBrowser_Test < TestCase
    properties
        simulation
    end

    methods
        function this = PlotBrowser_Test(name)
            this = this@TestCase(name);
        end
    end

    methods
        function setUp(this)
            this.simulation = edu.stanford.covert.cell.sim.SimulationFixture.load([], true);
        end

        function tearDown(~)
            close('all');
        end
    end

    methods
        function test(this)
            import edu.stanford.covert.cell.sim.analysis.PlotBrowser;
            browser = PlotBrowser(this.simulation);

            %test all plots
            messages = {};
            for i = 1:numel(browser.plotNames)
                try
                    browser.selectPlots(browser.plotNames{i});
                catch exception
                    messages{end+1} = sprintf('%20s: %s', browser.plotFunctions{i}, exception.message); %#ok<AGROW>
                end
            end
            assert(isempty(messages), ...
                sprintf('The following plot methods failed:\n%s',strjoin(sprintf('\n'),messages{:})));

            %close browser
            browser.close();
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%% Start the RL GUI
clear all;
addpath(genpath('./'))

%% Create grid world

gw_options = [];
gw_options.dims = [15,15];
grid_world = Grid_world(gw_options);

%% Create reward

r_options               = [];
r_options.x_goal        = [10,10];
r_options.x_bad         = rl_reward_firepit( 'central_pit' );

r_options.type           = 'discrete';
r_options.r_goal         =   10; % reward for reaching the goal.
r_options.r_pit          =   -10; % reward for falling into a pit.
r_options.r_step         =     -1; % reward for each time step.

r_options.grid_world    = grid_world;

reward                  = Reward_gridworld(r_options);

% Discrete reward function

reward_f                = @(s,x)reward_grid_world_goal(s,x,reward.s_goal,reward.s_bad,reward.r_goal,reward.r_pit,reward.r_step);

% Continuous reward function

%beta                    = 0.1;                          %  (s,x,beta,       s_goal,          x_goal  ,        x_bad,           r_bad,        r_goal)
%reward_fc               = @(s,x)reward_grid_world_continuous(s,x,beta, reward.s_goal,  reward.x_goal  ,  reward.s_bad,  reward.r_pit,  reward.r_goal);

%% Select agents starting position.

agent_pos               = [2,2];

%% Open GUI
global speed;
speed = 0.5;
hOut    = rl_demo_gui(grid_world,reward,agent_pos,reward_f);
handles = hOut{2};

author: epfl-lasa