text:
%----------------------------------------------------------------
% Depth Map Super-Resolution by Deep Multi-Scale Guidance, ECCV16
%                     written by T.-W. HUI
%---------------------------------------------------------------
clear all;
close all;
clc;

%% Settings
dataDir = '/path/to/training_data/';
saveDir = '/path/to/h5/';

chunksz = 128;

scale = 2;
size_dataY = 40;
size_dataD = 20;
size_label = 40;
border_dataY = 1;
border_label = 1;
stride = 22;
data_aug = true;
index_val = [1 20 28 58 64 66 69 73 75 79];

%% Initialization
dataY = zeros(size_dataY, size_dataY, 1, 1);
dataD = zeros(size_dataD, size_dataD, 1, 1);
label = zeros(size_label, size_label, 1, 1);

size_subIm = size_dataD*scale;

pad_dataY = (size_subIm - size_dataY)/2;
pad_label = (size_subIm - size_label)/2;
assert(pad_dataY >= 0, 'pad_dataY cannot be an negative number!');
assert(pad_label >= 0, 'pad_label cannot be an negative number!');

count = 0;

%% Generate training  data
load([dataDir '/RGB']);
load([dataDir '/Df']);

index_train = setdiff(1:numel(Df), index_val);
fprintf('%d images is used for training.\n', numel(index_train));

RGB = RGB(index_train);
Df  = Df(index_train);

h = ones(3,3)/9;

for k = 1 : size(Df,1)

    fprintf('Processing %d of %d, id: %d, ', k, size(Df,1), index_train(k));

    im_Dh = modcrop(Df{k}, scale);
    im_M  = im_Dh > 0;
    im_Dh = normalizeIm(im_Dh, [0 1]);
    im_Dl = imresize(im_Dh, 1/scale, 'bicubic');
    im_Dl_LF = imfilter(im_Dl, h, 'symmetric');
    im_label = im_Dh - imresize(im_Dl_LF, scale, 'bicubic');
    [H, W] = size(im_label);

    I = RGB{k};
    I = rgb2ycbcr(I);
    Y = double(I(:,:,1));
    im_Y = modcrop(Y, scale);
    im_Y = normalize_cleanIm(im_Y);
    im_Y_LF = imfilter(im_Y, h, 'symmetric');
    im_inputY = normalize_cleanIm(im_Y - im_Y_LF);

    im_grad_Dh = imgradient(im_Dh, 'sobel');

    stride_ = stride;

    if index_train(k) >= 75
        stride_ = floor(2.0*stride_);
    elseif index_train(k) >= 65
        stride_ = floor(1.5*stride_);
    end

    for x = 1 : stride_ : H-size_subIm+1
        for y = 1 :stride_ : W-size_subIm+1

            subIm_M = im_M(x : x+size_subIm-1, y : y+size_subIm-1);

            if sum(sum(subIm_M == 0)) == 0

                subim_dataY = im_inputY(x+pad_dataY : x+pad_dataY+size_dataY-1, y+pad_dataY : y+pad_dataY+size_dataY-1);

                subim_dataD = imresize(im_Dh(x : x+size_subIm-1, y : y+size_subIm-1), 1/scale, 'bicubic');
                subim_dataD = subim_dataD - imfilter(subim_dataD, h, 'symmetric');

                subim_label = im_label(x+pad_label : x+pad_label+size_label-1, y+pad_label : y+pad_label+size_label-1);

                count = count + 1;
                dataY(:, :, 1, count) = subim_dataY;
                dataD(:, :, 1, count) = subim_dataD;
                label(:, :, 1, count) = subim_label;

            end

        end
    end

end

fprintf('---------------------------------------\n')

fprintf('Total: %d patches.\n', count);
num_patches = count;

% Data augmntation
if data_aug
    num_patches = num_patches*2;
    dataY_tmp = dataY; clear dataY;
    dataD_tmp = dataD; clear dataD;
    label_tmp = label; clear label;
    dataY = cat(4, dataY_tmp, rot90(dataY_tmp,1));
    dataD = cat(4, dataD_tmp, rot90(dataD_tmp,1));
    label = cat(4, label_tmp, rot90(label_tmp,1));
end

fprintf('Number of patches used: %d\n', num_patches);

% Remove left-most column and bottom row of dataY and label
dataY = dataY(1:end-border_dataY,1:end-border_dataY,:,:);
label = label(1:end-border_label,1:end-border_label,:,:);

% Re-arrange order
order = randperm(num_patches);
dataY = dataY(:,:,:,order);
dataD = dataD(:,:,:,order);
label = label(:,:,:,order);

%% writing to HDF5
created_flag = false;
totalct = 0;

for batchno = 1 : floor(num_patches/chunksz)
    last_read = (batchno - 1)*chunksz;
    batchdata = dataY(:, :, :, last_read+1 : last_read+chunksz);

    startloc = struct('dat', [1,1,1,totalct+1], 'lab', [1,1,1,totalct+1]);
    curr_dat_sz = store2hdf5_v2(savepathY, '-Y', batchdata, [], [], ~created_flag, startloc, chunksz);
    created_flag = true;
    totalct = curr_dat_sz(end);
end

h5disp(savepathY);

created_flag = false;
totalct = 0;

for batchno = 1 : floor(num_patches/chunksz)
    last_read = (batchno - 1)*chunksz;
    batchdata = dataD(:, :, :, last_read+1 : last_read+chunksz);
    batchlabs = label(:, :, :, last_read+1 : last_read+chunksz);

    startloc = struct('dat', [1,1,1,totalct+1], 'lab', [1,1,1,totalct+1]);
    curr_dat_sz = store2hdf5_v2(savepathD, '-D', batchdata, [], batchlabs, ~created_flag, startloc, chunksz);
    created_flag = true;
    totalct = curr_dat_sz(end);
end

h5disp(savepathD);

author: twhui

repo: MSG-Net

save_path: github-repos/MATLAB/twhui-MSG-Net


text:
%% OPTI MEX Install Script

% This file enables me to quickly rebuild all solvers - it is not intended
% to be called by the user!

% My build platform:
% - Windows 8 x64
% - Visual Studio 2012
% - Intel Compiler XE (C++ & FORTRAN)
% - Intel Math Kernel Library

% YOU MUST BE IN THE BASE DIRECTORY OF THE OPTI TOOLBOX!

%% -- Intel C++ Solvers --
clc
fprintf(2,'Please Specify Intel C++ as your Compiler with VS2012 Linker (Option 3, Don''t Look for Installed Compilers)...\n\n');
mex -setup
%% CLP
opti_CLP_Install

%% -- VS2012 Solvers --
clc
fprintf(2,'Please Specify Visual Studio 2012 as your Compiler...\n\n');
mex -setup
%% BONMIN
opti_BONMIN_Install

%% CBC
opti_CBC_Install

%% CSDP
opti_CSDP_Install

%% DSDP
opti_DSDP_Install

%% FILTERSD
opti_FILTERSD_Install

%% GLPK
opti_GLPK_Install

%% HYBRJ
opti_HYBRJ_Install

%% IPOPT
opti_IPOPT_Install

%% L-BFGS-B
opti_LBFGSB_Install

%% LEVMAR
opti_LEVMAR_Install

%% LIPSOL
opti_LIPSOL_Install

%% LMDER
opti_LMDER_Install

%% LP_SOLVE
opti_LPSOLVE_Install

%% M1QN3
opti_M1QN3_Install

%% MKL TR NLS
opti_MKLTRNLS_Install

%% MUMPS
opti_MUMPS_Install

%% NL2SOL
opti_NL2SOL_Install

%% NLOPT
opti_NLOPT_Install

%% NOMAD
opti_NOMAD_Install

%% OOQP
opti_OOQP_Install

%% PSWARM
opti_PSWARM_Install

%% QSOPT
opti_QSOPT_Install

%% SCIP
opti_SCIP_Install

%% -- Utilities --
%% MKL JAC
opti_MKLJAC_Install

%% CoinUtils
opti_COINUTILS_Install

%% AMPL
opti_AMPL_Install

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function dcfplot(devkeyname,evalkeyname,devscrfilename,evalscrfilename,outfilename,plot_title,xmin,xmax,ymin,ymax,prior)
% Makes a Norm_DCF plot of the dev and eval scores for a system.
% Inputs:
%   devkeyname: The name of the file containing the Key for
%     the dev scores.
%   evalkeyname: The name of the file containing the Key for
%     the eval scores.
%   devscrfilename: The name of the file containing the Scores
%     for the dev trials.
%   evalscrfilename: The name of the file containing the
%     Scores the eval trials.
%   outfilename: The name for the PDF file that the plot will be
%     written in.
%   plot_title: A string for the plot title. (optional)
%   xmin, xmax, ymin, ymax: The boundaries of the plot. (optional)
%   prior: The effective target prior. (optional)

assert(isa(devkeyname,'char'))
assert(isa(evalkeyname,'char'))
assert(isa(devscrfilename,'char'))
assert(isa(evalscrfilename,'char'))
assert(isa(outfilename,'char'))

if ~exist('plot_title','var') || isempty(plot_title)
    plot_title = '';
end

if ~exist('xmin','var')
    xmin = -10;
    xmax = 0;
    ymin = 0;
    ymax = 1.2;
    prior = 0.001;
end

[dev_tar,dev_non] = get_tar_non_scores(devscrfilename,devkeyname);
[eval_tar,eval_non] = get_tar_non_scores(evalscrfilename,evalkeyname);

close all
plot_obj = Norm_DCF_Plot([xmin,xmax,ymin,ymax],plot_title);
plot_obj.set_system(dev_tar,dev_non,'dev')
plot_obj.plot_operating_point(logit(prior),'m--','new DCF point')
plot_obj.plot_curves([0 0 0 1 1 1 0 0],{{'b--'},{'g--'},{'r--'}})
plot_obj.set_system(eval_tar,eval_non,'eval')
plot_obj.plot_curves([0 0 1 1 1 1 0 1],{{'r','LineWidth',2},{'b'},{'g'},{'r'},{'k*'}})
plot_obj.display_legend()
plot_obj.save_as_pdf(outfilename)
end

function [tar,non] = get_tar_non_scores(scrfilename,keyname)
key = Key.read(keyname);
scr = Scores.read(scrfilename);
[tar,non] = scr.get_tar_non(key);
end

author: nesl


text:
function [dat, opt] = proc_pca(dat, varargin)
%PROC_PCA - Principal Component Analysis
% [dat_pca, pca_opt] = proc_pca(dat, <opt>)
%
%Synopsis for training PCA:
% [DAT_PCA_TRAIN, PCA_OPT] = proc_pca(DAT_TRAIN, <OPT>)
%
%Synopsis for applying PCA:
% DAT_PCA_TEST = proc_pca(DAT_TEST, PCA_OPT)
%
%Arguments:
% DAT_TRAIN     - data structure of continuous or epoched data
% OPT           - struct or property/value list of optional properties:
%  .whitening   - if 1, the output dimensions will all have unit variance
%                 (default 0)
% PCA_OPT       - a struct that contains: a bias vector, filters and field
%                 patterns of the sources
%
%Returns
% DAT_PCA_TRAIN,
% DAT_PCA_TEST  - updated data structure
% PCA_OPT       - a struct that contains: a bias vector, filters and field
%                 patterns of the sources
%

% Sven Daehne, 03.2011, sven.daehne@tu-berlin.de
% Matthias Schultze-Kraft, 12.2015, schultze-kraft@tu-berlin.de

props= {'whitening'         0       'BOOL'
        'filters'           []      'DOUBLE'
        'field_patterns'    []      'DOUBLE'
        'bias'              []      'DOUBLE'};

if nargin==0,
  dat = props; return
end

dat = misc_history(dat);
misc_checkType(dat, 'STRUCT(x clab y)');

opt = opt_proplistToStruct(varargin{:});
opt = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

[T,nChans,nEpos] = size(dat.x);

%% train PCA
if isempty(opt.filters) || isempty(opt.bias)
    % get the data matrix
    if ndims(dat.x)==3
        % since time structure does not matter, we can simply concatenate all
        % epochs to get one big data matrix
        X = permute(dat.x, [1,3,2]); % now channels are the last dimension
        X = reshape(X, [T*nEpos, nChans]);
    else
        X = dat.x;
    end
    % remove the mean here already
    b = mean(X,1);
    B = repmat(b, [T*nEpos, 1]);
    X = X - B;

    % perform PCA and compute whitening matrix
    C = cov(X);
    [V, D] = eig(C);
    [ev_sorted, sort_idx] = sort(diag(D), 'descend');
    V = V(:,sort_idx);
    D = diag(ev_sorted);

    if opt.whitening
        V = V * diag(diag(D).^-0.5);
    end

    opt.filters = V;
    opt.field_patterns = V';
    opt.eigenvalues = ev_sorted;
    opt.bias = b;

end

%% apply PCA
if not(length(opt.bias) == nChans)
    error('Dimension of bias must equal the number of channels!')
end
% make sure opt.bias is a row vector
if size(opt.bias, 1) > size(opt.bias,2)
    opt.bias = opt.bias';
end

% subtract bias, then apply filters
B = squeeze(repmat(opt.bias, [T,1,nEpos]));
dat.x = dat.x - B;
dat = proc_linearDerivation(dat, opt.filters);

author: bbci


text:
function [xo,Nout]=largestr(xi,p,varargin)
%LARGESTR   Keep fixed ratio of largest coefficients
%   Usage:  xo=largestr(x,p);
%           xo=largestr(x,p,mtype);
%           [xo,N]=largestr(...);
%
%   `largestr(x,p)` returns an array of the same size as *x* keeping
%   the fraction *p* of the coefficients. The coefficients with the largest
%   magnitude are kept.
%
%   `[xo,n]=largestr(xi,p)` additionally returns the number of coefficients
%   kept.
%
%   **Note:** If the function is used on coefficients coming from a
%   redundant transform or from a transform where the input signal was
%   padded, the coefficient array will be larger than the original input
%   signal. Therefore, the number of coefficients kept might be higher than
%   expected.
%
%   `largestr` takes the following flags at the end of the line of input
%   arguments:
%
%     'hard'    Perform hard thresholding. This is the default.
%
%     'wiener'  Perform empirical Wiener shrinkage. This is in between
%               soft and hard thresholding.
%
%     'soft'    Perform soft thresholding.
%
%     'full'    Returns the output as a full matrix. This is the default.
%
%     'sparse'  Returns the output as a sparse matrix.
%
%   **Note:** If soft- or Wiener thresholding is selected, one less
%   coefficient will actually be returned. This is caused by that
%   coefficient being set to zero.
%
%   See also:  largestn
%
%   References: ma98

%   AUTHOR : Peter L. Søndergaard
%   TESTING: OK
%   REFERENCE: OK

if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'thresh'};
[flags,keyvals]=ltfatarghelper({},definput,varargin);

if (prod(size(p))~=1 || ~isnumeric(p))
  error('p must be a scalar.');
end;

wascell=iscell(xi);

if wascell
  [xi,shape]=cell2vec(xi);
end;

% Determine the size of the array.
ss=numel(xi);

N=round(ss*p);

[xo,Nout]=largestn(xi,N,flags.outclass,flags.iofun);

if wascell
  xo=vec2cell(xo,shape);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function model = lleReconstruct(lleInfo, y)

% LLERECONSTRUCT Reconstruct an LLE form component parts.
% FORMAT
% DESC takes component parts of an LLE model and reconstructs the
% LLE model. The component parts are normally retrieved from a
% saved file.
% ARG lleInfo : the active set and other information stored in a structure.
% ARG y : the output target training data for the LLE.
% RETURN model : an LLE model structure that combines the component
% parts.
%
% SEEALSO : lleDeconstruct, lleCreate
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS

  model = lleInfo;
  model.Y = y;

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk,dxpdalpha] = project_points2(X,om,T,f,c,k,alpha)

%project_points.m
%
%[xp,dxpdom,dxpdT,dxpdf,dxpdc,dxpdk] = project_points2(X,om,T,f,c,k,alpha)
%
%Projects a 3D structure onto the image plane.
%
%INPUT: X: 3D structure in the world coordinate frame (3xN matrix for N points)
%       (om,T): Rigid motion parameters between world coordinate frame and camera reference frame
%               om: rotation vector (3x1 vector); T: translation vector (3x1 vector)
%       f: camera focal length in units of horizontal and vertical pixel units (2x1 vector)
%       c: principal point location in pixel units (2x1 vector)
%       k: Distortion coefficients (radial and tangential) (4x1 vector)
%       alpha: Skew coefficient between x and y pixel (alpha = 0 <=> square pixels)
%
%OUTPUT: xp: Projected pixel coordinates (2xN matrix for N points)
%        dxpdom: Derivative of xp with respect to om ((2N)x3 matrix)
%        dxpdT: Derivative of xp with respect to T ((2N)x3 matrix)
%        dxpdf: Derivative of xp with respect to f ((2N)x2 matrix if f is 2x1, or (2N)x1 matrix is f is a scalar)
%        dxpdc: Derivative of xp with respect to c ((2N)x2 matrix)
%        dxpdk: Derivative of xp with respect to k ((2N)x4 matrix)
%
%Definitions:
%Let P be a point in 3D of coordinates X in the world reference frame (stored in the matrix X)
%The coordinate vector of P in the camera reference frame is: Xc = R*X + T
%where R is the rotation matrix corresponding to the rotation vector om: R = rodrigues(om);
%call x, y and z the 3 coordinates of Xc: x = Xc(1); y = Xc(2); z = Xc(3);
%The pinehole projection coordinates of P is [a;b] where a=x/z and b=y/z.
%call r^2 = a^2 + b^2.
%The distorted point coordinates are: xd = [xx;yy] where:
%
%xx = a * (1 + kc(1)*r^2 + kc(2)*r^4 + kc(5)*r^6)      +      2*kc(3)*a*b + kc(4)*(r^2 + 2*a^2);
%yy = b * (1 + kc(1)*r^2 + kc(2)*r^4 + kc(5)*r^6)      +      kc(3)*(r^2 + 2*b^2) + 2*kc(4)*a*b;
%
%The left terms correspond to radial distortion (6th degree), the right terms correspond to tangential distortion
%
%Finally, convertion into pixel coordinates: The final pixel coordinates vector xp=[xxp;yyp] where:
%
%xxp = f(1)*(xx + alpha*yy) + c(1)
%yyp = f(2)*yy + c(2)
%
%
%NOTE: About 90 percent of the code takes care fo computing the Jacobian matrices
%
%
%Important function called within that program:
%
%rodrigues.m: Computes the rotation matrix corresponding to a rotation vector
%
%rigid_motion.m: Computes the rigid motion transformation of a given structure

if nargin < 7,
   alpha = 0;
   if nargin < 6,
      k = zeros(5,1);
      if nargin < 5,
         c = zeros(2,1);
         if nargin < 4,
            f = ones(2,1);
            if nargin < 3,
               T = zeros(3,1);
               if nargin < 2,
                  om = zeros(3,1);
                  if nargin < 1,
                     error('Need at least a 3D structure to project (in project_points.m)');
                     return;
                  end;
               end;
            end;
         end;
      end;
   end;
end;

[m,n] = size(X);

[Y,dYdom,dYdT] = rigid_motion(X,om,T);

inv_Z = 1./Y(3,:);

x = (Y(1:2,:) .* (ones(2,1) * inv_Z)) ;

bb = (-x(1,:) .* inv_Z)'*ones(1,3);
cc = (-x(2,:) .* inv_Z)'*ones(1,3);

dxdom = zeros(2*n,3);
dxdom(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(1:3:end,:) + bb .* dYdom(3:3:end,:);
dxdom(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdom(2:3:end,:) + cc .* dYdom(3:3:end,:);

dxdT = zeros(2*n,3);
dxdT(1:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(1:3:end,:) + bb .* dYdT(3:3:end,:);
dxdT(2:2:end,:) = ((inv_Z')*ones(1,3)) .* dYdT(2:3:end,:) + cc .* dYdT(3:3:end,:);

% Add distortion:

r2 = x(1,:).^2 + x(2,:).^2;

dr2dom = 2*((x(1,:)')*ones(1,3)) .* dxdom(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdom(2:2:end,:);
dr2dT = 2*((x(1,:)')*ones(1,3)) .* dxdT(1:2:end,:) + 2*((x(2,:)')*ones(1,3)) .* dxdT(2:2:end,:);

r4 = r2.^2;

dr4dom = 2*((r2')*ones(1,3)) .* dr2dom;
dr4dT = 2*((r2')*ones(1,3)) .* dr2dT;

r6 = r2.^3;

dr6dom = 3*((r2'.^2)*ones(1,3)) .* dr2dom;
dr6dT = 3*((r2'.^2)*ones(1,3)) .* dr2dT;

% Radial distortion:

cdist = 1 + k(1) * r2 + k(2) * r4 + k(5) * r6;

dcdistdom = k(1) * dr2dom + k(2) * dr4dom + k(5) * dr6dom;
dcdistdT = k(1) * dr2dT + k(2) * dr4dT + k(5) * dr6dT;
dcdistdk = [ r2' r4' zeros(n,2) r6'];

xd1 = x .* (ones(2,1)*cdist);

dxd1dom = zeros(2*n,3);
dxd1dom(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdom;
dxd1dom(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdom;
coeff = (reshape([cdist;cdist],2*n,1)*ones(1,3));
dxd1dom = dxd1dom + coeff.* dxdom;

dxd1dT = zeros(2*n,3);
dxd1dT(1:2:end,:) = (x(1,:)'*ones(1,3)) .* dcdistdT;
dxd1dT(2:2:end,:) = (x(2,:)'*ones(1,3)) .* dcdistdT;
dxd1dT = dxd1dT + coeff.* dxdT;

dxd1dk = zeros(2*n,5);
dxd1dk(1:2:end,:) = (x(1,:)'*ones(1,5)) .* dcdistdk;
dxd1dk(2:2:end,:) = (x(2,:)'*ones(1,5)) .* dcdistdk;

% tangential distortion:

a1 = 2.*x(1,:).*x(2,:);
a2 = r2 + 2*x(1,:).^2;
a3 = r2 + 2*x(2,:).^2;

delta_x = [k(3)*a1 + k(4)*a2 ;
   k(3) * a3 + k(4)*a1];

%ddelta_xdx = zeros(2*n,2*n);
aa = (2*k(3)*x(2,:)+6*k(4)*x(1,:))'*ones(1,3);
bb = (2*k(3)*x(1,:)+2*k(4)*x(2,:))'*ones(1,3);
cc = (6*k(3)*x(2,:)+2*k(4)*x(1,:))'*ones(1,3);

ddelta_xdom = zeros(2*n,3);
ddelta_xdom(1:2:end,:) = aa .* dxdom(1:2:end,:) + bb .* dxdom(2:2:end,:);
ddelta_xdom(2:2:end,:) = bb .* dxdom(1:2:end,:) + cc .* dxdom(2:2:end,:);

ddelta_xdT = zeros(2*n,3);
ddelta_xdT(1:2:end,:) = aa .* dxdT(1:2:end,:) + bb .* dxdT(2:2:end,:);
ddelta_xdT(2:2:end,:) = bb .* dxdT(1:2:end,:) + cc .* dxdT(2:2:end,:);

ddelta_xdk = zeros(2*n,5);
ddelta_xdk(1:2:end,3) = a1';
ddelta_xdk(1:2:end,4) = a2';
ddelta_xdk(2:2:end,3) = a3';
ddelta_xdk(2:2:end,4) = a1';

xd2 = xd1 + delta_x;

dxd2dom = dxd1dom + ddelta_xdom ;
dxd2dT = dxd1dT + ddelta_xdT;
dxd2dk = dxd1dk + ddelta_xdk ;

% Add Skew:

xd3 = [xd2(1,:) + alpha*xd2(2,:);xd2(2,:)];

% Compute: dxd3dom, dxd3dT, dxd3dk, dxd3dalpha

dxd3dom = zeros(2*n,3);
dxd3dom(1:2:2*n,:) = dxd2dom(1:2:2*n,:) + alpha*dxd2dom(2:2:2*n,:);
dxd3dom(2:2:2*n,:) = dxd2dom(2:2:2*n,:);
dxd3dT = zeros(2*n,3);
dxd3dT(1:2:2*n,:) = dxd2dT(1:2:2*n,:) + alpha*dxd2dT(2:2:2*n,:);
dxd3dT(2:2:2*n,:) = dxd2dT(2:2:2*n,:);
dxd3dk = zeros(2*n,5);
dxd3dk(1:2:2*n,:) = dxd2dk(1:2:2*n,:) + alpha*dxd2dk(2:2:2*n,:);
dxd3dk(2:2:2*n,:) = dxd2dk(2:2:2*n,:);
dxd3dalpha = zeros(2*n,1);
dxd3dalpha(1:2:2*n,:) = xd2(2,:)';

% Pixel coordinates:
if length(f)>1,
    xp = xd3 .* (f * ones(1,n))  +  c*ones(1,n);
    coeff = reshape(f*ones(1,n),2*n,1);
    dxpdom = (coeff*ones(1,3)) .* dxd3dom;
    dxpdT = (coeff*ones(1,3)) .* dxd3dT;
    dxpdk = (coeff*ones(1,5)) .* dxd3dk;
    dxpdalpha = (coeff) .* dxd3dalpha;
    dxpdf = zeros(2*n,2);
    dxpdf(1:2:end,1) = xd3(1,:)';
    dxpdf(2:2:end,2) = xd3(2,:)';
else
    xp = f * xd3 + c*ones(1,n);
    dxpdom = f  * dxd3dom;
    dxpdT = f * dxd3dT;
    dxpdk = f  * dxd3dk;
    dxpdalpha = f .* dxd3dalpha;
    dxpdf = xd3(:);
end;

dxpdc = zeros(2*n,2);
dxpdc(1:2:end,1) = ones(n,1);
dxpdc(2:2:end,2) = ones(n,1);

return;

% Test of the Jacobians:

n = 10;

X = 10*randn(3,n);
om = randn(3,1);
T = [10*randn(2,1);40];
f = 1000*rand(2,1);
c = 1000*randn(2,1);
k = 0.5*randn(5,1);
alpha = 0.01*randn(1,1);

[x,dxdom,dxdT,dxdf,dxdc,dxdk,dxdalpha] = project_points2(X,om,T,f,c,k,alpha);

% Test on om: OK

dom = 0.000000001 * norm(om)*randn(3,1);
om2 = om + dom;

[x2] = project_points2(X,om2,T,f,c,k,alpha);

x_pred = x + reshape(dxdom * dom,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on T: OK!!

dT = 0.0001 * norm(T)*randn(3,1);
T2 = T + dT;

[x2] = project_points2(X,om,T2,f,c,k,alpha);

x_pred = x + reshape(dxdT * dT,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on f: OK!!

df = 0.001 * norm(f)*randn(2,1);
f2 = f + df;

[x2] = project_points2(X,om,T,f2,c,k,alpha);

x_pred = x + reshape(dxdf * df,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on c: OK!!

dc = 0.01 * norm(c)*randn(2,1);
c2 = c + dc;

[x2] = project_points2(X,om,T,f,c2,k,alpha);

x_pred = x + reshape(dxdc * dc,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on k: OK!!

dk = 0.001 * norm(k)*randn(5,1);
k2 = k + dk;

[x2] = project_points2(X,om,T,f,c,k2,alpha);

x_pred = x + reshape(dxdk * dk,2,n);

norm(x2-x)/norm(x2 - x_pred)

% Test on alpha: OK!!

dalpha = 0.001 * norm(k)*randn(1,1);
alpha2 = alpha + dalpha;

[x2] = project_points2(X,om,T,f,c,k,alpha2);

x_pred = x + reshape(dxdalpha * dalpha,2,n);

norm(x2-x)/norm(x2 - x_pred)

author: strawlab

repo: MultiCamSelfCal


text:
function res = imRescale(img, varargin)
%IMRESCALE Rescale gray levels of image to get better dynamic
%
%   This function has been replaced by function 'imAdjustDynamic', and is
%   now deprecated.
%
%   RES = imRescale(IMG, [GMIN GMAX]);
%   all values below GMIN will be set to 0, all values greater than GMAX
%   will be set to 255, all values in between will be equally spaced
%   between 0 and 255.
%   The result is a 255 grayscale image.
%
%   RES = imRescale(IMG);
%   rescale using min and max values found in image.
%
%   See Also:
%   imLUT, imGrayscaleExtent, imadjust, mat2gray
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/03/2005.
%

%   HISTORY
%   25/03/2005 correct bug, due to computation with unsigned values
%   21/08/2009 extends to manage double images
%   05/11/2011 deprecated and replace by imAdjustDynamic

warning('matImage:deprecatedFunction', ...
    'function "imRescale" has been deprecated and replaced by "imAdjustDynamic"');

% process input arguments
if ~isempty(varargin)
    % use min and max values given as parameter
    var = varargin{1};
    min1 = double(var(1));
    max1 = double(var(end));
else
    % use min and max values computed from input image
    min1 = double(min(img(:)));
    max1 = double(max(img(:)));
end

% compute slope of linear transformation
a = double(255 / (max1 - min1));

% compute result image
% values below 0 or greater than 255 are automatically clipped when
% casting to uint8
res = uint8((img - min1) * a);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function val = IRget(options,name,value,flag)
%IRget Get options for IR Tools functions
%
% val = IRget(options,'name',value,flag)
%
% val = IRget(options,'name')
%   extracts the value of the named parameter from the options structure,
%   returning an empty matrix if the parameter value is not specified in
%   options.  It is sufficient to type only the leading characters that
%   uniquely identify the parameter, and case is ignored for parameter names.
%   [] is a valid options argument.
%
% val = IRget(options,'name',value)
%   extracts the named parameter as above, but returns VALUE if the named
%   parameter is not specified in options.
%
% val = IRget(options,'name',value,flag)
%   with flag = 'fast' no error checking is performed.
%
% See also: IRset, PRget, PRset

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD Licence. A separate license file should be provided as part
% of the package.

% Fast access with no error checking.
if (nargin == 4) && isequal(flag,'fast')
    val = IRgetfast(options,name,value);
    return
end

if nargin < 2
    error('Not enough input arguments.');
end
if nargin < 3
    value = [];
end

if ~isempty(options) && ~isa(options,'struct')
    error('First argument must be an options structure created with IRset.');
end

if isempty(options)
    val = value;
    return;
end
allfields = IRoption_fields;

Names = allfields;

name = deblank(name(:)'); % Force this to be a row vector.
j = find(strncmpi(name,Names,length(name)));
if isempty(j)             % If no matches.
    error(['Unrecognized property name ''%s''.  ' ...
           'See IRset for possibilities.'], name);
elseif length(j) > 1      % If more than one match.
    % Check for any exact matches (in case any names are subsets of others).
    k = find(strcmpi(name,Names));
    if length(k) == 1
        j = k;
    else
        error('Ambiguous property name ''%s'' ', name);
    end
end

if any(strcmp(Names,Names{j,:}))
    val = options.(Names{j,:});
    if isempty(val)
        val = value;
    end
else
    val = value;
end

%------------------------------------------------------------------
function value = IRgetfast(options,name,defaultopt)
%IRGETFAST- Get IR OPTIONS parameter with no error checking.
%   VAL = IRGETFAST(OPTIONS,FIELDNAME,DEFAULTOPTIONS) will get the
%   value of the FIELDNAME from OPTIONS with no error checking or
%   fieldname completion. If the value is [], it gets the value of the
%   FIELDNAME from DEFAULTOPTIONS, another OPTIONS structure which is
%   probably a subset of the options in OPTIONS.

if isempty(options)
     value = defaultopt.(name);
     return;
end
% We need to know if name is a valid field of options, but it is faster
% to use a try-catch than to test if the field exists and if the field name
% is correct.
try
    value = options.(name);
catch
    value = [];
end

if isempty(value)
    value = defaultopt.(name);
end

repo: IRtools


text:
classdef GradientVariationWithRadiusExperiment < handle

    properties (Access = private)
       iMesh
       gExperiment
    end

    properties (Access = private)
       nameCase
       inputFiles
       outputFolder
       levelSetParams
    end

    methods (Access = public)

        function obj = GradientVariationWithRadiusExperiment(s)
            obj.init(s);
        end

        function compute(obj)
            for im = 1:numel(obj.inputFiles)
                obj.iMesh = im;
                obj.createGradientVariationExperiment();
                obj.computeGradientVariationWithRadius();
            end
        end
    end

    methods (Access = private)

        function init(obj,cParams)
            obj.nameCase       = cParams.nameCase;
            obj.inputFiles     = cParams.inputFiles;
            obj.outputFolder   = cParams.outputFolder;
            obj.levelSetParams = cParams.levelSetParams;
        end

        function createGradientVariationExperiment(obj)
            s.inputFile = obj.inputFiles{obj.iMesh};
            s.iMesh     = obj.iMesh;
            s.levelSetParams = obj.levelSetParams;
            g = GradientVariationExperiment(s);
            obj.gExperiment = g;
        end

        function computeGradientVariationWithRadius(obj)
            s.mesh                 = obj.gExperiment.backgroundMesh;
            s.regularizedPerimeter = obj.gExperiment.regularizedPerimeter;
            s.inputFile            = obj.inputFiles{obj.iMesh};
            s.nameCase             = obj.nameCase;
            s.outputFolder         = obj.outputFolder;
            s.domainLength         = obj.gExperiment.domainLength();
            gComputer = GradientVariationWithRadiusComputer(s);
            gComputer.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function mkdirOptional(dirName)
% Create directory if it doesn't exist
if(~exist(dirName,'dir'))
    mkdir(dirName);
end
end

repo: CategoryShapes


text:
function varargout = spm_render_vol(varargin)
% surface render a memory mapped 8 bit image - a compiled routine
% FORMAT [REN, ZBUF, X, Y, Z] = spm_render_vol(V, A, [i j], [u n])
% V       -  is the memory mapped volume
% A       -  {4 x 4} affine transformation matrix
% [i j]   -  dimensions of REN
% [u n]   -  u is threhsold at which voxels are 'solid'
%            n is the number of nearest neighbours to use to determine the
%            surface orientation
% REN     -  is the rendered image
% ZBUF    -  distance from the view plane to the object's surface
% X, Y, Z -  are images containing the coordinates of the voxels on the
%            surface of the volume.
%_______________________________________________________________________
%
% [i j] defines the two dimensions of the output image. The coordinates
% in 3-D space of the voxels in this image are assumed to range from
% 1,1,0 to i,j,0.
%
% For each pixel in the volume, the coordinates (x,y,z & 1) are
% multiplied by the matrix A, to give the image coordinates that these
% voxels map to.
%
% The threshold at which voxels are assumed to be solid pertains to the
% 8-bit data i.e. {0 - 255}
%
% Illumination is assumed to be from the viewplane
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_render_vol.m 1143 2008-02-07 19:33:33Z spm $

%-This is merely the help file for the compiled routine
error('spm_render_vol.c not compiled - see Makefile')

author: spm


text:
function value = r8poly_value_horner ( m, c, x )

%*****************************************************************************80
%
%% R8POLY_VALUE_HORNER evaluates a polynomial using Horner's method.
%
%  Discussion:
%
%    The polynomial
%
%      p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m
%
%    is to be evaluated at the value X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the degree.
%
%    Input, real C(M+1), the polynomial coefficients.
%    C(I) is the coefficient of X^(I-1).
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the polynomial value.
%
  value = c(m+1);
  for i = m : -1 : 1
    value = value * x + c(i);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio]...
    = olmar2_run(fid, data, epsilon, alpha, tc, opts)
% This program simulates the OLMAR-2 algorithm
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%    = olmar2_run(fid, data, epsilon, alpha, tc, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% epsilon: mean reversion threshold
% alpha: trade off parameter for calculating moving average [0, 1]
% tc: transaction cost rate parameter
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%           = olmar2_run(fid, data, epsilon, alpha, tc, opts)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n, m] = size(data);

% Return variables
cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);

% Portfolio weights, starting with uniform portfolio
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);  % Last closing price adjusted portfolio
daily_portfolio = zeros(n, m);

% print file head
fprintf(fid, '-------------------------------------\n');
fprintf(fid, 'Parameters [epsilon:%.2f, alpha:%.2f, tc:%.4f]\n', ...
    epsilon, alpha, tc);
fprintf(fid, 'day\t Daily Return\t Total return\n');

fprintf(1, '-------------------------------------\n');
if(~opts.quiet_mode)
    fprintf(1, 'Parameters [epsilon:%.2f, alpha:%.2f, tc:%.4f]\n', ...
        epsilon, alpha, tc);
    fprintf(1, 'day\t Daily Return\t Total return\n');
end

data_phi = ones(1, m);

%% Trading
if (opts.progress)
	progress = waitbar(0,'Executing Algorithm...');
end
for t = 1:1:n,
    % Step 1: Receive stock price relatives
    if (t >= 2)
        [day_weight, data_phi] ...
            = olmar2_kernel(data(1:t-1, :), data_phi, day_weight, epsilon, alpha);
    end

    % Normalize the constraint, always useless
    day_weight = day_weight./sum(day_weight);
    daily_portfolio(t, :) = day_weight';

    if or((day_weight < -0.00001+zeros(size(day_weight))), (day_weight'*ones(m, 1)>1.00001))
        fprintf(1, 'mrpa_expert: t=%d, sum(day_weight)=%d, pause', t, day_weight'*ones(m, 1));
        pause;
    end

    % Step 2: Cal t's daily return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    % fprintf(1, '%d\t%.2f\t%.2f\t%.2f\n', t, day_weight(1), day_weight(2), daily_ret(t, 1));
    % Adjust weight(t, :) for the transaction cost issue
    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

    % Debug information
    % Time consuming part, other way?
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
    if (~opts.quiet_mode),
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
        end
    end
    if (opts.progress)
        if mod(t, 50) == 0
            waitbar((t/n));
        end
    end
end

% Debug Information
fprintf(fid, 'OLMAR-2(epsilon:%.2f, alpha:%.2f, tc:%.4f), Final return: %.2f\n', ...
    epsilon, alpha, tc, cum_ret);
fprintf(fid, '-------------------------------------\n');
fprintf(1, 'OLMAR-2(epsilon:%.2f, alpha:%.2f, tc:%.4f), Final return: %.2f\n', ...
    epsilon, alpha, tc, cum_ret);
fprintf(1, '-------------------------------------\n');
    if (opts.progress)
        close(progress);
    end
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function FillBetween(Xdata,YLower,YUpper,AreaColor,Alpha)

if size(Xdata) == size(YLower)
    fill([vertcat(Xdata,Xdata(end:-1:1))],[vertcat(YLower,YUpper(end:-1:1))],AreaColor,'FaceAlpha',Alpha);
else
    warning('MATLAB:InvlaidArgumentSize','The size of the input is not consistent');
end

author: Sable

repo: mcbench-benchmarks


text:
function chns = chnsCompute( I, varargin )
% Compute channel features at a single scale given an input image.
%
% Compute the channel features as described in:
%  P. Doll�r, Z. Tu, P. Perona and S. Belongie
%  "Integral Channel Features", BMVC 2009.
% Channel features have proven very effective in sliding window object
% detection, both in terms of *accuracy* and *speed*. Numerous feature
% types including histogram of gradients (hog) can be converted into
% channel features, and overall, channels are general and powerful.
%
% Given an input image I, a corresponding channel is a registered map of I,
% where the output pixels are computed from corresponding patches of input
% pixels (thus preserving overall image layout). A trivial channel is
% simply the input grayscale image, likewise for a color image each color
% channel can serve as a channel. Other channels can be computed using
% linear or non-linear transformations of I, various choices implemented
% here are described below. The only constraint is that channels must be
% translationally invariant (i.e. translating the input image or the
% resulting channels gives the same result). This allows for fast object
% detection, as the channels can be computed once on the entire image
% rather than separately for each overlapping detection window.
%
% Currently, three channel types are available by default (to date, these
% have proven the most effective for sliding window object detection):
%  (1) color channels (computed using rgbConvert.m)
%  (2) gradient magnitude (computed using gradientMag.m)
%  (3) quantized gradient channels (computed using gradientHist.m)
% For more information about each channel type, including the exact input
% parameters and their meanings, see the respective m-files which perform
% the actual computatons (chnsCompute is essentially a wrapper function).
% The converted color channels serve as input to gradientMag/gradientHist.
%
% Additionally, custom channels can be specified via an optional struct
% array "pCustom" which may have 0 or more custom channel definitions. Each
% custom channel is generated via a call to "chns=feval(hFunc,I,pFunc{:})".
% The color space of I is determined by pColor.colorSpace, use the setting
% colorSpace='orig' if the input image is not an 'rgb' image and should be
% left unchanged (e.g. if I has multiple channels). The input I will have
% type single and the output of hFunc should also have type single.
%
% "shrink" (which should be an integer) determines the amount to subsample
% the computed channels (in applications such as detection subsamping does
% not affect performance). The params for each channel type are described
% in detail in the respective function. In addition, each channel type has
% a param "enabled" that determines if the channel is computed. If
% chnsCompute() is called with no inputs, the output is the complete
% default params (pChns). Otherwise the outputs are the computed channels
% and additional meta-data (see below). The channels are computed at a
% single scale, for (fast) multi-scale channel computation see chnsPyramid.
%
% An emphasis has been placed on speed, with the code undergoing heavy
% optimization. Computing the full set of channels used in the BMVC09 paper
% referenced above on a 480x640 image runs over *100 fps* on a single core
% of a machine from 2011 (although runtime depends on input parameters).
%
% USAGE
%  pChns = chnsCompute()
%  chns = chnsCompute( I, pChns )
%
% INPUTS
%  I           - [hxwx3] input image (uint8 or single/double in [0,1])
%  pChns       - parameters (struct or name/value pairs)
%   .shrink       - [4] integer downsampling amount for channels
%   .pColor       - parameters for color space:
%     .enabled      - [1] if true enable color channels
%     .smooth       - [1] radius for image smoothing (using convTri)
%     .colorSpace   - ['luv'] choices are: 'gray', 'rgb', 'hsv', 'orig'
%   .pGradMag     - parameters for gradient magnitude:
%     .enabled      - [1] if true enable gradient magnitude channel
%     .colorChn     - [0] if>0 color channel to use for grad computation
%     .normRad      - [5] normalization radius for gradient
%     .normConst    - [.005] normalization constant for gradient
%     .full         - [0] if true compute angles in [0,2*pi) else in [0,pi)
%   .pGradHist    - parameters for gradient histograms:
%     .enabled      - [1] if true enable gradient histogram channels
%     .binSize      - [shrink] spatial bin size (defaults to shrink)
%     .nOrients     - [6] number of orientation channels
%     .softBin      - [0] if true use "soft" bilinear spatial binning
%     .useHog       - [0] if true perform 4-way hog normalization/clipping
%     .clipHog      - [.2] value at which to clip hog histogram bins
%   .pCustom      - parameters for custom channels (optional struct array):
%     .enabled      - [1] if true enable custom channel type
%     .name         - ['REQ'] custom channel type name
%     .hFunc        - ['REQ'] function handle for computing custom channels
%     .pFunc        - [{}] additional params for chns=hFunc(I,pFunc{:})
%     .padWith      - [0] how channel should be padded (e.g. 0,'replicate')
%   .complete     - [] if true does not check/set default vals in pChns
%
% OUTPUTS
%  chns       - output struct
%   .pChns      - exact input parameters used
%   .nTypes     - number of channel types
%   .data       - [nTypes x 1] cell [h/shrink x w/shrink x nChns] channels
%   .info       - [nTypes x 1] struct array
%     .name       - channel type name
%     .pChn       - exact input parameters for given channel type
%     .nChns      - number of channels for given channel type
%     .padWith    - how channel should be padded (0,'replicate')
%
% EXAMPLE - default channels
%  I=imResample(imread('peppers.png'),[480 640]); pChns=chnsCompute();
%  tic, for i=1:100, chns=chnsCompute(I,pChns); end; toc
%  figure(1); montage2(cat(3,chns.data{:}));
%
% EXAMPLE - default + custom channels
%  I=imResample(imread('peppers.png'),[480 640]); pChns=chnsCompute();
%  hFunc=@(I) 5*sqrt(max(0,max(convBox(I.^2,2)-convBox(I,2).^2,[],3)));
%  pChns.pCustom=struct('name','Std02','hFunc',hFunc); pChns.complete=0;
%  tic, chns=chnsCompute(I,pChns); toc
%  figure(1); im(chns.data{4});
%
% See also rgbConvert, gradientMag, gradientHist, chnsPyramid
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.23
% Copyright 2014 Piotr Dollar & Ron Appel.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get default parameters pChns
if(nargin==2), pChns=varargin{1}; else pChns=[]; end
if( ~isfield(pChns,'complete') || pChns.complete~=1 || isempty(I) )
  p=struct('enabled',{},'name',{},'hFunc',{},'pFunc',{},'padWith',{});
  pChns = getPrmDflt(varargin,{'shrink',4,'pColor',{},'pGradMag',{},...
    'pGradHist',{},'pCustom',p,'complete',1},1);
  pChns.pColor = getPrmDflt( pChns.pColor, {'enabled',1,...
    'smooth',1, 'colorSpace','luv'}, 1 );
  pChns.pGradMag = getPrmDflt( pChns.pGradMag, {'enabled',1,...
    'colorChn',0,'normRad',5,'normConst',.005,'full',0}, 1 );
  pChns.pGradHist = getPrmDflt( pChns.pGradHist, {'enabled',1,...
    'binSize',[],'nOrients',6,'softBin',0,'useHog',0,'clipHog',.2}, 1 );
  nc=length(pChns.pCustom); pc=cell(1,nc);
  for i=1:nc, pc{i} = getPrmDflt( pChns.pCustom(i), {'enabled',1,...
      'name','REQ','hFunc','REQ','pFunc',{},'padWith',0}, 1 ); end
  if( nc>0 ), pChns.pCustom=[pc{:}]; end
end
if(nargin==0), chns=pChns; return; end

% create output struct
info=struct('name',{},'pChn',{},'nChns',{},'padWith',{});
chns=struct('pChns',pChns,'nTypes',0,'data',{{}},'info',info);

% crop I so divisible by shrink and get target dimensions
shrink=pChns.shrink; [h,w,~]=size(I); cr=mod([h w],shrink);
if(any(cr)), h=h-cr(1); w=w-cr(2); I=I(1:h,1:w,:); end
h=h/shrink; w=w/shrink;

% compute color channels
p=pChns.pColor; nm='color channels';
I=rgbConvert(I,p.colorSpace); I=convTri(I,p.smooth);
if(p.enabled), chns=addChn(chns,I,nm,p,'replicate',h,w); end

% compute gradient magnitude channel
p=pChns.pGradMag; nm='gradient magnitude';
full=0; if(isfield(p,'full')), full=p.full; end
if( pChns.pGradHist.enabled )
  [M,O]=gradientMag(I,p.colorChn,p.normRad,p.normConst,full);
elseif( p.enabled )
  M=gradientMag(I,p.colorChn,p.normRad,p.normConst,full);
end
if(p.enabled), chns=addChn(chns,M,nm,p,0,h,w); end

% compute gradient histgoram channels
p=pChns.pGradHist; nm='gradient histogram';
if( p.enabled )
  binSize=p.binSize; if(isempty(binSize)), binSize=shrink; end
  H=gradientHist(M,O,binSize,p.nOrients,p.softBin,p.useHog,p.clipHog,full);
  chns=addChn(chns,H,nm,pChns.pGradHist,0,h,w);
end

% compute custom channels
p=pChns.pCustom;
for i=find( [p.enabled] )
  C=feval(p(i).hFunc,I,p(i).pFunc{:});
  chns=addChn(chns,C,p(i).name,p(i),p(i).padWith,h,w);
end

end

function chns = addChn( chns, data, name, pChn, padWith, h, w )
% Helper function to add a channel to chns.
[h1,w1,~]=size(data);
if(h1~=h || w1~=w), data=imResampleMex(data,h,w,1);
  assert(all(mod([h1 w1]./[h w],1)==0)); end
chns.data{end+1}=data; chns.nTypes=chns.nTypes+1;
chns.info(end+1)=struct('name',name,'pChn',pChn,...
  'nChns',size(data,3),'padWith',padWith);
end

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function [y_rvm] = rvr_predict(X,  model)
% RVR_PREDICT % Predicts labels and computes accuracy for new
% data with a learnt RVM model
%   input ----------------------------------------------------------------
%
%       o X        : (N x D), N  input data points of D dimensionality.
%
%
%
%   output ----------------------------------------------------------------
%
%       o y_rvm            : (N x 1), decision values
%
%
%% Predict labels given model and data

%%%%%% Test
weights    = model.weights;
kernel_    = model.kernel_;
width      = model.width;
bias       = model.bias;
RVs_idx    = model.RVs_idx;
RVs        = model.RVs;

% Compute RVM over test data and calculate error
PHI	= SB1_KernelFunction(X(:), RVs, kernel_, width);
y_rvm	= PHI*weights + bias;

end

author: epfl-lasa


text:
function [ x, w ] = rule_adjust ( a, b, c, d, n, x, w )

%*****************************************************************************80
%
%% RULE_ADJUST maps a quadrature rule from [A,B] to [C,D].
%
%  Discussion:
%
%    Most quadrature rules are defined on a special interval, like
%    [-1,1] or [0,1].  To integrate over an interval, the abscissas
%    and weights must be adjusted.  This can be done on the fly,
%    or by calling this routine.
%
%    If the weight function W(X) is not 1, then the W vector will
%    require further adjustment by the user.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the endpoints of the definition interval.
%
%    Input, real C, D, the endpoints of the integration interval.
%
%    Input, integer N, the number of abscissas and weights.
%
%    Input, real X(N), the abscissas.
%
%    Input, real W(N), the weights.
%
%    Output, real X(N), the adjusted abscissas.
%
%    Output, real W(N), the adjusted weights.
%
  x(1:n) = ( ( b - x(1:n)     ) * c   ...
           + (     x(1:n) - a ) * d ) ...
           / ( b          - a );

  w(1:n) = ( ( d - c ) / ( b - a ) ) * w(1:n);

  return
end


text:
%sim_coherenceOrder.m
%Dana Goerzen and Jamie Near, 2021.
%
% USAGE:
% out = sim_coherenceOrder(spinSys);
%
% DESCRIPTION:
% Creates an n x n matrix that represents the coherence order of each element
% of the density matrix for the given spin system For use in nulling signal
% incoherences during pulse sequence simulations
%
% INPUTS:
% spinSys     = spin system definition structure.
%
% OUTPUTS:
% out       = n x n coherence order matrix for spin system.

function out=sim_coherenceOrder(spinSys)
out=cell(length(spinSys),1);
for m=1:length(spinSys)
    Nspins=length(spinSys(m).J);

    p0=[0 -1;1 0];
    p=p0;
    for n=1:Nspins-1
        p=kron(ones(2),p)+kron(p0,ones(size(p)));
    end
    out{m,1}=p;
end
out=p;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function str = cvx_replicate_structure( str, sz )
%CVX_REPLICATE_STRUCTURE Apply matrix structure to N-D arrays.

szs  = size( str );
nmat = prod( sz );
omat = ones( 1, nmat );
[ r, c, v ] = find( str );
nelm = length( r );
oelm = ones( 1, nelm );
r = r( : );
nvec = ( 0 : nmat - 1 ) * szs( 1 );
r = r( :, omat ) + nvec( oelm, : );
c = c( : );
c = c( :, omat );
nvec = ( 0 : nmat - 1 ) * szs( 2 );
c = c( :, omat ) + nvec( oelm, : );
v = v( : );
v = v( :, omat );
str = sparse( r, c, v, nmat * szs( 1 ), nmat * szs( 2 ) );

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

repo: REST


text:
function g = kernDiagGradient(kern, x, covDiag)

% KERNDIAGGRADIENT Compute the gradient of the kernel's parameters for the diagonal.
% FORMAT
% DESC computes the gradient of functions of the diagonal of the
% kernel matrix with respect to the parameters of the kernel. The
% parameters' gradients are returned in the order given by the
% kernExtractParam command.
% ARG kern : the kernel structure for which the gradients are computed.
% ARG x : the input data for which the gradient is being computed.
% ARG factors : partial derivatives of the function of interest
% with respect to the diagonal elements of the kernel matrix.
% RETURN g : gradients of the relevant function with respect to each of the parameters. Ordering should match the ordering given in kernExtractParam.
%
% SEEALSO : kernDiagGradient, kernExtractParam, kernGradient

% KERN

fileName = [kern.type 'KernDiagGradient'];
if exist(fileName) == 2
  fhandle = str2func(fileName);
  g = fhandle(kern, x, covDiag);
else
  fhandle = str2func([kern.type 'KernGradient']);
  g = zeros(1, kern.nParams);
  for i = 1:size(x, 1)
    g = g ...
        + fhandle(kern, x(i, :), covDiag(i));
  end
end
% Check if parameters are being optimised in a transformed space.
factors = kernFactors(kern, 'gradfact');
g(factors.index) = g(factors.index).*factors.val;


text:
function generate_images (outputpath,amount,size)
selpath = outputpath;
noise = 0.0005; %noise 0 bis 0.05
z_move=10; %z_move
partAm=size^2*0.01;
Z=0.333; %sheet thickness
dt=3; %particle diameter
ddt=1; %particle diameter variation

for iiii = 0:amount-1
disp(['Progress: ' int2str(iiii/amount*100) ' %']);
	%% Generate random artificial particle images
	u = rand(1)*ones(size,size);
	v = rand(1)*ones(size,size);
	offsety=v; %zero displacement in y direction
	offsetx=u; %uniform x displacement increases from 0 to 99/10 pixels
	[x,y]=meshgrid(1:1:size);
	i=[];
	j=[];
	sizey=size;
	sizex=size;
	A=zeros(sizey,sizex);
	B=A;
	z0_pre=randn(partAm,1); %normal distributed sheet intensity
	randn('state', sum(100*clock)); %#ok<*RAND>
	z1_pre=randn(partAm,1); %normal distributed sheet intensity
	z0=z0_pre*(z_move/200+0.5)+z1_pre*(1-((z_move/200+0.5)));
	z1=z1_pre*(z_move/200+0.5)+z0_pre*(1-((z_move/200+0.5)));

	I0=255*exp(-(Z^2./(0.125*z0.^2))); %particle intensity
	I0(I0>255)=255;
	I0(I0<0)=0;

	I1=255*exp(-(Z^2./(0.125*z1.^2))); %particle intensity
	I1(I1>255)=255;
	I1(I1<0)=0;

	randn('state', sum(100*clock));
	d=randn(partAm,1)/2; %particle diameter distribution
	d=dt+d*ddt;
	d(d<0)=0;
	rand('state', sum(100*clock));
	x0=rand(partAm,1)*sizex;
	y0=rand(partAm,1)*sizey;
	rd = -8.0 ./ d.^2;

	xlimit1=floor(x0-d/2); %x min particle extent image1
	xlimit2=ceil(x0+d/2); %x max particle extent image1
	ylimit1=floor(y0-d/2); %y min particle extent image1
	ylimit2=ceil(y0+d/2); %y max particle extent image1
	xlimit2(xlimit2>sizex)=sizex;
	xlimit1(xlimit1<1)=1;
	ylimit2(ylimit2>sizey)=sizey;
	ylimit1(ylimit1<1)=1;

	%calculate particle extents for image2 (shifted image)
	x0integer=round(x0);
	x0integer(x0integer>sizex)=sizex;
	x0integer(x0integer<1)=1;
	y0integer=round(y0);
	y0integer(y0integer>sizey)=sizey;
	y0integer(y0integer<1)=1;

	xlimit3=zeros(partAm,1);
	xlimit4=xlimit3;
	ylimit3=xlimit3;
	ylimit4=xlimit3;
	for n=1:partAm
		xlimit3(n,1)=floor(x0(n)-d(n)/2-offsetx((y0integer(n)),(x0integer(n)))); %x min particle extent image2
		xlimit4(n,1)=ceil(x0(n)+d(n)/2-offsetx((y0integer(n)),(x0integer(n)))); %x max particle extent image2
		ylimit3(n,1)=floor(y0(n)-d(n)/2-offsety((y0integer(n)),(x0integer(n)))); %y min particle extent image2
		ylimit4(n,1)=ceil(y0(n)+d(n)/2-offsety((y0integer(n)),(x0integer(n)))); %y max particle extent image2
	end
	xlimit3(xlimit3<1)=1;
	xlimit4(xlimit4>sizex)=sizex;
	ylimit3(ylimit3<1)=1;
	ylimit4(ylimit4>sizey)=sizey;

	ctr=0;
	for n=1:partAm
		ctr=ctr+1;
		if ctr==10000
			ctr=0;
			%fprintf('.')
		end
		r = rd(n);
		for j=xlimit1(n):xlimit2(n)
			rj = (j-x0(n))^2;
			for i=ylimit1(n):ylimit2(n)
				A(i,j)=A(i,j)+I0(n)*exp((rj+(i-y0(n))^2)*r);
			end
		end
		for j=xlimit3(n):xlimit4(n)
			for i=ylimit3(n):ylimit4(n)
				B(i,j)=B(i,j)+I1(n)*exp((-(j-x0(n)+offsetx(i,j))^2-(i-y0(n)+offsety(i,j))^2)*-rd(n)); %place particle with gaussian intensity profile
			end
		end
	end

	%% Create random Background "glow"
	bg=im2bw(rand(size,size),0.99999);
	SE = strel('line',round(rand*90)+90,round(rand*180));
	bg2=imdilate(bg,SE);
	bg3 = imgaussfilt(double(bg2),40);
	bg3=bg3/max(max(bg3));
	bg3=bg3*0.1;
	bg=im2bw(rand(size,size),0.999997);
	SE = strel('disk',round(rand*100)+200,4);
	bg2=imdilate(bg,SE);
	bg4 = 0.12*imgaussfilt(double(bg2),100);
	bg3(isnan(bg3))=0;
	bg4(isnan(bg4))=0;

	A=A+bg3*rand*255+bg4*rand*255;
	B=B+bg3*rand*255+bg4*rand*255;

	A(A>255)=255;
	B(B>255)=255;
	A=imnoise(uint8(A),'gaussian',0,noise);
	B=imnoise(uint8(B),'gaussian',0,noise);

	A=uint8(A);
	B=uint8(B);

	%True velocities, save if desired
	x_real=x;
	y_real=y;
	u_real=u;
	v_real=v;

	imwrite(A,[selpath '\synth_image_' sprintf('%5.5d',iiii) '_A.jpg'],'Quality',97);
	imwrite(B,[selpath '\synth_image_' sprintf('%5.5d',iiii) '_B.jpg'],'Quality',97);

end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
%% Exercise with cosmo_correlation_measure
% Split-half correlation analysis using cosmo_correlation_measure
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Single subject analysis over all categories in VT

% Set the data path
config=cosmo_config();
study_path=fullfile(config.tutorial_data_path,'ak6');

subject_id='s01';
mask_label='vt_mask';
data_path=fullfile(study_path,subject_id); % data from subject s01

% Define data locations and load data from even and odd runs
mask_fn=fullfile(data_path, [mask_label '.nii']); % whole brain

% Load data for odd runs
data_odd_fn=fullfile(data_path,'glm_T_stats_odd.nii');
ds_odd=cosmo_fmri_dataset(data_odd_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',1);

% Load data for even runs
data_even_fn=fullfile(data_path,'glm_T_stats_even.nii');
ds_even=cosmo_fmri_dataset(data_even_fn,'mask',mask_fn,...
                            'targets',1:6,'chunks',2);

% Combine even and odd runs into single dataset, store in ds_odd_even
% Hint: use cosmo_stack
% >@@>
ds_odd_even=cosmo_stack({ds_odd, ds_even});
% <@@<

% remove constant features (due to liberal masking)
ds_odd_even=cosmo_remove_useless_data(ds_odd_even);

% display dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_odd_even);

% Part 1: compute correlation difference
% Hint: - apply cosmo_correlation_measure to ds_odd_even
%       - because correlation difference and split-half correlations
%         are the defaults, no second argument (options) is required
%         for cosmo_correlation_measure.
% >@@>
ds_corr=cosmo_correlation_measure(ds_odd_even);
% <@@<
fprintf('\nDataset output (correlation difference):\n');
cosmo_disp(ds_corr);

fprintf(['Average correlation difference between matching and '...
            'non-matching categories in %s for %s is %.3f\n'],...
            mask_label, subject_id, ds_corr.samples);

% Part 2: compute the raw Fisher-transformed correlation matrix,
% and store the results in 'c_raw'
%
% (Hint: use a struct 'args' with args.output='correlation' as second
% argument for cosmo_correlation_measure)
% >@@>
args=struct();
args.output='correlation';
c_raw=cosmo_correlation_measure(ds_odd_even,args);
% <@@<
fprintf('\nDataset output (Fisher-transformed correlations):\n');
cosmo_disp(c_raw)

% Because a measure returns .samples as a column vector, the
% confusion matrix is returned in a flattened form.
% The data can be put back in matrix form using cosmo_unflatten.
matrices=cosmo_unflatten(c_raw,1);
imagesc(matrices);
colorbar();

%% Group analysis over all categories in VT using correlation measure

% Compute correlation measure for each subject individually.
subject_ids={'s01','s02','s03','s04','s05','s06','s07','s08'};
nsubjects=numel(subject_ids);

mask_label='vt_mask';

% allocate a cell for the output of the measure for each subject
ds_corr_cell=cell(nsubjects,1);

% Apply correlation measure for each subject using
% cosmo_correlation_measure.
% - compute the correlation measure using cosmo_correlation_measure;
%   this should give a dataset structure with ds.samples of size 1x1
% - because statistics are computed later, add to the output:
%   * .sa.chunks : the subject number
%   * .sa.targets: 1 (for each subject)
% - store the result of the k-th subject in ds_corr_cell{k}

for subject_num=1:nsubjects
    subject_id=subject_ids{subject_num};

    % set path for this subject
    data_path=fullfile(study_path,subject_id);

    % Define data locations and load data from even and odd runs
    mask_fn=fullfile(data_path, [mask_label '.nii']); % whole brain

    data_odd_fn=fullfile(data_path,'glm_T_stats_odd.nii');
    ds_odd=cosmo_fmri_dataset(data_odd_fn,'mask',mask_fn,...
                                'targets',1:6,'chunks',1);

    data_even_fn=fullfile(data_path,'glm_T_stats_even.nii');
    ds_even=cosmo_fmri_dataset(data_even_fn,'mask',mask_fn,...
                                'targets',1:6,'chunks',2);

    % Combine even and odd runs
    ds_odd_even=cosmo_stack({ds_odd, ds_even});

    % remove constant features (due to liberal masking)
    ds_odd_even=cosmo_remove_useless_data(ds_odd_even);

    % apply correlation measure for this subject, and store
    % the result in 'ds_corr'
    % >@@>
    ds_corr=cosmo_correlation_measure(ds_odd_even);
    % <@@<
    % set targets and chunks for the output, so that cosmo_stat can be used
    % below
    ds_corr.sa.targets=1;
    ds_corr.sa.chunks=subject_num;

    % store result in the cell array 'ds_corr_cell' at the 'subject_num'-th
    % position
    % >@@>
    ds_corr_cell{subject_num}=ds_corr;
    % <@@<
end

% combine the data from all subjects (in 'ds_corr_cell') into one dataset
% using cosmo_stack, and assign the result to a variable 'ds_all'
% >@@>
ds_all=cosmo_stack(ds_corr_cell);
% <@@<

%% Compute group-level statistics

% run one-sample t-test again zero to compute the t-value and p-values
% against the null hypothesis of a mean of zero.
%
% hint: use matlab stats toolbox's ttest (if present), or
%       use cosmo_stat:
%       * for t-statistics, use 't' as second argument
%       * to convert a 't' to 'p' value, use 'p' as third argument
% >@@>
ds_t=cosmo_stat(ds_all,'t');     % t-test against zero
ds_p=cosmo_stat(ds_all,'t','p'); % convert to p-value

samples=ds_all.samples;

fprintf(['correlation difference in %s at group level: '...
           '%.3f +/- %.3f, %s=%.3f, p=%.5f (using cosmo_stat)\n'],...
            mask_label,mean(samples),std(samples),...
            ds_t.sa.stats{1},ds_t.samples,ds_p.samples);
% <@@<

% Using matlab's stat toolbox (if present)
if cosmo_check_external('@stats',false)
    % >@@>
    [h,p,ci,stats]=ttest(samples);
    fprintf(['correlation difference in %s at group level: '...
            '%.3f +/- %.3f, t_%d=%.3f, p=%.5f (using matlab stats '...
            'toolbox)\n'],...
            mask_label,mean(samples),std(samples),stats.df,stats.tstat,p);
    % <@@<
else
    fprintf('Matlab stats toolbox not found\n');
end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function imout=readspe(filename)
% imout=readspe(filename)
% Reads .spe file from the CSSTORM evaluation
fid = fopen(filename,'r','l');
header = fread(fid,2050,'*uint16');
imgs = fread(fid, inf, 'float32');
fclose(fid);
z_dim = double(header(724));
x_dim = double(header(22));
im=reshape(imgs,x_dim,x_dim,z_dim);
imnan=isnan(im);
imn=zeros(size(im));
imn(~imnan)=im(~imnan);
imout=rot90stack(fliplrstack(imn));

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function F = D2GaussFunctionRot(x,xdata)
%% x = [Amp, x0, wx, y0, wy, fi]
%[X,Y] = meshgrid(x,y)
%  xdata(:,:,1) = X
%  xdata(:,:,2) = Y
% Mrot = [cos(fi) -sin(fi); sin(fi) cos(fi)]
%%
xdatarot(:,:,1)= xdata(:,:,1)*cos(x(6)) - xdata(:,:,2)*sin(x(6));
xdatarot(:,:,2)= xdata(:,:,1)*sin(x(6)) + xdata(:,:,2)*cos(x(6));
x0rot = x(2)*cos(x(6)) - x(4)*sin(x(6));
y0rot = x(2)*sin(x(6)) + x(4)*cos(x(6));

F = x(1)*exp(   -((xdatarot(:,:,1)-x0rot).^2/(2*x(3)^2) + (xdatarot(:,:,2)-y0rot).^2/(2*x(5)^2) )    );

% figure(3)
% alpha(0)
% imagesc(F)
% colormap('gray')
% figure(gcf)%bring current figure to front
% drawnow
% beep
% pause %Wait for keystroke


text:
function linplus_test577 ( )

%*****************************************************************************80
%
%% TEST577 tests R8SS_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 9;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST577\n' );
  fprintf ( 1, '  For a symmetric skyline storage matrix,\n' );
  fprintf ( 1, '  R8SS_INDICATOR computes an indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  [ na, diag, a ] = r8ss_indicator ( n );

  r8ss_print ( n, na, diag, a, '  The R8SS indicator matrix:' );

  return
end


text:
function value = year_length_months_persian ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_PERSIAN returns the number of months in a Persian year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  value = 12;

  return
end


text:
classdef Postprocess < handle

    properties (Access = protected)
        outFileName
        resultsDir

        resultPrinter
        meshPrinter

        mshDataBase
        resDataBase
    end

     methods (Access = public)

        function obj = Postprocess(postCase,d)
            obj.init(d);
            obj.createOutputDirectory();
            obj.createMeshPrinter();
            obj.createResultPrinter(postCase);
        end

        function print(obj,iter,d)
            obj.printMeshFile(iter);
            obj.printResFile(iter,d)
        end

        function r = getResFile(obj)
            r = obj.outFileName;
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.outFileName = d.outFileName;
            obj.createResultsDirName();
            obj.computeDataBaseForMeshFile(d);
            obj.computeDataBaseForResFile(d);
        end

        function createMeshPrinter(obj)
            obj.meshPrinter = MeshPrinter();
        end

        function createResultPrinter(obj,postCase)
            obj.resultPrinter  = ResultsPrinter.create(postCase,obj.resDataBase);
        end

        function printResFile(obj,iter,d)
            obj.resultPrinter.print(iter,d);
        end

       	function printMeshFile(obj,iter)
            d = obj.mshDataBase;
            d.iter = iter;
            obj.meshPrinter.print(d);
        end

        function createOutputDirectory(obj)
            dir = obj.resultsDir;
            if ~exist(dir,'dir')
                mkdir(dir)
            end
        end

        function computeDataBaseForResFile(obj,dI)
            obj.resDataBase = dI;
            obj.resDataBase.resultsDir = obj.resultsDir;
        end

        function computeDataBaseForMeshFile(obj,dI)
            d.coordinates    = dI.coordinates;
            d.connectivities = dI.connectivities;
            d.outFileName    = dI.outFileName;
            d.npnod          = dI.npnod;
            d.pdim           = dI.pdim;
            d.nnode          = dI.nnode;
            d.nelem          = dI.nelem;
            d.ndim           = dI.ndim;
            d.etype          = dI.etype;
            d.resultsDir     = obj.resultsDir;
            obj.mshDataBase = d;
        end

        function createResultsDirName(obj)
           path = pwd;
           obj.resultsDir = fullfile(path,'Output',obj.outFileName);
        end

     end

end


text:
% rundesignsimMAP2
% Tor Wager, 8/29/01
% if your input var name is already a variable in the wkspace, uses that...
% otherwise creates one with zeros.

warning off
index = 1;
name = input('Enter the name of the output variable and files ','s');
% eval([name ' = zeros(10,10);'])

% -------------------------------------------------
% common user variables
% -------------------------------------------------
tickRes = [.1 .1];
ISIrange = [.1 16]
FCrange = [.01 .81];
GA.sizegenerations = 100;

myISI = ISIrange(1):tickRes(1):ISIrange(2);
myFC = FCrange(1):tickRes(2):FCrange(2);

ISIsize = length(myISI);
FCsize = length(myFC);

ISItick = 1:max(round(ISIsize./25),1):ISIsize;
ISIticklab = myISI(ISItick);
FCtick = 1:2:FCsize;
FCticklab = myFC(FCtick);

eval([' if ~(exist(''' name ''') == 1), disp(''Creating new variable''),' name ' = zeros(ISIsize,FCsize);, else,disp(''Using existing variable''),end'])

disp(['Starting ' num2str(prod([ISIsize FCsize])) ' surface points.'])

% ----------------------------------------------------------------
% * HRF GAMMA FUNCTION - spm 99
% ----------------------------------------------------------------

   HRF = spm_hrf(.1);						% SPM HRF sampled at .1 s
   HRF = HRF/ max(HRF);

% ----------------------------------------------------------------
% * make the map
% ----------------------------------------------------------------

for RR = 1:size(myISI,2)
	% set ISI
    GA.ISI = myISI(RR);

    % ----------------------------------------------------------------
    % * get smoothing matrix and autocorrelation matrix
    % ----------------------------------------------------------------
    numStim = ceil(GA.scanLength / (GA.ISI));
    numsamps = ceil(numStim*GA.ISI/GA.TR);
	disp('  ...getting S and svi')
	[S,Vi,svi] = getSmoothing(GA.HPlength,GA.LPsmooth,GA.TR,numsamps,GA.xc);
	if isempty(S), S = 1;, end

   	clear Vi

	for CC = 1:size(myFC,2)
		disp(['Starting model ' num2str(index)])

		% set rest frequency
		GA.freqConditions(end) = myFC(CC);
		GA.freqConditions(1:end-1) = (1 - myFC(CC)) / (size(GA.freqConditions,2)-1);

		disp(['FREQC and ISI are:	' num2str(GA.freqConditions) ' and ' num2str(GA.ISI)])

        eval(['mymat = ' name ';'])
        if mymat(RR,CC) == 0        % only do it if it's not done yet.

        % ----------------------------------------------------------------
        % * randomize initial set of organisms
        % ----------------------------------------------------------------
        unsortedList = [];
        numStimEachCond = ceil((numStim) * GA.freqConditions);
	    clear listMatrix

	    for i = 1:size(GA.freqConditions,2)
	    	% unsorted list of all stims in proper frequencies
   		    startIndex = size(unsortedList,1)+1;
		    % multiplied by 1.5 because lists were coming out too short.
   		    unsortedList(startIndex:(startIndex + 2*numStimEachCond(i)-1),1) = GA.conditions(i);
        end

        disp('  ...randomizing organism start state')
        for z = 1:GA.sizeGenerations % generate random ordering of x conditions
  	        stimList = getRandom(unsortedList);
            listMatrix(:,z) = stimList; % a row for each subject
        end

            [fitness,convar,xtxi,models] = testlist(listMatrix,GA,HRF,'svi',svi,'S',S);
            convar = 1./convar;
	        eval([name '(RR,CC) = mean(convar,2);']);
        end % if mymat...

		%eval(['save ' name num2str(index) ' ' name ' SIM'])
		index = index + 1;
	end
end

eval(['save ' name])
figure;
set(gcf,'Color','w')
eval(['surf(' name ')']); colormap(copper)
eval(['set(gca,''XTick'',1:size(' name ',2))'])
set(gca,'XTickLabel',myFC)
%eval(['set(gca,''YTick'',1:size(' name ',1))'])
%set(gca,'YTickLabel',myISI)
xlabel('Proportion of blank intervals','FontSize',14)
ylabel('ISI in s','FontSize',14)
title('Surface map of ISI vs. proportion of blank intervals','FontSize',18)
zlabel('Efficiency of contrast [1 -1]','FontSize',14)

set(gca,'YTick',ISItick)
set(gca,'YTickLabel',ISIticklab)

warning on

break
%for line plots
figure;plot(MAP1)
plot(MAP1,'LineWidth',2)
legend({'.01' '.21' '.41' '.61' '.81'})
set(gca,'XTick',1:2:40)
set(gca,'XTickLabels',.1:.4:7.9)
title('Efficiency as a function of ISI and jitter: Linear HRF function','FontSize',14)
grid on

eval(['saveas(gcf,''' name ''',''fig'')'])
eval(['saveas(gcf,''' name ''',''jpg'')'])

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function triangulation_orient ( prefix )

%*****************************************************************************80
%
%% MAIN is the main program for TRIANGULATION_ORIENT.
%
%  Discussion:
%
%    TRIANGULATION_ORIENT forces a triangulation to have positive orientation.
%
%    The user supplies a node file and a triangle file, containing
%    the coordinates of the nodes, and the indices of the nodes that
%    make up each triangle.  Either 3-node or 6-node triangles may
%    be used.
%
%    The program reads the data, and for each triangle, determines
%    whether the triangle has positive orientation.  This essentially
%    means that the vertices are listed in counter clockwise order.
%    If the vertices are listed in the wrong order, they are reordered.
%    The reordered triangle file is written out.
%
%    Note that for an order 6 triangulation, the vertices are listed
%    in the first three positions.
%
%  Usage:
%
%    triangulation_orient ( 'prefix' )
%
%    where 'prefix' is the common filename prefix:
%
%    * 'prefix'_nodes.txt contains the node coordinates,
%    * 'prefix'_elements.txt contains the element definitions.
%    * 'prefix'_orient_elements.txt will contain the oriented element definitions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_ORIENT\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Read a node dataset of NODE_NUM points in 2 dimensions.\n' );
  fprintf ( 1, '  Read an associated triangle file of TRIANGLE_NUM \n' );
  fprintf ( 1, '  triangles using 3 or 6 nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ensure that every triangle has positive orientation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Write the reoriented triangle file.\n' );
%
%  The command line argument is the common filename prefix.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORIENT:\n' );

    prefix = input ( ...
      'Please enter the filename prefix:' );

  end
%
%  Create the filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
  element_orient_filename = strcat ( prefix, '_orient_elements.txt' );
%
%  Read the node data.
%
  [ dim_num, node_num ] = r8mat_header_read (  node_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of points NODE_NUM = %d\n', node_num );

  node_xy(1:dim_num,1:node_num) = r8mat_data_read ( node_filename, ...
    dim_num, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', node_filename );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, dim_num, 5, ...
    '  First 5 nodes:' );
%
%  Read the element data.
%
  [ triangle_order, triangle_num ] = i4mat_header_read ( ...
    element_filename );

  if ( triangle_order ~= 3 && triangle_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORIENT - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 3-node or 6-node triangulation.\n' );
    error ( 'TRIANGULATION_L2Q - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', element_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Triangle order = %d\n', triangle_order );
  fprintf ( 1, '  Number of triangles TRIANGLE_NUM  = %d\n', triangle_num );

  triangle_node(1:triangle_order,1:triangle_num) = i4mat_data_read ( ...
    element_filename, triangle_order, triangle_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', element_filename );

  i4mat_transpose_print_some ( triangle_order, triangle_num, triangle_node, ...
    1, 1, triangle_order, 5, '  First 5 triangles:' );
%
%  Detect and correct 0-based indexing.
%
  triangle_node = mesh_base_one ( node_num, triangle_order, triangle_num, ...
    triangle_node );
%
%  Compute the area, and reorient if necessary.
%
  triangle_negative_num = 0;
  triangle_zero_num = 0;

  for triangle = 1 : triangle_num

    t3(1:2,1:3) = node_xy(1:2,triangle_node(1:3,triangle));

    area = triangle_area_2d ( t3 );

    if ( area < 0.0 )

      triangle_negative_num = triangle_negative_num + 1;

      node                      = triangle_node(2,triangle);
      triangle_node(2,triangle) = triangle_node(3,triangle);
      triangle_node(3,triangle) = node;

      if ( triangle_order == 6 )
        node                      = triangle_node(4,triangle);
        triangle_node(4,triangle) = triangle_node(6,triangle);
        triangle_node(6,triangle) = node;
      end
%
%  As a check, repeat the area calculation.
%  Now, we expect to get a positive result.
%
      t3(1:2,1:3) = node_xy(1:2,triangle_node(1:3,triangle));

      area = triangle_area_2d ( t3 );

      if ( area < 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'TRIANGULATION_ORIENT - Fatal error!\n' );
        fprintf ( 1, '  I thought I fixed the area but I did not!\n' );
        error ( 'TRIANGULATION_L2Q - Fatal error!' );
      end

    elseif ( area == 0.0 )

      triangle_zero_num = triangle_zero_num + 1;

    end

  end

  if ( 0 < triangle_zero_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORIENT - Warning!\n' );
    fprintf ( 1, '  You have %d triangles with\n', triangle_zero_num );
    fprintf ( 1, '  area equal to zero.\n' );
  end

  if ( 0 < triangle_negative_num )

    i4mat_write ( element_orient_filename, triangle_order, ...
      triangle_num, triangle_node );

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORIENT - Warning!\n' );
    fprintf ( 1, '  You have %d triangles with\n', triangle_negative_num );
    fprintf ( 1, '  negative area.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  We have reoriented these triangles to have positive\n' );
    fprintf ( 1,'  area, and written the new triangle data to \n' );
    fprintf ( 1, '  the triangle file "%s".\n', element_orient_filename );

    i4mat_transpose_print_some ( triangle_order, triangle_num, ...
      triangle_node, 1, 1, triangle_order, 5, '  First 5 triangles:' );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORIENT:\n' );
    fprintf ( 1, '  None of your triangles had negative area.\n' );
    fprintf ( 1, '  Therefore, no new triangle file was written.\n' );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_ORIENT\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function element_node = mesh_base_one ( node_num, element_order, ...
  element_num, element_node )

%*****************************************************************************80
%
%% MESH_BASE_ONE ensures that the element definition is one-based.
%
%  Discussion:
%
%    The ELEMENT_NODE array contains nodes indices that form elements.
%    The convention for node indexing might start at 0 or at 1.
%
%    Since a MATLAB program will naturally assume a 1-based indexing, it is
%    necessary to check a given element definition and, if it is actually
%    0-based, to convert it.
%
%    This function attempts to detect 0-based node indexing and correct it.
%
%    Thanks to Feifei Xu for pointing out that I was subtracting 1 when I
%    should have been adding 1!  29 November 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input/output, integer ELEMENT_NODE(ELEMENT_ORDE,ELEMENT_NUM), the element
%    definitions.
%
  node_min = min ( min ( element_node(1:element_order,1:element_num) ) );
  node_max = max ( max ( element_node(1:element_order,1:element_num) ) );

  if ( node_min == 0 && node_max == node_num - 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 0-based!\n' );
    fprintf ( 1, '  This will be converted to 1-based.\n' );
    element_node(1:element_order,1:element_num) = ...
      element_node(1:element_order,1:element_num) + 1;
  elseif ( node_min == 1 && node_max == node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 1-based!\n' );
    fprintf ( 1, '  No conversion is necessary.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE - Warning!\n' );
    fprintf ( 1, '  The element indexing is not of a recognized type.\n' );
    fprintf ( 1, '  NODE_MIN = %d\n', node_min );
    fprintf ( 1, '  NODE_MAX = %d\n', node_max );
    fprintf ( 1, '  NODE_NUM = %d\n', node_num );
  end

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function area = triangle_area_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_AREA_2D computes the area of a triangle in 2D.
%
%  Discussion:
%
%    If the triangle's vertices are given in counterclockwise order,
%    the area will be positive.  If the triangle's vertices are given
%    in clockwise order, the area will be negative!
%
%    An earlier version of this routine always returned the absolute
%    value of the computed area.  I am convinced now that that is
%    a less useful result!  For instance, by returning the signed
%    area of a triangle, it is possible to easily compute the area
%    of a nonconvex polygon as the sum of the (possibly negative)
%    areas of triangles formed by node 1 and successive pairs of vertices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real AREA, the area of the triangle.
%
  area = 0.5 * ( ...
      t(1,1) * ( t(2,2) - t(2,3) ) ...
    + t(1,2) * ( t(2,3) - t(2,1) ) ...
    + t(1,3) * ( t(2,1) - t(2,2) ) );

  return
end


text:
function glm = glmCreate(tr,nh,whiten)
%Create a GLM model for managing the GLM parameters
%
%   glm = glmCreate(tr,nh,whiten)
%
%

% It would be better to go and get the TR from the SESSION data.
if notDefined('tr'),     error('You must specify tr for the GLM'); end
if notDefined('nh'),     nh = 22; end
if notDefined('whiten'), whiten = 0; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize GLM model %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
glm.betas    = [];
glm.residual = [];
glm.stdevs   = [];
glm.sems     = [];
glm.nh       = nh;       % Number of time samples
glm.tr       = tr;       %
glm.C        = [];       % Covariance matrix
glm.voxDepNoise    = [];
glm.voxIndepNoise  = [];
glm.designMatrix   = [];
glm.dof            = []; % Degrees of freedom
glm.whiten         = whiten;  % Do not whiten the covariance

% glmSet() ...
if nh>1, glm.type = 'selective averaging';
else     glm.type = 'applied HRF';
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft