text:
function [cdf] = gp_kfcv_cdf(gp,x,y,varargin)
%GP_KFCV_CDF  K-fold cross validation to predict CDF for GP model
%
%  Description
%    [cdf] = GP_KFCV_CDF(GP, X, Y, OPTIONS)
%    Performs K-fold cross-validation for a GP model given input matrix X
%    and target vector Y.
%
%    OPTIONS is optional parameter-value pair
%      z          - optional observed quantity in triplet (x_i,y_i,z_i)
%                   Some likelihoods may use this. For example, in
%                   case of Poisson likelihood we have z_i=E_i,
%                   that is, expected value for ith case.
%      yt         - optional observed yt in test points.
%                   Default option for yt is yt=y.
%      inf_method - inference method. Possible methods are
%                    'MAP'      parameters optimized to MAP (default)
%                    'MCMC'     MCMC sampling using GP_MC
%                    'IA'       integration approximation using GP_IA
%                    'fixed'    parameters are fixed, it either use MAP
%                               or integration approximation, depending if
%                               GP is a single GP structure or a GP array
%                               (for example from GP_IA)
%      optimf     - function handle for an optimization function, which is
%                   assumed to have similar input and output arguments
%                   as usual fmin*-functions. Default is @fminscg.
%      opt        - options for the inference method. If 'MAP' is used
%                   use optimset to set options for optimization.
%                   Default options for optimization are 'GradObj'
%                   is 'on', 'LargeScale' is 'off', 'Display' is 'off'
%      k          - number of folds in CV, default k=10
%      rstream    - number of a random stream to be used for
%                   permuting the data befor division. This way
%                   same permutation can be obtained for different
%                   models. Default is 1. See doc RandStream for
%                   more information.
%      trindex    - k-fold CV training indices. A cell array with k
%                   fields each containing index vector for respective
%                   training set.
%      tstindex   - k-fold CV test indices. A cell array with k
%                   fields each containing index vector for
%                   respective test set.
%      display    - defines if messages are displayed.
%                   - 'iter' displays output at each iteration
%
%    The output argument is
%
%           cdf- Predictive Cumulative Distribiton function evaluated in y
%           (default) or yt.
%
%
%  See also
%    DEMO_MODELCOMPARISON1, DEMO_MODELCOMPARISON2
%
% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010-2011 Aki Vehtari
% Copyright (c) 2012 Ernesto Ulloa

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.FunctionName = 'GP_KFCV';
ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
ip.addParamValue('inf_method', 'MAP', @(x) ...
  ismember(x,{'MAP' 'LOO' 'KFCV' 'WAIC' 'WAICV' 'WAICG' 'MCMC' 'IA' 'fixed'}))
ip.addParamValue('optimf', @fminscg, @(x) isa(x,'function_handle'))
ip.addParamValue('opt', struct(), @isstruct)
ip.addParamValue('k', 10, @(x) isreal(x) && isscalar(x) && isfinite(x) && x>0)
ip.addParamValue('rstream', 1, @(x) isreal(x) && isscalar(x) && isfinite(x) && x>0)
ip.addParamValue('trindex', [], @(x) isempty(x) || iscell(x))
ip.addParamValue('tstindex', [], @(x) isempty(x) || iscell(x))
ip.addParamValue('display', 'on', @(x) islogical(x) || ...
  ismember(x,{'iter' 'fold'}))
ip.parse(gp, x, y, varargin{:});
z=ip.Results.z;
yt=ip.Results.yt;
inf_method=ip.Results.inf_method;
optimf=ip.Results.optimf;
opt=ip.Results.opt;
k=ip.Results.k;
rstream=ip.Results.rstream;
trindex=ip.Results.trindex;
tstindex=ip.Results.tstindex;
display = ip.Results.display;
if isequal(display,'fold');display='iter';end

[n,nin] = size(x);

gp_orig = gp;

if ismember(inf_method,{'MAP' 'LOO' 'KFCV' 'WAIC' 'WAICV' 'WAICG'})
  optdefault=struct('Display','off');
  opt=optimset(optdefault,opt);
end

if (isempty(trindex) && ~isempty(tstindex)) || (~isempty(trindex) && isempty(tstindex))
  error('gp_kfcv: If you give cross-validation indexes, you need to provide both trindex and tstindex.')
end

if isempty(trindex) || isempty(tstindex)
  [trindex, tstindex] = cvit(n, k, rstream);
end

% *** note: yt must can be: a scalar or a vector of size(y)x1 or even a matrix of
% size(y)xN

if ~isempty(yt)
  if size(yt,1)==size(y,1)
    yt=yt;
  elseif size(yt,1)==1
    yt=bsxfun(@times,ones(size(y)),yt);
  else
    error('size of yt does not match y nor is it a scalar');
  end
  ytflag=1;
else
  ytflag=0;
  yt=y;
end

cvws=[];
trw=[];
% loop over the crossvalidation sets

for i=1:length(trindex)
  if isempty(tstindex{i})
    continue
  end

  if isequal(display,'iter')
    fprintf(' The CV-fold number: %d/%d \n', i, k)
  end

  % Set the training and test sets for i'th cross-validation set
  xtr = x(trindex{i},:);
  ytr = y(trindex{i},:);
  yttr= yt(trindex{i},:);

  xtst = x(tstindex{i},:);
  ytst = y(tstindex{i},:);
  yttst = yt(tstindex{i},:);

  if ~isempty(z)
    ztr = z(trindex{i},:);
    zt = z;
    flagz=1;
    %* yt = y;
    %       opt_tst.zt = z;
    %       opt_tst.yt = y;
  else
    flagz=0;
    ztr = [];
    %*yt = y;
    zt = [];
    %       opt_tr = struct();
    %       opt_tst.yt = y;
  end

  gp = gp_orig;

  if iscell(gp)
    gptype=gp{1}.type;
  else
    gptype=gp.type;
  end
  tstind2 = [];

  switch gptype
    case {'FIC' 'CS+FIC'}
      tstind2 = trindex{i};
    case 'PIC'
      % Set the block indices for the cv set of data points. Variable
      % naming(e.g tstind2) because parfor loop.
      ntr = size(xtr,1);
      ntst = size(xtst,1);
      trind2 = [];
      for i1=1:length(gp.tr_index)
        tstind2{i1} = [];
        trind2{i1} = [];
        for j1 = 1:length(gp.tr_index{i1})
          indtmp = find( sum((xtr - repmat(x(gp.tr_index{i1}(j1),:),ntr,1)).^2,2) == 0 );
          if isempty( indtmp )
            indtmp = find( sum((xtst - repmat(x(gp.tr_index{i1}(j1),:),ntst,1)).^2,2) == 0 );
            tstind2{i1} = [tstind2{i1} indtmp];
          else
            trind2{i1} = [trind2{i1} indtmp];
          end
        end
      end
      if iscell(gp)
        for j=1:numel(gp)
          gp{j}.tr_index=trind2;
        end
      else
        gp.tr_index = trind2;
      end
  end

  % Conduct inference
  switch inf_method
    case 'MAP'
      if flagz
        gp=gp_optim(gp,xtr,ytr,'z',ztr(:,size(z,2)),'opt',opt, 'optimf', optimf);
        w=gp_pak(gp);
        cvws(i,:)=w;
      else
        gp=gp_optim(gp,xtr,ytr,'z',ztr,'opt',opt, 'optimf', optimf);
        w=gp_pak(gp);
        cvws(i,:)=w;
      end
    case {'LOO' 'KFCV' 'WAIC' 'WAICV' 'WAICG'}
      if ismember('optimf',ip.UsingDefaults)

        if flagz
          gp=gp_optim(gp,xtr,ytr,'z',ztr(:,size(z,2)),'opt',opt,'loss',inf_method);
        else
          gp=gp_optim(gp,xtr,ytr,'z',ztr,'opt',opt,'loss',inf_method);
        end

      else
        if flagz
          gp=gp_optim(gp,xtr,ytr,'z',ztr(:,size(z,2)),'opt',opt,'loss',inf_method, 'optimf', optimf);
        else
          gp=gp_optim(gp,xtr,ytr,'z',ztr,'opt',opt,'loss',inf_method, 'optimf', optimf);
        end
      end
      w=gp_pak(gp);
      cvws(i,:)=w;
    case 'MCMC'
      if numel(gp.jitterSigma2)>1
        gp=thin(gp,numel(gp.jitterSigma2)-1);
      end
      % Scaled mixture noise model is a special case
      % where we need to modify the noiseSigmas2 vector
      % to a right length
      if isequal(gp.lik.type, 'lik_smt')
        gp.lik.noiseSigmas2 = gp_orig.lik.noiseSigmas2(trindex{i});
        gp.lik.r = gp_orig.lik.r(trindex{i});
        gp.lik.U = gp_orig.lik.U(trindex{i});
        gp.lik.ndata = length(trindex{i});
      end
      % Pick latent values for the training set in this fold
      if isfield(gp,'latentValues')
        if (~isfield(gp.lik, 'nondiagW') || ismember(gp.lik.type, {'Softmax', 'Multinom', ...
            'LGP', 'LGPC'}))
          latentValues=reshape(gp_orig.latentValues, size(y,1), size(y,2));
          gp.latentValues=reshape(latentValues(trindex{i},:), size(y,2)*length(trindex{i}), 1);
          % gp.latentValues=gp_orig.latentValues(trindex{i});
        else
          if ~isfield(gp.lik, 'xtime')
            nl=length(gp.comp_cf);
            gp.latentValues=gp_orig.latentValues(trindex{i}+(0:nl-1)*n);
          else
            ntime=size(gp.lik.xtime,1);
            gp.latentValues=gp_orig.latentValues([1:ntime, (ntime+trindex{i})]);
          end
        end
      end
      if flagz
        gp = gp_mc(gp, xtr, ytr, 'z', ztr(:,size(z,2)), opt);
      else
        gp = gp_mc(gp, xtr, ytr, 'z', ztr, opt);
      end
      nburnin = floor(length(gp.etr)/3);
      gp = thin(gp,nburnin);
    case 'IA'
      if flagz
        [gp,P_TH] = gp_ia(gp, xtr, ytr, 'z', ztr(:,size(z,2)), opt);
      else
        [gp,P_TH] = gp_ia(gp, xtr, ytr, 'z', ztr, opt);
      end
    case 'fixed'
      % nothing to do here
  end

  if iscell(gp)
    gplik=gp{1}.lik;
  else
    gplik=gp.lik;
  end

  switch inf_method
    case {'LOO' 'KFCV' 'WAIC' 'WAICV' 'WAICG' 'MAP'}
      if ~isfield(gplik.fh,'trcov') && isfield(gp.lik.fh,'predcdf')
        for it=1:size(yt,2)
          [Eft, Varft] = gp_pred(gp, xtr, ytr, x, 'tstind', tstind2, 'z', ztr(:,it), 'yt', yt(:,it), 'zt', zt(:,it));
          cdftemp{it}= gp.lik.fh.predcdf(gplik, Eft, Varft,yt(:,it));
        end
      elseif isfield(gplik.fh,'trcov') && isfield(gp.lik.fh,'predcdf')
        for it=1:size(yt,2)
          [Eft, Varft] = gp_pred(gp, xtr, ytr, x, 'tstind', tstind2,'yt', yt(:,it));
          cdftemp{it}= gp.lik.fh.predcdf(gplik, Eft, Varft,yt(:,it));
        end
      else
        error('This likelihood has not been implemented for this function')
      end
    case 'MCMC'
      nsamples=size(gp.etr,1);
      for i2=1:nsamples
        Gp=take_nth(gp,i2);
        gplik=Gp.lik;
        if ~isfield(gplik.fh,'trcov') && isfield(gp.lik.fh,'predcdf')
          for it=1:size(yt,2)
            [Eft, Varft] = gp_pred(Gp, xtr, ytr, x, 'tstind', tstind2, 'z', ztr(:,it), 'yt', yt(:,it), 'zt', zt(:,it));
            cdftemp{it,i2}= gp.lik.fh.predcdf(gplik, Eft, Varft,yt(:,it));
          end
        elseif isfield(gplik.fh,'trcov') && isfield(gp.lik.fh,'predcdf')
          for it=1:size(yt,2)
            [Eft, Varft] = gp_pred(gp, xtr, ytr, x, 'tstind', tstind2,'yt', yt(:,it));
            cdftemp{it,i2}= gp.lik.fh.predcdf(gplik, Eft, Varft,yt(:,it));
          end
        else
          error('This likelihood has not been implemented for this function')
        end
      end
      cdftemp = mat2cell(mean(cell2mat(cdftemp),2),repmat(1043,1,10),1);
    case 'IA'
      nsamples=length(gp);
      for i2=1:nsamples
        Gp=gp{i2};
        gplik=Gp.lik;
        if ~isfield(gplik.fh,'trcov') && isfield(gplik.fh,'predcdf')
          for it=1:size(yt,2)
            [Eft, Varft] = gp_pred(Gp, xtr, ytr, x, 'tstind', tstind2, 'z', ztr(:,it), 'yt', yt(:,it), 'zt', zt(:,it));
            cdftemp{it,i2}= gplik.fh.predcdf(gplik, Eft, Varft,yt(:,it));
          end
        elseif isfield(gplik.fh,'trcov') && isfield(gplik.fh,'predcdf')
          for it=1:size(yt,2)
            [Eft, Varft] = gp_pred(gp, xtr, ytr, x, 'tstind', tstind2,'yt', yt(:,it));
            cdftemp{it,i2}= gplik.fh.predcdf(gplik, Eft, Varft,yt(:,it));
          end
        else
          error('This likelihood has not been implemented for this function')
        end
      end
      cdftemp = mat2cell(sum(bsxfun(@times,cell2mat(cdftemp),P_TH'),2),repmat(1043,1,10),1);
  end

  for it=1:size(yt,2)
    cdf_cv(tstindex{i},it)=cdftemp{it}(tstindex{i},:);
  end

end

cdf=cdf_cv;
end


text:
function kernel_new = dsKernel(kernel, tsub)
%% downsample/upsample the convolution kernel
%% inputs:
%   kernel: struct variable with fields {'kernel_type', 'pars', 'nMax', 'lb', 'ub', 'bound_pars'}
%       kernel_type: string, convolution kernel type. now support {'exp',
%       'exp2', 'vector'}
%       pars: parameters for the selected kernel type
%       nMax: length of the kernel
%       lb:     lower bound for each parameter
%       ub:     upper bound for each parameter
%       bound_pars: logical variable, bound the parameters or not {1, 0}
%% outputs
%   kernel: struct variable

%% Author: Pengcheng Zhou, Carnegie Mellon University, 2016

kernel_new = kernel;
if nargin<2 || isempty(tsub)
    return;
end

%% kernel size
kernel_new.nMax = ceil(kernel.nMax/tsub);

%% kernel type
kernel_type = kernel.type;
if strcmpi(kernel_type, 'exp')
    % single exponential function: ~ exp(-t/tau)
    kernel_new.pars = kernel.pars/tsub;
elseif strcmpi(kernel_type, 'vector')
    % single vector
len_kernel = length(kernel.pars);
    kernel_new.pars = resample(kernel.pars, ceil(len_kernel/tsub), len_kernel);
else
    % differencing of two exponential function:
    % ~  exp(-t/tau_d)-exp(-t/tau_r)
    kernel_new.pars = kernel.pars/tsub;
end

%% lower and upper bounds for parameters
kernel_new.lb = kernel.lb / tsub;
kernel_new.ub = kernel.ub / tsub;

author: flatironinstitute

repo: CaImAn-MATLAB


text:
%
%       HDR Toolbox demo Fusion:
%	   1) Load "Bottles_Small.pfm" HDR image
%	   2) Apply Fusion TMO
%	   3) Show the image without gamma encoding, because it is working in
%	   gamma space
%	   4) Save the tone mapped image as PNG
%
%       Author: Francesco Banterle
%       Copyright 2012-16 (c)
%

clear all;

disp('1) Load the image Bottles_Small.pfm using hdrimread');
img = hdrimread('demos/Bottles_Small.hdr');

disp('2) Apply Fusion Operator by Mertens et al.');
imgTMO = MertensTMO(img);

disp('3) Show the image after fusion, note that there is no need of gamma correction!');
h = figure(1);
set(h,'Name','Exposure fusion by Mertens et al.');
GammaTMO(imgTMO, 1.0, 0.5, 1);

disp('4) Save the tone mapped image as a PNG.');
imwrite(imgTMO, 'demos/output/Bottles_Small_TMO_fusion.png');

author: banterle


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%
temp= sigmoid( X*theta);
for i=1:m
    if temp(i) >=0.5
    p(i)=1;
    else
    p(i)=0;
    end

end

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

end

author: zzlyw

repo: machine-learning-exercises


text:
function [GLRLM] = getGLRLM(ROIonly,levels)
% -------------------------------------------------------------------------
% function [GLRLM] = getGLRLM(ROIonly,levels)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the Gray-Level Run-Length Matrix (GLRLM) of the
% region of interest (ROI) of an input volume. The input volume is assumed
% to be isotropically resampled. Only one GLRLM is computed per scan,
% simultaneously adding up all possible run-lengths in the 13 directions of
% the 3D space. To account for discretization length differences, runs
% constructed from voxels separated by a distance of sqrt(3) increment the
% GLRLM by a value of sqrt(3), runs constructed from voxels separated by a
% distance of sqrt(2) increment the GLRLM by a value of sqrt(2), and runs
% constructed from voxels separated by a distance of 1 increment the GLRLM
% by a value of 1. This function uses other functions from Wei's GLRLM
% toolbox [2].
%
% --> This function is compatible with 2D analysis (language not adapted in the text)
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Galloway, M. M. (1975). Texture analysis using gray level run lengths.
%     Computer Graphics and Image Processing, 4(2), 172–179.
% [2] Wei's GLRLM toolbox: Xunkai Wei, Gray Level Run Length Matrix Toolbox
%     v1.0, Software,Beijing Aeronautical Technology Research Center, 2007.
%     <http://www.mathworks.com/matlabcentral/fileexchange/17482-gray-level-run-length-matrix-toolbox>
% -------------------------------------------------------------------------
% INPUTS:
% - ROIonly: Smallest box containing the ROI, with the imaging data ready
%            for texture analysis computations. Voxels outside the ROI are
%            set to NaNs.
% - levels: Vector containing the quantized gray-levels in the tumor region
%           (or reconstruction levels of quantization).
%
% ** 'ROIonly' and 'levels' should be outputs from 'prepareVolume.m' **
% -------------------------------------------------------------------------
% OUTPUTS:
% - GLRLM: Gray-Level Run-Length Matrix of 'ROIOnly'.
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Martin Vallieres <mart.vallieres@gmail.com>
% - Xunkai Wei <xunkai.wei@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2013
% - Revision: May 2015
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    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) 2007-2012, Xunkai Wei
%     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.
% -------------------------------------------------------------------------

% PRELIMINARY
nLevel = length(levels);
if nLevel > 100
    adjust = 10000;
else
    adjust = 1000;
end
levelTemp = max(levels)+1;
ROIonly(isnan(ROIonly)) = levelTemp; % Last row needs to be taken out of the GLRLM
levels = [levels,levelTemp];

% QUANTIZATION EFFECTS CORRECTION
% In case (for example) we initially wanted to have 64 levels, but due to
% quantization, only 60 resulted.
uniqueVol = round(levels*adjust)/adjust;
ROIonly=round(ROIonly*adjust)/adjust;
NL = length(levels) - 1;

%INITIALIZATION
sizeV = size(ROIonly);
numInit = ceil(max(sizeV)*sqrt(3)); % Max run length
GLRLM = zeros(NL+1,numInit);

% START COMPUTATION
% Directions [1,0,0], [0 1 0], [1 1 0] and [-1 1 0] : 2D directions
% (x:right-left, y:top-bottom, z:3rd dimension)
if numel(size(ROIonly)) == 3
    nComp = sizeV(3); % We can add-up the GLRLMs taken separately in every image in the x-y plane
else
    nComp = 1;
end
for i = 1:nComp
    image = ROIonly(:,:,i);
    uniqueIm = unique(image);
    NLtemp = length(uniqueIm);
    indexRow = zeros(NLtemp,1);
    temp = image;
    for j = 1:NLtemp
        indexRow(j) = find(uniqueIm(j)==uniqueVol);
        image(temp==uniqueIm(j)) = j;
    end

    % [1,0,0]
    GLRLMtemp = rle_0(image,NLtemp);
    nRun = size(GLRLMtemp,2);
    GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun); % Cumulative addition into the GLRLM

    % [0 1 0]
    GLRLMtemp = rle_0(image',NLtemp);
    nRun = size(GLRLMtemp,2);
    GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun); % Cumulative addition into the GLRLM

    % [1 1 0]
    seq = zigzag(image);
    GLRLMtemp = rle_45(seq,NLtemp);
    nRun = size(GLRLMtemp,2);
    GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(2); % Cumulative addition into the GLRLM

    % [-1 1 0]
    seq = zigzag(fliplr(image));
    GLRLMtemp = rle_45(seq,NLtemp);
    nRun = size(GLRLMtemp,2);
    GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(2); % Cumulative addition into the GLRLM
end

if numel(size(ROIonly)) == 3 % 3D DIRECTIONS
    % Directions [0,0,1], [1 0 1] and [-1 0 1]
    % (x:right-left, y:top-bottom, z:3rd dimension)
    nComp = sizeV(1); % We can add-up the GLRLMs taken separately in every image in the x-z plane
    image = zeros(sizeV(3),sizeV(2));
    for i = 1:nComp
        for j = 1:sizeV(3)
            image(j,1:end) = ROIonly(i,1:end,j);
        end
        uniqueIm = unique(image);
        NLtemp = length(uniqueIm);
        indexRow = zeros(NLtemp,1);
        temp = image;
        for j=1:NLtemp
            indexRow(j) = find(uniqueIm(j)==uniqueVol);
            image(temp==uniqueIm(j)) = j;
        end

        % [0,0,1]
        GLRLMtemp = rle_0(image',NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun); % Cumulative addition into the GLRLM

        % [1 0 1]
        seq = zigzag(image);
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(2); % Cumulative addition into the GLRLM

        % [-1 0 1]
        seq = zigzag(fliplr(image));
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(2); % Cumulative addition into the GLRLM
    end

    % Directions [0,1,1] and [0 -1 1]
    % (x:right-left, y:top-bottom, z:3rd dimension)
    nComp = sizeV(2); % We can add-up the GLRLMs taken separately in every image in the y-z plane
    image = zeros(sizeV(1),sizeV(3));
    for i = 1:nComp
        for j = 1:sizeV(3)
            image(1:end,j) = ROIonly(1:end,i,j);
        end
        uniqueIm = unique(image);
        NLtemp = length(uniqueIm);
        indexRow = zeros(NLtemp,1);
        temp = image;
        for j = 1:NLtemp
            indexRow(j) = find(uniqueIm(j)==uniqueVol);
            image(temp==uniqueIm(j)) = j;
        end

        % [0,1,1]
        seq = zigzag(image);
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(2); % Cumulative addition into the GLRLM

        % [0 -1 1]
        seq = zigzag(fliplr(image));
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(2); % Cumulative addition into the GLRLM
    end

    % Four corners: [1,1,1], [-1,1,1], [-1,1,-1], [1,1,-1]
    % (x:right-left, y:top-bottom, z:3rd dimension)
    image = zeros(sizeV(3),sizeV(2));
    temp = rand(sizeV(3),sizeV(2));
    diagTemp = spdiags(temp);
    szDiag = size(diagTemp);
    diagMat1 = zeros(szDiag(1),szDiag(2),sizeV(1));
    diagMat2 = zeros(size(diagTemp,1),size(diagTemp,2),sizeV(1));
    for i = 1:sizeV(1)
        for j = 1:sizeV(3)
            image(j,1:end) = ROIonly(i,1:end,j);
        end
        try
            diagMat1(:,:,i)=spdiags(image);
        catch
            % Add a column at the beginning to prevent errors
            temp=spdiags(image);
            numberDiff=abs(size(temp,2)-size(diagMat1,2));
            if mod(numberDiff,2) % Odd difference number
                temp=padarray(temp,[0,(numberDiff+1)/2,0],0);
                diagMat1(:,:,i)=temp(:,1:end-1);
            else
                diagMat1(:,:,i)=padarray(temp,[0,numberDiff/2,0],0);
            end
        end
        try
            diagMat2(:,:,i)=spdiags(fliplr(image));
        catch
            % Add a column at the beginning to prevent errors
            temp = spdiags(fliplr(image));
            numberDiff = abs(size(temp,2)-size(diagMat2,2));
            if mod(numberDiff,2) % Odd difference number
                temp = padarray(temp,[0,(numberDiff+1)/2,0],0);
                diagMat2(:,:,i) = temp(:,1:end-1);
            else
                diagMat2(:,:,i) = padarray(temp,[0,numberDiff/2,0],0);
            end
        end
    end
    for j = 1:szDiag(2)
        index = (diagMat1(:,j,1)~=0);
        nTemp = sum(index);
        image1 = zeros(sizeV(1),nTemp);
        image2 = zeros(sizeV(1),nTemp);
        for k = 1:sizeV(1)
            image1(k,1:nTemp) = diagMat1(index(1:end),j,k)';
            image2(k,1:nTemp) = diagMat2(index(1:end),j,k)';
        end

        % 2 first corners
        uniqueIm = unique(image1);
        NLtemp = length(uniqueIm);
        indexRow = zeros(NLtemp,1);
        temp = image1;
        for i = 1:NLtemp
            indexRow(i) = find(uniqueIm(i)==uniqueVol);
            image1(temp==uniqueIm(i)) = i;
        end
        seq = zigzag(image1);
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(3); % Cumulative addition into the GLRLM
        seq = zigzag(fliplr(image1));
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(3); % Cumulative addition into the GLRLM

        % 2 last corners
        uniqueIm = unique(image2);
        NLtemp = length(uniqueIm);
        indexRow = zeros(NLtemp,1);
        temp = image2;
        for i = 1:NLtemp
            indexRow(i) = find(uniqueIm(i)==uniqueVol);
            image2(temp==uniqueIm(i)) = i;
        end
        seq = zigzag(image2);
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(3); % Cumulative addition into the GLRLM
        seq = zigzag(fliplr(image2));
        GLRLMtemp = rle_45(seq,NLtemp);
        nRun = size(GLRLMtemp,2);
        GLRLM(indexRow(1:NLtemp),1:nRun) = GLRLM(indexRow(1:NLtemp),1:nRun) + GLRLMtemp(1:NLtemp,1:nRun)*sqrt(3); % Cumulative addition into the GLRLM
    end
end

% REMOVE UNECESSARY COLUMNS
GLRLM(end,:) = [];
stop = find(sum(GLRLM),1,'last');
GLRLM(:,(stop+1):end) = [];

end

author: mvallieres

repo: radiomics


text:
function [f, logscale] = fwd1(engine, ev, t)
% Forwards pass for slice 1.

bnet = bnet_from_engine(engine);
ss = bnet.nnodes_per_slice;

CPDpot = cell(1,ss);
for n=1:ss
  fam = family(bnet.dag, n, 1);
  e = bnet.equiv_class(n, 1);
  CPDpot{n} = convert_to_pot(bnet.CPD{e}, engine.pot_type, fam(:), ev);
end
f.t = t;
f.evidence = ev;

pots = CPDpot;
slice1 = 1:ss;
CPDclqs = engine.clq_ass_to_node1(slice1);

[f.clpot, f.seppot] =  init_pot(engine.jtree_engine1, CPDclqs, CPDpot, engine.pot_type, engine.observed1);
[f.clpot, f.seppot] = collect_evidence(engine.jtree_engine1, f.clpot, f.seppot);
for c=1:length(f.clpot)
  [f.clpot{c}, ll(c)] = normalize_pot(f.clpot{c});
end
logscale = ll(engine.root1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function value = r8_chi_pdf ( df, rval )

%*****************************************************************************80
%
%% R8_CHI_PDF evaluates the PDF of a chi-squared distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, real DF, the degrees of freedom.
%    0.0 < DF.
%
%    Input, real RVAL, the point where the PDF is evaluated.
%
%    Output, real VALUE, the value of the PDF at RVAL.
%
  if ( df <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_CHI_PDF - Fatal error!\n' );
    fprintf ( 1, '  Degrees of freedom must be positive.\n' );
    error ( 'R8_CHI_PDF - Fatal error!' );
  end if

  if ( rval <= 0.0 )

    value = 0.0;

  else

    temp2 = df * 0.5;

    temp1 = ( temp2 - 1.0 ) * log ( rval ) - 0.5 * rval ...
      - temp2 * log ( 2.0 ) - r8_gamma_log ( temp2 );

    value = exp ( temp1 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y2 = discrim_plot(discf,y,varargin)
% y2 = discrim_plot(discf,y,[new figure],[column to do partial corr plot of])
%
% Plots correlation or partial correlation
% Color codes by median split of y
%
% see also:
% cluster_discrim
% cluster_discrim_montage

dofig = 1; dopr = 0;
if length(varargin) > 0, dofig = varargin{1};,end
if length(varargin) > 1, dopr = varargin{2};,end

y2 = mediansplit(y);

if dopr
    [discf(:,1),y,r,p,rrob,prob] = partialcor(discf,y,dopr);
    fprintf(1,'\nCalculating partial correlation.\n')
end

if dofig, figure('Color','w');,end

plot_correlation_samefig(discf(:,1),y,[],'ko',0,1);
wh = find(y2>0); wh2 = find(y2<0);
hold on; plot(discf(wh2,1),y(wh2),'bo','MarkerSize',10,'LineWidth',2);
hold on; plot(discf(wh,1),y(wh),'ro','MarkerSize',10,'LineWidth',2);
xlabel('Discriminant function'); ylabel('Behavior');

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [quality probs] = biqi(im)

%========================================================================
%
% -----------COPYRIGHT NOTICE STARTS WITH THIS LINE------------
% Copyright (c) 2009 The University of Texas at Austin
% All rights reserved.
%
% Permission is hereby granted, without written agreement and without license or royalty fees, to use, copy,
% modify, and distribute this code (the source files) and its documentation for
% any purpose, provided that the copyright notice in its entirety appear in all copies of this code, and the
% original source of this code, Laboratory for Image and Video Engineering (LIVE, http://live.ece.utexas.edu)
% and Center for Perceptual Systems (CPS, http://www.cps.utexas.edu) at the University of Texas at Austin (UT Austin,
% http://www.utexas.edu), is acknowledged in any publication that reports research using this code. The research
% is to be cited in the bibliography as:
%
% 1. A. K. Moorthy and A. C. Bovik, "A Modular Framework for Constructing Blind
% Universal Quality Indices", submitted to IEEE Signal Processing Letters (2009).
%
% 2. A. K. Moorthy and A. C. Bovik, "BIQI Software Release",
% URL: http://live.ece.utexas.edu/research/quality/biqi.zip, 2009.
%
% IN NO EVENT SHALL THE UNIVERSITY OF TEXAS AT AUSTIN BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL,
% OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF THIS DATABASE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF TEXAS
% AT AUSTIN HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% THE UNIVERSITY OF TEXAS AT AUSTIN SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE DATABASE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS,
% AND THE UNIVERSITY OF TEXAS AT AUSTIN HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
%
% -----------COPYRIGHT NOTICE ENDS WITH THIS LINE------------%

%Author  : Anush Krishna Moorthy
%Version : 1.0
%
%The authors are with the Laboratory for Image and Video Engineering
%(LIVE), Department of Electrical and Computer Engineering, The
%University of Texas at Austin, Austin, TX.
%
%Kindly report any suggestions or corrections to anushmoorthy@gmail.com
%
%========================================================================
%
% This is a demonstration of the Blind Image Quality Index (BIQI) .
% It is an implementation of the BIQI in the reference.
% The algorithm is described in:
% A. K. Moorthy and A. C. Bovik, "A Modular Framework for Constructing Blind
% Universal Quality Indices", submitted to IEEE Signal Processing Letters (2009).
%
%You can change this program as you like and use it anywhere, but please
%refer to its original source (cite our paper and our web page at
% http://live.ece.utexas.edu/research/quality/biqi.zip).
%
%Input : A test 8bits/pixel grayscale image loaded in a 2-D array
%Output: A quality score of the image. The score typically has a value
%        between 0 and 100 (0 represents the best quality, 100 the worst).
%
%Usage:
%
%1. Load the image, for example
%
%   image = rgb2gray(imread('testimage.jpg'));
%
%2. Call this function to calculate the quality score:
%
%    quality = biqi(image)
%
% Dependencies:
% MATLAB Wavelet  Toolbox
% You may need the MATLAB Image Processing Toolbox
% Binaries: svm-train, svm-scale (from LibSVM) - provided with release
% Other m files: jpeg_quality_score.m (provided with release)
% Data files: range2, range2_wn, range2_blur, range2_jp2k, model_89,
% model_89_wn, model_89_blur, model_89_jp2k, rang2_ff model_ff
%========================================================================

import biqi.*

if(size(im,3)~=1)
    im = rgb2gray(im);
end

%% First compute statistics
num_scales = 3; % somethings are hardcoded for this...please be careful when changing.
gam = 0.2:0.001:10;
r_gam = gamma(1./gam).*gamma(3./gam)./(gamma(2./gam)).^2;

[C S] = wavedec2(im,num_scales,'db9');
for p = 1:num_scales
    [horz_temp,vert_temp,diag_temp] = detcoef2('all',C,S,p) ;
    horz(p) = {[horz_temp(:)]};
    diag(p) = {[diag_temp(:)]};
    vert(p) = {[vert_temp(:)]};

    h_horz_curr  = cell2mat(horz(p));
    h_vert_curr  = cell2mat(vert(p));
    h_diag_curr  = cell2mat(diag(p));

    mu_horz(p) = mean(h_horz_curr);
    sigma_sq_horz(p)  = mean((h_horz_curr-mu_horz(p)).^2);
    E_horz = mean(abs(h_horz_curr-mu_horz(p)));
    rho_horz = sigma_sq_horz(p)/E_horz^2;
    [min_difference, array_position] = min(abs(rho_horz - r_gam));
    gam_horz(p) = gam(array_position);

    mu_vert(p) = mean(h_vert_curr);
    sigma_sq_vert(p)  = mean((h_vert_curr-mu_vert(p)).^2);
    E_vert = mean(abs(h_vert_curr-mu_vert(p)));
    rho_vert = sigma_sq_vert(p)/E_vert^2;
    [min_difference, array_position] = min(abs(rho_vert - r_gam));
    gam_vert(p) = gam(array_position);

    mu_diag(p) = mean(h_diag_curr);
    sigma_sq_diag(p)  = mean((h_diag_curr-mu_diag(p)).^2);
    E_diag = mean(abs(h_diag_curr-mu_diag(p)));
    rho_diag = sigma_sq_diag(p)/E_diag^2;
    [min_difference, array_position] = min(abs(rho_diag - r_gam));
    gam_diag(p) = gam(array_position);
end
rep_vec = [mu_horz mu_vert mu_diag sigma_sq_horz sigma_sq_vert sigma_sq_diag gam_horz gam_vert gam_diag];
rep_vec(:,1:9) = []; % remove the means...
%% Now classify

fid = fopen('test_ind.txt','w')
for j = 1:size(rep_vec,1)
    fprintf(fid,'%d ',j);
    for k = 1:size(rep_vec,2)
        fprintf(fid,'%d:%f ',k,rep_vec(j,k));
    end
    fprintf(fid,'\n');
end
fclose(fid);

system(['/usr/local/Cellar/libsvm/3.22/bin/svm-scale -r +biqi/range2 test_ind.txt >> test_ind_scaled']);
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-predict -b 1 test_ind_scaled +biqi/model_89 output_89']);
delete test_ind.txt test_ind_scaled

%% Quality along each dimension

% Write out SVM compatible

fid = fopen('test_ind.txt','w');
for j = 1:size(rep_vec,1)
    fprintf(fid,'%f ',j);
    for k = 1:size(rep_vec,2)
        fprintf(fid,'%d:%f ',k,rep_vec(j,k));
    end
    fprintf(fid,'\n');
end
fclose(fid);

% Jp2k quality
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-scale -r +biqi/range2_jp2k test_ind.txt >> test_ind_scaled']);
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-predict  -b 1 test_ind_scaled +biqi/model_89_jp2k output_blur']);
load output_blur
jp2k_score = output_blur;
delete output_blur test_ind_scaled

% JPEG quality
jpeg_score  = jpeg_quality_score(im);

% WN quality
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-scale -r +biqi/range2_wn test_ind.txt >> test_ind_scaled']);
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-predict -b 1 test_ind_scaled +biqi/model_89_wn output_blur']);
load output_blur
wn_score = output_blur;
delete output_blur test_ind_scaled

% Blur quality
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-scale -r +biqi/range2_blur test_ind.txt >> test_ind_scaled']);
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-predict  -b 1 test_ind_scaled +biqi/model_89_blur output_blur']);
load output_blur
blur_score = output_blur;
delete output_blur test_ind_scaled

% FF quality
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-scale -r +biqi/range2_ff test_ind.txt >> test_ind_scaled']);
system(['/usr/local/Cellar/libsvm/3.22/bin/svm-predict  -b 1 test_ind_scaled +biqi/model_89_ff output_blur']);
load output_blur
ff_score = output_blur;

delete output_blur
delete test_ind.txt test_ind_scaled

%% Final pooling

% figure out probabilities
fid = fopen('output_89','r');
fgetl(fid);
C = textscan(fid,'%f %f %f %f %f %f');
output = [C{1} C{2} C{3} C{4} C{5} C{6}];
fclose(fid);
probs = output(:,2:end);
scores  = [jp2k_score jpeg_score wn_score blur_score ff_score];
quality = sum(probs.*scores,2);
delete output_89
clc

author: dsoellinger


text:
function p21_default ( dim_num )

%*****************************************************************************80
%
%% P21_DEFAULT sets default values for problem 21.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p21_i4 ( 'D', '*', i4 );

  r8 = 0.0;
  p21_r8 ( 'D', '*', r8 );

  i4vec = [];
  p21_i4vec ( 'D', '*', dim_num, i4vec );

  return
end


text:
% Gives estimate of T0-normalized LF parameters given an Rd value
%
% Description
%  The Rd shape parameter [1] of the Liljencrants-Fant (LF) glottal model [2]
%  expresses a regression of the original shape parameters {te,tp,ta}. This
%  function gives the shape parameters {te,tp,ta} which correspond to a given Rd
%  value.
%
% Input
%  Rd  : The Rd shape parameter to convert.
%
% Output
%  te  : Glottal shape parameter, see [1]p.6 (assuming T0=1 ! (T0-normalized))
%  tp  : Glottal shape parameter, see [1]p.6 (assuming T0=1 ! (T0-normalized))
%  ta  : Glottal shape parameter, see [1]p.6 (assuming T0=1 ! (T0-normalized))
%
% See Also
%  gfm_spec_lf.m
%
% References
%  [1] G. Fant, "The LF-model revisited. Transformations and frequency domain
%      analysis", STL-QPSR 36(2-3):119-156, 1995.
%  [2] G. Fant, J. Liljencrants and Q. Lin, "A four-parameter model of glottal
%      flow", STL-QPSR, vol. 4, pp. 1-13, 1985.
%
% Copyright (c) 2011 University of Crete - Computer Science Department
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file 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 Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [te tp ta] = Rd2tetpta(Rd)

	Rap = (-1+4.8.*Rd)/100;                             % [1](2)
	Rkp = (22.4+11.8.*Rd)/100;                          % [1](3)
	Rgp = 1./(4*((0.11.*Rd./(1/2+1.2.*Rkp))-Rap)./Rkp); % [1] indirectly (4)

	tp = 1./(2.*Rgp); % [1]p.121
	te = tp.*(Rkp+1); % [1]p.121
	ta = Rap;         % [1]p.121

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function cirpva()
    % unclear purpose, having to do with a circle selection
%   This subroutine "circle"  selects the Ni closest earthquakes
%   around a interactively selected point.  Resets ZG.newcat and ZG.newt2
%   Operates on "primeCatalog".
 % turned into function by Celso G Reyes 2017

ZG=ZmapGlobal.Data; % used by get_zmap_globals

%  Input Ni:
%
report_this_filefun();
ZG=ZmapGlobal.Data;

delete(findobj('Tag','plos1'));

axes(h1)
%zoom off

titStr ='Selecting EQ in Circles                         ';
messtext= ...
    ['                                                '
    '  Please use the LEFT mouse button              '
    ' to select the center point.                    '
    ' The "ni" events nearest to this point          '
    ' will be selected and displayed in the map.     '];

msg.dbdisp(messtext, titStr);

% Input center of circle with mouse
%
[xa0,ya0]  = ginput(1);

stri1 = [ 'Circle: ' num2str(xa0,5) '; ' num2str(ya0,4)];
stri = stri1;
pause(0.1)

if met == 'ni'
    % take first ni and sort by time
    [ZG.newt2, max_rad] = ZG.primeCatalog.selectClosestEvents(ya0, xa0, [], ni);
    messtext = ['Radius of selected Circle:' num2str(maxrad)  ' km' ];
    disp(messtext)
elseif  met == 'ra'
    ZG.newt2 = ZG.primeCatalog.selectRadius(ya0, xa0, ra,'kilometer');
    messtext = ['Number of selected events: ' num2str(ZG.newt2.Count)  ];
    disp(messtext)
elseif met == 'ti'
    global t1 t2 t3 t4
    ZG.newt2 = copy(ZG.primeCatalog);
    lt =  ZG.newt2.Date >= t1 &  ZG.newt2.Date <t2 ;
    bdiff(ZG.newt2.subset(lt));
    ZG.hold_state=true;
    lt =  ZG.newt2.Date >= t3 &  ZG.newt2.Date <t4 ;
    bdiff(ZG.newt2.subset(lt));

end
R2 = ra;

%
% plot Ni clostest events on map as 'x':

set(gca,'NextPlot','add')
plot(ZG.newt2.Longitude,ZG.newt2.Latitude,'xk','Tag','plos1');

% plot circle containing events as circle
x = -pi-0.1:0.1:pi;
pl = plot(xa0+sin(x)*R2/(cosd(ya0)*111), ya0+cos(x)*R2/(cosd(ya0)*111),'k')

set(gcf,'Pointer','arrow')

%
newcat = ZG.newt2;                   % resets ZG.newcat and ZG.newt2

        ZG=ZmapGlobal;
        ctp=CumTimePlot(ZG.newt2);
        ctp.plot();
MyPvalClass.pvalcat

end

author: CelsoReyes


text:
function cvx_optval = norm_nuc( X ) %#ok

%NORM_NUC   Internal cvx version.

error( nargchk( 1, 1, nargin ) );
if ndims( X ) > 2,
    error( 'norm_nuc is not defined for N-D arrays.' );
elseif ~cvx_isaffine( X ),
    error( 'Input must be affine.' );
end

%
% Construct problem
%

[ m, n ] = size( X ); %#ok
cvx_begin sdp
    variable W1(m,m) symmetric
    variable W2(n,n) symmetric
    minimize(0.5*(trace(W1)+trace(W2)));
    [W1,X;X',W2] >= 0; %#ok
cvx_end

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

repo: REST


text:
function test_issue1387

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_read_event ft_read_header ft_preprocessing

%%
% load data and parse the events using '5*nanmedian' threshold on 8 channels
fileName       = dccnpath('/home/common/matlab/fieldtrip/data/test/issue1387/Test_4.edf');
bitChannels    = 3:10;
header         = ft_read_header(fileName);
labels         = header.label(bitChannels);
data_events	   = ft_read_event(fileName,'header',header,'detectflank','up','chanindx',bitChannels,'threshold','5*nanmedian');

%%
% load the raw data
cfg            = [];
cfg.dataset    = fileName;
cfg.continuous = 'yes';
cfg.channel    = 'all';
data           = ft_preprocessing(cfg);

%%
% parse out the event information for the 1st bit channel
idx            = cellfun(@(x)strcmpi(x,labels{1}),{data_events.type});
event.label	   = labels{1};
event.idx      = find(idx==1);
event.evnt     = data_events(event.idx);
event.samples  = [event.evnt.sample];
event.times	   = data.time{1}(event.samples);

%%
% ASSERT: There were 78 TTLs on the first channel, make sure this is the case, error otherwise.
assert(length(event.times)==78,'Number of events should be 78!')

%%
% Plot the events on top of the raw data for visual inspection
tm             = data.time{1};
ch             = data.trial{1}(3,:);
baseline       = median(ch(1:100));
ch             = ch - baseline;
ch             = ch / max(ch);
figure;
plot(tm,ch,'k-');
hold on;
y = repmat(0.5, [1 length(event.times)]);
plot(event.times,y,'ro','MarkerSize',6)
xlim([10 50]);
ylim([-0.05 1.05]);
title('TTL channel 1 should have 78 event markers aligned to the upstates')
xlabel('Time (s)')
ylabel('Normalised Amplitude')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip