text:
% 'region' is a class of objects that contains information on groups of
% voxels ("regions") defined in various ways (contiguous voxels above a
% threshold in an analysis, based on atlases, etc.)
%
% the region class replaces the "clusters" structure in the scnlab toolbox,
% and all or almost all of the "clusters" functions should work for regions
% structures as well.  The advantage of 'region' is that the class
% definition can help constrain the use and provide additional error
% checking, etc.
%
% Defining a region and initializing:
%
% *Usage:*
%   - r = region(obj1, [obj2], [keywords] )
%     - obj1: [fmri_data/statistic_image object to define regions]
%             Can also be char array of image filename
%     - obj2: Optional: fmri_data/statistic_image object to extract data from
%     - keywords: Optional: 'unique_mask_values' or 'contiguous_regions'
%     - 'noverbose' : suppress verbose output
%
%   - cl = region;  % generates an empty structure
%                   % You can add fields yourself if you want to, but best to
%                   define based on existing file name or image_vector object
%
% Define regions based on continuous voxel values (in this example image,
% there is only one set of contiguous voxels, so 1 region...)
%
%   - mask_image = which('brainmask.nii');
%   - cl = region(mask_image);
%
% *Inputs:*
%
% There are two ways to define which voxels are grouped into a 'region',
% which becomes an element of the region variable cl.
%
% enter 'contiguous_regions' -> group by contiguous blobs
%
% or    'unique_mask_values' -> group by unique values in mask .dat field
%
% *Note:* 'contiguous_regions' uses contiguity/clustering information stored
% in mask.volInfo.cluster, which may not have veridical contiguity info if
% you have manipulated it or incorporated anatomical information in working with
% the mask object, or if you have borrowed the volInfo structure from
% another source in creating it.
%
% *Examples:*
%
% Define regions based on continuous values in an anatomical mask:
%   - mask_image = which('atlas_labels_combined.img');
%   - cl = region(mask_image, 'unique_mask_values');
%
% Resample mask_image to space of data_comb fmri_data object, and extract
% averages.  Space is defined by data_comb.
%   - cl = extract_roi_averages(data_comb, mask_image, 'unique_mask_values');
%
% Define regions based on unique voxels values in mask_image, and extract
% data stored in data_comb object, resampled to space of mask_image.
% Space is defined by mask_image:
%   - cl = region(mask_image, data_comb, 'unique_mask_values');
%
% Reslice mask to space of functional images and define regions based on
% mask values in the functional space (good for extracting data, etc.)
%   - mask_image = which('anat_lbpa_thal.img');
%   - mask = fmri_mask_image(mask_image);
%   - mask = resample_to_image_space(mask, image_names(1, :));
%   - cl = region(mask);
%
% *Methods*
%
% Try typing methods(cl)
%
% methods for 'regions' include:
%
% Visualization methods:
%     montage, orthviews, surf, etc.
%
% *Programmers' Notes:*
%
% 8/3/2015 : Tor Wager: Fixed bug when applying region to thresholded
% statistic_image object.  Did not consider thresholding.
%
% 5/24/2017: Tor and Phil Kragel: first attempt to make compatible with
% use of enforce_variable_types method.
%
% 7/2018 : Tor - check for empty mask and skip data extraction if so

classdef region

    properties
        title = 'Untitled';     % Title of region object
        shorttitle = 'region';
        descrip1 = 'Region: a group of voxels.'
        descrip2 ='Some methods: orthviews, montage, surface, extract_data, table';

        XYZ
        XYZmm
        val
        val_descrip = 'Description of values for each voxel in val field.';

        Z
        Z_descrip = 'Legacy values for each voxel; Z-scores or other max stat';

        threshold % legacy, for compatibility
        voxSize
        M
        dim
        numVox
        numpeaks

        center
        mm_center

        timeseries
        contrastdata
        dat
        all_data

        source_images

        custom_info1
        custom_info1_descrip

        custom_info2
        custom_info2_descrip

    end % properties

    methods

        % class constructor method
        % takes maskinput in one of several forms:
        % 1) mask object
        % 2) mask image file
        % 3) mask vector data (???)
        %
        % takes one of two string inputs for how to define regions (see below)
        % 'contiguous_regions' (default)
        % 'unique_mask_values'

        function obj = region(maskinput, varargin)

            % initialize empty
            % return if no additional args

            obj.title = 'Untitled';
            obj.shorttitle = 'region';
            obj.descrip1 = 'Region: a group of voxels.';
            obj.descrip2 ='Some methods: orthviews, montage, surface, extract_data, table';
            obj.XYZ = [];
            obj.XYZmm = [];
            obj.val = [];
            obj.Z = [];
            obj.val_descrip = 'Values for each voxel, usually max stat.';
            obj.Z_descrip = 'Legacy values for each voxel; Z-scores or other max stat';
            obj.threshold = []; % legacy, for compatibility
            obj.voxSize = [];
            obj.M = [];
            obj.dim = [];
            obj.numVox = 0;
            obj.numpeaks = NaN;
            obj.center = [];
            obj.mm_center = [];
            obj.timeseries = [];
            obj.contrastdata = [];
            obj.dat = [];
            obj.all_data = [];
            obj.source_images = char([]);
            obj.custom_info1 = [];
            obj.custom_info1_descrip = 'Mask image name for region definition';
            obj.custom_info2 = [];
            obj.custom_info2_descrip = char([]);

            if nargin == 0
                return
            end

            % ---------------------------------
            % define mask object
            % ---------------------------------

            % Load mask into object if needed, apply .sig field for
            % statistic_image objects

            [mask, maskData] = prep_mask_object(maskinput);

            % Check for an incompatibility. Don't try to fix here, force
            % proper intended use.
            if isa(maskinput, 'atlas') && isa(maskinput.dat, 'int32')

                %atlas_obj.dat = double(atlas_obj.dat);
                error('Do not use region() for atlas objects. Use atlas2region instead.');

            end

            % ---------------------------------
            % define what to average over
            % ---------------------------------

            average_over = 'contiguous_regions'; %'contiguous_regions'  or 'unique_mask_values';
            doverbose = true;

            for varg = 1:length(varargin)
                if ischar(varargin{varg})
                    switch varargin{varg}

                        % reserved keywords
                        case 'contiguous_regions', average_over = 'contiguous_regions';
                        case 'unique_mask_values', average_over = 'unique_mask_values';

                        case 'noverbose'
                            doverbose = false;

                        otherwise
                            disp('region class constructor: Illegal string value for average_over.');
                            fprintf('You entered ''%s''\n Valid values are %s or %s\n', varargin{varg}, '''contiguous_regions''', '''unique_mask_values''');
                            error('Exiting');
                    end
                elseif isa(varargin{varg}, 'image_vector')
                        % data object to extract from at end

                        dataobj = varargin{varg};
                        varargin{varg} = [];

                        % Note: If you have manipulated an image_vector (e.g., fmri_data,
                        % statistic_image) object and eliminated some voxels, in order to create
                        % regions, contiguous voxels are automatically reparsed into regions using
                        % fmri_data.reparse_contiguous below.

                        dataobj = reparse_contiguous(dataobj, 'nonempty');
                else
                    % unknown input
                end
            end

            % If extracting data, we'll recreate the regions in
            % extract_roi_averages, so do that here and then exit.
            % Otherwise, continue to region definition.

            if exist('dataobj', 'var')
                % extract data
                if doverbose
                    disp('> Found image data, extracting region averages.');
                end

                dataobj = replace_empty(dataobj); % may need to do this to get voxels to line up

                cs = compare_space(dataobj, mask);

                % if empty, skip
                isemptymask = isempty(mask.dat) || all(mask.dat(:) == 0);
                if isemptymask
                    if doverbose
                        disp('No in-region voxels from which to extract data.');
                    end
                    return
                end

                if cs == 3
                    disp('Spaces for data object and mask object line up, but voxel numbers do not. Check.');
                    if doverbose, disp('> Resampling to mask space first.'); end
                    dataobj = resample_space(dataobj, mask); % resample data to mask space
                elseif cs
                    if doverbose, disp('> Resampling to mask space first.'); end
                    dataobj = resample_space(dataobj, mask); % resample data to mask space
                end

                if doverbose
                    obj = extract_roi_averages(dataobj, mask, average_over);
                else
                    obj = extract_roi_averages(dataobj, mask, average_over, 'noverbose');
                end

                return
            end

            % ---------------------------------
            % get which values to save later
            % ---------------------------------

            maskValues = maskData;  % values to save in .val field
            maskZ = maskData;       % values to save in .Z field later

            val_descrip = 'Input mask image values for each voxel.';
            Z_descrip = 'Input mask image values for each voxel.';

            if isa(mask, 'statistic_image')

                 val_descrip = 'Statistic effect value (.dat) for each voxel.';

                 mask = replace_empty(mask);            % after this p and dat now have all in-mask voxels

                 if ~isempty(mask.dat)
                     % Use values from image. These can be more precise and
                     % we don't want to convert them to Z-values if we
                     % don't have to.

                     maskZ = mask.dat;
                     Z_descrip = mask.type;

                 elseif ~isempty(mask.p)

                     if doverbose('Converting .dat data values to Z-values based on P-values'); end

                     Z_descrip = 'Z-score for each voxel.';
                     mask.p(mask.p==0) = eps;               % not to have Inf values in Z field
                     maskZ = sign(mask.dat) .* norminv(1 - mask.p ./ 2); % Z-score based on p-value, assuming two-tailed p-vals.

                 end

            end

            % ---------------------------------
            % get unique values for voxel grouping code
            % ---------------------------------

            switch average_over

                % Define integer codes for sets of voxels to average over.

                case 'unique_mask_values'

                    maskData = round(maskData);
                    u = unique(maskData)'; u(u == 0) = [];
                    nregions = length(u);

                    if doverbose
                        fprintf('Grouping voxels with unique mask values, assuming integer-valued mask: %3.0f regions\n', nregions);
                    end

                case 'contiguous_regions'

                    isinmask = maskData ~= 0 & ~isnan(maskData);

                    % re-make cluster ID for in-mask voxels
                    mask.volInfo.cluster(isinmask) = spm_clusters(mask.volInfo.xyzlist(isinmask, :)')';

                    % Old, no longer needed with newer SPM
                    %                     if sum(isinmask) < 50000
                    %                         mask.volInfo.cluster(isinmask) = spm_clusters(mask.volInfo.xyzlist(isinmask, :)')';
                    %                     else
                    %                         % don't, and print warning
                    %                         mask.volInfo.cluster(isinmask) = ones(sum(isinmask), 1);
                    %                         disp('Warning: spm_cluster will not parse clusters for masks with > 50000 voxels.');
                    %                     end

                    u = unique(mask.volInfo.cluster(isinmask)); u(u == 0) = [];
                    maskData(isinmask) = mask.volInfo.cluster(isinmask);
                    nregions = length(u);

                    if doverbose
                        fprintf('Grouping contiguous voxels: %3.0f regions\n', nregions);
                    end

                otherwise
                    error('This should never happen.');
            end

            % ---------------------------------
            % Now define the regions
            % ---------------------------------

            obj(1:nregions) = obj;  % fill in all empty fields

            for i = 1:nregions
                imgvec = maskData == u(i);

                obj(i).title = sprintf('Region %3.0f', i);
                obj(i).shorttitle = sprintf('Region%03d', i);
                obj(i).XYZ = double(mask.volInfo.xyzlist(imgvec,:))';

                myXYZ = double(obj(i).XYZ);                         % 5/24/17 tor&phil: recast as double. future: could recast XYZmm as int16
                obj(i).XYZmm = voxel2mm(myXYZ, mask.volInfo.mat);

                obj(i).val = maskValues(imgvec); %ones(1,size(obj(i).XYZ,2));
                obj(i).Z = maskZ(imgvec)'; %ones(1,size(obj(i).XYZ,2));

                obj(i).val_descrip = val_descrip;
                obj(i).Z_descrip = Z_descrip;

                obj(i).voxSize = abs(diag(mask.volInfo.mat(1:3, 1:3)));
                obj(i).M = mask.volInfo.mat;
                obj(i).dim = mask.volInfo.dim;
                obj(i).numVox = size(obj(i).XYZ, 2);

                obj(i).center = center_of_mass(myXYZ, double(obj(i).Z));
                obj(i).mm_center = center_of_mass(obj(i).XYZmm, double(obj(i).Z));

                if ~isempty(mask.volInfo.fname)
                    [~, ff, ee] = fileparts(mask.volInfo.fname);
                    obj(i).custom_info1 = [ff ee];
                    obj(i).custom_info1_descrip = 'Mask image name for region definition';
                else
                    obj(i).custom_info1 = 'No filename to associate';
                    obj(i).custom_info1_descrip = '';
                end

            end

            if nregions == 0
                % obligatory things even for empty regions
                obj(1).voxSize = abs(diag(mask.volInfo.mat(1:3, 1:3)));
                obj(1).M = mask.volInfo.mat;
                obj(1).dim = mask.volInfo.dim;

                obj(1).val_descrip = val_descrip;
                obj(1).Z_descrip = Z_descrip;

            end

        end % class constructor function

    end % methods

end % classdef

% Load mask into object if needed, apply .sig field for
% statistic_image objects
% * NOTE: some functionality may be redundant here; to-do is to refactor
% and run unit test

function  [mask, maskData] = prep_mask_object(maskinput)

if isa(maskinput, 'char') % string file name
    mask = fmri_mask_image(maskinput);

elseif isa(maskinput, 'image_vector')
    % case {'fmri_data', 'fmri_mask_image', 'statistic_image', 'image_vector'}

    % special for thresholded stats images: use threshold

    % if sig field exists, use sig voxels only
    if isa(maskinput, 'statistic_image') && ~isempty(maskinput.sig)

        for img = 1:size(maskinput.dat, 2)

            maskinput.dat(~maskinput.sig(:, img), img) = 0;

        end

    end

    mask = maskinput;

elseif isa(maskinput, 'atlas')
    error('Use atlas2region to convert an atlas object to a region object.');

else
    error('region class constructor: unknown mask input type.')
end

% Mask data can already be reduced to those indexed by
% wh_inmask or not.
% In addition, voxels/images with empty values may be removed.
% So insert those first.
% Note: .dat can sometimes have 2+ cols, so use only first one
mask = replace_empty(mask);

if size(mask.dat, 2) > 1
    disp('Warning: Mask has multiple images, will use first only.');
end

if size(mask.dat, 1) == mask.volInfo.n_inmask
    maskData = mask.dat(:, 1);

elseif size(mask.dat, 1) == mask.volInfo.nvox
    % We have a full-length vector
    maskData = mask.dat(mask.volInfo.wh_inmask(:, 1));
else
    error('Illegal size for mask.dat, because it does not match its volInfo structure.')
end

% If statistic_image, we need to consider thresholding
% so apply .sig field
if isa(mask, 'statistic_image')

    % old... we could delete these.. but just in case, I just comment them out for now (Wani).
    % if size(mask.dat, 1) == mask.volInfo.nvox
    %    error('statistic_image objects should not have data vector (.dat) the length of full image space.')
    % end

    maskData = maskData .* mask.sig(:, 1);

end

% may need to reparse contiguous voxels in the mask.
mask = reparse_contiguous(mask, 'nonempty');

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [smaIn, eIn, hHat, sHat, smaOut, eOut, oHat] = computeMultiFlyByParameters(vInfIn, vInfOut, gmu)
%computeFlyByParameters Summary of this function goes here
%   Detailed explanation goes here

    normVInfIn = sqrt(sum(abs(vInfIn).^2,1));
    normVinfOut = sqrt(sum(abs(vInfOut).^2,1));

    sHat = bsxfun(@rdivide,vInfIn,normVInfIn);
    energyIn = normVInfIn.^2./2;
    smaIn = -gmu./(2*energyIn);

    oHat = bsxfun(@rdivide,vInfOut,normVinfOut);
    energyOut = normVinfOut.^2./2;
    smaOut = -gmu./(2*energyOut);

    crossSO = cross(sHat, oHat);
    normCrossSO = sqrt(sum(abs(crossSO).^2,1));
    hHat = bsxfun(@rdivide,crossSO,normCrossSO);
    delta = acos(dot(sHat,oHat)); %may need the dim argument

    e1Eqn = @(e2) 1+(smaOut./smaIn).*(e2-1);
    f = @(e2) delta - asin(1./e1Eqn(e2)) - asin(1./e2);
    lb = 1.0 * ones(1,length(delta));
    ub = 1000000000 * ones(1,length(delta));

    [eOut, ~, ~] = bisection(f,lb,ub, [], 1E-8);
    eIn = e1Eqn(eOut);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function stats = rsa_regression(obj,design,study,varargin)
% Representational similarity-based analysis, including inferences about a stimulus/task model
% Constructs a rep. dissim. matrix (RDM) based on spatial covariance across images.
% Takes a stimulus/experimental design (design), which is a set of binary regressors specifying groups of images.
% Include a second grouping variable (study) of integers coding for blocks of images to resample within.
% The function regresses the RSA matrix on the design and returns bootstrapped
% statistics on each design regressor.
%   This method can be used to analyze generalizability across constructs and studies.
% e.g., see reference to Kragel et al. 2018 below. In this case, this method tests whether
% patterns of activity in an fmri data object are generalizable across
% different groupings specified in the matrix 'design'. 'study' is an integer
% valued vector grouping images according to study. This script does not
% currently implement repeated measures designs (only include one image per
% subject). See Kragel et al. (2018) Nature Neuroscience for details.
%
% Usage:
% ::
%
%    stats = rsa_regression(obj,design);
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2018 Phil Kragel
%
%     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/>.
% ..
%
% :Inputs:
%
%   **obj:**
%        An image object with one or more images loaded
%
%   **design:**
%        A binary matrix where each column corresponds to a grouping
%        variable and each row corresponds to an observation.
%
% :Optional inputs:
%
%	**'average_Euclidean'**
%		Use Euclidean distance of the response averaged across all voxels (univariate distance)
%	**'euclidean'**
%		Use Euclidean distance computed using all voxels
%	**'cosine'**
%		Compute distance using cosine similarity computed using all voxelst
%	**'correlation'**
%		Compute distance using correlation distance (1 - Pearson R) computed using all voxels (default)
% :Outputs:
%   **stats:**
%    Structure including:
%      gen_index  - sample estimate of generalization indices
%      bs_gen_index - bootstrap distribution for each generalization index
%      ste  - standard error of generalization index
%      Z  - Z score for each parameter
%      p  - p value for each parameter
%      sig - significance test at FDR q < .05
%
% :Examples:
% Example 1:
% ----------------------------------------------------------------------
% Download images from Kragel et al. 2018 Nature Neuroscience
% Kragel, P. A., Kano, M., Van Oudenhove, L., Ly, H. G., Dupont, P., Rubio, A., ? Wager, T. D. (2018). Generalizable representations of pain, cognitive control, and negative emotion in medial frontal cortex. Nature Neuroscience, 21(2), 283?289. doi:10.1038/s41593-017-0051-7
% 270 subject-level images systematically sampled from 18 studies across 3 domains
% [data_obj, names] = load_image_set('kragel18_alldata');
%
% Amygdala=select_atlas_subset(load_atlas('Canlab2018'),{'Amy'});
% masked_dat=apply_mask(data_obj,Amygdala);
%
% dsgn = condf2indic(ceil(data_obj.metadata_table.Studynumber/6));
% study = data_obj.metadata_table.Studynumber;
% stats = rsa_regression(masked_dat,dsgn,study);
%
% ..
%    Programmers' notes:
% List dates and changes here, and author of changes
%    01/08/2018, created function, Phil Kragel
% ..

%% convert design matrix into dissimilarity matrix
modelRDM=[];
for i=1:size(design,2) %for each specified grouping
   modelRDM(:,i)=pdist(design(:,i),'seuclidean'); %#ok<*AGROW> %compute squared euclidean distance
   modelRDM(:,i)=100000*modelRDM(:,i)/sum(modelRDM(:,i)); %normalize and scale
end

%supress distance warning - 0 value distances get removed...
warning('off','stats:pdist:ConstantPoints')

if any(strcmp([varargin{:}],'average_Euclidean'))
    brainRDM=pdist(nanmean(obj.dat)'); %distance is based on region average
elseif any(strcmp([varargin{:}],'euclidean'))
    brainRDM=pdist(obj.dat'); %euclidean

elseif any(strcmp([varargin{:}],'cosine'))
    brainRDM=pdist(obj.dat','cosine'); %euclidean

else %default to correlation
    brainRDM=pdist(obj.dat','correlation');
end
brainRDM(brainRDM<.00001)=NaN;
gen_index= glmfit([ones(length(modelRDM),1) double(modelRDM)],brainRDM','normal','constant','off');

num_it=1000;
bs_gen_index=zeros(num_it,size(gen_index,1));
parfor it=1:num_it
bs_gen_index(it,:) = random_resample_within_study(modelRDM,obj,study,varargin);
end

%compute stats from bootstrap distribution (normal approx)
b_ste = squeeze(nanstd(bs_gen_index));
b_mean = squeeze(nanmean(bs_gen_index));
b_ste(b_ste == 0) = Inf;
b_Z = b_mean ./ b_ste;
b_P = 2 * (1 - normcdf(abs(b_Z)));

% assign stats to output
stats.gen_index=gen_index;
stats.Z=b_Z;
stats.p=b_P;
stats.sig=b_P<FDR(b_P,.05);
stats.ste=b_ste;
stats.bs_gen_index=bs_gen_index;
stats.RDM=squareform(brainRDM);

end % main function

function beta = random_resample_within_study(modelRDM,obj,study,varargin)

bs_inds=zeros(size(obj.dat,2),1); %initialize

for i=1:max(study) %for each study

    study_inds=find(study==i); %find images that belong to this study
    bs_inds(study_inds)=study_inds(randi([1,length(study_inds)],1,length(study_inds))); %randomly resample with replacement

end
%supress distance warning - 0 value distances get removed...
warning('off','stats:pdist:ConstantPoints')

resampled_data=obj.dat(:,bs_inds)'; %random subsample and transpose for correlation

if any(strcmp([varargin{:}],'average_Euclidean'))
    brainRDM=pdist(nanmean(resampled_data')'); %distance is based on region average
elseif any(strcmp([varargin{:}],'euclidean'))
    brainRDM=pdist(resampled_data); %euclidean

elseif any(strcmp([varargin{:}],'cosine'))
    brainRDM=pdist(resampled_data,'cosine'); %euclidean

else %default to correlation
    brainRDM=pdist(resampled_data,'correlation');
end

brainRDM(brainRDM<.00001)=NaN; %exclude pairs that have exactly 0 distance

beta = glmfit([ones(length(modelRDM),1) double(modelRDM)],brainRDM','normal','constant','off'); %get parameter estimates using ols

end


text:
%%*******************************************************************
%% randinfsdp.m : creates random infeasible SDP problems with various block
%%                diagonal structures.
%%
%% [blk,Avec,C,b,X0,y0,Z0] =
%%      randinfsdp(dense_blk,sparse_blk,diag_blk,m,infeas,solve);
%%
%% E.g.
%%      randinfsdp([32 20],[10 5],100,10,infeas,solve);
%%
%%  dense_blk : for generating dense blocks, where
%%              dense_blk(i) is the dimension of the ith block
%%  sparse_blk: for generating a sparse block of small subblocks, where
%%              sparse_blk(i) is the size of the ith subblock.
%%  diag_blk: for generating a column vector of length specified by
%%            diag_blk (this corresponds to a diagonal block).
%%
%%  infeas = 1 if want primal infeasible pair of problems
%%         = 2 if want dual infeasible pair of problems
%%
%%  solve = 0 just to initialize
%%        = 1 if want to solve the problem.
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%******************************************************************

   function  [blk,Avec,C,b,X0,y0,Z0] = ...
                randinfsdp(dense_blk,sparse_blk,diag_blk,m,infeas,solve);

   if nargin < 6; solve = 0; end;
   if nargin < 5; error(' insufficient number of inputs '); end;
   if all(infeas-[1 2]);
      error(' infeas must be 1 or 2 ');
   end

   blk = [];
   if ~isempty(dense_blk);
      for k = 1:length(dense_blk);
          blk{k,1} = 's';  blk{k,2} = dense_blk(k);
      end;
   end;
   if ~isempty(sparse_blk);
      if size(sparse_blk,1) > size(sparse_blk,2); sparse_blk = sparse_blk'; end;
      tmp = size(blk,1);
      blk{tmp+1,1} = 's';
      blk{tmp+1,2} = sparse_blk;
   end;
   if ~isempty(diag_blk);
      tmp = size(blk,1);
      blk{tmp+1,1} = 'l';
      blk{tmp+1,2} = diag_blk;
   end;
   N = size(blk,1);

   if (infeas == 1),
%%
%% primal infeasible
%%
%% generate infeasibility certificate Zi and Z0
%%
   tmp_sp = sparse(sum(sparse_blk));
   for L = 1:2,
       T = [];
       if ~isempty(dense_blk);
          for k = 1:length(dense_blk);
              n = dense_blk(k);
              tmp = randn(n);
              tmp = tmp*tmp';
              tmp = 0.5*(tmp + tmp');
              %mineig = min(real(eig(tmp)));
              %if (mineig < 0); tmp = tmp - 1.1*mineig*eye(n); end;
              T{k,1} = tmp;
          end;
       end;
       if ~isempty(sparse_blk);
          for k = 1:length(sparse_blk);
              n = sparse_blk(k);
              pos = [sum(sparse_blk(1:k-1))+1 : sum(sparse_blk(1:k))];
              tmp = randn(n);
              tmp = tmp*tmp';
              tmp = 0.5*(tmp + tmp');
              %mineig = min(real(eig(tmp)));
              %if (mineig < 0); tmp = tmp - 1.1*mineig*eye(n); end;
              tmp_sp(pos,pos) = sparse(tmp);
          end;
          T{size(T,1)+1,1} = tmp_sp;
       end;
       if ~isempty(diag_blk);
          tmp = randn(diag_blk,1);
          tmp = tmp.*tmp;
          %mineig = min(tmp);
          %if (mineig < 0); tmp = tmp - 1.1*mineig*ones(diag_blk,1); end;
          T{size(T,1)+1,1} = tmp;
       end;
       if (L == 1);  Zi = T; else; Z0 = T; end;
   end;
%%
%% set up the matrices Ak and b
%%
   A  = cell(N,m);
   Ak_sp = sparse(sum(sparse_blk),sum(sparse_blk));
   for k = 1:m;
       Ak = [];
       if ~isempty(dense_blk);
          for j = 1:length(dense_blk);
              tmp = randn(dense_blk(j)); tmp = 0.5*(tmp+tmp');
              Ak{j,1} = tmp;
          end;
       end;
       if ~isempty(sparse_blk);
          for j = 1:length(sparse_blk);
              n = sparse_blk(j);
              tmp = randn(n); tmp = 0.5*(tmp+tmp');
              pos = [sum(sparse_blk(1:j-1))+1 : sum(sparse_blk(1:j))];
              Ak_sp(pos,pos) = tmp;
          end;
          Ak{size(Ak,1)+1,1} = Ak_sp;
       end;
       if ~isempty(diag_blk);
          Ak{size(Ak,1)+1,1} = randn(diag_blk,1);
       end;
       A(:,k) = Ak;
   end;

   y = ones(m,1);
   SAZm = ops(ops(Zi,'+',Asum(blk,A,y)),'/',m);
   SAZm = ops(ops(SAZm,'+',ops(SAZm,'transpose')),'*',0.5);
   yi  = randn(m,1);
   for k = 1:m,
      Ak = A(:,k);
      Ak = ops(ops(Ak,'-',SAZm),'/',yi(k));
      Ak = ops(ops(Ak,'+',ops(Ak,'transpose')),'*',0.5);
      A(:,k) = Ak;
   end;
   b = randn(m,1);
   if (b'*yi < 0), b = -b; end;
   y0  = randn(m,1);
   C = ops(Z0,'+',Asum(blk,A,y0));
   C = ops(ops(C,'+',ops(C,'transpose')),'*',0.5);
%%
%% (yi,Zi) is a primal infeasibility certificate
%%
   elseif (infeas == 2),
%%
%% dual infeasible
%%
%% generate infeasibility certificate Xi, positive definite X0, and C
%%
   tmp_sp_Xi  = sparse(sum(sparse_blk));
   tmp_sp_X0 = sparse(sum(sparse_blk));
   tmp_sp_C   = sparse(sum(sparse_blk));
   Xi = [];
   X0 = [];
   C  = [];
       if ~isempty(dense_blk);
          for k = 1:length(dense_blk);
              n = dense_blk(k);
              tmp = randn(n);
              tmp = tmp*tmp';
              tmp = 0.5*(tmp + tmp');
              %mineig = min(real(eig(tmp)));
              %if (mineig < 0); tmp = tmp - 1.1*mineig*eye(n); end;
              Xi{k,1} = tmp;
              tmp = randn(n);
              tmp = tmp*tmp';
              tmp = 0.5*(tmp + tmp');
              %mineig = min(real(eig(tmp)));
              %if (mineig < 0); tmp = tmp - 1.1*mineig*eye(n); end;
              X0{k,1} = tmp;
              tmp = randn(n); tmp = 0.5*(tmp + tmp');
              C{k,1} = tmp;
          end;
       end;
       if ~isempty(sparse_blk);
          for k = 1:length(sparse_blk);
              n = sparse_blk(k);
              pos = [sum(sparse_blk(1:k-1))+1 : sum(sparse_blk(1:k))];
              tmp = randn(n);
              tmp = tmp*tmp';
              tmp = 0.5*(tmp + tmp');
              %mineig = min(real(eig(tmp)));
              %if (mineig < 0); tmp = tmp - 1.1*mineig*eye(n); end;
              tmp_sp_Xi(pos,pos) = sparse(tmp);
              tmp = randn(n);
              tmp = tmp*tmp';
              tmp = 0.5*(tmp + tmp');
              %mineig = min(real(eig(tmp)));
              %if (mineig < 0); tmp = tmp - 1.1*mineig*eye(n); end;
              tmp_sp_X0(pos,pos) = sparse(tmp);
              tmp = randn(n);  tmp = 0.5*(tmp + tmp');
              tmp_sp_C(pos,pos) = sparse(tmp);
          end;
          Xi{size(Xi,1)+1,1} = tmp_sp_Xi;
          X0{size(X0,1)+1,1} = tmp_sp_X0;
          C{size(C,1)+1,1}   = tmp_sp_C;
       end;
       if ~isempty(diag_blk);
          n = diag_blk;
          tmp = randn(n,1);
          tmp = tmp.*tmp;
          %mineig = min(tmp);
          %if (mineig < 0); tmp = tmp - 1.1*mineig*ones(n,1); end;
          Xi{size(X0,1)+1,1} = tmp;
          tmp = randn(n,1);
          tmp = tmp.*tmp;
          %mineig = min(tmp);
          %if (mineig < 0); tmp = tmp - 1.1*mineig*ones(n,1); end;
          X0{size(X0,1)+1,1} = tmp;
          tmp = randn(n,1);
          C{size(C,1)+1,1} = tmp;
       end;
%%
%% set up the matrices Ak and b
%%
   trXX = blktrace(blk,Xi,Xi);
   A  = cell(N,m);
   Ak_sp = sparse(sum(sparse_blk),sum(sparse_blk));
   b = ones(m,1);
   for k = 1:m;
       Ak = [];
       if ~isempty(dense_blk);
          for j = 1:length(dense_blk);
              n = dense_blk(j);
              Aj  = randn(n);
              Ak{j,1} = (Aj + Aj')/2;
          end;
       end;
       if ~isempty(sparse_blk);
          for j = 1:length(sparse_blk);
              pos = [sum(sparse_blk(1:j-1))+1 : sum(sparse_blk(1:j))];
              n = sparse_blk(j);
              Aj  = randn(n);
              Ak_sp(pos,pos) = (Aj + Aj')/2;
          end;
          Ak{size(Ak,1)+1,1} = Ak_sp;
       end;
       if ~isempty(diag_blk);
          n = diag_blk;
          Ak{size(Ak,1)+1,1} = randn(n,1);
       end;
       trAX = blktrace(blk,Ak,Xi);
       Ak = ops(Ak,'-',ops(Xi,'*',(trAX/trXX)));
       Ak = ops(ops(Ak,'+',ops(Ak,'transpose')),'*',0.5);
       A(:,k) = Ak;
       b(k) = blktrace(blk,Ak,X0);
   end;
   if (blktrace(blk,C,Xi) > 0), C = ops(C,'*',-1); end;
%%
%% Xi is a dual infeasibility certificate
%%
   end;

%%
%% infeasible initial iterate
%%
   Avec = svec(blk,A,ones(size(blk,1),1));
   [X0,y0,Z0] = infeaspt(blk,Avec,C,b);
%%
   if solve;
      [obj,X,y,Z] = sqlp(blk,Avec,C,b,[],X0,y0,Z0);
   end;

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

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function nvar = p15_nvar ( option )

%*****************************************************************************80
%
%% P15_NVAR sets the number of variables for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, integer NVAR, the number of variables.
%
  nvar = 7;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ra=lpcar2ra(ar)
%LPCAR2RA Convert ar filter to inverse filter autocorrelation coefs. RA=(AR)

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

[nf,p1]=size(ar);
ra=zeros(nf,p1);
for i=1:p1
   ra(:,i)=sum(ar(:,1:p1+1-i).*ar(:,i:p1),2);
end

author: decouples


text:
% greBasedSegment.m
%
% APA, 3/20/2017

% directory containing all files
dirName = 'H:\Public\Aditya\mimExtensions\CERR_files_Sandra_contours_PC';
dirName = 'H:\Public\Aditya\mimExtensions\CERR_files_Sandra_contours_CT';
dirName = 'H:\Public\Aditya\mimExtensions\Atlas_Sanne\PC_cerr';
% dirName = 'H:\Public\Aditya\mimExtensions\Atlas_Sanne\CT_cropped_cerr';

% directory for writing the registered files (must have \ or / as last character)
registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_ROBINSON^HEATH_35487047\';
registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_ROBINSON^HEATH_35487047_CT\';
registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_MT160_PC\';
% registeredDir = 'H:\Public\Aditya\mimExtensions\registered_to_MT160_CT\';

dirS = dir(dirName);
dirS(1:2) = [];

% base scan file name
indBase = 3; %9 for Sandra's atlas, 3 for Sanne's
baseScan = fullfile(dirName,dirS(indBase).name);

% moving scan file names
indV = 1:length(dirS);
indV(indBase) = [];
movScanC = fullfile(dirName,{dirS(indV).name});

% registration callback
strNameToWarp = 'Parotid_L_SvD';
registerToAtlas(baseScan,movScanC,registeredDir,strNameToWarp)

% combine using the STAPLE and the GRE metric
regDirS = dir(registeredDir);
regDirS(1:2) = [];
regFilesC = strcat(registeredDir,{regDirS.name});
indBase = [3];
regFilesC(indBase) = [];
structNum = 2;
doseNum = 1;
scanNum = 1;
doseAllM = [];
strAllM = logical([]);
for i = 1:length(regFilesC)
    planC = loadPlanC(regFilesC{i},tempdir);
    indexS = planC{end};
    % Calculate the GRE metric
    baseScanNum = 1;
    movScanNum = 2;
    planC = calculateGRE(baseScanNum,movScanNum,planC);
    dose3M = getDoseOnCT(doseNum, scanNum, 'uniform', planC);
    str3M = getUniformStr(structNum,planC);
    strAllM(:,i) = str3M(:);
    doseAllM(:,i) = dose3M(:) .* str3M(:);
end

siz = size(str3M);

% STAPLE
numIter = 50;
confidence = 0.8;
numObservers = size(strAllM,2);
p = ones(1,numObservers)*0.999;
q = p;
[W,p,q] = staple(strAllM,confidence,p,q);
stapleStr3M = reshape(W > confidence,siz);
isUniform = 1;
scanNum = 1;
maskToCERRStructure(stapleStr3M,isUniform,scanNum,'STAPLE_80_pct_conf')

% Smooth contour
structNum = 2;
for slc = 1:length(planC{indexS.structures}(structNum).contour)
    for seg = 1:length(planC{indexS.structures}(structNum).contour(slc).segments)
        ptsM = planC{indexS.structures}(structNum).contour(slc).segments(seg).points;
        if isempty(ptsM)
            continue;
        end
        numPts = size(ptsM,1);
        intrvl = ceil(numPts*0.2/10);
        pts1M = spcrv(ptsM(1:intrvl:end,1:2)',3,100)';
        pts1M(:,3) = ptsM(1,3)*pts1M(:,1).^0;
        pts1M(end+1,:) = pts1M(1,:);
        planC{indexS.structures}(structNum).contour(slc).segments(seg).points = pts1M;
    end
end
reRasterAndUniformize

% GRE map
atlasGreV = sum(doseAllM,1) ./ sum(strAllM,1);
indToUse = atlasGreV < prctile(atlasGreV,30);
indToUse = 1:12;
weightedSegM = bsxfun(@times, strAllM(:,indToUse), (1./atlasGreV(indToUse)).^5);
weightedSegM(weightedSegM == 0) = NaN;
Wv = nansum(weightedSegM,2); % voxels weighted by GRE per registration

% Wv = nanmean(weightedSegM,2) < prctile(atlasGreV,50);
% % Wv = strAllM(:,2);
% segM = reshape(Wv,siz);
% maskToCERRStructure(segM,1,1,'GRE Weighted Majority')

% Combine GRE for each voxel
% gama = 1;
% invDoseAllM = 1./(doseAllM(:,indToUse)+eps);
% indZeroV = doseAllM(:,indToUse) > 0;
% invThrV = nanmean(doseAllM(:,indToUse),2);
% thrM = 1 ./ invThrV;
% thrM(invThrV < eps) = 0;
% thrV = mean(thrM,2)+1e10;
% % thr = mean(invDoseAllM(indZeroV)); % global thr
% % indZeroV = indZeroV & invDoseAllM > thrV;
% indZeroV = indZeroV & bsxfun(@le,invDoseAllM', thrV')';
% invDoseAllM(~indZeroV) = NaN;
% Wv = nansum(invDoseAllM.^gama , 2);
% % numMembersV = sum(strAllM,2);
% % Wv = Wv ./ numMembersV;

weightM = reshape(Wv,siz);
weightM(isnan(weightM)) = 0;
showIMDose(weightM,'ConsensusGRE',1);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function msi_file_get_index(fid, keyword, names)
%	MSI_FILE_GET_INDEX(fid, keyword)  Skip to a keyword in a file FID
%	and return the single valuea one based index array from a zero
%	based index array
%

frewind(fid);
test_line = sprintf('%s\t%%d', keyword);
index=0;
while feof(fid) == 0
    	line = fgetl(fid);
	matches = findstr(line, keyword);
	num = length(matches);
	if num > 0
	    % eat up keyword
	    [chopped,remainder] = strtok(line);
	    while(any(remainder));
		[chopped,remainder] = strtok(remainder, ',');
		names(index) = sscanf(chopped, '%s');
                index = index + 1;
	    end
	end
end

author: PatternRecognition

repo: OpenBMI


text:
function prob = appendMosekSDPdata(F_struc,K,prob)

prob.bardim = K.s;
prob.barc.subj = [];
prob.barc.subk = [];
prob.barc.subl = [];
prob.barc.val = [];
prob.bara.subi = [];
prob.bara.subj = [];
prob.bara.subk = [];
prob.bara.subl = [];
prob.bara.val = [];

C = F_struc(:,1);
A = -F_struc(:,2:end);

% -- Faster fix by Shahar
tops = [1 cumsum(K.s.^2)+1];
top = 1+K.f+K.l+sum(K.q);
[ii,jj,kk] = find(A(top:top + sum(K.s.^2)-1,:));
allcon = floor(interp1(tops,1:length(tops),ii,'linear'));
all_iilocal = ii-tops(allcon)'+1;
a = all_iilocal;
b = K.s(allcon);
allcol = ceil(a(:)./b(:))';
allrow = a(:)' - (allcol-1).*b(:)';
allvar = jj;
allval = kk;
% sort (for backward compatibility?)
[~,ind_sort] = sort(allcon);
allcon = allcon(ind_sort)';
allcol = allcol(ind_sort)';
allrow = allrow(ind_sort)';
allvar = allvar(ind_sort)';
allval = allval(ind_sort)';
% --

keep = find(allrow >= allcol);
allcol = allcol(keep);
allrow = allrow(keep);
allcon = allcon(keep);
allvar = allvar(keep);
allval = allval(keep);
prob.bara.subi = [prob.bara.subi allvar];
prob.bara.subj = [prob.bara.subj allcon];
prob.bara.subk = [prob.bara.subk allrow];
prob.bara.subl = [prob.bara.subl allcol];
prob.bara.val = [prob.bara.val allval];

for j = 1:length(K.s)
    n = K.s(j);
    Ci = C(top:top+n^2-1);
    Ci = tril(reshape(Ci,n,n));
    [k,l,val] = find(Ci);
    prob.barc.subj = [prob.barc.subj j*ones(1,length(k))];
    prob.barc.subk = [prob.barc.subk k(:)'];
    prob.barc.subl = [prob.barc.subl l(:)'];
    prob.barc.val = [prob.barc.val val(:)'];
    top = top + n^2;
end

author: shamilmamedov


text:
function MRSIStruct = setVoxelSize(MRSIStruct, dimension, value)
    switch(lower(dimension))
        case('x')
            MRSIStruct.voxelSize.x = value;
        case('y')
            MRSIStruct.voxelSize.y = value;
        case('z')
            MRSIStruct.voxelSize.z = value;
    end
end

author: CIC-methods

repo: FID-A

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


text:
function varargout = circleArcToPolyline(arc, N)
%CIRCLEARCTOPOLYLINE Convert a circle arc into a series of points
%
%   P = circleArcToPolyline(ARC, N);
%   convert the circle ARC into a series of N points.
%   ARC is given in the format: [XC YC R THETA1 DTHETA]
%   where XC and YC define the center of the circle, R its radius, THETA1
%   is the start of the arc and DTHETA is the angle extent of the arc. Both
%   angles are given in degrees.
%   N is the number of vertices of the resulting polyline, default is 65.
%
%   The result is a N-by-2 array containing coordinates of the N points.
%
%   [X Y] = circleArcToPolyline(ARC, N);
%   Return the result in two separate arrays with N lines and 1 column.
%
%
%   See also:
%   circles2d, circleToPolygon, drawCircle, drawPolygon
%
%
% ---------
% author : David Legland
% created the 22/05/2006.
% Copyright 2010 INRA - Cepia Software Platform.
%

% HISTORY
% 2011-03-30 use angles in degrees, add default value for N
% 2011-12-09 rename to circleArcToPolyline

% default value for N
if nargin < 2
    N = 65;
end

% vector of positions
t0 = deg2rad(arc(4));
t1 = t0 + deg2rad(arc(5));
t = linspace(t0, t1, N)';

% compute coordinates of vertices
x = arc(1) + arc(3) * cos(t);
y = arc(2) + arc(3) * sin(t);

% format output
if nargout <= 1
    varargout = {[x y]};
else
    varargout = {x, y};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [meg] = read_ctf_dat(filename)

% READ_CTF_DAT reads MEG data from an ascii format CTF file
%
% meg = read_ctf_dat(filename)
%
% returns a structure with the following fields:
%   meg.data        Nchans x Ntime
%   meg.time        1xNtime in miliseconds
%   meg.trigger     1xNtime with trigger values
%   meg.label       1xNchans cell-array with channel labels (string)

% Copyright (C) 2002, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(filename, 'r');

% read the sample number
line = fgetl(fid);
[tok, rem] = strtok(line, ':');
meg.sample = str2num(rem(2:end));

% read the time of each sample and convert to miliseconds
line = fgetl(fid);
[tok, rem] = strtok(line, ':');
meg.time = 1000*str2num(rem(2:end));

% read the trigger channel
line = fgetl(fid);
[tok, rem] = strtok(line, ':');
meg.trigger = str2num(rem(2:end));

% read the rest of the data
meg.data = [];
meg.label = {};
chan = 0;
while (1)
  line = fgetl(fid);
  if ~isempty(line) && line==-1
    % reached end of file
    break
  end
  [tok, rem] = strtok(line, ':');
  if ~isempty(rem)
    chan = chan + 1;
    meg.data(chan, :) = str2num(rem(2:end));
    meg.label{chan} = fliplr(deblank(fliplr(deblank(tok))));
  end
end

% convert to fT (?)
meg.data = meg.data * 1e15;

% apparently multiple copies of the data can be stored in the file
% representing the average, the standard deviation, etc.
% keep only the first part of the data -> average
tmp = find(diff(meg.time)<0);
if ~isempty(tmp)
  meg.data = meg.data(:,1:tmp(1));
  meg.time = meg.time(1:tmp(1));
  meg.trigger = meg.trigger(1:tmp(1));
  meg.sample = meg.sample(1:tmp(1));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ amin, amax ] = r8vec_range_2 ( n, a, amin, amax )

%*****************************************************************************80
%
%% R8VEC_RANGE_2 updates a range to include a new array.
%
%  Discussion:
%
%    Given a range AMIN to AMAX, and an array A, the routine will
%    decrease AMIN if necessary, or increase AMAX if necessary, so that
%    every entry of A is between AMIN and AMAX.
%
%    However, AMIN will not be increased, nor AMAX decreased.
%
%    This routine may be used to compute the maximum and minimum of a
%    collection of arrays one at a time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array.
%
%    Input, real AMIN, AMAX, the current legal range of values for A.
%
%    Output, real AMIN, AMAX, unchanged, or else "widened" so that all entries
%    of A are within the range.
%
  amax = max ( amax, max ( a(1:n) ) );
  amin = min ( amin, min ( a(1:n) ) );

  return
end


text:
classdef nncrop_wrapper < nntest
  methods (Test)

    function basicShrink(test)
      batchSize = 10 ;
      x1 = test.randn([10 10 3 batchSize]) ;
      x2 = test.randn([5 5 3 batchSize]) ;
      crop = [1 1] ;
      y = vl_nncrop_wrapper(x1, x2, crop) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      dzdx1 = vl_nncrop_wrapper(x1, x2, crop, dzdy) ;
      test.der(@(x1) vl_nncrop_wrapper(x1, x2, crop), ...
                               x1, dzdy, dzdx1, 1e-3*test.range) ;
    end
  end
end

repo: MFT


text:
%  Figure 10.27      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%  fig10_27.m is a script to generate Fig. 10.27, the
%  step response of the collocated
%  design for the satellite with PD compensation
clf;
m=[1, .1]; k0=[0, .091] ; d0=[0, .0036]; k1=[0, .4];
[f,g,h,j]=twomass(m,k0,d0); [f1,g,h,j] = twomass(m,k1,d0);
h1=[0, 0, 1, 0];

nc1=0.25*[2, 1];
dc1=[1/40, 1];

[ac,bc,cc,dc]=tf2ss(nc1, dc1);
[aol,bol,col,dol]= series(ac, bc,cc,dc,f,g,h1,j);
[acl]=aol-bol*col;
ccl2=[0*cc h]
[aol1,bol1,col1,dol1] = series(ac,bc,cc,dc,f1,g,h1,j);
acl1=aol1-bol1*col1;
t=0:.3:50;
sys1=ss(acl,bol,col,dol);
step(sys1,t);
hold on;
grid;
sys2=ss(acl1,bol1,col1,dol1);
step(sys2,t);
title('Closed-loop step response for D_5(s)G_{co}(s).')

author: Sable

repo: mcbench-benchmarks


text:
function [cerrFeatS,IBSIfeatS,pctDiffS] = compareRadiomicsWithIBSI1OrigImgWithInterp
% Function to compare radiomics features computed using CERR against
% the IBSI benchmark for configuration C.

%% 1. Calc. features for configuration 'C' using CERR
cerrFeatS = calcIBSIPhantomRadiomics('C');

%% Get IBSI benchmarks
ibsiConfigCResult = fullfile(fileparts(fileparts(getCERRPath)),...
            'Unit_Testing/data_for_cerr_tests/IBSI1_CT_phantom/IBSI_results_configC.mat');
temp = load(ibsiConfigCResult);
IBSIfeatS = temp.IBSIfeatS;

%% Compare results for each feature class
pctDiffS = struct();
%Shape features
shapeFeatC = fieldnames(IBSIfeatS.shapeS);
for n = 1:length(shapeFeatC)
    ibsiVal = IBSIfeatS.shapeS.(shapeFeatC{n});
    cerrVal = cerrFeatS.shapeS.(shapeFeatC{n});
    if ismember(shapeFeatC{n},{'volume','filledVolume'})
        cerrVal = cerrVal*1000; %cm^3 to mm^3
    elseif strcmp(shapeFeatC{n},'surfArea')
        cerrVal = cerrVal*100; %cm^2 to mm^2
    elseif ismember(shapeFeatC{n},{'majorAxis','minorAxis','leastAxis'}) || ...
            contains(shapeFeatC{n},'Diameter')
         cerrVal = cerrVal*10; %cm to mm
    elseif strcmp(shapeFeatC{n},'surfToVolRatio')
        cerrVal = cerrVal/10; %cm^-1 to mm^-1
    end
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.Shape.(shapeFeatC{n}) = pctDiff;
end

%IVH features
IVHfeatC = fieldnames(IBSIfeatS.Original.ivhFeaturesS);
for n = 1:length(IVHfeatC)
    ibsiVal = IBSIfeatS.Original.ivhFeaturesS.(IVHfeatC{n});
    cerrVal = cerrFeatS.Original.ivhFeaturesS.(IVHfeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.IVH.(IVHfeatC{n}) = pctDiff;
end

%First-order intensity features
firstOrdFeatC = fieldnames(IBSIfeatS.Original.firstOrderS);
for n = 1:length(firstOrdFeatC)
    ibsiVal = IBSIfeatS.Original.firstOrderS.(firstOrdFeatC{n});
    cerrVal = cerrFeatS.Original.firstOrderS.(firstOrdFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.FirstOrder.(firstOrdFeatC{n}) = pctDiff;
end

%GLCM features
%Avg
glcmFeatC = fieldnames(IBSIfeatS.Original.glcmFeatS.AvgS);
for n = 1:length(glcmFeatC)
    ibsiVal = IBSIfeatS.Original.glcmFeatS.AvgS.(glcmFeatC{n});
    cerrVal = cerrFeatS.Original.glcmFeatS.AvgS.(glcmFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.GLCM.Avg.(glcmFeatC{n}) = pctDiff;
end
%Merge
glcmFeatC = fieldnames(IBSIfeatS.Original.glcmFeatS.AvgS);
for n = 1:length(glcmFeatC)
    ibsiVal = IBSIfeatS.Original.glcmFeatS.AvgS.(glcmFeatC{n});
    cerrVal = cerrFeatS.Original.glcmFeatS.AvgS.(glcmFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.GLCM.Merge.(glcmFeatC{n}) = pctDiff;
end

%RLM features
%Avg
rlmFeatC = fieldnames(IBSIfeatS.Original.rlmFeatS.AvgS);
for n = 1:length(rlmFeatC)
    ibsiVal = IBSIfeatS.Original.rlmFeatS.AvgS.(rlmFeatC{n});
    cerrVal = cerrFeatS.Original.rlmFeatS.AvgS.(rlmFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.RLM.Avg.(rlmFeatC{n}) = pctDiff;
end
%Merge
rlmFeatC = fieldnames(IBSIfeatS.Original.rlmFeatS.CombS);
for n = 1:length(rlmFeatC)
    ibsiVal = IBSIfeatS.Original.rlmFeatS.CombS.(rlmFeatC{n});
    cerrVal = cerrFeatS.Original.rlmFeatS.CombS.(rlmFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.RLM.Merge.(rlmFeatC{n}) = pctDiff;
end

%NGTDM features
ngtdmFeatC = fieldnames(IBSIfeatS.Original.ngtdmFeatS);
for n = 1:length(ngtdmFeatC)
    ibsiVal = IBSIfeatS.Original.ngtdmFeatS.(ngtdmFeatC{n});
    cerrVal = cerrFeatS.Original.ngtdmFeatS.(ngtdmFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.NGTDM.(ngtdmFeatC{n}) = pctDiff;
end

%SZM features
szmFeatC = fieldnames(IBSIfeatS.Original.szmFeatS);
for n = 1:length(szmFeatC)
    ibsiVal = IBSIfeatS.Original.szmFeatS.(szmFeatC{n});
    cerrVal = cerrFeatS.Original.szmFeatS.(szmFeatC{n});
    pctDiff =(cerrVal-ibsiVal)*100/ibsiVal;
    pctDiffS.GLSZM.(szmFeatC{n}) = pctDiff;
end


text:
classdef TaizhouDataset < Datasets.CDDataset
    properties
        sampled logical=true;
    end
    methods
        function obj = TaizhouDataset(data_path, sampled)
            obj@Datasets.CDDataset(data_path);

            obj.refStr = 'Taizhou';

            if exist('sampled', 'var')
                obj.sampled = sampled;
            end

            obj.loaders.t1 = @Datasets.Loaders.tiffLoader;
            obj.loaders.t2 = @Datasets.Loaders.tiffLoader;

            if obj.sampled
                obj.loaders.ref = {@Datasets.Loaders.bmpLoader, @Datasets.Loaders.bmpLoader};
                obj.refList = {{fullfile(obj.dataPath, 'TaizhouChange_blackWhite.bmp'), ...
                    fullfile(obj.dataPath, 'TaizhouChange_blackWhite_unchange.bmp')}...
                    };
            else
                obj.loaders.ref = @Datasets.Loaders.bmpLoader;
                obj.refList = {fullfile(obj.dataPath, 'TaizhouChange_blackWhite.bmp')};
            end
        end
        function initFileSys(obj)
            obj.t1List = {fullfile(obj.dataPath, '2000TM.tif')};
            obj.t2List = {fullfile(obj.dataPath, '2003TM.tif')};
        end
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
function varargout = dirod(varargin)
% Date sorted folder listing
%
% USAGE:
%   dirod
%   dirod FILENAMEMASK
%   dirod('FILENAMEMASK')
%   LIST = dirod('FILENAMEMASK')
%
% INPUTS:
%   FILENAMEMASK - A file name mask, e.g. *, *.m or d*.m
%
% OUTPUTS:
%   LIST         - A directory structure ordered by date
%
% COMMENTS:
%   Similar to the windows command dir /od

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 6/7/2010

if isempty(varargin)
    arg = '*';
elseif length(varargin)==1
    arg = varargin{1};
else
    error('Too many input arguments.  0 or 1 input only');
end

dirData = dir(arg);
[~,index]=sort([dirData.datenum]);
dirData = dirData(index);

fileNameLen = length(dirData(1).name);
for i=2:length(dirData)
    fileNameLen = max(fileNameLen,length(dirData(i).name));
end

if nargout==0
    disp(' ');
    for i=1:length(dirData)
        disp([dirData(i).name repmat(' ',1,fileNameLen-length(dirData(i).name)+2) dirData(i).date]);
    end
    disp(' ');
else
    varargout{1} = dirData;
end

author: bashtage

repo: mfe-toolbox


text:
function lik = lik_negbinztr(varargin)
%LIK_NEGBINZTR  Create a zero-truncated Negative-binomial likelihood structure
%
%  Description
%    LIK = LIK_NEGBINZTR('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates zero-truncated Negative-binomial likelihood structure
%    in which the named parameters have the specified values. Any
%    unspecified parameters are set to default values.
%
%    Zero-truncated Negative-binomial can be used as a part of Hurdle model.
%
%    LIK = LIK_NEGBINZTR(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood structure with the named parameters
%    altered with the specified values.
%
%    Parameters for zero-truncated Negative-binomial likelihood [default]
%      disper       - dispersion parameter r [10]
%      disper_prior - prior for disper [prior_logunif]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                  __ n
%      p(y|f, z) = || i=1 [ (r/(r+mu_i))^r * gamma(r+y_i)
%                           / ( gamma(r)*gamma(y_i+1) )
%                           * (mu/(r+mu_i))^y_i / (1-(r/(r+mu_i))^r)]
%
%    where mu_i = z_i*exp(f_i) and r is the dispersion parameter. z
%    is a vector of expected mean and f the latent value vector
%    whose components are transformed to relative risk exp(f_i).
%    The last term (1-(r/(r+mu_i))^r) normalizes the truncated
%    distribution.
%
%    When using the zero-truncated Negbin likelihood you can give the
%    vector z as an extra parameter to each function that requires
%    also y. For example, you can call gp_optim as follows:
%      gp_optim(gp, x, y, 'z', z)
%    If z is not given or it is empty, then z_i=1 is used.
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%
% Copyright (c) 2007-2010 Jarno Vanhatalo & Jouni Hartikainen
% Copyright (c) 2010-2011 Aki Vehtari

% 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 = 'LIK_NEGBINZTR';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('disper',10, @(x) isscalar(x) && x>0);
  ip.addParamValue('disper_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Negbinztr';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Negbinztr')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('disper',ip.UsingDefaults)
    lik.disper = ip.Results.disper;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('disper_prior',ip.UsingDefaults)
    lik.p.disper=ip.Results.disper_prior;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_negbinztr_pak;
    lik.fh.unpak = @lik_negbinztr_unpak;
    lik.fh.lp = @lik_negbinztr_lp;
    lik.fh.lpg = @lik_negbinztr_lpg;
    lik.fh.ll = @lik_negbinztr_ll;
    lik.fh.llg = @lik_negbinztr_llg;
    lik.fh.llg2 = @lik_negbinztr_llg2;
    lik.fh.llg3 = @lik_negbinztr_llg3;
    lik.fh.tiltedMoments = @lik_negbinztr_tiltedMoments;
    lik.fh.tiltedMoments2 = @lik_negbinztr_tiltedMoments2;
    lik.fh.siteDeriv = @lik_negbinztr_siteDeriv;
    lik.fh.siteDeriv2 = @lik_negbinztr_siteDeriv2;
    lik.fh.upfact = @lik_negbinztr_upfact;
    lik.fh.predy = @lik_negbinztr_predy;
    lik.fh.predprcty = @lik_negbinztr_predprcty;
    lik.fh.invlink = @lik_negbinztr_invlink;
    lik.fh.recappend = @lik_negbinztr_recappend;
  end

end

function [w,s,h] = lik_negbinztr_pak(lik)
%LIK_NEGBINZTR_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_NEGBINZTR_PAK(LIK) takes a likelihood structure LIK and
%    combines the parameters into a single row vector W. This is a
%    mandatory subfunction used for example in energy and gradient
%    computations.
%
%       w = log(lik.disper)
%
%   See also
%   LIK_NEGBINZTR_UNPAK, GP_PAK

  w=[];s={};h=[];
  if ~isempty(lik.p.disper)
    w = log(lik.disper);
    s = [s; 'log(negbinztr.disper)'];
    h = [h 0];
    [wh,sh,hh] = lik.p.disper.fh.pak(lik.p.disper);
    w = [w wh];
    s = [s; sh];
    h = [h hh];
  end
end

function [lik, w] = lik_negbinztr_unpak(lik, w)
%LIK_NEGBINZTR_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    [LIK, W] = LIK_NEGBINZTR_UNPAK(W, LIK) takes a likelihood
%    structure LIK and extracts the parameters from the vector W
%    to the LIK structure. This is a mandatory subfunction used
%    for example in energy and gradient computations.
%
%   Assignment is inverse of
%       w = log(lik.disper)
%
%   See also
%   LIK_NEGBINZTR_PAK, GP_UNPAK

  if ~isempty(lik.p.disper)
    lik.disper = exp(w(1));
    w = w(2:end);
    [p, w] = lik.p.disper.fh.unpak(lik.p.disper, w);
    lik.p.disper = p;
  end
end

function lp = lik_negbinztr_lp(lik, varargin)
%LIK_NEGBINZTR_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_NEGBINZTR_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This subfunction
%    is needed when there are likelihood parameters.
%
%  See also
%    LIK_NEGBINZTR_LLG, LIK_NEGBINZTR_LLG3, LIK_NEGBINZTR_LLG2, GPLA_E

% If prior for dispersion parameter, add its contribution
  lp=0;
  if ~isempty(lik.p.disper)
    lp = lik.p.disper.fh.lp(lik.disper, lik.p.disper) +log(lik.disper);
  end

end

function lpg = lik_negbinztr_lpg(lik)
%LIK_NEGBINZTR_LPG  d log(prior)/dth of the likelihood
%                parameters th
%
%  Description
%    E = LIK_NEGBINZTR_LPG(LIK) takes a likelihood structure LIK and
%    returns d log(p(th))/dth, where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_NEGBINZTR_LLG, LIK_NEGBINZTR_LLG3, LIK_NEGBINZTR_LLG2, GPLA_G

  lpg=[];
  if ~isempty(lik.p.disper)
    % Evaluate the gprior with respect to disper
    ggs = lik.p.disper.fh.lpg(lik.disper, lik.p.disper);
    lpg = ggs(1).*lik.disper + 1;
    if length(ggs) > 1
      lpg = [lpg ggs(2:end)];
    end
  end
end

function ll = lik_negbinztr_ll(lik, y, f, z)
%LIK_NEGBINZTR_LL  Log likelihood
%
%  Description
%    LL = LIK_NEGBINZTR_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z, and
%    latent values F. Returns the log likelihood, log p(y|f,z).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also
%    LIK_NEGBINZTR_LLG, LIK_NEGBINZTR_LLG3, LIK_NEGBINZTR_LLG2, GPLA_E

  if numel(z)==0
    z=1;
  end

  r = lik.disper;
  mu = exp(f).*z;
  lp0=r.*(log(r) - log(r+mu));
  ll = sum(r.*(log(r) - log(r+mu)) + gammaln(r+y) - gammaln(r) - gammaln(y+1) + y.*(log(mu) - log(r+mu)) -log(1-exp(lp0)));
end

function llg = lik_negbinztr_llg(lik, y, f, param, z)
%LIK_NEGBINZTR_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_NEGBINZTR_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z and
%    latent values F. Returns the gradient of the log likelihood
%    with respect to PARAM. At the moment PARAM can be 'param' or
%    'latent'. This subfunction is needed when using Laplace
%    approximation or MCMC for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_NEGBINZTR_LL, LIK_NEGBINZTR_LLG2, LIK_NEGBINZTR_LLG3, GPLA_E

  if numel(z)==0
    z=1;
  end

  mu = exp(f).*z;
  r = lik.disper;
  switch param
    case 'param'
      % Derivative using the psi function
      llg = sum(1 + log(r./(r+mu)) - (r+y)./(r+mu) + psi(r + y) - psi(r));
      % add gradient of the normalization due to the truncation
      lp0=r.*(log(r) - log(r+mu));
      llg=llg-sum(1./(1 - exp(-lp0)).*(log(r./(mu + r)) - r./(mu + r) + 1));
      % correction for the log transformation
      llg = llg.*lik.disper;

% $$$             % Derivative using sum formulation
% $$$             llg = 0;
% $$$             for i1 = 1:length(y)
% $$$                 llg = llg + log(r/(r+mu(i1))) + 1 - (r+y(i1))/(r+mu(i1));
% $$$                 for i2 = 0:y(i1)-1
% $$$                     llg = llg + 1 / (i2 + r);
% $$$                 end
% $$$             end
% $$$             % correction for the log transformation
% $$$             llg = llg.*lik.disper;
    case 'latent'
      llg = y - (r+y).*mu./(r+mu);
      % add gradient of the normalization due to the truncation
      lp0=r.*(log(r) - log(r+mu));
      llg = llg -(1./(1-exp(-lp0)).*-r./(mu + r).*mu);
  end
end

function llg2 = lik_negbinztr_llg2(lik, y, f, param, z)
%LIK_NEGBINZTR_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_NEGBINZTR_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z, and
%    latent values F. Returns the Hessian of the log likelihood
%    with respect to PARAM. At the moment PARAM can be only
%    'latent'. LLG2 is a vector with diagonal elements of the
%    Hessian matrix (off diagonals are zero). This subfunction
%    is needed when using Laplace approximation or EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_NEGBINZTR_LL, LIK_NEGBINZTR_LLG, LIK_NEGBINZTR_LLG3, GPLA_E

  if numel(z)==0
    z=1;
  end

  mu = exp(f).*z;
  r = lik.disper;
  switch param
    case 'param'

    case 'latent'
      llg2 = - mu.*(r.^2 + y.*r)./(r+mu).^2;
      % add gradient of the normalization due to the truncation
      lp0=r.*(log(r) - log(r+mu));
      llg2=llg2+...
           (r.^2 + r.^2.*exp(-lp0).*(mu-1))./((mu + r).^2.*(exp(-lp0)-1).^2).*mu;
    case 'latent+param'
      llg2 = (y.*mu - mu.^2)./(r+mu).^2;
      % add gradient of the normalization due to the truncation
      lp0=r.*(log(r) - log(r+mu));
      llg2=llg2+(exp(lp0)./(exp(lp0) - 1).^2 .* (log(r) - log(mu + r) - r.*(1./(mu + r) - 1./r)) .* (-r./(mu + r)) -1./(1 - exp(-lp0)).*-mu./(mu + r).^2).*mu;

      % correction due to the log transformation
      llg2 = llg2.*lik.disper;

  end
end

function llg3 = lik_negbinztr_llg3(lik, y, f, param, z)
%LIK_NEGBINZTR_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_NEGBINZTR_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, incedence counts Y, expected counts Z and
%    latent values F and returns the third gradients of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. LLG3 is a vector with third gradients.  This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_NEGBINZTR_LL, LIK_NEGBINZTR_LLG, LIK_NEGBINZTR_LLG2, GPLA_E, GPLA_G

  if numel(z)==0
    z=1;
  end

  mu = exp(f).*z;
  r = lik.disper;
  switch param
    case 'param'

    case 'latent'
      llg3 = - mu.*(r.^2 + y.*r)./(r + mu).^2 + 2.*mu.^2.*(r.^2 + y.*r)./(r + mu).^3;
      % add gradient of the normalization due to the truncation
      lp0=r.*(log(r) - log(r+mu));
      llg3=llg3+ ...
           (exp(lp0).*(r.^2.*(r + r.*exp(2.*lp0)) + mu.^2.*r.^3 - mu.*r.^2.*(3.*r + exp(2.*lp0) + 1)) + exp(2.*lp0).*(mu.^2.*r.^3 - 2.*r.^3 + mu.*r.^2.*(3.*r + 2)))./((exp(lp0) - 1).^3.*(mu + r).^3).*mu;
    case 'latent2+param'
      llg3 = mu.*(y.*r - 2.*r.*mu - mu.*y)./(r+mu).^3;
      % add gradient of the normalization due to the truncation
      lp0=r.*(log(r) - log(r+mu));
      ip0=exp(-lp0);
      llg3=llg3+ ...
           (mu.*(2.*r + 2.*r.*ip0.*(mu - 1) + r.^2.*ip0.*(mu - 1).*(log(mu + r) - log(r) + r./(mu + r) - 1)))./((mu + r).^2.*(ip0 - 1).^2) - (2.*mu.*(r.^2 + r.^2.*ip0.*(mu - 1)))./((mu + r).^3.*(ip0 - 1).^2) - (2.*mu.*ip0.*(r.^2 + r.^2.*ip0.*(mu - 1)).*(log(mu + r) - log(r) + r./(mu + r) - 1))./((mu + r).^2.*(ip0 - 1).^3);
      % correction due to the log transformation
      llg3 = llg3.*lik.disper;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_negbinztr_tiltedMoments(lik, y, i1, sigma2_i, myy_i, z)
%LIK_NEGBINZTR_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_NEGBINZTR_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, incedence counts
%    Y, expected counts Z, index I and cavity variance S2 and
%    mean MYY. Returns the zeroth moment M_0, mean M_1 and
%    variance M_2 of the posterior marginal (see Rasmussen and
%    Williams (2006): Gaussian processes for Machine Learning,
%    page 55). This subfunction is needed when using EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

  if numel(z)==0
    avgE = ones(size(i1));
  else
    avgE = z(i1);
  end

  yy = y(i1);
  r = lik.disper;
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    if isscalar(sigma2_i)
      sigma2ii = sigma2_i;
    else
      sigma2ii = sigma2_i(i);
    end

    % get a function handle of an unnormalized tilted distribution
    % (likelihood * cavity = Negative-binomial * Gaussian)
    % and useful integration limits
    [tf,minf,maxf]=init_negbinztr_norm(yy(i),myy_i(i),sigma2ii,avgE(i),r,1);

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    if isnan(m_0)
      logM_0=NaN;
      return
    end
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigma2ii
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      if sigm2hati1(i) >= sigma2ii
        warning('lik_negbinztr_tilted_moments: sigm2hati1 >= sigm2_i');
        %sigm2hati1=sigm2_i-1e-9;
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [logM_0, m_1, sigm2hati1] = lik_negbinztr_tiltedMoments2(lik, y, i1, sigm2_i, myy_i, z, eta)
%LIK_NEGBINZTR_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_NEGBINZTR_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, incedence counts
%    Y, expected counts Z, index I and cavity variance S2 and
%    mean MYY. Returns the zeroth moment M_0, mean M_1 and
%    variance M_2 of the posterior marginal (see Rasmussen and
%    Williams (2006): Gaussian processes for Machine Learning,
%    page 55). This subfunction is needed when using EP for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

  if numel(z)==0
    avgE = ones(size(i1));
  else
    avgE = z(i1);
  end

  yy = y(i1);
  r = lik.disper;
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    % get a function handle of an unnormalized tilted distribution
    % (likelihood * cavity = Negative-binomial * Gaussian)
    % and useful integration limits
    [tf,minf,maxf]=init_negbinztr_norm(yy(i),myy_i(i),sigm2_i(i),avgE(i),r,eta);

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigm2_i(i)
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      if sigm2hati1(i) >= sigm2_i(i)
        warning('lik_negbinztr_tilted_moments: sigm2hati1 >= sigm2_i');
        %sigm2hati1=sigm2_i-1e-9;
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [g_i] = lik_negbinztr_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_NEGBINZTR_SITEDERIV  Evaluate the expectation of the gradient
%                      of the log likelihood term with respect
%                      to the likelihood parameters for EP
%
%  Description [M_0, M_1, M2] =
%    LIK_NEGBINZTR_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIK, incedence counts Y, expected
%    counts Z, index I and cavity variance S2 and mean MYY.
%    Returns E_f [d log p(y_i|f_i) /d a], where a is the
%    likelihood parameter and the expectation is over the
%    marginal posterior. This term is needed when evaluating the
%    gradients of the marginal likelihood estimate Z_EP with
%    respect to the likelihood parameters (see Seeger (2008):
%    Expectation propagation for exponential families). This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods and there are likelihood parameters.
%
%  See also
%    GPEP_G

  if numel(z)==0
    avgE = ones(size(i1));
  else
    avgE = z(i1);
  end

  yy = y(i1);
  r = lik.disper;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Negative-binomial * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_negbinztr_norm(yy,myy_i,sigm2_i,avgE,r, 1);
  % additionally get function handle for the derivative
  td = @deriv;

  % Integrate with quadgk
  [m_0, fhncnt] = quadgk(tf, minf, maxf);
  [g_i, fhncnt] = quadgk(@(f) td(f).*tf(f)./m_0, minf, maxf);
  g_i = g_i.*r;

  function g = deriv(f)
    mu = avgE.*exp(f);
    % Derivative using the psi function
    g = 1 + log(r./(r+mu)) - (r+yy)./(r+mu) + psi(r + yy) - psi(r);
    lp0=r.*(log(r) - log(r+mu));
    g = g -(1./(1 - exp(-lp0)).*(log(r./(mu + r)) - r./(mu + r) + 1));
  end
end

function [g_i] = lik_negbinztr_siteDeriv2(lik, y, i1, sigm2_i, myy_i, z, eta, lnZhat)
%LIK_NEGBINZTR_SITEDERIV  Evaluate the expectation of the gradient
%                      of the log likelihood term with respect
%                      to the likelihood parameters for EP
%
%  Description [M_0, M_1, M2] =
%    LIK_NEGBINZTR_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIK, incedence counts Y, expected
%    counts Z, index I and cavity variance S2 and mean MYY.
%    Returns E_f [d log p(y_i|f_i) /d a], where a is the
%    likelihood parameter and the expectation is over the
%    marginal posterior. This term is needed when evaluating the
%    gradients of the marginal likelihood estimate Z_EP with
%    respect to the likelihood parameters (see Seeger (2008):
%    Expectation propagation for exponential families). This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods and there are likelihood parameters.
%
%  See also
%    GPEP_G

  if numel(z)==0
    avgE = ones(size(i1));
  else
    avgE = z(i1);
  end

  yy = y(i1);
  r = lik.disper;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Negative-binomial * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_negbinztr_norm(yy,myy_i,sigm2_i,avgE,r, eta);
  % additionally get function handle for the derivative
  td = @deriv;

  % Integrate with quadgk
  [m_0, fhncnt] = quadgk(tf, minf, maxf);
  [g_i, fhncnt] = quadgk(@(f) td(f).*tf(f)./m_0, minf, maxf);
  g_i = g_i.*r;

  function g = deriv(f)
    mu = avgE.*exp(f);
    % Derivative using the psi function
    g = 1 + log(r./(r+mu)) - (r+yy)./(r+mu) + psi(r + yy) - psi(r);
    lp0=r.*(log(r) - log(r+mu));
    g = g -(1./(1 - exp(-lp0)).*(log(r./(mu + r)) - r./(mu + r) + 1));
%     g = eta.*g;
  end
end

function upfact = lik_negbinztr_upfact(gp, y, mu, ll, z)
  r = gp.lik.disper;
  sll = sqrt(ll);

  fh_e = @(f) negbinztr_pdf(y, exp(f).*z', r).*norm_pdf(f, mu, sll);
  EE = quadgk(fh_e, max(mu-6*sll,-30), min(mu+6*sll,30));

  fm = @(f) f.*negbinztr_pdf(y, exp(f).*z', r).*norm_pdf(f, mu, sll)./EE;
  mm  = quadgk(fm, max(mu-6*sll,-30), min(mu+6*sll,30));

  fV = @(f) (f - mm).^2.*negbinztr_pdf(y, exp(f).*z', r).*norm_pdf(f, mu, sll)./EE;
  Varp = quadgk(fV, max(mu-6*sll,-30), min(mu+6*sll,30));

  upfact = -(Varp - ll)./ll^2;
end

function [lpy, Ey, Vary] = lik_negbinztr_predy(lik, Ef, Varf, yt, zt)
%LIK_NEGBINZTR_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_NEGBINZTR_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the incedence counts YT, expected counts ZT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_NEGBINZTR_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  if numel(zt)==0
    zt=ones(size(Ef));
  end

  avgE = zt;
  r = lik.disper;

  lpy = zeros(size(Ef));
  Ey = zeros(size(Ef));
  EVary = zeros(size(Ef));
  VarEy = zeros(size(Ef));

  if nargout > 1
      % Evaluate Ey and Vary
      for i1=1:length(Ef)
        %%% With quadrature
        myy_i = Ef(i1);
        sigm_i = sqrt(Varf(i1));
        minf=myy_i-6*sigm_i;
        maxf=myy_i+6*sigm_i;

        F = @(f) exp(log(avgE(i1))+f+norm_lpdf(f,myy_i,sigm_i));
        Ey(i1) = quadgk(F,minf,maxf);

        F2 = @(f) exp(log(avgE(i1).*exp(f)+((avgE(i1).*exp(f)).^2/r))+norm_lpdf(f,myy_i,sigm_i));
        EVary(i1) = quadgk(F2,minf,maxf);

        F3 = @(f) exp(2*log(avgE(i1))+2*f+norm_lpdf(f,myy_i,sigm_i));
        VarEy(i1) = quadgk(F3,minf,maxf) - Ey(i1).^2;
      end
      Vary = EVary + VarEy;
  end

  % Evaluate the posterior predictive densities of the given observations
  lpy = zeros(length(yt),1);
  if (size(Ef,2) > 1) && (size(Ef,2) > 1) && size(yt,2) == 1
    % Approximate integral with sum of grid points when using corrected
    % marginal posterior pf
    for i1=1:length(yt)
      py = arrayfun(@(f) exp(lik.fh.ll(lik, yt(i1), f, zt(i1))), Ef(i1,:));
      pf = Varf(i1,:)./sum(Varf(i1,:));
      lpy(i1) = log(sum(py.*pf));
    end
  else
    for i1=1:length(yt)
      % get a function handle of the likelihood times posterior
      % (likelihood * posterior = Negative-binomial * Gaussian)
      % and useful integration limits
      [pdf,minf,maxf]=init_negbinztr_norm(...
        yt(i1),Ef(i1),Varf(i1),avgE(i1),r);
      % integrate over the f to get posterior predictive distribution
      lpy(i1) = log(quadgk(pdf, minf, maxf));
    end
  end
end

function [df,minf,maxf] = init_negbinztr_norm(yy,myy_i,sigm2_i,avgE,r, eta)
%INIT_NEGBINZTR_NORM
%
%  Description
%    Return function handle to a function evaluating
%    Negative-Binomial * Gaussian which is used for evaluating
%    (likelihood * cavity) or (likelihood * posterior) Return
%    also useful limits for integration. This is private function
%    for lik_negbinztr. This subfunction is needed by subfunctions
%    tiltedMoments, siteDeriv and predy.
%
%  See also
%    LIK_NEGBINZTR_TILTEDMOMENTS, LIK_NEGBINZTR_SITEDERIV,
%    LIK_NEGBINZTR_PREDY

% avoid repetitive evaluation of constant part
  ldconst = eta*(-gammaln(r)-gammaln(yy+1)+gammaln(r+yy))...
            - log(sigm2_i)/2 - log(2*pi)/2;
  % Create function handle for the function to be integrated
  df = @negbinztr_norm;
  % use log to avoid underflow, and derivates for faster search
  ld = @log_negbinztr_norm;
  ldg = @log_negbinztr_norm_g;
  ldg2 = @log_negbinztr_norm_g2;

  % Set the limits for integration
  % Negative-binomial likelihood is log-concave so the negbinztr_norm
  % function is unimodal, which makes things easier
  if yy==0
    % with yy==0, the mode of the likelihood is not defined
    % use the mode of the Gaussian (cavity or posterior) as a first guess
    modef = myy_i;
  else
    % use precision weighted mean of the Gaussian approximation
    % of the Negative-Binomial likelihood and Gaussian
    mu=log(yy/avgE);
    s2=(yy+r)./(yy.*r);
    modef = (myy_i/sigm2_i + mu/s2)/(1/sigm2_i + 1/s2);
  end
  % find the mode of the integrand using Newton iterations
  % few iterations is enough, since the first guess in the right direction
  niter=4;       % number of Newton iterations
  mindelta=1e-6; % tolerance in stopping Newton iterations
  for ni=1:niter
    g=ldg(modef);
    h=ldg2(modef);
    delta=-g/h;
    modef=modef+delta;
    if abs(delta)<mindelta
      break
    end
  end
  % integrand limits based on Gaussian approximation at mode
  modes=sqrt(-1/h);
  minf=modef-8*modes;
  maxf=modef+8*modes;
  modeld=ld(modef);
  iter=0;
  % check that density at end points is low enough
  lddiff=20; % min difference in log-density between mode and end-points
  minld=ld(minf);
  step=1;
  while minld<(modeld-lddiff) && minf<modef;
    % sometimes minf is too small
    minf=minf+step*modes;
    minld=ld(minf);
  end
  while minld>(modeld-lddiff)
    minf=minf-step*modes;
    minld=ld(minf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_negbinztr -> init_negbinztr_norm: ' ...
             'integration interval minimun not found ' ...
             'even after looking hard!'])
    end
  end
  maxld=ld(maxf);
  iter=0;
  step=1;
  while maxld>(modeld-lddiff)
    maxf=maxf+step*modes;
    maxld=ld(maxf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_negbinztr -> init_negbinztr_norm: ' ...
             'integration interval maximum not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = negbinztr_norm(f)
  % Negative-binomial * Gaussian
    mu = avgE.*exp(f);
    lp0=r.*(log(r) - log(r+mu));
    if lp0==0
      % exp(lp0)->1, that is, almost all the mass is in the zero part
      % approximate if yy=1, and give up if yy>1
      if yy==1
        integrand = exp(eta*(log(avgE)+f)...
                        -0.5*(f-myy_i).^2./sigm2_i -log(sigm2_i)/2 -log(2*pi)/2);
      else
        integrand = 0;
      end
    else
      integrand = exp(ldconst ...
                      +eta*(yy.*(log(mu)-log(r+mu))+r.*(log(r)-log(r+mu))) ...
                      -0.5*(f-myy_i).^2./sigm2_i ...
                      -eta*log(1-exp(lp0)));
    end
  end

  function log_int = log_negbinztr_norm(f)
  % log(Negative-binomial * Gaussian)
  % log_negbinztr_norm is used to avoid underflow when searching
  % integration interval
    mu = avgE.*exp(f);
    lp0=r.*(log(r) - log(r+mu));
    if lp0==0
      % exp(lp0)->1, that is, almost all the mass is in the zero part
      % approximate if yy=1, and give up if yy>1
      if yy==1
        log_int = eta*(log(avgE)+f) ...
                  -0.5*(f-myy_i).^2./sigm2_i - log(sigm2_i)/2 - log(2*pi)/2;
      else
        log_int=-Inf;
      end
    else
      log_int = ldconst...
                +eta*(yy.*(log(mu)-log(r+mu))+r.*(log(r)-log(r+mu)) -gammaln(r)-gammaln(yy+1)+gammaln(r+yy)) ...
                -0.5*(f-myy_i).^2./sigm2_i ...
                -eta*log(1-exp(lp0));
    end
  end

  function g = log_negbinztr_norm_g(f)
  % d/df log(Negative-binomial * Gaussian)
  % derivative of log_negbinztr_norm
    mu = avgE.*exp(f);
    lp0=r.*(log(r) - log(r+mu));
    if lp0==0
      % exp(lp0)->1, that is, almost all the mass is in the zero part
      % approximate if yy=1, and give up if yy>1
      g = eta+(myy_i - f)./sigm2_i;
    else
      g = -eta*(r.*(mu - yy))./(mu.*(mu + r)).*mu ...
          + (myy_i - f)./sigm2_i ...
          -eta*(1/(1 - exp(-lp0))*-r/(mu + r)*mu);
    end
  end

  function g2 = log_negbinztr_norm_g2(f)
  % d^2/df^2 log(Negative-binomial * Gaussian)
  % second derivate of log_negbinztr_norm
    mu = avgE.*exp(f);
    lp0=r.*(log(r) - log(r+mu));
    if lp0==0
      % exp(lp0)->1, that is, almost all the mass is in the zero part
      % approximate if yy=1, and give up if yy>1
      g2 = -1/sigm2_i;
    else
      g2 = -eta*(r*(r + yy))/(mu + r)^2.*mu ...
           -1/sigm2_i ...
           + eta*(r^2 + r^2*exp(-lp0)*(mu - 1))/((mu + r)^2*(exp(-lp0) - 1)^2)*mu;
    end
  end

end

function prctys = lik_negbinztr_predprcty(lik, Ef, Varf, zt, prcty)
%LIK_BINOMIAL_PREDPRCTY  Returns the percentiled of predictive density of y
%
%  Description
%    PRCTY = LIK_BINOMIAL_PREDPRCTY(LIK, EF, VARF YT, ZT)
%    Returns percentiles of the predictive density PY of YT, that is
%    This requires also the succes counts YT, numbers of trials ZT.
%    This subfunction is needed when using function gp_predprcty.
%
%  See also
%    GP_PREDPCTY

  if numel(zt)==0
    zt=ones(size(Ef));
  end

  opt=optimset('TolX',1e-7,'Display','off');
  nt=size(Ef,1);
  prctys = zeros(nt,numel(prcty));
  prcty=prcty/100;
  r = lik.disper;
  mu = zt.*exp(Ef);
  for i1=1:nt
    ci = sqrt(Varf(i1));
    for i2=1:numel(prcty)
      minf = floor(fminbnd(@(b) (quadgk(@(y) llvec(lik,y,Ef(i1)-1.96*ci,zt(i1)), 0, b)-prcty(i2)).^2,nbininv(prcty(i2), r, r./(r+zt(i1).*exp(Ef(i1))))-5,nbininv(prcty(i2), r, r./(r+zt(i1).*exp(Ef(i1))))+5,opt));
      if minf<0
        minf = 0;
      end
      maxf = floor(fminbnd(@(b) (quadgk(@(y) llvec(lik,y,Ef(i1)+1.96*ci,zt(i1)), 0, b)-prcty(i2)).^2,nbininv(prcty(i2), r, r./(r+zt(i1).*exp(Ef(i1))))-5,nbininv(prcty(i2), r, r./(r+zt(i1).*exp(Ef(i1))))+5,opt));
      if maxf<0
        maxf = 0;
      end
%       j=0;
%       figure;
%       for a=-2:0.1:50
%         j=j+1;
%         testi(j) = (quadgk(@(f) quadgk(@(y) llvec(lik,y,Ef(i1),zt(i1)), 0, a).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2;
%       end
%       plot(-2:0.1:50,testi)
%       if minf<maxf
%         set(gca,'XTick',[minf maxf])
%       else
%         set(gca,'XTick',[minf minf+1])
%       end
%       hold on;
%       a=floor(fminbnd(@(a) (quadgk(@(f) quadgk(@(y) llvec(lik,y,Ef(i1),zt(i1)), 0, a) ...
%              .*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2, minf, maxf,opt));
      a=floor(fminbnd(@(a) (quadgk(@(f) sum(llvec(lik,0:1:a,Ef(i1),zt(i1))) ...
             .*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4)-prcty(i2)).^2, minf, maxf,opt));
      if quadgk(@(f) sum(llvec(lik,0:1:a,Ef(i1),zt(i1))).*norm_pdf(f,Ef(i1),ci),Ef(i1)-6*ci,Ef(i1)+6*ci,'AbsTol',1e-4) < prcty(i2)
        a=a+1;
      end
      prctys(i1,i2)=a;
    end
  end

  function expll = llvec(lik,yt,f,z)
    % Compute vector of likelihoods of single predictions
    n = length(yt);
    if n>0
      for i=1:n
        expll(i) = exp(lik.fh.ll(lik, yt(i), f, z));
      end
    else
      expll = 0;
    end
  end

end

function mu = lik_negbinztr_invlink(lik, f, z)
%LIK_NEGBINZTR_INVLINK  Returns values of inverse link function
%
%  Description
%    MU = LIK_NEGBINZTR_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values MU of inverse link function.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_NEGBINZTR_LL, LIK_NEGBINZTR_PREDY

  if numel(z)==0
    z=1;
  end
  mu = bsxfun(@times,z,exp(f));
end

function reclik = lik_negbinztr_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = GPCF_NEGBINZTR_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK.  This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    % Initialize the record
    reclik.type = 'Negbinztr';

    % Initialize parameter
    reclik.disper = [];

    % Set the function handles
    reclik.fh.pak = @lik_negbinztr_pak;
    reclik.fh.unpak = @lik_negbinztr_unpak;
    reclik.fh.lp = @lik_negbinztr_lp;
    reclik.fh.lpg = @lik_negbinztr_lpg;
    reclik.fh.ll = @lik_negbinztr_ll;
    reclik.fh.llg = @lik_negbinztr_llg;
    reclik.fh.llg2 = @lik_negbinztr_llg2;
    reclik.fh.llg3 = @lik_negbinztr_llg3;
    reclik.fh.tiltedMoments = @lik_negbinztr_tiltedMoments;
    reclik.fh.tiltedMoments2 = @lik_negbinztr_tiltedMoments2;
    reclik.fh.siteDeriv = @lik_negbinztr_siteDeriv;
    reclik.fh.siteDeriv2 = @lik_negbinztr_siteDeriv2;
    reclik.fh.predy = @lik_negbinztr_predy;
    reclik.fh.predprcty = @lik_negbinztr_predprcty;
    reclik.fh.invlink = @lik_negbinztr_invlink;
    reclik.fh.recappend = @lik_negbinztr_recappend;
    reclik.p=[];
    reclik.p.disper=[];
    if ~isempty(ri.p.disper)
      reclik.p.disper = ri.p.disper;
    end
  else

    % Append to the record
    reclik.disper(ri,:)=lik.disper;
    if ~isempty(lik.p.disper)
      reclik.p.disper = lik.p.disper.fh.recappend(reclik.p.disper, ri, lik.p.disper);
    end
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
x=[1 2 3 4 5 6 7 8 9]
pause
find(x<6)
pause
y=[1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
pause
find(y<6)'
pause
[l,c,v]=find(y<6);
pause
l'
pause
c'
pause
v'


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

% Useful variables
[m, n] = size(X);

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu = (1 / m) * sum(X);
sigma2 = (1 / m) * sum((X - mu) .^ 2);

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

end

repo: machine-learning


text:
% in LMM / HGLMM, need to multiply by 10 to avoid numerical issue
LMM_HGLMM_MULT_FACTOR = 10;

global init_flag;
if isempty(init_flag)

    % location of this script
    pathstr = fileparts( mfilename('fullpath') );

    addpath([pathstr '/../utilities']);
    addpath([pathstr '/../FastICA_25']);
    addpath([pathstr '/../cca']);

    if isunix
        addpath([pathstr '/LMM_linux']);
        addpath([pathstr '/HGLMM_linux']);
    else
        addpath([pathstr '/LMM_win']);
        addpath([pathstr '/HGLMM_win']);
    end

    run([pathstr '/../vlfeat-0.9.18/toolbox/vl_setup']);

    clear pathstr;

    init_flag = true;
end

author: layumi

repo: Image-Text-Embedding


text:
function imdb = getMLIVEDatabase_train(Dir)

Scores1 = load([Dir, '\Part 1\Scores.mat']);
imdb.dmos1 = Scores1.DMOSscores;
Scores2 = load([Dir, '\Part 2\Scores.mat']);
imdb.dmos2 = Scores2.DMOSscores;

Imagelists1 = load([Dir, '\Part 1\Imagelists.mat']);
refnames1 = Imagelists1.refimgs;
distimgs1 = Imagelists1.distimgs;
ref4dist1 = Imagelists1.ref4dist;
Imagelists2 = load([Dir, '\Part 2\Imagelists.mat']);
refnames2 = Imagelists2.refimgs;
distimgs2 = Imagelists2.distimgs;
ref4dist2 = Imagelists2.ref4dist;

imdb.dmos = [imdb.dmos1,imdb.dmos2];

imdb.dir_dist1 = cell(1,225);
for i = 1:225
    imdb.dir_dist1{1,i} = fullfile('Part 1\blurjpeg',distimgs1{i,1});
end
imdb.dir_dist2 = cell(1,225);
for i = 1:225
    imdb.dir_dist2{1,i} = fullfile('Part 2\blurnoise',distimgs2{i,1});
end

imdb.imgDir =  cat(2,imdb.dir_dist1,imdb.dir_dist2);
imdb.dataset = 'MLIVE';
imdb.filenum = 450;

imdb.dmosall = imdb.dmos;
imdb.refnum = 15;
for i = 1:15
    imdb.refname{i} = num2str(i);
end
for i = 1:225
    imdb.refnames_all{i} = num2str(ref4dist1(i));
end
for i = 226:450
    imdb.refnames_all{i} = num2str(ref4dist2(i-225));
end

sel = randperm(15);
% sel = 1:15;
train_path = [];
train_dmos = [];
for i = 1:12
    train_sel = strcmpi(imdb.refname(sel(i)),imdb.refnames_all);
    train_sel = find(train_sel == 1);
    train_path = [train_path, imdb.imgDir(train_sel)];
    train_dmos = [train_dmos,imdb.dmos(train_sel)];
end

test_path = [];
test_dmos = [];
for i = 13:15
    test_sel = strcmpi(imdb.refname(sel(i)),imdb.refnames_all);
    test_sel = find(test_sel == 1);
    test_path = [test_path, imdb.imgDir(test_sel)];
    test_dmos = [test_dmos,imdb.dmos(test_sel)];
end

imdb.images.id = 1:450 ;
imdb.images.set = [ones(1,size(train_path,2)),2*ones(1,size(test_path,2))];
imdb.images.label = [train_dmos,test_dmos];
imdb.classes.description = {'MLIVE'};
imdb.images.name = [train_path,test_path] ;
imdb.imageDir = Dir;

repo: DBCNN


text:
% LOADODS: Loads data from an open document spreadsheet (ods) file into a cell
%          array
%
% SYNTAX: data = loadods(filename,options)
%
% INPUTS:
%   filename:   String representing the path and name of the ods file
%   options:    Structural variable containing optional arguments
%
% OUTPUTS:
%   data:       MxN cell array containing data from the spreadsheet
%
% OPTIONS:
%   sheet_name: Name of a specific sheet to load [defaults to the first sheet]
%   blank:      Value to store for blank spreadsheet cells [NaN]
%
% VERSION: 1.0
%

%Antonio Quarta - sgheppy88@gmal.com
%
% based on Alex Marten loadods.m function

function data = loadods(filename,opt1,opt2)
	% Set default options
            if nargin<2, sheetn = '';bchar=nan; end
            if nargin==2
                if  ~ischar(opt1)
                    sheetn = opt1;
                    bchar=nan;
                else
                    sheetn='';
                    bchar=opt2;
                end
            end
            if nargin==3
                sheetn = opt1;
                bchar=opt2;
            end

	% Check for the existence of the ods file
	if exist(filename)~=2
		error('The file does not exist');
	end

	% Create a temporary directory to unzip the ods file into
	if ~mkdir(tempdir,filename)
		error('Permission error');
	end
	dir_temp = strcat(tempdir,filename);

	% Unzip the contents of the ods file
	unzip(filename,dir_temp);

	% Load the XML file containing the spreadsheet data
	try
 		XMLfile = xmlread(strcat(dir_temp,'/content.xml'));
	catch
	   error('Unable to read the spreadsheet data');
	end

	% Parse down to the <office:spreadsheet> node

	nodes = XMLfile.getChildNodes;
        nodes2=nodes;
        i=0;
        trovato=0;
        while (~trovato && i < nodes.getLength )
            if (nodes.item(i).getNodeName ==  'office:document-content')
                node=nodes.item(i);
                nodes = node.getChildNodes;
                i=0;
            end
            if (nodes.item(i).getNodeName ==  'office:body')
                node=nodes.item(i);
                nodes = node.getChildNodes;
                i=0;
            end
            if (nodes.item(i).getNodeName ==  'office:spreadsheet')
                node=nodes.item(i);
                nodes = node.getChildNodes;
                trovato=1;
            end
            i=i+1;
        end

	% Find the requested sheet by name or default to the first sheet
	if ~isempty(sheetn)
		numSheets = nodes.getLength;
		for count = 1:numSheets
			sheet = nodes.item(count-1);
			if strcmp(get_attribute(sheet,'table:name'),sheetn)
				break
			end
		end
	else
            i=0;
            trovato=0;
            while (~trovato && i < nodes.getLength)
                if( nodes.item(i).getNodeName == 'table:table')
                    sheet = nodes.item(i);

                trovato=1;
               end
               i=i+1;
            end
	end

	% Get the begin of rows
 	nodes = sheet.getChildNodes;
	num_nodes = nodes.getLength;
	num_cols = 0;
	        for count = 1:num_nodes
		node = nodes.item(count-1);
            	if strcmp(char(node.getNodeName),'table:table-row')
                    count = count-2;
			break
		end
	end

	% Get the number of rows
	num_rows = num_nodes-count-1;

	% Initialize memory for the data
	data = cell(1,1);

	% Extract the data for the sheet
        rn=0;
        numrep=0;
	for row_num = 1:num_rows
            if ( ~strcmp(nodes.item(count+row_num).getNodeName,'#text') ...
                 && exist_validvalue( nodes.item(count+row_num)) )
                rn=rn+1;
                row = nodes.item(count+row_num);
                cols = row.getChildNodes;
                col_num = 0;
                num_items = cols.getLength-1;

                for item_num = 0:num_items
                    if ( ~strcmp(cols.item(item_num).getNodeName,'#text') ...
                         && exist_attribute(cols.item(item_num), ...
                                            'office:value') )
                        col = cols.item(item_num);
                        if (exist_attribute(col, ...
                                            'table:number-columns-repeated'))
                            numrep = str2num(get_attribute(col,['table:' ...
                                                'number-columns-repeated']));
                        end
                        value_type = get_attribute(col,'office:value-type');
                        if strcmp(value_type,'string')
                            temp = col.getChildNodes;
                            temp = temp.item(0);
                            temp = temp.getChildNodes;
                            temp = temp.item(0);
                            if any(strcmp(methods(temp),'getData'))
                                    value = char(temp.getData);
                            else
                                value = options.blank;
                            end
                        elseif strcmp(value_type,'float')
                            value = str2num(get_attribute(col,'office:value'));
                        else
                            value = bchar;
                        end
                        if numrep
                            for i=1:numrep
                            col_num = col_num+1;
                            data{rn,col_num} = value;
                            end
                            numrep=0;
                        else
                            col_num = col_num+1;
                            data{rn,col_num} = value;
                        end
                    end
                end
            end

 end

	% Remove the temporary files
	if ~rmdir(dir_temp,'s')
		warning('Temporary files could not be removed');
	end

end

% Returns the value of the attribute_name in node
function attribute = get_attribute(node,attribute_name)
	attribute = [];
	if node.hasAttributes
		attributes = node.getAttributes;
		num_attributes = attributes.getLength;
		for count = 1:num_attributes
			item = attributes.item(count-1);
		 	if strcmp(char(item.getName),attribute_name)
				attribute = char(item.getValue);
			end
		end
	end
end
function attribute = exist_attribute(node,attribute_name)
    if node.hasAttributes
        attributes=node.getAttributes;
        num_attributes=attributes.getLength;
        for count = 1:num_attributes
            item= attributes.item(count-1);
            if(strcmp(char(item.getName),attribute_name))
                attribute=1;
                return ;
            end
        end
        attribute=0;
        return ;
    end
    attribute=0;
    return ;
end
% It will return 1 if row contains any cells with a valid value.
% node is valid row
function result = exist_validvalue( node )
    cols = node.getChildNodes;
    col_num = 0;
    num_items = cols.getLength-1;
    for item_num = 0:num_items
        if (  ~strcmp(cols.item(item_num).getNodeName,'#text') ...
                         && exist_attribute(cols.item(item_num), ...
                                            'office:value') )
            result = 1 ;
            return ;
        end
    end
    result = 0;
    return ;
end