text:
function [X meaneffective sigmaeffective] = TruncatedGaussian(sigma, range, varargin)
% function X = TruncatedGaussian(sigma, range)
%          X = TruncatedGaussian(sigma, range, n)
%
% Purpose: generate a pseudo-random vector X of size n, X are drawn from
% the truncated Gaussian distribution in a RANGE braket; and satisfies
% std(X)=sigma.
% RANGE is of the form [left,right] defining the braket where X belongs.
% For a scalar input RANGE, the braket is [-RANGE,RANGE].
%
% X = TruncatedGaussian(..., 'double') or
% X = TruncatedGaussian(..., 'single') return an array X of of the
%    specified class.
%
% If input SIGMA is negative, X will be forced to have the same "shape" of
% distribution function than the unbounded Gaussian with standard deviation
% -SIGMA: N(0,-SIGMA). It is similar to calling RANDN and throw away values
% ouside RANGE. In this case, the standard deviation of the truncated
% Gaussian will be different than -SIGMA. The *effective* mean and
% the effective standard deviation can be obtained by calling:
%   [X meaneffective sigmaeffective] = TruncatedGaussian(...)
%
% Example:
%
% sigma=2;
% range=[-3 5]
%
% [X meaneff sigmaeff] = TruncatedGaussian(sigma, range, [1 1e6]);
%
% stdX=std(X);
% fprintf('mean(X)=%g, estimated=%g\n',meaneff, mean(X))
% fprintf('sigma=%g, effective=%g, estimated=%g\n', sigma, sigmaeff, stdX)
% hist(X,64)
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Last update: 19/April/2009

% We keep track this variables so as to avoid calling fzero if
% TruncatedGaussian is called succesively with the same sigma and range
persistent PREVSIGMA PREVRANGE PREVSIGMAC

% shape preserving?
shapeflag = (sigma<0);

% Force inputs to be double class
range = double(range);
if isscalar(range)
    % make sure it's positive
    range=abs(range);
    range=[-range range];
else
    range=sort(range); % right order
end
sigma = abs(double(sigma));

n=varargin;

if shapeflag
    % Prevent the same pdf as with the normal distribution N(0,sigma)
    sigmac = sigma;
else
    if diff(range)^2<12*sigma^2 % This imposes a limit of sigma wrt range
        warning('TruncatedGaussian:RangeSigmaIncompatible', ...
                'TruncatedGaussian: range and sigma are incompatible\n');
        sigmac = Inf;
    elseif isequal([sigma range], [PREVSIGMA PREVRANGE]) % See line #80
        sigmac = PREVSIGMAC; % Do not need to call fzero
    else
        % Search for "sigmac" such that the truncated Gaussian having
        % sigmac in the formula of its pdf gives a standard deviation
        % equal to sigma
        [sigmac res flag] = fzero(@scz,sigma,[],...
                                   sigma^2,range(1),range(2));
        sigmac = abs(sigmac); % Force it to be positive
        if flag<0 % Someting is wrong
            error('TruncatedGaussian:fzerofailled', ...
                  'Could not estimate sigmac\n');
        end
        % Backup the solution
        [PREVSIGMA PREVRANGE PREVSIGMAC] = deal(sigma,range,sigmac);
    end
end

% Compute effective standard deviation
meaneffective=meantrunc(range(1), range(2), sigmac);
sigmaeffective=stdtrunc(range(1), range(2), sigmac);

% Inverse of the cdf functions
if isinf(sigmac)
    % Uniform distribution to maximize the standard deviation within the
    % range. It is like a Gaussian with infinity standard deviation
    if any(strcmpi(n,'single'))
        range = single(range);
    end
    cdfinv = @(y) range(1)+y*diff(range);
else
    c = sqrt(2)*sigmac;
    e = erf(range/c);
    if any(strcmpi(n,'single'))
        % cdfinv will be single class
        c = single(c);
        e = single(e);
    end
    cdfinv = @(y) c*erfinv(e(1)+diff(e)*y);
end

% Generate random variable
X = cdfinv(rand(n{:}));
% Clip to prevent some nasty numerical issues with of erfinv function
% when argument gets close to +/-1
X = max(min(X,range(2)),range(1));

return

end % TruncatedGaussian

function m=meantrunc(lower, upper, s)
% Compute the mean of a trunctated gaussian distribution
if isinf(s)
    m = (upper+lower)/2;
else
    a = (lower/sqrt(2))./s;
    b = (upper/sqrt(2))./s;
    corr = sqrt(2/pi)*(-exp(-b.^2)+exp(-a.^2))./(erf(b)-erf(a));
    m = s.*corr;
end
end % vartrunc

function v=vartrunc(lower, upper, s)
% Compute the variance of a trunctated gaussian distribution
if isinf(s)
    v = (upper-lower)^2/12;
else
    a = (lower/sqrt(2))./s;
    b = (upper/sqrt(2))./s;
    if isinf(a)
        ea=0;
    else
        ea = a.*exp(-a.^2);
    end
    if isinf(b)
        eb = 0;
    else
        eb = b.*exp(-b.^2);
    end
    corr = 1 - (2/sqrt(pi))*(eb-ea)./(erf(b)-erf(a));
    v = s.^2.*corr;
end
end % vartrunc

function stdt=stdtrunc(lower, upper, s)
% Standard deviation of a trunctated gaussian distribution
stdt = sqrt(vartrunc(lower, upper, s)-meantrunc(lower, upper, s).^2);
end % stdtrunc

function res=scz(sc, targetsigma2, lower, upper)
% Gateway for fzero, aim the standard deviation to a target value
res = vartrunc(lower, upper, sc) - targetsigma2 - ...
      meantrunc(lower, upper, sc).^2;
end % scz

% End of file TruncatedGaussian.m


text:
%% im2patch
% Below is a demonstration of the features of the |im2patch| function

%%
clear; close all; clc;

%% Syntax
% |[F,V,C]=im2patch(M,IND,ptype);|

%% Description
% This function generates patch data for 3D images. The patch data is only
% generated for the voxels specified by the indices (logic or linear
% indices). The patch data is created according to the patch type desired
% (e.g. voxel type of slice type).

%% Examples

%%
% Plot settings
cMap=parula(250);
faceAlpha1=1;
faceAlpha2=0.5;
edgeColor1='none';
edgeColor2='none';
fontSize=15;

%% Example: Introduction to using |im2patch| for voxel plotting
% The voxel for which patch data is to be specified can be defined by
% supplying linear indices or a logic array.

% Simulating an image
M=rand(3,5,7);
% Example supplying linear indices, here all voxels
indPatch=1:numel(M);

%%
% Creating patch data for voxel display
[F,V,C]=im2patch(M,indPatch,'v');

cFigure;
title('patch type: v');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C,'k',faceAlpha2);
colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%%
% Study the size of the face and vertex arrays to confirm that the patch
% type |'v'| creates vertices and faces for each voxel. Shared vertices and
% faces are not removed. This option requires the most memory.
disp(num2str(size(F)));
disp(num2str(size(V)));

%%
% Creating patch data for voxel display with shared vertices and faces
% removed
[F,V,C]=im2patch(M,indPatch,'vu');

cFigure;
title('patch type: vu');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C,'k',faceAlpha2);
colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%%
% Study the size of the face and vertex arrays to confirm that the patch
% type |'vu'| ensures that shared faces and vertices are not shared. Each
% vertex and face is therefore unique. This saves memory with respect to
% using the |'v'| patch type for voxel display. The figures look identical
% except that shared faces may appear less dark when transparency is on
% since now one one face is used. Color information is shared too
% (averaged).
disp(num2str(size(F)));
disp(num2str(size(V)));

%%
% Creating patch data for voxel display with only non-shared faces and
% vertices
[F,V,C]=im2patch(M,indPatch,'vb');

cFigure;
title('patch type: vb');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C,'k',faceAlpha2);
colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%%
% Study the size of the face and vertex arrays to confirm that the patch
% type |'vb'| helps to plot only non-shared vertices and faces.
% In the case of an enclosed shape filled with voxels thus only the
% boundary faces are displayed. This path type appears no different than
% the other is transparency is not on however it is much lighter on memory
% than the above path types.
disp(num2str(size(F)));
disp(num2str(size(V)));

%% Example: Introduction to using |im2patch| for slice plotting

cFigure;
title('patch type: si, sj, sk');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;

%Setting up indices for I direction slices
S=round(size(M,1)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(S,:,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'si'); %Creating patch data for y mid-voxel slices
gpatch(F,V,C);

%Setting up indices for J direction slices
S=round(size(M,2)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,S,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sj'); %Creating patch data for x mid-voxel slices
gpatch(F,V,C);

%Setting up indices for Z direction slices
S=round(size(M,3)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,:,S)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sk'); %Creating patch data for z mid-voxel slices
gpatch(F,V,C);

colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;

drawnow;

%%
% The path type s*u are simular to s* but use shared vertices

%% Example: Comparison to standard MATLAB |imagesc| function and the patch type MATLAB functions |slice| and |pcolor|
% The comparison is what motivates the choice of coordinate system for
% im2patch i.e. that it is meant to aid in the visualization of image data
% as is expected of image data.

cFigure;
subplot(2,2,1);
title('MATLAB imagesc function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
imagesc(M(:,:,S));
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

subplot(2,2,2);
title('MATLAB slice function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
slice(M,[],[],S);
axisGeom(gca,fontSize);
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

subplot(2,2,3);
title('im2patch function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C);
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

subplot(2,2,4);
title('MATLAB pcolor function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
pcolor(M(:,:,S));
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

drawnow;

%%
% Note that the image size is wrong for the |slice| and |pcolor| commands.
% This is because intensities appear to be defined on voxels vertices (and
% are reinterpolated onto faces) for these functions instead of voxels
% centres as should be the case for image data.

%% Example: Creating and plotting combined voxel and slice patch data

%%
% Simulating 3D image
[X,Y,Z]=meshgrid(linspace(-4.77,4.77,25));
phi=(1+sqrt(5))/2;
M=1/6*(2 - (cos(X + phi*Y) + cos(X - phi*Y) + cos(Y + phi*Z) + cos(Y - phi*Z) + cos(Z - phi*X) + cos(Z + phi*X)));

%%
% Creating and plotting patch data. Last example illustrates the use of a
% specific logic description (i.e. a mask) for the voxels of interest

cFigure;
title('Combined voxel and slice plotting');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;

% Setting up indices for I direction slices
S=round(size(M,1)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(S,:,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'si'); %Creating patch data for y mid-voxel slices
gpatch(F,V,C,'none',faceAlpha2);

% Setting up indices for J direction slices
S=round(size(M,2)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,S,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sj'); %Creating patch data for x mid-voxel slices
gpatch(F,V,C,'none',faceAlpha2);

% Setting up indices for K direction slices
S=round(size(M,3)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,:,S)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sk'); %Creating patch data for z mid-voxel slices
gpatch(F,V,C,'none',faceAlpha2);

% Setting up indices for voxels to plot
L_mask=M>-0.2 & M<0;
[F,V,C]=im2patch(M,L_mask,'vb'); %Creating patch data for selection of high voxels

gpatch(F,V,C);
colormap(cMap); colorbar; caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Example: Shrinking patch data through combination with |scalePatch| function

% Using the function |scalePatch| patch data can be shrunk to aid
% visualisation (can be a tool to avoid memory costly transparency for
% instance).

[Ev,Vv,Cv]=im2patch(M,L_mask,'h'); %This creates a hexahedral element for each voxel
[Evs,Vvs]=scalePatch(Ev,Vv,0.5); %Apply voxel element scaling
[Fvs,Cvs]=element2patch(Evs,Cv); %Convert to quad faces for plotting

%%
% Plotting the voxels
cFigure;
title('Scaled (shrunk) patch data');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(Fvs,Vvs,Cvs);
colormap(cMap); colorbar; caxis([min(Cvs(:)) max(Cvs(:))]);
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Example: Medical image data and coordinate manipulation due to voxel size

% Get a 3D image
load mri;
M=squeeze(D); %example image data set
v=2./[1,1,.4]; %example voxel size

%%
% The voxels to display can be specified as a list (vector) of voxels
% numbers (linear indices) or using a mask (logic array).

%Defining row, column and slice indicices for slice patching
sliceIndexI=round(size(M,1)/2); %(close to) middle row
sliceIndexJ=round(size(M,2)/2); %(close to) middle column
sliceIndexK=round(size(M,3)/2); %(close to) middle slice

%Defining "masks" i.e. logic arrays with ones for voxels of interest
logicSliceI=false(size(M));
logicSliceI(sliceIndexI,:,:)=1;
logicSliceI=logicSliceI & M>0;

logicSliceJ=false(size(M));
logicSliceJ(:,sliceIndexJ,:)=1;
logicSliceJ=logicSliceJ & M>0;

logicSliceK=false(size(M));
logicSliceK(:,:,sliceIndexK)=1;
logicSliceK=logicSliceK & M>0;

%Defining voxel indices for voxels of interest
T_low=min(M(:))+((max(M(:))-min(M(:)))/10); %Threshold example
logicVoxels=(M>T_low);
logicVoxels(:,1:sliceIndexJ,:)=0;

%%
% Creating patch data
% The patch data consists of a matrix array defining the faces, a matrix
% array defining the vertices and a vector for the colour data. The
% vertices are based on the image coordinates however they are formatted
% as: [X(:) Y(:) Z(:)]. X relates to columns, Y to rows and Z to slices.
% Use a function like |im2cart| , or |im2mrcart| to convert image to cartesian
% coordinates.

[Fv,Vv,Cv]=im2patch(M,logicVoxels,'vb');
[Fx,Vx,Cx]=im2patch(M,logicSliceJ,'sj');
[Fy,Vy,Cy]=im2patch(M,logicSliceI,'si');
[Fz,Vz,Cz]=im2patch(M,logicSliceK,'sk');

% Convert image coordinates to cartesian coordinates
[Vv(:,1),Vv(:,2),Vv(:,3)]=im2cart(Vv(:,2),Vv(:,1),Vv(:,3),v);
[Vx(:,1),Vx(:,2),Vx(:,3)]=im2cart(Vx(:,2),Vx(:,1),Vx(:,3),v);
[Vy(:,1),Vy(:,2),Vy(:,3)]=im2cart(Vy(:,2),Vy(:,1),Vy(:,3),v);
[Vz(:,1),Vz(:,2),Vz(:,3)]=im2cart(Vz(:,2),Vz(:,1),Vz(:,3),v);

%%
%
cFigure;
title('MRI visualisation, slices and voxels in cartesian coordinates with aid of voxel size');
xlabel('X (mm)');ylabel('Y (mm)'); zlabel('Z (mm)'); hold on;
gpatch(Fv,Vv,Cv,edgeColor1,faceAlpha2);
gpatch(Fx,Vx,Cx,edgeColor1,faceAlpha1);
gpatch(Fy,Vy,Cy,edgeColor1,faceAlpha1);
gpatch(Fz,Vz,Cz,edgeColor1,faceAlpha1);
colormap(gray(250)); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Example: Combining colormap and RGB driven patch colours
% N.B. The figure renderer might have to be set to OPENGL.

%%
% Convert voxels colouring to RGB type, here a simple conversion to a gray
% scale description is used
Cv=(Cv*ones(1,3))./max(Cv(:));

%%
% Plotting the voxels

cFigure;
title('MRI visualisation, slices and voxels, colormap and RGB driven respectively');
xlabel('X (mm)');ylabel('Y (mm)'); zlabel('Z (mm)'); hold on;

%RGB driven
gpatch(Fv,Vv,Cv,edgeColor1,faceAlpha2);

%Colormap driven
gpatch(Fx,Vx,Cx,edgeColor1,faceAlpha1);
gpatch(Fy,Vy,Cy,edgeColor1,faceAlpha1);
gpatch(Fz,Vz,Cz,edgeColor1,faceAlpha1);

colormap(cMap); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function pix = pixelize(datain_all, intens, box, nx, ny, namefile, showplot)

% pix = pixelize(datain_all, intens, box, nx, ny, namefile,showplot)
% pixelize PALM data
% datain_all - input data, #points-by-2 matrix, each row correspond to xy
% coordinate of the datapoint
% intens - intensities of each center (stored in a0_phot)
% box = [xlim1, xlim2, ylim1, ylim2] - selected ROI
% nx - number of pixels in x direction
% ny - number of pixels in y direction
% namefile - (optional) name of the TXT file with matrix pix
if ~exist ('showplot', 'var')
    showplot = 1;
end

sx = (box(2)-box(1))/nx; %size of the pixels
sy = (box(4)-box(3))/ny;

[datain, indexGood] = ROIdata (datain_all, box(1), box(2), box(3),  box(4), showplot);
intensGood=intens(indexGood);
pixvec = ceil ([(datain(:,1) - box(1))/sx , (datain(:,2) - box(3))/sy ]);

% indrem_min = find (or(pixvec(:,1) < 0, pixvec(:,2) < 0));
% indrem_max = find (or(pixvec(:,1) > nx, pixvec(:,2) > ny));

pix = zeros(nx, ny);
ind = sub2ind(size(pix), pixvec(:,1), pixvec(:,2));
% pix(ind) = intens;
for ii=1:length(ind)
    pix(ind(ii))=pix(ind(ii))+intensGood(ii);
end

% pix = flipud(pix'); % to follow the matlab notation for imagesc(pix)...
pix = pix'; % to follow the matlab notation for imagesc(pix)...

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

if showplot
    figure
    imagesc(pix)
    set (gca, 'DataAspectRatio',[sy sx 1]);
end

if exist('namefile', 'var')
    if ~isempty(namefile)
        fid = fopen([namefile '.txt'], 'wt'); % Open for writing
        for i=1:size(pix,1)
            fprintf(fid, '%d ', pix(i,:));
            fprintf(fid, '\n');
        end
        fclose(fid);
    end
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function v = getfield(s,field)
%GETFIELD Get mutable structure field contents.
% f = getfield(s,'field') returns the contents of the specified field.
% This is equivalent to the syntax f = s.field.
%
% See also SETFIELD, RMFIELD, ISFIELD, FIELDNAMES.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

v = s.obj.get(field);
if isempty(v)
  error(sprintf('Reference to non-existent field ''%s''.',field));
end
v = fromJava(v);

author: andrewssobral

repo: lrslibrary


text:
function y = real( x )

%   Disciplined convex/geometric programming information for REAL:
%      REAL(X) may be freely applied to any CVX expression. However,
%      since REAL(X)=X for all real expressions (including convex,
%      concave, log-convex, and log-concave), it is only useful in
%      the complex affine case.

y = cvx( x.size_, real( x.basis_ ) );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function res = calc_pot_values(w)
%
%calculate potential values for non obstacle and foal locations in given world
%
r = size(w,1);
c = size(w,2);
for t = size(w,3):-1:1
    for x = 1:r
        for y = 1:c
            if w(x,y,t) == -1
                n = w(x,y,t+1);
                if x ~= 1
                    n = n + w(x-1,y,t+1);
                else
                    n = n + 1;% boundary condition
                end
                if x ~= r
                    n = n + w(x+1,y,t+1);
                else
                    n = n + 1;
                end
                if y ~= 1
                    n = n + w(x,y-1,t+1);
                else
                    n = n + 1;
                end
                if y ~= c
                    n = n + w(x,y+1,t+1);
                else
                    n = n + 1;
                end
                w(x,y,t) = n / 5;% take average of neighbors
            end
        end
    end
end
res = w;
end

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_nufft( pref )

if ( nargin == 0 )
    pref = chebfunpref();
end

% Choose a tolerance:
tol = eps;
rng(0)

% Test NUFFT-I
count = 1;
for N = 10.^(0:4)
    omega = (0:N-1)';
    omega = omega + 1.2*rand(N,1)/N;
    c = rand(N,1) + 1i*rand(N,1);
    exact = nudft1( c, omega );
    fast = chebfun.nufft( c, omega, 1 );
    pass(count) = norm( exact - fast, inf ) < 300*N*tol*norm(c,1);
    count = count + 1;
end

% Test on random inputs:
N = 50;
omega = N*randn(1,N);
c = rand(N,1) + 1i*rand(N,1);
F = exp(-2*pi*1i*((0:N-1)/N).'*omega);
f = chebfun.nufft(c, omega.', 1);
pass(count) = ( norm( f - F*c ) < 300*N*tol*norm(c,1) );
count = count + 1;

% Test NUFFT-II:
for N = 10.^(0:4)
    x = linspace(0,1,N+1)'; x(end) = [];
    x = x + 1.2/N;
    c = rand(N,1) + 1i*rand(N,1);
    exact = nudft2( c, x );
    fast = chebfun.nufft( c, x );
    pass(count) = norm( exact - fast, inf ) < 300*N*tol*norm(c,1);
    count = count + 1;
end

% Test on random inputs:
N = 50;
x = 100*rand(N,1);
c = rand(N,1) + 1i*rand(N,1);
F = exp(-2*pi*1i*x*(0:N-1));
f = chebfun.nufft(c, x);
pass(count) = ( norm( f - F*c ) < 300*N*tol*norm(c,1) );
count = count + 1;

% Test NUFFT-III
for N = 10.^(0:4)
    omega = (0:N-1)';
    omega = omega + 1.2*rand(N,1)/N;
    x = linspace(0,1,N+1)'; x(end) = [];
    x = x + 1.2/N;
    c = rand(N,1) + 1i*rand(N,1);
    exact = nudft3( c, x, omega );
    fast = chebfun.nufft( c, x, omega, 3 );
    pass(count) = norm( exact - fast, inf ) < 10000*N*tol*norm(c,1);
    count = count + 1;
end

% Test NUFFT-II on nonsquare inputs:
n = 101;
m = 3;
x = linspace(0,1,m+1)'; x(end) = [];
x = x + 1.01/m;
c = rand(n,1);
exact = nudft2( c, x );
fast = chebfun.nufft( c, x );
pass(count) = norm( exact - fast, inf ) < 100*tol*norm(c,1);
count = count + 1;

end

function f = nudft1( c, omega)

f = zeros(size(omega,1),1);
for j = 1:numel(f)
    f(j) = exp(-2*pi*1i*(j-1)/size(omega,1)*omega.')*c;
end
end

function f = nudft2( c, x)

f = zeros(size(x,1),1);
omega = 0:size(c,1)-1;
for j = 1:numel(f)
    f(j) = exp(-2*pi*1i*x(j)*omega)*c;
end
end

function f = nudft3( c, x, omega)

f = zeros(size(omega,1),1);
for j = 1:numel(f)
    f(j) = exp(-2*pi*1i*x(j)*omega.')*c;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% test for computations of nn

name = 'swissroll';
n = 3000;
options.sampling = 'rand';
[X,col] = load_points_set( name, n, options );
nbr_nn = 50;

options.use_nntools = 0;
[D1,nn_list1] = compute_nn_distance(X,nbr_nn, options);
options.use_nntools = 1;
[D2,nn_list2] = compute_nn_distance(X,nbr_nn, options);

disp( ['Error (should be 0)=' num2str(norm(D1-D2) + norm(nn_list1-nn_list2),2) '.'] );

L = nn_list1(1,:);
col(L) = 2;

plot_scattered(X,col);

author: gpeyre

repo: matlab-toolboxes


text:
function imResult = blendMode_Lighten(A, B, offsetW, offsetH)
%% Lighten blending mode: compares the color information for each pixel of
%   the base and the blend color and applies the lighter color as the
%   result.
%   Any pixels in the base image that are darker than the blend color are
%   replaced, and pixels that are lighter are left unchanged. No part of
%   the image will become darker.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_Lighten));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

ind = B > A;
C = abs(ind - 1) .* A + ind .* B;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end


text:
function [ a, b ] = p16_lim ( )

%*****************************************************************************80
%
%% P16_LIM returns the integration limits for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rx = rx3ViewToggle(rx, nToggles, views);
% Toggle between the interpolated and reference slices for the 3 view axes
% in the interpolated 3-view figure.
%
%  rx = rx3ViewToggle(rx, [nToggles=3], [views=1:3]);
%
%
% ras, 05/05/2009.
if ~exist('rx', 'var') | isempty(rx)
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

if notDefined('nToggles')
    % # of times to alternate b/w interp and ref images
    nToggles = 3;
end

if notDefined('views')
    % which of the 3 view windows to toggle
    views = 1:3;
end

%% get position from UI controls
for i = 1:3
    loc(i) = str2num( get(rx.ui.interpLoc(i), 'String') );
end

for ori = views
    % get interpolated image
    interp = rxInterpSlice(rx, loc(ori), ori);

    % get reference image
    switch ori
        case 1, ref = squeeze( rx.ref(loc(ori),:,:) );
        case 2, ref = squeeze( rx.ref(:,loc(ori),:) );
        case 3, ref = rx.ref(:,:,loc(ori));
    end

    % adjust contrast, brightness, etc.
    [interp ref] = rxGetComparisonImages(rx, interp, ref);

    if ishandle(rx.ui.interp3ViewAxes(ori))
        axes(rx.ui.interp3ViewAxes(ori))
    else
        figure('Name', ['mrRx: ' mfilename]);
    end

    rng = [min(interp(:)) max(interp(:))];

    % get cur axis limits and keep them, in case user has zoomed image
    x = get(gca, 'xlim');
    y = get(gca, 'ylim');

    for n = 1:nToggles
        imshow(ref, rng);
        xlim(x); ylim(y);
        pause(0.15);
        imshow(interp, rng);
        xlim(x); ylim(y);
        pause(0.15);
    end

    % if n toggles is 0, just switch the view to the reference image
    if nToggles < 1,
        imshow(ref, rng);
        xlim(x); ylim(y);
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_PA_10_6DOF(robot, T)
%   Solves the inverse kinematic problem for the MITSUBISHI PA-10 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_PA_10_6DOF returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('MISUBISHI', 'pa-10');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic( robot, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q = inversekinematic_pa_10_6DOF(robot, T)

%initialize, eight possible solutions
q=zeros(6,8);

%theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
%alpha = eval(robot.DH.alpha);

%L1=d(1);
%L2=a(2);
%L3=d(4);
L6=d(6);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Calculamos la posici�n de la mu�eca. W es el tercer vector en el extremo
%del brazo
W = T(1:3,3);
%Pm: posici�n de la mu�eca
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);

q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);

q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

%leave only the real part of the solutions
q=real(q);

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));
q(3,:) = normalize(q(3,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    qtemp = solve_spherical_wrist_pa_10(robot, q(:,i), T, 1); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    qtemp = solve_spherical_wrist_pa_10(robot, q(:,i), T, -1); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=a(2);
L3=d(4);

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(p1(2), p1(1));
gamma = (acos((L2^2+r^2-L3^2)/(2*r*L2)));
if ~isreal(gamma)
    disp('WARNING:inversekinematic_pa_10: the point is not reachable for this configuration, imaginary solutions');
    %gamma = real(gamma);
end
%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = beta + gamma-pi/2; %elbow up
q2(2) = beta - gamma-pi/2; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
%theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
%alpha = eval(robot.DH.alpha);

%See geometry
L2=a(2);
L3=d(4);

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = acos((L2^2 + L3^2 - r^2)/(2*L2*L3));

if ~isreal(eta)
   disp('WARNING:inversekinematic_pa_10: the point is not reachable for this configuration, imaginary solutions');
   %eta = real(eta);
end
%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = -pi + eta;
q3(2) = -eta + pi;

% Solve the special case of this spherical wrist
% For wrists that whose reference systems have been placed as in the
% ABB IRB 140--> use solve_spherical_wrist2
% For wrists with the same orientation as in the KUKA KR30_jet
%--> use solve_spherical_wrist
function q = solve_spherical_wrist_pa_10(robot, q, T, wrist)

% T is the noa matrix defining the position/orientation of the end
% effector's reference system
vx6=T(1:3,1);
vz5=T(1:3,3); % The vector a z6=T(1:3,3) is coincident with z5

% Obtain the position and orientation of the system 3
% using the already computed joints q1, q2 and q3
T01=dh(robot, q, 1);
T12=dh(robot, q, 2);
T23=dh(robot, q, 3);
T03=T01*T12*T23;

vx3=T03(1:3,1);
vy3=T03(1:3,2);
vz3=T03(1:3,3);

% find z4 normal to the plane formed by z3 and a
vz4=cross(vz3, vz5);	% end effector's vector a: T(1:3,3)

% in case of degenerate solution,
% when vz3 and vz6 are parallel--> then z4=0 0 0, choose q(4)=0 as solution
if norm(vz4) <= 0.00000001
    if wrist == 1 %wrist up
        q(4)=0;
    else
        q(4)=-pi; %wrist down
    end
else
    %this is the normal and most frequent solution
    cosq4=wrist*dot(vy3,vz4);
    sinq4=wrist*dot(-vx3,vz4);
    q(4)=atan2(sinq4, cosq4);
end
%propagate the value of q(4) to compute the system 4
T34=dh(robot, q, 4);
T04=T03*T34;
vx4=T04(1:3,1);
vy4=T04(1:3,2);

% solve for q5
cosq5=dot(-vy4,vz5);
sinq5=dot(vx4,vz5);
q(5)=atan2(sinq5, cosq5);

%propagate now q(5) to compute T05
T45=dh(robot, q, 5);
T05=T04*T45;
vx5=T05(1:3,1);
vy5=T05(1:3,2);

% solve for q6
cosq6=dot(vx6,vx5);
sinq6=dot(vx6,vy5);
q(6)=atan2(sinq6, cosq6);

repo: ARTE


text:
function b = isodd(x)
%ISODD  True for odd numbers.
%
%   ISODD(X) returns 1's where the elements of X are odd numbers and 0's where
%   they are not.
%
%   See also ISINT, ISEVEN.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-04-12 14:28:40 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 1, nargin));
   if ~isnumeric(x)
      error('Argument must be a numeric array.');
   end

   cls = class(x);                      % class of input argument
   if isempty(x)
      b = feval(cls, x);                % return empty array of same class
   else
      switch cls
         case 'double'
            b = mod(x, 2) == 1;
         case 'single'
            % "mod" is not defined for class "single"; so convert input to
            % double, compare, and convert back
            b = single(mod(double(x), 2) == 1);
         case {'uint8', 'uint16', 'uint32', 'uint64'}
            b = bitand(x, 1);
         case {'int8', 'int16', 'int32', 'int64'}
            error('Not implemented for classes int8, int16, int32, and int64.');
         otherwise
            error('Argument is of unrecognized class.');
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
x = uiuc_sample;
x = x(1:256, 1:256);

%% compute roto-translaiton scattering
filt_opt.J = 4;
filt_rot_opt.null = 1;
scat_opt.oversampling = 0;
Wop = wavelet_factory_3d(size(x), filt_opt, filt_rot_opt, scat_opt);

tic;
S_rt = scat(x, Wop);
toc;

%% arrange scattering layer in lexicographic order
var_y{1}.name = 'j';
var_y{1}.index = 1;
var_y{2}.name = 'j';
var_y{2}.index = 2;
var_x{1}.name = 'k2';
var_x{1}.index = 1;
var_x{2}.name = 'theta2';
var_x{2}.index = 1;

big_img = image_scat_layer_order(S_rt{3},var_x, var_y,0);
figure(2);
imagesc(big_img);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [ node_coord, node_att, node_marker ] = triangle_node_data_read ( ...
  node_filename, node_num, node_dim, node_att_num, node_marker_num )

%*****************************************************************************80
%
%% TRIANGLE_NODE_DATA_READ reads the data from a node file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
% Parameters:
%
%   Input, string NODE_FILENAME, the name of the node file.
%
%   Input, integer NODE_NUM, the number of nodes.
%
%   Input, integer NODE_DIM, the spatial dimension.
%
%   Input, integer NODE_ATT_NUM, number of node attributes
%   listed on each node record.
%
%   Input, integer NODE_MARKER_NUM, 1 if every node record
%   includes a final boundary marker value.
%
%   Output, real NODE_COORD(NODE_DIM,NODE_NUM), the nodal
%   coordinates.
%
%   Output, real NODE_ATT(NODE_ATT_NUM,NODE_NUM), the nodal
%   attributes.
%
%   Output, integer NODE_MARKER(NODE_MARKER_NUM,NODE_NUM), the
%   node markers.
%
  node = 0;

  input = fopen ( node_filename, 'rt' );

  if ( input < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_NODE_DATA_READ - Fatal error!\n' );
    fprintf ( 1, '  Unable to open file.\n' );
    error ( 'TRIANGLE_NODE_DATA_READ - Fatal error!' );
  end
%
%  Create format.
%
  format = '%d';
  for i = 1 : node_dim
    format = strcat ( format, '%g' );
  end
  for i = 1 : node_att_num
    format = strcat ( format, '%g' );
  end
  for i = 1 : node_marker_num
    format = strcat ( format, '%d' );
  end
%
%  Create arrays.
%
  node_coord = zeros ( node_dim, node_num );
  node_att = zeros ( node_att_num, node_num );
  node_marker = zeros ( node_marker_num, node_num );
%
%  Read values and place in the arrays.
%
  while ( 1 )

    text = fgets ( input );

    if ( text == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TRIANGLE_NODE_DATA_READ - Fatal error!\n' );
      fprintf ( 1, '  Unexpected end of file while reading.\n' );
      error ( 'TRIANGLE_NODE_DATA_READ - Fatal error!' );
    end

    if ( s_len_trim ( text ) == 0 )
      continue
    end

    if ( text(1) == '#' )
      continue
    end
%
%  Ignore the dimension line.
%
    if ( node == 0 )

    else

      [ value, count ] = sscanf ( text, format );

      if ( count < 1 + node_dim + node_att_num + node_marker_num )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'TRIANGLE_NODE_DATA_READ - Fatal error!\n' );
        fprintf ( 1, '  Line did not have enough data.\n' );
        error ( 'TRIANGLE_NODE_DATA_READ - Fatal error!' );
      end

      j = 1;

      for i = 1 : node_dim
        j = j + 1;
        node_coord(i,node) = value(j);
      end

      for i = 1 : node_att_num
        j = j + 1;
        node_att(i,node) = value(j);
      end

      for i = 1 : node_marker_num
        j = j + 1;
        node_marker(i,node) = value(j);
      end

    end

    node = node + 1;

    if ( node_num < node )
      break
    end

  end

  fclose ( input );

  return
end


text:
function PaV = calc_Pa(H0,betaV,dkV)
%Actuarial probability of complication using Cox model

PaV = 1 - exp(-H0 .* exp(betaV.*dkV));

%Test-1:
% tol = 10^-5;
% expectedp10 = 0.1;
% H0=0.01;
% B = 0.05;
% d10 = 47.096;
% p10 = calc_Pa(H0,B,d10);
% assertAlmostEqual(p10,expectedp10,tol);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [in] = au2in(au)
% Convert length from astronomical units to inches.
% Chad A. Greene 2012
in = au*5889679948464;


text:
function rs = acf(s, len)

%tstoolbox/@signal/acf
%   Syntax:
%     * acf(s, len)
%
%   Input arguments:
%     * len -length of the fft (optional)
%
%   Autocorrelation function for real scalar signals, using fft (of length
%   len). If len is ommited a default value is calculated. The maximum of
%   the calculated length is 128.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,2);

if (ndim(s) > 1) | (~isreal(data(s)))
	help(mfilename)
	return
end

if (nargin < 2)
	if dlens(s,1) > 256
		len = 128;
	else
		len = nextpow2(dlens(s,1)/4);
	end
end

c = acf(s.core, len);
rs = signal(c, s);	% special constructor calling syntax for working routines
a = getaxis(s, 1);
dl = delta(a);
a = setfirst(a, 0);
rs = setaxis(rs, 1, a);
rs = setyunit(rs, unit);		% acf values are scalars without unit
rs = addhistory(rs, 'Autocorrelation function');
rs = setlabel(rs, 'Autocorrelation function');
rs = addcommandlines(rs, 's = acf(s', len);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef PatternSearchUseCacheEnum < matlab.mixin.SetGet
    %PatternSearchUseCacheEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        UseCache('Use Cache',true)
        DoNotUseCache('Do Not Use Cache',false);
    end

    properties
        name char = '';
        optionVal(1,1) logical = true
    end

    methods
        function obj = PatternSearchUseCacheEnum(name, optionVal)
            obj.name = name;
            obj.optionVal = optionVal;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('PatternSearchUseCacheEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('PatternSearchUseCacheEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('PatternSearchUseCacheEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Semi-Soft GoDec (Zhou and Tao 2011)
% process_video('RPCA', 'SSGoDec', 'dataset/cctv.avi', 'output/demo_SSGoDec.avi');
rank = 1;
tau = 8;
power = 0;
L = SSGoDec(M,rank,tau,power);
S = M - L;


text:
function theta = angle3Points(varargin)
%ANGLE3POINTS Compute oriented angle made by 3 points.
%
%   ALPHA = angle3Points(P1, P2, P3);
%   Computes the angle between the points P1, P2 and P3.
%   Pi are either [1*2] arrays, or [N*2] arrays, in this case ALPHA is a
%   [N*1] array. The angle computed is the directed angle between line
%   (P2P1) and line (P2P3).
%   Result is always given in radians, between 0 and 2*pi.
%
%   See Also:
%   points2d, angles2d, angle2points
%
%
%   ---------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% created the 23/02/2004.
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY :
%   25/09/2005 : enable single parameter

if length(varargin)==3
    p1 = varargin{1};
    p2 = varargin{2};
    p3 = varargin{3};
elseif length(varargin)==1
    var = varargin{1};
    p1 = var(1,:);
    p2 = var(2,:);
    p3 = var(3,:);
end

% angle line (P2 P1)
theta = lineAngle(createLine(p2, p1), createLine(p2, p3));


text:
function DCM = spm_dcm_phase(DCM)
% Estimate parameters of a DCM model of phase-coupled responses
% FORMAT DCM = spm_dcm_phase(DCM)
%
% DCM
%    name: name string
%       M:  Forward model
%              M.dipfit - leadfield specification
%       xY: data   [1x1 struct]
%       xU: design [1x1 struct]
%
%   Sname: cell of source name strings
%
%   Connection constraints:
%
%       A: {[nr x nr double] }
%       B: {[nr x nr double]}   for GUI specification
%                               (Nfourier=1 & only sine terms)
%   or
%
%       As: [nr x nr x Nfourier]
%       Ac: [nr x nr x Nfourier]
%       Bs: [nr x nr x Nfourier]
%       Bc: [nr x nr x Nfourier]   for script specification
%
%   options.type         - 'ECD'
%__________________________________________________________________________
% Copyright (C) 2009 Wellcome Trust Centre for Neuroimaging
%
% Will Penny
% $Id: spm_dcm_phase.m 4884 2012-09-03 13:33:17Z guillaume $

% check options
%==========================================================================
clear spm_erp_L

% Filename and options
%--------------------------------------------------------------------------
try, DCM.name;                  catch, DCM.name           = 'DCM_PHASE'; end
try, DCM.options.Nmodes;        catch, DCM.options.Nmodes = 4;         end
try, h     = DCM.options.h;     catch, h                  = 1;         end
try, onset = DCM.options.onset; catch, onset              = 80;        end

% Data and spatial model
%==========================================================================
DCM    = spm_dcm_erp_dipfit(DCM, 1);

% Get data if not gotten already
%if ~isfield(DCM.xY,'source')
    DCM  = spm_dcm_phase_data(DCM);
%end

disp('Estimating Model ...');

xY     = DCM.xY;
xU     = DCM.xU;
xU.dt  = xY.dt;

if isfield(DCM,'A')
    % If using UI, copy model structures into As and Bs matrices
    DCM.As=DCM.A{1};
    DCM.Bs{1}=DCM.B{1};
end

% dimensions
%--------------------------------------------------------------------------
Nt     = length(xY.y);                  % number of trials
Nr     = size(DCM.As,1);                % number of sources
nx     = Nr;                            % number of states

Ns     = size(xY.y{1},1);               % number of samples
xU.u   = zeros(Ns,1);

Nu     = size(xU.X,2);                  % number of trial-specific effects

% Set error covariance components - one per region
%--------------------------------------------------------------------------
xY.Q   = spm_Ce(Ns*Nt*ones(1,Nr));

% Offsets
xY.X0  = sparse(Ns,0);

% Inputs
%==========================================================================

% trial-specific effects
%--------------------------------------------------------------------------
try
    if length(DCM.Bs) ~= Nu;
        warndlg({'Please ensure number of trial specific effects', ...
                 'encoded by DCM.xU.X & DCM.B are the same'})
    end
catch
    DCM.Bs = {};
end

% model specification and nonlinear system identification
%==========================================================================
M      = DCM.M;

M.freq = mean(DCM.options.Fdcm);
M.fb = 0.5*(DCM.options.Fdcm(2)-DCM.options.Fdcm(1));

try, M = rmfield(M,'g');  end
try, M = rmfield(M,'FS'); end

if ~isfield(M,'dipfit')
    M.dipfit=[];
end

% prior moments
%--------------------------------------------------------------------------
[pE,gE,pC,gC] = spm_phase_priors(DCM,M.fb,M.dipfit);

% likelihood model
%--------------------------------------------------------------------------
M.IS  = 'spm_gen_phase';
M.f   = 'spm_fx_phase';
M.G   = 'spm_lx_phase';

M.pE  = pE;
M.pC  = pC;
M.gE  = gE;
M.gC  = gC;
M.n   = nx;
M.l   = Nr;
M.ns  = Ns;

% Don't plot progress
try, M.nograph; catch, M.nograph=0; end

% Initial state
try
    M.x;
catch
    M.x=zeros(Nr,1);
end

% Set up initial phases
for n=1:length(DCM.xY.y)
    xx=DCM.xY.y{n}(1,:);
    M.trial{n}.x = double(xx(:)');
end

% keyboard
% myx=spm_gen_phase(pE,M,xU);

% EM: inversion
%--------------------------------------------------------------------------
[Qp,Qg,Cp,Cg,Ce,F] = spm_nlsi_N(M,xU,xY);

% Data ID
%==========================================================================
if isfield(M,'FS')
    try
        ID  = spm_data_id(feval(M.FS,xY.y,M));
    catch
        ID  = spm_data_id(feval(M.FS,xY.y));
    end
else
    ID  = spm_data_id(xY.y);
end

% Bayesian inference {threshold = prior} NB Prior on A,B  and C = exp(0) = 1
%==========================================================================
warning('off','SPM:negativeVariance');
dp  = spm_vec(Qp) - spm_vec(pE);
Pp  = spm_unvec(1 - spm_Ncdf(0,abs(dp),diag(Cp)),Qp);
warning('on','SPM:negativeVariance');

% neuronal and sensor responses (x and y)
%--------------------------------------------------------------------------
x   = feval(M.IS,Qp,M,xU);        % prediction (source space)
L   = feval(M.G, Qg,M);           % get gain matrix

% trial-specific responses (in mode, channel and source space)
%--------------------------------------------------------------------------
for i = 1:Nt
    s    = x{i};                  % prediction (source space)
    y{i} = s*L';                  % prediction (sensor space)
    %r  = xY.y{i} - y;            % residuals  (sensor space)

end

% Change back design matrix to user specified
xU.X=DCM.xU.oldX;

% store estimates in DCM
%--------------------------------------------------------------------------
DCM.M  = M;                    % model specification
DCM.xY = xY;                   % data structure
DCM.xU = xU;                   % input structure
DCM.Ep = Qp;                   % conditional expectation f(x,u,p)
DCM.Cp = Cp;                   % conditional covariances G(g)
DCM.Eg = Qg;                   % conditional expectation
DCM.Cg = Cg;                   % conditional covariances
DCM.Pp = Pp;                   % conditional probability
DCM.Ce = Ce;                   % conditional error covariance
DCM.F  = F;                    % Laplace log evidence
DCM.ID = ID;                   % data ID
DCM.y  = y;                    % Model predictions

% and save
%--------------------------------------------------------------------------
save(DCM.name, 'DCM', spm_get_defaults('mat.format'));
assignin('base','DCM',DCM)

author: spm


text:
function x = gendatav(pc,R0,nobs)
%function x = gendatav(pc,R0,nobs)
%
%  Simulates nobs observations of a vector AR-process
%  with partial correlations pc and covariance matrix R0.
%
%  See also FILTERV, RANDN.

%  S. de Waele, FEB 2001.

if ~isstatv(pc), error('Partial correlations non-stationairy!'), end
order = kingsize(pc,3)-1;
dim = size(pc,1);
I = eye(dim);
[rc,rcb] = pc2rcv(pc,R0);
[Pf,Pb] = pc2resv(pc,R0);
par = zeros(dim,dim,order+1);
parb = zeros(dim,dim,order+1);

x = zeros(dim,1,nobs); %The signal in matrix notation
%First observation
innovation = randncov(Pf(:,:,1));
x(:,:,1) = innovation;

par(:,:,1) = I;
parb(:,:,1)= I;
if order,
	par(:,:,2) = rc(:,:,2);
	parb(:,:,2)= rcb(:,:,2);
	par_o  = par;
	parb_o = parb;
end

%observation 2 to order
for obs = 2:order,
   innovation = randncov(Pf(:,:,obs));
   prediction = -prodsumv(flipdim(par(:,:,2:obs),3),x(:,:,1:obs-1));
   x(:,:,obs) = prediction+innovation;

   p=obs; %The next observations requires the AR(obs)-model
   %par(:,:,2:p) =  par_o(:,:,2:p) +flipdim(filterv(rc(:,:,p+1),1,parb_o(:,:,2:p)),3);
   par(:,:,2:p) =  par_o(:,:,2:p) +flipdim(armafilterv(parb_o(:,:,2:p),1,rc(:,:,p+1)),3);
   par(:,:,p+1)= rc(:,:,p+1);
   parb(:,:,2:p) =  parb_o(:,:,2:p) +flipdim(armafilterv(par_o(:,:,2:p),1,rcb(:,:,p+1)),3);
   parb(:,:,p+1)= rcb(:,:,p+1);

   par_o  = par;
   parb_o = parb;

end %for obs = 2:order,

%observations order+1 to nobs
innovation = zeros(dim,1,nobs-order);
for i = 1:nobs-order,
   innovation(:,:,i) = randncov(Pf(:,:,order+1));
end
x(:,:,order+1:end) = armafilterv(innovation,par,I,x(:,:,1:order));


text:
function pass = test_intops(pref)
% Test integral operators
% Toby Driscoll 28 May 2009
% Nick Hale 6 Jan 2012

% NOTE: Taken from V4 test chebop_intops.

if ( nargin == 0 )
    pref = cheboppref();
end
tol = 1e-10;

% From http://en.wikipedia.org/wiki/Integro-differential_equation
d = [0 5];
x = chebfun('x', d);
N = chebop(@(u) diff(u) + 2*u + 5*cumsum(u), d);
N.bc = @(x,u) u(0);
u = mldivide(N, 1, pref);
u_exact = 0.5*exp(-x).*sin(2*x);
err(1) = norm(u - u_exact);

% Fredholm
d = [0, 1];
x = chebfun(@(x) x, d);
K = @(x,y) sin(2*pi*(x-y));
A = chebop(@(x,u) u + fred(K,u), d);
u = x.*exp(x);
f = A*u;
err(2) = norm(u - mldivide(A, f, pref));

% Volterra
d = [0, pi];
x = chebfun(@(x) x, d);
K = @(x,y) x.*y;
A = chebop(@(x,u) u - volt(K,u), d);
f = x.^2.*cos(x) + (1-x).*sin(x);
u = mldivide(A, f, pref);
err(3) = norm( u - sin(x) );
err(4) = norm( A*u - f );

pass = err < tol;

%%
% From #2122
K = @(x,t) exp(t.*(x-t));
N = chebop(@(x,y) diff(y) + y - x.*(1+2*x).*volt(K, y) - 1 - 2*x);
N.lbc = 1;
y = N\0;
pass(5) = norm(N.op(y)) + abs(y(0)-1);

%%
% From #2179
K = @(x,y) exp(-(x-y));
A = chebop(@(x,u) diff(u) + fred(K,u));
A.lbc = 0;
u = A\1;
pass(6) = norm(A*u - 1, inf) < tol;

%%
% From https://groups.google.com/forum/#!topic/chebfun-users/0dpsggp1RyA .
% Also, see #2210.

L = chebop(@(u) diff(u)+sum(u), [0 1]);
L.lbc = 1;
u = L\0;
pass(7) = norm(u - (1- 2*chebfun('x',[0 1])/3), inf) < tol;

%%
% From #2185
A = chebop(@(x,u) diff(u,2) + sin(2*pi*x)*u + sum(u).*u);
A.bc = 'periodic';
u = A\1;
pass(8) = norm(A*u-1) < tol;

f = chebfun(@(x) sin(5*pi*x));
K = @(x,y) cos(-pi*(x-y));
A = chebop(@(x,u) 1/25*diff(u,2) + fred(K, u));
A.bc = 'periodic';
u = A\f;
pass(9) = norm(A*u-f) < tol;

%%
% From #2360.
K = @(x,y) 1+0*x;
N = chebop(@(x,u) - u + exp(x) + 1/3*x*(1-exp(3*x)) + x*volt(K, u^3), [0, 1]);
u = N\0;
pass(10) = norm(u-chebfun(@exp, [0 1])) < tol;

end


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Copyright 2015-2018 Jonas Koennemanm, Giovanni Licitra
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function [solution,times,problem] = vanderpol

  problem = ocl.Problem(10, @varsfun, @daefun, @pathcosts, 'N', 30);

  % intial state bounds
  problem.setInitialBounds('x',     0);
  problem.setInitialBounds('y',     1);

  % Get and set initial guess
  initialGuess = problem.getInitialGuess();
  initialGuess.states.x.set(-0.2);

  % Run solver to obtain solution
  [solution,times] = problem.solve(initialGuess);

  % plot solution
  figure
  hold on
  plot(times.states.value,solution.states.x.value,'-.','LineWidth',2)
  plot(times.states.value,solution.states.y.value,'--k','LineWidth',2)
  stairs(times.controls.value,solution.controls.F.value,'r','LineWidth',2)
  xlabel('time')
  legend({'x','y','u'})

  snapnow;
end

function varsfun(svh)
  % Scalar x:  -0.25 <= x <= inf
  % Scalar y: unbounded
  svh.addState('x', 'lb', -0.25, 'ub', inf);
  svh.addState('y');

  % Scalar u: -1 <= F <= 1
  svh.addControl('F', 'lb', -1, 'ub', 1);
end

function daefun(daeh,x,~,u,~)
  daeh.setODE('x', (1-x.y^2)*x.x - x.y + u.F);
  daeh.setODE('y', x.x);
end

function pathcosts(ch,x,~,u,~)
  ch.add( x.x^2 );
  ch.add( x.y^2 );
  ch.add( u.F^2 );
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
classdef CompressingVademecumData < handle

    properties (Access = private)
        qV
        rhoV
        xiV
        phiV
        cellDataSameSign
        cellDataDifferentSign
        phiIndeces
    end

    properties (Access = private)
        compressedFileName
        nMx
        nMy
        nPhiQuarter
        nPhiT
    end

    methods (Access = public)

        function obj = CompressingVademecumData()
            obj.init();
            obj.initVariables();
            obj.loadVademecumDataBase();
            obj.compute();
            obj.saveData();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.compressedFileName = 'OptimalSuperEllipses';
            obj.nMx = 20;
            obj.nMy = 20;
            obj.nPhiQuarter = 10;
            obj.nPhiT = obj.computeNphiT();
        end

        function nPT = computeNphiT(obj)
            n0  = obj.nPhiQuarter;
            n1  = 2*(n0 -1)+1;
            nPT = 2*(n1 -1)+1;
        end

        function initVariables(obj)
            nmx = obj.nMx;
            nmy = obj.nMy;
            np = obj.nPhiT;
            obj.qV   = zeros(nmx,nmy,np);
            obj.rhoV = zeros(nmx,nmy,np);
            obj.xiV  = zeros(nmx,nmy,np);
            obj.phiV = zeros(nmx,nmy,np);
        end

        function loadVademecumDataBase(obj)
            s.nMx = obj.nMx;
            s.nMy = obj.nMy;
            s.nPhi = obj.nPhiQuarter;
            v = VademecumDataLoader(s);
            obj.cellDataSameSign = v.cellDataSameSign;
            obj.cellDataDifferentSign = v.cellDataDifferentSign;
        end

        function compute(obj)
            for iMx = 1:obj.nMx
                for iMy = 1:obj.nMy
                    for iPhi = 1:nPhi
                        obj.computePhiIndices(iPhi)
                        obj.computeXiValues(iMx,iMy,iPhi);
                        obj.computeRhoValues(iMx,iMy,iPhi);
                        obj.computePhiValues(iMx,iMy,iPhi);
                        obj.computeQValues(iMx,iMy,iPhi);
                        obj.plotFigures(iMx,iMy,iPhi);
                    end
                end
            end
        end

        function computePhiIndices(obj,iPhi)
            iphi1 = iPhi;
            iphi2 = 2*nPhi - iPhi;
            iphi3 = 2*(nPhi -1) + iPhi;
            iphi4 = 2*(2*nPhi -1) - iPhi;
            obj.phiIndeces = [iphi1,iphi2,iphi3,iphi4];
        end

        function computeXiValues(obj,iMx,iMy,iPhi)
             [cS,cD] = obj.obtainCellData(obj,iMx,iMy,iPhi);
             xiQ = [cS.xi cS.xi cD.xi cD.xi];
             ind = obj.phiIndeces;
             obj.xiV(iMx,iMy,ind) = xiQ;
        end

        function computeRhoValues(obj,iMx,iMy,iPhi)
             [cS,cD] = obj.obtainCellData(obj,iMx,iMy,iPhi);
             rhoQ = [cS.rho cS.rho cD.rho cD.rho];
             ind = obj.phiIndeces;
             obj.rhoV(iMx,iMy,ind) = rhoQ;
        end

        function computePhiValues(obj,iMx,iMy,iPhi)
            [cS,cD] = obj.obtainCellData(obj,iMx,iMy,iPhi);
            qQ = [cS.phi pi/2-cS.phi cD.phi 3*pi/2-cD.phi];
            ind = obj.phiIndeces;
            obj.qV(iMx,iMy,ind) = qQ;
        end

        function computeQValues(obj,iMx,iMy,iPhi)
            [cS1,cD1] = obj.obtainCellData(obj,iMx,iMy,iPhi);
            [cS2,cD2] = obj.obtainCellData(obj,iMy,iMx,iPhi);
            qQ = [cS1.q cS2.q cD1.q cD2.q];
            ind = obj.phiIndeces;
            obj.qV(iMx,iMy,ind) = qQ;
        end

        function [cS,cD] = obtainCellData(obj,iMx,iMy,iPhi)
           cS = obj.cellDataSameSign{iMx}{iMx,iMy,iPhi};
           cD = obj.cellDataDifferentSign{iMx}{iMx,iMy,iPhi};
        end

        function plotFigures(obj,iMx,iMy,iPhi)
            [cS1,cD1] = obj.obtainCellData(obj,iMx,iMy,iPhi);
            [cS2,cD2] = obj.obtainCellData(obj,iMy,iMx,iPhi);
            obj.plotFigure(cS1,1);
            obj.plotFigure(cD1,2);
            obj.plotFigure(cS2,3);
            obj.plotFigure(cD2,4);
        end

       function saveData(obj)
            q = obj.qV;
            rho = obj.rhoV;
            xi  = obj.xiV;
            phi = obj.phiV;
            save(obj.compressedFileName,'q','rho','xi','phi')
        end

    end

    methods (Access = private, Static)

        function [xiP,rhoP,phiP,qP] = loadData(cP)
            xiP  = cP.xi;
            rhoP = cP.rho;
            phiP = cP.phi;
            qP   = cP.q;
        end

        function plotFigure(cellData,nFigure)
            figure(nFigure);
            clf
            cellData.mesh.plot;
            drawnow
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function imResult = blendMode_ColorBurn(A, B, offsetW, offsetH)
%% Color Burn blending mode: divides the inverted bottom layer by the top
%   layer, and then inverts the result. This darkens the top layer
%   increasing the contrast to reflect the color of the bottom layer.
%   The darker the bottom layer, the more its color is used. Blending with
%   white produces no difference.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_ColorBurn));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

ind = B == 0;
indCompl = abs(ind - 1);
tmp = max(0, 1 - (1 - A) ./ B);
C = ind .* B + indCompl .* tmp;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end


text:
function [err, Qd, dfx, dmat2] = QRDecom(dmat, cmat)
%
%   [err,] = QRDecom.m ()
%
%Purpose:
%
%
%
%Input Parameters:
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Gang Chen
%     Date : Tue Mar 23 16:09:23 EST 2004
%     SSCC/NIMH/ National Institutes of Health, Bethesda MD 20892

%Define the function name for easy referencing
FuncName = 'QRDecom.m';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

% Tolerance for computing rank from diag(R) after QR decomposition
%[nrows,ncols] = size(dmat);

% Find the null space of the constraints matrix
[Qc,Rc,Ec] = qr(cmat');
pc = Rrank(Rc);
Qc0 = Qc(:,pc+1:end);

% Do qr decomposition on design matrix projected to null space
Dproj = dmat*Qc0;
[Qd,Rd,Ed] = qr(Dproj,0);
dfx = Rrank(Rd);
Qd = Qd(:,1:dfx);
%Rd = Rd(1:dfx,1:dfx);

% Return reduced design matrix if requested
if nargout>3
   dmat2 = Qd' * dmat;
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function workbook=loadWorkbook(fileName,createEmpty)
% loadWorkbook
%   Loads an Excel file into a Workbook object using the Java library Apache POI
%
%   fileName    name of the Excel file. If it doesn't exist it will be
%               created
%   createEmpty true if an empty workbook should be created if the file
%               didn't exist (opt, default false)
%
%   workbook    Workbook object representing the Excel file
%
%   Usage: workbook=loadWorkbook(fileName,createEmpty)

if nargin<2
    createEmpty=false;
end

%Check if the user has MATLAB Text Analytics Toolbox installed, as it comes
%with its own conflicting version of the required Apache POI files
if exist('vaderSentimentScores.m')== 2
    error(['MATLAB Text Analytics Toolbox found. This should be uninstalled ' ...
           'if you want to read/write Excel files. See RAVEN GitHub Issues '...
           'page for instructions.'])
end

%Adds the required classes to the static Java path if not already added
addJavaPaths();

%Import required classes from Apache POI
import org.apache.poi.ss.usermodel.*;
import org.apache.poi.ss.util.*;
import java.io.FileInputStream;
import org.apache.poi.hssf.usermodel.*;
import org.apache.poi.xssf.usermodel.*;

%Check if the file exists
if ~isfile(fileName)
    if createEmpty==false
        EM='The Excel file could not be found';
        dispEM(EM);
    else
        %Create an empty workbook
        [~,~,I]=fileparts(fileName);
        if strcmpi(I,'.xls')
            workbook=HSSFWorkbook();
        else
            if strcmpi(I,'.xlsx')
                workbook=XSSFWorkbook();
            else
                EM='The file name must end in .xls or .xlsx';
                dispEM(EM);
            end
        end
    end
else
    %Opens the workbook. The input stream is needed since it will otherwise
    %keep the file open
    is=FileInputStream(getFullPath(fileName));
    workbook=WorkbookFactory.create(is);
    is.close();
end
end

author: SysBioChalmers

repo: RAVEN


text:
classdef PropertySetter < Algorithm

    properties (Access = public)
        node;
        property;
    end

    methods (Access = public)
        function obj = PropertySetter()
            obj.name = 'propertySetter';
            obj.inputPort = DataType.kAny;
            obj.outputPort = DataType.kNull;
        end

        function computedSignal = compute(obj,data)
            obj.node.(obj.property) = data;
            computedSignal = [];
        end

        function str = toString(obj)
            str = sprintf('%s_%s_%s',obj.name,obj.node.toString(),obj.property);
        end

        function editableProperties = getEditableProperties(obj)
            editableProperties = Property('property',obj.property);
        end

        function metrics = computeMetrics(~,input)
            metrics = Metric(1,1,0);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function stroud_test207 ( )

%*****************************************************************************80
%
%% STROUD_TEST207 tests the Stroud EN_R2 rules on monomials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 January 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'STROUD_TEST207\n' );
  fprintf ( 1, '  Demonstrate the use of Stroud rules for the region\n' );
  fprintf ( 1, '  EN_R2, that is, all of N-dimensional space, with the\n' );
  fprintf ( 1, '  weight function W(X) = exp ( - X1^2 - X2^2 ... -XN^2 )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use the formulas to integrate various monomials of\n' );
  fprintf ( 1, '  the form X1^ALPHA1 * X2^ALPHA2 * ... XN^ALPHAN\n' );
  fprintf ( 1, '  and compare to the exact integral.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The precision of each formula is known, and we only use\n' );
  fprintf ( 1, '  a formula if its precision indicates it should be able to\n' );
  fprintf ( 1, '  produce an exact result.\n' );

  for n = 1 : 7

    alpha = zeros ( n, 1 );
    en_r2_test ( n, alpha );

    alpha = zeros ( n, 1 );
    alpha(1) = 2;
    en_r2_test ( n, alpha );

    alpha = zeros ( n, 1 );
    alpha(2) = 4;
    en_r2_test ( n, alpha );

    alpha = zeros ( n, 1 );
    i = mod ( 3 - 1, n ) + 1;
    alpha(i) = 6;
    en_r2_test ( n, alpha );

    alpha = zeros ( n, 1 );
    alpha(1) = 2;
    alpha(2) = 4;
    en_r2_test ( n, alpha );

    alpha = zeros ( n, 1 );
    i = mod ( 4 - 1, n ) + 1;
    alpha(i) = 8;
    en_r2_test ( n, alpha );

    alpha = zeros ( n, 1 );
    i = mod ( 5 - 1, n ) + 1;
    alpha(i) = 10;
    en_r2_test ( n, alpha );

    alpha = 1 : n;
    en_r2_test ( n, alpha );

    alpha(1:n) = 2;
    en_r2_test ( n, alpha );

  end

  return
end


text:
function test_ft_write_mri

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_read_mri ft_write_mri

mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));
mri = ft_convert_units(mri, 'mm');

%%

dataformat = {
  'analyze'
  'analyze_hdr'
  'analyze_img'
  'analyze_old'
  'freesurfer_mgz'
  'mgh'
  'mgz'
  'nifti'
  'nifti2'
  'nifti_gz'
  'nifti_img'
  'nifti_spm'
  'vmp'
  'vmr'
  };

%%

p = tempname;

for i=1:numel(dataformat)
  filename = fullfile(p, ['Subject01_' dataformat{i}]);
  % pass it as separate fields
  ft_write_mri(filename, mri.anatomy, 'dataformat', dataformat{i}, 'transform', mri.transform, 'unit', mri.unit, 'coordsys', mri.coordsys);
end

% remove the temporary directory with all files inside it
rmdir(p, 's');

%%

p = tempname;

for i=1:numel(dataformat)
  filename = fullfile(p, ['Subject01_' dataformat{i}]);
  % pass it as a structure
  ft_write_mri(filename, mri, 'dataformat', dataformat{i});
end

% remove the temporary directory with all files inside it
rmdir(p, 's');

%%

p = tempname;

% scale between -1 and 1, we will skip the unsigned formats
dat = double(mri.anatomy);
dat = dat./max(abs(dat(:)));

mri.anatomy = double(dat);                         filename = fullfile(p, 'Subject01_double'); ft_write_mri(filename, mri, 'dataformat', 'analyze_old');
mri.anatomy = single(dat);                         filename = fullfile(p, 'Subject01_single'); ft_write_mri(filename, mri, 'dataformat', 'analyze_old');
mri.anatomy = int32(dat.*double(intmax('int32'))); filename = fullfile(p, 'Subject01_int32');  ft_write_mri(filename, mri, 'dataformat', 'analyze_old');
mri.anatomy = int16(dat.*double(intmax('int16'))); filename = fullfile(p, 'Subject01_int16');  ft_write_mri(filename, mri, 'dataformat', 'analyze_old');

mri.anatomy = double(dat);                         filename = fullfile(p, 'Subject01_double'); ft_write_mri(filename, mri, 'dataformat', 'nifti');
mri.anatomy = single(dat);                         filename = fullfile(p, 'Subject01_single'); ft_write_mri(filename, mri, 'dataformat', 'nifti');
mri.anatomy = int32(dat.*double(intmax('int32'))); filename = fullfile(p, 'Subject01_int32');  ft_write_mri(filename, mri, 'dataformat', 'nifti');
mri.anatomy = int16(dat.*double(intmax('int16'))); filename = fullfile(p, 'Subject01_int16');  ft_write_mri(filename, mri, 'dataformat', 'nifti');
mri.anatomy = uint8(dat.*double(intmax('int8')));  filename = fullfile(p, 'Subject01_int8');   ft_write_mri(filename, mri, 'dataformat', 'nifti');

% remove the temporary directory with all files inside it
rmdir(p, 's');


text:
function vector = huff2norm(zipped,info)
%HUFF2NORM   Huffman codification (decoder)
%   For vectors, HUFF2NORM(X,INFO) returns a decoded vector from a Huffman coded version X
%   with code words given by INFO.
%   For matrices, X(:) is used as input.
%
%   Input X must be of uint8 type, while the output is a uint8 array.
%
%   For more details about how is INFO structured, see NORM2HUFF help.
%
%   See also NORM2HUFF

%   $Author: Giuseppe Ridino' $
%   $Revision: 1.0 $  $Date: 10-May-2004 15:03:04 $

% ensure to handle uint8 input vector
if ~isa(zipped,'uint8'),
	error('input argument must be a uint8 vector')
end

% create the 01 sequence
len = length(zipped);
string = repmat(uint8(0),1,len.*8);
bitindex = 1:8;
for index = 1:len,
	string(bitindex+8.*(index-1)) = uint8(bitget(zipped(index),bitindex));
end

% adjust string
string = logical(string(:)'); % make a row of it
len = length(string);
string((len-info.pad+1):end) = []; % remove 0 padding
len = length(string);

% build output
weights = 2.^(0:51);
vector = repmat(uint8(0),1,info.length);
vectorindex = 1;
codeindex = 1;
code = 0;
for index = 1:len,
	code = bitset(code,codeindex,string(index));
	codeindex = codeindex+1;
	byte = decode(bitset(code,codeindex),info);
	if byte>0, % a code has been found
		vector(vectorindex) = byte-1;
		codeindex = 1;
		code = 0;
		vectorindex = vectorindex+1;
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function byte = decode(code,info)
byte = info.huffcodes(code);


text:
function [success] = downloadGithubRepositories(varargin)
	% Downloads Github repositories repositories.
	% Biafra Ahanonu
	% started: 2019.01.14 [10:23:05]
	% changelog
		% 2020.04.03 [11:15:32] - Allow inputs to use getOptions. Also allow force updating of git repository.
		% 2020.04.03 [14:02:33] - Save downloaded compressed files (e.g. zips) to a sub-folder.
		% 2020.06.28 [13:08:16] - Final implementation of force update, to bring to most current version of all git directories.
		% 2021.01.22 [13:18:25] - Update to allow regexp backup to find name of downloaded Github repo folder after unzipping, e.g. in cases where a release or non-master branch is downloaded. - IGNORE
		% 2021.02.01 [15:19:40] - Update `_external_programs` to call ciapkg.getDirExternalPrograms() to standardize call across all functions.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.01.27 [08:49:09] - If multiple repos requested but one already downloaded, previously exited before downloading the other repos.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	%========================
	% 1 = force update of the git repository, 0 = skip if already downloaded
	options.forceUpdate = 0;
	% Str: directory of external download path.
	options.signalExtractionDir = ciapkg.getDirExternalPrograms();
	options.gitNameDisp = {'NoRMCorre'};
	options.gitRepos = {'https://github.com/flatironinstitute/NoRMCorre/archive/master.zip'};
	options.outputDir = {'normcorre'};
	options.gitName = {'NoRMCorre-master'};

	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		signalExtractionDir = options.signalExtractionDir;

		gitNameDisp = options.gitNameDisp;
		gitRepos = options.gitRepos;
		outputDir = options.outputDir;
		gitName = options.gitName;
		nRepos = length(outputDir);

		% If forcing an update, make sure to remove all external programs from the path
		if options.forceUpdate==1
			ciapkg.download.updatePathCleanup(options.signalExtractionDir);
			pause(0.02)
		end

		for gitNo = 1:nRepos
			display(repmat('=',1,7))
			fprintf('%s\n',gitNameDisp{gitNo});

			outDirPath = [signalExtractionDir filesep outputDir{gitNo}];

			oldDir = [signalExtractionDir filesep gitName{gitNo}];
			newDir = [signalExtractionDir filesep outputDir{gitNo}];

			% If the new directory exist
			if exist(newDir,'dir')&options.forceUpdate==0
				fprintf('Package already exists (and no forced update), skipping: %s.\n',newDir)
				success = 1;
				continue;
			end

			% Make download directory
			rawSavePathDownload = [signalExtractionDir filesep '_downloads'];
			if ~exist(rawSavePathDownload,'dir');mkdir(rawSavePathDownload);fprintf('Made folder: %s',rawSavePathDownload);end
			% Download save path
			rawSavePathDownload = [rawSavePathDownload filesep outputDir{gitNo} '.zip'];

			% Remove directory and downloaded zip file if forcing an update.
			if options.forceUpdate==1
				if exist(outDirPath,'dir')
					fprintf('Removing directory to allow updates: %s\n', outDirPath);
					status = 0;
					tryCount = 1;
					% Try to delete several times else print out warning
					while status==0&tryCount<4
						[status,msg,msgID] = rmdir(outDirPath,'s');
						tryCount = tryCount + 1;
						pause(0.1);
					end
					if status==0
						warning('	Something went wrong removing the directory, check for proper permissions, etc.')
						warning(['	' msg])
						disp(['	' msgID])
					else
						disp('	Successfully removed directory!')
					end
				end
				if exist(rawSavePathDownload,'file')==2
					fprintf('Removing zip file to allow update: %s\n',rawSavePathDownload)
					delete(rawSavePathDownload)
				end
			else

			end

			% Download git repo zip
			if exist(rawSavePathDownload,'file')~=2
				fprintf('Downloading %s file to %s\n',gitRepos{gitNo},rawSavePathDownload)
				websave(rawSavePathDownload,gitRepos{gitNo});
			else
				fprintf('Already downloaded %s\n',rawSavePathDownload)
			end

			% To ensure do not have mixed extraction directories, still check directory is there before extracting even in case of forced update.
			if exist(outDirPath,'dir')
				fprintf('Already extracted %s\n', outDirPath);
				continue;
			end

			% Unzip the repo file
			fprintf('Unzipping file %s\n',rawSavePathDownload)
			filenames = unzip(rawSavePathDownload,signalExtractionDir);
			% cellfun(@disp,filenames,'UniformOutput',false)

			% Rename to proper folder for calciumImagingAnalysis
			fprintf('Renaming %s to %s \n',[signalExtractionDir filesep gitName{gitNo}],[signalExtractionDir filesep outputDir{gitNo}])
			if strcmp(oldDir,newDir)==1
				disp('Same directory, ignore name change!')
			else
				movefile(oldDir,newDir);
			end

			% fprintf('\n\n')
		end

		success = 1;
	catch err
		success = 0;
		disp('Check internet connection or download and unzip files manually, see:')
		try
			cellfun(@disp,gitRepos,'UniformOutput',false)
		catch err2
			display(repmat('@',1,7))
			disp(getReport(err2,'extended','hyperlinks','on'));
			display(repmat('@',1,7))
		end
		fprintf('\n\n')
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function tf=mocov_is_absolute_path(fn)
% indicate whether a path is absolute
%
% tf=mocov_is_absolute_path(fn)
%
% Input:
%   fn          filename or path
%
% Output:
%   tf          true if fn is an absolute path, false otherwise
%

    n=numel(fn);
    if ispc()
        tf=n>=2 && fn(2)==':';
    else
        tf=n>=1 && (fn(1)=='/' || fn(1)=='~');
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function disp(obj)
% Disp a NIFTI-1 object
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: disp.m 1143 2008-02-07 19:33:33Z spm $

sz = size(obj);
fprintf('NIFTI object: ');
if length(sz)>4,
    fprintf('%d-D\n',length(sz));
else
    for i=1:(length(sz)-1),
        fprintf('%d-by-',sz(i));
    end;
    fprintf('%d\n',sz(end));
end;
if prod(sz)==1,
    display(structn(obj))
end;
return;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
% Author: James A. Douthwaite

classdef agent_2D_formation_VO < agent_2D_VO & agent_formation
%% INITIALISE THE AGENT SPECIFIC PARAMETERS
    properties
        % PROPERTITIES UNIQUE TO AGENT CLASS
    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function [this] = agent_2D_formation_VO(varargin)
            % CALL THE SUPERCLASS CONSTRUCTOR
            this = this@agent_2D_VO(varargin);

            this.maxSamples = 1;

            % //////////////////// SENSOR PARAMETERS //////////////////////
%             [obj.SENSORS] = obj.GetDefaultSensorParameters();       % Default sensing
            [this.SENSORS] = this.GetCustomSensorParameters();       % Experimental sensing
            % /////////////////////////////////////////////////////////////

            % //////////////// Check for user overrides ///////////////////
            [this] = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Main
        function [this] = main(this,ENV,varargin)
            % INPUTS:
            % obj      - The agent object
            % TIME     - The current time structure
            % varargin - Cell array of inputs
            % OUTPUTS:
            % obj      - The updated object

            % PLOT AGENT FIGURE
            visualiseProblem = 0;
            visualiseAgent = 1;
            if this.objectID == visualiseAgent && visualiseProblem == 1
                overHandle = figure(100 + this.objectID);
                hold on; grid on;
                xlabel('x_{m}'); ylabel('y_{m}'); zlabel('z_{m}');
            end

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            % UPDATE THE AGENT WITH THE NEW ENVIRONMENTAL INFORMATION
            [this,obstacleSet,agentSet,~] = this.GetAgentUpdate(ENV,varargin{1});

            % ////////////////// FORMATION CONTROLLER /////////////////////
            % We wish to conduct formation control with the other agents
            % in the simulation only.
            L = 0;
            if ~isempty(agentSet)
                % PASS AGENT SET TO FORMATION CONTROLLER
                [heading,speed,L] = this.formationControl_distance(agentSet);          % Get the force vector

                desiredVelocity = heading*speed;
            end

            % ////////////////// OBSTACLE AVOIDANCE ///////////////////////
            % Modify the desired velocity with the augmented avoidance velocity.
%             avoidanceSet = [obstacleSet;agentSet];
            algorithm_start = tic; algorithm_indicator = 0;
%             if ~isempty(avoidanceSet) && avoidanceEnabled
%                 algorithm_indicator = 1;
%                 % GET THE UPDATED DESIRED VELOCITY
%                 [desiredHeadingVector,desiredSpeed] = this.GetAvoidanceCorrection(...
%                     dt,...
%                     desiredVelocity,...
%                     visualiseProblem);
%                 desiredVelocity = desiredHeadingVector*desiredSpeed;
%             end
            algorithm_dt = toc(algorithm_start);

            % //////////////////// AGENT CONTROL //////////////////////////
            [this] = this.Controller(ENV.dt,desiredVelocity);

            % ////////////// RECORD THE AGENT-SIDE DATA ///////////////////
            this = this.writeAgentData(ENV,algorithm_indicator,algorithm_dt);
            this.DATA.inputNames = {'$v_x$ (m/s)','$v_y$ (m/s)','$\dot{\psi}$ (rad/s)'};
            this.DATA.inputs(1:length(this.DATA.inputNames),ENV.currentStep) = this.localState(4:6);         % Record the control inputs
            this.DATA.lypanov(ENV.currentStep) = L;                        % Record the lypanov value
        end
    end
end

repo: OpenMAS


text:
function Fmenu = spm_Menu(action, varargin)
% SPM Menu window
% FORMAT Fmenu = spm_Menu('Create',Vis)
% Create the SPM Menu window (tag property set to 'Menu')
% Vis    - visibility: {'on'} or 'off'
% Fmenu  - handle of figure created
%
% FORMAT Fmenu = spm_Menu('Switch',Modality)
% Switch the SPM Menu window to the specified modality
%
% FORMAT spm_Menu('Close')
% Close the SPM Menu window
%__________________________________________________________________________
% Copyright (C) 2018-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_Menu.m 7755 2019-12-16 13:19:28Z spm $

if nargin < 1, action = 'Create'; end

switch lower(action)
    case 'create'
        spm_Menu('Close');
        Fmenu = create_menu(varargin{1:end});

    case 'switch'
        Fmenu = spm_figure('FindWin','Menu');
        [Modality,ModNum] = spm('CheckModality',varargin{1});
        if ~isempty(Fmenu)
            if strcmpi(Modality,'PET')
                set(findobj(Fmenu, '-regexp', 'Tag', 'FMRI'), 'Visible', 'off');
                set(findobj(Fmenu, '-regexp', 'Tag', 'EEG'),  'Visible', 'off');
                set(findobj(Fmenu, '-regexp', 'Tag', 'PET'),  'Visible', 'on' );
            elseif strcmpi(Modality,'FMRI')
                set(findobj(Fmenu, '-regexp', 'Tag', 'EEG'),  'Visible', 'off');
                set(findobj(Fmenu, '-regexp', 'Tag', 'PET'),  'Visible', 'off');
                set(findobj(Fmenu, '-regexp', 'Tag', 'FMRI'), 'Visible', 'on' );
            else
                set(findobj(Fmenu, '-regexp', 'Tag', 'PET'),  'Visible', 'off');
                set(findobj(Fmenu, '-regexp', 'Tag', 'FMRI'), 'Visible', 'off');
                set(findobj(Fmenu, '-regexp', 'Tag', 'EEG'),  'Visible', 'on' );
            end
            set(findobj(Fmenu,'Tag','Modality'),'Value',ModNum,'UserData',ModNum);
        else
            warning('SPM Menu window not found.');
        end

    case 'close'
        Fmenu = spm_figure('FindWin','Menu');
        delete(Fmenu);

    otherwise
        error('Unknown action ''%s''.',action);
end

%==========================================================================
%-function Fmenu = create_menu(Vis)
%==========================================================================
function Fmenu = create_menu(Vis)

if ~nargin, Vis = 'on'; end

%-Get size and scalings
%--------------------------------------------------------------------------
PF    = spm_platform('fonts');             %-Font names (for this platform)
FS    = spm('FontSizes');                  %-Scaled font sizes

WS    = spm('WinScale');                   %-Window scaling factors
Rect0 = spm('WinSize','0',1);
RectM = spm('WinSize','Menu',1).*WS;       %-Menu window rectangle: 400x445
Pos   = [Rect0(1) Rect0(2) 0 0] + RectM;

%-Create Menu window
%--------------------------------------------------------------------------
if ~isempty(which('openfig'))
    Fmenu = openfig(fullfile(spm('Dir'),'spm_Menu.fig'),'new','invisible');
else
    fprintf('\nFunction openfig.m is missing. Download it from:\n');
    fprintf(' https://savannah.gnu.org/bugs/download.php?file_id=33499\n');
    Fmenu = [];
end

set(Fmenu,'Name',[spm('Version') ': Menu']);
set(Fmenu,'Units','pixels', 'Position',Pos);
set(Fmenu,'Color',[1 1 1]*.8);

% Fmenu = figure('Tag','Menu',...
%     'Name',[spm('Version') ': Menu'],...
%     'IntegerHandle','off',...
%     'NumberTitle','off',...
%     'Units','pixels',...
%     'Position',Pos,...
%     'Resize','on',...
%     'Pointer','arrow',...
%     'Color',[1 1 1]*.8,...
%     'MenuBar','none',...
%     'DockControls','off',...
%     'CloseRequestFcn','spm(''Quit'')',...
%     'HandleVisibility','off',...
%     'DefaultUicontrolFontName',PF.helvetica,...
%     'Visible','off');

%-Set SPM colour
%--------------------------------------------------------------------------
set(findobj(Fmenu,'Tag', 'frame'),'BackgroundColor',spm('colour'));
if ispc && strcmpi(spm_check_version,'matlab')
    try
        %M   = getframe(Fmenu);
        %col = double(M.cdata(floor(size(M.cdata,1)/3),floor(size(M.cdata,2)/2),:));
        col = [204 204 204];
        set(findobj(Fmenu,'UserData','LABEL'),'BackgroundColor',col/255);
    end
elseif ismac
    set(findobj(Fmenu,'UserData','LABEL'),'Visible','off','Tag','');
end

%-Set Utils
%--------------------------------------------------------------------------
set(findobj(Fmenu,'Tag', 'Utils'), 'String',{'Utils...',...
    'CD',...
    'PWD',...
    'Run M-file',...
    'Load MAT-file',...
    'Save MAT-file',...
    'Delete files',...
    'Show SPM',...
    'Show MATLAB'});
set(findobj(Fmenu,'Tag', 'Utils'), 'UserData',{...
    ['spm(''FnBanner'',''CD'');' ...
     'cd(spm_select(1,''dir'',''Select new working directory''));' ...
     'spm(''alert"'',{''New working directory:'',[''    '',pwd]},''CD'',1);'],...
    ['spm(''FnBanner'',''PWD'');' ...
     'spm(''alert"'',{''Present working directory:'',[''    '',pwd]},''PWD'',1);'],...
    ['spm(''FnBanner'',''Run M-file'');' ...
     'spm(''Run'');' ...
     'fprintf(''%-40s: %30s\n'',''Completed'',spm(''time''));'],...
    ['spm(''FnBanner'',''Load MAT-file'');' ...
     'load(spm_select(1,''mat'',''Select MAT-file''));'],...
    ['spm(''FnBanner'',''Save MAT-file'');' ...
     'save(spm_input(''Output filename'',1,''s''), spm_get_defaults(''mat.format''));'],...
    ['spm(''FnBanner'',''Delete files'');' ...
     'spm(''Delete'');'],...
    ['spm(''FnBanner'',''Show SPM'');' ...
     'spm(''Show'');'],...
    ['spm(''FnBanner'',''Show Command Window'');' ...
     'commandwindow;']});

%-Set Toolboxes
%--------------------------------------------------------------------------
xTB       = spm('tbs');
if ~isempty(xTB)
    set(findobj(Fmenu,'Tag', 'Toolbox'),'String',{'Toolbox:' xTB.name });
    set(findobj(Fmenu,'Tag', 'Toolbox'),'UserData',xTB);
else
    set(findobj(Fmenu,'Tag', 'Toolbox'),'Enable','off');
end
set(Fmenu,'Visible',Vis);


text:
classdef cDPEA < ALGORITHM
% <multi> <real/integer/label/binary/permutation> <constrained>
% Constrained dual-population evolutionary algorithm

%------------------------------- Reference --------------------------------
% M. Ming, A. Trivedi, R. Wang, and D. Srinivasan, A dual-population based
% evolutionary algorithm for constrained multi-objective optimization, IEEE
% Transactions on Evolutionary Computation, 2021, 25(4): 739-753.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Mengjun Ming

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population1 = Problem.Initialization();
            Population2 = Population1;
            alpha       = 2./(1+exp(1).^(-Problem.FE*10/Problem.maxFE))-1;
            para        = ceil(Problem.maxFE/Problem.N)/2 - ceil(Problem.FE/Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population1)
                % Re-rank
                Population1 = Population1(randperm(Problem.N));
                Population2 = Population2(randperm(Problem.N));
                Lia   = ismember(Population2.objs,Population1.objs, 'rows');
                gamma = 1-sum(Lia)/Problem.N;
                [Population1,FrontNo1] = EnvironmentalSelection(Population1,Problem.N,alpha);
                [Population2,FrontNo2] = EnvironmentalSelection_noCon(Population2,Problem.N,alpha,gamma,para);

                % Offspring Reproduction
                Population_all   = [Population1,Population2];
                RankSolution_all = [FrontNo1,FrontNo2];
                MatingPool = TournamentSelection(2,2*Problem.N,RankSolution_all);
                Offspring  = OperatorGAhalf(Problem,Population_all(MatingPool));

                % Environmental Selection
                alpha = 2./(1+exp(1).^(-Problem.FE*10/Problem.maxFE)) - 1;
                para  = ceil(Problem.maxFE/Problem.N)/2 - ceil(Problem.FE/Problem.N);
                [Population1,~] = EnvironmentalSelection([Population1,Offspring],Problem.N,alpha);
                [Population2,~] = EnvironmentalSelection_noCon([Population2,Offspring],Problem.N,alpha,gamma,para);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [Xnorm] = L12norm(X)
% ||X||_{1,2} = sum_i||X^i||_2
Xnorm = sum(sqrt(sum(X.^2,2)));

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR