text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)

% This file shows how to perform training on your own dataset.

% 1. Training demo on PASCAL VOC 2012 dataset

% We here demonstrate the training on the original PASCAL VOC 2012 dataset for semantic segmentation.
% Before running this demo file for training, you need to
% download the PASCAL VOC 2012 dataset and place it at:
% ../datasets/voc2012_trainval

% 2. Training on other dataset:

% other datasets can be handled in a similar way as the VOC2012 dataset.
% Basically the users need to provide a file to describe your dataset:
% my_gen_ds_info_[your dataset name].m
% please refer to my_gen_ds_info_voc.m for details

% 3. Training instructions:

% Stage-1: using a fixed learning rate for training until the performance
% get stable on the validation set. You need to manually stop the running of this training code.

% Stage-2: chose one cached model, initialize from this model, and using a
% lower learning rate (e.g., multiplied by 0.1) to perform further training.

% For example, using the following setting for loading an existed model and
% using a lower learning rate:
% run_config.trained_model_path='../cache_data/voc2012_trainval/model_20161219094311_example/model_cache/epoch_70';
% run_config.learning_rate=5e-5;

% Please refer to the following demo file for perform further training with lower learning rate:
% demo_refinenet_train_reduce_learning_rate.m

% 4. Verify your trained model with performance evaluation:

% This demo file shows the training on the original PASCAL VOC 2012 dataset for semantic segmentation.
% This dataset consists of 1464 training images and 1449 validation images.
% For this dataset, the training will take around 3 days using a Titan-X card,
% including the training with a decreased learning rate.

% After training, the intersection over union (IoU) score of the trained model is supposed to be around 0.691.

% If you cannot get similar result, please verify whether you are using our modified version of MatConvNet,
% for which you need to replace the modified files in the original MatConvNet and compile after replacement.
% Details can be found at ./my_matconvnet_resnet/README.txt

function demo_refinenet_train()

rng('shuffle');

addpath('./my_utils');
dir_matConvNet='../libs/matconvnet/matlab';
run(fullfile(dir_matConvNet, 'vl_setupnn.m'));

run_config=[];
ds_config=[];

run_config.use_gpu=true;
% run_config.use_gpu=false;
run_config.gpu_idx=1;

% use a random model name:
model_name=['model_' datestr(now, 'YYYYmmDDHHMMSS')];

% or specify a model_name for training, if the cached snapshot files existed, it will resume training
% model_name='model_20161219094311_example';

% update the model name when reducing learning rate for training, e.g.,
% model_name='model_20161219094311_example_epoch70_low_learn_rate';

ds_name='voc2012_trainval';
gen_ds_info_fn=@my_gen_ds_info_voc;

% control the size of input images
run_config.input_img_short_edge_min=450;
run_config.input_img_short_edge_max=1100;

run_config.input_img_scale=1;

ds_config.use_dummy_gt=false;
run_config.use_dummy_gt=ds_config.use_dummy_gt;
ds_config.use_custom_data=false;

ds_config.ds_name=ds_name;
ds_config.gen_ds_info_fn=gen_ds_info_fn;
ds_config.ds_info_cache_dir=fullfile('../datasets', ds_name);

run_config.root_cache_dir=fullfile('../cache_data', ds_name, model_name);
mkdir_notexist(run_config.root_cache_dir);

run_config.model_name=model_name;

diary_dir=run_config.root_cache_dir;
mkdir_notexist(diary_dir);
diary(fullfile(diary_dir, 'output.txt'));
diary on

run_dir_name=fileparts(mfilename('fullpath'));
[~, run_dir_name]=fileparts(run_dir_name);
run_config.run_dir_name=run_dir_name;
run_config.run_file_name=mfilename();

fprintf('\n\n\n=====================================================================\n');
disp('gen ds_info');

disp('ds_config:');
disp(ds_config);

ds_info=gen_dataset_info(ds_config);

my_diary_flush();

fprintf('\n\n\n=====================================================================\n');
disp('run network');

run_config.run_evaonly=false;

% settings for training:

run_config.trained_model_path=[];
run_config.learning_rate=5e-4;

% init from a cached model
% e.g., can be used for resuming network training with a lower learning rate by
% initilizing from a cached model.
% the following is an example for loading a cached model and using a lower learning
% to continue the training

% please refer to the example script: demo_refinenet_train_reduce_learning_rate.m

% run_config.trained_model_path='../cache_data/voc2012_trainval/model_20161219094311_example/model_cache/epoch_70';
% run_config.learning_rate=5e-5;

run_config.epoch_num=600;

% turn on this option to cache all data into memory, if it's possible
% run_config.cache_data_mem=true;
run_config.cache_data_mem=false;

% random crop training:
run_config.crop_box_size=400;

% for cityscape, using a larger crop:
% run_config.crop_box_size=600;

% evaluate step: do evaluation every 5 epochs, can be set to 10:
run_config.eva_run_step=5;
run_config.snapshot_step=1;

% choose ImageNet pre-trained resnet:
run_config.init_resnet_layer_num=50;
% run_config.init_resnet_layer_num=101;
% run_config.init_resnet_layer_num=152;

% generate network
run_config.gen_network_fn=@gen_network_main;

run_config.gen_net_opts_fn=@gen_net_opts_model_type1;
% run_config.gen_net_opts_fn=@gen_net_opts_model_type2;

% uncomment the following for debug:
% select a subset for both training and evaluation.
% ds_info.train_idxes=ds_info.train_idxes(1:5);
% ds_info.test_idxes=ds_info.train_idxes;
% run_config.snapshot_step=inf;

train_opts=run_config.gen_net_opts_fn(run_config, ds_info.class_info);

disp('train_opts:');
disp(train_opts);

disp('train_opts.eva_param:');
disp(train_opts.eva_param);

my_diary_flush();

imdb=my_gen_imdb(train_opts, ds_info);

data_norm_info=[];
data_norm_info.image_mean=128;

imdb.ref.data_norm_info=data_norm_info;

if run_config.use_gpu
	gpu_num=gpuDeviceCount;
	if gpu_num>=1
		gpuDevice(run_config.gpu_idx);
    else
        error('no gpu found!');
	end
end

[net_config, net_exp_info]=prepare_running_model(train_opts);

% net_config can be changed here before running the network

my_net_tool(train_opts, imdb, net_config, net_exp_info);

fprintf('\n\n--------------------------------------------------\n\n');
disp('results are saved in:');
disp(run_config.root_cache_dir);

my_diary_flush();
diary off

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function out = grandAverage_prototye(cell_struct, bool_arith)
%% Description
% cell of averaged data
numData = length(cell_struct);

if nargin == 1
    bool_arith = false;
end

for i = 1:numData
    dat = cell_struct{i};

    if isequal(dat.class{1,2}, 'sgnr^2')
        dat.x = atanh(sqrt(abs(dat.x)).*sign(dat.x));
    end
    cell_struct{i} = dat;
end
[time, trials, channels] = size(cell_struct{1}.x);
dat_av = zeros(time, trials, channels);

ncls = size(cell_struct{1}.class, 1);
for cls = 1:ncls
    sW = 0;
    swV = 0;
    for v = 1:numData
        if bool_arith %% options-> Arithmetic or WeightedMean
            W = 1./cell_struct{v}.se(:, cls,:).^2;
        else
            % Arithmetic
            W = 1;
        end
        sW = sW + W;
        swV = swV + W.^2.*cell_struct{v}.se(:, cls, :).^2;
        dat_av(:, cls, :) = dat_av(:, cls,:) + W.*cell_struct{v}.x(:,cls,:);
    end
    dat_av(:, cls,:) = dat_av(:, cls,:)./sW;
    se(:, cls,:) = sqrt(swV)./sW;
end

if isequal(dat.class{1,2}, 'sgnr^2')
    dat_av = tanh(dat_av).*abs(tanh(dat_av));
end

out = cell_struct{1};
out.x = dat_av;
out.se = se;
end

author: PatternRecognition

repo: OpenBMI


text:
function [ m, n, nz_num, base ] = r8cc_read_size ( col_file, row_file )

%*****************************************************************************80
%
%% R8CC_READ_SIZE reads the sizes of a R8CC sparse matrix from a file.
%
%  Discussion:
%
%    The value of M is "guessed" to be the largest value that occurs in
%    the ROW file.  However, if a row index of 0 is encountered, then
%    the value of M is incremented by 1.
%
%    The value of N is the number of records in the COL file minus 1.
%
%    The value of NZ_NUM is simply the number of records in the ROW file.
%
%    The value of BASE is 0 or 1, depending on whether the program
%    "guesses" that the row and column indices are 0-based or 1-based.
%    Although the first entry of the COL array might be used as evidence,
%    this program makes its determination based on whether it encounters
%    a 0 index in the ROW file.
%
%    The R8CC format is the double precision sparse compressed column
%    format.  Associated with this format, we have an M by N matrix
%    with NZ_NUM nonzero entries.  We construct the column pointer
%    vector COL of length N+1, such that entries of column J will be
%    stored in positions COL(J) through COL(J+1)-1.  This indexing
%    refers to both the ROW and A vectors, which store the row indices
%    and the values of the nonzero entries.  The entries of the
%    ROW vector corresponding to each column are assumed to be
%    ascending sorted.
%
%    The R8CC format is equivalent to the MATLAB "sparse" format,
%    and the Harwell Boeing "real unsymmetric assembled" (RUA) format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992
%
%  Parameters:
%
%    Input, string COL_FILE, ROW_FILE, the names of the
%    column and row files that describe the structure of the matrix.
%
%    Output, integer M, N, the inferred number of rows and columns
%    in the sparse matrix.
%
%    Output, integer NZ_NUM, the number of nonzero entries in the
%    sparse matrix.
%
%    Output, integer BASE, is 0 if the row indexing is believed
%    to be 0-based, and 1 if the row-index is believed to be
%    1-based.  In uncertain cases, BASE = 1 is the default.
%

%
%  Default values.
%
  m = -1;
  n = -1;
  nz_num = -1;
  base = -1;
%
%  Check the COL file first.
%
  input_unit = fopen ( col_file );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8CC_READ_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', col_file );
    error ( 'R8CC_READ_SIZE - Fatal error!' );
    return;
  end

  n = -1;

  while ( 1 )

    [ col, count ] = fscanf ( input_unit, '%d', 1 );

    if ( count ~= 1 )
      break;
    end

    n = n + 1;

  end

  fclose ( input_unit );
%
%  Check the ROW file.
%
  input_unit = fopen ( row_file );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8CC_READ_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', row_file );
    error ( 'R8CC_READ_SIZE - Fatal error!' );
    return;
  end

  base = 1;
  m = 0;
  nz_num = 0;

  while ( 1 )

    [ row, count ] = fscanf ( input_unit, '%d', 1 );

    if ( count ~= 1 )
      break;
    end

    nz_num = nz_num + 1;
    m = max ( m, row );
    if ( row == 0 )
      base = 0;
    end

  end

  fclose ( input_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ADEM_observe
% This demo illustrates action-observation using synthetic writing under
% active inference. It shows how expectations about hidden states can be
% both cause and consequence of observed action (of self and others
% respectively). We first illustrate the generation of behaviour using a
% Lotka-Volterra form stable heteroclinic orbit. We then reproduce the
% same forces on the agent's arm but switching off the precision of
% proprioceptive inputs. This can be seen as attending selectively to
% visual inputs. The resulting inference calls upon the same hidden-states
% and implicit predictions (in a generalised or dynamic sense). These
% simulations can be regarded as simulations of mirror neuron responses.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_observe.m 7679 2019-10-24 15:54:07Z spm $

% hidden causes and states
%==========================================================================
% x    - hidden states
%   x.x(1) - joint angle
%   x.x(2) - joint angle
%   x.x(3) - angular velocity
%   x.x(4) - angular velocity
%
%   x.a(1) - attraction (location 1)
%   x.a(2) - attraction (location 2)
%   x.a(3) - attraction (location 3)
%    ...
%
% v    - causal states
%   v(1) - not used
%
%--------------------------------------------------------------------------

% parameters mapping from (unstable) point attractors to visual space
%--------------------------------------------------------------------------
P = [1.0  1.0;
     1.1  1.2;
     1.0  0.4;
     1.0  1.0;
     1.4  0.9;
     0.9  1.0]';
n = size(P,2);                                % number of attractors

% Recognition model (linear for expediency)
%==========================================================================
M(1).E.s = 1/2;                               % smoothness
M(1).E.n = 4;                                 % order of
M(1).E.d = 2;                                 % generalised motion

% level 1: Displacement dynamics and mapping to sensory/proprioception
%--------------------------------------------------------------------------
M(1).f   = 'spm_fx_dem_observe';              % plant dynamics
M(1).g   = 'spm_gx_dem_write';                % prediction

M(1).x.x = [pi/2; pi/2; 0; 0];                % physical states
M(1).x.a = sparse(1,1,3,n,1) - 4;             % attractor states
M(1).V   = exp(4);                            % error precision
M(1).W   = exp(8);                            % error precision
M(1).pE  = P;

% level 2: not used
%--------------------------------------------------------------------------
M(2).v  = 0;                                  % inputs
M(2).V  = exp(8);

% generative model
%==========================================================================

% first level
%--------------------------------------------------------------------------
G(1).f  = 'spm_fx_adem_write';
G(1).g  = 'spm_gx_adem_write';
G(1).x  = [pi/2; pi/2; 0; 0];                 % physical states
G(1).V  = exp(16);                            % error precision
G(1).W  = exp(16);                            % error precision
G(1).U  = sparse(1:4,1:4,exp(4),8,8);         % restriction

% second level
%--------------------------------------------------------------------------
G(2).v  = [0; 0];                             % exogenous forces
G(2).a  = [0; 0];                             % action forces
G(2).V  = exp(16);

% generate and invert
%==========================================================================
N       = 256;                                % length of data sequence
t       = (1:N)*8;
DEM.G   = G;
DEM.M   = M;
DEM.C   = sparse(2,N);
DEM     = spm_ADEM(DEM);

% overlay true values
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');
spm_DEM_qU(DEM.qU,DEM.pU)

subplot(2,2,3)
spm_dem_reach_movie(DEM)
title('click on finger for movie','FontSize',16)

% Action-observation
%==========================================================================
ADEM          = DEM;
ADEM.C        = DEM.qU.a{2};
ADEM.M(1).V   = exp([-8 -8 -8 -8  8  8  8  8]);  % remove proprioception
ADEM.G(1).U   = 0;                               % remove motor control
ADEM.M(1).x.x = DEM.M(1).x.x;
ADEM.M(1).x.a = DEM.M(1).x.a*0;
ADEM          = spm_ADEM(ADEM);

spm_figure('GetWin','Figure 2');
spm_DEM_qU(ADEM.qU,ADEM.pU)

subplot(2,2,3)
spm_dem_reach_movie(ADEM)
title('click on finger for movie','FontSize',16)

% Hidden (intended) states lead actual states
%==========================================================================
spm_figure('GetWin','Figure 3');

n  = 64;
p  = exp(DEM.qU.x{1}(5:end,:));
qx = P(end,:)*(p*diag(1./sum(p)));
px = DEM.pU.v{1}(end,:);

subplot(2,1,1)
plot(t,px,t,qx,'-.')
legend({'attained','intended'})
title('intended and attained position','FontSize',16)
xlabel('time')
ylabel('vertical displacement')
box off

qx = qx - mean(qx);
px = px - mean(px);

subplot(2,2,3)
xc = xcorr(qx,px,n,'coef');
plot([0:(n + n)] - n,xc,[0 0],[-1 1],'k-.')
title('action','FontSize',16)
xlabel('time')
ylabel('cross-correlation')
axis square tight

% repeat for inferred hidden (intended) states
%--------------------------------------------------------------------------
p  = exp(ADEM.qU.x{1}(5:end,:));
qx = P(end,:)*(p*diag(1./sum(p)));
px = ADEM.pU.v{1}(end,:);

qx = qx - mean(qx);
px = px - mean(px);

subplot(2,2,4)
xca = xcorr(qx,px,n,'coef');
plot((0:(n + n)) - n,xca,(0:(n + n)) - n,xc,':',[0 0],[-1 1],'k-.')
title('observation','FontSize',16)
xlabel('time')
ylabel('cross-correlation')
axis square tight

% coherence analysis
%==========================================================================
spm_figure('GetWin','Figure 4');

T  = 16;
xe = DEM.qU.w{1}(5:end,T:end);      % prediction error on attractor states
ve = DEM.qU.z{1}(4:8,T:end);        % prediction error on visual states
qx = full(sum(xe));                 % synthetic LFP
qv = full(sum(ve));                 % synthetic LFP
px = DEM.pU.v{1}(1,T:end);          % peripheral (plant motion)

% coherence analysis
%--------------------------------------------------------------------------
qx = spm_detrend(qx',2);
px = spm_detrend(px',2);

[Cqp,Hz] = mscohere(qx,px,64,48,[],1/0.008);
[Cxv,Hz] = mscohere(qx,qv,64,48,[],1/0.008);

subplot(2,1,1)
plot(t(T:end),xe,'r',t(T:end),ve,'b');  hold on
plot(t(T:end),px/8,'g','LineWidth',4); hold off
title('central and peripheral responses','FontSize',16)
xlabel('time (ms)')
ylabel('activity (au)')
axis tight

subplot(2,2,3)
plot(Hz,Cqp)
title('central-peripheral coherence','FontSize',16)
xlabel('frequency (Hz)')
ylabel('coherence')
axis square
axis([0 32 0 1])

subplot(2,2,4)
plot(Hz,Cxv)
title('central-central coherence','FontSize',16)
xlabel('frequency (Hz)')
ylabel('coherence')
axis square
axis([0 32 0 1])

% functional correlates
%==========================================================================
spm_figure('GetWin','Figure 5');

% under action
%--------------------------------------------------------------------------
x    = DEM.pU.v{1}(7,:);
y    = DEM.pU.v{1}(8,:);

u    = DEM.qU.x{1}(8,:);
i    = find(u > 2);

subplot(2,2,1)
plot(x,y,'LineWidth',4,'Color',[1 1 1]*.8), hold on
plot(x(i),y(i),'r.','MarkerSize',24),       hold off
title('action','FontSize',16)
xlabel('position (x)')
ylabel('position (y)')
axis square

x  = x + ([1:N] - N)/N;

subplot(2,2,3)
plot(x,y,'LineWidth',4,'Color',[1 1 1]*.8), hold on
plot(x(i),y(i),'r.','MarkerSize',24),       hold off
title('action','FontSize',16)
xlabel('position (x)')
ylabel('position (y)')
axis square ij

% and observation
%--------------------------------------------------------------------------
x    = ADEM.pU.v{1}(7,:);
y    = ADEM.pU.v{1}(8,:);

u    = ADEM.qU.x{1}(8,:);
i    = find(u > 2);

subplot(2,2,2)
plot(x,y,'LineWidth',4,'Color',[1 1 1]*.8), hold on
plot(x(i),y(i),'r.','MarkerSize',24),       hold off
title('observation','FontSize',16)
xlabel('position (x)')
ylabel('position (y)')
axis square

x  = x + ([1:N] - N)/N;

subplot(2,2,4)
plot(x,y,'LineWidth',4,'Color',[1 1 1]*.8), hold on
plot(x(i),y(i),'r.','MarkerSize',24),       hold off
title('observation','FontSize',16)
xlabel('position (x)')
ylabel('position (y)')
axis square ij

% Correlations between action and observation responses
%==========================================================================
spm_figure('GetWin','Figure 6');

qx =  DEM.qU.x{1};
px = ADEM.qU.x{1};

subplot(2,1,1)
bar3(spm_en(px',0)'*spm_en(qx',0))
title('correlations','FontSize',16)
xlabel('hidden unit (observation)')
ylabel('hidden unit (action)')
grid off

% last hidden state
%--------------------------------------------------------------------------
qx =  DEM.qU.x{1}(end,:);
px = ADEM.qU.x{1}(end,:);

subplot(2,2,3)
plot(px,qx,'.')
title('correlations','FontSize',16)
xlabel('hidden unit (observation)')
ylabel('hidden unit (action)')
box off
axis square

qx = spm_en(qx',0);
px = spm_en(px',0);

subplot(2,2,4)
xc = xcorr(qx,px,n,'coef');
plot([0:(n + n)] - n,xc,[0 0],[-1 1],'k-.')
title('action','FontSize',16)
xlabel('time')
ylabel('cross-correlation')
axis square tight

% Simulated ERP
%==========================================================================

% simulate deviant
%--------------------------------------------------------------------------
T          = 132;
C          = DEM.qU.a{2};
C(:,T:end) = -C(:,T:end);
NDEM       = ADEM;
NDEM.C     = C;
NDEM       = spm_ADEM(NDEM);

% Graphics
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 7');
t   = ([1:N] - T)*8;
i   = find(t > -256 & t < 512);

subplot(3,2,1)
spm_dem_reach_movie(ADEM)
title('observation','FontSize',16)

subplot(3,2,2)
spm_dem_reach_movie(NDEM)
title('violation','FontSize',16)

% Proprioceptive prediction error
%--------------------------------------------------------------------------
subplot(3,2,4)
plot(t(i),NDEM.qU.z{1}(1:4,i)'), hold on
plot(t(i),ADEM.qU.z{1}(1:4,i)',':'), hold off
title('proprioceptive error','FontSize',16)
xlabel('time')
ylabel('prediction error')
axis square
a = axis;

subplot(3,2,3)
plot(t(i),ADEM.qU.z{1}(1:4,i)')
title('proprioceptive error','FontSize',16)
xlabel('time')
ylabel('prediction error')
axis square
axis(a);

% Hidden state prediction error
%--------------------------------------------------------------------------
subplot(3,2,6)
plot(t(i),NDEM.qU.w{1}(:,i)'), hold on
plot(t(i),ADEM.qU.w{1}(:,i)',':'), hold off
title('error on hidden states','FontSize',16)
xlabel('time')
ylabel('prediction error')
axis square
a = axis;

subplot(3,2,5)
plot(t(i),ADEM.qU.w{1}(:,i)')
title('error on hidden states','FontSize',16)
xlabel('time')
ylabel('prediction error')
axis square
axis(a);

author: spm


text:
function [segclassifier] = mcmcTrainSegmentationClassifier2(features, labels, weights, maxdata, classparams)

if exist('classparams') && ~isempty(classparams)
    nnodes = classparams(1);
    ntrees = classparams(2);
    stopval = classparams(3);
else
    ntrees = 20;
    nnodes = 8;
    stopval = 0;
end

if ~exist('maxdata') || isempty(maxdata)
    maxdata = 25000;
end

[data, lab, w] = formatData(features, labels, weights, maxdata);

segclassifier = train_boosted_dt_2c(data, [], lab, ntrees, nnodes, stopval, w);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [data, lab, w] = formatData(features, labels, weights, maxdata)
% concatenate data
%disp(num2str(maxdata))
nimages = numel(features);

[tmp, nvars] = size(features{1});

% count segments
tseg = 0;
nseg = 0;
for f = 1:nimages
    nseg = nseg + sum((labels{f}~=0) & (weights{f}>0));
    tseg = tseg + numel(labels{f});
end
%disp(num2str([nseg tseg]))

data = zeros(nseg, nvars);
lab = zeros(nseg, 1);
w = zeros(nseg, 1);

% concatenate data
vc = 0;
for f = 1:nimages
    ind = find((labels{f}~=0) & (weights{f}>0));
    data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
    lab(vc+1:vc+numel(ind)) = (labels{f}(ind)>0)*2-1;
    w(vc+1:vc+numel(ind)) = weights{f}(ind); % weight according to area in image
    vc = vc + numel(ind);
end

if nseg > maxdata
    rind = randperm(nseg);
    rind = rind(1:maxdata);
    data = data(rind, :);
    lab = lab(rind);
    w = w(rind);
end

w = w / sum(w);

author: Cloud-CV

repo: object-proposals


text:
function f1 = p09_f1 ( x )

%*****************************************************************************80
%
%% P09_F1 evaluates the first derivative for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  f1 = 2.0 * x ...
    - 10.0 * cos ( x * x - 3.0 * x + 2.0 ) ...
    * ( 2.0 * x - 3.0 );

  return
end


text:
function [zeitgeberTime relativeTime absTime] = bz_getZeitgeberTime(basename,downSampleFactor)
% assumes pwd is basepath for now

%% get timestamps for recording starts
parts = strsplit(basename, '_');
date = parts{end-1};

if strcmp(parts{end},'merge') % if multiple sessions are merged, find the start times of each
    fileList = dir('*');

    for f = 1:length(fileList)
        if fileList(f).isdir
        for p = 1:length(parts)
            if ~isempty(strfind(fileList(f).name,parts{p}))
               matches(f,p) = 1;
            else
               matches(f,p) = 0;
            end
        end
        end
    end
    idx = find(sum(matches')==length(parts)-1);  % exlude '_merge' string
    for ind = 1:length(idx)
        temp = strsplit(fileList(idx(ind)).name,'_');
        recordingStart{ind} = temp{end};
    end
else
    recordingStart{1} = parts{end};
end

for ind = 1:length(recordingStart)
timeSeconds(ind) = int32(str2num(recordingStart{ind}(1:2)) * 60 * 60 + ...
                   str2num(recordingStart{ind}(3:4)) * 60 + ...
                   str2num(recordingStart{ind}(5:6)));
end

%% now get relative time from time.dat file
fileinfo = dir('time.dat');
num_samples = fileinfo.bytes/4; % int32 = 4 bytes
fid = fopen('time.dat', 'r');
t = fread(fid, num_samples, 'int32=>int32')';
t = downsample(t,downSampleFactor);
fclose(fid);

[a b] = min(t);  %% fixes output for time.dat files that 'wrap' at the int32 limit
if a < 0
    t = double(t)./20000;
    for i=b:length(t)
    t(i) = t(i-1)+1;
    end
else
    t = single(t)./20000;
end
% t = t / 20000;

idx = find(t==0);
idx(length(idx)+1) = length(t); % add end ts

for ind = 1:length(idx)-1
    zeitgeberTime(idx(ind):idx(ind+1)) = single(t(idx(ind):idx(ind+1)))/20000 + single(timeSeconds(ind));
end

zeitgeberTime = single(wrap(double(zeitgeberTime)/86400*pi*2,2))/ (2*pi) * 86400; % wrap it up

relativeTime = t;  % recording time

% zeitgeberTime = downsample(zeitgeberTime,downSampleFactor);
% relativeTime = downsample(relativeTime,downSampleFactor

% below is error prone if gaps in relative time, but WAY faster
absTime = linspace(datenum(2000 + str2num(date(1:2)),str2num(date(3:4)),str2num(date(5:6)),...
        str2num(parts{end}(1:2)),str2num(parts{end}(3:4)),double(relativeTime(1))),...
        datenum(2000 + str2num(date(1:2)),str2num(date(3:4)),str2num(date(5:6)),...
        str2num(parts{end}(1:2)),str2num(parts{end}(3:4)),double(relativeTime(end))),length(relativeTime));

% switch to this for accuracy
% for ind = 1:length(relativeTime)
%     absTime(ind) = datenum(2000 + str2num(date(1:2)),str2num(date(3:4)),str2num(date(5:6)),...
%         str2num(parts{end}(1:2)),str2num(parts{end}(3:4)),double(relativeTime(ind)));
% end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
data_masterdir = GetCurrentDataDir();

range_fish = [5,6,7];
M_ClusGroup = [2,2,2,2];
M_Cluster = [1,1,1,1];

% M_fish_set = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2];

%%
for i = 1:length(range_fish),
    i_fish = range_fish(i);
    disp(i_fish);

    LoadFullFish(hfig,i_fish);

    absIX = getappdata(hfig,'absIX');

    %% Cluster indexing
    i_ClusGroup = M_ClusGroup(i);
    i_Cluster = M_Cluster(i);
    Cluster = VAR(i_fish).ClusGroup{i_ClusGroup};
    numK = Cluster(i_Cluster).numK;
    gIX = Cluster(i_Cluster).gIX;

    cIX_abs = Cluster(i_Cluster).cIX_abs; % convert absolute index to index used for this dataset
    [~,cIX] = ismember(cIX_abs,absIX);
    setappdata(hfig,'cIX',cIX);
    toc

    %%
    periods = getappdata(hfig,'periods');
    if length(periods)>1,
        setappdata(hfig,'stimrange',1:length(periods));
    else
        setappdata(hfig,'stimrange',1);
    end
    UpdateTimeIndex(hfig); % set M_0
    M_0 = getappdata(hfig,'M_0');

    %%
    f.UpdateIndices(hfig,cIX,gIX,numK);

    %%
    isWkmeans = 1;
    AutoClustering(cIX,gIX,absIX,i_fish,M_0,isWkmeans)
end

author: xiuyechen

repo: FishExplorer


text:
function [values, components, full_data_objects, l2norms] = extract_gray_white_csf(obj,varargin)
% Extracts mean values (values) and top 5 component scores (components)
% from each of gray, white, and CSF masks.
%
% - Images must be in standard MNI space for this to apply.
% - Uses canonical masks in CANlab tools:'gray_matter_mask.img' 'canonical_white_matter.img' 'canonical_ventricles.img'
%
% This currently uses the images
%       'gray_matter_mask.img' 'canonical_white_matter.img' 'canonical_ventricles.img'
%       These images are based on the SPM8 a priori tissue probability
%       maps, but they have been cleaned up and made symmetrical and/or eroded
%       so that the white and CSF compartments are unlikely to contain very
%       much gray matter.  The gray compartment is currently more
%       inclusive. The potential value of this is that signal in the CSF/white
%       compartments may be removed from images prior to/during analysis
%
% :Usage:
% ::
%
%     [values, components] = extract_gray_white_csf(obj,options)
%
% :Inputs:
%
%   **obj:**
%        an image_vector (e.g., fmri_data) object
%
% :Options:
%
%   **'eval':**
%       A function handle to use for computing summary statistics of each
%       tissue class. Must accept exactly 1 argument and handle 'nan'
%       gracefully. e.g. '@(x1)(nanvar(x1))'
%
% :Outputs:
%
%   **values:**
%        mean gray matter, white, CSF. If 'eval' option is passed,
%        specified function will be used in place of mean.
%
%   **components:**
%        first 5 components from each tissue class, observation x 5
%
%   **full_data_objects:**
%        Masked data objects for {gray white CSF}
%
%   **l2norms:**
%        Length-adjusted L2 norms (divided by sqrt(nvox))
%
%   **masks:**
%        optional input for using different masks; it should provide gray,
%        white and ventricle images in order.
% ..
%    Tor Wager, July 21, 2015
% ..

% Programmers' notes:
% Jan 2017:  Issue with vector lengths if obj has removed images, fixed (Tor)
%
% Feb 2017: change to mask. Gray matter is now sparse, old gray_matter_mask.img thresholded at .5.
%
% Sept 2017: added option for custom function evaluation in place of mean.

fxn = @(x1)(nanmean(x1,1));
masks = {'gray_matter_mask_sparse.img' 'canonical_white_matter.img' 'canonical_ventricles.img'};

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'eval'
                fxn = varargin{i+1};
            case 'masks'
                masks = varargin{i+1};
        end
    end
end

numcomps = 5;

obj = remove_empty(obj);  % return only non-empty values
nimgs = size(obj.dat, 2); % - sum(obj.removed_images);

values = NaN * zeros(nimgs, length(masks));
components = cell(1, length(masks));
full_data_objects = [];

for i = 1:length(masks)

    maskname = which(masks{i});

    if isempty(maskname)
       fprintf('Image %s cannot be found on path.\n', masks{i});
       error('Exiting');
    else
        fprintf('Extracting from %s.\n', masks{i});
    end

    masked_obj = apply_mask(obj, maskname);

    % get all values, if requested
    if nargout > 2

        full_data_objects{i} = remove_empty(masked_obj);

    end

    % get means
    masked_obj.dat(masked_obj.dat == 0) = NaN;

    myvalues = fxn(masked_obj.dat);
    myvalues = myvalues(:);

    % may need to insert omitted - no, return in reduced space
%     if length(masked_obj.removed_images) > 1
%         myvalues = naninsert(masked_obj.removed_images, myvalues);
%     end

    values(:, i) = myvalues;

    % get components
    if nargout > 1

        % NaNs will mess this up - remove voxel-wise
        [wasnan, dataforpca] = nanremove(masked_obj.dat);
        if any(wasnan), fprintf('Removing %3.0f voxels with one or more NaNs\n', sum(wasnan)); end

        [~, components{i}] = pca(dataforpca', 'Economy', true, 'NumComponents', numcomps);

        % may need to insert omitted
%         if length(masked_obj.removed_images) > 1
%             components{i} = naninsert(masked_obj.removed_images, components{i});
%         end

    end

    if nargout > 3

        l2norms(:, i) = getnorms(masked_obj);

    end

end % end masks

end % function

function n = getnorms(masked_obj)

% Vector L2 norm / sqrt(length) of vector

% divide by this value to normalize image

normfun = @(x) sum(x .^ 2) .^ .5;

%masked_obj = remove_empty(masked_obj);

x = masked_obj.dat;

%nv = size(x, 1);

for i = 1:size(x, 2)

    % remove nans, 0s
    xx = x(:, i);
    xx(xx == 0) = [];
    xx(isnan(xx)) = [];

    % divide by sqrt(length) so number of elements will not change scaling

    n(i) = normfun(xx) ./ sqrt(length(xx));

end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function r8mat_normal_01_test ( )

%*****************************************************************************80
%
%% R8MAT_NORMAL_01_TEST tests R8MAT_NORMAL_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_NORMAL_01_TEST\n' );
  fprintf ( 1, '  R8MAT_NORMAL_01 computes a matrix of Normal 01 values.\n' );

  m = 5;
  n = 4;
  seed = 123456789;
  [ r, seed ] = r8mat_normal_01 ( m, n, seed );

  r8mat_print ( m, n, r, '  The matrix:' );

  return
end


text:
function C = interleave_rows(A,B)
  % Takes two nxm matrices and creates a 2nxm matrix such that even rows are
  % from the first matrix and odd rows are from the second.
  %
  % C = interleave_rows(A,B)
  %
  % Input:
  %   A  First nxm matrix, top row of output will be same as top of this matrix
  %      along with every other row
  %   B  Second nxm matrix, bottom row of output will be same as bottom of this
  %      matrix along with every other row
  %
  % Output:
  %   C  2nxm matrix, alternating rows
  %
  %
  if(size(A)~=size(B))
    error('Input matrices must be same shape');
  end
  C = reshape([A(:) B(:)]',size(A,1)*2, size(A,2));
end

author: alecjacobson

repo: gptoolbox


text:
function [doseBinsV, volsHistV] = loadDVHMatrix(DVHNum, planC)
%"loadDVHMatrix"
%   Return the doseBinsV and volsHistV with the binWidth shift included.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.
%
%Usage:
%   function [doseBinsV, volsHistV] = loadDVHMatrix(DVHNum, planC)

indexS = planC{end};

%Calculate width of bins, use to find middle of each bin.
binWidthsV = diff(planC{indexS.DVH}(DVHNum).DVHMatrix(:,1));
lastBinWidth = binWidthsV(end);
binWidthsV(end+1,1) = lastBinWidth;

%Extract dose bin values, adding half to binwidth to get middle.
doseBinsV = planC{indexS.DVH}(DVHNum).DVHMatrix(:,1) + binWidthsV/2;
volsHistV = planC{indexS.DVH}(DVHNum).DVHMatrix(:,2);
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function MASKSTATS = canlab_maskstats(msks,imgs,varargin)
% Produces comparison of pattern mask and images
% e.g., look at NPS pattern expression in set of beta images
%
% :Usage:
% ::
%
%    MASKSTATS = canlab_maskstats(maskfiles,imgfiles,[options])

% :Inputs:
%
%   **maskfiles:**
%      string or cellstring of mask filenames or fmri_data object

%   **imgfiles:**
%      string or cellstring of image filenames or fmri_data object
%
% :Optional Inputs:
%
%   **ts:**
%        timeseries treatment: each string in imgfiles is assumed to be
%        a 4D file. Data will be returned with one column per time series and
%        one volume per row. If not all timeseries are same length, all will
%        be NaN-padded to the length of the longest timeseries.
%
%        Note: does not work with imgfiles input as fmri_data object
%
%   **dir2cell:**
%        will sort stats into cells based on directory containing the imgfile
%        they belong to such that each cell contains one directory's worth of
%        stats which is a vector with a value for each imgfile.
%        :Examples:
%        ::
%
%             # Input: a list of single trial betas for a set of subjects
%             b = filenames('sub*/*heat_trials*.img');
%             ms = canlab_maskstats('nps',b,'dot_product','dir2cell');
%
%        Output: includes the set of cells that are input to a mediation
%        analysis.
%
%   **keepzeros:**
%        don't remove zeros from imgfiles before taking measurements
%
%   **keepzerosmask:**
%        don't remove zeros from maskfiles before taking measurements
%
%   **single:**
%        leave data as single (DEFAULT: convert to double)
%
%   **trinarize:**
%        trinarize maskfile (set values larger than 0 to 1 and values less than
%        zero to -1)
%
%   **noreshape:**
%        don't attempt to reshape results according to imgfiles array
%
%   **nobin:**
%        don't binarize mask before extracting mean or std
%
% :Note: ts, dir2cell, and noreshape are mutually exclusive options
%
% **Built-In Masks**
% The following strings can be given as the maskfile argument to
% call up built-in mask files:
%
%   **nps:**
%        weights_NSF_grouppred_cvpcr.img
%
%   **nps_thresh:**
%        weights_NSF_grouppred_cvpcr_FDR05.img
%
%   **nps_thresh_smooth:**
%        weights_NSF_grouppred_cvpcr_FDR05_smoothed_fwhm05.img
%
% :Measure Options:
%
%   **all:**
%        add: mean, dot_product, centered_dot_product,
%             cosine_similarity, and correlation
%
%   **mean:** (DEFAULT)
%      apply binarized mask to images and return means
%      mean(img .* abs(bin(mask)))
%
%   **std:**
%        apply binarized mask to images and return standard deviations
%        std(img .* abs(bin(mask)))
%
%   **dot_product:**
%        dot(mask, img)
%
%   **cosine_similarity:**
%        dot(mask, img) / (norm(mask) * norm(img))
%
%   **correlation:**
%        corr(mask, img)
%
%   **centered_dot_product:**
%        dot(mask-mean(mask), img-mean(img))
%
%
% :Details:
%   - imgfiles are spatially resampled to maskfiles
%   - voxels with zeros in maskfile are removed
%   - in-mask voxels with zeros in imgfiles will generate warnings
%
% ..
%    AUTHOR: Luka Ruzic 2013
% ..

OP = {}; % struct to contain desired "operations" (e.g., dot product, mean, etc)
%% parse arguments
ALLOPS = false; % adds all implemented operations
DO_TRINARIZE = false; % signed binarize option
DO_ZERO2NAN_DATA = true; % treat zeros as non-data in data input
DO_ZERO2NAN_MASK = true; % treat zeros as non-data in mask input
DO_RESHAPE = true; % attempt to reshape output data according to shape of input data
DO_BIN = true; % binarize mask
DO_TS = false; % timeseries mode
DO_CELLS = false; % cell array output mode

i=1;
while i<=numel(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case {'binmean' 'mean' 'norm' ...
                  'dot_product' 'centered_dot_product' 'cosine_similarity' 'correlation'}
                OP{end+1} = varargin{i}; %#ok
            case 'all'
                ALLOPS = true;
            case 'keepzeros'
                DO_ZERO2NAN_DATA = false;
            case 'keepzerosmask'
                DO_ZERO2NAN_MASK = false;
            case 'trinarize'
                DO_TRINARIZE = true;
            case 'nobin'
                DO_BIN = false;
            case 'noreshape'
                DO_RESHAPE = false;
                DO_CELLS = false;
                DO_TS = false;
            case 'ts'
                DO_RESHAPE = true;
                DO_TS = true;
                DO_CELLS = false;
            case 'dir2cell'
                DO_RESHAPE = true;
                DO_CELLS = true;
                DO_TS = false;
            otherwise
                error('Unrecognized argument %s',varargin{i})
        end
    elseif iscellstr(varargin{i})
        for j=1:numel(varargin{i})
            switch varargin{i}{j}
                case {'binmean' 'mean' 'std' 'nonbinmean' 'nonbinstd' 'norm' ...
                        'dot_product' 'centered_dot_product' 'cosine_similarity' 'correlation'}
                    OP{end+1} = varargin{i}{j}; %#ok
                case 'all'
                    ALLOPS = true;
                otherwise
                    error('Unrecognized argument %s',varargin{i}{j})
            end
        end
    else
        disp(varargin{i})
        error('Above argument unrecognized')
    end
    i=i+1;
end

if ALLOPS, OP = [OP 'mean' 'std' 'dot_product' 'centered_dot_product' 'cosine_similarity' 'correlation']; end

if isempty(OP), OP = {'mean'}; end

%% error checking
if isempty(msks), error('Must provide one or more maskfiles'); end
if ischar(msks), msks = cellstr(msks); end
if ~iscellstr(msks) && ~isa(msks,'fmri_data'), error('maskfiles must be a string, cell array of strings, or fmri_data object');end

if isempty(imgs), error('Must provide one or more imgfiles'); end
if ischar(imgs), imgs = cellstr(imgs); end
if ~iscellstr(imgs) && ~isa(imgs,'fmri_data'), error('imgfiles must be string, cell array of strings, or fmri_data object'); end

MASKSTATS(numel(msks)) = struct; % initializing

%% prepare mask files
if ~isa(msks,'fmri_data')
    for m = 1:numel(msks)
        switch msks{m}
            case 'nps'
                msks{m} = which('weights_NSF_grouppred_cvpcr.img');
            case 'nps_thresh'
                msks{m} = which('weights_NSF_grouppred_cvpcr_FDR05.img');
            case 'nps_thresh_smooth'
                msks{m} = which('weights_NSF_grouppred_cvpcr_FDR05_smoothed_fwhm05.img');
            otherwise
                if ~exist(msks{m},'file')
                    error('No such file: %s',msks{m})
                end
        end
    end
end

%% prepare data
% fprintf('LOADING IMAGE FILES\n')
if isa(imgs,'fmri_data')
    imgdata = imgs;
    imgdata = imgdata.replace_empty;
    if DO_TS
        warning('''ts'' mode is incompatible with fmri_data input: turning off'); %#ok
        DO_TS = false;
    end
    imgfiles = cellstr(imgdata.fullpath);
else
    imgfiles = imgs;
    evalc('imgdata = fmri_data(imgfiles);');
end

imgdata.dat = double(imgdata.dat); % convert singles to doubles
if DO_ZERO2NAN_DATA, imgdata.dat(imgdata.dat==0) = NaN; end  % remove zeros if desired

if DO_TS
    for i = 1:numel(imgfiles)
        ni = nifti(imgfiles{i});
        tsdim(i) = ni.dat.dim(4); %#ok
    end
end

if DO_CELLS
%     [ignore1 ignore2 x] = unique(regexprep(imgfiles,'/[^/]*$','')); %#ok
    [ignore1 ignore2 x] = unique(regexprep(cellstr(imgdata.fullpath),'/[^/]*$','')); %#ok
    reshaperows = histc(x,1:max(x));
end

%% produce comparison
if isa(msks,'fmri_data')
    allmaskdata = msks;
    allmaskdata = allmaskdata.replace_empty;
    nmasks = size(allmaskdata.dat,2);
else
    nmasks = numel(msks);
end
for m = 1:nmasks
    clear imgdat maskdat imgdatcent maskdatcent imgdatnorm maskdatnorm

    % save filenames
    % imgs (modify cell array according to mode)
    if DO_CELLS
        for i = 1:numel(reshaperows)
            st = sum(reshaperows(1:i-1))+1;
            MASKSTATS(m).imgfiles{i} = imgfiles(st:st+reshaperows(i)-1);
        end
    elseif DO_TS
        MASKSTATS(m).imgfiles = imgfiles';
    else
        MASKSTATS(m).imgfiles = imgfiles;
    end

    % load and prep mask data

    if isa(msks,'fmri_data')
        maskdata = allmaskdata;
        maskdata.dat(:,[1:size(allmaskdata.dat,2)]~=m) = [];
        MASKSTATS(m).maskfile = deblank(allmaskdata.fullpath(m,:));
    else
        evalc('maskdata = fmri_data(msks{m});');
        MASKSTATS(m).maskfile = msks{m};
    end
    fprintf('GETTING DATA FOR MASK: %s\n',MASKSTATS(m).maskfile);
    maskdata.dat = double(maskdata.dat); % convert singles to doubles
    %if DO_ZERO2NAN_MASK, maskdata.dat(maskdata.dat == 0) = NaN; end % SHOULD THIS HAPPEN HERE?   remove zeros
    maskdata = maskdata.resample_space(imgdata); % resample
    if DO_ZERO2NAN_MASK, maskdata.dat(maskdata.dat == 0) = NaN; end % remove zeros

    % select in-mask voxels
    wh_inmask = ~isnan(maskdata.dat);
    maskdat = maskdata.dat(wh_inmask);
    imgdat = imgdata.dat(wh_inmask,:);

    % trinarize
    if DO_TRINARIZE
        maskdat(maskdat<0) = -1;
        maskdat(maskdat>0) = 1;
    end

    % warn about in-mask zeros
    if DO_ZERO2NAN_DATA
        inmaskzeros = isnan(imgdat);
        if any(inmaskzeros)
            fprintf('WARNING! There are zeros in your data within this mask!\n')
            %fprintf('Consider carefully whether this impairs your ability to compare data points.\n')
            %fprintf('  These are being ignored as non-data but your results will be impacted!\n')
            for i = find(sum(inmaskzeros))
                fprintf(' %7d in-mask zeros in %s\n',sum(inmaskzeros(:,i)),deblank(imgdata.fullpath(i,:)));
            end
        end
    end

    % perform operation(s)
    for i = 1:numel(OP)
        switch OP{i}
            case 'mean'
                if DO_BIN
                    MASKSTATS(m).stats.mean = nanmean(imgdat)';
                else
                    MASKSTATS(m).stats.nonbinmean = nanmean(bsxfun(@times,imgdat,maskdat))';
                end

            case 'std'
                if DO_BIN
                    MASKSTATS(m).stats.std = nanstd(imgdat)';
                else
                    MASKSTATS(m).stats.nonbinstd = nanstd(bsxfun(@times,imgdat,maskdat))';
                end

            case 'norm'
                MASKSTATS(m).stats.norm = (nansum(imgdat .^ 2) .^ .5)';

            case 'dot_product'
                MASKSTATS(m).stats.dot_product = nansum(bsxfun(@times,imgdat,maskdat))';

            case 'centered_dot_product'
                imgdatcent = bsxfun(@minus,imgdat,nanmean(imgdat));
                maskdatcent = bsxfun(@minus,maskdat,nanmean(maskdat));
                MASKSTATS(m).stats.centered_dot_product = nansum(bsxfun(@times,imgdatcent,maskdatcent))';

            case 'cosine_similarity'
                imgdatnorm = nansum(imgdat .^ 2) .^ .5;
                maskdatnorm = nansum(maskdat .^ 2) .^ .5;
                MASKSTATS(m).stats.cosine_similarity = (nansum(bsxfun(@times,imgdat,maskdat)) ./ (imgdatnorm .* maskdatnorm))';

            case 'correlation'
                MASKSTATS(m).stats.correlation = corr(imgdat,maskdat,'rows','pairwise');

            otherwise
                error('Comparison type (%s) unrecognized',OP{i})
        end

        % reshape according to mode
        if DO_RESHAPE
            if DO_TS
                if numel(unique(tsdim)) == 1
                    tslen = tsdim(1);
                else
                    % not all timeseries are same length
                    % assume all start at same time and pad with NaNs
                    tslen = max(tsdim);
                    for j = 1:numel(tsdim)
                        cutind = ((j-1)*tslen) + tsdim(j);
                        MASKSTATS(m).stats.(OP{i}) = [...
                            MASKSTATS(m).stats.(OP{i})(1:cutind); ...
                            nan(tslen-tsdim(j),1); ...
                            MASKSTATS(m).stats.(OP{i})(cutind+1:end)];
                    end
                    % alternative coding option: vector -> cells -> NaN-pad the cells -> cell2mat %
                end
                MASKSTATS(m).stats.(OP{i}) = reshape(MASKSTATS(m).stats.(OP{i}),[tslen numel(imgfiles)]);
            elseif DO_CELLS
                MASKSTATS(m).stats.(OP{i}) = mat2cell(MASKSTATS(m).stats.(OP{i}),reshaperows,1)';
            else
                try MASKSTATS(m).stats.(OP{i}) = reshape(MASKSTATS(m).stats.(OP{i}),size(imgfiles)); end %#ok
            end
        end
    end
end

end


text:
% Fig. 5.23  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
% script to generate Fig. 5.23
n=1;
d=[1 1 0];
n1=[1 2 ];
d1=conv(d,[1 20]);
d2=conv(d,[1 10]);
 sys1=tf(n1,d);
rlocus(sys1)
hold on
sys2=tf(n1,d1);
rlocus(sys2,':')
sys3=tf(n1,d2);
rlocus(sys3,'--')
axis([-6 2 -3 3])
title('Fig. 5.23 Root locus with PD or lead compensation')
z=0:.1:.9;
 wn=1:1:6;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function y = idct2(x)
  t = idct(x);
  t = idct(t');
  y = t';

author: thomaskuestner


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bra.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function f = bra(x,y)
% Branin's function
function f = bra(x,y)
a=1;
b=5.1/(4*pi*pi);
c=5/pi;
d=6;
h=10;
ff=1/(8*pi);
if nargin == 1
  x1 = x(1);
  x2 = x(2);
else
  x1 = x;
  x2 = y;
end
f=a.*(x2-b.*x1.^2+c.*x1-d).^2+h.*(1-ff).*cos(x1)+h;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [cha, chp] = convexhullProps(grains)
% Area of the convex hull of grains
% Note: this is only reasonable for smooth and sufficiently large grains
%
%
% Syntax
%   [cha, chp] = convexhullArea(grains)
%
% Input
%  grains - @grain2d
%
% Output
%  cha - area of convex hull
%  chp - perimeter of covex hull
%
%

p = zeros(size(grains));

% store this in local variables for speed reasons
X = grains.V(:,1);
Y = grains.V(:,2);

poly = grains.poly;
% remove inclusions
incl = grains.inclusionId;
for i = find(incl>0).'
  poly{i} = poly{i}(1:end-incl(i));
end

% compute convex hull perimeters
for id = 1:length(grains)

  % extract coordinates
  xGrain = X(poly{id});
  yGrain = Y(poly{id});

  % compute convex hull
  ixy = convhull(xGrain,yGrain);

  % area
  cha(id) = polySgnArea(xGrain(ixy),yGrain(ixy));

  % perimeter
  chp(id) = sum(sqrt(...
    (xGrain(ixy(1:end-1)) - xGrain(ixy(2:end))).^2 + ...
    (yGrain(ixy(1:end-1)) - yGrain(ixy(2:end))).^2));

end

author: mtex-toolbox

repo: mtex

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


text:
function b=BlockWise(img,level,flag)

R_mean = mean(mean(img(:,:,1)));
G_mean = mean(mean(img(:,:,2)));
B_mean = mean(mean(img(:,:,3)));
W = size(img,1);
H = size(img,2);

switch flag
    case 5
        block = repmat([R_mean G_mean B_mean]',[1,32,32]);
        block = permute(block,[3,2,1]);
    case 4
        block = repmat([R_mean G_mean B_mean]',[1,32,32]);
        block = permute(block,[3,2,1]);
    case 3
        block = repmat([R_mean+30 G_mean B_mean]',[1,32,32]);
        block = permute(block,[3,2,1]);
    case 2
        block = repmat([R_mean+50 G_mean B_mean]',[1,32,32]);
        block = permute(block,[3,2,1]);
    case 1
        block = repmat([0 G_mean B_mean]',[1,32,32]);
        block = permute(block,[3,2,1]);
end

block = uint8(block);

for i=1:level
    r_W = randi([1 W-32],1);
    r_H = randi([2 H-32],1);
    img(r_W:r_W+31,r_H:r_H+31,:) = block;
end

b= img;

end

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


text:
name_list = {'bunny' 'david50kf' 'hand' 'lion-head' 'david_head' 'aphro' 'gargoyle' 'screwdriver' 'elephant-50kv'};

for itname = 1:length(name_list)
    name = name_list{itname};
    test_levelset_mesh;
end

author: gpeyre

repo: matlab-toolboxes


text:
function F = ismember(x,F)
% Internal class for constraint list

F = lmi(F);
F = ismember(x,F);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function figbolden(pos,pv,m)
%FIGBOLDEN embolden, resize and recolour the current figure =(POS,PV,M)
%
% Inputs: pos = [xmin ymin width height] gives the lower left corner position and the window size in pixels
%               [width height] leaves the lower left corner alone
%               [width] has a standard aspect ratio of 4:3
%               [-width/height] leaves the area unchanged but fixes the aspect ratio
%         pv is a cell name containing attribute-value pairs.
%            default = {'FontName' 'Arial'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8}
%            Note that the "listfonts" MATLAB command will list the available fonts
%         m is a mode string:
%                'l' lists the changes made
%                'd' append default pv settings (e.g. use [] for second argument)
%                'c' change default colours to improve contrast on a white background
%                    g->[0,0.7,0],c->[0,0.7,0.7],y->[0.83,0.83,0]
%                'x' suppresses all changes
%
% Bug: gives an error message if log axes have been used

%      Copyright (C) Mike Brookes 2003
%      Version: $Id: figbolden.m 9186 2016-12-13 08:51:37Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ps={'Title' 'XLabel' 'YLabel' 'Children'};
if nargin<3
    m='';
end
if nargin<2 || any(m=='d')
    pv={'FontName' 'Arial'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8};
end
pp={'Symbol';'Wingdings'};      % protected fonts
if nargin<1
    pos=[];
end
mlist=any(m=='l');  % list changes
mnotx=~any(m=='x'); % do changes
scsz=get(0,'screensize');
if length(pos)
    po=get(gcf,'position');
    if length(pos)>2            % position is specified
        po(1:2)=pos(1:2);
        pos(1:2)=[];      % remove xmin,ymin
    end
    if length(pos)>1
        po(3:4)=pos(1:2);
    else
        if pos(1)>0
            po(3:4)=[1 0.75]*pos(1);
        else
            po(3:4)=[-pos(1) 1]*sqrt(-po(3)*po(4)/pos(1)); % preserve area
        end
    end
    set(gcf,'position',po);
end
if any(m=='c')
    cc='gcy';
    cv=[0,0.7,0;0,0.7,0.7;0.83,0.83,0];
    for ic=1:3
        hlist=findobj(gcf,'color',cc(ic));
        if length(hlist)
            if mnotx
                set(hlist,'color',cv(ic,:));
            end
            if mlist
            for ih=1:length(hlist)
                 fprintf(['change %f Color: ''%c'' -> [%g %g %g]\n'],hlist(ih),cc(ic),cv(ic,:));
            end
            end
        end
    end
end
hlist=get(gcf,'children');
while length(hlist)
    pl=get(hlist(1));
    %fprintf('list length = %d, handle = %f\n',length(hlist),hlist(1));
    for i=1:size(pv,1)
        if isfield(pl,pv{i,1})
            if i>1 || all(~strcmpi(get(hlist(1),pv{i,1}),pp))
                pval=get(hlist(1),pv{i,1});
                if ~all(size(pval)==size(pv{i,2})) || ~all(pval(:) == pv{i,2}(:))
                    if mnotx
                        set(hlist(1),pv{i,1},pv{i,2})
                    end
                    if mlist
                        if ischar(pval)
                            fch='s';
                        else
                            fch='g';
                        end
                        fprintf(['change %f %s: %' fch ' -> %' fch '\n'],hlist(1),pv{i,1},pval,pv{i,2});
                    end
                end
            end
        end
    end
    for i=1:length(ps)
        if isfield(pl,ps{i})
            hlist=[hlist; get(hlist(1),ps{i})];
            %fprintf('add %f:%s\n',hlist(1),ps{i});
        end
    end
    hlist(1)=[];
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [params, names] = ggwhiteKernExtractParam(kern)

% GGWHITEKERNEXTRACTPARAM Extract parameters from the GG WHITE kernel structure.
% FORMAT
% DESC Extract parameters from the gaussian
%	gaussian white kernel structure into a vector of parameters for optimisation.
% RETURN param : vector of parameters extracted from the kernel. If the
%	   field 'transforms' is not empty in the kernel matrix, the
%	   parameters will be transformed before optimisation (for example
%	   positive only parameters could be logged before being returned).
% ARG kern : the kernel structure containing the parameters to be
%	   extracted.
%
% FORMAT
% DESC extract parameters and
%	their names from the gaussian gaussian white kernel structure.
% RETURN param : vector of parameters extracted from the kernel. If the
%	   field 'transforms' is not empty in the kernel matrix, the
%	   parameters will be transformed before optimisation (for example
%	   positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing parameter names.
% ARG kern : the kernel structure containing the parameters to be
%	   extracted.
%
% SEEALSO : ggwhiteKernParamInit, ggwhiteKernExpandParam, kernExtractParam, scg,
% conjgrad
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A Alvarez, 2009

% KERN

params = [kern.precisionG' kern.sigma2Noise kern.variance];

if nargout > 1
    if kern.isArd
        ynames = cell(1,kern.inputDimension);
        for i=1:kern.inputDimension,
            ynames{i}=['inverse width output ' num2str(i) '.'];
        end
        names = {ynames{:}, 'variance latent', 'sensitivity'};
    else
        names = {'inverse width output 1.' , 'variance latent', 'sensitivity'};
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [f] = ma_MinimizeMissionDuration(stateLog, eventID, celBodyData, maData)
%ma_maxTotalPropMass Summary of this function goes here
%   Detailed explanation goes here

    f = stateLog(end,1) - stateLog(1,1);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%go_calib_optim
%
%Main calibration function. Computes the intrinsic andextrinsic parameters.
%Runs as a script.
%
%INPUT: x_1,x_2,x_3,...: Feature locations on the images
%       X_1,X_2,X_3,...: Corresponding grid coordinates
%
%OUTPUT: fc: Camera focal length
%        cc: Principal point coordinates
%        alpha_c: Skew coefficient
%        kc: Distortion coefficients
%        KK: The camera matrix (containing fc and cc)
%        omc_1,omc_2,omc_3,...: 3D rotation vectors attached to the grid positions in space
%        Tc_1,Tc_2,Tc_3,...: 3D translation vectors attached to the grid positions in space
%        Rc_1,Rc_2,Rc_3,...: 3D rotation matrices corresponding to the omc vectors
%
%Method: Minimizes the pixel reprojection error in the least squares sense over the intrinsic
%        camera parameters, and the extrinsic parameters (3D locations of the grids in space)
%
%Note: If the intrinsic camera parameters (fc, cc, kc) do not exist before, they are initialized through
%      the function init_intrinsic_param.m. Otherwise, the variables in memory are used as initial guesses.
%
%Note: The row vector active_images consists of zeros and ones. To deactivate an image, set the
%      corresponding entry in the active_images vector to zero.
%
%VERY IMPORTANT: This function works for 2D and 3D calibration rigs, except for init_intrinsic_param.m
%that is so far implemented to work only with 2D rigs.
%In the future, a more general function will be there.
%For now, if using a 3D calibration rig, set quick_init to 1 for an easy initialization of the focal length

if ~exist('n_ima'),
   data_calib_no_read; % Load the images
   click_calib_no_read; % Extract the corners
end;

check_active_images;

check_extracted_images;

check_active_images;

desactivated_images = [];

recompute_extrinsic = (length(ind_active) < 100); % if there are too many images, do not spend time recomputing the extrinsic parameters twice..

%%% MAIN OPTIMIZATION CALL!!!!! (look into this function for the details of implementation)
go_calib_optim_iter;

if ~isempty(desactivated_images),

   param_list_save = param_list;

   fprintf(1,'\nNew optimization including the images that have been deactivated during the previous optimization.\n');
   active_images(desactivated_images) = ones(1,length(desactivated_images));
   desactivated_images = [];

   go_calib_optim_iter;

   if ~isempty(desactivated_images),
      fprintf(1,['List of images left desactivated: ' num2str(desactivated_images) '\n' ] );
   end;

   param_list = [param_list_save(:,1:end-1) param_list];

end;

%%%%%%%%%%%%%%%%%%%% GRAPHICAL OUTPUT %%%%%%%%%%%%%%%%%%%%%%%%

%graphout_calib;

repo: ekfmonoslam


text:
function test_suite = test_it
  initTestSuite;
end

function test_entropy
    data = [0 0 0 1 1 1 2 2 2 3 3 3];
    h = SPX_IT.entropy(data);
    assertElementsAlmostEqual(h, 2.0);
end

author: indigits

repo: sparse-plex


text:
function [U, S, V, details]  = lansvd(A, varargin)
    % We need to parse the options one by one
    params = inputParser;
    params.addRequired('A');
    % Number of singular values to compute
    params.addParameter('k', NaN);
    % Threshold for singular values
    params.addParameter('lambda', NaN);
    % Desired level of orthogonality
    params.addParameter('delta', -1);
    % Desired level of orthogonality after reorthogonalization
    params.addParameter('eta', -1);
    % Tolerance for iterated Gram-Schmidt procedure
    params.addParameter('gamma', -1);
    % Flag for classic/modified Gram-Schmidt
    params.addParameter('cgs', true);
    % Flag for extended local reorthogonalization
    params.addParameter('elr', false);
    % Verbosity level
    params.addParameter('verbosity', 0);
    % Maximum number of iterations
    params.addParameter('max_iters', -1);
    % Tolerance for convergence of singular values
    params.addParameter('tolerance', -1);
    % Initial vector specified by user
    params.addParameter('p0', []);
    % Parse the results
    params.parse(A, varargin{:});
    results = params.Results;
    options = struct;
    if(isobject(A))
        A = double(A);
    end
    if ~isnan(results.k)
        options.k = results.k;
    end
    if ~isnan(results.lambda)
        options.lambda = results.lambda;
    end
    if ~isfield(options, 'k') && ~isfield(options, 'lambda')
        % By default we compute 6 singular values
        options.k = 6;
    end
    options.delta  = results.delta;
    options.eta  = results.eta;
    options.gamma = results.gamma;
    options.cgs = results.cgs;
    options.elr = results.elr;
    options.verbosity = results.verbosity;
    options.max_iters = results.max_iters;
    options.tolerance = results.tolerance;
    if ~isempty(results.p0)
        options.p0 = results.p0;
    end
    if isstruct(A)
        M = A.M;
        N = A.N;
    elseif isnumeric(A)
        [M, N] = size(A);
    end
    details = struct;
    if min(M, N) < 1
        % A is empty
        U = []; S = []; V = []; ;
        if nargout == 1; U = zeros(0, 1); end;
        return;
    end
    if M*N == 1
        % This is the case of a singleton
        U = [1]; V = [1];
        a = A(1,1);
        S = [a];
        if nargout == 1; U = S; end;
        return;
    end
    if min(M, N) == 1
        % either a column vector or a row vector
        if isstruct(A)
            % function objects
            if N == 1
                % Single column
                A = A.A(1);
            else
                % Single row
                A = (A.At(1))';
            end
        end
        [U, S, V] = svd(full(A));
        if nargout == 1; U = S; end;
    end
    [U, S, V, alpha, beta, p, details] = mex_lansvd(A, options);
    % number of Lanczos vectors computed
    k_done = size(alpha, 1);
    p_norm = beta(k_done+1);
    if isfield(options, 'k')
        k = options.k;
    end
    if isfield(options, 'lambda')
        k = find(S <= options.lambda,1) - 1;
    end
    % Let's compute all the singular vectors if requested by caller
    if nargout>2 % computation of Ritz vectors
        % Form the k+1 x k bidiagonal matrix from alpha and beta
        B = spdiags([alpha(1:k_done) beta(2:k_done+1)], [0, -1], k_done+1, k_done);
        % Compute singular vectors
        [P,S,Q] = svd(full(B),0);
        % Singular values
        S = diag(S);
        % Keep the relevant K columns
        if size(Q,2)~=k
            Q = Q(:,1:k);
            P = P(:,1:k);
        end
        % Compute and normalize Ritz vectors (overwrites U and V to save memory).
        if p_norm~=0
            U = U*P(1:k_done,:) + (p/p_norm)*P(k_done+1,:);
        else
            U = U*P(1:k_done,:);
        end
        V = V*Q;
        % Make sure that the requested Ritz vectors are normalized
        for i=1:k
            nq = norm(V(:,i));
            if isfinite(nq) & nq~=0 & nq~=1
              V(:,i) = V(:,i)/nq;
            end
            nq = norm(U(:,i));
            if isfinite(nq) & nq~=0 & nq~=1
              U(:,i) = U(:,i)/nq;
            end
        end
    end % computation of Ritz vectors
    % Pick out desired part the spectrum
    S = S(1:k);
    if nargout == 1
        U = S;
    end
    if nargout >= 4
        details.alpha = alpha;
        details.beta = beta;
        details.p = p;
        details.k_done = k_done;
    end
end


text:
% POP_BIOSIG - import data files into EEGLAB using BIOSIG toolbox
%
% Usage:
%   >> OUTEEG = pop_biosig; % pop up window
%   >> OUTEEG = pop_biosig( filename, channels, type);
%
% Inputs:
%   filename - [string] file name
%
% Optional inputs:
%   'channels'   - [integer array] list of channel indices
%   'blockrange' - [min max] integer range of data blocks to import, in seconds.
%                  Entering [0 3] will import the first three blocks of data.
%                  Default is empty -> import all data blocks.
%   'ref'        - [integer] channel index or index(s) for the reference.
%                  Reference channels are not removed from the data,
%                  allowing easy re-referencing. If more than one
%                  channel, data are referenced to the average of the
%                  indexed channels. WARNING! Biosemi Active II data
%                  are recorded reference-free, but LOSE 40 dB of SNR
%                  if no reference is used!. If you do not know which
%                  channel to use, pick one and then re-reference after
%                  the channel locations are read in. {default: none}
%   'rmeventchan' - ['on'|'off'] remove event channel after event
%                  extraction. Default is 'on'.
%
% Outputs:
%   OUTEEG   - EEGLAB data structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Oct. 29, 2003-
%
% Note: BIOSIG toolbox must be installed. Download BIOSIG at
%       http://biosig.sourceforge.net
%       Contact a.schloegl@ieee.org for troubleshooting using BIOSIG.

% Copyright (C) 2003 Arnaud Delorme, SCCN, INC, UCSD, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG, command] = my_pop_biosig(filename, varargin);
EEG = [];
command = '';

if nargin < 1
	% ask user
	[filename, filepath] = uigetfile('*.*', 'Choose an BDF file -- pop_biosig()'); %%% this is incorrect in original version!!!!!!!!!!!!!!
    drawnow;

	if filename == 0 return; end
	filename = [filepath filename];

    % open file to get infos
    % ----------------------
    disp('Reading data file header...');
    dat = sopen(filename);

    % special BIOSEMI
    % ---------------
    if strcmpi(dat.TYPE, 'BDF')
        disp('We highly recommend that you choose a reference channel IF these are Biosemi data');
        disp('(e.g., a mastoid or other channel). Otherwise the data will lose 40 dB of SNR!');
    end
    uilist = { { 'style' 'text' 'String' 'Channel list (default all):' } ...
                 { 'style' 'edit' 'string' '' } ...
                 { 'style' 'text' 'String' [ 'Data range (in seconds) to read (default all [0 ' int2str(dat.NRec) '])' ] } ...
                 { 'style' 'edit' 'string' '' } ...
                 { 'style' 'text' 'String' 'Extract event - cannot be unset (set=yes)' } ...
                 { 'style' 'checkbox' 'string' '' 'value' 1 'enable' 'on' } {} ...
                 { 'style' 'text' 'String' 'Import continuous data (set=yes)' 'value' 1} ...
                 { 'style' 'checkbox' 'string' '' 'value' 0 } {} ...
                 { 'style' 'text' 'String' 'Reference chan(s) indices - required for BIOSEMI' } ...
                 { 'style' 'edit' 'string' '' } };
    geom = { [3 1] [3 1] [3 0.35 0.5] [3 0.35 0.5] [3 1] };
    result = inputgui( geom, uilist, 'pophelp(''pop_biosig'')', ...
                                 'Load data using BIOSIG -- pop_biosig()');
    if length(result) == 0 return; end

    % decode GUI params
    % -----------------
    options = {};
    if ~isempty(result{1}), options = { options{:} 'channels'   eval( [ '[' result{1} ']' ] ) }; end
    if ~isempty(result{2}), options = { options{:} 'blockrange' eval( [ '[' result{2} ']' ] ) }; end
    if length(result) > 2
        if ~isempty(result{4}), options = { options{:} 'ref'        eval( [ '[' result{4} ']' ] ) }; end
        if ~result{3},          options = { options{:} 'rmeventchan' 'off' }; end
    end
else
    options = varargin;
end

% decode input parameters
% -----------------------
g = finputcheck( options, { 'blockrange'  'integer' [0 Inf]    [];
                            'channels'    'integer' [0 Inf]    [];
                            'ref'         'integer' [0 Inf]    [];
                            'rmeventchan' 'string'  { 'on';'off' } 'on' }, 'pop_biosig');
if ischar(g), error(g); end

% import data
% -----------
EEG = eeg_emptyset;
if ~isempty(g.channels)
     dat = sopen(filename, 'r', g.channels,'OVERFLOWDETECTION:OFF');
else dat = sopen(filename, 'r', 0,'OVERFLOWDETECTION:OFF');
end
fprintf('Reading data in %s format...\n', dat.TYPE);

if ~isempty(g.blockrange)
    newblockrange    = g.blockrange;
    newblockrange(2) = min(newblockrange(2), dat.NRec);
    newblockrange    = newblockrange*dat.Dur;
    DAT=sread(dat, newblockrange(2)-newblockrange(1), newblockrange(1))';
else
    DAT=sread(dat, Inf)';% this isn't transposed in original!!!!!!!!
    newblockrange    = [];
end
dat = sclose(dat);

% convert to seconds for sread
% ----------------------------
EEG.nbchan          = size(DAT,1);
EEG.srate           = dat.SampleRate(1);
EEG.data            = DAT;
clear DAT;
% $$$ try  % why would you do the following???????  JO
% $$$     EEG.data            = EEG.data';
% $$$ catch,
% $$$     pack;
% $$$     EEG.data            = EEG.data';
% $$$ end;
EEG.setname 		= sprintf('%s file', dat.TYPE);
EEG.comments        = [ 'Original file: ' filename ];
EEG.xmin            = 0;
if strcmpi(dat.TYPE, 'BDF') || strcmpi(dat.TYPE, 'EDF')
    EEG.trials   = 1;
    EEG.pnts     = size(EEG.data,2);
else
    EEG.trials   = dat.NRec;
    EEG.pnts     = size(EEG.data,2)/dat.NRec;
end
if isfield(dat, 'Label') && ~isempty(dat.Label)
    EEG.chanlocs = struct('labels', cellstr(char(dat.Label)));
end
EEG = eeg_checkset(EEG);

% extract events % this part I totally revamped to work...  JO
% --------------
disp('Extracting events from last EEG channel...');
EEG.event = [];

% $$$ startval = mode(EEG.data(end,:)); % my code
% $$$ for p = 2:size(EEG.data,2)-1
% $$$     [codeout] = code(EEG.data(end,p));
% $$$     if EEG.data(end,p) > EEG.data(end,p-1) & EEG.data(end,p) >= EEG.data(end,p+1)
% $$$         EEG.event(end+1).latency =  p;
% $$$         EEG.event(end).type = bitand(double(EEG.data(end,p)-startval),255);
% $$$     end
% $$$ end

% lastout = mod(EEG.data(end,1),256);newevs = []; % andrey's code 8 bits
% codeout = mod(EEG.data(end,2),256);
% for p = 2:size(EEG.data,2)-1
%     nextcode = mod(EEG.data(end,p+1),256);
%     if codeout > lastout & codeout >= nextcode
%         newevs = [newevs codeout];
%         EEG.event(end+1).latency =  p;
%         EEG.event(end).type = codeout;
%     end
%     lastout = codeout;
%     codeout = nextcode;
% end

%lastout = mod(EEG.data(end,1),256*256);newevs = []; % andrey's code 16 bits
%codeout = mod(EEG.data(end,2),256*256);
%for p = 2:size(EEG.data,2)-1
%    nextcode = mod(EEG.data(end,p+1),256*256);
%    if (codeout > lastout) & (codeout >= nextcode)
%        newevs = [newevs codeout];
%        EEG.event(end+1).latency =  p;
%        EEG.event(end).type = codeout;
%    end
%    lastout = codeout;
%    codeout = nextcode;
%end

% Modifieded by Andrey (Aug.5,2008) to detect all non-zero codes:
thiscode = 0;lastcode=0;
for p = 1:size(EEG.data,2)-1
    prevcode = thiscode;
    thiscode = mod(EEG.data(end,p),256*256);   % andrey's code - 16 bits
    if (thiscode ~= 0) && (thiscode~=prevcode) && (thiscode~=lastcode) % fix to avoid repeated codes (per Ying's demand)
        EEG.event(end+1).latency =  p;
        EEG.event(end).type = thiscode;
        lastcode = thiscode;
    end
end

if strcmpi(g.rmeventchan, 'on')
    EEG.data(dat.BDF.Status.Channel,:) = [];
    EEG.nbchan = size(EEG.data,1);
    if ~isempty(EEG.chanlocs)
        EEG.chanlocs(dat.BDF.Status.Channel,:) = [];
    end
end
EEG = eeg_checkset(EEG, 'eventconsistency');

% $$$ if ~isempty(dat.EVENT)
% $$$     if isfield(dat, 'out') % Alois fix for event interval does not work
% $$$         if isfield(dat.out, 'EVENT')
% $$$             dat.EVENT = dat.out.EVENT;
% $$$         end
% $$$     end
% $$$     if ~isempty(newblockrange)
% $$$         interval(1) = newblockrange(1) * dat.SampleRate(1) + 1;
% $$$         interval(2) = newblockrange(2) * dat.SampleRate(1);
% $$$     else interval = [];
% $$$     end
% $$$     EEG.event = biosig2eeglabevent(dat.EVENT, interval); % Toby's fix
% $$$     if strcmpi(g.rmeventchan, 'on') & strcmpi(dat.TYPE, 'BDF') & isfield(dat, 'BDF')
% $$$         disp('Removing event channel...');
% $$$         EEG.data(dat.BDF.Status.Channel,:) = [];
% $$$         EEG.nbchan = size(EEG.data,1);
% $$$         if ~isempty(EEG.chanlocs)
% $$$             EEG.chanlocs(dat.BDF.Status.Channel,:) = [];
% $$$         end
% $$$     end
% $$$     EEG = eeg_checkset(EEG, 'eventconsistency');
% $$$ else
% $$$     disp('Warning: no event found. Events might be embedded in a data channel.');
% $$$     disp('         To extract events, use menu File > Import Event Info > From data channel');
% $$$ end

% rerefencing
% -----------
if ~isempty(g.ref)
    disp('Re-referencing...');
    EEG.data = EEG.data - repmat(mean(EEG.data(g.ref,:),1), [size(EEG.data,1) 1]);
    if length(g.ref) == size(EEG.data,1)
        EEG.ref  = 'average';
    end
    if length(g.ref) == 1
        disp([ 'Warning: channel ' int2str(g.ref) ' is now zeroed (but still present in the data)' ]);
    else
        disp([ 'Warning: data matrix rank has decreased through re-referencing' ]);
    end
end

% convert data to single if necessary
% -----------------------------------
EEG = eeg_checkset(EEG,'makeur');   % Make EEG.urevent field
EEG = eeg_checkset(EEG);

% history
% -------
if isempty(options)
    command = sprintf('EEG = my_pop_biosig(''%s'');', filename);
else
    command = sprintf('EEG = my_pop_biosig(''%s'', %s);', filename, vararg2str(options));
end;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%
%
%
%  Test the data density plot
%
%
x = randn(2048, 1);
y = randn(2048, 1);
x(1:512) = x(1:512) + 2.75;
x(1537:2048) = x(1537:2048) + 2.75;
y(1025:2048) = y(1025:2048) + 2.75;

% On scatter plot you probably can't see the data density
scatter(x, y);
% On data density plot the structure should be visible
DataDensityPlot(x, y, 32);


text:
function r8mat_ref_test ( )

%*****************************************************************************80
%
%% R8MAT_REF_TEST tests R8MAT_REF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2008
%
%  Author:
%
%    John Burkardt
%
  m = 4;
  n = 7;

  a = [ ...
    1.0, -2.0, 3.0, -1.0; ...
    3.0, -6.0, 9.0, -3.0; ...
    0.0,  0.0, 0.0,  0.0; ...
    2.0, -2.0, 0.0,  1.0; ...
    6.0, -8.0, 6.0,  0.0; ...
    3.0,  3.0, 6.0,  9.0; ...
    1.0,  1.0, 2.0,  3.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_REF_TEST\n' );
  fprintf ( 1, '  R8MAT_REF computes the row echelon form of a matrix.\n' );

  r8mat_print ( m, n, a, '  Input A:' );

  a = r8mat_ref ( m, n, a );

  r8mat_print ( m, n, a, '  REF form:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dt = EulerDT2D(Q, gamma)

% function dt = EulerDT2D(Q, gamma)
% purpose: compute the time step dt for the compressible Euler equations

Globals2D;

rho = Q(:,:,1); rhou = Q(:,:,2); rhov = Q(:,:,3); Ener = Q(:,:,4);
rho = rho(vmapM); rhou = rhou(vmapM); rhov = rhov(vmapM); Ener = Ener(vmapM);

u = rhou./rho; v = rhov./rho;
p = (gamma-1.0)*(Ener - rho.*(u.^2+v.^2)/2); c = sqrt(abs(gamma*p./rho));

dt = 1/max( ((N+1)^2)*.5*Fscale(:).*(sqrt ( u(:).^2 + v(:).^2 ) + c(:)));

rhoprange = [min(min(rho)), max(max(rho)), min(min(p)), max(max(p))]
return

repo: Aero-matlab


text:
function vl_testnn(varargin)
%VL_TESTNN Run MatConvNet test suite
% VL_TESTNN('option', value, ...) takes the following options:
%  `cpu`:: true
%    Run the CPU tests.
%
%  `gpu`:: false
%    Run the GPU tests.
%
%  `single`:: true
%    Perform tests in single precision.
%
%  `double`:: false
%    Perform tests in double precision.
%
%  `command`:: `'nn'`
%    Run only tests which name starts with the specified substring.
%    E.g. `vl_testnn('command', 'nnloss') would run only the nnloss tests.
%
%  `break`:: false
%    Stop tests in case of error.
%
%  `tapFile`:: ''
%    Output the test results to a file. If a specified file does
%    exist it is overwritten.
%
%  This function uses the Matlab unit testing framework which was
%  introduced in Matlab R2013a (v8.1).

% Copyright (C) 2015-16 Andrea Vedaldi, Karel Lenc.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.cpu = true ;
opts.gpu = false ;
opts.single = true ;
opts.double = false ;
opts.command = 'nn' ;
opts.break = false ;
opts.tapFile = '';
opts = vl_argparse(opts, varargin) ;

import matlab.unittest.constraints.* ;
import matlab.unittest.selectors.* ;
import matlab.unittest.plugins.TAPPlugin;
import matlab.unittest.plugins.ToFile;

% Choose which tests to run
sel = HasName(StartsWithSubstring(opts.command)) ;
if ~opts.gpu
  sel = sel & ~HasName(ContainsSubstring('device=gpu')) ;
end
if ~opts.cpu
  sel = sel & ~HasName(ContainsSubstring('device=cpu')) ;
end
if ~opts.double
  sel = sel & ~HasName(ContainsSubstring('dataType=double')) ;
end
if ~opts.single
  sel = sel & ~HasName(ContainsSubstring('dataType=single')) ;
end

% Run tests
root = fileparts(mfilename('fullpath')) ;
suite = matlab.unittest.TestSuite.fromFolder(fullfile(root, 'suite'), sel) ;
runner = matlab.unittest.TestRunner.withTextOutput('Verbosity',3);
if opts.break
  runner.addPlugin(matlab.unittest.plugins.StopOnFailuresPlugin) ;
end
if ~isempty(opts.tapFile)
  if exist(opts.tapFile, 'file')
    delete(opts.tapFile);
  end
  runner.addPlugin(TAPPlugin.producingOriginalFormat(ToFile(opts.tapFile)));
end
result = runner.run(suite);
display(result)

author: guosheng

repo: refinenet


text:
classdef StructuredMesh < handle

    properties (Access = public)
        nx
        ny
        x
        y
        mesh
    end

    properties (Access = private)
       xv
       yv
    end

    methods (Access = public)

        function obj = StructuredMesh(cParams)
            obj.init(cParams);
            s.coord  = obj.createCoordinates();
            s.connec = obj.createConnectivities();
            obj.createMesh(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            [obj.x,obj.y] = meshgrid(cParams.x,cParams.y);
            obj.nx = length(obj.x);
            obj.ny = length(obj.y);
        end

        function coord = createCoordinates(obj)
            coord(:,1) = reshape(obj.x',1,[]);
            coord(:,2) = reshape(obj.y',1,[]);
        end

        function connec = createConnectivities(obj)
            Nx = obj.nx;
            Ny = obj.ny;
            sw = obj.createVertex(1:Nx-1,1:Ny-1);
            se = obj.createVertex(2:Nx,1:Ny-1);
            ne = obj.createVertex(2:Nx,2:Ny);
            nw = obj.createVertex(1:Nx-1,2:Ny);
            connec = [sw se ne nw];
        end

        function v = createVertex(obj,xv,yv)
            Ny = obj.ny;
            v = bsxfun(@(x,y) x + Ny*(y-1),xv',yv);
            v = v(:);
        end

        function createMesh(obj,s)
            s.kFace = 0;
            obj.mesh = Mesh(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function res = vg_method_train(X, Y, args)

% Variational Garrote: performs linear regression with L0-norm penalty (Spike and Slab model)
%   See file test_vg.m for an example of use
%
% required parameters (n input dimension, p samples)
%   X     : n x p (training set, input)
%   Y     : 1 x s (training set, output)
%
% args: optional parameters (default)
%   method      : method for optimization 'dual' or 'regression' for fixed gamma ('dual')
%   maxiter     : maximum number of iterations for optimization for fixed gamma (1e4)
%   max_sum_m   : increases gamma values until sum(m)=max_sum_m  (n/2)
%   beta_max    : increases gamma values until beta=beta_max (1e3)
%   n_gamma     : number of gamma values to scan (50)
%   dmmin       : convergence threshold for mean field error (1e-12)
%   valset      : part of the training set used for validation (0.1*p)

%----------------
% REQUIRED PARAMS

n = size(X,1);

%----------------
% OPTIONAL PARAMS

% method for optimization {dual or regression} for fixed gamma
if ~isfield(args, 'method') method='dual'; else method=args.method; end

% maximum number of iterations for optimization for fixed gamma
if ~isfield(args, 'maxiter') maxiter=1e4; else maxiter=args.maxiter; end

% increases gamma values until sum(m)=max_sum_m
if ~isfield(args, 'max_sum_m') max_sum_m=n/2; else max_sum_m=args.max_sum_m; end

% increases gamma values until beta=beta_max
if ~isfield(args, 'beta_max') beta_max=1e3; else beta_max=args.beta_max; end

% number of gamma values to scan
if ~isfield(args, 'n_gamma') n_gamma=50; else n_gamma=args.n_gamma; end

% convergence threshold for mean field error
if ~isfield(args, 'dmmin') dmmin=1e-12; else dmmin=args.dmmin; end

% part of the training set used for validation (default 0.1*p)
if ~isfield(args, 'valset') valset=ceil(.25*n); else valset=args.valset; end

% randomly split training and validation datasets
dataok = false;
nits = 1;
p = size(X,2)-valset;
pv = valset;
while ~dataok && nits < 10
    it = randperm(size(X,2));
    xv = X(:,it(1:valset));
    yv = Y(:,it(1:valset));
    x = X(:,it(valset+1:end));
    y = Y(:,it(valset+1:end));

    % normalize training
    x=x-mean(x,2)*ones(1,p);
    dx=sqrt(1/p*sum(x.^2,2));
    x=x./(dx*ones(1,p));
    y=y-mean(y);

    % normalize validation
    xv=xv-mean(xv,2)*ones(1,pv);
    dxv=sqrt(1/pv*sum(xv.^2,2));
    xv=xv./(dxv*ones(1,pv));
    yv=yv-mean(yv);

    dataok = ~any(isnan(x(:))) && ~any(isnan(xv(:)));
    nits=nits+1;
end
if nits==10
    error('VG Error: Increase training set size');
end

%----------------
% compute garrote solution for range of gammas
% first from gamma_min to gamma_max and then in
% a second pass from gamma_max to gamma_min.

% C is input data covariance matrix.
if strcmp(method, 'regression')
    if n<=1500,
        C=x*x'/p;
    end;
end

% b is input output covariance
b=x*y'/p;

% sigma is output variance
sigmay=y*y'/p;

% set gamma range (min, max and step size)
delta=1e-8;
[b2sort,isort]=sort(b.^2,'descend');
bsort=b(isort);
gamma_min=log(delta*sigmay/p/max(abs(b)));
eps_gamma=0.001;
gamma_max=eps_gamma*gamma_min;
gamma_all =linspace(gamma_min,gamma_max,n_gamma);

% initial step size of mean field update
eta0=1; %e-2;
% initial step size for change in w in dual.m
eta_w0=0.02;

% input data variance
chi_ii=1/p*sum(x.^2,2);
if sum(abs(chi_ii-1)>1e-10),
    fprintf('input design matrix is not normalized\n');
    pause
end;

lg=n_gamma;
kl_all=inf(lg,2);
v_all=inf(lg,2,n);
m_all=inf(lg,2,n);
beta_all=inf(lg,2);
v_mf_all=inf(lg,n);
m_mf_all=inf(lg,n);
iter_all=inf(lg,2);
beta_mf_all=inf(1,lg);
error_mf_all=inf(1,lg);
errorv_mf_all=inf(1,lg);
m=zeros(1,n);

% the estimated inverse noise variance beta is initialized as the
% output variance
beta=1/sigmay;
i=0;

% for gamma is gamma_min to gamma_max, or when some criteria are
% satisfied
while (beta<beta_max)&&(i<n_gamma)&&(sum(m)<max_sum_m),
	i=i+1;
	gamma=gamma_all(i);
    eval(method);
	v_all(i,1,:)=v;
	m_all(i,1,:)=m;
	beta_all(i,1)=beta;
	iter_all(i,1)=iter;
	kl_all(i,1)=kl1;
	fprintf('gamma = %f beta = %f sum(m) = %f iter = %d kl = %f\n',gamma,beta,sum(m),iter,kl1);
end;

if beta>=beta_max
 fprintf('-----------------------------------------------------------\n');
 fprintf('beta > beta_max (%.3f > %.3f)\n', beta, beta_max);
 if i>1
    m = squeeze(m_all(i-1,1,:))';
 end
end
if sum(m)>=max_sum_m
 fprintf('-----------------------------------------------------------\n');
 fprintf('sum(m) > max_sum_m (%.3f > %.3f)\n', sum(m), max_sum_m);
end

% for gamma is current gamma decreasing to gamma_min
imax=i-1;
for i=imax:-1:1,
	gamma=gamma_all(i);
    eval(method);
	v_all(i,2,:)=v;
	m_all(i,2,:)=m;
	beta_all(i,2)=beta;
	iter_all(i,2)=iter;
	kl_all(i,2)=kl1;
	fprintf('gamma = %f beta = %f sum(m) = %f iter = %d kl = %f\n',gamma,beta,sum(m),iter,kl1);
end;

% select for each gamma from these two solutions the one with lowest KL
[klmin, imin]=min(kl_all,[],2);
for i=1:imax,
	v_mf_all(i,:)=squeeze(v_all(i,imin(i),:))';
	m_mf_all(i,:)=squeeze(m_all(i,imin(i),:))';
	beta_mf_all(i)=beta_all(i,imin(i));
	error_mf_all(i)=1/p*sum((y-v_mf_all(i,:)*x).^2,2);
	errorv_mf_all(i)=1/pv*sum((yv-v_mf_all(i,:)*xv).^2,2);
end;

% select the gamma that optimizes the validation error (errorv_mf_all)
[minerrorv i]=min(errorv_mf_all(1:imax));

res.gamma_mf=gamma_all(i);
res.v_mf=v_mf_all(i,:);
res.m_mf=m_mf_all(i,:);
res.n_mf1=sum(res.m_mf>0.5);

res.error_mf=error_mf_all(i);
res.errorv_mf=errorv_mf_all(i);
res.beta_mf=beta_mf_all(i);

if (res.beta_mf==beta_max)
	fprintf('beta_max too small: beta_mf %6.4f, beta_max %6.4f\n', beta_mf(iruns),beta_max);
	pause
end;
if (res.gamma_mf==gamma_min)
	fprintf('gamma at minimum range boundary\n');
end;
if (res.gamma_mf==gamma_max)
	fprintf('gamma at maxium range boundary\n');
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cp, vs, i1, iN] = params_legacy_cost(om, name, idx)
%PARAMS_LEGACY_COST  Returns cost parameters for legacy user-defined costs.
%   CP = OM.PARAMS_LEGACY_COST()
%   CP = OM.PARAMS_LEGACY_COST(NAME)
%   CP = OM.PARAMS_LEGACY_COST(NAME, IDX_LIST)
%   [CP, VS] = OM.PARAMS_LEGACY_COST(...)
%   [CP, VS, I1, IN] = OM.PARAMS_LEGACY_COST(...)
%
%   With no input parameters, it assembles and returns the parameters
%   for the aggregate legacy user-defined cost from all legacy cost sets
%   added using ADD_LEGACY_COST. The values of these parameters are cached
%   for subsequent calls. The parameters are contained in the struct CP,
%   described below.
%
%   If a NAME is provided then it simply returns parameter struct CP for the
%   corresponding named set. Likewise for indexed named sets specified
%   by NAME and IDX_LIST.
%
%   An optional 2nd output argument VS indicates the variable sets used by
%   this cost set. The size of CP.N will be consistent with VS.
%
%   If NAME is provided, optional 3rd and 4th output arguments I1 and IN
%   indicate the starting and ending row indices of the corresponding
%   cost set in the full aggregate cost matrix CP.N.
%
%   Let X refer to the vector formed by combining the corresponding varsets VS,
%   and F_U(X, CP) be the cost at X corresponding to the cost parameters
%   contained in CP, where CP is a struct with the following fields:
%       N      - nw x nx sparse matrix (optional, identity matrix by default)
%       Cw     - nw x 1 vector
%       H      - nw x nw sparse matrix (optional, all zeros by default)
%       dd, mm - nw x 1 vectors (optional, all ones by default)
%       rh, kk - nw x 1 vectors (optional, all zeros by default)
%
%   These parameters are used as follows to compute F_U(X, CP)
%
%       R  = N*X - rh
%
%               /  kk(i),  R(i) < -kk(i)
%       K(i) = <   0,     -kk(i) <= R(i) <= kk(i)
%               \ -kk(i),  R(i) > kk(i)
%
%       RR = R + K
%
%       U(i) =  /  0, -kk(i) <= R(i) <= kk(i)
%               \  1, otherwise
%
%       DDL(i) = /  1, dd(i) = 1
%                \  0, otherwise
%
%       DDQ(i) = /  1, dd(i) = 2
%                \  0, otherwise
%
%       Dl = diag(mm) * diag(U) * diag(DDL)
%       Dq = diag(mm) * diag(U) * diag(DDQ)
%
%       w = (Dl + Dq * diag(RR)) * RR
%
%       F_U(X, CP) = 1/2 * w'*H*w + Cw'*w
%
%   See also OPT_MODEL, ADD_LEGACY_COST, EVAL_LEGACY_COST.

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

if nargin > 1       %% individual set
    if nargin < 3
        idx = {};
    end
    if isempty(idx)                 %% name, no index provided
        if numel(om.cost.idx.i1.(name)) == 1    %% simple named set
            N  = om.cost.data.N.( name);
            Cw = om.cost.data.Cw.(name);
            [nw, nx] = size(N);
            e1 = ones(nw, 1);
            e0 = zeros(nw, 1);
            cp = struct(...
                'N',    N, ...
                'Cw',   Cw, ...
                'H',    sparse(nw, nw), ... %% default => no quadratic term
                'dd',   e1, ...             %% default => linear
                'rh',   e0, ...             %% default => no shift
                'kk',   e0, ...             %% default => no dead zone
                'mm',   e1 );               %% default => no scaling
            if isfield(om.cost.data.H, name) && ~isempty(om.cost.data.H.(name))
                cp.H = om.cost.data.H.(name);
            end
            if isfield(om.cost.data.dd, name) && ~isempty(om.cost.data.dd.(name))
                cp.dd = om.cost.data.dd.(name);
            end
            if isfield(om.cost.data.rh, name) && ~isempty(om.cost.data.rh.(name))
                cp.rh = om.cost.data.rh.(name);
            end
            if isfield(om.cost.data.kk, name) && ~isempty(om.cost.data.kk.(name))
                cp.kk = om.cost.data.kk.(name);
            end
            if isfield(om.cost.data.mm, name) && ~isempty(om.cost.data.mm.(name))
                cp.mm = om.cost.data.mm.(name);
            end
            if nargout > 1
                vs = om.cost.data.vs.(name);
                if nargout > 3
                    i1 = om.cost.idx.i1.(name);     %% starting row index
                    iN = om.cost.idx.iN.(name);     %% ending row index
                end
            end
        else                                    %% indexing required
            error('@opt_model/params_legacy_cost: legacy cost set ''%s'' requires an IDX_LIST arg', name);
        end
    else                            %% indexed named set
        % (calls to substruct() are relatively expensive ...
        % s = substruct('.', name, '{}', idx);
        % ... so replace it with these more efficient lines)
        sc = struct('type', {'.', '{}'}, 'subs', {name, idx});
        N  = subsref(om.cost.data.N,  sc);
        Cw = subsref(om.cost.data.Cw, sc);
        [nw, nx] = size(N);
        e1 = ones(nw, 1);
        e0 = zeros(nw, 1);
        cp = struct(...
            'N',    N, ...
            'Cw',   Cw, ...
            'H',    sparse(nw, nw), ... %% default => no quadratic term
            'dd',   e1, ...             %% default => linear
            'rh',   e0, ...             %% default => no shift
            'kk',   e0, ...             %% default => no dead zone
            'mm',   e1 );               %% default => no scaling
        if isfield(om.cost.data.H, name) && ~isempty(subsref(om.cost.data.H, sc))
            cp.H = subsref(om.cost.data.H,  sc);
        end
        if isfield(om.cost.data.dd, name) && ~isempty(subsref(om.cost.data.dd, sc))
            cp.dd = subsref(om.cost.data.dd, sc);
        end
        if isfield(om.cost.data.rh, name) && ~isempty(subsref(om.cost.data.rh, sc))
            cp.rh = subsref(om.cost.data.rh, sc);
        end
        if isfield(om.cost.data.kk, name) && ~isempty(subsref(om.cost.data.kk, sc))
            cp.kk = subsref(om.cost.data.kk, sc);
        end
        if isfield(om.cost.data.mm, name) && ~isempty(subsref(om.cost.data.mm, sc))
            cp.mm = subsref(om.cost.data.mm, sc);
        end
        if nargout > 1
            vs = subsref(om.cost.data.vs, sc);
            if nargout > 3
                sn = sc; sn(2).type = '()';         %% num array field
                i1 = subsref(om.cost.idx.i1, sn);   %% starting row index
                iN = subsref(om.cost.idx.iN, sn);   %% ending row index
            end
        end
    end
else                %% aggregate
    cache = om.cost.params;
    if isempty(cache)       %% build the aggregate
        nx = om.var.N;          %% number of variables
        nw = om.cost.N;         %% number of cost rows
        Nt = sparse(nx, nw);    %% use N transpose for speed
        Cw = zeros(nw, 1);
        H = sparse(nw, nw);     %% default => no quadratic term
        dd = ones(nw, 1);       %% default => linear
        rh = Cw;                %% default => no shift
        kk = Cw;                %% default => no dead zone
        mm = dd;                %% default => no scaling

        %% fill in each piece
        for k = 1:om.cost.NS
            name = om.cost.order(k).name;
            idx  = om.cost.order(k).idx;
            [cp, vs, i1, iN] = om.params_legacy_cost(name, idx);
            [mk, nk] = size(cp.N);      %% size of Nk
            if mk
                Nkt_full = sparse(nx, nw);
                if isempty(vs)
                    if nk == nx     %% full size
                        Nkt_full(:, i1:iN) = cp.N';
                    else            %% vars added since adding this cost set
                        Nk_all_cols = sparse(mk, nx);
                        Nk_all_cols(:, 1:nk) = cp.N;
                        Nkt_full(:, i1:iN) = Nk_all_cols';
                    end
                else
                    jj = om.varsets_idx(vs);    %% indices for var set
                    Nk_all_cols = sparse(mk, nx);
                    Nk_all_cols(:, jj) = cp.N;
                    Nkt_full(:, i1:iN) = Nk_all_cols';
                end
                Nt = Nt + Nkt_full;
                Cw(i1:iN) = cp.Cw;
                H_all_cols = sparse(mk, nw);
                H_all_cols(:, i1:iN) = cp.H';
                H(:, i1:iN) = H_all_cols';
                dd(i1:iN) = cp.dd;
                rh(i1:iN) = cp.rh;
                kk(i1:iN) = cp.kk;
                mm(i1:iN) = cp.mm;
            end
        end

        %% cache aggregated parameters
        om.cost.params = struct( ...
            'N', Nt', 'Cw', Cw, 'H', H, 'dd', dd, 'rh', rh, 'kk', kk, 'mm', mm );
        cp = om.cost.params;
    else                    %% return cached values
        cp = cache;
    end
    if nargout > 1
        vs = {};
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function au_qmr_ilu_s_d(p)
%
%  Deletes the global data that has been generated by 'au_qmr_ilu_s_i'
%  or 'au_qmr_ilu_s'.
%
%  Calling sequence:
%
%    au_qmr_ilu_s_d
%
%  Input:
%
%    p         vector containing the ADI shift parameters.
%
%
%  LYAPACK 1.0 (Thilo Penzl, August 1999)

l = length(p);

                     % Delete them regardless if they have been
                     % created at all.
for i = 1:l

  eval(lp_e( 'global LP_L',i,' LP_U',i ));

  eval(lp_e( 'clear global LP_L',i,' LP_U',i ));

end

global LP_P
clear global LP_P

author: Sable

repo: mcbench-benchmarks


text:
%  Figure 3.7      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.7
% fig3_07.m
% Satellite Double-Pulse Response Example 3.19
clf;
dI=1/5000;
numG=dI;
denG=[1 0 0];
sys=tf(numG,denG);
t=0:0.01:10;
% pulse input
u2=[zeros(1,500), 25*ones(1,10), zeros(1,100), -25*ones(1,10), zeros(1,381)];
[y2]=lsim(sys,u2,t);
figure();
plot(t,u2);
axis([0, 10, -26, 26]);
xlabel('Time (sec)');
ylabel('Thrust Fc');
title('Fig. 3.7(a): Thrust input')
% grid
nicegrid
pause;
% conversion to degrees
ff=180/pi;
y2=ff*y2;
figure();
plot(t,y2);
xlabel('Time (sec)');
ylabel('\theta (deg)');
title('Fig. 3.7(b): Satellite attitude')
% grid
nicegrid