text:
function prob = sedumi2opti(prob)
%SEDUMI2OPTI  Converts a SEDUMI Problem (At, b, C, K) to OPTI Format
% prob = sedumi2opti(prob)

%   Copyright (C) 2013 Jonathan Currie (I2C2)

if(~isfield(prob,'sdcone') || isempty(prob.sdcone) || ~isstruct(prob.sdcone))
    error('This function expects a prob.sdcone to contain a structure');
end
sdcone = prob.sdcone;
if((~isfield(sdcone,'A') && ~isfield(sdcone,'At')) || ~isfield(sdcone,'b') || (~isfield(sdcone,'C') && ~isfield(sdcone,'c')) || ~isfield(sdcone,'K'))
    error('This function expects prob.sdcone to contain the fields At (or A), b, C and K');
end
if(~isstruct(sdcone.K))
    error('SeDuMi prob.sdcone.K must be a structure!');
end
if(~isfield(sdcone.K,'s') && ~isfield(sdcone.K,'l'))
    error('SeDuMi prob.sdcone.K must contain the fields l or s (or both)');
end
if(isfield(sdcone.K,'f')), error('SeDuMi free variables are not supported in this interface - use SeDuMi directly to solve!'); end
if(isfield(sdcone.K,'q')), error('SeDuMi SOCP constraints are not supported in this interface - use SeDuMi directly to solve!'); end
if(isfield(sdcone.K,'r')), error('SeDuMi rotated cone constraints are not supported in this interface - use SeDuMi directly to solve!'); end
if(isfield(sdcone.K,'scomplex') || isfield(sdcone.K,'ycomplex') || isfield(sdcone.K,'xcomplex'))
    error('Complex variables are not supported in this interface - use SeDuMi to solve!');
end

%Read SeDuMi problem, convert to OPTI standard primal form
% MINIMIZE   dot(f,x) SUCH THAT  SUM A_i*x_i - C >= 0.
top = 1; dims = 0;
haveSDP = false; haveLP = false;
K = sdcone.K;
if(isfield(sdcone,'C'))
    c = sdcone.C;
else
    c = sdcone.c;
end
if(isfield(sdcone,'A'))
    A = sdcone.A;
else
    A = sdcone.At;
end
b = full(sdcone.b);

%Check c, b is a vector
if(size(c,2) > 1 && size(c,1) > 1)
    error('SeDuMi c field must be a column vector!');
end
if(size(b,2) > 1 && size(b,1) > 1)
    error('SeDuMi b field must be a column vector!');
end

%Flip A based on b dimensions
if(size(A,2) ~= length(sdcone.b))
    if(size(A,1) == length(sdcone.b))
        A = A';
    else
        error('SeDuMi A (At) appears to  have incorrect dimensions.\nExpected one dimension of A to be the same as length as %s.','b');
    end
end
if(size(c,2) > size(c,1))
    c = c';
end
%Check length of C against A
if(length(c) ~= size(A,1))
    error('SeDuMi c does not have the same number of elements as rows in A (At)');
end
%Check dims entered in K
if(isfield(K,'s') && ~isempty(K.s) && K.s(1) > 0)
    haveSDP = true;
    dims = dims + sum(K.s.^2);
end
if(isfield(K,'l') && ~isempty(K.l) && K.l(1) > 0)
    if(length(K.l) > 1)
        error('The SeDuMi K.l field must contain a scalar, the number of linear constraints');
    end
    haveLP = true;
    dims = dims + K.l;
end
%Check dims against rows
if(dims ~= size(A,1))
    error('SeDuMi A (At) does not have correct dimensions as per the specifications in K. Expected %d rows.',dims);
end

%Extract Objective
prob.f = full(-b);

%Extract Linear Constraints
if(haveLP)
    n = K.l;
    if(isfield(prob,'A') && ~isempty(prob.A)) %concatenate
        prob.A = [prob.A; A(top:top+n-1,:)];
        prob.b = full([prob.b; prob.c(top:top+n-1,:)]);
    else
        prob.A = A(top:top+n-1,:);
        prob.b = full(c(top:top+n-1,:));
    end
    top = top+n;
end
%Extract Semidefinite Constraints
if (haveSDP)
    prob.sdcone = cell(length(K.s),1);
    for i = 1:length(K.s)
        n = K.s(i);
        prob.sdcone{i} = [-c(top:top+n^2-1,:) -A(top:top+n^2-1,:)];
        top = top+n^2;
    end
else
    prob.sdcone = [];
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%HOMLINE Homogeneous line from two points
%
% L = HOMLINE(X1, Y1, X2, Y2) is a vector (3x1) which describes a line in
% homogeneous form that contains the two Euclidean points (X1,Y1) and (X2,Y2).
%
% Homogeneous points X (3x1) on the line must satisfy L'*X = 0.
%
% See also PLOT_HOMLINE.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% TODO, probably should be part of a HomLine class.

function l = homline(x1, y1, x2, y2)

    l = cross([x1 y1 1], [x2 y2 1]);

    % normalize so that the result of x*l' is the pixel distance
    % from the line
    l = l / norm(l(1:2));

author: Allopart

repo: rbpf-gmapping


text:
% WAVELET_FACTORY_1D Create wavelet cascade
%
% Usage
%    [Wop, filters] = WAVELET_FACTORY_1D(N)
%
%    [Wop, filters] = WAVELET_FACTORY_1D(N, filt_opt)
%
%    [Wop, filters] = WAVELET_FACTORY_1D(N, filt_opt, scat_opt)
%
% Input
%    N (int): The size of the signals to be transformed.
%    filt_opt (struct): The filter options, same as FILTER_BANK.
%    scat_opt (struct): The scattering options, same as WAVELET_LAYER_1D.
%
% Output
%    Wop (cell): A cell array of wavelet layer transforms needed for the
%       scattering transform.
%    filters (cell): A cell array of the filters used in defining the
%       wavelets.
%
% Description
%    In order to calculate the scattering coefficients of a signal, a set of
%    filter banks need to be defined first. Given these filter banks, wavelet
%    transforms can be defined on scattering layers. To obtain the former,
%    WAVELET_FACTORY_1D calls the FILTER_BANK function using the parameters
%    provided in filt_opt, the result of which is returned as filters.
%    Each filter bank is then used to create a layer operator, using the
%    function WAVELET_LAYER_1D. This function takes a layer, a filter bank
%    and a set of parameters as input. WAVELET_FACTORY_1D fixes the filter
%    bank as to the element of filters corresponding to the layer order, while
%    scat_opt is used as the parameters. The result is a cell array Wop,
%    of layer operators which take one layer as input, and return two layers
%    as output, A and V corresponding to the "average" and "variation" of the
%    input layer.
%
% See also
%    WAVELET_LAYER_1D, WAVELET_1D, FILTER_BANK

function [Wop, filters] = wavelet_factory_1dave(N, filt_opt, scat_opt)
	if nargin < 2
		filters = filter_bank(N);
    else
        filters = filter_bank(N, filt_opt);
    end

	if nargin < 3
		scat_opt.M = 2; % M is the scattering order
    else
        scat_opt = fill_struct(scat_opt, 'M', 2);
    end

    Wop = cell(1,scat_opt.M);
	for m = 0:scat_opt.M
		filt_ind = min(numel(filters), m+1);
		Wop{m+1} = @(X)(wavelet_layer_1dave(X, filters{filt_ind}, scat_opt));
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function freqplot(axs_hand,mode)
% FREQPLOT Frequency response plotting.
%          FREQPLOT handles the plotting of all graphs involved with
%          frequency.

% Author: Craig Borghesani
% Date: 8/8/94
% Revised: 10/18/94, 11/16/94, 12/10/94
% Copyright (c) 1999, Prentice-Hall

% mode represents all the options for frequency response.
% mode = 13, switching systems
% mode = 14, switching pages
% mode = 15, changing layouts

% if no frequency plots on current page, dont do anything
if ~length(axs_hand), return; end

% obtain handle information
f = gcf;
ui_data = get(f,'userdata');
ui_han = ui_data{1};

% current plant response
cur_sys   = get(ui_han(30),'userdata');
oth_sys   = [1:(cur_sys-1),(cur_sys+1):3];
plant_mat = get(ui_han(3+cur_sys),'userdata');
for_mat   = get(ui_han(6+cur_sys),'userdata');
bac_mat   = get(ui_han(9+cur_sys),'userdata');
sys_res   = get(ui_han(cur_sys),'userdata');
font_size = get(ui_han(42),'userdata');
stat_bar  = get(ui_han(43),'userdata');
cur_axs   = get(ui_han(32),'userdata');
w         = sys_res(1,:);
sys_cplx  = sys_res(2,:);

sys_axes  = ui_han(12);
show_all  = ui_han(25);
open_loop = ui_han(60);
clos_loop = ui_han(61);
exact     = ui_han(62);
straight  = ui_han(63);
margins   = ui_han(64);
range1    = 60:61;
range2    = 62:64;
fontname = 'times';

% extract data stored in axis handles
axs_data = [];
for k = 1:length(axs_hand),
 axs_data = [axs_data;get(axs_hand(k),'userdata')];
end
nyq = find(axs_data(:,1)==1);
bodm = find(axs_data(:,1)==2);
bodp = find(axs_data(:,1)==3);
nic = find(axs_data(:,1)==4);

if length(mode) & any(mode == [1,2]),

 if mode == 1, % open loop
  if strcmp(get(open_loop,'checked'),'on'),
   set(open_loop,'checked','off');
  else
   set(open_loop,'checked','on');
  end
 end

 if mode == 2, % closed loop
  if strcmp(get(clos_loop,'checked'),'on'),
   set(clos_loop,'checked','off');
  else
   set(clos_loop,'checked','on');
  end
 end

 ct = 2;
 for k = range2,
  if strcmp(get(ui_han(k),'checked'),'on') & ...
     strcmp(get(ui_han(k),'enable'),'on'),
   mode(ct) = k-(range1(1)-1);
  else
   mode(ct) = 0;
  end
  ct = ct + 1;
 end

elseif length(mode) & any(mode == [3,4,5]),

 ct = 2;
 for k = range1,
  if strcmp(get(ui_han(k),'checked'),'on') & ...
     strcmp(get(ui_han(k),'enable'),'on'),
   mode(ct) = k-(range1(1)-1);
  else
   mode(ct) = 0;
  end
  ct = ct + 1;
 end

end

if length(mode) < 2,
 if ~length(mode) | any(mode == [12:15]),

% determine various states of frequency response environments
  if ~length(mode), mode = 0; end
  ct = 2;
  for k = [range1,range2],
   if strcmp(get(ui_han(k),'checked'),'on') & ...
      strcmp(get(ui_han(k),'enable'),'on'),
    mode(ct) = k-(range1(1)-1);
   else
    mode(ct) = 0;
   end
   ct = ct + 1;
  end

 end
end

open_flag = strcmp(get(open_loop,'checked'),'on');
clos_flag = strcmp(get(clos_loop,'checked'),'on');
show_flag = strcmp(get(show_all,'label'),'Show All Systems');
hide_flag = strcmp(get(show_all,'label'),'Hide Other Systems');

if clos_flag,
 bac_cplx = termcplx(bac_mat,w);
 tmp_cplx = sys_cplx./bac_cplx;
 clos_cplx = tmp_cplx./(1+sys_cplx);
end

mode_len = 3;
if length(mode) == 3 & hide_flag,
 mode = [mode,14];
 mode_len = 4;
end

if any(mode==3), % exact plot

 if strcmp(get(exact,'checked'),'off') | length(mode)>mode_len,
  set(exact,'checked','on');

  if length(nyq), % nyquist plot

% update current system frequency response
   if all(mode ~= 12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(nyq,1+oth_sys),'vis','off');
     set(axs_data(nyq,4+oth_sys),'vis','off');
    elseif any(mode == 13),
     if open_flag,
      set(axs_data(nyq,1+oth_sys),'color','r');
      set(axs_data(nyq,4+oth_sys),'color','b');
     end
%     if clos_flag,
%      set(axs_data(nyq,4+oth_sys),'color','b');
%     end
    end

    if open_flag & any(mode == 1),
     set(stat_bar,'string','Computing exact open-loop frequency response');
     rl = real(sys_cplx);  im = imag(sys_cplx);
     set(axs_data(nyq,1+cur_sys),'xdata',rl,'ydata',im,...
           'color','g','vis','on');
     set(axs_data(nyq,4+cur_sys),'xdata',rl,'ydata',-im,...
           'color','m','vis','on');
    elseif any(mode == 1),
     set(stat_bar,'string','Removing exact open-loop frequency response');
     set(axs_data(nyq,1+cur_sys),'vis','off');
     set(axs_data(nyq,4+cur_sys),'vis','off');
    end

%    if clos_flag,
%     ers = 'norm';
%     if ~open_flag, ers = 'xor'; end
%     rl = real(clos_cplx);  im = imag(clos_cplx);
%     set(axs_data(nyq,4+cur_sys),'xdata',rl,'ydata',im,...
%           'color','m','vis','on');
%    elseif any(mode == 2),
%     set(axs_data(nyq,4+cur_sys),'vis','off');
%    end

   end

% switching pages or requesting to veiw all systems
   if (any(mode == 12) & show_flag) | (any(mode == 14) & hide_flag) | ...
      (any(mode == 15) & hide_flag),

    for k = 1:2,
     oth_res = get(ui_han(oth_sys(k)),'userdata');
     oth_cplx = oth_res(2,:);

%     if clos_flag,
%      oth_bac = get(ui_han(9+oth_sys(k)),'userdata');
%      bac_cplx = termcplx(oth_bac,w);
%      tmp_cplx = oth_cplx./bac_cplx;
%      oth_clos = tmp_cplx./(1+oth_cplx);
%     end

     if open_flag,
      rl=real(oth_cplx); im=imag(oth_cplx);
      set(axs_data(nyq,1+oth_sys(k)),'xdata',rl,'ydata',im,...
                'color','r','vis','on','erase','norm');
      set(axs_data(nyq,4+oth_sys(k)),'xdata',rl,'ydata',-im,...
            'color','b','vis','on');
     end

%     if clos_flag,
%      rl=real(oth_clos); im=imag(oth_clos);
%      set(axs_data(nyq,4+oth_sys(k)),'xdata',rl,'ydata',im,...
%                'color','b','vis','on','erase','norm');
%     end
    end

   elseif any(mode == 12) | any(mode == 14) | any(mode == 15),

    set(axs_data(nyq,1+oth_sys),'vis','off');
    set(axs_data(nyq,4+oth_sys),'vis','off');

   end

   if any(mode == 15), % setting axis limits to full after page layout
    pageview(1,axs_hand(nyq));
   end

  end

  if length(bodm), % bode (magnitude) plot

% update current system frequency response
   if all(mode~=12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(bodm,1+oth_sys),'vis','off');
     set(axs_data(bodm,4+oth_sys),'vis','off');
    elseif any(mode == 13),
     if open_flag,
      set(axs_data(bodm,1+oth_sys),'color','r');
     end
     if clos_flag,
      set(axs_data(bodm,4+oth_sys),'color','b');
     end
    end

    if open_flag & any(mode == 1),
     set(stat_bar,'string','Computing exact open-loop frequency response');
     wdb = w;
     db = 20*log10(abs(sys_cplx));
     set(axs_data(bodm,1+cur_sys),'xdata',wdb,'ydata',db,...
           'color','g','vis','on');
    elseif any(mode == 1),
     set(stat_bar,'string','Removing exact open-loop frequency response');
     set(axs_data(bodm,1+cur_sys),'vis','off');
    end

    if clos_flag & any(mode == 2),
     set(stat_bar,'string','Computing exact closed-loop frequency response');
     ers = 'norm';
     if ~open_flag, ers = 'norm'; end
     wdb = w;
     db = 20*log10(abs(clos_cplx));
     set(axs_data(bodm,4+cur_sys),'xdata',wdb,'ydata',db,...
           'color','m','vis','on');
    elseif any(mode == 2),
     set(stat_bar,'string','Removing exact closed-loop frequency response');
     set(axs_data(bodm,4+cur_sys),'vis','off');
    end

   end

% switching pages or requesting to veiw all systems
   if (any(mode == 12) & show_flag) | (any(mode == 14) & hide_flag) | ...
      (any(mode == 15) & hide_flag),

    for k = 1:2,
     oth_res = get(ui_han(oth_sys(k)),'userdata');
     oth_w = oth_res(1,:);
     oth_cplx = oth_res(2,:);

     if clos_flag,
      oth_bac = get(ui_han(9+oth_sys(k)),'userdata');
      bac_cplx = termcplx(oth_bac,w);
      tmp_cplx = oth_cplx./bac_cplx;
      oth_clos = tmp_cplx./(1+oth_cplx);
     end

     if open_flag,
      wdb = oth_w;
      db = 20*log10(abs(oth_cplx));
      set(axs_data(bodm,1+oth_sys(k)),'xdata',wdb,'ydata',db,'color','r',...
                                    'vis','on','erase','norm');
     end

     if clos_flag,
      wdb = oth_w;
      db = 20*log10(abs(oth_clos));
      set(axs_data(bodm,1+oth_sys(k)),'xdata',wdb,'ydata',db,'color','b',...
                                    'vis','on','erase','norm');
     end
    end

   elseif any(mode == 12) | any(mode == 14) | any(mode == 15),

    set(axs_data(bodm,1+oth_sys),'vis','off');
    set(axs_data(bodm,4+oth_sys),'vis','off');

   end

   if any(mode == 15), % setting axis limits to full after page change
    pageview(1,axs_hand(bodm));
   end

  end

  if length(bodp), % bode (phase) plot

   axsp=[get(axs_hand(bodp),'xlim'),get(axs_hand(bodp),'ylim')];

% update current system frequency response
   if all(mode~=12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(bodp,1+oth_sys),'vis','off');
     set(axs_data(bodp,4+oth_sys),'vis','off');
    elseif any(mode == 13),
     if open_flag,
      set(axs_data(bodp,1+oth_sys),'color','r');
     end
     if clos_flag,
      set(axs_data(bodp,4+oth_sys),'color','b');
     end
    end

    wph = w;
    if open_flag & any(mode == 1),
     set(stat_bar,'string','Computing exact open-loop frequency response');
     ph = phasecor(sys_cplx,[-360,0]);

     brk=find(abs(diff(ph))>170);
     t=1; pht=[]; wpht=[];
     for k=brk,
      pht=[pht,ph(t:k),NaN];
      wpht=[wpht,wph(t:k),NaN];
      t=k+1;
     end
     pht=[pht,ph(t:length(ph))];
     wpht=[wpht,wph(t:length(wph))];

% update current system frequency response
     set(axs_data(bodp,1+cur_sys),'xdata',wpht,'ydata',pht,...
          'color','g','vis','on');
    elseif any(mode == 1),
     set(stat_bar,'string','Removing exact open-loop frequency response');
     set(axs_data(bodp,1+cur_sys),'vis','off');
    end

    if clos_flag & any(mode == 2),
     set(stat_bar,'string','Computing exact closed-loop frequency response');
     ers = 'norm';
     if ~open_flag, ers = 'norm'; end

     ph = phasecor(clos_cplx,[-360,0]);

     brk=find(abs(diff(ph))>170);
     t=1; pht=[]; wpht=[];
     for k=brk,
      pht=[pht,ph(t:k),NaN];
      wpht=[wpht,wph(t:k),NaN];
      t=k+1;
     end
     pht=[pht,ph(t:length(ph))];
     wpht=[wpht,wph(t:length(wph))];

% update current system frequency response
     set(axs_data(bodp,4+cur_sys),'xdata',wpht,'ydata',pht,...
          'color','m','vis','on');
    elseif any(mode == 2),
     set(stat_bar,'string','Removing exact closed-loop frequency response');
     set(axs_data(bodp,4+cur_sys),'vis','off');
    end

   end

% switching pages or requesting to veiw all systems
   if (any(mode == 12) & show_flag) | (any(mode == 14) & hide_flag) | ...
      (any(mode == 15) & hide_flag),

    for k = 1:2,
     oth_res = get(ui_han(oth_sys(k)),'userdata');
     oth_w = oth_res(1,:);
     oth_cplx = oth_res(2,:);

     if clos_flag,
      oth_bac = get(ui_han(9+oth_sys(k)),'userdata');
      bac_cplx = termcplx(oth_bac,w);
      tmp_cplx = oth_cplx./bac_cplx;
      oth_clos = tmp_cplx./(1+oth_cplx);
     end

     wph = oth_w;
     if open_flag,
      ph = phasecor(oth_cplx,[-360,0]);

      brk=find(abs(diff(ph))>170);
      t=1; pht=[]; wpht=[];
      for k2=brk,
       pht=[pht,ph(t:k2),NaN];
       wpht=[wpht,wph(t:k2),NaN];
       t=k2+1;
      end
      pht=[pht,ph(t:length(ph))];
      wpht=[wpht,wph(t:length(wph))];

      set(axs_data(bodp,1+oth_sys(k)),'xdata',wpht,'ydata',pht,'color','r',...
                                    'vis','on','erase','norm');
     end

     if clos_flag,
      ph = phasecor(oth_clos,[-360,0]);

      brk=find(abs(diff(ph))>170);
      t=1; pht=[]; wpht=[];
      for k2=brk,
       pht=[pht,ph(t:k2),NaN];
       wpht=[wpht,wph(t:k2),NaN];
       t=k2+1;
      end
      pht=[pht,ph(t:length(ph))];
      wpht=[wpht,wph(t:length(wph))];

      set(axs_data(bodp,4+oth_sys(k)),'xdata',wpht,'ydata',pht,'color','b',...
                                    'vis','on','erase','norm');
     end

    end
   elseif any(mode == 12) | any(mode == 14) | any(mode == 15),

    set(axs_data(bodp,1+oth_sys),'vis','off');
    set(axs_data(bodp,4+oth_sys),'vis','off');

   end

   if any(mode == 15), % setting axis limits to full after page change
    pageview(1,axs_hand(bodp));
   end

  end

  if length(nic), % nichols plot

   axs = [get(axs_hand(nic),'xlim'),get(axs_hand(nic),'ylim')];

% update current system frequency response
   if all(mode~=12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(nic,1+oth_sys),'vis','off');
%     set(axs_data(nic,4+oth_sys),'vis','off');
    elseif any(mode == 13),
     if open_flag,
      set(axs_data(nic,1+oth_sys),'color','r');
     end
%     if clos_flag,
%      set(axs_data(nic,4+oth_sys),'color','b');
%     end
    end

    if open_flag & any(mode == 1),
     set(stat_bar,'string','Computing exact open-loop frequency response');
     db = 20*log10(abs(sys_cplx));
     ph = phasecor(sys_cplx,[-360,0]);
     brk=find(abs(diff(ph))>170);
     t=1; pht=[]; dbt=[];
     for k=brk,
      pht=[pht,ph(t:k),NaN];
      dbt=[dbt,db(t:k),NaN];
      t=k+1;
     end
     pht=[pht,ph(t:length(ph))];
     dbt=[dbt,db(t:length(db))];

% update current system frequency response
     set(axs_data(nic,1+cur_sys),'xdata',pht,'ydata',dbt,...
           'color','g','vis','on');

    elseif any(mode == 1),
     set(stat_bar,'string','Removing exact open-loop frequency response');
     set(axs_data(nic,1+cur_sys),'vis','off');
    end

%    if clos_flag,
%     ers = 'norm';
%     if ~open_flag, ers = 'norm'; end
%
%     db = 20*log10(abs(clos_cplx));
%     ph = phasecor(clos_cplx,[-360,0]);
%
%     brk=find(abs(diff(ph))>170);
%     t=1; pht=[]; dbt=[];
%     for k=brk,
%      pht=[pht,ph(t:k),NaN];
%      dbt=[dbt,db(t:k),NaN];
%      t=k+1;
%     end
%     pht=[pht,ph(t:length(ph))];
%     dbt=[dbt,db(t:length(db))];

% update current system frequency response
%     set(axs_data(nic,4+cur_sys),'xdata',pht,'ydata',dbt,...
%           'color','m','vis','on');
%    elseif any(mode == 2),
%     set(axs_data(nic,4+cur_sys),'vis','off');
%    end

   end

% switching pages or requesting to veiw all systems
   if (any(mode == 12) & show_flag) | (any(mode == 14) & hide_flag) | ...
      (any(mode == 15) & hide_flag),

    for k = 1:2,
     oth_res = get(ui_han(oth_sys(k)),'userdata');
     oth_w = oth_res(1,:);
     oth_cplx = oth_res(2,:);

%     if clos_flag,
%      oth_bac = get(ui_han(9+oth_sys(k)),'userdata');
%      bac_cplx = termcplx(oth_bac,w);
%      tmp_cplx = oth_cplx./bac_cplx;
%      oth_clos = tmp_cplx./(1+oth_cplx);
%     end

     if open_flag,
      db = 20*log10(abs(oth_cplx));
      ph = phasecor(oth_cplx,[-360,0]);

      brk=find(abs(diff(ph))>170);
      t=1; pht=[]; dbt=[];
      for k2=brk,
       pht=[pht,ph(t:k2),NaN];
       dbt=[dbt,db(t:k2),NaN];
       t=k2+1;
      end
      pht=[pht,ph(t:length(ph))];
      dbt=[dbt,db(t:length(db))];

      set(axs_data(nic,1+oth_sys(k)),'xdata',pht,'ydata',dbt,'color','r',...
                                    'vis','on','erase','norm');
     end

%     if clos_flag,
%      db = 20*log10(abs(oth_clos));
%      ph = phasecor(oth_clos,[-360,0]);
%
%      brk=find(abs(diff(ph))>170);
%      t=1; pht=[]; dbt=[];
%      for k2=brk,
%       pht=[pht,ph(t:k2),NaN];
%       dbt=[dbt,db(t:k2),NaN];
%       t=k2+1;
%      end
%      pht=[pht,ph(t:length(ph))];
%      dbt=[dbt,db(t:length(db))];
%
%      set(axs_data(nic,4+oth_sys(k)),'xdata',pht,'ydata',dbt,'color','b',...
%                                    'vis','on','erase','norm');
%     end

    end
   elseif any(mode == 12) | any(mode == 14) | any(mode == 15),

    set(axs_data(nic,1+oth_sys),'vis','off');
%    set(axs_data(nic,4+oth_sys),'vis','off');

   end

   if any(mode == 15), % setting axis limits to full after page change
    pageview(1,axs_hand(nic));
   end

  end
 else

  set(stat_bar,'string','Removing exact frequency response');
  set(axs_data(:,2:7),'vis','off');
  set(exact,'checked','off');

 end
end

if any(mode == 4), % straight line approximation of bode

 if strcmp(get(straight,'checked'),'off') | length(mode)>mode_len,

  if length(bodm) | length(bodp),
   set(straight,'checked','on');
   if clos_flag,
    [nump,denp] = termextr(plant_mat);
    [numg,deng] = termextr(for_mat);
    [numh,denh] = termextr(bac_mat);
    numcl = conv(conv(nump,numg),denh);
    num = conv(conv(nump,numg),numh);
    den = conv(conv(denp,deng),denh);
    ln = length(num); ld = length(den);
    num = [zeros(1,ld-ln),num];
    den = [zeros(1,ln-ld),den];
    dencl = den + num;
    clos_mat = termpars(numcl,dencl);
   end
  end
  term_mat = termjoin(plant_mat,for_mat,bac_mat);

  if length(bodm), % bode (magnitude)

   if all(mode~=12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(bodm,7+oth_sys),'vis','off');
     set(axs_data(bodm,10+oth_sys),'vis','off');
    elseif any(mode == 13),
     if open_flag,
      set(axs_data(bodm,7+oth_sys),'color','c');
     end
     if clos_flag,
      set(axs_data(bodm,10+oth_sys),'color','c');
     end
    end

    if open_flag & any(mode == 1),
     set(stat_bar,'string','Computing open-loop straight-line approximation');
     [frq_mag,st_mag]=straitln(term_mat,w,1);
     set(axs_data(bodm,7+cur_sys),'xdata',frq_mag,'ydata',st_mag,...
            'color','r','vis','on');
    elseif any(mode==1),
     set(stat_bar,'string','Removing open-loop straight-line approximation');
     set(axs_data(bodm,7+cur_sys),'vis','off');
    end

    if clos_flag & any(mode == 2),
     set(stat_bar,'string','Computing closed-loop straight-line approximation');
     [frq_mag,st_mag]=straitln(clos_mat,w,1);
     set(axs_data(bodm,10+cur_sys),'xdata',frq_mag,'ydata',st_mag,...
            'color','r','vis','on');
    elseif any(mode==2),
     set(stat_bar,'string','Removing closed-loop straight-line approximation');
     set(axs_data(bodm,10+cur_sys),'vis','off');
    end

   end

% switching pages or requesting to veiw all systems
   if (any(mode == 12) & show_flag) | (any(mode == 14) & hide_flag),

    for k = 1:2,
     oth_res = get(ui_han(oth_sys(k)),'userdata');
     oth_w = oth_res(1,:);

     if open_flag,
      oth_plant = get(ui_han(3+oth_sys(k)),'userdata');
      oth_for = get(ui_han(6+oth_sys(k)),'userdata');
      oth_bac = get(ui_han(9+oth_sys(k)),'userdata');
      oth_open = termjoin(oth_plant,oth_for,oth_bac);

      [frq_mag,st_mag]=straitln(oth_open,oth_w,1);

      set(axs_data(bodm,7+oth_sys(k)),'xdata',frq_mag,'ydata',st_mag,...
              'color','c','vis','on');
     end

     if clos_flag,
      [nump,denp] = termextr(oth_plant);
      [numg,deng] = termextr(oth_for);
      [numh,denh] = termextr(oth_bac);
      numcl = conv(conv(nump,numg),denh);
      num = conv(conv(nump,numg),numh);
      den = conv(conv(denp,deng),denh);
      ln = length(num); ld = length(den);
      num = [zeros(1,ld-ln),num];
      den = [zeros(1,ln-ld),den];
      dencl = den + num;
      oth_clos = termpars(numcl,dencl);

      [frq_mag,st_mag]=straitln(oth_clos,oth_w,1);

      set(axs_data(bodm,7+oth_sys(k)),'xdata',frq_mag,'ydata',st_mag,...
              'color','c','vis','on');
     end
    end

   elseif any(mode == 12) | any(mode == 14),

    set(axs_data(bodm,7+oth_sys),'vis','off');
    set(axs_data(bodm,10+oth_sys),'vis','off');

   end

  end

  if length(bodp), % bode (phase)

   if all(mode~=12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(bodp,7+oth_sys),'vis','off');
     set(axs_data(bodp,10+oth_sys),'vis','off');
    elseif any(mode == 13),
     if open_flag,
      set(axs_data(bodp,7+oth_sys),'color','c');
     end
     if clos_flag,
      set(axs_data(bodp,10+oth_sys),'color','c');
     end
    end

    if open_flag & any(mode == 1),
     set(stat_bar,'string','Computing open-loop straight-line phase approximation');
     [frq_mag,st_mag]=straitln(term_mat,w,2);
     set(axs_data(bodp,7+cur_sys),'xdata',frq_mag,'ydata',st_mag,...
            'color','r','vis','on');
    elseif any(mode==1),
     set(stat_bar,'string','Removing open-loop straight-line phase approximation');
     set(axs_data(bodp,7+cur_sys),'vis','off');
    end

    if clos_flag & any(mode == 2),
     set(stat_bar,'string','Computing closed-loop straight-line phase approximation');
     [frq_mag,st_mag]=straitln(clos_mat,w,2);
     set(axs_data(bodp,10+cur_sys),'xdata',frq_mag,'ydata',st_mag,...
            'color','r','vis','on');
    elseif any(mode==2),
     set(stat_bar,'string','Removing closed-loop straight-line phase approximation');
     set(axs_data(bodp,10+cur_sys),'vis','off');
    end

   end

% switching pages or requesting to veiw all systems
   if (any(mode == 12) & show_flag) | (any(mode == 14) & hide_flag),

    for k = 1:2,
     oth_res = get(ui_han(oth_sys(k)),'userdata');
     oth_w = oth_res(1,:);

     if open_flag,
      oth_plant = get(ui_han(3+oth_sys(k)),'userdata');
      oth_for = get(ui_han(6+oth_sys(k)),'userdata');
      oth_bac = get(ui_han(9+oth_sys(k)),'userdata');
      oth_open = termjoin(oth_plant,oth_for,oth_bac);

      [frq_mag,st_mag]=straitln(oth_open,oth_w,2);

      set(axs_data(bodp,7+oth_sys(k)),'xdata',frq_mag,'ydata',st_mag,...
              'color','c','vis','on');
     end

     if clos_flag,
      [nump,denp] = termextr(oth_plant);
      [numg,deng] = termextr(oth_for);
      [numh,denh] = termextr(oth_bac);
      numcl = conv(conv(nump,numg),denh);
      num = conv(conv(nump,numg),numh);
      den = conv(conv(denp,deng),denh);
      ln = length(num); ld = length(den);
      num = [zeros(1,ld-ln),num];
      den = [zeros(1,ln-ld),den];
      dencl = den + num;
      oth_clos = termpars(numcl,dencl);

      [frq_mag,st_mag]=straitln(oth_clos,oth_w,2);

      set(axs_data(bodp,7+oth_sys(k)),'xdata',frq_mag,'ydata',st_mag,...
              'color','c','vis','on');
     end
    end

   elseif any(mode == 12) | any(mode == 14),

    set(axs_data(bodp,7+oth_sys),'vis','off');
    set(axs_data(bodp,10+oth_sys),'vis','off');

   end

  end

 else

  set(stat_bar,'string','Removing straight-line approximation');
  set(axs_data([bodm,bodp],8:13),'vis','off');
  set(straight,'checked','off');

 end

end

if any(mode == 5), % margins

 if (strcmp(get(margins,'checked'),'off') | length(mode)>mode_len) & open_flag,
  set(margins,'checked','on');

  set(stat_bar,'string','Displaying margins');

  if all(mode ~= 12),
   db = 20*log10(abs(sys_cplx));
   deg = phase4(sys_cplx)*180/pi;
   [Gm,Wcg,Pm,Wcp] = termmarg(db,deg,w);
  end

% switching pages or requesting to veiw all systems
  if (any(mode == 12) & show_flag) | ...
     (any(mode == 14) & hide_flag),

   for k = 1:2,
    oth_res = get(ui_han(oth_sys(k)),'userdata');
    oth_w = oth_res(1,:);
    oth_cplx = oth_res(2,:);
    db = 20*log10(abs(oth_cplx));
    deg = phase4(oth_cplx)*180/pi;
    [oGm,oWcg,oPm,oWcp] = termmarg(db,deg,oth_w);
    oth_Gm(k) = oGm; oth_Wcg(k) = oWcg;
    oth_Pm(k) = oPm; oth_Wcp(k) = oWcp;
   end
  end

  if length(nyq), % nyquist gain and phase margins

   if all(mode ~= 12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(nyq,7+oth_sys),'vis','off');
     set(axs_data(nyq,10+oth_sys),'vis','off');
    elseif any(mode == 13),
     set(axs_data(nyq,7+oth_sys),'color','c');
     set(axs_data(nyq,10+oth_sys),'color','c');
    end

    if ~isnan(Wcg),
     set(axs_data(nyq,7+cur_sys),'xdata',[-1,-1/Gm],'ydata',[0,0],...
            'color','k','vis','on');
    else
     set(axs_data(nyq,7+cur_sys),'xdata',0,'ydata',0,'vis','off');
    end

    if ~isnan(Wcp),
     set(axs_data(nyq,10+cur_sys),...
               'xdata',[0,real(exp(i*(-180+Pm)*pi/180))],...
               'ydata',[0,imag(exp(i*(-180+Pm)*pi/180))],...
               'color','k','vis','on');
    else
     set(axs_data(nyq,10+cur_sys),'xdata',0,'ydata',0,'vis','off');
    end

   end

   if (any(mode == 12) & show_flag) | ...
      (any(mode == 14) & hide_flag),

    for k = 1:2,
     if ~isnan(oth_Wcg(k)),
      set(axs_data(nyq,7+oth_sys(k)),'xdata',[-1,-1/oth_Gm(k)],'ydata',[0,0],...
            'color','c','vis','on');
     else
      set(axs_data(nyq,7+oth_sys(k)),'xdata',0,'ydata',0,'vis','off');
     end

     if ~isnan(oth_Wcp(k)),
      set(axs_data(nyq,10+oth_sys(k)),...
                'xdata',[0,real(exp(i*(-180+oth_Pm(k))*pi/180))],...
                'ydata',[0,imag(exp(i*(-180+oth_Pm(k))*pi/180))],...
                'color','c','vis','on');
     else
      set(axs_data(nyq,10+oth_sys(k)),'xdata',0,'ydata',0,'vis','off');
     end

    end

   elseif any(mode == 12) | any(mode == 14),

    set(axs_data(nyq,7+oth_sys),'vis','off');
    set(axs_data(nyq,10+oth_sys),'vis','off');

   end
  end

  if length(bodm), % bode (magnitude)

   if all(mode ~= 12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(bodm,13+oth_sys),'vis','off');
    elseif any(mode == 13),
     set(axs_data(bodm,13+oth_sys),'color','c');
    end

    if ~isnan(Wcg),
     set(axs_data(bodm,13+cur_sys),'xdata',[Wcg,Wcg],'ydata',[0,-20*log10(Gm)],...
      'color','k','vis','on');
    else
     set(axs_data(bodm,13+cur_sys),'xdata',0,'ydata',0,'vis','off');
    end

   end

   if (any(mode == 12) & show_flag) | ...
      (any(mode == 14) & hide_flag),

    for k = 1:2,
     if ~isnan(oth_Wcg),
      set(axs_data(bodm,13+oth_sys(k)),'xdata',[oth_Wcg(k),oth_Wcg(k)],...
                  'ydata',[0,-20*log10(oth_Gm(k))],...
                  'color','c','vis','on');
     else
      set(axs_data(bodm,13+oth_sys(k)),'xdata',0,'ydata',0,'vis','off');
     end

    end

   elseif any(mode == 12) | any(mode == 14),

    set(axs_data(bodm,13+oth_sys),'vis','off');

   end

  end

  if length(bodp), % bode (phase)

   if all(mode ~= 12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(bodp,13+oth_sys),'vis','off');
    elseif any(mode == 13),
     set(axs_data(bodp,13+oth_sys),'color','c');
    end

    if ~isnan(Wcp),
     set(axs_data(bodp,13+cur_sys),'xdata',[Wcp,Wcp],'ydata',[-180,-180+Pm],...
            'color','k','vis','on');
    else
     set(axs_data(bodp,13+cur_sys),'xdata',0,'ydata',0,'vis','off');
    end

   end

   if (any(mode == 12) & show_flag) | ...
      (any(mode == 14) & hide_flag),

    for k = 1:2,
     if ~isnan(oth_Wcp(k)),
      set(axs_data(bodp,13+oth_sys(k)),...
          'xdata',[oth_Wcp(k),oth_Wcp(k)],'ydata',[-180,-180+oth_Pm(k)],...
          'color','c','vis','on');
     else
      set(axs_data(bodp,13+oth_sys(k)),'xdata',0,'ydata',0,'vis','off');
     end
    end

   elseif any(mode == 12) | any(mode == 14),

    set(axs_data(bodp,13+oth_sys),'vis','off');

   end

  end

  if length(nic), % nichols gain and phase margins

   if all(mode~=12),

% switching systems
    if any(mode == 13) & show_flag,
     set(axs_data(nic,7+oth_sys),'vis','off');
     set(axs_data(nic,10+oth_sys),'vis','off');
    elseif any(mode == 13),
     set(axs_data(nic,7+oth_sys),'color','c');
     set(axs_data(nic,10+oth_sys),'color','c');
    end

    if ~isnan(Wcg),
     set(axs_data(nic,7+cur_sys),'xdata',[-180,-180],...
            'ydata',[0,-20*log10(Gm)],'color','k','vis','on');
    else
     set(axs_data(nic,7+cur_sys),'xdata',0,'ydata',0,'vis','off');
    end

    if ~isnan(Wcp),
     set(axs_data(nic,10+cur_sys),'xdata',[-180,-180+Pm],...
            'ydata',[0,0],'color','k','vis','on');
    else
     set(axs_data(nic,10+cur_sys),'vis','off');
    end

   end

   if (any(mode == 12) & show_flag) | ...
      (any(mode == 14) & hide_flag),

    for k = 1:2,
     if ~isnan(oth_Wcg(k)),
      set(axs_data(nic,7+oth_sys(k)),'xdata',[-180,-180],...
            'ydata',[0,-20*log10(oth_Gm(k))],'color','c','vis','on');
     else
      set(axs_data(nic,7+oth_sys(k)),'xdata',0,'ydata',0,'vis','off');
     end

     if ~isnan(oth_Wcp(k)),
      set(axs_data(nic,10+oth_sys(k)),'xdata',[-180,-180+oth_Pm(k)],...
             'ydata',[0,0],'color','c','vis','on');
     else
      set(axs_data(nic,10+oth_sys(k)),'xdata',0,'ydata',0','vis','off');
     end
    end

   elseif any(mode == 12) | any(mode == 14),

    set(axs_data(nic,7+oth_sys),'vis','off');
    set(axs_data(nic,10+oth_sys),'vis','off');

   end

  end

 else

  if open_flag,
   set(stat_bar,'string','Removing margins');
   set(margins,'checked','off');
  end
  set(axs_data(nyq,8:13),'vis','off');
  set(axs_data(bodm,14:16),'vis','off');
  set(axs_data(bodp,14:16),'vis','off');
  set(axs_data(nic,8:13),'vis','off');

 end
end

if any(mode == 12) & show_flag,
 set(show_all,'label','Hide Other Systems');
elseif any(mode == 12),
 set(show_all,'label','Show All Systems');
end

author: Sable

repo: mcbench-benchmarks


text:
% pop_specparams() - Set plotting and statistics parameters for computing
%                    STUDY component spectra.
% Usage:
%       >> STUDY = pop_specparams(STUDY, 'key', 'val');
%
% Inputs:
%   STUDY        - EEGLAB STUDY set
%
% Plot options:
%   'topofreq'   - [real] Plot Spectrum scalp maps at one specific freq. (Hz).
%                  A frequency range [min max] may also be defined (the
%                  spectrum is then averaged over the interval) {default: []}
%   'freqrange'  - [min max] spectral frequency range (in Hz) to plot.
%                  {default: whole frequency range} .
%   'ylim'       - [mindB maxdB] spectral plotting limits in dB
%                  {default: from data}
%   'plotgroups' - ['together'|'apart'] 'together' -> plot subject groups
%                  on the same figure in different colors, else ('apart') on
%                  different figures {default: 'apart'}
%   'plotconditions' - ['together'|'apart'] 'together' -> plot conditions
%                  on the same figure in different colors, else ('apart')
%                  on different figures. Note: keywords 'plotgroups' and
%                  'plotconditions' cannot both be set to 'together'.
%                  {default: 'apart'}
%   'subtractsubjectmean' - ['on'|'off'] subtract individual subject mean
%                  from each spectrum before plotting and computing
%                  statistics. Default is 'off'.
%   'averagechan' - ['on'|'off'] average data channels when several are
%                  selected.
%
% See also: std_specplot()
%
% Authors: Arnaud Delorme, CERCO, CNRS, 2006-

% Copyright (C) Arnaud Delorme, CERCO
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [ STUDY, com ] = pop_specparams(STUDY, varargin);

STUDY = default_params(STUDY);
TMPSTUDY = STUDY;
com = '';
if isempty(varargin)

    enablecond         = fastif(length(STUDY.design(STUDY.currentdesign).variable(1).value)>1, 'on', 'off');
    enablegroup        = fastif(length(STUDY.design(STUDY.currentdesign).variable(2).value)>1, 'on', 'off');
    detachplots        = fastif(strcmpi(STUDY.etc.specparams.detachplots,'on'), 1, 0);
    plotconditions     = fastif(strcmpi(STUDY.etc.specparams.plotconditions, 'together'), 1, 0);
    plotgroups         = fastif(strcmpi(STUDY.etc.specparams.plotgroups,'together'), 1, 0);
    submean            = fastif(strcmpi(STUDY.etc.specparams.subtractsubjectmean,'on'), 1, 0);
    radio_averagechan  = fastif(strcmpi(STUDY.etc.specparams.averagechan,'on'), 1, 0);
    radio_scalptopo    = fastif(isempty(STUDY.etc.specparams.topofreq), 0, 1);
    if radio_scalptopo, radio_averagechan = 0; end;
    if radio_scalptopo+radio_averagechan == 0, radio_scalparray = 1; else radio_scalparray = 0; end;

    cb_radio = [ 'set(findobj(gcbf, ''userdata'', ''radio''), ''value'', 0);' ...
                 'set(gcbo, ''value'', 1);' ...
                 'set(findobj(gcbf, ''tag'', ''topofreq''), ''string'', '''');' ];
    cb_edit  = [ 'set(findobj(gcbf, ''userdata'', ''radio''), ''value'', 0);' ...
                 'set(findobj(gcbf, ''tag'', ''scalptopotext''), ''value'', 1);' ];

    uilist = { ...
           {'style' 'text'       'string' 'Spectrum plotting options' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'text'       'string' 'Frequency [low_Hz high_Hz]' } ...
           {'style' 'edit'       'string' num2str(STUDY.etc.specparams.freqrange) 'tag' 'freqrange' } ...
        {} {'style' 'text'       'string' 'Plot limits [low high]'} ...
           {'style' 'edit'       'string' num2str(STUDY.etc.specparams.ylim) 'tag' 'ylim' } ...
        {} {'style' 'checkbox'   'string' 'Subtract individual subject mean spectrum' 'value' submean 'tag' 'submean' } ...
        {} ...
           {'style' 'text'       'string' 'Spectrum plotting format' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'checkbox'   'string' 'Plot first variable on the same panel' 'value' plotconditions 'enable' enablecond  'tag' 'plotconditions' } ...
        {} {'style' 'checkbox'   'string' 'Plot second variable on the same panel' 'value' plotgroups 'enable' enablegroup 'tag' 'plotgroups' } ...
        {} {'style' 'checkbox'   'string' 'Detach plots' 'value' detachplots 'enable' 'on' 'tag' 'detachtag' } ...
        {} ...
           {'style' 'text'       'string' 'Multiple channels selection' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'radio'      'string' 'Plot channels in scalp array'    'value' radio_scalparray 'tag' 'scalparray'       'userdata' 'radio' 'callback' cb_radio} { } ...
        {} {'style' 'radio'      'string'  'Plot topography at freq. (Hz)' 'value' radio_scalptopo  'tag' 'scalptopotext' 'userdata' 'radio' 'callback' cb_radio} ...
           {'style' 'edit'       'string' num2str(STUDY.etc.specparams.topofreq) 'tag' 'topofreq' 'callback' cb_edit } ...
        {} {'style' 'radio'      'string' 'Average selected channels' 'value' radio_averagechan 'tag' 'averagechan' 'userdata' 'radio' 'callback' cb_radio} { } };
    cbline = [0.07 1.1];
    otherline = [ 0.07 0.6 .3];
    chanline  = [ 0.07 0.8 0.3];
    geometry = { 1 otherline otherline cbline 1 1 cbline cbline cbline 1 1 chanline chanline chanline };
    geomvert = [1.2 1 1 1 0.5 1.2 1 1 1 0.5 1.2 1 1 1 ];

    % component plotting
    % ------------------
    if isnan(STUDY.etc.specparams.topofreq)
        geometry(end-4:end) = [];
        geomvert(end-4:end) = [];
        uilist(end-10:end) = [];
    end;

    [out_param userdat tmp res] = inputgui( 'geometry' , geometry, 'uilist', uilist, 'geomvert', geomvert, ...
                                            'title', 'Spectrum plotting options -- pop_specparams()');
    if isempty(res), return; end;

    % decode inputs
    % -------------
    %if res.plotgroups & res.plotconditions, warndlg2('Both conditions and group cannot be plotted on the same panel'); return; end;
    if res.submean   , res.submean    = 'on'; else res.submean    = 'off'; end;
    if res.plotgroups, res.plotgroups = 'together'; else res.plotgroups = 'apart'; end;
    if res.plotconditions , res.plotconditions  = 'together'; else res.plotconditions  = 'apart'; end;
    if res.detachtag, res.detachtag = 'on'; else res.detachtag = 'off'; end;
    if ~isfield(res, 'topofreq'), res.topofreq = STUDY.etc.specparams.topofreq;
    else res.topofreq  = str2num( res.topofreq );
    end;
    if ~isfield(res, 'averagechan'), res.averagechan = STUDY.etc.specparams.averagechan;
    elseif res.averagechan, res.averagechan = 'on'; else res.averagechan = 'off';
    end;
    res.freqrange = str2num( res.freqrange );
    res.ylim      = str2num( res.ylim );

    % build command call
    % ------------------
    options = {};
    if ~strcmpi( res.plotgroups, STUDY.etc.specparams.plotgroups), options = { options{:} 'plotgroups' res.plotgroups }; end;
    if ~strcmpi( res.plotconditions , STUDY.etc.specparams.plotconditions ), options = { options{:} 'plotconditions'  res.plotconditions  }; end;
    if ~strcmpi( res.detachtag, STUDY.etc.specparams.detachplots), options = { options{:} 'detachplots' res.detachtag}; end;
    if ~strcmpi( res.submean   , STUDY.etc.specparams.subtractsubjectmean ), options = { options{:} 'subtractsubjectmean'  res.submean  }; end;
    if ~isequal(res.topofreq, STUDY.etc.specparams.topofreq),       options = { options{:} 'topofreq' res.topofreq }; end;
    if ~isequal(res.ylim, STUDY.etc.specparams.ylim),               options = { options{:} 'ylim' res.ylim      }; end;
    if ~isequal(res.freqrange, STUDY.etc.specparams.freqrange),     options = { options{:} 'freqrange' res.freqrange }; end;
    if ~isequal(res.averagechan, STUDY.etc.specparams.averagechan), options = { options{:} 'averagechan' res.averagechan }; end;
    if ~isempty(options)
        STUDY = pop_specparams(STUDY, options{:});
        com = sprintf('STUDY = pop_specparams(STUDY, %s);', vararg2str( options ));
    end;
else

    if strcmpi(varargin{1}, 'default')
        STUDY = default_params(STUDY);
    else
        for index = 1:2:length(varargin)
            if ~isempty(strmatch(varargin{index}, fieldnames(STUDY.etc.specparams), 'exact'))
                STUDY.etc.specparams = setfield(STUDY.etc.specparams, varargin{index}, varargin{index+1});
            end;
        end;
    end;
end;

% scan clusters and channels to remove specdata info if freqrange has changed
% ----------------------------------------------------------
if ~isequal(STUDY.etc.specparams.freqrange, TMPSTUDY.etc.specparams.freqrange) | ...
    ~isequal(STUDY.etc.specparams.subtractsubjectmean, TMPSTUDY.etc.specparams.subtractsubjectmean)
    rmfields = { 'specdata' 'specfreqs' };
    for iField = 1:length(rmfields)
        if isfield(STUDY.cluster, rmfields{iField})
            STUDY.cluster = rmfield(STUDY.cluster, rmfields{iField});
        end;
        if isfield(STUDY.changrp, rmfields{iField})
            STUDY.changrp = rmfield(STUDY.changrp, rmfields{iField});
        end;
    end;
end;

function STUDY = default_params(STUDY)
    if ~isfield(STUDY.etc, 'specparams'), STUDY.etc.specparams = []; end;
    if ~isfield(STUDY.etc.specparams, 'topofreq'),             STUDY.etc.specparams.topofreq = []; end;
    if ~isfield(STUDY.etc.specparams, 'freqrange'),            STUDY.etc.specparams.freqrange = []; end;
    if ~isfield(STUDY.etc.specparams, 'ylim'     ),            STUDY.etc.specparams.ylim      = []; end;
    if ~isfield(STUDY.etc.specparams, 'subtractsubjectmean' ), STUDY.etc.specparams.subtractsubjectmean  = 'off'; end;
    if ~isfield(STUDY.etc.specparams, 'plotgroups'),           STUDY.etc.specparams.plotgroups = 'apart'; end;
    if ~isfield(STUDY.etc.specparams, 'plotconditions'),       STUDY.etc.specparams.plotconditions  = 'apart'; end;
    if ~isfield(STUDY.etc.specparams, 'averagechan') ,         STUDY.etc.specparams.averagechan  = 'off'; end;
    if ~isfield(STUDY.etc.specparams, 'detachplots') ,         STUDY.etc.specparams.detachplots  = 'on'; end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [im, H] = repeatability_scale(im, scale_factor)
  im = imresize(im, scale_factor);
  H = eye(3);
  H(1,1) = scale_factor;
  H(2,2) = scale_factor;
  % because indices are 1 based:
  H(1,3) = -scale_factor + 1;
  H(2,3) = -scale_factor + 1;
end

author: hosang

repo: detection-proposals


text:
% /*----------------------------------------------------------------------
%     This file contains a simulation of the cart and pole dynamic system and
%  a procedure for learning to balance the pole.

%
clear all;
SEED=1;
rng(SEED);
SHOW_ANIMATION=1;%%if you want to visualize, set this to 1.

%%%%initialize the network
net=net_init_pole();

addpath(genpath('../CoreModules'));

SHOW_ANIMATION_EVERY_N=100; %every n trials
GAMMA   = 0.9;       % Discount factor for critic.
ACTIONS = 2;

MAX_FAILURES  =  5000;      % Termination criterion.
MAX_STEPS   =     100000;

MaxSteps=[];
Vs=[];
failures=0;
success=0;

opts.use_gpu=0;%don't use gpu for this application, since it will be slow
opts.parameters.mom =0.9;
opts.parameters.lr =1e-1;
opts.parameters.weightDecay=1e-3;
opts.parameters.clip=1e-1;

% Turning on the double buffering to plot the cart and pole
if SHOW_ANIMATION
    h = figure;
    set(h,'doublebuffer','on')
end

% Iterate through the action-learn loop.
while (failures < MAX_FAILURES)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%  RESET STARTS %%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Starting state is (0 0 0 0)
    x         = 0;       % cart position, meters
    x_dot     = 0;       % cart velocity
    theta     = 0;       % pole angle, radians
    theta_dot = 0.0;     % pole angular velocity

    state=[x;x_dot;theta;theta_dot];
    valid=is_valid_state(x,x_dot,theta,theta_dot);

    Inputs=zeros(4,MAX_STEPS);
    opts.dzdy =zeros(ACTIONS,MAX_STEPS);
    acts=zeros(1,MAX_STEPS);

    res(1).x=state;
    [ net,res,opts] = net_ff(net,res,opts);
    P=softmax(res(end).x,[]);
    %Choose action randomly according to the current policy.
    act=1+(rand(1)>P(1));
    samples=1;
    P(act)=P(act)-1;%der of softmaxlogloss
    opts.dzdy(:,samples)=P;
    Inputs(:,samples)=state;
    acts(samples)=act;

    %%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%  RESET  END %%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%

    failed=0;
    while samples < MAX_STEPS && failed==0

        if SHOW_ANIMATION&&(mod(failures,SHOW_ANIMATION_EVERY_N)==0||success)
            plot_Cart_Pole(x,theta)
            if success && samples>1000
                break;
            end
        end

        %Apply action to the simulated cart-pole
        [x,x_dot,theta,theta_dot]=Cart_Pole(act-1,x,x_dot,theta,theta_dot);

        state=[x;x_dot;theta;theta_dot];
        valid=is_valid_state(x,x_dot,theta,theta_dot);

        if valid<0
            %Failure occurred
            failed = 1;
            failures=failures+1;
            MaxSteps=[MaxSteps;samples];
            disp(['Trial was ' int2str(failures) ' steps '  num2str(samples)]);
            %Reinforcement upon failure is -1.
            discounted_r=-exp(2*log(GAMMA).*[samples:-1:1]);

            V=mean(-exp(log(GAMMA).*[samples:-1:1]));%E[V(s)]
            Vs=[Vs;V];
        else
            %Not a failure.r=0.
            failed = 0;
            res(1).x=state;
            [net,res,opts] = net_ff(net,res,opts);
            P=softmax(res(end).x,[]);
            %Choose action randomly according to the current policy.

            act=1+(rand(1)>P(1));
            samples=samples+1;
            acts(samples)=act;
            %P;
            %acts(1:samples)
            P(act)=P(act)-1;%der of softmaxlogloss
            opts.dzdy(:,samples)=P;
            Inputs(:,samples)=state;

        end

        if failed

            %%%%%%%%%%%%%%%%%%%%%
            opts.dzdy=discounted_r.*opts.dzdy(:,1:samples);

            %%This is redundant, just to recalculate some intermediate values
            %%and put them into the right place.
            res(1).x=Inputs(:,1:samples);
            [net,res,opts] = net_ff(net,res,opts);

            %%
            [ net,res,opts ] = net_bp( net,res,opts );
            [ net,res,opts ] = adam( net,res,opts );

        end

    end

    if success && samples>1000
        break;
    end
    if samples>=MAX_STEPS
        success=1;
        if SHOW_ANIMATION==0
            break;
        else
            continue;
        end
    end

  end

if (failures == MAX_FAILURES)
    disp(['Pole not balanced. Stopping after ' int2str(failures) ' failures ' ]);
else
    disp(['Pole balanced successfully for at least ' int2str(MAX_STEPS) ' steps ' ]);
end
close all;
figure;subplot(1,2,1);plot(Vs);title('Values');
subplot(1,2,2);plot(MaxSteps);title('Steps');

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
clear all; close all; clear classes; clc;

%% Set flags.
inspect_only = false;

%% Create shapes.
a = 420;  % lattice constant
t = 1;  % slab thickness
r = 0.29*a;  % hole radius

ad = 25;  % divider for a
td = 10;  % divider for t
dd = 10;  % divider for d = 2*r

mx = 20.5;
my = 7.5;
slab_yn = Box([-mx*a mx*a; -my*a -0.5*a; 0 t], [a/ad, a/ad, t]);
slab_yp = Box([-mx*a mx*a; 0.5*a my*a; 0 t], [a/ad, a/ad, t]);

rod = CircularCylinder(Axis.z, t, [0 0 t/2], r, [2*r/dd, 2*r/dd, t]);

%% Solve the system.
gray = [0.5 0.5 0.5];  % [r g b]
src_loc = 2*a;
[E, H, obj_array, src_array, J] = maxwell_run(...
	'OSC', 1e-9, 2500, ...
	'DOM', {'vacuum', 'none', 1}, [-mx*a mx*a; -my*a my*a; 0 t], [a/ad a/ad t], BC.p, [10*a 2*a 0], ...
	'OBJ', ...
		{'Palik/Si', gray}, Box([-mx*a mx*a; -a/2 a/2; 0 t]), ...
	'SRCJ', PointSrc(Axis.z, [src_loc, 0, 0.5]), ...
	inspect_only);

%% Visualize the solution.
if ~inspect_only
	figure;
	clear opts
% 	opts.withgrid = true;
	opts.withobjsrc = true;
	opts.withabs = true;
	opts.withpml = false;
	opts.phase = pi/2;
	figure(1);
	vis2d(E{Axis.z}, Axis.z, 0.5, obj_array, src_array, opts)
	%%
% 	figure(2);
% 	vis2d(H{Axis.y}, Axis.z, 0.5, obj_array, src_array, opts)

	%%
	flux_loc = 3*a;
	power_right = powerflux_patch(E, H, Axis.x, src_loc + flux_loc);
	power_left = -powerflux_patch(E, H, Axis.x, src_loc - flux_loc);
	fprintf('power:\n');
	fprintf('right = %s\n', num2str(power_right));
	fprintf('left = %s\n', num2str(power_left));
	fprintf('error = %s%%\n',num2str((power_left-power_right)/power_right*100));

	%%
	Sx = poynting(Axis.x, E{Axis.y}, E{Axis.z}, H{Axis.y}, H{Axis.z}, Axis.x, 0);
	[array, l] = Sx.data_original;
	figure(3);
	plot(l{2}, abs(array))
	mx*a - 10*a
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [L,yhat,st] = mci_pb_like (P,M,U,Y)
% Log-likelihood for Preece-Baines model
% FORMAT [L,yhat,st] = mci_pb_like (P,M,U,Y)
%
% P         parameters
% M,U,Y     as usual
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_pb_like.m 6548 2015-09-11 12:39:47Z will $

% Status flag (only used for dynamic systems)
st=[];

T=length(Y);
yhat = mci_pb_gen (P,M,U);
E=sum(sum((Y-yhat).^2));

L = M.logdet_Ce - 0.5*T*log(2*pi);
L = L - 0.5*M.iCe*E;

author: spm


text:
function [H,S,D]=sobi(X,n,p)
% SOBI - Second Order Blind Identification (SOBI) by joint diagonalization of
%          correlation  matrices. THIS CODE ASSUMES TEMPORALLY CORRELATED SIGNALS,
%          and uses correlations across times in performing the signal separation.
%          Thus, estimated time delayed covariance matrices must be nonsingular
%          for at least some time delays.
% Usage:
%         >> winv = sobi(data);
%         >> [winv,act] = sobi(data,n,p);
% Inputs:
%   data - data matrix of size [m,N] ELSE of size [m,N,t] where
%                m is the number of sensors,
%                N is the  number of samples,
%                t is the  number of trials (avoid epoch boundaries)
%         n - number of sources {Default: n=m}
%         p - number of correlation matrices to be diagonalized
%             {Default: min(100, N/3)} Note that for non-ideal data,
%             the authors strongly recommend using at least 100 time delays.
%
% Outputs:
%   winv - Matrix of size [m,n], an estimate of the *mixing* matrix. Its
%          columns are the component scalp maps. NOTE: This is the inverse
%          of the usual ICA unmixing weight matrix. Sphering (pre-whitening),
%          used in the algorithm, is incorporated into winv. i.e.,
%
%             >> icaweights = pinv(winv); icasphere = eye(m);
%
%   act  - matrix of dimension [n,N] an estimate of the source activities
%
%             >> data            = winv            * act;
%                [size m,N]        [size m,n]        [size n,N]
%             >> act = pinv(winv) * data;
%
% Authors:  A. Belouchrani and A. Cichocki (references: See function body)
% Note:     Adapted by Arnaud Delorme and Scott Makeig to process data epochs by
%           computing covariances while respecting epoch boundaries.

% REFERENCES:
% A. Belouchrani, K. Abed-Meraim, J.-F. Cardoso, and E. Moulines, ``Second-order
%  blind separation of temporally correlated sources,'' in Proc. Int. Conf. on
%  Digital Sig. Proc., (Cyprus), pp. 346--351, 1993.
%
%  A. Belouchrani and K. Abed-Meraim, ``Separation aveugle au second ordre de
%  sources correlees,'' in  Proc. Gretsi, (Juan-les-pins),
%  pp. 309--312, 1993.
%
%  A. Belouchrani, and A. Cichocki,
%  Robust whitening procedure in blind source separation context,
%  Electronics Letters, Vol. 36, No. 24, 2000, pp. 2050-2053.
%
%  A. Cichocki and S. Amari,
%  Adaptive Blind Signal and Image Processing, Wiley,  2003.

% Authors note: For non-ideal data, use at least p=100 the time-delayed covariance matrices.
DEFAULT_LAGS = 100;

[m,N,ntrials]=size(X);

if nargin<1 || nargin > 3

    help sobi

elseif nargin==1

    n=m; % Source detection (hum...)
    p=min(DEFAULT_LAGS,ceil(N/3)); % Number of time delayed correlation matrices to be diagonalized

elseif nargin==2

    p=min(DEFAULT_LAGS,ceil(N/3)); % Default number of correlation matrices to be diagonalized
    % Use < DEFAULT_LAGS delays if necessary for short data epochs
end

%
% Make the data zero mean
%
X(:,:)=X(:,:)-kron(mean(X(:,:)')',ones(1,N*ntrials));

%
% Pre-whiten the data based directly on SVD
%
[UU,S,VV]=svd(X(:,:)',0);
Q= pinv(S)*VV';
X(:,:)=Q*X(:,:);

% Alternate whitening code
% Rx=(X*X')/T;
% if m<n, % assumes white noise
%   [U,D]=eig(Rx);
%   [puiss,k]=sort(diag(D));
%   ibl= sqrt(puiss(n-m+1:n)-mean(puiss(1:n-m)));
%    bl = ones(m,1) ./ ibl ;
%   BL=diag(bl)*U(1:n,k(n-m+1:n))';
%   IBL=U(1:n,k(n-m+1:n))*diag(ibl);
% else    % assumes no noise
%    IBL=sqrtm(Rx);
%    Q=inv(IBL);
% end
% X=Q*X;

%
% Estimate the correlation matrices
%
k=1;
pm=p*m; % for convenience
for u=1:m:pm
    k=k+1;
    for t = 1:ntrials
        if t == 1
            Rxp=X(:,k:N,t)*X(:,1:N-k+1,t)'/(N-k+1)/ntrials;
        else
            Rxp=Rxp+X(:,k:N,t)*X(:,1:N-k+1,t)'/(N-k+1)/ntrials;
        end
    end
    M(:,u:u+m-1)=norm(Rxp,'fro')*Rxp;  % Frobenius norm =
end                                    % sqrt(sum(diag(Rxp'*Rxp)))

%
% Perform joint diagonalization
%
epsil=1/sqrt(N)/100;
encore=1;
V=eye(m);
step_n=0;
while encore
    encore=0;
    for p=1:m-1
        for q=p+1:m
            % Perform Givens rotation
            g=[   M(p,p:m:pm)-M(q,q:m:pm)  ;
                M(p,q:m:pm)+M(q,p:m:pm)  ;
                i*(M(q,p:m:pm)-M(p,q:m:pm)) ];
            [vcp,D] = eig(real(g*g'));
            [la,K]=sort(diag(D));
            angles=vcp(:,K(3));
            angles=sign(angles(1))*angles;
            c=sqrt(0.5+angles(1)/2);
            sr=0.5*(angles(2)-j*angles(3))/c;
            sc=conj(sr);
            oui = abs(sr)>epsil ;
            encore=encore | oui ;
            if oui  % Update the M and V matrices
                colp=M(:,p:m:pm);
                colq=M(:,q:m:pm);
                M(:,p:m:pm)=c*colp+sr*colq;
                M(:,q:m:pm)=c*colq-sc*colp;
                rowp=M(p,:);
                rowq=M(q,:);
                M(p,:)=c*rowp+sc*rowq;
                M(q,:)=c*rowq-sr*rowp;
                temp=V(:,p);
                V(:,p)=c*V(:,p)+sr*V(:,q);
                V(:,q)=c*V(:,q)-sc*temp;
            end %% if
        end %% q loop
    end %% p loop
    step_n=step_n+1;
    fprintf('%d step\n',step_n);
end %% while

%
% Estimate the mixing matrix
%
H = pinv(Q)*V;

%
% Estimate the source activities
%
if nargout>1
    S=V'*X(:,:); % estimated source activities
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%CROP_BORDERS Crop the borders of an image or stack of images
%
%   [B, vA, vB, bb_rel] = crop_borders(A, bcol, [padding])
%
%IN:
%   A - HxWxCxN stack of images.
%   bcol - Cx1 background colour vector.
%   padding - scalar indicating how much padding to have in relation to
%             the cropped-image-size (0<=padding<=1). Default: 0
%
%OUT:
%   B - JxKxCxN cropped stack of images.
%   vA     - coordinates in A that contain the cropped image
%   vB     - coordinates in B where the cropped version of A is placed
%   bb_rel - relative bounding box (used for eps-cropping)

% 06/03/15: Improved image cropping thanks to Oscar Hartogensis

function [A, vA, vB, bb_rel] = crop_borders(A, bcol, padding)
if nargin < 3
    padding = 0;
end
[h, w, c, n] = size(A);
if isscalar(bcol)
    bcol = bcol(ones(c, 1));
end
bail = false;
for l = 1:w
    for a = 1:c
        if ~all(col(A(:,l,a,:)) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
bcol = A(ceil(end/2),w,:,1);
bail = false;
for r = w:-1:l
    for a = 1:c
        if ~all(col(A(:,r,a,:)) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
bcol = A(1,ceil(end/2),:,1);
bail = false;
for t = 1:h
    for a = 1:c
        if ~all(col(A(t,:,a,:)) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
bcol = A(h,ceil(end/2),:,1);
bail = false;
for b = h:-1:t
    for a = 1:c
        if ~all(col(A(b,:,a,:)) == bcol(a))
            bail = true;
            break;
        end
    end
    if bail
        break;
    end
end
% Crop the background, leaving one boundary pixel to avoid bleeding on resize
%v = [max(t-padding, 1) min(b+padding, h) max(l-padding, 1) min(r+padding, w)];
%A = A(v(1):v(2),v(3):v(4),:,:);
if padding == 0  % no padding
    padding = 1;
elseif abs(padding) < 1  % pad value is a relative fraction of image size
    padding = sign(padding)*round(mean([b-t r-l])*abs(padding)); % ADJUST PADDING
else  % pad value is in units of 1/72" points
    padding = round(padding);  % fix cases of non-integer pad value
end
if padding > 0  % extra padding
    % Create an empty image, containing the background color, that has the
    % cropped image size plus the padded border
    B = repmat(bcol,(b-t)+1+padding*2,(r-l)+1+padding*2);
    % vA - coordinates in A that contain the cropped image
    vA = [t b l r];
    % vB - coordinates in B where the cropped version of A will be placed
    vB = [padding+1, (b-t)+1+padding, padding+1, (r-l)+1+padding];
    % Place the original image in the empty image
    B(vB(1):vB(2), vB(3):vB(4), :) = A(vA(1):vA(2), vA(3):vA(4), :);
    A = B;
else  % extra cropping
    vA = [t-padding b+padding l-padding r+padding];
    A = A(vA(1):vA(2), vA(3):vA(4), :);
    vB = [NaN NaN NaN NaN];
end
% For EPS cropping, determine the relative BoundingBox - bb_rel
bb_rel = [l-1 h-b-1 r+1 h-t+1]./[w h w h];
end

function A = col(A)
A = A(:);
end

author: marianux

repo: ecg-kit

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


text:
function [exploration,sws,rem] = BrainStates(s,t,f,q,emg,varargin)

%BrainStates - Determine brain state using LFP, EMG and movement.
%
%  USAGE
%
%    [exploration,sws,rem] = BrainStates(s,t,f,q,emg,<options>)
%
%    s              spectrogram
%    t              time bins for spectrogram
%    f              frequency bins for spectrogram
%    q              'quiescence' vector, obtained using QuietPeriods (see NOTE)
%    emg            electromyogram (pass [] if missing)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'nClusters'   number of clusters for K-means clustering (default = 2)
%     'method'      see below (default = 'hippocampus')
%     'nComponents' number of principal components (for 'pca' method only)
%                   (default = automatically computed to account for 85% of
%                   the variance)
%     'show'        plot K-means clustering in feature space ('kmeans'),
%                   clusters on spectrogram ('clusters') or both ('all')
%                   (default = 'none')
%    =========================================================================
%
%  METHODS
%
%    Option 'method' can take several possible values:
%
%     'pca'         PCA on the spectrogram
%     'hippocampus' theta / delta ratio
%     'amygdala'    gamma / low ratio
%     'cortex'      heuristic ratios described in Gervasoni et al. (2004)
%
%  OUTPUT
%
%    exploration    at each time t, whether the rat was exploring
%    sws            at each time t, whether the rat was in slow wave sleep
%    rem            at each time t, whether the rat was in rapid eye movement sleep
%
%  NOTE
%
%    Quiescence q is a time series of boolean values specifying for each
%    timestamp whether the animal was quiet (as determined by QuietPeriods).
%
%  SEE
%
%    See also QuietPeriods, SpectrogramBands.

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

% Defaults
nClusters = 2;
window = 5*1250;
show = 'none';
nComponents = 0;
method = 'hippocampus';

% Check number of parameters
if nargin < 5,
	error('Incorrect number of parameters (type ''help BrainStates'' for details).');
end

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

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i) ' is not a property (type ''help BrainStates'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'nclusters',
			nClusters = varargin{i+1};
			if ~isiscalar(nClusters,'>0'),
				error('Incorrect value for property ''nClusters'' (type ''help BrainStates'' for details).');
			end
		case 'method',
			method = lower(varargin{i+1});
			% 'direct' and 'ratios' are for backward compatibility (deprecated)
			if ~isstring_FMAT(method,'pca','hippocampus','cortex','amygdala','direct','ratios'),
				error('Incorrect value for property ''method'' (type ''help BrainStates'' for details).');
			end
		case 'ncomponents',
			nComponents = varargin{i+1};
			if ~isiscalar(nComponents,'>0'),
				error('Incorrect value for property ''nComponents'' (type ''help BrainStates'' for details).');
			end
		case 'show',
			show = lower(varargin{i+1});
			if ~isstring_FMAT(show,'kmeans','clusters','all','none'),
				error('Incorrect value for property ''show'' (type ''help BrainStates'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help BrainStates'' for details).']);
	end
end

% Compute and interpolate EMG and quiescence at spectrogram timestamps
emg0 = [];
if ~isempty(emg),
	emg0 = Interpolate(emg,t,'trim','off');
	emg0 = Smooth(abs(emg0(:,2)),5);
end
q0 = [];
if all(q(:,2)==0|q(:,2)==1),
	q0 = Interpolate(double(q),t,'trim','off');
	q0 = round(q0(:,2));
else
	q0 = Interpolate(q,t,'trim','off');
	q0 = q0(:,2);
end

% No movement info => undetermined
usable = ~isnan(q0);

% Exploration corresponds to movement periods
exploration = zeros(size(t));
exploration(usable) = ~q0(usable);

% Sleep/rest corresponds to quiescence
sleep = zeros(size(t));
sleep(usable) = q0(usable);
sleep = logical(sleep);

% Get ratios
bands = SpectrogramBands(s,f);

% Determine features for automatic data clustering
switch(method),
	case 'pca',
		% Compute PCA on spectrogram and reduce dimensionality
		S = s(:,sleep)';
		S(:,f>30) = 0;
		[eigenvectors,projected,lambda] = princomp(S,'econ');
		if nComponents == 0,
			nComponents = find(cumsum(lambda)/sum(lambda)>0.85);
			nComponents = nComponents(1);
		end
		eigenvectors = eigenvectors(:,1:nComponents);
		lambda = lambda(1:nComponents);
		projected = projected(:,1:nComponents);
		features = projected;
	case {'hippocampus','direct'},
		% Use theta/delta ratio
		features = bands.ratios.hippocampus(sleep);
	case {'cortex','ratios'},
		% Use heuristic ratios of Gervasoni et al. (2004)
		features = bands.ratios.cortex(sleep);
	case 'amygdala',
		% Use gamma / low ratio
		features = bands.ratios.amygdala(sleep);
end

% Add EMG to feature list
if ~isempty(emg0),
	features = [features emg0(sleep)];
end

% Normalize columns (we need this because K-means clustering uses Euclidean
% distances)
features = zscore(features);

% Cluster (K-means)
cluster = kmeans(features,nClusters);

% Identify clusters
sws = logical(zeros(size(t)));
rem = logical(zeros(size(t)));
% 1) Measure mean ratio for each state
switch(method),
    case 'amygdala',
        ratio = bands.ratios.amygdala(sleep);
    otherwise,
        ratio = bands.ratios.hippocampus(sleep);
end
for i = 1:nClusters,
	meanRatio(i) = mean(ratio(cluster==i));
end
% 2) REM is the quiet state with highest mean ratio
[~,i] = sort(meanRatio(:),1,'descend');
h = i(1);
rem(sleep) = cluster==h;
highest = meanRatio(h);
% 3) SWS is the quiet state with lowest mean ratio
[~,i] = sort(meanRatio(:),1,'ascend');
l = i(1);
sws(sleep) = cluster==l;
lowest = meanRatio(l);
% 4) REM must have a ratio at least twice as high as SWS
if highest < 2*lowest,
	sws(sleep) = 1;
	rem(sleep) = 0;
end

% Show K-means clustering in feature space
if strcmp(show,'kmeans') | strcmp(show,'all'),
	figure;
	% Determine number N of subplots (depends on number of principal components kept; clip at 15)
	n = size(features,2);
	N = min([n*(n-1)/2 21]);
	if N == 0, N = 1; end
	m0 = round(sqrt(N));
	n0 = ceil(N/m0);
	colors = jet(nComponents);
	if strcmp(method,'pca'),
		% Plot principal components
		subplot(m0+1,1,1);hold on;
		l = 'h = legend(';
		for i = 1:nComponents,
			plot(f,eigenvectors(:,i),'color',colors(i,:));
			l = [l '''' int2str(i) ''','];
		end
		l = [l(1:end-1) ',''location'',''northoutside'',''orientation'',''horizontal'');'];
		eval(l);
		set(h,'color',[1 1 1]*.7);
		shift = 1;
	else
		shift = 0;
	end
	% List pairs of principal component IDs to display in successive 2D projections
	colors = jet(nClusters);
	if N > 1,
		feature1 = [];
		feature2 = [];
		i = 1;
		while length(feature1) < N,
			feature1 = [feature1 repmat(i,1,n-i)];
			feature2 = [feature2 i+1:n];
			i = i+1;
		end
		% Plot 2D projections of data in PC space, colored by cluster ID
		for i = 1:N,
			subplot(m0+shift,n0,i+shift*n0);hold on;
			for j = 1:size(features,1),
				plot(features(j,feature1(i)),features(j,feature2(i)),'.','color',colors(cluster(j),:));
			end
			xlabel([int2str(feature1(i)) ' vs ' int2str(feature2(i))]);
		end
	else
		subplot(shift+1,1,shift+1);hold on;
		r = rand(size(features));
		for j = 1:size(features,1),
			plot(r(j),features(j),'.','color',colors(cluster(j),:));
		end
	end
end

%%%%%PLOT%%%%%%

% Redefine ratios for plotting
switch(method)
  case 'amygdala'
    ratio=bands.ratios.amygdala;
  case {'hippocampus','pca'}
    ratio=bands.ratios.hippocampus;
  case 'cortex'
    ratio=bands.ratios.cortex;
end

% Show clusters on spectrogram
if strcmp(show,'clusters') | strcmp(show,'all'),
	figure;
	% 1) Plot spectrogram and ratio
	subplot(2,1,1);hold on;
	PlotColorMap(log(s),1,'cutoffs',[0 12],'x',t,'y',f);
	plot(t,ratio*10,'k');
	l = ['h = legend(''' method ' ratio'',''q'''];
	ylim([0 30]);
	yLim = ylim;
	dy = yLim(2) - yLim(1);
	y0 = yLim(1);
	% ... quiescence...
	plot(t,5*ZeroToOne(q0)+y0+0.5*dy,'b');
	% ... and EMG
	if ~isempty(emg0),
		plot(t,ZeroToOne(emg0)*5+y0+0.65*dy,'r');
		l = [l ',''emg'''];
	end
	l = [l ');'];
	eval(l);
	set(h,'color',[1 1 1]);
	xlim([t(1) t(end)]);
	% 2) Plot spectrogram and clusters
	subplot(2,1,2);hold on;
	PlotColorMap(log(s),1,'cutoffs',[0 12],'x',t,'y',f);
	ylim([0 30]);
	l = 'h = legend(';
	yLim = ylim;
	dy = yLim(2) - yLim(1);
	y0 = yLim(1);
	% Plot clustering output
	c = nan(size(t));
	c(sleep) = cluster-1;
	plot(t,5*ZeroToOne(c)+y0+0.35*dy,'k');
	l = [l '''cluster'',''location'',''north'',''orientation'',''horizontal'');'];
	eval(l);
	% Show SWS/REM
	if ~isempty(sws),
		z = nan(size(t));
		z(sws) = t(sws);
		plot(z,15*~isnan(z),'b','linewidth',5);
	end
	if ~isempty(rem),
		z = nan(size(t));
		z(rem) = t(rem);
		plot(z,15*~isnan(z),'r','linewidth',5);
	end
	xlim([t(1) t(end)]);
	set(h,'color',[1 1 1]);
end


text:
function handle = findfig(figname)
    % function handle = findfig(FIGNAME)
    % Returns the handle of the foremost figure
    % with name FIGNAME, or 0 if not found
    %
    % Keith Rogers 9/26/94

    allfigs = get(groot,'children');
    handle = 0;
    for thisfig = (allfigs')
        if strcmp(get(thisfig,'name'),figname)
            handle = thisfig;
            break
        end
    end

author: CelsoReyes


text:
function [cp,tm,b,t,wh,oocvec] = change_point(t,tt,Mode,tthresh)
% :Usage:
% ::
%
%     [cp,tm,b,t,wh,ooc_vector] = change_point(t,tt,Mode,tthresh)
%
%     [cp,baseline mean,max t value,time of max t] = change_point(Zm,sterr,tt,tthresh)
%
% :Inputs:
%
%   **tt:**
%        is number of baseline timepoints
%
%   **t:**
%        is group t-value timeseries
%
%
%   **Mode: 'thresh' or 'time':**
%        if 'thresh': tthresh = threshold for significant t-values.
%
%        change_point finds the first significant supra-threshold t-value and
%        looks back in time to estimate the change point.
%
%        if Mode = 'time'
%
%        change_point finds the estimated change-point for process determined
%        to be out of control at time tthresh
%
% :Outputs:
%
%   **wh:**
%        indices of out of control points
%
%   **ooc_vector:**
%        indicator vector for ooc points
%

    switch Mode
        case 'thresh'
            % find t-values above threshold
            % see timeseries_mc_pvalue to get threshold

            oocvec = abs(t) > abs(tthresh);
            wh = find(oocvec);
            wh(wh <= tt) = [];      % eliminate values in baseline period

            if isempty(wh),
                b = [];
            else
                b = wh(1);
            end

        case 'time'
            b = tthresh;
            wh = b;
            oocvec = [];
        otherwise
            error('Mode must be thresh or time');
    end

    % tm = tmax, b = time of max t
    %[tm, b] = max(abs(t(tt+1:end)));              % Calculate maximum absolute t-value
    % b = b+tt;                                     % max t time

    tm = t(b);                                      % put sign back in

    % zero-crossing

    if ~isempty(b)
        [a,cp] = max((sign(t(1:b)) ~= sign(tm)).*(1:b));  % change point, last zero crossing
    else
        cp = [];
    end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function cmap= cmap_posneg(m)
%CMAP_POSNEG - Colormap going from cyan/blue over white to red/yellow
%
%Synposis:
% MAP= cmap_rainbow(M)
%
%Input:
% M: Size of the colormap (number of entries). Default value: Same size
%    as current colormap
%
%Output:
% MAP: A colormap matrix of size [M 3]
%
%Example:
% clf;
% colormap(cmap_posneg(65));
% imagesc(toeplitz(1:65)); colorbar;
%
%See also COLORMAP, HSV2RGB, CMAP_HSVFADE

% 01-2005 Benjamin Blankertz

if nargin<1 | isempty(m),
  m= size(get(gcf,'colormap'),1);
end

if mod(m,2)==0,
  warning('m should be odd -> using m+1');
  m= m+1;
end

mb6= floor((m-1)/6);
mq= m-1 - 6*mb6;
m1= mb6;
m2= mb6 + (mq>=2);
m3= mb6 + (mq>=4);
map1= cmap_hsvFade(m1+1, 0, 1, [0.5 1]);
map2= cmap_hsvFade(m2+1, [0 1/6], 1, 1);
map3= cmap_hsvFade(m3+1, 1/6, [1 0], 1);
map4= cmap_hsvFade(m3+1, 3/6, [0 1], 1);
map5= cmap_hsvFade(m2+1, [3/6 4/6], 1, 1);
map6= cmap_hsvFade(m1+1, 4/6, 1, [1 0.5]);
cmap= [map1; map2(2:end,:); map3(2:end,:); ...
       map4(2:end,:); map5(2:end,:); map6(2:end,:)];
cmap= flipud(cmap);

author: bbci


text:
% Convert an edge list to an adjacency list.
%
% INPUTS: edge list, mx3, m - number of edges
% OUTPUTS: adjacency list
%
% Note: Information about edge weights (if any) is lost.
% GB: last updated, September 25, 2012

function adjL = edgeL2adjL(el)

nodes = unique([el(:,1)' el(:,2)']);
adjL=cell(numel(nodes),1);

for e=1:size(el,1); adjL{el(e,1)}=[adjL{el(e,1)},el(e,2)]; end

author: aeolianine

repo: octave-networks-toolbox


text:
function bti_test01 ( )

%*****************************************************************************80
%
%% BTI_TEST01 tests BURGERS_TIME_INVISCID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTI_TEST01\n' );
  fprintf ( 1, '  Test BURGERS_TIME_INVISCID\n' );

  method = 1;
  nx = 81;
  nt = 200;
  t_max = 2.0;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Method: 1, Upwind nonconservative.\n' );
  fprintf ( 1, '  Initial condition: expansion\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_inviscid ( method, @ic_expansion, nx, nt, t_max, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 1 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers, inviscid, initial expansion, upwind nonconservative' )

  filename = 'bti_test01.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_pull1427

% MEM 6gb
% WALLTIME 1:00:00
% DEPENDENCY ft_prepare_mesh ft_prepare_headmodel

%%

% This function creates an hexahedral and tetrahedral volumetric mesh from a
% three-compartments volume to be passed as input to ft_prepare_headmodel with
% the method 'simbio'.

segprob = [];
segprob.brain = false(10,10,10); segprob.brain(4:7,4:7,4:7) = true;
segprob.skull = false(10,10,10); segprob.skull(3:8,3:8,3:8) = true;
segprob.scalp = false(10,10,10); segprob.scalp(2:9,2:9,2:9) = true;
% segprob.air = true(10,10,10);
% segprob.air(2:9,2:9,2:9) = false;
%segprob.brain = false(11,11,11); segprob.brain(4:8,4:8,4:8) = true;
%segprob.skull = false(11,11,11); segprob.skull(3:9,3:9,3:9) = true;
%segprob.scalp = false(11,11,11); segprob.scalp(2:10,2:10,2:10) = true;
segprob.dim = size(segprob.brain);
segprob.unit = 'cm';
segprob.coordsys = 'ctf';
segprob.transform = eye(4);
segprob.transform(1,4) = -0.5;
segprob.transform(2,4) = -0.5;
segprob.transform(3,4) = -0.5;

% it is more difficult to visualize a probabilistic segmentation than an indexed one
segindx = ft_datatype_segmentation(segprob, 'segmentationstyle', 'indexed');

cfg = [];
cfg.funparameter = 'tissue';
cfg.method = 'ortho';
cfg.location = [5 5 5]; % this is the center of the volume, in this plot it will be rounded off to the nearest voxel
ft_sourceplot(cfg, segindx);

% determine the range of the bounding box
[X, Y, Z] = ndgrid(1:segprob.dim(1), 1:segprob.dim(2), 1:segprob.dim(3));
voxpos = ft_warp_apply(segprob.transform, [X(:) Y(:) Z(:)]);
minmaxpos(:,1) = min(voxpos) - 0.5;
minmaxpos(:,2) = max(voxpos) + 0.5;

%%

cfg = [];
cfg.shift = 0.3;
cfg.method = 'hexahedral';
mesh_vol_hex1 = ft_prepare_mesh(cfg, segprob);

figure
ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
hold on
ft_plot_mesh(mesh_vol_hex1, 'surfaceonly', false, 'facecolor', 'none', 'edgecolor', 'm');
view(120, 30)

cfg = [];
cfg.method ='simbio';
cfg.conductivity = [0.33, 0.43, 0.53]; % order follows mesh.tissuelabel
ft_prepare_headmodel(cfg, mesh_vol_hex1)

%%

% make a volume with 2x the resolution
cfg = [];
cfg.resolution = 1/2;
cfg.xrange     = minmaxpos(1,:); % this is more robust than the dim, in case the origin is not inside the volume
cfg.yrange     = minmaxpos(2,:);
cfg.zrange     = minmaxpos(3,:);
cfg.method     = 'nearest';
segprob2 = ft_volumereslice(cfg, segprob);
segindx2 = ft_volumereslice(cfg, segindx);

cfg = [];
cfg.shift = 0.3;
cfg.method = 'hexahedral';
mesh_vol_hex2 = ft_prepare_mesh(cfg, segprob2);

% this one looks OK in terms of alignment
figure
% ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
ft_plot_ortho(segindx2.tissue, 'transform', segindx2.transform, 'location', [5 5 5], 'style', 'intersect'); % this one looks OK in terms of alignment
hold on
ft_plot_mesh(mesh_vol_hex2, 'surfaceonly', false, 'facecolor', 'none', 'edgecolor', 'm');
view(120, 30)

cfg = [];
cfg.method ='simbio';
cfg.conductivity = [0.33, 0.43, 0.53]; % order follows mesh.tissuelabel
ft_prepare_headmodel(cfg, mesh_vol_hex1)

%%

% make a volume with 3x the resolution
cfg = [];
cfg.resolution = 1/3;
cfg.xrange     = minmaxpos(1,:); % this is more robust than the dim, in case the origin is not inside the volume
cfg.yrange     = minmaxpos(2,:);
cfg.zrange     = minmaxpos(3,:);
cfg.method     = 'nearest';
segprob3 = ft_volumereslice(cfg, segprob);
segindx3 = ft_volumereslice(cfg, segindx);

cfg = [];
cfg.shift = 0.3;
cfg.method = 'hexahedral';
mesh_vol_hex3 = ft_prepare_mesh(cfg, segprob3);

figure
% ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
% ft_plot_ortho(segindx2.tissue, 'transform', segindx2.transform, 'location', [5 5 5], 'style', 'intersect');
ft_plot_ortho(segindx3.tissue, 'transform', segindx3.transform, 'location', [5 5 5], 'style', 'intersect'); % this one looks OK in terms of alignment
hold on
ft_plot_mesh(mesh_vol_hex3, 'surfaceonly', false, 'facecolor', 'none', 'edgecolor', 'm');
view(120, 30)

cfg = [];
cfg.method ='simbio';
cfg.conductivity = [0.33, 0.43, 0.53]; % order follows mesh.tissuelabel
ft_prepare_headmodel(cfg, mesh_vol_hex3)

%%

cfg = [];
cfg.method = 'tetrahedral';
mesh_vol_tet = ft_prepare_mesh(cfg, segprob);

figure
ft_plot_ortho(segindx.tissue, 'transform', segindx.transform, 'location', [5 5 5], 'style', 'intersect');
hold on
ft_plot_mesh(mesh_vol_tet, 'surfaceonly', false, 'facecolor', 'none', 'edgecolor', 'm');
view(120, 30)

cfg = [];
cfg.method ='simbio';
cfg.conductivity = [0.33, 0.43, 0.53]; % order follows mesh.tissuelabel
ft_prepare_headmodel(cfg, mesh_vol_tet)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Transcription
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 8/6/2011
classdef Transcription
    %plotting
    methods (Static = true)
        function plotRNA_Polymerase_States(sim, axesHandle, time, ~)
            import edu.stanford.covert.util.ConstantUtil;

            rnaPol = sim.state('RNAPolymerase');

            maxState = max(rnaPol.transcripts.transcriptionUnitLengths);
            if ndims(rnaPol.states) == 3
                colors = zeros([size(rnaPol.states) 3]);
            else
                colors = zeros([size(rnaPol.states) 1 3]);
            end
            for i = 1:size(rnaPol.states, 1)
                for k = 1:size(rnaPol.states, 3)
                    switch rnaPol.states(i, 1, k)
                        case rnaPol.notExistValue
                            colors(i,1,k,:) = [1 0 0];
                        case rnaPol.freeValue
                            colors(i,1,k,:) = [0 1 0];
                        case rnaPol.nonSpecificallyBoundValue
                            colors(i,1,k,:) = [0 0 1];
                        otherwise
                            colors(i,1,k,:) = repmat(max(0, rnaPol.states(i,1,k)) / maxState, 1, 3);
                    end
                end
            end

            colors = reshape(colors, [], size(rnaPol.states, 3), 3);
            image(colors, 'Parent', axesHandle);

            xlabel(axesHandle, 'Time (h)', 'fontSize', 8);
            ylabel(axesHandle, 'RNA Polymerase', 'fontSize', 8);

            xTick = get(axesHandle, 'XTick');
            xTickLabel = round((xTick - 0.5) * max(time) / ConstantUtil.secondsPerHour / max(xTick));
            set(axesHandle, 'XTickLabel', xTickLabel);
        end

        function plotRNA_Polymerase_State_Occupancies(sim, axesHandle, time, ~)
            rnaPol = sim.state('RNAPolymerase');

            labels = cell(4,1);
            labels{rnaPol.activelyTranscribingIndex} = 'Actively Transcribing';
            labels{rnaPol.specificallyBoundIndex} = 'Specifically Bound';
            labels{rnaPol.nonSpecificallyBoundIndex} = 'Non-specifically Bound';
            labels{rnaPol.freeIndex} = 'Free';

            plotHandles = plot(axesHandle, time / 3600, permute(rnaPol.stateOccupancies, [3 1 2]));
            legend(plotHandles, labels);
            xlabel(axesHandle, 'Time (h)', 'fontSize', 8);
            ylabel(axesHandle, 'Occupancy', 'fontSize', 8);
            xlim(axesHandle, [0 max(1, time(end))] / 3600);
        end

        function plotRNA_Polymerases(sim, axesHandle, time, ~)
            rnaPol = sim.state('RNAPolymerase');

            rnaPolymeraseGeneStates = rnaPol.states;
            rnaPolymeraseGenomeStates = zeros(size(rnaPolymeraseGeneStates));
            for i = 1:size(rnaPolymeraseGeneStates,1)
                for k = 1:size(rnaPolymeraseGeneStates,3)
                    switch rnaPolymeraseGeneStates(i, 1, k)
                        case rnaPol.specificallyBoundValue
                            rnaPolymeraseGenomeStates(i, 1, k) = rnaPol.transcripts.transcriptionUnitFivePrimeCoordinates(rnaPol.transcripts.boundTranscriptionUnits(i,1,k));
                        case rnaPol.activelyTranscribingValue
                            rnaPolymeraseGenomeStates(i, 1, k) = rnaPol.transcripts.transcriptionUnitFivePrimeCoordinates(rnaPol.transcripts.boundTranscriptionUnits(i,1,k))+...
                                (2*rnaPol.transcripts.transcriptionUnitDirections(rnaPol.transcripts.boundTranscriptionUnits(i,1,k))-1)*(rnaPolymeraseGeneStates(i,1,k)-1);
                    end
                end
            end

            plot(axesHandle, time/3600, permute(rnaPolymeraseGenomeStates, [3 1 2]), '.');
            xlabel(axesHandle, 'Time (h)', 'fontSize', 8);
            ylabel(axesHandle, 'Genome Coordinate', 'fontSize', 8);
            xlim(axesHandle, [0 max(1, time(end))] / 3600);
            ylim(axesHandle, [rnaPol.freeValue-1 rnaPol.transcripts.genomeLength+1]);
        end

        function plotRNA_Polymerase_Bound_Transcription_Units(sim, axesHandle, time, ~)
            rnaPol = sim.state('RNAPolymerase');

            numTimePoints = size(rnaPol.states,3);

            polymeraseTUs = rnaPol.transcripts.boundTranscriptionUnits;
            boundTUs = zeros(length(rnaPol.transcripts.transcriptionUnitLengths), 1, numTimePoints);
            for i = 1:size(polymeraseTUs,1)
                for k = 1:size(polymeraseTUs,3)
                    if polymeraseTUs(i,1,k) > 0
                        boundTUs(polymeraseTUs(i,1,k),1,k) = ...
                            boundTUs(polymeraseTUs(i,1,k),1,k) + 1;
                    end
                end
            end

            boundTUs = permute(boundTUs, [1 3 2]);
            colors = repmat(boundTUs/(max(max(boundTUs))),[1 1 3]);

            image(colors,'Parent',axesHandle);
            xlabel(axesHandle,'Time (h)','fontSize',8);
            ylabel(axesHandle,'Transcription Unit','fontSize',8);

            xTick = get(axesHandle,'XTick');
            xTickLabel = round((xTick-0.5)*max(time)/3600/max(xTick));
            set(axesHandle,'XTickLabel',xTickLabel);
        end

        function plotNTPs(sim, axesHandle, time, ~)
            p = sim.process('Transcription');

            plotHandles = plot(axesHandle, time / 3600, permute(p.substrates(p.substrateIndexs_ntp, :, :), [3 1 2]));
            legend(plotHandles, p.substrateWholeCellModelIDs(p.substrateIndexs_ntp));
            xlabel(axesHandle, 'Time (h)', 'fontSize', 8);
            ylabel(axesHandle, 'NTPs', 'fontSize', 8);
            xlim(axesHandle, [0 max(1, time(end))] / 3600);
        end

        function plotNMPs(sim, axesHandle, time, ~)
            p = sim.process('Transcription');

            plotHandles = plot(axesHandle, time / 3600, permute(p.substrates(p.substrateIndexs_nmp, :, :), [3 1 2]));
            legend(plotHandles, p.substrateWholeCellModelIDs(p.substrateIndexs_nmp));
            xlabel(axesHandle,'Time (h)', 'fontSize', 8);
            ylabel(axesHandle,'NMPs', 'fontSize', 8);
            xlim(axesHandle, [0 max(1, time(end))] / 3600);
        end

        %plot translation limits
        function plotTranscription_Limits(sim, axesHandle, time, ~)
            p = sim.process('Transcription');

            plotHandles = plot(axesHandle, time / 3600, permute(p.substrates(p.substrateIndexs_ntp), [3 1 2]));
            legend(plotHandles, p.substrateWholeCellModelIDs(p.substrateIndexs_ntp));
            xlabel(axesHandle, 'Time (h)', 'fontSize', 8);
            ylabel(axesHandle, 'NTP', 'fontSize', 8);
            xlim(axesHandle, [0 max(1, time(end))] / 3600);
        end

        function plotBound_Transcription_Factors(sim, axesHandle, time, ~)
            p = sim.process('Transcription');

            transcriptionFactors      = p.enzymes(p.enzymeIndexs_transcriptionFactors, :, :);
            boundTranscriptionFactors = p.boundEnzymes(p.enzymeIndexs_transcriptionFactors, :, :);

            plotHandles = plot(axesHandle, time / 3600, permute(boundTranscriptionFactors./(transcriptionFactors+boundTranscriptionFactors), [3 1 2]));
            legend(plotHandles, p.enzymeNames(p.enzymeIndexs_transcriptionFactors));
            xlabel(axesHandle, 'Time (h)', 'fontSize', 8);
            ylabel(axesHandle, 'Fraction Bound', 'fontSize', 8);
            xlim(axesHandle, [0 max(1, time(end))] / 3600);
            ylim(axesHandle, [0 1]);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [p1,p2,P1_l,P2_l] = idpLin2idpPnts(l)

% IDPLIN2IDPPNTS  IDP line to two IDP points conversion.
%   [P1,P2] = IDPLIN2IDPPNTS(L) extracts the two endpoints of the IDP line
%   L in the form of two IDP points with the same anchor of that of the
%   line.
%
%   [p1,p2,P1_l,P2_l] = IDPLIN2IDPPNTS(...) returns the Jacobians wrt L.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

p1 = l(1:6,:);
p2 = l([1:3 7:9],:);

if nargout > 2

    P1_l = [eye(6) zeros(6,3)];
    P2_l = [eye(3) zeros(3,6) ; zeros(3,6) eye(3)];

end

return

%% jac

syms x y z p1 y1 r1 p2 y2 r2 real
l = [x y z p1 y1 r1 p2 y2 r2]';

[p1,p2,P1_l,P2_l] = idpLin2idpPnts(l);

simplify(P1_l - jacobian(p1,l))
simplify(P2_l - jacobian(p2,l))

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [ won bet happy tt ] = auction( w, iter, e, verbose )
% Subroutine AUCTION is a implemetation of the auction Algorithm
% presented in Turing's Invisible Hand. Although well discribed the
% pseudo-algorithm presented there contained an error which was corrected
% by line 62-65 replacing
%           if ((w(ii,j) - p(j)) >= 0)
%
% [ won bet happiness ] = auction( w, iter, e, verbose )
% Inputs:
%   w (g x b)
%     the rows is the bidders
%     the cols is the goods
%     the values are every bidder i's percieved value of good j
%   iter
%     iterations iterations to do
%   e
%     e is the epsilon - the minumum bet increase (default e = 1/(ngoods+1)
%     )
%
% Copyright (C) 2012 - Pieter V. Reyneke, South Africa
%
% Open Source licence, use whereever you like but at own risk, keep this
% complete copyright notice in the code and please send me updates and
% report errors to pieter.reyneke@gmail.com. Will give recognition, if
% requested, after any valid comment leading to an update was recieved.
%
% Author: PV Reyneke

if (nargin<1)
    w = rand(3)*10;
end

if (nargin<2)
    iter = Inf;
end

if (nargin<4)
    verbose = 1;
end

p   = zeros(1,size(w,2));     % reset all goods' current bets (top - tracks)
won = zeros(1,size(w,2));     % reset cuurent winners for each good
q   = 1:size(w,1);            % create queue for bidders (left - obs)

if (nargin<3)
    e = 1/(length(p)+1);
end

if (verbose)
    display('The percieved value of goods (top) by each bidder (left):');
    display([num2str(w, '%3.0f') ]);
end

i = find(q>0,1);              % take first non-committed bidder
ii = q(i);
tdead = 0;
tic
while (~isempty(ii) && (iter > 0))
    q(q==ii) = 0;             % out of queue
    [v, j] = max(w(ii,:) - p); % search for good w greatest value
                              % above current bid for bidder i
    if ( ((w(ii,j) - (p(j) + e)) >= 0) && ...
         ...                  % if our current bidder i values good j
         ...                  % more than previous bidder plus epsilon
         (won(j) ~= ii) )     % and the bot is not on himself already
        if (won(j)>0)
            k = find(q==0,1); % find first open bidder space
            q(k:(end-1)) = q((k+1):end);
                              % shift left and put
            q(end) = won(j);  % previous highest bidder back in the queue
        end
        won(j) = ii;          % assign current bidder i to good j
        p(j) = p(j) + e;      % and save the new bet
    end

    if (verbose)
        tsss = toc;
        display(['Bets(in epsilons)  : ' num2str(p/e, '%3.0f') ]);
        display(['Bidder             : ' num2str(won, '%3.0f') ]);
        display(['Left todo          : ' num2str(q, '%3.0f') ]);
        display(['-----------------------------------------' ]);
        tdead = toc - tsss;
    end
    i = find(q>0,1);          % take next non-committed bidder
    ii = q(i);
    iter = iter-1;
end
bet = p;
tt = toc;

% measuring each bidder's happiness (the invariant of above loop)
happy = won;
iown = find(won > 0);
for (ii = won(iown))
    iii = won(ii);
    if (~isempty(iii))
        jjj = find(iii==won,1);
        if (~isempty(jjj))
            happy(iii) = (sum(e + w(iii,jjj) - p < w(iii,:) - p) == 0);
        end
    end
end
iown = find(won > 0);
for (ii = won(iown))
    iii = won(ii);
    happy(iii) = happy(iii) | (sum(w(iii,:) - p > 0));
end

j = 1;
cost = 0;
for (i=won)
    cost = cost + w(i,j);
    j = j+1;
end

if (verbose)
    display(['Assignment Cost    : ' num2str(cost, '%5.3f') ]);
    display(['Happiness          : ' num2str(happy, '%3.0f') ]);
    display(['Calc took          : ' num2str((tt - tdead) * 1000, '%5.3f') ' ms' ]);
end


text:
function f = toFunctionIn(disc, coeffs)
%TOFUNCTION   Convert discrete values of a TRIGSPEC to a CHEBFUN.
%   F = TOFUNCTIONIN(DISC, COEFFS) converts the coeffs returned by TRIGSPEC
%   to a CHEBFUN.
%
% See also TOVALUES, TOFUNCTIONOUT.

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

% Get the domain:
dom = disc.domain;

% Create a TRIGTECH object from the coefficients:
tech = trigtech({[], coeffs});

% Create a BNDFUN from the TRIGTECH:
fun{1} = bndfun(tech, struct('domain', dom));

% Create a CHEBFUN from the BNDFUN:
f = chebfun(fun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear, clf
Nfading=1024;  % IFFT size for the Npath x Nfading fading matrix
Nos=8;         % Fading oversampling factor
Npath=2;       % Number of paths
N=10000;
FadingType= 'class';
fm= [100 10];  % Doppler frequency
subplot(211)
[FadingMatrix,tf]= FWGN_ff(Npath,fm,Nfading,Nos,FadingType);
plot([1:Nfading]*tf,10*log10(abs(FadingMatrix(1,:))),'k:')
hold on
plot([1:Nfading]*tf,10*log10(abs(FadingMatrix(2,:))),'k-')
title('Modified FWGN in Frequency Domain')
xlabel('time[s]')
ylabel('Magnitude[dB]')
set(gca,'fontsize',9)
legend('Path 1, f_m=100Hz','Path 2, f_m=10Hz')
axis([0 0.5 -20 5])
subplot(212)
[FadingMatrix,tf]= FWGN_tf(Npath,fm,N,Nfading,Nos,FadingType);
plot([1:N]*tf,10*log10(abs(FadingMatrix(1,:))),'k:')
hold on
plot([1:N]*tf,10*log10(abs(FadingMatrix(2,:))),'k-')
title('Modified FWGN in Time Domain');
title(['Modified FWGN in Time Domain, Nfading=', num2str(Nfading),',Nos=',num2str(Nos),',T_s=',num2str(tf),'s']);
xlabel('time[s]')
ylabel('Magnitude[dB]')
set(gca,'fontsize',9)
legend('Path 1, f_m=100Hz','Path 2, f_m=10Hz')
axis([0 0.5 -20 5])

author: LyricYang