text:
function bvalfit()
    % bvalfit Calculates Freq-Mag functions (b-value) for two time-segments
    %
    %   Calculates Freq-Mag functions (b-value) for two time-segments
    %   finds best fit to the foreground for a modified background
    %   assuming a change in time of the following types:
    %   Mnew = Mold + d     , i.e. Simple magnitude shift
    %   Mnew = c*Mold + d   , i.e. Mag stretch plus shift
    %   Nnew = fac*Nold     , i.e. Rate change (N = number of events)
    %                                      R. Zuniga IGF-UNAM/GI-UAF  6/94
    %                                      Rev. 4/2001

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    %TODO DELETE THIS -> Z TOOL REMOVE COMPARE TWO RATES(FIT)

    global p

    report_this_filefun();

    % This is the info window text
    %
    ttlStr='Comparing Seismicity rates ';
    hlpStr1map= ...
        ['                                                '
        ' To be Implemented                              '
        '                                                '];

    if ic == 0
        format short;
        fac = 1.0;

        bvfig = figure;
        set(bvfig,'Units','normalized','NumberTitle','off','Name','b-value curves');
        set(bvfig,'pos',[ 0.435  0.3 0.5 0.5])

        if isempty(ZG.newcat)
            ZG.newcat = a;
        end
        maxmag = max(ZG.newcat.Magnitude);
        mima = min(ZG.newcat.Magnitude);
        if mima > 0
            mima = 0 ;
        end
        [t0b, teb] = bounds(ZG.newcat.Date) ;
        n = ZG.newcat.Count;
        tdiff = round(teb - t0b);

        % number of mag units
        nmagu = (maxmag-mima*10)+1;

        td12 = t2p(1) - t1p(1);
        td34 = t4p(1) - t3p(1);

        l = ZG.newcat.Date > t1p(1) & ZG.newcat.Date < t2p(1) ;
        backg =  ZG.newcat.subset(l);
        [bval,xt2] = hist(backg(:,6),(mima:0.1:maxmag));
        bval = bval/td12;                      % normalization
        bvalsum = cumsum(bval);                        % N for M <=
        bvalsum3 = cumsum(bval(end:-1:1));    % N for M >= (counted backwards)
        magsteps_desc = (maxmag:-0.1:mima);
        [cumux, xt] = hist(ZG.newcat.Date(l),t1p(1):days(ZG.bin_dur):t2p(1));

        l = ZG.newcat.Date > t3p(1) & ZG.newcat.Date < t4p(1) ;
        foreg = ZG.newcat.subset(l);
        bval2 = histogram(foreg(:,6),(mima:0.1:maxmag));
        bval2 = bval2/td34;                     % normallization
        bvalsum2 = cumsum(bval2);
        bvalsum4 = cumsum(bval2(end:-1:1));
        [cumux2, xt] = hist(ZG.newcat.Date(l),t3p(1):days(ZG.bin_dur):t4p(1));
        mean1 = mean(cumux);
        mean2 = mean(cumux2);
        var1 = cov(cumux);
        var2 = cov(cumux2);
        zscore = (mean1 - mean2)/(sqrt(var1/length(cumux)+var2/length(cumux2)));

        backg_be = log10(bvalsum);
        backg_ab = log10(bvalsum3);
        foreg_be = log10(bvalsum2);
        foreg_ab = log10(bvalsum4);

        orient landscape
        rect = [0.2,  0.2, 0.70, 0.70];           % plot Freq-Mag curves
        axes('position',rect)
        semilogy(magsteps_desc,bvalsum3,'om')
        set(gca,'NextPlot','add')
        semilogy(magsteps_desc,bvalsum4,'xb')

        te1 = max([bvalsum  bvalsum2 bvalsum4 bvalsum3]);
        te1 = te1 - 0.2*te1;
        title([file1 '   o: ' num2str(t1p(1)) ' - ' num2str(t2p(1)) '     x: ' num2str(t3p(1)) ' - '  num2str(t4p(1)) ],'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')

        xlabel('Magnitude','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
        ylabel('Cum. Number -normalized','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
        %  find b-values;
        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.s,'Linewidth',1.2)

        figure(mess);
        clf;
        cla;
        set(gcf,'Name','Magnitude selection ');
        set(gca,'visible','off');
        txt5 = text('Position',[.01 0.99 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','Please select two magnitudes to be used');
        txt1 = text('Position',[.01 0.84 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','in the calculation of straight line fit i.e.');
        txt2 = text('Position',[.01 0.66 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','b value of BACKGROUND (o)');

        figure(bvfig);
        seti = uicontrol(...
            'Units','normal','Position',[.4 .01 .2 .05],...
            'String','Select Mag1 ');

        pause(1)

        M1b = ginput(1);
        tx1 = text( M1b(1),M1b(2),['M1'] );
        set(seti,'String','Select Mag2');

        pause(0.1)

        M2b = ginput(1);
        tx2 = text( M2b(1),M2b(2),['M2'] );

        pause(0.1)
        delete(seti)

        ll = magsteps_desc > M1b(1) & magsteps_desc < M2b(1);
        x = magsteps_desc(ll);
        y = backg_ab(ll);
        p  = polyfit(x,y,1);                  % fit a line to background
        f = polyval(p,x);
        f = 10.^f;
        set(gca,'NextPlot','add')
        semilogy(x,f,'r')                         % plot linear fit to backg
        r = corrcoef(x,y);
        r = r(1,2);
        std_backg = std(y - polyval(p,x));      % standard deviation of fit

        figure(mess);
        clf;
        cla;
        set(gcf,'Name','Magnitude selection ');
        set(gca,'visible','off');
        txt5 = text(...
            'Position',[.01 0.99 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','Please select two magnitudes to be used');
        txt1 = text('Position',[.01 0.84 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','in the calculation of straight line fit i.e.');
        txt2 = text(...
            'Position',[.01 0.66 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','b value of FOREGROUND (x)');

        figure(bvfig);
        seti = uicontrol('Units','normal',...
            'Position',[.4 .01 .2 .05],'String','Select Mag1 ');

        pause(1)

        M1f = [];
        M1f = ginput(1);
        tx3 = text( M1f(1),M1f(2),['M1'] )
        set(seti','String','Select Mag2');

        pause(0.1)

        M2f = [];
        M2f = ginput(1);
        tx4 = text( M2f(1),M2f(2),['M2'] )

        pause(0.1)
        delete(seti)

        l = magsteps_desc > M1f(1) & magsteps_desc < M2f(1);
        x = magsteps_desc(l);
        y = foreg_ab(l);
        pp = polyfit(x,y,1);
        % fit a line to foreground
        f = polyval(pp,x);
        f = 10.^f;
        semilogy(x,f,'r')                   % plot fit to foreg
        rr = corrcoef(x,y);
        rr = rr(1,2);
        std_foreg = std(y - polyval(pp,x));      % standard deviation of fit

        figure(mess);
        clf
        set(gca,'visible','off')
        set(gcf,'Units','normalized','pos',[ 0.03  0.1 0.4 0.7])
        set(gcf,'Name','Compare Results');
        orient tall
        te = text(0.,0.99, ['   Catalogue : ' file1]) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        stri = [ 'Background (o):   ' num2str(t1p(1)) '  to  ' num2str(t2p(1)) ];
        te = text(0.01,0.93, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        aa_ = p(2) *1000.0;
        aa_ = round(aa_);
        aa_ = aa_/1000.0;
        bb = p(1) *1000.0;
        bb = round(bb);
        bb = bb/1000.0;          % round to 0.001
        stri = [' Log N = ' num2str(aa_)  num2str(bb) '*M ' ];
        te = text(0.01,0.88, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        stri = [ 'Foreground (x):   ' num2str(t3p(1)) '  to  ' num2str(t4p(1)) ];
        te = text(0.01,0.83, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        aa_ = pp(2) *1000.0;
        aa_ = round(aa_);
        aa_ = aa_/1000.0;
        bb = pp(1) *1000.0;
        bb = round(bb);
        bb = bb/1000.0;          % round to 0.001
        stri = [' Log N = ' num2str(aa_) num2str(bb) '*M '];
        te = text(0.01,0.78, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        disp([' Correlation coefficient for background = ', num2str(r) ]);
        disp([' Correlation coefficient for foreground = ', num2str(rr) ]);
        %  find simple shift
        % first find Mmin ( M for which the background relation
        % departs from straight line by more than std )
        ld = abs(backg_ab - polyval(p,magsteps_desc)) <= std_backg;
        [min_backg, ldb] = min(magsteps_desc(ld));        % Mmin of background
        n1 = backg_ab(ld);
        n1 = n1(ldb);                           % Cum number for Mmin background
        magi = (n1 - pp(2))/pp(1)  % magi is intercept of n1 with foreground linear fit
        dM = magi - min_backg;        % magnitude shift
        ld = abs(foreg_ab - polyval(pp,magsteps_desc)) <= std_foreg;
        [min_foreg, ldf] = min(magsteps_desc(ld));        % min_foreg is Mmin of foreground
        disp([' Mmin for background = ', num2str(min_backg) ]);
        disp([' Mmin for foreground = ', num2str(min_foreg) ]);
        stri = [ 'Minimum magnitude for Background = ' num2str(min_backg) ];
        te = text(0.01,0.73, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        stri = [ 'Minimum magnitude for Foreground = ' num2str(min_foreg) ];
        te = text(0.01,0.68, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        stri = ['Z score between both rates: '];
        te = text(0.01,0.63, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        stri = [' Z = ' num2str(zscore) ];
        te = text(0.01,0.58, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s);
        dM = (round(dM *10.0))/10;     % round to 0.1
        backg_new = [backg(:,1:5), backg(:,6)+dM, backg(:,7)];    %  add shift

        [bvalN,xt2] = hist(backg_new(:,6),(mima:0.1:maxmag));
        bvalN = bvalN/td12;                               % normalize
        bvalsumN = cumsum(bvalN);
        bvalsum3N = cumsum(bvalN(end:-1:1));
        backg_beN = log10(bvalsumN);
        backg_abN = log10(bvalsum3N);

        res =  (sum((bvalN - bval2).^2)/length(bval2))^0.5 ; % residual in histograms
        %%disp(['Average residual of simple shift = ', num2str(res)]);

        figure(mess);
        stri = [ 'Suggested single magnitude shift (d):']
        te = text(0.01,0.50, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold');
        stri = ['Mx = Mo + (', num2str(dM),')']
        te = text(0.01,0.45, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold');
        %  compute magnitude stretch and shift
        pause(0.1)

        mf = p(1)/pp(1);            % factor is calculated from ratio of b values
        mf = (round(mf *100.0))/100.0;    % round to 0.01
        dM = -mf*(pp(2) - p(2))/p(1);   %  find shift by diff of zero ordinates
        dM = (round(dM *100.0))/100.0;    % round to 0.01
        stri = [ 'Linear Mag correction (stretch, c, and shift, d):' ];
        te = text(0.01,0.38, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold');
        stri = [ 'Mx = ',num2str(mf), '* Mo + (', num2str(dM),')' ];
        te = text(0.01,0.33, stri) ;
        set(te,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold');
        set(gca,'NextPlot','add')
    end   % if ic

    if ic == 0 | ic == 2

        figure(bvfig);
        if ic == 2, clf, end
        bvalsumN = [ ];
        bvalsum3N = [ ];
        % Modify Magnitudes
        backg_new = [backg(:,1:5), (mf*backg(:,6))+dM, backg(:,7)];

        [bvalN,xt2] = hist(backg_new(:,6),(mima:0.1:maxmag));
        bvalN = bvalN/td12;                              % normalize
        bvalsumN = cumsum(bvalN);
        bvalsum3N = cumsum(bvalN(end:-1:1));
        backg_beN = log10(bvalsumN);
        backg_abN = log10(bvalsum3N);

        % residual in histograms
        res =  (sum((bvalN - bval2).^2)/length(bval2))^0.5 ;

        % find rate increase_decrease
        if ic ==0 |ic ==1
            rat = [ ];
            rat = bval2/bvalN;            % by mean of ratios
            l = 1 - (isnan(rat) + isinf(rat));
            fac1 = mean(rat(l));
            fac1 = fac1 *100.0;
            fac1 = round(fac1);
            fac1 = fac1/100.0;               % round to 0.01
            fac = fac1;
        end     % if ic

        ind = 0;                      %  find minimum magnitude for the rate change
        resm = [ ];
        if ic ==0
            fac = 1.0 ;
        end
        bvalN = bvalN*fac ;    % apply rate to all data

        bvalsumN = cumsum(bvalN);
        bvalsum3N = cumsum(bvalN(end:-1:1));
        backg_beN = log10(bvalsumN);
        backg_abN = log10(bvalsum3N);

        magi = magi *10.0;
        magi = round(magi);
        magi = magi/10.0;               % round to 0.1
        % residual in histograms
        res =  (sum((bvalN - bval2).^2)/length(bval2))^0.5 ;

        figure(mess);
        if ic == 0 | ic == 1
            stri = [ 'Suggested rate change (Nx = fac*No): \newline fac = ' num2str(fac1) ];
            te = text(0.01,0.27, stri) ;
            set(te,'FontSize',ZmapGlobal.Data.fontsz.s,'Visible','on');
        end   % if ic
        magis = maxmag;

        uicontrol('Units','normal','Position',[.88 .9 .11 .06],'String','Print  ', 'callback',@callbackfun_001)
        uicontrol('Units','normal','Position',[.88 .80 .11 .06],'String','Close  ', 'callback',@callbackfun_002)

        freq_field1=uicontrol('Style','edit',...
            'Position',[.30 .16 .13 .07],...
            'Units','normalized','String',num2str(dM),...
            'callback',@callbackfun_003);

        freq_field2=uicontrol('Style','edit',...
            'Position',[.75 .16 .13 .07],...
            'Units','normalized','String',num2str(mf),...
            'callback',@callbackfun_004);

        freq_field3=uicontrol('Style','edit',...
            'Position',[.30 .05 .13 .07],...
            'Units','normalized','String',num2str(fac),...
            'callback',@callbackfun_005);

        txt1 = text(...
            'Color',[0 0 0 ],...
            'Position',[.01 0.11 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.s ,...
            'String','Shift (d)');

        txt2 = text(...
            'Color',[0 0 0 ],...
            'Position',[.44 0.11 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.s ,...
            'String','  Stretch factor (c)');

        txt3 = text(...
            'Color',[0 0 0 ],...
            'Position',[.01 0.0 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.s ,...
            'String','Rate factor');

        go_button=uicontrol('Style','Pushbutton',...
            'Position',[.52 .01 .10 .07 ],...
            'Units','normalized',...
            'callback',@callbackfun_006,...
            'String','Go');
        res = bval2 - bvalN ;

        close(bvfig)
        % Find out if figure already exists and plot results of fit
        %
        bvfig=findobj('Type','Figure','-and','Name','Compare and fit two rates');

        ms3 = 5;

        %if isempty(bvfig)
        bvfig= figure_w_normalized_uicontrolunits( ...
            'Name','Compare and fit two rates',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'Visible','on', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)+200));

        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Print ',...
            'callback',@callbackfun_007)

        uicontrol('Units','normal',...
            'Position',[.0 .75 .08 .06],'String','Close ',...
            'callback',@callbackfun_008)

        uicontrol('Units','normal',...
            'Position',[.0 .85 .08 .06],'String','Info ',...
            'callback',@callbackfun_009)
        axis off
        %end % if figure exits

        %%figure(bvfig);
        delete(findobj(bvfig,'Type','axes'));
        % plot b-value plot
        %
        orient tall
        set(gcf,'PaperPosition',[2 1 5.5 7.5])
        rect = [0.20,  0.7, 0.70, 0.25];           % plot Freq-Mag curves
        axes('position',rect)
        set(gca,'NextPlot','add')
        figure(bvfig);
        set(gca,'NextPlot','add')
        pl = semilogy(magsteps_desc,bvalsum4,'xb');
        set(gca,'Yscale','log')
        set(gca,'NextPlot','add')
        set(pl,'MarkerSize',ms3)
        semilogy(magsteps_desc,bvalsum4,'-.b')
        pl = semilogy(magsteps_desc,bvalsum3N,'om');
        set(pl,'MarkerSize',ms3)
        semilogy(magsteps_desc,bvalsum3N,'m')
        te1 = max([bvalsum  bvalsum2 bvalsum4 bvalsum3]);
        te1 = te1 - 0.2*te1;

        ylabel('Cum. rate/year','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
        str = [ '   o: ' num2str(t1p(1),6) ' - ' num2str(t2p(1),4) '     x: ' num2str(t3p(1),6) ' - '  num2str(t4p(1),6) ];

        title(str,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
        %  find b-values;
        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.s,'Linewidth',1.0)
        p1 = gca;

        % Plot histogram
        %

        rect = [0.20,  0.40 0.70, 0.25];
        axes('position',rect)
        pl = plot(xt2,bvalN,'om');
        set(pl,'MarkerSize',ms3,'LineWidth',1.0)
        set(gca,'NextPlot','add')
        pl = plot(xt2,bval2,'xb');
        set(pl,'MarkerSize',ms3,'LineWidth',1.0)
        pl = plot(xt2,bval2,'-.b');
        set(pl,'MarkerSize',ms3,'LineWidth',1.0)
        pl = plot(xt2,bvalN,'m');
        set(pl,'MarkerSize',ms3,'LineWidth',1.0)
        disp([' Summation: ' num2str(sum(bval-bval2))])
        v = axis;
        xlabel('Magnitude ','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
        ylabel('rate/year','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.s,'Linewidth',1.0)

        uic = uicontrol('Units','normal','Position',[.35 .15 .30 .07],'String','Magnitude Signature? ', 'callback',@callbackfun_010);

    end   % if ic

    clear rat bvalNN mean1 mean2 ld l ll txt1 txt2 txt3 txt4 M1b M2b M1f M2f tx1 tx2 tx3 tx4;
    ic = 0;
    format;

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        printdlg;
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));

    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dM=str2double(freq_field1.String);
        freq_field1.String=num2str(dM);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        mf=str2double(freq_field2.String);
        freq_field2.String=num2str(mf);
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        fac=str2double(freq_field3.String);
        freq_field3.String=num2str(fac);
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ic = 2;
        bvalfit;
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        printdlg;
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        f1=gcf;
        f2=gpf;
        set(f1,'Visible','off');
    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        zmaphelp(ttlStr,hlpStr1map,hlpStr2map,hlpStr3map);
    end

    function callbackfun_010(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        delete(uic);
        synsig;
    end

end

author: CelsoReyes


text:
function bcdf = beta_cdf (x, a, b)
%  BETACDF Beta cumulative distribution function.
%     P = BETACDF(X,A,B) returns the beta cumulative distribution
%     function with parameters A and B at the values in X.

% Copyright (c) 1995-1997, 2005-2007 Kurt Hornik

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.
if (nargin ~= 3)
  error('must provide 3 parameters')
end

if (~isscalar (a) || ~isscalar(b))
  if ((size(a,1) ~= size(b,1)) || (size(b,1) ~= size(x,1)))
    error ('betainv: x, a and b must be of common size or scalars');
  end
end

[n,nin] = size(x);
bcdf = zeros(n,nin);

k = find ((a < 0) || (b < 0) || isnan (x));
if (any(k))
  bcdf(k) = NaN;
end

k = find((x >= 1) && (a > 0) && (b > 0));
if (any(k))
  bcdf(k) = 1;
end

k = find ((x > 0) && (x < 1) && (a > 0) && (b > 0));
if (any (k))
  if (isscalar(a) && isscalar(b))
    bcdf(k) = betainc(x(k), a, b);
  else
    bcdf(k) = betainc(x(k), a(k), b(k));
  end
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cdf = student_noncentral_cdf ( x, idf, d )

%*****************************************************************************80
%
%% STUDENT_NONCENTRAL_CDF evaluates the noncentral Student T CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Algorithm AS 5,
%    Applied Statistics,
%    Volume 17, 1968, page 193.
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, integer IDF, the number of degrees of freedom.
%
%    Input, real D, the noncentrality parameter.
%
%    Output, real CDF, the value of the CDF.
%
  a_max = 100;
  emin = 12.5;

  f = idf;

  if ( idf == 1 )

    a = x / sqrt ( f );
    b = f / ( f + x * x );
    drb = d * sqrt ( b );

    cdf2 = normal_01_cdf ( drb );
    cdf = 1.0 - cdf2 + 2.0 * tfn ( drb, a );

  elseif ( idf <= a_max )

    a = x / sqrt ( f );
    b = f / ( f + x * x );
    drb = d * sqrt ( b );
    sum2 = 0.0;

    fmkm2 = 0.0;
    if ( abs ( drb ) < emin )
      cdf2 = normal_01_cdf ( a * drb );
      fmkm2 = a * sqrt ( b ) * exp ( -0.5 * drb * drb ) * cdf2 / sqrt ( 2.0 * pi );
    end

    fmkm1 = b * d * a * fmkm2;
    if ( abs ( d ) < emin )
      fmkm1 = fmkm1 + 0.5 * b * a * exp ( - 0.5 * d * d ) / pi;
    end

    if ( mod ( idf, 2 ) == 0 )
      sum2 = fmkm2;
    else
      sum2 = fmkm1;
    end

    ak = 1.0;

    for k = 2 : 2 : idf - 2

      fk = k;

      fmkm2 = b * ( d * a * ak * fmkm1 + fmkm2 ) * ( fk - 1.0 ) / fk;

      ak = 1.0 / ( ak * ( fk - 1.0 ) );
      fmkm1 = b * ( d * a * ak * fmkm2 + fmkm1 ) * fk / ( fk + 1.0 );

      if ( mod ( idf, 2 ) == 0 )
        sum2 = sum2 + fmkm2;
      else
        sum2 = sum2 + fmkm1;
      end

      ak = 1.0 / ( ak * fk );

    end

    if ( mod ( idf, 2 ) == 0 )
      cdf2 = normal_01_cdf ( d );
      cdf = 1.0 - cdf2 + sum2 * sqrt ( 2.0 * pi );
    else
      cdf2 = normal_01_cdf ( drb );
      cdf = 1.0 - cdf2 + 2.0 * ( sum2 + tfn ( drb, a ) );
    end
%
%  Normal approximation.
%
  else

    a = sqrt ( 0.5 * f ) * exp ( gammaln ( 0.5 * ( f - 1.0 ) ) ...
      - gammaln ( 0.5 * f ) ) * d;

    temp = ( x - a ) / sqrt ( f * ( 1.0 + d * d ) / ( f - 2.0 ) - a * a );

    cdf2 = normal_01_cdf ( temp );
    cdf = cdf2;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = get(f, propName)
%GET   GET method for the CHEBOP2 class.
%   P = GET(N, PROP) returns the property P specified in the string PROP from
%   the CHEBOP2 N. Valid entries for the string PROP are:
%       'DOMAIN'         - The domain of definition of N.
%       'OP'             - The partial differential operator of N.
%       'LBC'            - The left boundary constraints of N.
%       'RBC'            - The right boundary constraints of N.
%       'UBC'            - The top boundary constraints of N.
%       'DBC'            - The bottom boundary constraints of N.
%
% The following are also supported:
%       'DIM'            - The size of the system.
%       'SCALE'          - The relative scale of the solution.
%       'COEFFS'         - The variable coefficients of N.
%       'XORDER'         - The order of differentiation in the x-direction.
%       'YORDER'         - The order of differentiation in the y-direction.
%       'U', 'S', 'V'    - The low rank structure of N.

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

if ( numel(f) > 1 )
    out = cell(numel(f));
    for k = 1:numel(f)
        out{k} = get(f(k), propName);
    end
    return
end

switch propName
    case fieldnames(f)
        % Allow access to any of the properties of F via GET:
        out = f.(propName);
    otherwise
        error('CHEBFUN:CHEBOP2:get:propName', ...
            [propName,' is not a valid chebop2 property.'])
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%
% ENVTOPO - Plot the envelope of a multichannel data epoch, plus envelopes and
%             scalp maps of specified or largest-contributing components. If a 3-D
%             input matrix, operates on the mean of the data epochs. Click on
%             individual axes to examine them in detail. The black lines represent
%             the max and min values across all channels at each time point.
%             The blue shading represents the max and min contributions of
%             the selected components tothose channels. The paired colored lines
%             represent the max and min contributions of the indicated component
%             across all channels.
% Usage:
%             >> envtopo(data,weights,'chanlocs',file_or_struct);
%             >> [cmpvarorder,cmpvars,cmpframes,cmptimes,cmpsplotted,sortvar] ...
%                                     = envtopo(data, weights, 'key1', val1, ...);
% Inputs:
%  data        = single data epoch (chans,frames), else it a 3-D epoched data matrix
%                (chans,frames,epochs) ->  processes the mean data epoch.
%  weights     = linear decomposition (unmixing) weight matrix. The whole matrix
%                should be passed to the function here.
%                    Ex: (EEG.icaweights*EEG.icasphere)
% Required keyword:
%  'chanlocs'  = [string] channel location filename or EEG.chanlocs structure.
%                For more information, see >> topoplot example
%
% Optional inputs:
%  'compnums'  = [integer array] vector of component indices to use in the
%                  calculations and to select plotted components from
%                  {default|[]: all}
%  'compsplot' = [integer] the number of largest contributing components to plot.
%                  compnums in the latter case is restricted in size by the
%                  internal variable MAXTOPOS (20) {default|[] -> 7}
%  'subcomps'  = [integer vector] indices of comps. to remove from the whole data
%                  before plotting. 0 -> Remove none. [] -> If 'compnums'
%                  also listed, remove *all* except 'compnums' {default: 0}
%  'limits'    = 0 or [minms maxms] or [minms maxms minuV maxuV]. Specify
%                  start/end plot (x) limits (in ms) and min/max y-axis limits
%                  (in uV). If 0, or if both minmx & maxms == 0 -> use latencies
%                  from 'timerange' (else 0:frames-1). If both minuV and
%                  maxuV == 0 -> use data uV limits {default: 0}
%  'timerange' = start and end input data latencies (in ms) {default: from 'limits'
%                  if any}. Note: Does NOT select a portion of the input data,
%                  just makes time labels.
%  'limcontrib' = [minms maxms]  time range (in ms) in which to rank component
%                  contribution (boundaries shown with thin dotted lines)
%                  {default|[]|[0 0] -> plotting limits}
%  'sortvar'   = ['mp'|'pv'|'pp'|'rp'] {default:'mp'}
%                  'mp', sort components by maximum mean back-projected power
%                  in the 'limcontrib' time range: mp(comp) = max(Mean(back_proj.^2));
%                  where back_proj = comp_map * comp_activation(t) for t in 'limcontrib'
%                  'pv', sort components by percent variance accounted for (EEG_PVAF)
%                    pvaf(comp) = 100-100*mean(var(data - back_proj))/mean(var(data));
%                  'pp', sort components by percent power accounted for (ppaf)
%                    ppaf(comp) = 100-100*Mean((data - back_proj).^2)/Mean(data.^2);
%                  'rp', sort components by relative power
%                    rp(comp) = 100*Mean(back_proj.^2)/Mean(data.^2);
%  'title'     = [string] plot title {default|[] -> none}
%  'plotchans' = [integer array] data channels to use in computing contributions and
%                  envelopes, and also for making scalp topo plots
%                  {default|[] -> all}, by calling TOPOPLOT.
%  'voffsets'  = [float array] vertical line extensions above the data max to
%                  disentangle plot lines (left->right heads, values in y-axis units)
%                  {default|[] -> none}
%  'colors'    = [string] filename of file containing colors for envelopes, 3 chars
%                  per line, (. = blank). First color should be "w.." (white)
%                  Else, 'bold' -> plot default colors in thick lines.
%                  {default|[] -> standard Matlab color order}
%  'fillcomp'  = int_vector>0 -> fill the numbered component envelope(s) with
%                  solid color. Ex: [1] or [1 5] {default|[]|0 -> no fill}
% 'plotproj'   = [1|0] Plot channel projections if only one component projection is
%                   plotted. {default: 0}
% 'fillcolor'  = Three elements RGB vector of the color to fill the component envelope
%                   if 'fillcomp' is set. {default|[.815 .94 1] %light blue}
%  'vert'      = vector of times (in ms) at which to plot vertical dashed lines
%                  {default|[] -> none}
%  'icawinv'   = [float array] inverse weight matrix. Normally computed by inverting
%                  the weights*sphere matrix (Note: If some components have been
%                  removed, the pseudo-inverse does not represent component maps
%                  accurately).
%  'icaact'    = [float array] component activations. {default: computed from the
%                  input weight matrix}
%  'envmode'   = ['avg'|'rms'] compute the average envelope or the root mean square
%                  envelope {default: 'avg'}
%  'sumenv'    = ['on'|'off'|'fill'] 'fill' -> show the filled envelope of the summed
%                  projections of the selected components; 'on' -> show the envelope
%                  only {default: 'fill'}
%  'actscale'  = ['on'|'off'] scale component scalp maps by maximum component activity
%                  in the designated (limcontrib) interval. 'off' -> scale scalp maps
%                  individually using +/- max(abs(map value)) {default: 'off'}
%  'dispmaps'  = ['on'|'off'] display component numbers and scalp maps {default: 'on'}
%  'topoplotkey','val' = optional additional TOPOPLOT arguments {default: none}
%  'axisoff'   = [real] percent of the figure y-dimension to be covered by the axis
%                 for the envelopes. Values are restricted to the range [0.60 1]{default: 0.6}
%
% 'limcontribweight' = [real] lineweight of limonctrib vertical lines.{default: 1.2}
% 'vertweight'       = [real] lineweight of specified vertical lines {default: 2}
%
% Outputs:
%
%  cmpvarorder = component numbers in decreasing order of max variance in data
%  cmpvars     = ('sortvar') max power for all components tested, sorted
%                   from highest to lowest. See 'sortvar' 'mp'
%  cmpframes   = frames of comvars for each component plotted
%  cmptimes    = times of compvars for each component plotted
%  cmpsplotted = indices of components plotted.
%                       unsorted_compvars = compvars(compsplotted);
%  sortvar      = The computed data used to sort components with.
%                 See 'sortvar' option above.
%
% Notes:
%  To label maps with other than component numbers, put four-char strings into
%  a local (pwd) file named 'envtopo.labels' (using . = space) in time-order
%  of their projection maxima
%
% Authors: Scott Makeig & Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/1998
%
% See also: TIMTOPO AXCOPY

% Copyright (C) 3-10-98 from timtopo.m Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.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.

% Edit History:
% 3-18-98 fixed bug in LineStyle for fifth component, topoplot maxproj with
%         correct orientations, give specified component number labels -sm
% 4-28-98 plot largest components, ranked by max projected variance -sm
% 4-30-98 fixed bug found in ICADEMO -sm
% 5-08-98 fixed bug found by mw () -sm
% 5-23-98 made vert. line styles for comps 6 & 11 correct -sm
% 5-30-98 added 'envtopo.labels' option -sm
% 5-31-98 implemented plotchans arg -sm
% 7-13-98 gcf->gca to allow plotting in subplots -sm
% 9-18-98 worked more to get plotting in subplot to work -- no luck yet! -sm
% 2-22-99 draw oblique line to max env value if data clipped at max proj -sm
% 2-22-99 added colorfile -sm
% 4-17-99 added support for drawing in subplots -t-pj
% 10-29-99 new versions restores search through all components for max 7 and adds
%          return variables (>7 if specified. Max of 7 comp envs still plotted. -sm
% 11-17-99 debugged new version -sm
% 12-27-99 improved help msg, moved new version to distribution -sm
% 01-21-00 added 'bold' option for colorfile arg -sm
% 02-28-00 added fill_comp_env arg -sm
% 03-16-00 added AXCOPY -sm & tpj
% 05-02-00 added vert option -sm
% 05-30-00 added option to show "envelope" of only 1 channel -sm
% 09-07-00 added [-n] option for compnums, added BOLD_COLORS as default -sm
% 12-19-00 updated ICAPROJ args -sm
% 12-22-00 trying 'axis square' for topoplots -sm
% 02-02-01 fixed bug in printing component 6 env line styles -sm
% 04-11-01 added [] default option for args -sm
% 01-25-02 reformated help & license, added links -ad
% 03-15-02 added readlocs and the use of eloc input structure -ad
% 03-16-02 added all topoplot options -ad

function [compvarorder,compvars,compframes,comptimes,compsplotted,sortvar] = envtopo(data,weights,varargin)

% icadefs;    % read toolbox defaults

sortvar = [];
all_bold = 0;
BOLD_COLORS = 1;    % 1 = use solid lines for first 5 components plotted
                    % 0 = use std lines according to component rank only
% FILL_COMP_ENV = 0;  % default no fill
%FILLCOLOR   = [.815 .94 1]; % use lighter blue for better env visibility
MAXTOPOS = 20;      % max topoplots to plot
% VERTWEIGHT = 2.0;  % lineweight of specified vertical lines
% LIMCONTRIBWEIGHT = 1.2; % lineweight of limonctrib vertical lines
MAX_FRAMES = 10000; % maximum frames to plot
MAXENVPLOTCHANS = 10;
xmin = 0; xmax = 0;

if nargin < 2
   help envtopo
   return
end

if nargin <= 2 || ischar(varargin{1})
	% 'key' 'val' sequences
	fieldlist = { 'chanlocs'      ''         []                       [] ;
				  'title'         'string'   []                       '';
				  'limits'        'real'     []                       0;
				  'timerange'     'real'     []                       [];
				  'plotchans'     'integer'  [1:size(data,1)]         [] ;
				  'icawinv'       'real'     []                       pinv(weights) ;
				  'icaact'        'real'     []                       [] ;
				  'voffsets'      'real'     []                       [] ;
				  'vert'          'real'     []                       [] ;
				  'fillcomp'      'integer'  []                       0 ;
				  'figure'        'integer'  []                       [] ;
				  'colorfile'     'string'   []                       '' ;
				  'colors'        'string'   []                       '' ;
				  'compnums'      'integer'  []                       [];
				  'compsplot'     'integer'  []                       7;
				  'subcomps'      'integer'  []                       0;
				  'envmode'       'string'   {'avg','rms'}            'avg';
				  'dispmaps'      'string'   {'on','off'}             'on';
				  'pvaf'          'string'   {'mp','mv','on','rp','rv','pv','pvaf','pp','off',''} '';
				  'sortvar'       'string'   {'mp','mv','rp','rv','pv','pvaf','pp'} 'mp';
				  'actscale'      'string'   {'on','off'}             'off';
				  'limcontrib'    'real'     []                       0;
				  'topoarg'       'real'     []                       0;
				  'sumenv'        'string'   {'on','off','fill'}      'fill';
                  'axisoff'       'real'     [0.6 1]               0.6;
                  'fillcolor'     'real'     []                       [.815 .94 1];
                  'limcontribweight'  'real' []                       1.2;
                  'plotproj'          'real' [0 1]                    0;
                  'vertweight'        'real' []                       2};

	% Note: Above, previous 'pvaf' arguments 'on' -> 'pv', 'off' -> 'rv'
	%       for backwards compatibility 11/2004 -sm

	[g varargin] = finputcheck( varargin, fieldlist, 'envtopo', 'ignore');
	if ischar(g), error(g); end
    if g.plotproj && strcmp(g.sumenv, 'fill'), g.sumenv = 'on'; end

else % dprecated - old style input args
	if nargin > 3,    g.chanlocs = varargin{1};
	else              g.chanlocs = [];
	end
        if ischar(varargin{2}), help envtopo; return; end
	if nargin > 4,	  g.limits = varargin{2};
	else              g.limits = 0; % [];
	end
        if ischar(varargin{3}), help envtopo; return; end
	if nargin > 5,    g.compnums = varargin{3};
	else              g.compnums = [];
	end
        if ~ischar(varargin{4}), help envtopo; return; end
	if nargin > 6,    g.title = varargin{4};
	else              g.title = '';
	end
        if ischar(varargin{5}), help envtopo; return; end
	if nargin > 7,    g.plotchans = varargin{5};
	else              g.plotchans = [];
	end
        if ischar(varargin{6}), help envtopo; return; end
	if nargin > 8,    g.voffsets = varargin{6};
	else              g.voffsets = [];
	end
        if ischar(varargin{7}), help envtopo; return; end
	if nargin > 9,    g.colorfile = varargin{7};
	else              g.colorfile = '';
	                  g.colors = '';
	end
        if ischar(varargin{8}), help envtopo; return; end
	if nargin > 10,   g.fillcomp = varargin{8};
	else              g.fillcom = 0;
	end
        if ischar(varargin{9}), help envtopo; return; end
	if nargin > 11,   g.vert = varargin{9};
	else              g.vert = [];
	end
    if nargin > 12, varargin =varargin(10:end); end
    g.sumenv = 'on';
    g.sortvar = 'mp';
    g.pvaf = [];
    g.timerange = [];
    g.icaact = [];
    g.limcontrib = 0;
    g.icawinv = pinv(weights);
    g.subcomps = 0;
    g.envmode = 'avg';
    g.dispmaps = 'on';
end
if ~isempty(g.figure), figure(g.figure); end;        % remember the current figure (for Matlab 7.0.0 bug)

if ~isempty(g.pvaf)
	g.sortvar = g.pvaf; % leave deprecated g.pvaf behind.
end

if strcmpi(g.sortvar,'on') || strcmpi(g.sortvar,'pvaf') || strcmpi(g.sortvar,'mv')
   g.sortvar = 'mp'; % update deprecated flags
end
if strcmpi(g.sortvar,'off') || strcmp(g.sortvar,'rv')
   g.sortvar = 'rp';
end

%
% Check input flags and arguments
%
if ndims(data) == 3
    data = mean(data,3); % average the data if 3-D
end
[chans,frames] = size(data);

if frames > MAX_FRAMES
   error('number of data frames to plot too large!');
end

if ischar(g.chanlocs)
    g.chanlocs = readlocs(g.chanlocs);  % read channel location information
    if length(g.chanlocs) ~= chans
     fprintf(...
      'envtopo(): locations for the %d data channels not in the channel location file.\n', ...
        chans);
        return
    end
end
% Checking dimension of chanlocs
if length(g.chanlocs) > size(data,1)
    fprintf(2,['\n envtopo warning: Dimensions of data to plot and channel location are not consistent\n' ...
        '                 As a result, the plots generated can be WRONG\n'...
        '                 If you are providing the path to your channel location file. Make sure\n'...
        '                 to load the file first, and to provide as an input only the channels that\n'...
        '                 correspond to the data to be plotted, otherwise just provide a valid chanlocs\n ']);
elseif length(g.chanlocs) > size(data,1)
        fprintf(2,['\n envtopo error: Dimensions of data to plot and channel location are not consistent\n' ...
        '                 Check the dimension of the channel location provided\n'...
        '                 Aborting plot ...\n']);
    exit;
end

if ~isempty(g.colors)
    g.colorfile = g.colors; % retain old usage 'colorfile' for 'colors' -sm 4/04
end

if ~isempty(g.vert)
   g.vert = g.vert/1000; % convert from ms to s
end
%
%%%%%% Collect information about the gca, then delete it %%%%%%%%%%%%%
%
uraxes = gca; % the original figure or subplot axes
pos=get(uraxes,'Position');
axcolor = get(uraxes,'Color');
delete(gca)

%
%%% Convert g.timerange, g.limits and g.limcontrib to sec from ms %%%%
%
g.timerange = g.timerange/1000;   % the time range of the input data
g.limits(1) = g.limits(1)/1000;   % the time range to plot
if length(g.limits) == 1   % make g.limits at least of length 2
    g.limits(1) = 0; g.limits(2) = 0;
else
    g.limits(2) = g.limits(2)/1000;  %
end
g.limcontrib = g.limcontrib/1000; % the time range in which to select largest components

%
%%%%%%%%%%%% Collect time range information %%%%%%%%%%%%%%%%%%%%%%%%%%
%
if length(g.limits) == 3 || length(g.limits) > 4 % if g.limits wrong length
   fprintf('envtopo: limits should be 0, [minms maxms], or [minms maxms minuV maxuV].\n');
end

xunitframes = 0; % flag plotting if xmin & xmax are in frames instead of sec
if ~isempty(g.timerange)   % if 'timerange' given
    if g.limits(1)==0 && g.limits(2)==0
         g.limits(1) = min(g.timerange); % if no time 'limits
         g.limits(2) = max(g.timerange); % plot whole data epoch
    end
else % if no 'timerange' given
    if g.limits(1)==0 && g.limits(2)==0 % if no time limits as well,
         fprintf('\nNOTE: No time limits given: using 0 to %d frames\n',frames-1);
         g.limits(1) = 0;
         g.limits(2) = frames-1;
         xunitframes     = 1; % mark frames as time unit instead of sec
    end
end

if isempty(g.timerange)
  xmin = g.limits(1); % (xmin, xmax) are data limits in sec
  xmax = g.limits(2);
else
  xmin = g.timerange(1); % (xmin, xmax) are data limits in sec
  xmax = g.timerange(2);
end

pmin = g.limits(1); % plot min and max sec
if pmin < xmin
   pmin = xmin;     % don't allow plotting beyond the data limits
end
pmax = g.limits(2);
if pmax > xmax
   pmax = xmax;
end

dt = (xmax-xmin)/(frames-1);  % sampling interval in sec
times=xmin*ones(1,frames)+dt*(0:frames-1); % time points in sec

%
%%%%%%%%%%%%%%% Find limits of the component selection window %%%%%%%%%
%
if any(g.limcontrib ~= 0)
    if xunitframes
       g.limcontrib = g.limcontrib*1000; % if no time limits, interpret
    end                                  % limcontrib as frames
    if g.limcontrib(1)<xmin
          g.limcontrib(1) = xmin;
    end
    if g.limcontrib(2)>xmax
          g.limcontrib(2) = xmax;
    end
    srate = (frames-1)/(xmax-xmin);
    limframe1  = round((g.limcontrib(1)-xmin)*srate)+1;
    limframe2  = round((g.limcontrib(2)-xmin)*srate)+1;
    g.vert(end+1) =  g.limcontrib(1);
    g.vert(end+1) =  g.limcontrib(2);
else
    limframe1 = 1;
    limframe2 = frames;
end

%
%%%%%%%%%%%%%%%%%%%%% Read line color information %%%%%%%%%%%%%%%%%%%%%
%
ENVCOLORS = strvcat('w..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..','m..','c..','r..','g..','b..');

if isempty(g.colorfile)
    g.colorfile = ENVCOLORS; % use default color order above
elseif ~ischar(g.colorfile)
      error('Color file name must be a string.');
end
if strcmpi(g.colorfile,'bold')
       all_bold = 1;
       g.colorfile = ENVCOLORS; % default colors
end
if exist(g.colorfile) == 2  % if an existing file
        cid = fopen(g.colorfile,'r');
        if cid <3,
            error('cannot open color file');
        else
            colors = fscanf(cid,'%s',[3 MAXENVPLOTCHANS]);
            colors = colors';
        end
else
        colors = g.colorfile;
end
[r c] = size(colors);
      for i=1:r
        for j=1:c
          if colors(i,j)=='.',
            if j==1
              error('Color file should have color letter in 1st column.');
            elseif j==2
              colors(i,j)='-';
            elseif j>2
              colors(i,j)=' ';
            end
          end
        end
      end
colors(1,1) = 'k'; % make sure 1st color (for data envelope) is black

%
%%%%%%%%%%%%%%%% Check other input variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[wtcomps,wchans] = size(weights);
if wchans ~= chans
     error('Sizes of weights and data do not agree');
end
if wtcomps ~= chans
   fprintf('Number of components not the same as number of channels.\n');
   fprintf('  - component scalp maps and time courses may not be correct.\n');
end

if isempty(g.voffsets) || ( size(g.voffsets) == [1,1] && g.voffsets(1) == 0 )
    g.voffsets = zeros(1,MAXTOPOS);
end
if isempty(g.plotchans) || g.plotchans(1)==0
    g.plotchans = 1:chans;
end
if max(g.plotchans) > chans || min(g.plotchans) < 1
    error('invalid ''plotchan'' index');
end

if g.compsplot < 0
   g.compsplot = abs(g.compsplot);
end

if g.compnums < 0 % legacy syntax
   g.compsplot = abs(g.compnums);
   g.compnums = [];
end
if isempty(g.compnums) || g.compnums(1) == 0
    g.compnums = 1:wtcomps; % by default, select from all components
end

if g.compsplot > MAXTOPOS
    fprintf('Can only plot a maximum of %d components.\n',MAXTOPOS);
    return
else
    MAXTOPOS = g.compsplot;
end

if max(g.compnums) > wtcomps || min(g.compnums)< 1
    error('Keyword ''compnums'' out of range (1 to %d)', wtcomps);
end

g.subcomps = abs(g.subcomps); % don't pass negative channel numbers
if max(g.subcomps) > wtcomps
    error('Keyword ''subcomps'' argument out of bounds');
end

%
%%%%%%%%%%%%%%% Subtract components from data if requested %%%%%%%%%%%%%
%

ncomps = length(g.compnums);

if isempty(g.subcomps) % remove all but compnums
     g.subcomps = 1:wtcomps;
     g.subcomps(g.compnums) = [];
else
  % g.subcomps 0    -> subtract no comps
  %            list -> subtract subcomps list
  if min(g.subcomps) < 1
      if length(g.subcomps) > 1
         error('Keyword ''subcomps'' argument incorrect.');
      end
      g.subcomps = [];   % if subcomps contains a 0 (or <1), don't remove components
  elseif max(g.subcomps) > wtcomps
      error('Keyword ''subcomps'' argument out of bounds.');
  end
end

g.icaact = weights*data;
if ~isempty(g.subcomps)
  fprintf('Subtracting requested components from plotting data: ');
  for k = 1:length(g.subcomps)
      fprintf('%d ',g.subcomps(k));
      if ~rem(k,32)
         fprintf('\n');
      end
  end
  fprintf('\n');
  g.icaact(g.subcomps,:) = zeros(length(g.subcomps),size(data,2));
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Process components %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

for i=1:ncomps-1
    if g.compnums(i) == 0
	       fprintf('Removing component number 0 from compnums.\n');
	       g.compnums(i)=[];
    elseif g.compnums(i)>wtcomps
	       fprintf('compnums(%d) > number of comps (%d)?\n',i,wtcomps);
	       return
    end
    for j=i+1:ncomps
	    if g.compnums(i)==g.compnums(j)
	       fprintf('Removing repeated component number (%d) in compnums.\n',g.compnums(i));
	       g.compnums(j)=[];
            end
    end
end

limitset = 0;
if isempty(g.limits)
  g.limits = 0;
end
if length(g.limits)>1
  limitset = 1;
end

%
%%%%%%%%%%%%%%% Compute plotframes and envdata %%%%%%%%%%%%%%%%%%%%%
%
ntopos = length(g.compnums);
if ntopos > MAXTOPOS
      ntopos = MAXTOPOS; % limit the number of topoplots to display
end

if max(g.compnums) > wtcomps || min(g.compnums)< 1
      fprintf(...
       'envtopo(): one or more compnums out of range (1,%d).\n',wtcomps);
      return
end

plotframes = ones(ncomps);

% toby 2.16.2006: maxproj will now contain all channel info, in case
%                 plotgrid is called in topoplot.
%                 NOT maxproj = zeros(length(g.plotchans),ncomps);
maxproj = zeros(chans,ncomps);

%
% first, plot the data envelope
%
envdata = zeros(2,frames*(ncomps+1));
envdata(:,1:frames) = envelope(g.icawinv(g.plotchans,:)*g.icaact, g.envmode);

fprintf('Data epoch is from %.0f ms to %.0f ms.\n',1000*xmin,1000*xmax);
fprintf('Plotting data from %.0f ms to %.0f ms.\n',1000*xmin,1000*xmax);
fprintf('Comparing maximum projections for components:  \n');
        if ncomps>32
           fprintf('\n');
        end

compvars = zeros(1,ncomps);
mapsigns = zeros(1,ncomps);

%
% Compute frames to plot
%
sampint  = (xmax-xmin)/(frames-1);     % sampling interval in sec
times    = xmin:sampint:xmax;          % make vector of data time values

[v minf] = min(abs(times-pmin));
[v maxf] = min(abs(times-pmax));
pframes  = minf:maxf;         % frames to plot
ptimes   = times(pframes);    % times to plot
if limframe1 < minf
   limframe1 = minf;
end
if limframe2 > maxf
   limframe2 = maxf;
end

%
%%%%%%%%%%%%%% find max variances and their frame indices %%%%%%%%%%%
%

if strcmp(g.sortvar,'pv') %Changed -Jean
	% Variance of the data in the interval, for calculating sortvar.
	vardat = mean(var(data(g.plotchans,limframe1:limframe2),1));
else
	% Compute data rms for sortvar
	powdat = mean(mean(data(g.plotchans,limframe1:limframe2).^2));
end

for c = 1:ncomps

      if isempty(g.icaact) % make the back-projection of component c

          % Changed to include all channels in computation for use in
          % topoplot, particularly with plotgrid option. toby 2.16.2006
          proj = g.icawinv(:,g.compnums(c))*weights(g.compnums(c),:)*data;
      else
          proj = g.icawinv(:,g.compnums(c))*g.icaact(g.compnums(c),:);
      end;

      % save the comp envelope for plotting component waveforms
      envdata(:,c*frames+1:(c+1)*frames) = envelope(proj(g.plotchans,:), g.envmode);

      % Find the frame (timepoint) of largest rms component value
      % and the relative value to those channels defined by plotchans.
      if length(g.plotchans) > 1
        [maxval,maxi] = max(mean((proj(g.plotchans,limframe1:limframe2)).^2));
      else % g.plotchans == 1 --> find absmax value
        [maxval,maxi] = max((abs(proj(g.plotchans,limframe1:limframe2))));
      end
      maxi = maxi+limframe1-1;

      % plotframes and compvars are needed for plotting the lines indicating
      % the timepoint a topoplot refers to.
      plotframes(c) = maxi;
      compvars(c)   = maxval;       % find value of max variance for comp c
      maxproj(:,c)  = proj(:,maxi); % maxproj now contains all channels, to handle
                                    % the 'plotchans'/topoplot'gridplot' conflict.
									% Toby 2.17.2006
      %
      %%%%%% Calculate sortvar, used to sort the components %%%%%%%%%%%
      %
      if strcmpi(g.sortvar,'mp')  % Maximum Power of backproj
          sortvar(c) = maxval;
          fprintf('IC%1.0f maximum mean power of back-projection: %g\n',c,sortvar(c));

      elseif strcmpi(g.sortvar, 'pv')   % Percent Variance
          % toby 2.19.2006: Changed to be consistent with eeg_pvaf().
          sortvar(c) = 100-100*mean(var(data(g.plotchans,limframe1:limframe2)...
                        - proj(g.plotchans,limframe1:limframe2),1))/vardat;
          fprintf('IC%1.0f percent variance accounted for(pvaf): %2.2f%%\n',c,sortvar(c));

      elseif strcmpi(g.sortvar,'pp')    % Percent Power
          sortvar(c) = 100-100*mean(mean((data(g.plotchans,limframe1:limframe2)...
                        - proj(g.plotchans,limframe1:limframe2)).^2))/powdat;
          fprintf('IC%1.0f percent power accounted for(ppaf): %2.2f%%\n',c,sortvar(c));

      elseif strcmpi(g.sortvar,'rp')    % Relative Power
          sortvar(c) = 100*mean(mean((proj(g.plotchans,limframe1:limframe2)).^2))/powdat;
          fprintf('IC%1.0f relative power of back-projection: %g\n',c,sortvar(c));
      else
          error('''sortvar'' argument unknown');
      end

end % component c
fprintf('\n');

%
%%%%%%%%%%%%%%% Compute component selection criterion %%%%%%%%%%%%%%%%%%%%%%%%%%
%

% compute sortvar
if ~xunitframes
  fprintf('  in the interval %3.0f ms to %3.0f ms.\n',...
					1000*times(limframe1),1000*times(limframe2));
end

[sortsortvar spx] = sort(sortvar);
sortsortvar = sortsortvar(end:-1:1);
spx      = spx(end:-1:1);
npercol = ceil(ncomps/3);

%
%%%%%%%%%%%%%%%%%%%%%%%%% Sort the components %%%%%%%%%%%%%%%%%%%%%%%%%%%
%

[tmp,compx]  = sort(sortvar');           % sort compnums on sortvar (as defined by input
                                      % 'sortvar', default is 'mp').
compx        = compx(ncomps:-1:1);    % reverse order of sort
compvars     = compvars(ncomps:-1:1)';% reverse order of sort (output var)
compvarorder = g.compnums(compx);     % actual component numbers (output var)
plotframes   = plotframes(compx);     % plotted comps have these max frames
compframes   = plotframes';           % frame of max variance in each comp (output var)
comptimes    = times(plotframes(compx)); % time of max variance in each comp (output var)
compsplotted = compvarorder(1:ntopos); % (output var)

%
%%%%%%%%%%%%%%%%%%%%%%%% Reduce to ntopos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

[plotframes,ifx] = sort(plotframes(1:ntopos));% sort plotframes on their temporal order
plottimes  = times(plotframes);       % convert to times in ms
compx      = compx(ifx);              % indices into compnums, in plotting order
maporder   = g.compnums(compx);       % comp. numbers, in plotting order (l->r)
maxproj    = maxproj(:,compx);        % maps in plotting order

vlen = length(g.voffsets); % extend voffsets if necessary
while vlen< ntopos
  g.voffsets = [g.voffsets g.voffsets(vlen)]; % repeat last offset given
  vlen=vlen+1;
end

head_sep = 1.2;
topowidth = pos(3)/(ntopos+(ntopos-1)/5); % width of each topoplot
if topowidth > 0.20    % adjust for maximum height
    topowidth = 0.2;
end

if rem(ntopos,2) == 1  % odd number of topos
   topoleft = pos(3)/2 - (floor(ntopos/2)*head_sep + 0.5)*topowidth;
else % even number of topos
   topoleft = pos(3)/2 - (floor(ntopos/2)*head_sep)*topowidth;
end

%
%%%%%%%%%%%%%%%%%%%% Print times and frames of comp maxes %%%%%%%%%%%%%%
%

% fprintf('\n');
fprintf('Plotting envelopes of %d component projections.\n',ntopos);
if length(g.plotchans) ~= chans
  fprintf('Envelopes computed from %d specified data channels.\n',...
      length(g.plotchans));
end

fprintf('Topo maps will show components: ');
for t=1:ntopos
  fprintf('%4d  ',maporder(t));
end

fprintf('\n');
if ~xunitframes
  fprintf('    with max var at times (ms): ');
  for t=1:ntopos
    fprintf('%4.0f  ',1000*plottimes(t));
  end
  fprintf('\n');
end

fprintf('                  epoch frames: ');
for t=1:ntopos
   fprintf('%4d  ',limframe1-1+plotframes(t));
end
fprintf('\n');

fprintf('    Component sortvar in interval:  ');
for t=1:ntopos
     fprintf('%4.2f ',sortvar(t));
end
fprintf('\n');

sumproj = zeros(size(data(g.plotchans,:))); % toby 2.21.2006 REDUNDANT CALCULATIONS!
for n = 1:ntopos
  if isempty(g.icaact)
      sumproj = sumproj + ...
		g.icawinv(g.plotchans,maporder(n))*weights(maporder(n),:)*data;
  else
      sumproj = sumproj + g.icawinv(g.plotchans,maporder(n))*g.icaact(maporder(n),:);
       % updated -sm 11/04
  end; % Note: sumproj here has only g.plotchans
end
rmsproj = mean(mean((data(g.plotchans,limframe1:limframe2).^2))); % find data rms in interval

if strcmpi(g.sortvar,'rp')
 	sumppaf = mean(mean(sumproj(:,limframe1:limframe2).^2));
    sumppaf = 100*sumppaf/rmsproj;
    ot   = 'rp';
elseif strcmpi(g.sortvar,'pv')
    sumppaf = mean(var((data(g.plotchans,limframe1:limframe2) ...
                                  - sumproj(:,limframe1:limframe2)),1));
    sumppaf = 100-100*sumppaf/vardat;
    ot   = 'pvaf';
else
   sumppaf = mean(mean((data(g.plotchans,limframe1:limframe2) ...
                                  - sumproj(:,limframe1:limframe2)).^2));
    sumppaf = 100-100*sumppaf/rmsproj;
    ot   = 'ppaf';
end

if ~xunitframes
   fprintf('    Summed component ''%s'' in interval [%4g %4g] ms: %4.2f%%\n',...
					ot, 1000*times(limframe1),1000*times(limframe2), sumppaf);
end

%
% Collect user-supplied Y-axis information
% Edited and moved here from 'Collect Y-axis information' section below -Jean
%
if length(g.limits) == 4
     if g.limits(3)~=0 || g.limits(4)~=0 % collect plotting limits from 'limits'
	 ymin = g.limits(3);
	 ymax = g.limits(4);
         ylimset = 1;
     end
else
  ylimset = 0; % flag whether hard limits have been set by the user
  ymin = min(min(g.icawinv(g.plotchans,:)*g.icaact(:,pframes))); % begin by setting limits from plotted data
  ymax = max(max(g.icawinv(g.plotchans,:)*g.icaact(:,pframes)));
end

fprintf('    Plot limits (sec, sec, uV, uV) [%g,%g,%g,%g]\n\n',xmin,xmax, ymin,ymax);

%
%%%%%%%%%%%%%%%%%%%%% Plot the data envelopes %%%%%%%%%%%%%%%%%%%%%%%%%
%
BACKCOLOR = [0.7 0.7 0.7];
FONTSIZE=12;
FONTSIZEMED=10;
FONTSIZESMALL=8;
newaxes=axes('position',pos);
axis off
set(newaxes,'FontSize',FONTSIZE,'FontWeight','Bold','Visible','off','Color',BACKCOLOR);
delete(newaxes) %XXX

% site the plot at bottom of the current axes
axe = axes('Position',[pos(1) pos(2) pos(3) g.axisoff*pos(4)],...
           'FontSize',FONTSIZE,'FontWeight','Bold');

g.limits = get(axe,'Ylim');
set(axe,'GridLineStyle',':')
set(axe,'Xgrid','off')
set(axe,'Ygrid','on')
axes(axe)
set(axe,'Color',axcolor);

%
%%%%%%%%%%%%%%%%% Plot the envelope of the summed selected components %%%%%%%%%%%%%%%%%
%
if BOLD_COLORS==1
  mapcolors = 1:ntopos+1;
else
  mapcolors = [1 maporder+1];
end

if strcmpi(g.sumenv,'on')  || strcmpi(g.sumenv,'fill') %%%%%%%% if 'sunvenv' %%%%%%%%%
 sumenv = envelope(sumproj(:,:), g.envmode);
 if ~ylimset && max(sumenv(:)) > ymax, ymax = max(sumenv(1,:)); end
 if ~ylimset && min(sumenv(:)) < ymin, ymin = min(sumenv(2,:)); end
 if strcmpi(g.sumenv,'fill')
    %
    % Plot the summed projection filled
    %
    mins = matsel(sumenv,frames,0,2,0);
    p=fill([times times(frames:-1:1)],...
        [matsel(sumenv,frames,0,1,0) mins(frames:-1:1)],g.fillcolor);
    set(p,'EdgeColor',g.fillcolor);
    hold on
    %
    % Overplot the data envelope so it is not covered by the fill()'d component
    %
    p=plot(times,matsel(envdata,frames,0,1,1),colors(mapcolors(1),1));% plot the max
    set(p,'LineWidth',2);                % component order (if BOLD_COLORS==0)
    p=plot(times,matsel(envdata,frames,0,2,1),colors(mapcolors(1),1));% plot the min
    set(p,'LineWidth',2);                % component order (if BOLD_COLORS==0)

 else % if no 'fill'
     if ntopos ==1 && g.plotproj
         hold on; pproj = plot(times,sumproj);
         set(pproj,'Tag','line_allprojections');
     end
    tmp = matsel(sumenv,frames,0,2,0);
    p=plot(times,tmp);% plot the min
    hold on
    set(p,'color',g.fillcolor,'Tag','line_envelope_2');
    set(p,'linewidth',3);
    p=plot(times,matsel(sumenv,frames,0,1,0));% plot the max
    set(p,'linewidth',3);
    set(p,'color',g.fillcolor,'Tag','line_envelope_2');
 end

 set(p,'Tag','patch_envelope');
end

if strcmpi(g.sortvar,'rp')
	t = text(double(xmin+0.1*(xmax-xmin)), ...
         double(ymin+0.1*(ymax-ymin)), ...
         ['rp ' num2str(sumppaf,'%4.2f') '%']);
elseif strcmpi(g.sortvar,'pv')
    	t = text(double(xmin+0.1*(xmax-xmin)), ...
         double(ymin+0.1*(ymax-ymin)), ...
         ['pvaf ' num2str(sumppaf,'%4.2f') '%']);
else
	t = text(double(xmin+0.1*(xmax-xmin)), ...
         double(ymin+0.1*(ymax-ymin)), ...
         ['ppaf ' num2str(sumppaf,'%4.2f') '%']);
end
set(t,'fontsize',FONTSIZESMALL,'fontweight','bold','Tag','text_inaxes')

%
% %%%%%%%%%%%%%%%%%%%%%%%% Plot the computed component envelopes %%%%%%%%%%%%%%%%%%
%
    envx = [1;compx+1];
    for c = 1:ntopos+1
        curenv = matsel(envdata,frames,0,1,envx(c));
        if ~ylimset && max(curenv) > ymax, ymax = max(curenv); end
        p=plot(times,curenv,colors(mapcolors(c),1),'Tag',['line_envelope_' num2str(c)]);% plot the max
        set(gca,'FontSize',FONTSIZESMALL,'FontWeight','Bold')
        if c==1                                % Note: use colors in original
            set(p,'LineWidth',2);              %       component order (if BOLD_COLORS==0)
        else
            set(p,'LineWidth',1);
        end
        if all_bold > 0
                set(p,'LineStyle','-','LineWidth',3);
        elseif mapcolors(c)>15                            % thin/dot 16th-> comp. envs.
                set(p,'LineStyle',':','LineWidth',1);
        elseif mapcolors(c)>10                            %
                set(p,'LineStyle',':','LineWidth',2);
        elseif mapcolors(c)>6                             % dot 6th-> comp. envs.
                set(p,'LineStyle',':','LineWidth',3);
        elseif mapcolors(c)>1
            set(p,'LineStyle',colors(mapcolors(c),2),'LineWidth',1);
            if colors(mapcolors(c),2) == ':'
                set(l1,'LineWidth',2);  % embolden dotted env lines
            end
        end
        hold on
        curenv = matsel(envdata,frames,0,2,envx(c));
        if ~ylimset && min(curenv) < ymin, ymin = min(curenv); end
        p=plot(times,curenv,colors(mapcolors(c),1),'Tag',['line_envelope_' num2str(c)]);% plot the min

        if c==1
            set(p,'LineWidth',2);
        else
            set(p,'LineWidth',1);
        end
        if all_bold > 0
                set(p,'LineStyle','-','LineWidth',3);
        elseif mapcolors(c)>15                            % thin/dot 11th-> comp. envs.
                set(p,'LineStyle',':','LineWidth',1);
        elseif mapcolors(c)>10
                set(p,'LineStyle',':','LineWidth',2);
        elseif mapcolors(c)>6                             % dot 6th-> comp. envs.
                set(p,'LineStyle',':','LineWidth',3);
        elseif mapcolors(c)>1
            set(p,'LineStyle',colors(mapcolors(c),2),'LineWidth',1);
            if colors(mapcolors(c),2) == ':'
                set(l1,'LineWidth',2);  % embolden dotted env lines
            end
        end
        if c==1 && ~isempty(g.vert)
            for v=1:length(g.vert)
                vl=plot([g.vert(v) g.vert(v)], [-1e10 1e10],'k--','Tag','line_vertline'); % plot specified vertical lines
                if any(g.limcontrib ~= 0) && v>= length(g.vert)-1;
                    set(vl,'linewidth',g.limcontribweight);
                    set(vl,'linestyle',':');
                else
                    set(vl,'linewidth',g.vertweight);
                    set(vl,'linestyle','--');
                end
            end
        end
        if g.limits(1) <= 0 && g.limits(2) >= 0    % plot vertical line at time zero
                vl=plot([0 0], [-1e10 1e10],'k');
                    set(vl,'linewidth',2,'Tag','line_vertatzero');
        end

        %
        % plot the n-th component filled
        %
        if g.fillcomp(1)>0 && find(g.fillcomp==c-1)
            fprintf('filling the envelope of component %d\n',c-1);
            mins = matsel(envdata,frames,0,2,envx(c));
            p=fill([times times(frames:-1:1)],...
                   [matsel(envdata,frames,0,1,envx(c)) mins(frames:-1:1)],...
                   colors(mapcolors(c),1));
            %
            % Overplot the data envlope again so it is not covered by the fill()'d component
            %
            p=plot(times,matsel(envdata,frames,0,1,envx(1)),colors(mapcolors(1),1));% plot the max
            set(p,'LineWidth',2);                % component order (if BOLD_COLORS==0)
            p=plot(times,matsel(envdata,frames,0,2,envx(1)),colors(mapcolors(1),1));% plot the min
            set(p,'LineWidth',2);                % component order (if BOLD_COLORS==0)
        end
    end  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%%%%%%%%%%%%%%%%%%%%%%% Extend y limits by 5% %%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~ylimset
  datarange = ymax-ymin;
  ymin = ymin-0.05*datarange;
  ymax = ymax+0.05*datarange;
end
axis([pmin pmax ymin ymax]);

%
%%%%%%%%%%%%%%%%%%%%%% Label axes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
set(axe,'Color',axcolor);
if ~xunitframes
   l= xlabel('Time (s)');
else % xunitframes == 1
   l= xlabel('Data (time points)');
end
set(l,'FontSize',FONTSIZEMED,'FontWeight','Bold');
if strcmpi(g.envmode, 'avg')
    l=ylabel('Potential (\muV)');
else
    l=ylabel('RMS of \muV');
end;
set(l,'FontSize',FONTSIZEMED,'FontWeight','Bold');
%
%%%%%%%%%%%%%% Draw maps and oblique/vertical lines %%%%%%%%%%%%%%%%%%%%%
%
% axall = axes('Units','Normalized','Position',pos,...
axall = axes('Position',pos,...
    'Visible','Off','Fontsize',FONTSIZE); % whole-figure invisible axes
axes(axall)
set(axall,'Color',axcolor);
axis([0 1 0 1])

width  = xmax-xmin;
pwidth  = pmax-pmin;
height = ymax-ymin;

if strcmpi(g.dispmaps, 'on')
    for t=1:ntopos % draw oblique lines from max env vals (or plot top)
                 % to map bases, in left to right order
        %
        %%%%%%%%%%%%%%%%%%% draw oblique lines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        if BOLD_COLORS==1
            linestyles = 1:ntopos;
        else
            linestyles = maporder;
        end
        axes(axall)
        axis([0 1 0 1]);
        set(axall,'Visible','off');
        maxenv = matsel(envdata,frames,plotframes(t),1,compx(t)+1);
        % max env val
        data_y = g.axisoff*(g.voffsets(t)+maxenv-ymin)/height;
        if (data_y > pos(2)+g.axisoff*pos(4))
            data_y = pos(2)+g.axisoff*pos(4);
        end
        l1 = plot([(plottimes(t)-pmin)/pwidth  ...
                   topoleft + 1/pos(3)*(t-1)*1.2*topowidth + (topowidth*g.axisoff)],...
                  [data_y (g.axisoff + 0.08)], ...
                  colors(linestyles(t)+1),'Tag',['line_topoconnect_' num2str(t)]); % 0.68 is bottom of topo maps
        if all_bold > 0
                set(l1,'LineStyle','-','LineWidth',3);
        elseif linestyles(t)>15                        % thin/dot 11th-> comp. envs.
                set(l1,'LineStyle',':','LineWidth',1);
        elseif linestyles(t)>10
                set(l1,'LineStyle',':','LineWidth',2);
        elseif linestyles(t)>5                     % dot 6th-> comp. envs.
                set(l1,'LineStyle',':','LineWidth',3);
        elseif linestyles(t)>1
            set(l1,'LineStyle',colors(linestyles(t)+1,2),'LineWidth',1);
            if colors(linestyles(t)+1,2) == ':'
                set(l1,'LineStyle',colors(linestyles(t)+1,2),'LineWidth',2);
            end
        end
        hold on
        %
        %%%%%%%%%%%%%%%%%%%% add specified vertical lines %%%%%%%%%%%%%%%%%%%%%%%%%
        %
        if g.voffsets(t) > 0
            l2 = plot([(plottimes(t)-xmin)/width  ...
                       (plottimes(t)-xmin)/width],...
                      [g.axisoff*(maxenv-ymin)/height ...
                       g.axisoff*(g.voffsets(t)+maxenv-ymin)/height],...
                      colors(linestyles(t)+1));
            if all_bold > 0
                    set(l2,'LineStyle','-','LineWidth',3);
            elseif linestyles(t)>15                      % thin/dot 11th-> comp. envs.
                    set(l2,'LineStyle',':','LineWidth',1);
            elseif linestyles(t)>10
                    set(l2,'LineStyle',':','LineWidth',2);
            elseif linestyles(t)>5                   % dot 6th-> comp. envs.
                    set(l2,'LineStyle',':','LineWidth',3);
            else
                set(l1,'LineStyle',colors(linestyles(t)+1,2),'LineWidth',1);
                if colors(linestyles(t)+1,2) == ':'
                    set(l1,'LineWidth',2);
                end
            end
        end
        set(gca,'Visible','off');
        axis([0 1 0 1]);
    end % t
end; % if g.dispmaps == on

%
%%%%%%%%%%%%%%%%%%%%%%%%% Plot the topoplots %%%%%%%%%%%%%%%%%%%%%%%%%%
%

if strcmpi(g.dispmaps, 'on')

    % common scale for colors
    % -----------------------
    if strcmpi(g.actscale, 'on')
        maxvolt = 0;
        for n=1:ntopos
            maxvolt = max(max(abs(maxproj(:,n))), maxvolt);
        end
    end

    [tmp tmpsort] = sort(maporder);
    [tmp tmpsort] = sort(tmpsort);

    for t=1:ntopos % left to right order  (maporder)
                   % axt = axes('Units','Normalized','Position',...
%         axt = axes('Units','Normalized','Position',...
%                    [pos(3)*topoleft+pos(1)+(t-1)*head_sep*topowidth pos(2)+0.66*pos(4) ...
%                     topowidth topowidth*head_sep]);
        axt = axes('Units','Normalized','Position',...
                   [pos(3)*topoleft+pos(1)+(t-1)*head_sep*topowidth  pos(2)+(g.axisoff+0.06)*pos(4) ...
                    topowidth topowidth*head_sep],...
                    'Tag',['axes_topo_' num2str(t)]);
        axes(axt)                             % topoplot axes
        cla

        if ~isempty(g.chanlocs)  % plot the component scalp maps
            if ~isempty(varargin)
                topoplot(maxproj(g.plotchans,t),g.chanlocs(g.plotchans), varargin{:});
            else  % if no varargin specified
                topoplot(maxproj(g.plotchans,t),g.chanlocs(g.plotchans),...
							'style','both','emarkersize',3);
            end
            axis square
            set(gca, 'userdata', ...
   ['text(-0.6, -0.6, ''' g.sortvar ': ' sprintf('%6.2f', sortvar(tmpsort(t))) ''');']);
        else
			axis off;
        end

        %
        %%%%%%%%%%%%% Scale colors %%%%%%%%%%%%%%%%%%%%%%%%%
        %
        if strcmpi(g.actscale, 'on')
            caxis([-maxvolt maxvolt]);
        end
        %
        %%%%%%%%%%%%%%%%%%%%%%%% label components %%%%%%%%%%%%%%%%%%%%%%%
        %
        if t==1
            chid = fopen('envtopo.labels','r');
            if chid <3,
                numlabels = 1;
            else
                fprintf('Will label scalp maps with labels from file %s\n',...
					'envtopo.labels');
                compnames = fscanf(chid,'%s',[4 MAXPLOTDATACHANS]);
                compnames = compnames';
                [r c] = size(compnames);
                for i=1:r
                    for j=1:c
                        if compnames(i,j)=='.',
                            compnames(i,j)=' ';
                        end
                    end
                end
                numlabels=0;
            end
        end
        if numlabels == 1
            complabel = int2str(maporder(t));        % label comp. numbers
        else
            complabel = compnames(t,:);              % use labels in file
        end
        text(0.00,0.80,['IC ' complabel],'FontSize',FONTSIZEMED,...
             'FontWeight','Bold','HorizontalAlignment','Center');
        % axt = axes('Units','Normalized','Position',[0 0 1 1],...
        axt = axes('Position',[0 0 1 1],...
                   'Visible','Off','Fontsize',FONTSIZE);
        set(axt,'Color',axcolor);           % topoplot axes
        drawnow
    end
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%% Plot a colorbar %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % axt = axes('Units','Normalized','Position',[.88 .58 .03 .10]);
%     axt = axes('Position',[pos(1)+pos(3)*1.015 pos(2)+0.6055*pos(4) ...
% 				pos(3)*.02 pos(4)*0.09]);
   axt = axes('Position',[pos(1)+pos(3)*1.015 pos(2)+(g.axisoff)*pos(4) ...
				pos(3)*.02 pos(4)*0.09]);
    if strcmpi(g.actscale, 'on')
        h=cbar(axt, [1:64],[-maxvolt maxvolt],3);
    else
        h=cbar(axt);                        % colorbar axes
        set(h,'Ytick',[]);

        axes(axall)
        set(axall,'Color',axcolor);
        tmp = text(0.50,1.05,g.title,'FontSize',FONTSIZE,...
						'HorizontalAlignment','Center',...
						'FontWeight','Bold');
        set(tmp, 'interpreter', 'none');
%         text(1,0.68,'+','FontSize',FONTSIZE,'HorizontalAlignment','Center');
        text(1,g.axisoff+0.08,'+','FontSize',FONTSIZE,'HorizontalAlignment','Center');
        % text(1,0.637,'0','FontSize',12,'HorizontalAlignment','Center',...
		%	'verticalalignment','middle');
%         text(1,0.61,'-','FontSize',FONTSIZE,'HorizontalAlignment','Center');
        text(1,g.axisoff+0.01,'-','FontSize',FONTSIZE,'HorizontalAlignment','Center');
    end
    axes(axall)
    set(axall,'layer','top'); % bring component lines to top

end
%
%%%%%%%%%%%%%%%%%%%%%%%%% turn on axcopy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
axcopy(gcf, ...
  'if ~isempty(get(gca,''''userdata'''')), eval(get(gca,''''userdata''''));end;');

return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function envdata = envelope(data, envmode)  % also in release as env()
  if nargin < 2
      envmode = 'avg';
  end
  if strcmpi(envmode, 'rms'); % return rms of pos and neg vals at each time point
      warning off;
      datnaeg = (data < 0).*data; % zero out positive values
      dataneg = -sqrt(sum(dataneg.^2,1) ./ sum(negflag,1));

      datapos = (data > 0).*data; % zero out negative values
      datapos = sqrt(sum(datapos.^2,1) ./ sum(posflag,1));

      envdata = [datapos;dataneg];
      warning on;
  else  % find max and min value at each time point
      if size(data,1)>1
          maxdata = max(data); % max at each time point
          mindata = min(data); % min at each time point
          envdata = [maxdata;mindata];
      else
          maxdata = max([data;data]); % max at each time point
          mindata = min([data;data]); % min at each time point
          envdata = [maxdata;mindata];
      end
  end

return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function Script_HOG_SVM_train_dyn()

% Change to your downloaded location
addpath('C:\liblinear\matlab')

%% load shared definitions and AU data
shared_defs;

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-9:0.5:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svm_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svm_test = @svm_test_linear;

pca_loc = '../../pca_generation/generic_face_rigid.mat';

%%
for a=1:numel(aus)

    au = aus(a);

    rest_aus = setdiff(all_aus, au);

    % load the training and testing data for the current fold
    [train_samples, train_labels, valid_samples, valid_labels, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic_dynamic(train_recs, devel_recs, au, rest_aus, SEMAINE_dir, hog_data_dir);

    train_samples = sparse(train_samples);
    valid_samples = sparse(valid_samples);

    %% Cross-validate here
    [ best_params, ~ ] = validate_grid_search_no_par(svm_train, svm_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);

    model = svm_train(train_labels, train_samples, best_params);

    [prediction, a, actual_vals] = predict(valid_labels, valid_samples, model);

    % Go from raw data to the prediction
    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;

    name = sprintf('models/AU_%d_dyn.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    write_lin_dyn_svm(name, means, svs, b, pos_lbl, neg_lbl);

    name = sprintf('results_SEMAINE_devel/AU_%d_dyn.mat', au);

    tp = sum(valid_labels == 1 & prediction == 1);
    fp = sum(valid_labels == 0 & prediction == 1);
    fn = sum(valid_labels == 1 & prediction == 0);
    tn = sum(valid_labels == 0 & prediction == 0);

    precision = tp/(tp+fp);
    recall = tp/(tp+fn);

    f1 = 2 * precision * recall / (precision + recall);

    save(name, 'model', 'f1', 'precision', 'recall', 'best_params', 'valid_labels', 'prediction');

end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
% -----------------------------------------------------------
% |Nagi Hatoum, Candlestick chart "ploter" function, 5/29/03 |
% |Edited and expanded by Jelle C.J. Veraa, 3/28/2009        |
% -----------------------------------------------------------
%
% Important: date plotting subroutine needs 'tlabel' by Carlos Adrian
%   Vargas Aguilera
%   Available at the MATLAB Central File Exchange site, File ID: #19314
%   http://www.mathworks.com/matlabcentral/fileexchange/19314
%
% function cndlV2(O,H,L,C)
%          cndlV2(O,H,L,C,date)
%          cndlV2(O,H,L,C,date,colorUp,colorDown,colorLine)
%          cndlV2(OHLC)
%          cndlV2(OHLC,date)
%          cndlV2(OHLC,date,colorUp,colorDown,colorLine)
%
% To use colors without a date, put '0' as date
%     i.e. cndlV2(OHLC,0,'b','r','k');
%
% required inputs: column vectors of O(pen), H(igh), L(ow)
% and C(lose) prices of commodity.
% Alternative: OHLC matrix of size [rowsx4]
%
% optional inputs [default]:
% - date: serial date number (make with 'datenum') [no dates, index#]
% - colorUp: Color for up candle                   ['w']
% - colorDown: Color for down candle               ['k']
% - colorLine: Color for lines                     ['k']
%
% Note: identical inputs as required for barChartPlot except for colors

function cndlV2(varargin)
% See if we have [OHLC] or seperate vectors and retrieve our
% required variables (Feel free to make this code more pretty ;-)
isMat = size(varargin{1},2);
indexShift = 0;
useDate = 0;

if isMat == 4,
    O = varargin{1}(:,1);
    H = varargin{1}(:,2);
    L = varargin{1}(:,3);
    C = varargin{1}(:,4);
else
    O = varargin{1};
    H = varargin{2};
    L = varargin{3};
    C = varargin{4};
    indexShift = 3;
end
if nargin+isMat < 7,
    colorDown = 'k';
    colorUp = 'w';
    colorLine = 'k';
else
    colorUp = varargin{3+indexShift};
    colorDown = varargin{4+indexShift};
    colorLine = varargin{5+indexShift};
end
if nargin+isMat < 6,
    date = (1:length(O))';
else
    if varargin{2+indexShift} ~= 0
        date = varargin{2+indexShift};
        useDate = 1;
    else
        date = (1:length(O))';
    end
end

% w = Width of body, change multiplier to draw body thicker or thinner
% the 'min' ensures no errors on weekends ('time gap Fri. Mon.' > wanted
% spacing)
w=.3*min([(date(2)-date(1)) (date(3)-date(2))]);
%%%%%%%%%%%Find up and down days%%%%%%%%%%%%%%%%%%%
d=C-O;
l=length(d);
hold on
%%%%%%%%draw line from Low to High%%%%%%%%%%%%%%%%%
for i=1:l
   line([date(i) date(i)],[L(i) H(i)],'Color',colorLine)
end
%%%%%%%%%%draw white (or user defined) body (down day)%%%%%%%%%%%%%%%%%
n=find(d<0);
for i=1:length(n)
    x=[date(n(i))-w date(n(i))-w date(n(i))+w date(n(i))+w date(n(i))-w];
    y=[O(n(i)) C(n(i)) C(n(i)) O(n(i)) O(n(i))];
    fill(x,y,colorDown)
end
%%%%%%%%%%draw black (or user defined) body(up day)%%%%%%%%%%%%%%%%%%%
n=find(d>=0);
for i=1:length(n)
    x=[date(n(i))-w date(n(i))-w date(n(i))+w date(n(i))+w date(n(i))-w];
    y=[O(n(i)) C(n(i)) C(n(i)) O(n(i)) O(n(i))];
    fill(x,y,colorUp)
end

if (nargin+isMat > 5) && useDate,
    tlabel('x');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold off

author: Sable

repo: mcbench-benchmarks


text:
function pp = sphere_imp_point_project_3d ( r, center, p )

%*****************************************************************************80
%
%% SPHERE_IMP_POINT_PROJECT_3D projects a point onto an implicit sphere in 3D.
%
%  Discussion:
%
%    An implicit sphere in 3D satisfies the equation:
%
%      sum ( ( P(1:DIM_NUM) - CENTER(1:DIM_NUM) )**2 ) = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real CENTER(3), the center of the sphere.
%
%    Input, real P(3), a point.
%
%    Output, real PP(3), the projected point.
%
  dim_num = 3;

  if ( r == 0.0 )

    pp(1:dim_num) = center(1:dim_num);

  elseif ( p(1:dim_num) == center(1:dim_num) )

    pp(1:dim_num) = center(1:dim_num) + r / sqrt ( dim_num );

  else

    norm = sqrt ( sum ( ( p(1:dim_num) - center(1:dim_num) ).^2 ) );

    pp(1:dim_num) = center(1:dim_num) + r * ( p(1:dim_num) - center(1:dim_num) ) / norm;

  end

  return
end


text:
function [Y] = LaplacianScore(X, W)
%	Usage:
%	[Y] = LaplacianScore(X, W)
%
%	X: Rows of vectors of data points
%	W: The affinity matrix.
%	Y: Vector of (1-LaplacianScore) for each feature.
%      The features with larger y are more important.
%
%    Examples:
%
%       fea = rand(50,70);
%       options = [];
%       options.Metric = 'Cosine';
%       options.NeighborMode = 'KNN';
%       options.k = 5;
%       options.WeightMode = 'Cosine';
%       W = constructW(fea,options);
%
%       LaplacianScore = LaplacianScore(fea,W);
%       [junk, index] = sort(-LaplacianScore);
%
%       newfea = fea(:,index);
%       %the features in newfea will be sorted based on their importance.
%
%	Type "LaplacianScore" for a self-demo.
%
% See also constructW
%
%Reference:
%
%   Xiaofei He, Deng Cai and Partha Niyogi, "Laplacian Score for Feature Selection".
%   Advances in Neural Information Processing Systems 18 (NIPS 2005),
%   Vancouver, Canada, 2005.
%
%   Deng Cai, 2004/08

if nargin == 0, selfdemo; return; end

[nSmp,nFea] = size(X);

if size(W,1) ~= nSmp
    error('W is error');
end

D = full(sum(W,2));
L = W;

allone = ones(nSmp,1);

tmp1 = D'*X;

D = sparse(1:nSmp,1:nSmp,D,nSmp,nSmp);

DPrime = sum((X'*D)'.*X)-tmp1.*tmp1/sum(diag(D));
LPrime = sum((X'*L)'.*X)-tmp1.*tmp1/sum(diag(D));

DPrime(find(DPrime < 1e-12)) = 10000;

Y = LPrime./DPrime;
Y = Y';
Y = full(Y);

%---------------------------------------------------
function selfdemo
% ====== Self demo using IRIS dataset
% ====== 1. Plot IRIS data after LDA for dimension reduction to 2D
load iris.dat

feaNorm = mynorm(iris(:,1:4),2);
fea = iris(:,1:4) ./ repmat(max(1e-10,feaNorm),1,4);

options = [];
options.Metric = 'Cosine';
options.NeighborMode = 'KNN';
options.WeightMode = 'Cosine';
options.k = 3;

W = constructW(fea,options);

[LaplacianScore] = feval(mfilename,iris(:,1:4),W);
[junk, index] = sort(-LaplacianScore);

index1 = find(iris(:,5)==1);
index2 = find(iris(:,5)==2);
index3 = find(iris(:,5)==3);
figure;
plot(iris(index1, index(1)), iris(index1, index(2)), '*', ...
     iris(index2, index(1)), iris(index2, index(2)), 'o', ...
     iris(index3, index(1)), iris(index3, index(2)), 'x');
legend('Class 1', 'Class 2', 'Class 3');
title('IRIS data onto the first and second feature (Laplacian Score)');
axis equal; axis tight;

figure;
plot(iris(index1, index(3)), iris(index1, index(4)), '*', ...
     iris(index2, index(3)), iris(index2, index(4)), 'o', ...
     iris(index3, index(3)), iris(index3, index(4)), 'x');
legend('Class 1', 'Class 2', 'Class 3');
title('IRIS data onto the third and fourth feature (Laplacian Score)');
axis equal; axis tight;

disp('Laplacian Score:');
for i = 1:length(LaplacianScore)
    disp(num2str(LaplacianScore(i)));
end

author: ZJULearning

repo: MatlabFunc


text:
function [binId, cumArea] = hist(grains,varargin)
% plot a grain size histogram
%
% Syntax
%   hist(grains)
%   hist(grains,n) % specify the number of bins
%
%   % use an arbitrary property for the histogramm
%   hist(grains,grains.equivalentRadius)
%
%   binId = hist(grains) % returns the bin for each grain
%
% Input
%   grains - @grain2d
%   n      - number of bins, default ist 15
%

[mtexFig,isNew] = newMtexFigure(varargin{:});
mtexFig.keepAspectRatio = false;

% exract area and maybe an aditional property
area = grains.area;
if nargin>1 && isnumeric(varargin{1}) && length(varargin{1}) == length(grains)
  prop = varargin{1};
  varargin(1) = [];
else
  prop = area;
end

if ~isempty(varargin) && isnumeric(varargin{1})
  nbins = varargin{1};
else
  nbins = 15;
end

idList = grains.indexedPhasesId;

bins = linspace(0,max(prop)+eps,nbins);

cumArea = zeros(numel(bins)-1,numel(idList));
binId = zeros(length(grains),1);

% loop through all phases
for id = 1:numel(idList)

  phaseId = grains.phaseId==idList(id);
  phaseArea = area(phaseId);
  phaseProp = prop(phaseId);

  % find for each area the binId
  [~,~,binIdLocal] = histcounts(phaseProp,bins);
  binId(phaseId) = binIdLocal;

  % compute the sum of areas belonging to the same bin
  for i = 1:numel(bins)-1
    cumArea(i,id) = sum(phaseArea(binIdLocal == i)) ./ sum(area);
  end

end

% plot the result as a bar plot
binCenter = 0.5*(bins(1:end-1)+bins(2:end));
binWidth = 1 + 0.5*size(cumArea,2)>1;
bar(binCenter,100*cumArea,'BarWidth',binWidth,'parent',mtexFig.gca)
xlim(mtexFig.gca,[bins(1),bins(end)])
if all(prop == area)
  title(mtexFig.gca,'grain size distribution')
  xlabel(mtexFig.gca,'grain area');
end
ylabel(mtexFig.gca,'relative area (%)')

min = grains.mineralList(idList);
legend(min{:})

if isNew, mtexFig.drawNow(varargin{:});end

if nargout == 0, clear binId; end

author: mtex-toolbox

repo: mtex

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


text:
%%*****************************************************************
%% HSDNTdirfun: compute (dX,dZ), given dy, for the NT direction.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

    function [par,dX,dy,dZ] = HSDNTdirfun(blk,At,par,Rd,EinvRc,xx);

    global solve_ok

    dX = cell(size(blk,1),1); dZ = cell(size(blk,1),1); dy = [];
    if (any(isnan(xx)) | any(isinf(xx)))
       solve_ok = 0;
       fprintf('\n  HSDNTdirfun: solution contains NaN or inf.');
       return;
    end
%%
    m = par.m;
    dy2 = xx(1:m+2);
%%
    for p=1:size(blk,1)
       pblk = blk(p,:);
       if strcmp(pblk{1},'l')
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),[],[],dy2));
          tmp   = par.dd{p}.*dZ{p};
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'q')
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),[],[],dy2));
          tmp = par.dd{p}.*dZ{p} + qops(pblk,qops(pblk,dZ{p},par.ee{p},1),par.ee{p},3);
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'s')
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),par.permA(p,:),par.isspAy(p),dy2));
          tmp   = Prod3(pblk,par.W{p},dZ{p},par.W{p},1);
          dX{p} = EinvRc{p}-tmp;
       end
    end
    dy  = dy2(1:m);
    par.dtau = dy2(m+1);
    par.dtheta = dy2(m+2);
    par.dkap = (par.mu./par.tau - par.kap) - par.kap*(par.dtau/par.tau);
%%*******************************************************************

author: yu-jiang

repo: radpbook

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


text:
function [flowx,flowy] = specific_wind(x,y,nel)
%constant_wind   Reference problem 3.3 convective wind
%   [flowx,flowy] = specific_wind(x,y,nel);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          nel        number of elements
%   specifies constant wind at angle theta to vertical
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      theta= -pi/6;
      flowx =  sin(theta)*ones(nel,1);
      flowy =  cos(theta)*ones(nel,1);
      return


text:
function anatIm = recomputeAnatImage(vw,displayMode,slice)
% anatIm = recomputeAnatImage(vw,[displayMode],[slice]);
%
% Return the anat image, clipped and or scaled
% according to the (non-hidden) vw's ui settings,
% for the specified display mode. displayMode
% defaults to the currently-selected mode (e.g.,
% 'anat' 'amp' 'map' etc). slice defaults to the
% current slice.
%
%
%
% ras 01/05. Added b/c I want to migrate to using
% brightness/contrast/possibly gamma, but also
% want back compatibility.
if isequal(viewGet(vw, 'Name'),'hidden')
    % hidden views: just return the anat
    % image w/o contrast adjustment
    %TODO: put slice error handling here
    anatIm = viewGet(view, 'anatomycurrentslice', slice);
    return
end

ui = viewGet(vw,'ui');

if ieNotDefined('slice')
	% Get curSlice from ui
	slice = viewGet(vw, 'curSlice');
end

if ieNotDefined('displayMode')
    displayMode = ui.displayMode;
end

modeStr = sprintf('%sMode',displayMode);
numGrays = ui.(modeStr).numGrays;

% Get anatomy image from vw (non-scaled)
anatIm = cropCurAnatSlice(vw,slice);

if isfield(ui,'brightness')
    % adjust img brightness/contrast
    brightness = get(ui.brightness.sliderHandle,'Value');
    contrast = get(ui.contrast.sliderHandle,'Value');

    % unlike the normal way contrast/brightness work,
    % I've found it's better to have 'contrast'
    % just change the upper bound of the anatClip,
    % and 'brightness' just shift the median value
    % up and down a bit:
    minVal = double(min(anatIm(:)));
	maxVal = (1-contrast)*double(max(anatIm(:)));
    % removed lines that turned off/on warnings...
	anatIm = (rescale2(double(anatIm),[minVal maxVal],[1 numGrays]));

    % brighten
    brightDelta = brightness - 0.5;
    if brightDelta ~= 0 % slowwww....
        anatIm = brighten(anatIm,brightDelta);
        anatIm = rescale2(anatIm,[],[1 numGrays]);
    end
else
    % do it the old way, using the
    % anat clip slider values

	% Get anatClip from sliders
	anatClip = getAnatClip(vw);

	% Rescale anatIm to [1:numGrays], anatClip determines the range
	% of anatomy values that gets mapped to the available grayscales.
	% If anatClip=[0,1] then there is no clipping and the entire
	% range of anatomy values is scaled to the range of available gray
	% scales.
    minVal = double(min(anatIm(:)));
	maxVal = double(max(anatIm(:)));
	anatClipMin = min(anatClip)*(maxVal-minVal) + minVal;
	anatClipMax = max(anatClip)*(maxVal-minVal) + minVal;
	warning off;
	anatIm = (rescale2(double(anatIm),[anatClipMin,anatClipMax],[1,numGrays]));
	warning backtrace;
end

return

% This way seemed sensible, but
% the images didn't look as great as
% I'd hoped:
%     % Will rescale the values in anatIm to be
%     % distributed with a median determined by
%     % the brightness param, and a range of values
%     % determined by contrast -- but it will all
%     % fit in the range 1:numGrays:
%     mu = round(brightness*numGrays);
%     sigma = (1-contrast) * (numGrays/2-1) + 1;
%     anatIm = rescale2(anatIm,[],[-sigma sigma]) + mu;
%
%     % need a 2nd clip step -- could set this in the
%     % rescale2 step, but it's the same # of cycles to
%     % do it this way, I think, and more legible:
%     anatIm(anatIm <= 1) = 1;
%     anatIm(anatIm >= numGrays) = numGrays;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function mean = normal_ms_mean ( mu, sigma )

%*****************************************************************************80
%
%% NORMAL_MS_MEAN returns the mean of the Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, SIGMA, the parameters of the PDF.
%    0.0 < SIGMA.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = mu;

  return
end


text:
function answer = is_function_called (line, func)
    % number of occurrences of each pattern:
    n1 = length( regexp(line, ['\<' func '\>']) );
    n2 = length( regexp(line, ['''.*\<' func '\>.*''']) );

    answer = n1 > n2;
end

%!test
%! myassert(is_function_called ('test', 'test'));
%! myassert(is_function_called ('abc test', 'test'));
%! myassert(is_function_called ('test abc', 'test'));
%! myassert(is_function_called ('abc test abc', 'test'));
%! myassert(~is_function_called ('abctest', 'test'));
%! myassert(~is_function_called ('testabc', 'test'));
%! myassert(~is_function_called ('abc ''test'' abc', 'test'));
%! myassert(~is_function_called ('''abc test abc''', 'test'));
%! myassert(is_function_called ('''abc'' test abc', 'test'));
%! myassert(is_function_called ('abc test ''abc''', 'test'));

%!test
%! myassert(is_function_called(...
%!     'myassert (s.identifier, ''myassert:error'');', 'myassert'));