text:
function [XYZworld,faces, RGBcolor, tsdf_value, tsdf_weight] = TSDFmeshSupposedToBeFasterButNot(data, frameIdtarget,frameIds,cameraRtC2W, is_color)

if ~exist('is_color','var')
    is_color = true;
end

cameraRtC2Wframes = cameraRtC2W(:,:,frameIds);
frameIdtargetID = find(frameIds==frameIdtarget);

mapRt = transformCameraRt(cameraRtC2Wframes(:,:,frameIdtargetID));

for i=1:length(frameIds)
    cameraRtC2Wframes(:,:,i) = concatenateRts(mapRt, cameraRtC2Wframes(:,:,i));
end

%% tsdf

voxel.unit = 0.01; % Kevin: 4mm = 0.004 meter. Kinect cannot go better than 3mm
voxel.mu_grid = 10; % used to be 4
voxel.size_grid = [512; 512; 1024]; % [512; 512; 512];

voxel.range(1,1) = - voxel.size_grid(1) * voxel.unit / 2;
voxel.range(1,2) = voxel.range(1,1) + (voxel.size_grid(1)-1) * voxel.unit;

voxel.range(2,1) = - voxel.size_grid(2) * voxel.unit / 2;
voxel.range(2,2) = voxel.range(2,1) + (voxel.size_grid(2)-1) * voxel.unit;

voxel.range(3,1) = -0.5; % - voxel.size_grid(3) * voxel.unit / 2;
voxel.range(3,2) = voxel.range(3,1) + (voxel.size_grid(3)-1) * voxel.unit;

voxel.mu = voxel.mu_grid * voxel.unit;

fprintf('memory = %f GB\n',  prod(voxel.size_grid) * 4 / (1024*1024*1024));
fprintf('space = %.2f m x %.2f m x %.2f m ', voxel.size_grid(1) * voxel.unit, voxel.size_grid(2) * voxel.unit, voxel.size_grid(3) * voxel.unit);
fprintf('= [%.2f,%.2f] x [%.2f,%.2f] x [%.2f,%.2f]\n',voxel.range(1,1),voxel.range(1,2),voxel.range(2,1),voxel.range(2,2),voxel.range(3,1),voxel.range(3,2));

%tsdf_value  = -ones([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)],'single');
tsdf_value  =  ones([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)],'single');
tsdf_weight = zeros([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)],'single');

if is_color
    tsdf_color  = zeros([voxel.size_grid(1),voxel.size_grid(2), voxel.size_grid(3)], 'uint32');
end

f = data.K(1,1);

ViewFrustumC = [...
    0 -data.K(1,3) -data.K(1,3)  data.K(1,3)  data.K(1,3);
    0 -data.K(2,3)  data.K(2,3)  data.K(2,3) -data.K(2,3);
    0    f    f    f    f];
ViewFrustumC = ViewFrustumC/f * 8; % 8 meter is the furthest depth of kinect

% precompute
[pX,pY]=meshgrid(1:640,1:480);

raycastingDirectionC = [pX(:)'-data.K(1,3); pY(:)'-data.K(2,3); f*ones(1,640*480)]; % clipping at 8 meter is the furthest depth of kinect
raycastingDirectionC = raycastingDirectionC ./ repmat(sqrt(sum(raycastingDirectionC.^2,1)),3,1);

%% get the grid there

underConsideration = true(1,prod(voxel.size_grid));

%% speed up by eliminate free space

% if you want to debug fast, enough some frames only
firstFrame = true;
%for i=1:33:length(frameIds)
for i=1:length(frameIds)

    fprintf('======================================== frame %d ========================================\n',i);

    camRtC2W = cameraRtC2Wframes(:,:,i);

    IMGcam = readImage(data.image{frameIds(i)});
    XYZcam = depth2XYZcamera(data.K, depthRead(data.depth{frameIds(i)}));

    ViewFrustumW = transformRT(ViewFrustumC, camRtC2W);

    range2test = [min(ViewFrustumW, [], 2) max(ViewFrustumW, [], 2)];

    % choose a bounding box to contain viewing frustum
    rangeGrid = (range2test - voxel.range(:,[1 1])) / voxel.unit + 1;
    rangeGrid(:,1) = max(1,floor(rangeGrid(:,1)));
    rangeGrid(:,2) = min(ceil (rangeGrid(:,2)),voxel.size_grid);
    rangeGrid = int32(rangeGrid);

    % get the grid there
    disp('meshgrid');
    tic;
    [Y,X,Z]=meshgrid(rangeGrid(1,1):rangeGrid(1,2),rangeGrid(2,1):rangeGrid(2,2),rangeGrid(3,1):rangeGrid(3,2)); % strange matlab syntax
    toc;
    disp('index generation');
    tic;
    X = X(:)'; Y = Y(:)'; Z = Z(:)';
    gridIndex = sub2ind(voxel.size_grid',X,Y,Z);
    toc;

    % speed up
    if firstFrame
        firstFrame = false;
    else
        disp('select voxels worth consideration');
        tic;
        underConsiderationHere = underConsideration(gridIndex);
        X = X(underConsiderationHere);
        Y = Y(underConsiderationHere);
        Z = Z(underConsiderationHere);
        gridIndex = gridIndex(underConsiderationHere);
        toc;
    end

    fprintf('%d voxels to consider\n',length(gridIndex));

    disp('world coordinate');
    tic;
    gridCoordinateW = [single(X)*voxel.unit + voxel.range(1,1); single(Y)*voxel.unit + voxel.range(2,1); single(Z)*voxel.unit + voxel.range(3,1)];
    clear X Y Z;
    toc;

    % transform the grid
    disp('transform');
    tic;
    gridCoordinateC = transformRT(gridCoordinateW, camRtC2W, true);
    toc;

    % select: in front of camera
    disp('select: in front of camera');
    tic;
    isValid = gridCoordinateC(3,:)>0; % find(gridCoordinateC(3,:)>0);
    gridCoordinateC = gridCoordinateC(:,isValid);
    gridIndex = gridIndex(isValid);
    toc;

    % select: project
    disp('select: project');
    tic;
    px = round(data.K(1,1)*(gridCoordinateC(1,:)./gridCoordinateC(3,:)) + data.K(1,3));
    py = round(data.K(2,2)*(gridCoordinateC(2,:)./gridCoordinateC(3,:)) + data.K(2,3));
    isValid = (1<=px & px <= 640 & 1<=py & py<= 480);
    gridCoordinateC = gridCoordinateC(:,isValid);
    gridIndex = gridIndex(isValid);
    py = py(isValid);
    px = px(isValid);
    toc;

    % select: valid depth
    disp('select: valid depth');
    tic;
    ind = sub2ind([480 640],py,px);
    isValid = XYZcam(ind+640*480*3)~=0;
    gridCoordinateC = gridCoordinateC(:,isValid);
    gridIndex = gridIndex(isValid);
    ind = ind(isValid);
    toc;

    % compare distance between measurement and the grid
    disp('compare distance between measurement and the grid');
    %{

            const float diff = (depth[px] - cameraX.z) * sqrt(1+sq(pos.x/pos.z) + sq(pos.y/pos.z));
            if(diff > -mu){
                const float sdf = fminf(1.f, diff/mu);
                float2 data = vol[pix];
                data.x = clamp((data.y*data.x + sdf)/(data.y + 1), -1.f, 1.f);
                data.y = fminf(data.y+1, maxweight);
                vol.set(pix, data);
            }

    %}
    tic;
    eta = (XYZcam(ind+640*480*2)- gridCoordinateC(3,:)) .* ((1+ (gridCoordinateC(1,:)./gridCoordinateC(3,:)).^2 + (gridCoordinateC(2,:)./gridCoordinateC(3,:)).^2 ).^0.5);
    toc;

    %tic;
    %XYZmeasure = [XYZcam(ind);XYZcam(ind+640*480);XYZcam(ind+640*480*2)];
    %Dmeasure = sqrt(sum(XYZmeasure.^2,1));
    %Dtran    = sqrt(sum(gridCoordinateC.^2,1));
    %eta = Dmeasure-Dtran;
    %toc;

    %
    %disp('speed up');
    %isValid = eta>voxel.mu*2;
    %underConsideration(gridIndex(isValid)) = false;

    % select: > - mu
    disp('select: > - mu');
    tic;
    isValid = eta>-voxel.mu;
    eta = eta(isValid);
    gridIndex = gridIndex(isValid);
    ind = ind(isValid);
    toc;

    new_value = min(1,eta/voxel.mu);

    disp('read write tsdf');
    tic;
    old_weight = tsdf_weight(gridIndex);
    new_weight = old_weight + 1;
    tsdf_weight (gridIndex)= new_weight;
    tsdf_value (gridIndex) = (tsdf_value(gridIndex).*old_weight +new_value)./new_weight;
    toc;

    % update speed up vector
    disp('updating speed up vector'); tic;
    underConsideration(gridIndex(eta>voxel.mu*2)) = false;
    toc;

    if is_color
        disp('coloring');
        tic;
        isValid = eta<=voxel.mu;
        eta = eta(isValid);
        gridIndex = gridIndex(isValid);
        ind = ind(isValid);

        IMGmeasure = [IMGcam(ind);IMGcam(ind+640*480);IMGcam(ind+640*480*2)]*255;

        colorArray = reshape(typecast(tsdf_color(gridIndex),'uint8'),4,[]);

        old_weight = colorArray(4,:);
        old_weight(double(old_weight)>=255) = 254;
        new_weight = old_weight + 1;

        newArray = uint8(round((double(colorArray(1:3,:)) .* repmat(double(old_weight),3,1) + double(IMGmeasure)) ./ repmat(double(new_weight),3,1)));
        newArray(4,:) = new_weight;

        tsdf_color(gridIndex) = typecast(newArray(:),'uint32');
        toc;
    end

    % visualizing the voxel
    %{
    figure(1)
    min_weight = min(tsdf_weight(:));
    max_weight = max(tsdf_weight(:));
    for ii=1:voxel.size_grid(3)
        if min(min(min(tsdf_value(:,:,ii)))) ~= max(max(max(tsdf_value(:,:,ii))))
            imagesc(tsdf_weight(:,:,ii)',[min_weight max_weight]); axis equal; axis tight; xlabel('x'); ylabel('y');
            title(['frame ' num2str(ii) ' = depth ' num2str((ii-1)*voxel.unit+voxel.range(3,1)) ' meter']);
            pause(0.05);
        end
    end

    figure(2)
    for ii=1:voxel.size_grid(3)
        if min(min(min(tsdf_value(:,:,ii)))) ~= max(max(max(tsdf_value(:,:,ii))))
            imagesc(tsdf_value(:,:,ii)',[-1 1]); axis equal; axis tight; xlabel('x'); ylabel('y');
            title(['frame ' num2str(ii) ' = depth ' num2str((ii-1)*voxel.unit+voxel.range(3,1)) ' meter']);
            pause(0.05);
        end
    end
    %}
end

%% meshing
disp('isosurfacing...');
tic;
fv = isosurface(tsdf_value,0);
toc;
% visualizaiton
%{
figure(3);
p = patch(fv);
p.FaceColor = 'red';
p.EdgeColor = 'none';
daspect([1,1,1])
view(3); axis tight
camlight
lighting gouraud
%}

%% remove bad points that are artifact of TSDF

XYZ = fv.vertices';

faces = fv.faces';

XYZ = XYZ([2 1 3],:);

XYZint = round(XYZ);

selected = 1< XYZint(1,:) & XYZint(1,:) < voxel.size_grid(1) & 1< XYZint(2,:) & XYZint(2,:) < voxel.size_grid(2) & 1< XYZint(3,:) & XYZint(3,:) < voxel.size_grid(3);
XYZint = XYZint(:,selected);
XYZ = XYZ(:,selected);
faces = faces(:,selected(faces(1,:)) & selected(faces(2,:)) & selected(faces(3,:)));
indNew = cumsum(double(selected));
faces = indNew(faces);

gridIndex = sub2ind(voxel.size_grid',XYZint(1,:),XYZint(2,:),XYZint(3,:));
valIndex = tsdf_value(gridIndex);
selected = -0.15 < valIndex & valIndex<0.15;
XYZint = XYZint(:,selected);
XYZ = XYZ(:,selected);
faces = faces(:,selected(faces(1,:)) & selected(faces(2,:)) & selected(faces(3,:)));
indNew = cumsum(double(selected));
faces = indNew(faces);
XYZworld = (XYZ-1) *  voxel.unit + repmat(voxel.range(:,1),1,size(XYZ,2));

%% coloring
if is_color
    RGBcolor = zeros(3,size(XYZworld,2));
    disp('');
    tic
    for i=1:size(XYZworld,2)
        RGBcolor(:,i) = interpolateTrilinearyColor(XYZ(1,i),XYZ(2,i),XYZ(3,i),tsdf_color);
    end
    toc;
else
    RGBcolor = [];
end

% visualizaiton

figure(4)
clf;
plot3(XYZworld(1,:),XYZworld(2,:),XYZworld(3,:),'.r');
XYZcam = depth2XYZcamera(data.K, depthRead(data.depth{frameIdtarget}));
XYZcam = reshape(XYZcam,[],4)';
XYZcam = XYZcam(:,XYZcam(4,:)~=0);
XYZcam = XYZcam(1:3,:);
hold on
plot3(XYZcam(1,:),XYZcam(2,:),XYZcam(3,:),'.b');
axis equal;
axis tight;
grid on;
xlabel('x');
ylabel('y');
zlabel('z');

figure(5);
clf
fvNew.vertices = XYZworld';
fvNew.faces = faces';
%fvNew.faces = fvNew.faces(randsample(size(fvNew.faces,1),10000),:);
p = patch(fvNew,'FaceColor','none','EdgeColor','red');
daspect([1,1,1])
view(3); axis tight


text:
function [mri, hdr, cpersist] = read_ctf_mri4(filename);

% READ_CTF_MRI reads header and imnage data from CTF format MRI file
%
% [mri, hdr] = read_ctf_mri(filename)
%
% See also READ_CTF_MEG4, READ_CTF_RES4

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

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

if fid<=0
  error(sprintf('could not open MRI file: %s\n', filename));
end

[cpersist] = read_cpersist(fid);

% turn warnings off
ws = warning('off');

% general header information
hdr.identifierString        = get_value(cpersist, '_CTFMRI_VERSION');      % CTF_MRI_FORMAT VER 4.1
hdr.imageSize               = get_value(cpersist, '_CTFMRI_SIZE');         % 256
hdr.dataSize                = get_value(cpersist, '_CTFMRI_DATASIZE');     % 1 or 2(bytes)
hdr.orthogonalFlag          = get_value(cpersist, '_CTFMRI_ORTHOGONALFLAG');   % if set then image is orthogonal
hdr.interpolatedFlag        = get_value(cpersist, '_CTFMRI_INTERPOLATEDFLAG'); % if set than image was interpolated
hdr.comment                 = get_value(cpersist, '_CTFMRI_COMMENT');

hdr.Image.modality          = get_value(cpersist, '_SERIES_MODALITY');
hdr.Image.manufacturerName  = get_value(cpersist, '_EQUIP_MANUFACTURER');
hdr.Image.instituteName     = get_value(cpersist, '_EQUIP_INSTITUTION');
hdr.Image.imagedNucleus     = get_value(cpersist, '_MRIMAGE_IMAGEDNUCLEUS');
hdr.Image.FieldStrength     = get_value(cpersist, '_MRIMAGE_FIELDSTRENGTH');
hdr.Image.EchoTime          = get_value(cpersist, '_MRIMAGE_ECHOTIME');
hdr.Image.RepetitionTime    = get_value(cpersist, '_MRIMAGE_REPETITIONTIME');
hdr.Image.InversionTime     = get_value(cpersist, '_MRIMAGE_INVERSIONTIME');
hdr.Image.FlipAngle         = get_value(cpersist, '_MRIMAGE_FLIPANGLE');

% euler angles to align MR to head coordinate system(angles in degrees !)
rotation = split_nvalue(get_value(cpersist, '_CTFMRI_ROTATE'));
hdr.rotate_coronal  = rotation(1);
hdr.rotate_sagittal = rotation(2);
hdr.rotate_axial    = rotation(3);

hdr.transformMatrix = split_nvalue(get_value(cpersist, '_CTFMRI_TRANSFORMMATRIX'));
hdr.transformMatrix = reshape(hdr.transformMatrix, 4, 4)';

mmPerPixel = split_nvalue(get_value(cpersist, '_CTFMRI_MMPERPIXEL'));
hdr.mmPerPixel_sagittal = mmPerPixel(1);
hdr.mmPerPixel_coronal  = mmPerPixel(2);
hdr.mmPerPixel_axial    = mmPerPixel(3);

% HeadModel_Info specific header items
hmNasion = split_nvalue(get_value(cpersist, '_HDM_NASION'));
hdr.HeadModel.Nasion_Sag = hmNasion(1);
hdr.HeadModel.Nasion_Cor = hmNasion(2);
hdr.HeadModel.Nasion_Axi = hmNasion(3);

hmLeftEar = split_nvalue(get_value(cpersist, '_HDM_LEFTEAR'));
hdr.HeadModel.LeftEar_Sag = hmLeftEar(1);
hdr.HeadModel.LeftEar_Cor = hmLeftEar(2);
hdr.HeadModel.LeftEar_Axi = hmLeftEar(3);

hmRightEar = split_nvalue(get_value(cpersist, '_HDM_RIGHTEAR'));
hdr.HeadModel.RightEar_Sag = hmRightEar(1);
hdr.HeadModel.RightEar_Cor = hmRightEar(2);
hdr.HeadModel.RightEar_Axi = hmRightEar(3);

hmSphere = split_nvalue(get_value(cpersist, '_HDM_DEFAULTSPHERE'));
hdr.HeadModel.defaultSphereX = hmSphere(1);
hdr.HeadModel.defaultSphereY = hmSphere(2);
hdr.HeadModel.defaultSphereZ = hmSphere(3);
hdr.HeadModel.defaultSphereRadius = hmSphere(4);

hmOrigin = split_nvalue(get_value(cpersist, '_HDM_HEADORIGIN'));
hdr.headOrigin_sagittal = hmOrigin(1);
hdr.headOrigin_coronal = hmOrigin(2);
hdr.headOrigin_axial = hmOrigin(3);

%fread(fid,204,'char'); % unused, padding to 1028 bytes

% revert to previous warning state
warning(ws);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ THE IMAGE DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mri = zeros(256, 256, 256);

for slice = 1:256
  name = sprintf('_CTFMRI_SLICE_DATA#%.5d', slice);
  offset = get_value(cpersist, name);

  fseek(fid, offset, 'bof');

  if(hdr.dataSize == 1)
    slicedata = uint8(fread(fid, [256 256], 'uint8'));
  elseif(hdr.dataSize == 2)
    slicedata = uint16(fread(fid, [256 256], 'uint16'));
  else
    error('Unknown datasize in CTF MRI file');
  end;

  mri(:, :, slice) = slicedata;
end;

%mri = reshape(mri, [256 256 256]);
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO POST-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Data is stored in PIR order (i.e. fastest changing direction goes from
% anterior to Posterior, then from superior to Inferior, finally from left
% to Right) assuming subject position was not too oblique and proper
% conversion to .mri format.  (Does this depend on hdr.imageOrientation?)
% On the other hand, the transformation matrix and fiducials are in RPI
% order so we must reorient the image data to match them.
mri = permute(mri, [3 1 2]);

transformMatrix = hdr.transformMatrix;

% Construct minimal transformation matrix if fiducials were not defined.
% Bring data into same orientation (head coordinates are ALS) at least.
if all(transformMatrix == 0)
  transformMatrix(1, 2) = -1;
  transformMatrix(2, 1) = -1;
  transformMatrix(3, 3) = -1;
  transformMatrix(4, 4) = 1;
end

% determine location of fiducials in MRI voxel coordinates
% flip the fiducials in voxel coordinates to correspond to the previous flip along left-right
hdr.fiducial.mri.nas = [hdr.HeadModel.Nasion_Sag hdr.HeadModel.Nasion_Cor hdr.HeadModel.Nasion_Axi];
hdr.fiducial.mri.lpa = [hdr.HeadModel.LeftEar_Sag hdr.HeadModel.LeftEar_Cor hdr.HeadModel.LeftEar_Axi];
hdr.fiducial.mri.rpa = [hdr.HeadModel.RightEar_Sag hdr.HeadModel.RightEar_Cor hdr.HeadModel.RightEar_Axi];

% Reorient the image data, the transformation matrix and the fiducials
% along the left-right direction.
% This may have been done only for visualization?  It can probably be
% "turned off" without problem.
if false
  mri = flipdim(mri, 1);
  flip = [-1 0 0 hdr.imageSize+1
           0 1 0 0
           0 0 1 0
           0 0 0 1    ];
  transformMatrix = flip*transformMatrix;
  hdr.fiducial.mri.nas = [hdr.fiducial.mri.nas, 1] * flip(1:3, :)';
  hdr.fiducial.mri.lpa = [hdr.fiducial.mri.lpa, 1] * flip(1:3, :)';
  hdr.fiducial.mri.rpa = [hdr.fiducial.mri.rpa, 1] * flip(1:3, :)';
end

% re-compute the homogeneous transformation matrices (apply voxel scaling)
scale = eye(4);
scale(1,1) = hdr.mmPerPixel_sagittal;
scale(2,2) = hdr.mmPerPixel_coronal;
scale(3,3) = hdr.mmPerPixel_axial;
hdr.transformHead2MRI = transformMatrix*inv(scale);
hdr.transformMRI2Head = scale*inv(transformMatrix);

% compute location of fiducials in MRI and HEAD coordinates
hdr.fiducial.head.nas = warp_apply(hdr.transformMRI2Head, hdr.fiducial.mri.nas, 'homogenous');
hdr.fiducial.head.lpa = warp_apply(hdr.transformMRI2Head, hdr.fiducial.mri.lpa, 'homogenous');
hdr.fiducial.head.rpa = warp_apply(hdr.transformMRI2Head, hdr.fiducial.mri.rpa, 'homogenous');

%
% Reads a series of delimited numbers from a string
%
% @param input string Delimited string to process
% @param delim string The delimiter (default: \\)
%
% @return values matrix Array containing the numbers found
%
  function [values] = split_nvalue(input, delim)
    if(nargin < 2), delim = '\\'; end;

    remain = input;
    values = [];

    while(numel(remain > 0))
      [value, remain] = strtok(remain, delim);
      values(end + 1) = str2num(value);
    end;
  end

%
% Reads a value from the CPersist structure
%
% @param cpersist struct-array The CPersist structure
% @param key string The name of the parameter
%
% @return value mixed The value of the named parameter
%
  function [value] = get_value(cpersist, key)
    idx = find(strcmp({cpersist.key}, key));

    if(numel(idx) < 1), error('Specified key does not exist.'); end;
    if(numel(idx) > 1), error('Specified key is not unique.'); end;

    value = cpersist(idx).value;
  end

%
% Processes the CTF CPersist structure into a struct-array
%
% @param fid numeric File handle from which to read the CPersist structure
% @return cpersist struct-array
%
  function [cpersist] = read_cpersist(fid)
    magic = char(fread(fid, 4, 'char'))';

    if(~strcmp(magic, 'WS1_')), error('Invalid CPersist header'); end;

    cpersist = struct('key', {}, 'value', {});

    while(~feof(fid))
      % Read label
      lsize = fread(fid, 1, 'int32');
      ltext = char(fread(fid, lsize, 'char'))';

      % Last label in file is always EndOfParameters
      if(strcmp(ltext, 'EndOfParameters')), return; end;

      % Read value
      vtype = fread(fid, 1, 'int32');
      value = read_cpersist_value(fid, vtype);

      cpersist(end + 1).key = ltext;
      cpersist(end).value = value;
    end
  end

%
% Reads a single value of type (vtype) from fid
%
% @param fid numeric  The file to read the value from
% @param vtype numeric The type of value to read
%
% @return value mixed The read value
%
  function [value] = read_cpersist_value(fid, vtype)
    switch vtype
      case 3
        vsize = fread(fid, 1, 'int32');
        value = ftell(fid);
        fseek(fid, vsize, 'cof');
      case 4
        value = fread(fid, 1, 'double');
      case 5
        value = fread(fid, 1, 'int32');
      case 6
        value = fread(fid, 1, 'int16');
      case 10
        vsize = fread(fid, 1, 'int32');
        value = char(fread(fid, vsize, 'char'))';
      otherwise
        error(['Unsupported valuetype (' num2str(vtype) ') found in CPersist object']);
        return
    end
  end

end

author: PatternRecognition

repo: OpenBMI


text:
function r = rank(varargin)
% rank = 1

r = 1;

author: mtex-toolbox

repo: mtex

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


text:
% Defines the class currentSourceViewer for visualization of EEG inverse solutions on a realistic head model.
% This class is part of MoBILAB software.
% For more details visit:  https://code.google.com/p/mobilab/
%
% Author: Alejandro Ojeda, SCCN/INC/UCSD, Jan-2012

classdef currentSourceViewer < handle
    properties
        hFigure
        hAxes
        streamObj
        hLabels
        hSensors
        hScalp
        hCortex
        hVector
        dcmHandle
        sourceMagnitud
        sourceOrientation
        scalpData
        pointer
        clim
    end
    methods
        function obj = currentSourceViewer(streamObj,J,V,figureTitle,channelLabels)
            if nargin < 3, V = [];end
            if nargin < 4, figureTitle = '';end
            if nargin < 5
                channelLabels = cell(size(V,1),1);
                for it=1:length(channelLabels), channelLabels{it} = num2str(it);end
            end
            if isempty(channelLabels)
                channelLabels = cell(size(V,1),1);
                for it=1:length(channelLabels), channelLabels{it} = num2str(it);end
            end
            obj.streamObj = streamObj;
            load(obj.streamObj.surfaces);
            color = [0.93 0.96 1];

            path = fileparts(which('currentSourceViewer'));
            loc = strfind(path,'mobilab');
            if ~isempty(loc)
                path = path(1:loc+length('mobilab')-1);
            end
            path = fullfile(path,'skin');

            try labelsOn  = imread([path filesep 'labelsOn.png']);
                labelsOff = imread([path filesep 'labelsOff.png']);
                sensorsOn = imread([path filesep 'sensorsOn.png']);
                sensorsOff = imread([path filesep 'sensorsOff.png']);
                scalpOn = imread([path filesep 'scalpOn.png']);
                scalpOff = imread([path filesep 'scalpOff.png']);
                vectorOn = imread([path filesep 'vectorOn.png']);
                vectorOff = imread([path filesep 'vectorOff.png']);
                prev = imread([path filesep '32px-Gnome-media-seek-backward.svg.png']);
                next = imread([path filesep '32px-Gnome-media-seek-forward.svg.png']);
            catch ME
                ME.rethrow;
            end
            if isa(streamObj,'struct'), visible = 'off';else visible = 'on';end
            obj.hFigure = figure('Menubar','figure','ToolBar','figure','renderer','opengl','Visible',visible,'Color',color,'Name',figureTitle);
            position = get(obj.hFigure,'Position');
            set(obj.hFigure,'Position',[position(1:2) 1.06*position(3:4)]);
            obj.hAxes = axes('Parent',obj.hFigure);
            toolbarHandle = findall(obj.hFigure,'Type','uitoolbar');

            hcb(1) = uitoggletool(toolbarHandle,'CData',labelsOff,'Separator','on','HandleVisibility','off','TooltipString','Labels On/Off','userData',{labelsOn,labelsOff},'State','off');
            set(hcb(1),'OnCallback',@(src,event)rePaint(obj,hcb(1),'labelsOn'),'OffCallback',@(src, event)rePaint(obj,hcb(1),'labelsOff'));

            hcb(2) = uitoggletool(toolbarHandle,'CData',sensorsOff,'Separator','off','HandleVisibility','off','TooltipString','Sensors On/Off','userData',{sensorsOn,sensorsOff},'State','off');
            set(hcb(2),'OnCallback',@(src,event)rePaint(obj,hcb(2),'sensorsOn'),'OffCallback',@(src, event)rePaint(obj,hcb(2),'sensorsOff'));

            hcb(3) = uitoggletool(toolbarHandle,'CData',scalpOff,'Separator','off','HandleVisibility','off','TooltipString','Scalp On/Off','userData',{scalpOn,scalpOff},'State','off');
            set(hcb(3),'OnCallback',@(src,event)rePaint(obj,hcb(3),'scalpOn'),'OffCallback',@(src, event)rePaint(obj,hcb(3),'scalpOff'));

            hcb(4) = uitoggletool(toolbarHandle,'CData',vectorOff,'Separator','off','HandleVisibility','off','TooltipString','Vectors On/Off','userData',{vectorOn,vectorOff},'State','off');
            set(hcb(4),'OnCallback',@(src,event)rePaint(obj,hcb(4),'vectorOn'),'OffCallback',@(src, event)rePaint(obj,hcb(4),'vectorOff'));

            uipushtool(toolbarHandle,'CData',prev,'Separator','off','HandleVisibility','off','TooltipString','Previous','ClickedCallback',@obj.prev);
            uipushtool(toolbarHandle,'CData',next,'Separator','off','HandleVisibility','off','TooltipString','Next','ClickedCallback',@obj.next);
            set(obj.hFigure,'WindowScrollWheelFcn',@(src, event)mouseMove(obj,[], event));

            obj.dcmHandle = datacursormode(obj.hFigure);
            obj.dcmHandle.SnapToDataVertex = 'off';
            set(obj.dcmHandle,'UpdateFcn',@(src,event)showLabel(obj,event));
            obj.dcmHandle.Enable = 'off';
            hold(obj.hAxes,'on');

            obj.hSensors = scatter3(obj.hAxes,obj.streamObj.channelSpace(:,1),obj.streamObj.channelSpace(:,2),...
                obj.streamObj.channelSpace(:,3),'filled','MarkerEdgeColor','k','MarkerFaceColor','y');
            set(obj.hSensors,'Visible','off');

            N = length(channelLabels);
            k = 1.1;
            obj.hLabels = zeros(N,1);
            for it=1:N, obj.hLabels(it) = text('Position',k*obj.streamObj.channelSpace(it,:),'String',channelLabels{it},'Parent',obj.hAxes);end
            set(obj.hLabels,'Visible','off');

            if size(J,1) == 3*size(surfData(end).vertices,1)
                J = reshape(J,[size(J,1)/3 3 size(J,2)]);
                Jm = squeeze(sqrt(sum(J.^2,2)));
                normals = J;
            else
                Jm = J;
                normals = geometricTools.getSurfaceNormals(surfData(end).vertices,surfData(end).faces,false);
            end
            obj.sourceMagnitud = Jm;
            obj.sourceOrientation = J;
            obj.pointer = 1;

            % vectors
            obj.hVector = quiver3(surfData(end).vertices(:,1),surfData(end).vertices(:,2),surfData(end).vertices(:,3),normals(:,1,1),normals(:,2,1),normals(:,3,1),2);
            set(obj.hVector,'Color','k','Visible','off');

            % cortex
            obj.hCortex = patch('vertices',surfData(end).vertices,'faces',surfData(end).faces,'FaceVertexCData',obj.sourceMagnitud(:,1),...
                'FaceColor','interp','FaceLighting','phong','LineStyle','none','FaceAlpha',1,'SpecularColorReflectance',0,...
                'SpecularExponent',50,'SpecularStrength',0.5,'Parent',obj.hAxes);
            camlight(0,180)
            camlight(0,0)

            % scalp
            if isempty(V)
                skinColor = [1,.75,.65];
                obj.hScalp = patch('vertices',surfData(1).vertices,'faces',surfData(1).faces,'facecolor',skinColor,...
                    'facelighting','phong','LineStyle','none','FaceAlpha',.85,'Parent',obj.hAxes,'Visible','off');
                obj.scalpData = [];
            else
                W = geometricTools.localGaussianInterpolator(obj.streamObj.channelSpace,surfData(1).vertices,32);
                obj.scalpData = W*V;
                obj.hScalp = patch('vertices',surfData(1).vertices,'faces',surfData(1).faces,'FaceVertexCData',obj.scalpData(:,1),...
                    'FaceColor','interp','FaceLighting','phong','LineStyle','none','FaceAlpha',0.85,'SpecularColorReflectance',0,...
                    'SpecularExponent',50,'SpecularStrength',0.5,'Parent',obj.hAxes,'Visible','off');
            end
            view(obj.hAxes,[90 0]);

            disp('Calibrating the color scale...')
            mxsrc = obj.getRobustLimits(obj.sourceMagnitud(:),10);
            mxscp = obj.getRobustLimits(obj.scalpData,1);
            obj.clim = struct('source',[-mxsrc mxsrc],'scalp',[-mxscp mxscp]);
            disp('Done.')

            colorbar
            % box on;
            hold(obj.hAxes,'off');
            axis(obj.hAxes,'equal','vis3d');
            axis(obj.hAxes,'off')
            set(obj.hAxes,'Clim',obj.clim.source);
            try
                colormap(bipolar(512, 0.99))
            catch
                warning('Bipolar colormap is missing, fallback with jet.')
            end
            if isprop(streamObj,'name'), objName = [streamObj.name ': '];else objName = '';end
            set(obj.hFigure,'Visible',visible,'userData',obj,'Name',[objName num2str(1) '/' num2str(size(obj.sourceMagnitud,2))]);
            rotate3d
            drawnow
        end
        function [mx,mn] = getRobustLimits(obj,vect,th)
            if isempty(vect)
                mx = 1;
                mn = -1;
                return
            end
            samples = unidrnd(numel(vect),round(0.75*numel(vect)),20);
            prc = prctile(vect(samples),[th 100-th]);
            mn = median(prc(1,:));
            mx = median(prc(2,:));
            if mx == mn && mx == 0
                mx = max(vect(:));
                mn = min(vect(:));
            end
        end
       %%
        function rePaint(obj,hObject,opt)
            CData = get(hObject,'userData');
            if isempty(strfind(opt,'Off'))
                set(hObject,'CData',CData{2});
            else
                set(hObject,'CData',CData{1});
            end
            switch opt
                case 'labelsOn'
                    set(obj.hLabels,'Visible','on');
                case 'labelsOff'
                    set(obj.hLabels,'Visible','off');
                case 'sensorsOn'
                    set(obj.hSensors,'Visible','on');
                case 'sensorsOff'
                    set(obj.hSensors,'Visible','off');
                case 'scalpOn'
                    set(obj.hCortex,'FaceAlpha',0.15);
                    if strcmp(get(obj.hVector,'Visible'),'on')
                        set(obj.hScalp,'Visible','on','FaceAlpha',0.65);
                    else
                        set(obj.hScalp,'Visible','on','FaceAlpha',0.85);
                    end
                    set(get(obj.hScalp,'Parent'),'Clim',obj.clim.scalp);
                case 'scalpOff'
                    set(obj.hScalp,'Visible','off');
                    set(obj.hCortex,'Visible','on','FaceAlpha',1);
                    set(get(obj.hCortex,'Parent'),'Clim',obj.clim.source);
                case 'vectorOn'
                    set(obj.hVector,'Visible','on');
                    set(obj.hCortex,'FaceAlpha',0.75);
                case 'vectorOff'
                    set(obj.hVector,'Visible','off');
                    set(obj.hCortex,'FaceAlpha',1);
            end
        end
        %%
        function output_txt = showLabel(obj,event_obj)
            persistent DT
            if strcmp(obj.dcmHandle.Enable,'off'),return;end
            if strcmp(get(obj.hVector,'Visible'),'on')
                set(obj.hVector,'Visible','off');
                set(obj.hCortex,'FaceAlpha',1);
            end
            if isempty(DT)
                load(obj.streamObj.surfaces);
                vertices = surfData(end).vertices;
                DT = DelaunayTri(vertices(:,1),vertices(:,2),vertices(:,3));
            end
            pos = get(event_obj,'Position');
            loc = nearestNeighbor(DT, pos);
            output_txt = obj.streamObj.atlas.label{obj.streamObj.atlas.colorTable(loc)};
            drawnow
            %updateCursor(obj.dcmHandle,pos);
        end
        %%
        function prev(obj,~,~)
            obj.pointer = obj.pointer-1;
            if obj.pointer < 1, obj.pointer = 1;end
            val = obj.sourceMagnitud(:,obj.pointer);
            set(obj.hCortex,'FaceVertexCData',val);
            if isprop(obj.streamObj,'name'), objName = [obj.streamObj.name ': '];else objName = '';end
            set(obj.hFigure,'Name',[objName num2str(obj.pointer) '/' num2str(size(obj.sourceMagnitud,2))]);
            if isempty(obj.scalpData), drawnow;return;end
            val = obj.scalpData(:,obj.pointer);
            set(obj.hScalp,'FaceVertexCData',val);
            try %#ok
                set(obj.hVector,'UData',obj.sourceOrientation(:,1,obj.pointer),'VData',obj.sourceOrientation(:,2,obj.pointer),'WData',obj.sourceOrientation(:,3,obj.pointer));
            end
        end
        function next(obj,~,~)
            obj.pointer = obj.pointer+1;
            n = size(obj.sourceMagnitud,2);
            if obj.pointer > n, obj.pointer = n;end
            val = obj.sourceMagnitud(:,obj.pointer);
            set(obj.hCortex,'FaceVertexCData',val);
            if isprop(obj.streamObj,'name'), objName = [obj.streamObj.name ': '];else objName = '';end
            set(obj.hFigure,'Name',[objName num2str(obj.pointer) '/' num2str(size(obj.sourceMagnitud,2))]);
            if isempty(obj.scalpData), drawnow;return;end
            val = obj.scalpData(:,obj.pointer);
            set(obj.hScalp,'FaceVertexCData',val);
            try %#ok
                set(obj.hVector,'UData',obj.sourceOrientation(:,1,obj.pointer),'VData',obj.sourceOrientation(:,2,obj.pointer),'WData',obj.sourceOrientation(:,3,obj.pointer));
            end
            drawnow
        end
        function mouseMove(obj,~,eventObj)
            obj.pointer = obj.pointer - eventObj.VerticalScrollCount;%*eventObj.VerticalScrollAmount;
            if obj.pointer < 1, obj.pointer = 1;end
            if obj.pointer > size(obj.sourceMagnitud,2), obj.pointer = size(obj.sourceMagnitud,2);end
            val = obj.sourceMagnitud(:,obj.pointer);
            if isprop(obj.streamObj,'name'), objName = [obj.streamObj.name ': '];else objName = '';end
            set(obj.hFigure,'Name',[objName num2str(obj.pointer) '/' num2str(size(obj.sourceMagnitud,2))]);
            set(obj.hCortex,'FaceVertexCData',val);
            if isempty(obj.scalpData), drawnow;return;end
            val = obj.scalpData(:,obj.pointer);
            set(obj.hScalp,'FaceVertexCData',val);
        end
    end
end

repo: REST


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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%CLASSSIZES Get sizes of classes in a dataset
%
%   [N,LABLIST] = CLASSSIZES(A)
%
% INPUT
%   A        Dataset
%
% OUTPUT
%   N        Vector of the class cardinalities
%   LABLIST  Label list
%
% DESCRIPTION
% Returns class cardinalities of the dataset A. The order of the classes
% is identical to the list of class labels returned in LABLIST.

author: marianux

repo: ecg-kit

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


text:
%% New features in MatlabBGL version 4.0
% MatlabBGL 3.0 was only released as a beta test.  Version 4.0 is a
% full release with quite a few new features.

%% Improved reweighted graphs
% In MatlabBGL 3.0, reweighted graphs were a pain to use.   Now they are
% simple!  We combine a structural matrix with a weight matrix.  As(i,j)=1
% if there is an edge between vertex i and j and A(i,j)=wij where wij is
% the weight of the edge between i and j.

As = cycle_graph(6,'directed',0); % compute a 6 node cycle graph
A = As; % set all the weights to be one initially
A(2,3) = 0; A(3,2) = 0; % make one edge have zero weight
fprintf('Edges\n');
full(As)
fprintf('Weights\n');
full(A)

%%
% Note that As is given as the graph in the following call, not A!
[d pred] = shortest_paths(As,1,'edge_weight',edge_weight_vector(As,A));
d(3) % distance from vertex 1 to vertex 3 should be just 1!

%% Graph layout algorithms
% Sometimes, it's really nice to see a picture of your graph.  The BGL
% implements a few graph layout algorithms and so these are now in
% MatlabBGL 4.0!

%%
G = grid_graph(6,5);
X = kamada_kawai_spring_layout(G);
gplot(G,X,'.-');

%%
G = grid_graph(6,5);
X = fruchterman_reingold_force_directed_layout(G);
gplot(G,X,'.-');

%%
G = grid_graph(6,5);
X = gursoy_atun_layout(G);
gplot(G,X,'.-');

%% Planar graph algorithms
% The Boost Graph Library received a new suite of planar graph algorithms.
% These are now in MatlabBGL too.

%%
% A grid in the xy plane is a planar graph.
G = grid_graph(6,5);
is_planar = boyer_myrvold_planarity_test(G)

%%
% Recall that K_5 (the clique on 5 vertices) is not a planar graph.  Let's
% see what happens.
K5 = clique_graph(5);
is_planar = test_planar_graph(K5) % helpful wrapper

%%
% We can also draw planar graphs
G = grid_graph(6,5);
X = chrobak_payne_straight_line_drawing(G);
gplot(G,X,'.-'); % it looks a little different!

%% New option syntax
% You probably noticed that the "struct" command that permeated MatlabBGL
% calls before is gone in these examples.  We've moved to a new option
% syntax that gives you the _choice_ between the MatlabBGL struct style
% arguments and a list of key-value pairs

%%
% We'll look at spanning trees on the clique graph with 5 vertices.
% Using Prim's algorithm, the spanning tree we get depends on the root.  We
% always get a star graph rooted at the vertex we pick as the root.
G = clique_graph(5);

%%
% Old style
full(mst(G,struct('root',5,'algname','prim')))

%%
% New style

%%
% Just to make sure it works
full(mst(G,'root',1,'algname','prim'))

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%ARDEM	Demonstrates modules of the ARfit package.

%   Revised: 30-Dec-99 Tapio Schneider

format short
format compact

echo on
clc

%  ARfit is a collection of Matlab modules for the modeling of
%  multivariate time series with autoregressive (AR) models. ARfit
%  contains modules for estimating parameters of AR models from given
%  time series data; for checking the adequacy of an estimated
%  AR model; for analyzing eigenmodes of an estimated AR model; and
%  for simulating AR processes.
%
%  This demo illustrates the use of ARfit with a bivariate AR(2)
%  process
%
%       v(k,:)' = w' + A1*v(k-1,:)' + A2*v(k-2,:)' + eta(k,:)',
%
%  where the row vectors eta(k,:) are independent and identically
%  distributed Gaussian noise vectors with zero mean and covariance
%  matrix C. The kth row v(k,:) of the 2-column matrix v represents an
%  observation of the process at instant k. The intercept vector w is
%  included to allow for a nonzero mean of the AR(p) process.

pause   	% Press any key to continue.

clc

%  Let us simulate observations from a bivariate AR(2) process,
%  choosing the parameters

w = [ 0.25 ;  0.1 ];

% for the intercept vector,

A1 = [ 0.4   1.2;   0.3   0.7 ];

%  and

A2 = [ 0.35 -0.3;  -0.4  -0.5 ];

%  for the AR coefficient matrices, and

C = [ 1.00  0.50;   0.50  1.50 ];

%  for the noise covariance matrix. The two 2x2 matrices A1 and A2 are
%  assembled into a single 2x4 coefficient matrix:

A = [ A1 A2 ];

%  We use the module ARSIM to simulate 200 observations of this AR
%  process:

v = arsim(w, A, C, 200);

pause   	% Press any key to continue.

clc

%  Suppose that we have no information about how the time series
%  data v are generated, but we want to try to fit an AR model to the
%  time series. That is, we must estimate the AR model parameters,
%  including the model order. Assuming that the correct model order
%  lies between 1 and 5, we use the module ARFIT to determine the
%  optimum model order using Schwarz's Bayesian Criterion (SBC):

pmin = 1;
pmax = 5;
[west, Aest, Cest, SBC, FPE, th] = arfit(v, pmin, pmax);

%  The output arguments west, Aest, and Cest of ARFIT are the
%  estimates of the intercept vector w, of the coefficient matrix A,
%  and of the noise covariance matrix C. (The matrix th will be needed
%  later in the computation of confidence intervals.) These parameters
%  are estimated for a model of order popt, where the optimum model
%  order popt is chosen as the minimizer of an approximation to
%  Schwarz's Bayesian Criterion. The selected order popt in our
%  example is:

m    = 2;                 % dimension of the state space
popt = size(Aest, 2) / m;

disp(['popt = ', num2str(popt)])

pause   	% Press any key to continue.

clc

%  Besides the parameter estimates for the selected model, ARFIT has
%  also returned approximations to Schwarz's Bayesian Criterion SBC
%  and to Akaike's Final Prediction Error FPE, each for models of
%  order pmin:pmax. In this demo, the model order was chosen as the
%  minimizer of SBC. Here are the SBCs for the fitted models of order
%  1,...,5:

disp(SBC)

%  To see if using Akaike's Final Prediction Error as a criterion to
%  select the model order would have resulted in the same optimum
%  model order, compare the FPEs:

disp(FPE)

%  Employing FPE as the order selection criterion, the optimum model
%  order would have been chosen as the minimizer of FPE. The values of
%  the order selection criteria are approximations in that in
%  evaluating an order selection criterion for a model of order p <
%  pmax, the first pmax-p observations of the time series are ignored.

pause   	% Press any key to continue.

clc

%  Next it is necessary to check whether the fitted model is adequate
%  to represent the time series v. A necessary condition for model
%  adequacy is the uncorrelatedness of the residuals. The module ARRES

[siglev,res] = arres(west,Aest,v);

%  returns the time series of residuals res as well as the
%  significance level siglev with which a modified Li-McLeod
%  portmanteau test rejects the null hypothesis that the residuals are
%  uncorrelated. A model passes this test if, say, siglev > 0.05.  In
%  our example, the significance level of the modified Li-McLeod
%  portmanteau statistic is

disp(siglev);

%  (If siglev is persistently smaller than about 0.05, even over
%  several runs of this demo, then there is most likely a problem with
%  the random number generator that ARSIM used in the simulation of
%  v.)

pause   	% Press any key to continue.

clc

if exist('xcorr')    % If the Signal Processing Toolbox is installed, plot
                     % autocorrelation function of residuals ...

%  Using ACF, one can also plot the autocorrelation function of, say,
%  the first component of the time series of residuals:

acf(res(:,1));

%  95% of the autocorrelations for lag > 0 should lie between the
%  dashdotted confidence limits for the autocorrelations of an IID
%  process of the same length as the time series of residuals res.

%  Since uncorrelatedness of the residuals is only a necessary
%  condition for model adequacy, further diagnostic tests should, in
%  practice, be performed. However, we shall go on to the estimation
%  of confidence intervals.

pause   	% Press any key to continue.
end

clc

%  Being reasonably confident that the model adequately represents the
%  data, we compute confidence intervals for the AR parameters with
%  the module ARCONF:

[Aerr, werr] = arconf(Aest, Cest, west, th);

%  The output arguments Aerr and werr contain margins of error such
%  that (Aest +/- Aerr) and (west +/- werr) are approximate 95%
%  confidence intervals for the individual AR coefficients and for the
%  components of the intercept vector w. Here is the estimated
%  intercept vector with margins of error for the individual
%  components in the second column:

disp([west werr])

%  For comparison, the `true' vector as used in the simulation:

disp(w)

pause   	% Press any key to display the other parameter estimates.

clc

%  The estimated coefficient matrix:

disp(Aest)

%  with margins of error for its elements:

disp(Aerr)

%  For comparison, the `true' AR coefficients:

disp(A)

pause   	% Press any key to continue.

echo off
%  Compute `true' eigenmodes from model parameters:
%  Eigenvectors and eigenvalues of corresponding AR(1) coefficient
%  matrix:
[Strue,EvTrue] = eig([A; eye(2) zeros(2)]);
EvTrue         = diag(EvTrue)';   % `true' eigenvalues
Strue          = adjph(Strue);    % adjust phase of eigenmodes

clc
echo on

%  Finally, ARMODE computes the eigendecomposition of the fitted AR
%  model:

[S, Serr, per, tau, exctn] = armode(Aest, Cest, th);

%  The columns of S are the estimated eigenmodes:

disp(S)

%  with margins of error Serr:

disp(Serr)

%  The intervals (S +/- Serr) are approximate 95% confidence intervals
%  for the individual components of the eigenmodes.  Compare the
%  estimated eigenmodes above with the eigenmodes obtained from the
%  `true' model parameters:

disp(Strue(3:4,:))

%  (Note that the estimated modes can be a permutation of the `true'
%  modes. The sign of the modes is also ambiguous.)

pause   	% Press any key to continue.

echo off
pertrue = 2*pi./abs(angle(EvTrue)); % `true' periods

clc

echo on
%  Associated with the eigenmodes are the following oscillation periods:

disp(per)

%  The second row contains margins of error for the periods such that
%  (per(1,k) +/- per(2,k)) is an approximate 95% confidence interval
%  for the period of eigenmode S(:,k). [Note that for a purely
%  relaxatory eigenmode, the period is infinite (Inf).] Compare the
%  above estimated periods with the `true' periods:

disp(pertrue)

pause   	% Press any key to get the damping time scales.

echo off
tautrue = -1./log(abs(EvTrue)); % `true' damping time scales

clc

echo on

%  The damping times associated with each eigenmode are:

disp(tau)

%  with margins of error again in the second row. For comparison, the
%  `true' damping times:

disp(tautrue)

pause   	% Press any key to get the excitation of each eigenmode.
echo off
%  Compute `true' excitation of eigenmodes from the designed parameters:
p  = 2;              % true model order

invStr = inv(Strue); % inverse of matrix with eigenvectors as columns

% covariance matrix of corresponding decoupled AR(1) system
CovDcpld = invStr*[C zeros(2,(p-1)*2); zeros((p-1)*2, p*2)]*invStr';

% diagonal of that covariance matrix
DgCovDcpld = real(diag(CovDcpld))';

% excitation
TrueExctn = DgCovDcpld(1:2*p)./(1-abs(EvTrue).^2);

% normalize excitation
TrueExctn = TrueExctn./sum(TrueExctn);

clc

echo on
%  ARMODE has also returned the excitations, measures of the relative
%  dynamical importance of the eigenmodes:

disp(exctn)

%  Compare the estimated excitations with the `true' excitations
%  computed from the parameters used in the simulation:

disp(TrueExctn)

echo off
disp('End')

author: Sable

repo: mcbench-benchmarks


text:
close all;
jobNames = {'SM+zDD+AnnealLin', 'SM+zDD', 'SM+cDD+AnnealLin', 'SM+cDD'};
jobIDs = [2134873 2134874 2134875 2134876 ];

for jj = jobIDs
    findBestAndWorstMCMCRunByLogPr(jj, 25);
end

plotLogPrVsTime( jobIDs(1:2), 1:25, jobNames(1:2) );
ylim( [-5.55 -5.34]*1e4 ); xlim( [0 35000]);
set(gcf,'Name', 'LogPr_zDD' );

plotLogPrVsTime( jobIDs(3:4), 1:25, jobNames(3:4) );
ylim( [-5.55 -5.34]*1e4 ); xlim( [0 35000]);
set(gcf,'Name', 'LogPr_cDD' );

plotLogPrVsTime( jobIDs(1:4), 1:25, jobNames(1:4) );
ylim( [-5.55 -5.34]*1e4 ); xlim( [0 35000]);
set(gcf,'Name', 'LogPr_CompareAnnealingInference' );

for jj = 1:length( jobNames )
    plotStateSeq(  jobIDs(jj), 'best' );
    set(gcf, 'Name', sprintf('StateSeq_best_%s', jobNames{jj})  );
end

plotHammDistVsTime( jobIDs(1:4), 1:25, jobNames(1:4) );

plotHammDistVsTime( jobIDs([1:2 4]), 1:25, jobNames([1:2 4]) );

author: michaelchughes

repo: NPBayesHMM


text:
% test for local fourier transform

n = 128;
name = 'barb';
M = crop(load_image(name),n);

w = 15; q = 7;
options.window_type = 'sin';
MF = perform_windowed_fourier_transform(M,w,q,n, options);
M1 = perform_windowed_fourier_transform(MF,w,q,n, options);

e = sum(M(:).^2);
ef = sum(abs(MF(:)).^2);
disp( ['Energy conservation error (should be 0) : ' num2str( (e-ef)/e*100 ) '%' ]);

author: gpeyre

repo: matlab-toolboxes


text:
function [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,options)
% spike inference using a constrained deconvolution approach:
%      min      sum(sp)
%    c,sp,b,c1
%      subject to: sp >= 0
%                   b >= 0
%                  G*c = sp
%                   c1 >= 0
%           ||y-b-c - c_in|| <= sn*sqrt(T)

%   Variables:
%   y:      raw fluorescence data (vector of length(T))
%   c:      denoised calcium concentration (Tx1 vector)
%   b:      baseline concentration (scalar)
%  c1:      initial concentration (scalar)
%   g:      discrete time constant(s) (scalar or 2x1 vector)
%  sn:      noise standard deviation (scalar)
%  sp:      spike vector (Tx1 vector)

%   USAGE:
%   [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,OPTIONS)
%   The parameters b,cin,g,sn can be given or else are estimated from the data

%   OPTIONS: (stuct for specifying options)
%         p: order for AR model, used when g is not given (default 2)
%    method: methods for performing spike inference
%   available methods: 'dual' uses dual ascent
%                       'cvx' uses the cvx package available from cvxr.com (default)
%                      'lars' uses the least regression algorithm
%                     'spgl1' uses the spgl1 package available from
%                     math.ucdavis.edu/~mpf/spgl1/  (usually fastest)
%   bas_nonneg:   flag for setting the baseline lower bound. if 1, then b >= 0 else b >= min(y)
%   noise_range:  frequency range over which the noise power is estimated. Default [Fs/4,Fs/2]
%   noise_method: method to average the PSD in order to obtain a robust noise level estimate
%   lags:         number of extra autocovariance lags to be considered when estimating the time constants
%   resparse:     number of times that the solution is resparsened (default 0). Currently available only with methods 'cvx', 'spgl'
%   fudge_factor: scaling constant to reduce bias in the time constant estimation (default 1 - no scaling)

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

defoptions.p = 2;
defoptions.method = 'cvx';
defoptions.bas_nonneg = 1;              % nonnegativity option for baseline estimation
defoptions.noise_range = [0.25,0.5];    % frequency range over which to estimate the noise
defoptions.noise_method = 'logmexp';    % method for which to estimate the noise level
defoptions.lags = 5;                    % number of extra lags when computing the AR coefficients
defoptions.resparse = 0;                % number of times to re-sparse solution
defoptions.fudge_factor = 1;            % fudge factor for time constants

if nargin < 6
    options = defoptions;
    if nargin < 5
        sn = [];
        if nargin < 4
            g = [];
            if nargin < 3
                c1 = [];
                if nargin < 2
                    b = [];
                end
            end
        end
    end
end

if ~isfield(options,'p');  options.p = defoptions.p;  end
if ~isfield(options,'method'); options.method = defoptions.method; end
if ~isfield(options,'bas_nonneg'); options.bas_nonneg = defoptions.bas_nonneg; end
if ~isfield(options,'noise_range'); options.noise_range = defoptions.noise_range; end
if ~isfield(options,'noise_method'); options.noise_method = defoptions.noise_method; end
if ~isfield(options,'lags'); options.lags = defoptions.lags; end
if ~isfield(options,'resparse'); options.resparse = defoptions.resparse; end
if ~isfield(options,'fudge_factor'); options.fudge_factor = defoptions.fudge_factor; end

method = options.method;
if isempty(b);
    bas_est = 1;
else
    bas_est = 0;
end
if isempty(c1)
    c1_est = 1;
else
    c1_est = 0;
end

y = y(:);
T = length(y);
y_full = y;
mis_data = isnan(y);
E = speye(T);
E(mis_data,:) = [];

if any(mis_data)
    y_full(mis_data) = interp1(find(~mis_data),y(~mis_data),find(mis_data));
end

if isempty(sn)
    sn = GetSn(y_full,options.noise_range,options.noise_method);
end
if isempty(g)
    g = estimate_time_constants(y_full,options.p,sn,options.lags);
    while max(abs(roots([1,-g(:)']))>1) && options.p < 5
        warning('No stable AR(%i) model found. Checking for AR(%i) model \n',options.p,options.p+1);
        options.p = options.p + 1;
        g = estimate_time_constants(y,options.p,sn,options.lags);
    end
    if options.p == 5
        g = 0;
    end
    %fprintf('Stable AR(%i) model found \n',options.p);
    % re-adjust time constant values
    rg = roots([1;-g(:)]);
    if ~isreal(rg); rg = real(rg) + .001*randn(size(rg)); end
    rg(rg>1) = 0.95 + 0.001*randn(size(rg(rg>1)));
    rg(rg<0) = 0.15 + 0.001*randn(size(rg(rg<0)));
    pg = poly(options.fudge_factor*rg);
    g = -pg(2:end);
end
if options.bas_nonneg  % lower bound for baseline
    b_lb = 0;
else
    b_lb = min(y);
end

if strcmpi(method,'dual'); method = 'dual';
elseif strcmpi(method,'cvx'); method = 'cvx';
elseif strcmpi(method,'lars'); method = 'lars';
elseif strcmpi(method,'spgl1'); method = 'spgl1';
else fprintf('Invalid choice of method. Using CVX \n'); method = 'cvx';
end

if strcmpi(method,'dual') && any(mis_data)
    warning('Dual method does not support missing data. Switching to CVX');
    method = 'cvx';
end

if options.resparse > 0 && (strcmpi(method,'dual') || strcmpi(method,'lars'))
    warning('Resparsening is not supported with chosen method. Switching to CVX');
    method = 'cvx';
end

pathCell = regexp(path, pathsep, 'split');
g = g(:);
G = spdiags(ones(T,1)*[-g(end:-1:1)',1],-length(g):0,T,T);
gd = max(roots([1,-g']));  % decay time constant for initial concentration
gd_vec = gd.^((0:T-1)');

switch method
    case 'dual'
         v = G'*ones(T,1);
        thr = sn*sqrt(T-sum(mis_data));
        if bas_est; b = 0; end
        if c1_est; c1 = 0; end
        myfun = @(Ald) lagrangian_temporal_gradient(Ald,thr^2,y(~mis_data)-b-c1*gd_vec(~mis_data),bas_est,c1_est);
        c = [G\max(G*y,0);zeros(bas_est);zeros(c1_est)];
        options_dual = optimset('GradObj','On','Display','Off','Algorithm','interior-point','TolX',1e-8);
        ld_in = 10;
        [ld,~,flag] = fmincon(myfun,ld_in,[],[],[],[],0,[],[],options_dual);
        if (flag == -2) || (flag == -3)
            warning('Problem seems unbounded or infeasible. Try a different method.');
        end
        if bas_est; b = c(T+bas_est); end
        if c1_est; c1 = c(end); end
        c = c(1:T);
        sp = G*c;
    case 'cvx'
        onPath = ~isempty(which('cvx_begin'));
        if onPath
            c = zeros(T,1+options.resparse);
            sp = zeros(T,1+options.resparse);
            bas = zeros(1+options.resparse,1);
            cin = zeros(1+options.resparse,1);
            w_ = ones(T,1);
            for rep = 1:options.resparse+1
                [c(:,rep),bas(rep),cin(rep)] = cvx_foopsi(y,b,c1,sn,b_lb,g,w_,~mis_data);
                sp(:,rep) = G*c(:,rep);
                w_ = 1./(max(sp(:,rep),0) + 1e-8);
            end
            sp(sp<1e-6) = 0;
            c = G\sp;
            b = bas;
            c1 = cin;
        else
            error('CVX does not appear to be on the MATLAB path. It can be downloaded from cvxr.com \n');
        end
    case 'lars'
         Ginv = E*[full(G\speye(T)),ones(T,bas_est),gd_vec*ones(1,c1_est)];
         if bas_est; b = 0; end
         if c1_est; c1 = 0; end
         [~, ~, spikes, ~, ~] = lars_regression_noise(y(~mis_data)-b_lb*bas_est - b - c1*gd_vec(~mis_data), Ginv, 1, sn^2*(T-sum(mis_data)));
         sp = spikes(1:T);
         b = (spikes(T+bas_est)+b_lb)*bas_est + b*(1-bas_est);
         c1 = spikes(end)*c1_est + c1*(1-c1_est);
         c = G\sp;
    case 'spgl1'
        onPath = ~isempty(which('spgl1'));
        if onPath
            Gx = @(x,mode) G_inv_mat(x,mode,T,g,gd_vec,bas_est,c1_est,E);
            c = zeros(T,1+options.resparse);
            sp = zeros(T,1+options.resparse);
            bas = zeros(1+options.resparse,1);
            cin = zeros(1+options.resparse,1);
            w_ = ones(T,1);
            for rep = 1:options.resparse+1
                if bas_est; b = 0; w_ = [w_;1e-10]; end
                if c1_est; c1 = 0; w_ = [w_;1e-10]; end
                options_spgl = spgSetParms('project',@NormL1NN_project ,'primal_norm', @NormL1NN_primal,'dual_norm',@NormL1NN_dual,'verbosity',0,'weights',w_);
                [spikes,r,~,~] = spg_bpdn( Gx, y(~mis_data)-b_lb*bas_est - (1-bas_est)*b-(1-c1_est)*c1*gd_vec(~mis_data), sn*sqrt(T-sum(mis_data)),options_spgl);
                c(:,rep) = G\spikes(1:T); %Gx([spikes(1:T);0],1);                                  %% calcium signal
                bas(rep) = b*(1-bas_est) + bas_est*spikes(T+bas_est)+b_lb*bas_est;       %% baseline
                cin(rep) = c1*(1-c1_est) + c1_est*spikes(end);
                sp(:,rep) = spikes(1:T);                                           %% spiking signal
                w_ = 1./(spikes(1:T)+1e-8);
            end
            b = bas;
            c1 = cin;
            %sn = norm(r)/sqrt(T);
        else
            error('SPGL1 does not appear to be on the MATLAB path. It can be downloaded from math.ucdavis.edu/~mpf/spgl1 \n');
        end
end

    function sn = GetSn(Y,range_ff,method)
        % estimate noise level with a power spectral density method
        L=length(Y);
%         if ~isempty(which('pmtm'))
%             [psd_Y,ff] = pmtm(Y,5/2,1000,1);
%         end
        if ~isempty(which('pwelch'));
            [psd_Y,ff]=pwelch(Y,round(L/8),[],1000,1);
        else
            xdft = fft(Y);
            xdft = xdft(:,1:round(L/2)+1);
            psd_Y = (1/L) * abs(xdft).^2;
            ff = 0:1/L:1/2;
            psd_Y(2:end-1) = 2*psd_Y(2:end-1);
        end
        ind=ff>range_ff(1);
        ind(ff>range_ff(2))=0;
        switch method
            case 'mean'
                sn=sqrt(mean(psd_Y(ind)/2));
            case 'median'
                sn=sqrt(median(psd_Y(ind)/2));
            case 'logmexp'
                sn = sqrt(exp(mean(log(psd_Y(ind)/2))));
        end
    end

    function g = estimate_time_constants(y,p,sn,lags)
        % estimate time constants from the sample autocovariance function

        lags = lags + p;
        if ~isempty(which('xcov')) %signal processing toolbox
            xc = xcov(y,lags,'biased');
        else
            ynormed = (y - mean(y));
            xc = nan(lags + 1, 1);
            for k = 0:lags
                xc(k + 1) = ynormed(1 + k:end)' * ynormed(1:end - k);
            end
            xc = [flipud(xc(2:end)); xc] / numel(y);
        end
        xc = xc(:);
        A = toeplitz(xc(lags+(1:lags)),xc(lags+(1:p))) - sn^2*eye(lags,p);
        g = pinv(A)*xc(lags+2:end);
    end

    function [f,grad] = lagrangian_temporal_gradient(Al,thr,y_raw,bas_flag,c1_flag)
        options_qp = optimset('Display','Off','Algorithm','interior-point-convex');
        H = [speye(T),ones(T,bas_flag),gd_vec*ones(1,c1_flag);...
            ones(bas_flag,T),T*ones(bas_flag),(1-gd^T)/(1-gd)*ones(c1_flag,bas_flag);...
            (gd_vec*ones(1,c1_flag))',(1-gd^T)/(1-gd)*ones(c1_flag,bas_flag),(1-gd^(2*T))/(1-gd^2)*ones(c1_flag,c1_flag)];
        Ay = [y_raw;sum(y_raw)*ones(bas_flag);gd_vec'*y_raw*ones(c1_flag)];
        c = quadprog(2*Al(1)*H,[v;zeros(bas_flag+c1_flag,1)]-2*Al(1)*Ay,[-G,sparse(T,bas_flag+c1_flag);sparse(bas_flag+c1_flag,T),-speye(bas_flag+c1_flag)]...
            ,[sparse(T,1);-b_lb*ones(bas_flag);zeros(c1_flag)],[],[],[],[],c,options_qp);
        f = v'*c(1:T);
        grad = [sum((c(1:T)-y_raw + c(T+bas_flag)*bas_flag + c(end)*gd_vec*c1_flag).^2)-thr];
        f = f + Al(:)'*grad;
    end

    function b = G_inv_mat(x,mode,NT,gs,gd_vec,bas_flag,c1_flag,Emat)
        if mode == 1
            b = filter(1,[1;-gs(:)],x(1:NT)) + bas_flag*x(NT+bas_flag) + c1_flag*gd_vec*x(end);
            b = Emat*b;
           %b = G\x(1:NT) + x(NT+bas_flag)*bas_flag + x(end)*c1_flag;
        elseif mode == 2
            x = Emat'*x;
            b = [flipud(filter(1,[1;-gs(:)],flipud(x)));ones(bas_flag,1)*sum(x);ones(c1_flag,1)*(gd_vec'*x)];
           %b = [G'\x;ones(bas_flag,1)*sum(x);ones(c1_flag,1)*(gd_vec'*x)] ;
        end
    end

end

author: zhoupc


text:
function s = tt_intvec2str(v)
%TT_VEC2STR Print integer vector to a string with brackets.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isempty(v)
    s = sprintf('[]');
    return;
end

s = ['[ ' sprintf('%d ',v(1:end)) ']'];

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function com = calcCoM()
global uLINK

M  = TotalMass(1);
MC = calcMC(1);
com = MC / M;

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function [ElementRnCnPn]=FindPRC_3D_QPOTTS(ElementNo,state)

e=ElementNo;
planes=size(state,3);
columns=size(state,2);
rows=size(state,1);
r=rows;
c=columns;

for k=1:planes
    if e<=(k*r*c) & e>((k-1)*r*c)
        pn=k;break
    end
end

for j=1:columns
    if e<=((pn-1)*r*c+j*r) & e>((pn-1)*r*c+(j-1)*r)
        cn=j;break
    end
end

for i=1:rows
    if e==(pn-1)*r*c+(cn-1)*r+i
        rn=i;break
    end
end

ElementRnCnPn=[rn,cn,pn];

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34985-monte-carlo-simulation-of-three-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 3D square-lattice - microstructure/FindPRC_3D_QPOTTS.m


text:
classdef PrincipalDirectionComputerIn3D < PrincipalDirectionComputer

    methods (Access = public, Static)

        function obj = PrincipalDirectionComputerIn3D(cParams)
            obj.ndim = 3;
            obj.init(cParams)
        end

    end

    methods (Access = public)

        function compute(obj,tensor)
            s1  = squeeze(tensor(1,1,:));
            s2  = squeeze(tensor(1,2,:));
            s3  = squeeze(tensor(1,3,:));
            s12 = squeeze(tensor(1,4,:));
            s13 = squeeze(tensor(1,5,:));
            s23 = squeeze(tensor(1,6,:));
            eG = obj.eigenComputer;
            for i = 1:obj.ndim
                for j = 1:obj.ndim
                    d = eG.eigenVectorFunction{i,j}(s1,s12,s13,s2,s23,s3);
                    obj.direction(i,j,:) = real(d);
                    obj.assertSmallImaginaryValue(d);
                end
            end
        end

    end

    methods (Access = private, Static)

        function assertSmallImaginaryValue(d)
            imagD = imag(d);
%             assert(norm(imagD(:)) < 1e6) % To be fixed
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function display_result_progress(cv, results_to_display, display_progress)

 % Displays the current mean and stdev when a single resample result is left out for different result types
 %  to give an idea whether the algorithm has converged to a stable decoding accuracy
 %  (these values could be used for stopping criteria, rather than running a fixed number of resampling runs)

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.edu)
%
%     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/>.

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

i = 1;
num_training_times = 0;
while num_training_times(1) == 0
    if i ~= 4
        num_training_times = size(results_to_display{i}, 3);
        num_test_times = size(results_to_display{i},4);
    else
        num_training_times = size(results_to_display{i}, 4);
        num_test_times = size(results_to_display{i}, 5);
    end
    i = i + 1;
end

 % default behavior if there is only 1 training point and display_progress.training_time_to_display_results == -1, is to just use this one training time for the results
 if (display_progress.training_time_to_display_results == -1)  && (num_training_times == 1)
     display_progress.training_time_to_display_results = 1;
 end

 % display_progress.training_time_to_display_results should be set if the num_training_times ~= num_test_times
 if (display_progress.training_time_to_display_results == -1) && (num_training_times~= num_test_times) && (num_test_times ~= 1)

     warning(['If the number of training times does not equal the number of test times, and one wants to display the resampled CV progress, ' ...
                'then display_progress.training_time_to_display_results should be set.  Using the default value of display_progress.training_time_to_display_results = 1, ' ...
                'to display the results here.'])

     display_progress.training_time_to_display_results = 1;

 end

current_results = [];

 for iResult = 1:length(results_to_display)

     iResample = size(results_to_display{1}, 1);

     if iResample == 0
         continue;
     end

    % if a particular training time that has been specified to display the results
    % note: if only saving results when training and testing at the same time period, then the displayed results will only be for one time point
    if display_progress.training_time_to_display_results > 0

        if iResult ~= 4
           curr_bootstrap_progress_results = squeeze(mean(results_to_display{iResult}(:, :, display_progress.training_time_to_display_results, :), 2));
        else
           curr_bootstrap_progress_results = squeeze(mean(mean(results_to_display{iResult}(:, :, :, display_progress.training_time_to_display_results, :), 3), 2));
        end

        if iResample == 1
           curr_bootstrap_progress_results = curr_bootstrap_progress_results';
        end

    else

        % if only saving the results for training and testing at the same time
        if ((length(size(results_to_display{iResult})) == 3) && (iResult ~=4)) || ((length(size(results_to_display{iResult})) == 4)  && (iResult == 4))

            if iResult ~=4
                 curr_bootstrap_progress_results = squeeze(mean(results_to_display{iResult}, 2));
            else
                 curr_bootstrap_progress_results = squeeze(mean(mean(results_to_display{iResult}, 3), 2));
            end

            if iResample == 1
               curr_bootstrap_progress_results = curr_bootstrap_progress_results';
            end

        else

             for iDisplayResultTrainingTime = 1:num_training_times

                 if iResult ~= 4
                    curr_bootstrap_progress_results(:, iDisplayResultTrainingTime) = squeeze(mean(results_to_display{iResult}(:, :, iDisplayResultTrainingTime, iDisplayResultTrainingTime), 2));
                 else
                    curr_bootstrap_progress_results(:, iDisplayResultTrainingTime)  = squeeze(mean(mean(results_to_display{iResult}(:, :, :, iDisplayResultTrainingTime, iDisplayResultTrainingTime), 3), 2));
                 end
             end

        end

    end

   if iResample == 1
       current_results = [current_results; curr_bootstrap_progress_results; NaN .* ones(size(curr_bootstrap_progress_results))];   % can't show cofficients of variation if only 1 bootstrap has been run (so skip the following part)
   else

       curr_total_mean = mean(curr_bootstrap_progress_results, 1);
       curr_one_bootstrap_left_out_mean = ((repmat(curr_total_mean, [iResample 1]) - (curr_bootstrap_progress_results./iResample)) .* (iResample./(iResample -1)))';   % means over bootstraps with the ith bootstrap result left out
       curr_stdev_one_bootstrap_left_out = std(curr_one_bootstrap_left_out_mean');  % stdev over results when one bootstrap sample is left out

       current_results = [current_results; curr_total_mean; curr_stdev_one_bootstrap_left_out; NaN .* ones(size(curr_total_mean))];

   end

 end

fprintf('\n')
disp(current_results)
fprintf('\n')

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


%  Script file: lsqfit.m
%
%  Purpose:
%    To perform a least-squares fit of an input data set
%    to a straight line, and print out the resulting slope
%    and intercept values.  The input data for this fit
%    comes from a user-specified input data file.
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/24/07    S. J. Chapman        Original code
%
% Define variables:
%   count     -- number of values read
%   filename  -- Input file name
%   fid       -- File id
%   msg       -- Open error message
%   n         -- Number of input data pairs (x,y)
%   slope     -- Slope of the line
%   sum_x     -- Sum of all input X values
%   sum_x2    -- Sum of all input X values squared
%   sum_xy    -- Sum of all input X*Y values
%   sum_y     -- Sum of all input Y values
%   x         -- An input X value
%   x_bar     -- Average X value
%   y         -- An input Y value
%   y_bar     -- Average Y value
%   y_int     -- Y-axis intercept of the line

% Initialize sums
n = 0; sum_x = 0; sum_y = 0; sum_x2 = 0; sum_xy = 0;

% Prompt user and get the name of the input file.
disp('This program performs a least-squares fit of an');
disp('input data set to a straight line. Enter the name');
disp('of the file containing the input (x,y) pairs:  ' );
filename = input(' ','s');

% Open the input file
[fid,msg] = fopen(filename,'rt');

% Check to see if the open failed.
if fid < 0

   % There was an error--tell user.
   disp(msg);

else

   % File opened successfully. Read the (x,y) pairs from
   % the input file.  Get first (x,y) pair before the
   % loop starts.
   [in,count] = fscanf(fid,'%g %g',2);

   while ~feof(fid)
      x = in(1);
      y = in(2);
      n      = n + 1;                  %
      sum_x  = sum_x + x;              % Calculate
      sum_y  = sum_y + y;              %   statistics
      sum_x2 = sum_x2 + x.^2;          %
      sum_xy = sum_xy + x * y;         %

      % Get next (x,y) pair
      [in,count] = fscanf(fid,'%f',[1 2]);

   end

   % Close the file
   fclose(fid);

   % Now calculate the slope and intercept.
   x_bar = sum_x / n;
   y_bar = sum_y / n;
   slope = (sum_xy - sum_x*y_bar) / (sum_x2 - sum_x*x_bar);
   y_int = y_bar - slope * x_bar;

   % Tell user.
   fprintf('Regression coefficients for the least-squares line:\n');
   fprintf('   Slope (m)     = %12.3f\n',slope);
   fprintf('   Intercept (b) = %12.3f\n',y_int);
   fprintf('   No of points  = %12d\n',n);

end


text:
filename='test2d_micro';
%filename = 'MicroQuad';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'chomog_alphabeta'};
weights=[1];
constraint = {'volumeConstraint'};
constraint_case = 'EQUALITY';
%incrementFactor = 1;
%designVariable = 'Density';
designVariable = 'LevelSet';
filterType = 'P1';
fracRadius = 0.51;
%optimizer = 'IPOPT';
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';

%optimizerUnconstrained = 'PROJECTED GRADIENT';
%line_search_initiator = 'INCREASING LAST STEP';
%incrementFactor = 1.95;

optimizerUnconstrained = 'SLERP';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target=1;
optimality_final = 0.5*1e-3;
constr_final =1e-12;

Vfrac_initial = 0.8;
optimality_initial = 0.5*1e-3;
constr_initial = 1e-12;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 1 0]';
micro.beta = [1 1 0]';

%nsteps = 10;

% For all tests
plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 200;


text:
function area = triangle_area ( vert1, vert2, vert3 )

%*****************************************************************************80
%
%% TRIANGLE_AREA returns the area of a triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    26 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real VERT1(2), VERT2(2), VERT3(2), the vertices of
%    the triangle.
%
%    Output, real AREA, the area of the triangle.
%
  area = 0.5 * ...
    ( ...
        ( vert2(1) - vert1(1) ) * ( vert3(2) - vert1(2) ) ...
      - ( vert3(1) - vert1(1) ) * ( vert2(2) - vert1(2) ) ...
    );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Js = JacobianSpace_sym(Slist, thetalist)
% *** CHAPTER 5: VELOCITY KINEMATICS AND STATICS ***
% Takes Slist: The joint screw axes in the space frame when the manipulator
%              is at the home position, in the format of a matrix with the
%              screw axes as the columns,
%       thetalist: A list of joint coordinates.
% Returns the corresponding space Jacobian (6xn real numbers).
% Example Input:
%
% clear; clc;
% Slist = [[0; 0; 1;   0; 0.2; 0.2], ...
%        [1; 0; 0;   2;   0;   3], ...
%        [0; 1; 0;   0;   2;   1], ...
%        [1; 0; 0; 0.2; 0.3; 0.4]];
% thetalist = [0.2; 1.1; 0.1; 1.2];
% Js = JacobianSpace(Slist, thetalist)
%
% Output:
% Js =
%         0    0.9801   -0.0901    0.9575
%         0    0.1987    0.4446    0.2849
%    1.0000         0    0.8912   -0.0453
%         0    1.9522   -2.2164   -0.5116
%    0.2000    0.4365   -2.4371    2.7754
%    0.2000    2.9603    3.2357    2.2251

Js = sym(Slist);
T = eye(4);
for i = 2: length(thetalist)
    T = T * simplify(MatrixExp6_sym(VecTose3(Slist(:, i - 1)),thetalist(i - 1)));
	Js(:, i) = Adjoint(T) * Slist(:, i);
end
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function spm_MDP_urn
% Demo for active inference with the urn task
%__________________________________________________________________________
%
% This demonstration uses the Urn or Beads Task to illustrate how choice
% behaviour can be simulated using active inference - in the context of
% Markov decision processes. In the urn task, a succession of draws
% from one of two urns are made and the agent has to decide whether the
% balls are being drawn from an urn with predominantly red or green balls.
% We model this in terms of a state-space with four dimensions: number of
% balls drawn (n),  number of green balls drawn (k), choice (undecided, red
% or green)and the true (hidden) state of the urn (red or green). With
% this relatively simple state-space, the utility of any hidden state is
% simply quantified by the log-odds ratio of making a correct
% decision. From binomial theory this is (2k - n)*log(p/(1 - p)), where p
% is the proportion of red or green balls. Having defined the utility
% function of states, we can then use the MDP formulation of active
% inference using variational Bayes to simulate choice behaviour.
%
% This routine first provides an illustration of a game in which a decision
% is delayed until the last draw to look at inferences during successive
% draws - with a special focus on precision. The illustration here shows
% a decrease in precision when an unexpected (green ball) is drawn during a
% sequence of red balls.
%
% We then characterise changes in choice probability (and latency to the
% decision) in terms of its dependency on threshold criteria (on the odds
% ratio) and hyperpriors about precision (alpha or the scale parameter of a
% standard gamma distribution). The routine concludes with an illustration
% of how to estimate model parameters using the likelihood of observed
% (simulated) choices.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_urn.m 6592 2015-11-06 16:20:48Z guillaume $

% set up and preliminaries
%==========================================================================
rng('default')

T     = 8;                         % number of trials in each game
Pa    = .85;                       % proportion of red or green goals

% hidden (and initial) states (a red urn with no drawers)
%--------------------------------------------------------------------------
D          = zeros(T,T,3,2);       % #balls - 1 x #green - 1 x #u x #urns
D(1,1,1,1) = 1;
S          = spm_vec(D);
D(1,1,1,2) = 1;
D          = spm_vec(D);

% likelihood - everything is seen  apart from the hidden states of the urn
%--------------------------------------------------------------------------
A     = kron([1 1],eye(3*T*T));

% transition probabilities (B{1} - wait; B{2} - red; B{3} - green)
%--------------------------------------------------------------------------
for i = 1:6
    Bn{i,i} = kron(spm_speye(T,T),spm_speye(T,T));
end
Br      = spm_speye(T,T,0)*Pa + spm_speye(T,T,-1)*(1 - Pa);
Bg      = spm_speye(T,T,0)*(1 - Pa) + spm_speye(T,T,-1)*Pa;
Bn{1,1} = kron(Br,spm_speye(T,T,-1));
Bn{4,4} = kron(Bg,spm_speye(T,T,-1));

B{1}    = spm_cat(Bn);
B{2}    = kron(eye(2),kron([0 0 0;1 1 0;0 0 1],eye(T*T)));
B{3}    = kron(eye(2),kron([0 0 0;0 1 0;1 0 1],eye(T*T)));

% priors over final state (exp(utility))
%--------------------------------------------------------------------------
L     = zeros(T,T,3,2);
W     = zeros(T,T,3,2);
rho   = log(Pa/(1 - Pa));
for n = 0:(T - 1)
    W(n + 1,:,:,:) = - n/8;
    for k = 0:n
        L(n + 1,k + 1,2,1) = (n - 2*k)*rho;
        L(n + 1,k + 1,3,2) = (2*k - n)*rho;
    end
end
C     = spm_softmax(spm_vec((L > 3) + W)*4);

% allowable policies (one decision before the game ends)
%--------------------------------------------------------------------------
V     = [ones(T,T) + eye(T,T), ones(T,T) + 2*eye(T,T)];

% MDP Structure
%==========================================================================
MDP.T = T;                          % process depth (the horizon)
MDP.S = S;                          % initial state
MDP.A = A;                          % likelihood
MDP.B = B;                          % transition probabilities (priors)
MDP.C = C;                          % terminal cost probabilities (priors)
MDP.D = D;                          % prior over initial states
MDP.V = V;                          % allowable policies

% Solve - an example game (with high offer at t = 10)
%==========================================================================
spm_figure('GetWin','Figure 1'); clf

% create a sequence of draws (outcomes) - with and oddball on the 4th
%--------------------------------------------------------------------------
k     = cumsum([0 0 0 1 0 0 0 0]);
a     = [1 1 1 1 1 1 2];
o     = zeros(1,length(a) - 1);
for n = 1:length(o)
    o(n) = sub2ind([T,T,3],n,k(n) + 1,1);
end

% Active inference (precluding a decision until the last trial)
%--------------------------------------------------------------------------
MDP.o    = o;
MDP.a    = a;
MDP.plot = gcf;
MDP.N    = 8;
MDP      = spm_MDP_game(MDP);

% plot convergence and precision
%--------------------------------------------------------------------------
subplot(4,2,7)
plot(MDP.d)
xlabel('Latency (updates)','FontSize',12)
ylabel('Precision of beliefs','FontSize',12)
title('Expected precision','FontSize',16)
spm_axis tight

% deconvolve to simulate dopamine responses
%--------------------------------------------------------------------------
nd  = length(MDP.d);
K   = tril(toeplitz(exp(-((1:nd) - 1)'/8)));

subplot(4,2,8)
plot(pinv(K)*MDP.d(:)), hold on
xlabel('Latency (updates)','FontSize',12)
ylabel('Precision of beliefs','FontSize',12)
title('Simulated dopamine responses','FontSize',16)
spm_axis tight

% Illustrate dependency on parameters
%==========================================================================
spm_figure('GetWin','Figure 2'); clf

% create a sequence of draws (outcomes) - no oddballs
%--------------------------------------------------------------------------
k     = cumsum([0 0 0 0 0 0 0 0]);
o     = zeros(1,length(k));
for n = 1:length(k)
    o(n) = sub2ind([T,T,3],n,k(n) + 1,1);
end

% probability distribution over time: P(1,:) is no action
%--------------------------------------------------------------------------
PrT      = @(P) [1 cumprod(P(1,1:end - 1))].*(1 - P(1,:));
MDP.plot = 0;                        % plot convergence
MDP.N    = 4;                        % number of variational iterations
MDP.a    = ones(1,T);                % and action
MDP.o    = o;                        % and outcomes

% beliefs about final state - decision threshold (log likelihood)
%--------------------------------------------------------------------------
DP    = MDP;
PF    = [];
DF    = [];
p     = linspace(0,8,8);
for i = 1:length(p)
    DP.C    = spm_softmax(spm_vec((L > p(i)) + W)*4);
    DP      = spm_MDP_game(DP);
    PF(i,:) = 1 - DP.P(1,:);
    DF(i,:) = PrT(DP.P);
end

% probability of accepting
%--------------------------------------------------------------------------
subplot(2,2,1)
plot(PF')
xlabel('Latency (trials)','FontSize',12)
ylabel('Probability of deciding','FontSize',12)
title('Increasing decision threshold','FontSize',16)
axis square xy

% distribution of acceptance latencies
%--------------------------------------------------------------------------
subplot(2,2,2)
plot(DF')
xlabel('Latency (trials)','FontSize',12)
ylabel('Latency of decision','FontSize',12)
title('Latency of decision','FontSize',16)
axis square xy

% Hyperpriors - prior precision (alpha)
%--------------------------------------------------------------------------
DP    = MDP;
PF    = [];
DF    = [];
p     = linspace(2,16,8);
for i = 1:length(p)
    DP.alpha  = p(i);
    DP      = spm_MDP_game(DP);
    PF(i,:) = 1 - DP.P(1,:);
    DF(i,:) = PrT(DP.P);
end

% probability of accepting
%--------------------------------------------------------------------------
subplot(2,2,3)
plot(PF')
xlabel('Latency (trials)','FontSize',12)
ylabel('Probability of deciding','FontSize',12)
title('Increasing prior precision','FontSize',16)
axis square xy

% distribution of acceptance latencies
%--------------------------------------------------------------------------
subplot(2,2,4)
plot(DF')
xlabel('Latency (trials)','FontSize',12)
ylabel('Latency of decision','FontSize',12)
title('Latency of decision','FontSize',16)
axis square xy

% Simulate multiple trials and try to infer prior precision (alpha)
%==========================================================================
spm_figure('GetWin','Figure 3'); clf

% Simulate multiple trials and record likelihood
%--------------------------------------------------------------------------
MDP.plot = 0;
MDP.N    = 4;

p     = linspace(2,16,8);
for t = 1:8

    % run game
    %----------------------------------------------------------------------
    MDP.s = [];
    MDP.o = [];
    MDP.a = [];
    MDP   = spm_MDP_game(MDP);

    % place outcomes in DP
    %----------------------------------------------------------------------
    DP    = MDP;
    [a j] = find(MDP.U);
    [o j] = find(MDP.O);
    DP.o  = o;
    DP.a  = a;
    y(t)  = find(a > 1);

    % get log-likelihood for different parameter values
    %----------------------------------------------------------------------
    for i = 1:length(p);
        DP.alpha = p(i);
        DP       = spm_MDP_game(DP);
        LL(i,t)  = sum(log(DP.P(find(DP.U))));
    end

end

% approximate the MAP with the ML and use the Laplace assumption
%--------------------------------------------------------------------------
Lp    = sum(LL,2);
dp    = mean(diff(p,1));
dLdpp = diff(Lp,2)/(dp^2);
[l i] = max(Lp(2:end - 1) + (dLdpp < 0)*1024);
Cp    = inv(-dLdpp(i));
Ep    = p(i + 1);

% plot responses
%--------------------------------------------------------------------------
subplot(2,2,1)
hist(y,1:T)
xlabel('Latency','FontSize',12)
ylabel('Probability','FontSize',12)
title('Distribution of responses','FontSize',16)
axis square

% plot log likelihood over trials
%--------------------------------------------------------------------------
subplot(2,2,2)
plot(p,LL - min(LL(:)))
xlabel('Parameter','FontSize',12)
ylabel('Probability','FontSize',12)
title('Log-likelihood over games','FontSize',16)
axis square

% plot likelihood
%--------------------------------------------------------------------------
subplot(2,2,3)
plot(p,Lp - min(Lp))
xlabel('Parameter','FontSize',12)
ylabel('Log-likelihood','FontSize',12)
title('Log-likelihood','FontSize',16)
axis square

% plot posterior
%--------------------------------------------------------------------------
subplot(2,2,4)
tp  = 8;
pp  = linspace(p(1),p(end),64);
plot(pp,spm_Npdf(pp,Ep,Cp)), hold on
plot([tp tp],[0 1/8],':'),   hold off
xlabel('Parameter','FontSize',12)
ylabel('Probability','FontSize',12)
title('Posterior probability','FontSize',16)
axis square

return

author: spm


text:
function[]=xlin
% XLIN   Sets x-axis scale to linear.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 J.M. Lilly --- type 'help jlab_license' for details

set(gca,'xscale','linear')

%Update for compatibility with ZOOL
if isappdata(gcf,'xlimitslock')
   if getappdata(gcf,'xlimitslock')
      h=gca;
      hall=axeshandles;
      for i=1:length(hall)
         axes(hall(i))
         set(gca,'xscale','linear')
      end
      axes(h)
   end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function a=randNestedParenth(n)
%%RANDNESTEDPARENTH Obtain a random character string of randomly nested
%                   parentheses. Such parantheses can represent random
%                   binary trees.
%
%INPUTS: n The number of open-closed parenthesis pairs.
%
%OUTPUTS: a A (2*n)X1 character string of n parenthesis pairs, with the
%           open parenthesis always coming before the close one.
%
%This function implements Algorithm W of Chapter 7.2.1.6 of [1].
%
%EXAMPLE:
% a=randNestedParenth(5).'
%We transpose it to make it easier to read.
%
%REFERENCES
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

p=n;
q=n;
m=1;
%Allocate space for the random set of nested parentheses.
a=repmat(' ',2*n,1);

while(q~=0)
    while(1)
        upperBound=(q+p)*(q-p+1);
        X=randi(upperBound-1);

        if(X<(q+1)*(q-p))
            q=q-1;
            a(m)=')';
            m=m+1;
            break;
        end
        p=p-1;
        a(m)='(';
        m=m+1;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function name = removestopwords(name, stopwords)
% remove stop of words. The returned sequence is also parsed into words:
%
% removestopwords('carFrontal')
% ans =
%   'car'
%
% You can also pass a cell array of strings with a list of words to be
% removed.
%
% If called without input arguments, it returs a cell array with the list
% of stopwords given by default.
%
% stopwords = removestopwords;
%
% You can change the list of stopwords by modifiying the file
% stopwords.txt.  (one stopword per line).
%
% Also, this function can be used to remove the stopwords from the LabelMe
% index. This might remove also useful information. For instance, the words
% 'blue', 'frontal', ... will disapear.
%
% D = removestopwords(D)
%

if nargin < 2
    % load list of stop words
    fid = fopen('stopwords.txt');
    C = textscan(fid,'%s');
    fclose(fid);
    stopwords = C{1};
end

if nargin > 0
    if isstruct(name)
        for i = 1:length(name)
            if isfield(name(i).annotation, 'object')
                N = length(name(i).annotation.object);
                for j = 1:N
                    nameobj = name(i).annotation.object(j).name;

                    % we avoid recursive call for speed.
                    w = lower(getwords(nameobj));
                    k = ismember(w, stopwords);
                    w = w(~k);
                    nameobj = strtrim(sprintf('%s ', w{:}));

                    name(i).annotation.object(j).name = nameobj;
                end
            end
        end
    else
        w = lower(getwords(name));
        j = ismember(w, stopwords);
        w = w(~j);
        name = strtrim(sprintf('%s ', w{:}));
    end
else
    name = stopwords;
end

author: CSAILVision

repo: LabelMeToolbox


text:
function [Q] = CouetteBC2D(xin, yin, nxin, nyin, mapI, mapO, mapW, mapC, Q, time);

% function [Q] = CouetteBC2D(xin, yin, nxin, nyin, mapI, mapO, mapW, mapC, Q, time);
% Purpose: evaluate solution for Couette flow

% Couette flow (mach .2 at inner cylinder)
gamma = 1.4;

% extract conserved variables
rho = Q(:,:,1); rhou = Q(:,:,2); rhov = Q(:,:,3); Ener = Q(:,:,4);

mapB = [mapI;mapO];

rad = sqrt(xin(mapB).^2 + yin(mapB).^2);
theta = atan2(yin(mapB), xin(mapB));

utheta = (-rad + 16./rad)/75;
p = 1 + (1./(75^2))*( (rad.^2)/2 - 32*log(rad) - 128./rad.^2 );

rho (mapB) = 1;
rhou(mapB) = -sin(theta).*utheta;
rhov(mapB) =  cos(theta).*utheta;
Ener(mapB) = p/(gamma-1) + 0.5*(rhou(mapB).^2 + rhov(mapB).^2)./rho(mapB);

% pack modified conserved variables
Q(:,:,1) = rho; Q(:,:,2) = rhou; Q(:,:,3) = rhov; Q(:,:,4) = Ener;
return

repo: Aero-matlab


text:
classdef GemImageSliceParameters < CoreBaseClass
    % GemImageSliceParameters. Parameters for extracting a 2D slice from a
    % 3D volume
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (SetObservable)
        Orientation = GemImageOrientation.XZ
        SliceNumber = [1, 1, 1] % The currently shown slice in 3 dimensions
        UpdateLock
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function stop = tfocs_stop( x, nonsmoothF, optTol )
% tfocs_stop : TFOCS stopping condition
%
%   $Revision: 0.1.2 $  $Date: 2012/09/15 $
%
  global subprob_Dg_y subprob_optim

  [ ~, x_prox ]   = nonsmoothF( x - subprob_Dg_y ,1);
    subprob_optim = norm( x_prox - x ,'inf');
    stop          = subprob_optim <= optTol;

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function dUpsilon = lfmvvGradientUpsilonMatrix(gamma, sigma2, t1, ...
    t2, mode, upsilon)

% LFMVVGRADIENTUPSILONMATRIX Gradient upsilon matrix vel. vel.
% FORMAT
% DESC computes the gradient of a portion of the LFMVV kernel.
% ARG gamma : Gamma value for system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG upsilon : precomputation of the upsilon matrix.
% ARG mode : operation mode, according to the derivative (mode 0,
% derivative wrt t1, mode 1 derivative wrt t2)
% RETURN upsilon : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010
%
% SEEALSO : lfmvvComputeUpsilonMatrix.m

% KERN

sigma = sqrt(sigma2);

if nargin<6
    upsilon = lfmvpComputeUpsilonMatrix(gamma, sigma2, t1, t2, mode);
    if nargin <5
        mode =0;
    end
end

if mode ==0
    dUpsilon = upsilon + gamma*lfmvpGradientUpsilonMatrix(gamma, sigma2, t1, t2, mode) ...
        - (2/(sqrt(pi)*sigma))*((1-gamma*t1).*exp(-gamma*t1))*(exp(-(t2.^2)/sigma2)).';
else
    dUpsilon = -upsilon - gamma*lfmvpGradientUpsilonMatrix(gamma, sigma2, t1, t2, mode);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function calendar_nyt_test04 ( )

%*****************************************************************************80
%
%% CALENDAR_NYT_TEST04 compares JED_TO_NYT and JED_TO_NYT_ISSUE_IDEAL.
%
%  Discussion:
%
%    JED_TO_NYT returns the issue number printed on the New York Times.
%    There were many "mistakes" and "accidents" and "corrections" in this system.
%
%    JED_TO_NYT_ISSUE_IDEAL returns an "ideal" issue number which keeps track
%    of every issue, in a sensible, usable way.
%
%  Volume  Issue   D M         Y
%  ------  -----  -- --------- ----
%       1      1  18 September 1851
%       1     92   3 January   1852
%       2    404   3 January   1853
%       2    416  17 January   1853
%       3    856  15 June      1854
%
%       4   1210   4 August    1855
%       5   1259   1 October   1855
%       5   1491  28 June      1856
%       6   1706   9 March     1857
%       7   2155  17 August    1858
%
%       8   2421  23 June      1859
%       9   2586   4 January   1860
%      10   2897   3 January   1861
%      10   3000   1 May       1861
%      12   3432  23 September 1862
%
%      13   3794  20 November  1863
%      14   4130  19 December  1864
%      14   4230  16 April     1865
%      15   4576  26 May       1866
%      17   5034  13 November  1867
%
%      17   5250  22 July      1868
%      20   6189  22 July      1871
%      47  14499   5 February  1898
%      47  15000   7 February  1898
%      52  16579  24 February  1903
%
%      53  16909  15 March     1904
%      54  17292  29 May       1905
%      55  17561  22 February  1906
%      57  18164  18 October   1907
%      58  18856   9 September 1909
%
%      66  21619   3 April     1917
%      74  24651  22 July      1925
%      89  29827  23 September 1939
%      89  30000  14 March     1940
%      94  31881   8 May       1945
%
%      94  31980  15 August    1945
%     113  38864  20 June      1964
%     114  39317  16 September 1965
%     117  40076  15 October   1967
%     118  40721  21 July      1969
%
%     120  41418  18 June      1971
%     127  44027   9 August    1978
%     128  44028   6 November  1978
%     141  48939  17 April     1992
%     144  50000  14 March     1995
%
%     147  50939   8 October   1997
%     149  51753  31 December  1999
%     149  51254   1 January   2000
%     150  51599  11 December  2000
%     154  53108  28 January   2005
%
%     157  54136  22 November  2007
%     157  54267  01 April     2008
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   04 February 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 52;
  issue_test = [ ...
        1,    92,   404,   416,   856, ...
     1210,  1259,  1491,  1706,  2155, ...
     2421,  2586,  2897,  3000,  3432, ...
     3794,  4130,  4230,  4576,  5034, ...
     5250,  6189, 14499, 15000, 16579, ...
    16909, 17292, 17561, 18164, 18856, ...
    21619, 24651, 29827, 30000, 31881, ...
    31980, 38864, 39317, 40076, 40721, ...
    41418, 44027, 44028, 48939, 50000, ...
    50939, 51753, 51254, 51599, 53108, ...
    54136, 54267 ];
  jed_test = [ ...
   2397383.50, 2397490.50, 2397856.50, 2397870.50, 2398384.50, ...
   2398799.50, 2398857.50, 2399128.50, 2399382.50, 2399908.50, ...
   2400218.50, 2400413.50, 2400778.50, 2400896.50, 2401406.50, ...
   2401829.50, 2402224.50, 2402342.50, 2402747.50, 2403283.50, ...
   2403535.50, 2404630.50, 2414325.50, 2414327.50, 2416169.50, ...
   2416554.50, 2416994.50, 2417263.50, 2417866.50, 2418558.50, ...
   2421321.50, 2424353.50, 2429529.50, 2429702.50, 2431583.50, ...
   2431682.50, 2438566.50, 2439019.50, 2439778.50, 2440423.50, ...
   2441120.50, 2443729.50, 2443818.50, 2448729.50, 2449790.50, ...
   2450729.50, 2451543.50, 2451544.50, 2451889.50, 2453398.50, ...
   2454426.50, 2454557.50 ];
  volume_test = [ ...
       1,    1,    2,    2,    3, ...
       4,    5,    5,    6,    7, ...
       8,    9,   10,   10,   12, ...
      13,   14,   14,   15,   17, ...
      17,   20,   47,   47,   52, ...
      53,   54,   55,   57,   58, ...
      66,   74,   89,   89,   94, ...
      94,  113,  114,  117,  118, ...
     120,  127,  128,  141,  144, ...
     147,  149,  149,  150,  154, ...
     157,  157 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALENDAR_NYT_TEST04\n' );
  fprintf ( 1, '  For the New York Times issue date:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED -> NYT1 by historical record.\n' );
  fprintf ( 1, '  JED -> NYT2 by "JED_TO_NYT"\n' );
  fprintf ( 1, '  JED -> NYT3 by "JED_TO_NYT_ISSUE_IDEAL"\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       JED                 Date             Volume     Issue     Issue     Issue\n' );
  fprintf ( 1, '                                                     (lookup)  (compute)  (Ideal)\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    jed = jed_test(test);

    [ y, m, d, f ] = jed_to_ymdf_common ( jed );
    s = ymdf_to_s_common ( y, m, d, f );

    issue1 = issue_test(test);
    volume1 = volume_test(test);

    [ volume2, issue2 ] = jed_to_nyt ( jed );

    issue3 = jed_to_nyt_issue_ideal ( jed );

    fprintf ( 1, '  %11.2f  %25s  %8d  %8d  %8d  %8d\n', ...
      jed, s, volume1, issue1, issue2, issue3 );

  end

  return
end