text:
function [nev] = read_neuralynx_nev(filename, varargin);

% READ_NEURALYNX_NEV reads the event information from the *.nev file in a
% Neuralynx dataset directory
%
% Use as
%   nev = read_neuralynx_hdr(datadir, ...)
%   nev = read_neuralynx_hdr(eventfile, ...)
%
% Optional input arguments should be specified in key-value pairs and may include
%   implementation  should be 1, 2 or 3 (default = 3)
%   value           number or list of numbers
%   mintimestamp    number
%   maxtimestamp    number
%   minnumber       number
%   maxnumber       number
%
% The output structure contains all events and timestamps.

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

% get the optional input arguments
flt_value        = keyval('value', varargin);
flt_mintimestamp = keyval('mintimestamp', varargin);
flt_maxtimestamp = keyval('maxtimestamp', varargin);
flt_minnumber    = keyval('minnumber', varargin);
flt_maxnumber    = keyval('maxnumber', varargin);
implementation   = keyval('implementation', varargin); if isempty(implementation), implementation = 3; end

if ft_filetype(filename, 'neuralynx_ds')
  % replace the directory name by the filename
  if     exist(fullfile(filename, 'Events.Nev'))
    filename = fullfile(filename, 'Events.Nev');
  elseif exist(fullfile(filename, 'Events.nev'))
    filename = fullfile(filename, 'Events.nev');
  elseif exist(fullfile(filename, 'events.Nev'))
    filename = fullfile(filename, 'events.Nev');
  elseif exist(fullfile(filename, 'events.nev'))
    filename = fullfile(filename, 'events.nev');
  end
end

% The file starts with a 16*1024 bytes header in ascii, followed by a
% number of records (c.f. trials).
%
% The format of an event record is
%   int16 PktId
%   int16 PktDataSize
%   int64 TimeStamp
%   int16 EventId
%   int16 TTLValue
%   int16 CRC
%   int32 Dummy
%   int32 Extra[0]
%   ...
%   int32 Extra[7]
%   char EventString[0]
%   ...
%   char EventString[127]
% PktId is usually 0x1002.
% PktDataSize is random data.
% Dummy is random data.
% CRC may contain random data.
% Extra is user-defined data.
% TTLValue is the value sent to the computer on a parallel input port.

fid = fopen(filename, 'rb', 'ieee-le');

headersize = 16384;
offset     = headersize;
if ~isempty(flt_minnumber)
  offset = offset + (flt_minnumber-1)*184;
end

nev = [];

if implementation==1
  if ~isempty(flt_maxnumber)
    warning('filtering on maximum number not yet implemneted');
  end
  % this is the slow way of reading it
  % it also does not allow filtering
  fseek(fid, offset, 'bof');
  while ~feof(fid)
    nev(end+1).PktStart     = fread(fid, 1, 'int16');
    nev(end  ).PktId        = fread(fid, 1, 'int16');
    nev(end  ).PktDataSize  = fread(fid, 1, 'int16');
    nev(end  ).TimeStamp    = fread(fid, 1, 'uint64');
    nev(end  ).EventId      = fread(fid, 1, 'int16');
    nev(end  ).TTLValue     = fread(fid, 1, 'uint16');
    nev(end  ).CRC          = fread(fid, 1, 'int16');
    nev(end  ).Dummy        = fread(fid, 1, 'int32');
    nev(end  ).Extra        = fread(fid, 8, 'int32');
    nev(end  ).EventString  = fread(fid, 128, 'char');
  end
end

if implementation==2
  if ~isempty(flt_maxnumber)
    warning('filtering on maximum number not yet implemneted');
  end
  % this is a faster way of reading it and it is still using the automatic type conversion from Matlab
  fp = offset;
  fseek(fid, fp+ 0, 'bof'); PktStart       = fread(fid, inf, 'uint16', 184-2);
  num = length(PktStart);
  fseek(fid, fp+ 2, 'bof'); PktId          = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+ 4, 'bof'); PktDataSize    = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+ 6, 'bof'); TimeStamp      = fread(fid, num, 'uint64=>uint64', 184-8);
  fseek(fid, fp+14, 'bof'); EventId        = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+16, 'bof'); TTLValue       = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+18, 'bof'); CRC            = fread(fid, num, 'uint16', 184-2);
  fseek(fid, fp+22, 'bof'); Dummy          = fread(fid, num, 'int32' , 184-4);
  % read each of the individual extra int32 values and concatenate them
  fseek(fid, fp+24+0*4, 'bof'); Extra1     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+1*4, 'bof'); Extra2     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+2*4, 'bof'); Extra3     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+3*4, 'bof'); Extra4     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+4*4, 'bof'); Extra5     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+5*4, 'bof'); Extra6     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+6*4, 'bof'); Extra7     = fread(fid, num, 'int32', 184-4);
  fseek(fid, fp+24+7*4, 'bof'); Extra8     = fread(fid, num, 'int32', 184-4);
  Extra = [Extra1 Extra2 Extra3 Extra4 Extra5 Extra6 Extra7 Extra8];
  % read the complete data excluding header as char and cut out the piece with the EventString content
  fseek(fid, fp, 'bof'); EventString        = fread(fid, [184 num], 'char');
  EventString = char(EventString(57:184,:)');
end

if implementation==3
  % this is an even faster way of reading it
  if isempty(flt_minnumber)
    flt_minnumber = 1;
  end
  if isempty(flt_maxnumber)
    flt_maxnumber = inf;
  end
  if fseek(fid, offset, 'bof')~=0
    error(ferror(fid));
  end
  buf = fread(fid, (flt_maxnumber-flt_minnumber+1)*184, 'uint8=>uint8');
  [PktStart, PktId , PktDataSize , TimeStamp , EventId , TTLValue , CRC , Dummy , Extra1 , Extra2 , Extra3 , Extra4 , Extra5 , Extra6 , Extra7 , Extra8 , EventString] = ...
    cstructdecode(buf, 'int16', 'int16', 'int16', 'uint64', 'int16', 'uint16', 'int16', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'int32', 'char128');
  % the cstructdecode function does not respect the byte order for other-endian numbers
  if bigendian
    PktStart     = swapbytes(PktStart);
    PktId        = swapbytes(PktId);
    PktDataSize  = swapbytes(PktDataSize);
    TimeStamp    = swapbytes(TimeStamp);
    EventId      = swapbytes(EventId);
    TTLValue     = swapbytes(TTLValue);
  % CRC          = swapbytes(CRC);
  % Dummy        = swapbytes(Dummy);
  % Extra1       = swapbytes(Extra1);
  % Extra2       = swapbytes(Extra2);
  % Extra3       = swapbytes(Extra3);
  % Extra4       = swapbytes(Extra4);
  % Extra5       = swapbytes(Extra5);
  % Extra6       = swapbytes(Extra6);
  % Extra7       = swapbytes(Extra7);
  % Extra8       = swapbytes(Extra8);
  end
end

fclose(fid);

if implementation~=1
  % make a selection of events
  sel = true(size(TTLValue));
  if ~isempty(flt_mintimestamp)
    sel = sel & (TimeStamp>=flt_mintimestamp);
  end
  if ~isempty(flt_maxtimestamp)
    sel = sel & (TimeStamp<=flt_maxtimestamp);
  end
  if ~isempty(flt_value)
    tmp = false(size(sel));
    for i=1:length(flt_value)
      tmp = tmp | (TTLValue==flt_value(i));
    end
    sel = sel & tmp;
  end

  % restructure the data into a struct array, by first making cell-arrays out of it ...
  numsel = sum(sel);
  PktStart      = mat2cell(PktStart(sel)     , ones(1,numsel), 1);
  PktId         = mat2cell(PktId(sel)        , ones(1,numsel), 1);
  PktDataSize   = mat2cell(PktDataSize(sel)  , ones(1,numsel), 1);
  TimeStamp     = mat2cell(TimeStamp(sel)    , ones(1,numsel), 1);
  EventId       = mat2cell(EventId(sel)      , ones(1,numsel), 1);
  TTLValue      = mat2cell(TTLValue(sel)     , ones(1,numsel), 1);
  % CRC           = mat2cell(CRC(sel)          , ones(1,numsel), 1);
  % Dummy         = mat2cell(Dummy(sel)        , ones(1,numsel), 1);
  % Extra1        = mat2cell(Extra1(sel)       , ones(1,numsel), 1);
  % Extra2        = mat2cell(Extra2(sel)       , ones(1,numsel), 1);
  % Extra3        = mat2cell(Extra3(sel)       , ones(1,numsel), 1);
  % Extra4        = mat2cell(Extra4(sel)       , ones(1,numsel), 1);
  % Extra5        = mat2cell(Extra5(sel)       , ones(1,numsel), 1);
  % Extra6        = mat2cell(Extra6(sel)       , ones(1,numsel), 1);
  % Extra7        = mat2cell(Extra7(sel)       , ones(1,numsel), 1);
  % Extra8        = mat2cell(Extra8(sel)       , ones(1,numsel), 1);
  EventString   = mat2cell(EventString(sel,:), ones(1,numsel), 128);
  EventNumber   = mat2cell(find(sel) + flt_minnumber - 1, ones(1,numsel), 1); % this helps in the external filtering for BCI
  % ... and then convert the cell-arrays into a single structure
  nev = struct(...
    'PktStart'    , PktStart     , ...
    'PktId'       , PktId        , ...
    'PktDataSize' , PktDataSize  , ...
    'TimeStamp'   , TimeStamp    , ...
    'EventId'     , EventId      , ...
    'TTLValue'    , TTLValue     , ...
    'EventString' , EventString  , ...
    'EventNumber' , EventNumber);
end

% remove null values and convert to strings
for i=1:length(nev)
  nev(i).EventString = nev(i).EventString(find(nev(i).EventString));
  nev(i).EventString = char(nev(i).EventString(:)');
end

author: PatternRecognition

repo: OpenBMI


text:
function theta = TLDA_initialize_SAE(numK, numM, numC, TrainData, TestData)
   %% Initialize parameters randomly based on layer sizes.

    train_x = [TrainData TestData]';
    %%  ex1 train a 10 hidden unit SDAE and use it to initialize a FFNN
    %  Setup and train a stacked denoising autoencoder (SDAE)
    rand('state',0)
    sae = saesetup([numM numK]);
    sae.ae{1}.activation_function       = 'sigm';
    sae.ae{1}.learningRate              = 100;
    sae.ae{1}.inputZeroMaskedFraction   = 0;
    opts.numepochs =   1;
    opts.batchsize = size(train_x,1);
    sae = saetrain(sae, train_x, opts);
    % visualize(sae.ae{1}.W{1}(:,2:end)')

    W1 = sae.ae{1}.W{1}(:,2:end);
    b1 = sae.ae{1}.W{1}(:,1:1);
    W11 = sae.ae{1}.W{2}(:,2:end);
    b11 = sae.ae{1}.W{2}(:,1:1);

    data=[TrainData TestData];
    hiddeninputs = W1 * data + b1 * ones(1, size(data,2)); % hiddensize * numpatches
    hiddenvalues = sigmoid( hiddeninputs ); % hiddensize * numpatches

    train_x = hiddenvalues';
    %%  ex1 train a 10 hidden unit SDAE and use it to initialize a FFNN
    %  Setup and train a stacked denoising autoencoder (SDAE)
    rand('state',0)
    sae = saesetup([numK numC]);
    sae.ae{1}.activation_function       = 'sigm';
    sae.ae{1}.learningRate              = 100;
    sae.ae{1}.inputZeroMaskedFraction   = 0;
    opts.numepochs =   1;
    opts.batchsize = size(train_x,1);
    sae = saetrain(sae, train_x, opts);
    % visualize(sae.ae{1}.W{1}(:,2:end)')

    W2 = sae.ae{1}.W{1}(:,2:end);
    b2 = sae.ae{1}.W{1}(:,1:1);
    W22 = sae.ae{1}.W{2}(:,2:end);
    b22 = sae.ae{1}.W{2}(:,1:1);

	% Convert weights and bias gradients to the vector form.
	theta = [W1(:) ; W2(:) ; W22(:) ; W11(:) ; b1(:) ; b2(:) ; b22(:) ; b11(:)];
end

function sigm = sigmoid(x)

    sigm = 1 ./ (1 + exp(-x));
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function r = taylorinit(v,k)
%TAYLORINIT   Initialization of taylor variable
%
%  x = taylorinit(v,k)
%
%The dependent variable x is identified and initialized with v, and Taylor
%coefficients up to k>0 are calculated. For example,
%
%  u = taylorinit(-3,5)
%
%initializes "u" to be a Taylor variable with value -3 and Taylor coefficients
%up to order k=5, i.e. [-3,1,0,0,0,0]. The default for k is 4. After initialization,
%the order of Taylor coefficients can be obtained by
%
%  k = taylororder
%
%In our example, the result will be k=5. The evaluation
%
%  f = inline('x*sin(x)-exp(x^2)')
%  y = f(u)
%
%computes the Taylor coefficients of f up to order 5. Note the first component of y
%is f(u), followed by f'(u), f''(u)/2!, f'''(u)/3!, etc. For example,
%
%  t0 = y{0}
%
%is the 0-th Taylor coefficients, namely f(u), and
%
%  sum( y{0:k} .* e.^(0:k) )
%
%is approximately equal to f(u+e) up to O(e^(k+1)). The subsequent statement
%
%  v = taylor( intval('3.14159_') )
%
%initializes v to be the constant with interval value infsup(3.14158,3.14160).
%
%Taylor expansions are implemented for univariate functions. However, it is
%sometimes useful to evaluate a function at an array of points. The statement
%
%  w = taylorinit((-3:1)',5)
%
%initializes w to be a vector of dependent variables with values -3,-2,-1,0,1. All
%operations with u are executed pointwise, so y=f(w) is executed as w.*sin(w)-exp(w.^2).
%In this case
%
%  y(3)
%
%is the Taylor expansion of f at -1 up to order 5.
%
%For other examples, see demotaylor.
%

% written  05/21/09     S.M. Rump
% modified 09/15/10     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  setappdata(0,'INTLAB_TAYLOR_END',0);

  if nargin==0
    setappdata(0,'INTLAB_TAYLOR_ORDER',0);
    return
  end

  if ~isa(v,'double') & ~isa(v,'intval')
    error('invalid initialization of dependent Taylor variables')
  end
  if nargin==1
    setappdata(0,'INTLAB_TAYLOR_ORDER',4);
  else
    setappdata(0,'INTLAB_TAYLOR_ORDER',k);
  end
  dummy.init = v;
  r = taylor( dummy , 'taylorinit' );

repo: OpenMAS


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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [tSeries] = tSeriesRemoveOutliers(tSeries, baseImage)
%
%

if ~exist('baseImage')
    baseImage = tSeries(1,:,:);
end

nFrames = size(tSeries,1);
nVoxels = size(tSeries,2);
nSlices = size(tSeries,3);

indexes = [];
backup = tSeries;

for curSlice = 1:nSlices

	LSE_ref = tSeries(:,:,curSlice) - repmat(baseImage(:,:,curSlice),[nFrames 1]);
	MSE_ref = sqrt(sum(LSE_ref.^2,2)/(nVoxels));

	mn = mean(MSE_ref);
	sd = std(MSE_ref);

	frame = 1;

	while (frame <= nFrames)

        frameOutlier = frame;
        while (frameOutlier <= nFrames) & ((MSE_ref(frameOutlier) - mn) > 2*sd)
            frameOutlier = frameOutlier + 1;
        end

        if frameOutlier ~= frame
            prev = frame - 1;
            next = frameOutlier;

            if frame == 1
                for i = frame:(frameOutlier - 1)
                    tSeries(i,:) = tSeries(next,:);
                end

            elseif frameOutlier > nFrames
                for i = frame:(frameOutlier - 1)
                    tSeries(i,:) = tSeries(prev,:);
                end

            else
                for i = frame:(frameOutlier - 1)
                    tSeries(i,:) = (tSeries(next,:) - tSeries(prev,:))/(next - prev)*(i - prev) + tSeries(prev,:);
                end

            end

            if frameOutlier - frame > length(indexes)
                indexes = [indexes zeros(1,frameOutlier - frame + length(indexes))];
            end
            indexes(frameOutlier - frame) = indexes(frameOutlier - frame) + 1;

        end

        frame = max(frame + 1,frameOutlier);
	end
end

if sum(indexes) ~= 0

    qstring = 'Here are the number of spike artefacts detected:\n';
	for index = 1:length(indexes)
        if indexes(index) ~= 0
            if indexes(index) == 1
                occurence = 'occurence';
            else
                occurence = 'occurences';
            end

            qstring = [qstring sprintf('\n     - %d consecutive artefacts: %d %s ',index,indexes(index)),occurence];
        end
	end

	qstring = [qstring sprintf('\n\nWould you like to remove them?')];

	button = questdlg(sprintf(qstring),'Remove Outliers','Yes','No','default');

else
    button = 'No';
end

if strcmp(button,'Yes')
   return
else
   tSeries = backup;
   return
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [colormap]=cbrewer(cname, ncol, interp_method)
% CBREWER - This function produces a colorbrewer table (rgb data) for a
% given type, name and number of colors of the colorbrewer tables.
% For more information on 'colorbrewer', please visit
% http://colorbrewer2.org/
%
% [colormap]=cbrewer(cname, ncol, interp_method)
%
% INPUT:
%  cname  name of colortable. One of the following
%      Sequential tables:
%        'Blues'
%        'BuGn'
%        'BuPu'
%        'GnBu'
%        'Greens'
%        'Greys'
%        'Oranges'
%        'OrRd'
%        'PuBu'
%        'PuBuGn'
%        'PuRd'
%        'Purples'
%        'RdPu'
%        'Reds'
%        'YlGn'
%        'YlGnBu'
%        'YlOrBr'
%        'YlOrRd'
%      Divergent tables:
%        'BrBG'
%        'PiYG'
%        'PRGn'
%        'PuOr'
%        'RdBu'
%        'RdGy'
%        'RdYlBu'
%        'RdYlGn'};
%      Qualitative tables:
%        'Accent'
%        'Dark2'
%        'Paired'
%        'Pastel1'
%        'Pastel2'
%        'Set1'
%        'Set2'
%        'Set3'
%   ncol  number of color in the table.
%   interp_method  if the table need to be interpolated, what method
%                  should be used for interp1? Default='pchip' (aka 'cubic')
%
% A note on the number of colors: Based on the original data, there is
% only a certain number of colors available for each type and name of
% colortable. When 'ncol' is larger then the maximum number of colors
% originally given, an interpolation routine is called (interp1) to produce
% the "extended" colormaps.
%
% See also: plot_brewer_cmap
%

%
% The tables were generated from an MS-Excel file provided on the website
% http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer_updates.html
%
% This product includes color specifications and designs developed by
% Cynthia Brewer (http://colorbrewer.org/).
%
%
% Author: Charles Robert
% email: tannoudji@hotmail.com
% Date: 06.12.2011

% load colorbrewer data
load('colorbrewer.mat')
% initialise the colormap is there are any problems
colormap=[];
if (~exist('interp_method', 'var'))
    interp_method='pchip';
end

ctype_names={'div', 'seq', 'qual'};
ctype = [];
for try_ctype = ctype_names
  try_ctype = try_ctype{1};
  if isfield(colorbrewer.(try_ctype),cname)
    ctype = try_ctype;
    break;
  end
end
if isempty(ctype)
  error(['Could not find colormap named: ' cname]);
end

if (ncol>length(colorbrewer.(ctype).(cname)))
  cbrew_init=colorbrewer.(ctype).(cname){length(colorbrewer.(ctype).(cname))};
  colormap=interpolate_cbrewer(cbrew_init, interp_method, ncol);
  colormap=colormap./255;
  return
end

if isempty(colorbrewer.(ctype).(cname){ncol})

  while isempty(colorbrewer.(ctype).(cname){ncol})
    ncol=ncol+1;
  end
  warning( ...
    ['The minimum number of colors for table ' cname ...
     ' is ' num2str(ncol) '.'])
end

colormap=(colorbrewer.(ctype).(cname){ncol})./255;

end
% Copyright (c) 2011, Charles Robert
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%

author: alecjacobson

repo: gptoolbox


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%

real_theta = theta(2:end, :);

regularization = lambda * sum(real_theta .^ 2) / (2 * m);
prediction = sum(((X * theta) .- y) .^ 2) / (2 * m);
J = prediction + regularization;

mask = ones(size(theta));
mask(1) = 0;

gradient0 = sum(((X * theta) .- y) .* X) / m;
offset = lambda * (theta .* mask) / m;
gradient = gradient0 .+ offset';
grad = gradient;

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

grad = grad(:);

end

author: ecmadao

repo: Coding-Guide


text:
% demo for Elogsig

gridM = -10:10;
gridV = 0:10:100;
Nmc = 1e4;
for i=1:length(gridM)
    m = gridM(i);
    for j=1:length(gridV)
        V = gridV(j);
        X = m + sqrt(V).*randn(Nmc,1);
        lsx1 = log(VBA_sigmoid(X));
        lsx2 = log(1-VBA_sigmoid(X));
        Els1(i,j) = mean(lsx1);
        Els2(i,j) = mean(lsx2);
        Els01(i,j) = VBA_Elogsig(m,V);
        Els02(i,j) = VBA_Elogsig(-m,V);
    end
end

figure,plot(VBA_vec(exp(Els1)),VBA_vec(exp(Els01)),'.')
figure,plot(VBA_vec(exp(Els2)),VBA_vec(exp(Els02)),'.')

author: MBB-team

repo: VBA-toolbox


text:
function [exo_irf_record,exo_irf_estimates]=irfexo(beta_gibbs,It,Bu,IRFperiods,IRFband,n,m,p,k,prior)

% create the cell aray that will store the values from the simulations
exo_irf_record=cell(n,m);

% deal with shocks in turn
for ii=1:m

   % step 1: repeat the simulation process a number of times equal to the number of Gibbs iterations
   for kk=1:It-Bu

   % step 3: draw beta from its posterior distribution
   beta=beta_gibbs(:,kk);
   % reshape to obtain B
   B=reshape(beta,k,n);
   % recover C, the matrix of coefficients on exogenous
   C=B(end-m+1:end,:);
   % create a matrix of zeros of dimension p*n
   Y=zeros(p,n);
   %  step 2: set the value of the last row, column i, equal to 1
   Y(p,:)=C(ii,:);

      % step 4: for each iteration kk, repeat the algorithm for periods T+1 to T+h
      for jj=1:IRFperiods-1
% % %           if prior==61
% % %          %use the function lagx to obtain the matrix X; retain only the last row
% % %    X=bear.lagx(Y,p-1);
% % %    X=X(end,:);
% % %
% % %
% % %     else

      % use the function lagx to obtain a matrix temp, containing the endogenous regressors
      temp=bear.lagx(Y,p-1);

      % define the vector X
      X=[temp(end,:) zeros(1,m)];
% % %           end
      % obtain the predicted value for T+jj
      yp=X*B;

      % concatenate yp at the top of Y
      Y=[Y;yp];

      % repeat until values are obtained for T+h
      end

   % step 5: record the results from current iteration in cell irf_record
      % loop over variables
      for jj=1:n
      % consider column jj of matrix 'Y' and trim the (p-1) initial periods: what remains is the series of IRFs for period T to period T+h-1, for variable jj
      temp=Y(p:end,jj);
      % record these values in the corresponding matrix of irf_record
      exo_irf_record{jj,ii}(kk,:)=temp';
      end

   % then go for next iteration
   end

% conduct the same process with shocks in other variables
end

% create first the cell that will contain the IRF estimates
exo_irf_estimates=cell(n,m);

% for the response of each variable to each shock, and each IRF period, compute the median, lower and upper bound from the Gibbs sampler records
% consider variables in turn
for ii=1:n
   % consider shocks in turn
   for jj=1:m
      % consider IRF periods in turn
      for kk=1:IRFperiods
      % compute first the lower bound
      exo_irf_estimates{ii,jj}(1,kk)=quantile(exo_irf_record{ii,jj}(:,kk),(1-IRFband)/2);
      % then compute the median
      exo_irf_estimates{ii,jj}(2,kk)=quantile(exo_irf_record{ii,jj}(:,kk),0.5);
      % finally compute the upper bound
      exo_irf_estimates{ii,jj}(3,kk)=quantile(exo_irf_record{ii,jj}(:,kk),1-(1-IRFband)/2);
      end
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [U,Out] = RecPF(m,n,aTV,aL1,picks,B,PsiT,Psi,opts,varargin)
% [U,Out] = RecPF(m,n,aTV,aL1,picks,B,PsiT,Psi,opts,varargin)
%
% RecPF solves the TVL1-L2 model:
%
%   min aTV*TV(u) + aL1*|PsiT*U|_1 + 0.5|Fp*U - B|_2^2
%
% Inputs:
%
%  m, n     -- size of image
%  aTV, aL1 -- regularization parameters in the model
%  picks    -- sample positions in Fourier domain
%  B        -- measurment vector
%  PsiT     -- sparsifying basis, PsiT*U is the sparse representation of U
%  Psi      -- inverse of PsiT, Psi*x = U reconstructs the image
%  opts      --- contains parameters for algorithm
%              * opts.mit_inn: maxium inner iteration number for each beta {default 30}
%              * opts.mit_out: maxium outer iteration number {default 10}
%              * opts.tol_inn: inner iteration error tolerance {1.e-3}
%                (this option is effective only when "opts.stc = 2", see below)
%              * opts.tol_rel_inn: tolerance of relative change in an inner iteration {5.e-2}
%              * opts.tol_rel_out: tolerance of relative change in an outer iteration {1.e-2}
%                (the above two options are effective only when "opts.stc = 1", see below)
%              * opts.beta0:   initial penalty parameter {2^5}
%              * opts.beta_max: final penalty parameter  {2^15}
%              * opts.beta_rate: increase rate of beta
%              * opts.idisp: 0 or nonzero, display inner iteration info or not
%              * opts.recordf: 0 or 1, keep (or not) history of function values, TV and fidelity;
%              * opts.U0: starting poing {default "a least squares solution"}
%              * opts.wbarflag: 1 (on) or 0 (off), controls wait bar.
%              * opts.stc: 1 (stopping criterion based on Relative Change)
%                          or 2 (stopping criterion based on optimality conditions) {default 1}
%              * opts.TVtype: 1 or 2, corresponding to anisotropic/isotropic TV {default 2}
%  varargin{1}  -- a m x n matrix contains local weights of TV. Specifically, the local weigthed
%                  TV is discretized as:
%                     TV(u) = sum_i w_i ||D_i u||.
%                  If all w_i == 1, then it is just the isotropic discritization of normal TV.
%                  We require all w_i > 0.
%  varargin{2}  -- true image (used to compute relative errors when it is present)
%
% Outputs:
%     U   --- reconsctructed image
%     Out --- a structrue contains
%             * Out.iter: total iteration number
%             * Out.Inner: inner iteration numbers for each beta
%             * Out.ftrue: function values at each iteration
%             * Out.Fvalvscpu: function values decrease with respect to CPU time
%             * Out.TVhist: total variation changes with iteration
%             * Out.FIDhist: fidelity history
%              {the above four subfields are present only when opts.recordf = 1}
%             * Out.ITvsERR: relative error vs. iteration number
%             * Out.CPUvsERR: CPU time vs. relative error
%              {the above two subfields are present only when the original image is present}

%
% Yin Zhang, 11-10-2007
% CAAM, Rice University, Copyright (2007)
%
% Junfeng Yang, last modified on Feb. 2, 2009
%

global Ux Uy FU PsiTU
global Wx Wy
global Z PsiZ
global beta
global remain

if aTV == 0 && aL1 == 0; error('No regularization' ); end
if aTV <  0 || aL1 <  0; error('Regularization < 0'); end
if ~exist('opts','var'); opts = []; end

[mit_inn,mit_out,tol_inn,tol_rel_inn,tol_rel_out,beta0,beta_max, ...
    beta_rate,TVtype,idisp,recordf,U0,wbarflag,stc] = setopts(opts);

[Nomin1,Denom1,Denom2] = getC(picks,B,aTV,m,n);

remain = setdiff(1:m*n,picks);
if isempty(U0); U = zeros(m,n); else U = U0; clear U0; end

if ~isempty(varargin)
    Weights = varargin{1};
    [mm,nn] = size(Weights);
    if mm~=m || nn~= n
        Weights = ones(m,n);
    elseif find(Weights < 0,1)
        error('negative weights are not allowed!');
    end
else
    Weights = ones(m,n);
end
if length(varargin) == 2
    I = varargin{2};
    nrmI = norm(I(:));  rer = norm(U(:)-I(:))/nrmI;
    RER = zeros(500,2); RER(1,:) = [0,rer];
    CPUvsRER = zeros(500,2); CPUvsRER(1,:) = [0,rer];
end

if aL1 > 0; PsiTU = PsiT(U); end

if aTV > 0
    Ux = [diff(U,1,2), U(:,1) - U(:,n)];
    Uy = [diff(U,1,1); U(1,:) - U(m,:)];
end

if recordf
    fval_true = funcval(U,aTV,aL1,B,picks);
    FVAL = zeros(500,2);
    FVAL(1,:) = [0,fval_true];
    Fvalvscpu = zeros(500,2);
    Fvalvscpu(1,:) = [0,fval_true];
    TVhist = zeros(500,1);
    FIDhist = zeros(500,1);
end

% initialization
Inn = zeros(100,1);
beta = beta0;
totalIter = 0;
Outiter = 0;
stopc = 0;
t0 = cputime;

if wbarflag == 1
    wbar = waitbar(0, 'RecPF is running, please wait ...');
end

%% Main loop
while ~stopc
    Outiter = Outiter + 1;
    if stc == 1; Uo = U; end

    Denom = Denom1;
    if aTV > 0; Denom = Denom + (aTV*beta)*Denom2; end
    if aL1 > 0; Denom = Denom + aL1*beta; end

    stopci = 0; Inniter = 0;
    while ~stopci

        % ================================
        %  Begin Alternating Minimization
        % ----------------
        %   W-subprolem
        % ----------------
        if aTV > 0;
            switch TVtype
                case 1;   % anisotropic TV
                    Wx = sign(Ux).* max(abs(Ux)-Weights./beta,0);
                    Wy = sign(Uy).* max(abs(Uy)-Weights./beta,0);
                case 2;   % isotropic TV
                    V = sqrt(Ux.^2 + Uy.^2);
                    S = max(V - Weights./beta, 0);
                    S = S ./ max(V,eps);
                    Wx = S.*Ux; Wy = S.*Uy;
                    clear V S
                otherwise; error('TVtype must be 1 or 2');
            end
        end

        % ----------------
        %   Z-subprolem
        % ----------------
        if aL1 > 0;
            Z = sign(PsiTU).*max(abs(PsiTU)-1/beta,0);
            PsiZ = Psi(Z);
        end

        % ----------------
        %   U-subprolem
        % ----------------
        if stc == 1; Ui = U; end
        rhs = 0;
        if aTV > 0
            rhs = [Wx(:,end) - Wx(:, 1), -diff(Wx,1,2)];
            rhs = rhs + [Wy(end,:) - Wy(1, :); -diff(Wy,1,1)];
            rhs = (aTV*beta)*rhs;
        end
        if aL1 > 0
            rhs = rhs + (aL1*beta)*PsiZ;
        end

        Nomin = Nomin1 + fft2(rhs);

        FU = Nomin./Denom;
        U = real(ifft2(FU));
        Inniter = Inniter + 1;
        totalIter = totalIter + 1;
        %
        %  End Alternating Minimization
        % ================================

        % -----------------------------------------
        % check if inner stopping criterion is met
        %
        if stc == 1
            if aTV > 0
                Ux = [diff(U,1,2), U(:,1) - U(:,n)];
                Uy = [diff(U,1,1); U(1,:) - U(m,:)];
            end
            if aL1 > 0; PsiTU = PsiT(U); end
            chg_inn = norm(Ui(:) - U(:))/norm(U(:));
            stopci = ((chg_inn < tol_rel_inn) || (Inniter >= mit_inn));
        elseif stc == 2
            [res_wn,res_wz,wzr,res_Zn,res_Zz,Zzr,res_u] = checkopt(U,aTV,aL1,PsiT,picks,B,0,Weights);
            res = [res_wn,res_wz,res_Zn,res_Zz,res_u];
            stopci = (max(res) < tol_inn || (Inniter >= mit_inn));
        else
            error('please specify a stopping criterion.');
        end

        if recordf && totalIter < 500
            [fval_true,tvu,fidu] = funcval(U,aTV,aL1,B,picks);
            TVhist(totalIter) = tvu;
            FIDhist(totalIter) = fidu;
            FVAL(totalIter+1,:) = [totalIter,fval_true];
            Fvalvscpu(totalIter+1,:) = [cputime - t0, fval_true];
        end
        if exist('I','var') && totalIter < 500
            rer = norm(U(:)-I(:))/nrmI;
            RER(totalIter+1,:) = [totalIter,rer];
            CPUvsRER(totalIter+1,:) = [cputime - t0, rer];
        end
        if idisp && stc == 2
            fprintf('Iter: %d, res_ wn: %4.1e, wz: %4.1e, Zn: %4.1e, Zz: %4.1e, u: %4.1e, wzr %2.0f%%, Zzr %2.0f%%\n', ...
                totalIter, res_wn,res_wz,res_Zn,res_Zz,res_u,100*wzr,100*Zzr);
        elseif idisp && stc == 1
            fprintf('Iter: %d, chg_inn %4.2f\n',totalIter,chg_inn);
        end

    end % inner
    Inn(Outiter) = Inniter;

    % ------------------------------------------
    % check if outer stopping criterion is met
    %
    beta = beta_rate*beta;
    if stc == 1
        chg_out = norm(U(:)-Uo(:))/norm(U(:));
        stopc = ((chg_out < tol_rel_out) || (Outiter >= mit_out));
    elseif stc == 2
        stopc = ((beta > beta_max) || (Outiter >= mit_out));
    else
        error('please specify a stopping criterion.');
    end
    if wbarflag == 1
        waitbar(log2(beta)/(log2(beta_max)+1), wbar)
    end

end % outer
Out.iter = totalIter;
Out.Inner = Inn(1:Outiter);
if wbarflag == 1
    close(wbar);
end
if recordf
    Out.ftrue = FVAL(1:totalIter+1,:);
    Out.Fvalvscpu = Fvalvscpu(1:totalIter+1,:);
    Out.TVhist = TVhist(1:totalIter);
    Out.FIDhist = FIDhist(1:totalIter);
end
if exist('I','var')
    Out.ITvsRER  = RER(1:totalIter+1,:);
    Out.CPUvsRER = CPUvsRER(1:totalIter+1,:);
end

%% ----------------- SUBFUNCTION ---------------------------
function [res_wn,res_wz,wzr,res_Zn,res_Zz,Zzr,res_u] = ...
    checkopt(U,aTV,aL1,PsiT,picks,f,flag,Weights)
%
% This function checks the optimality of
%                   (              beta                 )
% min   aTV * sum_i (||(W_i)||_2 + ---- ||D_iu - W_i||^2)
%                   (               2                   )
%
%                   (         beta             )
%     + aL1 *       ( |Z|_1 + ---- |Z-PsiT*u|^2 )
%                   (           2              )
%
%     + .5 |F_p*u - f|^2
%
% where W_i = (Wx_i;Wy_i);
%
% U     --- current point
% f     --- measurment vector
% flag  --- 0 (do not check res_u), otherwise check res_u;
%           Generally, res_u << 0, so setting flag = 0 is fine.
% Weights -- local weights of TV
%

% Junfeng Yang, Aug. 08, 2008

global Ux Uy FU PsiTU
global Wx Wy
global beta
global remain
global Z PsiZ

[m,n] = size(U);
mn = m*n;

% res_wn, res_wz, wzr
if aTV > 0
    Ux = [diff(U,1,2), U(:,1) - U(:,n)];
    Uy = [diff(U,1,1); U(1,:) - U(m,:)];

    Dx = Wx - Ux;
    Dy = Wy - Uy;

    V = sqrt(Wx.^2 + Wy.^2);
    Iz = find(V == 0);

    if isempty(Iz)
        In = 1:mn;
    else
        In = find(V ~= 0);
        V(Iz) = 1;
    end

    V = V*beta;

    RWx = Dx + Weights.*Wx./V;
    RWy = Dy + Weights.*Wy./V;
    S = sqrt(RWx.^2 + RWy.^2);

    wzr = length(Iz)/mn;
    if isempty(In)
        res_wn = 0;
    else
        res_wn = max(S(In));
    end

    if isempty(Iz);
        res_wz = -1;
    else
        res_wz = max(S(Iz) - Weights(Iz)/beta);
    end

else
    res_wn = 0;
    res_wz = 0;
    wzr = 0;
end

% res_Zn, res_Zz, Zzr
if aL1 > 0
    PsiTU = PsiT(U);
    absPsiTU = abs(PsiTU);

    ZIn = find(Z ~= 0);
    if isempty(ZIn)
        res_Zn = 0;
    else
        R3 = (1/beta)*sign(Z(ZIn)) + Z(ZIn) - PsiTU(ZIn);
        res_Zn = max(abs(R3));
    end

    ZIz = find(Z == 0);
    if isempty(ZIz)
        res_Zz = -1;
    else
        R4 = absPsiTU(ZIz) - 1/beta;
        res_Zz = max(R4);
    end

    Zzr = length(ZIz)/m/n;

else
    PsiTU = 0;
    res_Zn = 0;
    res_Zz = -1;
    Zzr = 0;
end

% res_u
res_u = 0;
if flag == 0
    return;
else
    % res_u
    FU(picks) = -f + FU(picks);
    FU(remain) = 0;
    Ru = ifft2(FU);
    Ru = real(Ru);

    if aTV > 0
        Dxx = [Dx(:,1) - Dx(:,n), diff(Dx,1,2)];
        Dyy = [Dy(1,:) - Dy(m,:); diff(Dy,1,1)];
        Ru = Ru + (beta*aTV)*(Dxx + Dyy);
    end

    if aL1 > 0
        Ru = Ru + (aL1*beta)*(U - PsiZ);
    end
    res_u = max(abs(Ru(:)));
end
%% ------------- SUBFUNCTION ---------------
function [Nomin1,Denom1,Denom2] = getC(picks,B,aTV,m,n)

% compute fixed quantities
Nomin1 = zeros(m,n);
Nomin1(picks) = B;
Denom1 = zeros(m,n);
Denom1(picks) = 1;
Denom2 = 0;
if aTV > 0
    Denom2 = abs(psf2otf([1,-1],[m,n])).^2 + abs(psf2otf([1;-1],[m,n])).^2;
end
%% ----------  SUBFUNCTION ---------------
function [mit_inn,mit_out,tol_inn,tol_rel_inn,tol_rel_out,beta0,beta_max, ...
    beta_rate,TVtype,idisp,recordf,U0,wbarflag,stc] = setopts(opts)

% define default option fields
mit_inn = 30;
mit_out = 10;
tol_inn  = 1.e-3;
tol_rel_inn = 5.e-3;
tol_rel_out = 1.e-2;
beta0 = 2^5;
beta_max = 2^15;
beta_rate = 2;
TVtype = 2;
idisp = 0;
recordf = 0;
U0 = [];
wbarflag = 0;
stc = 1;

% change to specified option fields if exist
if ~isempty(opts);
    if ~isa(opts,'struct'); error('L1pfi: opts not a struct'); end
    if isfield(opts,'mit_inn'); mit_inn = opts.mit_inn; end
    if isfield(opts,'mit_out'); mit_out = opts.mit_out; end
    if isfield(opts,'tol_inn'); tol_inn = opts.tol_inn; end
    if isfield(opts,'tol_rel_inn'); tol_rel_inn = opts.tol_rel_inn; end
    if isfield(opts,'tol_rel_out'); tol_rel_out = opts.tol_rel_out; end
    if isfield(opts,'beta0');   beta0 = opts.beta0; end
    if isfield(opts,'beta_max'); beta_max = opts.beta_max; end
    if isfield(opts,'beta_rate'); beta_rate = opts.beta_rate; end
    if isfield(opts,'TVtype'); TVtype = opts.TVtype; end
    if isfield(opts,'idisp'); idisp = opts.idisp; end
    if isfield(opts,'recordf'); recordf = opts.recordf; end
    if isfield(opts,'U0'); U0 = opts.U0; end
    if isfield(opts,'wbarflag'); wbarflag = opts.wbarflag; end
    if isfield(opts,'stc'); stc = opts.stc; end
end

repo: REST


text:
function geometry_test004 ( )

%*****************************************************************************80
%
%% TEST004 tests R8_ACOS;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  test_num = 9;

  test_x = [ 5.0, 1.2, 1.0, 0.9, 0.5, 0.0, -0.9, -1.0, -1.01 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST004\n' );
  fprintf ( 1, '  R8_ACOS computes an angle with a given cosine;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X, R8_ACOS(X), (Degrees)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : test_num

    x = test_x(i);

    temp1 = r8_acos ( x );
    temp2 = radians_to_degrees ( temp1 );

    fprintf ( 1, '  %8f  %12f  %12f\n', x, temp1, temp2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = uminus(varargin)
%UMINUS   Unary minus for a SPHEREFUN.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear all; close all;
I=imread('circuit.tif');
I=im2double(I);
BW=edge(I, 'canny');
[H, Theta, Rho]=hough(BW, 'RhoResolution', 0.5, 'ThetaResolution', 0.5);
figure;
set(0,'defaultFigurePosition',[100,100,1000,500]);
set(0,'defaultFigureColor',[1 1 1])
subplot(121);
imshow(BW);
subplot(122);
imshow(imadjust(mat2gray(H)));
axis normal;
hold on;
colormap hot;

author: HuangCongQing


text:
function figure_num = circle_segment_test06 ( figure_num )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_TEST06 samples using CIRCLE_SEGMENT_SAMPLE_FROM_HEIGHT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2013
%
%  Author:
%
%    John Burkardt
%
  if ( nargin < 1 )
    figure_num = 0;
  end

  n = 100;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIRCLE_SEGMENT_TEST06\n' );
  fprintf ( 1, '  CIRCLE_SEGMENT_SAMPLE_FROM_HEIGHT samples a circle segment.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot %d points from several segments.\n', n );
  fprintf ( 1, '\n' );

  r = 1.0;
  theta = pi;

  for i = 0 : 3

    h = circle_segment_height_from_angle ( r, theta );

    thetah = theta / 2.0;

    an = 0.5 * pi + linspace ( -thetah, +thetah, 51 );
    x = r * cos ( an );
    y = r * sin ( an );
    x = [ x, x(1) ];
    y = [ y, y(1) ];

    [ xs, ys, seed ] = circle_segment_sample_from_height ( r, h, n, seed );

    figure_num = figure_num + 1;
    figure ( figure_num );
    clf
    hold on
    plot ( x, y, 'k-' );
    plot ( xs, ys, 'b.', 'Markersize', 20 );
    grid on
    axis equal
    xlabel ( '<--- X --->' );
    ylabel ( '<--- Y --->' );
    title ( sprintf ( 'Segment with THETA = %g, R = %g H = %g\n', theta, r, h ) );
    hold off

    filename = sprintf ( 'sample_t%g.png', theta );
    print ( '-dpng', filename );
    fprintf ( 1, '  Created graphics file "%s".\n', filename );

    theta = theta / 2.0;

  end

  return
end


text:
function M = getmassmat3(elem2dof,volume,Dlambda,elemType,K)
%% GETMASSMAT Get mass matrix of the finite element space
% M = getmassmat(Dlambda,elem2dof,elemType,K) get mass matrix of the finite element
% space specified by elemType.
%
% The elemType can be:
%
% - "RT0": The lowest order Raviart-Thomas element
% - "ND1': The lowest order Nedelec element
%

if ~exist('K','var'), K = []; end
NT = size(elem2dof,1);

%% ND1: the lowest order edge element
if strcmp(elemType,'ND1')
    NE = double(max(elem2dof(:)));
    DiDj = zeros(NT,4,4);
    for i = 1:4
        for j = i:4
            DiDj(:,i,j) = dot(Dlambda(:,:,i),Dlambda(:,:,j),2);
            DiDj(:,j,i) = DiDj(:,i,j);
        end
    end
    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
    ii = zeros(21*NT,1); jj = zeros(21*NT,1); sM = zeros(21*NT,1);
    index = 0;
    for i = 1:6
        for j = i:6
            ii(index+1:index+NT) = double(elem2dof(:,i));
            jj(index+1:index+NT) = double(elem2dof(:,j));
            % mass matrix
            % locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
            i1 = locEdge(i,1); i2 = locEdge(i,2);
            j1 = locEdge(j,1); j2 = locEdge(j,2);
            Mij = 1/20*volume.*( (1+(i1==j1))*DiDj(:,i2,j2) ...
                               - (1+(i1==j2))*DiDj(:,i2,j1) ...
                               - (1+(i2==j1))*DiDj(:,i1,j2) ...
                               + (1+(i2==j2))*DiDj(:,i1,j1));
            sM(index+1:index+NT) = Mij;
            index = index + NT;
        end
    end
    diagIdx = (ii == jj);   upperIdx = ~diagIdx;
    M = sparse(ii(diagIdx),jj(diagIdx),sM(diagIdx),NE,NE);
    MU = sparse(ii(upperIdx),jj(upperIdx),sM(upperIdx),NE,NE);
    M = M + MU + MU';
end

%% RT0: the lowest order face element
if strcmp(elemType,'RT0')
    NF = double(max(elem2dof(:)));
    localFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; % normal vector is not oriented
    M = sparse(NF,NF);
    for i = 1:4
        for j = i:4
            % local to global index map
            ii = double(elem2dof(:,i));
            jj = double(elem2dof(:,j));
            i1 = localFace(i,1); i2 = localFace(i,2); i3 = localFace(i,3);
            j1 = localFace(j,1); j2 = localFace(j,2); j3 = localFace(j,3);
            % computation of mass matrix --- (phi_i, phi_j)
            Mij = 1/20*volume*4.*( ...
                  (1+(i1==j1))*dot(mycross(Dlambda(:,:,i2),Dlambda(:,:,i3),2), ...
                                   mycross(Dlambda(:,:,j2),Dlambda(:,:,j3),2),2)...
                 +(1+(i1==j2))*dot(mycross(Dlambda(:,:,i2),Dlambda(:,:,i3),2), ...
                                   mycross(Dlambda(:,:,j3),Dlambda(:,:,j1),2),2)...
                 +(1+(i1==j3))*dot(mycross(Dlambda(:,:,i2),Dlambda(:,:,i3),2), ...
                                   mycross(Dlambda(:,:,j1),Dlambda(:,:,j2),2),2)...
                 +(1+(i2==j1))*dot(mycross(Dlambda(:,:,i3),Dlambda(:,:,i1),2), ...
                                   mycross(Dlambda(:,:,j2),Dlambda(:,:,j3),2),2)...
                 +(1+(i2==j2))*dot(mycross(Dlambda(:,:,i3),Dlambda(:,:,i1),2), ...
                                   mycross(Dlambda(:,:,j3),Dlambda(:,:,j1),2),2)...
                 +(1+(i2==j3))*dot(mycross(Dlambda(:,:,i3),Dlambda(:,:,i1),2), ...
                                   mycross(Dlambda(:,:,j1),Dlambda(:,:,j2),2),2)...
                 +(1+(i3==j1))*dot(mycross(Dlambda(:,:,i1),Dlambda(:,:,i2),2), ...
                                   mycross(Dlambda(:,:,j2),Dlambda(:,:,j3),2),2)...
                 +(1+(i3==j2))*dot(mycross(Dlambda(:,:,i1),Dlambda(:,:,i2),2), ...
                                   mycross(Dlambda(:,:,j3),Dlambda(:,:,j1),2),2)...
                 +(1+(i3==j3))*dot(mycross(Dlambda(:,:,i1),Dlambda(:,:,i2),2), ...
                                   mycross(Dlambda(:,:,j1),Dlambda(:,:,j2),2),2));
            if ~isempty(K)
                Mij = Mij./K;
            end
            if (j==i)
                M = M + sparse(ii,jj,Mij,NF,NF);
            else
                M = M + sparse([ii;jj],[jj;ii],[Mij; Mij],NF,NF);
            end
        end
    end
end

repo: Aero-matlab


text:
function [ctf] = ctf_read_res4(folder,VERBOSE,COEFS);

% ctf_read_res4 - Read a CTF .res4 file
%
% ctf = ctf_read_res4( [folder], [verbose], [coefs])
%
% This function reads the resource information from a CTF .ds folder.  This
% resource information must be read before reading the .meg4 data file.
% All input arguments are optional.
%
% INPUTS
%
%   folder: the .ds directory containing the data to be read.  With
%           no input, it will prompt with a gui folder locator
%           (called by ctf_folder).
%
%   verbose: If verbose = 1, display 'ctf.setup' structure (default)
%            If verbose = 0, do not display 'ctf.setup' structure
%
%   coefs: an option to read the MEG sensor and reference coefficients,
%          which give the weights for calculation of synthetic 2nd or 3rd
%          order gradiometers.
%          If coefs = 1, read the sensor coefficients
%          If coefs = 0, do not read the sensor coefficients (this is the
%          default because it is assumed that data preprocessing with
%          CTF tools has already applied AND saved a synthetic
%          gradiometer transformation for the meg4 data file).  The
%          ctf.sensor.info(x).grad_order_no will indicate the gradient
%          order of the data (this only applies to MEG sensors, other
%          channels have zero values).
%
% OUTPUTS
%
%   ctf.folder - path of the .ds folder
%
%   ctf.res4.file - data file path/name
%   ctf.res4.header - data format header
%
%   ctf.setup - a header structure consisting of date, time, run name, run
%   title, subject, run description, operator, number of channels, number
%   of samples, sample rate, number of trials, duration, pretrigger_samples,
%   sensor filename, head zeroing, and number of filters.
%
%   ctf.sensor.index - a sensor structure consisting of EEG sensors, MEG
%   sensors, reference sensors, and other sensors.
%
%   ctf.sensor.info - a structure with gain and offset information consisting
%   of proper gain, Q gain, io gain, io offset, and index.
%
% The proper gain is the channel-wide constant quotient of a value in raw
% units and a value in user units. (For EEG this is one because there is
% no conversion; i.e., the values are in volts, for MEG there is conversion
% between phi0 and teslas).
% proper gain = raw value / user value
%
% The Q-gain is used to convert the 32-bit integers produced by the
% electronics to the raw values stored in the data file.
% raw value = 32 bit value / Q gain
%
% The IO-gain is essentially a correction factor that is applied to the
% raw values before they are stored to disk.
%
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%      <                                                      >
%      <                    DISCLAIMER:                       >
%      <                                                      >
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. >
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   >
%      <                    OFFICIAL USE.                     >
%      <                                                      >
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:27 $

% Copyright (C) 2003  Darren L. Weber
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% Modified: 11/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - modified from NIH code readresfile.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.1 $';
fprintf('\nCTF_READ_RES4 [v %s]\n',ver(11:15)); tic;

if ~exist('folder','var'),
    ctf = ctf_folder;
else
    ctf = ctf_folder(folder);
end
%[folderPath,folderName,folderExt] = fileparts(ctf.folder);
ctf.res4.file = findres4file(ctf.folder);

if ~exist('VERBOSE','var'), VERBOSE = 1; end
if ~exist('COEFS','var'), COEFS = 0; end

%----------------------------------------------------------------
% open the data file

[fid,message] = fopen(ctf.res4.file,'rb','ieee-be.l64');
if fid < 0, error('cannot open .res4 file'); end

%-------------------------------------------------------------
% READ HEADER

fseek(fid,0,-1);
ctf.res4.header = char(fread(fid,8,'char'))';

% check for the right format
if strmatch('MEG41RS',ctf.res4.header),
    % OK, we can handle this format
else
    msg = sprintf('This function is designed to read MEG41RS format.\nIt may not read "%s" format correctly',ctf.res4.header);
    warning(msg);
end

%-------------------------------------------------------------
% READ SETUP

%---DATE/TIME
fseek(fid, 778,-1);
ctf.setup.time = char(fread(fid,255,'char'))';
fseek(fid,1033,-1);
ctf.setup.date = char(fread(fid,255,'char'))';

%---NUMBER OF SAMPLES
fseek(fid,1288,-1);
ctf.setup.number_samples = (fread(fid,1,'int32')');

%---NUMBER OF CHANNELS
fseek(fid,1292,-1);
ctf.setup.number_channels = (fread(fid,1,'int16')');

%---SAMPLE RATE
fseek(fid,1296,-1);
ctf.setup.sample_rate = fread(fid,1,'double')';
ctf.setup.sample_msec = 1000 / ctf.setup.sample_rate;
ctf.setup.sample_sec  =    1 / ctf.setup.sample_rate;

%---NUMBER OF TRIALS
fseek(fid,1312,-1);
ctf.setup.number_trials = fread(fid,1,'int16')';
fseek(fid, 776,-1);
ctf.setup.number_trials_averaged = fread(fid,1,'int16');

%---DURATION
fseek(fid,1304,-1);
ctf.setup.duration_total = fread(fid,1,'double');
ctf.setup.duration_trial = ctf.setup.duration_total / ctf.setup.number_trials;

%---PRE_TRIG POINTS
fseek(fid,1316,-1);
ctf.setup.pretrigger_samples = fread(fid,1,'int32');
ctf.setup.pretrigger_msec = (ctf.setup.pretrigger_samples / ctf.setup.sample_rate) * 1000;

%---HEAD ZEROING
fseek(fid,1348,-1);
h_zero = fread(fid,1,'int32')';
no_yes = {'no','yes'};
ctf.setup.head_zero = no_yes{h_zero+1};

%---RUN NAME
fseek(fid,1360,-1);
ctf.setup.run_name = char(fread(fid,32,'char'))';

%---RUN TITLE
fseek(fid,1392,-1);
ctf.setup.run_title = char(fread(fid,256,'char'))';

%---SUBJECT
fseek(fid,1712,-1);
ctf.setup.subject = char(fread(fid,32,'char'))';

%---OPERATOR
fseek(fid,1744,-1);
ctf.setup.operator = char(fread(fid,32,'char'))';

%---SENSOR FILE NAME
fseek(fid,1776,-1);
ctf.setup.sensor_file_name = char(fread(fid,60,'char'))';

%---RUN DESCRIPTION & FILTERS
fseek(fid,1836,-1);
run_size = fread(fid,1,'int32');
fseek(fid,1844,-1);
ctf.setup.run_description = char(fread(fid,run_size,'char'));

ctf.setup.number_filters = fread(fid,1,'int16');

for i = 1:ctf.setup.number_filters,
    ctf.setup.filters(i).freq     = fread(fid,1,'double');
    ctf.setup.filters(i).class    = fread(fid,1,'int32');
    ctf.setup.filters(i).type     = fread(fid,1,'int32');
    ctf.setup.filters(i).numparam = fread(fid,1,'int16');
    ctf.setup.filters(i).params   = fread(fid,ctf.setup.filters(i).numparam,'double');
end

if(COEFS)
    b = 1846 + run_size;
    if(ctf.setup.number_filters == 0)
        np = 0;
    else
        np = sum([ctf.setup.filters.numparam]);
        warning('3rd gradient + hardware filter parameters not fully tested! let''s see what happens... :)');
    end
    nf = ctf.setup.number_filters;
    f = ( nf * 18 ) + ( np * 8 );
    offset = b + f + ctf.setup.number_channels * 1360;
end

%-------------------------------------------------------------
% CREATE TIME ARRAYS

% the time arrays must be based on increments of the sample_msec
ctf.setup.time_msec = [0:ctf.setup.number_samples - 1]' * ctf.setup.sample_msec;
ctf.setup.time_msec = ctf.setup.time_msec - ctf.setup.pretrigger_msec;

% adjust the sample point closest to zero so that it is zero, if it
% is reasonably close to zero, say within 3 decimal places for msec timing
zero_index = find(abs(ctf.setup.time_msec) == min(abs(ctf.setup.time_msec)));
zero_value = ctf.setup.time_msec(zero_index);
if (-0.0001 < zero_value) & (zero_value < 0.0001),
    ctf.setup.time_msec(zero_index) = 0;
end

ctf.setup.start_msec = ctf.setup.time_msec(1);
ctf.setup.end_msec   = ctf.setup.time_msec(end);

ctf.setup.time_sec  = ctf.setup.time_msec / 1000;
ctf.setup.start_sec = ctf.setup.time_sec(1);
ctf.setup.end_sec   = ctf.setup.time_sec(end);

%-------------------------------------------------------------
% PRINT SETUP
if VERBOSE,
    ctf_print_setup(ctf);
end

%-------------------------------------------------------------
% READ SENSOR INFORMATION

ctf.sensor.info = struct(...
    'proper_gain',[],...
    'q_gain',[],...
    'io_gain',[],...
    'io_offset',[],...
    'index',[],...
    'extra',[],...
    'label',[],...
    'grad_order_no',[]);

% read channel names
for chan = 1:ctf.setup.number_channels,
    temp = fread(fid,32,'char');
    ctf.sensor.info(chan).label = parse_sensor_label(temp);
end

for chan = 1:ctf.setup.number_channels,

    % The proper gain is the channel-wide constant quotient of a value in raw
    % units and a value in user units. (For EEG this is one because there is
    % no conversion; i.e., the values are in volts, for MEG there is conversion
    % between phi0 and teslas).
    % proper gain = raw value / user value
    %
    % The Q-gain is used to convert the 32-bit integers produced by the
    % electronics to the raw values stored in the data file.
    % raw value = 32 bit value / Q gain
    %
    % The IO-gain is essentially a correction factor that is applied to the
    % raw values before they are stored to disk.

    %ftell(fid);

    ctf.sensor.info(chan).index = fread(fid,1,'int16');
    ctf.sensor.info(chan).extra = fread(fid,1,'int16');
    id = fread(fid,1,'int32')+1;
    ctf.sensor.info(chan).proper_gain = fread(fid,1,'double');
    ctf.sensor.info(chan).q_gain = fread(fid,1,'double');
    ctf.sensor.info(chan).io_gain = fread(fid,1,'double');
    ctf.sensor.info(chan).io_offset = fread(fid,1,'double');
    fread(fid,1,'int16');
    ctf.sensor.info(chan).grad_order_no = fread(fid,1,'int16');
    fread(fid,1,'int32');

    %fseek(fid,ftell(fid)+6,0);

    for pos = 1:8,
        ctf.sensor.info(chan).coil(pos).position.x = fread(fid,1,'double');
        ctf.sensor.info(chan).coil(pos).position.y = fread(fid,1,'double');
        ctf.sensor.info(chan).coil(pos).position.z = fread(fid,1,'double');
        fread(fid,1,'double');
        ctf.sensor.info(chan).coil(pos).orient.x = fread(fid,1,'double');
        ctf.sensor.info(chan).coil(pos).orient.y = fread(fid,1,'double');
        ctf.sensor.info(chan).coil(pos).orient.z = fread(fid,1,'double');
        fread(fid,1,'double');
        fread(fid,1,'int16');
        fread(fid,1,'int32');
        fread(fid,1,'int16');
        fread(fid,1,'double');

        %fseek(fid,ftell(fid)+56,0);
        %fseek(fid,ftell(fid)-80,0);
    end

    for pos = 1:8,
        ctf.sensor.info(chan).hcoil(pos).position.x = fread(fid,1,'double');
        ctf.sensor.info(chan).hcoil(pos).position.y = fread(fid,1,'double');
        ctf.sensor.info(chan).hcoil(pos).position.z = fread(fid,1,'double');
        fread(fid,1,'double');
        ctf.sensor.info(chan).hcoil(pos).orient.x = fread(fid,1,'double');
        ctf.sensor.info(chan).hcoil(pos).orient.y = fread(fid,1,'double');
        ctf.sensor.info(chan).hcoil(pos).orient.z = fread(fid,1,'double');
        fread(fid,1,'double');
        fread(fid,1,'int16');
        fread(fid,1,'int32');
        fread(fid,1,'int16');
        fread(fid,1,'double');

        %fseek(fid,ftell(fid)+56,0);
        %fseek(fid,ftell(fid)+80,0);
    end
    %fseek(fid,ftell(fid)+1288,-1);
end

%-------------------------------------------------------------
% Find channel types and define channel sets, see the
% System Administrators .pdf, 'Channel Sets Configuration'

ctf = ctf_channel_sets(ctf);

%-------------------------------------------------------------
% Channel coordinates, in centimeters, in subject head space

for chan = 1:ctf.setup.number_channels,

    switch ctf.sensor.info(chan).index,

        case {0,1,5},
            %0=Reference Magnetometers
            %1=Reference Gradiometers
            %5=MEG Channels

            coord = [ctf.sensor.info(chan).hcoil(1:2).position];
            ctf.sensor.info(chan).location = [coord.x; coord.y; coord.z];

            orient = [ctf.sensor.info(chan).hcoil(1).orient];
            ctf.sensor.info(chan).orientation = [orient.x; orient.y; orient.z];

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % This ensures that the orientation of the sensor is away from the
            % center of a sphere.  It uses the sign of the dot product between
            % the orientation vector and the location vector.
            tmp = ctf.sensor.info(chan).orientation' * ctf.sensor.info(chan).location;
            tmp = sign(tmp(1));
            ctf.sensor.info(chan).orientation = tmp * ctf.sensor.info(chan).orientation;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        case 9,
            %EEG Channels

            coord = [ctf.sensor.info(chan).hcoil(1:2).position];
            ctf.sensor.info(chan).location = [coord.x; coord.y; coord.z];
            ctf.sensor.info(chan).orientation = [];

    end
end

%-------------------------------------------------------------------------
% assign sensor locations and orientations
% simplify information in the ctf.sensor.info struct

ctf.sensor.location = zeros(3,ctf.setup.number_channels);
ctf.sensor.orientation = zeros(3,ctf.setup.number_channels);

ctf.sensor.label = [];
ctf.sensor.location = [];
ctf.sensor.orientation = [];

for c = 1:ctf.setup.number_channels,

    % All channels have a label
    ctf.sensor.label{1,c} = strtok(ctf.sensor.info(c).label,'-');

    % All channels have a location
    % EEG channels do not have any orientation
    if length(ctf.sensor.type.meg_ref) > 1, ind2 = 2; else ind2 = 1; end;

    switch ctf.sensor.info(c).index,

        case {ctf.sensor.type.meg_sens, ctf.sensor.type.meg_ref(1), ctf.sensor.type.meg_ref(ind2)},
        % modification above based on bug 438 of EEGLAB - Arnaud Delorme, August 2007
            %0=Reference Channels,
            %1=More Reference Channels,
            %5=MEG Channels

            % MEG channels are radial gradiometers, so they have an inner (1) and
            % an outer (2) location - it might be better to take the average of
            % their locations
            if ~isempty(ctf.sensor.info(c).location),
                ctf.sensor.location(:,c) = ctf.sensor.info(c).location(:,1);
            end

            if ~isempty(ctf.sensor.info(c).orientation),
                ctf.sensor.orientation(:,c) = ctf.sensor.info(c).orientation(:,1);
            end

        case ctf.sensor.type.eeg_sens,
            %9=EEG Channels

            if ~isempty(ctf.sensor.info(c).location),
                ctf.sensor.location(:,c) = ctf.sensor.info(c).location(:,1);
            end

    end
end

%%%%%% Coefficient reading appended by SSD %%%%%%

if(COEFS)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % NUMBER OF COEFFICIENTS, byte offset = b+f+nc*1360, byte size = 1

    fseek(fid,offset,-1);
    ctf.res4.numberCoefficients = fread(fid,1,'int16');

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % SENSOR COEFFICIENT RECORDS, byte offset = b+f+nc*1360+2, byte size = 1992

    if VERBOSE & ctf.res4.numberCoefficients,
        fprintf('...reading %d coefficients\n',ctf.res4.numberCoefficients);
    end

    SENSOR_LABEL  = 31;
    MAX_NUM_COEFS = 50;
    MAX_BALANCING = MAX_NUM_COEFS;

    hexadef = {'00000000','47314252','47324252','47334252','47324f49','47334f49'};
    strdef = {'NOGRAD','G1BR','G2BR','G3BR','G2OI','G3OI'};

    for i = 1:ctf.res4.numberCoefficients,

        % read the sensor name (channel label)
        temp = fread(fid,[1,32],'char');
        sensorName = parse_sensor_label(temp);
        sensorIndex = strmatch( sensorName, {ctf.sensor.info.label} );

        % read the coefficient type
        coefType = fread(fid,1,'bit32');

        padding = fread(fid,1,'int32'); % not sure why this is needed???

        % read the coefficient record
        numberCoefs = fread(fid,1,'int16');

        if numberCoefs > MAX_NUM_COEFS,
            msg = sprintf('numberCoefs > MAX_NUM_COEFS\n');
            warning(msg);
        end

        sensor_list = char(fread(fid,[SENSOR_LABEL,MAX_BALANCING],'uchar')');
        % clean-up the sensor_list
        sensor_list = sensor_list(1:numberCoefs,:);
        for j=1:numberCoefs,
            temp = strtok(sensor_list(j,:),char(0));
            temp = strtok(temp,'-');
            % check if this sensor is a reference
            refLabels = {ctf.sensor.info(ctf.sensor.index.meg_ref).label};
            refIndex = strmatch(temp,refLabels);
            if refIndex,
                % ensure this one has the same label
                temp = refLabels{refIndex};
            end

            new_sensor_list(j) = refIndex;
        end
        sensor_list = ctf.sensor.index.meg_ref(:,new_sensor_list)';

        coefs_list = fread(fid,MAX_BALANCING,'double');
        % clean-up the coefs_list
        coefs_list = coefs_list(1:numberCoefs,:)';

        % allocate the coefficient parameters into the ctf struct
        ctf.res4.sensorCoef(i).sensorName = sensorName;
        ctf.res4.sensorCoef(i).coefType   = coefType;
        ctf.res4.sensorCoef(i).coefRec.numberCoefs = numberCoefs;
        ctf.res4.sensorCoef(i).coefRec.sensor_list = sensor_list;
        ctf.res4.sensorCoef(i).coefRec.coefs_list  = coefs_list;

        % DLW:
        % This is a brainstorm variable, note the use of coefType
        % Not clear why this is checked and allocated as such
        coefType = find( hex2dec(hexadef) == coefType );
        if coefType,
            CoefInfo{sensorIndex,coefType-1}.numberCoefs = numberCoefs;
            CoefInfo{sensorIndex,coefType-1}.sensor_list = sensor_list;
            CoefInfo{sensorIndex,coefType-1}.coefs = coefs_list;
        end
    end

    % Channel Gains
    gain_chan = zeros(size(ctf.setup.number_channels,1),1);
    gain_chan(ctf.sensor.index.meg_sens) = ([ctf.sensor.info(ctf.sensor.index.meg_sens).proper_gain]'.*[ctf.sensor.info(ctf.sensor.index.meg_sens).q_gain]');
    gain_chan(ctf.sensor.index.meg_ref) = ([ctf.sensor.info(ctf.sensor.index.meg_ref).proper_gain]'.*[ctf.sensor.info(ctf.sensor.index.meg_ref).q_gain]');
    %  gain_chan(ieegsens) = 1./([SensorRes(ieegsens).qGain]'*1e-6);
    %     gain_chan(ieegsens) = 1./([SensorRes(ieegsens).qGain]');
    %     gain_chan(iothersens) = ([SensorRes(iothersens).qGain]'); % Don't know exactly which gain to apply here

    % Calculus of the matrix for nth-order gradient correction
    % Coefficients for unused reference channels are weigthed by zeros in
    % the correction matrix.
    Gcoef = zeros(length(ctf.sensor.index.meg_sens),length(min(ctf.sensor.index.meg_ref):max(ctf.sensor.index.meg_ref)));
    grad_order_no = 3*ones(306,1);
    for k = 1:length(ctf.sensor.index.meg_sens)

        % Reference coils for channel k
        if grad_order_no(ctf.sensor.index.meg_sens(k)) == 0
            %Data is saved as RAW
            %Save 3rd order gradient sensor-list for subsequent correction if requested later by the user
            [refs] = (CoefInfo{ctf.sensor.index.meg_sens(k),3}.sensor_list);
            Gcoef(k,refs-min(ctf.sensor.index.meg_ref)+1) = CoefInfo{ctf.sensor.index.meg_sens(k),3}.coefs ...
                .* gain_chan(refs)'/gain_chan(ctf.sensor.index.meg_sens(k));
        else
            [refs] = (CoefInfo{ctf.sensor.index.meg_sens(k),grad_order_no(ctf.sensor.index.meg_sens(k))}.sensor_list);
            Gcoef(k,refs-min(ctf.sensor.index.meg_ref)+1) = CoefInfo{ctf.sensor.index.meg_sens(k),grad_order_no(ctf.sensor.index.meg_sens(k))}.coefs ...
                .* gain_chan(refs)/gain_chan(ctf.sensor.index.meg_sens(k));
        end
        ctf.sensor.info(ctf.sensor.index.meg_sens(k)).Gcoef = Gcoef(k,:);
    end
end   %% end COEF block

fclose(fid);

t = toc; fprintf('...done (%6.2f sec)\n\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -------------------------------------------------------
function res4name = findres4file( folder )

% find file name if truncated or with uppercase extension
% added by Arnaud Delorme June 15, 2004

res4name = dir([ folder filesep '*.res4' ]);
if isempty(res4name)
    res4name = dir([ folder filesep '*.RES4' ]);
end

if isempty(res4name)
    error('No file with extension .res4 or .RES4 in selected folder');
else
    res4name = [ folder filesep res4name(1).name ];
end;
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -------------------------------------------------------
function sensorName = parse_sensor_label(temp)

% sensorName = parse_sensor_label(temp)
% parse sensor label names

temp(temp>127) = 0;
temp(temp<0) = 0;
temp = strtok(temp,char(0));
temp = strtok(temp,'-');
sensorName = char(temp)';

return


text:
function binsize_str = binsizelabel(binsize)
    binsize_str = '';
    if binsize == 1.0/1440
        binsize_str = 'per minute';
    elseif binsize == 1.0/24
        binsize_str = 'per hour';
    elseif binsize == 1.0
        binsize_str = 'per day';
    elseif binsize == 7.0
        binsize_str = 'per week';
    elseif binsize >= 28 & binsize <=31
        binsize_str = 'per month';
    elseif binsize >= 365 & binsize <= 366
        binsize_str = 'per year';
    end

    %disp(sprintf('binsize label',binsize_str));
end

author: geoscience-community-codes

repo: GISMO


text:
% bsb_test.m
x=[-0.5;-0.4];
beta=0.5;
c=100;bsb(x,beta,c)
text(-0.5,-0.4,'(-0.5,-0.4)')


text:
function scaleDoseROE(hObj,hEvt,hFig)
% Compute TCP/NTCP at scaled fraction size/number.
%
% AI 05/12/21

ud = guidata(hFig);
protS = ud.Protocols;

%Get selected scale
userScale = get(hObj,'Value');
xScale = userScale;

%Get scale & clear any previous markers
switch(ud.plotMode)
    case {1,2}
        y1PlotAxis = ud.handle.modelsAxis(1); %NTCP axis
        y2PlotAxis = [];
        maxDelFrx = round(max([protS.numFractions])/2); %rounded
    case 3
        y1PlotAxis = ud.handle.modelsAxis(2); %NTCP axis
        y2PlotAxis = ud.handle.modelsAxis(3); %TCP/BED axis
        fxSizScaleV = linspace(0.5,1.5,99);
        xIdx = abs(fxSizScaleV-userScale) < eps;
        if isempty(xIdx)
            return %Invalid scale factor entered
        end
    case 4
        y1PlotAxis = ud.handle.modelsAxis(4); %NTCP axis
        y2PlotAxis = ud.handle.modelsAxis(5); %TCP/BED axis
        maxDelFrx = round(max([protS.numFractions])/2);
end
hScaled_y1 = findall(y1PlotAxis,'type','line','LineStyle','-.');
delete(hScaled_y1);
if ~isempty(y2PlotAxis)
    hScaled_y2 = findall(y2PlotAxis,'type','line','LineStyle','-.');
    delete(hScaled_y2);
end

%Clear previous readouts
if isfield(ud,'scaleDisp')
    set(ud.scaleDisp,'String','');
end
y1ReadoutsH = findobj(y1PlotAxis,'tag','NTCPreadout');
y2ReadoutsH = findobj(y2PlotAxis,'tag','TCPBEDreadout');
yReadoutsH = [y1ReadoutsH;y2ReadoutsH];
for nLabel = length(yReadoutsH):-1:1
    delete(yReadoutsH(nLabel));
end

if isfield(ud,'y1Disp')
    ud.y1Disp = gobjects;
end
if isfield(ud,'y2Disp')
     ud.y2Disp = gobjects;
end

xLmtV = get(y1PlotAxis,'xLim');
hDisp_y1 = text(xLmtV(1),0,'','Parent',y1PlotAxis,...
    'FontSize',8,'Color',[.3 .3 .3],'HitTest','on',...
    'PickableParts','visible');
if ~isempty(y2PlotAxis)
    hDisp_y2 = text(xLmtV(2),0,'','Parent',y2PlotAxis,...
        'FontSize',8,'Color',[.3 .3 .3],'HitTest','on',...
        'PickableParts','visible');
end

%Set color order
colorM = ud.plotColorOrderM;

%Scale plots as selected
modNum = 0;
y1 = 0;
y2 = 0;
for l = 1:numel(ud.Protocols)

    nMod = length(ud.Protocols(l).model);
    if l == ud.foreground
        pColorM = [colorM,ones(size(colorM,1),1)];
    else
        wt = 0.4;
        pColorM = [colorM,repmat(wt,size(colorM,1),1)];
    end

    %Get plan no.
    planNum = ud.planNum;

    %Loop over models
    for k = 1:nMod
        modNum = modNum+1;

        % Get params
        modelsC = ud.Protocols(l).model;
        paramsS = modelsC{k}.parameters;

        % Get struct
        strNum = modelsC{k}.strNum;
        paramsS.structNum = strNum;

        % Get plan
        paramsS.planNum = planNum;
        paramsS.numFractions.val = ud.Protocols(l).numFractions;
        paramsS.frxSize.val = ud.Protocols(l).totalDose/ud.Protocols(l).numFractions;
        if isfield(modelsC{k},'abRatio')
            paramsS.abRatio.val = modelsC{k}.abRatio;
        end

        % Get dose bins
        if isfield(modelsC{k},'dv')
            dose0C = modelsC{k}.dv{1};
            vol0C = modelsC{k}.dv{2};

            %Scale
            if ud.plotMode==3
                scdoseC = cellfun(@(x) x*userScale,dose0C,'un',0);
                paramsS.frxSize.val = userScale*paramsS.frxSize.val;
            else
                nFProtocol = paramsS.numFractions.val;
                scNumFrx = userScale + nFProtocol;
                nfrxV = linspace(-maxDelFrx,maxDelFrx,99);
                [~,xIdx] = min(abs(nfrxV-userScale));
                if isempty(xIdx)
                    return %Invalid scale factor entered
                end
                paramsS.numFractions.val = scNumFrx;
                scdoseC = cellfun(@(x) x*scNumFrx/nFProtocol,dose0C,'un',0);
            end
            %Apply fractionation correction where required
            eqScaledDoseC = frxCorrectROE(modelsC{k},strNum,paramsS.numFractions.val,scdoseC);

            % Pass as vector if nStr==1
            if numel(strNum) == 1
                vol0C = vol0C{1};
                eqScaledDoseC = eqScaledDoseC{1};
            end

            % Compute probability
            cpNew = feval(modelsC{k}.function,paramsS,eqScaledDoseC,vol0C);

            % Set plot color
            clrIdx = mod(k,size(pColorM,1))+1;

            if strcmpi(modelsC{k}.type,'NTCP') %y1 axis
                currAx = 'y1';
                loc = hObj.Min;
                hplotAx = y1PlotAxis;
                y1 = y1+1;
                count = y1;
                hDisp_y1(count) = text(xLmtV(1),0,'','Parent',y1PlotAxis,...
                   'FontSize',8,'Color',[0 0 0],'BackgroundColor',[1 1 1],...
                   'EdgeColor',pColorM(clrIdx,:),'LineWidth',1.5,...
                   'FontWeight','Bold','Tag','NTCPreadout');
                txtPos = xLmtV(1) - 0.15*abs(xLmtV(1));
                if ud.plotMode==1 || ud.plotMode==2
                    xScale = ud.NTCPCurve(k).XData(xIdx);
                end
                skip=0;
            else %y2 axis
                currAx = 'y2';
                if ud.plotMode==1 || ud.plotMode==2
                    %Skip
                    skip = 1;
                else
                    loc = hObj.Max;
                    hplotAx = y2PlotAxis;
                    y2 = y2+1;
                    count = y2;
                    hDisp_y2(count) = text(xLmtV(2),0,'','Parent',y2PlotAxis,...
                        'FontSize',8,'Color',[0 0 0],'BackgroundColor',...
                        [1 1 1],'EdgeColor',pColorM(clrIdx,:),'LineWidth',...
                        1.5,'FontWeight','Bold','Tag','TCPBEDreadout');
                    txtPos = xLmtV(2)+.02;
                    skip=0;
                end
            end

            if ~skip %Error here: TO DO! Check!

                plot([xScale xScale],[0 cpNew],'Color',pColorM(clrIdx,:),...
                    'LineStyle','-.','linewidth',1.5,'parent',hplotAx);
                plot([loc xScale],[cpNew cpNew],'Color',pColorM(clrIdx,:),...
                    'LineStyle','-.','linewidth',1.5,'parent',hplotAx);

                if strcmp(currAx,'y1')
                    set(hDisp_y1(count),'Position',[txtPos,cpNew],'String',sprintf('%.3f',...
                        cpNew),'Edge',pColorM(clrIdx,:),'LineWidth',2,...
                        'FontWeight','Bold');
                    %Make labels draggable
                    draggable(hDisp_y1(count), "v", [0.01 0.1]);
                else
                    set(hDisp_y2(count),'Position',[txtPos,cpNew],'String',sprintf('%.3f',...
                        cpNew),'Edge',pColorM(clrIdx,:),'LineWidth',2,...
                        'FontWeight','Bold');
                    %Make labels draggable
                    draggable(hDisp_y2(count), "v", [0.01 0.1]);
                end

            end

        end
    end
end

scaleVal = sprintf('%.3f',xScale);
hXDisp = text(xScale,-.03,scaleVal,'Parent',y1PlotAxis,...
    'FontSize',8,'Color',[.3 .3 .3]);
ud.scaleDisp = hXDisp;
ud.y1Disp = hDisp_y1;
if ~isempty(y2PlotAxis)
    ud.y2Disp = hDisp_y2;
end

guidata(hFig,ud);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function T=operatormatrix(Op)
%OPERATORMATRIX  Matrix representation of an operator
%   Usage: T=operatormatrix(Op);
%
%   `T=operatormatrix(Op)` returns the matrix representation *T* of the
%   operator *Op*. The operator object *Op* must have been created using
%   |operatornew|.
%
%   See also: operatornew, operator, operatoreigs

if nargin<1
  error('%s: Too few input parameters.',upper(mfilename));
end;

if ~isstruct(Op)
  error('%s: First argument must be a operator definition structure.',upper(mfilename));
end;

T=operator(Op,eye(Op.L));

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% dl = gp_dlogpdf(y_invK_dK_invK_y, trace_invK_dK)

% dl = gp_dlogpdf(dK, invK, invK_y)

function dl = gp_dlogpdf(y_invK_dK_invK_y, trace_invK_dK)
dl = 0.5*y_invK_dK_invK_y - 0.5*trace_invK_dK;
end
% $$$ function dl = gp_dlogpdf(dK, invK, invK_y)
% $$$ dl = 0.5*(invK_y'*dK*invK_y) - 0.5*traceprod(invK,dK);
% $$$ end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function cmatrix = get_condition_assignments(obj)
% Condition assignments
%   - Indicator matrix coding for which columns in X belong to the same
%     modeled condition, and are part of the same HRF fit
%   - There is one set of columns for each condition modeled, and one set of
%     columns for each parametric modulator of each condition
%   - Because parametric modulators may not exist for all conditions, we need
%     to build this dynamically for modulators.
%
% Design matrix build (which calls method get_session_X) builds columns in
% this order:
%
% All within Session:
% Regressors of interest, basis functions within conditions
% Parametric modulators, basis functions within conditions
% Covariates of no interest
%
% Then:
% Baselines (session/run intercepts)
%
% This method is called automatically in the build method.

switch obj.build_method
    case 'Separate sessions'
        nsess = length(obj.Sess);

    case 'Concatenated sessions'
        nsess = 1;

    otherwise
        error('Unknown build method for fmri_model object. See help for valid strings.');

end

for s = 1:nsess
    nconds(s) = length(obj.Sess(s).U);
end
if any(diff(nconds)), error('Variable numbers of conditions; this should have been caught in build.'); end

% cell of indicator matrices for each session separately
all_indic = cell(1, nsess);

for s = 1:nsess

    % regs of interest
    % -----------------------------------------------------------
    sessindic = {};

    for i = 1:nconds(s)

        % Allow for variable number of bfs.
        nbf = size(obj.xBF(i).bf, 2);

        % indicate a set of regressors based on basis set
        sessindic{i} = true(nbf, 1);

    end

    % concatenate into block diagonal
    sessindic = blkdiag(sessindic{:});

    % Modulators
    % -----------------------------------------------------------
    pmindic = {};

    for i = 1:nconds(s)
        is_pm = isfield(obj.Sess(s).U(i), 'P') && ~isempty(obj.Sess(s).U(i).P) && isfield(obj.Sess(s).U(i).P, 'P') && ~isempty(obj.Sess(s).U(i).P.P);

        if is_pm
            % Allow for variable number of bfs.
            nbf = size(obj.xBF(i).bf, 2);

            pmindic{end + 1} = true(nbf, 1);
        end
    end
    if ~isempty(pmindic)
        pmindic = blkdiag(pmindic{:});
    else
        pmindic = [];
    end

    % add them together in block diag form
    all_indic{s} = blkdiag(sessindic, pmindic);

end % session

cmatrix = logical(blkdiag(all_indic{:}));

ncols = size(obj.xX.X, 2);

if size(cmatrix, 1) > ncols
    error('More conditions of interest than rows in obj.xX.X!! Must build xX before running this.');
end

% add rows of zeros to match total number of columns
z = false(ncols - size(cmatrix, 1), size(cmatrix, 2));

cmatrix = [cmatrix; z];

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function maxFig(varargin)

%% Parse input

switch nargin
    case 0
        hf=gcf;
    case 1
        hf=varargin{1};
end

if isempty(hf)
    hf=gcf; %Take the current figure window handle if the handle is empty
end

%% Turn off warning (JavaFrame will become obsolete, this will lead to error in the future)
% warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame');
%
% %% Maximize using JaveFrame
% try
%     pause(0.01);
%     jFrame=get(hf,'JavaFrame');
% %     set(jFrame,'Maximized',1);
%     jFrame.setMaximized(1);
% catch exception
%     rethrow(exception)
% end

%%

jFrame=get(hf,'JavaFrame');
drawnow;  %Next line throws error without pause or drawnow here
jFrame.fHG2Client.setMaximized(true);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% MEEG time-lock searchlight
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm where a participant saw
% images of six object categories.
%
%
% The code presented here can be adapted for other MEEG analyses, but
% there please note:
% * the current examples do not perform baseline corrections or signal
%   normalizations, which may reduce discriminatory power.
%
% Note: running this code requires FieldTrip.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% get timelock data in CoSMoMVPA format

% set configuration
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'meg_obj6');

% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

% load preprocessed data
data_fn=fullfile(data_path,'meg_obj6_s00.mat');
data_tl=load(data_fn);

% Show data_tl
% >@@>
cosmo_disp(data_tl);
% <@@<

%%

% convert to cosmomvpa struct, using cosmo_meeg_dataset
% >@@>
ds=cosmo_meeg_dataset(data_tl);
% <@@<

% show the dataset
% >@@>
cosmo_disp(ds);
% <@@<

%%

% set the targets in ds.sa.targets (trial condition)
% Hint: use the first column from ds.sa.trialinfo
% >@@>
ds.sa.targets=ds.sa.trialinfo(:,1); % 6 categories
% <@@<

% set the chunks in ds.sa.chunks (independent measurements)
% all trials are here considered to be independent, so the chunks
% must all have a different value
% >@@>
nsamples=size(ds.samples,1);
ds.sa.chunks=(1:nsamples)';
% <@@<

% in addition give a label to each trial
index2label={'body','car','face','flower','insect','scene'};
ds.sa.labels=cellfun(@(x)index2label(x),num2cell(ds.sa.targets));

% just to check everything is ok
cosmo_check_dataset(ds);

%% Count number of channels, time points and trials
% >@@>
fprintf('There are %d channels, %d time points and %d trials\n',...
        numel(unique(ds.fa.chan)),numel(unique(ds.fa.time)),...
        size(ds.samples,1));
% <@@<
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Part I: compute difference between faces and scenes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For each time point and sensor; then visualize the results
% for the magnetometer (meg_axial) sensors.

% slice 'ds' twice to get 'ds_face' and 'ds_scene', each with only trials
% from the face and scene categories

% >@@>
ds_face=cosmo_slice(ds,cosmo_match(ds.sa.labels,'face'));
ds_scene=cosmo_slice(ds,cosmo_match(ds.sa.labels,'scene'));
% <@@<

% prepare dataset for output
ds_faceVSscene=cosmo_slice(ds_face,1);
ds_faceVSscene.sa=struct(); % destroy sample attributes

% Compute difference between average of faces versus average of scenes;
% store the result in the samples field of ds_faceVSscene
% >@@>
ds_faceVSscene.samples=mean(ds_face.samples)-mean(ds_scene.samples);
% <@@<

% Convert ds_faceVSscene to a fieldtrip structure and convert
ft_faceVSscene=cosmo_map2meeg(ds_faceVSscene);
%%
% Use FieldTrip to visualize the face versus house contrast
chantype='meg_axial';
layout=cosmo_meeg_find_layout(ds_faceVSscene,'chantype',chantype);

figure();
cfg=struct();
cfg.interactive='yes';
cfg.zlim=[-1 1];
cfg.layout=layout;

% show figure with plots for each sensor
ft_multiplotER(cfg, ft_faceVSscene);
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part 2: run searchlight over time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set MVPA parameters
fprintf('The input has feature dimensions %s\n', ...
                cosmo_strjoin(ds.a.fdim.labels,', '));

% only select relevant time period and sensors
sensor_posterior_axial={'MEG1631', 'MEG1641', 'MEG1731', 'MEG1841', ...
                        'MEG1911', 'MEG1921', 'MEG1941', 'MEG2231', ...
                        'MEG2311', 'MEG2321', 'MEG2341', 'MEG2431', ...
                        'MEG2441', 'MEG2511', 'MEG2531'};

% define the mask
msk=cosmo_dim_match(ds,'time',@(t) t>=-.1 & t<=.4,...
                        'chan',sensor_posterior_axial);

% first slice the dataset, then use cosmo_dim_prune to avoid using
% non-selected data. Assign the result to ds_sel
%%%% >>> Your code here <<< %%%%
ds_sel_orig=cosmo_slice(ds,msk,2);
cosmo_disp(ds_sel_orig);

ft_orig=cosmo_map2meeg(ds_sel_orig);
%%
ds_sel=cosmo_dim_prune(ds_sel_orig);
cosmo_disp(ds_sel);
ft_sel=cosmo_map2meeg(ds_sel);

%%

% define the neighborhood for time with a time radius of 2 time points,
% and assign to time_nbrhood.
% Hint: use cosmo_interval_neighborhood,

time_nbrhood=cosmo_interval_neighborhood(ds_sel,'time','radius',2);
cosmo_disp(time_nbrhood)

%%

%%%% >>> Your code here <<< %%%%

% Define the measure to be cosmo_crossvalidation_measure,
% and assign to measure
measure=@cosmo_crossvalidation_measure;

%%%% >>> Your code here <<< %%%%

%nsamples=numel(ds_sel.samples);
%rp=cosmo_randperm(nsamples);
%ds_sel.sa.targets=ds_sel.sa.targets(rp);

partitions=cosmo_independent_samples_partitioner(ds_sel,...
                                    'fold_count',5,...
                                    'test_ratio',0.2);

cosmo_disp(partitions)
%%
% Define the partitioning scheme using
% cosmo_independent_samples_partitioner, and assign to partitions.
% Use 'fold_count',5 to use 5 folds,
% and use 'test_ratio',.2 to use 20% of the data for testing (and 80% for
% training) in each fold.

%%%% >>> Your code here <<< %%%%

% Use the LDA classifier and the partitions just defined,
measure_args=struct();
measure_args.partitions=partitions;
measure_args.classifier=@cosmo_classify_lda;

ds_sl=cosmo_searchlight(ds_sel,time_nbrhood,measure,measure_args);

plot(ds_sl.a.fdim.values{1},ds_sl.samples)
xlabel('time');
ylabel('classification accuracy');

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Part III: channel-time searchlight
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select limited time period
msk=cosmo_dim_match(ds,'time',@(t) t>=-.1 & t<=.4);

% first slice the dataset using the mask and assign to 'ds_sel'
ds_sel=cosmo_slice(ds,msk,2);

% Now use cosmo_dim_prune to avoid attempts to using
% non-selected data in the searchlight
ds_sel=cosmo_dim_prune(ds_sel);

% Set the seachlight parameters
chan_count=10; % 10 center channels in each searchlight
time_radius=2; % 2*2+1=5 time bins
chan_type='meg_combined_from_planar';

% define the neighborhood for each dimensions
% First, set 'chan_nbrhood' using cosmo_meeg_chan_neighborhood,
% and use the 'chantype' and 'count' parameters to set the channel type
% and the number of sensors in each searchlight.
chan_nbrhood=cosmo_meeg_chan_neighborhood(ds_sel,'chantype',chan_type,...
                                    'count',chan_count);

cosmo_disp(chan_nbrhood);
%%

%%%% >>> Your code here <<< %%%%

% Second, set 'time_nbrhood' using cosmo_interval_neighborhood,
% using the 'time' dimension

time_nbrhood=cosmo_interval_neighborhood(ds_sel,'time','radius',2);

%%%% >>> Your code here <<< %%%%

% cross neighborhoods for chan-time searchlight
% Hint: use cosmo_cross_neighborhood, and use chan_nbrhood and time_nbrhood
% (in that order) in a cell as the second argument
nbrhood=cosmo_cross_neighborhood(ds_sel,{chan_nbrhood,time_nbrhood});

%%%% >>> Your code here <<< %%%%

% print how many neighbors features have on average
nbrhood_nfeatures=cellfun(@numel,nbrhood.neighbors);
fprintf('Features have on average %.1f +/- %.1f neighbors\n', ...
            mean(nbrhood_nfeatures), std(nbrhood_nfeatures));

%%

% set the 'measure' variable to a function handle to the
% split-half correlation measure
measure=@cosmo_correlation_measure;

% Define the partitioning scheme using
% cosmo_independent_samples_partitioner.
% Use 'fold_count',1 to use 1 folds,
% and use 'test_ratio',.5 to use 50% of the data for testing (and 50% for
% training) in the single fold.
partitions=cosmo_independent_samples_partitioner(ds_sel,'fold_count',1,...
                                            'test_ratio',0.5);

%%%% >>> Your code here <<< %%%%
measure_args=struct();
measure_args.partitions=partitions;

%% run searchlight
% run the searchlight using the parameters above, and assign the result
% to a varibale 'ds_sl'

ds_sl=cosmo_searchlight(ds_sel,nbrhood,measure,measure_args);

%%%% >>> Your code here <<< %%%%

%% visualize timeseries results

% deduce layout from output
layout=cosmo_meeg_find_layout(ds_sl);
fprintf('The output uses layout %s\n', layout.name);

% map ds_sl to a FieldTrip structure. Assign the result to 'sl_ft'
sl_ft=cosmo_map2meeg(ds_sl);

%%%% >>> Your code here <<< %%%%

figure();
cfg = [];
cfg.interactive = 'yes';
cfg.zlim=[-1 1];
cfg.layout       = layout;

% show figure with fisher-transformed correlations for each sensor
ft_multiplotER(cfg, sl_ft);

%% visualize topology results
% show figure with topology for 100 before to 400ms after stimulus onset
% in bins of 50 ms
figure();
cfg.xlim=-0.1:0.05:0.4;
ft_topoplotER(cfg, sl_ft);

%% Show citation information
cosmo_check_external('-cite');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA