text:
function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p,plt,f)
% function [dz,vdz,Adz]=two_group_test_coherence(J1c1,J2c1,J1c2,J2c2,p,plt,f)
% Test the null hypothesis (H0) that data sets J1c1,J2c1,J1c2,J2c2 in
% two conditions c1,c2 have equal population coherence
%
% Usage:
% [dz,vdz,Adz]=two_sample_test_coherence(J1c1,J2c1,J1c2,J2c2,p)
%
% Inputs:
% J1c1   tapered fourier transform of dataset 1 in condition 1
% J2c1   tapered fourier transform of dataset 1 in condition 1
% J1c2   tapered fourier transform of dataset 1 in condition 2
% J2c2   tapered fourier transform of dataset 1 in condition 2
% p      p value for test (default: 0.05)
% plt    'y' for plot and 'n' for no plot
% f      frequencies (useful for plotting)
%
%
% Dimensions: J1c1,J2c2: frequencies x number of samples in condition 1
%              J1c2,J2c2: frequencies x number of samples in condition 2
%              number of samples = number of trials x number of tapers
% Outputs:
% dz    test statistic (will be distributed as N(0,1) under H0
% vdz   Arvesen estimate of the variance of dz
% Adz   1/0 for accept/reject null hypothesis of equal population
%       coherences based dz ~ N(0,1)
%
% Note: all outputs are functions of frequency
%
% References: Arvesen, Jackkknifing U-statistics, Annals of Mathematical
% Statisitics, vol 40, no. 6, pg 2076-2100 (1969)

if nargin < 4; error('Need four sets of Fourier transforms'); end;
if nargin < 6 || isempty(plt); plt='n'; end;

%
% Test for matching dimensionalities
%
if size(J1c1)~=size(J2c1) | size(J1c2)~=size(J2c2) | size(J1c1,1)~=size(J1c2,1);
    error('Need matching dimensionalities for the Fourier transforms: Check the help file for correct dimensionalities');
else;
    m1=size(J1c1,2); % number of samples, condition 1
    m2=size(J1c2,2); % number of samples, condition 2
    dof1=2*m1; % number of degrees of freedom in the first condition estimates
    dof2=2*m2; % number of degrees of freedom in the second condition estimates
end;
if nargin < 7 || isempty(f); f=size(J1c1,1); end;
if nargin < 5 || isempty(p); p=0.05; end; % set the default p value

%
% Compute the individual condition spectra, coherences
%
S12c1=conj(J1c1).*J2c1; % individual sample cross-spectrum, condition 1
S12c2=conj(J1c2).*J2c2; % individual sample cross-spectrum, condition 2
S1c1=conj(J1c1).*J1c1; % individual sample spectrum, data 1, condition 1
S2c1=conj(J2c1).*J2c1; % individual sample spectrum, data 2, condition 1
S1c2=conj(J1c2).*J1c2; % individual sample spectrum, data 1, condition 2
S2c2=conj(J2c2).*J2c2; % individual sample spectrum, data 2, condition 2

Sm12c1=squeeze(mean(S12c1,2)); % mean cross spectrum, condition 1
Sm12c2=squeeze(mean(S12c2,2)); % mean cross spectrum, condition 2
Sm1c1=squeeze(mean(S1c1,2)); % mean spectrum, data 1, condition 1
Sm2c1=squeeze(mean(S2c1,2)); % mean spectrum, data 2, condition 1
Sm1c2=squeeze(mean(S1c2,2)); % mean spectrum, data 1, condition 1
Sm2c2=squeeze(mean(S2c2,2)); % mean spectrum, data 2, condition 1

Cm12c1=abs(Sm12c1./sqrt(Sm1c1.*Sm2c1)); % mean coherence, condition 1
Cm12c2=abs(Sm12c2./sqrt(Sm1c2.*Sm2c2)); % mean coherence, condition 2

Ccm12c1=Cm12c1; % mean coherence saved for output
Ccm12c2=Cm12c2; % mean coherence saved for output
%
% Compute the statistic dz, and the probability of observing the value dz
% given an N(0,1) distribution i.e. under the null hypothesis
%
z1=atanh(Cm12c1)-1/(dof1-2); % Bias-corrected Fisher z, condition 1
z2=atanh(Cm12c2)-1/(dof2-2); % Bias-corrected Fisher z, condition 2
dz=(z1-z2)/sqrt(1/(dof1-2)+1/(dof2-2)); % z statistic
%
% The remaining portion of the program computes Jackknife estimates of the mean (mdz) and variance (vdz) of dz
%
samples1=[1:m1];
samples2=[1:m2];
%
% Leave one out of one sample
%
for i=1:m1;
    ikeep=setdiff(samples1,i); % all samples except i
    Sm12c1=squeeze(mean(S12c1(:,ikeep),2)); % 1 drop mean cross-spectrum, condition 1
    Sm1c1=squeeze(mean(S1c1(:,ikeep),2)); % 1 drop mean spectrum, data 1, condition 1
    Sm2c1=squeeze(mean(S2c1(:,ikeep),2)); % 1 drop mean spectrum, data 2, condition 1
    Cm12c1(:,i)=abs(Sm12c1./sqrt(Sm1c1.*Sm2c1)); % 1 drop coherence, condition 1
    z1i(:,i)=atanh(Cm12c1(:,i))-1/(dof1-4); % 1 drop, bias-corrected Fisher z, condition 1
    dz1i(:,i)=(z1i(:,i)-z2)/sqrt(1/(dof1-4)+1/(dof2-2)); % 1 drop, z statistic, condition 1
    ps1(:,i)=m1*dz-(m1-1)*dz1i(:,i);
%      ps1(:,i)=dof1*dz-(dof1-2)*dz1i(:,i);
end;
ps1m=mean(ps1,2);
for j=1:m2;
    jkeep=setdiff(samples2,j); % all samples except j
    Sm12c2=squeeze(mean(S12c2(:,jkeep),2)); % 1 drop mean cross-spectrum, condition 2
    Sm1c2=squeeze(mean(S1c2(:,jkeep),2)); % 1 drop mean spectrum, data 1, condition 2
    Sm2c2=squeeze(mean(S2c2(:,jkeep),2)); % 1 drop mean spectrum, data 2, condition 2
    Cm12c2(:,j)=abs(Sm12c2./sqrt(Sm1c2.*Sm2c2)); % 1 drop coherence, condition 2
    z2j(:,j)=atanh(Cm12c2(:,j))-1/(dof2-4); % 1 drop, bias-corrected Fisher z, condition 2
    dz2j(:,j)=(z1-z2j(:,j))/sqrt(1/(dof1-2)+1/(dof2-4)); % 1 drop, z statistic, condition 2
    ps2(:,j)=m2*dz-(m2-1)*dz2j(:,j);
%      ps2(:,j)=dof2*dz-(dof2-2)*dz2j(:,j);
end;

%
% Leave one out, both samples
% and pseudo values
% for i=1:m1;
%     for j=1:m2;
%         dzij(:,i,j)=(z1i(:,i)-z2j(:,j))/sqrt(1/(dof1-4)+1/(dof2-4));
%         dzpseudoval(:,i,j)=m1*m2*dz-(m1-1)*m2*dz1i(:,i)-m1*(m2-1)*dz2j(:,j)+(m1-1)*(m2-1)*dzij(:,i,j);
% %         dzpseudoval(:,i,j)=dof1*dof2*dz-(dof1-2)*dof2*dz1i(:,i)-dof1*(dof2-2)*dz2j(:,j)+(dof1-2)*(dof2-2)*dzij(:,i,j);
%     end;
% end;
% dzah=sum(sum(dzpseudoval,3),2)/(m1*m2);
ps2m=mean(ps2,2);
% dzar=(sum(ps1,2)+sum(ps2,2))/(m1+m2);
vdz=sum((ps1-ps1m(:,ones(1,m1))).*(ps1-ps1m(:,ones(1,m1))),2)/(m1*(m1-1))+sum((ps2-ps2m(:,ones(1,m2))).*(ps2-ps2m(:,ones(1,m2))),2)/(m2*(m2-1));
% vdzah=sum(sum((dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))).*(dzpseudoval-dzah(:,ones(1,m1),ones(1,m2))),3),2)/(m1*m2);
%
% Test whether H0 is accepted at the specified p value
%
Adz=zeros(size(dz));
x=norminv([p/2 1-p/2],0,1);
indx=find(dz>=x(1) & dz<=x(2));
Adz(indx)=1;

if strcmp(plt,'y');
    if isempty(f) || nargin < 6;
        f=linspace(0,1,length(dz));
    end;
    %
    % Compute the coherences
    %
    S121=mean(conj(J1c1).*J2c1,2);
    S122=mean(conj(J1c2).*J2c2,2);
    S111=mean(conj(J1c1).*J1c1,2);
    S221=mean(conj(J2c1).*J2c1,2);
    S112=mean(conj(J1c2).*J1c2,2);
    S222=mean(conj(J2c2).*J2c2,2);
    C121=abs(S121)./sqrt(S111.*S221);
    C122=abs(S122)./sqrt(S112.*S222);
    %
    % Plot the coherence
    %
    subplot(311);
    plot(f,C121,f,C122); legend('Data 1','Data 2');
    set(gca,'FontName','Times New Roman','Fontsize', 16);
    ylabel('Coherence');
    title('Two group test for coherence');
    subplot(312);
    plot(f,dz);
    set(gca,'FontName','Times New Roman','Fontsize', 16);
    ylabel('Test statistic');
    conf=norminv(1-p/2,0,1);
    line(get(gca,'xlim'),[conf conf]);
    line(get(gca,'xlim'),[-conf -conf]);
    subplot(313);
    plot(f,vdz);
    set(gca,'FontName','Times New Roman','Fontsize', 16);
    xlabel('frequency'); ylabel('Jackknifed variance');
end;
% Adzar=zeros(size(dzar));
% indx=find(dzar>=x(1) & dzar<=x(2));
% Adzar(indx)=1;
%
% Adzah=zeros(size(dzah));
% indx=find(dzah>=x(1) & dzah<=x(2));
% Adzah(indx)=1;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function TetraFile = fem_hexa2tetra(HexaFile)
% FEM_HEXA2TETRA: Converts hexahedral mesh to tetrahedral mesh
%
% USAGE: TetraFile = fem_hexa2tetra(HexaFile)

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

% Load file
bst_progress('start', 'Convert FEM mesh', ['Loading file "' HexaFile '"...']);
HexaFile = file_fullpath(HexaFile);
FemMat = load(HexaFile);
% Already tetrahedral
if (size(FemMat.Elements,2) == 4)
    disp(['BST> Warning: Mesh is already tetrahedral: ' HexaFile])
    TetraFile = HexaFile;
    return;
end

% Convert to tetrahedral
bst_progress('text', 'Converting to tetrahedral...');
[tetraElem, tetraNode, tetraLabel] = hex2tet(FemMat.Elements, FemMat.Vertices, FemMat.Tissue, 4);
% Update output structure
FemMat.Vertices = tetraNode;
FemMat.Elements = tetraElem(:, [2 1 3 4]);
FemMat.Tissue = tetraLabel;
FemMat.Comment = sprintf('FEM %dV (hexa2tetra, %d layers)', length(FemMat.Vertices), length(FemMat.TissueLabels));
% Add history
FemMat = bst_history('add', FemMat, 'fem_hexa2tetra', 'Converted to tetrahedral');

% Output filename
[fPath, fBase, fExt] = bst_fileparts(HexaFile);
TetraFile = file_unique(bst_fullfile(fPath, [fBase, '_hexa', fExt]));
% Get subject
[sSubject, iSubject] = bst_get('SurfaceFile', HexaFile);
% Save new surface in Brainstorm format
bst_progress('text', 'Saving tetra mesh...');
bst_save(TetraFile, FemMat, 'v7');
db_add_surface(iSubject, TetraFile, FemMat.Comment);

% Close progress bar
bst_progress('stop');

author: brainstorm-tools


text:
function [h,ax] = plot3d(v,data,varargin)
% plot spherical data
%
% Syntax
%   plot3d(v,data)
%
% Input
%
% See also
% savefigure

% -------------------- GET OPTIONS ----------------------------------------

% where to plot
if check_option(varargin,'parent')
  ax = get_option(varargin,'parent');
else
  ax = gca;
end

% scale and shift if required
v = v .* get_option(varargin,'scale',1);
v = v + get_option(varargin,'shift',0);

% plot
%v = v .* reshape(data,size(v));
h = surf(v.x,v.y,v.z,reshape(double(data),size(v,1),size(v,2),[]),'parent',ax,...
  'edgeColor','none');

% colormap
if numel(v) == numel(data), mtexColorMap(ax,getMTEXpref('defaultColorMap')); end

if ~ishold

  axis(ax,'equal','vis3d','off');

  set(ax,'XDir','rev','YDir','rev',...
    'XLim',[-1,1],'YLim',[-1,1],'ZLim',[-1,1]);
end

author: mtex-toolbox

repo: mtex

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


text:
function [heatMaps, prediction] = applyModel_human3d(test_image, param, rectangle, modelID)

caffe.reset_all();

%% check to use click mode or not
click = param.click;

%% select model, 1 for pose estimation, 2 for person detection
model = param.model;
model = model(modelID);
boxsize = model.boxsize;
np = model.np;

%%
oriImg = imread(test_image);
makeFigure = 0;
octave = param.octave;
starting_range = param.starting_range;
ending_range = param.ending_range;
assert(starting_range <= ending_range, 'starting ratio should <= ending ratio');
assert(octave>=1, 'octave should >= 1');

starting_scale = boxsize/(size(oriImg,1)*ending_range);
ending_scale = boxsize/(size(oriImg,1)*starting_range);
multiplier = 2.^(log2(starting_scale):(1/octave):log2(ending_scale));

% set the center and roughly scale range (overwrite the config!) according to the rectangle
    x_start = max(rectangle(1), 1);
    x_end = min(rectangle(1)+rectangle(3), size(oriImg,2));
    y_start = max(rectangle(2), 1);
    y_end = min(rectangle(2)+rectangle(4), size(oriImg,1));

    middle_range = (y_end - y_start) / size(oriImg,1);
    starting_range = middle_range - 0.25;
    ending_range = middle_range + 0.25;

    center = [(x_start + x_end)/2, (y_start + y_end)/2];

    starting_scale = boxsize/(size(oriImg,1)*ending_range);
    ending_scale = boxsize/(size(oriImg,1)*starting_range);
    multiplier = 2.^(log2(starting_scale):(1/octave):log2(ending_scale));

% data container for each scale
score = cell(1);;
peakValue = zeros(length(multiplier), np+1);
pad = cell(1, length(multiplier));
ori_size = cell(1, length(multiplier));

for m = 1
    scale = multiplier(m);
    imageToTest = imresize(oriImg, scale);
    ori_size{m} = size(imageToTest);
    center_s = center * scale;
    [imageToTest, pad{m}] = padAround(imageToTest, boxsize, center_s, model.padValue);
    imageToTest = preprocess(imageToTest, 0.5);
    if(m==1 || ~click)
        caffe.reset_all();
        system(sprintf('sed -i "4s/.*/input_dim: %d/" %s', size(imageToTest,2), model.deployFile));
        system(sprintf('sed -i "5s/.*/input_dim: %d/" %s', size(imageToTest,1), model.deployFile));
        net = caffe.Net(model.deployFile, model.caffemodel, 'test');
    end

    score{m} = applyDNN(imageToTest, net);

    pool_time = size(imageToTest,1) / size(score{m},1);
    score{m} = imresize(score{m}, pool_time);
    score{m} = resizeIntoScaledImg(score{m}, pad{m});
    score{m} = imresize(score{m}, [size(oriImg,2) size(oriImg,1)]);

    if(makeFigure)
        title(sprintf('Current Scale: %f, TOTAL: %f', multiplier(m), sum(peakValue(m,1:np))));
    end
end

%% make heatmaps into the size of scaled image according to pad

    final_score = zeros(size(score{1,1}));
    for m = 1:size(score,2)
        final_score = final_score + score{m};
    end
    heatMaps = permute(final_score, [2 1 3]);

    % generate prediction
    prediction = zeros(np,2);
    for j = 1:np
        [prediction(j,1), prediction(j,2)] = findMaximum(final_score(:,:,j));
    end

function img_out = preprocess(img, mean)
    img_out = double(img)/256;
    img_out = double(img_out) - mean;
    img_out = permute(img_out, [2 1 3]);
    img_out = img_out(:,:,[3 2 1]);
    boxsize = 368;
    centerMapCell = produceCenterLabelMap([boxsize boxsize], boxsize/2, boxsize/2);
    img_out(:,:,4) = centerMapCell{1};

function scores = applyDNN(images, net)
    input_data = {single(images)};
    % do forward pass to get scores
    % scores are now Width x Height x Channels x Num
    s_vec = net.forward(input_data);
    scores = s_vec{1}; % note this score is transposed

function [img_padded, pad] = padAround(img, boxsize, center, padValue)
    center = round(center);
    h = size(img, 1);
    w = size(img, 2);
    pad(1) = boxsize/2 - center(2); % up
    pad(3) = boxsize/2 - (h-center(2)); % down
    pad(2) = boxsize/2 - center(1); % left
    pad(4) = boxsize/2 - (w-center(1)); % right

    pad_up = repmat(img(1,:,:), [pad(1) 1 1])*0 + padValue;
    img_padded = [pad_up; img];
    pad_left = repmat(img_padded(:,1,:), [1 pad(2) 1])*0 + padValue;
    img_padded = [pad_left img_padded];
    pad_down = repmat(img_padded(end,:,:), [pad(3) 1 1])*0 + padValue;
    img_padded = [img_padded; pad_down];
    pad_right = repmat(img_padded(:,end,:), [1 pad(4) 1])*0 + padValue;
    img_padded = [img_padded pad_right];

    center = center + [max(0,pad(2)) max(0,pad(1))];

    img_padded = img_padded(center(2)-(boxsize/2-1):center(2)+boxsize/2, center(1)-(boxsize/2-1):center(1)+boxsize/2, :); %cropping if needed

function [img_padded, pad] = padRightDownCorner(img, padValue)
    h = size(img, 1);
    w = size(img, 2);

    pad(1) = 76; % up
    pad(2) = 76; % left
    pad(3) = 76 + 8 - mod((152+h), 8); % down
    pad(4) = 76 + 8 - mod((152+w), 8); % right

    img_padded = img;
    pad_up = repmat(img_padded(1,:,:), [pad(1) 1 1])*0 + padValue;
    img_padded = [pad_up; img_padded];
    pad_left = repmat(img_padded(:,1,:), [1 pad(2) 1])*0 + padValue;
    img_padded = [pad_left img_padded];
    pad_down = repmat(img_padded(end,:,:), [pad(3) 1 1])*0 + padValue;
    img_padded = [img_padded; pad_down];
    pad_right = repmat(img_padded(:,end,:), [1 pad(4) 1])*0 + padValue;
    img_padded = [img_padded pad_right];

function [x,y] = findMaximum(map)
    [~,i] = max(map(:));
    [x,y] = ind2sub(size(map), i);

function score = resizeIntoScaledImg(score, pad)
    np = size(score,3)-1;
    score = permute(score, [2 1 3]);
    if(pad(1) < 0)
        padup = cat(3, zeros(-pad(1), size(score,2), np), zeros(-pad(1), size(score,2), 1));
        score = [padup; score]; % pad up
    else
        score(1:pad(1),:,:) = []; % crop up
    end

    if(pad(2) < 0)
        padleft = cat(3, zeros(size(score,1), -pad(2), np), zeros(size(score,1), -pad(2), 1));
        score = [padleft score]; % pad left
    else
        score(:,1:pad(2),:) = []; % crop left
    end

    if(pad(3) < 0)
        paddown = cat(3, zeros(-pad(3), size(score,2), np), zeros(-pad(3), size(score,2), 1));
        score = [score; paddown]; % pad down
    else
        score(end-pad(3)+1:end, :, :) = []; % crop down
    end

    if(pad(4) < 0)
        padright = cat(3, zeros(size(score,1), -pad(4), np), zeros(size(score,1), -pad(4), 1));
        score = [score padright]; % pad right
    else
        score(:,end-pad(4)+1:end, :) = []; % crop right
    end
    score = permute(score, [2 1 3]);

function label = produceCenterLabelMap(im_size, x, y) %this function is only for center map in testing
    sigma = 21;
    %label{1} = zeros(im_size(1), im_size(2));
    [X,Y] = meshgrid(1:im_size(1), 1:im_size(2));
    X = X - x;
    Y = Y - y;
    D2 = X.^2 + Y.^2;
    Exponent = D2 ./ 2.0 ./ sigma ./ sigma;
    label{1} = exp(-Exponent);

author: flyawaychase


text:
function newmodel = addMultipleMetabolites(model, metIDs, varargin)
% Simultaneously add multiple metabolites. Allows the use of the field names of the model struct to specify additional properties of the metabolites added.
%
% USAGE:
%
%    newmodel = addMultipleMetabolites(model, metIDs, varargin)
%
% INPUTS:
%    model:         The model to add the Metabolite batch to.
%    metIDs:        The metabolite IDs to add. No duplicate IDs may be provided
%                   and no ID may be present in the supplied model.
%
% OPTIONAL INPUTS:
%    varargin:      fieldName, Value pairs with additional properties for the
%                   added metabolites. In addition, a 'printLevel' flag can be used (default: 0).
%                   The given values fields will be set according to the values.
%                   Only fields associated with mets defined in the COBRA definitions (except `S`)
%                   or fields already in the model may be used.
%                   Examples for this use would be:
%
%                     * 'metNames',{'Glucose';'Pyruvate'}
%                     * 'metCharges',[0;-2]
%
% OUTPUTS:
%
%    newmodel:     The model structure with the additional metabolites.
%
% EXAMPLE:
%
%    % To add metabolites, with charges, formulas and KEGG ids:
%    model = addMultipleMetabolites(model,{'A','b','c'},'metCharges', [ -1 1 0], 'metFormulas', {'C','CO2','H2OKOPF'}, 'metKEGGID',{'C000012','C000023','C000055'})
%

if numel(unique(metIDs)) < numel(metIDs)
    % check, if there are either duplicate metabolite IDS to be added OR if
    % any metabolite is already in the model.
    error('Duplicate Metabolite ID in the given IDs detected.');
end

if checkIDsForTypeExist(model,metIDs,'mets')
    % check, if there are either duplicate metabolite IDS to be added OR if
    % any metabolite is already in the model.
    [tf,dups] = checkIDsForTypeExist(model,metIDs,'mets');
    if any(ismember(model.mets,dups))
        pres = ismember(model.mets,dups);
        error('The following Metabolite ID(s) are already present in the model:\n%s',strjoin(model.mets(pres),'\n'));
    else
        error('The following Metabolite ID(s) are already present Ids of constraints in the model:\n%s',strjoin(dups,'\n'));
    end
end

nMets = numel(model.mets);

%We extract those fields which are either associated with the mets field
%(2nd or 3rd column contains 'mets' in the definitions), and we also look up which fields are in the model and
%associated with the mets field (from the sizes)
fieldDefs = getDefinedFieldProperties();
fieldDefs = fieldDefs(cellfun(@(x) strcmp(x,'mets'), fieldDefs(:,2)) | cellfun(@(x) strcmp(x,'mets'), fieldDefs(:,3)));
modelMetFields = getModelFieldsForType(model,'mets');

% extract additional parameters which are not fields (currently only
% printLevel
printLevel = 0;

for parameter = 1:2:numel(varargin)
    if strcmp(varargin{parameter},'printLevel')
        printLevel = varargin{parameter+1};
    end
end

%Then we add the ids.
model.mets = [model.mets;columnVector(metIDs)];
if printLevel > 0
    fprintf('Adding the following Metabolites to the model:\n%s\n',strjoin(metIDs,'\n'));
end

%Now, add the the data from the additional supplied fields, and check that
%they are in sync
for field = 1:2:numel(varargin)
    %If the field does not exist and is not defined, we will not add that
    %information, as we don't know how the field should look.
    cfield = varargin{field};
    if strcmp('S',cfield) || (~any(ismember(fieldDefs(:,1),cfield)) && ~any(ismember(modelMetFields,cfield)))
        if printLevel > 0 && ~strcmp('printLevel',cfield)
            warning('Field %s is excluded',cfield);
        end
        continue;
    end
    %Now add the field data.
    if ~isfield(model,cfield)
        %If its not yet in the model, create it
        %according to its defaults and replace the final elements.
        if printLevel > 2
            fprintf('Creating model field %s\n',cfield);
        end
        model = createEmptyFields(model,cfield);
        model.(cfield)((end-numel(varargin{field+1})+1):end) = columnVector(varargin{field+1});
    else
        %Or just extend it with the supplied data.0
        model.(cfield) = [model.(cfield);columnVector(varargin{field+1})];
    end

end

%Extend the remaining fields, to keep the fields in sync
newmodel = extendModelFieldsForType(model,'mets','originalSize',nMets);

author: opencobra

repo: cobratoolbox


text:
% dipfit_gridsearch() - do initial batch-like dipole scan and fit to all
%                       data components and return a dipole model with a
%                       single dipole for each component.
%
% Usage:
%  >> EEGOUT = dipfit_gridsearch( EEGIN, varargin)
%
% Inputs:
%    ...
%
% Optional inputs:
%   'component' - vector with integers, ICA components to scan
%   'xgrid'     - vector with floats, grid positions along x-axis
%   'ygrid'     - vector with floats, grid positions along y-axis
%   'zgrid'     - vector with floats, grid positions along z-axis
%
% Output:
%    ...
%
% Author: Robert Oostenveld, SMI/FCDC, Nijmegen 2003, load/save by
%         Arnaud Delorme
%         Thanks to Nicolas Robitaille for his help on the CTF MEG
%         implementation

% SMI, University Aalborg, Denmark http://www.smi.auc.dk/
% FC Donders Centre, University Nijmegen, the Netherlands http://www.fcdonders.kun.nl

% Copyright (C) 2003 Robert Oostenveld, SMI/FCDC roberto@smi.auc.dk
%
% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [EEGOUT] = dipfit_gridsearch(EEG, varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert the optional arguments into a configuration structure that can be
% understood by FIELDTRIPs dipolefitting function
if nargin>2
  cfg = struct(varargin{:});
else
  help dipfit_gridsearch
  return
end

% specify the FieldTrip DIPOLEFITTING configuration
cfg.model      = 'moving';
cfg.gridsearch = 'yes';
cfg.nonlinear  = 'no';
% add some additional settings from EEGLAB to the configuration
tmpchanlocs    = EEG.chanlocs;
cfg.channel    = { tmpchanlocs(EEG.dipfit.chansel).labels };
if isfield(EEG.dipfit, 'vol')
    cfg.vol        = EEG.dipfit.vol;
elseif isfield(EEG.dipfit, 'hdmfile')
    cfg.hdmfile    = EEG.dipfit.hdmfile;
else
    error('no head model in EEG.dipfit')
end
if isfield(EEG.dipfit, 'elecfile') & ~isempty(EEG.dipfit.elecfile)
    cfg.elecfile = EEG.dipfit.elecfile;
end
if isfield(EEG.dipfit, 'gradfile') & ~isempty(EEG.dipfit.gradfile)
    cfg.gradfile = EEG.dipfit.gradfile;
end

% convert the EEGLAB data structure into a structure that looks as if it
% was computed using FIELDTRIPs componentanalysis function
comp = eeglab2fieldtrip(EEG, 'componentanalysis', 'dipfit');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Added code to handle CTF data with multipleSphere head model           %
%  This code is copy-pasted in dipfit_gridSearch, dipfit_nonlinear        %
%  The flag .isMultiSphere is used by dipplot                             %
%  Nicolas Robitaille, January 2007.                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Do some trick to force fieldtrip to use the multiple sphere model
if strcmpi(EEG.dipfit.coordformat, 'CTF')
   cfg = rmfield(cfg, 'channel');
   comp = rmfield(comp, 'elec');
   cfg.gradfile = EEG.dipfit.chanfile;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(cfg, 'component')
  % default is to scan all components
  cfg.component = 1:size(comp.topo,2);
end

% for each component scan the whole brain with dipoles using FIELDTRIPs
% dipolefitting function
source = ft_dipolefitting(cfg, comp);

% reformat the output dipole sources into EEGLABs data structure
for i=1:length(cfg.component)
  EEG.dipfit.model(cfg.component(i)).posxyz = source.dip(i).pos;
  EEG.dipfit.model(cfg.component(i)).momxyz = reshape(source.dip(i).mom, 3, length(source.dip(i).mom)/3)';
  EEG.dipfit.model(cfg.component(i)).rv     = source.dip(i).rv;
end

EEGOUT = EEG;

author: PatternRecognition

repo: OpenBMI


text:
function x=make_odd(x)
if(~mod(x,2)), x=x-1; end

author: alexandrebarachant


text:
function[varargout]=mspec(varargin)
% MSPEC  Multitaper power and cross spectra.
%   _______________________________________________________________________
%
%   *|* mspec.png --- Figure illustrating the multitaper spectrum.
%   Type 'jhelp mspec' to view this image. *|*
%   _______________________________________________________________________
%
%   MSPEC implements spectral and cross-spectral analysis using the multi-
%   taper method for real or complex-valued data, and along any dimension.
%
%   MSPEC is to be run after calling SLEPTAP to compute the multitapers.
%
%   Confidence intervals can be computed by calling MCONF.
%   _______________________________________________________________________
%
%   Real-valued time series
%
%   [F,S]=MSPEC(X,PSI) returns the power spectrum of the real-valued time
%   series X at positive frequencies using data tapers PSI.  The spectrum
%   at negative frequencies, which is not returned, is identical.
%
%   X may be an array with an arbitrary number of dimensions and with time
%   along it first dimensions.  Let M be its number of rows, M=SIZE(X,1).
%
%   PSI is a then matrix of K data tapers having M rows and K columns.
%
%   F is an array of frequencies with FlOOR(M/2)+1 rows, while the spectral
%   matrix S has FlOOR(M/2)+1 rather than M rows, and the same size as X
%   along all of its other dimensions.
%
%   The spectrum can also be computed with time oriented along a different
%   dimension, as described below.
%
%   By default, MSPEC removes the mean from each time series before
%   computing the spectra. This is suppressed by MSPEC(...,'nodemean').
%   ______________________________________________________________________
%
%   Cross-spectra of real-valued data
%
%   [F,SXX,SYY,SXY]=MSPEC(X,Y,PSI) computes the cross-spectrum of two
%   real-valued time series or sets of time series.  Here SXX and SYY are
%   the one-sided spectra of X and Y, while SXY is their cross spectrum.
%
%   See TWOSPECPLOT for plotting SXX and SYY simultaneously.
%   ______________________________________________________________________
%
%   Rotary spectra of complex-valued data
%
%   [F,SPP,SNN,SPN]=MSPEC(Z,PSI) where Z is complex-valued computes the so-
%   called "rotary spectra". Here SPP and SNN are the positively-rotating
%   and negatively rotating spectra, and SPN is the rotary cross spectrum.
%   the one-sided spectra of X and Y, while SPN is their cross spectrum.
%
%   Note that the rotary spectra are defined such that SXX+SYY=SPP+SNN.
%
%   The rotary spectra SPP and SNN are normalized such that the sum of SPP
%   over all frequencies plus that of SNN approximates the variance of Z.
%
%   See TWOSPECPLOT for plotting SPP and SNN simultaneously.
%   ______________________________________________________________________
%
%   Periodogram
%
%   MPSEC can be used to form the naive spectral estimator, known as the
%   periodogram. Although this is not generally a good way to estimate the
%   spectrum, it can be useful as a comparision.
%
%   MSPEC(X,[]) or MSPEC(X,Y,[]) with PSI empty uses the default, or boxcar
%   taper, normalized to unit energy. This returns the periodogram.
%   ______________________________________________________________________
%
%   Sample rate
%
%   [F,S]=MSPEC(DT,...) specifies the sample interval to be used in the
%   calculation of the frequency array F. DT defaults to unity.
%
%   Spectral values depend linearly upon the sample rate in order that the
%   integral of the spectra over frequency approximate the variance.
%   ______________________________________________________________________
%
%   Spectra along arbitary dimension
%
%   MSPEC(...,DIM) computes the spectrum with time oriented along dimension
%   DIM, with the default behavior corresponding to DIM=1.
%
%   Let M be the length of the input X, Z, or X and Y along dimension DIM.
%   PSI is again a matrix of K data tapers having M rows and K columns.
%
%   F will be again an array of frequencies with FlOOR(M/2)+1 rows, while
%   the output spectral matrices will the same size as the input arrays,
%   but with FlOOR(M/2)+1 rather than M elements along dimension DIM.
%   ______________________________________________________________________
%
%   Normalizations
%
%   By default, MSPEC uses *radian* frequency as in cos(f t).  Optionally
%   MSPEC(,...,'cyclic') will use *cyclic* frequency, as in cos(2 pi f t).
%
%   MSPEC is normalized to approximately recover the time series variance.
%   For the MSPEC periodogram, this recovery is exact, although the
%   expressions are complicated somewhat by the use of one-sided spectra.
%
%   For simplicity, the normalizations will be explained for the case of a
%   single time series with M oriented along rows, that is, with DIM=1.
%
%   Real-valued data
%
%   [F,S]=MSPEC(DT,X,[]) where X is a real-valued time series of length M
%   recovers the variance of X, STD(X,1).^2, as follows:
%
%     2*(1/2/pi)*(F(2)-F(1))*SUM(S(2:end))               -- M odd
%     2*(1/2/pi)*(F(2)-F(1))*(SUM(S(2:end-1))+S(end)/2)  -- M even
%
%   where the initial factor of two accounts for the fact that the spectrum
%   at negative frequencies is the same as that at positive frequencies.
%
%   Note that the zero frequency is omitted in the summation, and for even
%   time series length, the power at the Nyquist S(end) must be divided by
%   two to avoid double-counting by the one-sided spectrum.  The "1" in the
%   argument of STD forces STD to use an N rather than N-1 normalization.
%
%   Complex-valued data
%
%   [F,SPP,SNN]=MSPEC(DT,Z,[]) where Z is a complex-valued time series of
%   length M recovers the variance of Z, STD(Z,1).^2, as follows:
%
%     (1/2/pi)*(F(2)-F(1))*(SUM(SPP(2:end))+SUM(SNN(2:end)))   -- M odd
%     (1/2/pi)*(F(2)-F(1))*(SUM(SPP(2:end))+SUM(SNN(2:end-1))) -- M even
%
%   Again the modification for even M prevents the power at the Nyquist
%   from being double-counted.  This modification is necessary because the
%   negative rotary spectrum duplicates the Nyquist when M is even.
%   ______________________________________________________________________
%
%   Cross-spectra of complex-valued data
%
%   To compute the cross-spectra of two complex-valued time series or sets
%   of time series Z1 and Z2, run MSPEC repeatedly.
%
%   [F,SP1P1,SP2P2,SP1P2]=MSPEC(Z1,Z2,PSI);
%   [F,SN1N1,SN2N2,SN1N2]=MSPEC(CONJ(Z1),CONJ(Z2),PSI);
%
%   The first call returns the spectra and cross-spectra of Z1 and Z2 at
%   positive frequencies, while the second returns their spectra and the
%   *conjugate* of the cross-spectrum at negative frequencies.  Finally
%
%   [F,SP1P1,SN2N2,SP1N2]=MSPEC(Z1,CONJ(Z2),PSI);
%   [F,SN1N1,SP2P2,SN1P2]=MSPEC(CONJ(Z1),Z2,PSI);
%
%   returns the so-called outer or complementary cross-spectra.
%   ______________________________________________________________________
%
%   Adaptive spectra
%
%   MSPEC(...,LAMBDA,'adaptive'), where LAMBDA contains the eigenvalues of
%   the tapers as computed by SLEPTAP, alternately uses the "adaptive"
%   multitaper method of Thomson (1982).
%
%   This implementation follows that of Park et al. (1987a), JGR.
%
%   For cross-spectra or for rotary spectra, the weights appearing in the
%   adaptive spectra are derived for the total spectrum of each signal
%   compoment, that is for SXX+SYY or SPP+SNN as appropriate.  Then the
%   separate spectra and co-spectra are computed using identical weights.
%   ______________________________________________________________________
%
%   Cell array input / output
%
%   MSPEC generates cell array output given cell array input.
%
%   Let's say one has P different time series, X1, X2,..., XP.  Put these
%   into a cell array X{1}=X1, X{2}=X2, ..., X{P}=XP, and then use
%   "[psi,lambda]=sleptap(cellength(x))" to make a cell array of tapers.
%
%   [F,S]=MSPEC(X,PSI) then returns cell arrays F and S corresponding
%   to the Fourier frequencies and spectra of the P arrays.
%
%   The other argument forms given above also work.  In particular,
%   specifiying the sample time through MPSEC(DT,...) works, with DT either
%   a scalar or an array of the same length as the cell array X.
%
%   The spectra can then be plotted with CELLPLOT(F,S), or TWOSPECPLOT for
%   a pair of output spectra.
%   ______________________________________________________________________
%
%   Parallelization
%
%   MSPEC(..., 'parallel') when the input fields X, X and Y, or Z are cell
%   arrays, parallelizes the spectral estimation by looping over the cells
%   with a PARFOR loop.  This requires Matlab's Parallel Computing Toolbox.
%   ______________________________________________________________________
%
%   Example
%
%   The example at the top of this help file shows clockwise (left) and
%   counterclockwise (right) rotary spectra from moored current meter
%   measurements of the ocean currents in the Labrador Sea.
%
%   The periodogram is in gray, and blue and red are multitaper spectra
%   with P=4 and P=32, respectively.  The local Coriolis frequency is
%   marked with a dashed line.  Tidal and inertial peaks are apparent.
%
%   The main point of this figure is to show that increasing P increases
%   the degree of frequency-domain smoothing.
%   ______________________________________________________________________
%
%   'mspec --t' runs some tests.
%   'mspec --f' generates the above sample figure from Bravo mooring data.
%
%   See also:  SLEPTAP, MCONF, HERMFUN, MSVD, TWOSPECPLOT.
%
%   Usage   [f,s]=mspec(x,psi);
%           [f,s]=mspec(dt,x,psi);
%           [f,s]=mspec(dt,x,psi,dim);
%           [f,spp,snn,spn]=mspec(z,psi);
%           [f,sxx,syy,sxy]=mspec(x,y,psi);
%           [f,sxx,syy,sxy]=mspec(x,y,psi,dim);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2020 J.M. Lilly --- type 'help jlab_license' for details

%   The Cartesian and rotary spectra are then related by a unitary
%   transformation.  For details see
%
%       Lilly and Olhede (2010).  Bivariate instantaneous frequency and
%           bandwidth.  IEEE Trans. Sig. Proc.

if strcmpi(varargin{1},'--t')
    mspec_test;return
end
if strcmpi(varargin{1},'--f')
    type makefigs_mspec;
    makefigs_mspec;
    return
end

%Sort out input arguments
dim=1;
deltat=1;
lambda=1;  %This means use the average multitaper spectrum
detrendstr='demean';
normstr='rad';
cores='serial';

for i=1:4
    if ischar(varargin{end})
        if ~isempty(strfind(varargin{end},'ada'))
            lambda=varargin{end-1};
            varargin=varargin(1:end-2);
            if length(lambda)~=numel(lambda)
                if lambda~=floor(lambda)
                    error('Looks like you forgot to input LAMBDA with the adaptive algorithm.')
                end
            end
        elseif ~isempty(strfind(varargin{end},'dem'))||~isempty(strfind(varargin{end},'nod'))
             detrendstr=varargin{end};
             varargin=varargin(1:end-1);
        elseif strcmpi(varargin{end}(1:3),'par')||strcmpi(varargin{end}(1:3),'ser')
             cores=varargin{end};
             varargin=varargin(1:end-1);
        else
             normstr=varargin{end};
             varargin=varargin(1:end-1);
        end
    end
end

if strcmpi(cores(1:3),'par')
    if exist('parpool')~=2
        disp('Sorry, parallel algorithm requires the Parallel Computing Toolbox.')
        disp('Defaulting to the standard algorithm.')
        cores='serial';
    end
end

if isscalar(varargin{end})&&~iscell(varargin{end})
    dim=varargin{end};
    varargin=varargin(1:end-1);
end

if isscalar(varargin{1})
  deltat=varargin{1};
  varargin=varargin(2:end);
else
  if length(varargin)>2
      bool1=~iscell(varargin{2})&&(numel(varargin{1})==size(varargin{2},2));
      bool2= iscell(varargin{2})&&(numel(varargin{1})==length(varargin{2}));
      if bool1||bool2
          deltat=varargin{1};
          deltat=deltat(:)';
          varargin=varargin(2:end);
      end
  end
end

x=varargin{1};
if length(varargin)==1
    error('Taper not specified.')
end
psi=varargin{end};
na=length(varargin);

y=[];
if na==3
    y=varargin{2};
end

if iscell(x)||iscell(y)
    %All of this is just to handle cell array input
    if isempty(y)
        y=cell(size(x));
    elseif isempty(x)
        x=cell(size(y));
    end
    if ~iscell(psi)
        psio=psi;
        psi=cell(size(x));
        for i=1:length(x)
            psi{i}=psio;
        end
    end
    if size(lambda,2)==1
        lambda=vrep(lambda,length(x),2);
    end
    if length(deltat)==1
        deltat=deltat+zeros(size(x));
    end

    if strcmpi(cores(1:3),'ser')
        for i=1:length(x)
            disp(['SLEPTAP computing spectra for time series #' int2str(i) ' of ' int2str(length(x)) '.'])
            %iscell(deltat),iscell(x),iscell(y),iscell(psi)
            %vsize(deltat(i),x{i},y{i},psi{i},lambda(:,i))
            cellout{i}=mspec_one(deltat(i),x{i},y{i},psi{i},lambda(:,i),detrendstr,normstr,dim);
        end
    elseif strcmpi(cores(1:3),'par')
        %Exactly the same but with a parfor
        parfor i=1:length(x)
            disp(['SLEPTAP computing spectra for time series #' int2str(i) ' of ' int2str(length(x)) '.'])
            %iscell(deltat),iscell(x),iscell(y),iscell(psi)
            %vsize(deltat(i),x{i},y{i},psi{i},lambda(:,i))
            cellout{i}=mspec_one(deltat(i),x{i},y{i},psi{i},lambda(:,i),detrendstr,normstr,dim);
        end
    end

    for i=1:length(x)
        for j=1:length(cellout{i})
            varargout{j}{i,1}=cellout{i}{j};
        end
    end
else
    varargout=mspec_one(deltat,x,y,psi,lambda,detrendstr,normstr,dim);
end

function[cellout]=mspec_one(dt,x,y,psi,lambda,detrendstr,normstr,dim)

if ~isscalar(dt)
    dt=vrep(dt,length(fourier(size(x,1))),1);
end

if ~isreal(x)&&isempty(y)
    y=conj(x);
end

if strcmpi(detrendstr(1:3),'dem')
    x=x-vrep(vmean(x,dim),size(x,dim),dim);
    if ~isempty(y)
       y=y-vrep(vmean(y,dim),size(y,dim),dim);
    end
end

if isempty(psi)
    psi=frac(1,sqrt(size(x,1)))+zeros(size(x(:,1)));
end

[f,mmatx,mmaty]=mtrans(x,y,psi,dim);

N=lnsd(x)+1;
if isempty(y) %One time series
     if length(lambda)==1
         cellout{2}=avgspec(mmatx,mmatx,N).*dt;
     else
         var=vrep(squared(vstd(x,dim)),floor(size(x,dim)/2)+1,dim);
         %Variance same size as original input field
         cellout{2}=adaptspec(abs(mmatx).^2,lambda,var,N).*dt;
     end
     cellout{3}=zeros(size(cellout{2}));
     cellout{4}=zeros(size(cellout{2}));
else         %Two time series
     if length(lambda)==1
        cellout{2}=avgspec(mmatx,mmatx,N).*dt;
        cellout{3}=avgspec(mmaty,mmaty,N).*dt;
        cellout{4}=avgspec(mmatx,mmaty,N).*dt;
        cellout{4}(isnan(cellout{4}))=nan+1i*nan;%cross-spectrum should have complex nans
     else
        %For two time series one should do the adaptive spectra on both
        %with the same coefficients
        var=vrep(squared(vstd(x,dim))+squared(vstd(y,dim)),floor(size(x,dim)/2)+1,dim);

        [~,dk]=adaptspec(abs(mmatx).^2+abs(mmaty).^2,lambda,var,N);
        cellout{2}=frac(sum(dk.^2.*abs(mmatx).^2.*dt,N),sum(abs(dk).^2,N));
        cellout{3}=frac(sum(dk.^2.*abs(mmaty).^2.*dt,N),sum(abs(dk).^2,N));
        cellout{4}=frac(sum(dk.^2.*mmatx.*conj(mmaty).*dt,N),sum(abs(dk).^2,N));
        cellout{4}(isnan(cellout{4}))=nan+1i*nan;%cross-spectrum should have complex nans
     end
end

%Corrections for zero component, and for Nyquist with even and odd length
%Both zero are Nyquist are shared for even length time series, so divide both by two.
%Only zero is shared for odd length, since the Nyquist does not appear.
%This is best visualized by drawing N equally spaced points on the unit circle.

% for i=2:length(cellout)
%     cellout{i}(1,:)=cellout{i}(1,:)./2;
%     if iseven(size(x,1))
%        cellout{i}(end,:)=cellout{i}(end,:)./2;
%     end
% end

if contains(normstr,'cyc')
    f=f/2/pi;
    %for i=2:length(cellout)
    %    cellout{i}=cellout{i}/2/pi;
    %end
end

if isscalar(dt)
    cellout{1}=f./dt;
else
    cellout{1}=vrep(f,size(dt,2),2)./dt;
end

function[S]=avgspec(mmat1,mmat2,N)
eigspec=mmat1.*conj(mmat2);
S=mean(eigspec,N);

function[s,dk]=adaptspec(eigspec,lambda,var,N)

sold=vindex(eigspec,1,N);
s=frac(1,2)*(vindex(eigspec,1,N)+vindex(eigspec,2,N));
%Start just with first two

tol=1e-4;

var=vrep(var,length(lambda),N);
lambdamat=zeros(size(var));
%vsize(var,lambdamat,eigspec)
for i=1:length(lambda)
    lambdamat=vindexinto(lambdamat,lambda(i),i,N);
end
%figure,plot(lambdamat)
bkmat=var.*(1-lambdamat);  %Outer product;

% sold=s;
% smat=vrep(s,size(eigspec,N),N);
% dk=(smat.*real(sqrt(lambdamat)))./(smat.*lambdamat+bkmat);  %Outer products
% s=frac(sum(dk.^2.*eigspec,N),sum(abs(dk).^2,N));

i=0;
%maxmax(abs(s-sold)./sold)
while anyany(abs(s-sold)./sold>tol)&&i<20
    i=i+1;
	sold=s;
    smat=vrep(s,size(eigspec,N),N);
	dk=(smat.*real(sqrt(lambdamat)))./(smat.*lambdamat+bkmat);  %Outer products
	s=frac(sum(dk.^2.*eigspec,N),sum(abs(dk).^2,N));
 %   maxmax(abs(s-sold)./sold)
end

if i~=20
   disp(['Adaptive spectral estimate took ' int2str(i) ' iterations.'])
else
   disp(['Adaptive spectral loop terminated at ' int2str(i) ' iterations.'])
end

function[varargout]=mtrans(varargin)
% MTRANS  Multitaper "eigentransform" computation.
%
%   [F,W]=MTRANS(X,PSI,DIM) returns the multitaper "eigentransform" matrix
%   for use in multitaper spectral estimates or eigenspectral SVD analysis.
%
%       X  --  M x N matrix containing N length M time series
%     PSI  --  M x K matrix of K data tapers
%       W  --  [M/2] x K x  N eigentransform matrix (for real X)
%              [M/2] x K x 2N eigentransform matrix (for complex X)
%
%   In the above, [M/2] means M/2 if M is even, and (M-1)/2 is M is odd.
%
%   F is the angular Fourier frequency, in radians.
%
%   [F,W1,W2,...,WN]=MTRANS(X1,X2,...,XN,PSI,DIM) also works, where X1,...,
%   XN are all M x N matrices.
%
%   MTRANS(X,[]) with PSI empty uses the default taper, so that the square
%   of W corresponds to the periodogram.
%
%   See also: SLEPTAP, HERMFUN, MSPEC, MSVD.
%
%   Usage:  [f,w]=mtrans(x,psi,dim);
%           [f,wx,wy]=mtrans(x,y,psi,dim);
%           [f,wx,wy,wz]=mtrans(x,y,z,psi,dim);
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2015 J.M. Lilly --- type 'help jlab_license' for details

%Sort out input arguments
dim=varargin{end};
psi=varargin{end-1};
x=varargin(1:end-2);

%Remove empties
clear bool
for i=1:length(x)
    bool(i)=isempty(x{i});
end
x=x(~bool);

%Size check
bool=false(4,length(x));
for i=1:4
    for j=1:length(x)
       bool(i,j)=aresame(size(x{j},i),size(x{1},i));
    end
end
if ~allall(bool)
   error('All input arguments should be the same size')
end

%This complicated code makes a vector that permutes PSI such that it has
%its first dimension along DIM and its last along dimension LNSD(X{1})+1.
n=3;
permutevec=zeros(lnsd(x{1})+1,1);
for i=1:lnsd(x{1})
    if i==dim
        permutevec(i)=1;
    else
        permutevec(i)=n;
        n=n+1;
    end
end
permutevec(end)=2;

if isempty(psi)
    psimat=[];
else
    psimat=permute(psi,permutevec);
    for i=1:lnsd(x{1})
        if i~=dim
            psimat=vrep(psimat,size(x{1},i),i);
        end
    end
end

f=fourier(size(x{1},dim));

index=1:length(f);
varargout{1}=f;
varargout{2}=[];
varargout{3}=[];

for i=1:size(x,2)
    Nnans=length(find(~isfinite(x{i})));
    if Nnans>0
        disp(['MSPEC finding non-finite data values.  Spectrum will be undefined.'])
    end

    N=lnsd(x{i})+1;
    xmat=vrep(x{i},size(psimat,N),N);
    mmat=fft(psimat.*xmat,[],dim);

    mmat(isnan(mmat))=nan;%This is to prevent imaginary NaNs
    %    varargout{i+1}=mmat(index,:,:);
    varargout{i+1}=vindex(mmat,index,dim);
end

function[]=mspec_test

tol=1e-10;
load bravo94
cv=bravo94.rcm.cv(:,2:end);
[psi,lambda]=sleptap(length(cv),8);

[~,spp,snn,spn]=mspec(cv,psi,lambda,'adaptive');
[~,spp2,snn2,spn2]=mspec(conj(cv)',psi,2,lambda,'adaptive');
bool=aresame(spp,spp2',1e-9)&&aresame(snn,snn2',1e-9)&&aresame(spn,conj(spn2)',1e-9);
reporttest('MSPEC transposed input for adaptive spectrum',bool)

[~,sxx,syy,sxy]=mspec(real(cv),imag(cv),psi);
[~,spp,snn,spn]=mspec(cv,psi);

[~,spp2,snn2,spn2]=mspec(conj(cv)',psi,2);
bool=aresame(spp,spp2',1e-9)&&aresame(snn,snn2',1e-9)&&aresame(spn,conj(spn2)',1e-9);
reporttest('MSPEC transposed input for average spectrum',bool)

S(1,1,:,:)=sxx;
S(2,2,:,:)=syy;
S(1,2,:,:)=sxy;
S(2,1,:,:)=conj(sxy);

SZ(1,1,:,:)=spp;
SZ(2,2,:,:)=snn;
SZ(1,2,:,:)=spn;
SZ(2,1,:,:)=conj(spn);

T=sqrt(2)*vrep(vrep(tmat,size(S,3),3),size(S,4),4);

SZ2=matmult(matmult(T,S,1),conj(permute(T,[2 1 3 4])),1);

reporttest('MSPEC for (x,y) and (z,z^*) are orthogonal transforms with matrix SQRT(2)*T',aresame(SZ,SZ2,1e-10))

tol=1e-10;
x=bravo94.rcm.cv(:,3);
xo=vfilt(x,24,'nonans');
%num=yf2num(bravo.rcm.yearf);
%t=num-yf2num(floor(bravo.rcm.yearf(1)));
t=yearfrac(bravo94.rcm.num);

psi=sleptap(length(x),8);

p0=vsum(abs(real(x)-vmean(real(x),1)).^2,1)./length(x);
[f,sp]=mspec(real(x),psi);
p1=2*frac(1,2*pi)*(vsum(sp,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for real Bravo, unit sample rate',abs(p1-p0)./p0<4/100);

[f,sp]=mspec(t(2)-t(1),real(x),psi);
p1=2*frac(1,2*pi)*(vsum(sp,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for real Bravo, non-unit sample rate',abs(p1-p0)./p0<4/100);

[f,sp]=mspec(t(2)-t(1),real(x),[]);
p1=2*frac(1,2*pi)*(vsum(sp,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for real Bravo, non-unit sample rate periodogram',abs(p1-p0)./p0<4/100);

[f,sp]=mspec(t(2)-t(1),real(x),psi,'cyc');
p1=2*(vsum(sp,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for real Bravo, non-unit sample rate, cyclic frequency',abs(p1-p0)./p0<4/100);

p0=vsum(abs(x-vmean(x,1)).^2,1)./length(x);
[f,sp,sn]=mspec(x,psi);
p1=frac(1,2*pi)*(vsum(sp,1)+vsum(sn,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for complex Bravo, unit sample rate',abs(p1-p0)./p0<4/100);

[f,sp,sn]=mspec(t(2)-t(1),x,psi);
p1=frac(1,2*pi)*(vsum(sp,1)+vsum(sn,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for complex Bravo, non-unit sample rate',abs(p1-p0)./p0<4/100);

[f,sx,sy]=mspec(real(x),imag(x),psi);
p1=2*frac(1,2*pi)*(vsum(sx+sy,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for bivariate Bravo, unit sample rate',abs(p1-p0)./p0<4/100);

[f,sx,sy]=mspec(t(2)-t(1),real(x),imag(x),psi);
p1=2*frac(1,2*pi)*(vsum(sx+sy,1)).*(f(2)-f(1));
reporttest('MSPEC satisfies Parseval''s theorem to within 4% for bivariate Bravo, non-unit sample rate',abs(p1-p0)./p0<4/100);

reporttest('MSPEC 2*(SXX+SYY)=SPP+SNN for Bravo, non-unit sample rate',aresame(2*sx+2*sy,sp+sn,tol));

N=1001;
cv=randn(N,1)+sqrt(-1)*randn(N,1);
cv=cv-vmean(cv,1);
[~,spp,snn]=mspec(cv,[]);
reporttest('MSPEC periodogram matches variance exactly for zero-mean complex signal and odd length', aresame(sum(spp)./N+sum(snn(2:end))./N,sum(abs(cv).^2)./N,1e-10))

N=1001;
cv=randn(N,1)+sqrt(-1)*randn(N,1)+17;
[~,spp,snn]=mspec(cv,[],'nodemean');
reporttest('MSPEC periodogram matches variance exactly for non-zero-mean complex signal and odd length', aresame(sum(spp)./N+sum(snn(2:end))./N,sum(abs(cv).^2)./N,1e-10))

N=1000;
cv=randn(N,1)+sqrt(-1)*randn(N,1);
cv=cv-vmean(cv,1);
[~,spp,snn]=mspec(cv,[]);
reporttest('MSPEC periodogram matches variance exactly for zero-mean complex signal and even length', aresame(sum(spp)./N+sum(snn(2:end-1))./N,sum(abs(cv).^2)./N,1e-10))

N=1000;
cv=randn(N,1)+sqrt(-1)*randn(N,1)+17;
[~,spp,snn]=mspec(cv,[],'nodemean');
reporttest('MSPEC periodogram matches variance exactly for non-zero-mean complex signal and even length', aresame(sum(spp)./N+sum(snn(2:end-1))./N,sum(abs(cv).^2)./N,1e-10))

N=1000;
dt=3600;
cv=randn(N,1)+sqrt(-1)*randn(N,1)+17;
[~,spp,snn]=mspec(dt,cv,[],'nodemean');
reporttest('MSPEC periodogram matches variance exactly for non-zero-mean complex signal and even length, and non-unit sample rate', aresame(sum(spp)./(dt*N)+sum(snn(2:end-1))./(dt*N),sum(abs(cv).^2)./N,1e-10))

N=1000;
cv=randn(N,1)+sqrt(-1)*randn(N,1)+17;
tic;[~,spp,snn]=mspec(cv,[],'nodemean');toc
tic;S=squared(fft(cv));toc;
b1=aresame(S(1:length(spp)),spp*N,1e-6);
b2=aresame(flipud(S(length(spp)+1:end)),snn(2:end-1)*N,1e-6);

reporttest('MSPEC recovers periodogram for even N', b1&&b2)

N=1001;
cv=randn(N,1)+sqrt(-1)*randn(N,1)+17;
tic;[~,spp,snn]=mspec(cv,[],'nodemean');toc
tic;S=squared(fft(cv));toc;
b1=aresame(S(1:length(spp)),spp*N,1e-6);
b2=aresame(flipud(S(length(spp)+1:end)),snn(2:end)*N,1e-6);
reporttest('MSPEC recovers periodogram for odd N', b1&&b2)

mspec_test_frequency;

function[]=mspec_test_frequency

T=10000;
fo=100./T;
x=cos([1:T]'*2*pi*fo);
psi=sleptap(T,1,1);
[f,sp,sn]=mspec(x+sqrt(-1)./1e10,psi);

[~,jp]=max(sp);
[~,jn]=max(sn);
bool=aresame(frac(1,2*pi)*f(jp),fo,1e-12)&&aresame(frac(1,2*pi)*f(jn),fo,1e-12);
reporttest('MSPEC frequency matches expected exactly, even number of points',bool);

T=10000-1;
fo=100/T;
x=cos([1:T]'*2*pi*fo);
psi=sleptap(T,1,1);
[f,sp,sn]=mspec(x+sqrt(-1)./1e10,psi);

[~,jp]=max(sp);
[~,jn]=max(sn);
bool=aresame(frac(1,2*pi)*f(jp),fo,1e-12)&&aresame(frac(1,2*pi)*f(jn),fo,1e-12);
reporttest('MSPEC frequency matches expected exactly, odd number of points',bool);

% [f,Cuv]=mspec(real(cv),imag(x),psi);
% [f,Suu]=mspec(real(cv),psi);
% [f,Svv]=mspec(imag(cv),psi);
%
% gammauv=Cuv;
% for i=1:size(Suu,2)
%   gammauv(:,i)=Cuv(:,i)./sqrt(Suu(:,i).*Svv(:,3));
% end
% figure,
%
%
% plot(f,abs(gammauv)),xlog,yoffset 1
% title('Cross-spectrum of u(t) at each depth vs. u(t) at #3')

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [y]=bivnor(a,b,rho)
%function [y]=bivnor(a,b,rho)
%
% This function give an approximation of
% cumulative bivariate normal probabilities
%
%                       1                    x^2 - 2*r*x*y + y^2
%    f(x,y,r)= ---------------------- exp(- ----------------------)
%               2*pi*(1-r^2)^(1/2)              2 - 2r^2
%
%    bivnor(a,b,ro)=Int(a..infinity) Int(b..infinity) f(x,y,r) dx dy
%
% Note : Mex compile the file bivnor.c before using it.
%
% Based on Fortran code in Commun. ACM oct 1973 p638 Algo 462
% Translated to C by Ajay Shah (ajayshah@usc.edu)
% Sligtly modified for Matlab compatibility by Moranvil william (moranviw@onid.orst.edu)
%
% 2004 William Moranvil (moranviw@onid.orst.edu)
%

author: Sable

repo: mcbench-benchmarks


text:
clear all;clc

a = [1 0 1 1 0 1 0 0 1 0 1];%---Bit stream
A = 1;  %---Amplitude
Tb = 10e-6; %---Bit time period

%---Get Unipolar NRZ
[x t] = LineEncoder('uninrz',a,Tb,A);
figure(1);
subplot(5,1,1);plot(t,x,'LineWidth',2);grid on;
title(['Unipolar NRZ line coding for ', num2str(a)]);
xlabel('Time in sec');
ylabel('Amplitude');
axis([0,max(t),min(x)-A,max(x)+A]);

%---Get Unipolar RZ
[x t] = LineEncoder('unirz',a,Tb,A);
subplot(5,1,2);plot(t,x,'LineWidth',2);grid on;
title(['Unipolar RZ line coding for ', num2str(a)]);
xlabel('Time in sec');
ylabel('Amplitude');
axis([0,max(t),min(x)-A,max(x)+A]);

%---Get Polar RZ
[x t] = LineEncoder('polrz',a,Tb,A);
subplot(5,1,3);plot(t,x,'LineWidth',2);grid on;
title(['Polar RZ line coding for ', num2str(a)]);
xlabel('Time in sec');
ylabel('Amplitude');
axis([0,max(t),min(x)-A,max(x)+A]);

%---Get Polar NRZ
[x t] = LineEncoder('polnrz',a,Tb,A);
subplot(5,1,4);plot(t,x,'LineWidth',2);grid on;
title(['Polar NRZ line coding for ', num2str(a)]);
xlabel('Time in sec');
ylabel('Amplitude');
axis([0,max(t),min(x)-A,max(x)+A]);

%---Get Manchester
[x t] = LineEncoder('manchester',a,Tb,A);
subplot(5,1,5);plot(t,x,'LineWidth',2);grid on;
title(['Manchester line coding for ', num2str(a)]);
xlabel('Time in sec');
ylabel('Amplitude');
axis([0,max(t),min(x)-A,max(x)+A]);


text:
close all;
clear all;
clc;

data_file_path = 'bin/mmv_phase_transition_noiseless_s_8.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'mmv_noiseless_s_8';
options.subtitle = 'Noiseless, S = 8';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function C = encode(p,R)
%
% Compute vector of transmit costs for each stage of the KD-tree
%

N = getNpts(p);  bw = getBW(p,1);
if (size(p.bandwidth,2) > 2*N)
  error('Encoding of variable bandwidths not yet supported...');
end;
if (any( abs(getWeights(p) - 1/N) > 2*eps ))
  error('Encoding of variable weights not yet supported...');
end;
C = zeros(1,2*N); C(1) = 2*getDim(p)*R;		% Root node...

for i=1:N-1,            % calc costs of splitting each node; 1/2 cost to each side
 if (~isLeaf(i,N))
  mu0 = p.means(:,i);				% get means of parent, children
  mu1 = p.means(:, double(p.leftch(i))+1);
  mu2 = p.means(:, double(p.rightch(i))+1);
  sig0= sqrt( p.bandwidth(:,i) );		% and bw's of parent, children
  sig1= sqrt( p.bandwidth(:,double(p.leftch(i))+1) );
  sig2= sqrt( p.bandwidth(:,double(p.rightch(i))+1) );
  sig0UB = sqrt(max( sig0.^2 - bw.^2, 2^-R));	% don't round down too far...

  [tmp,splitDim] = max( sig0UB );		% which dimension are we splitting on...
  %[tmp,splitDim2] = max( abs(mu0-mu1) + abs(mu0-mu2) );
  %if (splitDim ~= splitDim2) warning('Disagreement?'); end;

  % COMPUTE COST OF SENDING MEANS
  muMax = max(mu1,mu2); sigDiff = sig0.^2 - (mu1.^2 + mu2.^2 - mu0.^2);
  for j=splitDim,
    costMu = gauss( muMax(j), mu0(j), sig0UB(j).^2 ,1,R );
  end;
  for j = [1:splitDim-1,splitDim+1:getDim(p)];
    costMu = costMu + gauss2( mu1(j), mu0(j), 1*sig0UB(j).^2 ,1,R );
  end;

  % COMPUTE COST OF SENDING VARIANCES
  if (isLeaf(p.leftch(i),N) && isLeaf(p.rightch(i),N)), costSig = 0; costMu = 0;
  elseif (isLeaf(p.leftch(i),N) || isLeaf(p.rightch(i),N)), costSig = 0;
  else
    for j=splitDim,
      costSig = gauss2( sig1(j).^2 , sig0(j).^2/2, sig0(j).^2/4, 1, R);
    end;
    for j = [1:splitDim-1,splitDim+1:getDim(p)]
      costSig = costSig + gauss2( sig1(j).^2 , sig0(j).^2, sig0(j).^2/2, 1, R);
    end;
  end;

  %[costMu,costSig]
  C(double(p.leftch(i))+1) = .5 * (C(i) + costMu + costSig);
  C(double(p.rightch(i))+1) = .5 * (C(i) + costMu + costSig);
 end;
end;

function v = gauss(x,mu,sig2,sigOut2,R)		% SIMPLE 1-SIDED GAUSSIAN
  v = R-log2(  2*1/sqrt(2*pi*sig2) .* exp(-(x-mu).^2 ./(2*sig2)) );

%function v = gauss(x,mu,sig2,sigOut2,R)	% 1-SIDED W/ OUTLIER PROCESS
%  v = R-log2(  2*1/sqrt(2*pi*sig2) .* exp(-(x-mu).^2 ./(2*sig2)) ); + ...
%               1/sqrt(2*pi*sig2) .* exp(-(x-mu).^2 ./(2*sigOut2)) );

function v = gauss2(x,mu,sig2,sigOut2,R)	% SIMPLE 2-SIDED GAUSSIAN
  v = R-log2(  1/sqrt(2*pi*sig2) .* exp(-(x-mu).^2 ./(2*sig2)) );

%function v = gauss2(x,mu,sig2,sigOut2,R)	% 2-SIDED W/ OUTLIER PROCESS
%  v = R-log2(  1/sqrt(2*pi*sig2) .* exp(-(x-mu).^2 ./(2*sig2))  + ...
%               1/sqrt(2*pi*sig2) .* exp(-(x-mu).^2 ./(2*sigOut2)) );

function b = isLeaf(ind,N)
  ind = double(ind);
  if (ind <= 0 || ind > 2*N) b = 0;
  elseif (ind <= N-1) b = 0;
  else b = 1;
  end;

author: ShapeNet

repo: RenderForCNN


text:
function test_bug1911

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_databrowser

% When ft_movieplotER is called within ft_databrowser by right-clicking on a
% segment of data and selecting ft_movieplotER the movieplot is opened in the
% ft_databrowser window. Only a quarter of the topoplot is visible in this case.
%
% I have only experienced this problem in Windows.

cd(dccnpath('/home/common/matlab/fieldtrip/data/test'));
load bug1911.mat

cfg = [];
cfg.layout = 'easycapM10.lay';
ft_databrowser(cfg,timelock);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function img=headPhantom(varargin)
%HEADXPHANTOM returns the head phantom
%
%  IMG=HEADXPHANTOM() returns 128^3 image IMG
%
%  IMG=HEADXPHANTOM(SZ) returns a SZ^3 image IMG if SZ is scalar, or a [SZ(1)
%   SZ(2) SZ(3] image, of SZ is a vector.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%                     and
%                     https://www.mathworks.com/matlabcentral/fileexchange/view_license?file_info_id=35548
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Kyungsang Kim, modified by Ander Biguri
%--------------------------------------------------------------------------

% Deal with inputs
if nargin==0
    sz=[128,128,128];
end
if nargin==1
    if ~isnumeric(varargin{1})
        error('TIGRE:headphantom:invalid input','input is not numerical');
    end
    nsz=max(size(varargin{1}));
    if nsz==2 || nsz>3
        error('TIGRE:headphantom:invalid input','input is not 1x1 or 1x3');
    end
    if nsz==1
        sz=[varargin{1},varargin{1},varargin{1}];
    else
        sz=varargin{1};
    end
end
% load data
curr_path=mfilename('fullpath');
data_path=curr_path(1:end-length('/MATLAB/Test_data/MRheadbrain/headPhantom'));
data_path=[data_path '/Common/data/'];
data=load([data_path 'head.mat']);
img=data.img;

% interpolate data to get desired size
[y, x, z]=...
   ndgrid(linspace(1,size(img,1),sz(1)),...
          linspace(1,size(img,2),sz(2)),...
          linspace(1,size(img,3),sz(3)));

imOut=interp3(img,x,y,z,'nearest');
% out!
img=imOut;

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function [Operations,MasterOperations] = TS_LinkOperationsWithMasters(Operations,MasterOperations)
% TS_LinkOperationsWithMasters    Link Operations with MasterOperations using Label field

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

numOps = height(Operations);
numMops = height(MasterOperations);

% ------------------------------------------------------------------------------
% Match operations to a master ID
% ------------------------------------------------------------------------------
for i = 1:numOps
    theMasterMatch = strcmp(Operations.Label{i},MasterOperations.Label);
    if sum(theMasterMatch)==0
        error('No master match for operation: %s',Operations.Name{i});
    end
    Operations.MasterID(i) = MasterOperations.ID(theMasterMatch);
end

% No longer need the label field from the Operations table
Operations(:,'Label') = [];

%-------------------------------------------------------------------------------
% Check that all master operations are required
%-------------------------------------------------------------------------------
mastersNeeded = ismember(MasterOperations.ID,Operations.MasterID);
if ~all(mastersNeeded)
    warning(['%u/%u master operations are not used by the %u operations' ...
                         ' and will be removed.'],...
                        sum(~mastersNeeded),numMops,numOps);
    MasterOperations = MasterOperations(mastersNeeded,:);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% GETMSWTFEAT Gets the Multiscale Wavelet Transform features, these
% include: Energy, Variance, Standard Deviation, and Waveform Length
% feat = getmswtfeat(x,winsize,wininc,SF)
% ------------------------------------------------------------------
% The signals in x are divided into multiple windows of size
% "winsize" and the windows are spaced "wininc" apart.
% Inputs
% ------
%    x: 		columns of signals
%    winsize:	window size (length of x)
%    wininc:	spacing of the windows (winsize)
%    SF:        sampling frequency (Not used in the current implementation, but I left you some options down there)
% Outputs
% -------
%    feat:     WT features organized as [Energy, Variance, Waveform Length, Entropy]

% Example
% -------
% feat = getmswtfeat(rand(1024,1),128,32,32)
% Assuming here rand(1024,1) (this can be any one dimensional signal,
% for example EEG or EMG) is a one dimensional signal sampled at 32
% for 32 seconds only. Utilizing a window size of 128 at 32 increments,
% features are extracted from the wavelet tree.
% I assumed 10 decomposition levels (J=10) below in the code.
% For a full tree at 10 levels you should get 11 features
% as we have decided to extract 4 types of features then we get 11 x 4 =44
% features.
% =========================================================================
% Multiscale Wavelet Transform feature extraction code by Dr. Rami Khushaba
% Research Fellow - Faculty of Engineering and IT
% University of Technology, Sydney.
% Email: Rami.Khushaba@uts.edu.au
% URL: www.rami-khushaba.com (Matlab Code Section)
% last modified 29/08/2012
% last modified 09/02/2013

function feat = getmswtfeat(x,winsize,wininc,SF)

if nargin < 4
    if nargin < 3
        if nargin < 2
            winsize = size(x,1);
        end
        wininc = winsize;
    end
    error('Please provide the sampling frequency of this signal')
end

datawin = ones(winsize,1);
datasize = size(x,1);
Nsignals = size(x,2);

%% Chop the signal according to a sliding window approach
numwin = floor((datasize - winsize)/wininc)+1;
% allocate memory
feat = zeros(winsize,numwin);
st = 1;
en = winsize;
for i = 1:numwin
    curwin = x(st:en,:).*repmat(datawin,1,Nsignals);
    feat(1:winsize,i) = detrend(curwin);

    st = st + wininc;
    en = en + wininc;
end
%% ---------------- Various options for J -----------------------
% Note I put SF above in the inputs because you can use SF to determine the
% best decompisition level J, however for simplicity here I put it J=10;
J=10;% Number of decomposition levels which can also be set using
% or J=wmaxlev(winsize,'Sym5');
% or J=(log(SF/2)/log(2))-1;
%% Multisignal one-dimensional wavelet transform decomposition
dec = mdwtdec('col',feat,J,'db4');
% Proceed with Multisignal 1-D decomposition energy distribution

if isequal(dec.dirDec,'c')
    dim = 1;
end
[cfs,longs] = wdec2cl(dec,'all');
level = length(longs)-2;

if dim==1
    cfs = cfs';
    longs = longs';
end
numOfSIGs  = size(cfs,1);
num_CFS_TOT = size(cfs,2);
absCFS   = abs(cfs);
absCFS0   = (cfs);
cfs_POW2 = absCFS.^2;
Energy  = sum(cfs_POW2,2);
percentENER = 0*ones(size(cfs_POW2));
notZER = (Energy>0);
percentENER(notZER,:) = 100*cfs_POW2(notZER,:)./Energy(notZER,ones(1,num_CFS_TOT));

%% or try this version below and tell us which  one is the best on your data
% percentENER(notZER,:) = cfs_POW2(notZER,:);

%% Pre-define and allocate memory
tab_ENER = zeros(numOfSIGs,level+1);
tab_VAR = zeros(numOfSIGs,level+1);
% tab_STD = zeros(numOfSIGs,level+1);
tab_WL = zeros(numOfSIGs,level+1);
tab_entropy = zeros(numOfSIGs,level+1);

%% Feature extraction section
st = 1;
for k=1:level+1
    nbCFS = longs(k);
    en  = st+nbCFS-1;
    tab_ENER(:,k) = mean(percentENER(:,st:en),2);%.*sum(abs(diff(absCFS0(:,st:en)')'),2); % energy per waveform length
    tab_VAR(:,k) = var(percentENER(:,st:en),0,2); % variance of coefficients
    %     tab_STD(:,k) = std(percentENER(:,st:en),[],2); % standard deviation of coefficients
    tab_WL(:,k) = sum(abs(diff(percentENER(:,st:en)').^2))'; % waveform length
    percentENER(:,st:en) = percentENER(:,st:en)./repmat(sum(percentENER(:,st:en),2),1,size(percentENER(:,st:en),2));
    tab_entropy(:,k) = -sum(percentENER(:,st:en).*log(percentENER(:,st:en)),2)./size(percentENER(:,st:en),2);
    st = en + 1;
end
feat =([log1p([tab_ENER tab_VAR  tab_WL]) tab_entropy]);


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

X = [ones(m, 1) X];
temp = sigmoid(X * Theta1');
temp = [ones(m, 1) temp];
r = sigmoid(temp * Theta2');
for i = 1:m
    [~, temp] = max(r(i, :));
    p(i) = temp;
end

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

end

repo: machine-learning


text:
%% overfitting #1
data = importdata('BCPF_fiber_rmselb_m30_r5.csv');
rmse = data(1:10,:);
lb = data(11:20,:);

color = [55,126,184]/255;
gray = [211,211,211]/255;
orange = [ 0.91,0.41,0.17];
blue = [0,0.4470,0.7410];
light_blue = [0.5843,0.8157,0.9882];
light_red =  [255,204,204]/255;

figure('DefaultAxesFontSize',16);
ax = gca;

yyaxis left;
xbound = [1,1,1:200,200,200,fliplr(1:200)];
ybound = [max(rmse(:,1)),min(rmse(:,1)),min(rmse),min(rmse(:,end)),max(rmse(:,end)),fliplr(max(rmse))];
hold on;
s = fill(xbound,ybound,color,'EdgeColor',gray,'linewidth',0.5);
alpha(s,0.10);
hold on;
plot(rmse','-','color',light_blue,'linewidth',0.5);
hold on;
h1 = plot(mean(rmse),'-','color',blue,'linewidth',2,'DisplayName','rmse');

box on;
grid on;
ax.XTick = [0,25,50,100,200];
ax.XTickLabel = [0,25,50,100,200];
ylim([4.60,5.60]);
ax.YTick = [4.60,4.70,4.80,4.90,5.10,5.60];
ax.YTickLabel = ['4.60';'4.70';'4.80';'4.90';'5.10';'5.60'];
ax.GridLineStyle = '-.';
xlabel('Epoch');
ylabel('RMSE (km/h)');

yyaxis right;
xbound = [1,1,1:200,200,200,fliplr(1:200)];
ybound = [max(lb(:,1)),min(lb(:,1)),min(lb),min(lb(:,end)),max(lb(:,end)),fliplr(max(lb))];
hold on;
s = fill(xbound,ybound,gray,'EdgeColor',gray,'linewidth',0.5);
alpha(s,0.10);
hold on;
plot(lb','-','color',light_red,'linewidth',0.5);
hold on;
h2 = plot(mean(lb),'-','color',orange,'linewidth',2,'DisplayName','lb');
ylim([-1.24e7,-1.19e7]);
ylabel('Lower bound');

legend([h1,h2],{'Average of RMSE';'Average of lower bound'},'location','north');

set(gcf, 'PaperSize', [6 4.5]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 6 4.5]);
saveas(gcf,'overfitting_ms30_r5','pdf');

%% overfitting #2
data = importdata('BCPF_fiber_rmselb_m30_r10.csv');
rmse = data(1:10,:);
lb = data(11:20,:);

figure('DefaultAxesFontSize',16);
ax = gca;

yyaxis left;
xbound = [1,1,1:200,200,200,fliplr(1:200)];
ybound = [max(rmse(:,1)),min(rmse(:,1)),min(rmse),min(rmse(:,end)),max(rmse(:,end)),fliplr(max(rmse))];
hold on;
s = fill(xbound,ybound,color,'EdgeColor',gray,'linewidth',0.5);
alpha(s,0.10);
hold on;
plot(rmse','-','color',light_blue,'linewidth',0.5);
hold on;
h1 = plot(mean(rmse),'-','color',blue,'linewidth',2,'DisplayName','rmse');

box on;
grid on;
ax.XTick = [0,25,50,100,200];
ax.XTickLabel = [0,25,50,100,200];
ylim([4.35,8.60]);
ax.YTick = [4.35,4.60,4.90,5.20,6.00,8.60];
ax.YTickLabel = ['4.35';'4.60';'4.90';'5.20';'6.00';'8.60'];
ax.GridLineStyle = '-.';
xlabel('Epoch');
ylabel('RMSE (km/h)');

yyaxis right;
xbound = [1,1,1:200,200,200,fliplr(1:200)];
ybound = [max(lb(:,1)),min(lb(:,1)),min(lb),min(lb(:,end)),max(lb(:,end)),fliplr(max(lb))];
hold on;
s = fill(xbound,ybound,gray,'EdgeColor',gray,'linewidth',0.5);
alpha(s,0.10);
hold on;
plot(lb','-','color',light_red,'linewidth',0.5);
hold on;
h2 = plot(mean(lb),'-','color',orange,'linewidth',2,'DisplayName','lb');
ylim([-1.205e7,-1.19e7]);
ylabel('Lower bound');

legend([h1,h2],{'Average of RMSE';'Average of lower bound'},'location','north');

set(gcf, 'PaperSize', [6 4.5]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 6 4.5]);
saveas(gcf,'overfitting_ms30_r10','pdf');

author: xinychen

repo: academic-drawing


text:
function out = SY_Trend(y)
% SY_Trend  Quantifies various measures of trend in a time series.
%
%---INPUT:
% y, the input time series.
%
%---OUTPUTS:
% Linearly detrends the time series using detrend, and returns the ratio of
% standard deviations before and after the linear detrending. If a strong linear
% trend is present in the time series, this operation should output a low value.
%
% Also fits a line and gives parameters from that fit, as well as statistics on
% a cumulative sum of the time series.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

if ~BF_iszscored(y)
    warning('The input time series should be z-scored')
end

N = length(y);

% Ratio of std before and after linear detrending:
out.stdRatio = std(detrend(y)) / std(y);

% Linear fit:
[out.gradient,out.intercept] = LinearFit(1:N,y);

% Stats on the cumulative sum:
yC = cumsum(y);
out.meanYC = mean(yC);
out.stdYC = std(yC);
[out.gradientYC,out.interceptYC] = LinearFit(1:N,yC);

% Mean cumsum in first and second half of the time series:
out.meanYC12 = mean(yC(1:floor(N/2)));
out.meanYC22 = mean(yC(floor(N/2)+1:end));

% ------------------------------------------------------------------------------
function [m,b] = LinearFit(xData,yData)
    if size(xData,1) ~= size(yData,1);
        yData = yData';
    end
    coeff = polyfit(xData,yData,1);
    m = coeff(1); b = coeff(2);
end
% ------------------------------------------------------------------------------

end


text:
function [vtx,elt,col,data] = mshReadMsh(filename)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshReadMsh.m                                  |
%|    #    |   VERSION    : 0.42                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 31.12.2018                                    |
%| ( === ) |   SYNOPSIS   : Read .msh files (particle, edge, triangular   |
%|  `---'  |                and tetrahedral                               |
%+========================================================================+

% Open
fid = fopen(filename,'r');
if (fid==-1)
    error('mshReadMsh.m : cant open the file');
end

% Read file using keywords
str = fgets(fid);
while ~(str==-1)
    if contains(str,'$Nodes')
        vtx = readNodes(fid);
    elseif contains(str,'$Elements')
        [elt,col] = readElements(fid);
    elseif contains(str,'$NodeData')
        data = readNodesData(fid,size(vtx,1));
    end
    str = fgets(fid);
end

% Close file
fclose(fid);

% Only keep higher elements (tetra > triangle > edge > particles)
ord = sum(elt>0,2);
dim = max(ord);
elt = elt(ord==dim,1:dim);
col = col(ord==dim);
end

function vtx = readNodes(fid)
% Nodes
Nvtx = str2double(fgets(fid));
vtx  = zeros(Nvtx,3);
for i = 1:Nvtx
    tmp      = str2num(fgets(fid));
    vtx(i,:) = tmp(2:4);
end

% Verify nodes ending
str = fgets(fid);
if ~contains(str,'$EndNodes')
    error('mshReadMsh.m : unavailable case');
end
end

function [elt,col] = readElements(fid)
% Initialize
Nelt = str2double(fgets(fid));
elt  = zeros(Nelt,4);
col  = zeros(Nelt,1);

% Elements (up to tetra)
for i = 1:Nelt
    tmp    = str2num(fgets(fid));
    col(i) = tmp(4);
    if (tmp(2) == 15)    % particles
        elt(i,1)   = tmp(6);
    elseif (tmp(2) == 1) % segment
        elt(i,1:2) = tmp(6:7);
    elseif (tmp(2) == 2) % triangles
        elt(i,1:3) = tmp(6:8);
    elseif (tmp(2) == 4) % tetrahedra
        elt(i,1:4) = tmp(6:9);
    end
end

% Verify elements ending
str = fgets(fid);
if ~contains(str,'$EndElements')
    error('mshReadMsh.m : unavailable case');
end
end

function data = readNodesData(fid,Nvtx)
% Go to Node Data
str = fgets(fid);
while ~contains(str,num2str(Nvtx))
    str = fgets(fid);
end

% Initialization with first row
tmp       = str2num(fgets(fid));
data      = zeros(Nvtx,length(tmp)-1);
data(1,:) = tmp(2:end);

% Nodes data
for i = 2:Nvtx
    tmp       = str2num(fgets(fid));
    data(i,:) = tmp(2:end);
end

% Verify node data ending
str = fgets(fid);
if ~contains(str,'$EndNodeData')
    error('mshReadMsh.m : unavailable case');
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function img = atlasMapRgn2FullImage(img,imgRegion,X,Y,fillvalue)
%
%   img = atlasMapRgn2FullImage(img,imgRegion,X,Y,[fillvalue = -1])
%
% Author: Wandell
% Purpose:
%    We use atlasCreate (via CP2TFORM and IMTRANFORM) to transform the standard
%    atlas into the coordinate frame of the data. This produces a set of
%    standard atlas values that fall within an imgRegion inside the full
%    image of the data (img). The coordinates of imgRegion are at the
%    locations between X(1):X(2) and Y(1):Y(2) in the destination image,
%    img.  Because of the limitations on the affine transformations,
%    sometimes the data in imgRegion are   invalid,  having a fill value
%    (default of -1).  This routine copies the valid data in imgRegion into
%    the full image.
%
% Examples:
%    atlasE = atlasMapRgn2FullImage(atlasE,atlasET,X,Y);
%

if ~exist('fillvalue','var') , fillvalue = -1; end

xCoords = round(X(1):X(2));
yCoords = round(Y(1):Y(2));
for ii=1:length(xCoords)
    for jj=1:length(yCoords)
        if imgRegion(jj,ii) ~= fillvalue
            img(yCoords(jj),xCoords(ii)) = imgRegion(jj,ii);
        end
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function result = analyzeHardStop(debug, ground_truth, tFailure_s, PerformanceSpec)

% Author: Alexander Wischnewski     Date: 24-02-2019
%
% Description:
%   scans the debug data and ground truth whether a hard stop
%   was performed within the given specifications
%
% Input:
%   debug:          Debug structure of the vehicle
%   ground_truth:   Physical ground truth data from simulation
%   tFailure_s:     Timestamp where the failure occured
%   PerformanceSpec Structure with fields:
%       AccxMin_mps2        Minimum mean negative acceleration applied until vehicle
%                           stops
%       tDetectionMax_s     Maximum time to detect the hard stop
%
% Outputs:
%   result:       1 if a valid emergency stop was detected, 0 otherwise

result = 1;

disp(' ');
disp('------------------------------');
disp('Check Hard Stop: ');
disp('------------------------------');
% find timestamp where vehicle has gone to hard stop
idxDetect = find(uint16(debug.debug_mloc_statemachine_debug_TUMVehicleState.Data) == 60, 1);
if isempty(idxDetect)
    disp('The vehicle did not detect the hard stop');
    result = 0;
    return
else
    tFailuredDetected_s = debug.debug_mloc_statemachine_debug_TUMVehicleState.Time(idxDetect);
    disp(['The vehicle took ' num2str(tFailuredDetected_s - tFailure_s) ' s to detect the failure.']);
end
% check if fault detection time was within spec
if(tFailuredDetected_s <= 0 && tFailuredDetected_s >= PerformanceSpec.tDetectionMax_s)
    disp('The vehicle did not detect the emergency within the specified time');
    result = 0;
    return
end
% find point where failure was detected
idxDetect_gt = find(ground_truth.SimRealState_ax_mps2.Time > tFailuredDetected_s, 1);
% find timestamp where velocity has gone to new stopped
% for very low velocites, the stop controller does not behave precisly in
% sim, therefore its tested against roughly 5kph.
idxStop = find(ground_truth.SimRealState_vx_mps.Data(idxDetect_gt:end) < 2, 1) + idxDetect_gt;
if isempty(idxStop)
    disp('The vehicle did not stop.');
    result = 0;
    return
else
    tStop_s = ground_truth.SimRealState_vx_mps.Time(idxStop);
    disp(['The vehicle took ' num2str(tStop_s - tFailure_s) ' s to stop.']);
end
% check if control performance was ok
meanAccX = mean(ground_truth.SimRealState_ax_mps2.Data(idxDetect_gt:idxStop));
disp(['The mean longitudinal accleeration was: ' num2str(meanAccX)]);
if(meanAccX > PerformanceSpec.AccxMin_mps2)
    disp('The vehicle did not stay within the specified deceleration bounds.');
    result = 0;
    return
end

disp(' ');
disp('All checks passed.');

author: TUMFTM


text:
%% autoEncoder class
%
% Copyright (c) 2018-present, Mahmoud Afifi
% York University, Canada
% mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
% All rights reserved.
%
% Please cite the following work if this program is used:
% Mahmoud Afifi, Brian Price, Scott Cohen, and Michael S. Brown,
% "When color constancy goes wrong: Correcting improperly white-balanced
% images", CVPR 2019.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

classdef autoEnc
    properties
        weights
        bias
    end
    methods
        function feature = encode(obj,hist)
            feature = (obj.sigmf(obj.weights * reshape(hist,[],1) + obj.bias,[1,0]))';
        end
        function y = sigmf(obj, x, params)
            a = cast(params(1),'like',x);
            c = cast(params(2),'like',x);
            y = 1./(1 + exp(-a*(x-c)));

        end
    end
end

author: mahmoudnafifi


text:
function ic = ic(icfh, x_struct,p_struct,x,p, userdata)
% initial condition function
icHandler = ocl.Constraint(userdata);
x = ocl.Variable.create(x_struct,x);
p = ocl.Variable.create(p_struct,p);
icfh(icHandler,x,p)
ic = icHandler.values;
assert(all(icHandler.lowerBounds==0) && all(icHandler.upperBounds==0),...
  'In initial condition are only equality constraints allowed.');

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [testFileList]=getTestFiles(testSet)

% function [testFileList]=getTestFiles(testSet)
% ------------------------------------------------------------------------
% This function obtains the list of test files from the documentation
% folder. The test files are either the HELP_* files (testSet='help'), the
% DEMO_* files (testSet='demo') or all (testSet='all').
%
% Change log:
% 2021/08/05 Created
% ------------------------------------------------------------------------

%%
testFileList=fullfile(fileparts(fileparts(mfilename('fullpath'))),'docs');

%Get list of M-files in documentation folder
allFiles_publish = dir(fullfile(testFileList,'*.m'));
allFiles_publish={allFiles_publish(1:end).name};
allFiles_publish=sort(allFiles_publish(:));

logicExclude=~gcontains(allFiles_publish,'obj_DEMO_'); %Skip objective functions
logicDemo=gcontains(allFiles_publish,'DEMO_'); %Logic for demo files
logicHelp=gcontains(allFiles_publish,'HELP_'); %Logic for help files

%Get list of files, either all or just help or demo files
switch lower(testSet) %Lower so both HELP and help etc work
    case 'all'
        testFileList=allFiles_publish((logicHelp | logicDemo) & logicExclude);
    case 'demo'
        testFileList=allFiles_publish(logicDemo & logicExclude);
    case 'help'
        testFileList=allFiles_publish(logicHelp & logicExclude);
    otherwise
        error([testMode,' is not a valid option. Use: all, demo, or help']);
end

end
%%
% _*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:
function trans = createScaling3d(varargin)
%CREATESCALING3D Create the 4x4 matrix of a 3D scaling.
%
%   TRANS = createScaling3d(S);
%   returns the scaling transform corresponding to a scaling factor S in
%   each direction. S can be a scalar, or a 1-by-3 vector containing the
%   scaling factor in each direction.
%
%   TRANS = createScaling3d(SX, SY, SZ);
%   returns the scaling transform corresponding to a different scaling
%   factor in each direction.
%
%   The returned matrix has the form :
%   [SX  0  0  0]
%   [ 0 SY  0  0]
%   [ 0  0 SZ  0]
%   [ 0  0  0  0]
%
%   See also
%   transforms3d, transformPoint3d, transformVector3d, createTranslation3d,
%   createRotationOx, createRotationOy, createRotationOz
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2006-04-20
% Copyright 2006-2022 INRA - TPV URPOI - BIA IMASTE

%% default arguments
sx = 1;
sy = 1;
sz = 1;
center = [0 0 0];

%% process input parameters
if nargin == 1
    % only one argument -> scaling factor
    [sx, sy, sz]= parseScalingFactors(varargin{1});

elseif nargin == 2
    % 2 arguments, giving center and uniform scaling
    center = varargin{1};
    [sx, sy, sz]= parseScalingFactors(varargin{2});

elseif nargin == 3
    % 3 arguments, giving scaling in each direction
    sx = varargin{1};
    sy = varargin{2};
    sz = varargin{3};

elseif nargin == 4
    % 4 arguments, giving center and scaling in each direction
    center = varargin{1};
    sx = varargin{2};
    sy = varargin{3};
    sz = varargin{4};
end

%% create the scaling matrix
trans = [...
    sx 0 0 center(1)*(1-sx);...
    0 sy 0 center(2)*(1-sy);...
    0 0 sz center(3)*(1-sz);...
    0 0 0 1];

%% Helper function
function [sx, sy, sz] = parseScalingFactors(var)

if length(var)==1
    % same scaling factor in each direction
    sx = var;
    sy = var;
    sz = var;
elseif length(var)==3
    % scaling is a vector, giving different scaling in each direction
    sx = var(1);
    sy = var(2);
    sz = var(3);
else
    error('wrong size for first parameter of "createScaling3d"');
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [summaryTable, booleanTable, problemTypes, solvers] = getAvailableGAMSSolvers
% This function return the GAMS solvers which are available in your system
% according to the license.
%
% USAGE:
%
%         [summaryTable, booleanTable, problemTypes, solvers] = getAvailableGAMSSolvers
%
% OUTPUTS:
%    summaryTable    Type: cell array
%                    Description: matrix summarizing which problem can be
%                    solved by which solver. A "yes" in position (i,j)
%                    means that the GAMS solver "j" is available in your
%                    system to solve a problem of type "i". Please note
%                    that headers are included in this summaryTable, so the
%                    first column and the first row have labels for solvers
%                    and problem types, respectively.
%    booleanTable    Type: double matrix
%                    Description: matrix summarizing which problem can be
%                    solved by which solver. A "1" in position (i,j) means
%                    that the GAMS solver "j" is available in your system
%                    to solve a problem of type "i"
%    problemTypes    Type: cell array of strings
%                    Description: list of problem types that can be solved
%                    in GAMS. E.g: LP, MIP, etc
%    solvers         Type: cell array of strings
%                    Description: list of solvers available in GAMS
%
% .. Author: - Sebastiian Mendoza, May 30th 2017, Center for Mathematical Modeling, University of Chile, snmendoz@uc.cl

gamsPath = which('gams');  % verify that gams is installed
if isempty(gamsPath)
    error('GAMS is not installed or GAMS path has not been added to MATLAB path.');
end

%verify that licememo.gms in the path of the system
licememoFullPath = which('licememo.gms');
if isempty(licememoFullPath)
    error('licememo.gms in not in MATLAB path.');
end

%copy licememo.gms to the current folder
current = pwd;
if ~strcmp(regexprep(licememoFullPath, {'\\licememo.gms', '/licememo.gms'}, {'', ''}), current)
    copyfile(licememoFullPath);
end

%run licememo.gms
run = system('gams licememo');

%if licememo.gms ran correctly.
if run == 0
    % read the content of the file generated by licememo.gms to find the
    % GAMS solvers available
    fileID = fopen('licememo.txt', 'r');
    summaryTable = cell(2, 2);
    booleanTable = zeros(120, 30);
    tline = fgetl(fileID);
    solvers = {};
    while ischar(tline)
        if strfind(tline,'Solver/Modeltype Matrix:')
            tline = fgetl(fileID); if isempty(tline); tline = fgetl(fileID); end;
            problemTypes = {};
            aux = regexp(tline, ' ', 'split');
            for i = 1:length(aux)
                if ~isempty(aux{i})
                    problemTypes(end + 1) = aux(i);
                end
            end
            nProblemType = length(problemTypes);

            nSolvers = 0;
            tline = fgetl(fileID);
            while ischar(tline)
                solver = regexprep(tline(1:12), ' ', '');
                nSolvers = nSolvers + 1;
                solvers{nSolvers} = solver;
                for i = 1:length(problemTypes)
                    if length(tline) < 12 + 7 * (i - 1)
                        continue;
                    elseif length(tline)>12 + 7 * (i - 1) && length(tline)< 12 + 7 * (i) && ~isempty(regexprep(tline(12 + 7 * (i - 1):end), ' ', ''))
                        summaryTable{nSolvers,i} = 'yes';
                        booleanTable(nSolvers,i) = 1;
                        continue;
                    elseif ~isempty(regexprep(tline(12 + 7 * (i - 1):12 + 7 * i), ' ', ''))
                        summaryTable{nSolvers,i} = 'yes';
                        booleanTable(nSolvers,i) = 1;
                    end
                end
                tline = fgetl(fileID);
            end
            break;
        end
        tline = fgetl(fileID);
    end

    booleanTable = booleanTable(1:nSolvers, 1:nProblemType);
    summaryTable = [{''}, problemTypes; [solvers', summaryTable]];
    fclose(fileID);
    %delete autogenerated file
    delete('licememo.txt')
else

    summaryTable ={};
    booleanTable = 0;
    problemTypes = {};
    solvers = {};
end
%delete autogenerated file and licememo.gms
delete('licememo.gms')
delete('licememo.lst')

end

author: opencobra


text:
function [niiC, scanFileNameC] = scan2nii(planC,scanNumV,tmpDirPath,reorientFlag)

if ~exist('tmpDirPath','var') || isempty(tmpDirPath) || ~exist(tmpDirPath,'dir')
    tmpDirPath = fullfile(getCERRPath, 'ImageRegistration', 'tmpFiles');
end

if ~exist('reorientFlag','var') || isempty(reorientFlag)
    reorientFlag = 1;
end

for i = 1:numel(scanNumV)
    scanNum = scanNumV(i);
    [scanUniqName, ~] = genScanUniqName(planC,scanNum);
    [affineMat,scan3M_RAS,voxel_size] = getPlanCAffineMat(planC, scanNum, 1);
    qOffset = affineMat(1:3,end)';
    scanFileName = fullfile(tmpDirPath, ['scan_' num2str(scanNumV(i)) '_' scanUniqName '.nii']);
    niiC{i} = vol2nii(scan3M_RAS,affineMat,qOffset,voxel_size,scanFileName);
    scanFileNameC{i} = scanFileName;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function test_pull574_reref

% MEM 2gb
% WALLTIME 01:00:00
% DEPENDENCY ft_preprocessing preproc ft_preproc_reference

% test the ft_preprocessing for re-referencing to rest
% Li Dong $ 2020.3.10, (Lidong@uestc.edu.cn) and J.M.Schoffelen

load(dccnpath('/home/common/matlab/fieldtrip/data/test/pull574.mat')); % load data and leadfield calculated by FieldTrip

% -------------------------------------------------------------------------
% A snippet of code to calculate leadfield using FieldTrip
% cfg                        = [];
% cfg.elec                   = elec_aligned; % aligned eletrode positions
% cfg.grid                   = grid; % dipole positions
% cfg.headmodel              = vol;  % headmodel created by 'ft_prepare_headmodel'
% cfg.grid.resolution        = 6; % number for automatic sourcemodel generation
% cfg.normalize       = 'yes';
% [lf]        = ft_prepare_leadfield(cfg); % calculate leadfield
% -----------
% An full workflow of creating leadfield using real headmodel (BEM and FEM)
% can be seen in:
% http://www.fieldtriptoolbox.org/workshop/baci2017/forwardproblem/

% A full example script for leadfield based on FEM headmodel can be seen in
% http://www.fieldtriptoolbox.org/development/project/example_fem/
% -------------------------------------------------------------------------
% rest re-referencing
cfg           = [];
cfg.reref     = 'yes';
cfg.refmethod = 'rest';     %  if select 'rest','leadfield' is required.
cfg.leadfield = lf;
%              The leadfield can be a matrix (channels X sources)
%              which is calculated by using the forward theory, based on
%              the electrode montage, head model and equivalent source
%              model. It can also be the output of ft_prepare_leadfield.m
%              (e.g. lf.leadfield) based on real head modal using FieldTrip.

% cfg.refchannel = data.label([1:3,5:60],1); % use first 60 channels
cfg.refchannel = {'all'};    % vector with indices of the selected channels
%                                  %  (re-referenced channels), or 'all'.
data_eeg_rest  = ft_preprocessing(cfg,data);
% -------------------------------------------------------------------------
% median re-referencing
cfg = [];
cfg.reref       = 'yes';
cfg.refmethod   = 'median';    % median
cfg.refchannel  = {'all'};     % use first 60 channels

data_eeg_median = ft_preprocessing(cfg,data);

% plot the results

figure;
k = 1; % kth channel
data1 = data_eeg_rest.trial{1,1};
data2 = data_eeg_median.trial{1,1};

plot(data1(k,1:100),'-b');
hold on;
plot(data2(k,1:100),'-r');
plot(data1(k,1:100)-data2(k,1:100),'-.k')
hold off;grid on;
legend('rest','median','difference');


text:
function linpack_s_test15 ( )

%*****************************************************************************80
%
%% TEST15 tests SPBFA and SPBSL.
%
%  Discussion:
%
%    SPBFA and SPBSL are for a positive definite symmetric band matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  lda = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST15\n' );
  fprintf ( 1, '  For a positive definite symmetric band matrix,\n' );
  fprintf ( 1, '  SPBFA computes the LU factors.\n' );
  fprintf ( 1, '  SPBSL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to matrix A and right hand side B.
%
%  The problem is just an enlarged version of the
%  problem for N = 5, which is:
%
%  Matrix A is ( 2 -1  0  0  0)    right hand side B is  (1)
%              (-1  2 -1  0  0)                          (0)
%              ( 0 -1  2 -1  0)                          (0)
%              ( 0  0 -1  2 -1)                          (0)
%              ( 0  0  0 -1  2)                          (1)
%
%
%  solution is   (1)
%                (1)
%                (1)
%                (1)
%                (1)
%
%  Set the right hand side.
%
  b(1) =     1.0;
  b(2:n-1) = 0.0;
  b(n) =     1.0;
%
%  Set the number of nonzero diagonals.
%
  m = 1;
%
%  Set the value of the subdiagonal and diagonal.
%
  a(1,1:n) = -1.0;
  a(2,1:n) =  2.0;
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, info ] = spbfa ( a, lda, n, m );

  if ( info ~= 0 )
    fprintf ( 1, '  Error!  SPBFA returns INFO = %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  b = spbsl ( a, lda, n, m, b );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last 5 entries of the solution:\n' );
  fprintf ( 1, '  (All should be 1):\n' );
  fprintf ( 1, '\n' );

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

  return
end

author: johannesgerer

repo: jburkardt-m