text:
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESCENT(X, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters
    theta1=theta(1,1)-alpha/m*sum(X*theta-y);
    theta2=theta(2,1)-alpha/m*sum((X*theta-y) .* X(:,2));
    theta(1,1)=theta1;
    theta(2,1)=theta2;

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %

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

    % Save the cost J in every iteration
    J_history(iter) = computeCost(X, y, theta);

end

end

author: loserChen

repo: Coursera-MachineLearning


text:
function view_tensors(D, HA, delta, dim_order)
% VIEW_TENSORS    Renders diffusion tensors in 3D, colour coded by HA
%
%
% Inputs:
%
%   D is the diffusion tensor from fit_DT in 2 dimensions only (i.e. a 3D
%   volume with [x, y, DT])
%
%   HA is the helix angle (in 2D). HA values of NaN are coloured grey.
%
%   delta is the spacing between tensors (default 1.5E-3)
%
%   dim_order can be used (for example) when the tensors are referenced to
%   sagittal slices, whereas you want to display an axial slice. (default
%   [1 2 3])

% Author: Darryl McClymont <darryl.mcclymont@gmail.com>
% Copyright © 2014-2015 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% This function is an extension of plotDTI in the fanDTasia toolbox

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

if nargin < 4
    dim_order = [1 2 3]; % dimension
end
if nargin < 3
    delta = 1.5E-3; % spacing between tensors
end

% convert from 6x1 to 3x3 tensor representation
DT_33 = zeros([3 3 size(D,1) size(D,2)]);

DT_33(dim_order(1),dim_order(1),:,:) = D(:,:,1);
DT_33(dim_order(1),dim_order(2),:,:) = D(:,:,2);
DT_33(dim_order(1),dim_order(3),:,:) = D(:,:,3);
DT_33(dim_order(2),dim_order(1),:,:) = D(:,:,2);
DT_33(dim_order(2),dim_order(2),:,:) = D(:,:,4);
DT_33(dim_order(2),dim_order(3),:,:) = D(:,:,5);
DT_33(dim_order(3),dim_order(1),:,:) = D(:,:,3);
DT_33(dim_order(3),dim_order(2),:,:) = D(:,:,5);
DT_33(dim_order(3),dim_order(3),:,:) = D(:,:,6);

sz=size(DT_33);
if length(sz)==2
    nx=1;ny=1;
elseif length(sz)==3
    nx=sz(3);ny=1;
elseif length(sz)==4
    nx=sz(3);ny=sz(4);
end

ha = HA / 90;

hold on;
for i=1:nx
    for j=1:ny

        % don't waste time on nulled voxels
        if DT_33(1,1,i,j) == 0
            continue
        end

        % Jet colourmap
        if isnan(ha(i,j))
            C = [0.8 0.8 0.8];
        else
            C = [red(ha(i,j)), green(ha(i,j)), blue(ha(i,j))];
        end

        % get the eigenvectors of the tensor
        [v,l]=eig(round(DT_33(:,:,i,j)*1E7)/1E7);

        % generate the ellipsoid
        [X,Y,Z]=ellipsoid(0,0,0,l(1,1),l(2,2),l(3,3),10);
        sz=size(X);
        for x=1:sz(1)
            for y=1:sz(2)
                A=[X(x,y) Y(x,y) Z(x,y)]';
                A=v*A;
                X(x,y)=A(1);
                Y(x,y)=A(2);
                Z(x,y)=A(3);
            end
        end
        X=X+(i-1)*delta*2;
        Y=Y+(j-1)*delta*2;

        % display the ellipsoid
        surf(real(X),real(Y),real(Z), 'faceColor', C, 'EdgeAlpha', 0, 'EdgeColor', 'none');

    end
end

axis equal
view([0 90]);
set(gca,'GridLineStyle','none')
set(gca,'XTick',[])
set(gca,'YTick',[])
set(gca,'ZTick',[])

view([-1.2 0.5 0.8])
lighting phong
light('Position',[0 0 1],'Style','infinite','Color',[ 0.8 0.8 0.8]);

end

% these functions are to get the 'jet' colourmap (thanks stackexchange)
function val_out = interpolate( val,  y0,  x0,  y1,  x1 )
    val_out = (val-x0).*(y1-y0)./(x1-x0) + y0;
end

function val_out = base( val )
    if ( val <= -0.75 )
        val_out = 0;
    elseif ( val <= -0.25 )
         val_out = interpolate( val, 0.0, -0.75, 1.0, -0.25 );
    elseif ( val <= 0.25 )
        val_out =  1.0;
    elseif ( val <= 0.75 )
        val_out = interpolate( val, 1.0, 0.25, 0.0, 0.75 );
    else
        val_out = 0.0;
    end
end

function v = red( gray )
    v = base( gray - 0.5 );
end

function v = green( gray )
    v = base( gray );
end

function v = blue( gray )
    v = base( gray + 0.5 );
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


close all; clear all

imgname = 'trace-without-GP';
outfile = fullfile(pwd,'/simresults/20-01-15-out-GP-without-GP.mat');
kplot = 558;

imgname = 'trace-with-GP';
outfile = fullfile(pwd,'/simresults/20-01-15-out-GP-with-GP-optimized.mat');
kplot = 412

load(outfile)

[trackdata, x0, th0, w] = RaceTrack.loadTrack_02();
track = RaceTrack(trackdata, x0, th0, w);

trackAnim = SingleTrackAnimation(track, out.mu_x_pred_opt, out.var_x_pred_opt, out.u_pred_opt, out.x_ref);
trackAnim.initTrackAnimation();
drawnow;

% k = find(~isnan(out.xhat(1,:)), 1, 'last' ) - 1;

trackAnim.mu_x_pred_opt  = out.mu_x_pred_opt;
trackAnim.var_x_pred_opt = out.var_x_pred_opt;
trackAnim.u_pred_opt     = out.u_pred_opt;
trackAnim.x_ref          = out.x_ref;
trackAnim.updateTrackAnimation(kplot); % 558

delete(trackAnim.h_car)
delete(trackAnim.h_x_ref)
delete(trackAnim.h_mu_x_pred_opt)
cellfun( @(h) delete(h), trackAnim.h_var_x_pred_opt )
legend('off')
title('')
fp.savefig(imgname,'format','epsc')

%% CHECK RELAXED BARRIER FUNCTION PLOT

x = -0.5:0.001:0.5;
% Relaxied barrier function for (x<=lambda)
gamma = 10000;
lambda = -0.2;
y = 0.5*(sqrt((4+gamma*(lambda-x).^2)/gamma) - (lambda-x));
figure('Color','w'); hold on; %grid on;
plot(x,y,'LineWidth',2)
ylim([0,0.1])
xlim([-0.4,0])
set(gca,'YTickLabel',[]);
set(gca,'XTickLabel',[]);
%%%%
lambda = -0.2;
x = -0.6:0.001:(lambda-eps);
y = 0.1*-log(lambda-x);
figure('Color','w'); hold on; %grid on;
plot(x,y,'LineWidth',2)
xlim([-0.6,0])


text:
function cascade_data(model, data_year, pca)

% cascade_data(model, data_year, pca)
%
% Compute score statistics for filters and deformation models
% using the model and training data associated with the model
% and dataset year 'data_year'.  If PCA is given a value > 0, then
% the score statistics are computed using a PCA projection of the
% model and training data.
%
% The score statistics are written in to the file
% class_year_cascade_data_pcaK_data_year.inf, which is used
% by cascade_model.m.
%
% model      object detector
% data_year  dataset year as a string (e.g., '2007')
% pca        number of PCA components to project onto (if pca > 0)

setVOCyear = model.year;
globals;
pascal_init;

% if using PCA, project the model
if pca > 0
  load('pca.mat');
  [ignore, model] = projectmodel(model, coeff, pca);
end

% files used by scorestat.cc
pca_str = num2str(pca);
class_year = [model.class '_' model.year];
detfile = [cscdir class_year '_cascade_data_det_' data_year '.inf'];
inffile = [cscdir class_year '_cascade_data_pca' pca_str '_' data_year '.inf'];

% load training images
pos = pascal_data(model.class, true, data_year);
if pca == 0
  % file for recording the detection location and component of each detection
  detfid = fopen(detfile, 'w');
  validinfo = [];
else
  % read the detection location and component used for each detection by the
  % non-PCA model
  pattern = ['%n%n%n' repmat('%n', 1, model.numblocks)];
  detfid = fopen(detfile);
  D = textscan(detfid, '%d %d %d %d %d', 'Delimiter', ' ');
  fclose(detfid);
  validinfo = [D{2} D{3} D{4} D{5}];
end

% reserve an int at the beginning of the file
inffid = fopen(inffile, 'wb');
fwrite(inffid, 0, 'int32');
% compute detections on training images
num = process(model, pos, validinfo, pca, detfid, inffid);
% rewind and write sample size
frewind(inffid);
fwrite(inffid, num, 'int32');
fclose(inffid);
if pca == 0
  fclose(detfid);
end
fprintf('Wrote score statistics for %d examples.\n', num);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Write detection info file and detection feature vector file.
function num = process(model, pos, validinfo, pca, detfid, inffid)

numpos = length(pos);
pixels = model.minsize * model.sbin;
minsize = prod(pixels);

num = 0;
batchsize = 16;
% process positive examples in parallel batches
for j = 1:batchsize:numpos
  % do batches of detections in parallel
  thisbatchsize = batchsize - max(0, (j+batchsize-1) - numpos);
  % data for this batch
  pardata = {};
  parfor k = 1:thisbatchsize
    i = j+k-1;
    fprintf('%s %s: cascade data (PCA=%d): %d/%d', procid(), model.class, pca, i, numpos);
    bbox = [pos(i).x1 pos(i).y1 pos(i).x2 pos(i).y2];
    valid = struct('c',0,'x',0,'y',0,'l',0);
    % skip small examples
    if (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1) < minsize
      pardata{k}.det = [];
      if pca == 0
        pardata{k}.detinfo = valid;
      end
      fprintf(' (too small)\n');
      continue;
    end
    % get example
    im = imreadx(pos(i));
    [im, bbox] = croppos(im, bbox);
    pyra = featpyramid(im, model);
    if pca > 0
      pyra = projectpyramid(model, pyra);
      valid.c = validinfo(i,1);
      valid.x = validinfo(i,2);
      valid.y = validinfo(i,3);
      valid.l = validinfo(i,4);
    end
    [det, boxes, info] = gdetectvalid(pyra, model, 0, bbox, 0.5, valid);
    pardata{k}.detinfo = getdetinfo(model, info);
    pardata{k}.det = det;
    pardata{k}.boxes = boxes;
    pardata{k}.info = info;
    pardata{k}.pyra = pyra;

    if ~isempty(det)
      fprintf(' (%f)\n', det(end));
      num = num + 1;
      %showboxes(im, det);
    else
      fprintf(' (skip)\n');
    end
  end
  % write data to disk sequentially
  for k = 1:thisbatchsize
    i = j+k-1;
    if ~isempty(pardata{k}.det)
      writescores(pardata{k}.pyra, model, pardata{k}.info, inffid);
    end
    if pca == 0
      % record the detection location and component for this detection
      detinfo = pardata{k}.detinfo;
      fprintf(detfid, '%d %d %d %d %d\n', ...
              i, detinfo.c, detinfo.x, detinfo.y, detinfo.l);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Retrieve the (x,y,scale) location of the detection in feature
% pyramid coordinates.
function d = getdetinfo(model, info)

d.c = 0;
d.x = 0;
d.y = 0;
d.l = 0;

if ~isempty(info)
  % See: gdetectwrite.m and getdetections.cc
  DET_IND = 2;    % rule index (i.e. component number)
  DET_X   = 3;    % x coord (filter and deformation)
  DET_Y   = 4;    % y coord (filter and deformation)
  DET_L   = 5;    % level (filter)

  ruleind = info(DET_IND, model.start, 1);
  d.c = ruleind;
  d.x = info(DET_X, model.start, 1);
  d.y = info(DET_Y, model.start, 1);
  d.l = info(DET_L, model.start, 1);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function blas0_test02 ( )

%*****************************************************************************80
%
%% BLAS0_TEST02 tests R8_ABS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2014
%
%  Author:
%
%    John Burkardt
%
  r8_hi = 5.0;
  r8_lo = -3.0;
  test_num = 10;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLAS0_TEST02\n' );
  fprintf ( 1, '  R8_ABS returns the absolute value of an R8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X         R8_ABS(X)\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    [ r8, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
    r8_absolute = r8_abs ( r8 );
    fprintf ( 1, '  %10.6f  %10.6f\n', r8, r8_absolute );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% OFDM_basic.m
clear all
NgType=1; % NgType=1/2 for cyclic prefix/zero padding
if NgType==1
    nt='CP';
elseif NgType==2
    nt='ZP';
end
Ch=0;  % Ch=0/1 for AWGN/multipath channel
if Ch==0
    chType='AWGN';
    Target_neb=100;
else
    chType='CH';
    Target_neb=500;
end
figure(Ch+1), clf
PowerdB=[0 -8 -17 -21 -25]; % Channel tap power profile 'dB'
Delay=[0 3 5 6 8];          % Channel delay 'sample'
Power=10.^(PowerdB/10);     % Channel tap power profile 'linear scale'
Ntap=length(PowerdB);       % Chanel tap number
Lch=Delay(end)+1;           % Channel length
Nbps=4;
M=2^Nbps;                   % Modulation order=2/4/6 for QPSK/16QAM/64QAM
Nfft=64;                    % FFT size
Ng=Nfft/4;                  % Guard interval length
Nsym=Nfft+Ng;               % Symbol duration
Nvc=Nfft/4;                 % Nvc=0: no virtual carrier
Nused=Nfft-Nvc;
EbN0=[0:5:20];    % EbN0
N_iter=1e5;       % Number of iterations for each EbN0
Nframe=3;         % Number of symbols per frame
sigPow=0;         % Signal power initialization
file_name=['OFDM_BER_' chType '_' nt '_' 'GL' num2str(Ng) '.dat'];
fid=fopen(file_name, 'w+');
norms=[1 sqrt(2) 0 sqrt(10) 0 sqrt(42)];     % BPSK 4-QAM 16-QAM
for i=0:length(EbN0)
   randn('state',0);
   rand('state',0);
   %Ber2=ber(); % BER initialization
   Neb=0;
   Ntb=0; % Initialize the number of error/total bits
   for m=1:N_iter
      % Tx______________________________________________________________
      X= randint(1,Nused*Nframe,M); % bit: integer vector
      Xmod= qammod(X,M,0,'gray')/norms(Nbps);
      if NgType~=2
          x_GI=zeros(1,Nframe*Nsym);
       elseif NgType==2
           x_GI= zeros(1,Nframe*Nsym+Ng);
        % Extend an OFDM symbol by Ng zeros
      end
      kk1=[1:Nused/2];
      kk2=[Nused/2+1:Nused];
      kk3=1:Nfft;
      kk4=1:Nsym;
      for k=1:Nframe
         if Nvc~=0
             X_shift= [0 Xmod(kk2) zeros(1,Nvc-1) Xmod(kk1)];
         else
             X_shift= [Xmod(kk2) Xmod(kk1)];
         end
         x= ifft(X_shift);
         x_GI(kk4)= guard_interval(Ng,Nfft,NgType,x);
         kk1=kk1+Nused;
         kk2=kk2+Nused;
         kk3=kk3+Nfft;
         kk4=kk4+Nsym;
      end
      if Ch==0
          y= x_GI;  % No channel
      else  % Multipath fading channel
        channel=(randn(1,Ntap)+j*randn(1,Ntap)).*sqrt(Power/2);
        h=zeros(1,Lch);
        h(Delay+1)=channel; % cir: channel impulse response
        y = conv(x_GI,h);
      end
      if i==0 % Only to measure the signal power for adding AWGN noise
        y1=y(1:Nframe*Nsym);
        sigPow = sigPow + y1*y1';
        continue;
      end
      % Add AWGN noise________________________________________________
      snr = EbN0(i)+10*log10(Nbps*(Nused/Nfft)); % SNR vs. Eb/N0
      noise_mag = sqrt((10.^(-snr/10))*sigPow/2);
      y_GI = y + noise_mag*(randn(size(y))+j*randn(size(y)));
      % Rx_____________________________________________________________
      kk1=(NgType==2)*Ng+[1:Nsym];
      kk2=1:Nfft;
      kk3=1:Nused;
      kk4=Nused/2+Nvc+1:Nfft;
      kk5=(Nvc~=0)+[1:Nused/2];
      if Ch==1
         H= fft([h zeros(1,Nfft-Lch)]); % Channel frequency response
         H_shift(kk3)= [H(kk4) H(kk5)];
      end
      for k=1:Nframe
         Y(kk2)= fft(remove_GI(Ng,Nsym,NgType,y_GI(kk1)));
         Y_shift=[Y(kk4) Y(kk5)];
         if Ch==0
             Xmod_r(kk3) = Y_shift;
         else
             Xmod_r(kk3)= Y_shift./H_shift;  % Equalizer - channel compensation
         end
         kk1=kk1+Nsym;
         kk2=kk2+Nfft;
         kk3=kk3+Nused;
         kk4=kk4+Nfft;
         kk5=kk5+Nfft;
      end
      X_r=qamdemod(Xmod_r*norms(Nbps),M,0,'gray');
      Neb=Neb+sum(sum(de2bi(X_r,Nbps)~=de2bi(X,Nbps)));
      Ntb=Ntb+Nused*Nframe*Nbps;  %[Ber,Neb,Ntb]=ber(bit_Rx,bit,Nbps);
      if Neb>Target_neb
          break;
      end
   end
   if i==0
     sigPow= sigPow/Nsym/Nframe/N_iter;
     fprintf('Signal power= %11.3e\n', sigPow);
     fprintf(fid,'%%Signal power= %11.3e\n%%EbN0[dB]       BER\n', sigPow);
    else
     Ber = Neb/Ntb;
     fprintf('EbN0=%3d[dB], BER=%4d/%8d =%11.3e\n', EbN0(i), Neb,Ntb,Ber)
     fprintf(fid, '%d\t%11.3e\n', EbN0(i), Ber);
     if Ber<1e-6
         break;
     end
   end
end
if (fid~=0)
    fclose(fid);
end
disp('Simulation is finished');
plot_ber(file_name,Nbps);

author: LyricYang


text:
function ChannelMat = in_channel_curry_rs3(ChannelFile)
% IN_CHANNEL_CURRY_RS3:  Read 3D cartesian positions for a set of electrodes from Curry .rs3 file.
%
% USAGE:  ChannelMat = in_channel_curry_res3(ChannelFile)
%
% INPUTS:
%     - ChannelFile : Full path to the file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2011

% Open file
fid = fopen(ChannelFile, 'r');
if (fid == -1)
    error('Cannot open file.');
end
% Initialize indices structure
iChannel = 1;
curBlock = '';
ChannelMat = db_template('channelmat');
ChannelMat.Comment = 'Curry rs3';

% Read file line by line
while 1
    % Read line
    read_line = fgetl(fid);
    % Empty line: go to next line
    if isempty(read_line)
        continue
    end
    % End of file: stop reading
    if (read_line(1) == -1)
        break
    end

    % Check if beginning/end of block

    if ~isempty(regexp(read_line, 'SENSORS\w* START_LIST'))
        curBlock = 'pos';
        iChannel = 1;
    elseif ~isempty(regexp(read_line, 'LABELS\w* START_LIST'))
        curBlock = 'labels';
        iChannel = 1;
    elseif ~isempty(regexp(read_line, 'SENSORS\w* END_LIST')) || ~isempty(regexp(read_line, 'LABELS\w* END_LIST'))
        curBlock = '';
    else
        switch (curBlock)
            case 'labels'
                ChannelMat.Channel(iChannel).Name = read_line;
                iChannel = iChannel + 1;
            case 'pos'
                xyz = sscanf(read_line, '%f %f %f') / 1000;
                ChannelMat.Channel(iChannel).Type    = 'EEG';
                ChannelMat.Channel(iChannel).Loc     = [-xyz(2); xyz(1); xyz(3)];
                ChannelMat.Channel(iChannel).Orient  = [];
                ChannelMat.Channel(iChannel).Comment = '';
                ChannelMat.Channel(iChannel).Weight  = 1;
                iChannel = iChannel + 1;
        end
    end
end
% Close file
fclose(fid);

author: brainstorm-tools


text:
if ~isempty(strfind(fsub,'.mat'))
    dat = load(fsub);
    if isfield(dat,'X')
        X = dat.X;
    else
        X = getfield(dat,char(fieldnames(dat)));
    end
    clear dat
elseif ~isempty(strfind(fsub,'.txt'))
    X = dlmread(fsub);
else % no extension - assumed to be an SPM file
    try % Note that we are not taking care of bad samples here
        D = spm_eeg_load(fsub);
        X = permute(D(:,:,:),[2 3 1]); clear D
        X = reshape(X,[size(X,1)*size(X,2) size(X,3)]);
    catch
        if ~exist('spm_eeg_load','file')
            error('spm_eeg_load() function not found - is SPM in the path?')
        else
            error('Incorrect data format - input files must be .mat, .txt or an SPM file');
        end
    end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)

%function [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)
%
%% READ_RICOH_DATA reads continuous or averaged MEG data
%% generated by the RICOH MEG system and software,
%% and allows the data to be used in FieldTrip.
%%
%% Use as
%%   [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)
%%
%% This is a wrapper function around the function getRData
%%
%% See also READ_RICOH_HEADER, READ_RICOH_EVENT

if ~ft_hastoolbox('ricoh_meg_reader')
    ft_error('cannot determine whether Ricoh toolbox is present');
end

hdr = hdr.orig; % use the original Ricoh header, not the FieldTrip header

% default is to select all channels
if nargin<5
  chanindx = 1:hdr.channel_count;
end

handles = definehandles;

switch hdr.acq_type
  case handles.AcqTypeEvokedAve
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getRData(filename, start_sample, sample_length);

  case handles.AcqTypeContinuousRaw
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getRData(filename, start_sample, sample_length);

  % Unlike Yokogawa system, "AcqTypeEvokedRaw" is not supported for Ricoh system.

  otherwise
    ft_error('unknown data type');
end

if size(dat,1)~=hdr.channel_count
  ft_error('could not read all channels');
elseif size(dat,2)~=(endsample-begsample+1)
  ft_error('could not read all samples');
end

% select only the desired channels
dat = dat(chanindx,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: spm


text:
% bpExampleScript

ctol = 0.000001;  % convergence tolerance
T = 0.1;       % annealing temperature
maxiter = Inf; % max iter for maxBeliefPropBethe

% factor for P(x1 | x2)
factor2var{1} = [1 2]; % which variables the factor involves
factors{1} = [0.9 0.5 ; ...
              0.1 0.5]; % P(x1=1|x2=1) = 0.9, P(x1=2|x2=1)=0.1
% factor for P(x2 | x3)
factor2var{2} = [2 3];
factors{2} = [0.25 0.01 ; ...
              0.75 0.99];
% factor for P(x3)
factor2var{3} = [3];
factors{3} = [0.25 ; ...
              0.75];  % P(x1=1)=0.25, P(x1=2)=0.75

nvals = [2 ; 2 ; 2]; % number of states for each node (variable)

[vals_max, bel_max] = maxBeliefPropBethe(factors, factor2var, nvals, ctol, T, maxiter);
bel_max = reshape(cell2mat(bel_max), [nvals(1) numel(nvals)])';

bel_marg = marginalBeliefPropBethe(factors, factor2var, nvals, ctol);
bel_marg = reshape(cell2mat(bel_marg), [nvals(1) numel(nvals)])';

% marginal answer should be:
%   P(x3=1) = 0.25;
%   P(x2=1) = (0.25)*P(x3=1) + (0.01)*P(x3=2) =  0.07
%   P(x1=1) = (0.9)*P(x2=1) + 0.5*P(x2=2) = 0.528
bel_marg

% max solution should be:
%   x1 = 2;
%   x2 = 2;
%   x3 = 1 or 2; (equally likely)
bel_max

author: Cloud-CV

repo: object-proposals


text:
%%**********************************************************************
%% NTpred: Compute (dX,dy,dZ) for NT direction.
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************************

 function [par,dX,dy,dZ,coeff,L,hRd] = ...
          NTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);

    global schurfun schurfun_par
%%
%% compute NT scaling matrix
%%
    [par.W,par.G,par.sv,par.gamx,par.gamz,par.dd,par.ee,par.ff] = ...
     NTscaling(blk,X,Z,Zchol,invZchol);
%%
%% compute schur matrix
%%
    m = length(rp);
    schur = sparse(m,m);
    UU = []; EE = []; Afree = [];
    dX = cell(size(blk,1),1); dy = []; dZ = cell(size(blk,1),1);
%%
    for p = 1:size(blk,1)
       pblk = blk(p,:);
       if strcmp(pblk{1},'l')
          [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,par.dd);
       elseif strcmp(pblk{1},'q');
          [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,par.dd,par.ee);
       elseif strcmp(pblk{1},'s')
          if isempty(schurfun{p})
             schur = schurmat_sblk(blk,At,par,schur,p,par.W);
          elseif isstr(schurfun{p})
             schurtmp = sparse(m,m);
             if ~isempty(par.permZ{p})
                Wp = par.W{p}(par.permZ{p},par.permZ{p});
             else
                Wp = par.W{p};
             end
             eval(['schurtmp = ',schurfun{p},'(Wp,Wp,schurfun_par(p,:));']);
             schur = schur + schurtmp;
          end
       elseif strcmp(pblk{1},'u')
          Afree = [Afree, At{p}'];
       end
    end
%%
%% compute rhs
%%
    [rhs,EinvRc,hRd] = NTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu);
%%
%% solve linear system
%%
    [xx,coeff,L] = linsysolve(par,schur,UU,Afree,EE,rhs);
%%
%% compute (dX,dZ)
%%
    [dX,dy,dZ] = NTdirfun(blk,At,par,Rd,EinvRc,xx,m);
%%**********************************************************************

repo: REST


text:
function [names, ids] = getMultiSpeciesModelId(modelJoint, nameTagsModels, nameTagHost, metTagRe, rxnTagRe, compCom, compHost)
% Get names and IDs for metabolites and exchange reactions in the [u] and [b] space
%
% USAGE:
%
%    [names, ids] = getMultiSpeciesModelId(modelJoint, nameTagsModels, nameTagHost, metTagRe, rxnTagRe, compCom, compHost)
%
% INPUTS:
%    modelJoint:       COBRA multi-organism model
%    nameTagsModels:   cell array of tags for species to identify the respective
%                      reactions and metabolites from `modelJoint.rxns` and `.mets`
%
% OPTIONAL INPUTS:
%    nameTagHost:      string of tag for the host model if exist. Input [] if no host is present
%    metTagRe:         a regular expression to identify the tag in nameTagsModels from `modelJoint.mets`.
%                      Use '%s' for the tag in `nameTagsModels` for each species.
%                      Default '^%s', the beginning of the string. In this case, if the tag for species 1 is
%                      'SP1', then mets with id 'SP1glc-D[e]' will be identified as belonging to species 1.
%                      E.g., 'met[compartment_SP1]' where 'SP1' is the tag for species 1, input '\[[^_]+_%s\]$'.
%    rxnTagRe:         a regular expression to idenify the tag in `nameTagsModels` from `modelJoint.rxns`.
%                      Default '^%s', the beginning of the string
%    compCom:          compartment Id for the community exchange compartment. Default 'u'
%    compHost:         compartment Id for the exchange compartment specific to the host. Default 'b'
%
% OUTPUTS:
%    names:            structure of reaction/metabolite name (.rxns/.mets) with the following fields:
%
%                        * spAbbr - species abbreviation (= `nameTags`). Host put at the end
%                        * EXcom - #met[u]-by-1 cell. Exchange reactions for community metabolites `met[u]`
%                        * EXhost - #met[b]-by-1 cell. Exchange reactions for host-specific exchange metabolites `met[b]`
%                        * EXsp - #met[u]-by-#species cell. `EXsp(i,j)` is the species-community exchange reactions
%                          for the i-th met[u] and the j-th species
%                        * Mcom - #met[u]-by-1 cell. All community metabolites `met[u]`
%                        * Mhost - #met[u]-by-1 cell. Host-specific exchange metabolites `met[b]`
%                        * Msp - #met[u]-by-#species cell. `Msp(i,j)` is the metabolite in `met[e]` of the `j`-th species
%                          participating in reaction `EXsp(i,j)`
%                        * rxnSps - #rxns-by-1 cell. `rxnSps(i)` would be `spAbbr(j)` if the `i`-th reaction belongs to the `j`-th species
%                        * metSps - #mets-by-1 cell. `metSps(i)` would be `spAbbr(j)` if the i-th metabolite belongs to the `j`-th species
%
%    ids:              structure of reaction/metabolite index with the same fields as names except without `spAbbr`

if nargin < 7 || isempty(compHost)
    compHost = 'b';
end
if nargin < 6 || isempty(compCom)
    compCom = 'u';
end
compHost = regexprep(compHost, '^\[(.*)\]$', '$1');
compCom = regexprep(compCom, '^\[(.*)\]$', '$1');
% regular expression for identifying species
if nargin < 5 || isempty(rxnTagRe)
    rxnTagRe = '^%s';
end
if nargin < 4 || isempty(metTagRe)
    metTagRe = '^%s';
end
metCompartment = getCompartment(modelJoint.mets);
% mets in the common exchange space ([u]) or the biomass (use biomass[c] to be consistent with createMultipleSpeciesModel)
metCom = strcmp(metCompartment, compCom) | strcmp(modelJoint.mets, 'biomass[c]');
% mets in exchange space accessible only by the host ([b])
metHost = strcmp(metCompartment, compHost);
hostExist = any(metHost) && nargin >= 3 && ~isempty(nameTagHost);   % use && to avoid error if nameTagHost not given
if hostExist && iscell(nameTagHost)
    nameTagHost = nameTagHost{1};
end
if ischar(nameTagsModels)
    nameTagsModels = {nameTagsModels};
end
% get the regular expression of identifiers for species-specific mets and rxns
nSp = numel(nameTagsModels);
[metTagReSp, rxnTagReSp] = deal(repmat({''}, nSp + hostExist, 1));
for jSp = 1:nSp
    metTagReSp{jSp} = strrep(metTagRe, '%s', nameTagsModels{jSp});
    rxnTagReSp{jSp} = strrep(rxnTagRe, '%s', nameTagsModels{jSp});
end
metHost = find(metHost);
if hostExist
    % add regular expression for host
    metTagReSp{nSp + 1} = strrep(metTagRe, '%s', nameTagHost);
    rxnTagReSp{nSp + 1} = strrep(rxnTagRe, '%s', nameTagHost);
    % try to sort metHost to have the same order as metCom by comparing met names
    metHostName = regexprep(modelJoint.mets(metHost), metTagReSp{nSp + 1}, '');
    metHostName = strrep(metHostName, ['[' compHost ']'], '');
    metComName = strrep(modelJoint.mets(metCom), ['[' compCom ']'], '');
    [yn, id] = ismember(metComName, metHostName);
    if all(yn) && numel(union(metComName, metHostName)) == numel(metComName)
        % one-to-one mapping exists. Reorder
        metHost = metHost(id);
    end
end

% get system exchange rxns
rxnEX = find(sum(modelJoint.S ~= 0, 1) == 1);
% exchange rxns for mets in [u], in the same order as metCom
[EXcom, ~] = find(modelJoint.S(metCom, rxnEX)');
EXcom = rxnEX(EXcom);
% exchange rxns for mets in [b], in the same order as metHost
[EXhost, ~] = find(modelJoint.S(metHost, rxnEX)');
EXhost = rxnEX(EXhost);

% Get the species number for each met and rxn. 0 for those in [u]
[nMets, nRxns] = size(modelJoint.S);
metSps = zeros(nMets, 1);
rxnSps = zeros(nRxns, 1);
for jSp = 1:nSp
    metSps(~cellfun(@isempty, regexp(modelJoint.mets, metTagReSp{jSp}, 'once')) & ~metCom) = jSp;
    rxnSps(~cellfun(@isempty, regexp(modelJoint.rxns, rxnTagReSp{jSp}, 'once'))) = jSp;
    %     metSps(strncmp(modelJoint.mets, nameTagsModels{j}, numel(nameTagsModels{j})) & ~metCom) = j;
    %     rxnSps(strncmp(modelJoint.rxns, nameTagsModels{j}, numel(nameTagsModels{j}))) = j;
    rxnSps(EXcom) = 0;  % in case some mets or exchange rxns in [u] have prefix equal to the name tags
    % check if the model is properly compartmentalized.
    if nnz(modelJoint.S(metSps == jSp, rxnSps ~= jSp)) > 0 ...
            || nnz(modelJoint.S(metSps ~= jSp & metSps > 0, rxnSps == jSp)) > 0
        error('The model for species %s is not correctly compartmentalized. Check the name tags.');
    end
end
if hostExist
    metSps(~cellfun(@isempty, regexp(modelJoint.mets, metTagReSp{nSp + 1}, 'once')) & ~metCom) = nSp + 1;
    rxnSps(~cellfun(@isempty, regexp(modelJoint.rxns, rxnTagReSp{nSp + 1}, 'once'))) = nSp + 1;
    %     metSps(strncmp(modelJoint.mets, nameTagHost, numel(nameTagHost))) = nSp + 1;
    %     rxnSps(strncmp(modelJoint.rxns, nameTagHost, numel(nameTagHost))) = nSp + 1;
end

% species-community exchange rxns (between members' [e] and [u])
EXeu = any(modelJoint.S(metCom, :), 1);
EXeu(EXcom) = false;  % exclude system input/output exchange rxns
EXeu = find(EXeu);
% EXsp is a #met[u] x #species matrix with (i,j)-entry being the index for exchange reactionn
% for the i-th met in [u] between [u] and [e] of the j-th organism.
% Msp similar but the index for the exchanged metabolite in [e] of the organism
[EXsp, Msp] = deal(zeros(sum(metCom), nSp + hostExist));
for j = 1:numel(EXeu)
    mJcom = find(modelJoint.S(:, EXeu(j)) ~= 0 & metCom);  % met in [u]
    mJsp = find(modelJoint.S(:, EXeu(j)) ~= 0 & ~metCom);  % met in [e]
    EXsp(sum(metCom(1:mJcom)), rxnSps(EXeu(j))) = EXeu(j);
    Msp(sum(metCom(1:mJcom)), rxnSps(EXeu(j))) = mJsp;
end
% store all rxns/mets and their ids
names.spAbbr = nameTagsModels(:);  % species abbreviation
if hostExist
    names.spAbbr = [names.spAbbr; {nameTagHost}];
end
names.spName = names.spAbbr;
names.EXcom = modelJoint.rxns(EXcom);
ids.EXcom = EXcom(:);
names.EXhost = modelJoint.rxns(EXhost);
ids.EXhost = EXhost(:);
names.EXsp = repmat({''}, size(EXsp, 1), size(EXsp, 2));
names.EXsp(EXsp ~= 0) = modelJoint.rxns(EXsp(EXsp ~= 0));
ids.EXsp = EXsp;
names.Mcom = modelJoint.mets(metCom);
ids.Mcom = find(metCom);
names.Mhost = modelJoint.mets(metHost);
ids.Mhost = metHost;
names.Msp = repmat({''}, size(EXsp, 1), size(EXsp, 2));
names.Msp(Msp ~= 0) = modelJoint.mets(Msp(Msp ~= 0));
ids.Msp = Msp;
names.rxnSps = cell(nRxns, 1);
names.rxnSps(rxnSps > 0) = names.spAbbr(rxnSps(rxnSps > 0));
names.rxnSps(rxnSps == 0) = {'com'};
ids.rxnSps = rxnSps;
names.metSps = cell(nMets, 1);
names.metSps(metSps > 0) = names.spAbbr(metSps(metSps > 0));
names.metSps(metSps == 0) = {'com'};
ids.metSps = metSps;

author: opencobra

repo: cobratoolbox


text:
function [val]=gy(x,quant)

% indicator function for quantiles values
% returns 1 if the value 'x' is lower than or equal to the threshold value 'quant'
% returns 0 otherwise

val=x<=quant;

author: european-central-bank

repo: BEAR-toolbox


text:
function [varargout]=axis_inc(varargin)

% function [h]=axis_inc(s,ha)
% ------------------------------------------------------------------------
% This function widens the axis limits for the axis with the handle |hs|
% using the scale factor |s| (e.g. if s=2 the axis width is doubled). The
% scaling parameter |s| can be a single scalar or a vector such that each
% axis direction is scaled differently.
%
% Change log:
% 2019/06/26 Added variable input/output handling
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        s=varargin{1};
        ha=gca;
    case 2
        s=varargin{1};
        ha=varargin{2};
end

if numel(s)==1
    s=s*ones(1,3);
end

%% Set new axis limits

%Get axis limits
limSet={'XLim','YLim','ZLim'};

for q=1:1:numel(limSet)
    try
        limNow=get(ha,limSet{q});
        w=limNow(2)-limNow(1);
        if w<=eps(0)
            w=1; %Force 1 if too small
        end
        limInc=w*(s(q)-1)/2; %Axis limit increment
        limNew=limNow+[-limInc limInc];
        set(ha,limSet{q},limNew);
    catch

    end
end

%% Collect output
varargout{1}=ha;

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [net,stats] = cnn_train_dag(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-15 Andrea Vedaldi.
% 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).

% Todo: save momentum with checkpointing (a waste?)

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;
opts.profile = false ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

state.getBatch = getBatch ;
evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end
stats = [] ;

% setup GPUs
numGpus = numel(opts.gpus) ;
if numGpus > 1
  if isempty(gcp('nocreate')),
    parpool('local',numGpus) ;
    spmd, gpuDevice(opts.gpus(labindex)), end
  end
  if exist(opts.memoryMapFile)
    delete(opts.memoryMapFile) ;
  end
elseif numGpus == 1
  gpuDevice(opts.gpus)
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('resuming by loading epoch %d\n', start) ;
  [net, stats] = loadState(modelPath(start)) ;
end

for epoch=start+1:opts.numEpochs

  % train one epoch
  state.epoch = epoch ;
  state.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  state.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  state.val = opts.val ;
  state.imdb = imdb ;

  if numGpus <= 1
    stats.train(epoch) = process_epoch(net, state, opts, 'train') ;
    stats.val(epoch) = process_epoch(net, state, opts, 'val') ;
  else
    savedNet = net.saveobj() ;
    spmd
      net_ = dagnn.DagNN.loadobj(savedNet) ;
      stats_.train = process_epoch(net_, state, opts, 'train') ;
      stats_.val = process_epoch(net_, state, opts, 'val') ;
      if labindex == 1, savedNet_ = net_.saveobj() ; end
    end
    net = dagnn.DagNN.loadobj(savedNet_{1}) ;
    stats__ = accumulateStats(stats_) ;
    stats.train(epoch) = stats__.train ;
    stats.val(epoch) = stats__.val ;
    clear net_ stats_ stats__ savedNet_ ;
  end

  if ~evaluateMode
    saveState(modelPath(epoch), net, stats) ;
  end

  figure(1) ; clf ;
  plots = setdiff(...
    cat(2,...
        fieldnames(stats.train)', ...
        fieldnames(stats.val)'), {'num', 'time'}) ;
  for p = plots
    p = char(p) ;
    values = zeros(0, epoch) ;
    leg = {} ;
    for f = {'train', 'val'}
      f = char(f) ;
      if isfield(stats.(f), p)
        tmp = [stats.(f).(p)] ;
        values(end+1,:) = tmp(1,:)' ;
        leg{end+1} = f ;
      end
    end
    subplot(1,numel(plots),find(strcmp(p,plots))) ;
    plot(1:epoch, values','o-') ;
    xlabel('epoch') ;
    title(p) ;
    legend(leg{:}) ;
    grid on ;
  end
  drawnow ;
  print(1, modelFigPath, '-dpdf') ;
end

% -------------------------------------------------------------------------
function stats = process_epoch(net, state, opts, mode)
% -------------------------------------------------------------------------

if strcmp(mode,'train')
  state.momentum = num2cell(zeros(1, numel(net.params))) ;
end

numGpus = numel(opts.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  if strcmp(mode,'train')
    sate.momentum = cellfun(@gpuArray,state.momentum,'UniformOutput',false) ;
  end
end
if numGpus > 1
  mmap = map_gradients(opts.memoryMapFile, net, numGpus) ;
else
  mmap = [] ;
end

stats.time = 0 ;
stats.num = 0 ;
subset = state.(mode) ;
start = tic ;
num = 0 ;

for t=1:opts.batchSize:numel(subset)
  batchSize = min(opts.batchSize, numel(subset) - t + 1) ;

  for s=1:opts.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+opts.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    inputs = state.getBatch(state.imdb, batch) ;

    if opts.prefetch
      if s == opts.numSubBatches
        batchStart = t + (labindex-1) + opts.batchSize ;
        batchEnd = min(t+2*opts.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
      state.getBatch(state.imdb, nextBatch) ;
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (s ~= 1) ;
      net.eval(inputs, opts.derOutputs) ;
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % extract learning stats
  stats = opts.extractStatsFn(net) ;

  % accumulate gradient
  if strcmp(mode, 'train')
    if ~isempty(mmap)
      write_gradients(mmap, net) ;
      labBarrier() ;
    end
    state = accumulate_gradients(state, net, opts, batchSize, mmap) ;
  end

  % print learning statistics
  time = toc(start) ;
  stats.num = num ;
  stats.time = toc(start) ;

  fprintf('%s: epoch %02d: %3d/%3d: %.1f Hz', ...
    mode, ...
    state.epoch, ...
    fix(t/opts.batchSize)+1, ceil(numel(subset)/opts.batchSize), ...
    stats.num/stats.time * max(numGpus, 1)) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s:', f) ;
    fprintf(' %.3f', stats.(f)) ;
  end
  fprintf('\n') ;
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulate_gradients(state, net, opts, batchSize, mmap)
% -------------------------------------------------------------------------
for p=1:numel(net.params)

  % bring in gradients from other GPUs if any
  if ~isempty(mmap)
    numGpus = numel(mmap.Data) ;
    tmp = zeros(size(mmap.Data(labindex).(net.params(p).name)), 'single') ;
    for g = setdiff(1:numGpus, labindex)
      tmp = tmp + mmap.Data(g).(net.params(p).name) ;
    end
    net.params(p).der = net.params(p).der + tmp ;
  else
    numGpus = 1 ;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = ...
          (1 - thisLR) * net.params(p).value + ...
          (thisLR/batchSize) * net.params(p).der ;

    case 'gradient'
      thisDecay = opts.weightDecay * net.params(p).weightDecay ;
      thisLR = state.learningRate * net.params(p).learningRate ;
      state.momentum{p} = opts.momentum * state.momentum{p} ...
        - thisDecay * net.params(p).value ...
        - (1 / batchSize) * net.params(p).der ;
      net.params(p).value = net.params(p).value + thisLR * state.momentum{p} ;

    case 'otherwise'
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function mmap = map_gradients(fname, net, numGpus)
% -------------------------------------------------------------------------
format = {} ;
for i=1:numel(net.params)
  format(end+1,1:3) = {'single', size(net.params(i).value), net.params(i).name} ;
end
format(end+1,1:3) = {'double', [3 1], 'errors'} ;
if ~exist(fname) && (labindex == 1)
  f = fopen(fname,'wb') ;
  for g=1:numGpus
    for i=1:size(format,1)
      fwrite(f,zeros(format{i,2},format{i,1}),format{i,1}) ;
    end
  end
  fclose(f) ;
end
labBarrier() ;
mmap = memmapfile(fname, 'Format', format, 'Repeat', numGpus, 'Writable', true) ;

% -------------------------------------------------------------------------
function write_gradients(mmap, net)
% -------------------------------------------------------------------------
for i=1:numel(net.params)
  mmap.Data(labindex).(net.params(i).name) = gather(net.params(i).der) ;
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

stats = struct() ;

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;

      if g == 1
        stats.(s).(f) = 0 ;
      end
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
stats = struct() ;
for i = 1:numel(sel)
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net, stats)
% -------------------------------------------------------------------------
net_ = net ;
net = net_.saveobj() ;
save(fileName, 'net', 'stats') ;

% -------------------------------------------------------------------------
function [net, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function yout2 = ola(yin,fragsize,overlap)

noover  = 1 - overlap;
Move    = fragsize * noover;
Nframe  = size(yin,1);
NFFT    = size(yin,2);
yout    = zeros(1, Nframe*Move+NFFT);

for n = 1:Nframe
    framestart  = 1 + (n-1)*Move;
    frameend    = fragsize + (n-1)*Move;
    yout(framestart:frameend) = yout(framestart:frameend) + yin(n,1:fragsize);
end

yout2 = yout(1:Nframe*Move+fragsize*overlap);

author: XiaoxiangGao


text:
function sedumi_binary_error()
%Throw an error indicating the sedumi binaries are not available

if exist('OCTAVE_VERSION','builtin'),
    software = 'Octave';
else
    software = 'Matlab';
end

fs = filesep;
mpath = mfilename('fullpath');
temp = strfind( mpath, fs );
mpath = mpath( 1 : temp(end) - 1 );

error( 'SeDuMi:NoBinaries', ...
    [ 'The SeDuMi binaries for your platform are not installed.\n', ...
      'To rectify, run the following commands from the %s command line:\n\n', ...
      '   cd %s\n', ...
      '   install_sedumi\n\n', ...
      'For more information, read %s%sInstall.txt.' ], ...
    software, mpath, mpath, fs );

author: yu-jiang

repo: radpbook

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


text:
%QUESTION NO:4

%For the following 2-class problem determine the decision boundaries
%obtained by LMS and perceptron learning laws.
% Class C1 : [-2 2]', [-2 3]', [-1 1]', [-1 4]', [0 0]', [0 1]', [0 2]',
%            [0 3]' and [1 1]'
% Class C2 : [ 1 0]', [2 1]', [3 -1]', [3 1]', [3 2]', [4 -2]', [4 1]',
%            [5 -1]' and [5 0]'

clear;
inp=[-2 -2 -1 -1 0 0 0 0 1 1 2 3 3 3 4 4 5 5;2 3 1 4 0 1 2 3 1 0 1 -1 1 2 -2 1 -1 0];
out=[1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0];
choice=input('1: Perceptron Learning Law\n2: LMS Learning Law\n Enter your choice :');
switch choice
    case 1
            network=newp([-2 5;-2 4],1);
            network=init(network);
            y=sim(network,inp);
            figure,plot(inp,out,inp,y,'o'),title('Before Training');
            axis([-10 20 -2.0 2.0]);
            network.trainParam.epochs = 20;
            network=train(network,inp,out);
            y=sim(network,inp);
            figure,plot(inp,out,inp,y,'o'),title('After Training');
            axis([-10 20 -2.0 2.0]);
            display('Final weight vector and bias values : \n');
            Weights=network.iw{1};
            Bias=network.b{1};
            Weights
            Bias
            Actual_Desired=[y' out'];
            Actual_Desired
    case 2
            network=newlin([-2 5;-2 4],1);
            network=init(network);
            y=sim(network,inp);
            network=adapt(network,inp,out);
            y=sim(network,inp);
            display('Final weight vector and bias values : \n');
            Weights=network.iw{1};
            Bias=network.b{1};
            Weights
            Bias
            Actual_Desired=[y' out'];
            Actual_Desired
    otherwise
            error('Wrong Choice');
end

author: Sable

repo: mcbench-benchmarks


text:
function [numRows numColumns area] = BoxSize(bbox)
% [numRows numColumns Surface] = BoxSize(bbox)
%
% Retrieves number of rows, columns, and surface area from bounding box
%
% bbox:         4 x N Bounding box as [rowBegin colBegin rowEnd colEnd]
%
% numRows:      Number of rows of boxes
% numColumns:   Number of columns of boxes
% area:         Area of boxes
%
%     Jasper Uijlings - 2013

% Box is empty
if isempty(bbox)
    numRows = 0;
    numColumns = 0;
    area = 0;
    return
end

numRows = bbox(:,3) - bbox(:,1) + 1;
numColumns = bbox(:,4) - bbox(:,2) + 1;
area = numRows .* numColumns;

author: zouchuhang

repo: LayoutNet


text:
report_this_filefun(mfilename('fullpath'));

niv = 0.6:0.2:4;
bz2 = [];

tdiff = round((teb - t0b)*365/par1);
wai = waitbar(0,' Please Wait ...  ');
set(wai,'NumberTitle','off','Name','Makegrid  -Percent done');;
ni = str2double(prmptdlg('Number of events in each window?','100'));
na = str2double(prmptdlg('Number of random samples drawn ?','30'));
for iwl = 0.6:0.2:4
    iwl = iwl*365/par1
    zr = [];
    for i = 1:na;
        l = ceil(rand([ni 1])*a.Count);
        [cumu, xt] = hist(a(l,3),(t0b:par1/365:teb));
        for j = 2:tdiff-iwl
            cu = [cumu(1:j-1) cumu(j+iwl+1:length(cumu))];
            mean1 = mean(cu);
            mean2 = mean(cumu(j:j+iwl));
            var1 = cov(cu);
            var2 = cov(cumu(j:j+iwl));
            as(j) = (mean1 - mean2)/(sqrt(var1/(length(cumu)-iwl)+var2/iwl));
        end     % for j
        zr = [zr as];
    end
    bz2 = [bz2 ; zr];
    [len,len2] = size(bz2);
    waitbar(len/length(niv))
end
close(wai)
figure
pl =plot(niv,prctile2(bz2',50),'b');
set(pl,'LineWidth',2.0)
hold on
pl=plot(niv,prctile2(bz2',99),'b');
set(pl,'LineWidth',2.0)
pl=plot(niv,max(bz2'),'b-.');
set(pl,'LineWidth',2.0)
pl=plot(niv,prctile2(bz2',1),'b');
set(pl,'LineWidth',2.0)
pl=plot(niv,min(bz2'),'b-.');
set(pl,'LineWidth',2.0)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
grid
xlabel('Windowlength in [years]')
ylabel('Range of z')
title(['ni  =  ' num2str(ni) 'events, ' num2str(na) ' random samples'])

matdraw

author: CelsoReyes


text:
classdef prtFeatSel < prtAction
    % prtFeatSel

    properties (SetAccess = protected)

        isSupervised = true   %True
        isCrossValidateValid = true; % False
    end

    methods
        function obj = prtFeatSel()
            % As an action subclass we must set the properties to reflect
            % our dataset requirements
            obj.classTrain = 'prtDataSetStandard';
            obj.classRun = 'prtDataSetStandard';
            obj.classRunRetained = true;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function imgExt = getImgExt(seqFolder)
% return image extension including the dot
% e.g.   '.jpg'

imgExt = '.jpg';

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
% function Figure2A()

%% Stragegy implementation with Historical closing odds
% Comparison of the returns of our strategy vs a random bet strategy
% at clossing odds

% comment the next 3 lines if using Matlab
% warning("off")
% pkg load statistics
% pkg load nan

% clear all
dbstop if error
addpath('./aux_files/')
addpath('./aux_plot/')
addpath('./strategies/')

%% Parameters
dat_dir = '../data/';
file_name = 'closing_odds.csv';
bet = 50; % money on each bet
marg = 0.05; % margin odds above the mean.
nSamps = 2000; % number of returns to calculate (with replacement) for the random strategy
rand('seed',1) % use always the same seed to get same results
runStrategies = 1; % 1: run both strategies, 0: load results from disk

%% Run strategies

fid = fopen([dat_dir file_name], 'r');
% 1. match_table_id: unique identifier of the game
% 2. league of the game
% 3. match date
% 4. home team
% 5. 90-minute score of home team
% 6. away team
% 7. 90-minute score of away team
% 8. average closing odds home win
% 9. average closing odds draw
% 10. average closing odds away win
% 11. maximum offered closing odds home win
% 12. maximum offered closing odds draw
% 13. maximum offered closing odds away win
% 14. name of bookmaker offering maximum closing odds for home win
% 15. name of bookmaker offering maximum closing odds for draw
% 16. name of bookmaker offering maximum closing odds for away win
% 17. number of available closing odds for home win
% 18. number of available closing odds for draw
% 19. number of available closing odds for away win
C = textscan(fid, '%s %s %s %s %f %s %f %f %f %f %f %f %f %s %s %s %f %f %f', 'delimiter', ',');
fclose(fid);

dat = [C{5} C{7} C{8} C{9} C{10} C{11} C{12} C{13} C{17} C{18} C{19}];

if runStrategies

    %% Implement Our Strategy
    s1 = beatTheBookie(dat, bet, marg);
    s1.name = 'BeatTheBookies';
    fprintf('Finished running "beatTheBookie"\n')

    % Proportion of Home, Draw or away for games selected by our strategy.
    s1.pHome = mean(s1.ids == 1);
    s1.pDraw = mean(s1.ids == 2);
    s1.pAway = mean(s1.ids == 3);

    %% Implement Random bet strategy
    nGamesStrategy = length(s1.money) -1; % number of games that were selected for the historical betting
    s2 = randomBetStrategy(dat, nSamps, nGamesStrategy, bet , s1);
    s2.name = 'RandomStrategy';
    fprintf('Finished running "random Strategy"\n')

    save([dat_dir 'returns_HistoricalClosingOdds'], 's1', 's2', 'bet')

else
    %% Or load pre-calculated results from disk
    load([dat_dir 'returns_HistoricalClosingOdds.mat']);

end

%% Mean closing odds and Expected accuracy

% Compute descriptive stats
mS1 = mean(s1.mean_odds);
mS2 = mean(s2.mean_odds(:));
stdS1 = std(s1.mean_odds);
stdS2 = std(s2.mean_odds(:));

% These are are the intercepts obtained in the regression analysis of
% Figure 1 (see Figure1.m)
offsets = [-0.034, -0.057, -0.037];

% Calculate Expected Accuracy of our strategy
s1_prob = [ ( 1 ./ s1.mean_odds(s1.ids==1 )) + offsets(1) , (1 ./ s1.mean_odds(s1.ids==2)) + offsets(2), ...
    (1 ./ s1.mean_odds(s1.ids==3)) + offsets(3)];
s1_accuracy = mean(s1.accuracy);
s1_expectedAccuracy = mean(s1_prob);

% Calculate Expected Accuracy of Random bet Strategy
for m = 1 : size(s2.mean_odds, 1)

    oddsHome = s2.mean_odds(m, s2.ids(m,:)==1);
    oddsDraw = s2.mean_odds(m, s2.ids(m,:)==2);
    oddsAway = s2.mean_odds(m, s2.ids(m,:)==3);
    s2_prob = [ ( 1 ./ oddsHome) + offsets(1), (1 ./ oddsDraw) + offsets(2), ...
            (1 ./ oddsAway) + offsets(3)];
    s2_expectedAccuracies(m) = mean(s2_prob);

end
s2_accuracy = mean(mean(s2.accuracy,2));
s2_expectedAccuracy = mean(s2_expectedAccuracies);

randomStrategyMean = nanmean(s2.money(:,end));
randomStrategyStd = nanstd(s2.money(:,end));

delta_sigma = (s1.money(end) - randomStrategyMean) / randomStrategyStd; % distance to the mean in standard deviations

p = normcdf(s1.money(end),randomStrategyMean,randomStrategyStd);
% percentage of z values expected to lie above zσ.  CI = (−zσ, zσ)
prop = (1 - p);
fraction = 1 / prop; % expressed as fraction

clc
fprintf('Mean odds of our strategy: %2.3f (STD=%2.3f) \nMean Odds Random Bet Strategy: %2.3f (STD= %2.3f) \n', ...
    mS1, stdS1, mS2, stdS2);

fprintf('Beat The Bookie statistics:\n');
fprintf('# of bets: %2.0f \n Return: %2.4f\n Profit: %2.0f\n Expected Accuracy: %2.1f\n Accuracy: %2.2f \n',length(s1.money)-1, ...
    s1.money(end)/((length(s1.money)-1) * bet),s1.money(end), s1_expectedAccuracy * 100, s1_accuracy * 100);

fprintf('Random bet strategy statistics:\n');
fprintf('# of bets: %2.0f \n Return: %2.4f\n Profit: %2.0f\n STD: %2.4f\n Expected Accuracy: %2.1f\n Accuracy: %2.2f \n',length(s2.money), ...
    randomStrategyMean/((length(s2.money)-1)*bet), randomStrategyMean, randomStrategyStd, s2_expectedAccuracy * 100, s2_accuracy * 100);

%% Figure 2A: Compare "Beat the bookie" with the Random Bet Strategy
f1 = figure(1); clf;
set(gcf, 'Position', [0 0 1200 800], 'InvertHardCopy', 'on', 'PaperPositionMode', 'auto')
hold on

% Random strategy
p3 = plot(mean(s2.money), 'r', 'LineWidth', 3);
p1 = plot(s2.money', '-r', 'LineWidth', 3);
for m = 1 : length(p1)
    p1(m).Color(4) = 0.01;
end

% Beat the bookie
p2= plot(s1.money, 'b', 'LineWidth', 3);
p2.Color(4) = 0.8;

xlabel('Game Number')
ylabel('Returns [U$D]')
fontSize = 16;
set(gca, 'FontSize', fontSize)
legend([p2, p3], 'Our Strategy', 'Random Bet Strategy', 'Location', 'SouthWest')
legend boxoff

% Change color of line back to black without affecting the legend
p3 = plot(mean(s2.money), 'k', 'LineWidth', 3);
p3.Color(4) = 0.7;

set(gca, 'YTick', -200000:50000:150000, 'YTickLabel', {-200000 -150000 -100000 -50000 0 50000 100000 150000})
set(gca, 'XTick', 0:10000:60000, 'XTickLabel', {0 10000 20000 30000 40000 50000 60000})

xlim([0 80000])
ylim([-180000 125000])
set(gca, 'FontSize', fontSize)

% Draw curly brace
drawbrace([length(s1.money)-1 s1.money(end)], [length(s1.money)-1 randomStrategyMean], 20, 'Color', 'k', 'LineWidth', 2);

tit = sprintf('%2.2f', delta_sigma);
ht = text(62000, 15000, [tit ' \sigma']);
set(ht,'Rotation',270)
set(ht,'FontSize',20)

% Draw histogram inset
axes('position', [0.71 0.20 0.22 0.3]);

final_returns = s2.money(:,end);
[counts,bins] = hist(final_returns, 30); %# get counts and bin locations
h = barh(bins,counts);
h.FaceColor = [1.0 0 0];
set(gca,'visible','off');

print(f1, '-dpng', '../figures/Figure2A.png')
print(f1, '-depsc', '../figures/Figure2A.eps')

repo: BeatTheBookie


text:
function s=get_name(a)
s=[get_name(a.algorithm)];
eval_name

n=a.func; if length(n)>6 n=[n(1:6) '..']; end;
s=[s  '(' n ')'];

if ~isempty(a.p1) s=[s ' p1=' dis(a.p1)]; end
if ~isempty(a.p2) s=[s ' p2=' dis(a.p2)]; end
if ~isempty(a.p3) s=[s ' p3=' dis(a.p3)]; end
if ~isempty(a.p4) s=[s ' p4=' dis(a.p4)]; end

function s=dis(s)

  if isa(s,'double')
    s=num2str(s);
  else
    if isa(s,'algorithm')
      s=get_name(s);
    end
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [nlines] = dicomrt_nASCIIlines(filename)
% dicomrt_nASCIIlines(filename)
%
% Returns the number of lines in an ASCII file.
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

fid=fopen(filename);
nlines = 0;
while 1
    if feof(fid)
        break;
    end
    nlines = nlines+1;
    fgetl(fid);
end
fclose(fid);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
clear, clc;

% This is an example for running the function LogisticC
%
%  Problem:
%
%  min  f(x,c) = - weight_i * log (p_i) + 1/2 * rsL2 * ||x||_2^2
%  s.t. \|x\|_1 <=z
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  weight_i denotes the weight for the i-th sample
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
% [3]  Jun Liu, Jianhui Chen, and Jieping Ye,
%      Large-Scale Sparse Logistic Regression, KDD, 2009.
%
%% ------------   History --------------------
%
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1;

m=1000;  n=1000;    % The data matrix is of size m x n

randNum=1;          % a random number

% ---------------------- generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);         % the data matrix

y=[ones(n/2,1);...
    -ones(n/2, 1)];  % the response

z=40;                % the radius of the L1 ball

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=40;    % maximum number of iterations

% Normalization
opts.nFlag=0;       % without normalization

% Regularization
opts.rsL2=0;        % the squared two norm term

% Group Property
opts.sWeight=[1,1]; % set the weight for positive and negative samples

%----------------------- Run the code LogisticC -----------------------
fprintf('\n lFlag=0 \n');
opts.lFlag=0;       % Nemirovski's line search
tic;
[x1, c1,funVal1, ValueL1]= LogisticC(A, y, z, opts);
toc;

fprintf('\n lFlag=1 \n');
opts.lFlag=1;       % adaptive line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x2, c2, funVal2, ValueL2]= LogisticC(A, y, z, opts);
toc;

figure;
plot(funVal1,'-r');
hold on;
plot(funVal2,'--b');
legend('lFlag=0', 'lFlag=1');
xlabel('Iteration (i)');
ylabel('The objective function value');

% --------------------- compute the pathwise solutions ----------------
opts.fName='LogisticC';      % set the function name to 'LogisticC'
Z=[10, 20, 30, 40];          % set the parameters

% run the function pathSolutionLogistic
fprintf('\n Compute the pathwise solutions, please wait...');
[X,C]=pathSolutionLogistic(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
function rsamobj = medfilt1(rsamobj, nsamples_to_average_over)
% MEDFILT1 Median filter for RSAM objects
% For removing outliers due to transient events, or telemetry spikes
    for c = 1:numel(rsamobj)
        if nargin==1
            nsamples_to_average_over = 3;
        end
        rsamobj(c).data = medfilt1(rsamobj(c).data, nsamples_to_average_over);
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function NR = esvm_show_transfer_figure(I, models, topboxes, overlays, ...
                                        current_rank, corr)
% Show a figure with the detections of the exemplar svm model
% NOTE(TJM): this function needs cleanup, but works in the pipeline
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

topboxes = topboxes(1:min(2,size(topboxes,1)),:);

%if enabled, try to transfer objects too
%add_one = 1;
add_one = 0;

figure(1)
clf
if add_one == 1
  ha = tight_subplot(2, 2, .05, ...
                     .1, .01);
  NR = [2 2];
else
  ha = tight_subplot(1, 2, .05, ...
                     .01, .01);
  %hb = tight_subplot(1, 2, .05, ...
  %                   .1, .01);
  NR = [2 2];

end

PADDER = 100;
Ipad = pad_image(I,PADDER);

%Cut out detection chunks
for q = size(topboxes,1):-1:1
  if q>size(topboxes,1)
    continue
  end

  bb = topboxes(q,1:4);

  bb = bb + PADDER;
  r1 = [round(bb(2)) round(bb(4))];
  r2 = [round(bb(1)) round(bb(3))];
  r1 = cap_range(r1,1,size(Ipad,1));
  r2 = cap_range(r2,1,size(Ipad,2));

  mask = zeros(size(Ipad,1),size(Ipad,2));
  mask(r1(1):r1(2), r2(1):r2(2))=1;
  mask = logical(mask);
  chunks{q} = maskCropper(Ipad, mask);
end

%%NOTE: below is a loop, but we only process the first element
%since we only show the top association
N = min(1,size(topboxes,1));
for i = 1:N

  mid = topboxes(i,6);
  flip = topboxes(i,7);
  hogpic = HOGpicture(models{topboxes(i,6)}.model.w);
  hogpic = jettify(hogpic);
  if flip == 1
    hogpic = flip_image(hogpic);
  end

  Iex = esvm_get_exemplar_icon(models, mid, topboxes(i,7));

  Iex1 = imresize(chunks{i},[size(Iex,1) size(Iex,2)]);
  Iex1 = max(0.0,min(1.0,Iex1));

  Iex1 = imresize(hogpic,[size(Iex1,1) size(Iex,2)], 'nearest');
  Iex1 = max(0.0,min(1.0,Iex1));

  PPP = round(size(Iex,1)*.2);
  Iex2 = Iex1(1:PPP,:,:)*0+1;
  Ishow = cat(1,Iex1,Iex2,Iex);

  SHOW_SEG = 1;
  if SHOW_SEG == 1

    Ia = (overlays{1}.mini_overlay.I.* ...
              repmat(overlays{1}.mini_overlay.alphamask,[1 1 3]));
    Ib = Iex;
    alphamap = overlays{1}.mini_overlay.alphamask;
    alphamap(alphamap>0) = .5;
    alphamap(alphamap==0) = .7;
    alphamap = repmat(alphamap,[1 1 3]);
    Itotal = Ia.*alphamap + Ib.*(1-alphamap);

    Imeta = Itotal;
    offset = [0 size(Ishow,1)+size(Iex2,1)];

    %Ishow = cat(1,Ishow,Iex2,Imeta);

  elseif strcmp(models{mid}.cls,'bus')
    m1 = repmat(double(overlays{1}.mini_overlay.alphamask>0)*.8,[1 ...
                    1 3]);
    m2 = repmat(double(overlays{1}.mini_overlay.alphamask>0)*.2,[1 ...
                    1 3]);
    m2(m2(:)==0) = 1;

    Imeta = overlays{1}.mini_overlay.I.*m1 +...
            Iex.*m2;
    offset = [0 size(Ishow,1)+size(Iex2,1)]
    overlays{1}.mini_overlay.faces = ...
        cellfun2(@ ...
                 (x)bsxfun(@plus,x,offset),overlays{1}.mini_overlay.faces);

    %Ishow = cat(1,Ishow,Iex2,Imeta);
  end

  bb1 = [1 1 size(Iex,2) size(Iex,1)];
  bb2 = bb1;

  bb2([2 4]) = bb2([2 4]) + size(Iex,1) + PPP;
  axes(ha(1));
  imagesc(Ishow)

  axis image
  axis off

  bouter = [1 1 size(Ishow,2) size(Ishow,1)];
  bouter(1) = bouter(1) - round(size(Ishow,2)*.05);
  bouter(2) = bouter(2) - round(size(Ishow,1)*.05);
  bouter(3) = bouter(3) + round(size(Ishow,2)*.05);
  bouter(4) = bouter(4) + round(size(Ishow,1)*.05);
  plot_bbox(bouter, '', [0 0 0], [0 0 0], 0, [2 1])

  if isfield(overlays{1}.mini_overlay,'faces')
    plot_faces(overlays{1}.mini_overlay.faces);
  end

  axis image
  curcolor = [0 1 0];
  if corr == 0
    curcolor = [1 0 0 ];
  end

  plot_bbox(bb1, sprintf('Exemplar-SVM %s %d',models{mid}.cls,mid), ...
            curcolor, curcolor, 0, [2 1])

  %sprintf('Exemplar.%d',mid)

  plot_bbox(bb2, sprintf('Exemplar Image %d',mid), [1 1 0], [1 1 0], 0, [2 1])

  id_string = '';
  if isfield(models{abs(mid)},'curid')
    try
      oid = models{abs(mid)}.objectid;
    catch
      oid = -1;
    end
    id_string = sprintf('%s.%d',models{abs(mid)}.curid, ...
                        oid);
  else
    id_string = sprintf('%s@%s',models{abs(mid)}.models_name,models{abs(mid)}.cls);
  end
  lrstring='';
  if flip == 1
    lrstring = '.LR';
  end

  sprintf('Exemplar: %s%s',id_string,lrstring);

end

extraI = overlays{1}.I;

clipped_top = clip_to_image(topboxes(1,:),[1 1 size(extraI,2) ...
                    size(extraI,1)]);
clipped_top([1 2]) = clipped_top([1 2])+2;
clipped_top([3 4]) = clipped_top([3 4])-2;

axes(ha(2));

PPP = round(size(I,1)*.05);
I2 = I(1:PPP,:,:)*0+1;
Ishow = cat(1,I,I2,extraI);
%Ishow = extraI;

imagesc(Ishow)
axis image
axis off
%title(sprintf('Detection (Rank=%d, Score=%.3f)',current_rank, ...
%              topboxes(1,end)))

curcolor = [0 1 0];
if corr == 0
  curcolor = [1 0 0 ];
end

plot_bbox(clipped_top,sprintf('%s: %.3f',models{1}.cls,clipped_top(end)),...
          curcolor, curcolor,0,[2 1]);

if 0
  %%NOTE: not sure why this is not working?
%% find xform from exemplar onto new detection
start_box = [1 1 size(Iex,2) size(Iex,1)];
end_box = topboxes(1,1:4);
xform = find_xform(start_box,end_box);

faces2 = cellfun2(@(x)apply_xform(x,xform), ...
                  overlays{1}.mini_overlay.faces);

offset = [0 size(I,1)+size(I2,1)];
offset = [0 0];
faces2 = cellfun2(@ ...
                  (x)bsxfun(@plus,x,offset),faces2);
end
%plot_faces(faces2);

%axes(ha(3));
%imagesc(extraI);

%sprintf('Exemplar.%d',mid)
%orig green
clipped_top([2 4]) = clipped_top([2 4]) + size(I,1)+size(I2,1);

%dont do yellow box for segmentation
if 1 %%~strcmp(models{mid}.cls,'bus')
  plot_bbox(clipped_top,'Exemplar',[1 1 0],[1 1 0],0,[2 1]);
end
axis image
axis off
%title('Appearance Transfer')%sprintf('Exemplar Inpainting: %.3f',topboxes(1,end)))
%title(sprintf('Exemplar Inpainting: %.3f',topboxes(1,end)))

title('Exemplar-SVM detection','FontSize',12);

if add_one == 1
  %transfer objects here

  is_complement = sum(ismember(models{1}.cls,{'bicycle', ...
                    'motorbike','horse'}));

  if isfield(exemplar_overlay,'friendbb') & ...
        size(exemplar_overlay.friendbb,1)>0

    axes(ha(4));

    imagesc(I)
    axis image
    axis off

    c = jet(20);

    for q = 1:size(exemplar_overlay.friendbb,1)
      [aa,bb] = find(ismember(VOCopts.classes, ...
                              exemplar_overlay.friendclass{q}));
      curcolor = c(aa,:);
      %Yellow color always
      curbb = exemplar_overlay.friendbb(q,:);
      curbb = clip_to_image(curbb,[1 1 size(I,2) size(I,1)]);

      plot_bbox(curbb,...
                exemplar_overlay.friendclass{q},curcolor,curcolor,1,[2 ...
                    1]);

      [aa,bb] = find(ismember(VOCopts.classes, ...
                              models{1}.cls));
      curcolor = c(aa,:);
      curcolor = [1 1 0];
      curcolor = [0 1 0];
      if corr == 0
        curcolor = [1 0 0 ];
      end

      curbb = clip_to_image(topboxes(1,:),[1 1 size(I,2) size(I, ...
                                                  1)]);

      plot_bbox(curbb,models{1}.cls,curcolor,curcolor,0,[2 1]);
      title(sprintf('Detection: Predicted %s', ...
                    exemplar_overlay.friendclass{1}));

    end
  elseif ~is_complement && isfield(exemplar_overlay,'overlay2')
    %% SHOW TRANSFERRED BUSES
    axes(ha(4));
    imagesc(exemplar_overlay ...
            .overlay2.segI);
    title(exemplar_overlay.overlay2.title);
    axis image
    axis off

    %Faces are for transferring geometry
    if isfield(exemplar_overlay.overlay2,'faces')
      faces = exemplar_overlay.overlay2.faces;
      for q = 1:length(faces)
        if length(faces{q})>0
          f = faces{q};
          %%estimate Xform here:
          source = topboxes(1,:);
          target = models{source(6)}.gt_box;

          if topboxes(1,7) == 1
            target = flip_box(target,models{source(6)}.sizeI);
            f = flip_faces(f,models{source(6)}.sizeI);

            fprintf(1,'got here for a flip')
            %source = flip_box(source,size(I));
          end

          %%BUG: we should be using the new find_xform function...
          clear xs ys
          xs(:,1) = target([1 2])';
          xs(:,2) = target([3 2])';
          xs(:,3) = target([1 4])';
          xs(:,4) = target([3 4])';

          ys(:,1) = source([1 2])';
          ys(:,2) = source([3 2])';
          ys(:,3) = source([1 4])';
          ys(:,4) = source([3 4])';

          xs(3,:) = 1;
          ys(3,:) = 1;

          A=ys*pinv(xs);
          A(abs(A)<.000001) = 0;

          f2 = f';
          f2(3,:) = 1;

          f2 = A*f2;
          f2(3,:) = [];
          f = f2';

          hold on;
          plot([f(:,1); f(1,1)],[f(:,2); f(1,2)],'k','LineWidth', ...
               4);

        end
      end
    end

  else
    %% SHOW FRIENDS
    axes(ha(4))
    imagesc(I)
    axis image
    axis off

    titler = sprintf('Detection');
    if isfield(models{mid},'friendclass')
      titler = [titler sprintf(': No %s',models{mid}.friendclass)];
    end
    title(titler);

    [aa,bb] = find(ismember(VOCopts.classes, ...
                            models{1}.cls));
    c = jet(20);
    curcolor = c(aa,:);
    curcolor = [1 1 0];

    clipped_top = clip_to_image(topboxes(1,:),[1 1 size(extraI,2) ...
                    size(extraI,1)]);
    clipped_top([1 2]) = clipped_top([1 2])+2;
    clipped_top([3 4]) = clipped_top([3 4])-2;
    plot_bbox(clipped_top,models{1}.cls,curcolor,curcolor,0,[2 1]);
  end
end

author: quantombone

repo: exemplarsvm


text:
%--- help for ts/describe ---
%
%  Print the description of the time series object
%
%  Args:
%
%     db (ts object): time series to describe
%
%  Note:
%
%     The function prints
%
%        - mean
%        - standard deviation
%        - min
%        - 25th percentile
%        - 50th percentile
%        - 75th percentile
%        - max
%        - variable names
%        - variable descriptions
%
%

author: jmaih


text:
function E = keplerEq(M,e,eps)
% Function solves Kepler's equation M = E-e*sin(E)
% Input - Mean anomaly M [rad] , Eccentricity e and Epsilon
% Output  eccentric anomaly E [rad].
   	En  = M;
	Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
	while ( abs(Ens-En) > eps )
		En = Ens;
		Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
	end;
	E = Ens;


text:
function [models] = featureSelection(X,Y,maxOrder,nBoot,Info,imbalance,seed)
% -------------------------------------------------------------------------
% function [models] = featureSelection(X,Y,maxOrder,nBoot,Info,imbalance,seed)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection according to the 0.632+
% bootstrap methodology for an input matrix of features and and input
% outcome vector, and for multiple model orders as defined by the user.
% See ref. [1] for more details. This function uses logistic regression
% utilities from DREES <http://www.cerr.info/drees>, and a fast
% implementation of AUC calculation by Enric Junqué de Fortuny that is
% available at: <http://www.mathworks.com/matlabcentral/fileexchange/41258-faster-roc-auc>
%
% NOTE: This function now maximizes 0.5*AUC + 0.5*(1-abs(SENSITIVITY-SPECIFICITY))
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - X: Matrix of size [nInst X nFeat], specifying the numerical data of the
%      features of the input features, where 'nInst' refers to the number
%      of instances in X, and 'nFeat' to the number of features in X.
%      Each column is a different feature.
% - Y: Column vector of size [nInst X 1] specifying the outcome status
%      (1 or 0) for all instances.
% - maxOrder: Integer specifying the maximal model order to construct.
% - nBoot: Number of bootstrap samples to use.
% - Info: Cell of size [nFeat X 1] of strings specifying the name of each
%         feature in 'X'.
% - imbalance: String specifying the type of imbalance-adjustement strategy
%              employed. Either 'IABR' for imbalance-adjusted bootstrap
%              resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%              logistic regression (formal reference to come).
% - seed: (optional input). Numerical number to use as random generator
%         seed for bootstrapping experiments.
%         --> Ex: 54288
% -------------------------------------------------------------------------
% OUTPUTS:
% - models: Structure specifying the resulting multivariable models for the
%           input feature set in 'data', for each model order. Example for
%           order 4:
%        --> models.Order4.Data: Matrix of size [nInst X order], specifying
%                                the selected features, in order of
%                                selection.
%        --> models.Order4.Name: Cell specifying the names of the selected
%                                features, in order of selection.
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Martin Vallieres <mart.vallieres@gmail.com>
% - DREES development team <http://www.cerr.info/drees> (logistic regression)
% - Enric Junqué de Fortuny (fastAUC.cpp)
% -------------------------------------------------------------------------
% HISTORY:
% - Creation - May 2015
% - Revision I - July 2015: including imbalance-adjusted logistic regression
% - Revision II - July 2015: maximizing 0.5*AUC + 0.5*(1-abs(SENSITIVITY-SPECIFICITY))
% - Revision III - December 2016: Including initial seed as input (optional)
%                                 for reproducibility of bootstrapping experiments
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2016  Martin Vallieres
% --> Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
%
%    _______________________________________________________________
%
% --> Copyright (c) 2013, Enric Junqué de Fortuny
%     All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% 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 OWNER 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.
% -------------------------------------------------------------------------

% RANDOM NUMBER GENERATOR SEED
if nargin == 7
    rng(seed)
else
    if ~RandStream.getGlobalStream.Seed
        rng('shuffle')
    end
end

% INITIALIZATION
top = 1-1/exp(1);
low = 1/exp(1);
nFeat = size(X,2);
models = struct; % Final model structure

% GETTING BOOTSTRAP SAMPLES FOR ALL EXPERIMENTS + LOGISTIC REGRESSION TYPE
if strcmp(imbalance,'IABR')
    logisticRegression = @(x,y) applyStandardLR(x,y);
    [trainSets,testSets] = buildBootSet(Y,nBoot,'adjust'); % 'trainSets' is a matrix, 'testSets' is a cell
elseif strcmp(imbalance,'IALR')
    logisticRegression = @(x,y) applyEnsembleLR(x,y);
    [trainSets,testSets] = buildBootSet(Y,nBoot); % 'trainSets' is a matrix, 'testSets' is a cell
end

% FORWARD FEATURE SELECTION (for all different starters)
modelMat = zeros(nFeat,maxOrder);
metricMat = zeros(nFeat,maxOrder);
for i = 1:nFeat
     indLeft = 1:nFeat;
     % Order 1
     indLeft(i) = [];
     modelMat(i,1) = i;
     Xtrain = X(:,i);
     [Xtrain] = normalizeZeroOne(Xtrain);
     [coeff] = logisticRegression(Xtrain,Y);
     [resp] = responseLR(Xtrain,coeff);
     [aucData,sensData,specData,~] = calcPerformMetrics(resp,Y,0);
     aucTemp = 0; sensTemp = 0; specTemp = 0;
     for n = 1:nBoot
         Xtrain = X(trainSets(:,n),i);Xtest = X(testSets{n},i); Ytrain = Y(trainSets(:,n),1); Ytest = Y(testSets{n},1);
         [Xtrain,Xtest] = normalizeZeroOne(Xtrain,Xtest);
         [coeff] = logisticRegression(Xtrain,Ytrain);
         [resp] = responseLR(Xtest,coeff);
         [aucBoot,sensBoot,specBoot,~] = calcPerformMetrics(resp,Ytest,0);

         % FOR AUC
         alpha = top/(1-low*(aucData-aucBoot)/(aucData-0.5+eps));
         if alpha > 1
             alpha = 1;
         elseif alpha < top
             alpha = top;
         end
         if aucBoot < 0.5
             aucBoot = 0.5;
         end
         aucTemp = aucTemp + (1-alpha)*aucData + alpha*aucBoot;

         % For SENSITIVITY
         alpha = top/(1-low*(sensData-sensBoot)/(sensData+eps));
         if alpha < top
             alpha = top;
         end
         sensTemp = sensTemp + (1-alpha)*sensData+alpha*sensBoot;

         % For SPECIFICITY
         alpha = top/(1-low*(specData-specBoot)/(specData+eps));
         if alpha < top
             alpha = top;
         end
         specTemp = specTemp + (1-alpha)*specData+alpha*specBoot;

     end
     aucTemp = aucTemp/nBoot; sensTemp = sensTemp/nBoot; specTemp = specTemp/nBoot;
     metricMat(i,1) = 0.5*aucTemp + 0.5*(1-abs(sensTemp-specTemp));

     % Going for orders 2 to maxOrder
     for j = 2:maxOrder
         maxMetric = 0;
         for k = 1:(nFeat-j+1)
             indexModel = [modelMat(i,1:(j-1)),indLeft(k)];
             Xtrain = X(:,indexModel);
             [Xtrain] = normalizeZeroOne(Xtrain);
             [coeff] = logisticRegression(Xtrain,Y);
             [resp] = responseLR(Xtrain,coeff);
             [aucData,sensData,specData,~] = calcPerformMetrics(resp,Y,0);
             aucTemp = 0; sensTemp = 0; specTemp = 0;
             for n=1:nBoot
                 Xtrain = X(trainSets(:,n),indexModel); Xtest = X(testSets{n},indexModel); Ytrain = Y(trainSets(:,n),1); Ytest = Y(testSets{n},1);
                 [Xtrain,Xtest] = normalizeZeroOne(Xtrain,Xtest);
                 [coeff] = logisticRegression(Xtrain,Ytrain);
                 [resp] = responseLR(Xtest,coeff);
                 [aucBoot,sensBoot,specBoot,~] = calcPerformMetrics(resp,Ytest,0);

                 % FOR AUC
                 alpha = top/(1-low*(aucData-aucBoot)/(aucData-0.5+eps));
                 if alpha > 1
                     alpha = 1;
                 elseif alpha < top
                     alpha = top;
                 end
                 if aucBoot < 0.5
                     aucBoot = 0.5;
                 end
                 aucTemp = aucTemp + (1-alpha)*aucData + alpha*aucBoot;

                 % For SENSITIVITY
                 alpha = top/(1-low*(sensData-sensBoot)/(sensData+eps));
                 if alpha < top
                     alpha = top;
                 end
                 sensTemp = sensTemp + (1-alpha)*sensData+alpha*sensBoot;

                 % For SPECIFICITY
                 alpha = top/(1-low*(specData-specBoot)/(specData+eps));
                 if alpha < top
                     alpha = top;
                 end
                 specTemp = specTemp + (1-alpha)*specData+alpha*specBoot;

             end
             aucTemp = aucTemp/nBoot; sensTemp = sensTemp/nBoot; specTemp = specTemp/nBoot;
             metricTemp = 0.5*aucTemp + 0.5*(1-abs(sensTemp-specTemp));
             if metricTemp >= maxMetric
                 maxMetric = metricTemp;
                 index = indLeft(k);
             end
         end
         modelMat(i,j) = index;
         metricMat(i,j) = maxMetric;
         indLeft(find(indLeft==index)) = [];
     end
end

% OBTAINING MAXIMUM RESULTS FOR EVERY MODEL ORDER (maximum from all different starters)
[~,indMax] = max(metricMat);
for i = 1:maxOrder
    nameOrder = ['Order',num2str(i)];
    models.(nameOrder).Data = X(:,modelMat(indMax(i),1:i));
    models.(nameOrder).Name = cell(i,1);
    for j = 1:i
        models.(nameOrder).Name{j} = Info{modelMat(indMax(i),j)};
    end
end

end

author: mvallieres

repo: radiomics


text:
function sigma = dtiComputeImageNoise(dwRaw, bvals, brainMask, noiseCalcMethod)
% Compute image noise from a 4-d DWI image
%
% sigma = dtiComputeImageNoise(dwRaw, bvals, brainMask, [noiseCalcMethod])
%
% Takes in dwi data and returns sigma, the image noise
%
% dwRaw    - Nifti structure of the dwi data
% bvals    - A 1xN vector of bvalues for each volume in dwRaw
% brainMask- A 3-D matrix of binary values  denoting which voxels are
%            insided the brain. This is only necesary for noiseCalcMethod =
%            'b0'
% noiseCalcMethod- Whether to calculate the noise based on the corner of
%                  the image (default) or based on the standard deviation
%                  of the b=0 images ('b0').

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

if strcmp('corner', noiseCalcMethod)

    % According to Henkelman (1985), the expected signal variance (sigma) can
    % be computed as 1.5267 * SD of the background (thermal) noise.
    sz = size(dwRaw.data);
    x  = 10;
    y  = 10;
    z  = round(sz(3)/2);

    [x,y,z,s] = ndgrid(x-5:x+5, y-5:y:5, z-5:z+5, 1:sz(4));
    noiseInds = sub2ind(sz, x(:), y(:), z(:), s(:));
    sigma     = 1.5267 * std(double(dwRaw.data(noiseInds)));

elseif strcmp('b0', noiseCalcMethod)

    % Number of volumes in the dw dataset
    numVols   = size(dwRaw.data,4);
    % Get brainmask indices
    brainInds = find(brainMask);
    % preallocate a 2d array (with a 2nd dimension that is a singleton). The
    % first dimension is the number of volumes and the 3rd is each voxel
    % (within the brain mask).
    data      = zeros(numVols,1,length(brainInds));
    % Loop over the volumes and assign the voxels within the brain mask to data
    for ii=1:numVols
        tmp = double(dwRaw.data(:,:,:,ii));
        data(ii,1,:) = tmp(brainInds);
    end

    % Find which volumes ar b=0
    b0inds = find(bvals ==  0);
    n = length(b0inds);
    % Pull out the b=0 volumes
    dataB0 = squeeze(data(b0inds,1,:));
    % Calculate the median of the standard deviation. We do not think that
    % this needs to be rescaled. Henkelman et al. (1985) suggest that this
    % aproaches the true noise as the signal increases.
    sigma = median(std(dataB0,0,1));

    % std of a sample underestimates sigma (see http://nbviewer.ipython.org/4287207/)
    % This can be very big for small n (e.g., 20% for n=2)
    % We can compute the underestimation bias:
    bias = sigma * (1 - sqrt(2 / (n-1)) * (gamma(n / 2) / gamma((n-1) / 2)));

    % and correct for it:
    sigma = sigma + bias;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function t = fieldnames(obj)
% Fieldnames of a file-array object
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: fieldnames.m 7147 2017-08-03 14:07:01Z spm $

t = {...
    'fname'
    'dim'
    'dtype'
    'offset'
    'scl_slope'
    'scl_inter'
};

author: spm


text:
function out = equal(A,dim)
% check all elements of A to be equal

if nargin == 1
  out = isempty(A) || all(A == A(1));
  return;
end

%make dim first dimension
A = shiftdim(A,dim-1);
A = reshape(A,size(A,1),[]);

B = repmat(A(1,:),size(A,1),1);

out = all(A == B);

author: mtex-toolbox

repo: mtex

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


text:
function A = matrix(k, n)
%MATRIX  Test matrices accessed by number.
%        MATRIX(K, N) is the N-by-N instance of matrix number K in
%        a set of test matrices comprising those in MATLAB plus those
%        in the Matrix Computation Toolbox,
%        with all other parameters set to their default.
%        N.B. - Only those matrices which are full and take an arbitrary
%               dimension N are included.
%             - Some of these matrices are random.
%        MATRIX(K) is a string containing the name of the K'th matrix.
%        MATRIX(0) is the number of matrices, i.e. the upper limit for K.
%        Thus to set A to each N-by-N test matrix in turn use a loop like
%             for k=1:matrix(0)
%                 A = matrix(k, N);
%                 Aname = matrix(k);   % The name of the matrix
%             end
%        MATRIX(-1) returns the version number and date of the
%        Matrix Computation Toolbox.
%        MATRIX with no arguments lists the names and numbers of the M-files in the
%        collection.

%         References:
%         N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%            Second edition, Society for Industrial and Applied Mathematics,
%            Philadelphia, PA, 2002; sec. 20.5.

% Matrices from gallery:
matrices = [%
'cauchy  '; 'chebspec'; 'chebvand'; 'chow    ';
'circul  '; 'clement '; 'condex  ';
'cycol   '; 'dramadah'; 'fiedler ';
'forsythe'; 'frank   '; 'gearmat '; 'grcar   ';
'invhess '; 'invol   '; 'ipjfact '; 'jordbloc';
'kahan   '; 'kms     '; 'krylov  '; 'lehmer  ';
'lesp    '; 'lotkin  '; 'minij   '; 'moler   ';
'orthog  '; 'parter  '; 'pei     '; 'prolate ';
'randcolu'; 'randcorr'; 'rando   '; 'randsvd ';
'redheff '; 'riemann '; 'ris     '; 'smoke   ';
'toeppd  '; 'triw    ';];
n_gall = length(matrices);

% Other MATLAB matrices:
matrices = [matrices;
'hilb    '; 'invhilb '; 'magic   '; 'pascal  ';
'rand    '; 'randn   ';];
n_MATLAB = length(matrices);

% Matrices from Matrix Computation Toolbox:
matrices = [matrices;
'augment '; 'gfpp    '; 'magic   '; 'makejcf ';
'rschur  '; 'vand    '];
n_mats = length(matrices);

if nargin == 0

   rows = ceil(n_mats/5);
   temp = zeros(rows,5);
   temp(1:n_mats) = 1:n_mats;

   for i = 1:rows
      for j = 1:5
        if temp(i,j) == 0, continue, end
        fprintf(['%2.0f: ' sprintf('%s',matrices(temp(i,j),:)) '  '], ...
                temp(i,j))
      end
      fprintf('\n')
   end
   fprintf('Matrices 1 to %1.0f are from MATLAB\.', n_MATLAB)

elseif nargin == 1
   if k == 0
      A = length(matrices);
   elseif k > 0
      A = deblank(matrices(k,:));
   else
      % Version number and date of collection.
      A = 'Version 1.2, September 5 2002';
   end
else
   if k <= n_gall
      A = eval( ['gallery(''' deblank(matrices(k,:)) ''',n)'] );
   else
      A = eval( [deblank(matrices(k,:)) '(n)'] );
   end
end

author: opencobra

repo: cobratoolbox


text:
function qwgw_test05 ( )

%*****************************************************************************80
%
%% TEST05 tests QWGW for the generalized Laguerre weight.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2014
%
%  Author:
%
%    John Burkardt
%

%
%  The quadrature interval is [0,+oo).
%  Set the number of points.
%
  a = 0.0;
  n = 5;
%
%  Set the weight function parameter.
%
  alpha = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05:\n' );
  fprintf ( 1, '  Compute points and weights for Gauss quadrature\n' );
  fprintf ( 1, '  with the generalized Laguerre weight w(x) = x^alpha * exp(-x).\n' );
  fprintf ( 1, '  Order N = %d\n', n );
  fprintf ( 1, '  ALPHA = %g\n', alpha );
  fprintf ( 1, '  Interval = [0,+oo)\n' );
%
%  Set the recursion coefficients.
%
  aj = zeros ( n, 1 );
  bj = zeros ( n, 1 );

  for j = 1 : n
    aj(j) = alpha + 2 * j - 1;
  end

  for j = 1 : n - 1
    bj(j) = j * ( alpha + j );
  end
  bj(n) = 0.0;

  bj(1:n) = sqrt ( bj(1:n) );

  mu0 = gamma ( alpha + 1 )
%
%  Compute the points and weights.
%
  [ x, w ] = sgqf ( n, aj, bj, mu0 );

  r8vec_print ( n, x, '  Abscissas:' );
  r8vec_print ( n, w, '  Weights:' );

  return
end

author: johannesgerer

repo: jburkardt-m