text:
% Filter design
% FIR Filter Design via Spectral Factorization and Convex Optimization
% EE364 lecture, Filter design and equalization
%
%  fir_chebychev_design.m              - Chebychev design of an FIR filter given a desired H(w)
%  one_over_f_filter.m                 - Design a 1/f spectrum shaping (pink-noise) filter
%  equalizer_design.m                  - Equalizer design example
%  iir_mag_design_bandpass_max_atten.m - Maximize stopband attenuation of a bandpass IIR filter
%  fir_lin_phase_lowpass_max_atten.m   - Maximize stopband attenuation of a linear phase lowpass FIR filter
%  fir_mag_design_lowpass_max_atten.m  - Maximize stopband attenuation of a lowpass FIR filter (magnitude design)
%  iir_mag_design_lowpass_max_atten.m  - Maximize stopband attenuation of a lowpass IIR filter
%  fir_lin_phase_lowpass_min_order.m   - Minimize order of a linear phase lowpass FIR filter
%  fir_mag_design_lowpass_min_order.m  - Minimize order of a lowpass FIR filter (magnitude design)
%  fir_lin_phase_lowpass_min_ripple.m  - Minimize stopband ripple of a linear phase lowpass FIR filter
%  fir_lin_phase_lowpass_min_trans.m   - Minimize transition bandwidth of a linear phase lowpass FIR filter
%  spectral_fact.m                     - Spectral factorization using Kolmogorov 1939 approach.
help Contents

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [st, stack] = stt                                             % setting [directory, number of frames, map setting]

%% main directories
st.dr.nm   = '04';                                                     % dataset name
if     filesep == '/';  mdr = sprintf('%s','~');                       % linux
elseif filesep == '\';  mdr = sprintf('%s','D:');                      % windows
end
st.dr.pts  = [mdr filesep 'dataset' filesep 'sequences' filesep ...                                        % directory of points
              st.dr.nm filesep 'velodyne' filesep];
st.dr.pose = [mdr filesep 'dataset' filesep 'poses' filesep];                                              % directory of poses
st.dr.img  = [mdr filesep 'dataset' filesep 'color' filesep st.dr.nm filesep 'image_2' filesep];           % directory of images
st.dr.clb  = [mdr filesep 'dataset' filesep 'calibration' filesep st.dr.nm filesep];                       % directory of calibration
st.dr.rec  = [filesep 'media' filesep 'ali' filesep 'TOSHIBA EXT' filesep 'result' filesep 'mar' filesep]; % directory of record
%% setting
st.st.st   = 1;                                                        % start frames
st.st.tn   = size(dir(sprintf('%s*.bin',st.dr.pts)),1);                % number of frames
%% local grid options
st.vm.xf   = +40;                                                      % x direction and front (x: -35 ~ +55)
st.vm.xb   = -20;                                                      % x direction and behind
st.vm.yl   = +12;                                                       % left and right (y: -15 ~ +15)
st.vm.yr   = -12;
st.vm.zu   = 2.5;                                                      % z direction and up (z: 0 ~ +2.5)
st.vm.zd   = -1;                                                       % z direction and down
st.vm.nxf  = +30;                                                       % new x direction and front (x: -35 ~ +55)
st.vm.nxb  = -10;                                                       % new x direction and behind
st.vm.nyl  = +10;                                                       % new left and right (y: -15 ~ +15)
st.vm.nyr  = -10;
st.vm.bs   = 3;                                                        % blind spot radius
st.vm.tr   = 1;                                                        % minimum number of points to make a valid voxel
%% cell size
st.vx.x    = 0.2;   % 0.2
st.vx.y    = 0.2;   % 0.2
st.vx.ix   = (st.vm.xf - st.vm.xb) / st.vx.x;                          % mat size
st.vx.iy   = (st.vm.yl - st.vm.yr) / st.vx.y;
%% road detection options
st.rd.vr   = 0.02;                                                     % road blocks should have variance lower than a threshold
st.rd.mx   = 0.8;                                                      % road blocks should have height lower than a threshold
%% stack option
st.stc.sz  = 90;                                                       % stack size (history)
st.stc.in  = 60;                                                       % number of integrating measurements
stack.mat  = zeros(st.vx.ix, st.vx.iy, st.stc.sz);                     % initialize stack
stack.ind  = zeros(st.vx.ix, st.vx.iy, st.stc.sz);
stack.pts  = zeros(st.vx.ix* st.vx.iy, 3, st.stc.sz);                  % associated pts
stack.ptn  = zeros(st.vx.ix* st.vx.iy, 3, st.stc.sz);                  % associated ptn
%% foreground detection
st.fr.sz   = 2;                                                        % check neighbourhood (sensor error in x-y) 0.6
%% plot options
figure('units','normalized','outerposition',[0 0 1 1])  % figure()
st.pl.vl   = '.b';                                                     % velodyne point's shape and color
st.pl.cl   = 'g';
st.pl.az   = -25;                                                      % driverview (-90, 20) [azimuth, elevation]
st.pl.el   = 45;                                                       % globalview (-25, 45) topview(0, 90)
st.pl.x    = 0.5 * max(st.vm.nxf, abs(st.vm.nxb));                     % the coordinates size
st.pl.y    = 0.5 * st.vm.nyl;
st.pl.z    = st.vm.zu - st.vm.zd;
st.cr.bs   = 8;                                                        % coordinate bias to show on image
%% Colors
st.pl.cv   = zeros(15, 3);
st.pl.cv(1, :)   = [255, 0, 0] / 255;         % red
st.pl.cv(2, :)   = [0, 255, 0] / 255;         % green
st.pl.cv(3, :)   = [0, 0, 255] / 255;         % blue
st.pl.cv(4, :)   = [176, 23, 31] / 255;       % indian red
st.pl.cv(5, :)   = [30, 144, 255] / 255;      % dodger blue
st.pl.cv(6, :)   = [142, 142, 142] / 255;     % gray 56
st.pl.cv(7, :)   = [142, 142, 56] / 255;      % olive lab
st.pl.cv(8, :)   = [255, 215, 0] / 255;       % gold
st.pl.cv(9, :)   = [128, 0, 128] / 255;       % purple
st.pl.cv(10, :)  = [205, 198, 115] / 255;     % khaki 3
st.pl.cv(11, :)  = [255, 20, 147] / 255;      % deep pink
st.pl.cv(12, :)  = [255, 165, 0] / 255;       % orange
st.pl.cv(13, :)  = [210, 180, 140] / 255;     % tan
st.pl.cv(14, :)  = [139, 35, 35] / 255;       % brown 4
st.pl.cv(15, :)  = [139, 69, 19] / 255;       % saddle brown
st.pl.cv(16, :)  = [0, 0, 0] / 255;           % black
st.pl.cv(17, :)  = [255, 255, 255] / 255;     % white

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%demo_SingleRoomDRIR - Demonstrates the upgrade from SingleRoomDRIR to SingleRoomSRIR.
% The demo script downloads the file from the SOFA database,
% upgrades the conventions and saves it as a new file in the local database.

% #Author: Michael Mihocic: (23.12.2022)
%
% SOFA Toolbox - demo script
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% Load SOFA file
db='thk';
fn='DRIR_CR1_VSA_50RS_L.sofa';
Obj=SOFAload(['db://' fullfile('database',db,fn)]);

%% Upgrade convention
disp(['Old conventions: ' Obj.GLOBAL_SOFAConventions ' v' Obj.GLOBAL_SOFAConventionsVersion])
[Obj,modified] = SOFAupgradeConventions(Obj);
disp(['New conventions: ' Obj.GLOBAL_SOFAConventions ' v' Obj.GLOBAL_SOFAConventionsVersion])

%% save the SOFA file
SOFAfn=fullfile([SOFAdbPath,'\','sofatoolbox_test','\',fn]);
disp(['Saving:  ' SOFAfn]);
Obj=SOFAsave(SOFAfn, Obj);

author: sofacoustics

repo: SOFAtoolbox


text:
function [success,msg] = ConvertLarvaeReid2JAABA(varargin)

success = false;
msg = {};

MINDT = .001;

trxfns = {'x','y','a','b','theta','x_mm','y_mm','a_mm','b_mm','theta_mm','dt','timestamps'};
trxunits = {
  parseunits('px'),...
  parseunits('px'),...
  parseunits('px'),...
  parseunits('px'),...
  parseunits('rad'),...
  parseunits('mm'),...
  parseunits('mm'),...
  parseunits('mm'),...
  parseunits('mm'),...
  parseunits('rad'),...
  parseunits('sec'),...
  parseunits('sec')};

perframeunits = {...
  'timestamp','s'
  'frame','frame'
  'objectid','unit'
  'objectnumber','unit'
  'goodnumber','unit'
  'persistence','s'
  'area_mm','mm^2'
  'speed','mm/s'
  'angularspeed','deg/s'
  'skeletonlength','mm'
  'width','mm'
  'midlinelength','mm'
  'kink','deg'
  'bias','unit'
  'pathlength_mwt','mm'
  'curvature','deg'
  'changeindirection','unit'
  'x_mm','mm'
  'y_mm','mm'
  'velx','mm/s'
  'vely','mm/s'
  'orientation','deg'
  'crab','mm/s'
  'flux','unit'
  'headangle','deg'
  'headvecx','unit'
  'headvecy','unit'
  'tailvecx','unit'
  'tailvecy','unit'
  'tailx','px'
  'taily','px'
  'hcdist','mm'
  'headx','px'
  'heady','px'
  'heading','deg'
  'bodyangle','deg'
  'headspeed','mm/s'
  'tailspeed','mm/s'
  'centerspeed','mm/s'
  'bodyanglespeed','deg/s'
  'concheadspeed','unit/s'
  'dfromsource','mm'
  'source','unit'
  'Chead','unit'
  'Ctail','unit'
  'sourceangle','deg'
  'curvingrate','deg/s'
  'bearingangforcurvingrate','deg'
  'cdotc','unit/s'};

[inmoviefile,...
  blobsfile,datfiles,kinmatfile,eventmatfile,...
  expdir,moviefilestr,trxfilestr,perframedirstr,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  pxpermm,fps,overridefps,overridearena,...
  dosoftlink] = myparse(varargin,...
  'inmoviefile','',...
  'blobsfile','','datfiles',{},'kinmatfile','','eventmatfile','',...
  'expdir','','moviefilestr','movie.ufmf','trxfilestr','trx.mat','perframedirstr','perframe',...
  'arenatype','None','arenacenterx',0,'arenacentery',0,...
  'arenaradius',123,'arenawidth',123,'arenaheight',123,...
  'pxpermm',1,'fps',30,...
  'overridefps',false,'overridearena',false,...
  'dosoftlink',false);

% check that required inputs are given
if isempty(blobsfile),
  msg = 'Input blobs file is empty';
  return;
end
if ~iscell(blobsfile),
  blobsfile = {blobsfile};
end
for i = 1:numel(blobsfile),
  if ~exist(blobsfile{i},'file'),
    msg = sprintf('Input blobs file %s does not exist',blobsfile{i});
    return;
  end
end

for i = 1:numel(datfiles),
  if ~exist(datfiles{i},'file'),
    msg = sprintf('Input dat file %s does not exist',datfiles{i});
    return;
  end
end

if ~exist(kinmatfile,'file'),
  msg = sprintf('Input mat file %s does not exist',kinmatfile);
  return;
end

if ~isempty(eventmatfile) && ~exist(eventmatfile,'file'),
  msg = sprintf('Input mat file %s does not exist',eventmatfile);
  return;
end

% output file locations
moviefile = fullfile(expdir,moviefilestr);
trxfile = fullfile(expdir,trxfilestr);
perframedir = fullfile(expdir,perframedirstr);

%% read in the trx from the blobs file

trx = ReadMWTBlobFile(blobsfile);
nflies = numel(trx);

% make sure all the timestamps agree with each other
timestamps = unique([trx.timestamps]);
dts = diff(timestamps);
% sanity check to look for rounding errors
if any(dts < MINDT),
  msg = sprintf('Rounding error check failed: there are timestamps < %f apart',MINDT);
  return;
end

% add dt, fps
if ~overridefps,
  fps = 1/nanmedian(dts);
end
for i = 1:nflies,
  trx(i).dt = diff(trx(i).timestamps);
  trx(i).fps = fps;
end

trxids = [trx.id];

%% read in kin mat file

try
  load(kinmatfile,'kinData');
catch %#ok<CTCH>
  msg = sprintf('Error loading kinData from kin mat file %s',kinmatfile);
  return;
end

dataall = [];
for i = 1:numel(kinData), %#ok<USENS>
  dataall = structappend(dataall,kinData{i});
end
fns = setdiff(fieldnames(dataall),{'frame'});
for i = 1:numel(dataall),
  id = unique(dataall(i).objID);
  if numel(id) ~= 1,
    msg = sprintf('More than one objID found in kinData{%d}',i);
    return;
  end
  dataall(i).id = id; %#ok<AGROW>
  dataall(i).firstframe = min(dataall(i).frame)+1; %#ok<AGROW>
  dataall(i).endframe = max(dataall(i).frame)+1; %#ok<AGROW>
  dataall(i).nframes = dataall(i).endframe-dataall(i).firstframe+1; %#ok<AGROW>
  [ism,fidx] = ismember(dataall(i).firstframe-1:dataall(i).endframe-1,dataall(i).frame);
  if any(~ism),
    fprintf('%d frames with no data for target %d\n',nnz(~ism),id);
  end

  n = dataall(i).nframes;
  for j = 1:numel(fns),
    fn = fns{j};
    tmp = nan(1,n);
    tmp(ism) = dataall(i).(fn)(fidx);
    dataall(i).(fn) = tmp; %#ok<AGROW>
  end

end

%% read in the eventData

try
  load(eventmatfile,'eventData');
catch %#ok<CTCH>
  msg = sprintf('Error loading eventData from kin mat file %s',eventmatfile);
  return;
end

if numel(eventData) > numel(kinData), %#ok<USENS>
  msg = sprintf('Less trajectories for kinData file %s than eventData file %s',kinmatfile,eventmatfile);
  return;
end

% add eventData to dataall
fns = fieldnames(eventData{1});
for i = 1:numel(eventData),
  datacurr = eventData{i};
  n = dataall(i).nframes;
  for j = 1:numel(fns),
    fn = fns{j};
    if max(datacurr.(fn)) > n,
      msg = sprintf('Maximum index for %s bigger than number of frames in trajectory %d',fn,i);
      return;
    end
    tmp = zeros(1,n);
    tmp(datacurr.(fn)) = 1;
    dataall(i).(['event_',fn]) = tmp; %#ok<AGROW>
  end
end

%% match up ids

% find starts of videos

% which sequence does the kindata correspond to
dataids = [dataall.id];
isstart = [true,dataids(1:end-1)>dataids(2:end)];
idxstart = find(isstart);
idxend = [idxstart(2:end)-1,numel(dataids)];

maxmatches = -1;
for i = 1:numel(idxstart),
  idxcurr = idxstart(i):idxend(i);
  dataids_curr = [dataall(idxcurr).id];
  [ism,order] = ismember(dataids_curr,trxids);
  ismatch = [dataall(idxcurr(ism)).nframes] == [trx(order(ism)).nframes];
  fracmatches = nnz(ismatch) / nnz(ism);
  if fracmatches > maxmatches,
    bestmatch = i;
    maxmatches = fracmatches;
  end
end

msg{end+1} = sprintf('Choosing to match kindata sequence %d with %s, fraction of ids for which nframes matches = %f',bestmatch,blobsfile{1},maxmatches);
dataall = dataall(idxstart(bestmatch):idxend(bestmatch));
dataids = [dataall.id];

if numel(unique(dataids)) < numel(dataids),
  msg = 'Some repeat ids found in dataids';
  return;
end
badtrxids = ~ismember(trxids,dataids);
baddataids = ~ismember(dataids,trxids);
if any(badtrxids) || any(baddataids),
  msg{end+1} = sprintf('%d ids in blobs not in kindata, %d ids in kindata not in blobs',nnz(badtrxids),nnz(baddataids));
end
dataall(baddataids) = [];
dataids(baddataids) = [];
trx(badtrxids) = [];
trxids(badtrxids) = [];
[~,idx] = ismember(dataids,trxids);
dataall = dataall(idx);

badids = false(1,numel(trx));
for i = 1:numel(trx),
  if trx(i).firstframe ~= dataall(i).firstframe || ...
      trx(i).endframe ~= dataall(i).endframe,
    badids(i) = true;
  end
end
if any(badids),
  msg{end+1} = sprintf('%d id nframes do not match between kindata and blobs files for some ids, removing these',nnz(badids));
  dataall(badids) = [];
  trx(badids) = [];
end

if isempty(trx),
  msg{end+1} = 'After removing mismatches between kindata and blobs, no trajectories left.';
  return;
end

%% add trajectory features that are in kindata files

for i = 1:numel(trx),
  if isfield(dataall,'headx'),
    trx(i).headx = dataall(i).headx;
  end
  if isfield(dataall,'heady'),
    trx(i).heady = dataall(i).heady;
  end
  if isfield(dataall,'tailx'),
    trx(i).tailx = dataall(i).tailx;
  end
  if isfield(dataall,'taily'),
    trx(i).tailx = dataall(i).taily;
  end
  if isfield(dataall,'cx'),
    trx(i).x_mm = dataall(i).cx;
  end
  if isfield(dataall,'cy'),
    trx(i).y_mm = dataall(i).cy;
  end
end

%% convert to mm

if strcmpi(arenatype,'None'),
  arenacenterx = 0;
  arenacentery = 0;
end

if isfield(trx,'x_mm'),
  if ~overridearena,
    idx = [trx.x_mm] > 0;
    tmp = [trx.x]./[trx.x_mm];
    s = nanstd(tmp);
    if s < 1,
      tmp = nanmean(tmp(idx));
      if ~isnan(tmp),
        pxpermm = tmp;
      end
    end
  end
  dx = [trx.x_mm] - ([trx.x]-arenacenterx)/pxpermm;
  if nanstd(dx) > 1,
    msg = sprintf('pxpermm does not match the relationship between x and x_mm');
    return;
  end
  arenacenterx_mm = nanmean(dx);
else
  arenacenterx_mm = 0;
end

if isfield(trx,'y_mm'),
  dy = [trx.y_mm] - ([trx.y]-arenacentery)/pxpermm;
  if nanstd(dy) > 1,
    msg = sprintf('pxpermm does not match the relationship between y and y_mm');
    return;
  end
  arenacentery_mm = nanmean(dy);
else
  arenacentery_mm = 0;
end

for i = 1:numel(trx),
  trx(i).pxpermm = pxpermm;
  if ~isfield(trx,'x_mm'),
    trx(i).x_mm = (trx(i).x - arenacenterx) / pxpermm;
  end
  if ~isfield(trx,'y_mm'),
    trx(i).y_mm = (trx(i).y - arenacentery) / pxpermm;
  end
  trx(i).a_mm = trx(i).a / pxpermm;
  trx(i).b_mm = trx(i).b / pxpermm;
  trx(i).theta_mm = trx(i).theta;
  if isfield(trx,'xspine_mm'),
    trx(i).xspine = trx(i).xspine_mm * pxpermm;
  end
  if isfield(trx,'yspine_mm'),
    trx(i).yspine = trx(i).yspine_mm * pxpermm;
  end
end

%% over-ride timestamps

if overridefps,
  timestamps = timestamps(1)+(0:numel(timestamps)-1)/fps;
  for i = 1:numel(trx),
    trx(i).timestamps = timestamps(trx(i).firstframe:trx(i).endframe);
    trx(i).dt(:) = 1/fps;
  end
end

%% set landmark parameters

arenaradius_mm = arenaradius / pxpermm;
arenawidth_mm = arenawidth / pxpermm;
arenaheight_mm = arenaheight / pxpermm;
trx = SetLandmarkParameters(trx,arenatype,arenacenterx_mm,arenacentery_mm,...
  arenaradius_mm,arenawidth_mm,arenaheight_mm);

%% create directory

% create the experiment directory
if ~exist(expdir,'dir'),
  [success1,msg] = mkdir(expdir);
  if ~success1,
    return;
  end
end

%% create per-frame directory

if ~exist(perframedir,'dir'),
  [success1,msg] = mkdir(perframedir);
  if ~success1,
    return;
  end
end

%% save the per-frame data

perframefns = setdiff(fieldnames(dataall),{'frame','objID','objN','goodN','persists','id','firstframe','endframe','nframes','times'});

for i = 1:numel(perframefns),
  fn = perframefns{i};
  data = {dataall.(fn)}; %#ok<NASGU>
  j = find(strcmp(perframeunits(:,1),fn),1);
  if ~isempty(j),
    units = parseunits(perframeunits{j,2}); %#ok<NASGU>
  else
    units = parseunits('unit'); %#ok<NASGU>
  end
  outfile = fullfile(perframedir,[fn,'.mat']);
  try
    save(outfile,'data','units');
  catch ME,
    msg = getReport(ME);
    return;
  end
end

for i = 1:numel(trxfns),
  fn = trxfns{i};
  if isfield(trx,fn),
    data = {trx.(fn)}; %#ok<NASGU>
    units = trxunits{i}; %#ok<NASGU>
    %units = parseunits('unit'); %#ok<NASGU>
    outfile = fullfile(perframedir,[fn,'.mat']);
    try
      save(outfile,'data','units');
    catch ME,
      msg = getReport(ME);
      return;
    end
  end
end

%% create the trx file

try
  save(trxfile,'trx','timestamps');
catch ME,
  msg = getReport(ME);
  return;
end

if ~exist(trxfile,'file'),
  msg = sprintf('Failed to save trx to file %s',trxfile);
  return;
end

%% copy/soft-link movie

if ~isempty(inmoviefile),

  if strcmp(fullfile(inmoviefile),fullfile(moviefile)),
    fprintf('Input and out movie files are the same, not copying/linking.\n');
  else

    if dosoftlink,
      if isunix,
        cmd = sprintf('ln -s %s %s',inmoviefile,moviefile);
        unix(cmd);
        % test to make sure it worked
        [status,result] = unix(sprintf('readlink %s',moviefile));
        result = strtrim(result);
        if status ~= 0 || ~strcmp(result,inmoviefile),
          warndlg(sprintf('Failed to make soft link, copying %s to %s instead',inmoviefile,moviefile));
          dosoftlink = false;
        end
      elseif ispc,
        cmd = sprintf('mkshortcut.vbs /target:"%s" /shortcut:"%s"',inmoviefile,moviefile);
        fprintf('Making a Windows shortcut file at "%s" with target "%s"\n',inmoviefile,moviefile);
        system(cmd);
        % test to make sure that worked
        [equalmoviefile,didfind] = GetPCShortcutFileActualPath(moviefile);
        if ~didfind || ~strcmp(equalmoviefile,inmoviefile),
          warndlg(sprintf('Failed to make shortcut, copying %s to %s instead',inmoviefile,moviefile));
          dosoftlink = false;
        end
      else
        warndlg(sprintf('Unknown OS, not soft-linking movie file %s',inmoviefile));
        dosoftlink = false;
      end
    end

    if ~dosoftlink,
      [success1,msg1] = copyfile(inmoviefile,moviefile);
      if ~success1,
        msg = msg1;
        success = false;
        return;
      end
    end

  end
end

success = true;

author: kristinbranson

repo: JAABA


text:
function [outC,ptListC] = stackDLMaskFiles(outPath,outFmt,passedScanDim)
% stackDLMaskFiles.m Reads output mask files and returns 3D stack.
%--------------------------------------------------------------------------
%INPUTS:
% outPath       : Path to generated files
%                 Note: Assumes output filenames are of the form:
%                 prefix_slice# if  passedScanDim = '2D' and
%                 prefix_3D if passedScanDim = '3D'.
% outFmt        : Output format ('H5' or 'NRRD').
% passedScanDim : '2D' or '3D'.
%------------------------------------------------------------------------
% AI 6/29/21

switch outFmt

    case 'H5'
        [outC,ptListC] = stackHDF5Files(outPath,passedScanDim);

    case 'NRRD'

        dirS = dir(fullfile(outPath,'outputNRRD','*.nrrd'));
        fileNameC = {dirS.name};
        ptListC = unique(strtok(fileNameC,'_'));
        outC = cell(length(ptListC),1);

        for p = 1:length(ptListC)

            %Assumes 3D mask file
            fileName = fullfile(outPath,'outputNRRD',fileNameC{1});
            mask3M = nrrdread_opensrc(fileName);
            outC{p} = mask3M;
        end

    otherwise

        error('invalid output format %s',outFmt);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function Fitness = Fit(PopObj,PopCon)
% Calculate the fitness of each solution

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    N  = size(PopObj,1);
    CV = sum(max(0,PopCon),2);

    %% Detect the dominance relation between each two solutions
    Dominate = false(N);
    for i = 1 : N-1
        for j = i+1 : N
            if CV(i) < CV(j)
                Dominate(i,j) = true;
            elseif CV(i) > CV(j)
                Dominate(j,i) = true;
            else
                k = any(PopObj(i,:)<PopObj(j,:))-any(PopObj(i,:)>PopObj(j,:));
                if k == 1
                    Dominate(i,j) = true;
                elseif k == -1
                    Dominate(j,i) = true;
                end
            end
        end
    end

    %% Calculate S(i)
    S = sum(Dominate,2);

    %% Calculate R(i)：
    R = zeros(1,N);
    for i = 1 : N
        R(i) = sum(S(Dominate(:,i)));
    end

    %% Calculate D(i)
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Distance = sort(Distance,2);
    D = 1./(Distance(:,floor(sqrt(N)))+2);

    %% Calculate the fitnesses
    Fitness = R + D';
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function tq = gettrimeshquan( p, t)
% tq.ds  - edge length connecting node
%           2-3, 3-1, and 1-2 respectively
% tq.J   - 0.5*element area
% tq.vang - angle assocaited with node 1, 2, 3 respectively
%

edge = [2  3  1
        3  1  2] ;

ne = size( t, 1) ;

p1 = p(t(:, edge(1,:))',:)  ;
p2 = p(t(:, edge(2,:))',:)  ;

%
ds = reshape(sqrt(sum((p1 - p2).^2,2)), 3, ne )' ;

idxcm = [ 1 2 3
          2 3 1
          3 1 2 ] ;

% Internal angle
tq.vang = zeros(ne,3) ;
for i = 1: 3
   tq.vang(:,i) = acos((-ds(:,idxcm(i,1)).^2 +  ds(:,idxcm(i,2)).^2 + ...
           ds(:,idxcm(i,3)).^2)./(2*ds(:,idxcm(i,2)).*ds(:,idxcm(i,3)))) ;
end
%

% Length
tq.ds = ds ;

% Jacobian
xr = 0.5*(p(t(:,2),:) - p(t(:,1),:)) ;
xs = 0.5*(p(t(:,3),:) - p(t(:,1),:)) ;

tq.J = (xr(:,1).*xs(:,2) - xs(:,1).*xr(:,2)) ;

% idxtb = find( vang' < 25*pi/180 ) ;
% ietb = unique(ceil(idxtb/3)) ;

%
%  Bank, Randolph E.,
%      PLTMG: A Software Package for Solving Elliptic Partial Differential Equations,
%      User's Guide 6.0, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1990.
tq.qm = 4*sqrt(3)*(2*tq.J)./(sum(tq.ds.^2,2)) ;

author: CHLNDDEV


text:
%msre_newton_solver is a function.
%    [T, retcode, c, Q] = msre_newton_solver(obj, inner_probs, T0, x0, control_shocks, options, kron_method)
%

author: jmaih


text:
function varargout = funlogger(x,varargin)
%FUNLOGGER Record log of calls to a specific function.
%
%   FUNLOGGER([],FUN) initialize logger for function FUN. FUN must be either
%   a string with the function name or a function handle for function F(X).
%   F takes a numeric array as input and returns a single scalar as output.
%
%   FUNLOGGER([],FUN,N) stores up to N function calls (default N=1e4).
%
%   FUNLOGGER([],FUN,N,OUTSIZE) specifies that FUN returns a numeric array
%   with OUTSIZE elements (default OUTSIZE=1).
%
%   FVAL = FUNLOGGER(X) returns the value of the logged function at X,
%   where X is a numeric array.
%
%   [FVAL,FUNLOG] = FUNLOGGER(X) also returns the function log structure.
%
%   [...] = FUNLOGGER(X,ARG1,ARG2,...) passes a series of additional
%   arguments to the call of logged function.
%
%   FUNLOG = FUNLOGGER() returns the function log structure.
%
%
%   Author: Luigi Acerbi
%   Version: Sep/11/2015
%
persistent funlog;     % Record log of function calls

%% Case 1: No arguments ==> return function log
if nargin < 1
    varargout = {funlog};

    %if isstruct(x)  % Swapped variable order
    %    temp = probstruct;
    %    probstruct = x;
    %    x = temp;
    %end

%% Case 2: Empty X, no output, pass function name or handle ==> Re-initialize log struct
elseif isempty(x) && nargout == 0 && ...
        (ischar(varargin{1}) || isa(varargin{1},'function_handle'))
    funlog = [];
    fun = varargin{1};
    if nargin > 2; N = varargin{2}; else N = []; end
    if nargin > 3; outsize = varargin{3}; else outsize = []; end

    % Store function name and function handle
    if ischar(fun)
        funlog.FuncName = fun;
        funlog.FuncHandle = str2func(fun);
    else
        funlog.FuncName = func2str(fun);
        funlog.FuncHandle = fun;
    end

    % Number of stored function values
    if isempty(N); N = 1e4; end

    % Size of function output arguments
    if isempty(outsize); outsize = 1; end
    funlog.OutSize = outsize;

    funlog.Clock = tic;     % Time of initialization
    funlog.FuncCount = 0;
    funlog.X = [];
    funlog.Y = [];

    funlog.N = N;           % Size of X matrix, temporary
    funlog.last = 0;        % Last entry

%% Case 3: Pass X, get output ==> Standard function call
else
    if isempty(x) && nargin == 2 && isnumeric(varargin{1})
        % MCS call
        x = varargin{1}';
        varargin = [];
        mcs_flag = true;
    else
        mcs_flag = false;
    end

    % You need to initialize FUNLOGGER first
    if isempty(funlog)
        error('FUNLOGGER has not been initialized.');
    end

    % Call function
    func = funlog.FuncHandle;

    % if isfield(probstruct,'octstruct'); x = optimct(x,probstruct.octstruct,1); end

    % Check if need to pass probstruct
    try
        if nargin > 1 && ~mcs_flag
            fval = func(x,varargin{:});
        else
            fval = func(x);
        end
    catch funErr
        warning(['Error in executing the logged function ''' funlog.FuncName ''' with input:']);
        x
        rethrow(funErr);
        % fval = NaN(1,funlog.OutSize);
    end

    funlog.FuncCount = funlog.FuncCount + 1;

    % Record log
    if isempty(funlog.X)
        funlog.X = NaN(funlog.N,length(x));
        funlog.Y = NaN(funlog.N,funlog.OutSize);
    end

    funlog.last = max(1,mod(funlog.last+1, funlog.N+1));
    funlog.X(funlog.last,:) = x;
    try
        funlog.Y(funlog.last,:) = fval;
    catch
        pause
    end

    varargout{1} = fval;
    if nargout > 1; varargout{2} = funlog; end
end

end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [G, pixels, patches] = gsp_patch_graph(img, param)
%GSP_PATCH_GRAPH Create a graph by NN patches of an image
%   Usage :  G = compute_patch_graph( img );
%            [G, f] = gsp_patch_graph( img, param );
%            [G, f, patches] = gsp_patch_graph( img, param );
%
%   Input parameters:
%       img         : Input image (unknown pixels marked as NaN)
%       param       : Structure of optional parameters
%
%   Output parameters:
%       G           : Resulting graph
%       f           : Image signal
%       patches     : Patches
%
%   'compute_patch_graph( path , param )' creates a graph between pixels in
%   an image by connecting them using the euclidean distance between
%   patches around the pixels
%
%   Additional parameters
%   ---------------------
%
%   * *param.patch_size*      : int     the patch size in pixels (odd)
%   * *param.scale*           : float   to rescale the input image
%   * *param.rho*             : float   spatial constraint
%   * *use_incomplete_patch*  : boolean use incomplete patch for the graph construction
%   * *param.nnparam*         : struct  parameters for graph construction
%

% Author: Johan Paratte, Michael Defferrard, Nathanael Perraudin
% Date: November 2014
%

    if nargin < 2
    % Define parameters
        param = {};
    end

    if ~isfield(param, 'patch_size'), param.patch_size = 5; end
    if ~isfield(param, 'rho'), param.rho = 0.001; end
    if ~isfield(param, 'use_incomplete_patch'), param.use_incomplete_patch = 0; end
    if ~isfield(param, 'nnparam')
       param.nnparam = {};
       param.nnparam.center = 0;
       param.nnparam.resize = 0;
       param.nnparam.k = 10;
       param.nnparam.use_l1 = 0;
       param.nnparam.use_flann = 1;
    end

    if (mod(param.patch_size, 2) == 0)
        disp('Patch size must be odd, converting to closest odd number');
        param.patch_size = param.patch_size + 1;
    end

%     if length(size(img)) > 2
%        img = rgb2gray(img);
%     end

    %Extract patches
    [oheight, owidth,Nc] = size(img);

    psize = param.patch_size;

    %Parameters
    if isfield(param, 'scale')
        img = imresize(img, param.scale);
    end

    margin = floor(psize / 2);

    dim = psize*psize;
    %height = (oheight - 2*margin);
    %width = (owidth - 2*margin);

    %Expand the image to compensate the margin
    h = oheight;
    w = owidth;
    new_img = zeros(h + 2*margin, w + 2*margin,Nc);
    %Copy image
    new_img(1+margin:end-margin, 1+margin:end-margin,:) = img;
    %Four lines
    new_img(1:margin, 1+margin:end-margin,:) = img(margin:-1:1, :,:);
    new_img(1+margin:end-margin, 1:margin,:) = img(:, margin:-1:1,:);
    new_img(end-margin:end, 1+margin:end-margin,:) = img(end:-1:end-margin, :,:);
    new_img(1+margin:end-margin, end-margin:end,:) = img(:,end:-1:end-margin,:);
    %Four corners
    new_img(1:margin, 1:margin,:) = img(margin:-1:1, margin:-1:1,:);
    new_img(1:margin, end-margin:end,:) = img(margin:-1:1, end:-1:end-margin,:);
    new_img(end-margin:end, 1:margin,:) = img(end:-1:end-margin, margin:-1:1,:);
    new_img(end-margin:end, end:-1:end-margin,:) = img(end:-1:end-margin, end:-1:end-margin,:);

    % Signals on the graph.
    patches = zeros(h*w, dim*Nc+2);
    pixels  = double(reshape(img, h*w, Nc)); % zeros(h*w, 1);
    coords  = zeros(h*w, 2);  % gsp_plot_signal takes [x,y]
    coords(:,2) = repmat((h:-1:1).', w, 1);
    coords(:,1) = reshape(repmat((1:+1:w), h, 1), [], 1);

    % Pre-allocation.
    nUnknown = sum(sum(img(:,:,1)<0)) * 2;
    unknowns = zeros(nUnknown, 1);
    iUnknown = 0;

    count = 1;

    % For speed improvement, one of the for loop can probably be removed
    % here
    for w = 1:owidth
       for h = 1:oheight
            patches(count, 1:Nc*dim) = reshape(new_img(h:h+psize-1, w:w+psize-1,:), 1, Nc*dim);
            patches(count, Nc*dim+1:end) = [param.rho*w; param.rho*h];
            count = count + 1;
        end
    end

    fprintf('Compute graph with %d vertices\n', size(patches,1));

    if param.use_incomplete_patch

        G = gsp_rmse_mv_graph(patches, param.nnparam);
    else
        % Ignore patches which contain unknown pixel values.
        % We do not want them to be connected with known patches.
        clearedPatches = patches;

        for count = 1:size(clearedPatches,1)
            if sum(isnan(clearedPatches(count,:)))
                clearedPatches(count,:) = -1e3;
                % List of unknown patches.
                iUnknown = iUnknown + 1;
                unknowns(iUnknown) = count;
            end
        end

        G = gsp_nn_graph(clearedPatches, param.nnparam);

        % Disconnect unknown patches.
        unknowns = unknowns(1:iUnknown);
        G.W(unknowns,:) = 0;
        G.W(:,unknowns) = 0;

        % Update G.A, G.d and G.Ne.
        G = gsp_graph_default_parameters(G);
    end

    % These are coordinates in the high dimensional patch space. We prefer
    % to visualize the graph in the 2D image space.
%     patches = G.coords;
    G.coords = coords;
    G.plotting.limits = [1,w,1,h];

end

repo: gspbox


text:
% eeg_options() - eeglab option script
%
% Note: DO NOT EDIT, instead use pop_editoptions() or the menu
%       /File/Maximize memory in EEGLAB gui

% STUDY options (set these checkboxes if you intend to work with studies)
option_storedisk = 0 ; % If set, keep at most one dataset in memory. This allows processing hundreds of datasets within studies.
option_savetwofiles = 1 ; % If set, save not one but two files for each dataset (header and data). This allows faster data loading in studies.
% Memory options
option_single = 1 ; % If set, use single precision under Matlab 7.x. This saves RAM but can lead to rare numerical imprecisions.
option_memmapdata = 0 ; % If set, use memory mapped array under Matlab 7.x. This may slow down some computation.
% ICA options
option_computeica = 0 ; % If set, precompute ICA activations. This requires more RAM but allows faster plotting of component activations.
option_scaleicarms = 1 ; % If set, scale ICA component activities to RMS (Root Mean Square) in microvolt (recommended).
% Folder options
option_rememberfolder = 1 ; % If set, when browsing to open a new dataset assume the folder/directory of previous dataset.

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ n, a, seed ] = tree_rb_yule ( n, a, seed )

%*****************************************************************************80
%
%% TREE_RB_YULE adds two nodes to a rooted binary tree using the Yule model.
%
%  Discussion:
%
%    The Yule model is a simulation of how an evolutionary family tree
%    develops.  We start with a root node.  The internal nodes of the tree
%    are inactive and never change.  Each pendant or leaf node of the
%    tree represents a biological family that can spontaneously "fission",
%    developing two new distinct sub families.  In graphical terms, the node
%    becomes internal, with two new leaf nodes depending from it.
%
%    The tree is stored in inorder traversal form.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer N, the number of nodes in the input
%    tree.  On output, this number has been increased, usually by 2.
%
%    Input/output, integer A(*), the preorder traversal form
%    for the rooted binary tree.  The number of entries in A is N.
%
%    Input/output, integer SEED, a seed for the random number
%    generator.
%
  if ( n <= 0 )
    n = 1;
    a(1) = 0;
    return
  end
%
%  Count the expected number of leaves, which are the 0 values.
%
  nleaf = floor ( ( n + 1 ) / 2 );
%
%  Choose a random number between 1 and NLEAF.
%
  ileaf = i4_uniform_ab ( 1, nleaf, seed );
%
%  Locate leaf number ILEAF.
%
  j = 0;
  jleaf = 0;
  for i = 1 : n
    if ( a(i) == 0 )
      jleaf = jleaf + 1;
    end
    if ( jleaf == ileaf )
      j = i;
      break
    end
  end
%
%  Replace '0' by '100'
%
  a(n+2:-1:j+2) = a(n:-1:j);
  a(j) = 1;
  a(j+1) = 0;

  n = n + 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%DEMO_REGRESSION_ADDITIVE2  Regression demonstration with additive Gaussian
%                           process using linear, squared exponential and
%                           neural network covariance fucntions
%
%  Description
%    Gaussian process solutions in 2D regression problem using
%    constant, linear, squared exponential (sexp) and neural
%    network covariance functions, and with various additive
%    combinations of these four covariance functions. The noisy
%    observations y are assumed to satisfy
%
%         y = f + e,    where e ~ N(0, s^2)
%
%    where f is an unknown underlying function. A zero mean Gaussian
%    process prior is assumed for f
%
%         f ~ N(0, K),
%
%    where K is the covariance matrix whose elements are given by one of
%    the following six covariance function:
%
%    - constant + linear
%    - costant + sexp for 1. input + linear for 2. input
%    - sexp for 1. input + sexp for 2. input
%    - sexp
%    - neural network for 1. input + neural network for 2. input
%    - neural network
%
%    A prior is assumed for parameters of the covariance functions,
%    and the inference is done with a MAP estimate for parameter
%    values.
%
%    For more detailed discussion of  covariance functions, see e.g.
%
%    Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian
%    Processes for Machine Learning. The MIT Press.
%
%
%  See also
%    DEMO_REGRESSION1
%
% Copyright (c) 2010 Jaakko Riihimäki, Aki Vehtari

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

% REGRESSION TOY DATA
% Create an example regression data
x=rand(225,2)*2-1;
y=3*norm_cdf(2*x(:,2))+2*norm_cdf(4*x(:,1));
% add some noise
y=y+randn(size(y))*0.25;
y=y-mean(y);
[n, nin] = size(x);
% create equally spaced points to visualise the predictions:
[xt1,xt2]=meshgrid(-2:0.1:2,-2:0.1:2);
xt=[xt1(:) xt2(:)];
nxt=size(xt1,1);

% Assume a Student-t distribution for the GP parameters
pt = prior_t('nu', 4, 's2', 10);	% a prior structure

% Create a Gaussian noise model
lik = lik_gaussian('sigma2', 0.2^2);

% Set a small amount of jitter
jitter=1e-4;

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);

% CONSTANT + LINEAR COVARIANCE FUNCTION
disp('Constant + linear covariance function')
% constant covariance function
gpcf_c = gpcf_constant('constSigma2', 1, 'constSigma2_prior', pt);
% linear covariance function
gpcf_l = gpcf_linear('coeffSigma2_prior', pt);
gp = gp_set('lik', lik, 'cf', {gpcf_c gpcf_l}, 'jitterSigma2', jitter);

% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Compute predictions in a grid using the MAP estimate
Eft_map = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure
set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,nxt,nxt));
hold on
plot3(x(:,1), x(:,2), y, '*')
xlabel('x_1'), ylabel('x_2')
title('The predicted underlying function (constant + linear)');

% CONSTANT + SQUARED EXPONENTIAL COVARIANCE FUNCTION (W.R.T. THE
% FIRST INPUT DIMENSION) + LINEAR (W.R.T. THE SECOND INPUT
% DIMENSION)
fprintf(['Constant + squared exponential covariance function\n' ...
         '(w.r.t. the first input dimension) + linear (w.r.t.\n' ...
         'the second input dimension)\n'])

%Covariance function for the first input variable
gpcf_s1 = gpcf_sexp('selectedVariables', 1, 'lengthScale',0.5, ...
                    'lengthScale_prior', pt, 'magnSigma2', 0.15, ...
                    'magnSigma2_prior', pt);
% gpcf_s1 can be construted also as
%metric1 = metric_euclidean('components', {[1]}, 'lengthScale',[0.5], ...
%                           'lengthScale_prior', pt);
%gpcf_s1 = gpcf_sexp('magnSigma2', 0.15, 'magnSigma2_prior', pt, ...
%                    'metric', metric1);
%Covariance function for the second input variable
gpcf_l2 = gpcf_linear('selectedVariables', 2, 'coeffSigma2_prior', pt);
gp = gp_set('lik', lik, 'cf', {gpcf_c gpcf_s1 gpcf_l2}, 'jitterSigma2', jitter);

% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Compute predictions in a grid using the MAP estimate
Eft_map = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure
set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,nxt,nxt));
hold on
plot3(x(:,1), x(:,2), y, '*')
xlabel('x_1'), ylabel('x_2')
title('The predicted underlying function (sexp for 1. input + linear for 2. input )');

% ADDITIVE SQUARED EXPONENTIAL COVARIANCE FUNCTION
disp('Additive squared exponential covariance function')

% Covariance function for the second input variable
gpcf_s2 = gpcf_sexp('selectedVariables', 2,'lengthScale',0.5, ...
                    'lengthScale_prior', pt, 'magnSigma2', 0.15, ...
                    'magnSigma2_prior', pt);
% gpcf_s2 can be construted also as
%metric2 = metric_euclidean('components', {[2]},'lengthScale',[0.5], 'lengthScale_prior', pt);
%gpcf_s2 = gpcf_sexp('magnSigma2', 0.15, 'magnSigma2_prior', pt, 'metric', metric2);
gp = gp_set('lik', lik, 'cf', {gpcf_s1,gpcf_s2}, 'jitterSigma2', jitter);

% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Compute predictions in a grid using the MAP estimate
Eft_map = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure
set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,nxt,nxt));
hold on
plot3(x(:,1), x(:,2), y, '*')
xlabel('x_1'), ylabel('x_2')
title('The predicted underlying function (additive sexp)');

% SQUARED EXPONENTIAL COVARIANCE FUNCTION
disp('Squared exponential covariance function')

gpcf_s = gpcf_sexp('lengthScale', ones(1,nin), 'magnSigma2', 0.2^2, ...
                   'lengthScale_prior', pt, 'magnSigma2_prior', pt);
gp = gp_set('lik', lik, 'cf', {gpcf_s}, 'jitterSigma2', jitter);

% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Compute predictions in a grid using the MAP estimate
Eft_map = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure
set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,nxt,nxt));
hold on
plot3(x(:,1), x(:,2), y, '*')
xlabel('x_1'), ylabel('x_2')
title('The predicted underlying function (sexp)');

% ADDITIVE NEURAL NETWORK COVARIANCE FUNCTION
disp('Additive neural network covariance function');

gpcf_nn1 = gpcf_neuralnetwork('weightSigma2', 1, 'biasSigma2', 1, 'selectedVariables', [1], ...
                              'weightSigma2_prior', pt, 'biasSigma2_prior', pt);
gpcf_nn2 = gpcf_neuralnetwork('weightSigma2', 1, 'biasSigma2', 1, 'selectedVariables', [2], ...
                              'weightSigma2_prior', pt, 'biasSigma2_prior', pt);
gp = gp_set('lik', lik, 'cf', {gpcf_nn1,gpcf_nn2}, 'jitterSigma2', jitter);

% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Compute predictions in a grid using the MAP estimate
Eft_map = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure
set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,nxt,nxt));
hold on
plot3(x(:,1), x(:,2), y, '*')
xlabel('x_1'), ylabel('x_2')
title('The predicted underlying function (additive neural network)');

% NEURAL NETWORK COVARIANCE FUNCTION
disp('Neural network covariance function')

gpcf_nn = gpcf_neuralnetwork('weightSigma2', ones(1,nin), 'biasSigma2', 1, ...
                             'weightSigma2_prior', pt, 'biasSigma2_prior', pt);
gp = gp_set('lik', lik, 'cf', {gpcf_nn}, 'jitterSigma2', jitter);

% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);

% Compute predictions in a grid using the MAP estimate
Eft_map = gp_pred(gp, x, y, xt);

% Plot the prediction and data
figure
set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,nxt,nxt));
hold on
plot3(x(:,1), x(:,2), y, '*')
xlabel('x_1'), ylabel('x_2')
title('The predicted underlying function (neural network)');

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
%% Interpolating EBSD Data
%
%%
% In the section <EBSDDenoising.html Denoising> and <EBSDFilling.html
% Filling Missing Data> we have discussed how to work with noisy EBSD data
% the contained non indexed pixels. Hereby, we made the assumption that the
% grid before and after the operations is the same.
%
% In this section we explain how to interpolate an EBSD map at positions
% that do not belong to the grid. Lets us consider a simple example

mtexdata twins;

[grains, ebsd.grainId] = calcGrains(ebsd('indexed'));

% this command here is important :)
ebsd = ebsd.project2FundamentalRegion(grains);

plot(ebsd('indexed'),ebsd('indexed').orientations)

%%
% Now we can use the command <EBSD.interp.html |interp|> to interpolate the
% orientation at arbitrary coordinates |x| and |y|.

x = 30.5; y = 5.5;
e1 = interp(ebsd,x,y)

%%
% By default the command <EBSD.interp.html |interp|> performs inverse
% distance interpolation. This is different to

e2 = ebsd('xy',x,y)

%%
% which returns the nearest neighbour EBSD measurement. Lets have a look at
% the difference

angle(e1.orientations,e2.orientations)./degree

%% Change of the measurement grid
% The command <EBSD.interp.html |interp|> can be used to evaluate the EBSD
% map on a different grid, which might have higher or lower resolution or
% might even be rotated. Lets demonstrate this

% define a rotated coarse grid
omega = 5*degree;
[xmin, xmax, ymin, ymax] = ebsd.extent;
x = linspace(xmin-cos(omega)*ymax,xmax,100);
y = linspace(ymin-sin(omega)*xmax,ymax,50);
[x,y] = meshgrid(x,y);

xy = [cos(omega) -sin(omega); sin(omega) cos(omega) ] * [x(:),y(:)].';

% define the EBSD data set on this new grid
ebsdNewGrid = interp(ebsd,xy(1,:),xy(2,:))

% plot the regridded EBSD data set
plot(ebsdNewGrid('indexed'),ebsdNewGrid('indexed').orientations)

%%
% Note, that we have not rotated the EBSD data but only the grid. All
% orientations as well as the position of all grains remains unchanged.
%
% Another example is the change from a square to an hexagonal grid or vice
% versa. In this case the command <EBSD.interp.html |interp|> is
% implicitely called by the command <EBSD.gridify.html |gridify|>. In order
% to demonstrate this functionality we start by EBSD data on a hex grid

mtexdata ferrite silent

plot(ebsd,ebsd.orientations)

%%
% and resample the data on a square grid. To do so we first define a
% smaller square unit cell corresponding to the hexagonal unit cell

% define a square unit cell
hexUnitCell = abs(round(ebsd.unitCell,4));
minUnit = min(hexUnitCell(hexUnitCell>0));
squnitCell = minUnit * [-1 -1;-1 1; 1 1; 1 -1];

% use the square unit cell for gridify
ebsd = ebsd.gridify('unitCell',squnitCell);

plot(ebsd,ebsd.orientations)

author: mtex-toolbox

repo: mtex

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


text:
function thePartition=unrankPartition(rank,n)
%%UNRANKPARTITION Given the rank of a partiion of the integer n, when
%                 considering the partitions given in reverse lexicographic
%                 order, find the partition of the given rank. A partition
%                 of the integer n is a set of integers that sum to n.
%
%INPUTS: rank The integer rank of the partition. rank>0.
%           n The integer being partitioned.
%
%OUTPUTS: thePartition A numElementsX1  partition of n of the given rank.
%
%Let p(n,m) be the number of partitions of n with up to m parts. It is
%known that p(n,m)=p(n,m-1)+p(n-m,m). This relation is given in Chapter
%7.2.1.4 of [1]. The rank obtained here is derived using that recurrence
%relation.
%
%EXAMPLE:
%Here, we demonstrate that the unrankPartition function is consistent with
%the rankPartition function.
% n=20;
% [thePartition,theData]=getNextPartition(n);
% while(~isempty(thePartition))
%     assert(all(unrankPartition(rankPartition(thePartition),n)==thePartition))
%     [thePartition,theData]=getNextPartition(n,theData);
% end
%There should be no error indicating that the unranked partitions are all
%the same as the ranked partitions.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%September 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

thePartition=zeros(n,1);%We do not know how many parts to have...

rank=numberOfPartitions(n)-rank;

if(rank<0)
    error('The rank given is more than the total number of partitions of n')
end

curIdx=0;
nLeft=n;
while(rank>0)
    parts=1;
    numParts=numMPartitions(nLeft,parts,true);
    while(numParts<rank)
        parts=parts+1;
        numParts=numMPartitions(nLeft,parts,true);
    end
    if(numParts>rank)
        parts=parts-1;
        numParts=numMPartitions(nLeft,parts,true);
    end

    curIdx=curIdx+1;
    thePartition(curIdx)=parts+1;
    rank=rank-numParts;
    nLeft=nLeft-thePartition(curIdx);
end

while(nLeft>0)
    curIdx=curIdx+1;
    thePartition(curIdx)=1;
    nLeft=nLeft-1;
end

%Shrink to fit.
thePartition=thePartition(1:curIdx);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function skeleton_plot(cc, lab, varargin)
% SKELETON_PLOT  Plot a labelled segmentation and its skeleton
%
% SKELETON_PLOT(CC, LAB)
%
%   CC and LAB are the outputs from function skeleton_label().
%
%   This function plots the branches from the skeleton with the two end
%   points of each one (note than several branches can share the same end
%   point). If LAB is provided and not empty, the branches are plotted on
%   the labelled segmentation.
%
% SKELETON_PLOT(CC, LAB, ...)
%
%   It is possible to add more input arguments, that will be passed to the
%   plot functions. See 'help plot' for the available options.
%
% See also: skeleton_label.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2011 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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 3 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(1, Inf, nargin, 'struct'));
error(nargoutchk(0, 1, nargout, 'struct'));

% defaults
col = 'b';

% save hold status for later
HOLD = ishold;

% plot
if (nargin > 1 && ~isempty(lab))
    imagesc(lab)
    axis xy
    hold on
    col = 'w';
end

% loop every branch in the skeleton
for I = 1:cc.NumObjects
    % list of voxels in current branch
    br = cc.PixelIdxList{I};

    % plot skeleton branches
    if (length(cc.ImageSize) == 2)
        % convert image linear indices to image r, c, s indices
        [r, c] = ind2sub(cc.ImageSize, br);

        % plot branch in order
        plot(c, r, col, varargin{:})
        hold on

        % plot starting and end points
        plot(c([1 end]), r([1 end]), [col 'o'], varargin{:})
    elseif (length(cc.ImageSize) == 3)
        % convert image linear indices to image r, c, s indices
        [r, c, s] = ind2sub(cc.ImageSize, br);

        % plot branch in order
        plot3(c, r, s, 'w', varargin{:})
        hold on

        % plot starting and end points
        plot3(c([1 end]), r([1 end]), s([1 end]), [col 'o'], varargin{:})
    else
        error('Can only print 2D or 3D skeletons')
    end
end

% recover hold status
if HOLD
    hold on
else
    hold off
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [c1,c2] = simpleXover(p1,p2,bounds,Ops)
% Simple crossover takes two parents P1,P2 and performs simple single point
% crossover.
%
% function [c1,c2] = simpleXover(p1,p2,bounds,Ops)
% p1      - the first parent ( [solution string function value] )
% p2      - the second parent ( [solution string function value] )
% bounds  - the bounds matrix for the solution space
% Ops     - Options matrix for simple crossover [gen #SimpXovers].

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

numVar = size(p1,2)-1; 			% Get the number of variables
% Pick a cut point randomly from 1-number of vars
cPoint = round(rand * (numVar-2)) + 1;

c1 = [p1(1:cPoint) p2(cPoint+1:numVar+1)]; % Create the children
c2 = [p2(1:cPoint) p1(cPoint+1:numVar+1)];

author: HuangCongQing


text:
function BuildBCMaps3D()

% function BuildBCMaps3D
% Purpose: Build specialized nodal maps for various types of boundary conditions,
%           specified in BCType.

Globals3D;

% create label of face nodes with boundary types from BCType
bct    = BCType';
bnodes = ones(Nfp, 1)*bct(:)'; bnodes = bnodes(:);

% find location of boundary nodes in face and volume node lists
mapI = find(bnodes==In);           vmapI = vmapM(mapI);
mapO = find(bnodes==Out);	   vmapO = vmapM(mapO);
mapW = find(bnodes==Wall);         vmapW = vmapM(mapW);
mapF = find(bnodes==Far);          vmapF = vmapM(mapF);
mapC = find(bnodes==Cyl);          vmapC = vmapM(mapC);
mapD = find(bnodes==Dirichlet);    vmapD = vmapM(mapD);
mapN = find(bnodes==Neuman);       vmapN = vmapM(mapN);
mapS = find(bnodes==Slip);         vmapS = vmapM(mapS);
return;

repo: Aero-matlab


text:
function r=rot4z(theta,xscale,yscale,zscale)
    % r = rot4z(theta)
    %
    % rotz produces a 4x4 rotation matrix representing
    % a rotation by theta radians about the z axis.
    %
    %	Argument definitions:
    %
    %	theta = rotation angle in radians
    if nargin == 1
        c = cos(theta);
        s = sin(theta);
        r = [c -s  0 0;
            s  c  0 0;
            0  0  1 0;
            0  0  0 1];
    else
        c = cos(theta);
        s = sin(theta);
        r = [c -s  0 xscale;
            s  c  0 yscale;
            0  0  1 zscale;
            0 0 0  1];
    end

author: CelsoReyes


text:
function y = cvx_subsasgn( y, varargin )
temp.type = '()';
temp.subs = varargin(1:end-1);
y = subsasgn( y, temp, varargin{end} );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function varargout = trace(varargin)
% TRACE integral of a CHEBFUN2 along its diagonal
%
% TRACE(f) is the integral of function f(x,x).
%
% See also DIAG.

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

[varargout{1:nargout}] = trace@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = spm_dartel_jacobian(job)
% Generate Jacobian determinant fields
% FORMAT spm_dartel_jacobian(job)
% job.flowfields - Filenames of flowfields
% job.K          - 2^K timesteps are used
%
% Note that K needs to be reasonably large in order to obtain reasonable
% Jacobian determinant fields.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_dartel_jacobian.m 5506 2013-05-14 17:13:43Z john $

PU = job.flowfields;
K  = job.K;

spm_progress_bar('Init',numel(PU),'Creating Jacobian det fields','Number complete');
for i=1:numel(PU),
    NU = nifti(PU{i});
    [pth,nam,ext,num] = spm_fileparts(NU.dat.fname);
    [y,dt] = spm_dartel_integrate(NU.dat,[1 0], K);
    clear y

    NO = NU;
    NO.dat.fname=fullfile(pth,['jac_' nam(3:end) ext]);
    NO.dat.scl_slope = 1.0;
    NO.dat.scl_inter = 0.0;
    NO.dat.dtype     = 'float32-le';
    NO.dat.dim = NU.dat.dim(1:3);
    NO.mat  = NU.mat;
    NO.mat0 = NU.mat;
    NO.mat_intent  = 'Aligned';
    NO.mat0_intent = 'Aligned';
    NO.descrip = 'Dartel Jacobian';
    create(NO);
    NO.dat(:,:,:)=dt;
    spm_progress_bar('Set',i);
end;
spm_progress_bar('Clear');

PU = job.flowfields;
out.files = cell(numel(PU),1);
for i=1:numel(PU),
    [pth,nam,ext] = fileparts(PU{i});
    fname         = fullfile(pth,['jac_' nam(3:end) ext]);
    out.files{i}  = fname;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dark_mode(use_dark_mode)
  % DARK_MODE  Switch command and figures windows to dark background with light
  % text.
  %
  % dark_mode
  % dark_mode(flag)
  %
  % Optional input:
  %   flag  whether to switch to dark mode {false} or back to default (true)
  %
  % Examples:
  %   % Place this in startup.m to add a timer that checks every 5mins whether
  %   % its ≥6am or ≥6pm and toggles dark mode accordingly.
  %   dmt = timer('TimerFcn',@(x,y) dark_mode(mod(clock*[0;0;0;1;0;0]+6,24)<12),'Period',60*5,'ExecutionMode','fixedDelay');
  %   start(dmt);
  %
  %   % Detect if Mac OS X is dark and match
  %   dark_mode(system('defaults read -g AppleInterfaceStyle 2&>/dev/null')==0);
  %

  % Adapted from:
  % https://github.com/benhager/solarized-matlab/blob/master/setupSolarized.m

  if ~exist('use_dark_mode','var') || use_dark_mode
    sys = false;
    txc = [1 1 1];
    bgc = [30 30 30]/255;
    fbg = [1 1 1].*40/255;
    alc = [1 1 1];
    ac = [30 30 30]/255;
    str = [254.1188  106.3358  255.0000]/255;
    ustr = [203  75  22] / 255;
    hyp = [ 38 139 210] / 255;
    cmt = [ 54 139  54] / 255;
    scmd = [198 160  20] / 255;
    errs = [230 40 40]/255;
  else
    sys = true;
    txc = [0 0 0];
    bgc = [1 1 1];
    fbg = get(0,'factoryFigureColor');
    alc = 'k';
    ac = 'w';
    str = [160  32 240]/255;
    ustr = [178   0   0]/255;
    hyp = [0 0 1];
    cmt = [ 34 139  34] / 255;
    scmd = [178 140   0] / 255;
    errs = [230 0 0]/255;
  end
  kwd = hyp;
  warn = ustr;

  % Main colors
  com.mathworks.services.Prefs.setBooleanPref('ColorsUseSystem',sys);
  com.mathworks.services.Prefs.setColorPref('ColorsText',java.awt.Color(txc(1), txc(2), txc(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsText');
  com.mathworks.services.Prefs.setColorPref('ColorsBackground',java.awt.Color(bgc(1), bgc(2), bgc(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsBackground');
  com.mathworks.services.Prefs.setColorPref('Colors_M_Strings',java.awt.Color(str(1), str(2), str(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Strings');
  com.mathworks.services.Prefs.setColorPref('Colors_M_UnterminatedStrings',java.awt.Color(ustr(1), ustr(2), ustr(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_UnterminatedStrings');
  com.mathworks.services.Prefs.setColorPref('Colors_M_Keywords',java.awt.Color(kwd(1), kwd(2), kwd(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Keywords');
  com.mathworks.services.Prefs.setColorPref('Colors_M_Comments',java.awt.Color(cmt(1), cmt(2), cmt(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Comments');
  com.mathworks.services.Prefs.setColorPref('Colors_M_SystemCommands',java.awt.Color(scmd(1), scmd(2), scmd(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_SystemCommands');
  com.mathworks.services.Prefs.setColorPref('Colors_M_Errors',java.awt.Color(errs(1), errs(2), errs(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Errors');
  com.mathworks.services.Prefs.setColorPref('Colors_HTML_HTMLLinks',java.awt.Color(hyp(1), hyp(2), hyp(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_HTML_HTMLLinks');
  com.mathworks.services.Prefs.setColorPref('Colors_M_Warnings',java.awt.Color(warn(1), warn(2), warn(3)));
  com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Warnings');
  %com.mathworks.services.Prefs.setBooleanPref('ColorsUseMLintAutoFixBackground',afhb);
  %com.mathworks.services.Prefs.setColorPref('ColorsMLintAutoFixBackground',java.awt.Color(afh(1), afh(2), afh(3)));
  %com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsMLintAutoFixBackground');
  %com.mathworks.services.Prefs.setColorPref('Editor.NonlocalVariableHighlighting.TextColor',java.awt.Color(vwss(1), vwss(2), vwss(3)));
  %com.mathworks.services.ColorPrefs.notifyColorListeners('Editor.NonlocalVariableHighlighting.TextColor');
  %com.mathworks.services.Prefs.setBooleanPref('EditorCodepadHighVisible',hsb);
  %com.mathworks.services.Prefs.setColorPref('Editorhighlight-lines', java.awt.Color(hsbc(1), hsbc(2), hsbc(3)));
  %com.mathworks.services.Prefs.setBooleanPref('Editorhighlight-caret-row-boolean',hclb);
  %com.mathworks.services.Prefs.setColorPref('Editorhighlight-caret-row-boolean-color',java.awt.Color(hcl(1), hcl(2), hcl(3)));
  %com.mathworks.services.ColorPrefs.notifyColorListeners('Editorhighlight-caret-row-boolean-color');
  %com.mathworks.services.Prefs.setBooleanPref('EditorShowLineNumbers',slnb);
  %com.mathworks.services.Prefs.setBooleanPref('EditorRightTextLineVisible',rtlb);

  % Figure windows
  set(0,'defaultfigurecolor',fbg)
  arrayfun(@(f) set(f,'Color',fbg),get(groot, 'Children'));
  % axes labels
  set(groot,{'DefaultAxesXColor','DefaultAxesYColor','DefaultAxesZColor','DefaultAxesColor','DefaultTextColor','DefaultLineColor'},{alc,alc,alc,ac,alc,alc})
  arrayfun(@(a) set(a,{'Xcolor','YColor','ZColor','Color'},{alc,alc,alc,ac}),findall(0,'type','axes'));
  arrayfun(@(t) set(t,'Color',alc),findall(0,'type','text'));

end

author: alecjacobson

repo: gptoolbox


text:
function ouotput = MRCG_features(sig, sampFreq)
% This function computes MRCG features

beta = 1000 / sqrt( sum(sig .^ 2) / length(sig) );
sig = sig .* beta;
sig = reshape(sig, length(sig), 1);
g = gammatone(sig, 64, [50 8000], sampFreq); % Gammatone filterbank responses

cochlea1 = log10(cochleagram(g,sampFreq*0.025,sampFreq*0.010));
cochlea2 = log10(cochleagram(g,sampFreq*0.200,sampFreq*0.010));
cochlea1 = cochlea1(:,:);
cochlea2 = cochlea2(:,:);

cochlea3  = get_avg(cochlea1,5,5);
cochlea4  = get_avg(cochlea1,11,11);
all_cochleas = [cochlea1; cochlea2; cochlea3; cochlea4];

del = deltas(all_cochleas);
ddel = deltas(deltas(all_cochleas,5),5);

ouotput = [all_cochleas;del;ddel];

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%mmf writer
function [keyframes,nframes]=mmfwriter2(inputfile,outputfile,obj,firstframe,lastframe)

%obj=VideoReader(inputfile);
foutput=fopen(outputfile,'wb');
%load text
ftexto=fopen('/groups/visitors/home/albam/Desktop/mmfreader/text.text');

texto=[];
for i=1:15
linea=fgets(ftexto);
texto=[texto linea];
end

%write header text
fwrite(foutput,texto,'char')
% end of the header
fwrite(foutput,0,'char')
%id_code
d=hex2dec('a3d2d45d');
fwrite(foutput,d,'ulong')
%header size
fwrite(foutput,10240,'int')
%calculate keyframes
keyframes=fix((lastframe-firstframe+1)/180)+1;
%keyframe
fwrite(foutput,keyframes,'int')
%other recording parameters TODO
% zero padded to reach header size
for i=ftell(foutput)+1:10240
fwrite(foutput,0,'uint8');
end
%calculate nframes
%
nframes=zeros(1,keyframes);
nframes(1:end-1)=180;
nframes(end)=(lastframe-firstframe+1)-(keyframes-1)*180;

isfirst=1;
% for the stacks
for stack=1:keyframes

%Image Stacks
begining_header=ftell(foutput);
%idcode
d=hex2dec('bb67ca20');
fwrite(foutput,d,'ulong')
%header size
fwrite(foutput,512,'int')
%total size of stack in disk TODO
fwrite(foutput,0,'int')
%nframes in this stack default=180 but the last one can be smaller
fwrite(foutput,nframes(stack),'int');
% zero padded to reach header size
for i=ftell(foutput)-begining_header+1:512
fwrite(foutput,0,'uint8');
end

%Backgorund image
begining_header=ftell(foutput);
%header size
fwrite(foutput,112,'int')
%imageid
fwrite(foutput,0,'int')
%numchannels
fwrite(foutput,1,'int')
%alphachannel
fwrite(foutput,0,'int')
%depth
fwrite(foutput,8,'int')
%colormodel
fwrite(foutput,[0 0 0 0],'char')
%channel seq
fwrite(foutput,[0 0 0 0],'char')
%data order
fwrite(foutput,0,'int')
% origin
fwrite(foutput,0,'int')
%align should 4 or 8 TODO
fwrite(foutput,0,'int')
%width
fwrite(foutput,obj.width,'int')
%height
fwrite(foutput,obj.height,'int')
%roi 0=> full image
fwrite(foutput,0,'int')
%mask roi
fwrite(foutput,0,'int')
%im id
fwrite(foutput,0,'int')
%tile info
fwrite(foutput,0,'int')
%images size (heightxwidth
fwrite(foutput,obj.width*obj.height,'int')
%imagedata
fwrite(foutput,0,'int')
%widthstep
fwrite(foutput,obj.width,'int')
%bordermode
fwrite(foutput,0,'int')
%borderconst
fwrite(foutput,0,'int')
%imagedataorigin
fwrite(foutput,0,'int')
% zero padded to reach header size
for i=ftell(foutput)-begining_header+1:112
fwrite(foutput,0,'uint8');
end
s=0;
im=zeros(obj.width,obj.height,nframes(stack));
if isfirst==1
    for i=firstframe:firstframe+nframes(stack);
        s=s+1;
        imrgbtemp=read(obj,i);
        imrgbtemp=double(imrgbtemp);
        im(:,:,s)=(sum(imrgbtemp,3)/3)';
    end
    isfirst=0;
else
    for i=firstframe+sum(nframes(1:stack-1)):firstframe+sum(nframes(1:stack))
        s=s+1;
        imrgbtemp=read(obj,i);
        imrgbtemp=double(imrgbtemp);
        im(:,:,s)=(sum(imrgbtemp,3)/3)';
    end
end

%%creting mean backgfix
bkgim=median(im,3);
%writing background image
fwrite(foutput,reshape(bkgim',1,obj.width*obj.height))

%Extracting larvae out of the images

 frames=zeros(size(im));
 SE=strel('disk',8);
 bounding=cell(1,nframes(stack));

for i=1:nframes(stack)
   begining_header=ftell(foutput);
   %id code
   d=hex2dec('f80921af');
   fwrite(foutput,d,'ulong')
   %header size
   fwrite(foutput,1024,'int')
   %depth
   fwrite(foutput,8,'int')
   %nchannels
   fwrite(foutput,1,'int')

   %calculating numims
    frames(:,:,i)=abs(im(:,:,i)-bkgim);
    %threshold automatically set
    %level(i)=graythresh(frames(:,:,i));
     bw=zeros(size(frames(:,:,1)));
     bw(frames(:,:,i)*10>max(max(frames(:,:,i))))=1;
     bw2=bwareaopen(bw,7);
     dilated=imdilate(bw2,SE);
     cc=bwconncomp(dilated);

     %numims
     fwrite(foutput,cc.NumObjects,'int')
     % zero padded to reach header size
        for k=ftell(foutput)-begining_header+1:1024
        fwrite(foutput,0,'uint8');
        end
%     %Creating bounding box
     STATS=regionprops(cc,'BoundingBox');
      bounding{i}=zeros(cc.NumObjects,4);
     for j=1:cc.NumObjects
         bounding{i}(j,:)=STATS(j,1).BoundingBox;
         %temp=zeros(fix(bounding{i}(j,4)),fix(bounding{i}(j,3)));
         if fix(bounding{i}(j,2))==0 && fix(bounding{i}(j,1))==0
            temp=im(1:fix(bounding{i}(j,4)),1:fix(bounding{i}(j,3)),i);
        elseif fix(bounding{i}(j,2))==0
            temp=im(1:fix(bounding{i}(j,4)),fix(bounding{i}(j,1)):fix(bounding{i}(j,1))+fix(bounding{i}(j,3))-1,i);
        elseif fix(bounding{i}(j,1))==0
            temp=im(fix(bounding{i}(j,2)):fix(bounding{i}(j,2))+fix(bounding{i}(j,4))-1,1:fix(bounding{i}(j,3)),i);
        else
            temp=im(fix(bounding{i}(j,2)):fix(bounding{i}(j,2))+fix(bounding{i}(j,4))-1,fix(bounding{i}(j,1)):fix(bounding{i}(j,1))+fix(bounding{i}(j,3))-1,i);
         end
        fwrite(foutput,fix(bounding{i}(j,:)),'int')
        fwrite(foutput,reshape(temp',1,fix(bounding{i}(j,3))*fix(bounding{i}(j,4))))
      end

end

end
fclose(foutput)

author: kristinbranson

repo: JAABA


text:
% This function can read in all the information from the HTK model file
% with specifications of AURORA2 baseline system.

% AUTHOR: XIAO XIONG
% CREATED: 7 Jul, 2006
% LAST MODIFIED: 7 Jul, 2006
function model = readHTKmodel(file_name,digits_file)

% digits = {'sil','one','two','three','four','five','six','seven','eight','nine','zero','oh'};
if nargin < 2
    digits_file = 'E:\Research\Workspace\Projects\Normalization\AURORA2\LIB\words4';
end
if nargin < 1
    file_name = 'E:\Research\Workspace\Projects\Normalization\AURORA2\MODEL\MFCC_0DA\hmm20\models';
end
silence = {'sil','sp','silst'};
FID = fopen(file_name,'r');
digits = my_cat(digits_file);
digits2 = {};
for i=1:length(digits)
    if length(digits{i})>0
        digits2{end+1} = digits{i};
    end
end
digits = digits2;

% Read in the global parameters
fsearch('~o',FID);
tmp = textscan(FID,'%s %d %d',1);
tmp = fgetl(FID); tmp = fgetl(FID);
idx = regexp(tmp,'<');
Dimension = str2num(tmp(idx(1)+10:idx(2)-1));        % Dimension of the feature vector
idx = regexp(tmp,'<');
feature_type = tmp(idx(3)+1:idx(4)-2);         % type of the feature vector
cov_type = tmp(idx(4)+1:end-1);;             % type of the covariance matrix

% Read in the silst model
endOfFile = fsearch('~s "silst"',FID);
if endOfFile==0
    model.silst.NUMSTATES = 1;
    tmp = textscan(FID,'%s %d',1);
    if strcmp(tmp{1},'<MEAN>')
        model.silst.NUMMIXES = 1;
        fseek(FID,-7,'cof');
    else
        model.silst.NUMMIXES = tmp{2};
    end
    % model.('silst').NUMMIXES = tmp{2};
    for i=1:model.silst.NUMMIXES
        if model.silst.NUMMIXES>1
            tmp = textscan(FID,'%s %d',1);
            tmp = textscan(FID,'%f ',1);
            model.silst.weight(i) = tmp{1};
        else
            model.silst.weight(i) = 1;
        end

        tmp = textscan(FID,'%s %d',1);
        tmp = textscan(FID,'%f ',Dimension);
        model.silst.mean(:,i) = tmp{1};

        tmp = textscan(FID,'%s %d',1);
        tmp = textscan(FID,'%f ',Dimension);
        model.silst.var(:,i) = tmp{1};

        tmp = textscan(FID,'%s %f',1);
        model.silst.GCONST(:,i) = tmp{2};
    end
    model.silst.weight = model.silst.weight(:);
end

% Read in the sp model
frewind(FID);
endOfFile = fsearch(sprintf('~h "sp"'),FID);
if endOfFile==0
    fsearch(sprintf('<TRANSP> 3'),FID);
    model.sp = model.silst;
    model.sp.NUMSTATES = 1;
    for j=1:model.sp.NUMSTATES+2
        tmp = textscan(FID,'%f ',3);
        model.sp.TRANSP(j,:) = tmp{1};
    end
end

% Read in the digit models
for i=1:length(digits)
    frewind(FID);
    fsearch(sprintf('~h "%s"',digits{i}),FID);
    tmp = textscan(FID,'%s',1);
    tmp = textscan(FID,'%s %d',1);

    if strcmp(digits{i},'@')
        currDigit = 'AA';
    else
        currDigit = digits{i};
    end
    model.(currDigit).NUMSTATES = tmp{2}-2;

    for j=1:model.(currDigit).NUMSTATES
        tmp = textscan(FID,'%s %d',1);
        tmp = textscan(FID,'%s',1);
        if strcmp(tmp{1},'~s')  % this is a silst state
            tmp = textscan(FID,'%s',1);
            model.(currDigit).NUMMIXES(j) = model.silst.NUMMIXES;
            if model.silst.NUMMIXES == 1
                model.(currDigit).mean(:,1,j) = model.silst.mean;
                model.(currDigit).var(:,1,j) = model.silst.var;
                model.(currDigit).weight(1,j) = model.silst.weight;
            else
                model.(currDigit).mean(:,:,j) = model.silst.mean;
                model.(currDigit).var(:,:,j) = model.silst.var;
                model.(currDigit).weight(:,j) = model.silst.weight;
            end
            continue;
        elseif strcmp(tmp{1},'<MEAN>')
            model.(currDigit).NUMMIXES(j) = 1;
            fseek(FID,-7,'cof');
        else
            tmp = textscan(FID,'%d',1);
            model.(currDigit).NUMMIXES(j) = tmp{1};
        end
        for k=1:model.(currDigit).NUMMIXES(j)
            if model.(currDigit).NUMMIXES(j)>1
                tmp = textscan(FID,'%s %d',1);
                tmp = textscan(FID,'%f ',1);
                model.(currDigit).weight(k,j) = tmp{1};
            else
                model.(currDigit).weight(k,j) = 1;
            end
            tmp = textscan(FID,'%s %d',1);
            tmp = textscan(FID,'%f ',Dimension);
            model.(currDigit).mean(:,k,j) = tmp{1};

            tmp = textscan(FID,'%s %d',1);
            tmp = textscan(FID,'%f ',Dimension);
            model.(currDigit).var(:,k,j) = tmp{1};

            tmp = textscan(FID,'%s %f',1);
            model.(currDigit).GCONST(k,j) = tmp{2};
        end
    end
    % Read in the trainsit matrix
    tmp = textscan(FID,'%s %d',1);
    for j=1:model.(currDigit).NUMSTATES+2
        tmp = textscan(FID,'%f ',model.(currDigit).NUMSTATES+2);
        model.(currDigit).TRANSP(j,:) = tmp{1};
    end
end
fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function op = prox_dualize( dualProx, NEG )
%PROX_DUALIZE   Define a proximity function by its dual
%    OP = PROX_DUALIZE( dualOp ) returns an operator implementing the
%    dual of the function dualProx. You can verify they are duals
%    via test_proxPair( dualOp, OP ).
%
%    OP = PROX_DUALIZE( dualOp, 'neg' )
%    OP = PROX_DUALIZE( dualOp, 'negative' )
%       will return the scaled dual of dualOp; that is,
%       dualOp(x) and OP(-x) are duals.
%       The negative is useful because this is the version TFOCS
%       expects for the SCD formulation.
%       For 1-homogenous functions (e.g. norms), this has no effect,
%       since ||x|| = ||-x||.
%
% Warning: if you can calculate the dual function explicitly,
%   it is likely more computationally efficient to do so, rather
%   than rely on this code. This code requires some tricks, some
%   of which can sometimes be expensive; also, it will call dualOp
%   so it is at least as slow as dualOp; and it may require high
%   precision from dualOp. This code can break down if dualOp
%   is not numerically stable.

if nargin < 2, NEG = ''; end
if strcmpi(NEG,'neg') || strcmpi(NEG,'negative')
    op = @(varargin)dualize(dualProx,-1,varargin{:} );
else
    op = @(varargin)dualize(dualProx,1,varargin{:} );
end

function [ v, x ] = dualize( dualProx, scale, x, t )
vec     = @(x) x(:);
myDot   = @(x,y) x(:)'*y(:);
if scale == -1
    x   = -x;
end
switch nargin,
    case 3
        if nargout == 2,
            error( 'This function is not differentiable.'  );
        else
            % This case is a bit tricky...
            %   If the function is non-differentiable, then standard exact
            %   penalty function results tell us that for a sufficiently
            %   small stepsize, we can remove the effect of smoothing.
            %   In other cases, we don't have an exact value, but we hope this
            %   is a reasonable approximation.

%             t       = 1e-15;
%             [~,x2]  = dualProx( x/t, 1/t );
%             v       = myDot(x,x2) - dualProx( x2 );

            % However, some functions break down when 1/t is huge
            % So we will slowly decrease it
            vOld    = Inf;
            t       = 1e-5;
            ok      = false;
            iter    = 0;
            while ~ok && t > eps
                [~,x2]  = dualProx( x/t, 1/t );
                v       = myDot(x,x2) - dualProx( x2 );
                if abs(v-vOld)/max( 1e-10, abs(v) ) < 1e-4
                    % due to exact penalty, we expect that
                    %   for t < t_cutoff, v=vOld up to machine accuracy
                    ok = true;
                else
                    t   = t/10;
                    vOld = v;
                    iter = iter + 1;
                    %fprintf('%d and v is %.2e\n', iter, v );
                end
            end

        end

    case 4
        % This is exact.
        [ignore,x2]  = dualProx( x/t, 1/t );
        x1      = x - t*x2; % Moreau's identity, equation (8.1) in the user guide
        v       = myDot(x1,x2) - dualProx( x2 );

        % If we think it is an indicator function, then round down to zero:
        if abs(v) < 100*eps,
            v  = 0;
        end
        x   = scale*x1;

    otherwise,
        error( 'Not enough arguments.' );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
%TestSuiteInDir Test suite requiring temporary directory change
%   The TestSuiteInDir class defines a test suite that has to be run by first
%   changing to a specified directory.
%
%   The setUp method adds the starting directory to the path and then uses cd to
%   change into the specified directory.  The tearDown method restores the
%   original path and directory.
%
%   TestSuiteInDir methods:
%       TestSuiteInDir  - Constructor
%       gatherTestCases - Add test cases found in the target directory
%
%   See also TestSuite

%   Steven L. Eddins
%   Copyright 2009 The MathWorks, Inc.

classdef TestSuiteInDir < TestSuite & TestComponentInDir

    methods
        function self = TestSuiteInDir(testDirectory)
            %TestCaseInDir Constructor
            %   TestCaseInDir(testName, testDirectory) constructs a test case
            %   using the specified name and located in the specified directory.
            self = self@TestComponentInDir(testDirectory);

            if strcmp(testDirectory, '.')
                self.Name = pwd;
                self.Location = pwd;
            else
                [pathstr, name] = fileparts(testDirectory);
                self.Name = name;
                self.Location = testDirectory;
            end
        end

        function gatherTestCases(self)
            %gatherTestCases Add test cases found in the target directory
            %   suite.gatherTestCases() automaticall finds all the test cases in
            %   the directory specified in the constructor call and adds them to
            %   the suite.
            current_dir = pwd;
            c = onCleanup(@() cd(current_dir));

            cd(self.TestDirectory);
            tmp = TestSuite.fromPwd();
            self.TestComponents = tmp.TestComponents;
        end
    end
end

author: ricedsp


function dataout=modeWPE_NLMS_BSS(nummics, numrefs, datain)
%
% Perform wpe, nlms, and bss.
% nummics:              no. of mic channels
% numrefs:              no. of reference channels
% datain:               input data
% dataout:              output data
%

config;

%% perform wpe
addpath('wpe_v1.33');

mic=datain(:, 1:2);
ref=datain(:, 3);

cfgs='wpe_v1.33/settings/local.m';
early=wpe(mic, cfgs);

%% perform nlms
addpath('Speex-AEC-matlab-master');

%    Usage:
%
%       speex_mdf_out = speex_mdf(Fs, u, d, filter_length, frame_size, dbg_var_name);
%
%       Fs                  sample rate
%       u                   speaker signal, column vector in range [-1; 1]
%       d                   microphone signal, column vector in range [-1; 1]
%       filter_length       typically 250ms, i.e. 4096 @ 16k FS
%                           must be a power of 2
%       frame_size          typically 8ms, i.e. 128 @ 16k Fs
%                           must be a power of 2
%       dbg_var_name        internal state variable name to trace.
%                           Default: 'st.leak_estimate'.
%
%    Jonathan Rouach <jonr@waves.com>
%
filter_length=stftshift*AEC_FLEN;
frame_size=stftshift;
nearend1=speex_mdf(fs, ref, early(:, 1), filter_length, frame_size);
nearend2=speex_mdf(fs, ref, early(:, 2), filter_length, frame_size);
nearend=[nearend1.e, nearend2.e];

%% perform stft
addpath('stft2');

M=nummics;
N=M;

Xtf=cell(M, 1);
for m=1:M
    Xtf{m}=stft(nearend(:, m), stftshift, fftsize, false);
end
[K, T]=size(Xtf{1});

Ytf=cell(M, 1);
for m=1:M
    Ytf{m}=zeros(K, T);
end

%% params go to config

%% space for bss
% the weighted correlation matrices
C1=cell(K, 1);
C2=cell(K, 1);
for k=1:K
    C1{k}=STABLE_EPS*eye(M, M);
    C2{k}=STABLE_EPS*eye(M, M);
end

% demixing matrices
Demix=cell(K, 1);
for k=1:K
    Demix{k}=eye(N, M);
end

%% perform bss
for tau=1:T
    %% perform bss
    Bssout=zeros(K, M);

    %
    % calculate nonlinearity
    %
    phi1=0;
    phi2=0;

    for k=1:K
        x=zeros(M, 1);
        for m=1:M
            x(m)=Xtf{m}(k, tau);
        end

        y=Demix{k}*x;
        % output data
        Bssout(k, :)=y.';

        phi1=phi1+abs(y(1))^2;
        phi2=phi2+abs(y(2))^2;
    end

    phi1=(1-BF_FORGET)*(phi1+VAR_BIAS)^((GAMMA-2)/2);
    phi2=(1-BF_FORGET)*(phi2+VAR_BIAS)^((GAMMA-2)/2);

    % update the demixing matrices
    for k=1:K
        %
        % accumulate the weighted correlation
        %
        x=zeros(M, 1);
        for m=1:M
            x(m)=Xtf{m}(k, tau);
        end

        C1{k}=BF_FORGET*C1{k}+phi1*(x*x');
        C2{k}=BF_FORGET*C2{k}+phi2*(x*x');

        %
        % solve gev problem
        %
        [Ev, Ed]=eig(C2{k}+BF_DIAGLOAD*eye(M, M), C1{k}+BF_DIAGLOAD*eye(M, M));
        if Ed(1, 1)>=Ed(2, 2)
            e1=Ev(:, 1);
            e2=Ev(:, 2);
        else
            e1=Ev(:, 2);
            e2=Ev(:, 1);
        end

        D=[e1'; e2'];

        %
        % solve the scaling ambiguity
        %
        A=inv(D);

        if abs(A(1, 1))>=abs(A(2, 1))
            a1=A(1, 1);
        else
            a1=A(2, 1);
        end

        if abs(A(2, 2))>=abs(A(1, 2))
            a2=A(2, 2);
        else
            a2=A(1, 2);
        end

        D=diag([a1; a2])*D;
        Demix{k}=D;
    end

    for m=1:M
        Ytf{m}(:, tau)=Bssout(:, m);
    end
end

%% perform istft and output signal
dataout=zeros(dataLength(T, stftshift, fftsize ), N);
for n=1:N
    dataout(:, n)=istft(Ytf{n}, stftshift, false);
end

end


text:
function [mach] = ftps2mach(ftps)
% Convert speed from feet per second to mach number (at STP!)
% Chad A. Greene 2012
mach = ftps*0.0008886297376093 ;

author: Sable

repo: mcbench-benchmarks


text:
function pageplot(mode,option)
% PAGEPLOT Update all graphs on current page.
%          PAGEPLOT manages the updating of all the currently viewed
%          axes objects.

% Author: Craig Borghesani
% Date: 10/18/94
% Revised: 11/13/94
% Copyright (c) 1999, Prentice-Hall

% obtain handle information
f = gcf;
ui_data = get(f,'userdata');
ui_han = ui_data{1};
sys_axes = ui_han(12);
cur_sys = get(ui_han(30),'userdata');
cur_page = get(ui_han(31),'userdata');
cur_axs = get(ui_han(32),'userdata');
font_size = get(ui_han(42),'userdata');
page_setup = get(ui_han(34),'userdata');
page_layout = get(ui_han(35),'userdata');

plant_mat = get(ui_han(3+cur_sys),'userdata');
for_mat = get(ui_han(6+cur_sys),'userdata');
bac_mat = get(ui_han(9+cur_sys),'userdata');
sys_res = get(ui_han(cur_sys),'userdata');

plot_type = page_setup(:,cur_page); plot_type(find(plot_type==0)) = [];
plot_hand = page_layout(:,cur_page); plot_hand(find(plot_hand==0)) = [];

if mode == 0, % update all visible plots
 freqplot(plot_hand(find(plot_type>0 & plot_type<5)),option);
 rootplot(plot_hand(find(plot_type==5)),option);
 gainplot(plot_hand(find(plot_type>5 & plot_type<8)),option)
 timeplot(plot_hand(find(plot_type==8)),option);

elseif mode == 1, % only frequency response plots
 freqplot(plot_hand(find(plot_type>0 & plot_type<5)),option);

elseif mode == 2, % only root locus plots
 rootplot(plot_hand(find(plot_type==5)),option);
 gainplot(plot_hand(find(plot_type>5 & plot_type<8)),option)

elseif mode == 3, % only time response plots
 timeplot(plot_hand(find(plot_type==8)),option);

elseif mode == 4, % only gain plots
 gainplot(plot_hand(find(plot_type>5 & plot_type<8)),option)

end

% display system information
% option == 13 - switching systems
% option == 14 - changing pages
% option == 15 - changing page layouts
% option == [] - changing transfer function

if ~length(option) | any(option == [13:15]),

 data = get(cur_axs,'userdata');
 environ = data(1);
 if all(environ~=[5:7]),
  systemin(environ, sys_axes, font_size, sys_res, plant_mat, for_mat, bac_mat);
 end

end


text:
BTB.MatDir= fullfile(BTB.DataDir, 'demoMat');
filename= fullfile('VPean_10_07_26', 'NIRS', 'real_movementVPean');

%% Load data
[cnt, mrk, mnt] = file_loadNIRSMatlab(filename, 'Signal','oxy');

%% define classes according to makers
stimDef= {1,2;'left','right'};
mrk = mrk_defineClasses(mrk, stimDef);

%% intervalls for display and segmentation, and fequencies for filter
ival_base=  [-1000 0];
ival_epo= [-1000 15000];
ival_scalps= [0:2500:12500];
clab=[];

%% segmentation and baseline correction
epo= proc_segmentation(cnt, mrk, ival_epo);
epo= proc_baseline(epo, ival_base);

%% r-values
epo_r= proc_rSquareSigned(epo);

%% display
fig_set(1)
H= grid_plot(epo, mnt, defopt_erps);
grid_addBars(epo_r, 'HScale',H.scale);

fig_set(2);
H= plot_scalpEvolution(epo, mnt, ival_scalps, ...
                       defopt_scalp_erp, ...
                       'ExtrapolateToMean', 1);

fig_set(4, 'Resize',[1 0.5]);
H= plot_scalpEvolution(epo_r, mnt, ival_scalps, ...
                       defopt_scalp_r);

author: bbci


text:
function range = p02_range ( )

%*****************************************************************************80
%
%% P02_RANGE returns an interval bounding the root for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = -10.0;
  range(2) = 100.0;

  return
end

author: johannesgerer


text:
function r=roteu2ro(m,t)
%ROTEU2QR converts a sequence of Euler angles to a real unit quaternion
% Inputs:
%
%     M(1,n)   a string of n characters from the set {'x','y','z'}
%              or, equivalently, a vector whose elements are 1, 2, or 3
%     T(n,1)   n rotation angles. A positive rotation is clockwise if
%              looking along the axis away from the origin.
%
% Outputs:
%
%     R(3,3)   Input rotation matrix
%              Plots a diagram if no output specified
%
% The string M specifies the axes about which the rotations are performed.
% You cannot have the same axis in adjacent positions and so there are 12
% possibilities. Common ones are "ZXZ" and "ZYX". A positive rotation is clockwise
% if looking along the axis away from the origin; thus a rotation of +pi/2
% around Z rotates [1 0 0]' to [0 1 0]'.
%
% Inverse conversion: If m has length 3 with adjacent characters distinct,
%                     then rotro2eu(m,roteu2ro(m,t))=t.
%
% Inverse rotation:   roteu2ro(m,t)*roteu2ro(fliplr(m),-fliplr(t))=eye(3)

%
%      Copyright (C) Mike Brookes 2007-2012
%      Version: $Id: roteu2ro.m 2171 2012-07-12 07:33:03Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=rotqr2ro(roteu2qr(m,t));


text:
%function Main_CNN_ImageNet_minimal()
% Minimalistic demonstration of how to run an ImageNet CNN model

% setup toolbox
addpath(genpath('../CoreModules'))
% download a pre-trained CNN from the web
if ~exist('imagenet-vgg-f.mat', 'file')
  fprintf('Downloading a model ... this may take a while\n') ;
  urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-f.mat', ...
    'imagenet-vgg-f.mat') ;
end
net = load('imagenet-vgg-f.mat') ;

% obtain and preprocess an image
im = imread('test_im.JPG') ;
im_ = single(im) ; % note: 255 range
im_ = imresize(im_, net.meta.normalization.imageSize(1:2)) ;
im_ = im_ - net.meta.normalization.averageImage ;

% run the CNN
opts=[];
opts.use_gpu=0;%unless you have a good gpu
opts.use_nntoolbox=0; %Requires Neural Network Toolbox to use it.

opts.training=0;
opts.use_corr=1;
res(1).x=im_;

if opts.use_gpu
    net=SwitchProcessor(net,'gpu');
end
tic;
[ net,res,opts ] = net_ff( net,res,opts );
toc;
% show the classification result
scores = squeeze(gather(res(end).x)) ;
[bestScore, best] = max(scores) ;
figure(1) ; clf ; imagesc(im) ;
title(sprintf('%s (%d), score %.3f',...
   net.meta.classes.description{best}, best, bestScore)) ;
drawnow;

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function cj = del_in_source(tdeps, cj)

% function cj = del_in_source(tdeps, cj)
% delete foreign target dependencies from own source dependencies.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: del_in_source.m 1716 2008-05-23 08:18:45Z volkmar $

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

for k = 1:numel(tdeps)
    sitem = subsref(cj, tdeps(k).src_exbranch); % Source item to deal with
    stind = false(1,numel(sitem.sdeps));
    for l = 1:numel(sitem.sdeps)
        stind(l) = ~(isequaltarget(tdeps(k), sitem.sdeps(l)) && ...
                     isequalsource(tdeps(k), sitem.sdeps(l)));
    end;
    sitem.sdeps = sitem.sdeps(stind);
    cj = subsasgn(cj, tdeps(k).src_exbranch,sitem);
end;

author: spm


text:
function [dKSpace, measPara, file] = fLoadkSpace(path)
%% Load kSpace data
%
% (c) Martin Schwartz and Thomas Kuestner
% -------------------------------------------------------------------------

    file='';
    warning('off','MATLAB:dispatcher:nameConflict');
    para = 0;
    % prepare m-files[
    currpath = fileparts(mfilename('fullpath'));

    % load ADCs
    if(ischar(path))
        if(isdir(path))
            [file, path] = uigetfile({'*.mat;*.dat', 'ADC measdata (*.mat, *.dat)'} ,'Select ADC measdata file', 'MultiSelect', 'off', path);
            if(isequal(file,0))
                error('CS_reconstruction(): User Termination');
            end
            [~,filename,ext] = fileparts(file);
        else
            [path, filename, ext] = fileparts(path);
        end
        lReadFromFile = true;
    else
        lReadFromFile = false;
    end

    % read inputs
    if(lReadFromFile)
        if(exist(fullfile(path,[filename,'.mat']),'file'))
            matfile = whos('-file', fullfile(path,[filename,'.mat']));
        end
        if(exist(fullfile(path,[filename,'.mat']),'file') && ~ismember('iLC', {matfile.name}))
            % load mat file to check if it contains iLC/iSP or something else
            data = load(fullfile(path,[filename,'.mat']));
            kSpace = data.kSpace;
            measPara.dim = zeros(1,5);
            measPara.LCall = zeros(1,4);
            dimNames = {'y (phase)', 'x (frequency)', 'z (phase)', 't (time)', 'channels'};
            LCallNames = {'slices', 'acquisitions', 'echos', 'repetitions'};
            dimDefault = ones(1,5);
            LCallDefault = [size(kSpace,1), size(kSpace,4), size(kSpace,5), size(kSpace,3)];
            tmp = size(kSpace{1});
            if(LCallDefault(1) > 1 && length(tmp) > 2)
                dimDefault(1:2) = tmp(1:2);
                dimDefault(4) = tmp(3);
            else
                dimDefault(1:length(tmp)) = tmp;
            end
            dimDefault(5) = size(kSpace,2);
            if(isfield(data,'dim') && isfield(data,'LCall'))
                measPara.dim = data.dim;
                measPara.LCall = data.LCall;
            elseif(isfield(data,'dim') && ~isfield(data,'LCall'))
                measPara.dim = data.dim;
                fprintf('Please specify amount of:\n');
                for i=1:length(measPara.LCall)
                    helper = input(sprintf('%s [%d]: ',LCallNames{i}, LCallDefault(i)),'s');
                    if(isempty(helper)), helper = LCallDefault(i); end;
                    measPara.LCall(i) = helper;
                end

            elseif(~isfield(data,'dim') && isfield(data,'LCall'))
                measPara.LCall = data.LCall;
                fprintf('Please specify dimensions:\n');
                for i=1:length(measPara.dim)
                    helper = input(sprintf('%s [%d]: ',dimNames{i}, dimDefault(i)),'s');
                    if(isempty(helper)), helper = dimDefault(i); end;
                    measPara.dim(i) = helper;
                end

            else
                fprintf('Please specify dimensions:\n');
                for i=1:length(measPara.dim)
                    helper = input(sprintf('%s [%d]: ',dimNames{i}, dimDefault(i)),'s');
                    if(isempty(helper)), helper = dimDefault(i); end;
                    measPara.dim(i) = helper;
                end
                fprintf('Please specify amount of:\n');
                for i=1:length(measPara.LCall)
                    helper = input(sprintf('%s [%d]: ',LCallNames{i}, LCallDefault(i)),'s');
                    if(isempty(helper)), helper = LCallDefault(i); end;
                    measPara.LCall(i) = helper;
                end
            end
            clear 'LCallDefault' 'dimDefault';
            loadedVars = {'dimension', 'oversampling', 'measPara', 'iLC', 'iLCPositions', 'drecksMDH'};
            for i=1:length(loadedVars)
                if(isfield(data,loadedVars{i}))
                    if(any(strcmp(loadedVars{i},{'dimension','oversampling'})))
                        eval(sprintf('measPara.%s = data.%s;', loadedVars{i}, loadedVars{i}));
                    else
                        eval(sprintf('%s = data.%s;', loadedVars{i}, loadedVars{i}));
                    end
                else
                    if(strcmp(loadedVars{i},'measPara'))
                        continue;
                    end
                    eval(sprintf('%s = [];', loadedVars{i}));
                end
            end
            clear 'data';
        else
            [dData, iLC, iEvalInfoMask, drecksMDH, iLCPositions] = fMeas_main(fullfile(path, [filename, ext]));
            measPara = [];
        end
    end

    %% get data dimension and loop counters
    espresso.reconType = 1;
    postproc.aniso = [];
    [measPara, ~, ~] = extractMeasPara( iLC, measPara, espresso, postproc, logical([0 0 0]), drecksMDH, iLCPositions );

    %% extract kSpace data from measurement file
    [kSpace, ~] = extractKSpace_main(dData, iLC, iEvalInfoMask, measPara);

    dKSpace = zeros([size(kSpace{1}),size(kSpace,2)]);
    for iI = 1:size(kSpace,2)
        dKSpace(:,:,:,iI) = kSpace{iI};
    end
    % y-x-z-c -> x-y-z-c
    dKSpace = permute(dKSpace, [2 1 3 4]);

end

author: thomaskuestner


text:
function  distSQ = find3DNeighbourDists2(mesh,scaling)
%
%   distSQ = find3DNeighbourDists2(mesh,scaling);
%
% Return a list of the (squared) distances between each node and its
% neighbours. If there are n connections in the mesh, there should be n
% entries in distSQ See also find2DNeighboutDists - does the same thing
% with the 2d mesh.
%
% ARW 031201 - Now takes a scaling argument (equiv of dimdist for
% mrManDist) The returned distances are now in millimeters.
%
% AUTHOR: WADE
% DATE : 062700

if notDefined('scaling'), scaling=[1 1 1]; end

spX = mesh.connectionMatrix;
spY = mesh.connectionMatrix;
spZ = mesh.connectionMatrix;
nVerts = length(mesh.connectionMatrix);

% Comments anyone?
xI = sparse((1:nVerts),(1:nVerts),mesh.vertices(1,:),nVerts,nVerts)*scaling(1);
yI = sparse((1:nVerts),(1:nVerts),mesh.vertices(2,:),nVerts,nVerts)*scaling(2);
zI = sparse((1:nVerts),(1:nVerts),mesh.vertices(3,:),nVerts,nVerts)*scaling(3);

spX = spX*xI;
spY = spY*yI;
spZ = spZ*zI;

% Cols of sp(X,Y,Z) are now ordinates of corresponding vertices
spX = spX - spX';
spY = spY - spY';
spZ = spZ - spZ';
distSQ = (spX.^2 + spY.^2 + spZ.^2);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft