text:
% PAC - compute phase-amplitude coupling (power of first input
%         correlation with phase of second). There is no graphical output
%         to this function.
%
% Usage:
%   >> pac(x,y,srate);
%   >> [coh,timesout,freqsout1,freqsout2,cohboot] ...
%                     = pac(x,y,srate,'key1', 'val1', 'key2', val2' ...);
% Inputs:
%    x       = [float array] 2-D data array of size (times,trials) or
%              3-D (1,times,trials)
%    y       = [float array] 2-D or 3-d data array
%    srate   = data sampling rate (Hz)
%
%    Most important optional inputs
%       'method'    = ['mod'|'corrsin'|'corrcos'|'latphase'] modulation
%                     method or correlation of amplitude with sine or cosine of
%                     angle (see ref). 'laphase' compute the phase
%                     histogram at a specific time and requires the
%                     'powerlat' option to be set.
%       'freqs'     = [min max] frequency limits. Default [minfreq 50],
%                     minfreq being determined by the number of data points,
%                     cycles and sampling frequency. Use 0 for minimum frequency
%                     to compute default minfreq. You may also enter an
%                     array of frequencies for the spectral decomposition
%                     (for FFT, closest computed frequency will be returned; use
%                     'padratio' to change FFT freq. resolution).
%       'freqs2'    = [float array] array of frequencies for the second
%                     argument. 'freqs' is used for the first argument.
%                     By default it is the same as 'freqs'.
%       'wavelet'   = 0  -> Use FFTs (with constant window length) { Default }
%                   = >0 -> Number of cycles in each analysis wavelet
%                   = [cycles expfactor] -> if 0 < expfactor < 1,  the number
%                     of wavelet cycles expands with frequency from cycles
%                     If expfactor = 1, no expansion; if = 0, constant
%                     window length (as in FFT)            {default wavelet: 0}
%                   = [cycles array] -> cycle for each frequency. Size of array
%                      must be the same as the number of frequencies
%                     {default cycles: 0}
%       'wavelet2'  = same as 'wavelet' for the second argument. Default is
%                     same as cycles. Note that if the lowest frequency for X
%                     and Y are different and cycle is [cycles expfactor], it
%                     may result in discrepancies in the number of cycles at
%                     the same frequencies for X and Y.
%       'ntimesout' = Number of output times (int<frames-winframes). Enter a
%                     negative value [-S] to subsample original time by S.
%       'timesout'  = Enter an array to obtain spectral decomposition at
%                     specific time values (note: algorithm find closest time
%                     point in data and this might result in an unevenly spaced
%                     time array). Overwrite 'ntimesout'. {def: automatic}
%       'powerlat'  = [float] latency in ms at which to compute phase
%                     histogram
%       'tlimits'   = [min max] time limits in ms.
%
%    Optional Detrending:
%       'detrend'   = ['on'|'off'], Linearly detrend each data epoch   {'off'}
%       'rmerp'     = ['on'|'off'], Remove epoch mean from data epochs {'off'}
%
%    Optional FFT/DFT Parameters:
%       'winsize'   = If cycles==0: data subwindow length (fastest, 2^n<frames);
%                     If cycles >0: *longest* window length to use. This
%                     determines the lowest output frequency. Note that this
%                     parameter is overwritten if the minimum frequency has been set
%                     manually and requires a longer time window {~frames/8}
%       'padratio'  = FFT-length/winframes (2^k)                    {2}
%                     Multiplies the number of output frequencies by dividing
%                     their spacing (standard FFT padding). When cycles~=0,
%                     frequency spacing is divided by padratio.
%       'nfreqs'    = number of output frequencies. For FFT, closest computed
%                     frequency will be returned. Overwrite 'padratio' effects
%                     for wavelets. Default: use 'padratio'.
%       'freqscale' = ['log'|'linear'] frequency scale. Default is 'linear'.
%                     Note that for obtaining 'log' spaced freqs using FFT,
%                     closest correspondent frequencies in the 'linear' space
%                     are returned.
%       'subitc'    = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence
%                    (ITC) from x and y. This computes the  'intrinsic' coherence
%                     x and y not arising from common synchronization to
%                     experimental events. See notes. {default: 'off'}
%       'itctype'   = ['coher'|'phasecoher'] For use with 'subitc', see TIMEF
%                     for more details {default: 'phasecoher'}.
%       'subwin'    = [min max] sub time window in ms (this windowing is
%                     performed after the spectral decomposition).
%
% Outputs:
%        pac         = Matrix (nfreqs1,nfreqs2,timesout) of coherence (complex).
%                      Use 20*log(abs(crossfcoh)) to visualize log spectral diffs.
%        timesout    = Vector of output times (window centers) (ms).
%        freqsout1   = Vector of frequency bin centers for first argument (Hz).
%        freqsout2   = Vector of frequency bin centers for second argument (Hz).
%        alltfX      = single trial spectral decomposition of X
%        alltfY      = single trial spectral decomposition of Y
%
% Author: Arnaud Delorme, SCCN/INC, UCSD 2005-
%
% Ref: Testing for Nested Oscilations (2008) J Neuro Methods 174(1):50-61
%
% See also: TIMEFREQ, CROSSF

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

function [crossfcoh, timesout1, freqs1, freqs2, crossfcohall, alltfX, alltfY] = pac(X, Y, srate, varargin);

if nargin < 1
    help pac;
    return;
end

% deal with 3-D inputs
% --------------------
if ndims(X) == 3, X = reshape(X, size(X,2), size(X,3)); end
if ndims(Y) == 3, Y = reshape(Y, size(Y,2), size(Y,3)); end
frame = size(X,2);

g = finputcheck(varargin, ...
                { 'alpha'         'real'     [0 0.2]                  [];
                  'baseboot'      'float'    []                       0;
                  'boottype'      'string'   {'times','trials','timestrials'}  'timestrials';
                  'detrend'       'string'   {'on','off'}              'off';
                  'freqs'         'real'     [0 Inf]                  [0 srate/2];
                  'freqs2'        'real'     [0 Inf]                  [];
                  'freqscale'     'string'   { 'linear','log' }       'linear';
                  'itctype'       'string'   {'phasecoher','phasecoher2','coher'}  'phasecoher';
                  'nfreqs'        'integer'  [0 Inf]                  [];
                  'lowmem'        'string'   {'on','off'}              'off';
                  'method'        'string'   { 'mod','corrsin','corrcos','latphase' }         'mod';
                  'naccu'         'integer'  [1 Inf]                   250;
                  'newfig'        'string'   {'on','off'}              'on';
                  'padratio'      'integer'  [1 Inf]                   2;
                  'rmerp'         'string'   {'on','off'}              'off';
                  'rboot'         'real'     []                        [];
                  'subitc'        'string'   {'on','off'}              'off';
                  'subwin'        'real'     []                        []; ...
                  'gammapowerlim' 'real'     []                        []; ...
                  'powerlim'      'real'     []                        []; ...
                  'powerlat'      'real'     []                        []; ...
                  'gammabase'     'real'     []                        []; ...
                  'timesout'      'real'     []                        []; ...
                  'ntimesout'     'integer'  []                        200; ...
                  'tlimits'       'real'     []                        [0 frame/srate];
                  'title'         'string'   []                        '';
                  'vert'          { 'real','cell' }  []                [];
                  'wavelet'       'real'     [0 Inf]                   0;
                  'wavelet2'      'real'     [0 Inf]                   [];
                  'winsize'       'integer'  [0 Inf]                   max(pow2(nextpow2(frame)-3),4) }, 'pac');

if ischar(g), error(g); end

% more defaults
% -------------
if isempty(g.wavelet2), g.wavelet2 = g.wavelet; end
if isempty(g.freqs2),   g.freqs2   = g.freqs;   end

% remove ERP if necessary
% -----------------------
X = squeeze(X);
Y = squeeze(Y);X = squeeze(X);
trials = size(X,2);
if strcmpi(g.rmerp, 'on')
    X = X - repmat(mean(X,2), [1 trials]);
    Y = Y - repmat(mean(Y,2), [1 trials]);
end

% perform timefreq decomposition
% ------------------------------
[alltfX freqs1 timesout1] = timefreq(X, srate, 'ntimesout',  g.ntimesout, 'timesout',  g.timesout,  'winsize',  g.winsize, ...
                                'tlimits', g.tlimits, 'detrend',   g.detrend,   'itctype',  g.itctype, ...
                                'subitc',  g.subitc,  'wavelet',   g.wavelet,   'padratio', g.padratio, ...
                                'freqs',   g.freqs,   'freqscale', g.freqscale, 'nfreqs',   g.nfreqs);
[alltfY freqs2 timesout2] = timefreq(Y, srate, 'ntimesout',  g.ntimesout, 'timesout',  g.timesout,  'winsize',  g.winsize, ...
                                'tlimits', g.tlimits, 'detrend',   g.detrend,   'itctype',  g.itctype, ...
                                'subitc',  g.subitc,  'wavelet',   g.wavelet2,  'padratio', g.padratio, ...
                                'freqs',   g.freqs2,  'freqscale', g.freqscale, 'nfreqs',   g.nfreqs);

% check time limits
% -----------------
if ~isempty(g.subwin)
    ind1      = find(timesout1 > g.subwin(1) & timesout1 < g.subwin(2));
    ind2      = find(timesout2 > g.subwin(1) & timesout2 < g.subwin(2));
    alltfX    = alltfX(:, ind1, :);
    alltfY    = alltfY(:, ind2, :);
    timesout1 = timesout1(ind1);
    timesout2 = timesout2(ind2);
end
if length(timesout1) ~= length(timesout2) || any( timesout1 ~= timesout2)
    disp('Warning: Time points are different for X and Y. Use ''timesout'' to specify common time points');
    [vals ind1 ind2 ] = intersect_bc(timesout1, timesout2);
    fprintf('Searching for common time points: %d found\n', length(vals));
    if length(vals) < 10, error('Less than 10 common data points'); end
    timesout1 = vals;
    timesout2 = vals;
    alltfX = alltfX(:, ind1, :);
    alltfY = alltfY(:, ind2, :);
end

% scan across frequency and time
% -------------------------------
%if isempty(g.alpha)
%    disp('Warning: if significance mask is not applied, result might be slightly')
%    disp('different (since angle is not made uniform and amplitude interpolated)')
%end

cohboot =[];
if ~strcmpi(g.method, 'latphase')
    for find1 = 1:length(freqs1)
        for find2 = 1:length(freqs2)
            for ti = 1:length(timesout1)

                % get data
                % --------
                tmpalltfx = squeeze(alltfX(find1,ti,:));
                tmpalltfy = squeeze(alltfY(find2,ti,:));

                %if ~isempty(g.alpha)
                %    tmpalltfy = angle(tmpalltfy);
                %    tmpalltfx = abs(  tmpalltfx);
                %    [ tmp cohboot(find1,find2,ti,:) newamp newangle ] = ...
                %        bootcircle(tmpalltfx, tmpalltfy, 'naccu', g.naccu);
                %    crossfcoh(find1,find2,ti) = sum ( newamp .* exp(j*newangle) );
                %else
                tmpalltfy = angle(tmpalltfy);
                tmpalltfx = abs(  tmpalltfx);
                if strcmpi(g.method, 'mod')
                    crossfcoh(find1,find2,ti) = sum( tmpalltfx .* exp(j*tmpalltfy) );
                elseif strcmpi(g.method, 'corrsin')
                    tmp = corrcoef( sin(tmpalltfy), tmpalltfx);
                    crossfcoh(find1,find2,ti) = tmp(2);
                else
                    tmp = corrcoef( cos(tmpalltfy), tmpalltfx);
                    crossfcoh(find1,find2,ti) = tmp(2);
                end
            end
        end
    end
elseif 1
    % this option computes power at a given latency
    % then computes the same as above (vectors)

    %if isempty(g.powerlat)
    %    error('You need to specify a latency for the ''powerlat'' option');
    %end

    gammapower = mean(10*log10(alltfX(:,:,:).*conj(alltfX)),1); % average all frequencies for power
    if isempty(g.gammapowerlim)
        g.gammapowerlim = [ min(gammapower(:)) max(gammapower(:)) ];
    end
    fprintf('Gamma power limits: %3.2f to %3.2f\n', g.gammapowerlim(1), g.gammapowerlim(2));
    power = 10*log10(alltfY(:,:,:).*conj(alltfY));
    if isempty(g.powerlim)
        for freq = 1:size(power,1)
            g.powerlim(freq,:) = [ min(power(freq,:)) max(power(freq,:)) ];
        end
    end
    for freq = 1:size(power,1)
        fprintf('Freq %d power limits: %3.2f to %3.2f\n', freqs2(freq), g.powerlim(freq,1), g.powerlim(freq,2));
    end

    % power plot
    %figure; plot(timesout2/1000, (mean(power(9,:,:),3)-mean(power(9,:)))/50);
    %hold on; plot(linspace(0, length(Y)/srate, length(Y)), mean(Y'), 'g');

    % phase with power
    % figure; plot(timesout2/1000, (mean(phaseangle(9,:,:),3)-mean(phaseangle(9,:)))/50);
    % hold on; plot(timesout1/1000, (mean(gammapower,3)-mean(gammapower(:)))/100, 'r');
    %figure; plot((mean(phaseangle(9,:,:),3)-mean(phaseangle(9,:)))/50+j*(mean(gammapower,3)-mean(gammapower(:)))/100, '.');

    matsize               = 32;
    matcenter             = (matsize-1)/2+1;
    matrixfinalgammapower = zeros(size(alltfY,1),size(alltfX,3),matsize,matsize);
    matrixfinalcount      = zeros(size(alltfY,1),size(alltfX,3),matsize,matsize);

    % get power indices
    if isempty(g.gammabase)
        g.gammabase = mean(gammapower(:));
    end
    fprintf('Gamma power average: %3.2f\n', g.gammabase);
    gammapoweradd  = gammapower-g.gammabase;
    gammapower     = floor((gammapower-g.gammapowerlim(1))/(g.gammapowerlim(2)-g.gammapowerlim(1))*(matsize-2))+1;
    phaseangle     = angle(alltfY);
    posx           = zeros(size(power));
    posy           = zeros(size(power));
    for freq = 1:length(freqs2)
        fprintf('Processing frequency %3.2f\n', freqs2(freq));
        power(freq,:,:) = (power(freq,:,:)-g.powerlim(freq,1))/(g.powerlim(freq,2)-g.powerlim(freq,1))*(matsize-3)/2+1;
        complexval      = power(freq,:,:).*exp(j*phaseangle(freq,:,:));
        posx(freq,:,:)  = round(real(complexval)+matcenter);
        posy(freq,:,:)  = round(imag(complexval)+matcenter);
        for trial = 1:size(alltfX,3) % scan trials
            for time = 1:size(alltfX,2)
                %matrixfinal(freq,posx(freq,time,trial),posy(freq,time,trial),gammapower(1,time,trial)) = ...
                %    matrixfinal(freq,posx(freq,time,trial),posy(freq,time,trial),gammapower(1,time,trial))+1;
                matrixfinalgammapower(freq,trial,posx(freq,time,trial),posy(freq,time,trial)) = ...
                    matrixfinalgammapower(freq,trial,posx(freq,time,trial),posy(freq,time,trial))+gammapoweradd(1,time,trial);
                matrixfinalcount(freq,trial,posx(freq,time,trial),posy(freq,time,trial)) = ...
                    matrixfinalcount(freq,trial,posx(freq,time,trial),posy(freq,time,trial))+1;
            end
        end
        %matrixfinal(freq,:,:,:) = convn(squeeze(matrixfinal(freq,:,:,:)), gs, 'same');
        %tmpmat = posx(index,:)+(posy(index,:)-1)*64+(gammapower(:)-1)*64*64;
        matrixfinalcount(freq, find(matrixfinalcount(freq,:) == 0)) = 1;
        matrixfinalgammapower(freq,:,:,:) = matrixfinalgammapower(freq,:,:,:)./matrixfinalcount(freq,:,:,:);
    end

    % average and smooth
    matrixfinalgammapowermean = squeeze(mean(matrixfinalgammapower,2));
    for freq = 1:length(freqs2)
        matrixfinalgammapowermean(freq,:,:) = conv2(squeeze(matrixfinalgammapowermean(freq,:,:)), gauss2d(5,5), 'same');
    end
    %matrixfinalgammapower = matrixfinalgammapower/size(alltfX,3)/size(alltfX,2);

    %vect = linspace(-pi,pi,50);
    %for f = 1:length(freqs2)
    %    crossfcoh(f,:) = hist(tmpalltfy(f,:), vect);
    %end

    % smoothing of output image
    % -------------------------
    %gs = gauss2d(6, 6, 6);
    %crossfcoh = convn(crossfcoh, gs, 'same');
    %freqs1    = freqs2;
    %timesout1 = linspace(-180, 180, size(crossfcoh,2));

    crossfcoh    = matrixfinalgammapowermean;
    crossfcohall = matrixfinalgammapower;

else
    % this option computes power at a given latency
    % then computes the same as above (vectors)

    %if isempty(g.powerlat)
    %    error('You need to specify a latency for the ''powerlat'' option');
    %end

    gammapower = mean(10*log10(alltfX(:,:,:).*conj(alltfX)),1); % average all frequencies for power
    if isempty(g.gammapowerlim)
        g.gammapowerlim = [ min(gammapower(:)) max(gammapower(:)) ];
    end
    power = 10*log10(alltfY(:,:,:).*conj(alltfY));
    if isempty(g.powerlim)
        for freq = 1:size(power,1)
            g.powerlim(freq,:) = [ min(power(freq,:)) max(power(freq,:)) ];
        end
    end

    % power plot
    %figure; plot(timesout2/1000, (mean(power(9,:,:),3)-mean(power(9,:)))/50);
    %hold on; plot(linspace(0, length(Y)/srate, length(Y)), mean(Y'), 'g');

    % phase with power
    % figure; plot(timesout2/1000, (mean(phaseangle(9,:,:),3)-mean(phaseangle(9,:)))/50);
    % hold on; plot(timesout1/1000, (mean(gammapower,3)-mean(gammapower(:)))/100, 'r');
    %figure; plot((mean(phaseangle(9,:,:),3)-mean(phaseangle(9,:)))/50+j*(mean(gammapower,3)-mean(gammapower(:)))/100, '.');

    matsize               = 64;
    matcenter             = (matsize-1)/2+1;
    matrixfinal           = zeros(size(alltfY,1),64,64,64);
    matrixfinalgammapower = zeros(size(alltfY,1),matsize,matsize);
    matrixfinalcount      = zeros(size(alltfY,1),matsize,matsize);

    % get power indices
    gammapoweradd  = gammapower-mean(gammapower(:));
    gammapower     = floor((gammapower-g.gammapowerlim(1))/(g.gammapowerlim(2)-g.gammapowerlim(1))*(matsize-1))+1;
    phaseangle     = angle(alltfY);
    posx           = zeros(size(power));
    posy           = zeros(size(power));
    gs             = gauss3d(6, 6, 6);
    for freq = 1:size(alltfY)
        fprintf('Processing frequency %3.2f\n', freqs2(freq));
        power(freq,:,:) = (power(freq,:,:)-g.powerlim(freq,1))/(g.powerlim(freq,2)-g.powerlim(freq,1))*(matsize-2)/2;
        complexval      = power(freq,:,:).*exp(j*phaseangle(freq,:,:));
        posx(freq,:,:)  = round(real(complexval)+matcenter);
        posy(freq,:,:)  = round(imag(complexval)+matcenter);
        for trial = 1:size(alltfX,3) % scan trials
            for time = 1:size(alltfX,2)
                %matrixfinal(freq,posx(freq,time,trial),posy(freq,time,trial),gammapower(1,time,trial)) = ...
                %    matrixfinal(freq,posx(freq,time,trial),posy(freq,time,trial),gammapower(1,time,trial))+1;
                matrixfinalgammapower(freq,posx(freq,time,trial),posy(freq,time,trial)) = ...
                    matrixfinalgammapower(freq,posx(freq,time,trial),posy(freq,time,trial))+gammapoweradd(1,time,trial);
                matrixfinalcount(freq,posx(freq,time,trial),posy(freq,time,trial)) = ...
                    matrixfinalcount(freq,posx(freq,time,trial),posy(freq,time,trial))+1;
            end
        end
        %matrixfinal(freq,:,:,:) = convn(squeeze(matrixfinal(freq,:,:,:)), gs, 'same');
        %tmpmat = posx(index,:)+(posy(index,:)-1)*64+(gammapower(:)-1)*64*64;
        matrixfinalcount(freq, find(matrixfinalcount(freq,:) == 0)) = 1;
        matrixfinalgammapower(freq,:,:) = matrixfinalgammapower(freq,:,:)./matrixfinalcount(freq, :,:);
        matrixfinalgammapower(freq,:,:) = conv2(squeeze(matrixfinalgammapower(freq,:,:)), gauss2d(5,5), 'same');
    end
    %matrixfinalgammapower = matrixfinalgammapower/size(alltfX,3)/size(alltfX,2);

    %vect = linspace(-pi,pi,50);
    %for f = 1:length(freqs2)
    %    crossfcoh(f,:) = hist(tmpalltfy(f,:), vect);
    %end

    % smoothing of output image
    % -------------------------
    %gs = gauss2d(6, 6, 6);
    %crossfcoh = convn(crossfcoh, gs, 'same');
    %freqs1    = freqs2;
    %timesout1 = linspace(-180, 180, size(crossfcoh,2));

    crossfcoh = matrixfinalgammapower;

end

% 7/31/2014 Ramon: crossfcohall sometimes does not exist depending on choice of input options
if ~exist('crossfcohall', 'var')
    crossfcohall = [];
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef MimModelCacheEntry < CoreBaseClass

    properties
        CurrentHash
        RemoteHash
    end

    methods
        function obj = MimModelCacheEntry(currentHash, remoteHash)
            obj.CurrentHash = currentHash;
            obj.RemoteHash = remoteHash;
        end

        function updateHashes(obj, currentHash, remoteHash)
            obj.CurrentHash = currentHash;
            obj.RemoteHash = remoteHash;
        end

        function changed = hasChanged(obj)
            changed = ~isequal(obj.CurrentHash, obj.RemoteHash);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef LIRCMOP14 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP with large infeasible regions

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Huang, Y. Fang, Y. You, J. Mo, C. Wei, and E.
% Goodman, An improved epsilon constraint-handling method in MOEA/D for
% CMOPs with large infeasible regions, Soft Computing, 2019, 23:
% 12491-12510.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            [popsize,variable_length] = size(X);
            sum1 = zeros(popsize,1);
            for j = 3 : variable_length
                sum1 = sum1+10*(X(:,j)-0.5).^2;
            end
            PopObj(:,1) = (1.7057+sum1).*cos(0.5*pi*X(:,1)).*cos(0.5*pi*X(:,2));
            PopObj(:,2) = (1.7057+sum1).*cos(0.5*pi*X(:,1)).*sin(0.5*pi*X(:,2));
            PopObj(:,3) = (1.7057+sum1).*sin(0.5*pi*X(:,1));
            gx          = PopObj(:,1).^2+PopObj(:,2).^2+PopObj(:,3).^2;
            PopCon(:,1) = (gx-9).*(4-gx);
            PopCon(:,2) = (gx-3.61).*(3.24-gx);
            PopCon(:,3) = (gx-3.0625).*(2.56-gx);
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R = sqrt(3.0625)*R./repmat(sqrt(sum(R.^2,2)),1,3);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            a = linspace(0,pi/2,10)';
            R = {sin(a)*cos(a')*sqrt(3.0625),sin(a)*sin(a')*sqrt(3.0625),cos(a)*ones(size(a'))*sqrt(3.0625)};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function model = standard_classifier_train(X,Y,T,options,binsize)
% Fit a classifier to the data given in X with labels Y of the type
% specified in the options - see below for details
%
% INPUT
% X: Brain data, (total time by regions) or (time by trials by regions)
% Y: Stimulus, (total time by q). This variable can either be entered as a
%       single vector of class assignments (ie q=1 and y takes values
%       1,2,...N=total number of classes); or can be entered as a matrix of
%       labels where each column contains 1 or 0 denoting whether that
%       class is active.
% T: Length of series
%       options             structure with the following subfields:
%           classifier      String indicating calssifer type;
%                   logistic    Logistic Regression classifier (default)
%                   SVM         Linear support vector machine classifier
%                   LDA         Linear Discriminant Analysis classifier
%           regularisation  String indicating type of regularisation;
%                   L1          Lasso regularisation
%                   L2          Ridge regularisation (default)
%                   ARD         Automatic relevance determination
%           lambda          Vector of regularisation penalty values to be
%                           tested
% binsize: NOT YET IMPLMENTED
%           how many consecutive time points will be used for each estiamtion.
%           By default, this is 1, such that one decoding model
%               is estimated using 1 time point
%
% OUTPUT
%
%
% Author: Cam Higgins, OHBA, University of Oxford (2019)

if nargin < 4 || isempty(options), options = struct(); end
if nargin < 5, binsize = 1; end

if ~all(T==T(1)), error('All elements of T must be equal for time-aligned classification'); end
N = length(T); ttrial = T(1);

if size(Y,1) < size(Y,2); Y = Y'; end
q = size(Y,2);
if size(Y,1) == length(T) % one value per trial
    responses = Y;
    Ystar = Y;
elseif length(Y(:)) ~= (ttrial*N*q)
    error('Incorrect dimensions in Y')
else
    responses = reshape(Y,[ttrial N q]);
    Ystar = reshape(repmat(responses(1,:,:),[ttrial,1,1]),[ttrial*N q]);
    responses = permute(responses(1,:,:),[2 3 1]); % N x q
    if ~all(Y(:)==Ystar(:))
        error('For time-aligned classification, the same stimulus must be presented for the entire trial');
    end
end

Ycopy = Y;
if size(Ycopy,1) == N
    Ycopy = repmat(reshape(Ycopy,[1 N q]),[ttrial 1 1]);
else
    Ycopy = reshape(Ycopy,[ttrial N q]);
end
% no demeaning by default if this is a classification problem
if ~isfield(options,'demeanstim'), options.demeanstim = 0; end

options.Nfeatures = 0;
options.K = 1;
[X,Y,T] = preproc4hmm(X,Y,T,options);
Q_star = size(Y,2); % for multinomial and LDA (mean inserted) models

p = size(X,2);

X = reshape(X,[ttrial N p]);
Y = reshape(Y,[ttrial N Q_star]);

if (binsize/2)==0
    warning(['binsize must be an odd number, setting to ' num2str(binsize+1)])
    binsize = binsize + 1;
end

if any(T~=ttrial), error('All trials must have the same length'); end
if nargin<5, binsize = 1; end

if isfield(options,'classifier')
    classifier = options.classifier; options = rmfield(options,'classifier');
    if ~(strcmp(classifier,'logistic') || strcmp(classifier,'SVM') ||...
            strcmp(classifier,'LDA') || strcmp(classifier,'KNN') || ...
            strcmp(classifier,'decisiontree') || strcmp(classifier,'SVM_rbf'))
        error('options.classifier can only take values "logistic", "SVM" or "LDA"');
    end
else, classifier = 'logistic';
end
model = struct();
model.classifier = classifier;
if isfield(options,'lambda')
    lambda = options.lambda; options = rmfield(options,'lambda');
else, lambda = 0.001;
end
if isfield(options,'regularisation')
    regtype = options.regularisation; options = rmfield(options,'regularisation');
    if ~(strcmp(regtype,'L1') || strcmp(regtype,'L2') || strcmp(regtype,'ARD'));
        error('options.regularisation can only take values "L1", "L2" or "ARD"');
    end
else, regtype = 'L2';
end

halfbin = floor(binsize/2);
%nwin = round(ttrial / binsize);
%binsize = floor(ttrial / nwin);
%RidgePen = lambda * eye(p);

%reshape to 3d vectors:
if strcmp(classifier,'logistic')
    % alternative: model = synchronousLogReg(X,Y,T,options)
    if length(unique(Ystar))<=2
        Y(Y==-1) = 0;
    end
    for t  = 1:ttrial
        if options.verbose
            if strcmp(regtype,'L1')
                fprintf(['\nRunning LASSO Regression analysis on data for t=',int2str(t),'th sample']);
            elseif strcmp(regtype,'L2')
                fprintf(['\nRunning RIDGE Regression analysis on data for t=',int2str(t),'th sample']);
            elseif strcmp(regtype,'ARD')
                error('ARD regularisation not implemented yet for logistic regression');
            end
        end
        for ilambda = 1:length(lambda)
            betas = zeros(p,Q_star);
            intercepts = zeros(1,Q_star);
            for c = 1:q
                if length(unique(Ystar))>2
                    vp = Y(t,:,c)~=0;
                    Ysample = Y(t,vp,c)'==1;
                else
                    vp = true(1,N);
                    Ysample = Y(t,vp,c)';
                end
                if strcmp(regtype,'L1')
                    [betas(:,c),fitInfo] = lassoglm(squeeze(X(t,vp,:)),Ysample, ...
                        'binomial','Alpha', 1, 'Lambda', lambda(ilambda),'Standardize', false);
                    intercepts(c) = fitInfo.Intercept;
                elseif strcmp(regtype,'L2')
                    if p>1
                        [betas(:,c),fitInfo] = lassoglm(squeeze(X(t,vp,:)),Ysample, ...
                            'binomial','Alpha', 0.00001, 'Lambda', lambda(ilambda),'Standardize', false);
                    else
                        [betas(:,c),fitInfo] = lassoglm(squeeze(X(t,vp,:))',Ysample, ...
                            'binomial','Alpha', 0.00001, 'Lambda', lambda(ilambda),'Standardize', false);
                    end
                    intercepts(c) = fitInfo.Intercept;
                end
            end
            model.betas(t,:,:,ilambda) = betas;
            model.intercepts(t,:,ilambda) = intercepts;
        end
        model.lambda = lambda;
    end
elseif strcmp(classifier,'SVM')
    model.SVM=cell(ttrial,Q_star);
    for t=1:ttrial
        if options.verbose;
            fprintf(['Inferring SVM (linear) classifiers for t = ',int2str(t),'th sample\n']);
        end
        for iStim=1:Q_star
            model.SVM{t,iStim} = fitclinear(squeeze(X(t,:,:)),squeeze(Y(t,:,iStim)));%,ones(sum(tselect(:,t))));
        end
    end
elseif strcmp(classifier,'SVM_rbf')
    model.SVM=cell(ttrial,Q_star);
    for t=1:ttrial
        if options.verbose;
            fprintf(['Inferring SVM (w radial basis function) classifiers for t = ',int2str(t),'th sample\n']);
        end
        for iStim=1:Q_star
            model.SVM{t,iStim} = fitcsvm(squeeze(X(t,:,:)),squeeze(Y(t,:,iStim)));%,ones(sum(tselect(:,t))));
        end
    end
elseif strcmp(classifier,'LDA')
    if isfield(options,'CVmethod')
        options_LDA = rmfield(options,'CVmethod');
    elseif isfield(options,'accuracyType')
        options_LDA = rmfield(options,'accuracyType');
    end
    options_LDA.useUnsupervisedGamma = 1; %this ensures gamma is fixed for standard decoding
    options_LDA.verbose = 1;
    options_LDA.useParallel = 0;
    options_LDA.standardise = 0;
    options_LDA.K = ttrial;
    options_LDA.Gamma = repmat(eye(ttrial),length(T),1);
    if isfield(options,'covtype');
        options_LDA.covtype = options.covtype;
    else
        options_LDA.covtype = 'full';
    end
    X_LDA = reshape(X,[ttrial*N p]);
    Y_LDA = reshape(Y,[ttrial*N Q_star]);
    if options.verbose
        fprintf(['Inferring LDA classifiers for all timepoints \n']);
    end
    [model,Gamma] = encodertrain(X_LDA,Y_LDA,T,options_LDA);
    model.Gamma = Gamma(1:ttrial,:);
    %model=rmfield(model,{'prior','Dir_alpha','Dir2d_alpha','P','Pi','features','K'});
    model.classifier='LDA';
elseif strcmp(classifier,'KNN')
    model.X_train = X;
    model.Y_train = Y;
elseif strcmp(classifier,'decisiontree')
    for t=1:ttrial
        Y_temp=zeros(N,1);
        for n=1:N;Y_temp(n) = find(squeeze(Y(t,n,:)));end
        model.decisiontree{t} = fitctree(squeeze(X(t,:,:)),Y_temp);
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function FOut = ReOrderVariables( FIn, NewVarOrder )

    %Check that new variable order is OK:
    [isMem MAP] = ismember( NewVarOrder, FIn.var );
    assert( length(FIn.var) == length(NewVarOrder), 'Different # variables!' );
    assert( all(isMem), 'Some of new variables not in factor!' );
    assert( length(NewVarOrder) == length(unique(NewVarOrder)), 'Some variables repeated!');

    %Set up variables and cardinality; preallocate
    FOut.var = NewVarOrder;
    FOut.card = FIn.card( MAP );
    FOut.val = ones(1,prod(FIn.card));

    for i=1:prod(FIn.card)
        %Loop through every value in FIn, and copy to FOut as approprate
        OldAssignment = IndexToAssignment(i, FIn.card);
        NewAssignment = OldAssignment( MAP );

        newIdx = AssignmentToIndex( NewAssignment, FOut.card );
        FOut.val( newIdx ) = FIn.val( i );
    end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function result=qn(x)

%QN is a scale estimator which does not require an auxiliary location estimate.
% Essentially it is the first quartile of all pairwise distances between
% two data points. Its definition is given by
%      qn(x)= c_n 2.2219{|x_i - x_j|; i<j}_(k)
% with c_n a small sample correction factor to make the qn unbiased at the
% normal distribution. It can resist 50% outliers.
% If x is a matrix, the scale estimate is computed on the columns of x. The
% result is then a row vector. If x is a column vector, the output is a scalar.
%
% Reference:
%   Rousseeuw P.J., Croux C. (1993),
%   "Alternatives to the median absolute deviation",
%   Journal of the American Statistical Association, 88, 1273-1283.
%
% Required input argument:
%    x: either a data matrix with n observations in rows, p variables in columns
%       or a column vector of length n.
%
% I/O: result=qn(x);
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% This is help on the C MEX-file "qn.dll". If full Matlab code is required,
% use the qnm.m function.
%
% Original Fortran Code by C. Croux, translated to MATLAB by S. Verboven.
% Last update : 13/07/2000

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function res = lt(a,b)
%LT           Implements  a < b  for gradients, compares only a.x and b.x
%

% written  10/16/98     S.M. Rump
% modified 12/18/02     S.M. Rump  complex comparison
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'gradient')
    if isreal(a) & isreal(b.x)
      res = ( a<b.x );
    else
      res = ( real(a)<real(b.x) ) & ( imag(a)<imag(b.x) );
    end
  elseif ~isa(b,'gradient')
    if isreal(a.x) & isreal(b)
      res = ( a.x<b );
    else
      res = ( real(a.x)<real(b) ) & ( imag(a.x)<imag(b) );
    end
  else
    if isreal(a.x) & isreal(b.x)
      res = ( a.x<b.x );
    else
      res = ( real(a.x)<real(b.x) ) & ( imag(a.x)<imag(b.x) );
    end
  end

repo: OpenMAS


text:
function [ bool ] = gsp_test_duality_coefficient( gcoeff,hcoeff,tol )
%GSP_TEST_DUALITY_COEFFICIENT Test if the coefficient are from dual filters
%   Usage:  bool = gsp_test_duality_coefficient( gcoeff,hcoeff );
%           bool = gsp_test_duality_coefficient( gcoeff,hcoeff,tol );
%
%   Input parameters:
%       gcoeff  : coefficient of the filter 1 (matrix $N$ x $M$ )
%       hcoeff  : coefficinet of the filter 2 (matrix $N$ x $M$ )
%       tol     : tolerance for the test (default 1e-5)
%
%   Ouput paramters:
%       bool    : boolean
%
%   This function test if two discrete filterbanks are dual. Each filter is
%   a column in the matrix *gcoeff* or *hcoeff*. $M$ is the number of
%   filters and $N$ the number of coefficients (size of the graph signal).
%

% Author: Nathanael Perraudin
% Date  : 13 july 2014
% testing: test_dual

if nargin<3
    tol = 1e-5;
end

v = sum(gcoeff.*hcoeff,2);
A = min(v);
B = max(v);
bool = ( (A-B) / ((A+B)/2) ) < tol;

end

repo: gspbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This scripts checks whether the required data is already available.
% If if so and possible, it also initializes the modules
% viewimage, imgModel, trafo, distance, regularizer for later use;
% see setup2DhandsData.m for an example
%==============================================================================
caller          = dbstack;        % identify the name of the calling function
caller          = caller(min(length(caller),2)).name;
expfile         = fullfile(FAIRpath,'temp',[caller,'.mat']);
expfileExists   = (exist(expfile,'file') == 2);

if ~expfileExists,
  return;         % matfile has to be generated
end;

fprintf('%s: load(%s); and initialize modules\n',...
  caller,expfile(length(FAIRpath)+1:end));
clear caller
load(expfile);

% initialize toolbox
if ~isempty(who(expfile,'viewPara')),    viewImage('reset',viewPara{:});  end;
if ~isempty(who(expfile,'imgPara')),     imgModel('reset',imgPara{:});    end;
if ~isempty(who(expfile,'traPara')),     trafo('reset',traPara{:});       end;
if ~isempty(who(expfile,'disPara')),     distance('reset',disPara{:});    end;
if ~isempty(who(expfile,'regPara')),     regularizer('reset',regPara{:}); end;
%==============================================================================


text:
function Bnew = pollBADS2N(B,x,gpstruct,LB,UB,optimState,options)
%POLLBADS2N Poll 2N fixed basis vectors (Generalized Pattern Search).

nvars = numel(x);

if isempty(B)
    Bnew = [eye(nvars); -eye(nvars)];
    Bnew = Bnew*optimState.C';

    if 0
        if size(gpstruct.x,1) > (2*nvars + nvars*nvars/4)
            try
                % H = hessian(@(xi_) gppred(xi_,gpstruct), x);
                H = fhess(@(xi_) gppred(xi_,gpstruct), x, [], 'step', optimState.searchmeshsize);
                H = nearestSPD(H,10);
                % eig(H)
                L = chol(H);
                % Bnew = solve_chol(L,Bnew')';
                Hinv = solve_chol(L,eye(nvars));
                Hinv = nearestSPD(Hinv,10);
                [V,D] = eig(Hinv);
                lD = real(log(diag(D)));
                lD(lD < max(lD) + 0.5*log(eps)) = max(lD) + 0.5*log(eps);
                while 1
                    delta = max(lD) - min(lD);
                    if delta < log(1e6); break; end
                    lD = 0.8*lD;
                end
                D = exp(0.5*lD);
                Bnew = bsxfun(@times, [V'; -V'], [D(:);D(:)]);
                % Hinv = inv(H);
                % Bnew = (chol(Hinv)*Bnew')';
            catch
                % Use diagonal matrix
            end
        end
    end

    % Global vector normalization
    D = sqrt(sum(Bnew(1:nvars,:).*Bnew(1:nvars,:),2));
    N = exp(log(D) - mean(log(D)))./D;
    Bnew = bsxfun(@times, Bnew, [N;N]);
    Bnew = bsxfun(@rdivide, Bnew, gpstruct.pollscale);

else
    Bnew = [];
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
% TERNARY.COORDS calculate rectangular coordinates of fractions on a ternary plot
%   [X, Y] = TERNARY.COORDS(FA, FB) returns the rectangular X and Y coordinates
%   for the point with a fraction defined by FA and FB.  It is assumed that
%   FA and FB are sensible fractions.
%
%   [X, Y] = TERNARY.COORDS(FA, FB, FC) returns the same.  FC is assumed to be
%   the remainder when subtracting FA and FB from 1.

% Author: Carl Sandrock 20050211

% Modifications

% Modifiers

function [x, y] = coords(fA, fB, fC)
if nargin < 3
    fC = 1 - (fA + fB);
end

y = fB*sin(deg2rad(60));
x = fA + y*cot(deg2rad(60));%fA+fb*cos(deg2rad(60))
return

author: Sable

repo: mcbench-benchmarks


text:
%% cellEig
% Below is a demonstration of the features of the |cellEig| function

%% Syntax
% |[V,D]=cellEig(C);|

%% Description
% Computes eigenvalues and eigenvectors for each matrix contained in the
% cell array C, i.e. [v,d]=eig(c) is executed for each cell entry. The
% output is two cell arrays, i.e. the cell V containing the eigenvectors
% and the cell D containing the eigenvalues.

%% Examples

%%
clear; close all; clc;

%% Example: Calculating eigenvalues for matrices contained in cells
% Creating example cell containing two matrices

M1=rand(3,3);
M1=M1*M1';
M2=rand(5,5);
M2=M2*M2';

C={M1,M2};
[V,D]=cellEig(C);

%%
% Contained in the output cells are the eigenvectors and eigenvalues of
% each of the matrices e.g. for the first
v1=V{1}
d1=D{1}

%%
% and the second entry
v2=V{2}
d2=D{2}

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function plotx2north
% set the default plot direction of the x-axis

setMTEXpref('xAxisDirection','north');

author: mtex-toolbox

repo: mtex

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


text:
function days = month_length_republican ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_REPUBLICAN returns the number of days in a Republican month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer DAYS, the number of days
%    in the month.
%

%
%  Copy the input.
%
  m2 = m;
  y2 = y;
%
%  Check the input.
%
  [ y2, m2, ierror ] = ym_check_republican ( y2, m2 );

  if ( ierror ~= 0 )
    days = 0;
    return
  end
%
%  Get the number of days in the month.
%
  if ( 1 <= m2 && m2 <= 12 )
    days = 30;
  else if ( m2 == 13 )
    if ( year_is_leap_republican ( y2 ) )
      days = 6;
    else
      days = 5;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
addpath(genpath('./dependence'))
addpath data/

slaveRun('localhost', 10000);

author: huashiyiqike

repo: LSTM-MATLAB


text:
function v = fdm_2d_vector(n0,f_str)
%
%  Generates a vector v which contains the values of a function f(x,y)
%  on an equidistant grid in the interior of the unit square. The grid
%  points are numbered consistently with those used in the function
%  'fdm_2d_matrix'.
%
%  This function is just used as an easy way to generate test problems
%  rather than to solve PDEs.
%
%  Calling sequence:
%
%    v = fdm_2d_vector( n0, f_str)
%
%  Input:
%
%    n0        number of inner grid points in each dimension;
%    f_str     string describing the function f in the space variables 'x'
%              and 'y', e.g., f_str = 'sin(x+2*y)+3'.
%
%  Output:
%
%    v         vector of order n = n0^2 containing the values of f(x,y).
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

% Input data not completely checked!

na = nargin;

if na~=2
  error('Wrong number of input parameters.');
end

h = 1.0/(n0+1);

n2 = n0*n0;

v = zeros(n2,1);

i = 0;

for iy = 1:n0
  y = iy*h;
  for ix = 1:n0
    x = ix*h;
    i = i+1;
    v(i) = eval(f_str);
  end
end


text:
function X = newtonraphson(Eqn_Str,Start_Point,Max_Iter)

% Function :
%        X = newtonraphson(fun_str,Start_Point)
%              Finds the root of an equation by NEWTON-RAPHSON METHOD.
%
%         x(n+1) = x(n) - (F(x(n))/diff(F(x(n)));
%
% INPUTS :
%        Eqn_Str     : The equation whose root has to be find.
%                      Eqn_Str Should be an string format.
%                      Example :
%                       if  F(x) = x^4 - x - 10
%                       the Eqn_Str = 'x^4 - x - 10';
%
%        Start_Point : Initial value of root X ==> x0
%
%        Max_Iter    : Maximum number of iterations
%
% OUTPUT :
%        X           : Estimeted root of Equation.
%--------------------------------------------------------------------------
% By :-
%      SANDEEP SOLANKI
%      rtm_sandeep@rediffmail.com
%--------------------------------------------------------------------------
if nargin < 2
    error('Not Enough Input Arguments');
end

if nargin < 3
    Max_Iter = 1000;
end

if ~ischar(Eqn_Str)
    error('Funtion Should Be an String');
end

fx = inline(Eqn_Str);
fxd = inline(diff(Eqn_Str));

x0 = Start_Point;
disp(['F(X) = ' Eqn_Str]);
disp(['X0 = ' num2str(x0)]);
% Iterating
for i = 1 : Max_Iter
    s1=sprintf('    Iteration : %1.0f',i);
    disp(s1);
    x1 = x0 - (fx(x0)/fxd(x0));

    s2=sprintf('              X(%0.0f) = %0.15f',i,x1);
    disp(s2);
    if x1 == x0
        disp('Terminating Process : Value of Root Repeated');
        break;
    else
        x0 = x1;
    end
end

if x1 == NaN
    error('Start Point is Not Correct');
end
X = x1;


text:
function data= bbci_load(bbci)
%BBCI_LOAD - Load data set
%
%Synopsis:
%  DATA= bbci_load(BBCI);

bbci= bbci_calibrate_setDefaults(bbci);
BC= bbci.calibrate;

% fullfile does not work for CELL, so we use strcat here:
if ~isempty(BC.folder) && ~fileutil_isAbsolutePath(BC.file),
    data.filename= strcat(BC.folder, filesep, BC.file);
else
    data.filename = BC.file;
end
filespec= strcat(data.filename,'.*');
if iscell(filespec),
    data.fileinfo= cellfun(@dir, filespec, 'UniformOutput',false);
else
    data.fileinfo= {dir(filespec)};
end
if isempty(BC.montage_fcn),
    [data.cnt, data.mrk, data.mnt]= ...
        BC.read_fcn(data.filename, BC.read_param{:});
else
    [data.cnt, data.mrk]= BC.read_fcn(data.filename, BC.read_param{:});
    data.mnt= BC.montage_fcn(data.cnt.clab, BC.montage_param{:});
end
if ~isempty(BC.marker_fcn),
    data.mrk= BC.marker_fcn(data.mrk, BC.marker_param{:});
end
data.isnew= 1;

author: bbci


text:
function program_03 ( )

%*****************************************************************************80
%
%% PROGRAM_03 demonstrates barycentric coordinates.
%
%  Discussion:
%
%    The program
%    * reads a triangle T (defined by three points),
%    * computes and prints the area;
%    * computes and prints the angles;
%    * reads a point P
%    * computes the signed areas of the subtriangles;
%    * computes the barycentric coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROGRAM_03 - Barycentric Coordinates\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define a triangle T:\n' );

  t_v1 = input ( '  Enter [ T.v1.x, T.v1.y]: ' );
  t_v2 = input ( '  Enter [ T.v2.x, T.v2.y]: ' );
  t_v3 = input ( '  Enter [ T.v3.x, T.v3.y]: ' );
%
%  Compute the area.
%
  t_area = 0.5 * ( t_v1(1) * ( t_v2(2) - t_v3(2) ) ...
                 + t_v2(1) * ( t_v3(2) - t_v1(2) ) ...
                 + t_v3(1) * ( t_v1(2) - t_v2(2) ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Triangle area = %f\n', t_area );
%
%  Check various points.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now we compute barycentric coordinates for points.\n' );
  fprintf ( 1, '\n' );

  while ( 1 )

    p = [];
    fprintf ( 1, '\n' );
    p = input ( '  Enter a point [ P.x, P.y], or RETURN to quit: ' );

    if ( isempty ( p ) )
      break
    end
%
%  Compute the areas of the three subtriangles.
%
    tpbc_area = 0.5 * ( p(1)    * ( t_v2(2) - t_v3(2) ) ...
                      + t_v2(1) * ( t_v3(2) - p(2)    ) ...
                      + t_v3(1) * ( p(2)    - t_v2(2) ) );

    tapc_area = 0.5 * ( t_v1(1) * ( p(2)    - t_v3(2) ) ...
                      + p(1)    * ( t_v3(2) - t_v1(2) ) ...
                      + t_v3(1) * ( t_v1(2) - p(2)    ) );

    tabp_area = 0.5 * ( t_v1(1) * ( t_v2(2) - p(2)    ) ...
                      + t_v2(1) * ( p(2)    - t_v1(2) ) ...
                      + p(1)    * ( t_v1(2) - t_v2(2) ) );

    xi1 = tpbc_area / t_area;
    xi2 = tapc_area / t_area;
    xi3 = tabp_area / t_area;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  PBC = %8.4f  XI1 = %8.4f\n', tpbc_area, xi1 );
    fprintf ( 1, '  APC = %8.4f  XI2 = %8.4f\n', tapc_area, xi2 );
    fprintf ( 1, '  APB = %8.4f  XI3 = %8.4f\n', tabp_area, xi3 );
    sum1 = tpbc_area + tapc_area + tabp_area;
    sum2 = xi1 + xi2 + xi3;
    fprintf ( 1, '  Sum = %8.4f  Sum = %8.4f\n', sum1, sum2 );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROGRAM_03\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

diff = R .* (X * Theta' - Y);
sum_squared = @(A) sum(sum(A .^ 2));

J = (sum_squared(diff) + lambda * (sum_squared(X)  + sum_squared(Theta))) /2;
X_grad = diff * Theta + lambda * X;
Theta_grad = diff' * X + lambda * Theta;

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

grad = [X_grad(:); Theta_grad(:)];

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
function z = triumtriu(x,y,K)
% z = triumtriu(x,y,K)
%
% TRIUMTRIU  Computes y = x * y
%   Both r and u should be upper triangular.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

Ks = K.s;
if isempty(Ks)
    z = [];
    return
end
Kq = Ks .* Ks;
nr = K.rsdpN;
nc = length(Ks);
N  = K.N - K.sblkstart(1) + 1;
z  = zeros(N,1);
xi = length(x) - N;
zi = 0;
for i = 1 : nc
    ki = Ks(i);
    qi = Kq(i);
    XX = x(xi+1:xi+qi);
    YY = y(xi+1:xi+qi);
    xi = xi + qi;
    if i > nr
        XX = XX + 1j * x(xi+1:xi+qi);
        YY = YY + 1j * y(xi+1:xi+qi);
        xi = xi + qi;
    end
    ZZ = triu(reshape(XX,ki,ki)) * triu(reshape(YY,ki,ki));
    ZZ = ZZ + triu(ZZ,1)';
    z(zi+1:zi+qi) = real(ZZ);
    zi = zi + qi;
    if i > nr
        z(zi+1:zi+qi) = imag(ZZ);
        zi = zi + qi;
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function filename=cosmo_make_temp_filename(prefix,suffix)
% give temporary filename that does not exist when this function is called
%
% filename=cosmo_make_temp_filename(prefix,suffix)
%
% Inputs:
%   prefix              optional, string with prefix for temporary file
%                       (default: '')
%   suffix              string [or cellstring] with suffix [suffixes] for
%                       temporary file [files]
%                       (default: '')
%
% Output:
%   filename            filename that does not exist when this function was
%                       called (in the current directory), starting with
%                       prefix, ending with suffix, and with a random infix
%                       string.
%                       If suffix is a cellstring, then filename is a
%                       cellstring with as many elements as suffix,
%                       each with the same prefix and random infix.
%
% Examples:
%   % generate random filename
%   fn=cosmo_make_temp_filename();
%
%   % generate random filename starting with 'foo'
%   fn=cosmo_make_temp_filename('foo');
%
%   % generate random filename starting with 'foo' and ending with '.bar'
%   fn=cosmo_make_temp_filename('foo','.bar');
%
%   % generate two random filename, both starting with 'foo', and ending
%   % with '.bar' and '.baz', respectively. The output fns is a cell with
%   % two strings
%   fns=cosmo_make_temp_filename('foo',{'.bar','.bar'});
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if nargin<2, suffix=''; end
    if nargin<1, prefix=''; end

    check_input(prefix,suffix);

    suffix_is_cell=iscell(suffix);

    if suffix_is_cell
        suffix_cell=suffix;
    else
        suffix_cell={suffix};
    end

    nsuffix=numel(suffix_cell);

    filename_cell=cell(nsuffix,1);
    while true
        infix=generate_random_infix();

        does_exist=false;

        for k=1:nsuffix
            fn=[prefix infix suffix_cell{k}];

            if exist(fn,'file')
                does_exist=true;
                break;
            end

            filename_cell{k}=fn;
        end

        if ~does_exist
            break;
        end
    end

    if suffix_is_cell
        filename=filename_cell;
    else
        filename=filename_cell{1};
    end

function infix=generate_random_infix(infix_length)
    if nargin<1
        infix_length=20;
    end

    char_val_min=double('a');
    char_val_max=double('z');
    char_val_range=char_val_max-char_val_min;

    rand_chars=char(rand(1,infix_length)*char_val_range+char_val_min);
    infix=['tmp_' rand_chars];

function check_input(prefix,suffix)
    if ~ischar(prefix)
        error('prefix must be a string');
    end

    if ~ischar(suffix) && ~iscellstr(suffix)
        error('suffix must be string or cell with strings');
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%% Quadrature
% Using quad to integrate f(x)

clear
clc

%% Evaluating quad
n=0;
for i=[2,5.5,9]
    n=n+1;
    F=@(x)((1./x)-1).^(i-1).*exp(1-(1./x))./(x.^2);
    [q(n,1) p(n,1)]=quad(F,0,1,1.0e-8);
end
[q p]

repo: Aero-matlab


text:
%%*******************************************************************
%%  Read in a problem in SeDuMi format and convert to DSDP format
%%
%%  [AC,b] = readsedumifile(fname)
%%
%%  Input: fname.mat = name of the file containing SDP data in
%%                     SeDuMi format.
%%
%%  Input: AC  = DSDP formatted block structure and constraint data
%%         b   = Dual objective vector
%%
%% DSDP: version 5.0
%% Copyright (c) 2004 by
%% S. Benson Y. Ye
%% Last modified: 12 Jan 04
%%******************************************************************

  function [AC,b]=readsedumifile(fname);
%%
%%  First, load the matlab file containing At, c, b, and K
%%
  K.q = [];
  K.l = [];
  A = 0;
  At = 0;
  compressed = 0;
  if exist([fname,'.mat.gz']);
     compressed = 1;
     unix(['gunzip ', fname,'.mat.gz']);
  elseif exist([fname,'.gz']);
     compressed = 2;
     unix(['gunzip ', fname,'.gz']);
  elseif exist([fname,'.mat.Z']);
     compressed = 3;
     unix(['uncompress ', fname,'.mat.Z']);
  elseif exist([fname,'.Z']);
     compressed = 4;
     unix(['uncompress ', fname,'.Z']);
  end
  if exist([fname,'.mat']) | exist(fname)
     eval(['load ', fname]);
  else
     error('** Problem not found, please specify the correct path.');
  end

  if (compressed == 1)
     unix(['gzip ', fname,'.mat']);
  elseif (compressed == 2)
     unix(['gzip ', fname]);
  elseif (compressed == 3)
     unix(['compress ', fname,'.mat']);
  elseif (compressed == 4)
     unix(['compress ', fname]);
  end
%%
  if (size(c,1) == 1), c = c'; end;
  if (size(b,1) == 1), b = b'; end;
  if (At == 0), At = A'; end;  clear A;
  [nn,mm] = size(At); if (max(size(c)) == 1); c = c*ones(nn,1); end;

  if ~isfield(K,'l'); K.l = 0; end
  if ~isfield(K,'q'); K.q = 0; end
  if ~isfield(K,'s'); K.s = 0; end
  if K.l == 0 | isempty(K.l); K.l = 0; end;
  if sum(K.q) == 0 | isempty(K.q); K.q = 0; end
  if sum(K.s) == 0 | isempty(K.s); K.s = 0; end
%%
%%
%%
%%  Then convert formats

  [AC b]=readsedumi(At,b,c,K);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function BrainNet=tHOFC(BOLD)
% Topographical similarity-based high-order FC network construction
%
% Input:
% BOLD. A cell array with size of N x 1, each cell is a matrix of BOLD signals (#time points x #ROIs) from one of N subject
%       Each subject may have different #time points but should have the same #ROIs
%
% Output:
% BrainNet. High-order brain networks (#ROIs x #ROIs x #Subjects)
%
% By Han Zhang, hanzhang@med.unc.edu
% IDEA lab https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, UNC Chapel Hill
%
% Cite: 1. Zhang, H., Chen, X., Shi, F., Li, G., Kim, M., Giannakopoulos, P., Haller, S., Shen, D., Jun 28, 2016.
%          Topographic Information based High-Order Functional Connectivity and its Application in Abnormality Detection
%          for Mild Cognitive Impairment, Journal of Alzheimer's Disease. In press. DOI: 10.3233/JAD-160092.
%       2. Zhang, H., Chen, X., Zhang, Y., Shen, D., Test-retest reliability of �high-order� functional connectivity
%          in young healthy adults. Frontiers in Neuroscience, In Press. DOI: 10.3389/fnins.2017.00439.

[nTime,nROI]=size(BOLD{1});
nSubj=length(BOLD);
BrainNet=zeros(nROI,nROI,nSubj,'single');

% Compute low-order FC network
pc=PC(BOLD);   % Pearson's correlation-based network
pc=atanh(pc);       % Fisher z transformation

% Compute topographical similarity which defines high-order FC
for ns=1:nSubj
    tempNet=zeros(nROI,nROI);
    for i=1:nROI
        for j=(i+1):nROI
            temp1=pc(:,i,ns);
            temp2=pc(:,j,ns);
            temp1([i,j])=[];                    % ignore self connections
            temp2([i,j])=[];                    % ignore self connections
            tempNet(i,j)=corr(temp1,temp2);     % only upper triangle
        end
    end
    BrainNet(:,:,ns)=tempNet+tempNet';          % form symmetric adjancency matrix (main diag are zeros)
end

author: zzstefan

repo: BrainNetClass


text:
function [bval] = function(fval,bounds,bits)
% function [bval] = f2b(fval,bounds,bits)
%
% Return the binary representation of the float number fval.
%
% fval   - the float representation of the number
% bval   - the binary representation of the number
% bounds - the bounds on the variables
% bits   - the number of bits to represent each variable

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

scale=(2.^bits-1)./ (bounds(:,2)-bounds(:,1))'; %The range of the variables
numV=size(bounds,1);
cs=[0 cumsum(bits)];
bval=[];
for i=1:numV
  fval(i)=(fval(i)-bounds(i,1)) * scale(i);
  bval=[bval rem(floor(fval(i)*pow2(1-bits(i):0)),2)];
end

author: Grootzz

repo: GA-BP

save_path: github-repos/MATLAB/Grootzz-GA-BP


text:
function [tfr,t,f]=tfrideal(iflaws,t,N,trace)
%TFRIDEAL Ideal TFR for given instantaneous frequency laws.
%	[TFR,T,F]=TFRIDEAL(IFLAWS,T,N,TRACE) generates the ideal
%	time-frequency representation corresponding to the
%	instantaneous frequency laws of the components of a signal.
%
%	IFLAWS : (M,P)-matrix where each column corresponds to
%		 the instantaneous frequency law of an (M,1)-signal,
%		 These P signals do not need to be present at the same time.
%		 The values of IFLAWS must be between 0 and 0.5.
%	T      : the time instant(s)      (default : 1:M).
%	N      : number of frequency bins (default : M).
%	TRACE  : if nonzero, the progression of the algorithm is shown
%                                         (default : 0).
%	TFR    : output time-frequency matrix, of size (N,length(t)).
%		 If nargout=0, a contour plot of TFR is automatically
%		 displayed on the screen.
%	F      : vector of normalized frequencies.
%
%	Example :
%         N=140; t=0:N-1; [x1,if1]=fmlin(N,0.05,0.3);
%         [x2,if2]=fmsin(70,0.35,0.45,60);
%         if2=[zeros(35,1)*NaN;if2;zeros(35,1)*NaN];
%         tfrideal([if1 if2]);
%
%	See also PLOTIFL, PLOTSID, and all the time-frequency
%        representations listed in the CONTENTS file (TFR*)

%	O. Lemoine, F. Auger - March, April 1996.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin==0),
 error('at least one parameter required');
end;

[ifrow,ifcol]=size(iflaws);

if (nargin==1),
 t=1:ifrow; N=ifrow; trace=0;
elseif (nargin==2),
 N=ifrow; trace=0;
elseif (nargin==3),
 trace=0;
end;

[trow,tcol]=size(t);
if (trow~=1),
 error('T must only have one row');
end;

tfr=zeros(N,tcol);

if any(any(iflaws>0.5)) | any(any(iflaws<0)),
 error('The values of IFLAWS must be between 0 and 0.5');
end

if trace, disp('Ideal time-frequency distribution'); end;

for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 ti= t(icol);
 for fi=1:ifcol,
  if isnan(iflaws(ti,fi)),
   tfr(fi,icol)=NaN;
  else
   tfr(round(iflaws(ti,fi)*2*(N-1))+1,icol)=1;
  end
 end
end;

if (nargout==0),
 f=(0:N-1)/(2*N);
 axes(gca); contour(t,f,tfr,1,'k');
 xlabel('Time'); ylabel('Normalized frequency');
 title('Ideal time-frequency representation');
 grid;
elseif (nargout==3),
 f=(0.5*(0:N-1)/N)';
end;

author: HeLiangHIT


text:
function [label, model, L] = vbgm(X, init, prior)
% Perform variational Bayesian inference for Gaussian mixture.
%   X: d x n data matrix
%   init: k (1 x 1) or label (1 x n, 1<=label(i)<=k) or center (d x k)
% Reference: Pattern Recognition and Machine Learning by Christopher M. Bishop (P.474)
% Written by Michael Chen (sth4nth@gmail.com).

fprintf('Variational Bayesian Gaussian mixture: running ... \n');
[d,n] = size(X);
if nargin < 3
    prior.alpha = 1;
    prior.kappa = 1;
    prior.m = mean(X,2);
    prior.v = d+1;
    prior.M = eye(d);   % M = inv(W)
end
tol = 1e-20;
maxiter = 2000;
L = -inf(1,maxiter);
converged = false;
t = 1;

model.R = initialization(X,init);
while  ~converged && t < maxiter
    t = t+1;
    model = vmax(X, model, prior);
    model = vexp(X, model);
    L(t) = vbound(X,model,prior)/n;
    converged = abs(L(t)-L(t-1)) < tol*abs(L(t));
end
L = L(2:t);
label = zeros(1,n);
[~,label(:)] = max(model.R,[],2);
[~,~,label] = unique(label);
if converged
    fprintf('Converged in %d steps.\n',t-1);
else
    fprintf('Not converged in %d steps.\n',maxiter);
end

function R = initialization(X, init)
[d,n] = size(X);
if length(init) == 1  % random initialization
    k = init;
    idx = randsample(n,k);
    m = X(:,idx);
    [~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
    [u,~,label] = unique(label);
    while k ~= length(u)
        idx = randsample(n,k);
        m = X(:,idx);
        [~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
        [u,~,label] = unique(label);
    end
    R = full(sparse(1:n,label,1,n,k,n));
elseif size(init,1) == 1 && size(init,2) == n  % initialize with labels
    label = init;
    k = max(label);
    R = full(sparse(1:n,label,1,n,k,n));
elseif size(init,1) == d  %initialize with only centers
    k = size(init,2);
    m = init;
    [~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2),[],1);
    R = full(sparse(1:n,label,1,n,k,n));
else
    error('ERROR: init is not valid.');
end
% Done
function model = vmax(X, model, prior)
alpha0 = prior.alpha;
kappa0 = prior.kappa;
m0 = prior.m;
v0 = prior.v;
M0 = prior.M;
R = model.R;

nk = sum(R,1); % 10.51
alpha = alpha0+nk; % 10.58
nxbar = X*R;
kappa = kappa0+nk; % 10.60
m = bsxfun(@times,bsxfun(@plus,kappa0*m0,nxbar),1./kappa); % 10.61
v = v0+nk; % 10.63

[d,k] = size(m);
M = zeros(d,d,k);
sqrtR = sqrt(R);

xbar = bsxfun(@times,nxbar,1./nk); % 10.52
xbarm0 = bsxfun(@minus,xbar,m0);
w = (kappa0*nk./(kappa0+nk));
for i = 1:k
    Xs = bsxfun(@times,bsxfun(@minus,X,xbar(:,i)),sqrtR(:,i)');
    xbarm0i = xbarm0(:,i);
    M(:,:,i) = M0+Xs*Xs'+w(i)*(xbarm0i*xbarm0i'); % 10.62
end

model.alpha = alpha;
model.kappa = kappa;
model.m = m;
model.v = v;
model.M = M; % Whishart: M = inv(W)
% Done
function model = vexp(X, model)
alpha = model.alpha; % Dirichlet
kappa = model.kappa;   % Gaussian
m = model.m;         % Gasusian
v = model.v;         % Whishart
M = model.M;         % Whishart: inv(W) = V'*V

n = size(X,2);
[d,k] = size(m);

logW = zeros(1,k);
EQ = zeros(n,k);
for i = 1:k
    U = chol(M(:,:,i));
    logW(i) = -2*sum(log(diag(U)));
    Q = (U'\bsxfun(@minus,X,m(:,i)));
    EQ(:,i) = d/kappa(i)+v(i)*dot(Q,Q,1);    % 10.64
end

ElogLambda = sum(psi(0,bsxfun(@minus,v+1,(1:d)')/2),1)+d*log(2)+logW; % 10.65
Elogpi = psi(0,alpha)-psi(0,sum(alpha)); % 10.66

logRho = (bsxfun(@minus,EQ,2*Elogpi+ElogLambda-d*log(2*pi)))/(-2); % 10.46
logR = bsxfun(@minus,logRho,logsumexp(logRho,2)); % 10.49
R = exp(logR);

model.logR = logR;
model.R = R;
% Done
function L = vbound(X, model, prior)
alpha0 = prior.alpha;
kappa0 = prior.kappa;
m0 = prior.m;
v0 = prior.v;
M0 = prior.M;

alpha = model.alpha; % Dirichlet
kappa = model.kappa;   % Gaussian
m = model.m;         % Gasusian
v = model.v;         % Whishart
M = model.M;         % Whishart: inv(W) = V'*V
R = model.R;
logR = model.logR;

[d,k] = size(m);
nk = sum(R,1); % 10.51

Elogpi = psi(0,alpha)-psi(0,sum(alpha));

Epz = dot(nk,Elogpi);
Eqz = dot(R(:),logR(:));
logCalpha0 = gammaln(k*alpha0)-k*gammaln(alpha0);
Eppi = logCalpha0+(alpha0-1)*sum(Elogpi);
logCalpha = gammaln(sum(alpha))-sum(gammaln(alpha));
Eqpi = dot(alpha-1,Elogpi)+logCalpha;
L = Epz-Eqz+Eppi-Eqpi;

U0 = chol(M0);
sqrtR = sqrt(R);
xbar = bsxfun(@times,X*R,1./nk); % 10.52

logW = zeros(1,k);
trSW = zeros(1,k);
trM0W = zeros(1,k);
xbarmWxbarm = zeros(1,k);
mm0Wmm0 = zeros(1,k);
for i = 1:k
    U = chol(M(:,:,i));
    logW(i) = -2*sum(log(diag(U)));

    Xs = bsxfun(@times,bsxfun(@minus,X,xbar(:,i)),sqrtR(:,i)');
    V = chol(Xs*Xs'/nk(i));
    Q = V/U;
    trSW(i) = dot(Q(:),Q(:));  % equivalent to tr(SW)=trace(S/M)
    Q = U0/U;
    trM0W(i) = dot(Q(:),Q(:));

    q = U'\(xbar(:,i)-m(:,i));
    xbarmWxbarm(i) = dot(q,q);
    q = U'\(m(:,i)-m0);
    mm0Wmm0(i) = dot(q,q);
end

ElogLambda = sum(psi(0,bsxfun(@minus,v+1,(1:d)')/2),1)+d*log(2)+logW; % 10.65
Epmu = sum(d*log(kappa0/(2*pi))+ElogLambda-d*kappa0./kappa-kappa0*(v.*mm0Wmm0))/2;
logB0 = v0*sum(log(diag(U0)))-0.5*v0*d*log(2)-logmvgamma(0.5*v0,d);
EpLambda = k*logB0+0.5*(v0-d-1)*sum(ElogLambda)-0.5*dot(v,trM0W);

Eqmu = 0.5*sum(ElogLambda+d*log(kappa/(2*pi)))-0.5*d*k;
logB =  -v.*(logW+d*log(2))/2-logmvgamma(0.5*v,d);
EqLambda = 0.5*sum((v-d-1).*ElogLambda-v*d)+sum(logB);

EpX = 0.5*dot(nk,ElogLambda-d./kappa-v.*trSW-v.*xbarmWxbarm-d*log(2*pi));

L = L+Epmu-Eqmu+EpLambda-EqLambda+EpX;


text:
function [lfp,indices] = GetLFP(channels,varargin)

%GetLFP - Get local field potentials.
%
%  Load local field potentials from disk (unlike spikes or positions, LFP data
%  is usually too large to keep in memory).
%
%  USAGE
%
%    [lfp,indices] = GetLFP(channels,<options>)
%
%    channels       list of channels to load (use keyword 'all' for all)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'restrict'    list of time intervals to read from the LFP file
%     'intervals'   same as 'restrict' (for backwards compatibility)
%     'select'      select channel by ID ('id', counted from 0 a la NeuroScope)
%                   or by number ('number', counted from 1 a la Matlab)
%                   (default = 'id')
%    =========================================================================
%
%  OUTPUT
%
%    lfp            list of (time,voltage1,...,voltageN) tuples
%    indices        for each tuple, the index of the interval it falls in
%
%  EXAMPLES
%
%    % channel ID 5 (= # 6), from 0 to 120 seconds
%    lfp = GetLFP(5,'restrict',[0 120]);
%    % same, plus from 240.2 to 265.23 seconds
%    lfp = GetLFP(5,'restrict',[0 120;240.2 265.23]);
%    % multiple channels
%    lfp = GetLFP([1 2 3 4 10 17],'restrict',[0 120]);
%    % channel # 3 (= ID 2), from 0 to 120 seconds
%    lfp = GetLFP(3,'restrict',[0 120],'select','number');

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

error('This function is now deprecated from /buzcode, try using bz_GetLFP.m')

global DATA;
if isempty(DATA),
	error('No session defined (did you forget to call SetCurrentSession? Type ''help Data'' for details).');
end

% Default values
intervals = [0 Inf];
select = 'id';

if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help GetLFP'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+1) ' is not a property (type ''help GetLFP'' for details).']);
  end
  switch(lower(varargin{i})),
    case {'intervals','restrict'},
      intervals = varargin{i+1};
      if ~isdmatrix(intervals) || size(intervals,2) ~= 2,
        error('Incorrect value for property ''intervals'' (type ''help GetLFP'' for details).');
      end
    case 'select',
      select = lower(varargin{i+1});
      if ~isstring_FMAT(select,'id','number'),
        error('Incorrect value for property ''select'' (type ''help GetLFP'' for details).');
      end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help GetLFP'' for details).']);
  end
end

filename = [DATA.session.path '/' DATA.session.basename '.lfp'];
if ~exist(filename,'file'),
    filename = [DATA.session.path '/' DATA.session.basename '.lfp'];
    if ~exist(filename,'file'),
        error(['File ''' filename ''' not found.']);
    end
end
nChannels = DATA.nChannels;
if isa(channels,'char') && strcmp(lower(channels),'all'),
	channels = (1:nChannels)-1;
end

if strcmp(select,'id'),
	channels = channels + 1;
end

nIntervals = size(intervals,1);
lfp = [];
indices = [];
for i = 1:nIntervals,
	duration = (intervals(i,2)-intervals(i,1));
	start = intervals(i,1);
	% Load data
	data = LoadBinary(filename,'duration',duration,'frequency',DATA.rates.lfp,'nchannels',nChannels,'start',start,'channels',channels);
    if isinteger(data)
        error('GetLFP does not work with integer data - use bz_GetLFP instead');
    end
	t = start:(1/DATA.rates.lfp):(start+(length(data)-1)/DATA.rates.lfp);t=t';
	lfp = [lfp ; t data];
	indices = [indices ; i*ones(size(t))];
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function v = haar_2d ( u )

%*****************************************************************************80
%
%% HAAR_2D computes the Haar transform of an array.
%
%  Discussion:
%
%    Thanks to Stephen Becker for pointing out that a previous version of
%    the haar_2d code was not inverted by haar_2d_inverse in cases where
%    M and N were not powers of 2, 05 March 2014.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real U(M,N), the vector to be transformed.
%
%    Output, real V(M,N), the transformed vector.
%
  [ m, n ] = size ( u );

  v = u;

  s = sqrt ( 2.0 );

  w = zeros ( m, n );
%
%  Determine K, the largest power of 2 such that K <= M.
%
  k = 1;
  while ( k * 2 <= m )
    k = k * 2;
  end
%
%  Transform all columns.
%
  while ( 1 < k )

    k = floor ( k / 2 );

    w(  1:  k,:) = ( v(1:2:2*k-1,:) + v(2:2:2*k,:) ) / s;
    w(k+1:k+k,:) = ( v(1:2:2*k-1,:) - v(2:2:2*k,:) ) / s;

    v(1:2*k,:) = w(1:2*k,:);

  end
%
%  Determine K, the largest power of 2 such that K <= N.
%
  k = 1;
  while ( k * 2 <= n )
    k = k * 2;
  end
%
%  Transform all rows.
%
  while ( 1 < k )

    k = floor ( k / 2 );

    w(:,  1:  k) = ( v(:,1:2:2*k-1) + v(:,2:2:2*k) ) / s;
    w(:,k+1:k+k) = ( v(:,1:2:2*k-1) - v(:,2:2:2*k) ) / s;

    v(:,1:2*k) = w(:,1:2*k);

  end

  return
end


text:
function Q = spm_MDP_VB_game(MDP)
% auxiliary plotting routine for spm_MDP_VB - multiple trials
% FORMAT Q = spm_MDP_VB_game(MDP)
%
% MDP.P(M,T)      - probability of emitting action 1,...,M at time 1,...,T
% MDP.X           - conditional expectations over hidden states
% MDP.R           - conditional expectations over policies
% MDP.O(O,T)      - a sparse matrix encoding outcomes at time 1,...,T
% MDP.S(N,T)      - a sparse matrix encoding states at time 1,...,T
% MDP.U(M,T)      - a sparse matrix encoding action at time 1,...,T
% MDP.W(1,T)      - posterior expectations of precision
%
% MDP.xn  = Xn    - simulated neuronal encoding of policies
% MDP.wn  = wn    - simulated neuronal encoding of precision
% MDP.da  = dn    - simulated dopamine responses (deconvolved)
% MDP.rt  = rt    - simulated dopamine responses (deconvolved)
%
% returns summary of performance:
%
%     Q.X  = x    - expected hidden states
%     Q.R  = u    - final policy expectations
%     Q.S  = s    - initial hidden states
%     Q.O  = o    - final outcomes
%     Q.p  = p    - performance
%     Q.q  = q    - reaction times
%
% please see spm_MDP_VB
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_VB_game.m 7766 2020-01-05 21:37:39Z karl $

% numbers of transitions, policies and states
%--------------------------------------------------------------------------
if iscell(MDP(1).X)
    Nf = numel(MDP(1).B);                 % number of hidden state factors
    Ng = numel(MDP(1).A);                 % number of outcome factors
else
    Nf = 1;
    Ng = 1;
end

% graphics
%==========================================================================
Nt    = numel(MDP);                      % number of trials
Ne    = MDP(1).T;                        % number of outcomes per trial
Np    = size(MDP(1).V,2);                % number of policies
for i = 1:Nt

    % assemble expectations of hidden states and outcomes
    %----------------------------------------------------------------------
    for j = 1:Ne
        for k = 1:Ne
            for f = 1:Nf
                try
                    x{f}{i,1}{k,j} = gradient(MDP(i).xn{f}(:,:,j,k)')';
                catch
                    x{f}{i,1}{k,j} = gradient(MDP(i).xn(:,:,j,k)')';
                end
            end
        end
    end

    % fill-in selected policies if missing
    %----------------------------------------------------------------------
    try
        MDP(i).v;
    catch
        MDP(i).v = [];
    end

    s(:,i) = MDP(i).s(:,1);
    o(:,i) = MDP(i).o(:,end);
    u{1,i} = MDP(i).R;
    v{1,i} = MDP(i).v;
    w(:,i) = MDP(i).dn,2;

    % assemble context learning
    %----------------------------------------------------------------------
    for f = 1:Nf
        try
            try
                D = MDP(i).d{f};
            catch
                D = MDP(i).D{f};
            end
        catch
            try
                D = MDP(i).d;
            catch
                D = MDP(i).D;
            end
        end
        d{f}(:,i) = D/sum(D);
    end

    % assemble performance
    %----------------------------------------------------------------------
    p(i)  = 0;
    for g = 1:Ng
        try
            U = spm_softmax(MDP(i).C{g});
        catch
            U = spm_softmax(MDP(i).C);
        end
        for t = 1:Ne
            p(i) = p(i) + log(U(MDP(i).o(g,t),t))/Ne;
        end
    end

    % reaction time
    %----------------------------------------------------------------------
    try
        q(i) = sum(MDP(i).rt(2:end));
    catch
        q(i) = 0;
    end

end

% assemble output structure if required
%--------------------------------------------------------------------------
if nargout
    Q.X  = x;            % expected hidden states
    Q.R  = u;            % final policy expectations
    Q.S  = s;            % inital hidden states
    Q.O  = o;            % final outcomes
    Q.p  = p;            % performance
    Q.q  = q;            % reaction times
    return
end

% Initial states and expected policies (habit in red)
%--------------------------------------------------------------------------
col   = {'r.','g.','b.','c.','m.','k.'};
u     = spm_cat(u);
v     = spm_cat(v);
t     = 1:Nt;
subplot(6,1,1)
if Nt < 64
    MarkerSize = 24;
else
    MarkerSize = 16;
end
image(t,1:Np,64*(1 - u)),  hold on
plot(linspace(1,Nt,numel(v)),v,'c.','MarkerSize',16)
for f = 1:Nf
    for i = 1:max(s(f,:))
        j = find(s(f,:) == i);
        plot(t(j),j - j + f - Nf,col{rem(i - 1,6)+ 1},'MarkerSize',MarkerSize)
    end
end
set(gca,'YLim',[-Nf,Np] + 1/2), box off
try
    E = spm_softmax(spm_cat({MDP.e}));
    plot(Np*(1 - E(end,:)),'r:')
end
title('Initial state and policy selection')
xlabel('Trial'),ylabel('Policy'), hold off

% Performance
%--------------------------------------------------------------------------
q     = q - mean(q);
q     = q/std(q);
subplot(6,1,2), bar(p,'k'),   hold on
plot(q,'.c','MarkerSize',16), hold on
plot(q,':c')
for g = 1:Ng
    for i = 1:max(o(g,:))
        j = find(o(g,:) == i);
        plot(t(j),j - j + 3 + g,col{rem(i - 1,6)+ 1},'MarkerSize',MarkerSize)
    end
end
title('Final outcome, performance and reaction times')
ylabel('Expected utility'), spm_axis tight, hold off, box off

% Initial states (context)
%--------------------------------------------------------------------------
subplot(6,1,3)
col   = {'r','b','g','c','m','k','r','b','g','c','m','k'};
for f = 1:Nf
    if Nf > 1
        plot(spm_cat(x{f}),col{f}), hold on
    else
        plot(spm_cat(x{f}))
    end
end
title('State estimation (ERPs)'), ylabel('Response'),
spm_axis tight, hold off, box off

% Precision (dopamine)
%--------------------------------------------------------------------------
subplot(6,1,4)
w   = spm_vec(w);
if Nt > 8
    fill([1 1:length(w) length(w)],[0; w.*(w > 0); 0],'k'), hold on
    fill([1 1:length(w) length(w)],[0; w.*(w < 0); 0],'k'), hold off
else
    bar(w,1.1,'k')
end
title('Precision (dopamine)')
ylabel('Precision','FontSize',12), spm_axis tight, box off
YLim = get(gca,'YLim'); YLim(1) = 0; set(gca,'YLim',YLim);
set(gca,'XTickLabel',{});

% learning - D
%--------------------------------------------------------------------------
for f = 1:Nf
    subplot(6*Nf,1,Nf*4 + f), image(64*(1 - d{f}))
    if f < 2
        title('Context Learning')
    end
    set(gca,'XTick',1:Nt);
    if f < Nf
        set(gca,'XTickLabel',{});
    end
    set(gca,'YTick',1);
    try
        set(gca,'YTickLabel',MDP(1).label.factor{f});
    end
    try
        set(gca,'YTickLabel',MDP(1).Bname{f});
    end

end
if isfield(MDP(1),'c')
    title('Learning (C and D)')
else
    return
end

% Habit learning
%--------------------------------------------------------------------------
k     = round(linspace(1,Nt,6));
for j = 1:length(k)
    h = MDP(k(j)).c;
    h = h*diag(1./sum(h));
    subplot(6,6,30 + j), image(64*(1 - h))
    axis image
end

author: spm


text:
function data = simulate_slideNeg(setup)

%This function simulates the stick sliding to the left

%Set up the integration algorithm
Tspan = setup.Tspan;
IC = [  setup.IC.th;
    setup.IC.x;
    setup.IC.dth;
    setup.IC.dx];
eventFunc = @(t,z)events_slideNeg(t,z,setup.P);
options = odeset(...
    'RelTol',setup.tol,...
    'AbsTol',setup.tol,...
    'Vectorized','on',...
    'MaxStep',setup.odeMaxStep,...
    'Events',eventFunc);
userfun = @(t,z)dynamics_slideNeg(t,z,setup.P);

%Run simulation
sol = feval(setup.solver,userfun,Tspan,IC,options);

%Format for post processing
tspan = [sol.x(1),sol.x(end)];
nTime = ceil(setup.dataFreq*diff(tspan));
t = linspace(tspan(1),tspan(2),nTime);
Z = deval(sol,t);
[~, C, E] = dynamics_slideNeg(t,Z,setup.P);

%Store in a nice format for plotting
data.time = t;
data.state.th = Z(1,:);
data.state.x = Z(2,:);
data.state.dth = Z(3,:);
data.state.dx = Z(4,:);
data.state.y = zeros(1,nTime) + setup.IC.y;
data.state.dy = zeros(1,nTime);
data.contact.h = C(1,:);
data.contact.v = C(2,:);
data.energy.potential = E(1,:);
data.energy.kinetic = E(2,:);
data.P = setup.P;

%Get transitions for finite state machine
data.phase = 'SLIDE_NEG';
if isempty(sol.ie)
    data.exit = 'TIMEOUT';
elseif sol.xe(end) ~= sol.x(end);
    %An event was triggered, but it was non-terminal
    data.exit = 'TIMEOUT';
else
    switch sol.ie(end)
        case 1
            data.exit = 'FALL_NEG';
        case 2
            data.exit = 'FALL_POS';
        case 3
            data.exit = 'STUCK';
        case 4
            data.exit = 'LIFT';
        case 5
            data.exit = 'LIFT';
        otherwise
            error('Invalid exit condition in simulate flight')
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [eparams, errors] = calibrateBdtClassifier(data, eclassifier, lab, ncv)
% [eparams, errors] = calibrateEdgeClassifier(efeatures, adjlist, imsegs,
% eclassifier, ncv)
lab(lab == -1 ) = 0;
nfeat = size(data, 1);
for k = 1:ncv
    if ncv > 1
        testind = [(k-1)*nfeat/ncv+1:k*nfeat/ncv];
        trainind = setdiff([1:nfeat], testind);
    else
        trainind = (1:nfeat);
    end
    edata{k} = data(trainind,:);
    elab{k} = lab(trainind);
    econf{k} = test_boosted_dt_mc(eclassifier, edata{k});
%    econf{k} = 1 ./ (1+exp(-econf{k}));
end

for k = 1:ncv
    disp(['iter: ' num2str(k)])
    if ncv>1
        traink = setdiff([1:ncv], k);
    else
        traink = k;
    end
    eparams{k} = fminunc(@(x) objective(x, cat(1, econf{traink}), cat(1, elab{traink})), [-1 0], optimset('TolFun', 0.001));
end

for k = 1:ncv
    econf{k} = 1 ./ (1+exp(eparams{k}(1)*econf{k}+eparams{k}(2)));
end

elab = cat(1, elab{:});
econf = cat(1, econf{:});

eerror = mean((econf>0.5)~=elab);

econf2 = 1-abs(elab-econf);

ind1 = find(elab==0);
ind2 = find(elab==1);
px = [0.025:0.05:0.975];
f1 = ksdensity(econf(ind1), px, 'support', [0 1]);
f2 = ksdensity(econf(ind2), px, 'support', [0 1]);
fc = ksdensity(econf2, px, 'support', [0 1]);
%fc = fc;

errors.err = eerror;
errors.pneg = f1;
errors.ppos = f2;
errors.conf = fc;
errors.px = px;

medFS = 18;
bigFS = 20;

% figure(1), hold on, plot(px, fc, 'y', 'LineWidth', 2);
% %axis([0 1 0 1])
% xlabel('Confidence in True Label', 'FontSize', medFS)
% ylabel('Frequency', 'FontSize', medFS)
% title('Same Label Confidence', 'FontSize', bigFS)
% set(gca, 'FontSize', medFS)
%
% figure(2), hold on, plot(px, f2 ./ (f1+f2), 'y', 'LineWidth', 2)
% hold on, plot(px, px, '--k')
% axis([0 1 0 1])
% xlabel('Estimated Probability', 'FontSize', medFS)
% ylabel('Empirical Probability', 'FontSize', medFS)
% %title('Same Label Confidence', 'FontSize', bigFS)
% set(gca, 'FontSize', medFS)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function err = objective(param, econf, elab)
econf = 1./(1+exp(param(1)*econf+param(2)));
px = [0.025:0.05:0.975];
f1 = ksdensity(econf(elab==0), px, 'support', [0 1+eps])+eps;
f2 = ksdensity(econf(elab==1), px, 'support', [0 1+eps])+eps;
f1 = f1 / sum(f1+f2);
f2 = f2 / sum(f1+f2);
err = sum((f1+f2).*(px - f2./(f1+f2)).^2);
disp(num2str([sum((f1+f2).*abs(px - f2./(f1+f2))) param]))

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function f = uplus(f)
%+   BALLFUN unary plus.
%   +F of a BALLFUN is F.
%
%   G = UPLUS(F) is called for the syntax '+F'.

% Copyright 2018 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun