text:
function score = IGDX(Population,POS)
% <min> <multi/many> <real/integer/label/binary/permutation> <large/none> <constrained/none> <expensive/none> <sparse/none> <multimodal> <dynamic/none>
% Inverted generational distance in the decision space

%------------------------------- Reference --------------------------------
% A. Zhou, Q. Zhang, and Y. Jin, Approximating the set of Pareto-optimal
% solutions in both the decision and objective spaces by an estimation of
% distribution algorithm, IEEE Transactions on Evolutionary Computation,
% 2009, 13(5): 1167-1189.
%------------------------------- 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".
%--------------------------------------------------------------------------

    PopDec = Population.decs;
    if size(PopDec,2) ~= size(POS,2)
        score = nan;
    else
        score = mean(min(pdist2(POS,PopDec),[],2));
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = animation(varargin)
% ANIMATION M-file for animation.fig
%      ANIMATION, by itself, creates a new ANIMATION or raises the existing
%      singleton*.
%
%      H = ANIMATION returns the handle to a new ANIMATION or the handle to
%      the existing singleton*.
%
%      ANIMATION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ANIMATION.M with the given input arguments.
%
%      ANIMATION('Property','Value',...) creates a new ANIMATION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before animation_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to animation_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help animation
% Last Modified by GUIDE v2.5 09-Jul-2013 19:10:43
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @animation_OpeningFcn, ...
                   'gui_OutputFcn',  @animation_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before animation is made visible.
function animation_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to animation (see VARARGIN)

    handles.portfolio = varargin{1};
    handles.windowAnimation = varargin{2};
    handles.status = 1;
    [r c] = size(handles.portfolio);
    handles.nframes = r;
    handles.c = c;
    guidata(hObject, handles);

    startAnimating(hObject, eventdata, handles, varargin);

    % Choose default command line output for temp
    handles.output = hObject;

    % Update handles structure
    guidata(hObject, handles);
% UIWAIT makes animation wait for user response (see UIRESUME)
% uiwait(handles.animationGui);

function startAnimating(hObject, eventdata, handles, varargin)
    c = handles.c;
    handles.status = 1;

    try
        for i = 1:handles.nframes
            start = max(1, i-handles.windowAnimation+1);
            finish = i;
            toPlot = mean(handles.portfolio(start:finish,:));
            bar(handles.animationGraph,toPlot);
            axis([0.5 (c+0.5) 0 1]);
            grid on;
            if ishandle(handles.animationGui)
                M(i) = getframe;
                set(handles.time,'String',strcat('Time Period: ', int2str(i)));
            else
                return;
            end
            if handles.status == 0
                return;
            end
        end
    catch exception
    end

    handles.output = hObject;

% --- Outputs from this function are returned to the command line.
function varargout = animation_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in closeAnimation.
function closeAnimation_Callback(hObject, eventdata, handles)
% hObject    handle to closeAnimation (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    handles.status = 0;
    guidata(hObject, handles);
    close(handles.animationGui);

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function ex4bvp
%EX4BVP  Example 4 of the BVP tutorial.
%   This is a nerve impulse model considered in Example 7.1 of
%   R. Seydel, From Equilibrium to Chaos, Elsevier, 1988.  The
%   differential equations
%
%      y1' = 3*(y1 + y2 - 1/3*y1^3 + lambda)
%      y2' = -(y1 - 0.7 + 0.8*y2)/3
%
%   are to be solved subject to periodic boundary conditions
%
%      y1(0) = y1(T)
%      y2(0) = y2(T)
%
%   This is an example of non-separated boundary conditions, meaning
%   that some conditions involve values of the solution at both ends
%   of the interval.  Periodic boundary conditions are a common source
%   of such problems.  BVP4C is unusual in that it accepts problems
%   with non-separated boundary conditions.
%
%   The parameter lambda has the value -1.3 in the text cited.  The
%   period T is unknown, which is to say that the length of the
%   interval is not known. Such problems require some preparation.
%   By scaling the independent variable to tau = t/T the problem is
%   posed on the fixed interval [0,1]. When this is done, T becomes
%   an unknown parameter, but BVP4C provides for unknown parameters.

% Copyright 2002, The MathWorks, Inc.

%  Periodic functions evaluated in EX4INIT are used as a guess for the
%  solution on a crude mesh of 5 equally spaced points.  A guess of 2*pi
%  is provided for the unknown parameter T.
solinit = bvpinit(linspace(0,1,5),@ex4init,2*pi);

options = bvpset('stats','on');

sol = bvp4c(@ex4ode,@ex4bc,solinit,options);
T = sol.parameters;

%  The independent variable needs to be rescaled to its original value
%  before plotting the solution.  The initial guess is also plotted.
figure
plot(T*sol.x,sol.y(1,:),T*solinit.x,solinit.y(1,:),'o')
legend('solution found','initial guess');
axis([0 T -2.2 2]);
title('Nerve impulse model');
ylabel('solution y_1(t)');
xlabel(['Periodic solution with period ',num2str(sol.parameters,4)]);

% --------------------------------------------------------------------------

function v = ex4init(x)
%EX4INIT  Guess function for Example 4 of the BVP tutorial.
%   V = EX4INIT(X) returns a column vector V that is a guess for y(x).
v = [ sin(2*pi*x)
      cos(2*pi*x)];

% --------------------------------------------------------------------------

function dydt = ex4ode(t,y,T);
%EX4ODE  ODE funcion for Example 4 of the BVP tutorial.
dydt = [ 3*T*(y(1) + y(2) - 1/3*(y(1)^3) - 1.3);
         (-1/3)*T*(y(1) - 0.7 + 0.8*y(2)) ];

% --------------------------------------------------------------------------

function res = ex4bc(ya,yb,T)
%EX4BC  Boundary conditions for Example 4 of the BVP tutorial.
res = [ya(1) - yb(1)
       ya(2) - yb(2)
       T*(-1/3)*(ya(1) - 0.7 + 0.8*ya(2)) - 1];

author: Sable

repo: mcbench-benchmarks


text:
function value = f20_s ( dim_num, x )

%*****************************************************************************80
%
%% F20_S evaluates a function of a vector used in defining P2(Q).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real X(DIM_NUM), the value of the argument.
%
%    Output, real VALUE, the value of F20_S(X).
%
  value = 1.0;
  for i = 1 : dim_num
    value = value * ( 1.0 + ( f2 ( x(i) ) - 1.0 ) );
  end

  value = value - 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

for i = 1:size(X, 2)
	mu(:,i) = mean(X_norm(:,i));
	sigma(:,i) = std(X_norm(:,i));
	X_norm(:,i) = (X_norm(:,i) - mu(:,i)) / sigma(:,i);
end

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

end

repo: MachineLearning


text:
function [c, d] = lpdec1(x, h, g, extmod)
% LPDEC1   One-level Laplacian pyramid decomposition
%
%	[c, d] = lpdec1(x, h, g)
%
% Input:
%   x:      input signal
%   h, g:   two biorthogonal 1-D lowpass filters
%   extmod: [optional] extension mode (default is 'per')
%
% Output:
%   c:      coarse signal at half size
%   d:      detail signal at full size
%
% See also:	LPREC1

if ~exist('extmod', 'var')
    extmod = 'per';
end

nd = ndims(x);

% Computer the coarse signal by filter and downsample
c = x;
for dim = 1:nd
    c = filtdn(c, h, dim, extmod, 0);
end

% Compute the detail signal by upsample, filter, and subtract
% Even size filter needs to be adjusted to obtain perfect reconstruction
adjust = mod(length(g) + 1, 2);

p = c;
for dim = 1:nd
    p = upfilt(p, g, dim, extmod, adjust);
end

d = x - p;


text:
  function atten = xray_filters(mtype, thickness, energy, varargin)
%|function atten = xray_filters(mtype, thickness, energy, [options])
%|
%| Compute X-ray photon survival probability as a function of energy
%| for various materials.
%| in
%|	mtype			'aluminum', 'copper', ...
%|				can be a cell array {L} for multiple filters
%|	thickness		in cm (can be an array [L] if mtype is cell)
%|	energy	[N,1]		in keV (vector)
%| option
%|	'units'	cm | mm		default: cm
%| out
%|	atten	[N,L]		unitless survival probabilities
%|
%| Copyright 2001-04-27, Jeff Fessler, University of Michigan

if nargin == 1 && streq(mtype, 'test'), xray_filters_test, return, end
if nargin < 3, ir_usage, end

arg.units = 'cm';
arg = vararg_pair(arg, varargin);

if iscell(mtype)
	LL = length(mtype);
	if length(thickness) ~= LL, error 'thickness / material mismatch', end
	atten = zeros(length(energy), LL);
	for ll=1:LL
		atten(:,ll) = xray_filters(mtype{ll}, thickness(ll), ...
			energy(:), 'units', arg.units);
	end
return
end

mass_atten = xray_read_atten(mtype, energy, 'units', arg.units); % [N,L]
density = xray_read_dens(mtype, 'units', arg.units); % [L,1]
atten = exp(-mass_atten .* thickness .* density);

function xray_filters_test
kev = [20:200]';
mtype = {'lead', 'copper', 'aluminum'};
t = xray_filters(mtype, [0.01 0.2 0.1], kev);
if im
	clf, semilogy(kev, t, '-o'), axis([20 200 10^-2 1]),
	ir_legend(mtype)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function check = zipf_check ( a )

%*****************************************************************************80
%
%% ZIPF_CHECK checks the parameter of the Zipf PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the PDF.
%    1.0 < A.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( a <= 1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ZIPF_CHECK - Fatal error!\n' );
    fprintf ( 1, '  A <= 1.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
function y = square_pos( x )

%SQUARE_POS    Square of positive part.
%   SQUARE_POS(X) is the square of the postive parts of the elements of X;
%   i.e., SQUARE_POS(X)=MAX(X,0).^2. X must be real.
%
%   Disciplined convex programming information:
%       SQUARE_POS(X) is convex and nondecreasing in X. Thus when used in
%       CVX expressions, X must be convex (or affine).

narginchk(1,1);
if ~isreal( x ),
    error( 'Argument must be real.' );
end

y = square( max( x, 0 ) );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [iCM, cCM, oCM, sm] = Itti_Saliency(img,varargin)

verbose=0;
pictures=0;

% Load image
if verbose
    fprintf('Loading %s\n',filename);
end
image=img;
if pictures
    ShowImage(1,image,'Image');
end
image=double(image);

% Extract luminance and color channels
if verbose
    fprintf('Extracting early channels\n');
end

if length(size(image)) == 3
	[iIm,rIm,gIm,bIm,yIm]=ExtractChannels(image);

   % if max(max(rIm))>0 && max(max(gIm))>0 && max(max(bIm))>0
     if max(max(rIm))>0 && max(max(gIm))>0
        % Create pyramids
        if verbose
            fprintf('Creating pyramids\n');
        end
        iPyr=GaussianPyramid(iIm);
        rPyr=GaussianPyramid(rIm);
        gPyr=GaussianPyramid(gIm);
        bPyr=GaussianPyramid(bIm);
        yPyr=GaussianPyramid(yIm);
        oPyr=OrientationPyramid(iPyr);

        % Create feature maps
        if verbose
            fprintf('Creating feature maps\n');
        end
        iFM=IntensityFeatureMap(iPyr);
        [rgFM byFM]=ColorFeatureMap(rPyr,gPyr,bPyr,yPyr);
        oFM=OrientationFeatureMap(oPyr);

        % Create conspicuity maps
        if verbose
            fprintf('Conspicuity maps\n');
        end
        [iCM cCM oCM]=ConspicuityMap(iFM,oFM, rgFM,byFM);
        if pictures
            ShowImage(2,iCM,'Intensity CM');
            ShowImage(3,cCM,'Color CM');
            ShowImage(4,oCM,'Orientation CM');
        end

        % Create saliency map
        if verbose
            fprintf('Saliency map\n');
        end
        sm=SaliencyMap(iCM,oCM,cCM);

   %     s=size(image);
   %     s=s(1:2);
   %    sm=imresize(sm,s,'bilinear');
   %    sm=sm/max(max(sm))*255;
        if pictures
            ShowImage(5,sm,'Saliency');
        end
    else
        iIm = rgb2hsv(image);
        iIm = iIm(:,:,3);
        iPyr=GaussianPyramid(iIm);
        oPyr=OrientationPyramid(iPyr);

        if verbose
            fprintf('Creating feature maps\n');
        end
        iFM=IntensityFeatureMap(iPyr);
        oFM=OrientationFeatureMap(oPyr);

        % Create conspicuity maps
        if verbose
            fprintf('Conspicuity maps\n');
        end
        [iCM cCM oCM]=ConspicuityMap(iFM,oFM);
        if pictures
            ShowImage(2,iCM,'Intensity CM');
            ShowImage(3,oCM,'Orientation CM');
        end
        if verbose
            fprintf('Saliency map\n');
        end
        sm=SaliencyMap(iCM,oCM);
        s=size(image);
        s=s(1:2);
        sm=imresize(sm,s,'bilinear');
        if pictures
            ShowImage(4,sm,'Saliency');
        end
    end
else
        iIm = image;
        iPyr=GaussianPyramid(iIm);
        oPyr=OrientationPyramid(iPyr);

        if verbose
            fprintf('Creating feature maps\n');
        end
        iFM=IntensityFeatureMap(iPyr);
        oFM=OrientationFeatureMap(oPyr);

        % Create conspicuity maps
        if verbose
            fprintf('Conspicuity maps\n');
        end
        [iCM cCM oCM]=ConspicuityMap(iFM,oFM);
        if pictures
            ShowImage(2,iCM,'Intensity CM');
            ShowImage(3,oCM,'Orientation CM');
        end
        if verbose
            fprintf('Saliency map\n');
        end
        sm=SaliencyMap(iCM,oCM);
        s=size(image);
        s=s(1:2);
        sm=imresize(sm,s,'bilinear');
        if pictures
            ShowImage(4,sm,'Saliency');
        end
end

% ------------------------------------------------------------------------
% ExtractChannels
% ------------------------------------------------------------------------

function [iIm,rIm,gIm,bIm,yIm]=ExtractChannels(image)

r = image(:,:,1);
g = image(:,:,2);
b = image(:,:,3);

iIm = (r+g+b)/3;

% iIm = rgb2hsv(image);
% iIm = iIm(:,:,3);

% Normalize (see Itti, Koch & Niebur, p. 1255)
normalizer = iIm;
maxIm = max(max(iIm));

zeroentry = find(normalizer == 0);
normalizer(zeroentry) = 1e-10;
t = iIm < maxIm/10;
r = r ./ normalizer;
g = g ./ normalizer;
b = b ./ normalizer;
r(t) = 0;
g(t) = 0;
b(t) = 0;

% Channels R,G,B,Y
rIm = r - (g + b)/2;
gIm = g - (r + b)/2;
bIm = b - (r + g)/2; % negtive
yIm = (r + g)/2 - abs(r - g)/2 - b;
rIm(rIm<0) = 0;
gIm(gIm<0) = 0;
bIm(bIm<0) = 0;
yIm(yIm<0) = 0;

% ------------------------------------------------------------------------
% GaussianPyramid
% ------------------------------------------------------------------------

function pyramid = GaussianPyramid(image)

pyramid{1} = image;

for level=2:9
    im = gausmooth(pyramid{level-1});
    s=ceil(size(pyramid{level-1})/2.0);
	pyramid{level} = imresize(im,s);
end

% ------------------------------------------------------------------------
% GaussianSmooth
% ------------------------------------------------------------------------
function im = gausmooth(im)

[m,n] = size(im);
GaussianDieOff = .0001;
pw = 1:30;
ssq = 2;
width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');
if isempty(width)
    width = 1;  % the user entered a really small sigma
end
t = (-width:width);
gau = exp(-(t.*t)/(2*ssq))/sum(exp(-(t.*t)/(2*ssq)));

im = imfilter(im, gau,'conv','replicate');   % run the filter accross rows
im = imfilter(im, gau','conv','replicate'); % and then accross columns
% im = im/max(max(im));

% ------------------------------------------------------------------------
% OrientationPyramid
% ------------------------------------------------------------------------

function oPyr = OrientationPyramid(iPyr)

gabor{1}=gabor_fn(1,0.5,0,2,0);
gabor{2}=gabor_fn(1,0.5,0,2,pi/4);
gabor{3}=gabor_fn(1,0.5,0,2,pi/2);
gabor{4}=gabor_fn(1,0.5,0,2,pi*3/4);

% gabor{1}=gabor_fn(1,1,0,2.333,0);
% gabor{2}=gabor_fn(1,1,0,2.333,pi/4);
% gabor{3}=gabor_fn(1,1,0,2.333,pi/2);
% gabor{4}=gabor_fn(1,1,0,2.333,pi*3/4);

for l=3:9
    for o=1:4
        oPyr{l,o}=imfilter(iPyr{l},gabor{o},'symmetric');
    end
end

% ------------------------------------------------------------------------
% GaborFilterBank (with complex Gabors)
% ------------------------------------------------------------------------

function gb=gabor_fn(bw,gamma,psi,lambda,theta)
% bw    = bandwidth, (1)
% gamma = aspect ratio, (0.5)
% psi   = phase shift, (0)
% lambda= wave length, (>=2)
% theta = angle in rad, [0 pi)

sigma = lambda/pi*sqrt(log(2)/2)*(2^bw+1)/(2^bw-1);
sigma_x = sigma;
sigma_y = sigma/gamma;

sz=fix(8*max(sigma_y,sigma_x));
if mod(sz,2)==0, sz=sz+1;end

% alternatively, use a fixed size
% sz = 60;

[x y]=meshgrid(-fix(sz/2):fix(sz/2),fix(sz/2):-1:fix(-sz/2));
% x (right +)
% y (up +)

% Rotation
x_theta=x*cos(theta)+y*sin(theta);
y_theta=-x*sin(theta)+y*cos(theta);

gb=exp(-0.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+psi);
%imshow(gb/2+0.5);

function gabor = GaborFilterBank

% low-pass filter for 5 x 5 kernels
l=[1 3 8 3 1]/16.0;
lpf = l'*l;
ssd = 1;

%4 orientations
% for o=1:4
%     angle=(o-1)*pi/4;
%     c=cos(angle);
%     s=sin(angle);
%     for x=1:5
%         for y=1:5
%             m=pi/2*(c*(x-3)+s*(y-3));
%             kernel(y,x)=cos(m)+i*sin(m);
%         end
%     end
%     gabor{o}=kernel .* lpf;
% end

for o = 1:4
    sz_x=fix(6*sqrt(ssd));
    sz_y=fix(6*sqrt(ssd));

    [x y]=meshgrid(-fix(sz_x/2):fix(sz_x/2),fix(-sz_y/2):fix(sz_y/2));

    % Rotation
    angle=(o-1)*pi/4;
    x_theta=x*cos(angle)+y*sin(angle);
    y_theta=-x*sin(angle)+y*cos(angle);

    match = find(x_theta==0);
    x_theta(match) = x_theta(match) + 1e-10;
    gabor{o}=exp(-.5*(x_theta.^2/ssd+y_theta.^2/ssd)).*cos(2*pi./x_theta);
end

% ------------------------------------------------------------------------
% IntensityFeatureMap
% ------------------------------------------------------------------------

function iFM = IntensityFeatureMap(iPyr)

for c=3:5
    for delta = 3:4
        s = c + delta;
        iFM{c,s} = abs(Subtract(iPyr{c}, iPyr{s}));
    end
end

% ------------------------------------------------------------------------
% ColorFeatureMap
function [rgFM, byFM] = ColorFeatureMap(rPyr, gPyr, bPyr, yPyr)

for c=3:5
    for delta = 3:4
        s = c + delta;
        rgFM{c,s} = abs(Subtract(rPyr{c}-gPyr{c}, gPyr{s}-rPyr{s}));
        byFM{c,s} = abs(Subtract(bPyr{c}-yPyr{c}, yPyr{s}-bPyr{s}));
    end
end

% ------------------------------------------------------------------------
% OrientationFeatureMap
% ------------------------------------------------------------------------

function oFM = OrientationFeatureMap(oPyr)

for c=3:5
    for delta = 3:4
        s = c + delta;
        for o=1:4
            oFM{c,s,o}=abs(Subtract(oPyr{c,o},oPyr{s,o}));
        end
    end
end

% ------------------------------------------------------------------------
% ConspicuityMap
% ------------------------------------------------------------------------

function [iCM,cCM,oCM] = ConspicuityMap(varargin)

iFM = varargin{1};
oFM = varargin{2};

dim=size(iFM{3,6});
iCM=zeros(dim);
for c=3:5
    for delta = 3:4
        s = c + delta;
%         weight=s*c;
        weight=1;
        iCM = Add(iCM,weight*Normalize(iFM{c,s}));
    end
end

oCM=zeros(dim);
for c=3:5
    for delta = 3:4
        s = c + delta;
%         weight=s*c;
        weight=1;
        for o=1:4
            oCM=Add(oCM,weight*Normalize(oFM{c,s,o}));
        end
    end
end

cCM=zeros(dim);
if length(varargin) == 4
    rgFM = varargin{3};
    byFM = varargin{4};
    for c=3:5
        for delta = 3:4
            s = c + delta;
    %         weight=s*c;
            weight=1;
            cCM = Add(cCM,weight*Normalize(rgFM{c,s}));
            cCM = Add(cCM,weight*Normalize(byFM{c,s}));
        end
    end
end

% ------------------------------------------------------------------------
% Normalize
% ------------------------------------------------------------------------

function normalized = Normalize(map)

% Normalize map to range [0..1]
minValue = min(min(map));
map = map-minValue;
maxValue = max(max(map));
if maxValue>0
    map = map/maxValue;
end

% Position of local maxima
lmax = LocalMaxima(map);

% Position of global maximum
gmax = (map==1.0);

% Local maxima excluding global maximum
lmax = lmax .* (gmax==0);

% Average of local maxima excluding global maximum
nmaxima=sum(sum(lmax));
if nmaxima>0
    m = sum(sum(map.*lmax))/nmaxima;
else
    m = 0;
end
normalized = map*(1.0-m)^2;

% ------------------------------------------------------------------------
% LocalMaxima
% ------------------------------------------------------------------------

function maxima = LocalMaxima(A)

nRows=size(A,1);
nCols=size(A,2);
% compare with bottom, top, left, right
maxima =           (A > [A(2:nRows, :);   zeros(1, nCols)]);
maxima = maxima .* (A > [zeros(1, nCols); A(1:nRows-1, :)]);
maxima = maxima .* (A > [zeros(nRows, 1), A(:, 1:nCols-1)]);
maxima = maxima .* (A > [A(:, 2:nCols),   zeros(nRows, 1)]);

% ------------------------------------------------------------------------
% Subtract
% ------------------------------------------------------------------------

function result = Subtract(im1, im2)

im2 = imresize(im2, size(im1), 'bilinear');
result = im1 - im2;

% ------------------------------------------------------------------------
% Add
% ------------------------------------------------------------------------

function result = Add(im1, im2)

im2 = imresize(im2, size(im1), 'bilinear');
result = im1 + im2;

% ------------------------------------------------------------------------
% SaliencyMap
% ------------------------------------------------------------------------

function sm=SaliencyMap(varargin)

iCM = varargin{1};
oCM = varargin{2};
if length(varargin) == 2
    sm=(Normalize(iCM)+Normalize(oCM))/2;
else
    cCM = varargin{3};
    sm=(Normalize(iCM)+Normalize(cCM)+Normalize(oCM))/3;
end

% ------------------------------------------------------------------------
% ShowImage
% ------------------------------------------------------------------------

function ShowImage(nFigure,image,fTitle)

% figure(nFigure);
figure;
imagesc(image);
if (size(image, 3) == 1)
    colormap('gray');
    colorbar;
end
axis image;
title(fTitle);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [fl re]=lines(im_texto)
% Divide text in lines
% im_texto->input image; fl->first line; re->remain line
% Example:
% im_texto=imread('TEST_3.jpg');
% [fl re]=lines(im_texto);
% subplot(3,1,1);imshow(im_texto);title('INPUT IMAGE')
% subplot(3,1,2);imshow(fl);title('FIRST LINE')
% subplot(3,1,3);imshow(re);title('REMAIN LINES')
im_texto=clip(im_texto);
num_filas=size(im_texto,1);
for s=1:num_filas
    if sum(im_texto(s,:))==0
        nm=im_texto(1:s-1, :); % First line matrix
        rm=im_texto(s:end, :);% Remain line matrix
        fl = clip(nm);
        re=clip(rm);
        %*-*-*Uncomment lines below to see the result*-*-*-*-
        %         subplot(2,1,1);imshow(fl);
        %         subplot(2,1,2);imshow(re);
        break
    else
        fl=im_texto;%Only one line.
        re=[ ];
    end
end

function img_out=clip(img_in)
[f c]=find(img_in);
img_out=img_in(min(f):max(f),min(c):max(c));%Crops image


text:
%% Solves the Krusell and Smith (1998)
% Solving Krusell-Smith model using perturbation method involves simple 5
% step process given by
%
% # Solve for Steady State
% # Linearize Model Equations
% # Solve out Static Constraint (optional reduce model)
% # Solve Linear System
% # Compute Impulse Response Functions
% # Check for Internal Consistency (if model reduction is used)
% where each step just requires about one function call.
%
% In this example file, the entire problem will be solved total of 3 times.
% First, KS model will be solved without any model reduction. Second, KS
% model will be solved with model reduction, to compare how well the model
% reduction works. Lastly, different parts of observables will be
% introduced to compare how the model reduction works with different
% requirements. We have updates to Internal Consistency check coming in the
% future.
%
% Estimated Runtime: 1.5 seconds
%
% REFERENCES:
%
% * Ahn, SeHyoun, Greg Kaplan, Benjamin Moll, Thomas Winberry, and
% Christian Wolf. "When Inequality Matters for Macro and Macro Matters for
% Inequality."
% * Krusell, Per, and Anthony A. Smith, Jr. "Income and wealth heterogeneity in the macroeconomy." Journal of political Economy 106.5 (1998): 867-896.
%
% REQUIRES:
%
% * auto diff toolbox: <https://github.com/sehyoun/MATLABAutoDiff>
% * phact toolbox: <https://github.com/gregkaplan/phact>
% * compute_steady_state.m
% * equilibrium_conditions.m
% * plot_steady_state.m
% * plot_IRFs.m
% Relevant files except for auto diff toolbox, can be found in
%    /examples/KrusellSmith folder of phact toolbox provided on github
%
%%

%% Setup the toolbox
% Just need to include folders containing the files in the path.

addpath('/home/sehyoun/Dropbox/0Packages/AutoDiff');
addpath('/home/sehyoun/Dropbox/0Packages/PHACT');

%% Set options for this example run
% example_case : 1   run full and compare with only value function reduction
%                2   run full and compare with only state space reduction
%                3   run full and compare with both reductions
example_case = 3;
reduceDist_hor = 50;        % m of K_m(A,b)

% initialize shocks for simulation
T = 200;
N = 2000;
vAggregateShock = zeros(1,N);
vAggregateShock(1,1) = 1;
%vAggregateShock = randn(1,N);     % Uncomment for one realization instead
                                  % of IRFs

%% Step 0: Set Parameters
% The script sets up parameters relevant for the model
set_parameters;

%% Step 1: Solve for Steady State
% Non-stochastic steady state can be found using any methods. In
%    particular, example codes can be found at
%    <http://www.princeton.edu/~moll/HACTproject.htm>.

tStart = tic;
fprintf('Computing steady state...\n')
global IfSS IbSS I0SS varsSS A

[rSS,wSS,KSS,ASS,uSS,cSS,VSS,gSS,dVUSS,dVfSS,dVbSS,IfSS,IbSS,I0SS] = ...
    compute_steady_state();

fprintf('Time to compute steady state: %.3g seconds\n\n\n',toc(tStart));

% Store steady state values
varsSS = zeros(nVars,1);
varsSS(1:2*I,1) = reshape(VSS,2*I,1);
ggSS = reshape(gSS,2*I,1);
varsSS(2*I+1:4*I-1,1) = ggSS(1:2*I-1);
varsSS(4*I,1) = 0;
varsSS(4*I+1,1) = KSS;
varsSS(4*I+2,1) = rSS;
varsSS(4*I+3,1) = wSS;
varsSS(4*I+4,1) = (KSS ^ aalpha) * (zAvg ^ (1 - aalpha));
CSS = sum(cSS(:) .* gSS(:) * da);
varsSS(4*I+5,1) = CSS;
varsSS(4*I+6,1) = ddelta * KSS;

% plot steady state results
plot_steady_state;

%% Step 2: Linearize Model Equations
% For computing derivatives, the codes written for solving for the
%    steady-state can be used almost verbatim using automatic
%    differentiation toolbox as long as only the functions supported by
%    automatic differentation are used. For list of supported functions and
%    documentation of relevant syntax check
%    <https://github.com/sehyoun/MATLABAutoDiff>. Example usage/syntax of
%    automatic differentiation can be found at
%    <http://sehyoun.com/EXAMPLE_AutoDiff_Syntax.html>
fprintf('Taking derivatives of equilibrium conditions...\n')
t0 = tic;

% Prepare automatic differentiation
vars = zeros(2*nVars+nEErrors+1,1);
vars = myAD(vars);

% Evaluate derivatives
derivativesIntermediate = equilibrium_conditions(vars);

% Extract out derivative values
derivs = getderivs(derivativesIntermediate);

tDerivs = toc(t0);
fprintf('...Done!\n')
fprintf('Time to compute derivatives: %2.4f seconds\n\n\n',tDerivs)
if tDerivs > 1
    warning('If you do not compile mex files for automatics differentiation, matrix vector multiplication will be slow');
    disp('Press any key to continue...');
    pause();
end

% Unpackage derivatives
mVarsDerivs = derivs(:,1:nVars);
mVarsDotDerivs = derivs(:,nVars+1:2*nVars);
mEErrorsDerivs = derivs(:,2*nVars+1:2*nVars+nEErrors);
mShocksDerivs = derivs(:,2*nVars+nEErrors+1);

%% Step 3 (full): Solve out Static Constraints and/or Reduce Models
% For the first run, a unreduced model will be solved

% rename derivative matrix to match notation in paper
g0 = mVarsDotDerivs;
g1 = -mVarsDerivs;
c = sparse(nVars,1);
psi = -mShocksDerivs;
pi = -mEErrorsDerivs;

% Solve out static constratins
[~,inv_state_red,g0,g1,constant,pi,psi] = clean_G0_sparse(g0,g1,c,pi,psi);

%% Step 4: Solve Linear System
t0 = tic;
fprintf('Solving linear system...\n')

[G1,~,impact,eu] = schur_solver(g0,g1,c,psi,pi,1,1,1);

fprintf('...Done!\n')
fprintf('Existence and uniqueness? %2.0f and %2.0f\n',eu);
fprintf('Time to solve full linear system: %2.4f seconds\n\n\n',toc(t0))

%% Step 5: Simulate Impulse Response Functions
t0 = tic;

[simulated,~] = simulate(G1,impact,T,N,vAggregateShock,'implicit',inv_state_red);%,4*I:4*I+6);

fprintf('...Done!\n')
fprintf('Time to simulate model: %2.4f seconds\n\n\n',toc(t0))

big_simul = simulated + varsSS;
simulated = simulated(4*I:4*I+6,:);

varsSS_small = varsSS(4*I:4*I+6,1);
% Add state-states back in to get values in levels
vAggregateTFP			= simulated(1,:) + varsSS_small(1);
vAggregateOutput		= simulated(5,:) + varsSS_small(5);
vAggregateConsumption	= simulated(6,:) + varsSS_small(6);
vAggregateInvestment 	= simulated(7,:) + varsSS_small(7);

% Compute log differences for plotting
vAggregateTFP_full		= vAggregateTFP;
vAggregateOutput_full	= log(vAggregateOutput) - log(varsSS_small(5));
vAggregateConsumption_full	= log(vAggregateConsumption) - log(varsSS_small(6));
vAggregateInvestment_full	= log(vAggregateInvestment) - log(varsSS_small(7));

%% Step 3 (reduced): Solve Out Static Constraints and/or Reduce Models

% Set relevant parameters for example cases
if (example_case == 1)
    reduceV = 1;
    reduceDistribution = 0;
elseif (example_case == 2)
    reduceV = 0;
    reduceDistribution = 1;
else
    reduceV = 1;
    reduceDistribution = 1;
end

% rename derivatives to match notation in paper
g0 = mVarsDotDerivs;
g1 = -mVarsDerivs;
c = sparse(nVars,1);
psi = -mShocksDerivs;
pi = -mEErrorsDerivs;

t0 = tic;
fprintf('Model Reduction ...\n')

%%
% *State space reduction using Krylov subspace method*
if reduceDistribution == 1
    % State space reduction
    [state_red,inv_state_red,n_g_red] = krylov_reduction(g0,g1,n_v,n_g,reduceDist_hor);
    [g1,psi,pi,c,g0] = change_basis(state_red,inv_state_red,g1,psi,pi,c,g0);
else
    % Clean G0
    [state_red,inv_state_red,g0,g1,c,pi,psi] = clean_G0_sparse(g0,g1,c,pi,psi);
    n_g_red = n_g;
end

%%
% *Value function reduction using spline inspired bases*
if reduceV == 1
    % Create knot points for spline (the knot points are not uniformly spaced)
    knots = linspace(amin,amax,n_knots-1)';
    knots = (amax-amin)/(2^c_power-1)*((knots-amin)/(amax-amin)+1).^c_power+amin-(amax-amin)/(2^c_power-1);
    % Function calls to create basis reduction
    [from_spline, to_spline] = oneDquad_spline(x,knots);
    [from_spline, to_spline] = extend_to_nd(from_spline,to_spline,n_prior,n_post);
    n_splined = size(from_spline,2);
    [from_spline, to_spline] = projection_for_subset(from_spline,to_spline,0,n_g_red);

    % Reduce the decision vector
    [g1,psi,~,c,g0] = change_basis(to_spline,from_spline,g1,psi,pi,c,g0);
    pi = to_spline * pi * from_spline(1:n_v,1:n_splined);
elseif reduceV == 0
    from_spline = speye(n_g_red + n_v);
    to_spline = speye(n_g_red + n_v);
    n_splined = n_v;
end
fprintf('...Done!\n')
fprintf('Time to reduce dimensionality: %2.4f seconds\n\n\n',toc(t0))

%% Step 4: Solve Linear System
t0 = tic;
fprintf('Solving reduced linear system...\n')

[G1,~,impact,eu,F] = schur_solver(g0,g1,c,psi,pi,1,1,1);

fprintf('...Done!\n')
fprintf('Existence and uniqueness? %2.0f and %2.0f\n',eu);
fprintf('Time to solve linear system: %2.4f seconds\n\n\n',toc(t0))

%% Step 5: Simulate Impulse Response Functions
fprintf('Simulating Model...\n')
t0 = tic;

trans_mat = inv_state_red*from_spline;
[simulated,vTime] = simulate(G1,impact,T,N,vAggregateShock,'implicit',trans_mat,4*I:4*I+6);

fprintf('...Done!\n')
fprintf('Time to simulate model: %2.4f seconds\n\n\n',toc(t0))

% Add state-states back in to get values in levels
vAggregateTFP = simulated(1,:) + varsSS_small(1);
vAggregateOuput = simulated(5,:) + varsSS_small(5);
vAggregateConsumption = simulated(6,:) + varsSS_small(6);
vAggregateInvestment = simulated(7,:) + varsSS_small(7);

% Compute log differences for plotting
vAggregateTFP_reduced = vAggregateTFP;
vAggregateOutput_reduced = log(vAggregateOutput) - log(varsSS_small(5));
vAggregateConsumption_reduced = log(vAggregateConsumption) - log(varsSS_small(6));
vAggregateInvestment_reduced = log(vAggregateInvestment) - log(varsSS_small(7));

%% Step 6: Internal Consistency Check
% For large problem, we are forced to take explicit update, so it can take
% awhile to run, but the good thing is that this only needs to be run at
% the end. We are currently working on this, so speed for this part might
% improve in the future. (For small problems, the speed is not an issue). A
% different consistency check is in works, and will be updated in the
% future.

g1 = -mVarsDerivs;
psi = -mShocksDerivs;
from_red = inv_state_red * from_spline;
to_red = to_spline * state_red;
[epsilon] = internal_consistency_check(G1,impact,n_g_red,from_red,to_red,g1,psi,F,n_v,n_g,1000,varsSS,1,0);

%% (optional) Step 7: Plot relevant values
% Plot impulse response functions
plot_IRFs;

%% State Space Reduction with More Observables
% It is possible to include other variables into observables. In this part,
%    different parts of the distribution will be included in the model
%    reduction to show how the Krylov subspace based model reduction
%    behaves. To that regard, the reduced model is solved while including
%    different parts of distribution as part of the observable. In this
%    case, a very strong requirement of knowing the probability density
%    value at different grid points of the distribution g is taken. In the
%    first case, 80th through 89th grid points were included in the parts
%    that we consider to be observable. For the second example, 85th
%    through 94th grid points were included in observable.

% rename derivative matrix to match notation in paper
g0 = mVarsDotDerivs;
g1 = -mVarsDerivs;
c = sparse(nVars,1);
psi = -mShocksDerivs;
pi = -mEErrorsDerivs;

% Solve full model
[~,inv_state_red,g0_full,g1_full,c_full,pi_full,psi_full] = clean_G0_sparse(g0,g1,c,pi,psi);
[G1,~,impact,eu] = schur_solver(g0_full,g1_full,c_full,psi_full,pi_full,1,1,1);
[simulated,vTime] = simulate(G1,impact,T,N,vAggregateShock,'implicit',inv_state_red);
simulated_full = simulated + varsSS;

% Solve reduced without extra observables
[state_red,inv_state_red,n_g_red] = krylov_reduction(g0,g1,n_v,n_g,20);
[g1_red,psi_red,pi_red,c_red,g0_red] = change_basis(state_red,inv_state_red,g1,psi,pi,c,g0);
[G1,~,impact,eu] = schur_solver(g0_red,g1_red,c_red,psi_red,pi_red,1,1,1);
[simulated,vTime] = simulate(G1,impact,T,N,vAggregateShock,'implicit',inv_state_red);
simulated_red = simulated + varsSS;

%%
% Solve reduced with extra observables for 80th through 89th grid point

[state_red,inv_state_red,n_g_red] = krylov_reduction(g0,g1,n_v,n_g,20,spdiags(ones(n_g,1),79,10,n_g));
[g1_80,psi_80,pi_80,c_80,g0_80] = change_basis(state_red,inv_state_red,g1,psi,pi,c,g0);
[G1,~,impact,eu] = schur_solver(g0_80,g1_80,c_80,psi_80,pi_80,1,1,1);
[simulated,vTime] = simulate(G1,impact,T,N,vAggregateShock,'implicit',inv_state_red);
simulated_obs = simulated + varsSS;
diff_red = (abs(simulated_red(n_v+1:end,:) - simulated_full(n_v+1:end,:)))./simulated_full(n_v+1:end,:);
diff_obs80 = (abs(simulated_obs(n_v+1:end,:) - simulated_full(n_v+1:end,:)))./simulated_full(n_v+1:end,:);

figure;
plot(70:99,log(diff_red(70:99,11)));
hold on;
plot(70:99,log(diff_obs80(70:99,11)));
grid on
legend('Location','northwest');
legend('without','with grid 80~89 included');
title(['Error Comparison at time ',num2str(vTime(11))],'interpreter','latex','fontsize',14)
ylabel('$\log$(relative error)','interpreter','latex')
xlabel('Distribution Location','interpreter','latex')

%%
% Solve reduced with extra observables for 85th through 94th grid point

[state_red,inv_state_red,n_g_red] = krylov_reduction(g0,g1,n_v,n_g,20,spdiags(ones(n_g,1),84,10,n_g));
[g1,psi,pi,c,g0] = change_basis(state_red,inv_state_red,g1,psi,pi,c,g0);
[G1,~,impact,eu] = schur_solver(g0,g1,c,psi,pi,1,1,1);
[simulated,vTime] = simulate(G1,impact,T,N,vAggregateShock,'implicit',inv_state_red);
simulated_obs = simulated + varsSS;
diff_obs = (abs(simulated_obs(n_v+1:end,:) - simulated_full(n_v+1:end,:)))./simulated_full(n_v+1:end,:);

figure;
plot(70:99,log(diff_red(70:99,11)));
hold on;
plot(70:99,log(diff_obs(70:99,11)));
grid on
legend('Location','northwest');
legend('without','with grid 85~94 included');
title(['Error Comparison at time ',num2str(vTime(11))],'interpreter','latex','fontsize',14)
ylabel('$\log$(relative error)','interpreter','latex')
xlabel('Distribution Location','interpreter','latex')

%%
% Note how the errors decreases only for the parts that we required to be
% observable while the irrelevant parts were reduced away. This is the
% magic of Krylov subspace based reduction in action. Though we could solve
% the full model for the simple Krusell-Smith case, this ability to reduce
% away irrelevant part will be essential for bigger problems.
%

%%
% Lastly, some of the functions have lengthy parameters, you can always
% type "help function_name;" to see documentation. All functions contain the
% function call syntax as the last line, and thoese can be copied and
% pasted into the program file. For example, to use the
% krylov_reduction function, you can call
help krylov_reduction;

%%
% and just copy the lastline of the documentation string. Alternately, you
% can call "doc function_name" to read the doc string in the MATLAB help
% browser, and "edit function_name" to see the function codes.

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function nx= sizeXrec(i,P)
% returns the number of k-ToM's hidden sates
% function nx = sizeXrec(i,P)
% Marie Devaine wrote this in November 2015 (comments: JD).
% The formula follows from an analysis of the recursive sequence.
% IN:
%   - i: k-ToM's sophistication level
%   - P: total nb of evol/obs params
% OUT:
%   - nx: number of k-ToM's hidden sates

if i==0
    nx = 2;
else
    i_1 = i-1;
    if i_1>0
        power2 = exp(((i_1-1):-1:0).*log(2));
        S = 2^i + power2*((1:i_1)*(3*P+2)-1)';
    else
        S = 2;
    end
    nx = S+ i*(3*P+2)-1;
end
end

author: MBB-team

repo: VBA-toolbox


text:
function [ n_data, a, x, fx ] = gamma_inc_values ( n_data )

%*****************************************************************************80
%
%% GAMMA_INC_VALUES returns some values of the incomplete Gamma function.
%
%  Discussion:
%
%    The (normalized) incomplete Gamma function P(A,X) is defined as:
%
%      PN(A,X) = 1/Gamma(A) * Integral ( 0 <= T <= X ) T**(A-1) * exp(-T) dT.
%
%    With this definition, for all A and X,
%
%      0 <= PN(A,X) <= 1
%
%    and
%
%      PN(A,INFINITY) = 1.0
%
%    In Mathematica, the function can be evaluated by:
%
%      1 - GammaRegularized[A,X]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, the parameter of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  a_vec = [ ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.10E+01, ...
     0.10E+01, ...
     0.10E+01, ...
     0.11E+01, ...
     0.11E+01, ...
     0.11E+01, ...
     0.20E+01, ...
     0.20E+01, ...
     0.20E+01, ...
     0.60E+01, ...
     0.60E+01, ...
     0.11E+02, ...
     0.26E+02, ...
     0.41E+02  ];

  fx_vec = [ ...
     0.7382350532339351E+00, ...
     0.9083579897300343E+00, ...
     0.9886559833621947E+00, ...
     0.3014646416966613E+00, ...
     0.7793286380801532E+00, ...
     0.9918490284064973E+00, ...
     0.9516258196404043E-01, ...
     0.6321205588285577E+00, ...
     0.9932620530009145E+00, ...
     0.7205974576054322E-01, ...
     0.5891809618706485E+00, ...
     0.9915368159845525E+00, ...
     0.1018582711118352E-01, ...
     0.4421745996289254E+00, ...
     0.9927049442755639E+00, ...
     0.4202103819530612E-01, ...
     0.9796589705830716E+00, ...
     0.9226039842296429E+00, ...
     0.4470785799755852E+00, ...
     0.7444549220718699E+00 ];

  x_vec = [ ...
     0.30E-01, ...
     0.30E+00, ...
     0.15E+01, ...
     0.75E-01, ...
     0.75E+00, ...
     0.35E+01, ...
     0.10E+00, ...
     0.10E+01, ...
     0.50E+01, ...
     0.10E+00, ...
     0.10E+01, ...
     0.50E+01, ...
     0.15E+00, ...
     0.15E+01, ...
     0.70E+01, ...
     0.25E+01, ...
     0.12E+02, ...
     0.16E+02, ...
     0.25E+02, ...
     0.45E+02 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [clusters,m,mnames] = cluster_FA(clusters,OPT)
% [clusters,m,mnames]  = cluster_FA(clusters,OPT)
% Gets independent components for canonical timeseries values from clusters
%
% Requires OPT structure with the following fields (examples shown here):
%   O.HP = 100;
%   O.TR = 1.5;
%   O.doHP = 1;
%   O.doLP = 0;
%   O.scanadjust = 1;
%   O.percent = 0;
% 	O.filtertype = 'spm';
% 	O.nruns = 2;
%   O.adjustmatrix = custom adjustment matrix to regress out (e.g., movement params)
%   O.trimts = 3;
%
% Recommended now: average, or FA.  ICA can sometimes reverse sign, so can
% PCA potentially.  Not sure about sign thing.  ICA and PCA enforce
% orthogonal canonical variates.
%
% Returns: clusters, with pcscore, ica, and varimax fields
%          average of voxels is stored in timeseries
%
% Maybe the best thing is to take voxels that load highly on ICs or PCs and
% return selective averages of those subsets - so gives subclusters.
%
% m is a matrix of rows = time (or observations), cols = timecourses within
% and across regions
%
% see also cluster_princomp.m

mnames = {};

for i = 1:length(clusters)

    m = clusters(i).all_data;
    m2 = roi_timeseries(m,OPT);

    close all

    clusters(i).pcscore = m2.pcscore;
    clusters(i).ica = m2.ica;
    clusters(i).varimax = m2.varimax;

    tmp = [clusters(i).title '_' num2str(i) '_' num2str(clusters(i).mm_center(1)) num2str(clusters(i).mm_center(2)) num2str(clusters(i).mm_center(3))];
    tmp(tmp==' ' | tmp=='.') = ['_'];

    for j = 1:size(clusters(i).varimax,2)
        mnames = [mnames {tmp}];
    end
end

m = cat(2,clusters.varimax);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function value = r8_besk0 ( x )

%*****************************************************************************80
%
%% R8_BESK0 evaluates the Bessel function K of order 0 of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the Bessel function K of order 0 of X.
%
  persistent bk0cs
  persistent ntk0
  persistent xmax
  persistent xsml

  if ( isempty ( ntk0 ) )

    bk0cs = [ ...
      -0.353273932339027687201140060063153E-01, ...
      +0.344289899924628486886344927529213, ...
      +0.359799365153615016265721303687231E-01, ...
      +0.126461541144692592338479508673447E-02, ...
      +0.228621210311945178608269830297585E-04, ...
      +0.253479107902614945730790013428354E-06, ...
      +0.190451637722020885897214059381366E-08, ...
      +0.103496952576336245851008317853089E-10, ...
      +0.425981614279108257652445327170133E-13, ...
      +0.137446543588075089694238325440000E-15, ...
      +0.357089652850837359099688597333333E-18, ...
      +0.763164366011643737667498666666666E-21, ...
      +0.136542498844078185908053333333333E-23, ...
      +0.207527526690666808319999999999999E-26, ...
      +0.271281421807298560000000000000000E-29, ...
      +0.308259388791466666666666666666666E-32 ]';

    ntk0 = r8_inits (bk0cs, 16, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r8_mach ( 3 ) );
    xmax = - log ( r8_mach ( 1 ) );
    xmax = xmax - 0.5 * xmax * log ( xmax ) / ( xmax + 0.5 );

  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_BESK0 = Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R8_BESK0 = Fatal error!' )
  elseif ( x <= xsml )
    y = 0.0;
    value = - log ( 0.5 * x ) * r8_besi0 ( x ) ...
      - 0.25 + r8_csevl ( 0.5 * y - 1.0, bk0cs, ntk0 );
  elseif ( x <= 2.0 )
    y = x * x;
    value = - log ( 0.5 * x ) * r8_besi0 ( x ) ...
      - 0.25 + r8_csevl ( 0.5 * y - 1.0, bk0cs, ntk0 );
  elseif ( x <= xmax )
    value = exp ( - x ) * r8_besk0e ( x );
  else
    value = 0.0;
  end

  return
end


text:
classdef KNN < Algorithm
    %KNN Basic k-nearest neighbors algorithm based on Euclidean distance

    properties
        description = 'k-nearest neighbors algorithm';
        % Parameters to optimize and default value
        parameters = struct('k', 5);
    end

    methods
        function obj = KNN(varargin)
            %KNN constructs an object of the class KNN. Default k is 5
            %
            %   OBJ = KNN('k', neighbours)
            %   builds KNN with NEIGHBOURS as number of neighbours to consider
            %   to label new patterns.
            obj.parseArgs(varargin);
        end

        function [projectedTrain, predictedTrain]= privfit( obj, train, param)
            if(nargin == 3)
                obj.parameters.k = param.k;
            end

            % save train data in the model structure
            obj.model.train = train;
            obj.model.parameters = obj.parameters;
            % Predict train labels
            [projectedTrain, predictedTrain] = predict(obj, train.patterns);
        end

        function [projected, predicted] = privpredict(obj, testPatterns)
            % Variables aliases
            x = obj.model.train.patterns;
            xlabel = obj.model.train.targets;
            k = obj.model.parameters.k;

            dist = pdist2(testPatterns,x);
            % indicies of nearest neighbors
            [~,nearest] = sort(dist,2);
            % k nearest
            nearest = nearest(:,1:k);
            % mode of k nearest
            val = xlabel(nearest);
            predicted = mode(val,2);

            % dummy value for projections
            projected = -1.*ones(length(testPatterns),1);
        end
    end
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
% Distortionless data hiding based on integer wavelet transform (Watermark Extraction)
% Guorong Xuan   Jiang Zhu   Jidong Chen   Shi, Y.Q.   Zhicheng Ni   Wei Su
% Department of Computer Science, Tongji University, Shanghai
%
% This paper appears in: Electronics Letters
% Publication Date: 5 Dec 2002
% Volume: 38,  Issue: 25
% On page(s): 1646- 1648
% ISSN: 0013-5194

% August 19, 2007, 3:38pm
% updated 13 December, 2008
% By: Asad (asad_82@yahoo.com)

clear all;
close all;

load WatermarkInfo.mat;
disp('------------------------ Extraction -------------------------------');

% read image and convert to gray scale if necessary
% a = imread('Watermarked Image.bmp');
% get the single channel
% a = a(:,:,1);
a = watermarkedImage;

% STEP 1: perfom integer wavelet decomposition
LS = liftwave('cdf2.2','Int2Int');
[CA,CH,CV,CD] = lwt2(double(a),LS);

% STEP 2: extract the embeded signal from 5th bit of CH, CV and CD
index = 1;
% preallocate memory to fasten things
embededSignal = zeros(size(CH,1)*size(CH,2)+size(CV,1)*size(CV,2)+size(CD,1)*size(CD,2),1);
for i=1:size(CH,1)
    for j=1:size(CH,2)
        % for constructing binary image using CH
        if CH(i,j) ~= -ERROR_NUM
            binSeq = dec2bin(abs(CH(i,j)),8);
            if binSeq(BITPLANE_NUMBER) == '1'
                embededSignal(index,1) = 1;
            else
                embededSignal(index,1) = 0;
            end
            index = index + 1;
        end
    end
end

% extract the embeded signal from 5th bit of CV
for i=1:size(CV,1)
    for j=1:size(CV,2)
        % for constructing binary image using CV
        if CV(i,j) ~= -ERROR_NUM
            binSeq = dec2bin(abs(CV(i,j)),8);
            if binSeq(BITPLANE_NUMBER) == '1'
                embededSignal(index,1) = 1;
            else
                embededSignal(index,1) = 0;
            end
            index = index + 1;
        end
    end
end

% extract the embeded signal from 5th bit of CD
for i=1:size(CD,1)
    for j=1:size(CD,2)
        % for constructing binary image using CD
        if CD(i,j) ~= -ERROR_NUM
            binSeq = dec2bin(abs(CD(i,j)),8);
            if binSeq(BITPLANE_NUMBER) == '1'
                embededSignal(index,1) = 1;
            else
                embededSignal(index,1) = 0;
            end
            index = index + 1;
        end
    end
end

% STEP 3: extract header information
% obtain count1 of CH for use in aritmatic decoding
binNum = '';
for i=1:8
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
count1(1,1) = num;

binNum = '';
for i=9:16
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
count1(1,2) = num;

% obtain count2 of CV for use in aritmatic decoding
binNum = '';
for i=17:24
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
count2(1,1) = num;

binNum = '';
for i=25:32
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
count2(1,2) = num;

% obtain count3 of CD for use in aritmatic decoding
binNum = '';
for i=33:40
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
count3(1,1) = num;

binNum = '';
for i=41:48
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
count3(1,2) = num;

% obtain length of compressed CH
binNum = '';
for i=49:64
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
lenCH5 = num;

% obtain length of compressed CV
binNum = '';
for i=65:80
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
lenCV5 = num;

% obtain length of compressed CD
binNum = '';
for i=81:96
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
lenCD5 = num;

% obtain length of watermark
binNum = '';
for i=97:128
    if embededSignal(i,1) == 1
        binNum = strcat(binNum,'1');
    else
        binNum = strcat(binNum,'0');
    end
end
num = bin2dec(binNum);
lenWatermark = num;

% STEP 4: Extract the respective signals CH, CV, CD and decode(uncompress) to get the original bit sequence
HL = 128;
% construct sequence 1(CH) and decode
seq11 = embededSignal(HL+1:(HL+lenCH5));
CH5 = arithdeco(seq11,count1,size(CH,1)*size(CH,2));

% construct sequence 2(CV) and decode
seq22 = embededSignal((HL+lenCH5)+1:(HL+lenCH5+lenCV5));
CV5 = arithdeco(seq22,count2,size(CV,1)*size(CV,2));

% construct sequence 3(CD) and decode
seq33 = embededSignal((HL+lenCH5+lenCV5)+1:(HL+lenCH5+lenCV5+lenCD5));
CD5 = arithdeco(seq33,count3,size(CD,1)*size(CD,2));

% additional step is to verify the embeded watermark so extract it (not compared to the original but can be)
watermarkE = embededSignal((HL+lenCH5+lenCV5+lenCD5)+1:(HL+lenCH5+lenCV5+lenCD5+lenWatermark));
if isequal(watermark,watermarkE')
    disp('Watermark correct');
else
    disp('Watermark Error')
end
watermarkE = reshape(watermarkE,WM_SIZE,WM_SIZE);
figure,imshow(watermarkE,[]),title('Retrieved Watermark')

% STEP 5: Restore the Image i.e. remove the watermark and insert the
% uncompressed 5th bit data back into the image
CH5 = reshape(CH5,size(CH,1),size(CH,1));
CV5 = reshape(CV5,size(CV,1),size(CV,1));
CD5 = reshape(CD5,size(CD,1),size(CD,1));

neg = 0;

for x=1:size(CH,1)
    for y=1:size(CH,2)
        if CH(x,y) ~= -ERROR_NUM
            % restore 5th bit of CH
            neg = 0;
            if CH(x,y) < 0
                neg = 1;
            end
            binSeq = dec2bin(abs(CH(x,y)),8);
            if CH5(x,y) == 2
                binSeq(BITPLANE_NUMBER) = '1';
            else
                binSeq(BITPLANE_NUMBER) = '0';
            end
            num = bin2dec(binSeq);
            if neg == 1
                CH(x,y) = num * -1;
            else
                CH(x,y) = num;
            end
        end
        % restore 5th bit of CV
        if CV(x,y) ~= -ERROR_NUM
            neg = 0;
            if CV(x,y) < 0
                neg = 1;
            end
            binSeq = dec2bin(abs(CV(x,y)),8);
            if CV5(x,y) == 2
                binSeq(BITPLANE_NUMBER) = '1';
            else
                binSeq(BITPLANE_NUMBER) = '0';
            end
            num = bin2dec(binSeq);
            if neg == 1
                CV(x,y) = num * -1;
            else
                CV(x,y) = num;
            end
        end

        % restore 5th bit of CD
        if CD(x,y) ~= -ERROR_NUM
            neg = 0;
            if CD(x,y) < 0
                neg = 1;
            end
            binSeq = dec2bin(abs(CD(x,y)),8);
            if CD5(x,y) == 2
                binSeq(BITPLANE_NUMBER) = '1';
            else
                binSeq(BITPLANE_NUMBER) = '0';
            end
            num = bin2dec(binSeq);
            if neg == 1
                CD(x,y) = num * -1;
            else
                CD(x,y) = num;
            end
        end
    end
end

%STEP 6: Take inverse integer wavelet transform to get the original
%distortionless image back
restoredImage = ilwt2(CA,CH,CV,CD,LS);
figure,imshow(restoredImage,[]),title('Restored Image (Distortionless)');

difference = abs(double(originalImage) - double(restoredImage));
figure,imshow(difference,[]),title('Difference Image');


text:
% Fig. 6.12   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num1=[10 10];
den1=[1 10];
w=logspace(-2,3,100);
[m1,p1]=bode(num1,den1,w);
num2=[10 -10];
[m2,p2]=bode(num2,den1,w);
figure(1)
loglog(w,m1,w,m2);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.12 Bode Plot for a NMP System (a) magnitude');
grid;
pause;
figure(2)
semilogx(w,p1,w,p2);
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('Fig 6.12 (b) phase');


text:
function calpak_test535 ( )

%*****************************************************************************80
%
%% CALPAK_TEST535 tests YEAR_IS_LEAP_GREEK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST535\n' );
  fprintf ( 1, '  For the Greek calendar:\n' );
  fprintf ( 1, '  YEAR_IS_LEAP_GREEK reports leap years.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Leap?\n' );
  fprintf ( 1, '\n' );

  for y = 1 : 10
    sy = y_to_s_greek ( y );
    fprintf ( 1, '  %10s  %d\n', sy, year_is_leap_greek ( y ) );
  end

  return
end


text:
%  Make NIfTI structure specified by an N-D matrix. Usually, N is 3 for
%  3D matrix [x y z], or 4 for 4D matrix with time series [x y z t].
%  Optional parameters can also be included, such as: voxel_size,
%  origin, datatype, and description.
%
%  Once the NIfTI structure is made, it can be saved into NIfTI file
%  using "save_nii" command (for more detail, type: help save_nii).
%
%  Usage: nii = make_nii(img, [voxel_size], [origin], [datatype], [description])
%
%  Where:
%
%	img:		Usually, img is a 3D matrix [x y z], or a 4D
%			matrix with time series [x y z t]. However,
%			NIfTI allows a maximum of 7D matrix. When the
%			image is in RGB format, make sure that the size
%			of 4th dimension is always 3 (i.e. [R G B]). In
%			that case, make sure that you must specify RGB
%			datatype, which is either 128 or 511.
%
%	voxel_size (optional):	Voxel size in millimeter for each
%				dimension. Default is [1 1 1].
%
%	origin (optional):	The AC origin. Default is [0 0 0].
%
%	datatype (optional):	Storage data type:
%		2 - uint8,  4 - int16,  8 - int32,  16 - float32,
%		32 - complex64,  64 - float64,  128 - RGB24,
%		256 - int8,  511 - RGB96,  512 - uint16,
%		768 - uint32,  1792 - complex128
%			Default will use the data type of 'img' matrix
%			For RGB image, you must specify it to either 128
%			or 511.
%
%	description (optional):	Description of data. Default is ''.
%
%  e.g.:
%     origin = [33 44 13]; datatype = 64;
%     nii = make_nii(img, [], origin, datatype);    % default voxel_size
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function nii = make_nii(varargin)

   nii.img = varargin{1};
   dims = size(nii.img);
   dims = [length(dims) dims ones(1,8)];
   dims = dims(1:8);

   voxel_size = [0 ones(1,7)];
   origin = zeros(1,5);
   descrip = '';

   switch class(nii.img)
      case 'uint8'
         datatype = 2;
      case 'int16'
         datatype = 4;
      case 'int32'
         datatype = 8;
      case 'single'
         if isreal(nii.img)
            datatype = 16;
         else
            datatype = 32;
         end
      case 'double'
         if isreal(nii.img)
            datatype = 64;
         else
            datatype = 1792;
         end
      case 'int8'
         datatype = 256;
      case 'uint16'
         datatype = 512;
      case 'uint32'
         datatype = 768;
      otherwise
         error('Datatype is not supported by make_nii.');
   end

   if nargin > 1 & ~isempty(varargin{2})
      voxel_size(2:4) = double(varargin{2});
   end

   if nargin > 2 & ~isempty(varargin{3})
      origin(1:3) = double(varargin{3});
   end

   if nargin > 3 & ~isempty(varargin{4})
      datatype = double(varargin{4});

      if datatype == 128 | datatype == 511
         dims(5) = [];
         dims(1) = dims(1) - 1;
         dims = [dims 1];
      end
   end

   if nargin > 4 & ~isempty(varargin{5})
      descrip = varargin{5};
   end

   if ndims(nii.img) > 7
      error('NIfTI only allows a maximum of 7 Dimension matrix.');
   end

   maxval = round(double(max(nii.img(:))));
   minval = round(double(min(nii.img(:))));

   nii.hdr = make_header(dims, voxel_size, origin, datatype, ...
	descrip, maxval, minval);

   switch nii.hdr.dime.datatype
   case 2
      nii.img = uint8(nii.img);
   case 4
      nii.img = int16(nii.img);
   case 8
      nii.img = int32(nii.img);
   case 16
      nii.img = single(nii.img);
   case 32
      nii.img = single(nii.img);
   case 64
      nii.img = double(nii.img);
   case 128
      nii.img = uint8(nii.img);
   case 256
      nii.img = int8(nii.img);
   case 511
      img = double(nii.img(:));
      img = single((img - min(img))/(max(img) - min(img)));
      nii.img = reshape(img, size(nii.img));
      nii.hdr.dime.glmax = double(max(img));
      nii.hdr.dime.glmin = double(min(img));
   case 512
      nii.img = uint16(nii.img);
   case 768
      nii.img = uint32(nii.img);
   case 1792
      nii.img = double(nii.img);
   otherwise
      error('Datatype is not supported by make_nii.');
   end

   return;					% make_nii

%---------------------------------------------------------------------
function hdr = make_header(dims, voxel_size, origin, datatype, ...
	descrip, maxval, minval)

   hdr.hk   = header_key;
   hdr.dime = image_dimension(dims, voxel_size, datatype, maxval, minval);
   hdr.hist = data_history(origin, descrip);

   return;					% make_header

%---------------------------------------------------------------------
function hk = header_key

    hk.sizeof_hdr       = 348;			% must be 348!
    hk.data_type        = '';
    hk.db_name          = '';
    hk.extents          = 0;
    hk.session_error    = 0;
    hk.regular          = 'r';
    hk.dim_info         = 0;

    return;					% header_key

%---------------------------------------------------------------------
function dime = image_dimension(dims, voxel_size, datatype, maxval, minval)

   dime.dim = dims;
   dime.intent_p1 = 0;
   dime.intent_p2 = 0;
   dime.intent_p3 = 0;
   dime.intent_code = 0;
   dime.datatype = datatype;

   switch dime.datatype
   case 2,
      dime.bitpix = 8;  precision = 'uint8';
   case 4,
      dime.bitpix = 16; precision = 'int16';
   case 8,
      dime.bitpix = 32; precision = 'int32';
   case 16,
      dime.bitpix = 32; precision = 'float32';
   case 32,
      dime.bitpix = 64; precision = 'float32';
   case 64,
      dime.bitpix = 64; precision = 'float64';
   case 128
      dime.bitpix = 24;  precision = 'uint8';
   case 256
      dime.bitpix = 8;  precision = 'int8';
   case 511
      dime.bitpix = 96;  precision = 'float32';
   case 512
      dime.bitpix = 16; precision = 'uint16';
   case 768
      dime.bitpix = 32; precision = 'uint32';
   case 1792,
      dime.bitpix = 128; precision = 'float64';
   otherwise
      error('Datatype is not supported by make_nii.');
   end

   dime.slice_start = 0;
   dime.pixdim = voxel_size;
   dime.vox_offset = 0;
   dime.scl_slope = 0;
   dime.scl_inter = 0;
   dime.slice_end = 0;
   dime.slice_code = 0;
   dime.xyzt_units = 0;
   dime.cal_max = 0;
   dime.cal_min = 0;
   dime.slice_duration = 0;
   dime.toffset = 0;
   dime.glmax = maxval;
   dime.glmin = minval;

   return;					% image_dimension

%---------------------------------------------------------------------
function hist = data_history(origin, descrip)

   hist.descrip = descrip;
   hist.aux_file = 'none';
   hist.qform_code = 0;
   hist.sform_code = 0;
   hist.quatern_b = 0;
   hist.quatern_c = 0;
   hist.quatern_d = 0;
   hist.qoffset_x = 0;
   hist.qoffset_y = 0;
   hist.qoffset_z = 0;
   hist.srow_x = zeros(1,4);
   hist.srow_y = zeros(1,4);
   hist.srow_z = zeros(1,4);
   hist.intent_name = '';
   hist.magic = '';
   hist.originator = origin;

   return;					% data_history

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


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

%% Syntax
% |[FT,VT,CT,CT_c]=foamWrap(F,V,C,cPar);|

%% Description
% Use |foamWrap| to generate a foam like structure on top of an input mesh

%% Examples

clear; close all; clc;

%%
% Plot Settings

fontSize=15;
faceAlpha=1;
edgeColor=0.1*ones(1,3);
edgeWidth=1;
cmap=gjet(250);

%%
% Create surface model

[F,V,~]=geoSphere(2,1); %Geodesic sphere
% [F,V]=parasaurolophus;
% [F,V]=cow;
% [F,V]=stanford_bunny;

[F,V,C,indIni]=triPolyDualRefine(F,V);

%%
cFigure;
gpatch(F,V,C);
colormap(cmap);
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%%

cPar.n=3;
cPar.dirFlip=1;
cPar.foamThickness=0.05; %Empty uses default which is mean edgelength based
cParSmooth.Method='HC';
cParSmooth.n=25;
cPar.cParSmooth=cParSmooth;

%%
L_remove=true(size(F,1),1);
[FT,VT,CT,CT_c]=foamWrap(F,V,C,cPar);

%%
cFigure; hold on;
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);

gpatch(FT,VT,CT_c,'none',1);

axisGeom(gca,fontSize);
view(0,58.25);
camlight headlight;
axis off;
colormap(gray(4)); icolorbar;
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:
% POP_TOPOCHANSEL - pop up a topographic interface to select channels
%
% Pops up a topographic interface to select multiple channels with the
% mouse. Click a polygon around the electrodes you wish to select. Right
% click to finish.
%
% Usage:
%   >> [chanlist] = pop_topochansel(chanlocs,selection);
%
% Inputs:
%   chanstruct     - channel structure. See READLOCS
%                    (optional if EEG structure with chanlocs in caller or
%                    base workspace)
%   selection      - currently selected electrodes.
%                    (optional)
%
%
% Output:
%   chanlist      - indices of selected channels
%   cellchannames - names of selected channel names in a cell array
%   strchannames  - names of selected channel names in a concatenated string
%                   (channel names are separated by space characters)
%

function [chanlist, cellchannames, strchannames] = pop_topochansel(chanlocs,select,varargin)
if nargin == 0 || isempty(chanlocs)
    try
        chanlocs = evalin('caller','chanlocs');
    catch
        try
            tmp = evalin('caller','EEG');
            chanlocs = tmp.chanlocs;
        catch
            error('tell me where the electrodes are'),
        end
    end
end
if nargin >= 2
    if ischar(select)|| iscellstr(select)
        select = chnb(select,{chanlocs.labels});
    end
end
g = finputcheck( varargin, ...
                 { 'labels'    'string'    {'on' 'off'}   'off'
                  'cellstrout' 'string'    {'on' 'off'}   'off'
                    },'ignore');

h0 = figure(23537);clf;
set(gcf,'numbertitle','off','name','Channel selection');
topoplot([],chanlocs,'electrodes','on', 'style','blank');

chi = get(gca,'children');
% chi = chi(1:numel(chanlocs));
todel = [];
for i = 1:numel(chi)
    try
        if numel(get(chi(i),'XData')) == sum(~emptycells({chanlocs.X}))
            continue % we found electrodes
        else
            todel(end+1) = i;
        end
    catch
        todel(end+1) = i;
    end
end
chi(todel) = [];
hold on
X = get(chi,'XData');
Y = get(chi,'YData');
delete(chi)
if strcmp(g.labels,'on')
    for i = 1:numel(X)
        text(X(i),Y(i),chanlocs(i).labels);
    end
end

if exist('select','var')
    plot(X(select),Y(select),'r.','markersize',20);
else
    select = [];
end
[dum dum chanlist] = lasso(X,Y);
if isempty(chanlist)
    disp('No new elecs selected. Keeping old selection')
    chanlist = select;
end
strchannames = '';
for i = 1:numel(chanlist)
    strchannames = [strchannames chanlocs(chanlist(i)).labels];
    if i ~= numel(chanlist)
        strchannames = [strchannames ' '];
    end
    cellchannames{i} = chanlocs(chanlist(i)).labels;
end
if strcmp(g.cellstrout,'on')
    chanlist = cellchannames;
end
if nargout == 0
    disp(strchannames)
    disp(chanlist)
    clear
end
% close(h0);

function [selx,sely,indexnr]=lasso(x,y)

% lasso -  enables the selection/encircling of (clusters of) events in a scatter plot by hand
%          using the mouse
%
% Input:    x,y                 - a set of points in 2 column vectors.
% Output:   selx,sely,indexnr   - a set of selected points in 3 column vectors
%
% Note:   After the scatter plot is given, selection by mouse is started after any key press.
%         This is done to be able to ZOOM or CHANGE AXES etc. in the representation before selection
%         by mouse.
%         Encircling is done by pressing subsequently the LEFT button mouse at the requested positions
%         in a scatter plot.
%         Closing the loop is done by a RIGHT button press.
%
% T.Rutten V2.0/9/2003
% downloaded and adapted from mathworks fileexchange by M. Chaumon 2011

plot(x,y,'ob')

las_x=[];
las_y=[];

c=1;

key=0;

while c==1
    try
        [a,b,c]=ginput(1);
        las_x=[las_x;a];las_y=[las_y;b];
        line(las_x,las_y)
    catch
        selx = [];sely = []; indexnr = [];
        return
    end
end

las_x(length(las_x)+1)=las_x(1);
las_y(length(las_y)+1)=las_y(1);

line(las_x,las_y)
pause(.2)

in=inpolygon(x,y,las_x,las_y);

ev_in=find(in>0);

selx=x(ev_in);
sely=y(ev_in);
plot(selx,sely,'r.','markersize',20);
drawnow

indexnr=ev_in;

function nb = chnb(channame, varargin)

% CHNB - return channel number corresponding to channel names in an EEG
%           structure
%
% Usage:
%   >> [nb] = chnb(channame);
%   >> [nb] = chnb(channame, labels);
%
% Input:
%   channame      - name of channels to search. Either a string with space
%                   separated channel names, or a cell array of strings.
%                   Note that regular expressions can be used to match
%                   several channels. See regexp.
%   labels        - channel labels as found in EEG.chanlocs.labels.
%
% Output:
%   nb            - channel numbers in the EEG structure found in the
%                   caller workspace (i.e. where the function is called
%                   from) or in the base workspace, if no EEG structure
%                   exists in the caller workspace.
%
error(nargchk(1,2,nargin));
if nargin == 2
    labels = varargin{1};
else

    try
        EEG = evalin('caller','EEG');
    catch
        try
            EEG = evalin('base','EEG');
        catch
            error('Could not find EEG structure');
        end
    end
    if not(isfield(EEG,'chanlocs'))
        error('No channel list found');
    end
    labels = {EEG.chanlocs.labels};
end
if ischar(channame)
    tmp = regexp(channame,'(\S*) ?','tokens');
    channame = {};
    for i = 1:numel(tmp)
        channame{i} = tmp{i}{1};
    end
    if isempty(channame)
        nb = [];
        return
    end
end

nb = regexpcell(labels,channame,'exactignorecase');

function idx = regexpcell(c,pat, cmds)

% idx = regexpcell(c,pat, cmds)
%
% Return indices idx of cells in c that match pattern(s) pat (regular expression).
% Pattern pat can be char or cellstr. In the later case regexpcell returns
% indexes of cells that match any pattern in pat.
%
% cmds is a string that can contain one or several of these commands:
% 'inv' return indexes that do not match the pattern.
% 'ignorecase' will use regexpi instead of regexp
% 'exact' performs an exact match (regular expression should match the whole strings in c).
% 'all' (default) returns all indices, including repeats (if several pat match a single cell in c).
% 'unique' will return unique sorted indices.
% 'intersect' will return only indices in c that match ALL the patterns in pat.
%
% v1 Maximilien Chaumon 01/05/09
% v1.1 Maximilien Chaumon 24/05/09 - added ignorecase
% v2 Maximilien Chaumon 02/03/2010 changed input method.
%       inv,ignorecase,exact,combine are replaced by cmds

error(nargchk(2,3,nargin))
if not(iscellstr(c))
    error('input c must be a cell array of strings');
end
if nargin == 2
    cmds = '';
end
if not(isempty(regexpi(cmds,'inv', 'once' )))
    inv = true;
else
    inv = false;
end
if not(isempty(regexpi(cmds,'ignorecase', 'once' )))
    ignorecase = true;
else
    ignorecase = false;
end
if not(isempty(regexpi(cmds,'exact', 'once' )))
    exact = true;
else
    exact = false;
end
if not(isempty(regexpi(cmds,'unique', 'once' )))
    combine = 2;
elseif not(isempty(regexpi(cmds,'intersect', 'once' )))
    combine = 3;
else
    combine = 1;
end

if ischar(pat)
    pat = cellstr(pat);
end

if exact
    for i_pat = 1:numel(pat)
        pat{i_pat} = ['^' pat{i_pat} '$'];
    end
end

for i_pat = 1:length(pat)
    if ignorecase
        trouv = regexpi(c,pat{i_pat}); % apply regexp on each pattern
    else
        trouv = regexp(c,pat{i_pat}); % apply regexp on each pattern
    end
    idx{i_pat} = [];
    for i = 1:numel(trouv)
        if not(isempty(trouv{i}))% if there is a match, store index
            idx{i_pat}(end+1) = i;
        end
    end
end
switch combine
    case 1
        idx = [idx{:}];
    case 2
        idx = unique([idx{:}]);
    case 3
        for i_pat = 2:length(pat)
            idx{1} = intersect(idx{1},idx{i_pat});
        end
        idx = idx{1};
end
if inv % if we want to invert result, then do so.
    others = 1:numel(trouv);
    others(idx) = [];
    idx = others;
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [tip1, tip2] = scimat_rv_crescent_tips(scimat, m)
% SCIMAT_RV_CRESCENT_TIPS  Extract the tips of the crescent-shaped curve
% in all slices of the Right Ventricle.
%
% [X1, X2] = scimat_rv_crescent_tips(SCIMAT, M)
%
%   X1, X2 are 3-colum matrices where each row are the real world
%   coordinates of each of the tips of the RV crescent shape.
%
%   SCIMAT is the struct with the RV segmentation mask (see "help scimat"
%   for details).
%
%   M is a 3-column matrix, where each row has the real world coordinates
%   of a centroid (typically, the central Left Ventricle's central curve
%   computed with scimat_centroids()). There must be as many centroids as
%   slices. Centroids with NaN coordinates will be skipped.
%
%   In the case of a RV slice that has no corresponding LV centroid, the
%   closest LV centroid will be used.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010,2014 University of Oxford
% Version: 0.2.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 2);
nargoutchk(0, 2);

% volume size
sz = [scimat.axis.size];

% if (size(m, 1) ~= sz(3))
%     error('There must be a centroid per slice in SCIMAT, even if the centroid in NaN')
% end
if (size(m, 2) ~= 3)
    error('M must be a 3-column matrix')
end

% init outputs
tip1 = nan(sz(3), 3);
tip2 = nan(sz(3), 3);

% iterate slices
for I = 1:sz(3)
    % extract slice
    im = scimat.data(:,:, I);

    % get linear indices of all pixels in current slice
    idx = find(im);

    % if RV slice is empty, skip it
    if isempty(idx)
        continue
    end

    % convert linear index to multiple subscripts
    [ir, ic] = ind2sub( sz(1:2), idx );

    % convert indices to real world coordinates and make colum vectors
    x = scimat_index2world([ir, ic, I+zeros(length(ir), 1) ], ...
        scimat);

    % compute a centroid for the slice
    xm = mean(x, 1);

    % compute the distance from the slice centroid to each axis centroid
    d = dmatrix(xm', m');

    % find closest axis centroid to the slice centroid
    [~, inn] = min(d);

    % center the slice points around the axis centroid
    x = x - m(inn * ones(size(x, 1), 1), :);

    % compute polar coordinates of the centered pixel coordinates
    [phi, th] = cart2pol(x(:,1), x(:,2));

    % assume that the point with the largest and smallest azimuth values
    % are the crescent tips
    [~, idx1] = max(phi);
    [~, idx2] = min(phi);

    % assign the tip points to the output variable, undoing the previous
    % centering
    tip1(I, :) = x(idx1, :) + m(inn, :);
    tip2(I, :) = x(idx2, :) + m(inn, :);

end


text:
function [inpts,volpts]=mrUndoAPoint(inpts,volpts)
%function [inpts,volpts]=mrUndoAPoint(inpts,volpts)
%
%Removes the last selected pair of alignment points from the list

%6/16/96	gmb	wrote it.

npoints = size(volpts,1);

inpts=inpts(1:npoints-1,:);
volpts=volpts(1:npoints-1,:);

npoints = size(volpts,1);
disp([num2str(npoints),' pairs of points selected.']);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function R = sw_ht_converse(a, n_vals, eps)
% Compute the converse sum rate (at the symmetrical rate point) from Han's point-to-point converse
% for the joint source distribution: [a 1/3*(1-a); 1/3*(1-a) 1/3*(1-a)].

% Input:
% a: parameter to define the joint source distribution, should be in range [0.25, 1)
% n_vals: a vector containing the blocklengths
% eps: target error probability
% Output:
% R: lower bounds of the rates for the blocklengths specified in n_vals
% Other:
% One can also choose to return
% R_sw_ht_hi: upper bounds of the rates for the blocklengths specified in n_vals
% R_sw_ht_lo: lower bounds of the rates for the blocklengths specified in n_vals

    % Compute the optimal log_c coefficients from binary hypothesis testing,
    % which serve as the start point for computing the parameters for
    % composite hypothesis testing.
    [~, log_c_binary] = p2p_ht_converse(a, n_vals, eps);

    % Compute log binomial coefficients
    log_b = binomial_coeff(max(n_vals));

    R_sw_ht_hi = zeros(n_vals,1);
    R_sw_ht_lo = zeros(n_vals,1);

    for ni = 1:length(n_vals)
        n = n_vals(ni);
        m = 0:n;
        log_Pr = m*log(a)+(n-m)*log(1/3*(1-a));  % log_Pr stores all the log joint probabilities in ascending order
        log_Pm = m*log(a+1/3*(1-a))+(n-m)*log(1/3*(1-a)+1/3*(1-a));  % log_Pm stores all the log marginal probabilities in ascending order

        % Initial points
        x = log_c_binary(n) - 0.1;
        y = log_c_binary(n) - 0.0001;

        [beta1, beta2] = my_bisect(n,eps,x,log_Pr,log_Pm,log_b);
        bx = [beta1 beta2];
        diffx = beta1 - beta2 > 0;
        [beta1, beta2] = my_bisect(n,eps,y,log_Pr,log_Pm,log_b);
        by = [beta1 beta2];
        diffy = beta1 - beta2 > 0;

        % For catching special cases: boolean 1 -> zero exists, 0 -> no zero
        disp(strcat('n =',32,num2str(n),',',32,num2str(diffx ~= diffy)));

        if diffx ~= diffy
            % Do bisection search if zero exists
            iter = 0;
            while abs(x-y) > 0.0001 && iter < 20
                c = (x+y)/2;
                [beta1, beta2] = my_bisect(n,eps,c,log_Pr,log_Pm,log_b);
                if (beta1 - beta2 > 0) == diffx
                   x = c;
                   bx = [beta1 beta2];
                else
                   y = c;
                   by = [beta1 beta2];
                end
                iter = iter + 1;
            end
            [M,I] = min([max(bx) max(by)]);  % Minimum upper bound
            R_sw_ht_hi(ni) = M;
            lo = [min(bx) min(by)];
            R_sw_ht_lo(ni) = lo(I);          % A lower bound
        else
            % Scan a range of possible parameters for the best bound
            if (n < 100)
                c_val = log_c_binary(n)-0.1:0.0001:log_c_binary(n)-0.0001;
            else
                c_val = log_c_binary(n)-0.005:0.0001:log_c_binary(n)-0.0001;
            end
            beta1_val = zeros(length(c_val),1);
            beta2_val = zeros(length(c_val),1);
            for ci = 1:length(c_val)
                c = c_val(ci);
                [beta1, beta2] = my_bisect(n,eps,c,log_Pr,log_Pm,log_b);
                beta1_val(ci) = beta1;
                beta2_val(ci) = beta2;
            end
            R_sw_ht_hi(ni) = min(max(beta1_val,beta2_val));  % Minimum upper bound
            R_sw_ht_lo(ni) = max(min(beta1_val,beta2_val));  % Maximum lower bound
        end
    end
    R = R_sw_ht_lo;
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function M = fixedrankMNquotientfactory(m, n, k)
% Manifold of m-by-n matrices of rank k with quotient geometry.
%
% function M = fixedrankMNquotientfactory(m, n, k)
%
% This follows the quotient geometry described in the following paper:
% P.-A. Absil, L. Amodei and G. Meyer,
% "Two Newton methods on the manifold of fixed-rank matrices endowed
%  with Riemannian quotient geometries", arXiv, 2012.
%
% Paper link: http://arxiv.org/abs/1209.0068
%
% A point X on the manifold is represented as a structure with two
% fields: M and N. The matrix M (mxk) is orthonormal, while the matrix N
% (nxk) is full-rank.
%
% Tangent vectors are represented as a structure with two fields (M, N).

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors:
% Change log:

M.name = @() sprintf('MN'' quotient manifold of %dx%d matrices of rank %d', m, n, k);

M.dim = @() (m+n-k)*k;

% Choice of the metric is motivated by the symmetry present in the
% space
M.inner = @(X, eta, zeta) eta.M(:).'*zeta.M(:) + eta.N(:).'*zeta.N(:);

M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

M.dist = @(x, y) error('fixedrankMNquotientfactory.dist not implemented yet.');

M.typicaldist = @() 10*k;

symm = @(X) .5*(X+X');
stiefel_proj = @(M, H) H - M*symm(M'*H);

M.egrad2rgrad = @egrad2rgrad;
    function eta = egrad2rgrad(X, eta)
        eta.M = stiefel_proj(X.M, eta.M);
    end

M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)

        % Directional derivative of the Riemannian gradient
        Hess.M = ehess.M - eta.M*symm(X.M'*egrad.M);
        Hess.M = stiefel_proj(X.M, Hess.M);

        Hess.N = ehess.N;

        % Projection onto the horizontal space
        Hess = M.proj(X, Hess);
    end

M.proj = @projection;
    function etaproj = projection(X, eta)

        % Start by projecting the vector from Rmp x Rnp to the tangent
        % space to the total space, that is, eta.M should be in the
        % tangent space to Stiefel at X.M and eta.N is arbitrary.
        eta.M = stiefel_proj(X.M, eta.M);

        % Now project from the tangent space to the horizontal space, that
        % is, take care of the quotient.

        % First solve a Sylvester equation (A symm., B skew-symm.)
        A = X.N'*X.N + eye(k);
        B = eta.M'*X.M + eta.N'*X.N;
        B = B-B';
        omega = lyap(A, -B);

        % And project along the vertical space to the horizontal space.
        etaproj.M = eta.M + X.M*omega;
        etaproj.N = eta.N + X.N*omega;

    end

M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        A = t*X.M'*eta.M;
        S = t^2*eta.M'*eta.M;
        Y.M = [X.M t*eta.M]*expm([A -S ; eye(k) A])*eye(2*k, k)*expm(-A);

        % re-orthonormalize (seems necessary from time to time)
        [Q R] = qr(Y.M, 0);
        Y.M = Q * diag(sign(diag(R)));

        Y.N = X.N + t*eta.N;

    end

% Factor M lives on the Stiefel manifold, hence we will reuse its
% random generator.
stiefelm = stiefelfactory(m, k);

M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        Y.M = uf(X.M + t*eta.M); % This is a valid retraction
        Y.N = X.N + t*eta.N;
    end

M.hash = @(X) ['z' hashmd5([X.M(:) ; X.N(:)])];

M.rand = @random;
    function X = random()
        X.M = stiefelm.rand();
        X.N = randn(n, k);
    end

M.randvec = @randomvec;
    function eta = randomvec(X)
        eta.M = randn(m, k);
        eta.N = randn(n, k);
        eta = projection(X, eta);
        nrm = M.norm(X, eta);
        eta.M = eta.M / nrm;
        eta.N = eta.N / nrm;
    end

M.lincomb = @lincomb;

M.zerovec = @(X) struct('M', zeros(m, k), 'N', zeros(n, k));

M.transp = @(x1, x2, d) projection(x2, d);

end

% Linear combination of tangent vectors
function d = lincomb(x, a1, d1, a2, d2) %#ok<INMSL>

if nargin == 3
    d.M = a1*d1.M;
    d.N = a1*d1.N;
elseif nargin == 5
    d.M = a1*d1.M + a2*d2.M;
    d.N = a1*d1.N + a2*d2.N;
else
    error('Bad use of fixedrankMNquotientfactory.lincomb.');
end

end

function A = uf(A)
[L, unused, R] = svd(A, 0);
A = L*R';
end

author: andrewssobral

repo: lrslibrary


text:
function [collectedEyeData, img, fixMap] = showAllEyeDataForImage2(filename, FixationDB, STIMULI, HASH)

if strcmp(filename(end-2:end), 'jpg')
    stimname = filename(1:end-4);
elseif strcmp(filename(end-2:end), 'mat')
    stimname = filename(1:end-4);
else
    disp('filename not .mat or .jpg');
end

if exist(fullfile(STIMULI, strcat(stimname, '.jpg')))
    img = imread(fullfile(STIMULI, strcat(stimname, '.jpg')));
else
    fprintf('image file not found\n')
    fullfile(STIMULI, strcat(stimname, '.jpg'))
    collectedEyeData = [];
    img = NaN;
    fixMap = NaN;
    return
end

% get the eyetracking data from the FixationDB
load(HASH);
index = hash.(stimname);
collectedEyeData = FixationDB(index).fixationData;
[users{1:length(collectedEyeData)}] = deal(collectedEyeData.user);
fixMap = zeros(size(img, 1), size(img, 2));

for u=1:length(users)

    appropFix = collectedEyeData(u).appropFix;

    % create new image with fixations and save
    for j=1:size(appropFix, 1)
        pt = appropFix(j, :);
        fixMap(pt(2), pt(1))=1;
    end

end

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function [Keff_p_medio, Keff_s_medio, Keff_n_medio] = interpolateElectrolyteConductivities(Keff_p,Keff_s,Keff_n,param)
%	interpolateElectrolyteConductivities interpolates electrolyte conductivities at the edges of control volumes using harmonic mean.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

%% Positive electrode mean conductivity
beta_p = 0.5;
Keff_p_medio = Keff_p(1:end-1).*Keff_p(2:end)./(beta_p*Keff_p(2:end)+(1-beta_p)*Keff_p(1:end-1));

% The last element of Keff_p_medio will be the harmonic mean of the
% elements at the interface positive-separator

beta_p_s = param.deltax_p*param.len_p/2 /(param.deltax_s*param.len_s/2+param.deltax_p*param.len_p/2);

Keff_p_s_interface = Keff_p(end)*Keff_s(1) / (beta_p_s*Keff_s(1) + (1-beta_p_s)*Keff_p(end));

Keff_p_medio = [Keff_p_medio;Keff_p_s_interface];

%% Separator mean conductivity
% Compute the harmonic mean values for the separator
beta_s = 0.5;
Keff_s_medio = Keff_s(1:end-1).*Keff_s(2:end)./(beta_s*Keff_s(2:end)+(1-beta_s)*Keff_s(1:end-1));

% The last element of Keff_s_medio will be the harmonic mean of the
% elements at the interface separator-negative

beta_s_n = param.deltax_s*param.len_s/2 /(param.deltax_s*param.len_s/2+param.deltax_n*param.len_n/2);

Keff_s_n_interface = Keff_s(end)*Keff_n(1) / (beta_s_n*Keff_n(1) + (1-beta_s_n)*Keff_s(end));

Keff_s_medio = [Keff_s_medio;Keff_s_n_interface];

%% Negative electrode mean conductivity
% Compute the harmonic mean values for the negative electrode

beta_n = 0.5;
Keff_n_medio = Keff_n(1:end-1).*Keff_n(2:end)./(beta_n*Keff_n(2:end)+(1-beta_n)*Keff_n(1:end-1));

Keff_n_medio = [Keff_n_medio;0]; % the cc interface is not used. The zero is only to match the dimensions
end

author: lionsimbatoolbox

repo: LIONSIMBA