text:
%MDL_OFFSET6 A minimalistic 6DOF robot arm with shoulder offset
%
% MDL_OFFSET6 is a script that creates the workspace variable off6 which
% describes the kinematic characteristics of a simple arm manipulator with
% a spherical wrist and a shoulder offset, using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%
% Notes::
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also SerialLink, mdl_simple6, mdl_puma560, mdl_twolink.

% MODEL: generic, 6DOF, standard_DH

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function r = mdl_offset6()

    % robot length values (metres)
    L1 = 1;
    L2 = 1;
    O1 = 0;
    O2 = 0.2;
    O3 = 0;

    % and build a serial link manipulator

    robot = SerialLink([
        Revolute('alpha', -pi/2, 'a', O1, 'd', 0)
        Revolute('alpha', 0,     'a', L1, 'd', O2)
        Revolute('alpha', pi/2,  'a', L2, 'd', O3)
        Revolute('alpha', pi/2, 'a', 0,  'd', 0)
        Revolute('alpha', -pi/2,  'a', 0,  'd', 0)
        Revolute('alpha', 0,     'a', 0,  'd', 0)
        ], ...
        'name', 'Offset6');

    % place the variables into the global workspace
    if nargout == 1
        r = robot;
    elseif nargout == 0
        assignin('base', 'off6', robot);
        assignin('base', 'qz', [0 0 0 0 0 0]); % zero angles, arm up
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function step_drawRF(theConf,Name_batch,layerID,epochNum)
showNum=10;
integrateEnergy=0.003; %0.005;
visualizationStep=50;

rng(1);

fileRoot='./draw_RF';
mkdir(fileRoot);
createImages(Name_batch,layerID,fileRoot,epochNum,showNum,integrateEnergy,theConf,visualizationStep);
end

function createImages(nameList,layerID,fileRoot,epochNum,showNum,integrateEnergy,theConf,visualizationStep)
num=numel(nameList);
for i=1:num
    filename=sprintf('./%s/%s/net-epoch-%d.mat',theConf.output.dir,nameList{i},epochNum);
    net=load(filename);
    theConf.output.dir=['./',fileRoot,'/'];
    theConf.data.Name_batch=nameList{i};
    generateAMTOriginalImages(net.net,layerID,showNum,integrateEnergy,theConf,visualizationStep);
    drawRawFMaskF(net.net,layerID,showNum,theConf,visualizationStep);
    computeRF(net.net,layerID,showNum,theConf,visualizationStep);
    drawRF(net.net,layerID,theConf,visualizationStep);
end
end

function generateAMTOriginalImages(net,layerID,showNum,integrateEnergy,theConf,visualizationStep)
batchSize=32;
samplingTimes=1000;

objset=readAnnotation(theConf.data.Name_batch,theConf);
if(strcmp(theConf.data.Name_batch,'cub200'))
    objset=objset(1:500);
end
imgNum=numel(objset);
root_ori=fullfile(theConf.output.dir,'img_ori');
mkdir(root_ori);
convnet=getConvNetPara(net);
net=vl_simplenn_tidy(net);
net=our_vl_simplenn_move(net,'gpu');
assert(layerID>0);
labelNum=size(net.layers{end}.class,3);
ih=net.meta.normalization.imageSize(1);
iw=net.meta.normalization.imageSize(2);
isFlip=zeros(1,imgNum);
targetLayer=convnet.targetLayers(layerID)-1;
upperLayer=convnet.targetLayers(layerID+1)-1;
a=load(sprintf('./mat/%s/imdb.mat',theConf.data.Name_batch),'meta');
for img=1:batchSize:imgNum
    imgList=img:min(img+batchSize-1,imgNum);
    batS=numel(imgList);
    net.layers{end}.class=gpuArray(ones(1,1,labelNum,batS));
    images=getImg(imgList,theConf,[ih,iw],objset,a.meta);
    res=our_vl_simplenn(net,gpuArray(images));
    clear images
    x=res(targetLayer+1).x;
    x=maskX(x,net,upperLayer);
    if(img==1)
        depth=size(x,3);
        depStat=zeros(depth,samplingTimes*imgNum,'single');
        depMax=zeros(depth,imgNum,'single');
    end
    for dep=1:depth
        tmp=x(:,:,dep,:);
        tmplen=numel(tmp);
        if(tmplen>0)
            tmp=reshape(tmp(randperm(tmplen)),[1,tmplen]);
            inputlen=numel(imgList)*samplingTimes;
            tmp=repmat(tmp,[1,ceil(inputlen/tmplen)]);
            depStat(dep,(img-1)*samplingTimes+(1:inputlen))=tmp(1:inputlen);
        end
    end
    depMax(:,imgList)=reshape(max(max(x,[],1),[],2),[depth,batS]);
end
[depStat,~]=sort(depStat,2,'descend');
validMap=zeros(depth,imgNum);
tarP=round(integrateEnergy*size(depStat,2));
tau=depStat(:,tarP);
clear depStat
for dep=1:depth
    validList=find(depMax(dep,:)>=tau(dep));
    validNum=numel(validList);
    if(validNum<showNum)
        [tttmp,validList]=sort(depMax(dep,:),'descend');
        validList=sort(validList(1:showNum));
        validNum=showNum;
        tau(dep)=tttmp(showNum);
    end
    validMap(dep,validList)=true;
    assert(validNum>=showNum);
    tmp=randperm(validNum);
    validMap(dep,validList(tmp(showNum+1:end)))=false;
end
count=zeros(depth,1);
validDepthVisualization=1:visualizationStep:depth;
imgList=1:imgNum;
for dep=1:depth
    if(ismember(dep,validDepthVisualization))
        validList=find(validMap(dep,imgList)>0);
        validNum=numel(validList);
        for id=1:validNum
            count(dep)=count(dep)+1;
            imgID=imgList(validList(id));
            obj=objset(imgID);
            isF=isFlip(imgID);
            Iori=getI(obj,theConf,[ih,iw],isF);
            filename=sprintf('%s/%s-%04d-%05d.jpg',root_ori,theConf.data.Name_batch,dep,count(dep));
            imwrite(uint8(Iori),filename);
        end
    end
end
filename=sprintf('%s/%s_tau.mat',root_ori,theConf.data.Name_batch);
save(filename,'tau');
end

function computeRF(net,layerID,showNum,theConf,visualizationStep)
batchSize=64;
stride=9; %6; %3;
occludeSize=33; %22; %11;
topOccNum=4; %15;

root_ori=fullfile(theConf.output.dir,'img_ori');
convnet=getConvNetPara(net);
net=vl_simplenn_tidy(net);
net=our_vl_simplenn_move(net,'gpu');
assert(layerID>0);
ih=net.meta.normalization.imageSize(1);
iw=net.meta.normalization.imageSize(2);
targetLayer=convnet.targetLayers(layerID)-1;
depth=size(net.layers{targetLayer}.weights{1},4);
occludeNum=ceil((ih-occludeSize+1)/stride)*ceil((iw-occludeSize+1)/stride);
occludeList=repmat(struct('sh',[],'sw',[]),[1,occludeNum]);
c=0;
for sh=1:stride:ih-occludeSize+1
    for sw=1:stride:iw-occludeSize+1
        c=c+1;
        occludeList(c).sh=sh;
        occludeList(c).sw=sw;
    end
end
a=load(sprintf('./mat/%s/imdb.mat',theConf.data.Name_batch),'meta');
for dep=1:visualizationStep:depth
    filename_out=sprintf('%s/%s-%04d_template.mat',root_ori,theConf.data.Name_batch,dep);
    if(exist(filename_out,'file'))
        disp(filename_out);
        continue;
    end
    template=zeros(ih*2,iw*2,showNum);
    for img=1:showNum
        filename=sprintf('%s/%s-%04d-%05d.jpg',root_ori,theConf.data.Name_batch,dep,img);
        I=imread(filename);
        [oriValue,~,m]=getValue(net,single(I)-single(a.meta.dataMean),dep,convnet,layerID);
        if(img==1)
            map=zeros(size(m,1),size(m,2),showNum);
        end
        map(:,:,img)=m;
        tarValue=zeros(1,occludeNum);
        pHW=zeros(2,occludeNum);
        for occStart=1:batchSize:occludeNum
            occList=occStart:min(occStart+batchSize,occludeNum);
            images=getImg_occlude(occludeList(occList),I,occludeSize,a.meta);
            [tarValue(occList),pHW(:,occList),~]=getValue(net,images,dep,convnet,layerID);
            clear images
        end
        [v,idx]=sort(tarValue,'ascend');
        for i=1:topOccNum
            id=idx(i);
            val=max(oriValue-v(i),0);
            hmin=occludeList(id).sh-pHW(1,id);
            hmax=occludeList(id).sh+occludeSize-1-pHW(1,id);
            wmin=occludeList(id).sw-pHW(2,id);
            wmax=occludeList(id).sw+occludeSize-1-pHW(2,id);
            hmin=max(round(hmin+ih),1);
            hmax=min(round(hmax+ih),ih*2);
            wmin=max(round(wmin+iw),1);
            wmax=min(round(wmax+iw),iw*2);
            template(hmin:hmax,wmin:wmax,img)=template(hmin:hmax,wmin:wmax,img)+val;
        end
    end
    save(filename_out,'template','map');
    disp(filename_out);
    clear template map
end
end

function drawRF(net,layerID,theConf,visualizationStep)
drawThreshold=0.5;

root_ori=fullfile(theConf.output.dir,'img_ori');
root_RF=fullfile(theConf.output.dir,'img_RF');
mkdir(root_RF);
filename=sprintf('%s/%s_tau.mat',root_ori,theConf.data.Name_batch);
tau=load(filename,'tau');
convnet=getConvNetPara(net);
assert(layerID>0);
ih=net.meta.normalization.imageSize(1);
iw=net.meta.normalization.imageSize(2);
targetLayer=convnet.targetLayers(layerID)-1;
depth=size(net.layers{targetLayer}.weights{1},4);
for dep=1:visualizationStep:depth
    filename=sprintf('%s/%s-%04d_template.mat',root_ori,theConf.data.Name_batch,dep);
    a=load(filename,'template','map');
    [xh,xw,showNum]=size(a.map);
    RFMask=getRFMask(a,ih,iw,xh,xw,layerID,convnet);
    for img=1:showNum
        filename=sprintf('%s/%s-%04d-%05d.jpg',root_ori,theConf.data.Name_batch,dep,img);
        I=imread(filename);
        assert((size(I,1)==ih)&&(size(I,2)==iw));
        map=a.map(:,:,img);
        theTau=tau.tau(dep);
        wei=getRegion(ih,iw,xh,xw,map,theTau,drawThreshold,RFMask);
        I_RF=drawEdge(wei,I);
        filename=sprintf('%s/%s-%04d-%05d.jpg',root_RF,theConf.data.Name_batch,dep,img);
        imwrite(uint8(I_RF),filename);
    end
end
end

function RFMask=getRFMask(a,ih,iw,xh,xw,layerID,convnet)
T=sum(a.template,3);
Gauss=fspecial('gaussian',50,30);
T=imfilter(T,Gauss,'replicate','same');
RFMask=zeros(ih,iw,xh*xw,'single');
c=0;
for W=1:xw
    for H=1:xh
        c=c+1;
        thePos=x2p_([H;W],layerID,convnet);
        hmin=ih-thePos(1)+1;
        wmin=iw-thePos(2)+1;
        RFMask(:,:,c)=T(hmin:hmin+ih-1,wmin:wmin+ih-1);
    end
end
RFMask=gpuArray(RFMask);
end

function region=getRegion(ih,iw,xh,xw,map,tau,drawThreshold,RFMask)
map=reshape(map,[1,1,xh*xw]);
tau=min(tau,max(map(:)));
map(map<tau)=0;
map=min(map.*(1.5/max(max(map(:)),0.00000000001)),1);
map(map>=tau)=1;
map=repmat(map,[ih,iw,1]);
wei=max(RFMask.*map,[],3);
region=single(gather(wei>=drawThreshold*max(wei(:))));
end

