text:
% PRTBRVMVB - PRT BRV Multivariate Binary Observation Model

classdef prtBrvMvb < prtBrv & prtBrvVbOnline & prtBrvVbMembershipModel & prtBrvVbOnlineMembershipModel

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Properties required by prtAction
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties (SetAccess = private)
        name = 'Multi-varite Binary Bayesian Random Variable';
        nameAbbreviation = 'BRVMVB';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrv
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        function self = estimateParameters(self, x)
            self = conjugateUpdate(self, self, x);
        end

        function y = predictivePdf(self, x)
            %%%% FIXME
            % The true predictive here is a product of beta-binomials
            % Since that isn't implemented yet we use the average
            % variational loglikelihood

            y = conjugateVariationalAverageLogLikelihood(self, x);
        end

        function val = getNumDimensions(self)
            val = length(self.model.countOfOnes);
        end

        function self = initialize(self, x)
            x = self.parseInputData(x);
            if ~self.model.isValid
                self.model = self.model.defaultParameters(size(x,2));
            end
        end

        % Optional methods
        %------------------------------------------------------------------
        function kld = conjugateKld(obj, priorObj)
            betaKlds = zeros(1,length(obj.model.countOfOnes));
            for iDim = 1:length(obj.model.countOfOnes)
                betaKlds(iDim) = prtRvUtilDirichletKld([obj.model.countOfZeros(iDim) obj.model.countOfOnes(iDim)],[priorObj.model.countOfZeros(iDim) priorObj.model.countOfOnes(iDim)]);
            end

            kld = sum(betaKlds);
        end

        function s = posteriorMeanStruct(obj)
            s.probabilities = obj.model.countOfOnes./(obj.model.countOfOnes + obj.model.countOfZeros);
        end

        function plotCollection(objs,colors)

            nComponents = length(objs);

            if nargin < 2
                colors = prtPlotUtilClassColors(nComponents);
            end

            countOfOnesMat = cell2mat(arrayfun(@(c)c.model.countOfOnes,objs,'uniformOutput',false));
            countOfZerosMat = cell2mat(arrayfun(@(c)c.model.countOfZeros,objs,'uniformOutput',false));

            probMat = countOfOnesMat ./ (countOfOnesMat + countOfZerosMat);

            h = plot(probMat');

            for iLine = 1:length(h)
                set(h(iLine),'color',colors(iLine,:));
            end
            title('Posterior Mean Probabilities')
            ylim([0 1])
            xlim([1 size(probMat,2)])

        end

        function val = plotLimits(self)
            val = [1 self.nDimensions 0 1];
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvVb
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function [self, training] = vbBatch(self,x)
            % Since we are purely conjugate we actually don't need vbBatch
            % However we must implement it.
            self = conjugateUpdate(self,x);
            training = struct([]);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvMembershipModel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        function [phiMat, priorObjs] = collectionInitialize(objs, priorObjs, X) % Vector of objects

            learningInitialMembershipFactor = objs(1).initFudgeFactor;

            [classMeans,kmMembership] = prtUtilKmeans(X,length(objs),'handleEmptyClusters','random','distanceMetricFn',@prtDistanceHamming,'maxIterations',100,'logicalMeans',true); %#ok<ASGLU>

            phiMat = zeros(size(X,1),length(objs));
            for iComp = 1:length(objs)
                phiMat(kmMembership == iComp, iComp) = learningInitialMembershipFactor;
            end
            phiMat(phiMat==0) = (1-learningInitialMembershipFactor)./(length(objs)-1);

            % We should normalize this just in case the
            % learningInitialMembershipFactor was set poorly
            phiMat = bsxfun(@rdivide,phiMat,sum(phiMat,2));

        end

        function obj = weightedConjugateUpdate(obj, priorObj, x, weights)

            x = obj.parseInputData(x);
            priorObj = priorObj.initialize(x);

            if nargin < 4 || isempty(weights)
                weights = ones(size(x,1),1);
            end

            obj.model.countOfOnes = sum(bsxfun(@times,x,weights)) + priorObj.model.countOfOnes;

            obj.model.countOfZeros = sum(bsxfun(@times,~x,weights)) + priorObj.model.countOfZeros;

        end

        function self = conjugateUpdate(self, prior, x)
            x = self.parseInputData(x);

            self = weightedConjugateUpdate(self, prior, x, ones(size(x,1),1));
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvVbMembershipModel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function y = conjugateVariationalAverageLogLikelihood(obj, x)

            [logProb0, logProb1] = obj.expectedValueLogProbabilities;

            y = sum(bsxfun(@times,x,logProb1) + bsxfun(@times,~x,logProb0),2);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvVbOnlineMembershipModel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function obj = vbOnlineInitialize(obj, x) %#ok<INUSD>

            randDraw = (rand(1,obj.nDimensions)>0.5);

            obj.model.countOfOnes = randDraw+0.5;
            obj.model.countOfZeros  = (1-randDraw) + 0.5;
        end

        function [self, training] = vbOnlineUpdate(self, priorObj, x, lambda, D, prevObj)
            x = self.parseInputData(x);
            [self, training] = vbOnlineWeightedUpdate(self, priorObj, x, ones(size(x,1),1), lambda, D, prevObj);
        end

        function [obj, training] = vbOnlineWeightedUpdate(obj, priorObj, x, weights, lambda, D, prevObj)
            x = obj.parseInputData(x);

            S = size(x,1);

            if nargin < 4 || isempty(weights)
                weights = ones(size(x,1),1);
            end

            obj.model.countOfOnes = (D/S*sum(bsxfun(@times,x,weights)) + priorObj.model.countOfOnes)*lambda + prevObj.model.countOfOnes*(1-lambda);
            obj.model.countOfZeros = (D/S*sum(bsxfun(@times,~x,weights)) + priorObj.model.countOfZeros)*lambda + prevObj.model.countOfZeros*(1-lambda);

            training = struct([]);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Properties for prtBrvMvb use
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties
        model = prtBrvMvbHierarchy;
    end

    properties (Hidden)
        initFudgeFactor = 1; % Between zero and 1, probably > 0.9
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvMvb use
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function self = prtBrvMvb(varargin)
            if nargin < 1
                return
            end
            self = constructorInputParse(self,varargin{:});
        end

        function [logProb0, logProb1] = expectedValueLogProbabilities(obj)
            psi1 = psi(obj.model.countOfOnes);
            psi0 = psi(obj.model.countOfZeros);

            psiSum = psi(obj.model.countOfOnes+obj.model.countOfZeros);

            logProb1 = psi1 - psiSum;
            logProb0 = psi0 - psiSum;
        end

        function model = modelDraw(obj)
            model.probabilities = zeros(1,length(obj.model.countOfOnes));
            for iDim = 1:length(obj.model.countOfOnes)
                model.probabilities(iDim) = prtRvUtilDirichletRnd([obj.model.countOfZeros(iDim) obj.model.countOfOnes(iDim)]);
            end
        end
    end

    methods (Hidden)
        function x = parseInputData(self,x) %#ok<MANU>
            if isnumeric(x) || islogical(x)
                return
            elseif prtUtilIsSubClass(class(x),'prtDataSetBase')
                x = x.getObservations();
            else
                error('prt:prtBrvMvb:parseInputData','prtBrvMvb requires a prtDataSet or a numeric 2-D matrix');
            end
        end
    end
end


text:
% This is used to set the path to visible images
%
% Please change the path in line 16 to your own path before running
%
% Author:Xingchen Zhang, Ping Ye, Gang Xiao
% Contact: xingchen.zhang@imperial.ac.uk
%
% If you use this code, please site the following paper:
%
% X. Zhang, P. Ye, G. Xiao. VIFB:A Visible and Infrared Benchmark. In
% Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern
% Recognition Workshops, 2020.

function imgs=configImgsIR

    path = 'Your own path\VIFB\input\';

    img_IR = {struct('name','carLight','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','carShadow','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','carWhite','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','elecbike','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','fight','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','kettle','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','labMan','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','man','path',strcat(path,'IR\'),'ext','jpg'),...
     struct('name','manCall','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','manCar','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','manlight','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','manWalking','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','manwithbag','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','nightCar','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','peopleshadow','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','running','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','snow','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','tricycle','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','walking','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','walking2','path',strcat(path,'IR\'),'ext','jpg'),...
    struct('name','walkingnight','path',strcat(path,'IR\'),'ext','jpg'),...
    };

imgs=[img_IR];

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [vertices, faces] = mergeMeshVertices(vertices, faces, vertexInds, varargin)
%MERGEMESHVERTICES Merge two vertices and removes eventual degenerated faces.
%
%   output = mergeMeshVertices(input)
%
%   Example
%   mergeMeshVertices
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2019-01-31,    using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019 INRA - Cepia Software Platform.

newPos = vertices(vertexInds(1), :);
if nargin > 3
    newPos = varargin{1};
end

vertices(vertexInds(1), :) = newPos;
vertices(vertexInds(2:end), :) = NaN;

% replace face-vertex indices by index of first vertex
faces(ismember(faces, vertexInds)) = vertexInds(1);

% need to check existence of degenerated faces with same vertex twice
nFaces = size(faces, 1);
dgnFaces = false(nFaces, 1);
dims = [1 2;1 3;2 3];
for i = 1:3
    dgnFaces = dgnFaces | faces(:,dims(i,1)) == faces(:,dims(i,2));
end

% remove degenerated faces
faces(dgnFaces, :) = [];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function h = specgram2(s, T, varargin)
   %SPECGRAM2 - plots spectrograms of waveforms with waveform along top
   %  h = specgram2(spectralobject, waveforms) generates a spectrogram from
   %  the waveform(s), overwriting the current figure.  The waveform will be
   %  displayed along the top of the spectrogram. The return value is a handle
   %  to the spectrogram, and is optional.
   %
   %  The spectrograms will be created in the same shape as the passed
   %  waveforms.  ie, if W is a 2x3 matrix of waveforms, then
   %  specgram2(spectralobject,W) will generate a 2x3 plot of spectra.
   %
   %  Many additional behaviors can be modified through the passing of
   %  additional parameters, as listed further below.  These parameters are
   %  always passed in pairs, as such:
   %
   %  specgram2(spectralobject, waveforms,'PARAM1',VALUE1,...,'PARAMn',VALUEn)
   %    Any number of these parameters may be passed to specgram2.
   %
   %  specgram2(..., 'axis', AXIS_HANDLE)
   %    Specify the axis AXIS_HANDLE within which the spectrogram will be
   %    generated.  The boundary of the axis becomes the boundary for the
   %    entire spectra plot.  For a matrix of waveforms, this area is
   %    subdivided into NxM subplots, where N and M are the size of the
   %    waveform matrix.
   %
   %  specgram2(..., 'xunit', XUNIT)
   %    Spedifies the x-unit scale to be used with the spectrogram.  The
   %    default unit is 'seconds'.
   %    valid xunits:
   %     'seconds','minutes','hours','days','doy' (day of year),and 'date'
   %
   %  specgram2(..., 'colormap', ALTERNATEMAP)
   %    Instead of using the default colormap, any colormap may be used.  An
   %    alternate way of setting the global map is by using the SETMAP
   %    function.  ALTERNATEMAP will either be a name (eg. grayscale) or an
   %    Nx3 numeric. Type HELP GRAPH3D to see additional useful colormaps.
   %
   %
   %  specgram2(..., 'colorbar', COLORBAR_OPTION)
   %    Generates a spectrogram from the waveform and uses a specific map
   %    valid COLORBAR_OPTION values: 'horiz' (default),'vert','none',
   %      'HORIZ' places a single colorbar below all plots
   %      'VERT' places a single colorbar to the right of all plots
   %      'NONE' supresses the colorbar placement
   %
   %  specgram2(..., 'yscale', YSCALE)
   %    Choosing 'log' Allows the y-axis to be generated on a log-frequency
   %    scale, with uneven vertical cell spacing.  The default value is
   %    'normal', and provides the standard spectrogram view.
   %    valid yscales: 'normal', 'log' (see NOTE below)
   %
   %    NOTE: In order to use the log scale, UIMAGESC needs to be available on
   %    the matlab path.  This routine was created by Frederic Moisy, and may
   %    be downloaded from the maltabcentral fileexchange (File ID: 11368).
   %    If this routine is not found,then the original spectrogram will be
   %    created.
   %
   %  specgram2(..., 'fontsize', FONTSIZE)
   %    Specify the font size for a spectrogram.  The default font size is 8.
   %
   %  specgram2(..., 'innerLabels', SHOWINNERLABELS)
   %    Suppress the labling of the inside graphs by setting SHOWINNERLABELS
   %    to false.  If this is false, then the frequency label only shows on
   %    the leftmost spectrograms, and the X-unit label only shows on the
   %    bottommost spectrograms.
   %
   %  Example 1:
   %    % The following plots a waveform using an alternate mapping, an xunit of
   %    % of 'hours', and with the y-axis plotted using a log scale.
   %    specgram2(spectralobject, waveform,...
   %      'colormap', alternateMap,'xunit','h','yscale','log')
   %
   %
   %  Example 2:
   %    % create an arbitrary subplot, and then plot multiple spectra
   %    a = subplot(3,2,1);
   %    specgram2(spectralobject,waves,'axis',a); % waves is an NxM waveform
   %
   %   See also SPECTRALOBJECT/SPECGRAM

   if ~isa(T,'TraceData')
            try
               T = SeismicTrace(T);
               disp('successfully converted to a SeismicTrace');
            catch er
      error('Spectralobject:specgram2:invalidArgument','Should work on a trace (ex. TraceData, SeismicTrace), not a %s',class(T));
            end
   end

   %% search for relevent property pairs passed as parameters
   p = parseSpecgramInputs(s, varargin);

   %% figure out exactly WHERE to plot the spectrogram(s)
   %find out area(axis) in which the spectrograms will be plotted
   clabel= 'Relative Amplitude  (dB)';

   if p.Results.axis == 0,
      clf;
      pos = get(gca,'position');
   else
      pos = get(p.Results.axis,'position');
   end
   if ~isempty(p.Results.position)
      pos = p.Results.position;
   end

   %% If there are multiple waveforms...
   % subdivide the axis and loop through specgram2 with individual waveforms.
   if numel(T) > 1
      if p.Results.axis== 0,
         myaxis = gca;
      end
      %create the colorbar if desired
      TraceSpectra.createcolorbar(s,p.Results.colorbar, clabel, p.Results.fontsize);
      h = TraceSpectra.subdivide_axes(myaxis, size(T));
      % remainingproperties = TraceSpectra.property2varargin(proplist);
      remainingproperties = TraceSpectra.buildParameterList( p.Unmatched);
      for n=1:numel(h)
         keepYlabel =  ~p.Results.innerlabels || (n <= size(h,1));
         keepXlabel = ~p.Results.innerlabels || (mod(n,size(h,2))==0);
         specgram2(s,T(n),...
            'xunit',p.Results.xunit,...
            'axis',h(n),...
            'fontsize',p.Results.fontsize,...
            'useXlabel',keepXlabel,...
            'useYlabel',keepYlabel,...
            'colorbar','none',...
            remainingproperties{:});
      end
      return
   end

   %% Plot the spectrogram with a wiggle on top and colorbar below

   %plot the wiggle
   ax_wiggle = subplot('position',wigglePosition(pos));
   plot(T,'xunit',p.Results.xunit,'autoscale',true,'fontsize',p.Results.fontsize);

   % make the axis tight, and keep axis info for later use with spectra
   axis(ax_wiggle,'tight');
   xAxisLims = get(ax_wiggle,'xlim');
   ticnos = get(ax_wiggle,'xtick');

   %plot the spectra
   ax_spectra = subplot('position',spectraPosition(pos));
   additionalParams = TraceSpectra.buildParameterList( p.Unmatched);
   specgram(s,T,...
      'xunit',p.Results.xunit,...
      'fontsize',p.Results.fontsize,...
      'yscale',p.Results.yscale,...
      'colorbar','none',...
      'axis',ax_spectra,...
      'suppressXlabel',p.Results.useXlabel,...
      'suppressYlabel',p.Results.useYlabel,...
      additionalParams{:});

   %make the axis match exactly with the waveform above
   set(ax_spectra,'xtick',ticnos);
   if ~strcmpi(p.Results.yscale,'log')
      % axis scaling doesn't work quite right at a log scale
      xlim(ax_spectra, xAxisLims);
   end

   title(''); %clear the title

   %create the colorbar if desired
   TraceSpectra.createcolorbar(s,p.Results.colorbar, clabel, p.Results.fontsize);
end

function wigPos = wigglePosition(pos)
   % wigglePosition   wiggle occupies the top 15% of the axis
   % pos = [ left, bottom, width, height ]
   % wavepos = [ left, bottom + height * 0.85,  width , height * 0.15] ;
   wigPos = pos .* [1, 1, 1, 0.15] + [0,  pos(4) * 0.85, 0, 0] ;
end

function specPos = spectraPosition(pos)
   %spectraPosition   spectra occupies the bottom 85% of the axis
   specPos = pos .* [1, 1, 1, 0.85 ];
end

author: geoscience-community-codes

repo: GISMO


text:
function [h, compUpAP, compUp] =  lfmComputeH3AP(gamma1_p, gamma1_m, sigma2, t1, ...
    t2, preFactor, mode)

% LFMCOMPUTEH3AP Helper function for computing part of the LFMAP kernel.
% FORMAT
% DESC computes a portion of the LFMAP kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG preFactor : precomputed constants.
% ARG mode: indicates the correct derivative.
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

% Evaluation of h

if nargout>1
    [compUpAP{1}, compUp{1}] = lfmapComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode);
    [compUpAP{2}, compUp{2}] = lfmapComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
    h = preFactor(1)*compUpAP{1} + preFactor(2)*compUpAP{2};
else
    h = preFactor(1)*lfmapComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode) ...
        + preFactor(2)*lfmapComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function objNew = copy(obj)
% COPY Create an indipendent copy of a point cloud object.
% ------------------------------------------------------------------------------
% DESCRIPTION/NOTES
% The pointCloud class is a subclass of the handle class, i.e. it is NOT a value
% class. More informations about the difference between handle and value classes
% can be found in [1]. To create an indipendent copy of a point cloud object,
% this method can be used.
% ------------------------------------------------------------------------------
% OUTPUT
% 1 [objNew]
%   Copy of point cloud object.
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 Copy an object and select a subset of points in the new object.
%   pcOrig = pointCloud('Lion.xyz');
%   pcCopy = pcOrig.copy;
%   % Attention: if instead of the previous line 'pcCopy = pcOrig;' is used,
%   % pcOrig and pcCopy share the same object, i.e. they are NOT indipendent.
%   pcOrig.select('RandomSampling', 10); % select a subset of points in pcOrig
%   pcOrig.info; % in pcOrig only a subset of points is selected
%   pcCopy.info; % in pcCopy still all points are selected
% ------------------------------------------------------------------------------
% REFERENCES
% [1] % https://de.mathworks.com/help/matlab/matlab_oop/comparing-handle-and-value-classes.html
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Reference: https://de.mathworks.com/matlabcentral/newsreader/view_thread/257925

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'COPY'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));

% Copy -------------------------------------------------------------------------

% Deactivate logging
originalLogLevel = msg('O', 'GetLogLevel');
msg('O', 'SetLogLevel', 'off');

% Instantiate new object of the same class
objNew = feval(class(obj), obj.X);

% Copy all non-hidden properties
p = properties(obj);
for i = 1:length(p)
    if ~any(strcmp(p{i}, {'cog' 'lim' 'noPoints'}))
        objNew.(p{i}) = obj.(p{i});
    end
end

% Reactivate logging
msg('O', 'SetLogLevel', originalLogLevel);

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

end

author: pglira


text:
function value = function_3d_name ( ifunc )

%*****************************************************************************80
%
%% FUNCTION_3D_NAME returns the name of the current 3D function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IFUNC, the index of the function.
%
%    Output, character FNAME(7), the name of the function.
%
  if ( ifunc == 1 )
    value = '      1';

  elseif ( ifunc == 2 )
    value = '      X';
  elseif ( ifunc == 3 )
    value = '      Y';
  elseif ( ifunc == 4 )
    value = '      Z';

  elseif ( ifunc == 5 )
    value = '    X*X';
  elseif ( ifunc == 6 )
    value = '    X*Y';
  elseif ( ifunc == 7 )
    value = '    X*Z';
  elseif ( ifunc == 8 )
    value = '    Y*Y';
  elseif ( ifunc == 9 )
    value = '    Y*Z';
  elseif ( ifunc == 10 )
    value = '    Z*Z';

  elseif ( ifunc == 11 )
    value = '    X^3';
  elseif ( ifunc == 12 )
    value = '  X*Y*Z';
  elseif ( ifunc == 13 )
    value = '  Z*Z*Z';

  elseif ( ifunc == 14 )
    value = '    X^4';
  elseif ( ifunc == 15 )
    value = 'X^2 Z^2';
  elseif ( ifunc == 16 )
    value = '    Z^4';

  elseif ( ifunc == 17 )
    value = '    X^5';
  elseif ( ifunc == 18 )
    value = '    X^6';
  elseif ( ifunc == 19 )
    value = '      R';
  elseif ( ifunc == 20 )
    value = ' SIN(X)';
  elseif ( ifunc == 21 )
    value = ' EXP(X)';
  elseif ( ifunc == 22 )
    value = '1/(1+R)';
  elseif ( ifunc == 23 )
    value = 'SQRT(R)';
  else
    value = '???????';
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [C,T,A]=naca4_3d(airfoil,y,nc)
% determine airfoil skin and camber coordinates from NACA designation
%% NacaCoord.m
% NACA 4-digit or 5-digit series airfoil
%
% outputs:  C = matrix of camber line coordinates, x,y,z
%               ordered trailing edge first
%           T = matrix of skin coordinates
%               ordered from trailing edge over top surface then
%               along bottom surface to trailing edge
%           A = area of airfoil cross section (normalized)
%
% NOTE: Chord of output airfoil is 1
%
% inputs:   airfoil = NACA 4- or 5-digit airfoil as string
%           y = spanwise location
%           nc = number of chordwise panels
%
% Ref:      http://www.aerospaceweb.org/question/airfoils/q0041.shtml

%           Ira H. Abbott and Albert E. Von Doenhoff, Theory of Wing
%           Sections, 1959, Dover Publications, NY.
%
% Coordinate system:
%   x: chord direction, origin at LE, positive is pointing toward TE
%   y: station along span, perpendicular to aircraft midplane
%   z: perpendicular to chord, positive is up

%This 'if' statement establishes the constants for either 4- or 5- digit
%airfoils.
if length(airfoil) == 4
    m = str2num(airfoil(1))/100;
    p = str2num(airfoil(2))/10;
    t = str2num(airfoil(3:4))/100;

elseif length(airfoil) == 5
    switch airfoil(1:3)
        case '210'
            m = 0.0580;
            k1 = 361.4;
        case '220'
            m = 0.1260;
            k1 = 51.64;
        case '230'
            m = 0.2025;
            k1 = 15.957;
        case '240'
            m = 0.2900;
            k1 = 6.643;
        case '250'
            m = 0.3910;
            k1 = 3.230;
        otherwise
    end
    t = str2num(airfoil(4:5))/100;
end

%% Constants
% nc          % number of line segments describing top surface
% (= number of line segments describing bottom surface)
npc = nc+1;      % number of points along mean camber line
nps = 2*nc+1;    % number of points along surface
dx=1/nc;        % x-direction increment (along chord)

%% Fill coordinate matrices with y coordinate
C = zeros(npc,3); % mean camber
C(:,2) = y;
T = zeros(nps,3); %skin
T(:,2) = y;

%% Mean Camber Line and Slope
% dzcdx is the derivative of the mean camber line
% theta is the local angle of the mean camber line
% x=1:-dx:0; % order x locations starting at trailing edge, x=1
x=0:dx:1; % order x locations starting at leading edge, x=0
if length(airfoil) == 4
    for i = 1:nc+1
        if x(i)<p    % leading edge
            zc(i) = m/(p^2)*(2*p*x(i)-x(i)^2);  %Eqn 6.4 Abbot
            dzcdx(i) = m/p^2*(2*p-2*x(i));
            theta(i) = atan(dzcdx(i));
        else            % trailing edge
            zc(i) = m/(1-p)^2*((1-2*p)+2*p*x(i)-x(i)^2);    %Eqn 6.4 Abbot
            dzcdx(i) = m/(1-p)^2*(2*p-2*x(i));
            theta(i) = atan(dzcdx(i));
        end
    end
elseif length(airfoil) == 5
    for i = 1:nc+1
        if x(i)<m    % leading edge
            zc(i) = 1/6*k1*(x(i)^3-3*m*x(i)^2+m^2*(3-m)*x(i));  %Eqn 6.6 Abbot
            dzcdx(i) = 1/6*k1*(3*x(i)^2-6*m*x(i)+m^2*(3-m));
            theta(i) = atan(dzcdx(i));
        else            % trailing edge
            zc(i) = 1/6*k1*m^3*(1-x(i));    %Eqn 6.6 Abbot
            dzcdx(i) = -1/6*k1*m^3;
            theta(i) = atan(dzcdx(i));
        end
    end
end

C(:,1)=-x;  % output mean camber line coordinates; negative x due to axis convention
C(:,3)=-zc; %negative z due to axis convention

%% Thickness distribution is the same for both 4- and 5-digit airfoils
for i=1:nc+1
    zt(i) = t/0.2*[0.2969*x(i)^0.5-0.1260*x(i)-0.3516*x(i)^2+0.2843*x(i)^3-0.1015*x(i)^4];  %Eqn 6.2 Abbott
end

%% Area
A=0;
for i=1:nc
    da=dx*(zt(i)+zt(i+1));  %Take the average of zt(i) and zt(i+1) and double it
    A=A+da;
end
A;

%% Upper skin
for i=1:nc+1
    xu(i) = x(i)-zt(i)*sin(theta(i));  %Eqn 6.1 Abbott
    zu(i) = zc(i)+zt(i)*cos(theta(i));  %Eqn 6.1 Abbott
    T(i,1) = -xu(i);  %negative x due to sign convention
    T(i,3) = -zu(i);  %negative z due to sign convention
end

%% lower skin
% exclude point nc+1 at leading edge
% because it is included as the last point of the upper surface

for i=1:nc
    xl(i) = x(i)+zt(i)*sin(theta(i));    %Eqn 6.1 Abbott
    zl(i) = zc(i)-zt(i)*cos(theta(i));    %Eqn 6.1 Abbott
end

% reverse order so LE to TE
for i=nc+2:nps
    T(i,1) = -xl(2*(nc+1)-i);%negative x due to sign convention
    T(i,3) = -zl(2*(nc+1)-i);%negative z due to sign convention
end
%figure
%plot(xu,zu,'o')
%hold on
%plot(xl,zl,'or')
%axis equal

author: Sable

repo: mcbench-benchmarks


text:
function wiener_iter(ns_file, es_file, iter_num)

%
%  Implements the basic iterative Wiener filtering algorithm [1].
%
%  Usage:  wiener_iter(noisyFile, outputFile,NumberOfIterations)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%         NumberOfIterations - number of iterations (recommended 2-4)
%
%  Example call:  wiener_iter('sp04_babble_sn10.wav','out_wien.wav',3);
%
%  References:
%   [1] Lim, J. and Oppenheim, A. V. (1978). All-pole modeling of degraded speech.
%       IEEE Trans. Acoust. , Speech, Signal Proc., ASSP-26(3), 197-210.
%
% Authors: Yi Hu and Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<3
   fprintf('Usage: wiener_iter(noisyfile.wav,outFile.wav,NumIter) \n\n');
   return;
end

NF_SABSENT= 6;
%this is the number of speech-absent frames to estimate the initial
%noise power spectrum

[nsdata, Fs, bits]= wavread( ns_file);	%nsdata is a column vector

nwind= floor( 20* Fs/ 1000);	%this corresponds to 20ms window
if rem( nwind, 2)~= 0 nwind= nwind+ 1; end	%made window length even
noverlap= nwind/ 2;	w= hanning( nwind);	rowindex= ( 1: nwind)';
pred_order= 12;	%LPC order is set to 12
FFTlen=2*nwind;

%we assume the first NF_SABSENT frames are speech absent, we use them to estimate the noise power spectrum
noisedata= nsdata( 1: nwind* NF_SABSENT);	noise_colindex= 1+ ( 0: NF_SABSENT- 1)* nwind;
noisematrixdata = zeros( nwind, NF_SABSENT);
noisematrixdata( :)= noisedata( ...
   rowindex( :, ones(1, NF_SABSENT))+ noise_colindex( ones( nwind, 1), :)- 1);
noisematrixdata= noisematrixdata.* w( :, ones( 1, NF_SABSENT)) ;	%WINDOWING NOISE DATA
noise_r0= sum( sum( noisematrixdata.^ 2)/ nwind)/ NF_SABSENT;	%noise energy
noise_ps= mean( (abs( fft( noisematrixdata,FFTlen))).^ 2, 2);
%NOTE!!!! it is a column vector

nslide= nwind- noverlap;

x= nsdata( nwind* NF_SABSENT+ 1: end);	%x is to-be-enhanced noisy speech
nx= length( x);	ncol= fix(( nx- noverlap)/ nslide);
colindex = 1 + (0: (ncol- 1))* nslide;
if nx< (nwind + colindex(ncol) - 1)
   x(nx+ 1: nwind+ colindex(ncol) - 1) = ...
      rand( nwind+ colindex( ncol)- 1- nx, 1)* (2^ (-15));   % zero-padding
end

es_old= zeros( noverlap, 1);
%es_old is actually the second half of the previous enhanced speech frame,
%it is used for overlap-add

img=sqrt(-1);
for k= 1: ncol

   y= x( colindex( k): colindex( k)+ nwind- 1);
   y= y.* w;	%WINDOWING NOISY SPEECH DATA

   y_spec= fft(y,FFTlen);	y_specmag= abs( y_spec);	y_specang= angle( y_spec);
   %they are the frequency spectrum, spectrum magnitude and spectrum phase, respectively

   y_ps= y_specmag.^ 2;	%power spectrum of noisy speech

   %we must have a set of initial LPC coefficients to use for iterative wiener algorithm,
   %we get it from original noisy speech file

   lpc_coeffs= (lpc(y, pred_order))';
   %now we get initial lpc coefficients of all current speech frame

   for m= 1: pred_order+ 1
      %exp_matrix(m, :)= exp(-i* (m- 1)* ((1: nwind)- 1)* 2* pi/ nwind);
      exp_matrix(m, :)= exp(-img* (m- 1)* ((1: FFTlen)- 1)* 2* pi/FFTlen);
   end

   x_old_spec=y_spec;
   for n=1:iter_num

      xx= 1./ (abs( exp_matrix'* lpc_coeffs).^ 2);

      lpc_energy= mean( xx);

      tmp= y_ps- noise_ps;

      g= max( mean( tmp)./ lpc_energy, 1e-16);

      tmp1= g.* xx;
      h_spec= tmp1./ (tmp1+ noise_ps);  % Wiener filter

      es_tmpspec= x_old_spec.* h_spec;
      es_tmp= real( ifft( es_tmpspec,FFTlen));

      x_old_spec = fft(es_tmp, FFTlen);
      % ----
      if n~= iter_num
         lpc_coeffs= lpc( es_tmp, pred_order)';
      end
   end

   es_data( colindex( k): colindex( k)+ nwind- 1)= [es_tmp( 1: noverlap)+ es_old;...
         es_tmp( noverlap+ 1: nwind)];
   %overlap-add
   es_old= es_tmp( nwind- noverlap+ 1: nwind);
end

wavwrite( es_data, Fs, bits, es_file);

author: jtkim-kaist

repo: Speech-enhancement


text:
function output = TestCopySynthRegressionR(n_test)
output = false;
original_speech_dir = '/Users/kawahara/Music/VCTK_CORPUS/VCTK-Corpus/wav48/';
if isOctave
  target_analysis_dir = '/Users/kawahara/m-file/STRAIGHTV40_007e/analysisDataO/';
  target_wave_dir = '/Users/kawahara/m-file/STRAIGHTV40_007e/waveDataO/';
else
  target_analysis_dir = '/Users/kawahara/m-file/STRAIGHTV40_007e/analysisDataR/';
  target_wave_dir = '/Users/kawahara/m-file/STRAIGHTV40_007e/waveDataR/';
end;
target_files = dir([target_wave_dir '*.wav']);
n_files = length(target_files);
selected_id = randi(n_files, n_test);
command1 = 'rand("seed", 12345);';
command2 = 'randn("seed", 12345);';
for ii = 1:n_test
  tmp_name = target_files(selected_id(ii)).name;
  [x, fs] = audioread([original_speech_dir tmp_name(1:4) '/' tmp_name]);
  disp([num2str(ii) ': ' tmp_name ' ' datestr(now)]);
  if isOctave
    eval(command1);
    eval(command2);
  else
    rng(12345); % initialize frozen random number
  end;
  f0raw = MulticueF0v14(x, fs);
  ap = exstraightAPind(x, fs, f0raw);
  n3sgram=exstraightspec(x, f0raw, fs);
  tmp_name_root = tmp_name(1:end - 4);
  if ~CheckAnalysisData(f0raw, ap, n3sgram, target_analysis_dir, tmp_name_root)
    disp(['Failed: ' tmp_name ' data is not similar.']);
  end;
  wave_pathname = [target_wave_dir tmp_name];
  [sy, fs] = audioread(wave_pathname);
  if isOctave
    eval(command1);
    eval(command2);
  else
    rng(12345); % initialize frozen random number
  end;
  y = exstraightsynth(f0raw,n3sgram,ap,fs);
  disp(['Relative error SD: ' num2str(100 * std(sy - y / max(abs(y)) * 0.9) / std(sy)) ' %']);
  if std(sy - y / max(abs(y)) * 0.9) / std(sy) > 10 ^ (-3)
    disp(['Failed! ' tmp_name ' copy synthesis test.']);
    keyboard
    return;
  end;
end;
disp(['Success! ' num2str(n_test) ' files are passed copy-synth regression.']);
output = true;
end

author: HidekiKawahara


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% SABR density script using different parameters

% base scenario
clear; clc;
t=0; T=1;
f = 0.03;
alpha = 0.075;
beta = 0.5;
nu = 0.1;
rho = 0;

sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);

lcoeff = 1; ucoeff = 1;
lowerbound = f-lcoeff*f;
upperbound = f + ucoeff*f;

xvals = lowerbound:.0005:upperbound;

zvals1d = density1d(xvals,sabrdensity,0,1);

% change rho
rho = -.9;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
rho = .9;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \rho')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\rho=-0.9','\rho=0','\rho=0.9')
hold off;
rho = 0.0;

% change nu
nu = 0.05;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
nu = .15;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \nu')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\nu=0.05','\nu=0.4','\nu=0.2')
hold off;
nu = 0.1;

% change beta
beta = 0.4;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
beta = 0.8;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \beta')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\beta=0.4','\beta=0.5','\beta=0.6')
hold off;
beta = 0.5;

% change alpha
alpha = 0.05;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dm = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dm)
alpha = .1;
sabrdensity = @(x,y) psabr3(t,T,f,x,alpha,y,beta,nu,rho);
zvals1dp = density1d(xvals,sabrdensity,0,1);
intdensity(xvals,zvals1dp)
figure('Color',[1 1 1]);
hold on;
plot(xvals,zvals1dm,'r','LineWidth',2); plot(xvals,zvals1d,'b','LineWidth',2); plot(xvals,zvals1dp,'g','LineWidth',2);
title('Plot of SABR densities for different \alpha')
xlabel('0\leq x \leq .1')
ylabel('density')
legend('\alpha=0.05','\alpha=0.075','\alpha=0.1')
hold off;
alpha = 0.075;


text:
function [F_xw,F_sprocedure] = filter_sprocedure(F_xw,w,uncertaintyModel,ops)

% Select the constraints to analyse
if any(is(F_xw,'elementwise'))
    F_lp = F_xw(find(is(F_xw,'elementwise')));
    % and save the rest for later analysis
    F_xw = F_xw(find(~is(F_xw,'elementwise')));

    p = sdpvar(F_lp);
    keep = ones(1,length(p));
    F_sprocedure = [];

    % We cannot use a duality based SOS-decomposition
    opsin = ops;
    ops.verbose = 0;
    ops.sos.model=2;

    % All decision variables (i.e. not variables in SOS)
    Parameters = [recover(setdiff(depends(p),depends(w)))];
    for i = 1:length(p)
        d = degree(p(i),w);
        if all(d<=2) & any(d==2)
            if opsin.verbose
                disp(' - Using exact S-procedure to eliminate uncertainty');
            end
            lambda = sdpvar(1);
            if isempty(uncertaintyModel{1}.r)
                g = uncertaintyModel{1}.g;
            else
                e = (w-uncertaintyModel{1}.center);
                g = uncertaintyModel{1}.r^2-e'*e;
            end
            Parameters = [Parameters;lambda];
            F_sprocedure = [F_sprocedure, sos(p(i)-lambda*g), lambda>=0];
            %F_sprocedure = [F_sprocedure, lambda>0, compilesos(sos(p(i)-lambda*g),[],ops,s)];
            keep(i) = 0;
        elseif 0%any(d>2)
            [multiplier, lambda] = polynomial(w,2);
            e = (w-uncertaintyModel{1}.center);
            g = uncertaintyModel{1}.r^2-e'*e;
            Parameters = [Parameters;lambda];
            F_sprocedure = [F_sprocedure,lambda>=0,sos(p(i)-multiplier*g), sos(multiplier)];
            %F_sprocedure = [F_sprocedure, lambda>0, compilesos(sos(p(i)-multiplier*g)+sos(multiplier),[],ops,s)];
            keep(i) = 0;
        end
    end
    if ~isempty(F_sprocedure)
        F_sprocedure = compilesos(F_sprocedure,[],ops,Parameters);
    end
    if any(keep)
        F_xw = [F_xw,p(find(keep))>=0];
    end
else
    F_sprocedure = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function val=calcRMSE(xTrue,xEst,is3D)
%%CALCRMSE Compute the scalar root-mean-squared error (RMSE) of estimates
%          compared to true values.
%
%INPUTS: xTrue The truth data. This is either an xDimXNumSamples matrix or
%              an xDimXNXnumSamples matrix. The latter formulation is
%              useful when the MSE over multiple Monte Carlo runs of an
%              entire length-N track is desired. In the first formulation,
%              we take N=1. Alternatively, if the same true value is used
%              for all numSamples, then xTrue can just be an xDimXN matrix.
%              Alternatively, if xTrue is the same for all numSamples and
%              all N, then just an xDimX1 matrix can be passed. N and
%              numSamples are inferred from xEst.
%         xEst An xDimXnumSamples set of estimates or an xDimXNXnumSamples
%              set of estimates (if values at N times are desired).
%         is3D An optional indicating that xEst is 3D. This is only used if
%              xEst is a matrix. In such an instance, there is an ambiguity
%              whether xEst is truly 2D, so N=1, or whether numSamples=1
%              and xEst is 3D with the third dimension being 1. If this
%              parameter is omitted or an empty matrix is passed, it is
%              assumed that N=1.
%
%OUTPUTS: val A 1XN vector of scalar RMSE values.
%
%The RMSE is discussed with relation to alternative error measured in [1]
%and is given in Equation 1 of [1].
%
%EXAMPLE:
%Here, we show that the RMSE of samples of a Gaussian random variable is
%related to the square root of the trace of a covariance matrix.
% R=[28,   4, 10;
%     4,  22, 16;
%     10, 16, 16];%The covariance matrix.
% xTrue=[10;-20;30];
% numRuns=100000;
% xEst=GaussianD.rand(numRuns,xTrue,R);
% rootTrace=sqrt(trace(R))
% val=calcRMSE(xTrue,xEst)
%One will see that the root-trace and the sample RMSE are close.
%
%REFERENCES:
%[1] X. R. Li and Z. Zhao, "Measures of performance for evaluation of
%    estimators and filters," in Proceedings of SPIE: Conference on Signal
%    and Data processing of Small Targets, vol. 4473, San Diego, CA, 29
%    Jul. 2001, pp. 530-541.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(is3D))
    is3D=false;
end

xDim=size(xEst,1);
if(ismatrix(xEst)&&is3D==false)
    N=1;
    numSamples=size(xEst,2);
    xEst=reshape(xEst,[xDim,1,numSamples]);

    if(size(xTrue,2)==1)
        %If the true values are the same for all samples.
        xTrue=repmat(xTrue,[1,1,numSamples]);
    else
        xTrue=reshape(xTrue,[xDim,1,numSamples]);
    end
else
    N=size(xEst,2);
    numSamples=size(xEst,3);

    if(ismatrix(xTrue))
        if(size(xTrue,2)==1)
            %If the true values are the same for all samples and for all N.
            xTrue=repmat(xTrue,[1,N,numSamples]);
        else
            %If the true values are the same for all samples.
            xTrue=repmat(xTrue,[1,1,numSamples]);
        end
    end
end

val=zeros(1,N);
for k=1:N
    val(k)=sqrt(sum(sum((xEst(:,k,:)-xTrue(:,k,:)).^2,3)/numSamples));
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function Accuracy=P300_Ncorrect(smt,filepath,subject_info,spellerText_on)

file=fullfile(filepath,subject_info.subject,subject_info.session, 'log');
% get log file
list=ls(sprintf('%s\\*_p300_param*', file));
load_file=fullfile(file,list(end,:));
load(load_file);

score=0;
for i=1:length(spellerText_on)
    if strcmp(SAVE{i}.char(1),spellerText_on(i))
        score=score+1;
    end
end

Accuracy=score/length(spellerText_on);
disp(Accuracy);

for i=1:length(spellerText_on)
    a(i)=SAVE{i}.char(1);
    disp(SAVE{i}.char(1));
end
end

author: PatternRecognition

repo: OpenBMI


text:
function z = cpf_tangent(V, lam, Ybus, Sbusb, Sbust, pv, pq, ...
                            zprv, Vprv, lamprv, parameterization, direction)
%CPF_TANGENT  Computes normalized tangent predictor for continuation power flow
%   Z = CPF_TANGENT(V, LAM, YBUS, SBUSB, SBUST, PV, PQ, ...
%                                 ZPRV, VPRV, LAMPRV, PARAMETERIZATION, DIRECTION)
%
%   Computes a normalized tangent predictor for the continuation power flow.
%
%   Inputs:
%       V : complex bus voltage vector at current solution
%       LAM : scalar lambda value at current solution
%       YBUS : complex bus admittance matrix
%       SBUSB : handle of function returning nb x 1 vector of complex
%               base case injections in p.u. and derivatives w.r.t. |V|
%       SBUST : handle of function returning nb x 1 vector of complex
%               target case injections in p.u. and derivatives w.r.t. |V|
%       PV : vector of indices of PV buses
%       PQ : vector of indices of PQ buses
%       ZPRV : normalized tangent prediction vector from previous step
%       VPRV : complex bus voltage vector at previous solution
%       LAMPRV : scalar lambda value at previous solution
%       PARAMETERIZATION : value of cpf.parameterization option.
%       DIRECTION: continuation direction (+1 for postive lambda
%                  increase, -1 otherwise)
%
%   Outputs:
%       Z : the normalized tangent prediction vector

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

%% sizes
nb = length(V);
npv = length(pv);
npq = length(pq);

Vm = abs(V);

%% compute Jacobian for the power flow equations
[dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V);
[dummy, neg_dSdb_dVm] = Sbusb(Vm);
[dummy, neg_dSdt_dVm] = Sbust(Vm);
dSbus_dVm = dSbus_dVm - neg_dSdb_dVm - lam * (neg_dSdt_dVm - neg_dSdb_dVm);

j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
j12 = real(dSbus_dVm([pv; pq], pq));
j21 = imag(dSbus_dVa(pq, [pv; pq]));
j22 = imag(dSbus_dVm(pq, pq));

J = [   j11 j12;
        j21 j22;    ];

Sxf = Sbust(Vm) - Sbusb(Vm);    %% "transfer" at current voltage level
dF_dlam = -[real(Sxf([pv; pq])); imag(Sxf(pq))];
[dP_dV, dP_dlam] = cpf_p_jac(parameterization, zprv, V, lam, Vprv, lamprv, pv, pq);

%% linear operator for computing the tangent predictor
J = [   J   dF_dlam;
      dP_dV dP_dlam ];

% J = [ J, dF_dlam;
%       z([pv; pq; nb+pq; 2*nb+1])'];

%% compute normalized tangent predictor
z = zeros(size(zprv));
s = zeros(npv+2*npq+1, 1);
s(end,1) = sign(direction);
z([pv; pq; nb+pq; 2*nb+1]) = J\s;   %% tangent vector
z = z/norm(z);                      %% normalize tangent predictor

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function model = rmGridFit_twoGaussiansDoGbetafixed(model,prediction,data,params,t)
% rmGridFit_twoGaussiansDoG - core of DoG fit
%
% model = rmGridFit_twoGaussiansDoG(model,prediction,data,params);
%
% Second gaussian can go negative only.
%
% 2008/01 SOD: split of from rmGridFit.

% input check
if nargin < 4,
    error('Not enough arguments');
end

% some variables we need
rssinf         = inf(size(data(1,:)),'single');
trends         = t.trends;
t_id           = t.dcid+1;

% stimulus sequence is needed because make second pRF on the fly
allstimimages = rmDecimate(params.analysis.allstimimages,...
    params.analysis.coarseDecimate);

%-----------------------------------
%--- fit different receptive fields profiles
%--- another loop --- and a slow one too
%-----------------------------------
tic; progress = 0;
for n=1:numel(params.analysis.x0),
    %-----------------------------------
    % progress monitor (10 dots) and time indicator
    %-----------------------------------
    if floor(n./numel(params.analysis.x0).*10)>progress,
        if progress==0,
            % print out estimated time left
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d hours.\t(%s)\n',...
                    mfilename, ceil(esttime./3600), datestr(now));
            else
                fprintf(1, '[%s]:Estimated processing time: %d minutes.\t(%s)\n',...
                    mfilename, ceil(esttime./60), datestr(now));
            end;
            fprintf(1,'[%s]:Grid (x,y,sigma) fit:',mfilename);drawnow;
        end;
        % progress monitor
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    %-----------------------------------
    %--- first make all second rf profiles
    %-----------------------------------
    sigmaNew = params.analysis.sigmaRatio.*params.analysis.sigmaMajor(n);
    % limit to sigmaRatioMaxVal
    sigmaNew = sigmaNew(sigmaNew<=params.analysis.sigmaRatioMaxVal);
    % add sigmaRatioInfVal which essentially is on/off
    sigmaNew = [sigmaNew; params.analysis.sigmaRatioInfVal]; %#ok<AGROW>
    %betaratio = params.analysis.betaRatioAlpha.*(params.analysis.sigmaMajor(n).^2./sigmaNew.^2);
    betaratio = (params.analysis.sigmaMajor(n)./sigmaNew).^params.analysis.betaRatioAlpha;

    % Now we make it: slightly different call for speed reasons.
    z = zeros(size(sigmaNew));
    tmprf   = rfGaussian2d(params.analysis.X - params.analysis.x0(n),...
        params.analysis.Y - params.analysis.y0(n),...
        sigmaNew,sigmaNew, ...
        false, z, z);
    prediction2 = allstimimages*tmprf;

    %-----------------------------------
    %--- try two rf profiles, yet another loop
    %-----------------------------------
    for sr = 1:numel(sigmaNew),
        %-----------------------------------
        % Now apply GLM, see *_oneGaussian for logic.
        % New rules for this one:
        % 1. first Gaussian has to be positive
        % 2. there should be a positive response in the
        % center at all times. This
        % assumes Gaussians are unscaled (see rfGaussian2d).
        %-----------------------------------
        predictionNew = prediction(:,n) - betaratio(sr,:).*prediction2(:,sr);
        X     = [predictionNew trends];
        b     = pinv(X)*data;
        rss   = rssinf;
        keep  = b(1,:)>0;
        rss(keep) = sum((data(:,keep)-X*b(:,keep)).^2);

        %-----------------------------------
        %--- store data with lower rss
        %-----------------------------------
        minRssIndex = rss < model.rss;

        % now update
        model.x0(minRssIndex)       = params.analysis.x0(n);
        model.y0(minRssIndex)       = params.analysis.y0(n);
        model.s(minRssIndex)        = params.analysis.sigmaMajor(n);
        model.s_major(minRssIndex)        = params.analysis.sigmaMajor(n);
        model.s_minor(minRssIndex)        = params.analysis.sigmaMajor(n);
        model.s_theta(minRssIndex)        = params.analysis.theta(n);
        model.rss(minRssIndex)      = rss(minRssIndex);
        model.b([1 t_id],minRssIndex)    = b(:,minRssIndex);
        model.s2(minRssIndex)       = sigmaNew(sr);
    end;
end;

% end time monitor
et  = toc;
if floor(esttime/3600)>0,
    fprintf(1,'Done[%d hours].\t(%s)\n', ceil(et/3600), datestr(now));
else
    fprintf(1,'Done[%d minutes].\t(%s)\n', ceil(et/60), datestr(now));
end;
drawnow;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [faces,vertices] = mesh_emse2mat3d(vertex,patch)

% mesh_emse2mat3d - Convert emse vertex/patches to matlab vertices/faces
%
% Useage: [faces,vertices] = mesh_emse2mat3d(vertex,patch)
%
%           vertex & patch are generated from emse files with
%           mesh_emse2matlab (see this file for more help).
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  12/98 Abbas Kouzani
%           09/01 Darren.Weber_at_radiology.ucsf.edu
%                 - converted to function, rather than script
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%S=sprintf('load %s',file);eval(S);

vertices = zeros(length(vertex),3);
for i=1:length(vertex),
    vertices(i,1) = vertex(i).location1;
    vertices(i,2) = vertex(i).location2;
    vertices(i,3) = vertex(i).location3;
end

faces = zeros(length(patch),3);
for i=1:length(patch),

    v = patch(i).vertex1;   flag = 0;
    for j=1:length(vertex),
        if isequal(v, vertex(j).address)
            faces(i,1) = j; flag = 1; break;
        end
    end
    if (flag==0) fprintf('...error in patch(%d).vertex1!\n',i); break;end

    v = patch(i).vertex2;   flag = 0;
    for j=1:length(vertex),
        if isequal(v, vertex(j).address)
            faces(i,2) = j; flag = 1; break;
        end
    end
    if (flag==0) fprintf('...error in patch(%d).vertex2!\n',i); break;end

    v = patch(i).vertex3;   flag = 0;
    for j=1:length(vertex),
        if isequal(v, vertex(j).address)
            faces(i,3) = j; flag = 1; break;
        end
    end
    if (flag==0) fprintf('...error in patch(%d).vertex3!\n',i); break;end
end

%disp('...Saving');
%GET1=input(' Enter output filename: ','s');
%S=sprintf('save %s faces vertices',GET1);
%eval(S);


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function varargout = About_MU(varargin)
% ABOUT_MU M-file for About_MU.fig
%      ABOUT_MU, by itself, creates a new ABOUT_MU or raises the existing
%      singleton*.
%
%      H = ABOUT_MU returns the handle to a new ABOUT_MU or the handle to
%      the existing singleton*.
%
%      ABOUT_MU('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ABOUT_MU.M with the given input arguments.
%
%      ABOUT_MU('Property','Value',...) creates a new ABOUT_MU or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before About_MU_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to About_MU_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help About_MU

% Last Modified by GUIDE v2.5 07-Nov-2011 13:21:51

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @About_MU_OpeningFcn, ...
                   'gui_OutputFcn',  @About_MU_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before About_MU is made visible.
function About_MU_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to About_MU (see VARARGIN)

% Choose default command line output for About_MU
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes About_MU wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = About_MU_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function Y = Wavedb1Phi(X,trans)
persistent s;

level = 2;

if ~trans;
    % Phi' * X
    [Y,s] = wavedec2(X,level,'db1');
else
    Y = waverec2(X(:),s,'db1');
end

repo: REST


text:
function out=mtimes(x,y)

precisionx=0;   precisiony=0;
xmp=isa(x,'mp');
if xmp
 precisionx=x(1).precision;
 ymp=isa(y,'mp');
else
 ymp=true;
end
if ymp
 precisiony=y(1).precision;
end
precision=max(precisionx,precisiony);
ex=numel(x);
ey=numel(y);

if ex==1 | ey==1
 if ~xmp,  x=mp(x);  end
 if ~ymp,  y=mp(y);  end
 if ex==1
  out_rval=cell(size(y));
  out_ival=out_rval;

  [xrval,xival]=getVals(x,1);
  %xrval=x(1).rval;  xival=x(1).ival;  if isempty(xival), xival='0'; end

  xReal=hasimag(xival);
  for ii=1:ey
   [yrval,yival]=getVals(y,ii);
   %yrval=y(ii).rval;  yival=y(ii).ival;  if isempty(yival), yival='0'; end
   if xReal | hasimag(yival)
    [out_rval{ii},out_ival{ii}]=mpfr_mulc(precision,xrval,xival,yrval,yival);
   else
    %if isempty(xrval),  'xxxxxxxxxxxx',kb, end
    out_rval{ii}=mpfr_mul(precision,xrval,yrval);
   end
  end
 else
  out_rval=cell(size(x));
  out_ival=out_rval;
  [yrval,yival]=getVals(y,1);
  yReal=hasimag(yival);
  for ii=1:ex
   [xrval,xival]=getVals(x,ii);
   if hasimag(xival) | yReal
    [out_rval{ii},out_ival{ii}]=mpfr_mulc(precision,xrval,xival,yrval,yival);
   else
    out_rval{ii}=mpfr_mul(precision,xrval,yrval);
   end
  end
 end
 out=class(struct('rval',out_rval,...
                  'ival',out_ival,...
                  'precision',precision),'mp');
else %%% out=x*y;
 sx=size(x);
 sy=size(y);
 if sx(2)~=sy(1)
  disp(['size of matrix 1 => ',num2str(sx(1)),'x',num2str(sx(2))]);
  disp(['size of matrix 2 => ',num2str(sy(1)),'x',num2str(sy(2))]);
  error('array size mismatch for mtimes')
 else
  out=mp(zeros(sx(1),sy(2)),precision);
  for ii=1:sx(1)
   for jj=1:sy(2)
    out(ii,jj)=sum(x(ii,1:end).*y(1:end,jj).');
   end
  end
 end
end

author: opencobra

repo: cobratoolbox


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

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

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %
	%h=X*theta;
	%delta=1/m*(sum((h-y)'*X));  %delta is only 1x1 dimension
	%delta=1/m*(sum(h-y)*sum(X)); %array correct but turns to NaN
	%delta=1/m*(h-y).*X;
	delta=1/m*(X'*X*theta-X'*y);
	theta=theta-alpha.*delta;
	%fprintf(theta);

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

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

end
J_history
end

author: yhyap

repo: machine-learning-coursera


text:
function ndx = array2ind( siz, in_array, ~ )
%ARRAY2IND Linear index from multiple subscripts.
%   ARRAY2IND is used to determine the equivalent single index
%   corresponding to a given set of subscript values.
%
%   IND = ARRAY2IND(SIZ,[I,J]) returns the linear index equivalent to the
%   row and column subscripts in the arrays I and J for a matrix of
%   size SIZ.
%
%   IND = ARRAY2IND(SIZ,[I1,I2,...,IN]) returns the linear index
%   equivalent to the N subscripts in the arrays I1,I2,...,IN for an
%   array of size SIZ.
%
%   I1,I2,...,IN must have the same size, and IND will have the same size
%   as I1,I2,...,IN. For an array A, if IND = ARRAY2IND(SIZE(A),I1,...,IN)),
%   then A(IND(k))=A(I1(k),...,IN(k)) for all k.
siz = double(siz);
if length(siz)<2
    error('Invalid Input Size');
end

if length(siz)~=length(in_array)
    error('Invalid Input Size');
end

%Compute linear indices
k = [1 cumprod(siz(1:end-1))];
ndx = 1;
s = size(in_array(1)); %For size comparison

for i = 1:length(siz),
    v = in_array(i);
    %%Input checking
    if ~isequal(s,size(v))
        %Verify sizes of subscripts
        error('SubscriptVectorSize');
    end
    if (any(v(:) < 1)) || (any(v(:) > siz(i)))
        %Verify subscripts are within range
        error('IndexOutOfRange');
    end
    ndx = ndx + (v-1)*k(i);
end

end


text:
function value = r8mat_amax ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_AMAX returns the maximum absolute value entry of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the M by N matrix.
%
%    Output, real VALUE, the maximum absolute value entry of A.
%
  value = max ( max ( abs ( a(1:m,1:n) ) ) );

  return
end


text:
% PRMF (Wang et al. 2012)
% process_video('RPCA', 'PRMF', 'dataset/demo.avi', 'output/demo_PRMF.avi');
X = normalize(M);
rk = 2;
lambdaU = 1;
lambdaV = 1;
tol = 1e-2;
[P, Q] = RPMF(X, rk, lambdaU, lambdaV, tol);
L = P * Q;
%S = abs(X - P * Q);
S = X - L;

author: andrewssobral

repo: lrslibrary


text:
function SV=gmm2sv(mix)
%--------------------------------------------------------------------------
%Usage: SuperVECTOR=gmm2sv(mix);
%       input : mix: GMM structure of Netlab
%       output: SuperVECTOR: GMM Supervector
%
% This program concatenates means of GMM to form SuperVector.
% Written by Md. Sahidullah (Graduate Student, IIT Kharagpur)
% Use this code works with Netlab toolbox
% http://www1.aston.ac.uk/eas/research/groups/ncrg/resources/netlab/
% If you have any query or suggestion please mail me sahidullahmd@gmail.com
temp=mix.centres';
SV=temp(:);
%--------------END OF CODE-------------------------------------------------


text:
function out = read_tiny_big_metadata(words_in,ind,fields)

%%% Routine to read in metadata for tiny images from giant binary file
%holding all metdata.

% valid field names:
% 1. keyword
% 2. filename
% 3. width
% 4. height
% 5. colors
% 6. date
% 7. engine
% 8. thumb_url
% 9. source_url
% 10. page
% 11. ind_page
% 12. ind_engine
% 13. ind_overall
% 14. label (1 = correct, 0 = incorrect, -1 = unlabelled)

%%% There are two modes of operation:
%
% 1. raw indexing mode
% --------------------
% Inputs:
% 1. words_in - 1 x nImages vector, each element being an integer
% in the range 1 to 73,777,893 (or whatever the total number of images in the dataset is).
% 2.  ind - 1 x nFields cell array holding requested fields (see
% list above) [optional - if ommitted, all fields will be returned].
% Outputs:
% 1. out - 1 x nImages structure holding the fields requested

% e.g. out = read_tiny_big_metadata([1:10],{'filename' 'label'});

% ++++++++++++++++++++++++++++++++++++++++++++++++++++
%
% 2. word based indexing.
% -----------------------
% Inputs:
% 1. word_in - 1 x nWords cell array, each entry holding one noun
% 2. ind - 1 x nWords cell array, each holding indices of images
% for the corresponding noun
% 3.  fields - 1 x nFields cell array holding requested fields (see
% list above)
% Outputs:
% 1. out - 1 x nImages structure holding the fields requested

% e.g. out = read_tiny_big_metadata('aardvark',[1:10],{'filename' 'label'});

% Requires: read_tiny_metadata_big_core.mex file

%%%% Hints regarding operation:
% The access time per image drops as the total number of images to be
% gathered increases. e.g. for 500,000 images load time is
% 1ms/image. For 100,000 images it is 2ms/image. So try and pass in big
% blocks of indices.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET PATHS BELOW TO tiny_images.bin AND tiny_index.mat

% path to giant binary file
binary_fname = '/mit/tiny/data/tiny_metadata.bin';
% indexing info .mat file
binary_data_fname = '/mit/tiny/data/tiny_index.mat';

% END
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TOTAL_NUM_IMGS = 79302017;

field_names = {'keyword' 'filename' 'width' 'height' 'colors' 'date' ...
    'engine' 'thumb_url' 'source_url' 'page' 'ind_page' ...
    'ind_engine' 'ind_overall' 'label'};

% if field actually a number
numeric_flag = [0 0 1 1 1 0 0 0 0 1 1 1 1 1];

if (nargin==1)
    RAW_INDEX_MODE = 1;
    ind = field_names;
elseif (nargin==2)
    RAW_INDEX_MODE = 1;
    if ~iscell(ind)
        ind = {ind};
    end
else % full 3 inputs
    RAW_INDEX_MODE = 0;

    if ~iscell(ind)
        % convert to cell arrays
        words_in = {words_in};
        ind = {ind};
    end

    if ~iscell(fields)
        fields = {fields};
    end

    nWords = length(words_in);

end

if RAW_INDEX_MODE

    %% just go and call binary with raw indices

    %% convert field names into numbers
    field_ind = [];
    for a=1:length(ind)
        field_ind = [ field_ind , strmatch(ind{a},field_names) ];
    end

    if isempty(field_ind)
        error('No valid field names');
    else
        nFields = length(field_ind);
    end

    %% check that all indices are >0 and <TOTAL_NUM_IMGS
    if ((sum(words_in<1)==0) & (sum(words_in>TOTAL_NUM_IMGS)==0))

        % now sort for speed
        [sorted_words,sort_ind] = sort(words_in);

        %% use MEX for sppedy reading in
        out = read_tiny_metadata_big_core(binary_fname, uint64(sorted_words), uint32(field_ind));
    else
        error('indices <0 or greater than total number of images');
    end

else % words_in and indices passed in

    %% convert field names into numbers
    field_ind = [];
    for a=1:length(fields)
        field_ind = [ field_ind , strmatch(fields{a},field_names) ];
    end

    if isempty(field_ind)
        error('No valid field names');
    else
        nFields = length(field_ind);
    end

    % load up index data
    load(binary_data_fname);

    % compute hashes for all words_in
    hashes = compute_hash_function(words_in,hash_key);

    img_offsets = [];

    for a=1:nWords

        % find word using hash
        match_ind = find(hashes(a) == hash);

        if any(match_ind)
            % check all ind are >=1 and <=total number of images

            if ((sum(ind{a}<1)==0) & (sum(ind{a}>num_imgs(match_ind))==0))
                word_offset = offset(match_ind);
                img_offsets = [img_offsets , uint64(double(word_offset) + [ind{a}])];
            else
                error('indices <0 or greater than number per image');
            end
        else
            error('word not found');
        end

        % now sort for speed
        [sorted_offsets,sort_ind] = sort(img_offsets);

        out = read_tiny_metadata_big_core( binary_fname, uint64(sorted_offsets), uint32(field_ind));

    end
end

% convert to string and trim whitespace
for a=1:nFields
    if numeric_flag(field_ind(a))
        tmp_out = int16(str2num(char(out{a}')));

        %%% hack workaround due to 'co' and 'in' problem
        if (field_ind(a)==14) & isempty(tmp_out)
            out{a}(find(out{a}(:)==99)) = 49;
            out{a}(find(out{a}(:)==105)) = 48;
            out{a}(find(out{a}(:)==110)) = 32;
            out{a}(find(out{a}(:)==111)) = 32;
            tmp_out = int16(fast_str2num(char(out{a}')));
        end

        out{a} = tmp_out;

    else
        out{a} = cellstr(char(out{a}'));
    end
    % sort back
    out{a}(sort_ind,:) = out{a};
end

% put into cell array
out = cell2struct(out,field_names(field_ind),2);

author: huashiyiqike

repo: LSTM-MATLAB


text:
function print_names = WindowFeatureParams2String(pff,wfs)

print_names = {};

if isempty(wfs),
  return;
end
iscellinput = iscell(wfs{1});
if ~iscellinput,
  wfs = {wfs};
end
n = numel(wfs);

print_names = cell(1,n);
for i = 1:n,
  print_names{i} = pff;
  for j = 1:numel(wfs{i}),
    if ischar(wfs{i}{j}),
      print_names{i} = [print_names{i},'_',wfs{i}{j}];
    else
      print_names{i} = [print_names{i},'_',num2str(wfs{i}{j})];
    end
  end
end

if ~iscellinput,
  print_names = print_names{1};
end

author: kristinbranson

repo: JAABA


text:
function [nlogL,nlogLvar,exitflag,output] = ibslike(fun,params,respMat,designMat,options,varargin)
%IBSLIKE Unbiased negative log-likelihood via inverse binomial sampling.
%   NLOGL = IBLLIKE(FUN,PARAMS,RESPMAT,DESIGNMAT) returns unbiased estimate
%   NLOGL of the negative of the log-likelihood for the simulated model
%   and data, calculated using inverse binomial sampling (IBS).
%   FUN is a function handle to a function that simulates the model's
%   responses (see below).
%   PARAMS is the parameter vector used to simulate the model's responses.
%   RESPMAT is a "response" data matrix, where each row correspond to one
%   observation or "trial" (e.g., a trial of a psychophysical experiment),
%   and each column represents a different response feature (e.g., the
%   subject's response and reported confidence level). Responses need to
%   belong to a finite set.
%   DESIGNMAT is an optional experimental design matrix, where each row
%   corresponds to one trial, and each column corresponds to a different
%   trial feature (such as condition, stimulus value, etc.).
%
%   FUN takes as input a vector of parameters PARAMS and an experimental
%   design matrix DMAT (one row per trial), and generates a matrix of
%   simulated model responses (one row per trial, corresponding to rows of
%   DMAT). DMAT is built by the algorithm out of rows of DESIGNMAT.
%
%   DESIGNMAT can be omitted or left empty, in which case FUN needs to
%   accept a parameter vector PARAMS and an array of trial numbers T, and
%   returns a matrix of simulated responses, where the i-th row contains
%   the simulated response for trial T(i) (the indices in T may repeat).
%
%   NLOGL = IBSLIKE(FUN,PARAMS,RESPMAT,DESIGNMAT,OPTIONS) uses options in
%   structure OPTIONS to replace default values. (To be explained...)
%
%   NLOGL = IBSLIKE(...,VARARGIN) additional arguments are passed to FUN.
%
%   [NLOGL,NLOGLVAR] = IBSLIKE(...) also returns an estimate NLOGVAR of the
%   variance of the log likelihood.
%
%   [NLOGL,NLOGLVAR,EXITFLAG] = IBSLIKE(...) returns an EXITFLAG that
%   describes the exit condition. Possible values of EXITFLAG and the
%   corresponding exit conditions are
%
%    2  IBS terminated after reaching the maximum runtime specified by the
%       user (the estimate can be arbitrarily biased).
%    1  IBS terminated after reaching the negative log-likelihood
%       threshold specified by the user (the estimate is biased).
%    0  Correct run of IBS; the estimate is unbiased.
%
%   [NLOGL,NLOGLVAR,EXITFLAG,OUTPUT] = IBSLIKE(...) returns a structure
%   OUTPUT with additional information about the sampling.
%
%   OPTIONS = IBSLIKE('defaults') returns a basic default OPTIONS structure.
%
%   EXITFLAG = IBSLIKE('test') runs some tests. Here EXITFLAG is 0 if
%   everything works correctly.
%
%   Test code on binomial sampling:
%      p = 0.7; Ntrials = 100;                  % Define binomial probability
%      fun = @(x,dmat) rand(size(dmat)) < x;    % Simulating function
%      rmat = fun(p,NaN(Ntrials,1));            % Generate responses
%      [nlogL,nlogLvar,pc,output] = ibslike(fun,p,rmat);
%      nlogL_true = -log(p)*sum(rmat == 1) - log(1-p)*sum(rmat == 0);
%      fprintf('Ground truth: %.4g, Estimate: %.4g ± %.4g.\n',nlogL_true,nlogL,sqrt(nlogLvar));
%
%   Reference:
%   van Opheusden*, B., Acerbi*, L. & Ma, W. J. (2020), "Unbiased and
%   efficient log-likelihood estimation with inverse binomial sampling".
%   (* equal contribution), PLoS Computational Biology 16(12): e1008483.
%   Link: https://doi.org/10.1371/journal.pcbi.1008483
%
%   See also @.

%--------------------------------------------------------------------------
% IBS: Inverse Binomial Sampling for unbiased log-likelihood estimation
% To be used under the terms of the MIT License
% (https://opensource.org/licenses/MIT).
%
%   Authors (copyright): Luigi Acerbi and Bas van Opheusden, 2020-2022
%   e-mail: luigi.acerbi@helsinki.fi, svo@princeton.edu
%   URL: http://luigiacerbi.com
%   Version: 0.96
%   Release date: Jan 21, 2021
%   Code repository: https://github.com/acerbilab/ibs
%--------------------------------------------------------------------------

if nargin < 4; designMat = []; end
if nargin < 5; options = []; end

t0 = tic;

% Default options
% defopts.Display       = 'off';        % Level of display on screen
defopts.Nreps           = 10;           % # independent log-likelihood estimates per trial
defopts.NegLogLikeThreshold = Inf;      % Stop sampling if estimated nLL is above threshold (incompatible with vectorized sampling)
defopts.Vectorized      = 'auto';       % Use vectorized sampling algorithm with acceleration
defopts.Acceleration    = 1.5;          % Acceleration factor for vectorized sampling
defopts.NsamplesPerCall = 0;            % # starting samples per trial per function call (0 = choose automatically)
defopts.MaxIter         = 1e5;          % Maximum number of iterations (per trial and estimate)
defopts.ReturnPositive  = false;        % If true, the first returned output is the *positive* log-likelihood
defopts.ReturnStd       = false;        % If true, the second returned output is the standard deviation of the estimate
defopts.MaxTime         = Inf;          % Maximum time for a IBS call (in seconds)
defopts.TrialWeights    = [];           % Vector of per-trial weights to yield a weighted sum of log-likelihood

%% If called with no arguments or with 'defaults', return default options
if nargout <= 1 && (nargin == 0 || (nargin == 1 && ischar(fun) && strcmpi(fun,'defaults')))
    if nargin < 1
        fprintf('Basic default options returned (type "help ibslike" for help).\n');
    end
    nlogL = defopts;
    return;
end

%% If called with the first argument as 'test', run test
if ischar(fun) && strcmpi(fun,'test')
    if nargin < 2; options = []; else; options = params; end
    figure;
    subplot(1,3,1);
    exitflag(1) = runtest1(options);
    subplot(1,3,2);
    exitflag(2) = runtest2(options);
    subplot(1,3,3);
    exitflag(3) = runtest3(options);
    nlogL = any(exitflag);
    return;
end

for f = fields(defopts)'
    if ~isfield(options,f{:}) || isempty(options.(f{:}))
        options.(f{:}) = defopts.(f{:});
    end
end

Ntrials = size(respMat,1);

% Add hard-coded options
options.MaxSamples = 1e4;                       % Maximum # of samples per function call
options.AccelerationThreshold = 0.1;            % Keep accelerating until threshold is passed (in s)
options.VectorizedThreshold  = 0.1;             % Max threshold for using vectorized algorithm (in s)
options.MaxMem = 1e6;                           % Maximum number of samples for vectorized implementation
options.MaxMem = max(min(Ntrials,1e4),10)*100;  % Maximum number of samples for vectorized implementation

% NSAMPLESPERCALL should be a scalar integer
if ~isnumeric(options.NsamplesPerCall) || ~isscalar(options.NsamplesPerCall)
    error('ibslike:NsamplesPerCall','OPTIONS.NsamplesPerCall should be a scalar integer.');
end

% ACCELERATION should be a scalar equal or greater than 1
if ~isnumeric(options.Acceleration) || ~isscalar(options.Acceleration) || ...
        options.Acceleration < 1
    error('ibslike:Acceleration','OPTIONS.Acceleration should be a scalar equal or greater than one.');
end

% NEGLOGLIKETHRESHOLD should be a scalar greater than 0 (or Inf)
if ~isnumeric(options.NegLogLikeThreshold) || ~isscalar(options.NegLogLikeThreshold) || ...
        options.NegLogLikeThreshold <= 0
    error('ibslike:NegLogLikeThreshold','OPTIONS.NegLogLikeThreshold should be a positive scalar (including Inf).');
end

% MAXTIME should be a positive scalar (including Inf)
if ~isnumeric(options.MaxTime) || ~isscalar(options.MaxTime) || ...
        options.MaxTime <= 0
    error('ibslike:MaxTime','OPTIONS.MaxTime should be a positive scalar (or Inf).');
end

% WEIGHT vector should be a scalar or same as number of trials
weights = options.TrialWeights(:);
if isempty(weights); weights = 1; end

if numel(weights) ~= 1 && numel(weights) ~= Ntrials
    error('ibslike:NumWeights','OPTIONS.TrialWeights should be empty, a scalar or an array of weights with as many elements as the number of trials.');
end

Trials = (1:Ntrials)';
funcCount = 0;

simdata = []; elapsed_time = [];

% Use vectorized or loop version?
if ischar(options.Vectorized) && options.Vectorized(1) == 'a'
    if options.Nreps == 1
        vectorized_flag = false;
    else
        % First full simulation to determine computation time
        fun_clock = tic;
        if isempty(designMat)    % Pass only trial indices
            simdata = fun(params,Trials(:),varargin{:});
        else                    % Pass full design matrix per trial
            simdata = fun(params,designMat(Trials(:),:),varargin{:});
        end
        elapsed_time = toc(fun_clock);
        vectorized_flag = elapsed_time < options.VectorizedThreshold;
        funcCount = 1;
    end
else
    vectorized_flag = logical(options.Vectorized);
    if options.Nreps == 1 && vectorized_flag
        vectorized_flag = false;
        warning('Vectorized IBS requires OPTIONS.Nreps > 1. Switching to non-vectorized algorithm.');
    end
end

if vectorized_flag
    [nlogL,K,Nreps,Ns,fc,exitflag] = ...
        vectorized_ibs_sampling(fun,params,respMat,designMat,simdata,elapsed_time,t0,options,varargin{:});
else
    [nlogL,K,Nreps,Ns,fc,exitflag] = ...
        loop_ibs_sampling(fun,params,respMat,designMat,simdata,elapsed_time,t0,options,varargin{:});
end
funcCount = funcCount + fc;

% Variance of estimate per trial
if nargout > 1
    K_max = max(max(K(:),1));
    Ktab = -(psi(1,1:K_max)' - psi(1,1));
    LLvar = Ktab(max(K,1));
    nlogLvar = sum(LLvar,2)./Nreps.^2;
end

% OUTPUT structure with additional information
if nargout > 2
    output.funcCount = funcCount;
    output.NsamplesPerTrial = Ns/Ntrials;
    output.nlogL_trials = nlogL;
    output.nlogLvar_trials = nlogLvar;
end

% Return negative log-likelihood and variance summed over trials
nlogL = sum(nlogL.*weights);
if options.ReturnPositive; nlogL = -nlogL; end
if nargout > 1
    nlogLvar = sum(nlogLvar.*(weights.^2));
    if options.ReturnStd    % Return standard deviation instead of variance
        nlogLvar = sqrt(nlogLvar);
    end
end

end

%--------------------------------------------------------------------------
function [nlogL,K,Nreps,Ns,fc,exitflag] = vectorized_ibs_sampling(fun,params,respMat,designMat,simdata0,elapsed_time0,t0,options,varargin)

Ntrials = size(respMat,1);
Trials = (1:Ntrials)';
Ns = 0;
fc = 0;
exitflag = 0;

Psi_tab = [];   % Empty PSI table

% Empty matrix of K values (samples-to-hit) for each repeat for each trial
K_mat = zeros([max(options.Nreps),Ntrials]);

% Matrix of rep counts
K_place0 = repmat((1:size(K_mat,1))',[1,Ntrials]);

% Current rep being sampled for each trial
Ridx = ones(1,Ntrials);

% Current vector of "open" K values per trial (not reached a "hit" yet)
K_open = zeros(1,Ntrials);

targetHits = options.Nreps(:)'.*ones(1,Ntrials);
MaxIter = options.MaxIter*max(options.Nreps);

% Starting samples
if options.NsamplesPerCall == 0
    samples_level = options.Nreps;
else
    samples_level = options.NsamplesPerCall;
end

for iter = 1:MaxIter
    % Pick trials that need more hits, sample multiple times
    T = Trials(Ridx <= targetHits);
    if isfinite(options.MaxTime) && toc(t0) > options.MaxTime
        T = [];
        exitflag = 2;
    end
    if isempty(T); break; end

    Ttrials = numel(T);    % Number of trials under consideration

    % With accelerated sampling, might request multiple samples at once
    Nsamples = min(options.MaxSamples,max(1,round(samples_level)));
    MaxSamples = ceil(options.MaxMem / Ttrials);
    Nsamples = min(Nsamples, MaxSamples);
    Tmat = repmat(T,[1,Nsamples]);

    % Simulate trials
    if iter == 1 && Nsamples == 1 && ~isempty(simdata0)
        simdata = simdata0;
        elapsed_time = elapsed_time0;
    else
        fun_clock = tic;
        if isempty(designMat)    % Pass only trial indices
            simdata = fun(params,Tmat(:),varargin{:});
            fc = fc + 1;
        else                    % Pass full design matrix per trial
            simdata = fun(params,designMat(Tmat(:),:),varargin{:});
            fc = fc + 1;
        end
        elapsed_time = toc(fun_clock);
    end

    % Check that the returned simulated data have the right size
    if size(simdata,1) ~= numel(Tmat)
        error('ibslike:SizeMismatch', ...
            'Number of rows of returned simulated data does not match the number of requested trials.');
    end

    Ns = Ns + Ttrials;

    % Accelerated sampling
    if options.Acceleration > 0 && elapsed_time < options.AccelerationThreshold
        samples_level = samples_level*options.Acceleration;
    end

    % Check new "hits"
    hits_temp = all(respMat(Tmat(:),:) == simdata,2);

    % Build matrix of new hits (sandwich with buffer of hits, then removed)
    hits_new = [ones(1,Ttrials);reshape(hits_temp,size(Tmat))';ones(1,Ttrials)];

    % Warning: from now on it's going to be incomprehensible
    % (all vectorized for speed)

    % Extract matrix of Ks from matrix of hits for this iteration
    h = size(hits_new,1);
    list = find(hits_new(:) == 1)-1;
    row = floor(list/h)+1;
    col = mod(list,h)+1;
    delta = diff([col;1]);
    remidx = delta <= 0;
    delta(remidx) = [];
    row(remidx) = [];
    indexcol = find(diff([0;row]));
    col = 1 + (1:numel(row))' - indexcol(row);
    K_iter = zeros(size(T,1),max(col));
    K_iter(row + (col-1)*size(K_iter,1)) = delta;

    % This is the comprehensible version that we want to get to:
    %
    %   for iTrial = 1:Ntrials
    %       index = find(hits_new(iTrial,:),targetHits(iTrial));
    %       K = diff([0 index]);
    %       logL(iTrial) = sum(Ktab(K))/numel(index);
    %   end

    % Add still-open K to first column
    K_iter(:,1) = K_open(T)' + K_iter(:,1);

    % Find last K position for each trial
    [~,idx_last] = min([K_iter,zeros(Ttrials,1)],[],2);
    idx_last = idx_last - 1;
    ii = sub2ind(size(K_iter),(1:Ttrials)',idx_last);

    % Subtract one hit from last K (it was added)
    K_iter(ii) = K_iter(ii) - 1;
    K_open(T) = K_iter(ii)';

    % For each trial, ignore entries of K_iter past max # of reps
    idx_mat = bsxfun(@plus,Ridx(T)',repmat(0:size(K_iter,2)-1,[Ttrials,1]));
    K_iter(idx_mat > (options.Nreps)) = 0;

    % Find last K position for each trial again
    [~,idx_last2] = min([K_iter,zeros(Ttrials,1)],[],2);
    idx_last2 = idx_last2 - 1;

    % Add current K to full K matrix
    K_iter_place = bsxfun(@ge,K_place0(:,1:Ttrials),Ridx(T)) & bsxfun(@le,K_place0(:,1:Ttrials),Ridx(T) + idx_last2'- 1);
    K_place = false(size(K_place0));
    K_place(:,T) = K_iter_place;
    Kt = K_iter';
    K_mat(K_place) = Kt(Kt > 0);
    Ridx(T) = Ridx(T) + idx_last' - 1;

    % Compute log-likelihood only if requested for thresholding
    if isfinite(options.NegLogLikeThreshold)
        Rmin = min(Ridx(T));    % Find repeat still ongoing
        if Rmin > size(K_mat,1); continue; end
        [LL_temp,Psi_tab] = get_LL_from_K(Psi_tab,K_mat(Rmin,:));
        nLL_temp = -sum(LL_temp,1);
        if nLL_temp > options.NegLogLikeThreshold
            idx_move = Ridx == Rmin;
            Ridx(idx_move) = Rmin+1;
            K_open(idx_move) = 0;
            exitflag = 1;
        end
    end
end

if ~isempty(T)
    error('ibslike:ConvergenceFail', ...
        'Maximum number of iterations or time limit reached and algorithm did not converge. Check FUN and DATA.');
end

% Log likelihood estimate per trial and run lengths K for each repetition
Nreps = sum(K_mat > 0,1)';
[LL_mat,Psi_tab] = get_LL_from_K(Psi_tab,K_mat);
nlogL = sum(-LL_mat',2)./Nreps;
K = K_mat';

end

%--------------------------------------------------------------------------
function [nlogL,K,Nreps,Ns,fc,exitflag] = loop_ibs_sampling(fun,params,respMat,designMat,simdata0,elapsed_time0,t0,options,varargin)

Ntrials = size(respMat,1);
Trials = (1:Ntrials)';
MaxIter = options.MaxIter;
exitflag = 0;

K = zeros(Ntrials,options.Nreps);
Ns = 0;
fc = 0;
Psi_tab = [];

for iRep = 1:options.Nreps

    offset = 1;
    hits = zeros(Ntrials,1);
    if isfinite(options.MaxTime) && toc(t0) > options.MaxTime
        exitflag = 2;
        break;
    end

    for iter = 1:MaxIter
        % Pick trials that need more hits, sample multiple times
        T = Trials(hits < 1);
        if isempty(T); break; end

        % Simulate trials
        if iter == 1 && iRep == 1 && ~isempty(simdata0)
            simdata = simdata0;
        elseif isempty(designMat)    % Pass only trial indices
            simdata = fun(params,T(:),varargin{:});
            fc = fc + 1;
        else                    % Pass full design matrix per trial
            simdata = fun(params,designMat(T(:),:),varargin{:});
            fc = fc + 1;
        end

        % Check that the returned simulated data have the right size
        if size(simdata,1) ~= numel(T)
            error('ibslike:SizeMismatch', ...
                'Number of rows of returned simulated data does not match the number of requested trials.');
        end

        Ns = Ns + numel(T); % Count samples
        hits_new = all(respMat(T(:),:) == simdata,2);
        hits(T) = hits(T) + hits_new;

        K(T(hits_new),iRep) = offset;
        offset = offset + 1;

        % Terminate if negative log likelihood is above a given threshold
        if isfinite(options.NegLogLikeThreshold)
            K(hits < 1,iRep) = offset;
            [LL_mat,Psi_tab] = get_LL_from_K(Psi_tab,K(:,iRep));
            nlogL_sum = -sum(LL_mat,1);
            if nlogL_sum > options.NegLogLikeThreshold
                T = [];
                exitflag = 1;
                break;
            end
        end

        % Terminate if above maximum allowed runtime
        if isfinite(options.MaxTime) && toc(t0) > options.MaxTime
            T = [];
            exitflag = 2;
            break;
        end

    end

    if ~isempty(T)
        error('ibslike:ConvergenceFail', ...
            'Maximum number of iterations reached and algorithm did not converge. Check FUN and DATA.');
    end
end

Nreps = sum(K > 0,2);
[LL_mat,Psi_tab] = get_LL_from_K(Psi_tab,K);
nlogL = sum(-LL_mat,2)./Nreps;

end

%--------------------------------------------------------------------------
function [LL_mat,Psi_tab] = get_LL_from_K(Psi_tab,K_mat)
%GET_LL_FROM_K Convert matrix of K values into log-likelihoods.

K_max = max(1,max(K_mat(:)));
if K_max > numel(Psi_tab)   % Fill digamma function table
    Psi_tab = [Psi_tab; (psi(1) - psi(numel(Psi_tab)+1:K_max)')];
end
LL_mat = Psi_tab(max(1,K_mat));

end

%--------------------------------------------------------------------------
function exitflag = runtest1(options)

Nreps = 1e3;
RMSE_tol = 2/sqrt(Nreps);

% Binomial probability model
p_model = exp(linspace(log(1e-3),log(1),10));
fun = @(x,dmat) rand(size(dmat)) < x;   % Simulating function
rmat = fun(1,NaN);

fprintf('\n');
fprintf('TEST 1: Using IBS to compute log(p) of Bernoulli distributions with %d repeats.\n',Nreps);
fprintf('We consider p = %s.\n',mat2str(p_model,3));

options.Nreps = Nreps;
options.NegLogLikeThreshold = Inf;

nlogL = zeros(1,numel(p_model));
nlogLvar = zeros(1,numel(p_model));
for iter = 1:numel(p_model)
    [nlogL(iter),nlogLvar(iter)] = ibslike(fun,p_model(iter),rmat,[],options);
end

% We expect the true value to be almost certainly (> 99.99%) in this range
LL_min = (-nlogL - 4*sqrt(nlogLvar));
LL_max = (-nlogL + 4*sqrt(nlogLvar));

exitflag = any(log(p_model) < LL_min) | any(log(p_model) > LL_max);

rmse = sqrt(mean((-nlogL - log(p_model)).^2));
fprintf('Average RMSE of log(p) estimates across p: %.4f.\n',rmse);

exitflag = exitflag | (rmse > RMSE_tol);

if exitflag
    fprintf('Test FAILED. Something might be wrong.\n');
else
    fprintf('Test PASSED. IBS estimates are calibrated and close to ground truth.\n');
end

% Plot figure
xx = log(p_model);
h(1) = plot(xx,xx,'k-','LineWidth',2); hold on;

yy = -nlogL;
xxerr = [xx, fliplr(xx)];
yyerr_down = yy - 1.96*sqrt(nlogLvar);
yyerr_up = yy + 1.96*sqrt(nlogLvar);
yyerr = [yyerr_down, fliplr(yyerr_up)];
fill(xxerr, yyerr,'b','FaceAlpha',0.5,'LineStyle','none'); hold on;
h(2) = plot(xx,yy,'b-','LineWidth',2); hold on;

box off;
set(gca,'TickDir','out');
set(gcf,'Color','w');
xlabel('True log(p)');
ylabel('Estimated log(p)')
%xlim([-5 5]);
hl = legend(h,'True log(p)','IBS estimate (95% CI)');
set(hl,'Location','NorthWest','Box','off');
title('IBS estimation test');

end

%--------------------------------------------------------------------------
function exitflag = runtest2(options)

% Binomial probability model
Ntrials = 100;
p_true = 0.9*rand() + 0.05;             % True probability
p_model = 0.9*rand() + 0.05;            % Model probability
fun = @(x,dmat) rand(size(dmat)) < x;   % Simulating function
Nexps = 2e3;

options.NegLogLikeThreshold = Inf;

fprintf('\n');
fprintf('TEST 2: Using IBS to compute the log-likelihood of a binomial distribution.\n');
fprintf('Parameters: p_true=%.2g, p_model=%.2g, %d trials per experiment.\n',p_true,p_model,Ntrials);
fprintf('The distribution of z-scores should approximate a standard normal distribution (mean 0, SD 1).\n');

zscores = zeros(1,Nexps);
for iter = 1:Nexps
    rmat = fun(p_true,NaN(Ntrials,1));            % Generate data
    [nlogL,nlogLvar] = ibslike(fun,p_model,rmat,[],options);
    nlogL_exact = -log(p_model)*sum(rmat == 1) - log(1-p_model)*sum(rmat == 0);
    zscores(iter) = (nlogL_exact - nlogL)/sqrt(nlogLvar);
end

edges = -4.75:0.5:4.75;
nz = histc(zscores,edges);
h(1) = bar(edges,nz,'histc');
hold on;
xx = linspace(-5,5,1e4);
h(2) = plot(xx,Nexps*exp(-xx.^2/2)/sqrt(2*pi)/2,'k-','LineWidth',2);

box off;
set(gca,'TickDir','out');
set(gcf,'Color','w');
xlabel('z-score');
ylabel('pdf')
xlim([-5 5]);
hl = legend(h,'z-scores histogram','expected pdf');
set(hl,'Location','NorthEast','Box','off');
title('Calibration test');

exitflag = abs(mean(zscores)) > 0.15 || abs(std(zscores) - 1) > 0.1;

fprintf('Distribution of z-scores (%d experiments). Mean: %.4g. Standard deviation: %.4g.\n',Nexps,mean(zscores),std(zscores));
if exitflag
    fprintf('Test FAILED. Something might be wrong.\n');
else
    fprintf('Test PASSED. We verified that IBS is unbiased (~zero mean) and calibrated (SD ~1).\n');
end

end

%--------------------------------------------------------------------------
function exitflag = runtest3(options)

Nreps = 100;
RMSE_tol = 4/sqrt(Nreps);

% Binomial probability model
p_model = exp(linspace(log(1e-3),log(0.1),10));
fun = @(x,dmat) rand(size(dmat)) < x;   % Simulating function
rmat = fun(1,NaN);
thresh = -log(0.01);
p_target = max(p_model,exp(-thresh));

fprintf('\n');
fprintf('TEST 3: Log-likelihood thresholding at log(p) = %.3f.\n',-thresh);
fprintf('Using IBS to compute thresholded log(p) of Bernoulli distributions with %d repeats.\n',Nreps);
fprintf('We consider p = %s.\n',mat2str(p_model,3));

options.Nreps = Nreps;
options.NegLogLikeThreshold = thresh;
options.Acceleration = 1;

nlogL = zeros(1,numel(p_model));
nlogLvar = zeros(1,numel(p_model));
for iter = 1:numel(p_model)
    [nlogL(iter),nlogLvar(iter)] = ibslike(fun,p_model(iter),rmat,[],options);
end

% We expect the true value to be almost certainly (> 99.99%) in this range
LL_min = (-nlogL - 4*sqrt(nlogLvar));
LL_max = (-nlogL + 4*sqrt(nlogLvar));

% We expect the estimates to be (almost) correct away from the threshold
idx = log(p_model) > -thresh*0.75;
exitflag = any(log(p_model(idx)) < LL_min(idx)) | any(log(p_model(idx)) > LL_max(idx));

% We expect the estimates to be above the true value below the threshold
LL_thresh = (-nlogL - sqrt(nlogLvar));
idx_below = log(p_model) < -thresh;
exitflag = exitflag | any(log(p_model(idx_below)) > LL_thresh(idx_below));
% exitflag = any(log(p_target) < LL_min) | any(log(p_target) > LL_max);

rmse = sqrt(mean((-nlogL(idx) - log(p_target(idx))).^2));
fprintf('Average RMSE of log(p) estimates across p: %.4f.\n',rmse);

exitflag = exitflag | (rmse > RMSE_tol);

if exitflag
    fprintf('Test FAILED. Something might be wrong.\n');
else
    fprintf('Test PASSED. IBS estimates are calibrated and close to (thresholded) ground truth.\n');
end

% Plot figure
xx = log(p_model);
h(1) = plot(xx,xx,'k-','LineWidth',2); hold on;

yy = -nlogL;
xxerr = [xx, fliplr(xx)];
yyerr_down = yy - 1.96*sqrt(nlogLvar);
yyerr_up = yy + 1.96*sqrt(nlogLvar);
yyerr = [yyerr_down, fliplr(yyerr_up)];
fill(xxerr, yyerr,'b','FaceAlpha',0.5,'LineStyle','none'); hold on;

h(2) = plot(xx,yy,'b-','LineWidth',2); hold on;
h(3) = plot([xx(1),xx(end)],-thresh*[1 1],'k:','LineWidth',2);

box off;
set(gca,'TickDir','out');
set(gcf,'Color','w');
xlabel('True log(p)');
ylabel('Estimated log(p) with thresholding')
%xlim([-5 5]);
hl = legend(h,'True log(p)','IBS estimate (95% CI)','Threshold');
set(hl,'Location','NorthWest','Box','off');
title('Thresholded IBS test');

end

%   TODO:
%   - Fix help and documentation
%   - Optimal allocation of estimates?

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
 function y = mtimes_block(ob, x, iblock, nblock)
%function y = mtimes_block(ob, x, iblock, nblock)
% y = G(i'th block) * x	or y = G'(i'th block) * x
% in either case the project data will be "small"
% iblock is 1,...,nblock

% support 'exists' option for seeing if this routine is available
if nargin == 2 & ischar(x) & streq(x, 'exists')
	y = 1;
	return
end

if nargin ~= 4, ir_usage, end

nb = ob.nb;
na = ob.na;

%
% forward projection
%
if ~ob.is_transpose
	if ob.apower ~= 1, error notdone, end

	% if needed, expand concise column
	if ob.is_masked
		idim = size(x);
		np = sum(ob.mask(:));
		nxy = numel(ob.mask);
		if idim(1) ~= np
			error 'size mismatch'
		end
		x = embed(x, ob.mask);
		x = reshape(x, nxy, idim(2));
	end

	y = double(wtfmex('dsc,proj', ob.arg', single(x), uint8(ob.mask), ...
		int32(iblock-1), int32(nblock), ...
		int32(ob.nthread), int32(ob.chat)));

	% fix: extract the relevent columns - should do in wtfmex?
	ia = iblock:nblock:na;
	nv = length(ia);

	if ndims(y) == 3			% [nb,na,nz]
		y = y(:,ia,:);			% [nb,nv,nz]

	elseif size(y,1) == nb*na		% [nba,nz]
		nz = size(y,2);
		y = reshape(y, nb, na, nz);	% [nb,na,nz]
		y = y(:,ia,:);			% [nb,nv,nz]
		y = reshape(y, nb*nv, nz);	% [nb*nv,nz]

	elseif size(y,1) == nb			% [nb,na]
		y = y(:,ia);			% [nb,nv]

	else
		error size
	end

%
% backprojection
%
else
	if ob.apower ~= 1, error notdone, end

	ia = iblock:nblock:na;
	nv = length(ia);

	if ndims(x) == 3		% [nb,nv,nz]
		error todo

	elseif size(x,1) == nb*nv	% [nb*nv,nz]
		nz = size(x,2);
		tmp = zeros(nb, na, nz);
		tmp(:,ia,:) = reshape(x, [nb nv nz]);
		x = reshape(tmp, [nb*na nz]);

	elseif size(x,1) == nb		% [nb,nv]
		error todo

	else
		error bug
	end

	y = double(wtfmex('dsc,back', ob.arg', single(x), uint8(ob.mask), ...
		int32(iblock-1), int32(nblock), ...
		int32(ob.nthread), int32(ob.chat)));

	if ob.is_masked
		if size(y,1) == numel(ob.mask)
			y = y(ob.mask,:);	% [nxy,nz] -> [np,nz]
		else
			y = y(ob.mask);		% [nx,ny] -> [np]
		end
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function J = computeCost(X, y, theta)
%COMPUTECOST Compute cost for linear regression
%   J = COMPUTECOST(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = 1 / (2 * m) * (X * theta - y)' * (X * theta - y);

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

end

repo: machine-learning