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

function pathS = MC_NIGGOU(S0,r,d,T,alpha,beta,delta,lambda,a,b,NTime,NSim,NBatches)
% discretization of Normal Inverse Gaussian process with
% Gamma Ornstein-Uhlenbeck clock

intNt = 10;                 % steps in between orginial grid points
allsteps = intNt * NTime;   % All Nts that have to be simulated
dT = T / NTime;        % Delta for time discretization
time = 0 : dT : T;    % dummy for martingale correction

pathS = zeros(NSim,NTime+1,NBatches); % output
lnS = zeros(NSim,NTime+1);            % used in batch
lnS(:,1) = log(S0*exp(-d*T));         % set S(0) dividend adjusted

% precompute constants
y0 = 1;
psiNIG = (-1i)*(-delta)*(sqrt(alpha^2-(beta+1)^2) ...
    -sqrt(alpha^2-beta^2));                       % char exp
phiGOU = 1i*psiNIG*y0/lambda*(1-exp(-lambda*time)) ...
         + lambda*a./(1i*psiNIG-lambda*b) ...
         .*(b*log(b./(b-1i*psiNIG/lambda ...
         *(1-exp(-lambda*time))))-1i*psiNIG*time);% char func
omegaT = -phiGOU;           % martingale correction
omegaT(1) = 0;              % martingale correction in 0

for l = 1 : NBatches        % batch loop
    yy = ones(NSim,allsteps+1);             % spot clock
    Np = poissrnd(a*lambda/allsteps*T,[NSim,allsteps]);   % Poissonians

    for k = 1 : NSim
        for j = 1 : allsteps                  % generating OU process
            if Np(k,j) > 0
                Ex = -log(rand(Np(k,j),1))/b; % RVs with exponential law
                U = exp(-lambda * T / allsteps * rand(Np(k,j),1));          % Uniforms
                yy(k,j+1) = (1-lambda*T/allsteps)*yy(k,j) + sum(Ex .* U);
         else
                yy(k,j+1) = (1-lambda*T/allsteps)*yy(k,j);
            end
        end
    end

    ZZ = T*cumsum(yy,2)/allsteps;       % integrated Gamma O-U process
    Y = zeros(NSim,NTime+1);            % intergrated time
    % compute the integrated time at discretization steps
    for m=2:NTime+1
        Y(:,m) = ZZ(:,(m-1)*intNt);
    end

    for m=2:NTime+1     % time loop
        a_par = Y(:,m)-Y(:,m-1);             % IG param
        b_par = delta*sqrt(alpha^2-beta^2);  % IG param
        theta = a_par/b_par;
        chi = a_par.^2;
        %Yvec = ones(NSim,1);
        %for n = 1:NSim
        %    Yvec(n) = randraw('ig', [theta(n), chi(n)], 1 );
        %end
        chisq1 = randn(NSim,1).^2;
        Yvec = theta + 0.5*theta./chi .* ( theta.*chisq1 - ...
            sqrt(4*theta.*chi.*chisq1 + theta.^2.*chisq1.^2) );
        Ind = find(rand(NSim,1) >= theta./(theta+Yvec));
        Yvec(Ind) = theta(Ind).^2./Yvec(Ind);   % subordinator

        Zvec = randn(NSim,1);       % Gaussian
        lnS(:,m) = lnS(:,m-1) + (r-d)*dT + omegaT(m)-omegaT(m-1) ...
                     + beta*delta^2*Yvec + delta*sqrt(Yvec).*Zvec;
    end
    pathS(:,:,l) = exp(lnS);   % spot paths
end


text:
% demonstrate use of hhist

figure;
x = sample([0.1 0.2 0.3 0.4], 1000);
hhist(x);

figure;
% verify density for product of indep normals
x = prod(randn(2,100000));
ts = linspace(-1,1,1000);
z = hhist(x, ts);
h = plot(ts, z, ts, besselk(0,abs(ts))*1/pi);
%set(h(2),'LineWidth',2)
axis_pct

author: andrewssobral

repo: lrslibrary


text:
function f = diff( f, varargin )
%DIFF   Derivative of a diskfun in Cartesian coordinates.
%
%   F = DIFF( F ) computes the first derivative of F with respect to x.
%
%   F = DIFF( F, DIM )  computes the first derivative of F. If DIM = 1, the
%   derivative is taken in the x-direction. If DIM = 2, the derivative
%   is taken in the y-direction.
%
%   F = DIFF( F, DIM, K) computes the kth derivatives of F in the variable
%   given by DIM.
%
% See also DISKFUN/LAPLACIAN, DISKFUN/DIFFX, DISKFUN/DIFFY

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( isempty(f) )
    f = diskfun;
    return
end

% Parse user inputs:
if ( nargin == 1 )
    dim = 1;
    K = 1;
elseif ( nargin == 2 )
    dim = varargin{1};
    K = 1;
else
    dim = varargin{1};
    K = varargin{2};
end

if ( dim ~= 1 && dim ~= 2  )
    error('DISKFUN:DIFF:DIM', 'Unrecognized coordinate dimension.');
end

if ( abs( K - round(K) ) > eps )
    error('DISKFUN:DIFF:DIFFORDER', 'Fractional derivatives not allowed.')
end
K = round( K );

% Implement higher derivatives as repeated (iterated) differentiation
for j=1:K
    f = onediff(f, dim);
end

end
% TODO: This code will not work for complex valued diskfuns, if we ever
% allow them.

function f = onediff(f, dim)

% Make sure f has as short a length as possible.
f=simplify(f);

% We are going to work at the tech level to make things faster.
[C, ~, R] = cdr( f );

% Do everything with even length columns since then no special
% modifications are required for dividing by rho series expansion.
parity = mod(length(C),2);
n = length(C)+parity;
m = length(R);

% The variable coefficients in the definitions of the derivatives means
% that the length of the columns and rows will increase by one wave number
% after taking the derivatives with respect to x and y.
m = m+2;
n = n+2;

% Matrices for multiplying by sin(theta), cos(theta) and 1/r in coefficient space.
Msinm = .5i*spdiags(ones(m,1)*[-1,1],[-1 1],m,m);
Mcosm = .5*spdiags(ones(m,1)*[1,1],[-1 1],m,m);
Mn = ultraS.multmat(n, [0;1], 0);

% Work at the tech level to make things faster.
ctechs = C.funs{1}.onefun;
rtechs = R.funs{1}.onefun;

% Alias to get the extra zeros in place.
ctechs.coeffs = ctechs.alias(ctechs.coeffs, n);
rtechs.coeffs = rtechs.alias(rtechs.coeffs, m);

% Compute the derivatives.
dCdr = diff(ctechs);
dRdth = diff(rtechs)/pi;

% 1/r*col
rinv = Mn \ ctechs.coeffs;

if (dim==1) %d/dx

    %CDR for -1/r.*sin(th).*d/dth
    C1 = rinv;
    R1 = -Msinm*dRdth.coeffs;

    %CDR for cos(th).*d/dr
    C2 = dCdr.coeffs;
    R2 = Mcosm*rtechs.coeffs;

else  %d/dy

    %CDR for 1/r.*cos(th).*d/dth
    C1 = rinv;
    R1 = Mcosm*dRdth.coeffs;

    %CDR for sin(th).*d/dr
    C2 = dCdr.coeffs;
    R2 = Msinm*rtechs.coeffs;
end

% Put pieces back together
f1 = f;
c1techs = chebtech2({' ',C1});
f1.cols.funs{1}.onefun = c1techs;
r1techs = real(trigtech({'',R1}));
f1.rows.funs{1}.onefun = r1techs;

% Parity changes
temp = f1.idxPlus;
f1.idxPlus = f1.idxMinus;
f1.idxMinus = temp;

f2 = f;
c2techs = chebtech2({'',C2});
f2.cols.funs{1}.onefun = c2techs;
r2techs = real(trigtech({'',R2}));
f2.rows.funs{1}.onefun = r2techs;

% Parity changes
temp = f2.idxPlus;
f2.idxPlus = f2.idxMinus;
f2.idxMinus = temp;

% Compression plus may not preserve the expansion properties we want.
% So we sample each piece add them together and construct a diskfun.
% TODO: Fix this so everything is done in coefficient space, like this
% f = f1 + f2;

% Constructor requires even m
m = m + mod(m, 2);

f = diskfun(sample(f1,m,n/2+1)+sample(f2,m,n/2+1));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function display(a)
% function display(A)
%
% DESCRIPTION
%   Displays a polynomial.
%
% INPUTS
%   A: polynomial
%
% OUTPUTS
%   NONE
%
% SYNTAX
%   display(A);

% 6/7/2002: PJS  Initial Coding
% 6/9/2002: PJS  Use char conversion and display matrices

%   Get String Representation of a
s = char(a);
sza = size(a);

%   Display the Polynomial
if ~isequal(get(0,'FormatSpacing'),'compact')
    disp(' ');
end

% Compute sizes of character arrays and max size
nchar = cellfun('size',s,2);
maxr = max(nchar,[],1);
numc = sum(maxr);
rootprops = fieldnames(get(0));
if any(strcmp(rootprops,'CommandWindowSize'))
    ws = get(0,'CommandWindowSize');
else
    % Older Versions of matlab don't allow access to the
    % window size.  Modify window size (ws) if polynomial
    % line breaks are not in the right spot.
    ws = 80;
end
maxchar=ws(1)-6-3*(sza(2)-1);

if isempty(a)
    disp([inputname(1) ' = ']);
    disp([' ' s{1}]);

elseif max(numc) < maxchar
    % Display as matrix if all rows can fit on screen
    disp([inputname(1) ' = ']);
    for i1 = 1:sza(1);
        if all(sza==[1 1])
            d = '  ';
        else
            d = '  [ ';
        end
        for i2 = 1:sza(2);
            d = [d blanks(maxr(i2)-nchar(i1,i2)) s{i1,i2}];
            if i2~=sza(2)
                d = [d ', '];
            elseif ~all(sza==[1 1])
                d = [d ']'];
            end
        end
        disp(d);
    end

else
    % Else display entry-by-entry going down columns
    for i2 = 1:sza(2);
        if min(sza) > 1
            disp(['----  Column ' int2str(i2) ' ----------'])
            disp(' ')
        end
        for i1 = 1:sza(1);
            % Display the (i1,i2) entry
            if ~all(sza==[1 1]);
                disp([inputname(1) '(' int2str(i1) ',' int2str(i2) ')  = ']);
            else
                disp([inputname(1) ' = ']);
            end

            % Break lines at +,-, or *
            sr = s{i1,i2};
            while ~isempty(sr) %length(sr)>0
                if length(sr) < (ws(1)-6)
                    disp(['  ' sr]);
                    sr = [];
                else
                    idx1 = sort([strfind(sr,'-') strfind(sr,'+') strfind(sr,'*')]);
                    idx1 = [setdiff(idx1,1) length(sr)];
                    %idx2 = max(find(idx1< (ws(1)-6) ));
                    idx2 = find(idx1< (ws(1)-6), 1, 'last' );
                    if isempty(idx2)
                        disp(['  ' sr(1:idx1(1)-1)]);
                        sr = sr(idx1(1):end);
                    else
                        disp(['  ' sr(1:idx1(idx2)-1)]);
                        sr = sr(idx1(idx2):end);
                    end
                end
            end
            if ~all([i1 i2]==sza)
                disp(' ');
            end
        end
    end

end

if ~isequal(get(0,'FormatSpacing'),'compact')
    disp(' ');
end

author: yu-jiang

repo: radpbook

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


text:
function x = RecoverNonlinearSolversolution(model,xout);

if isempty(model.nonlinearindicies)
    x = xout(:);
else
    x = zeros(length(model.c),1);
    for i = 1:length(model.linearindicies)
        x(model.linearindicies(i)) = xout(i);
    end
    x = x(1:length(model.c));
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% example_meshscript:
%
% quick script to collect a bunch of mesh images automagically, and save in
% the 'Images/' folder. This loops across different data types (which in
% this case represent different population receptive field models), and
% loops across maps produced by each pRF model (the variance explained, the
% polar angle, eccentricity, and pRF size), THEN loops across different
% view angles for each hemisphere, taking a snapshot and saving it as an
% PNG image with a stereotyped name. The images should then be easy to sort
% through. The final line calls a separate function (also included as an
% example) which then loads these images, and compsites the results across
% hemispheres into a PowerPoint file.
%
%
% ras, 080408.
rois = {'aFus_faces' 'pFus_faces' 'FBA_event' 'aEBA' ...
		'sEBA'  'TOS_event' 'PPA' ...
		'MT.mat' 'STS_faces_event'};
meshes = {'s3_left_new_football' 's3_right_new_football'};
dtNames = {'CheckerboardAverages' 'HouseAverages' 'FaceAverages' 'Averages' ...
		   'Odd' 'Even' 'FaceOdd' 'FaceEven' 'HouseOdd' 'HouseEven' ...
		   'CheckerboardOdd' 'CheckerboardEven'};
meshViews = { {'VentralZoom' 'Medial' 'LateralOccipital'} ...
			  {'VentralZoom' 'Medial' 'LateralOccipital'} };
fields = {'co' 'map' 'ph' 'amp'};
fieldNames = {'Variance Explained' 'Eccentricity' 'Polar Angle' 'pRF Diameter'};
lr = 'LR';  leftright = {'Left' 'Right'};  rightleft = {'right' 'left'};

%% open the view
mrvCleanWorkspace;
mrVista 3
% VOLUME{1} = prfLoadROIs(VOLUME{1}, rois);
for h = 1:2
	meshPath = fullfile('3DAnatomy', leftright{h}, '3DMeshes', meshes{h});
	VOLUME{1} = meshLoad(VOLUME{1}, meshPath, 1);
end

%% loop across data types, pRF models
for dt = dtNames
	VOLUME{1} = selectDataType(VOLUME{1}, dt{1});

	%% load the model / maps
	% (the 'prfLoadModel' function is a personal function: this calls the
	% 'rmSelect' function to load the most recent pRF model.)
% 	VOLUME{1} = prfLoadModel(VOLUME{1});
	VOLUME{1} = rmSelect(VOLUME{1}, 'mostrecent');

	%% loop across view fields
	for f = 1:length(fields)
		VOLUME{1} = setDisplayMode(VOLUME{1}, fields{f});

		% set viewing parameters appropriate to the field
		if isequal(fields{f}, 'co')
			VOLUME{1} = setCothresh(VOLUME{1}, 0.01); % minimal thresholding
			VOLUME{1}.ui.coMode.clipMode = [0 .6];
		else
			VOLUME{1} = setCothresh(VOLUME{1}, 0.1);
			VOLUME{1}.ui.coMode.clipMode = [0 .5];
		end

		%% loop across meshes (hemispheres)
		for h = 1:2
			if isequal(fields{f}, 'ph')
% 				VOLUME{1}.ui.phMode.cmap(129:end,:) = circshift(hsv(128), 32);
				VOLUME{1} = cmapPolarAngleRGB(VOLUME{1}, rightleft{h});
			else
				VOLUME{1} = setPhWindow(VOLUME{1}, [0 2*pi]);
			end
			VOLUME{1}.meshNum3d = h;
			meshColorOverlay(VOLUME{1});

			%% loop across mesh view settings
			for j = 1:3
				meshRetrieveSettings(VOLUME{1}.mesh{h}, meshViews{h}{j});
				img = mrmGet(VOLUME{1}.mesh{h}, 'screenshot') ./ 255;
				saveName = sprintf('%s %sH %s %s.png',  meshViews{h}{j}, ...
									 lr(h), fieldNames{f}, dt{1});
				savePath = fullfile(pwd, 'Images', saveName);
				imwrite(img, savePath, 'png');
				fprintf('Saved %s.\n', savePath);
			end
		end
	end
end

example_prfMakePPTFromImages;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function setup(varargin)
% SETUP  Add the required search paths to MATLAB
run matconvnet/matlab/vl_setupnn ;
addpath matconvnet/examples ;

opts.useGpu = false ;
opts.verbose = false ;
opts.enableImReadJPEG = false ;
opts = vl_argparse(opts, varargin) ;

try
  vl_nnconv(single(1),single(1),[]) ;
catch
  warning('VL_NNCONV() does not seem to be compiled. Trying to compile it now.') ;
  vl_compilenn('enableGpu', opts.useGpu, ...
               'enableImReadJPEG', opts.enableImReadJPEG, ...
               'verbose', opts.verbose) ;
end

if opts.useGpu
  try
    vl_nnconv(gpuArray(single(1)),gpuArray(single(1)),[]) ;
  catch
    vl_compilenn('enableGpu', opts.useGpu, ...
                 'enableImReadJPEG', opts.enableImReadJPEG, ...
                 'verbose', opts.verbose) ;
    warning('GPU support does not seem to be compiled in MatConvNet. Trying to compile it now') ;
  end
end

author: vedaldi

repo: practical-cnn


text:
function mori = parents(mori)
% variants of an orientation relationship
%
% Syntax
%
%   ori_parents = ori_child * inv(mori.parents)
%
% Input
%  mori - child to parent @orientation relationship
%  ori_child - child orientation
%
% Output
%  ori_parents - all possible parent @orientation
%
% Example
%   % parent symmetry
%   cs_fcc = crystalSymmetry('m-3m', [3.6599 3.6599 3.6599], 'mineral', 'Iron fcc');
%
%   % child symmetry
%   cs_bcc = crystalSymmetry('m-3m', [2.866 2.866 2.866], 'mineral', 'Iron bcc')
%
%   % define a bcc child orientation
%   ori_bcc = orientation.goss(cs_bcc)
%
%   % define Nishiyama Wassermann fcc to bcc orientation relation ship
%   NW = orientation.NishiyamaWassermann (cs_fcc,cs_bcc)
%
%   % compute a fcc parent orientation related to the bcc child orientation
%   ori_fcc = ori_bcc * NW
%
%   % compute all symmetrically possible parent orientations
%   ori_fcc = unique(ori_bcc.symmetrise * NW)
%
%   % same using the function parents
%   ori_fcc2 = ori_bcc * NW.parents
%
% See also
% orientation/variants
%

% store child symmetry
CS_child = mori.SS;

% symmetrise only with respect to child symmetry
mori = CS_child * mori;

% ignore all variants symmetrically equivalent
% with respect to the parent symmetry
mori.SS = crystalSymmetry('1');
mori = unique(mori);
mori.SS = CS_child;

author: mtex-toolbox

repo: mtex

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


text:
function test11results
%TEST11RESULTS analyze results from test11.m
% Example:
%   test11results
% See also test11, cholmod_test

% Copyright 2006-2007, Timothy A. Davis, University of Florida

load Results
index = UFget ;

c = E1(1:kkk) < 1 & T1(1:kkk) > 0 ;
m = E2(1:kkk) < 1 & T2(1:kkk) > 0 ;
cgood = find (c) ;	%#ok
mgood = find (m) ;	%#ok
good  = find (c | m) ;
bad = find (~(c|m)) ;

fl_per_lnz = FL(1:kkk) ./ LNZ(1:kkk) ;
speedup = T1(1:kkk) ./ T2(1:kkk) ;

[ignore ii] = sort (fl_per_lnz (good)) ;
good = good (ii) ;

fprintf ('MATLABtime CHOLMOD(time,flop,nnz(L)) speedup problem\n') ;
for k = good
    i = f (k) ;
%    fprintf ('%4d: t1 %10.2f t2 %10.2f fl %6.1e lnz %6.1e   %s/%s\n', ...
%	i, T1(k), T2(k), FL(k), LNZ(k), index.Group{i}, index.Name{i}) ;
     fprintf ('%10.4f %10.4f  %6.1e  %6.1e  %5.2f   %s/%s\n', ...
	T1(k), T2(k), FL(k), LNZ(k), speedup(k), ...
	index.Group{i}, index.Name{i}) ;
end

fprintf ('\nfailed in both:\n') ;
for k = bad
     i = f (k) ;
     fprintf ('%10.4f %10.4f  %6.1e  %6.1e  %5.2f   %s/%s\n', ...
	T1(k), T2(k), FL(k), LNZ(k), speedup(k), ...
	index.Group{i}, index.Name{i}) ;
end

% figure (3)
clf
loglog (fl_per_lnz (good), speedup (good), 'x') ;
axis ([1 4000 .1 50]) ;
xlabel ('Cholesky flop count / nnz(L)') ;
ylabel ('(MATLAB x=A\\b time) / (CHOLMOD time)') ;

drawnow

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ cimages ] = getSurfaceLabel(imdir,imagename,workspcdir, vpdata)

moveFolder = './roomModel/SpatialLayout_shrink/';

outimgdir = [workspcdir 'Images/'];
if ~exist(outimgdir,'dir')
    mkdir(outimgdir);
end
workspcdir = [workspcdir 'data/'];
if ~exist(workspcdir,'dir')
    mkdir(workspcdir);
end

% tempdir='./sptiallayouttempworkspace/';
%Compute vps
boxlayout=[];
surface_labels=[];

img=imread([imdir imagename]);
[h w kk]=size(img);
if ~exist('vpdata','var')
    [vp p All_lines]=getVP(imdir,imagename,0,workspcdir);
    VP=vp;
    if numel(VP)<6
        aa = zeros(h,w);
        indd = zeros(h,w);
        return;
    end

    vp=[VP(1) VP(2);VP(3) VP(4);VP(5) VP(6)];
    [vp P]=ordervp(vp,h,w,p);
    [vv linemem]=max(P,[],2);
    vpdata.vp=vp;
    vpdata.lines=All_lines;
    vpdata.linemem=linemem;
    vpdata.dim=[h w];
%     visvp(vpdata,img)
end

%Get segmentation and GC surface confidence maps
% sigma=num2str(0.8);
% k1=num2str(100);
% min1=num2str(100);
% inputim=[workspcdir imagename(1:end-4) '.ppm'];
% outputim=[workspcdir imagename(1:end-4) '.pnm'];
% [s w]=system(['./segment' ' ' sigma ' '  k1  ' ' min1 ' ' inputim ' ' outputim ]);

segext='pnm';
nsegments=[5 15 25 35 40 60 80 100];
% fn=['../Imsegs/' imagename(1:end-4) '.' segext];
fn=[moveFolder '/Imsegs/' imagename(1:end-4) '.' segext];
imseg = processSuperpixelImage(fn);

tic
imdata = mcmcComputeImageData(im2double(img), imseg);% made changes here
toc

load(fullfile([moveFolder '/LabelClassifier/'], 'Classifiers_gc.mat'));

spfeatures = mcmcGetAllSuperpixelData(imdir, imseg);
[efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imseg(1));

nsegments=[5 15 25 35 40 60 80 100];

pE{1} = test_boosted_dt_mc(eclassifier, efeatures{1});
pE{1} = 1 ./ (1+exp(ecal(1)*pE{1}+ecal(2)));
smaps{1} = msCreateMultipleSegmentations(pE{1}, adjlist{1}, ...
    imseg(1).nseg, nsegments);

for k = 1:numel(nsegments)
    if max(smaps{1}(:, k))>0
        segfeatures{1, k} = mcmcGetSegmentFeatures(imseg, ...
            spfeatures{1}, imdata, smaps{1}(:, k), (1:max(smaps{1}(:, k))));
    end
end

%Get surface label confidences initial from GC
normalize = 1;
pg=zeros(imseg.nseg,7);%7 labels
%get P(L/I)
pg = msTest(imseg, segfeatures, smaps, ...
    labelclassifier, segclassifier,normalize);

filename=fullfile(workspcdir, [imagename(1:end-4) '_lc_gc.mat' ]);
save(filename,'pg');

%visualize

%Compute intergral images for features
tic
[integData]=getIntegralimages([imagename],vpdata,imseg,500,workspcdir,imdir);
toc

%Get candidate layouts and their features
[polyg, Features] = getcandboxlayout( vpdata.vp,vpdata.dim(1),vpdata.dim(2),integData);

%Get initial estimate
Features1=Features;

% load ../LearntClassifiers/pf_i.mat
load([moveFolder '/LearntClassifiers/pf_i.mat']);

lay_score=[];
numL=size(Features,1);
% change features
if(numel(weights)==14)
    tmpf1=sum(Features(:,1:5).*Features(:,11:15),2);
    tmpf2=sum(Features(:,1:5).*Features(:,16:20),2);
    tmpf3=sum(Features(:,6:10).*Features(:,11:15),2);
    tmpf4=sum(Features(:,6:10).*Features(:,16:20),2);
    Features=[Features(:,1:10) tmpf1 tmpf2 tmpf3 tmpf4];

end
%evaluate
score=repmat(weights,[numL,1]).*Features;
score=sum(score,2);
[vv ii]=sort(score,'descend');

boxlayout.polyg=polyg;
boxlayout.init=[vv ii];

% load ../LearntClassifiers/pf_il.mat
load([moveFolder '/LearntClassifiers/pf_il.mat']);
Features=Features1;
if(numel(weights)==59)
    tmpf1=sum(Features(:,1:5).*Features(:,11:15),2);
    tmpf2=sum(Features(:,1:5).*Features(:,16:20),2);
    tmpf3=sum(Features(:,6:10).*Features(:,11:15),2);
    tmpf4=sum(Features(:,6:10).*Features(:,16:20),2);
    Features=[Features(:,1:10) tmpf1 tmpf2 tmpf3 tmpf4 Features(:,31:75) ];%Features(:,114:122)];

end
score=repmat(weights,[numL,1]).*Features;
score=sum(score,2);
[vv ii]=sort(score,'descend');

boxlayout.reestimated=[vv ii];

lay_scores=[vv ii];
save([workspcdir imagename(1:end-4) '_layres.mat'],'polyg','lay_scores');%,'avg_pg');

figure(101);
drawnow;
for lay=1:25
    layoutid=ii(lay);
    Polyg=[];
    for fie=1:5
        Polyg{fie}=[];
        if size(polyg{layoutid,fie})>0
            Polyg{fie}=polyg{layoutid,fie};
        end
    end

    tempimg=displayout(Polyg,w,h,img);
    subplot(5,5,lay);imshow(uint8(tempimg),[]);title(num2str(vv(lay)));
end
saveas(101,[outimgdir imagename(1:end-4) '_boxlayouts.png']);

Polyg=[];
for fie=1:5
    Polyg{fie}=[];
    if size(polyg{ii(1),fie})>0
        Polyg{fie}=polyg{ii(1),fie};
    end
end
ShowGTPolyg(img,Polyg,103);
saveas(103,[outimgdir imagename(1:end-4) '_boxlayout.png']);

%Re-Compute Surface labels (GC+box layout features)
% load(['../LabelClassifier/' 'Classifiers_stage2.mat']);
load([moveFolder '/LabelClassifier/' 'Classifiers_stage2.mat']);
tic
xspfields=[];%save per image
numSup=imseg.nseg;
for lay=1:numL %each layout
    xspf=[];
    for  supno=1:numSup
        tempbndy=imdata.tracedbndy{supno}{1};
        if size(tempbndy,1) > 30

            YY1=tempbndy(1:10:end,1);XX1=tempbndy(1:10:end,2);
        else

            YY1=tempbndy(:,1);XX1=tempbndy(:,2);
        end

        for fi=1:5  %each field
            xarea=0;
            if size(polyg{lay,fi},1)>0

                XX2=polyg{lay,fi}(:,1);
                YY2=polyg{lay,fi}(:,2);

                [in on]=inpolygon(XX1,YY1,[XX2;XX2(1)],[YY2;YY2(1)]);

                if numel(find(in==1))==length(in)
                    X0=XX1;Y0=YY1;
                    xarea=polyarea([X0;X0(1)],[Y0;Y0(1)]);
                elseif  numel(find(in==1))==0

                    X0=[];Y0=[];
                    xarea=0;
                else
                    xarea=polyintarea(XX1,YY1,XX2,YY2,0);

                end

            end
            xspf(supno,fi)=xarea;
        end

    end
    xspfields{lay}=xspf;
end
toc

nsp=imseg.nseg; %num of suppixels
smap = [1:nsp];
smaps{1} = smap(:);

clear segfeatures
for k = 1: 1%numel(nsegments)
    features = mcmcGetSegmentFeatures(imseg, ...
        spfeatures{1}, imdata, smaps{1}(:, k), (1:max(smaps{1}(:, k))));

    for lay=1:numL
        tempfeatures=features;
        [fieldfeatures]=segfieldsfeat_sup(imseg,imseg.nseg,smaps{1}(:,k),xspfields{lay});
        tempfeatures(:,95:100)=fieldfeatures(1:size(features,1),:);
        tempfeatures(:,101:106)=pg{1}(:,1:6);
        segfeatures{lay,k}  = tempfeatures;
    end
end

if numel(vv)>100
    A = [vv(1) 1;vv(100) 1];
else
    A = [vv(1) 1;vv(end) 1];
end
b = [0;log(49)];
params = inv(A)*b;

%     conf = 1./(1+exp(params(1)*vv(1:5)+params(2)));
conf = 1./(1+exp(params(1)*vv(1)+params(2)));
conf=conf./sum(conf);

avg_pg=zeros(imseg.nseg,7);%7 labels

tic
for j=1:1 %take best 5 scored layouts
    lay=ii(j);

    %get P(L/F,I)
    pg = msTest(imseg, segfeatures(lay, :), smaps, ...
        labelclassifier);%, segclassifier,normalize);
    avg_pg=avg_pg+pg{1}.*conf(j);

end
toc

filename=fullfile(workspcdir,[imagename(1:end-4) '_lc_st2.mat' ]);
save(filename,'avg_pg');

surface_labels.restimated={avg_pg};
surface_labels.init=pg;

rl=[255  255 0  0 0 255 0];
gl=[255 0 255  255 0 0 0];
bl=[0  0 255 0 255 255 0];
rl1=[255 202 132 255 0 255 0 191 0 0];
gl1=[236 255 112 255 191 0 255 21 0 0];
bl1=[139 112 255  255 255 0 0 133 255 0];
cimages = msPg2confidenceImages(imseg,{avg_pg});
[aa indd]=max(cimages{1}(:,:,1:6),[],3);

figure(102);
clear mask_color;
mask_r = rl(indd);
mask_g = gl(indd);
mask_b = bl(indd);
mask_color(:,:,1) = mask_r;
mask_color(:,:,2) = mask_g;
mask_color(:,:,3) = mask_b;

hsvmask=rgb2hsv(mask_color);
hsvmask(:,:,3)=aa*255;
%     hsvmask(:,:,2)=aa;
mask_color=hsv2rgb(hsvmask);

tempimg = double(img)*0.5 + mask_color*0.5;
%         tempimg =  mask_color;
imshow(uint8(tempimg));
saveas(102,[outimgdir imagename(1:end-4) '_surfacelabels.png']);

end

author: zouchuhang

repo: LayoutNet


text:
function y = s2y(s);

% Y = s2y(S)
%
% Scattering to Admittance transformation
%
% y = (I-s) * inv(I+s)
%
% for square matrices at multiple frequencies
%
% 27.09.2002

I = diag(ones(1, size(s,2)));

for i=1:size(s,3)
   y(:,:,i) = (I-s(:,:,i)) * inv(I+s(:,:,i));
end;

author: Sable

repo: mcbench-benchmarks


text:
function [l,d,perm] = mchol(A,mu)
% [l,d,perm] = mchol(A,mu)
% Compute the Gill-Murray modified LDL factorization of A,

if nargin < 2
    mu = 1e-12;
end

n = size(A,1);
l = eye(n);
d = zeros(n,1);
perm = 1:n;

for i = 1:n
    c(i,i) = A(i,i);
end

% Compute modification parameters
gamma = max(abs(diag(A)));
xi = max(max(abs(setdiag(A,0))));
delta = mu*max(gamma+xi,1);
if n > 1
    beta = sqrt(max([gamma xi/sqrt(n^2-1) mu]));
else
    beta = sqrt(max([gamma mu]));
end

for j = 1:n

    % Find q that results in Best Permutation with j
    [maxVal maxPos] = max(abs(diag(c(j:end,j:end))));
    q = maxPos+j-1;

    % Permute d,c,l,a
    d([j q]) = d([q j]);
    perm([j q]) = perm([q j]);
    c([j q],:) = c([q j],:);
    c(:,[j q]) = c(:,[q j]);
    l([j q],:) = l([q j],:);
    l(:,[j q]) = l(:,[q j]);
    A([j q],:) = A([q j],:);
    A(:,[j q]) = A(:,[q j]);

    for s = 1:j-1
        l(j,s) = c(j,s)/d(s);
    end
    for i = j+1:n
        c(i,j) = A(i,j) - sum(l(j,1:j-1).*c(i,1:j-1));
    end
    theta = 0;
    if j < n
        theta = max(abs(c(j+1:n,j)));
    end
    d(j) = max([abs(c(j,j)) (theta/beta)^2 delta]);
    if j < n
        for i = j+1:n
            c(i,i) = c(i,i) - (c(i,j)^2)/d(j);
        end
    end
end


text:
%%%% EVALUATES CONSTRIANT ERROR

% Course: Robotic Manipulation and Mobility
% Advisor: Dr. V. Krovi
%
% Homework Number: MIDTERM
%
% Names: Sourish Chakravarty
% 	Hrishi Lalit Shah

function [Cerr]= CONSTRAINT_ERROR_1(th1,x2,y2,th2)

global l1 lc1 l2 lc2

Cerr=[];

for i=1:length(th1)

    c1=cos(th1(i));
    c2=cos(th2(i));
    s1=sin(th1(i));
    s2=sin(th2(i));

    C= [x2(i)-l1*c1-lc2*c2;
        y2(i)-l1*s1-lc2*s2];
    Cerr=[Cerr;abs(C')];
end

return


text:
function dataset = voc0712plus_trainval_ss(dataset, usage, use_flip)
% Pascal voc 0712 trainval set with selective search
% set opts.imdb_train opts.roidb_train
% or set opts.imdb_test opts.roidb_train

% change to point to your devkit install
devkit2007                      = voc2007_devkit();
devkit2012                      = voc2012_devkit();

switch usage
    case {'train'}
        dataset.imdb_train    = {  imdb_from_voc(devkit2012, 'trainval', '2012', use_flip), ...
                                    imdb_from_voc(devkit2007, 'trainval', '2007', use_flip), ...
                                    imdb_from_voc(devkit2007, 'test', '2007', use_flip)};
        dataset.roidb_train   = cellfun(@(x) x.roidb_func(x, 'with_selective_search', true), dataset.imdb_train, 'UniformOutput', false);
    case {'test'}
        error('only supports one source test currently');
    otherwise
        error('usage = ''train'' or ''test''');
end

end

author: ShaoqingRen


text:
function [imageVol, lastfile] = GE_readVolume(baseFileName, volSize, depth, im_offset)
%
%GE_readVolume
%
% [imageVol, lastfile] = GE_readVolume(baseFileName, [nX nY nZ], depth, im_offset)
%
% reads the volume for passnum from the series which is stored
% starting in startDir and returns the name of the last file read
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%
% 2002.01.24 RFD: code needs comments! Apparently, this code assumes a particular dir
%   structure involving numbered 'run' directories. I simplified things by just assuming
%   all files are in the same dir.
%

% initialize some variables
nX = volSize(1);
nY = volSize(2);
nZ = volSize(3);
sliceSize = nX*nY;
imageVol = zeros(nX, nY, nZ);

for i = 1:nZ
    imNum = sprintf('.%.3d', i);
    imageFile = fullfile('', [baseFileName,imNum]);
    % Open the file
    [fid,message] = fopen(imageFile,'r','b');
    if (fid == -1)
        fprintf('Cannot Open %s.\n',imageFile);
        break
    end

    % Skip to the data
    fseek(fid,im_offset,-1);
    % Read the slice
    buffer = fread(fid,sliceSize,sprintf('int%d',depth));
    % append the slice to the imageSet
    imageVol(:,:,i) = reshape(buffer, nX, nY);

    % Close the file
    status = fclose(fid);

    % Set the lastfile
    lastfile = imageFile;
end

return


text:
function test_suite = test_vectorNorm
% One-line description here, please.
%   output = testVectorNorm(input)
%
%   Example
%   testVectorNorm
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testEuclidean(testCase) %#ok<*DEFNU>

v = [3 4];
norm = vectorNorm(v);
testCase.assertEqual(5, norm, 'AbsTol', .01);

function testEuclideanArray(testCase)

v = [3 4;4 3;6 8;5 12];
norm = vectorNorm(v);
testCase.assertEqual([5;5;10;13], norm, 'AbsTol', .01);

function testExplicitEuclideanArray(testCase)

v = [3 4;4 3;6 8;5 12];
norm = vectorNorm(v, 2);
testCase.assertEqual([5;5;10;13], norm, 'AbsTol', .01);

function testNorm1Array(testCase)

v = [3 4;4 3;6 8;5 12];
norm = vectorNorm(v, 1);
testCase.assertEqual([7;7;14;17], norm, 'AbsTol', .01);

function testNormInfArray(testCase)

v = [3 4;4 3;6 8;5 12];
norm = vectorNorm(v, inf);
testCase.assertEqual([4;4;8;12], norm, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% FMRIB1
%
% Files
%   decimate2           - Clone of decimate.m by the mathworks
%   eegplugin_fmrib     - eegplugin_fmrib() -       EEGLAB plugin from the Centre for
%   fastranc            - fastranc.m -  adaptive noise cancellation for fmrib_fastr.m
%   fmrib_fastr         - fmrib_fastr() -  Remove FMRI gradient artifacts from EEG data using
%   fmrib_pas           - fmrib_pas() - Remove pulse
%   fmrib_qrsdetect     - fmrib_qrsdetect() - Detect QRS peaks from ECG channel using combined
%   guifastr            - Application M-file for guifastr.fig
%   guipas              - Application M-file for guifastr.fig
%   pca_calc            - Principal Component Analysis
%   pop_fmrib_fastr     - pop_fmrib_fastr() -  GUI for fmrib_fastr: Remove gradient
%   pop_fmrib_pas       - pop_fmirb_pas() -  GUI for fmrib_pas: Remove hear-related pulse artifacts
%   pop_fmrib_qrsdetect - pop_fmrib_qrsdetect() -  GUI for fmrib_qrsdetect:
%   prcorr2             - - Compute correlation coefficient over all dimentions
%   qrscorrect          - qrscorrect() -  Correct for false positive and negative QRS peaks
%   trigcorrect         - trigcorrect() -   correct for missing FMRI slice triggers

author: PatternRecognition

repo: OpenBMI


text:
function tssk = t_mp_dm_converter_mpc2(quiet)
%T_MP_DM_CONVERTER_MPC2  Tests for MP.DM_CONVERTER_MPC2.

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

% define_constants;
if quiet
    verbose = 0;
else
    verbose = 1;
end

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case_ext';
mpc0 = loadcase(casefile);
mpc0.bus(end, MU_VMIN) = 0;         %% add columns
mpc0.gen(end, MU_QMIN) = 0;         %% add columns
mpc0.branch(end, MU_ANGMAX) = 0;    %% add columns
dmc = mp.dm_converter_mpc2().build();

t_begin(27, quiet);

t = 'dmc constructor : ';
dmc = mp.dm_converter_mpc2();
t_ok(isa(dmc, 'mp.dm_converter'), [t 'class']);
t_is(length(dmc.element_classes), 5, 12, [t '# of element_classes']);
t_ok(isempty(dmc.elements), [t 'elements empty']);

t = 'dmc.build() : ';
dmc.build();
t_is(length(dmc.elements), 5, 12, [t '# of elements']);
t_str_match(dmc.elements{1}.name, 'bus', [t 'element{1} is bus']);
t_str_match(dmc.elements{2}.name, 'gen', [t 'element{2} is gen']);
t_str_match(dmc.elements{3}.name, 'load', [t 'element{3} is load']);
t_str_match(dmc.elements{4}.name, 'branch', [t 'element{4} is branch']);
t_str_match(dmc.elements{5}.name, 'shunt', [t 'element{5} is shunt']);

t = 'dm constructor : ';
dm = mp.data_model_cpf();
t_ok(isa(dm, 'mp.data_model'), [t 'class']);
t_is(length(dm.element_classes), 5, 12, [t '# of element_classes']);
t_ok(isempty(dm.elements), [t 'elements empty']);

t = 'dm.build(mpc, dmc) : ';
dm.build(mpc0, dmc);
t_is(length(dm.elements), 5, 12, [t '# of elements']);
t_str_match(dm.elements{1}.name, 'bus', [t 'element{1} is bus']);
t_str_match(dm.elements{2}.name, 'gen', [t 'element{2} is gen']);
t_str_match(dm.elements{3}.name, 'load', [t 'element{3} is load']);
t_str_match(dm.elements{4}.name, 'branch', [t 'element{4} is branch']);
t_str_match(dm.elements{5}.name, 'shunt', [t 'element{5} is shunt']);

t = 'mpc = dmc.export(dm, mpc0) : ';
mpc = dmc.export(dm, mpc0);
t_ok(isequal(mpc, mpc0), [t 'mpc == mpc0']);

t = 'mpc = dmc.export(updated_dm, mpc0) : ';
mpc1 = mpc0;
mpc1.bus([2;3;6;9], VM) = [1.02;1.03;1.06;1.09];
dm.elements.bus.tab.vm([2;3;6;9]) = [1.02;1.03;1.06;1.09];
mpc1.bus([3;5;7], PD) = [30;50;70];
mpc1.bus([3;5;7], QD) = [18;30;42];
load0 = dm.elements.load.tab;
load = [load0; load0(1:2, :)];
k = [4;1;5];
load.uid(k) = [4;1;5];
load.bus(k) = [30;5;6];
load.source_uid(k) = [3;5;7];
load.pd(k) = [30;50;70];
load.qd(k) = [18;30;42];
dm.elements.load.tab = load;
dm.elements.load.rebuild(dm);
mpc = dmc.export(dm, mpc0);
t_ok(~isequal(mpc, mpc0), [t 'mpc ~= mpc0']);
t_ok( isequal(mpc, mpc1), [t 'mpc == updated_mpc']);

t = 'mpc = dmc.export(dm) : ';
mpc1 = struct( ...
    'version', mpc0.version, ...
    'baseMVA', mpc0.baseMVA, ...
    'bus', mpc0.bus, ...
    'gen', mpc0.gen, ...
    'branch', mpc0.branch, ...
    'gencost', mpc0.gencost, ...
    'bus_name', {mpc0.bus_name} );
dmc = mp.dm_converter_mpc2().build();
dm = mp.data_model_opf().build(mpc1, dmc);
mpc = dmc.export(dm);
mpc1.gen(:, MBASE) = 0;
mpc1.branch(7, BR_STATUS) = 0;
t_ok(isequal(mpc, mpc1), [t 'mpc == mpc0']);

t = 'mpc = dmce.export(dme, mpc, ''qcost'') : ';
%% define named indices into data matrices
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
dme = dm.elements.gen;
dmce = dmc.elements.gen;
dme.tab.cost_qg = dme.tab.cost_pg;
dme.tab.cost_qg.poly_coef = dme.tab.cost_qg.poly_coef ./ 5;
ng = size(mpc.gen, 1);
epc = mpc.gencost(1:ng, :);
eqc = epc; eqc([1;3], COST) = eqc([1;3], COST) / 5;
mpc = dmce.export(dme, mpc, {'cost_qg'});
t_is(size(mpc.gencost, 1), 2*ng, 12, [t 'size(gencost, 1)']);
t_is(mpc.gencost(1:ng,:), epc, 12, [t 'pcost']);
t_is(mpc.gencost(ng+1:2*ng,:), eqc, 12, [t 'qcost']);

t = 'mpc = dmce.export(dme, mpc, ''qcost'', ridx) : ';
ridx = [2;3];
dme.tab.cost_pg.poly_coef = dme.tab.cost_pg.poly_coef * 2;
dme.tab.cost_pg.pwl_qty = dme.tab.cost_pg.pwl_qty + 10;
dme.tab.cost_pg.pwl_cost = dme.tab.cost_pg.pwl_cost + 10;
dme.tab.cost_qg.poly_coef = dme.tab.cost_qg.poly_coef * 3;
dme.tab.cost_qg.pwl_qty = dme.tab.cost_qg.pwl_qty + 1;
dme.tab.cost_qg.pwl_cost = dme.tab.cost_qg.pwl_cost + 1;
epc(2, COST:end) = epc(2, COST:end) + 10;
epc(3, COST:end) = epc(3, COST:end) * 2;
eqc(2, COST:end) = eqc(2, COST:end) + 1;
eqc(3, COST:end) = eqc(3, COST:end) * 3;
mpc = dmce.export(dme, mpc, {'cost_pg', 'cost_qg'}, ridx);
t_is(mpc.gencost(1:ng,:), epc, 12, [t 'pcost']);
t_is(mpc.gencost(ng+1:2*ng,:), eqc, 12, [t 'qcost']);

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function r_size  = nwspgr_size ( type, dim, k, sym, compress )

%*****************************************************************************80
%
%% NWSPGR_SIZE determines the size of a sparse grid rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2014
%
%  Author:
%
%    Original MATLAB version by Florian Heiss, Viktor Winschel.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Florian Heiss, Viktor Winschel,
%    Likelihood approximation by numerical integration on sparse grids,
%    Journal of Econometrics,
%    Volume 144, 2008, pages 62-80.
%
%  Parameters:
%
%    Input, string TYPE, selects the 1D integration rule:
%    func:  any function name. Function must accept level l and return nodes N
%    and weights W for univariate quadrature rule with polynomial exactness
%    2*L-1 as [n, w] = feval(func,level)
%
%    Input, integer DIM, the dimension of the integration problem.
%
%    Input, integer K, the level.
%
%    Input, integer SYM.  If the rule is
%    symmetric, specifying SYM to be 1 will allow the code to run faster.
%    But this also requires that the user quadrature rule function only
%    return the nonnegative abscissas and their weights.
%
%    Input, integer COMPRESS,
%    0, do not compress the rule ( = do not merge duplicate points.)
%    1, compress the rule.
%
%    Output, integer R_SIZE, the "size" of the rule.
%
  sym = logical ( sym );
%
%  Retrieve the 1D rules of levels 1 to K.
%  Create cell arrays X1D and W1D containing the node and weight information.
%  The array N1D stores the length or order of each rule.
%
%  If the user indicates that it's a symmetric rule,
%  keep only the positive orthant.
%
  n1d = zeros ( 1, k );
  x1d = cell ( k, 1 );
  w1d = cell ( k, 1 );

  for level = 1 : k

    [ x, w ] = feval ( type, level );

    if ( sym )
      [ numnew, dummy ] = size ( x );
      [ x, sortvec ] = sortrows ( x );
      w = w ( sortvec );
      x = x((floor(numnew/2)+1):numnew,:);
      w = w((floor(numnew/2)+1):numnew,:);
    end

    n1d(level) = length ( w );
    x1d{level} = x;
    w1d{level} = w;

  end
%
%  Initialization.
%
  minq = max ( 0, k - dim );
  maxq = k - 1;
  nodes = [];
  weights = [];
%
%  Loop for max ( 0, K - DIM ) <= Q <= K - 1.
%
  for q = minq : maxq

    r = length ( weights );
%
%  BQ is the combinatorial coefficient applied to the component
%  product rules which have level Q.
%
    bq = ( -1 )^( maxq - q ) * nchoosek ( dim - 1, dim + q - k );
%
%  Compute the D-dimensional row vectors that sum to DIM+Q.
%
    is = get_seq ( dim, dim + q );
%
%  Preallocate new rows for nodes and weights.
%
    Rq = prod ( n1d(is), 2 );
    sRq = sum ( Rq );
    nodes   = [ nodes;   zeros(sRq,dim) ];
    weights = [ weights; zeros(sRq,1) ];
%
%  Generate each of the product rules indicated by IS, and
%  insert them into NODES and WEIGHTS.
%
    for j = 1 : size(is,1)
      midx = is(j,:);
      [ newn, neww ] = tensor_product ( x1d(midx), w1d(midx) );
      nodes((r+1):(r+Rq(j)),:) = newn;
      weights((r+1):(r+Rq(j))) = bq .* neww;
      r = r + Rq(j);
    end
%
%  Sort the nodes and merge repeated values.
%
    [ nodes, sortvec ] = sortrows ( nodes );
    weights = weights(sortvec);
    keep = 1;
    lastkeep = 1;

    for j = 2 : size(nodes,1)
      if ( compress )

        if ( nodes(j,:) == nodes(j-1,:) )
          weights(lastkeep) = weights(lastkeep) + weights(j);
        else
          lastkeep = j;
          keep = [ keep ; j ];
        end
      else
        lastkeep = j;
        keep = [ keep ; j ];
      end
    end

    nodes = nodes(keep,:);
    weights = weights(keep);

  end
%
%  The rule has been computed.
%  If we used symmetry, we now have to extend the rule.
%
  if ( sym )

    nr = length ( weights );
    m = x1d{1};

    for j = 1 : dim

      keep = zeros(nr,1);
      numnew = 0;

      for r = 1 : nr
        if ( nodes(r,j) ~= m )
          numnew = numnew + 1;
          keep(numnew) = r;
        end
      end

      if ( 0 < numnew )
        nodes = [nodes ; nodes(keep(1:numnew),:)];
        nodes(nr+1:nr+numnew,j) = 2*m - nodes(nr+1:nr+numnew,j);
        weights = [weights ; weights(keep(1:numnew))];
        nr = nr + numnew;
      end

    end

    [ nodes, sortvec ] = sortrows ( nodes );
    weights = weights(sortvec);

  end

  r_size = length ( weights );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y=ewmacovariance(R1,R2,n,s,lambda)
%
% implemented by ewmaestimatevar function
%
%% Written By Ali Najjar
%
if s==1
    y=R1(n-s)*R2(n-s);
else
    y=((1-lambda)*((R1(n-1))*(R2(n-1))))+(lambda*ewmacovariance(R1,R2,n-1,s-1,lambda));
end
end


text:
function [hy,model,Sigma] = ml_gpr(X,y,model,epsilon,rbf_var,method)
%ML_GP Gaussian Process Regression
%
%   Desired structure for ml_kcv ------------------------------------------
%
%  [hy,model] = f(X_train,labels_train,model) : used to train
%
%  hy         = f(X_test,[],model) : used to predict
%
%   input -----------------------------------------------------------------
%
%       o X         : (N x D), N input samples of dimension D.
%
%       o y         : (N x 1), N output predictions of dimension 1.
%
%       o model     : struct
%
%       o epsilon   : measurement noise
%
%       o rbf_var   : kernel width
%
%   output ----------------------------------------------------------------
%
%
%

Sigma = [];

if ~exist('method','var'), method='gpml'; end

if ~isempty(y)
% Train the model
    model.X_train = X;
    model.y_train = y;

    if strcmp(method,'gdc')

        [hy,Sigma] = gaussian_process(X,y,X,epsilon,rbf_var);

    elseif strcmp(method,'gpml')

        meanfunc     = {@meanZero};
        covfunc      = {@covSEiso};
        ell          = rbf_var;     % kernel width of RBF covariance function.
        sf           = 1;           % signal variance (not measurement noise)
        sn           = epsilon;      % measurement noise
        hyp          = [];
        hyp.cov      = log([ell; sf]);
        hyp.lik      = log(sn);
        [hy,Sigma]   = gp(hyp, @infExact, meanfunc, covfunc, @likGauss, X, y, X);

    else

        error(['No such method defined: ' method]);

    end

else
% Test the model

    if strcmp(method,'gdc')

        [hy,Sigma]    = gaussian_process(model.X_train,model.y_train,X,epsilon,rbf_var);

    elseif strcmp(method,'gpml')

        meanfunc     = {@meanZero};
        covfunc      = {@covSEiso};
        ell          = rbf_var;     % kernel width of RBF covariance function.
        sf           = 1;           % signal variance (not measurement noise)
        sn           = epsilon;      % measurement noise
        hyp          = [];
        hyp.cov      = log([ell; sf]);
        hyp.lik      = log(sn);

        [hy,Sigma]   = gp(hyp, @infExact, meanfunc, covfunc, @likGauss, model.X_train,model.y_train, X);

    else

        error(['No such method defined: ' method]);

    end

end

end

author: epfl-lasa


text:
function mff_setup

% MFF_SETUP ensures that the JVM is available and that the MFF java
% implementation is on the java classpath.

% this toolbox requires JAVA to be running and properly configured
if ~usejava('jvm')
  error('the EGI_MFF format requires that the java virtual machine is available, see http://www.fieldtriptoolbox.org/getting_started/egi')
end

% check that the java archive (jar file) is on the path
% the user can add it to the static path using classpath.txt
if ~(exist('com.egi.services.mff.api.MFFFactory', 'class')==8)
    jarpath = fullfile(fileparts(mfilename('fullpath')), 'java', 'M*.jar');
    jarlist = dir(jarpath);
    if size(jarlist, 1) == 1

        % Some MATLAB versions have a bug in javaclasspath that causes global
        % variables to be erased. Hence we keep a local copy and restore it afterwards.
        % See http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2133

        var = whos('global');
        glob = [var.global];
        name = {var.name};
        name = name(glob);
        for i=1:length(name)
        eval(sprintf('global %s', name{i}));
        eval(sprintf('localcopy.%s = %s;', name{i}, name{i}));
        end

        % try to add it to the dynamic classpath
        jarpath = fullfile(fileparts(jarpath), jarlist(1).name);
        warning('adding %s to your Java path', jarpath);
        javaclasspath(jarpath);

        % The following pertains to http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2133
        % also available as http://bit.ly/1o7QlRk, for which we cannot provide a solution.
        warning('In some versions of Matlab, javaclasspath clears all global variables due to a bug in MATLAB. Please execute %s.m immediately after starting MATLAB. See http://bit.ly/1o7QlRk', mfilename('fullpath'));

        % restore the global variables
        for i=1:length(name)
        eval(sprintf('global %s', name{i}));
        eval(sprintf('%s = localcopy.%s;', name{i}, name{i}));
        end
        clear localcopy

    else
        error('the EGI_MFF format requires one MFF .jar file on your classpath, see http://www.fieldtriptoolbox.org//getting_started/egi')
    end
end

% check once more, give an error if it still does not exist
if ~(exist('com.egi.services.mff.api.MFFFactory', 'class')==8)
  error('the EGI_MFF format requires one MFF .jar file on your classpath, see http://www.fieldtriptoolbox.org/getting_started/egi')
end


text:
function vw = makeAnnotationString(vw)
%
% function vw = makeAnnotationString(vw)
%
% djh, 3/2001
% ras, 11/2004: made it a uicontrol.
pos = [.15 .95 .5 .04];

bgcolor = get(gcf,'Color');

h = uicontrol('Style','text',...
              'Units','Normalized',...
              'Position',pos,...
              'FontName','Helvetica',...
              'FontSize',10,...
              'FontWeight','bold',...
              'BackgroundColor',bgcolor,...
              'ForegroundColor',[.1 0 0],...
              'HorizontalAlignment','center',...
              'String','');

vw.ui.annotationHandle = h;

return

% OLD:
%
% set(INPLANE{1}.ui.annotationHandle,'string','mumble');%
% set(VOLUME{1}.ui.annotationHandle,'string','mumble');%
% set(FLAT{1}.ui.annotationHandle,'string','mumble');
% annotationPos = [0.5 0.95 0.4 0.04];
% annotationAxis = subplot('position',annotationPos);
% axis off;
% text(0,0,'','FontSize',12,'HorizontalAlignment','center',...
%      'FontWeight','bold');
% annotationHandle = get(annotationAxis,'Children');
% vw.ui.annotationHandle = annotationHandle;
%
% % Return the current axes to the main image
% set(gcf,'CurrentAxes',vw.ui.mainAxisHandle);


text:
function [model] = bp_resume(model)
% Resume discriminative finetuning of CDBN.

rng('shuffle');
kernels;
debug = 0;

data_list = read_data_list(model.data_path, model.classnames, ....
    model.volume_size + 2 * model.pad_size, 'train', debug);

param = [];
param.epochs = 100;
param.lr = 0.01;
param.weight_decay = 5*10^-4;
param.momentum = 0.9;
param.batch_size = 32;
param.snapshot_iter = 10;
param.snapshot_name = 'bp_finetune_iter';
param.test_iter = 5;
batch_size = param.batch_size;

fprintf('Resume discriminative funetuning the CDBN\n');
fprintf('lr = %f, wd = %d, momentum = %f\n', param.lr, param.weight_decay, param.momentum);

% prepare data and label
[new_list, label] = balance_data(data_list, batch_size);
n = length(new_list);
batch_num = n / batch_size;
assert(batch_num == floor(batch_num));

% prepare model
numLayer = model.numLayer;
for iter = 1 : param.epochs
    loss_all = 0;
    shuffle_index = randperm(n);
    for b = 1 : batch_num
        batch_index = shuffle_index((b-1)*batch_size + 1 : b * batch_size);
        batch = read_batch(model, new_list(batch_index), false);
        batch_label = label(batch_index,:);
        [model, activation] = bp_forward(model, batch);
        [model, loss] = bp_backward(model, activation, batch_label, param.weight_decay);
        loss_all = loss_all + loss;
        model = bp_update(model, param);
    end
    loss_all = loss_all / batch_num;
    fprintf('iteration: %d, loss: %f\n', iter, loss_all);

    if mod(iter, param.snapshot_iter) == 0
        fprintf('snapshoting to %s_%d\n', param.snapshot_name, iter);
        snapshot_name = sprintf('%s_%d', param.snapshot_name, iter);
        save(snapshot_name, 'model');
    end

    if mod(iter, param.test_iter) == 0
        test_loss = bp_test(model);
        fprintf('test loss: %f\n', test_loss);
    end
end

for l = 2 : numLayer
    model.layers{l} = rmfield(model.layers{l},'grdw');
    model.layers{l} = rmfield(model.layers{l},'grdc');
    model.layers{l} = rmfield(model.layers{l},'histw');
    model.layers{l} = rmfield(model.layers{l},'histc');
end

save('bp_finetuned_model', 'model');

author: zhirongw


text:
%% Create artificial low-rank image
perHei = 10;
perWid = 20;
hei = 50*perHei;
wid1 = 5*perWid;  rect1 = ones(hei, wid1) * 75;
wid2 = 15*perWid; rect2 = ones(hei, wid2) * 184;
wid3 = 8*perWid;  rect3 = ones(hei, wid3) * 223;
wid4 = 18*perWid; rect4 = ones(hei, wid4) * 113;
wid5 = 10*perWid; rect5 = ones(hei, wid5) * 38;
A = uint8([rect1 rect2 rect3 rect4 rect5]);
D = imnoise(A, 'salt & pepper');
A = double(A);
D = double(D);
E = D - A;
%% Create low-rank matrix
% m = 100; % m = 100, 200, 400, 800
% r = 5;  % rank(A) = 5, 10, 20, 40
% U = randn(m,r);
% V = randn(m,r);
% A = U*V';
% E = (imnoise(zeros(size(A)),'salt & pepper',0.1) > 0) .* (rand(size(A))-0.5)*2*500;
% D = A + E;
%% Apply APG to recover A_hat and E_hat
[rows, cols] = size(D);
lambda = rows^(-1/2);
tic
[A_hat,E_hat,numIter] = proximal_gradient_rpca(D, lambda);
toc

norm(A_hat-A, 'fro') / norm(A, 'fro')

figure; imshow(A,[]); title('original A');
figure; imshow(A_hat,[]); title('recovered A');

author: YimianDai

path: github-repos/MATLAB/YimianDai-Image-Processing-Codes-for-Easier-Understanding/Image-Processing-Codes-for-Easier-Understanding-874302799e48852624bc3760b58b46bd9360f238/src/(NIPS 2009) Robust Principal Component Analysis Exact Recovery of Corrupted Low-Rank Matrices by Convex Optimization/Demo_APG.m


text:
function [g,geq,dg,dgeq,xevaled] = fmincon_con(x,model,xevaled)

global latest_xevaled
global latest_x_xevaled
% Early bail for linear problems
g = [];
geq = [];
dg = [];
dgeq = [];
if model.linearconstraints
    xevaled = [];
    return
end

if nargin<3
    if isequal(x,latest_x_xevaled)
        xevaled = latest_xevaled;
    else
        xevaled = zeros(1,length(model.c));
        xevaled(model.linearindicies) = x;
        xevaled = apply_recursive_evaluation(model,xevaled);
        latest_x_xevaled = x;
        latest_xevaled = xevaled;
    end
end

if model.nonlinearinequalities
    g = full(model.Anonlinineq*xevaled(:)-model.bnonlinineq);
end

if nnz(model.K.q) > 0
    top = 1;
    for i = 1:length(model.K.q)
        z = model.F_struc(top:top+model.K.q(i)-1,:)*[1;xevaled];
      %  g = [g;-(z(1)^2 - z(2:end)'*z(2:end))];
        g = [g;-(z(1) - sqrt(z(2:end)'*z(2:end)))];
        top = top + model.K.q(i);
    end
end

if model.nonlinearequalities
    geq = full(model.Anonlineq*xevaled(:)-model.bnonlineq);
end

dgAll_test = [];

if nargout == 2 || ~model.derivative_available
    return
elseif ~isempty(dgAll_test) & isempty(model.evalMap)
    dgAll = dgAll_test;
elseif isempty(model.evalMap) & (model.nonlinearinequalities==0) & (model.nonlinearequalities==0) & (model.nonlinearcones==0) & any(model.K.q)
    dg = computeConeDeriv(model,xevaled);
elseif isempty(model.evalMap) & (model.nonlinearinequalities | model.nonlinearequalities | model.nonlinearcones)
    n = length(model.c);
    linearindicies = model.linearindicies;
 %   xevaled = zeros(1,n);
 %   xevaled(linearindicies) = x;
    % FIXME: This should be vectorized

    news = model.fastdiff.news;
    allDerivemt = model.fastdiff.allDerivemt;
    c = model.fastdiff.c;

    if model.fastdiff.univariateDifferentiates
        zzz = c.*(x(model.fastdiff.univariateDiffMonom).^model.fastdiff.univariateDiffPower);
    else
        %  X = repmat(x(:)',length(c),1);
        O = ones(length(c),length(x));
        nz = find(allDerivemt);
        %  O(nz) = X(nz).^allDerivemt(nz);
        O(nz) = x(ceil(nz/length(c))).^allDerivemt(nz);
        zzz = c.*prod(O,2);
    end

    newdxx = model.fastdiff.newdxx;
    newdxx(model.fastdiff.linear_in_newdxx) = zzz;
    %newdxx = newdxx';

    if ~isempty(model.Anonlineq)
        dgeq = model.Anonlineq*newdxx;
    end
    if ~isempty(model.Anonlinineq)
        dg = model.Anonlinineq*newdxx;
    end

    if nnz(model.K.q)>0
        dg = [dg;computeConeDeriv(model,xevaled,newdxx);];
    end
else
    requested = model.fastdiff.requested;
    dx = apply_recursive_differentiation(model,xevaled,requested,model.Crecursivederivativeprecompute);
    conederiv = computeConeDeriv(model,xevaled,dx);
    if ~isempty(model.Anonlineq)
        dgeq = [model.Anonlineq*dx];
    end
    if ~isempty(model.Anonlinineq)
        dg = [model.Anonlinineq*dx];
    end
    if ~isempty(conederiv)
        dg = [dg;conederiv];
    end
end

if model.nonlinearequalities
    dgeq = dgeq';
end
if model.nonlinearinequalities | any(model.K.q)
    dg = dg';
end

function conederiv = computeConeDeriv(model,z,dzdx)
conederiv = [];
z = z(:);
if any(model.K.q)
    top = 1 + model.K.f + model.K.l;
    for i = 1:length(model.K.q)
        d = model.F_struc(top,1);
        c = model.F_struc(top,2:end)';
        b = model.F_struc(top+1:top+model.K.q(i)-1,1);
        A = model.F_struc(top+1:top+model.K.q(i)-1,2:end);

        if nargin == 2
            % No inner derivative
            A = A(:,model.linearindicies);
            c = c(model.linearindicies);
            % -c'*x - d + ||Ax+b||>=0
            e = A*z(model.linearindicies) + b;
            smoothed = sqrt(10^-10 + e'*e);
            conederiv = [conederiv;-c'+(A'*b + A'*A*z(model.linearindicies))'/smoothed];
           % conederiv = [conederiv;(2*A(:,model.linearindicies)'*(A(:,model.linearindicies)*z(model.linearindicies)+b)-2*c(model.linearindicies)*(c(model.linearindicies)'*z(model.linearindicies)+d))'];
        else

            % -c'*x - d + ||Ax+b||>=0
            e = A*z + b;
            smoothed = sqrt(10^-10 + e'*e);
            % conederiv = [conederiv;-c'+(dzdx'*A'*b + dzdx'*A'*A*z(model.linearindicies))'/smoothed];
            aux = z'*(A'*A-c*c')*dzdx+(b'*A-d*c')*dzdx;
            conederiv = [conederiv;(-dzdx'*c+aux'/smoothed)'];

            % inner derivative
            % aux = 2*z'*(A'*A-c*c')*dzdx+2*(b'*A-d*c')*dzdx;
            % conederiv = [conederiv;aux];

        end
        top = top + model.K.q(i);
    end
end

author: shamilmamedov


text:
function rtk=baserefpos(rtk,opt,obs,nav)

global glc

if opt.basepostype==glc.POSOPT_POS
    rtk.basepos=opt.basepos;
elseif opt.basepostype==glc.POSOPT_SPP
    rtk.basepos=avepos(rtk,opt,obs,nav);
elseif opt.basepostype==glc.POSOPT_RINEX
    if norm(obs.sta.pos)<=0
        error('Have no base position in renix header!!!\n');
    end
    [~,Cne]=xyz2blh(obs.sta.pos);
    dr=Cne'*obs.sta.del;
    rtk.basepos=obs.sta.pos+dr;
else
    error('Unable to obtain base position!!!\n');
end

if norm(rtk.basepos)==0
    error('The coordinates of the base station is zero!!!');
end

return

repo: GINav


text:
%% Define params
s_cam_definitions;

%% Load image for pathway xforms
ni_seg = niftiRead(fullfile(mask_dir,seg_filename));
% Calc xform from camino space to acpc space
xform = ni_seg.qto_xyz;
xform(1:3,1:3) = eye(3);
xform_cam_to_acpc = xform;
xform_cam_to_acpc(1:3,4) = xform_cam_to_acpc(1:3,4)+diag(ni_seg.qto_xyz(1:3,1:3))/2;
xform_img_to_acpc = ni_seg.qto_xyz;
fd_size = ni_seg.dim;

%% Fix image header
%filename = fullfile(camino_dir,'acm','t_mesd_gm_rh_exit_1cm_acm_sc.nii.gz');
filename = fullfile(camino_dir,'fa.nii.gz');
ni = niftiRead(filename);
dtiWriteNiftiWrapper(ni.data, xform_img_to_acpc, filename);

%% Load generated paths
% Load camino path file
tractsFile = fullfile(camino_dir,'tracts_clean.Bfloat');
fg_temp = mtrImportFibers(tractsFile);
fg = dtiNewFiberGroup;
fg.fibers = fg_temp.fibers;
%fg = dtiXformFiberCoords(fg, xform_cam_to_acpc);
mtrExportFibers(fg,fullfile(camino_dir,'tracts.pdb'));

%% Make fiber density map of all occ lobe
tractsFile = fullfile(camino_dir,'tracts_occ.Bfloat');
fg_temp = mtrImportFibers(tractsFile,eye(4));
fg = dtiNewFiberGroup;
fg.fibers = fg_temp.fibers;
fg = dtiXformFiberCoords(fg, xform_cam_to_acpc);
fd_img = dtiComputeFiberDensityNoGUI(fg, xform_img_to_acpc, fd_size);
dtiWriteNiftiWrapper(fd_img, xform_img_to_acpc, fullfile(camino_dir,'tracts_occ.nii.gz'));
% Intersect occ map with cortex
tracts_filename = 'tracts_occ.nii.gz';
ni_tracts = niftiRead(fullfile(camino_dir, tracts_filename));
ni_c = niftiRead(fullfile(mask_dir,'gm.nii.gz'));
ni_tracts.fname = fullfile(camino_dir,'tracts_occ_gm.nii.gz');
ni_tracts.data(ni_c.data==0) = 0;
writeFileNifti(ni_tracts);
% Intersect with retinotopic cortex
ni_tracts = niftiRead(fullfile(camino_dir, tracts_filename));
ni_c = niftiRead(fullfile(mask_dir,'atlas_mask.nii.gz'));
ni_tracts.fname = fullfile(camino_dir,'tracts_occ_atlas.nii.gz');
ni_tracts.data(ni_c.data==0) = 0;
writeFileNifti(ni_tracts);

%% Make fiber density map of exiting occ lobe
tractsFile = fullfile(camino_dir,'tracts_occ_exit.Bfloat');
fg_temp = mtrImportFibers(tractsFile,eye(4));
fg = dtiNewFiberGroup;
fg.fibers = fg_temp.fibers;
fg = dtiXformFiberCoords(fg, xform_cam_to_acpc);
fd_img = dtiComputeFiberDensityNoGUI(fg, xform_img_to_acpc, fd_size);
dtiWriteNiftiWrapper(fd_img, xform_img_to_acpc, fullfile(camino_dir,'tracts_occ_exit.nii.gz'));
% Intersect occ map with cortex
tracts_filename = 'tracts_occ_exit.nii.gz';
ni_tracts = niftiRead(fullfile(camino_dir, tracts_filename));
ni_c = niftiRead(fullfile(mask_dir,'gm.nii.gz'));
ni_tracts.fname = fullfile(camino_dir,'tracts_occ_exit_gm.nii.gz');
ni_tracts.data(ni_c.data==0) = 0;
writeFileNifti(ni_tracts);
% Intersect with retinotopic cortex
ni_tracts = niftiRead(fullfile(camino_dir, tracts_filename));
ni_c = niftiRead(fullfile(mask_dir,'atlas_mask.nii.gz'));
ni_tracts.fname = fullfile(camino_dir,'tracts_occ_exit_atlas.nii.gz');
ni_tracts.data(ni_c.data==0) = 0;
writeFileNifti(ni_tracts);

%% Tract surface measurements

%% Fix residual image
ni_b0 = niftiRead(fullfile(dti_dir,'b0.nii.gz'));
ni = niftiRead(fullfile(camino_dir,'wdt_rmsqr.nii.gz'));
ni_b0.fname =  fullfile(camino_dir,'wdt_rmsqr.nii.gz');
ni_b0.data = ni.data;
%ni_b0.data = ni.data(1:-1:end,:,:);
writeFileNifti(ni_b0);

%% Examine residuals
ni_mse = niftiRead(fullfile(camino_dir,'wdt_rmsqr.nii.gz'));
ni_std = niftiRead(fullfile(camino_dir,'wdt_std.nii.gz'));
ni_std.data(ni_std.data>100) = 100;
ni_mse.data(ni_mse.data>100) = 100;
showMontage(ni_std.data)
showMontage(ni_mse.data)
e_ratio = zeros(ni_std.dim);
e_ratio(ni_std.data>0) = ni_mse.data(ni_std.data>0) ./ ni_std.data(ni_std.data>0);
showMontage(e_ratio)

% Fix label image
filename_labels = fullfile(camino_dir,'cbs_labels_1_1_0.nii.gz');
ni = niftiRead(filename_labels);
ni_labels = ni_occ;
ni_labels.fname = filename_labels;
ni_labels.data = ni.data;
writeFileNifti(ni_labels);

filename_labels = fullfile(camino_dir,'cbs_labelcp_1_1_0.nii.gz');
ni = niftiRead(filename_labels);
ni_labels = ni_occ;
ni_labels.fname = filename_labels;
ni_labels.data = ni.data;
writeFileNifti(ni_labels);


text:
function clear_class(class_name)
  % CLEAR_CLASS  clears all variables of a specific class
  %
  % clear_class(class_name)
  %
  % Inputs:
  %   class_name  name of class to clear as string
  %
  % See also: clear
  %

  % Get a cell list of all variable names as strings
  s_v = who;
  s_v
  % loop over variabls
  for ii = 1:numel(s_v)
    % Get this variable as string
    s_s = s_v{ii};
    s_s
    s_d = eval(s_s);
    class(s_d)
    if strcmp(class(s_d),class_name)
      fprintf('clearing: %s...',s_s);
      clear(s_s);
    end
  end

end

author: alecjacobson

repo: gptoolbox


text:
%compute speed in the central-head direction
function [data,units]=compute_velch(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velch=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    % just a little faster
    velch{i} = trx(larva).vx_cm.*cos(trx(larva).centralheadang(1,1:end-1)) + ...
      trx(larva).vy_cm.*sin(trx(larva).centralheadang(1,1:end-1));
    % velch{1,i}=trx(larva).velmag_ctr.*(cos(trx(larva).velang).*cos(trx(larva).centralheadang(1,1:end-1))+sin(trx(larva).velang).*sin(trx(larva).centralheadang(1,1:end-1)));
end

units=parseunits('mm/s');
data=velch;

author: kristinbranson

repo: JAABA


text:
function [success] = viewAcceptedRejectedCellExtraction(inputImages,inputSignals,valid,inputMovie,varargin)
	% Plots cell image, peak transients, and image from activity in movie for accepted and rejected cells, useful for publications.
	% Biafra Ahanonu
	% started: 2014.01.03 [19:13:01]
	% inputs
		%
	% outputs
		%

	[success] = ciapkg.view.viewAcceptedRejectedCellExtraction(inputImages,inputSignals,valid,inputMovie,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function varargout = chopdeal(data,sizes)
% Chop a vector into multiple vectors of varying sizes.
% Outputs... = chopdeal(Data,Sizes)
%
% This function operates simiarly to mat2cell, except that it accepts only vectors and returns
% the results as multiple outputs rather than in single cell array. This makes the function
% equivalent to celldeal(mat2cell(Data,1,Sizes)).
%
% In:
%   Data : a vector to chop (single, double, or cell array).
%
%   Sizes : a double vector of sizes for the chop outputs;
%           the sizes must sum to the length of the data.
%
% Out:
%   Outputs... : one output for each size, each of size [1xSizes(k)].
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-01

% the below implementation is the MATLAB fallback for the mex function of same name.
if nargout ~= length(sizes)
    error('The number of output arguments must match the number of elements in sizes.'); end
varargout = cell(1,length(sizes));
p=0;
for s=1:length(sizes)
    varargout{s} = data(p+(1:sizes(s)));
    p = p+sizes(s);
end

repo: REST


text:
function v = niftiCompareRoiPair(roi1,roi1LabelVal,roi2,roi2LabelVal)
%
% function v = niftiCompareRoiPair(roi1,roi1LabelVal,roi2,roi2LabelVal)
%
% This function takes two nifti ROIs and computes volume and a measure of
% overlap (percent agreement) as a means of evaluating reliabitliy in the
% case of two raters drawing the same ROI. The roi.pixdim field is used to
% compute the volume appropriately.
%
% EXAMPLE USAGE:
% v = niftiCompareRoiPair('r1c.nii.gz',1,'r2c.nii.gz',3);
% v =
%         roi1Name: 'r1c.nii.gz'
%          roi1Val: 1
%         roi2Name: 'r2c.nii.gz'
%          roi2Val: 3
%            units: 'mm^3'
%              dim: [1 1 1]
%       roi1Volume: 5265
%       roi2Volume: 4996
%       overlapVol: 4868
%       overlayVol: 5393
%     percentAgree: 90.2652%
%
% HISTORY
% 2011.06.01 LMP: wrote it.
% 2011.06.03 LMP: now calls for a label value for each roi. This allows
%            files with multiple labels to be used. User must specify the
%            value of each. Now uses the dim of the roi itself to calculate
%            the volume.

%%

% Check input arguments and prompt for them if they don't exist
if(~exist('roi1','var') || isempty(roi1))
  dd = pwd; roi1 = mrvSelectFile([],'*.nii.gz','Select ROI 1 File',dd);
end

if(~exist('roi1LabelVal','var')) || isempty(roi1LabelVal)
    roi1LabelVal = str2double(inputdlg('Enter ROI1 Label Value','ROI 1 Label Value'));
end

if(~exist('roi2','var') || isempty(roi2))
  dd = pwd; roi2 = mrvSelectFile([],'*.nii.gz','Select ROI 2 File',dd);
end

if(~exist('roi2LabelVal','var')) || isempty(roi2LabelVal)
    roi2LabelVal = str2double(inputdlg('Enter ROI2 Label Value','ROI 2 Label Value'));
end

% Read in the nifti files
roi1 = niftiRead(roi1);
roi2 = niftiRead(roi2);
if roi1.pixdim==roi1.pixdim
    mmPerVoxel = roi1.pixdim; % Dimensions of the t1 image
else error('ROI dimensions must agree');
end

% Find the points in the data array that make up the ROI
r1Inds = find(roi1.data==roi1LabelVal); % points in roi1 where there is any label >0
r2Inds = find(roi2.data==roi2LabelVal); % points in roi2 where there is any label >0
a = r1Inds(:,1);
b = r2Inds(:,1);
c = ismember(a,b); % Find the overlap indices

% Get the inds of the overlay ROI - that is the ROI that would be created
% if you overlay and combine the two ROIs together (with the overlap only
% being counted once)
d = vertcat(a,b);
e = unique(d); % Overlay indices

% Populate the structure
v.roi1Name = roi1.fname;
v.roi1Val  = roi1LabelVal;
v.roi2Name = roi2.fname;
v.roi2Val  = roi2LabelVal;
v.units = 'mm^3';
v.dim = roi1.pixdim;

% Calculate the volume and percent agreement
v.roi1Volume = length(a)*prod(mmPerVoxel);
v.roi2Volume = length(b)*prod(mmPerVoxel);
v.overlapVol = sum(c)*prod(mmPerVoxel);
v.overlayVol = length(e)*prod(mmPerVoxel);
v.percentAgree = v.overlapVol/v.overlayVol*100;

% Print the results to the screen
fprintf('\n================================\n');
% fprintf('T1 Image        = %s\n',v.imgName);
fprintf('ROI1 Name       = %s\n',v.roi1Name);
fprintf('ROI1 Label Val  = %s\n',num2str(roi1LabelVal));
fprintf('ROI2 Name       = %s\n',v.roi2Name);
fprintf('ROI2 Label Val  = %s\n',num2str(roi2LabelVal));
fprintf('Units           = %s\n',v.units);
fprintf('Voxel Dim       = %s\n',num2str(v.dim));
fprintf('Roi1 Volume     = %s\n',num2str(v.roi1Volume));
fprintf('Roi2 Volume     = %s\n',num2str(v.roi2Volume));
fprintf('Overlap Volume  = %s\n',num2str(v.overlapVol));
fprintf('Overlay Volume  = %s\n',num2str(v.overlayVol));
fprintf('Percent Agree   = %s\n',num2str(v.percentAgree));
fprintf('================================\n');

return

author: vistalab


text:
tic;
neuron.options.nk = 1; %round(T/(60*neuron.Fs)); % number of knots for spline basis, the interval between knots is 180 seconds
nr_patch = str2double(get(edit_nr, 'string'));
nc_patch = str2double(get(edit_nc, 'string'));

patch_par = round([nr_patch, nc_patch]); %1;  % divide the optical field into 3 X 3 patches and do initialization patch by patch
K = str2double(get(edit_K, 'string'));  % maximum number of neurons to search within each patch. you can use [] to search the number automatically
neuron.options.bd = 0; % boundaries to be removed due to motion correction
[center, Cn, pnr] = neuron.initComponents_endoscope(Y, K, patch_par, debug_on, save_avi); %#ok<ASGLU>
figure;
imagesc(Cn);
hold on; plot(center(:, 2), center(:, 1), 'or');
colormap; axis off tight equal;

[~, srt] = sort(max(neuron.C, [], 2)./get_noise_fft(neuron.C), 'descend');
neuron.orderROIs(srt);
neuron_init = neuron.copy();

author: zhoupc


text:
classdef SOP_F16 < PROBLEM
% <single> <real> <expensive/none>
% Six-hump camel-back function

%------------------------------- Reference --------------------------------
% X. Yao, Y. Liu, and G. Lin, Evolutionary programming made faster, IEEE
% Transactions on Evolutionary Computation, 1999, 3(2): 82-102.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 1;
            obj.D = 2;
            obj.lower    = zeros(1,obj.D) - 5;
            obj.upper    = zeros(1,obj.D) + 5;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = 4*PopDec(:,1).^2-2.1*PopDec(:,1).^4+PopDec(:,1).^6/3+PopDec(:,1).*PopDec(:,2)-4*PopDec(:,2).^2+4*PopDec(:,2).^4;
        end
        %% Generate the minimum objective value
        function R = GetOptimum(obj,N)
            R = -1.0316;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function out=prox_sum_k_largest_abs(x,k,alpha)
%PROX_SUM_K_LARGEST_ABS computes the proximal operator of the function
%                     alpha*(sum of k largest absolute values of x(:))
%
%  Usage:
%  out = PROX_SUM_K_LARGEST_ABS(x,k,alpha)
%  ===========================================
%  INPUT:
%  x - point to be projected (vector/matrix)
%  k - positive integer
%  alpha - positive scalar
%  ===========================================
%  Assumptions:
%  k in {1,...,length(x(:))}
%  ===========================================
%  Output:
%  out - proximal operator at x

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

if (nargin < 3)
    error ('usage: prox_sum_k_largest_abs(x,k,alpha)') ;
end

if (alpha < 0)
    error('usage: prox_sum_k_largest_abs(x,k,alpha) - alpha should be positive')
end

eps = 1e-10 ; % default value
if ((k < 1) || ( k > length(x(:))) || (abs(round(k) - k)  > eps))
    error('usage: prox_sum_k_largest_abs(x,k,alpha) - k should be in {1,...,length(x(:))}')
end

out = x - alpha * proj_l1ball_box (x/alpha,ones(size(x)),k,ones(size(x))) ;

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function dat = N_class_loss(algo,dat)

  x = dat.X;
  y = dat.X;
%   [1:length(x)]';

  lss = x-y;
  lss(lss ~= 0) = 1;
  lss = sum(lss)/length(x);

  dat=data([get_name(dat) ' -> N_class_loss=' num2str(lss,4) ],[],lss);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function rbfhKernDisplay(kern, spacing)

% RBFHKERNDISPLAY Display parameters of the RBFH kernel.
% FORMAT
% DESC displays the parameters of the radial basis function heat
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : rbfhKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('RBFH inverse width time: %2.4f (length scale %2.4f)\n', ...
        kern.inverseWidthTime, 1/sqrt(kern.inverseWidthTime));
fprintf(spacing);
fprintf('RBFH inverse width space: %2.4f (length scale %2.4f)\n', ...
        kern.inverseWidthSpace, 1/sqrt(kern.inverseWidthSpace));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function colsum = r8col_sum ( m, n, a )

%*****************************************************************************80
%
%% R8COL_SUM sums the columns of an R8COL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the array to be examined.
%
%    Output, real COLSUM(N), the sums of the columns.
%
  for j = 1 : n
    colsum(j) = sum ( a(1:m,j) );
  end

  return
end


text:
function [y] = spm_gx_hdm(x,u,P,M)
% Simulated BOLD response to input.
% FORMAT [y] = spm_gx_hdm(x,u,P,M)
% y    - BOLD response (%)
% x    - state vector     (see spm_fx_fmri)
% P    - Parameter vector (see spm_fx_fmri)
%__________________________________________________________________________
%
% This function implements the BOLD signal model described in:
%
% Stephan KE, Weiskopf N, Drysdale PM, Robinson PA, Friston KJ (2007)
% Comparing hemodynamic models with DCM. NeuroImage 38: 387-401.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston & Klaas Enno Stephan
% $Id: spm_gx_hdm.m 6856 2016-08-10 17:55:05Z karl $

% biophysical constants for 1.5 T:
%==========================================================================

% hemodynamic parameters
%--------------------------------------------------------------------------
%   H(1) - signal decay                                   d(ds/dt)/ds)
%   H(2) - autoregulation                                 d(ds/dt)/df)
%   H(3) - transit time                                   (t0)
%   H(4) - exponent for Fout(v)                           (alpha)
%   H(5) - resting oxygen extraction                      (E0)
%   H(6) - ratio of intra- to extra-vascular components   (epsilon)
%--------------------------------------------------------------------------
if isstruct(P)
    H     = [0.64 0.32 2.00 0.32 0.4];
    for i = 1:numel(P.decay)
        H(6)   = P.epsilon;
        y(i,1) = spm_gx_hdm(x(i,:),u(i),H,M);
    end
    return
end

% echo time (seconds)
%--------------------------------------------------------------------------
try
    TE = M(1).TE;
catch
    TE = 0.04;
end

% resting venous volume
%--------------------------------------------------------------------------
V0    = 100*0.08;

% slope r0 of intravascular relaxation rate R_iv as a function of oxygen
% saturation Y:  R_iv = r0*[(1-Y)-(1-Y0)]
%--------------------------------------------------------------------------
r0    = 25;

% frequency offset at the outer surface of magnetized vessels
%--------------------------------------------------------------------------
nu0   = 40.3;

% region-specific resting oxygen extraction fractions
%--------------------------------------------------------------------------
E0    = P(5);

% region-specific ratios of intra- to extravascular components of
% the gradient echo signal (prior mean = 1, log-normally distributed
% scaling factor)
%--------------------------------------------------------------------------
epsi  = exp(P(6));

% coefficients in BOLD signal model
%--------------------------------------------------------------------------
k1    = 4.3.*nu0.*E0.*TE;
k2    = epsi.*r0.*E0.*TE;
k3    = 1 - epsi;

% exponentiation of hemodynamic state variables
%--------------------------------------------------------------------------
x     = exp(x);

% BOLD signal
%--------------------------------------------------------------------------
v     = x(3);
q     = x(4);
y     = V0*(k1.*(1 - q) + k2.*(1 - (q./v)) + k3.*(1 - v));

author: spm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compare batch opt of Hawkes with online/stochastic opt of Hawkes
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
% data simulation
options.N = 50; % the number of sequences
options.Nmax = 1000; % the maximum number of events per sequence
options.Tmax = 100; % the maximum size of time window
options.tstep = 0.2;% the step length for computing sup intensity
options.M = 50; % the number of steps
options.GenerationNum = 5; % the number of generations
D = 10; % the dimension of Hawkes processes

nTest = 5;
Nout = 40;

disp('Fast simulation of Hawkes processes with exponential kernel')
para.mu = rand(D,1)/D;
para.A = rand(D, D);
para.A = 0.65 * para.A./max(abs(eig(para.A)));
para.A = reshape(para.A, [D, 1, D]);
para.w = 1;
Seqs = SimulationFast_Thinning_ExpHP(para, options);

err1 = zeros(nTest, Nout);
err2 = zeros(nTest, Nout);

for n = 1:nTest

% initialize
model.A = rand(D,1,D)./(D^2);
model.mu = rand(D,1)./D;
model.kernel = 'exp';
model.w = 1;
model.landmark = 0;

disp('Learning HP by batch opt')
alg1.LowRank = 0;
alg1.Sparse = 0;
alg1.GroupSparse = 0;
alg1.outer = Nout;
alg1.rho = 0.1;
alg1.inner = 1;
alg1.thres = 1e-5;
alg1.Tmax = [];
alg1.storeErr = 1;
alg1.storeLL = 0;
alg1.truth = para;
model1 = Learning_MLE_Basis( Seqs, model, alg1 );

disp('Learning HP by stochastic opt')
alg2.LowRank = 0;
alg2.Sparse = 0;
alg2.GroupSparse = 0;
alg2.epoch = Nout;
alg2.rho = 0.1;
alg2.eventbatch = 20;
alg2.seqbatch = 10;
alg2.historyL = 20;
alg2.thres = 1e-5;
alg2.Tmax = [];
alg2.storeErr = 1;
alg2.storeLL = 0;
alg2.truth = para;
model2 = Learning_MLE_Basis_Stoc( Seqs, model, alg2 );

err1(n,:) = model1.err(:,3)';
err2(n,:) = model2.err(:,3)';

end

em1 = mean(err1);
ev1 = std(err1);
em2 = mean(err2);
ev2 = std(err2);

figure
hold on
errorbar(em1, ev1)
errorbar(em2, ev2)
%plot(1:Nout, model2.err(:,3), 'r-', 1:Nout, model1.err(:,3), 'b-');
legend('Batch HP', 'Stochastic HP')

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function [x,fval,exitflag,info] = mosekqcqp(H,f,A,rl,ru,Q,l,qrl,qru,lb,ub,x0,opts)
%MOSEKQCQP Solve a QCQP using MOSEK
%
%   min 0.5*x'*H*x + f'*x      subject to:     rl <= A*x <= ru
%    x                                         qrl <= x'Qx + l'x <= qru
%                                              lb <= x <= ub
%
%
%   x = mosekqcqp(H,f,A,rl,ru,Q,l,qrl,qru,lb,ub) solves a QCQP where H and f
%   are the objective matrix and vector respectively, A,rl,ru are the
%   linear constraints, Q,l,qrl,qru are the quadratic constraints (see below)
%   and lb,ub are the bounds.
%
%   x = mosekqcqp(H,...,ub,x0) uses x0 as the initial solution guess.
%
%   x = mosekqcqp(H,...,x0,opts) uses opts to pass mosekset options to the
%   MOSEK solver.
%
%   [x,fval,exitflag,info] = mosekqcqp(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   Quadratic Constraint Form:
%       - Single Quadratic Constraint
%           Q is matrix
%           l is a column vector
%           qrl is a scalar
%           qru is a scalar
%       - Multiple Quadratic Constraints
%           Q is a row cell array, each cell is a quadratic constraint Q
%           l is a matrix, each column for each quadratic constraint l
%           qrl is a column vector, each row for each quadratic constraint qrl
%           qru is a column vector, each row for each quadratic constraint qru

%   This function is based in parts on examples from the MOSEK Toolbox,
%   Copyright (c) 1998-2011 MOSEK ApS, Denmark.

%   Copyright (C) 2012 Jonathan Currie (IPL)

t = tic;

% Handle missing arguments
if nargin < 13, opts = mosekset('warnings','off'); else opts = mosekset(opts); end
if nargin < 12, x0 = []; end
if nargin < 11, ub = []; end
if nargin < 10, lb = []; end
if nargin < 9, qru = []; end
if nargin < 8, qrl = []; end
if nargin < 7, l = []; end
if nargin < 6, Q = []; end
if nargin < 5, error('You must supply at least 5 arguments to mosekqcqp'); end

%Build MOSEK Command
[cmd,param] = mosekBuild(opts);

%Create Problem Structure
prob = mosekProb(H,f,A,rl,ru,Q,l,qrl,qru,lb,ub,[],[],x0,opts);

%Call MOSEK
[rcode,res] = mosekopt(cmd,prob,param);

%Extract Results
[x,fval,exitflag,info] = mosekRes(prob,rcode,res,t);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function title = p06_title ( )

%*****************************************************************************80
%
%% P06_TITLE returns the title of problem p06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Discontinuous';

  return
end

author: johannesgerer


text:
function [V,mean_discrepancy]=mdp_eval_policy_TD_0(P,R,discount,policy,N)

% mdp_eval_policy_TD_0  Evaluation of the value function, using the TD(0) algorithm
%
% Arguments
% -------------------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA)  = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount  = discount rate in ]0; 1[
%   policy(S) = optimal policy
%   N(optional) = number of iterations to execute, default value: 10000.
%              It is an integer greater than the default value.
% Evaluation --------------------------------------------------------------
%   V(S)   = optimal value function.
%   mean_discrepancy(N/100) = vector of V discrepancy mean over 100 iterations
%             Then the length of this vector for the default value of N is 100.

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

% Arguments checking
if iscell(P); S = size(P{1},1); else S = size(P,1); end;
if (nargin < 4 || nargin > 5)
    disp('------------------------------------------')
    disp('The number of arguments must be 4 or 5')
    disp('------------------------------------------')
elseif (discount <= 0 || discount >= 1)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0,1[')
    disp('--------------------------------------------------------')
elseif size(policy,1)~=S || any(mod(policy,1)) || any(policy<1) || any(policy>S)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: policy must be a (1xS) vector with integer from 1 to S')
    disp('--------------------------------------------------------')
elseif (nargin == 5) && (N < 10000)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: N must be upper than 10000')
    disp('--------------------------------------------------------')
else

    % initialization of optional argument
    if nargin < 5; N = 10000; end;

    % Initializations
    s=randi([1,S]);      % Initial state choice
    V=zeros(S,1);
    mean_discrepancy=[]; % vector of mean V variations
    discrepancy=[];      % vector of V variation

    for n=1:N

        % Reinitialisation of trajectories every 100 transitions
	if (mod(n,100)==0); s = randi([1,S]); end;

        % Select an action: here action of the policy
        a = policy(s);

        % Simulate next state s_new
        p_s_new = rand(1);
        p = 0;
        s_new = 0;
        while ((p <= p_s_new) && (s_new < S))
            s_new = s_new+1;
            if iscell(P)
                p = p+P{a}(s,s_new);
            else
                p = p+P(s,s_new,a);
            end
        end;

        % Update V
        if iscell(R)
            r = R{a}(s,s_new);
        else
            if ndims(R) == 3
                r = R(s,s_new,a);
            else
                r = R(s,a);
            end
        end
        delta = r + discount*V(s_new)-V(s);
        dV = (1/sqrt(n+2))*delta;
        V(s) = V(s)+dV;

        % Update current state
        s = s_new;

        % Memorize mean V variations on each trajectory (100 transitions)
        discrepancy(mod(n,100)+1) = abs(dV);
        if (length(discrepancy)==100)
            mean_discrepancy = [mean_discrepancy mean(discrepancy)];
            discrepancy = [];
        end;
    end;

end;

author: Sable


text:
function lines = strflushleft(lines, offset)
%STRFLUSHLEFT Flush left each line of text.
%
%   STRFLUSHLEFT(LINES, OFFSET) where LINES is a cell array of lines and OFFSET
%   is a non-negative integer, does a left flush of the text in each line and
%   indents it to the specified offset.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 11:12:55 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 2, nargin));
   if nargin < 2
      offset = 0;
   end

   spc = ' ';
   tab = sprintf('\t');
   pad = spc(:,ones(offset,1));
   lines = lines(:);

   for i = 1:length(lines)
      k = find(lines{i} ~= spc & lines{i} ~= tab);
      if isempty(k)
         lines{i} = '';
      else
         lines{i} = [pad lines{i}(k(1):k(end))];
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function BED = calc_BED(paramS,varargin)
%BED = calc_BED(paramS)
%Ref: Comparison Between Mechanistic Radiobiological Modeling Vs.
%Fowler BED Equation in Evaluating Lung Cancer Radiotherapy Outcome
%for a Broad Range of Fractionation, J Jeong et al., AAPM 2017
%-----------------------------------------------------------------------
% INPUTS
% paramS : Parameter dictionary with fields:
%          d  - Fraction size
%          n  - No. fractions
%          T  - No. treatment days
%          alpha
%          abRatio
% Note : For a 3D dose distibution, calculate 3D BED by setting input
%        paramS.frxSize.val = doseArray3M/numFrx
%-----------------------------------------------------------------------
% AI 12/4/17
% AI 07/30/18 Updated to handle 3D dose distibution

%Define constants
Tk = paramS.Tk.val;         %Kick-off time of repopulation (days)
Tp =  paramS.Tp.val;        %Potential tumor doubling time (days)

alpha = paramS.alpha.val;
abRatio = paramS.abRatio.val;  %alpha/beta for tumor
d = paramS.frxSize.val;
n = paramS.numFractions.val;
if isfield(paramS,'treatmentDays')
    txDaysV = paramS.treatmentDays.val;
    %Check for numeric input
    if ~isnumeric(txDaysV)
        txDaysV = str2num(txDaysV);
    end
    %Otherwise assume function specified
    if isempty(txDaysV)
        txDaysV = eval([txDaysV,'(',num2str(n),')']);
    end
    T = txDaysV(end);
else
    % Default: Compute length of treatment assuming one fraction every
    % weekday with weekend breaks.
    T = floor(n/5)*7 + mod(n,5);
end

%Compute BED
BED = n.* d .* (1 + d./abRatio);
if T > Tk
  BED = BED -  log(2) * (T-Tk)/(alpha*Tp);
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function value = index0n ( n, i_min, i, i_max )

%*****************************************************************************80
%
%% INDEX0N indexes an N-dimensional array by columns, with zero base.
%
%  Discussion:
%
%    Entries of the array are indexed starting at entry
%      ( I_MIN(1), I_MIN(2),...,I_MIN(N) ),
%    and increasing the first index up to I_MAX(1),
%    then the second and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of indices.
%
%    Input, integer I_MIN(N), the minimum indices.
%
%    Input, integer I(N), the indices.
%
%    Input, integer I_MAX(N), for maximum indices.
%
%    Output, integer VALUE, the index of element I.
%
  index_min = 0;

  value = ( i(n) - i_min(n) );

  for j = n - 1 : -1 : 1
    value = value * ( i_max(j) + 1 - i_min(j) ) + ( i(j) - i_min(j) );
  end
  value = value + index_min;

  return
end


text:
function kern = rbfard2KernParamInit(kern)

% RBFARD2KERNPARAMINIT RBFARD2 kernel parameter initialisation.
% The automatic relevance determination version of the radial basis
% function kernel (RBFARD2) is a very smooth non-linear kernel and is a
% popular choice for generic use.
%
% k(x_i, x_j) = sigma2 * exp(-1/2 *(x_i - x_j)'*A*(x_i - x_j))
%
% The parameters are sigma2, the process variance (kern.variance), the
% diagonal matrix of input scales (kern.inputScales, constrained to be
% positive).
%
% SEEALSO : rbfKernParamInit
%
% FORMAT
% DESC initialises the automatic relevance determination radial basis function
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Michalis K. Titsias, 2009

% KERN

% This parameter is restricted positive.
kern.variance = 1;
kern.inputScales = 0.999*ones(1, kern.inputDimension);
kern.nParams = 1 + kern.inputDimension;

kern.transforms(1).index = [1:kern.nParams];
kern.transforms(1).type = optimiDefaultConstraint('positive');

kern.isStationary = true;


text:
function res = logical(a)
%LOGICAL      Like Matlab function "logical" for intval
%
%Call
%
%   L = logical(A)
%
%Same functionality as Matlab/logical for intval quantity A
%

% written  12/06/05     S.M. Rump
%

  wng = warning;
  warning off

  if a.complex                % complex input
    res = logical(real(a.mid)) | logical(imag(a.mid)) | logical(a.rad);
  else                        % real input
    res = logical(a.inf) | logical(a.sup);
  end

  warning(wng)

repo: OpenMAS


text:
function fiff_write_complex(fid,kind,data)
%
% fiff_write_complex(fid,kind,data)
%
% Writes a single-precision complex tag to a fif file
%
%     fid           An open fif file descriptor
%     kind          Tag kind
%     data          The data
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.1  2006/09/23 14:43:56  msh
%   Added routines for writing complex and double complex matrices.
%   Added routine for writing double-precision real matrix.
%
%

me='MNE:fiff_write_complex';

if nargin ~= 3
        error(me,'Incorrect number of arguments');
end

FIFFT_COMPLEX_FLOAT=20;
FIFFV_NEXT_SEQ=0;
nel=numel(data);
datasize=2*nel*4;
count = fwrite(fid,int32(kind),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFT_COMPLEX_FLOAT),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(datasize),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFV_NEXT_SEQ),'int32');
if count ~= 1
   error(me,'write failed');
end
for k = 1:nel
   count = fwrite(fid,real(data(k)),'single');
   if count ~= 1
      error(me,'write failed');
   end
   count = fwrite(fid,imag(data(k)),'single');
   if count ~= 1
      error(me,'write failed');
   end
end
return;

author: fieldtrip

repo: fieldtrip


text:
function a = kahan_inverse ( alpha, n )

%*****************************************************************************80
%
%% KAHAN_INVERSE returns the inverse of the Kahan matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar that defines A.  A typical
%    value is 1.2.  The "interesting" range of ALPHA is 0 < ALPHA < PI.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  ci = cos ( alpha );

  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = 1.0;
      elseif ( i == j - 1 )
        a(i,j) = ci;
      elseif ( i < j )
        a(i,j) = ci * ( 1.0 + ci )^(j-i-1);
      else
        a(i,j) = 0.0;
      end

    end
  end
%
%  Scale the columns.
%
  for j = 1 : n
    si = sin ( alpha )^j;
    a(1:n,j) = a(1:n,j) / si;
  end

  return
end


text:
% KEVALS - Returns the number of kernel evaluations
%
% Syntax: kernel_evals = kevals();
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function kernel_evals = kevals();

global kernel_evals;

if (isempty(kernel_evals))
   kernel_evals = 0;
end;


text:
function geometry_test013 ( )

%*****************************************************************************80
%
%% TEST013 tests CIRCLE_LUNE_CENTROID_2D, CIRCLE_SECTOR_CENTROID_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n_test = 12;

  r = 2.0;
  center(1:2,1) = [ 5.0; 3.0 ];
  theta1 = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST013\n' );
  fprintf ( 1, '  CIRCLE_LUNE_CENTROID_2D computes the centroid of a\n' );
  fprintf ( 1, '    circular lune, defined by joining the endpoints\n' );
  fprintf ( 1, '    of a circular arc.\n' );
  fprintf ( 1, '  CIRCLE_SECTOR_CENTROID_2D computes the centroid of a\n' );
  fprintf ( 1, '    circular sector, defined by joining the endpoints\n' );
  fprintf ( 1, '    of a circular arc to the center.\n' );

  circle_imp_print_2d ( r, center, '  The implicit circle:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first angle of our lune and sector is always 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                         Lune                       Sector\n' );
  fprintf ( 1, '  THETA2           X             Y             X             Y\n' );
  fprintf ( 1, '\n' );

  for i = 0 : n_test

    theta2 = i * 2.0 * pi / n_test;

    centroid1(1:2,1) = circle_lune_centroid_2d ( r, center, theta1, theta2 );

    centroid2(1:2,1) = circle_sector_centroid_2d ( r, center, theta1, theta2 );

    fprintf ( 1, '  %12f  %12f  %12f  %12f  %12f\n', ...
      theta2, centroid1(1:2,1), centroid2(1:2,1) );

  end

  return
end


text:
function [ icc, ccc, acc ] = cc_data_read ( prefix, ncc, n )

%*****************************************************************************80
%
%% CC_DATA_READ reads data about a sparse matrix in CC format.
%
%  Discussion:
%
%    Three files are presumed to exist:
%    * prefix_icc.txt contains NCC ICC values;
%    * prefix_ccc.txt contains N+1 CCC values;
%    * prefix_acc.txt contains NCC ACC values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, a common prefix for the filenames.
%
%    Input, integer NCC, the number of CC elements.
%
%    Input, integer N, the number of columns in the matrix.
%
%    Output, integer ICC(NCC), the CC rows.
%
%    Output, integer CCC(N+1), the compressed CC columns.
%
%    Output, real ACC(NCC), the CC values.
%
  filename_icc = strcat ( prefix, '_icc.txt' );
  icc = i4vec_data_read ( filename_icc, ncc );

  filename_ccc = strcat ( prefix, '_ccc.txt' );
  ccc = i4vec_data_read ( filename_ccc, n + 1 );

  filename_acc = strcat ( prefix, '_acc.txt' );
  acc = r8vec_data_read ( filename_acc, ncc );

  return
end