function [tarValue,pHW,map]=getValue(net,images,dep,convnet,layerID)
batchS=size(images,4);
net.layers{end}.class=gpuArray(ones(1,1,size(net.layers{end}.class,3),batchS));
res=our_vl_simplenn(net,gpuArray(images));
targetLayer=convnet.targetLayers(layerID)-1;
upperLayer=convnet.targetLayers(layerID+1)-1;
x=res(targetLayer+1).x;
x=maskX(x,net,upperLayer);
[xh,xw,~,~]=size(x);
xx=reshape(x(:,:,dep,:),[xh*xw,batchS]);
[~,idx]=max(xx,[],1);
pw=reshape(ceil(idx./xh),[1,batchS]);
ph=reshape(idx-(pw-1).*xh,[1,batchS]);
pHW=x2p_([ph;pw],layerID,convnet);
tarValue=reshape(sum(sum(x(:,:,dep,:),1),2),[batchS,1]);
map=reshape(gather(x(:,:,dep,:)),[xh,xw,batchS]);
end

function images=getImg_occlude(occludeList,I,occludeSize,meta)
batS=numel(occludeList);
images=repmat(single(I)-single(meta.dataMean),[1,1,1,batS]);
for i=1:batS
    hList=occludeList(i).sh:occludeList(i).sh+occludeSize-1;
    wList=occludeList(i).sw:occludeList(i).sw+occludeSize-1;
    images(hList,wList,:,i)=0;
end
end

function images=getImg(imgList,theConf,tarSize,objset,meta)
batS=numel(imgList);
ih=tarSize(1);
iw=tarSize(2);
images=zeros(ih,iw,3,batS,'single');
isFlip=false;
for i=1:batS
    I=getI(objset(imgList(i)),theConf,tarSize,isFlip);
    images(:,:,:,i)=single(I);
end
images=images-repmat(meta.dataMean,[1,1,1,batS]);
end

function [x,mask]=maskX(x,net,upperLayer)
x=max(x,0);
if(~strcmp(net.layers{upperLayer}.type,'conv_mask'))
    x=gather(x);
    return;
end
[h,w,depth,batchS]=size(x);
%[mu_x,mu_y,~]=getMu(x);
[mu_x,mu_y]=getMu_max(x);
posTempX=gpuArray(single(repmat(linspace(-1,1,h)',[1,w,depth,batchS])));
posTempY=gpuArray(single(repmat(linspace(-1,1,w),[h,1,depth,batchS])));
mask=abs(posTempX-repmat(mu_x,[h,w,1,1]));
mask=mask+abs(posTempY-repmat(mu_y,[h,w,1,1]));
mask=max(1-mask.*repmat(reshape(net.layers{upperLayer}.weights{3},[1,1,depth]),[h,w,1,batchS]),0);
mask(:,:,net.layers{upperLayer}.filter~=1,:)=1;
%assert(isempty(find(net.layers{upperLayer}.filter~=1)));
x=gather(x.*mask);
end

function [mu_x,mu_y]=getMu_max(x)
[h,w,depth,batchS]=size(x);
x=reshape(x,[h*w,depth,batchS]);
[~,p]=max(x,[],1);
p=reshape(p,[1,1,depth,batchS]);
mu_y=ceil(p./h);
mu_x=p-(mu_y-1).*h;
tmp=linspace(-1,1,h);
mu_x=reshape(tmp(mu_x),[1,1,depth,batchS]);
mu_y=reshape(tmp(mu_y),[1,1,depth,batchS]);
end

function drawRawFMaskF(net,layerID,showNum,theConf,visualizationStep)
root_ori=fullfile(theConf.output.dir,'img_ori');
root_maskF=fullfile(theConf.output.dir,'img_maskF');
mkdir(root_maskF);
convnet=getConvNetPara(net);
net=vl_simplenn_tidy(net);
net=our_vl_simplenn_move(net,'gpu');
assert(layerID>0);
ih=net.meta.normalization.imageSize(1);
targetLayer=convnet.targetLayers(layerID)-1;
depth=size(net.layers{targetLayer}.weights{1},4);
a=load(sprintf('./mat/%s/imdb.mat',theConf.data.Name_batch),'meta');
for dep=1:visualizationStep:depth
    for img=1:showNum
        filename=sprintf('%s/%s-%04d-%05d.jpg',root_ori,theConf.data.Name_batch,dep,img);
        I=imread(filename);
        [raw_x,mask,mask_x]=getMaskMap(net,single(I)-single(a.meta.dataMean),dep,convnet,layerID);
        filename=sprintf('%s/%s-%04d-%05d_rawX.jpg',root_maskF,theConf.data.Name_batch,dep,img);
        imagewrite(raw_x,filename);
        filename=sprintf('%s/%s-%04d-%05d_mask.jpg',root_maskF,theConf.data.Name_batch,dep,img);
        imagewrite(mask,filename);
        filename=sprintf('%s/%s-%04d-%05d_maskX.jpg',root_maskF,theConf.data.Name_batch,dep,img);
        imagewrite(mask_x,filename);
    end
    fprintf('%s-%04d.jpg\n',theConf.data.Name_batch,dep)
end
end

function [raw_x,mask,mask_x]=getMaskMap(net,images,dep,convnet,layerID)
net.layers{end}.class=gpuArray(ones(1,1,size(net.layers{end}.class,3),1));
res=our_vl_simplenn(net,gpuArray(images));
targetLayer=convnet.targetLayers(layerID)-1;
upperLayer=convnet.targetLayers(layerID+1)-1;
raw_x=res(targetLayer+1).x;
[mask_x,mask]=maskX(raw_x,net,upperLayer);
raw_x=gather(raw_x(:,:,dep));
mask=gather(mask(:,:,dep));
mask_x=gather(mask_x(:,:,dep));
end

function imagewrite(x,filename)
scaleW=16;

[xh,xw]=size(x);
x=uint8(x./max(x(:)).*255);
x=reshape(repmat(reshape(x,[1,xh,xw]),[scaleW,1,1]),[xh*scaleW,xw]);
x=reshape(repmat(reshape(x',[1,xw,xh*scaleW]),[scaleW,1,1]),[xw*scaleW,xh*scaleW]);
x=x';
imwrite(x,filename);
end

repo: interpretableCNN


text:
function [Acc]=parameter_sensitivity_test_loocv(All_Feat,nSubj,label,meth_Net,lambda_lasso)
% This function performs parameter sensitivity test using all the subject data with LOOCV, and
% it is run across all the combination of parameters.
% Input:
%       All_Feat: the constructed brain network using one combination of parameter;
%       nSubj: number of subjects;
%       label: the label for each subject; e.g., -1 for normal controls and 1 for patients;
%       lambda_lasso: lambda in the lasso feature selection;
%
% Written by Zhen Zhou, zzstefan@email.unc.edu
% IDEA lab, https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, University of North Carolina at Chapel Hill

e = 1:nSubj;
cpred = zeros(nSubj,1);
score = zeros(nSubj,1);

% LOOCV for testing
for i=1:nSubj
    Tst_ind = i;
    telabel = label(i);
    Trn_ind = e;
    Trn_ind(i) = [];
    trlabel = label;
    trlabel(i) = [];

    Feat = All_Feat;
    trFe = Feat(Trn_ind,:);
    teFe = Feat(Tst_ind,:);

    % Feature selection ag
    if ~strcmpi(meth_Net,'dHOFC')
        pval=0.05;  % generally use pval<0.05 as a threshold for feature selection
        [~,p]=ttest2(trFe(trlabel==-1,:),trFe(trlabel==1,:));
        trFe=trFe(:,p<pval);
        teFe=teFe(:,p<pval);

        midw=lasso(trFe,trlabel,'Lambda',lambda_lasso);  % parameter lambda for sparsity
        trFe=trFe(:,midw~=0);
        teFe=teFe(:,midw~=0);
    else
        midw=lasso(trFe,trlabel,'Lambda',lambda_lasso);  % parameter lambda for sparsity
        trFe=trFe(:,midw~=0);
        teFe=teFe(:,midw~=0);
    end

    % Feature normalization ag
    Mtr=mean(trFe);
    Str=std(trFe);
    trFe=trFe-repmat(Mtr,size(trFe,1),1);
    trFe=trFe./repmat(Str,size(trFe,1),1);
    teFe=teFe-Mtr;
    teFe=teFe./Str;

    % train SVM model ag
    classmodel=svmtrain(trlabel,trFe,'-t 0 -c 1 -q'); % linear SVM (require LIBSVM toolbox)
    % classify ag
    [cpred(i),~,score(i)]=svmpredict(telabel,teFe,classmodel,'-q');
end

Acc=100*sum(cpred==label)/nSubj;
%[AUC,SEN,SPE,F1]=perfeval(label,cpred,score);

author: zzstefan

repo: BrainNetClass


text:
function InsertHeader(varargin)
% InsertHeader - Adds a template to the top of current file in editor for documentation.
%
% Syntax: InsertHeader(varargin)
% varargin - optional input of full file name
%   Example 1:  InsertHeader adds header to top of currently selected file in
%   editor
%   Example 2: InsertHeader('InsertHeader') added the current header information to this
%   file. NOTE:  THIS CANNOT BE UNDONE!  THE FIRST TIME YOU USE THIS, PLEASE
%   HAVE A BACK UP COPY OF THE FILE IN CASE IT MAKES A MISTAKE.
%   Custom defined headers can be created in the subfunction defined at
%   the end of the file.  I wish to express thanks to the authors of tedit.m and
%   newfcn.m which helped me to create this routine.
%  See also:  tedit, newfcn
%
% Subfunctions: EditorCurrentFile, parseForSubFcns, insertHeaderTemplateFile

% AUTHOR    : Dan Kominsky
% Copyright 2011  Prime Photonics, LC.
%%

if nargin==1
	CurrentFile=varargin{1};
else
	% Call subfunction to determine the current file:
	[CurrentFile]=EditorCurrentFile;
end
% Ask to make sure I am going to add the header to the correct file
resp=questdlg(['Add Header to file: ' CurrentFile '?'],'Header Query');
if ~strcmpi(resp,'Yes'),return;end;

if verLessThan('matlab','7.12.0') % Matlab interface changed with 7.12.0
  old_style = true;
% Define the handle for the java commands:
edhandle=com.mathworks.mlservices.MLEditorServices;
% Save the current document so changes aren't lost
edhandle.saveDocument(CurrentFile);
[PATH,FILENAME]=fileparts(CurrentFile);

% Load the existing file
rawText=char(edhandle.builtinGetDocumentText(CurrentFile));

else
  old_style = false;
  % Define the handle for the java commands:
  editorObject = matlab.desktop.editor.getActive;
  % Save the current document so changes aren't lost
  editorObject.save;
  % Load the text:
  rawText=editorObject.Text;

  [PATH,FILENAME]=fileparts(CurrentFile);

end

% For some reason the previous command replaces (at least on windows) each
% carraige return with a carraige return/line feed combo.  The next line removes
% the double spacing.   It may not be necessary on other platforms!
% In some cases it might use just the char(10), so convert that to char(13);
currentText=regexprep(rawText,[char(13) char(10)],char(13));
currentText=regexprep(currentText,char(10),char(13));

% Check to find the true end of the first line.  If it is a continuation, keep
% going until the function definition is ended.
ellipsisTest=true;
firstLine='';
remainder=currentText;
while ellipsisTest
	[tmpfirstLine,remainder]=strtok(remainder,char(13));
	firstLine=[firstLine,char(13),tmpfirstLine];
	if sum(firstLine(end-4:end)==46)>=3
		ellipsisTest=true;
	else
		ellipsisTest=false;
	end
	pause(0.01)
end

funcLabel=any(regexp(firstLine,'function'));
syntaxStr={};
inputs={};
outputs={};
count=1;
if funcLabel
	syntaxStr=regexprep(firstLine,'function','');
	hasEquals=any(regexp(syntaxStr,'='));
	if hasEquals
	[beforeEq,afterEq]=strtok(strtrim(syntaxStr),'=');
	else
		beforeEq='';afterEq=strtrim(syntaxStr);
	end
	tmpoutputs=regexp(beforeEq,'(\w+)','tokens');
	for n=1:numel(tmpoutputs)
		outputs{n}=char(tmpoutputs{n});
	end
	tmpinputs=regexp(afterEq,'(\w+)','tokens');
	funcname=char(tmpinputs{1});
	for n=2:numel(tmpinputs)
		inputs{n-1}=char(tmpinputs{n});
	end
end

  subFcnName=parseForSubFcns(CurrentFile);

% insertHeaderTemplateFile contains the text to be added into the file
tmpl=insertHeaderTemplateFile(FILENAME,inputs,outputs,syntaxStr,subFcnName);

% Convert to a string
parseIn = sprintf('%s\n',tmpl{:});

% Concatenate the firstline (function definition), header, and body:
concatenatedText=strcat(firstLine,parseIn,remainder);
% Strip any extra leading carraige returns:
while (double(concatenatedText(1)))==13 || (double(concatenatedText(1)))==10
	concatenatedText(1)=[];
end

  if old_style
% Close the file:
edhandle.closeDocument(CurrentFile);
% Open and OVERWRITE the file
fid=fopen(CurrentFile,'wt');
fprintf(fid,'%s\n',concatenatedText);
fclose(fid);
% Reopen the file to the begining.
edhandle.openDocumentToLine(CurrentFile,1);
  else
    editorObject.save;
    fid=fopen(CurrentFile,'wt');
    fprintf(fid,'%s\n',concatenatedText);
    fclose(fid);
    editorObject.reload;
  end
end

% ------------------------
function [CurFile, varargout] = EditorCurrentFile

  if verLessThan('matlab','7.12.0') % Matlab interface changed with 7.12.0
% Define the handle for the set of java commands:
desktopHandle=com.mathworks.mlservices.MatlabDesktopServices.getDesktop;
% Determine the last selected document in the editor:
lastDocument=desktopHandle.getLastDocumentSelectedInGroup('Editor');
% Strip off the '*' which indicates that it has been modified.
CurFile=strtok(char(desktopHandle.getTitle(lastDocument)),'*');

if nargout>1
	varargout{1}=lastDocument;
end
  else
    editorObject = matlab.desktop.editor.getActive;
    CurFile = editorObject.Filename;
  end

end

function [functions varargout]=parseForSubFcns(myFunc)
  if exist(myFunc,'file')
    % We have a full path to the file
    str= mlintmex('-calls',myFunc);
  else
    % File is on the path, and we must find it:
    str= mlintmex('-calls',which(myFunc));
  end
  [a,b,c,d,e] = regexp(str,'([NS])\d* (\d+) (\d+) (\w+)\n');
  if isempty(e)
    functions = [];varargout{1}=[];
    return
  end
  if nargout < 1
    sprintf('Sub/Nest\t\tLine#\t\tCol#\t\tName')
    for iFcn = 1:numel(e)
      S=sprintf('%8s\t\t%6s\t\t%6s\t\t%s',e{iFcn}{1:4});
      disp(S)
    end
  else
    functions = [];
    for iFcn = 1:numel(e)
      S=sprintf('%s, ',e{iFcn}{4});
      functions = [functions S];
      cellList{iFcn} = e{iFcn}{4};
    end
    functions(end-1:end) = [];
    if nargout>1
      varargout{1} = cellList;
    end
  end
end

function tmpl=insertHeaderTemplateFile(FILENAME,varargin)
% insertHeaderTemplateFile Contains the header which is automatically inserted
%   tmpl = insertHeaderTemplateFile(FILENAME)
%
%   Example
%   insertHeaderTemplateFile('currentFile.m') returns the cell array of strings
%   which is inserted into file headers.
%
%   See also

% AUTHOR    : Dan Kominsky
  % Copyright 2006 Prime Photonics, LC.
%%
applyComment=@(x)['% ' x ' - Description'];
commaPad=@(x)[x{:}, ', '];

if nargin>1
	inputLines=cellfun(applyComment,varargin{1}','UniformOutput',false);
else
	inputLines='';
end
if nargin>2
	outputLines=cellfun(applyComment,varargin{2}','UniformOutput',false);
else
	outputLines='';
end
if nargin>3
  synt_str_raw = strtrim(varargin{3});
  synt_str_eols = regexprep(synt_str_raw,'\n','\n%');
  synt_str_adj = regexprep(synt_str_eols,'\r','\r%');
	syntaxStr=['% Syntax: ' synt_str_adj];
else
	syntaxStr='';
end
if nargin>4
    subFcns = ['% Subfunctions: ' varargin{4}];
else
	subFcns='% Subfunctions: none';
end
topLines={ ...
	''
  ['% ' FILENAME  ' - One line description of what the function or script performs (H1 line)']
	'% Optional file header info (to give more details about the function than in the H1 line)'};
midLines=[syntaxStr; inputLines;outputLines];
bottomLines={	'%   Example'
	'%  Line 1 of example '
	'%'};
authorLines={...
	'%   See also: '
	''
	'% AUTHOR    : Jebediah Obediah'
	['% Copyright ' datestr(now,'yyyy') '  XYZ Corp.']
	'%%'};

tmpl=[topLines;midLines;bottomLines;subFcns;authorLines];
end

author: Sable

repo: mcbench-benchmarks


text:
%PREX_COST PRTools example on cost matrices and rejection
%
% Prtools example code to show the use of cost matrices and how
% to introduce a reject class.

help prex_cost
randstate = randreset;
echo on
                % Generate a three class problem
  randn('state',1);
  rand('state',1);
  n = 30;
  class_labels = char('apple','pear','banana');
  a = [gendatb([n,n]);  gendatgauss(n,[-2 6])];
  laba = genlab([n n n],class_labels);
  a = setlabels(a,laba);
                % Compute a simple ldc
  w = ldc(a);
                % Scatterplot and classifier
  figure;
  gridsize(30);
  scatterd(a,'legend');
  plotc(w);

                % Define a classifier with a new cost matrix,
					 % which puts a high cost on misclassifying
					 % pears to apples
  cost = [0.0  1.0  1.0;
          9.0  0.0  1.0;
		    1.0  1.0  0.0];
  wc = w*classc*costm([],cost,class_labels);
  plotc(wc,'b');

          % Define a classifier with a cost matrix where
					 % an outlier class is introduced. For this an
					 % extra column in the cost matrix has to be defined.
					 % Furthermore, the class labels have to be supplied
					 % to give the new class a name.
  cost = [0.0  1.0  1.0  0.2;
          9.0  0.0  1.0  0.2;
          1.0  1.0  0.0  0.2];
  class_labels = char('apple','pear','banana','reject');
  wr = w*classc*costm([],cost,class_labels);
  plotc(wr,'--')

echo off
randreset(randstate);
disp(' ')
disp('   The black decision boundary shows the standard ldc classifier');
disp('   for this data. When the misclassification cost of a pear to an');
disp('   apple is increased, we obtain the blue classifier. When on top');
disp('   of that a rejection class is introduced, we get the blue dashed');
disp('   classifier. In that case, all objects between the dashed lines');
disp('   are rejected.');
fprintf('\n');
fprintf('  Cost of basic classifier  =  %4.2f\n',...
             a*w*testcost([],cost,class_labels));
fprintf('  Cost of cost classifier   =  %4.2f\n',...
             a*wc*testcost([],cost,class_labels));
fprintf('  Cost of reject classifier =  %4.2f\n',...
             a*wr*testcost([],cost,class_labels));

author: marianux

repo: ecg-kit

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


text:
function createHelpDemoDocumentation

%% Get toolbox location

filePath=mfilename('fullpath');
toolboxPath=fileparts(fileparts(filePath));
helpPath=fullfile(toolboxPath,'docs','html');

%% FIND HTML FILES

allFiles = dir(fullfile(helpPath,'*.html'));
allFiles={allFiles(1:end).name};
allFiles=sort(allFiles(:));

logicDemo=gcontains(allFiles,'DEMO_');

demoFiles=allFiles(logicDemo);
numDemoFiles=numel(demoFiles);

logicHelp=gcontains(allFiles,'HELP_');

helpFiles=allFiles(logicHelp);
numHelpFiles=numel(helpFiles);

logicStart=gcontains(allFiles,'GIBBON_product_page');

startFile=allFiles{logicStart};

%% ADJUST helptoc.XML
%TOP SECTION
T_top={};
T_top(1,1)={'<?xml version="1.0" encoding="utf-8"?>'};
T_top(2,1)={'<toc version="2.0">'};
T_top(3,1)={['    <tocitem target="',startFile,'">The GIBBON Toolbox']};

%GETTING STARTED SECTION
T_start=cell(3,1);
T_start{1}='        <tocitem target="GettingStarted.html" image="HelpIcon.GETTING_STARTED">Getting Started';
T_start{2}=['            <tocitem target="',startFile,'">',startFile(1:end-5),'</tocitem>',startFile(1:end-4)];
T_start{3}='        </tocitem>';

%FUNCTION HELP
T_help=cell(numHelpFiles+2,1);
T_help{1}='        <tocitem target="funclist.html" image="HelpIcon.FUNCTION">Functions';
for q=1:1:numHelpFiles
    currentFile=helpFiles{q};
    currentName=currentFile(1:end-5);
    T_help{q+1}=['            <tocitem target="',currentFile,'">',currentName(6:end),'</tocitem>',currentName];
end
T_help{end}='        </tocitem>';

%DEMO EXAMPLES
T_demo=cell(numDemoFiles+2,1);
T_demo{1}='        <tocitem target="gibbonExampes.html" image="HelpIcon.EXAMPLES">Examples';
for q=1:1:numDemoFiles
    currentFile=demoFiles{q};
    currentName=currentFile(1:end-5);
    T_demo{q+1}=['            <tocitem target="',currentFile,'">',currentName(6:end),'</tocitem>',currentName];
end
T_demo{end}='        </tocitem>';

%BOTTOM SECTION
T_bottom={};
T_bottom(1,1)={'    </tocitem>'};
T_bottom(2,1)={'</toc>'};

T=[T_top;T_start;T_help;T_demo;T_bottom];

saveName=fullfile(helpPath,'helptoc.xml');

cell2txtfile(saveName,T,0,0);

%% Add searchable help
addHelpSearch;

disp('You may need to restart MATLAB to allow for the help and documentation integration changes to take effect');

disp('Use gdoc command to open GIBBON documentation page in MATLAB, or use https://www.gibboncode.org/Documentation/');

%%
% _*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:
% mff_importcoordinates - import information from MFF 'coordinates.xml' file
%
% Usage:
%   [chanlocs, ref] = mff_importcoordinates(mffFile);
%
% Inputs:
%  mffFile - filename/foldername for the MFF file
%
% Output:
%  chanlocs - EEGLAB channel location structure
%  ref      - reference

% This file is part of mffmatlabio.
%
% mffmatlabio 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.
%
% mffmatlabio 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 mffmatlabio.  If not, see <https://www.gnu.org/licenses/>.

function [chanlocs, reference, rVal] = mff_importcoordinates(mffFile)

chanlocs = [];
rVal = true;

p = fileparts(which('mff_importsignal.m'));
warning('off', 'MATLAB:Java:DuplicateClass');
javaaddpath(fullfile(p, 'MFF-1.2.2-jar-with-dependencies.jar'));
warning('on', 'MATLAB:Java:DuplicateClass');

% Create an MFFFactory object.
mfffactorydelegate = javaObject('com.egi.services.mff.api.LocalMFFFactoryDelegate');
mfffactory = javaObject('com.egi.services.mff.api.MFFFactory', mfffactorydelegate);

% Create Signal object and read in event track file.
% Note that 3 is the value associated with the event track resource type.
sLayoutURI = fullfile(mffFile, 'coordinates.xml');
sensorLayoutType = javaObject('com.egi.services.mff.api.MFFResourceType', javaMethod('valueOf', 'com.egi.services.mff.api.MFFResourceType$MFFResourceTypes', 'kMFF_RT_Coordinates'));

cLayout = mfffactory.openResourceAtURI(sLayoutURI, sensorLayoutType);

listfid  = [];
ref      = [];
reference = {};
if ~isempty(cLayout)

    %try

    if cLayout.loadResource()

        %  The information about the Coordinate Layout Resource.
        fprintf('Coordinate acquisition Time: %s\n', char(cLayout.getAcquisitionTime()));
        fprintf('Coordinate acquisition Method: %s\n', char(cLayout.getAcquisitionMethod()));

        %  Check to see if we have a default subject.
        %  fprintf('Default Subject: %d\n', cLayout.getDefaultSubject());

        %  Get the sesnors.
        sLayout =  cLayout.getSensorLayout();

        if ~isempty(sLayout )

            %  Information about the sensor layout.

            fprintf('Coordinate Layout Name: %s\n', char(sLayout.getName()));

            %  If Montage
            if ~isempty(sLayout.getOriginalLayout())
                fprintf('Cordinate Original Layout Name: %s\n', char(sLayout.getOriginalLayout()) );
            end
            sensors = sLayout.getSensors();

            if ~isempty(sensors )

                for iSensor = 1:sensors.size
                    sensor = sensors.get(iSensor-1);
                    description = char(sensor.getName());
                    chanlocs(iSensor).labels = [ 'E' num2str(sensor.getNumber()) ];
                    chanlocs(iSensor).description = description;
                    chanlocs(iSensor).X = sensor.getX();
                    chanlocs(iSensor).Y = sensor.getY();
                    chanlocs(iSensor).Z = sensor.getZ();
                    chanlocs(iSensor).identifier = sensor.getIdentifier();
                    sensorType = sensor.getType();
                    switch sensorType
                        case 0, chanlocs(iSensor).type = 'EEG';
                        case 1, chanlocs(iSensor).type = 'EEG'; ref = [ref iSensor]; reference = { reference{:} chanlocs(iSensor).labels };
                        case 2, chanlocs(iSensor).type = 'FID'; listfid = [ listfid iSensor ];
                    end
                end

            end
        else
            fprintf('Error: Coordinate Layout is null\n');
            rVal = false;
        end

    else
        fprintf('Error: could not load the coordinates resource.\n');
        rVal = false;
    end

end

if ~isempty(ref)
    for iSensor = 1:length(chanlocs)
        if strcmpi(chanlocs(iSensor).type, 'EEG'), chanlocs(iSensor).ref = [ chanlocs(ref).labels ]; end
    end
end

if ~isempty(chanlocs)
    chanlocs = convertlocs(chanlocs, 'cart2all');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Figure 9.23       Feedback Control of Dynamic Systems, 5e
%                      Franklin, Powell, Emami
%
clf;
ka=10;
sim('fig9_22')
figure(1)
plot(uantib(:,1),uantib(:,2))
axis([0 10 -1.2 1.2])
hold on
figure(2)
plot(yantib(:,1),yantib(:,2))
hold on
ka=0;
sim('fig9_22')
figure(1)
plot(uantib(:,1),uantib(:,2))
title('Fig. 9.23(b) Control with and without antiwindup')
xlabel('Time (sec)');
ylabel('Control');
nicegrid;
gtext('With anti-windup')
gtext('Without anti-winduup')
figure(2)
plot(yantib(:,1),yantib(:,2))
xlabel('Time (sec)');
ylabel('Output');
title('Fig. 9.23(a) Output with and without antiwindup')
gtext('With anti-windup')
gtext('Without anti-winduup')
nicegrid;
hold off


text:
function test_nint_test ( )

%*****************************************************************************80
%
%% TEST_NINT_TEST tests the TEST_NINT library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_NINT_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TEST_NINT library.\n' );

  test_nint_test01 ( );
  test_nint_test02 ( );
  test_nint_test03 ( );
  test_nint_test04 ( );
  test_nint_test05 ( );
  test_nint_test06 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_NINT_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef EIMEGO < ALGORITHM
% <multi> <real/integer> <expensive>
% Expected improvement matrix based efficient global optimization
% InfillCriterionIndex --- 1 --- infill criterion index number

%------------------------------- Reference --------------------------------
% D. Zhan, Y. Cheng, and J. Liu, Expected improvement matrix-based infill
% criteria for expensive multiobjective optimization, IEEE Transactions on
% Evolutionary Computation, 2017, 21(6): 956-975.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Dawei Zhan

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            % 1 for Euclidean distance-based EIM criterion (default)
            % 2 for Maximin distance-based EIM criterion
            % 3 for Hypervolume-based EIM criterion
            InfillCriterionIndex = Algorithm.ParameterSet(1);

            %% Generate the initial design points
            % number of design variables
            D = Problem.D;
            % number of objective functions
            M = Problem.M;
            % number of initial design points
            N  = 11*D-1;
            % generate initial design points using Latin Hypercube sampling
            PopDec = repmat(Problem.upper-Problem.lower,N,1).*UniformPoint(N,D,'Latin')+repmat(Problem.lower,N,1);
            % calculate initial design points
            Population   = Problem.Evaluation(PopDec);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % scale the objective values to 0 and 1
                PopDec = Population.decs;
                PopObj = Population.objs;
                N  = size(PopDec,1);
                PopObjScaled = (PopObj-repmat(min(PopObj),N,1))./(repmat(max(PopObj),N,1)-repmat(min(PopObj),N,1));
                % bulid Kriging models for all the objective functions
                KrigingModel = cell(1,M);
                for i = 1 : M
                    KrigingModel{i}= dacefit(PopDec,PopObjScaled(:,i),'regpoly0','corrgauss',1*ones(1,D),0.001*ones(1,D),1000*ones(1,D));
                end
                % select one candidate with the maximum EIM value using GA
                PopDec     = InfillSamplingEIM(Problem,KrigingModel,PopObjScaled,InfillCriterionIndex);
                Population = [Population,Problem.Evaluation(PopDec)];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function S3G = equispacedSO3Grid(CS,SS,varargin)
% defines a equispaced grid in the orientation space
%
% Syntax
%   S3G = equispacedSO3Grid(CS,SS,'points',n)
%   S3G = equispacedSO3Grid(CS,'resolution',res)
%
%   % fill only a ball with radius of 20 degree
%   S3G = equispacedSO3Grid(CS,'maxAngle',20*degree)
%
% Input
%  CS  - @crystalSymmetry
%  SS  - @specimenSymmetry
%   n  - approximate number of points
%  res - resolution in radiant
%
% Output
%  S3G - @SO3Grid
%
% Options
%  maxAngle - radius of the ball to be filles
%  center - center of the ball
%
% See also
% equispacedS2Grid, SO3Grid/SO3Grid

% extract specimen symmetry if provided
if nargin == 1
  SS = specimenSymmetry('1');
elseif ~isa(SS,'symmetry')
  varargin = [{SS},varargin];
  SS = specimenSymmetry('1');
end

% may be we should populate only a ball
maxAngle = get_option(varargin,'maxAngle',2*pi);

if maxAngle < pi/2/CS.multiplicityZ
  S3G = localOrientationGrid(CS,SS,maxAngle,varargin{:});
  return
end

% get fundamental region
[maxAlpha,maxBeta,maxGamma] = fundamentalRegionEuler(CS,SS,'SO3Grid',varargin{:});
maxGamma = maxGamma/2; % we will consider the interval -maxGamma/2 .. maxGamma/2
if ~check_option(varargin,'center'), maxGamma = min(maxGamma,maxAngle);end

% determine resolution
if check_option(varargin,'points')

  points = get_option(varargin,'points');

  switch CS.LaueName  % special case: cubic symmetry
    case 'm-3'
      points = 3*points;
    case 'm-3m'
      points = 2*points;
  end

  % calculate number of subdivisions for the angles alpha,beta,gamma
  res = 2/(points/( maxBeta*maxGamma))^(1/3);

  if  maxAngle < pi * 2 && maxAngle < maxBeta
    res = res * maxAngle; % bug: does not work properly for all syms
  end

else

  res = get_option(varargin,'resolution',5*degree);

end

alphabeta = equispacedS2Grid('resolution',res,...
  'maxTheta',maxBeta,'minRho',0,'maxRho',maxAlpha,...
  no_center(res),'restrict2minmax');

ap2 = round(2*maxGamma/res);

[beta,alpha] = polar(alphabeta);

% calculate gamma shift
re = cos(beta).*cos(alpha) + cos(alpha);
im = -(cos(beta)+1).*sin(alpha);
dGamma = atan2(im,re);
dGamma = repmat(reshape(dGamma,1,[]),ap2,1);
gamma = -maxGamma + (0:ap2-1) * 2 * maxGamma / ap2;

% arrange alpha, beta, gamma
gamma  = dGamma+repmat(gamma.',1,length(alphabeta));
alpha = repmat(reshape(alpha,1,[]),ap2,1);
beta  = repmat(reshape(beta,1,[]),ap2,1);

ori = orientation.byEuler(alpha,beta,gamma,'ZYZ',CS,SS,varargin{:});

gamma = S1Grid(gamma,-maxGamma+dGamma(1,:),...
  maxGamma+dGamma(1,:),'periodic','matrix');

res = 2 * maxGamma / ap2;

% eliminiate 3 fold symmetry axis of cubic symmetries
% TODO: this should be done better!!
ind = fundamental_region(ori,CS,specimenSymmetry);

if nnz(ind) ~= 0
  % eliminate those rotations
  ori(ind) = [];

  % eliminate from index set
  gamma = subGrid(gamma,~ind);
  alphabeta  = subGrid(alphabeta,GridLength(gamma)>0);
  gamma(GridLength(gamma)==0) = [];

end

S3G = SO3Grid(ori,alphabeta,gamma,'resolution',res);

if check_option(varargin,'maxAngle')
  center = get_option(varargin,'center',quaternion.id);
  S3G = subGrid(S3G,center,maxAngle);
end

end

% ----------------------------------------------------------------
function s = no_center(res)

if mod(round(2*pi/res),2) == 0
  s = 'no_center';
else
  s = '';
end
end

% ---------------------------------------------------------------
function ind = fundamental_region(q,cs,ss)

if isempty(q), ind = []; return; end

c = {};

% eliminiate 3 fold symmetry axis of cubic symmetries
switch cs.LaueName

  case   {'m-3m','m-3'}

    c{end+1}.v = vector3d([1 1 1 1 -1 -1 -1 -1],[1 1 -1 -1 1 1 -1 -1],[1 -1 1 -1 1 -1 1 -1]);
    c{end}.h = sqrt(3)/3;

    if strcmp(cs.LaueName,'m-3m')
      c{end+1}.v = vector3d([1 -1 0 0 0 0],[0 0 1 -1 0 0],[0 0 0 0 1 -1]);
      c{end}.h = sqrt(2)-1;
    end
end

switch ss.LaueName
  case 'mmm'
   c{end+1}.v = vector3d([-1 0],[0 -1],[0 0]);
   c{end}.h = 0;
end

% find rotation not part of the fundamental region
rodrigues = Rodrigues(q); clear q;
ind = false(length(rodrigues),1);
for i = 1:length(c)
  for j = 1:length(c{i}.v)
    p = dot(rodrigues,1/norm(c{i}.v(j)) * c{i}.v(j));
    ind = ind | (p(:)>c{i}.h);
  end
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [feat,idxs_bbox_pair,cb1,cb2] = get_spatial_features_same_part_regr_img(locations,idxs_bbox_pair)

cb1 = locations(idxs_bbox_pair(:,1), :);
cb2 = locations(idxs_bbox_pair(:,2), :);

deltaX = cb2(:,1)-cb1(:,1);
deltaY = cb2(:,2)-cb1(:,2);

feat = cat(2, deltaX, deltaY);

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function price = FFTCOS_UpAndOut(n, Nex,H,Rb, L, c, cp, type, S0, t, r, q, ...
                                    strike, varargin)

% Function FFTCOS_UpAndOut calculates the price of a discretely
% monitored Up-and-Out Call or Put option by use of the Fourier-Cosine
% Series Expansion introduced by Oosterlee & Fang
%------------------------------------------------------------------------

% Nex := number of examination points
% H := Barrier
% Rb := Rebate

dt = t / Nex;                    % time interval
Ngrid = 2 ^ n;                   % Grid points
Nstrike = size(strike,1);        % number of strikes

x = double(log(S0 ./ strike));       % center
h = double(log(H ./ strike));

a = double(c(1) + x - L * sqrt(c(2) + sqrt(c(3))));   % lower trunc
b = double(c(1) + x + L * sqrt(c(2) + sqrt(c(3))));   % upper trunc

Grid_i = repmat((0:Ngrid-1)',1,Nstrike);    % Grid index

% Set up function handles
if cp == 1
    vk = @(x) calcv(Grid_i, x, h, a, b, cp, strike);
    if h >= 0
        V = vk(0);
    end
else
    vk = @(x) calcv(Grid_i, a, x, a, b, cp, strike);
    if h >= 0
        V = vk(0);
    else
        V = vk(h);
    end
end

cv = @(y) cvalue(a, h, a, b, Ngrid, y, type, dt, r, q, varargin{:});

aux = pi * Grid_i * diag(1./(b-a));
G = ((sin(aux * diag(b-a))-sin(aux*diag(h-a)))./aux);
G = exp(-r * dt) * 2 * Rb * [((b-h)./(b-a))';G(2:end,:)]* diag(1./(b-a));

V = V + G;
for m = Nex-1:-1:1              % backward induction
    V = cv(V) + G;
end

cfval = exp(feval(@CF, type,aux, dt,r,q,varargin{:}));

pF = cfval .* exp( 1i * aux * diag(x - a) );
pF(1,:) = 0.5*pF(1,:);
price = exp(-r * dt) * sum(real(pF) .* V) ;  % Option value at t_0

end


text:
function o = logBarrierRound(problem)
%Input: a structure P with the following fields
%  .Aeq
%  .beq
%  .lb
%  .ub
% describing the polytope {Aeq x = beq, lb <= x <= ub}
%Output:
% o - problem structure

A = problem.Aeq; b = problem.beq;
lb = problem.lb; ub = problem.ub;

f = ConvexProgram.LinearProgram(A, b, [], lb, ub);
f.normalize();
assert(f.feasible, 'The problem is not feasible.')

A = double(f.A);
b = double(f.b);
x = double(f.interior);
lb = double(f.barrier.lb);
ub = double(f.barrier.ub);

x0 = f.x0;
idx = f.idx;
scale = double(f.scale);

o = struct;
o.P = struct('A', A, 'b', b, 'lb', lb, 'ub', ub, 'x', x); % rounded polytope with feasible point x
o.T = struct('x0', x0, 'idx', idx, 'scale', scale); % used to recover sample in original space by x(idx) = x0(idx) + scale.*samples
end

author: opencobra

repo: cobratoolbox


text:
function mc = moc ( a, b, n, f )

%*****************************************************************************80
%
%% MOC estimates the modulus of continuity of a function over an interval.
%
%  Discussion;
%
%    The modulus of continuity function MC(T) for a function F(X) over an
%    interval [A,B] is defined as
%
%      MC(T) = max | F(X+DX) - F(X) | for 0 <= DX <= T, and X and X+DX in [A,B].
%
%    The modulus of continuity function is a monotone increasing function,
%    with MC(0) = 0.
%
%    This function estimates the modulus of continuity based on a discrete
%    set of data at N equally spaced points in the interval [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the left and right endpoints of the interval.
%
%    Input, integer N, the number of equally spaced sample points.
%
%    Input, function F(x), a handle to the function.
%
%    Output, real MC(N), the modulus of continuity function estimated at
%    0, H or less, 2*H or less, ..., (N-1)*H or less.
%

%  Compute the maximum difference with a separation DX of exactly 0*H, 1*H, 2*H,
%  ..., (N-1)*H.
%
  mc1 = moc1 ( a, b, n, f );
%
%  Compute the maximum difference with a separation DX of 0*H or less,
%  1*H or less, 2*H or less, ..., (N-1)*H or less.
%
  mc = zeros ( n, 1 );
  for i = 1 : n - 1
    mc(i+1) = max ( mc(i), mc1(i+1) );
  end

  return
end


text:
function [M_R] = vox2ras_rsolveAA(Vc_C, inPlaneRotation, varargin)
%%
%% NAME
%%
%%     vox2ras_rsolveAA.m (vox2ras_r{otation}solveA{uto}A{lign})
%%
%% AUTHOR
%%
%%	Rudolph Pienaar
%%
%% SYNOPSIS
%%
%%     [M_R] = vox2ras_rsolveAA(Vc_C, inPlaneRotation, ch_orientation)
%%
%% ARGUMENTS
%%
%%      Vc_C		in      column vector defining a direction cosine for
%%					a volume
%%	inPlaneRotation	in	scalar float defining the in plane rotation as
%%					read from the data meas.asc file
%%	ch_orientation	in	optional character string definining the plane of
%%					the main slab orientation
%%	M_R		out	vox2ras rotational candidate
%%
%% DESCRIPTION
%%
%%	"vox2ras_rsolveAA" attempts to find a candidate vox2ras rotation matrix using
%%	techniques described in autoaligncorrect.cpp.
%%
%%	The optional character string, ch_orientation, defaults to 's' (sagittal)
%%	and defines the orientation of the main slab. This string is one of:
%%
%%		's'	sagittal
%%		'c'	coronal
%%		't'	transverse
%%
%% PRECONDITIONS
%%
%%	o The vector C is read from a Siemens meas.asc file such that
%%		ci	= sSliceArray.asSlice[0].sNormal.dSag
%%		cj	= sSliceArray.asSlice[0].sNormal.dCor
%%		ck	= sSliceArray.asSlice[0].sNormal.dTra
%%
%% POSTCONDITIONS
%%
%%	o All returned matrices are 4x4.
%%	o Only the rotations of the vox2ras matrix are determined by this function.
%%		The center of k-space is not determined.
%%
%% SEE ALSO
%%
%%	vox2ras_rsolve	- determine the rotational component of a vox2ras matrix
%%				using Siemens reference orientations indirectly
%%	vox2ras_ksolve	- determine the k-space col in RAS of a vox2ras matrix
%%	vox2ras_dfmeas	- main function: determines the vox2ras matrix from a
%%			  Siemens meas.asc file.
%%
%% HISTORY
%%
%% 02 June 2004
%% o Initial design and coding.
%%

%
% vox2ras_rsolveAA.m
%
% Original Author: Rudolph Pienaar
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

M_R		= zeros(4, 4);
Vc_Cn		= Vc_C./(norm(Vc_C));
%Vc_Cn		= Vc_C;
ch_orientation	= 's';
if length(varargin)
	ch_orientation	= varargin{1}(1);
end

%% phase reference vector -
%%	adapted from Andre van der Kouwe's "autoaligncorrect.cpp"
Vc_P	= zeros(3, 1);
switch ch_orientation
    case 't'
    	Vc_P(1)	= 0;
	Vc_P(2)	=  Vc_Cn(3)*sqrt(1/(Vc_Cn(2)*Vc_Cn(2)+Vc_Cn(3)*Vc_Cn(3)));
	Vc_P(3)	= -Vc_Cn(2)*sqrt(1/(Vc_Cn(2)*Vc_Cn(2)+Vc_Cn(3)*Vc_Cn(3)));
    case 'c'
	Vc_P(1)	=  Vc_Cn(2)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
	Vc_P(2)	= -Vc_Cn(1)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
    	Vc_P(3)	= 0;
    case 's'
	Vc_P(1)	= -Vc_Cn(2)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
	Vc_P(2)	=  Vc_Cn(1)*sqrt(1/(Vc_Cn(1)*Vc_Cn(1)+Vc_Cn(2)*Vc_Cn(2)));
    	Vc_P(3)	= 0;
    otherwise
        fprintf(1, 'Unknown orientation parameter passed. Returning with dummy M_R');
	return;
end

%% The readout reference vector is the cross product of Vc_Cn and Vc_P
Vc_R	= cross(Vc_Cn, Vc_P);

M_R(1:3, 1)	= Vc_P;
M_R(1:3, 2)	= Vc_R;
M_R(1:3, 3)	= Vc_C;

%% The above calculated rotation matrices define an (x,y) plane
%%	given by the first two column vectors. These reference
%%	vectors need to be rotated by inPlaneRotation to
%%	reach the final rotational vox2ras.

theta_f		= inPlaneRotation;

M3_Mu	= [	 cos(theta_f)	 sin(theta_f)	0
		-sin(theta_f)	 cos(theta_f)	0
		 	0		0	1];
M3_R	= M_R(1:3,1:3)	* M3_Mu;
M_R(1:3,1:3)	= M3_R;

%% The MGH vox2ras matrix inverts the Readout column
M_R		= M_R * [ 1  0  0  0
			  0 -1  0  0
			  0  0  1  0
			  0  0  0  1];

%% All done!


text:
function [H] = rigidbody(f);

% RIGIDBODY creates the homogenous spatial transformation matrix
% for a 6 parameter rigid-body transformation
%
% Use as
%   [H] = rigidbody(f)
%
% The transformation vector f should contain the
%   x-shift
%   y-shift
%   z-shift
% followed by the
%   pitch (rotation around x-axis)
%   roll  (rotation around y-axis)
%   yaw   (rotation around z-axis)

% Copyright (C) 2000-2005, Robert Oostenveld
%
% 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

% $Log: rigidbody.m,v $
% Revision 1.1  2009/01/30 04:02:11  arno
% *** empty log message ***
%
% Revision 1.4  2006/04/13 10:37:38  roboos
% added a ; to the end of a line
%
% Revision 1.3  2005/08/15 08:15:32  roboos
% reimplemented the rotate function, which contained an error (the error is in the AIR technical reference)
% changed all functions to be dependent on the rotate, translate and scale function
% all functions now behave consistenly, which also means that they are not compleetly backward compatible w.r.t. the order of the rotations
%
% Revision 1.2  2004/05/19 09:57:07  roberto
% added GPL copyright statement, added CVS log item
%

% compute the homogenous transformation matrix for the translation
T = translate(f([1 2 3]));

% compute the homogenous transformation matrix for the rotation
R = rotate_new(f([4 5 6]));

% compute the homogenous transformation matrix for the combination
H = T*R;

repo: REST


text:
function value = givens_determinant ( n )

%*****************************************************************************80
%
%% GIVENS_DETERMINANT returns the determinant of the GIVENS matrix.
%
%  Discussion:
%
%    Since a formula for the eigenvalues is known, we compute the
%    determinant as the product of those values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real VALUE, the determinant.
%
  value = 1.0;

  for i = 1 : n
    angle = ( 2 * i - 1 ) * pi / ( 4 * n );
    value = value * 0.5 / ( cos ( angle ) )^2;
  end

  return
end


text:
function net=initWeight(net,opts)
num=numel(net.layers);
for layer=1:num
    if(strcmp(net.layers{layer}.type,'conv'))
        [h,w,in,out]=size(net.layers{layer}.weights{1});
        net.layers{layer}.weights={init_weight(opts, h, w, in, out, 'single'),ones(out, 1, 'single')*opts.initBias};
    end
end
end


text:
function test_pull1602

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_specest_irasa ft_freqanalysis

% simulate data
time = (1:2000)/1000;
for rpt = 1:20
  if false
    % We do not have that many concurrent licenses of the DSP Systems toolbox (https://nl.mathworks.com/products/dsp-system.html)
    % which was used in the initial version of this script. Therefore this test script does not use the DSP toolbox by default.
    dspobj = dsp.ColoredNoise('Color', 'pink', 'SamplesPerFrame', length(time));
    signal = dspobj()';
  else
    % integrate over time to create pink noise with 1/f power spectrum
    signal = cumsum(randn(1,length(time)));
  end

  % add line noise
  data.trial{1,rpt}     = signal + cos(2*pi*50*time) + cos(2*pi*100*time);
  data.time{1,rpt}      = time;
  data.label{1}         = 'chan';
  data.trialinfo(rpt,1) = rpt;
end

% using unfiltered data
cfg = [];
cfg.method  = 'irasa';
cfg.output  = 'original';
cfg.pad     = 'nextpow2';
freq = ft_freqanalysis(cfg, data);

cfg.output = 'fractal';
freqI = ft_freqanalysis(cfg, data);

% what happens if we use bandpassfiltered data?
cfg2            = [];
cfg2.bpfilter   = 'yes';
cfg2.bpfilttype = 'firws';
cfg2.bpfreq     = [60 150];
datafilt1 = ft_preprocessing(cfg2, data);

cfg2.bpfilter   = 'no';
cfg2.hpfilter   = 'yes';
cfg2.hpfreq     = 60;
cfg2.hpfilttype = 'firws';
datafilt2 = ft_preprocessing(cfg2, data);

cfg2.bpfilter   = 'no';
cfg2.hpfilter   = 'no';
cfg2.dftfilter  = 'yes'; % notch filter to filter out line noise
cfg2.dftfreq    = [50 100];
datafilt3 = ft_preprocessing(cfg2, data);

cfg.output = 'original';
freqfilt1  = ft_freqanalysis(cfg, datafilt1);
freqfilt2 = ft_freqanalysis(cfg, datafilt2);
freqfilt3 = ft_freqanalysis(cfg, datafilt3);

cfg.output = 'fractal';
freqfiltI = ft_freqanalysis(cfg, datafilt1);
freqfiltI2 = ft_freqanalysis(cfg, datafilt2);
freqfiltI3 = ft_freqanalysis(cfg, datafilt3);

allpow = [
  freq.powspctrm
  freqfilt1.powspctrm
  freqfilt2.powspctrm
  freqfilt3.powspctrm
  freqI.powspctrm
  freqfiltI.powspctrm
  freqfiltI2.powspctrm
  freqfiltI3.powspctrm
  ];

figure;
semilogy(freq.freq, allpow);
legend({
  'orig-unfilterd'
  'orig-bpfiltered'
  'orig-hpfiltered'
  'orig-dftfiltered'
  'frac-unfiltered'
  'frac-bpfiltered'
  'frac-hpfiltered'
  'frac-dftfiltered'
  });


text:
function y = daub18_transform ( n, x )

%*****************************************************************************80
%
%% DAUB18_TRANSFORM computes the DAUB18 transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%    N must be a power of 2 and at least 4.
%
%    Input, real X(N), the vector to be transformed.
%
%    Output, real Y(N), the transformed vector.
%
  c = [ ...
     3.807794736387834E-02; ...
     2.438346746125903E-01; ...
     6.048231236901111E-01; ...
     6.572880780513005E-01; ...
     1.331973858250075E-01; ...
    -2.932737832791749E-01; ...
    -9.684078322297646E-02; ...
     1.485407493381063E-01; ...
     3.072568147933337E-02; ...
    -6.763282906132997E-02; ...
     2.509471148314519E-04; ...
     2.236166212367909E-02; ...
    -4.723204757751397E-03; ...
    -4.281503682463429E-03; ...
     1.847646883056226E-03; ...
     2.303857635231959E-04; ...
    -2.519631889427101E-04; ...
     3.934732031627159E-05 ];
  p = 17;
  y(1:n,1) = x(1:n);
  m = n;
  q = floor ( ( p - 1 ) / 2 );

  while ( 4 <= m )

    i = 1;
    z(1:m,1) = 0.0;

    for j = 1 : 2 : m - 1

      mh = floor ( m / 2 );
      for k = 0 : 2 : p - 1
        j0 = i4_wrap ( j + k,     1, m );
        j1 = i4_wrap ( j + k + 1, 1, m );
        z(i,1)    = z(i,1)    + c(  k+1) * y(j0) + c(  k+2) * y(j1);
        z(i+mh,1) = z(i+mh,1) + c(p-k+1) * y(j0) - c(p-k  ) * y(j1);
      end

      i = i + 1;

    end

    y(1:m,1) = z(1:m);

    m = floor ( m / 2 );

  end

  return
end


text:
%% Load the data into MATLAB from a binary log file
% Usage: >> [datapoints, numpoints] = readdata('datafile.log')
% Header information format:
%           String "MWLOGV##"
%           Time/Date 4 bytes (time())
%           Number of Signals per record Logged 1 bytes (256 max)
%           Data Type of Signals Logged  1 bytes (1-10)
%           Number of bytes per record 2 (65535 max)
% Plot Data Example: plot([1:numpoints], datapoints(1,:), [1:numpoints], datapoints(2,:))
% MathWorks Pilot Engineering 2015
% Steve Kuznicki
function [datapts, numpts] = px4_read_binary_file(dataFile)
%%
datapts = 0;
numpts = 0;

if nargin == 0
    dataFile = 'data.bin';
end

fid = fopen(dataFile, 'r');
% load the header information
hdrToken = fread(fid, 8, 'char');
if strncmp(char(hdrToken),'MWLOGV',6) == true
    logTime = uint32(fread(fid, 1, 'uint32'));
    numflds = double(fread(fid, 1, 'uint8'));
    typefld = uint8(fread(fid, 1, 'uint8'));
    recSize = uint16(fread(fid, 1, 'uint16'));
    fieldTypeStr = get_elem_type(typefld);
    datapts = fread(fid, double([numflds, Inf]), fieldTypeStr);
    fclose(fid);
    numpts = size(datapts,2);
end

end

%% get the element type string
function [dtypeStr] = get_elem_type(dtype)
    switch(dtype)
        case 1
            dtypeStr = 'double';
        case 2
            dtypeStr = 'single';
        case 3
            dtypeStr = 'int32';
        case 4
            dtypeStr = 'uint32';
        case 5
            dtypeStr = 'int16';
        case 6
            dtypeStr = 'uint16';
        case 7
            dtypeStr = 'int8';
        case 8
            dtypeStr = 'uint8';
        case 9
            dtypeStr = 'logical';
        case 10
            dtypeStr = 'embedded.fi';
    end
end

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [isSessionInfo] = bz_isSessionInfo(SessionInfo)
% USAGE
% [isSessionInfo] = bz_isSessionInfo(SessionInfo)
%
% INPUT
%       SessionInfo   - struct with the following fields
%                   .channels
%                   .region
%                   .depth
%                   .nChannels
%                   .FileName
%                   .Date
%
%
% OUTPUT
%      logical true if struct meets event criteria, false if otherwise
%
% written by david tingley, 2017

if isfield(SessionInfo,'channels')
     if isvector(SessionInfo.channels)
         isSessionInfo = true;
     else
         warning('.channels field is not formatted correctly')
         isSessionInfo = false;
     end
end
if isfield(SessionInfo,'region')
     if iscell(SessionInfo.region)
         isSessionInfo = true;
     else
         warning('.region field is not formatted correctly')
         isSessionInfo = false;
     end
end
if isfield(SessionInfo,'depth')
     if isnumeric(SessionInfo.depth)
         isSessionInfo = true;
     else
         warning('.depth field is not formatted correctly')
         isSessionInfo = false;
     end
end
if isfield(SessionInfo,'spikeGroups')
     if isstruct(SessionInfo.spikeGroups) % should we go deeper here to check the struct?
         isSessionInfo = true;
     else
         warning('.spikeGroups field is not formatted correctly')
         isSessionInfo = false;
     end
end
if isfield(SessionInfo,'FileName')
     if ischar(SessionInfo.FileName)
         isSessionInfo = true;
     else
         warning('.FileName field is not formatted correctly')
         isSessionInfo = false;
     end
end
if isfield(SessionInfo,'Date')
     if ischar(SessionInfo.Date)
         isSessionInfo = true;
     else
         warning('.Date field is not formatted correctly')
         isSessionInfo = false;
     end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%%=====================================================================
%% $RCSfile: tps_test_script.m,v $
%% $Author: bjian $
%% $Date: 2008/11/24 08:59:03 $
%% $Revision: 1.1 $
%%=====================================================================

% preparation: set ctrl pts
ctrl_pts = rand(25,2);
[PP,kernel] = tps_set_ctrl_pts(ctrl_pts);
lambda = 0;

m = 10;
n = 10;
err = ones(m,n);
for i=1:m
    % step1: set landmarks;
    landmarks = rand(100,2);
    [U,Pm,Q1,Q2,R] = tps_set_landmarks(landmarks,ctrl_pts);

    for j=1:n
        % step2: set parameters;
        param0 = rand(25,2);
        [target, bending] = tps_warp(PP,kernel,U,Pm,param0);
        param1 = tps_compute_param(PP,kernel,U,Pm,Q1,Q2,R,lambda,target);
        err(i,j) = norm(param1-param0);
    end
end
max(err(:))


text:
function plotRegionOverlapScript

load('~/data/occlusion/labelme/results/result_LM_region_mean_covering.mat');

%% Plot by area
area_vals = [0.0025 0.005 0.01 0.2 0.4 1];

areas = cat(2, cover_occ1.area)';

% assign index for largest area_vals that area is less than or equal to
area_ind = zeros(size(areas));
for k = 1:numel(area_vals)
    area_ind = area_ind + (areas > area_vals(k));
end
area_ind = area_ind+1;

names = {'Pb+ucm', 'globalPb+ucm', 'occ_1', 'occ_{final}'};
area_hist = zeros(numel(area_vals), 4);
ovmax{1} = cat(1, cover_ucm1.maxov);
ovmax{2} = cat(1, cover_ucm2.maxov);
ovmax{3} = cat(1, cover_occ1.maxov);
ovmax{4} = cat(1, cover_occave.maxov);

area_hist = zeros(numel(area_vals), numel(ovmax));
for k = 1:numel(ovmax)
    for k2 = 1:numel(area_vals)
        area_hist(k2, k) = mean(ovmax{k}(area_ind==k2));
    end
end

na = numel(area_vals);
nl = numel(ovmax);
figure(1), hold off;
style = {'--', '--', '-', '-' ,'-'};
markers = {'+', 'x', 'o', '+', 'x'};
for k = 1:numel(ovmax)
    plot(area_hist(:, k), 'LineStyle', style{k}, 'MarkerSize', 12, 'Marker', markers{k}, 'LineWidth', mod(k,3)+2, 'Color', hsv2rgb(k/nl, 1, k/nl)); hold on;
end
xlabels = strtokAll(num2str(area_vals), ' ');
xlabels = [[{'0'} ; xlabels(1:end-1)] repmat({'-'}, na, 1)   xlabels(:) [repmat({'|'}, na-1, 1) ; {' '}]]';
axis([1 6 0 1]); set(gca, 'XTick', 1:na, 'XTickLabel', cat(2, xlabels{:}), 'FontSize', 16);
legend(names, 'Location', 'NorthWest')
xlabel('Region Area', 'FontSize', 18); ylabel('Average Overlap', 'FontSize', 18)
print -f1 -depsc ~/data/occlusion/labelme/figs/plots/overlap_comparison_by_area.eps

%% Cumulative overlap
figure(2), hold off,
ind = true(size(areas));
for k = 1:numel(ovmax)
    y = sort(ovmax{k}(ind), 'descend');
    plot((1:numel(y))/numel(y), y, 'LineStyle', style{k}, 'LineWidth', mod(k,3)+2, 'Color', hsv2rgb(k/nl, 1, k/nl)); hold on;
end
axis([0 1 0 1]); set(gca, 'XTick', 0:0.1:1, 'FontSize', 16);
legend(names, 'Location', 'NorthEast')
xlabel('Recall', 'FontSize', 18); ylabel('Overlap', 'FontSize', 18)
print -f2 -depsc ~/data/occlusion/labelme/figs/plots/overlap_overall_comparison.eps

%% Cumulative overlap by area
figure(3), hold off,
for k = 1:na-1
    ind = areas>=area_vals(k);
    y = sort(ovmax{end}(ind), 'descend');
    plot((1:numel(y))/numel(y),y, 'LineStyle', style{k}, 'LineWidth', mod(k,3)+2, 'Color', hsv2rgb(k/(na-1), 1, k/(na-1))); hold on;
end
axis([0 1 0 1]); set(gca, 'XTick', 0:0.1:1, 'FontSize', 16);
for k = 1:na-1, lstr{k} = [' >' num2str(area_vals(k))]; end
legend(lstr, 'Location', 'SouthWest')
xlabel('Recall', 'FontSize', 18); ylabel('Overlap', 'FontSize', 18);
print -f3 -depsc  ~/data/occlusion/labelme/figs/plots/overlap_occlusion.eps

author: Cloud-CV

repo: object-proposals


text:
%  name for the perturbation parameter entering the structural model
%  explicitly
%  if parameterize is true, return also the default (MAIH's perturbation)
%  values
%

author: jmaih


text:
function varargout = cvx_run_solver( sfunc, varargin )
global cvx___
settings_arg = varargin{end};
settings     = varargin{end-1};
inputs       = varargin(1:end-nargout-2);
dumpfile = '';
custom_on = false;
if isstruct( settings ),
    for f = fieldnames( settings )',
        sval = settings.(f{1});
        if isequal( f{1}, 'dumpfile' ),
            dumpfile = sval;
        else
            custom_on = true;
            inputs{settings_arg}.(f{1}) = sval;
        end
    end
end
if custom_on,
    fprintf( 'NOTE: custom settings have been set for this solver.\n' );
end
if ~isempty( dumpfile ),
    if ~ischar( dumpfile ) || size( dumpfile, 1 ) > 1,
        cvx_throw( 'Invalid filename for the dumpfile.' );
    elseif length(dumpfile) < 4 || ~strcmpi(dumpfile(end-3:end),'.mat'),
        dumpfile = [ dumpfile, '.mat' ];
    end
    fprintf( 'Saving output to: %s\n', dumpfile );
    fprintf( '------------------------------------------------------------\n');
    inp_names = cell(1,length(inputs));
    for k = 1 : length(inp_names),
        inp_names{1,k} = inputname(k+1);
    end
    dstruct = cell2struct( inputs, inp_names, 2 ); %#ok
    save( dumpfile, '-struct', 'dstruct' );
    diaryfile = [ dumpfile, '.txt' ];
    fid = fopen( diaryfile, 'w+' );
    if fid ~= 0,
        fclose( fid );
        diary( diaryfile );
    end
elseif custom_on,
    fprintf( '------------------------------------------------------------\n');
end
errmsg = [];
if cvx___.isoctave,
    fflush(1);
end
try
   [ varargout{1:nargout} ] = sfunc( inputs{:} );
catch errmsg
   [ varargout{1:nargout} ] = deal( [] );
end
if ~isempty( dumpfile ),
    if fid ~= 0,
        diary( 'off' );
        fid = fopen( diaryfile, 'r' );
        if fid ~= 0,
            output = fread( fid, Inf, '*char' )';
            fclose( fid );
            delete( diaryfile );
        end
    end
    if fid == 0,
        output = '<Could not save>';
    end
    otp_names = varargin(end-nargout-1:end-2);
    dstruct = cell2struct( [ inputs, varargout, output ], [ inp_names, otp_names, 'output' ], 2 ); %#ok
    save( dumpfile, '-struct', 'dstruct' );
end
if ~isempty( errmsg ),
    rethrow( errmsg );
end

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

author: yu-jiang

repo: radpbook

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


text:
function setup_WSDDN()
%SETUP_WSDDN Sets up WSDDN, by adding its folders to the Matlab path

root = fileparts(mfilename('fullpath')) ;
addpath(root, [root '/matlab'], [root '/pascal'], [root '/core']) ;
addpath([vl_rootnn '/examples/']) ;
addpath([vl_rootnn '/examples/imagenet/']) ;

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function [sliceNum, planC] = getNoseSlice(outerStrMask3M,planC,outerStrName)
% Automatically identify nose slice in H&N CT scans based on
% first non-zero row of mask.
%
% AI 10/1/19
%
%------------------------------------------------------------------------
% INPUT
% outerStrMask3M   : Mask of pt outline. Set to [] to use structure name
%                    instead.
% planC
% outerStrName     : Structure name corresponding to pt outline
%------------------------------------------------------------------------

if ~isempty(outerStrMask3M)
    mask3M = outerStrMask3M;
else
    %Get mask of outer structure
    indexS = planC{end};
    strC = {planC{indexS.structures}.structureName};
    strIdx = getMatchingIndex(outerStrName,strC,'exact');
    [mask3M, planC] = getStrMask(strIdx, planC);
end

startSliceIdx = 11;
endSliceIdx = min(size(mask3M,3),100);
supMask3M = double(mask3M(:,:,startSliceIdx:endSliceIdx));

%Identify first non-zero row
[sel,rowIdxM] = max(supMask3M, [], 1);
rowIdxM = squeeze(sel.*(rowIdxM)).';
rowIdxM(rowIdxM==0) = nan;
minRowV = nanmin(rowIdxM,[],2);

%Smooth
minRowV = movmean(minRowV,5,'omitnan');

%Compute difference & identify min
%[~,mins] = min(diff([NaN;minRowV]));
[~,mins] = findpeaks(-minRowV,'MinPeakWidth',2, 'MaxPeakWidth',40);
if ~isempty(mins)
    sliceNum = mins(1) + startSliceIdx -1;
else
    sliceNum=1; %default
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function kern = whitefixedKernParamInit(kern)

% WHITEFIXEDKERNPARAMINIT WHITEFIXED kernel parameter initialisation.
% The white noise kernel arises from assuming independent Gaussian
% noise for each point in the function. The variance of the noise is
% given by the kern.variance parameter. The fixed white noise kernel
% is a simple variant of the white kernel that doesn't allow the noise
% parameter to be optimised. It is useful when the level of noise is
% known a priori.
%
% This kernel is not intended to be used independently, it is provided
% so that it may be combined with other kernels in a compound
% kernel.
%
% SEEALSO : cmpndKernParamInit
%
% FORMAT
% DESC initialises the fixed parameter white noise
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Nathaniel J. King, 2006

% KERN

kern.variance = exp(-2);
kern.nParams = 0;

kern.isStationary = true;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
C= 2;
T= 10000;
cnt= struct('fs', 100);
cnt.x= 50*randn(T, C);
cnt.clab= str_cprintf('Ch%d', 1:C);

M= 100;
mrk= struct; %('fs', cnt.fs);
mrk.time= round(linspace(0, T/cnt.fs*1000, M+2));
mrk.time([1 end])= [];
%mrk.desc= str_cprintf('S%3d', ceil(rand(1,M)*10))';
mrk.desc= ceil(rand(1,M)*10);

% --- Setup a very simple system for (simulated) online processing
bbci= struct;
bbci.source.acquire_fcn= @bbci_acquire_offline;
bbci.source.acquire_param= {cnt, mrk};
bbci.source.log.output= 'screen';
data= bbci_recordSignals(bbci, '/tmp/rec_test');

[cnt_re, mrk_re]= file_readBV(data.source.record.filename);

isequal(cnt.clab, cnt_re.clab)
max(abs(cnt.x(:)-cnt_re.x(:)))

isequal(mrk.pos, mrk_re.pos(2:end))
isequal(mrk.desc, mrk_re.desc(2:end))

author: bbci


text:
% distance to wall
function [data,units] = compute_distnose2wall_animaldir_rect(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  xnose = trx(fly).x_mm + 2*trx(fly).a_mm.*cos(trx(fly).theta_mm);
  ynose = trx(fly).y_mm + 2*trx(fly).a_mm.*sin(trx(fly).theta_mm);
  nose = [xnose;ynose];
  body = [trx(fly).x_mm;trx(fly).y_mm ];
  tl = [trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly)];
  bl = [trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly)];
  tr = [trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly)];
  br = [trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly)];

  dtop = getDist(nose,body,tl,tr);
  dleft = getDist(nose,body,bl,tl);
  dright = getDist(nose,body,tr,br);
  dbottom = getDist(nose,body,br,bl);

  data{i} = min([dtop;dleft; dright; dbottom],[],1);

end
units = parseunits('mm');

function d = getDist(a,b,u,v)

x1=a(1,:);
x2=b(1,:);
x3=u(1);
x4=v(1);
y1=a(2,:);
y2=b(2,:);
y3=u(2);
y4=v(2);

ua = ((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))./((y4-y3)*(x2-x1)-(x4-x3)*(y2-y1));
x = x1 + ua.*(x2 - x1);
y = y1 + ua.*(y2 - y1);

side = sign((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3));
d = sqrt( (x-x1).^2 + (y-y1).^2);
d = side.*d;

% sign of ua tell us if the intersection point was away from the
% the body or not. ua>0 means that the intersection point was towards the
% tail rather than towards the head.
% side tells us whether the nose was beyond the wall or not.
% ua>0 & side>0 means that nose was inside the wall and the intersection
% point was behind the mice. In this case set the dist to 0.
d(ua>0 & side>0) = inf;

author: kristinbranson

repo: JAABA


text:
%s = rand('twister');
%rand('twister', s);

len = 200;              % number of 1st level observations (time points within subjects)
sub = 20;               % number of second-level units (subjects)

fixed_slope = 0.5;      % population fixed-effect slope
fixed_int = 0.3;        % pop. fixed-effect intercept

rand_slope = 0.5;         % 'true' random effect standard dev: slope
rand_int = 0.5;           % 'true' random effect standard dev: intercept

noise_sigma = 0.5;      % measurement error/unexplained std of 1st level observations
between_sigma = 0.5;    % measurement error/unexplained std of 2nd level covariate

% Generate simulated data
% 1) create predictors
x = zeros(len,sub);
y = x;
x(11:20,:) = 1;
x(111:120,:) = 1;

% 2) Create instances of effects for each subject
c = normrnd(fixed_slope, rand_slope, sub, 1);       % slope between-subjects variations
d = normrnd(fixed_int, rand_int, sub, 1);                  % intercept between-subjects variations

% 3) Create a between-ss covariate that is related to the slope but not the
% intercept
covt = scale(c + normrnd(0, between_sigma, sub, 1), 1);

% 4) Create a between-ss covariate that is related to the intercept but not the
% intercept
covti = scale(d + normrnd(0, between_sigma, sub, 1), 1);

corrcoef([c d covt covti])

% Add between-subjects error (random effects) and measurement noise
% (within-subjects error)

for i=1:sub
    y(:,i) = d(i) + c(i) .* x(:,i) + normrnd(0, noise_sigma, len, 1);
end

%% Run the model - without a 2nd-level covariate
out = igls(y, x);  % for igls
fprintf('\t\t%s\t%s\t%s\t\n', 'Effect', 'Intercept', 'Slope');
disp('True random-effect variances:');
fprintf('\t\t\t%3.3f\t%3.3f\n', [rand_slope rand_int]);

disp('Input random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', std([d c]))

disp('Est.  random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', sqrt(out.betastar)');

%% Run the model - with an unrelated 2nd-level covariate
out = igls(y, x, 'covariate', (1:20)');  % for igls
fprintf('\t\t%s\t%s\t%s\t\n', 'Effect', 'Intercept', 'Slope');
disp('True random-effect variances:');
fprintf('\t\t\t%3.3f\t%3.3f\n', [rand_slope rand_int]);

disp('Input random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', std([d c]))

disp('Est.  random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', sqrt(out.betastar)');

%% Run the model - with a 2nd-level covariate related to slope
out = igls(y, x, 'covariate', covt);  % for igls
fprintf('\t\t%s\t%s\t%s\t\n', 'Effect', 'Intercept', 'Slope');
disp('True random-effect variances:');
fprintf('\t\t\t%3.3f\t%3.3f\n', [rand_slope rand_int]);

disp('Input random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', std([d c]))

disp('Est.  random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', sqrt(out.betastar)');

%% Run the model - with a 2nd-level covariate related to intercept
out = igls(y, x, 'covariate', covti);  % for igls
fprintf('\t\t%s\t%s\t%s\t\n', 'Effect', 'Intercept', 'Slope');
disp('True random-effect variances:');
fprintf('\t\t\t%3.3f\t%3.3f\n', [rand_slope rand_int]);

disp('Input random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', std([d c]))

disp('Est.  random-effect variances: ');
fprintf('\t\t\t%3.3f\t%3.3f\n', sqrt(out.betastar)');

%%
clear *pvals* *fpr*

mytype = 'i';
iterations = 200;

len = 100; sub = 10;

fixpvals = zeros(iterations, 2);
randpvals = zeros(iterations, 2);
LRTrandpvals = zeros(iterations, 2);
chic = zeros(iterations, 1);
chid = zeros(iterations, 1);

c_fixed = 0;  % slope pop. average
d_fixed = 0;   % intercept pop. average
c_rand = 0;    % slope std across Ss
d_rand = 0;    % intercept std.
noise_std = 2.0;  % within-subjects noise std

x = zeros(len,sub);
y = x;
x(1:2:end,:) = 1;                   % create signal

% --------------------------------------------------------------------------
% This block: No true random effects, so test false positive rate for
% p_randvariance

verbstr = 'noverbose';

for i = 1:iterations

    c = normrnd(c_fixed,c_rand,sub,1);       % slope between-subjects variations
    d = normrnd(d_fixed,d_rand,sub,1);         % intercept between-subjects variations

    % Add between-subjects error (random effects) and measurement noise
    % (within-subjects error)
           y = zeros(len, sub);
    for s = 1:sub
        y(:,s) = d(s) + c(s) .* x(:,s) + normrnd(0,noise_std,len,1);
    end
    out = igls(y, x, 'type', mytype, verbstr);  % for igls

    fpr_converged(i) = out.isconverged;

    fpr_betastar(i, :) = out.betastar';

    fixpvals(i, :) = out.p;
    randpvals(i, :) = out.p_randvariance'; %[out.p_randvariance_d out.p_randvariance_c];

    LRTrandpvals(i, :) = out.pLRT_randvariance';

    if mod(i, 10) == 0, fprintf(1, '%3.0f ', i); end

    verbstr = 'noverbose';
end
fprintf('\n')

%pvals_fpr = pvals;
alph = .1;

fixfpr =  sum(fixpvals < alph) ./ iterations;
randfpr =  sum(randpvals < alph) ./ iterations;
LRTfpr =  sum(LRTrandpvals < alph) ./ iterations;

fprintf('\tfixed\t\trand fx\t\t\n');
fprintf('TPR/FPR @ alpha = %3.3f:\t%3.5f\t%3.5f\t%3.5f\t%3.5f\t%3.5f\t%3.5f\n', alph, fixfpr, randfpr, LRTfpr);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function value = grad(psi,co2)

N = psi.bandwidth;
A = psi.A(:);

k = 1:0.5:(N+1)/2;
k(2:2:end) = 0;
l = (N+1)/2:0.5:N;
l(end-1:-2:1)=0;
M = triu(hankel(k,l)-toeplitz([0,k(1:end-1)]));

Ahat = (-4) * (M*A(2:end));

Dpsi = SO3Kernel(Ahat);

% TODO: Interpolation error: Try
% psi = SO3DeLaValleePoussinKernel(90);
% psi2=SO3Kernel(psi.A);
% omega = 35*degree:0.01:pi;
% plot(omega/degree,psi.grad(cos(omega/2)) )
% hold on
% plot(omega/degree, psi2.grad(cos(omega/2)) )
% hold off

if nargin == 2
  value = Dpsi.eval(co2).*co2.*sqrt(1-co2.^2);
else
  value = SO3KernelHandle(@(co2) Dpsi.eval(co2).*co2.*sqrt(1-co2.^2));
end


text:
% This subroutine assigns creates a 3D grid with
% spacing dx,dy, dz (in degreees). The size will
% be selected interactiVELY. The pvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/98

report_this_filefun(mfilename('fullpath'));
global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dx = 0.1;
    dy = 0.1 ;
    dz = 5.00 ;
    ni = 300;
    R = 10000;

    def = {'0.1','0.1',num2str(dz),num2str(max(a.Depth)), num2str(min(a.Depth))};

    tit ='Three dimesional z-value analysis';
    prompt={...
        'Spacing in Longitude (dx in [deg])',...
        'Sapcing in Latitude  (dy in [deg])',...
        'Spacing in Depth    (dz in [km ])',...
        'Depth Range: deep limit [km] ',...
        'Depth Range: shallow limit',...
        };

    ni2 = inputdlg(prompt,tit,1,def);

    l = ni2{1}; dx= str2double(l);
    l = ni2{2}; dy= str2double(l);
    l = ni2{3}; dz= str2double(l);
    l = ni2{4}; z1= str2double(l);
    l = ni2{5}; z2= str2double(l);

    sel = 'ca'; zgrid3d

end   % if sel == 'in'

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'
    selgp3dB
    zvect=[z2:dz:z1];
    gz = zvect;
    itotal = length(t5);
    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    zvg = NaN(length(gx),length(gy),length(gz),300);
    ram  = NaN(length(gx),length(gy),length(gz),300);

    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name',' 3D gridding - percent done');;
    drawnow
    %
    %

    z0 = 0; x0 = 0; y0 = 0; dt = 1;
    % loop over all points
    for il =1:length(t5);

        x = t5(il,1);
        y = t5(il,2);
        z = t5(il,3);
        allcount = allcount + 1.;
        i2 = i2+1;

        % calculate distance from center point and sort wrt distance
        di = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2 + ((a.Depth - z)).^2 ) ;
        [s,is] = sort(di);

        l2 = find(is <= 300);

        %[cumu, xt] = hist(b(:,3),(t0b:(teb-t0b)/99:teb));

        zvg(t5(il,5),t5(il,6),t5(il,7),:) = is(1:300);
        ram(t5(il,5),t5(il,6),t5(il,7),:) = di(is(1:300));
        if rem(allcount,20) == 0;  waitbar(allcount/itotal) ;end
    end  % for xt5
    % save data
    %
    catSave3 =...
        [ 'zmap_message_center.set_info(''Save Grid'',''  '');think;',...
        '[file1,path1] = uiputfile(fullfile(hodi, ''eq_data'', ''*.mat''), ''Grid Datafile Name?'') ;',...
        ' sapa2 = [''save '' path1 file1 '' zvg ram gx gy gz dx dy dz  par1 tdiff t0b teb a main faults mainfault coastline yvect xvect tmpgri ll''];',...
        ' if length(file1) > 1, eval(sapa2),end , done']; eval(catSave3)

    close(wai)
    watchoff

    gz = -gz;
    zv2 = zvg;
    sel = 'no';
    lta_winy = 2;
    zv4 = zv2;
    tiz = 10;
    slm = 'new'; slicemapz;

end  % if cal

author: CelsoReyes


text:
%EPS2PDF  Convert an eps file to pdf format using ghostscript
%
% Examples:
%   eps2pdf source dest
%   eps2pdf(source, dest, crop)
%   eps2pdf(source, dest, crop, append)
%   eps2pdf(source, dest, crop, append, gray)
%   eps2pdf(source, dest, crop, append, gray, quality)
%
% This function converts an eps file to pdf format. The output can be
% optionally cropped and also converted to grayscale. If the output pdf
% file already exists then the eps file can optionally be appended as a new
% page on the end of the eps file. The level of bitmap compression can also
% optionally be set.
%
% This function requires that you have ghostscript installed on your
% system. Ghostscript can be downloaded from: http://www.ghostscript.com
%
%IN:
%   source - filename of the source eps file to convert. The filename is
%            assumed to already have the extension ".eps".
%   dest - filename of the destination pdf file. The filename is assumed to
%          already have the extension ".pdf".
%   crop - boolean indicating whether to crop the borders off the pdf.
%          Default: true.
%   append - boolean indicating whether the eps should be appended to the
%            end of the pdf as a new page (if the pdf exists already).
%            Default: false.
%   gray - boolean indicating whether the output pdf should be grayscale or
%          not. Default: false.
%   quality - scalar indicating the level of image bitmap quality to
%             output. A larger value gives a higher quality. quality > 100
%             gives lossless output. Default: ghostscript prepress default.

% Copyright (C) Oliver Woodford 2009-2011

% Suggestion of appending pdf files provided by Matt C at:
% http://www.mathworks.com/matlabcentral/fileexchange/23629

% Thank you to Fabio Viola for pointing out compression artifacts, leading
% to the quality setting.
% Thank you to Scott for pointing out the subsampling of very small images,
% which was fixed for lossless compression settings.

% 9/12/2011 Pass font path to ghostscript.

function eps2pdf(source, dest, crop, append, gray, quality)
% Intialise the options string for ghostscript
options = ['-q -dNOPAUSE -dBATCH -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress -sOutputFile="' dest '"'];
% Set crop option
if nargin < 3 || crop
    options = [options ' -dEPSCrop'];
end
% Set the font path
fp = font_path();
if ~isempty(fp)
    options = [options ' -sFONTPATH="' fp '"'];
end
% Set the grayscale option
if nargin > 4 && gray
    options = [options ' -sColorConversionStrategy=Gray -dProcessColorModel=/DeviceGray'];
end
% Set the bitmap quality
if nargin > 5 && ~isempty(quality)
    options = [options ' -dAutoFilterColorImages=false -dAutoFilterGrayImages=false'];
    if quality > 100
        options = [options ' -dColorImageFilter=/FlateEncode -dGrayImageFilter=/FlateEncode -c ".setpdfwrite << /ColorImageDownsampleThreshold 10 /GrayImageDownsampleThreshold 10 >> setdistillerparams"'];
    else
        options = [options ' -dColorImageFilter=/DCTEncode -dGrayImageFilter=/DCTEncode'];
        v = 1 + (quality < 80);
        quality = 1 - quality / 100;
        s = sprintf('<< /QFactor %.2f /Blend 1 /HSample [%d 1 1 %d] /VSample [%d 1 1 %d] >>', quality, v, v, v, v);
        options = sprintf('%s -c ".setpdfwrite << /ColorImageDict %s /GrayImageDict %s >> setdistillerparams"', options, s, s);
    end
end
% Check if the output file exists
if nargin > 3 && append && exist(dest, 'file') == 2
    % File exists - append current figure to the end
    tmp_nam = tempname;
    % Copy the file
    copyfile(dest, tmp_nam);
    % Add the output file names
    options = [options ' -f "' tmp_nam '" "' source '"'];
    try
        % Convert to pdf using ghostscript
        [status, message] = ghostscript(options);
    catch me
        % Delete the intermediate file
        delete(tmp_nam);
        rethrow(me);
    end
    % Delete the intermediate file
    delete(tmp_nam);
else
    % File doesn't exist or should be over-written
    % Add the output file names
    options = [options ' -f "' source '"'];
    % Convert to pdf using ghostscript
    [status, message] = ghostscript(options);
end
% Check for error
if status
    % Report error
    if isempty(message)
        error('Unable to generate pdf. Check destination directory is writable.');
    else
        error(message);
    end
end
return

% Function to return (and create, where necessary) the font path
function fp = font_path()
fp = user_string('gs_font_path');
if ~isempty(fp)
    return
end
% Create the path
% Start with the default path
fp = getenv('GS_FONTPATH');
% Add on the typical directories for a given OS
if ispc
    if ~isempty(fp)
        fp = [fp ';'];
    end
    fp = [fp getenv('WINDIR') filesep 'Fonts'];
else
    if ~isempty(fp)
        fp = [fp ':'];
    end
    fp = [fp '/usr/share/fonts:/usr/local/share/fonts:/usr/share/fonts/X11:/usr/local/share/fonts/X11:/usr/share/fonts/truetype:/usr/local/share/fonts/truetype'];
end
user_string('gs_font_path', fp);
return


text:
function vo = ch_qmulv(q, vi)
% 向量通过四元数做3D旋转
%
% Inputs: q - Qb2n
%            vi - 需要旋转的向量
% Output: vout - output vector, such that vout = q*vin*conjugation(q)
%
% See also  q2mat, qconj, qmul.

%     qi = [0; vi];
%     qo = ch_qmul(ch_qmul(q,qi),ch_qconj(q));
%     vo = qo(2:4,1);

    qo1 =              - q(2) * vi(1) - q(3) * vi(2) - q(4) * vi(3);
    qo2 = q(1) * vi(1)                + q(3) * vi(3) - q(4) * vi(2);
    qo3 = q(1) * vi(2)                + q(4) * vi(1) - q(2) * vi(3);
    qo4 = q(1) * vi(3)                + q(2) * vi(2) - q(3) * vi(1);
    vo = vi;
    vo(1) = -qo1 * q(2) + qo2 * q(1) - qo3 * q(4) + qo4 * q(3);
    vo(2) = -qo1 * q(3) + qo3 * q(1) - qo4 * q(2) + qo2 * q(4);
    vo(3) = -qo1 * q(4) + qo4 * q(1) - qo2 * q(3) + qo3 * q(2);

author: yandld


text:
function change(file)
[np,nf]=num(file);
fprintf('Number_of_point:%d\n',np);
fprintf('Number_of_face:%d\n',nf);
fid=fopen('OFF','w');
fprintf(fid,'OFF\n');
fprintf(fid,'%d %d 0\n',np,nf);
fid2=fopen(file);
H=0;
while 1
    tline=fgetl(fid2);
    H=H+1;
    if ~ischar(tline) , break , end
    if H > 1
        disp(tline)
        fprintf(fid,'%s\n',tline);
    end

end
fclose(fid2);
fclose(fid);
function [np,nf]=num(file)
fid=fopen(file);
C=0;
while 1
    tline=fgetl(fid);
    C=C+1;
    if C==1
        num=textscan(tline,'%s %d %d %d');
        break;
    end
     if ~ischar(tline) , break , end

end
fclose(fid);
np=num{1,2};
nf=num{1,3};

author: Sable

repo: mcbench-benchmarks