text:
classdef Balls<handle

%--------------------------------------------------------------------------
%Collision Ball Libraries
%Version 1.00
%Created by Stepen
%Created 20 March 2013
%Last modified 20 May 2013
%--------------------------------------------------------------------------
%Constructor Format:
%- Balls()
%  creates an empty Collision Ball object.
%- Balls(x,y,u,v)
%  creates a Collision Ball where x, y, u, and v respectively are column
%  array with same row length defining balls position and velocity.
%- Balls(x,y,u,v,m,r,c)
%  creates a Collision Ball where x, y, u, v, m, r, and c respectively are
%  column array (except for c is 3column array) with same row length
%  defining balls position, velocity, mass, radius, and color.
%Public Non-static Method:
%- addBall(x,y,u,v,m,r,c)
%  add a ball to Collision Ball object.
%- addRandomBall(x,y,u,v,m,r,c)
%  add a random ball to Collision Ball object.
%- removeBall(ballID)
%  remove selected ball from Collision Ball object.
%- isOccupied(x,y,r)
%  check whether a ball at x and y position with radius of r can be added
%  to Collision Ball ojbect without conflict.
%- moveBall(dt)
%  move all balls in Collision Ball object for given time interval.
%- draw()
%  plot all balls in Collision Ball object at current time.
%- play(dt)
%  move and animate all balls in Collision Ball object using given time
%  interval as animation refresh rate.
%Public Static Method:
%- isOutBound(x,y,radius)
%  check a ball at x and y position with radius of r is within Collision
%  Ball object's containment space.
%--------------------------------------------------------------------------

%CodeStart-----------------------------------------------------------------
%Declaring object properties
    properties(SetAccess=protected)
        n_ball
        mass
        radius
        color
        x
        y
        u
        v
        ANIMATIONSTAT=false;
    end
%Declaring constant
    properties(Constant=true)
        MASS_MAX=100;
        RADIUS_MAX=20;
        XMIN=-100;
        XMAX=100;
        YMIN=-100;
        YMAX=100;
        UMIN=-25;
        UMAX=25;
        VMIN=-25;
        VMAX=25;
    end
%Declaring constructor
    methods
        function this=Balls(varargin)
            %Initiating field
            this.n_ball=0;
            this.mass=zeros(0,1);
            this.radius=zeros(0,1);
            this.color=zeros(0,3);
            this.x=zeros(0,1);
            this.y=zeros(0,1);
            this.u=zeros(0,1);
            this.v=zeros(0,1);
            %Assigning input arguments to its appropriate field
            switch nargin
                case 7
                    for i=1:numel(varargin{1})
                        this.addBall(varargin{1}(i),varargin{2}(i),...
                                     varargin{3}(i),varargin{4}(i),...
                                     varargin{5}(i),varargin{6}(i),...
                                     varargin{7}(i,:))
                    end
                case 4
                    for i=1:numel(varargin{1})
                        this.addBall(varargin{1}(i),varargin{2}(i),...
                                     varargin{3}(i),varargin{4}(i),...
                                     10,10,[0,0,0])
                    end
                case 1
                    Balls.checksinglevarargin(varargin{1});
                    for i=1:varargin{1}
                        ballValidity=false;
                        while ~ballValidity
                            %Creating random ball
                            [xr,yr,ur,vr,mr,rr,cr]=Balls.generateball();
                            %Checking random ball validity
                            ballValidity=(~this.isOccupied(xr,yr,rr))&&...
                                         (~Balls.isOutBound(xr,yr,rr));
                            %Adding ball
                            if ballValidity
                                this.addBall(xr,yr,ur,vr,mr,rr,cr)
                            end
                        end
                    end
                case 0

                otherwise
                    error('Unexpected constructor format!');
            end
        end
    end
%Declaring public method
    methods
        %Declaring method to add a ball
        function addBall(this,x,y,u,v,mass,radius,color)
            %Checking input arguments
            switch nargin
                case 1
                    this.addRandomBall()
                    return;
                case 5
                    Balls.checkfourvarargins(x,u,v,y)
                    mass=10;
                    radius=10;
                    color=[0,0,255];
                case 8
                    Balls.checksevenvarargins(x,y,u,v,mass,radius,color);
                otherwise
                    error('Unexpected format of input argument!');
            end
            %Checking for validity
            if this.isOccupied(x,y,radius)
                error('New ball is conflicting with the existing balls!');
            end
            if Balls.isOutBound(x,y,radius)
                error('New ball is out of bounds!');
            end
            %Adding ball
            this.n_ball=this.n_ball+1;
            this.mass=[this.mass;mass];
            this.radius=[this.radius;radius];
            this.color=[this.color;color];
            this.x=[this.x;x];
            this.y=[this.y;y];
            this.u=[this.u;u];
            this.v=[this.v;v];
            %Updating drawing plot (if drawing already exist)
            mainfig=Balls.getfigurehandle();
            if ~isempty(mainfig)
                rectangle('Parent',Balls.getaxeshandle(),...
                          'FaceColor',color/255,...
                          'Curvature',[1,1],...
                          'Position',[x-radius,y-radius,...
                                      2*radius,2*radius]);
            end
        end
        %Declaring method to add a random ball
        function addRandomBall(this)
            %Initiating variable for while loop
            canBeAdded=false;
            %Generating random ball until it is valid
            while ~canBeAdded
                %Generating random position
                [xr,yr,ur,vr,mr,rr,cr]=Balls.generateball();
                %Checking random ball validity
                canBeAdded=(~this.isOccupied(xr,yr,rr))&&...
                           (~Balls.isOutBound(xr,yr,rr));
            end
            %Adding random ball
            this.addBall(xr,yr,ur,vr,mr,rr,cr);
        end
        %Declaring method to remove a ball
        function removeBall(this,ballID)
            %Checking input value
            Balls.checksinglevarargin(ballID);
            if ballID>this.n_ball
                error('Ball ID should be less that number of balls!')
            end
            %Deleting ball
            this.mass(ballID)=[];
            this.radius(ballID)=[];
            this.color(ballID,:)=[];
            this.x(ballID)=[];
            this.y(ballID)=[];
            this.u(ballID)=[];
            this.v(ballID)=[];
            this.n_ball=this.n_ball-1;
            %Updating ball drawing (if drawing already exist)
            mainfig=Balls.getfigurehandle();
            if ~isempty(mainfig)
                rectlist=findobj(mainfig,'Type','rect');
                delete(rectlist(this.n_ball+2-ballID))
            end
        end
        %Declaring method to check whether a position is available
        function stat=isOccupied(this,x,y,radius)
            %Initiating output variable
            stat=false;
            %Checking space for each ball
            for i=1:this.n_ball
                S=norm([this.x(i)-x,this.y(i)-y]);
                if S<(this.radius(i)+radius)
                    stat=true;
                    return
                end
            end
        end
        %Declaring function to move the balls
        function moveBall(this,dt)
            %Predicting ball position without considering collision
            [xf,yf]=predictposition(this,dt);
            %Finding collision
            stat=Balls.findcollision(xf,yf,this.radius);
            %while ~isempty(stat)
                %Correcting the speed of colliding balls
                for i=1:size(stat,1)
                    ID1=stat(i,1);
                    ID2=stat(i,2);
                    [this.u(ID1),this.v(ID1),...
                     this.u(ID2),this.v(ID2)]=Balls.collide(this.mass(ID1),...
                                                            this.x(ID1),...
                                                            this.y(ID1),...
                                                            this.u(ID1),...
                                                            this.v(ID1),...
                                                            this.mass(ID2),...
                                                            this.x(ID2),...
                                                            this.y(ID2),...
                                                            this.u(ID2),...
                                                            this.v(ID2));
                end
                %Finding out-of-bount balls
                hstat=(xf<Balls.XMIN+this.radius)|...
                      (xf>Balls.XMAX-this.radius);
                this.u(hstat)=-this.u(hstat);
                vstat=(yf<Balls.YMIN+this.radius)|...
                      (yf>Balls.YMAX-this.radius);
                this.v(vstat)=-this.v(vstat);
                %Limiting maximum speed
                this.u(this.u<Balls.UMIN)=Balls.UMIN;
                this.u(this.u>Balls.UMAX)=Balls.UMAX;
                this.v(this.v<Balls.VMIN)=Balls.VMIN;
                this.v(this.v>Balls.VMAX)=Balls.VMAX;
                %Recalculating ball position with updated velocity
                [xf,yf]=predictposition(this,dt);
                %Rechecking collision
                %stat=Balls.findcollision(xf,yf,this.radius);
                %Updating ball position
                this.x=xf;
                this.y=yf;
            %end
        end
        %Declaring function to draw the balls
        function draw(this)
            %Checking for existing figure
            mainfig=Balls.getfigurehandle();
            %Creating or getting figure handle
            if isempty(mainfig)
                %Creating brand new figure and axes
                WinSize=[1,1,800,800];
                ScrSize=get(0,'ScreenSize');
                mainfig=figure('Name','Elastic Collision Balls',...
                               'NumberTitle','off',...
                               'Menubar','none',...
                               'Units','pixels',...
                               'Position',[(ScrSize(3)-WinSize(3))/2,...
                                           (ScrSize(4)-WinSize(4))/2,...
                                           WinSize(3),WinSize(4)]);
                mainaxs=axes('Parent',mainfig,...
                             'Units','normalized',...
                             'Position',[0.1,0.1,0.8,0.8]);
            else
                %Getting axes handle
                mainaxs=Balls.getaxeshandle();
                %Deleting previous rectangle object
                delete(Balls.getrecthandle());
            end
            %Drawing white rectangle as field
            rectangle('Parent',mainaxs,...
                      'FaceColor','w',...
                      'Position',[Balls.XMIN,Balls.YMIN,...
                                  Balls.XMAX-Balls.XMIN,...
                                  Balls.YMAX-Balls.YMIN]);
            %Drawing balls as curve rectangle
            for i=1:this.n_ball
                rectangle('Parent',mainaxs,...
                          'FaceColor',this.color(i,:)/255,...
                          'Curvature',[1,1],...
                          'Position',[this.x(i)-...
                                      this.radius(i),...
                                      this.y(i)-...
                                      this.radius(i),...
                                      2*this.radius(i),...
                                      2*this.radius(i)]);
            end
            %Correcting axes
            axis(mainaxs,[Balls.XMIN,Balls.XMAX,Balls.YMIN,Balls.YMAX]);
            axis(mainaxs,'equal');
            axis(mainaxs,'off');
        end
        %Declaring function to animate the balls
        function play(this,dt)
            %Checking for existing figure
            mainfig=Balls.getfigurehandle();
            %Creating or getting figure handle
            if isempty(mainfig)
                %Creating brand new figure and axes
                WinSize=[1,1,800,800];
                ScrSize=get(0,'ScreenSize');
                mainfig=figure('Name','Elastic Collision Balls',...
                               'NumberTitle','off',...
                               'Menubar','none',...
                               'Units','pixels',...
                               'Position',[(ScrSize(3)-WinSize(3))/2,...
                                           (ScrSize(4)-WinSize(4))/2,...
                                           WinSize(3),WinSize(4)]);
                mainaxs=axes('Parent',mainfig,...
                             'Units','normalized',...
                             'Position',[0.1,0.1,0.8,0.8]);
            else
                %Getting axes handle
                mainaxs=Balls.getaxeshandle();
                %Deleting previous rectangle object
                delete(Balls.getrecthandle());
            end
            %Drawing rectangle for visualization field
            rectangle('Parent',mainaxs,...
                      'FaceColor','w',...
                      'Position',[Balls.XMIN,Balls.YMIN,...
                                  Balls.XMAX-Balls.XMIN,...
                                  Balls.YMAX-Balls.YMIN]);
            %Drawing balls as curved rectangle
            ballrect=zeros(this.n_ball,1);
            for i=1:this.n_ball
                ballrect(i)=rectangle('Parent',mainaxs,...
                                      'FaceColor',this.color(i,:)/255,...
                                      'Curvature',[1,1],...
                                      'Position',[this.x(i)-...
                                                  this.radius(i),...
                                                  this.y(i)-...
                                                  this.radius(i),...
                                                  2*this.radius(i),...
                                                  2*this.radius(i)]);
            end
            %Animating balls
            this.ANIMATIONSTAT=true;
            while this.ANIMATIONSTAT
                %Timing process for animation delay
                t=tic;
                %Moving balls
                moveBall(this,dt)
                %Redrawing balls
                for i=1:this.n_ball
                    try
                        set(ballrect(i),'Position',[this.x(i)-...
                                                    this.radius(i),...
                                                    this.y(i)-...
                                                    this.radius(i),...
                                                    2*this.radius(i),...
                                                    2*this.radius(i)])
                    catch
                        this.ANIMATIONSTAT=false;
                        return
                    end
                end
                %Correcting axes
                axis(mainaxs,[Balls.XMIN,Balls.XMAX,Balls.YMIN,Balls.YMAX]);
                axis(mainaxs,'equal');
                axis(mainaxs,'off');
                %Delaying animation
                tt=toc(t);
                pause(dt-tt)
                if (dt-tt)<0
                    disp('Interval time too short for animation!')
                    disp('Animation might not be displayed properly!')
                end
            end
        end
    end
%Declaring public static method
    methods(Static=true)
    %Declaring method to check whether a ball is out of bound
        function stat=isOutBound(x,y,radius)
            %Initiating output variable
            stat=false;
            %Checking space for ball position
            if (x<Balls.XMIN+radius)||(x>Balls.XMAX-radius)||...
               (y<Balls.YMIN+radius)||(y>Balls.YMAX-radius)
                stat=true;
                return;
            end
        end
    end
%Declaring private method
    methods(Access=protected)
        %Declaring function to move ball by neglecting collision
        function [x,y]=predictposition(this,dt)
            x=this.x+(this.u*dt)/2;
            y=this.y+(this.v*dt)/2;
        end
    end
%Declaring static private methods
    methods(Access=protected,...
            Static=true)
        %Declaring function to get Balls' figure
        function fig=getfigurehandle()
            %Initiating output variable
            fig=[];
            %Finding all figure objects
            list=findobj(0,'type','figure');
            %Checking each figure name
            for i =1:numel(list)
                if isequal(get(list(i),'Name'),'Elastic Collision Balls')
                    fig=list(i);
                    return;
                end
            end
        end
        %Declaring object to get Balls' axes
        function axs=getaxeshandle()
            axs=findobj(Balls.getfigurehandle(),'type','axes');
        end
        %Declaring function to get Balls' rectangle object
        function rectlist=getrecthandle()
            rectlist=findobj(Balls.getfigurehandle(),'type','rect');
        end
        %Declaring function to check input argument for constructor
        function checksinglevarargin(var1)
            %Checking varargin array size
            if numel(var1)~=1
                error('Input number of balls must be a scalar!');
            end
            %Checking varargin data type
            if (var1<=0)||(mod(var1,1)~=0)
                error('Input number of balls must be a positive integer!');
            end
        end
        function checkfourvarargins(var1,var2,var3,var4)
            %Checking varargin array size
            if (numel(var1)~=1)||(numel(var2)~=1)||...
               (numel(var3)~=1)||(numel(var4)~=1)
                error('Input x, y, u, and v must be a scalar!');
            end
            %Checking varargin data type
            if (~isreal(var1))||(~isreal(var2))||...
               (~isreal(var3))||(~isreal(var4))
                error('Input x, y, u, and v must be real numbers!');
            end
        end
        function checksevenvarargins(var1,var2,var3,var4,var5,var6,var7)
            %Checking varargin for x,y,u,v
            Balls.checkfourvarargins(var1,var2,var3,var4);
            %Checking varargin dimension size
            if (numel(var5)~=1)||(numel(var6)~=1)||(numel(var7)~=3)
                error(['Input mass, radius, and color must be a 2D',...
                       ' column array (3 column array for color)!']);
            end
            %Checking varargin data type
            if (~isreal(var5))||(var5<=0)%||(var5>Balls.MASS_MAX)
                error(['Input mass must be a positive real number',...
                       ' smaller than ',num2str(Balls.MASS_MAX),'!']);
            end
            if (~isreal(var6))||(var6<=0)||(var6>Balls.RADIUS_MAX)
                error(['Input radius must be a positive real number',...
                       ' smaller than ',num2str(Balls.RADIUS_MAX),'!']);
            end
            if (sum(var7<0)~=0)||(sum(var7>255)~=0)||...
               (sum(mod(var7,1)~=0)~=0)
                error(['Input color must be a positive integer',...
                       ' between 0 and 255!']);
            end
        end
        %Declaring function to generate random ball
        function [xr,yr,ur,vr,mr,rr,cr]=generateball()
            mr=rand*Balls.MASS_MAX;
            rr=rand*Balls.RADIUS_MAX;
            cr=floor(rand(1,3)*256);
            xr=Balls.XMIN+rr+...
                (rand*(Balls.XMAX-Balls.XMIN-(2*rr)));
            yr=Balls.YMIN+rr+...
                (rand*(Balls.YMAX-Balls.XMIN)-(2*rr));
            ur=Balls.UMIN+(rand*(Balls.UMAX-Balls.UMIN));
            vr=Balls.VMIN+(rand*(Balls.VMAX-Balls.VMIN));
        end
        %Declaring function to find colliding balls
        function stat=findcollision(x,y,r)
            %Preallocating array for output variable
            stat=zeros(0,2);
            %Finding collision
            for i=1:numel(x)-1
                if isempty(find(stat(:,2)==i, 1))
                    for j=i+1:numel(x)
                        S=norm([x(i)-x(j),y(i)-y(j)]);
                        if S<(r(i)+r(j))
                            stat=[stat;[i,j]];
                            %break
                        end
                    end
                end
            end
        end
        %Declaring function to resolve 2D collision
        function [uf1,vf1,uf2,vf2]=collide(m1,x1,y1,u1,v1,m2,x2,y2,u2,v2)
            %Finding normal angle between two balls
            theta=atan2(y2-y1,x2-x1);
            %Transforming balls velocity to normal coordinate
            VN1=(u1*cos(theta))+(v1*sin(theta));
            VT1=(-u1*sin(theta))+(v1*cos(theta));
            VN2=(u2*cos(theta))+(v2*sin(theta));
            VT2=(-u2*sin(theta))+(v2*cos(theta));
            %Resolving colllision for normal axis
            VN1_F=((VN1*(m1-m2))+(2*m2*VN2))/(m1+m2);
            VN2_F=((VN2*(m2-m1))+(2*m1*VN1))/(m1+m2);
            %Retransforming balls velocity to original coordinate
            uf1=(VN1_F*cos(theta))-(VT1*sin(theta));
            vf1=(VN1_F*sin(theta))+(VT1*cos(theta));
            uf2=(VN2_F*cos(theta))-(VT2*sin(theta));
            vf2=(VN2_F*sin(theta))+(VT2*cos(theta));
        end
    end
%CodeEnd-------------------------------------------------------------------

end

repo: Aero-matlab


text:
function aviname = MakeJAABAResultsMovie_PositiveClips(expdir,scorefn,varargin)

%%

textcolor = [0,0,0];
fontsize = 14;
trx_linewidth = 1;
target_linewidth = 2;
trxcolor = [0,0,0];
labelcolor = [.7,0,0];
trx_rad = 25;
compression = 'None';
allowedcompressions = {'Indeo3', 'Indeo5', 'Cinepak', 'MSVC', 'RLE', 'None','Uncompressed AVI','Motion JPEG AVI'};
useVideoWriter = exist('VideoWriter','file');
aviname = '';
fps = 10;
usemencoder = true;

[nzoomr,nzoomc,...
  figpos,...
  moviefilestr,trxfilestr,...
  hfig,...
  nframesout,...
  weight_id,weight_tmid,...
  weight_length,weight_sumscore,...
  weight_meanscore,...
  min_sumscore,...
  min_length,...
  pxwidthradius,...
  pxheightradius,...
  awidthradius,...
  aheightradius,...
  pxborder,...
  aborder,...
  nframesborder,...
  labelcolor,...
  useVideoWriter,...
  avifileTempDataFile,aviname,fps,...
  usemencoder,...
  behavior,...
  debug] ...
  = myparse(varargin,...
  'nzoomr',4,'nzoomc',5,...
  'figpos',[50,50,1500,1000],...
  'moviefilestr','movie.ufmf',...
  'trxfilestr','registered_trx.mat',...
  'hfig',1001,...
  'nframesout',500,...
  'weight_id',100,...
  'weight_tmid',1,...
  'weight_length',1,...
  'weight_sumscore',1,...
  'weight_meanscore',.01,...
  'min_sumscore',0,...
  'min_length',0,...
  'pxwidthradius',[],...
  'pxheightradius',[],...
  'awidthradius',3,...
  'aheightradius',3,...
  'pxborder',5,...
  'aborder',1,...
  'nframesborder',5,...
  'labelcolor',labelcolor,...
  'useVideoWriter',useVideoWriter,...
  'avifileTempDataFile','',...
  'aviname',aviname,...
  'fps',fps,...
  'usemencoder',usemencoder,...
  'behavior','',...
  'debug',false);

TMID = 1;
LENGTH = 2;
SUMSCORE = 3;
MEANSCORE = 4;

w = nan(4,1);
w(TMID) = weight_tmid;
w(LENGTH) = weight_length;
w(SUMSCORE) = weight_sumscore;
w(MEANSCORE) = weight_meanscore;

if ~ischar(compression),
  compression = '';
end
if ~ispc && ~useVideoWriter,
  compression = 'None';
end
if ~isempty(compression) && ~any(strcmpi(compression,allowedcompressions)),
  fprintf('Unknown compressor %s\n',compression);
  compression = '';
end

if isempty(aviname),
  [~,experiment_name] = fileparts(expdir);
  if isempty(behavior),
    tmp = regexprep(scorefn,'\.mat$','');
    tmp = regexprep(tmp,'[sS]cores_?','');
    behavior = regexprep(tmp,'[sS]core_?','');
  end
  aviname = fullfile(expdir,sprintf('JAABA%sExamples_%s.avi',behavior,experiment_name));
end

%% open files
moviename = fullfile(expdir,moviefilestr);
[readframe,nframes,fid] = get_readframe_fcn(moviename);
im = readframe(1);
[nr,nc,ncolors] = size(im);

trxname = fullfile(expdir,trxfilestr);
[trx,trxname,loadsucceeded,timestamps] = load_tracks(trxname);
if ~loadsucceeded,
  error('Could not load trx from file %s',trxname);
end
nids = length(trx);
trxendframes = [trx.endframe];
trxfirstframes = [trx.firstframe];
T0 = min(trxfirstframes);
T1 = max(trxendframes);

scorefilename = fullfile(expdir,scorefn);
tmp = load(scorefilename);
scores = cell(1,nids);
labels = cell(1,nids);
for j = 1:nids,
  scores{j} = tmp.allScores.scores{j}(trxfirstframes(j):trxendframes(j));
  labels{j} = tmp.allScores.postprocessed{j}(trxfirstframes(j):trxendframes(j))>0;
 end

%% compute statistics of bouts

X = [];
ids = [];
allt0s = [];
allt1s = [];
alllengths = [];

for i = 1:nids,
  idx = scores{i} > 0;
  [t0s,t1s] = get_interval_ends(idx);
  t1s = t1s - 1;
  sumscores = [0,cumsum(scores{i})];
  allt0s = [allt0s,t0s]; %#ok<AGROW>
  allt1s = [allt1s,t1s]; %#ok<AGROW>
  ids = [ids,i+zeros(size(t0s))]; %#ok<AGROW>
  xcurr = nan(4,numel(t0s));
  xcurr(TMID,:) = (t0s+t1s)/2;
  xcurr(LENGTH,:) = t1s-t0s+1;
  alllengths = [alllengths,xcurr(LENGTH,:)]; %#ok<AGROW>
  xcurr(SUMSCORE,:) = sumscores(t1s+1)-sumscores(t0s);
  xcurr(MEANSCORE,:) = xcurr(SUMSCORE,:) ./ xcurr(LENGTH,:);
  X = [X,xcurr]; %#ok<AGROW>
end
badidx = X(SUMSCORE,:) < min_sumscore | X(LENGTH,:) < min_length;
X(:,badidx) = [];
allt0s(badidx) = [];
allt1s(badidx) = [];
ids(badidx) = [];
nbouts = size(X,2);

% get the rank for each of these stats
for i = 1:size(X,1),
  [~,order] = sort(X(i,:));
  [~,order] = sort(order);
  X(i,:) = order / nbouts; %#ok<AGROW>
end

ischosen = false(1,nbouts);
did = ones(1,nbouts);
dx = ones(size(X));
nframesfilled = zeros([nzoomr,nzoomc]);
boutis = cell(nzoomr,nzoomc);
while true,

  % choose the most different example
  d = weight_id*did + sum(bsxfun(@times,w,dx),1);
  d(ischosen) = -1;
  maxd = max(d);
  idxcurr = find(d == maxd);
  if numel(idxcurr) > 1,
    i = randsample(idxcurr,1);
  else
    i = idxcurr;
  end

  if ischosen(i),
    break;
  end

  % update distances
  did(ids == ids(i)) = 0;
  dx = min(dx,abs(bsxfun(@minus,X,X(:,i))));

  % try to put this somewhere
  [nfilledcurr,j] = min(nframesfilled(:));

  if nfilledcurr >= nframesout,
    break;
  end

  t0 = max(1,allt0s(i)-nframesborder);
  t1 = min(trx(ids(i)).nframes,allt1s(i)+nframesborder);
  n = t1-t0+1;

  [r,c] = ind2sub([nzoomr,nzoomc],j);
  boutis{r,c}(end+1) = i;
  nframesfilled(j) = nframesfilled(j) + n;
  ischosen(i) = true;

end

nframesout = min(nframesout,max(nframesfilled(:)));

tsout = nan(nzoomr,nzoomc,nframesout);
idsout = nan(nzoomr,nzoomc,nframesout);
t1sout = nan(nzoomr,nzoomc,nframesout);

for r = 1:nzoomr,
  for c = 1:nzoomc,
    off = 0;
    for i = 1:numel(boutis{r,c}),
      j = boutis{r,c}(i);
      t0 = max(1,allt0s(j)-nframesborder) + trx(ids(j)).firstframe - 1;
      t1 = min(trx(ids(j)).nframes,allt1s(j)+nframesborder) + trx(ids(j)).firstframe - 1;
      n = t1-t0+1;
      idsout(r,c,off+1:off+n) = ids(j);
      tsout(r,c,off+1:off+n) = t0:t1;
      t1sout(r,c,off+1:off+n) = t1;
      off = off + n;
    end
  end
end

%% output window sizes

a = [trx.a];
meana = nanmedian(a)*4;
if isempty(pxwidthradius),
  pxwidthradius = awidthradius*meana;
end
if isempty(pxheightradius),
  pxheightradius = aheightradius*meana;
end
if isempty(pxborder),
  pxborder = aborder*meana;
end

axh = figpos(4)/nzoomr;
axw = figpos(3)/nzoomc;
pxheightradius1 = min((pxwidthradius)*axh / axw,(nr-1)/2);
pxwidthradius1 = min((pxheightradius)*axw / axh,(nc-1)/2);
pxheightradius = max(pxheightradius1,pxheightradius);
pxwidthradius = max(pxwidthradius1,pxwidthradius);
ax = MakeJAABAResultsMovie_ResetAxis(nc/2,nr/2,1,nr,nc,pxwidthradius,pxheightradius,pxborder);
axround = max(1,[min(nc,[floor(ax(1)),ceil(ax(2))]),...
  min(nr,[floor(ax(3)),ceil(ax(4))])]);
ax = repmat(ax(:),[1,nzoomr,nzoomc]);
axround = repmat(axround(:),[1,nzoomr,nzoomc]);

%% create the figure

figure(hfig);
clf;
set(hfig,'Position',figpos);

hax = createsubplots(nzoomr,nzoomc,0);
hax = reshape(hax,[nzoomr,nzoomc]);
hinfo = nan(nzoomr,nzoomc);
him = nan(nzoomr,nzoomc);
htrx = nan(nzoomr,nzoomc);
htarget = nan(nzoomr,nzoomc);
hlabel = nan(nzoomr,nzoomc);
for i = 1:nzoomr*nzoomc,
  him(i) = imagesc(0,'Parent',hax(i),[0,255]);
  axis(hax(i),'image');
  axis(hax(i),ax(:,i)');
  hold(hax(i),'on');
  htrx(i) = plot(hax(i),nan,nan,'.-','LineWidth',trx_linewidth,'color',trxcolor);
  hlabel(i) = plot(hax(i),nan,nan,'.','LineWidth',trx_linewidth,'color',labelcolor);
  htarget(i) = plot(hax(i),nan,nan,'-','LineWidth',target_linewidth,'color','k');
  hinfo(i) = text(ax(1,i),ax(4,i),'','Parent',hax(i),'Color',textcolor,...
    'FontUnits','pixels','FontSize',fontsize,...
    'HorizontalAlignment','left',...
    'VerticalAlignment','bottom');

  set(hax(i),'XTickLabel',{},'YTickLabel',{});
end
colormap gray;

if ~debug,
  set(hfig,'Visible','off');
end

%%

isfirstframe = true;

for j = 1:nframesout,

  if mod(j,10) == 0,
    fprintf('Frame %d / %d\n',j,nframesout);
  end
  for r = 1:nzoomr,
    for c = 1:nzoomc,
      t = tsout(r,c,j);
      if isnan(t),
        set(him(r,c),'CData',0);
        set(htrx(r,c),'XData',nan,'YData',nan);
        set(hlabel(r,c),'XData',nan,'YData',nan);
        set(htarget(r,c),'XData',nan,'YData',nan);
        set(hinfo(r,c),'String','');
        continue;
      end
      fly = idsout(r,c,j);
      i = t+trx(fly).off;
      im = readframe(t);

      i0 = max(1,min(trx(fly).nframes,i-trx_rad));
      i1 = max(1,min(trx(fly).nframes,i+trx_rad));
      tmpx = trx(fly).x(i0:i1);
      tmpy = trx(fly).y(i0:i1);

      set(htrx(r,c),'XData',tmpx,...
        'YData',tmpy);
      set(hlabel(r,c),'XData',tmpx(labels{fly}(i0:i1)),...
        'YData',tmpy(labels{fly}(i0:i1)));

      updatefly(htarget(r,c),trx(fly).x(i),trx(fly).y(i),...
        trx(fly).theta(i),trx(fly).a(i),trx(fly).b(i));
      if labels{fly}(i),
        set(htarget(r,c),'Color',labelcolor);
      else
        set(htarget(r,c),'Color',trxcolor);
      end

      outofbounds = trx(fly).x(i)-pxborder < ax(1,r,c) || ...
        trx(fly).x(i)+pxborder > ax(2,r,c) || ...
        trx(fly).y(i)-pxborder < ax(3,r,c) || ...
        trx(fly).y(i)+pxborder > ax(4,r,c);

      if j == 1 || tsout(r,c,j)-1 ~= tsout(r,c,j-1) || outofbounds,
        j1 = t1sout(r,c,j)+trx(fly).off;
        ax(:,r,c) = MakeJAABAResultsMovie_ResetAxis(trx(fly).x(i:j1),trx(fly).y(i:j1),1,nr,nc,pxwidthradius,pxheightradius,pxborder);
        axis(hax(r,c),ax(:,r,c)');
        set(hinfo(r,c),'Position',ax([1,4]));
        axround(:,r,c) = max(1,[min(nc,[floor(ax(1,r,c)),ceil(ax(2,r,c))]),...
          min(nr,[floor(ax(3,r,c)),ceil(ax(4,r,c))])]);
        set(hinfo(r,c),'Position',ax([1,4],r,c));
      end
      set(him(r,c),'CData',im(axround(3,r,c):axround(4,r,c),axround(1,r,c):axround(2,r,c)),...
        'XData',axround([1,2],r,c),'YData',axround([3,4],r,c));
      set(hinfo(r,c),'String',sprintf('Target %d, t = %fs',fly,timestamps(t)));
    end
  end

  if isfirstframe,
    if ~debug,

      if useVideoWriter,
        if strcmpi(compression,'None') || strcmpi(compression,'Uncompressed AVI'),
          profile = 'Uncompressed AVI';
        else
          profile = 'Motion JPEG AVI';
        end
        aviobj = VideoWriter(aviname,profile); %#ok<TNMLP>
        set(aviobj,'FrameRate',fps);
        if ~strcmpi(profile,'Uncompressed AVI'),
          set(aviobj,'Quality',100);
        end
        open(aviobj);
      else
        if isempty(avifileTempDataFile),
          aviobj = avifile(aviname,'fps',fps,'quality',100,'compression',compression);  %#ok<REMFF1>
        else
          aviobj = myavifile(aviname,'fps',fps,'quality',100,'compression',compression,...
            'TempDataFile',avifileTempDataFile);
          fprintf('Temporary data file for avi writing: %s\n',aviobj.TempDataFile);
        end
      end
    end

  end

  if ~debug,

    if isfirstframe,
      fr = getframe_invisible(hfig);
      [height,width,~] = size(fr);
      gfdata = getframe_initialize(hfig);
      [fr,height,width] = getframe_invisible_nocheck(gfdata,[height,width],false,false);
    else
      fr = getframe_invisible_nocheck(gfdata,[height,width],false);
    end

    if useVideoWriter,
%       fr = getframe(hfig);
%       height = size(fr.cdata,1);
%       width = size(fr.cdata,2);
      writeVideo(aviobj,fr);
    else
      aviobj = addframe(aviobj,hfig);
    end
    %set(hfig,'Visible','off');
  end

  pos = get(hfig,'Position');
  if any(pos(3:4)~=figpos(3:4)),
    pos(3:4) = figpos(3:4);
    set(hfig,'Position',pos);
  end

  if debug,
    drawnow;
  end

  isfirstframe = false;

end

%% clean up

if ~debug,
  if useVideoWriter,
    close(aviobj);
  else
    aviobj = close(aviobj); %#ok<NASGU>
  end
end

fclose(fid);

%% compress using mencoder

if ~debug && isunix && usemencoder,
  if ~useVideoWriter,
    width = aviobj.Width;
    height = aviobj.Height;
  end
  [path,base,ext] = fileparts(aviname);
  tmpfile = fullfile(path,[base,'_tmp',ext]);
  newheight = 4*ceil(height/4);
  newwidth = 4*ceil(width/4);
  % subtitles are upside down, so encode with subtitles and flip, then flip
  % again
  cmd = sprintf('mencoder %s -o %s -ovc xvid -xvidencopts fixed_quant=4 -vf scale=%d:%d',...
    aviname,tmpfile,newwidth,newheight);
  status = system(cmd);
  if status ~= 0,
    error('Failed to compress avi to %s',tmpfile);
  end
  unix(sprintf('mv %s %s',tmpfile,aviname));
end

author: kristinbranson

repo: JAABA


text:
function value = index_to_level_open ( dim_num, t, order, level_max )

%*****************************************************************************80
%
%% INDEX_TO_LEVEL_OPEN determines the level of a point given its index.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer T(DIM_NUM), the grid index of a point.
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, integer LEVEL_MAX, the level with respect to which the
%    index applies.
%
%    Output, integer VALUE, the first level on which
%    the point associated with the given index will appear.
%
  value = 0;

  for dim = 1 : dim_num

    s = round ( t(dim) );

    s = i4_modp ( t, order );

    if ( s == 0 )

      level = 0;

    else

      level = level_max;

      while ( mod ( s, 2 ) == 0 )
        s = floor ( s / 2 );
        level = level - 1;
      end

    end

    if ( level == 0 )
      level = 1;
    elseif ( level == 1 )
      level = 0;
    end

    value = value + level;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
%
% [node,face]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
%    or
% [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
%
% create the surface and (optionally) tetrahedral mesh of a 3D cylinder
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   c0, c1:  cylinder axis end points
%   r:   radius of the cylinder
%   tsize: maximum surface triangle size on the sphere
%   maxvol: maximu volume of the tetrahedral elements
%   ndiv: approximate the cylinder surface into ndiv flat pieces, if
%         ignored, ndiv is set to 20
%
% output:
%   node: node coordinates, 3 columns for x, y and z respectively
%   face: integer array with dimensions of NB x 3, each row represents
%         a surface mesh triangle
%   elem: (optional) integer array with dimensions of NE x 4, each row
%         represents a tetrahedron
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<3)
    error('you must at least provide c0, c1, and r');
end
if(r<=0 || all(c0==c1))
    error('invalid cylinder parameters');
end
c0=c0(:);
c1=c1(:);
len=sqrt(sum((c0-c1).*(c0-c1)));
% define the axial vector v0 and a perpendicular vector t
v0=c1-c0;

% calculate the cylinder end face nodes
if(nargin<6) ndiv=20; end

dt=2*pi/ndiv;
theta=dt:dt:2*pi;
cx=r*cos(theta);
cy=r*sin(theta);
p0=[cx(:) cy(:) zeros(ndiv,1)];
p1=[cx(:) cy(:) len*ones(ndiv,1)];
pp=[p0;p1];
no=rotatevec3d(pp,v0)+repmat(c0',size(pp,1),1);

count=1;
for i=1:ndiv-1
   fc{count}={[i i+ndiv i+ndiv+1 i+1],1}; count=count+1;
end
i=ndiv;
fc{count}={[i i+ndiv 1+ndiv 1],1}; count=count+1;
fc{count}={1:ndiv,2};count=count+1;  % bottom inner circle
fc{count}={1+ndiv:2*ndiv,3};count=count+1;  % top inner circle

if(nargin==3)
    tsize=len/10;
end
if(nargin<5)
    maxvol=tsize*tsize*tsize;
end
[node,elem,face]=surf2mesh(no,fc,min(no),max(no),1,maxvol,[0 0 1],[],0);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function nam = get_name(com,detailed)

nam = com.name;

%%%% SGE support
if isdeferred(com), nam = get_name(com.deferred, nam, com.is_data); end
%%%%

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = process_extract_values( varargin )
% PROCESS_EXTRACT_VALUES Extract values from multiple files and concatenate in a matrix file.
%
% USAGE:                   OutputFiles = process_extract_values('Run',     sProcess, sInputs)
%       [ConcatMat, FileType, matName] = process_extract_values('Extract', sProcess, sInputs, OPTIONS)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015-2022

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Extract values';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Extract';
    sProcess.Index       = 350;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Statistics?highlight=%28Extract+values%29#Histograms';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === GENERIC EXTRACT OPTIONS
    sProcess = DefineExtractOptions(sProcess);

    % === CONCATENATE
    sProcess.options.dim.Comment = {'Concatenate signals (dimension 1)', 'Concatenate time (dimension 2)'};
    sProcess.options.dim.Type    = 'radio';
    sProcess.options.dim.Value   = 2;
    % === OUTPUT COMMENT
    sProcess.options.Comment.Comment = 'Comment (empty=default): ';
    sProcess.options.Comment.Type    = 'text';
    sProcess.options.Comment.Value   = '';
end

%% ===== DEFINE EXTRACT OPTIONS =====
function sProcess = DefineExtractOptions(sProcess)
    % === SELECT: TIME WINDOW
    sProcess.options.timewindow.Comment    = 'Time window:';
    sProcess.options.timewindow.Type       = 'timewindow';
    sProcess.options.timewindow.Value      = [];
    sProcess.options.timewindow.InputTypes = {'data', 'results', 'timefreq', 'matrix'};
    % === SELECT: CHANNELS
    sProcess.options.sensortypes.Comment    = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type       = 'text';
    sProcess.options.sensortypes.Value      = '';
    sProcess.options.sensortypes.InputTypes = {'data'};
    % === SELECT: FREQUENCY RANGE
    sProcess.options.freqrange.Comment    = 'Frequency range: ';
    sProcess.options.freqrange.Type       = 'freqrange';
    sProcess.options.freqrange.Value      = [];
    sProcess.options.freqrange.InputTypes = {'timefreq'};
    % === SELECT: ROWS
    sProcess.options.rows.Comment    = 'Signals names or indices (empty=all): ';
    sProcess.options.rows.Type       = 'text';
    sProcess.options.rows.Value      = '';
    sProcess.options.rows.InputTypes = {'timefreq', 'matrix'};

    % === SCOUTS SELECTION ===
    sProcess.options.scoutsel.Comment    = 'Use scouts';
    sProcess.options.scoutsel.Type       = 'scout_confirm';
    sProcess.options.scoutsel.Value      = {};
    sProcess.options.scoutsel.InputTypes = {'results'};
    % === SCOUTS: FUNCTION
    sProcess.options.scoutfunc.Comment    = {'Mean', 'Max', 'PCA', 'Std', 'All', 'Scout function:'};
    sProcess.options.scoutfunc.Type       = 'radio_line';
    sProcess.options.scoutfunc.Value      = 1;
    sProcess.options.scoutfunc.InputTypes = {'results'};

    % === NORM XYZ
    sProcess.options.isnorm.Comment    = 'Compute absolute values (or norm for unconstrained sources)';
    sProcess.options.isnorm.Type       = 'checkbox';
    sProcess.options.isnorm.Value      = 0;
    sProcess.options.isnorm.InputTypes = {'results'};
    % === ABSOLUTE VALUE
    sProcess.options.isabs.Comment    = 'Compute absolute values';
    sProcess.options.isabs.Type       = 'checkbox';
    sProcess.options.isabs.Value      = 0;
    sProcess.options.isabs.InputTypes = {'data', 'timefreq', 'matrix'};

    % === AVERAGE: TIME
    sProcess.options.avgtime.Comment    = 'Average selected time window';
    sProcess.options.avgtime.Type       = 'checkbox';
    sProcess.options.avgtime.Value      = 0;
    sProcess.options.avgtime.InputTypes = {'data', 'results', 'timefreq', 'matrix'};
    % === AVERAGE: CHANNELS
    sProcess.options.avgrow.Comment    = 'Average selected signals';
    sProcess.options.avgrow.Type       = 'checkbox';
    sProcess.options.avgrow.Value      = 0;
    sProcess.options.avgrow.InputTypes = {'data', 'timefreq', 'matrix'};
    % === AVERAGE: FREQUENCY
    sProcess.options.avgfreq.Comment    = 'Average selected frequency band';
    sProcess.options.avgfreq.Type       = 'checkbox';
    sProcess.options.avgfreq.Value      = 0;
    sProcess.options.avgfreq.InputTypes = {'timefreq'};
    % === MATCH ROWS WITH NAMES
    sProcess.options.matchrows.Comment    = 'Match signals between files using their names';
    sProcess.options.matchrows.Type       = 'checkbox';
    sProcess.options.matchrows.Value      = 1;
    sProcess.options.matchrows.InputTypes = {'timefreq', 'matrix'};
end

%% ===== GET EXTRACT OPTIONS =====
function OPTIONS = GetExtractOptions(sProcess, sInputs)
    % Time window
    if isfield(sProcess.options, 'timewindow') && ~isempty(sProcess.options.timewindow) && ~isempty(sProcess.options.timewindow.Value) && iscell(sProcess.options.timewindow.Value)
        OPTIONS.TimeWindow = sProcess.options.timewindow.Value{1};
    else
        OPTIONS.TimeWindow = [];
    end
    % Sensor type
    if ismember(sInputs(1).FileType, {'data'}) && isfield(sProcess.options, 'sensortypes') && ~isempty(sProcess.options.sensortypes) && ~isempty(sProcess.options.sensortypes.Value)
        OPTIONS.SensorTypes = sProcess.options.sensortypes.Value;
    else
        OPTIONS.SensorTypes = [];
    end
    % Row indices
    if ismember(sInputs(1).FileType, {'results', 'timefreq', 'matrix'}) && isfield(sProcess.options, 'rows') && ~isempty(sProcess.options.rows) && ~isempty(sProcess.options.rows.Value)
        OPTIONS.Rows = sProcess.options.rows.Value;
    else
        OPTIONS.Rows = [];
    end
    % Freq indices
    if ismember(sInputs(1).FileType, {'timefreq'}) && isfield(sProcess.options, 'freqrange') && isfield(sProcess.options.freqrange, 'Value') && iscell(sProcess.options.freqrange.Value) && (length(sProcess.options.freqrange.Value) == 3) && (length(sProcess.options.freqrange.Value{1}) == 2)
        OPTIONS.FreqRange = sProcess.options.freqrange.Value{1};
    else
        OPTIONS.FreqRange = [];
    end
    % Scouts: Selection
    if isfield(sProcess.options, 'scoutsel') && isfield(sProcess.options.scoutsel, 'Value') && isfield(sProcess.options, 'scoutfunc') && isfield(sProcess.options.scoutfunc, 'Value')
        OPTIONS.ScoutSel = sProcess.options.scoutsel.Value;
        switch lower(sProcess.options.scoutfunc.Value)
            case {1, 'mean'}, OPTIONS.ScoutFunc = 'mean';
            case {2, 'max'},  OPTIONS.ScoutFunc = 'max';
            case {3, 'pca'},  OPTIONS.ScoutFunc = 'pca';
            case {4, 'std'},  OPTIONS.ScoutFunc = 'std';
            case {5, 'all'},  OPTIONS.ScoutFunc = 'all';
            otherwise,  bst_report('Error', sProcess, [], 'Invalid scout function.');  return;
        end
    else
        OPTIONS.ScoutSel = [];
    end
    % Absolute values / Norm
    OPTIONS.isAbsolute = 0;
    if isfield(sProcess.options, 'isabs') && isfield(sProcess.options.isabs, 'Value')
        OPTIONS.isAbsolute = sProcess.options.isabs.Value;
    end
    if isfield(sProcess.options, 'isnorm') && isfield(sProcess.options.isnorm, 'Value')
        OPTIONS.isAbsolute = sProcess.options.isnorm.Value;
    end
    % Averages
    if isfield(sProcess.options, 'avgtime') && isfield(sProcess.options.avgtime, 'Value')
        OPTIONS.isAvgTime = sProcess.options.avgtime.Value;
    else
        OPTIONS.isAvgTime = 0;
    end
    if isfield(sProcess.options, 'avgrow') && isfield(sProcess.options.avgrow, 'Value')
        OPTIONS.isAvgRow = sProcess.options.avgrow.Value;
    else
        OPTIONS.isAvgRow = 0;
    end
    if isfield(sProcess.options, 'avgfreq') && isfield(sProcess.options.avgfreq, 'Value')
        OPTIONS.isAvgFreq = sProcess.options.avgfreq.Value;
    else
        OPTIONS.isAvgFreq = 0;
    end
    % Match signals between files using their names
    if isfield(sProcess.options, 'matchrows') && isfield(sProcess.options.matchrows, 'Value') && ~isempty(sProcess.options.matchrows.Value)
        OPTIONS.isMatchRows = sProcess.options.matchrows.Value;
    else
        OPTIONS.isMatchRows = 1;
    end
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess)
    Comment = [sProcess.Comment, ': [', process_extract_time('GetTimeString', sProcess) ']'];
    if isfield(sProcess.options, 'freqrange') && isfield(sProcess.options.freqrange, 'Value') && iscell(sProcess.options.freqrange.Value) && (length(sProcess.options.freqrange.Value) == 3) && (length(sProcess.options.freqrange.Value{1}) == 2)
        FreqRange = sProcess.options.freqrange.Value{1};
        if (FreqRange(1) == FreqRange(2))
            Comment = [Comment, ' ' num2str(FreqRange(1)) 'Hz'];
        else
            Comment = [Comment, ' ' num2str(FreqRange(1)) '-' num2str(FreqRange(2)) 'Hz'];
        end
    end
    if isfield(sProcess.options, 'sensortypes') && isfield(sProcess.options.sensortypes, 'Value') && ~isempty(sProcess.options.sensortypes.Value)
        Comment = [Comment, ' ', sProcess.options.sensortypes.Value];
    end
    if isfield(sProcess.options, 'rows') && isfield(sProcess.options.rows, 'Value') && ~isempty(sProcess.options.rows.Value)
        Comment = [Comment, ' ', sProcess.options.rows.Value];
    end
    if isfield(sProcess.options, 'scoutsel') && isfield(sProcess.options.scoutsel, 'Value') && iscell(sProcess.options.scoutsel.Value) && (length(sProcess.options.scoutsel.Value) >= 2) && ~isempty(sProcess.options.scoutsel.Value{1,2})
        listScouts = sProcess.options.scoutsel.Value{1,2};
        scoutComment = sprintf(' %s', listScouts{:});
        if (length(scoutComment) > 20)
            Comment = [Comment, ' ', num2str(length(listScouts)), ' scouts'];
        else
            Comment = [Comment, scoutComment];
        end
    end
    if (isfield(sProcess.options, 'isabs') && isfield(sProcess.options.isabs, 'Value') && sProcess.options.isabs.Value) || ...
       (isfield(sProcess.options, 'isnorm') && isfield(sProcess.options.isnorm, 'Value') && sProcess.options.isnorm.Value)
        Comment = [Comment, ' abs'];
    end
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Initialize returned variable
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get generic options
    OPTIONS = GetExtractOptions(sProcess, sInputs);
    % Output options
    OPTIONS.Dim = sProcess.options.dim.Value;

    % ===== EXTRACT =====
    [newMat, FileType] = Extract(sProcess, sInputs, OPTIONS);
    if isempty(newMat)
        return;
    end

    % ===== SAVE FILE =====
    % Get output study
    [sStudy, iStudy, Comment] = bst_process('GetOutputStudy', sProcess, sInputs);
    % Comment: forced in the options
    if isfield(sProcess.options, 'Comment') && isfield(sProcess.options.Comment, 'Value') && ~isempty(sProcess.options.Comment.Value)
        newMat.Comment = sProcess.options.Comment.Value;
    % Comment: Process default
    else
        newMat.Comment = FormatComment(sProcess);
        if (length(sInputs) > 1)
            newMat.Comment = [newMat.Comment, sprintf(' (%d files)', length(sInputs))];
        end
    end
    % Timefreq: the file tag might be longer than the file type
    if strcmpi(FileType, 'timefreq')
        fileTag = bst_process('GetFileTag', sInputs(1).FileName);
    else
        fileTag = FileType;
    end
    % Output filename
    OutputFiles{1} = bst_process('GetNewFilename', bst_fileparts(sStudy.FileName), [fileTag '_concat']);
    % Save on disk
    bst_save(OutputFiles{1}, newMat, 'v6');
    % Register in database
    db_add_data(iStudy, OutputFiles{1}, newMat);
end

%% ===== EXTRACT =====
function [newMat, newFileType, matName] = Extract(sProcess, sInputs, OPTIONS)

    % ===== GET INPUTS =====
    % Intialize returned variables
    newMat = [];
    newFileType = [];
    OutValue = cell(1, length(sInputs));
    OutNames = {};
    OutEvents = [];
    OutTFmask = [];
    ChannelFile = [];
    ChannelFlag = [];
    % Default options
    defOPTIONS = struct(...
        'TimeWindow',  [], ...
        'SensorTypes', [], ...
        'Rows',        [], ...
        'FreqRange',   [], ...
        'ScoutSel',    [], ...
        'ScoutFunc',   [], ...
        'isAbsolute',  0, ...
        'isAvgTime',   0, ...
        'isAvgRow',    0, ...
        'isAvgFreq',   0, ...
        'isMatchRows', 0, ...
        'Dim',         2, ...  % Dimension on which to concatenate (0,1,2,3,4)
        'isBadZero',   0);     % If 1, all the bad values are set to zero
    OPTIONS = struct_copy_fields(OPTIONS, defOPTIONS, 0);

    % ===== CHECK INPUTS =====
    % Types of files in input
    inFileType = sInputs(1).FileType;
    % Is time modified with these options
    isTimeModified = ~OPTIONS.isAvgTime && ((OPTIONS.Dim == 1) || (length(sInputs) == 1));
    % Forbid concatenating in dimension 1 multiple source files
    if strcmpi(inFileType, 'results') && (OPTIONS.Dim == 1) && (length(sInputs) > 1) && isempty(OPTIONS.ScoutSel)
        bst_report('Error', sProcess, sInputs(1), 'Cannot concatenate full source files in dimension 1.');
        return;
    end
    % Do not accept PAC files
    if ~isempty(strfind(sInputs(1).FileName, '_pac'))|| ~isempty(strfind(sInputs(1).FileName, '_dpac'))
        bst_report('Error', sProcess, sInputs(1), 'PAC files are not supported yet. Ask on the forum if you need it.');
        return;
    end

    % ===== PREPARE HISTORY =====
    % History
    strHistory = ['Extract data: ' FormatComment(sProcess)];
    if OPTIONS.isAvgRow
        strHistory = [strHistory ' AvgRow'];
    end
    if OPTIONS.isAvgTime
        strHistory = [strHistory ' AvgTime'];
    end
    if OPTIONS.isAvgFreq
        strHistory = [strHistory ' AvgFreq'];
    end
    if (length(sInputs) > 1)
        strHistory = [strHistory sprintf(' (%d files)', length(sInputs))];
    end
    % History: Add all the file names
    historyMat.History = [];
    historyMat = bst_history('add', historyMat, 'extract', strHistory);
    if ismember(OPTIONS.Dim, [1 2 3 4])
        historyMat.History = repmat(historyMat.History, length(sInputs) + 1, 1);
        historyMat.History(2:end,3) = cellfun(@(c)cat(2, ' - File: ', c), {sInputs.FileName}, 'UniformOutput', 0);
    end

    % ===== GET COMMON ROWS =====
    % Get the list of common rows (only the files with named signals: matrix and timefreq)
    if OPTIONS.isMatchRows && ismember(inFileType, {'matrix','timefreq','pmatrix','ptimefreq'})
        bst_progress('text', 'Reading the signal names...');
        % Cannot use this option in combination with connectivity files
        if ~isempty(strfind(sInputs(1).FileName, '_connect1')) || ~isempty(strfind(sInputs(1).FileName, '_connectn'))
            bst_report('Error', sProcess, [], 'EXTRACT> Cannot use the option "Match signals" on connectivity files.');
            return;
        end
        % Identify the list of all the rows in all the files
        [DestRowNames, AllRowNames, iRowsSrc, iRowsDest, Messages] = process_stdrow('GetUniformRows', {sInputs.FileName}, 'all');
        if ~isempty(Messages)
            bst_report('Warning', sProcess, [], Messages);
        end
        if isempty(DestRowNames)
            bst_report('Warning', sProcess, [], [...
                'EXTRACT> Could not find a common list of rows: Trying to concatenate directly the matrices...' 10 ...
                'EXTRACT> To avoid this warning, uncheck the option "Match signals" in the process options.']);
        end
    else
        DestRowNames = [];
    end

    % ===== LOOP ON FILES =====
    DisplayUnits = [];
    for iInput = 1:length(sInputs)
        bst_progress('text', sprintf('Reading input files... [%d/%d]', iInput, length(sInputs)));

        % === LOAD FILE ===
        isAbsApplied = 0;
        % Extract scout
        if strcmpi(sInputs(1).FileType, 'results') && ~isempty(OPTIONS.ScoutSel)
            % Options for LoadInputFile()
            LoadOptions.LoadFull    = 1;    % Load full source results
            LoadOptions.IgnoreBad   = 1;    % Do not read bad segments
            LoadOptions.ProcessName = func2str(sProcess.Function);
            LoadOptions.TargetFunc  = OPTIONS.ScoutFunc;
            LoadOptions.isNorm      = OPTIONS.isAbsolute;
            % Load reference signal
            sLoaded = bst_process('LoadInputFile', sInputs(iInput).FileName, OPTIONS.ScoutSel, OPTIONS.TimeWindow, LoadOptions);
            if isempty(sLoaded) || isempty(sLoaded.Data)
                bst_report('Error', sProcess, sInputs(iInput), 'Could not extract scouts.');
                return;
            end
            % If there is only one time point: collapse the data that was read
            if ~isempty(OPTIONS.TimeWindow) && (OPTIONS.TimeWindow(1) == OPTIONS.TimeWindow(2))
                sLoaded.Data = sLoaded.Data(:,1,:,:);
                sLoaded.Time = OPTIONS.TimeWindow(1);
            end
            % Add components labels
            if (sLoaded.nComponents == 3)
                sLoaded.RowNames = sLoaded.RowNames(:)';
                sLoaded.RowNames = [cellfun(@(c)cat(2,c,'.1'), sLoaded.RowNames, 'UniformOutput', 0); ...
                                    cellfun(@(c)cat(2,c,'.2'), sLoaded.RowNames, 'UniformOutput', 0); ...
                                    cellfun(@(c)cat(2,c,'.3'), sLoaded.RowNames, 'UniformOutput', 0)];
                sLoaded.RowNames = sLoaded.RowNames(:);
            end
            % Convert to matrix structure
            FileMat = db_template('matrixmat');
            FileMat.Value       = sLoaded.Data;
            FileMat.Description = sLoaded.RowNames;
            FileMat.Comment     = sLoaded.Comment;
            FileMat.Time        = sLoaded.Time;
            FileMat.nAvg        = sLoaded.nAvg;
            FileMat.Leff        = sLoaded.Leff;
            FileMat.SurfaceFile = sLoaded.SurfaceFile;
            FileMat.Atlas       = sLoaded.Atlas;
            FileMat.DisplayUnits= sLoaded.DisplayUnits;
            % Interpret as matrix file in the rest of the function
            inFileType = 'matrix';
            MatValues  = FileMat.Value;
            %isAbsApplied = 1;
        % Regular file
        else
            % Load file
            [FileMat, matName] = in_bst(sInputs(iInput).FileName, OPTIONS.TimeWindow);
            % File could not be loaded
            if isempty(FileMat) || isempty(FileMat.(matName))
                bst_report('Error', sProcess, sInputs(iInput), 'Could not load anything from the input file. Check the requested time window.');
                return;
            end
            % Unconstrained source files: compute norm
            if strcmpi(inFileType, 'results') && OPTIONS.isAbsolute && (FileMat.nComponents > 1)
                FileMat = process_source_flat('Compute', FileMat, 'rms');
                isAbsApplied = 1;
            end
            % Get values to extract
            MatValues = FileMat.(matName);
        end
        % Make sure that units are the same for all the files
        if isfield(FileMat, 'DisplayUnits') && ~isempty(FileMat.DisplayUnits)
            if isempty(DisplayUnits)
                DisplayUnits = FileMat.DisplayUnits;
            elseif ~strcmpi(FileMat.DisplayUnits, DisplayUnits)
                DisplayUnits = 'Mixed units';
            end
        end

        % Get the signals descriptions
        switch (inFileType)
            case 'data'
                % Load channel file (only if new one)
                if ~isequal(ChannelFile, sInputs(iInput).ChannelFile)
                    ChannelFile = sInputs(iInput).ChannelFile;
                    ChannelMat = in_bst_channel(ChannelFile);
                end
                % Select sensors
                if ~isempty(OPTIONS.SensorTypes)
                    % Find selected channels
                    iChannels = channel_find(ChannelMat.Channel, OPTIONS.SensorTypes);
                    if isempty(iChannels)
                        bst_report('Error', sProcess, sInputs(iInput), 'Could not load anything from the input file. Check the sensor selection.');
                        return;
                    end
                    % Keep only selected channels
                    MatValues = MatValues(iChannels,:,:);
                    Description = {ChannelMat.Channel(iChannels).Name}';
                    FileMat.ChannelFlag = FileMat.ChannelFlag(iChannels);
                else
                    Description = {ChannelMat.Channel.Name}';
                end
                % Set the bad values to zero
                if OPTIONS.isBadZero
                    MatValues(FileMat.ChannelFlag == -1, :) = 0;
                end
                % Report good/bad channels
                if isempty(ChannelFlag)
                    ChannelFlag = FileMat.ChannelFlag;
                % When setting the bad channels to zero: set as good all the channels that are good in at least one file
                elseif OPTIONS.isBadZero
                    ChannelFlag(FileMat.ChannelFlag == 1) = 1;
                % Regular case: Set as bad all the channels that are bad in at least one file
                else
                    ChannelFlag(FileMat.ChannelFlag == -1) = -1;
                end
                Freqs = [];
                TFmask = [];
            case 'results'
                Description = [];
                Freqs = [];
                TFmask = [];
            case 'timefreq'
                Description = FileMat.RowNames;
                % Get file frequency vector
                if iscell(FileMat.Freqs)
                    BandBounds = process_tf_bands('GetBounds', FileMat.Freqs);
                    FreqVector = mean(BandBounds,2);
                else
                    FreqVector = FileMat.Freqs;
                end
                % Rounds the frequency vector, to have the same level of precision as the process (3 significant digits)
                FreqVector = round(FreqVector * 1000) / 1000;
                % Get TFmask
                if isfield(FileMat, 'TFmask') && ~isempty(FileMat.TFmask) % && ((length(sInputs) == 1) || (OPTIONS.Dim == 0))
                    TFmask = FileMat.TFmask;
                else
                    TFmask = [];
                end
                % Keep only selected frequencies
                if ~isempty(OPTIONS.FreqRange) && ~isempty(FreqVector) && ~isequal(FreqVector, 0)
                    iFreqs = find((FreqVector >= OPTIONS.FreqRange(1)) & (FreqVector <= OPTIONS.FreqRange(2)));
                    if isempty(iFreqs)
                        bst_report('Error', sProcess, sInputs(iInput), 'Invalid frequency range.');
                        return;
                    end
                    MatValues = MatValues(:,:,iFreqs);
                    % Keep only the selected frequencies
                    if iscell(FileMat.Freqs)
                        Freqs = FileMat.Freqs(iFreqs,:);
                    else
                        Freqs = FileMat.Freqs(iFreqs);
                    end
                    FreqVector = FreqVector(iFreqs);
                    % Report selection on TFmask
                    if ~isempty(TFmask)
                        TFmask = TFmask(iFreqs,:);
                    end
                else
                    Freqs = FileMat.Freqs;
                end
            case 'matrix'
                Description = FileMat.Description;
                Freqs = [];
                TFmask = [];
        end

        % === KEEP COMMON ROWS ===
        if ~isempty(DestRowNames) && ~isempty(Description) && ~isequal(Description, DestRowNames)
            tmpValues = zeros(length(DestRowNames), size(MatValues,2), size(MatValues,3));
            tmpValues(iRowsDest{iInput},:,:) = MatValues(iRowsSrc{iInput},:,:);
            MatValues = tmpValues;
            Description = DestRowNames;
        end

        % === SELECT ROWS ===
        if ~isempty(OPTIONS.Rows)
            % Cannot do any sensor selection on connectivity NxN results
            if isfield(FileMat, 'RefRowNames') && (length(FileMat.RefRowNames) > 1)
                bst_report('Error', sProcess, sInputs(iInput), 'Cannot select rows for connectivity [NxN] results.');
                return;
            end
            % Cannot select signals in sources results
            if ~isempty(Description) && isnumeric(Description)
                bst_report('Error', sProcess, sInputs(iInput), 'Cannot select rows in source maps.');
                return;
            end
            % Try to find row names in the file
            SelRowNames = strtrim(str_split(OPTIONS.Rows, ',;'));
            iRows = find(ismember(Description, SelRowNames));
            % Else, try selecting by index number
            if isempty(iRows)
                iRows = str2num(OPTIONS.Rows);
                iRows((iRows < 1) | (iRows > length(Description))) = [];
            end
            % Else, try selecting by channel type (for time-frequency computed on channels data)
            if isempty(iRows) && strcmpi(inFileType, 'timefreq') && isequal(FileMat.DataType, 'data') && ~isempty(sInputs(iInput).ChannelFile)
                % Load channel file (only if new one)
                if ~isequal(ChannelFile, sInputs(iInput).ChannelFile)
                    ChannelFile = sInputs(iInput).ChannelFile;
                    ChannelMat = in_bst_channel(ChannelFile);
                end
                % Select sensors
                iChanTF = channel_find(ChannelMat.Channel, OPTIONS.Rows);
                % If something was found: try finding the corresponding channel names in the TF file
                if ~isempty(iChanTF)
                    iRows = find(ismember(Description, {ChannelMat.Channel(iChanTF).Name}));
                end
            end
            % Invalid indices selected
            if (any(iRows > size(MatValues,1)))
                bst_report('Error', sProcess, sInputs(iInput), 'Invalid signals selection.');
                return;
            end
            % If there is something to select
            if ~isempty(iRows)
                MatValues = MatValues(iRows,:,:);
                if ~isempty(Description)
                    Description = Description(iRows);
                end
            end
        else
            iRows = [];
        end
        % Error: Nothing selected
        if isempty(MatValues)
            bst_report('Error', sProcess, sInputs(iInput), 'Could not load anything from the input file. Please check the signal selection.');
            return;
        end

        % === ABSOLUTE VALUES ===
        if OPTIONS.isAbsolute && ~isAbsApplied
            MatValues = abs(MatValues);
        end
        % === AVERAGE ROWS ===
        if OPTIONS.isAvgRow && (size(MatValues,1) > 1)
            % Cannot average rows on connectivity NxN results
            if isfield(FileMat, 'RefRowNames') && (length(FileMat.RefRowNames) > 1)
                bst_report('Error', sProcess, sInputs(iInput), 'Cannot average rows for connectivity [NxN] results.');
                return;
            end
            % Average rows
            MatValues = mean(MatValues, 1);
            Description = {'AVG'};
        end
        % === AVERAGE TIME ===
        if OPTIONS.isAvgTime && (size(MatValues,2) > 1)
            % Compute average in time
            MatValues = mean(MatValues, 2);
            % If we are concatenating multiple files in time dimension: keep only one time point for each file
            if (OPTIONS.Dim == 2) && (length(sInputs) > 1)
                FileMat.Time = FileMat.Time(1);
            % If we are not concatenating in time: duplicate the matrix in time so we keep a time segment
            else
                MatValues = [MatValues, MatValues];
                FileMat.Time = [FileMat.Time(1), FileMat.Time(end)];
            end
            % Discard edge effects map
            TFmask = [];
        end
        % === AVERAGE FREQUENCY ===
        if OPTIONS.isAvgFreq && (size(MatValues,3) > 1) && ~isempty(Freqs)
            MatValues = mean(MatValues, 3);
            Freqs = {'AVG', [num2str(FreqVector(1)), ', ' num2str(FreqVector(end))], 'mean'};
            % Edge effects map: good only if good in all the frequencies
            if ~isempty(TFmask)
                TFmask = all(TFmask, 1);
            end
        end
        % === TFMASK ===
        if ~isempty(TFmask)
            if isempty(OutTFmask)
                OutTFmask = TFmask;
            elseif isequal(size(OutTFmask), size(TFmask))
                OutTFmask = OutTFmask & TFmask;
            end
        end

        % === STORE SIGNALS ===
        % Add to load of loaded blocks
        OutValue{iInput} = MatValues;
        % Concatenating in dimension 1: Add the file name
        if (OPTIONS.Dim == 1) && ~isempty(Description)
            % Sources maps: Cannot concatenate in dimension 1
            if ~iscell(Description)
                % If there are more than one file: error
                if (length(sInputs) > 1)
                    bst_report('Error', sProcess, sInputs, 'Cannot concatenate source files in dimension #1.');
                    return;
                end
            % Other file types
            else
                Description = cellfun(@(c)cat(2, c, ' @ ', sInputs(iInput).FileName), Description, 'UniformOutput', 0);
            end
        % No concatenation: Update the modified signals list
        elseif (OPTIONS.Dim == 0) && ~isempty(Description)
            switch (inFileType)
                case 'data',      FileMat.Description = Description;
                case 'timefreq',  FileMat.RowNames    = Description;
                case 'matrix',    FileMat.Description = Description;
            end
        end
        % First signal
        if (iInput == 1)
            % Save signal names
            OutNames = Description(:);
            % Get time definition
            if (length(FileMat.Time) == 1)
                if (OPTIONS.Dim == 2) && (length(sInputs) > 1)
                    sfreq = 1;
                    tstart = 1;
                else
                    sfreq = 1000;
                    tstart = FileMat.Time(1);
                end
            elseif isfield(FileMat, 'TimeBands') && ~isempty(FileMat.TimeBands)
                sfreq = 1;
                tstart = 1;
                FileMat.TimeBands = [];
            else
                sfreq = 1/(FileMat.Time(2) - FileMat.Time(1));
                tstart = FileMat.Time(1);
            end
            % Save first file in memory: will be used as a base (or store all the files if not need of concatenation)
            if (OPTIONS.Dim == 0)
                LoadedMat = cell(1, length(sInputs));
            else
                LoadedMat = cell(1,1);
            end
            LoadedMat{1} = FileMat;
        else
            % Make sure the dimensions are ok
            switch (OPTIONS.Dim)
                case 1
                    isDimOk = (size(OutValue{1},2) == size(OutValue{iInput},2)) && (size(OutValue{1},3) == size(OutValue{iInput},3));
                    % Concatenate names of signals
                    OutNames = cat(1, OutNames, Description(:));
                case 2
                    isDimOk = (size(OutValue{1},1) == size(OutValue{iInput},1)) && (size(OutValue{1},3) == size(OutValue{iInput},3));
                case 3
                    error('Dimension 3 is not supported for concatenation.');
                case 4
                    isDimOk = (size(OutValue{1},1) == size(OutValue{iInput},1)) && (size(OutValue{1},2) == size(OutValue{iInput},2)) && (size(OutValue{1},3) == size(OutValue{iInput},3));
                case 0   % No check, just stacking all the loaded files
                    isDimOk = 1;
                    LoadedMat{iInput} = FileMat;
                otherwise
                    error('Concatenation dimension is not specified.');
            end
            % Error management
            if ~isDimOk
                bst_report('Error', sProcess, sInputs(iInput), sprintf('The dimensions of file #%d do not match file #1.', iInput));
                return;
            end
        end

        % === EVENTS ===
        if isfield(FileMat, 'Events') && ~isempty(FileMat.Events) && isTimeModified
            if isempty(OutEvents)
                OutEvents = FileMat.Events;
            else
                sFile.events = OutEvents;
                sFile.prop.sfreq = sfreq;
                sFile = import_events(sFile, [], FileMat.Events);
                OutEvents = sFile.events;
            end
        end
    end

    % ===== CONCATENATE SIGNALS ====
    bst_progress('text', 'Concatenating files...');
    % Concatenate data matrix
    if (OPTIONS.Dim ~= 0)
        OutValue = {cat(OPTIONS.Dim, OutValue{:})};
    end

    % ===== CREATE STRUCTURE =====
    % Determine the output file format
    switch (inFileType)
        case 'data'
            if OPTIONS.isAvgRow || ((OPTIONS.Dim == 1) && (length(sInputs) > 1)) || (length(ChannelMat.Channel) ~= size(OutValue{1},1))
                newFileType = 'matrix';
            else
                newFileType = inFileType;
            end
        case 'results'
            if OPTIONS.isAvgRow || ((OPTIONS.Dim == 1) && (length(sInputs) > 1)) || ~isempty(iRows)
                newFileType = 'matrix';
            else
                newFileType = inFileType;
            end
        case 'timefreq'
            % Simplify output (not available for statistics)
            % if (length(FreqVector) == 1) && (OPTIONS.Dim ~= 0) && (OPTIONS.Dim ~= 4)
            if (length(FreqVector) == 1) && (OPTIONS.Dim ~= 0) && (OPTIONS.Dim ~= 4) && (length(sInputs) > 1)   %% ADDED FOR EB TO PLOT TOPOGRAPHY OF EXTRACT FREQUENCY OF A SINGLE FILE
                % TF of data or matrix
                if iscell(OutNames)
                    newFileType = 'matrix';
                % TF of sources
                else
                    newFileType = 'results';
                end
            else
                newFileType = inFileType;
            end
        case 'matrix'
            newFileType = inFileType;
    end
    % If the input and output file types are the same: keep the loaded structure
    if isequal(newFileType, inFileType)
        newMat = LoadedMat;
    % Else: Create a new empty sources structure
    elseif strcmpi(newFileType, 'results')
        newMat = {db_template('resultsmat')};
        % Copy source information
        newMat{1}.HeadModelFile = LoadedMat{1}.HeadModelFile;
        newMat{1}.HeadModelType = LoadedMat{1}.HeadModelType;
        newMat{1}.SurfaceFile   = LoadedMat{1}.SurfaceFile;
        newMat{1}.nAvg          = LoadedMat{1}.nAvg;
        newMat{1}.Leff          = LoadedMat{1}.Leff;
    % Else: Create a new empty matrix structure
    else
        newMat = {db_template('matrixmat')};
        newMat{1}.nAvg = LoadedMat{1}.nAvg;
        newMat{1}.Leff = LoadedMat{1}.Leff;
        if (OPTIONS.Dim == 0)
            newMat = repmat(newMat, 1, length(LoadedMat));
            % Copy channel names
            if strcmpi(inFileType, 'data')
                for i = 1:length(LoadedMat)
                    newMat{i}.Description = LoadedMat{i}.Description;
                    newMat{i}.ChannelFlag = LoadedMat{i}.ChannelFlag;
                    newMat{i}.nAvg        = LoadedMat{i}.nAvg;
                    newMat{i}.Leff        = LoadedMat{i}.Leff;
                end
            end
        end
    end
    % Process all the structures in the same way
    for i = 1:length(LoadedMat)
        % Replicate value in time for display if necessary
        if ismember(newFileType, {'data', 'results', 'matrix'}) && (size(OutValue{i},2) == 1) && (OPTIONS.Dim ~= 0)
            OutValue{i} = [OutValue{i}, OutValue{i}];
        end
        % Copy the values in the appropriate fields
        switch (newFileType)
            case 'data'
                newMat{i}.F           = OutValue{i};
                newMat{i}.Events      = OutEvents;
                newMat{i}.ChannelFlag = ChannelFlag;
            case 'results'
                newMat{i}.ImageGridAmp = OutValue{i};
                % Disconnect from parent file if the time was altered
                if isTimeModified || (length(sInputs) > 1)
                    newMat{i}.DataFile = [];
                end
            case 'timefreq'
                newMat{i}.TF     = OutValue{i};
                newMat{i}.TFmask = OutTFmask;
                newMat{i}.Freqs  = Freqs;
                % If there is not concatenation: do not force the row names to be all the same
                if (OPTIONS.Dim ~= 0)
                    newMat{i}.RowNames = OutNames;
                end
                % Disconnect from parent file if the time was altered
                if (length(sInputs) > 1)
                    newMat{i}.DataFile = [];
                end
                % Change the type of data if necessary
                switch (newMat{i}.DataType)
                    case 'data'
                        if OPTIONS.isAvgRow || ((OPTIONS.Dim == 1) && (length(sInputs) > 1)) || ~isempty(iRows)
                            newMat{i}.DataType = 'matrix';
                        end
                    case 'results'
                        if OPTIONS.isAvgRow || ((OPTIONS.Dim == 1) && (length(sInputs) > 1)) || ~isempty(iRows)
                            newMat{i}.DataType = 'matrix';
                        end
                end
            case 'matrix'
                newMat{i}.Value  = OutValue{i};
                newMat{i}.Events = OutEvents;
                % If there is not concatenation: do not force the row names to be all the same
                if (OPTIONS.Dim ~= 0)
                    newMat{i}.Description = OutNames;
                end
        end
        % Final time vector
        newMat{i}.Time = (0:size(OutValue{i},2)-1) ./ sfreq + tstart;
        % Add units
        newMat{i}.DisplayUnits = DisplayUnits;
        % Add common history field
        newMat{i} = bst_history('add', newMat{i}, historyMat.History);
    end
    % Just return one structure
    if (OPTIONS.Dim ~= 0)
        newMat = newMat{1};
    end
    % Return fieldname
    switch (newFileType)
        case 'data',     matName = 'F';
        case 'results',  matName = 'ImageGridAmp';
        case 'timefreq', matName = 'TF';
        case 'matrix',   matName = 'Value';
    end
end

author: brainstorm-tools


text:
function Nimg = Gscale(img,levels,gsize,sigma);
%
%  Function to generate a gaussian-pyramid for the given input image
%
% Input:
%        img: input image-matrix grayscale
%        levels:  number of levels of the pyramid
%        gsize: size of the gaussian kernel [w h] ([5 5] normally provides a smooth output)
%        sigma:  sigma for gaussian kernel
% Output:
%        Nimg:  is a struct consisting of images from each level
%            :  Nimg.img;
% Usage:
%      im = imread('cameraman.tif');
%      Nimg = Gscale(im,3,[5 5],1.6);
%      i = 2; %select a level
%      figure; imshow(Nimg(i).img);
%
% Author: Pranam Janney                     Date: 24th July 2006 15:39
% Email: pranamjanney@yahoo.com
%
% Revised Version 1.0.1                     Date: 04th August 2006, 10:50
%

%guassian filter  with a sigma=1.6
g = fspecial('gaussian',gsize,sigma);

%pyramid
for i = 1:levels
    if i == 1
        im = imfilter(img,g,'conv');
        Nimg(i).img = im;
    else
        %perform guassian filtering
        im = imfilter(Nimg(i-1).img,g,'conv');
        %perform downsampling (horizontal)
        im1 = im(:,1:2:size(Nimg(i-1).img,2));
        %vertical
        im2 = im1(1:2:size(Nimg(i-1).img,1),:);
        %store it in a struct format
        Nimg(i).img = im2;
    end
end

%End

repo: ekfmonoslam


text:
function [detvtx,pth] = SStat_mass_FDR2(pval,maskvtx,rate)
% [detvtx,pth] = SStat_mass_FDR2(pval,maskvtx,rate)
%
% Two-stage FDR approach to achieve tighter control of the FDR. This
% procedure is more powerful than the original FDR procedure implemented in
% SStat_mass_FDR.
%
% Input
% pval: P-values.
% maskvtx: Mask's vertices (1-based). Default [] (all vertices included).
% rate: Expected FDR (between 0 and 1).
%
% Output
% detvtx: Detected vertices (1-based).
% pth: FDR threshold.
%
% $Revision: 1.1 $  $Date: 2015/01/06 17:03:59 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:03:59 $
%    $Revision: 1.1 $
% References: Benjamini, Y., Krieger, A.M., Yekutieli, D. (2006). Adaptive
% linear step-up procedures that control the false discovery rate.
% Biometrika, 93, 491-507.
%
if nargin < 3
    rate = 0.05;
    if nargin < 2
        maskvtx = [];
    end
end;
nv0 = length(pval);
if isempty(maskvtx)
   maskvtx = 1:nv0;
end;
p = pval(maskvtx);
nv = length(p);
%% First stage (m0 estimation)
q0 = rate/(1+rate);
pth0 = SStat_mass_FDR(p,q0);
detv0 = maskvtx(p <= pth0);
ndetv0 = length(detv0);
m0 = nv-ndetv0;
%% Second stage
if (ndetv0 ~= 0) && (ndetv0 ~= nv)
    pth = SStat_mass_FDR(p,q0*nv/m0);
    detvtx = maskvtx(p <= pth);
else
    detvtx = detv0;
    pth = pth0;
end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [grade,X,SR,improve,AppSet] = LocalSearch1(Problem,X,SR,improve,AppSet)
% Local Search 1 of MTS

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    if ~improve
        SR = SR / 2;
        if all(SR<1e-8)
            SR = (Problem.upper-Problem.lower).*(rand(1,length(SR))/10+0.4);
        end
    end
    improve = false;
    grade   = 0;
    for i = randperm(length(SR))
        old_X  = X;
        dec    = X.dec;
        dec(i) = dec(i) + SR(i)*(rand*2-1);
        X      = Problem.Evaluation(dec);
        [grade,improve,AppSet] = Grading(X,old_X,grade,improve,AppSet,Problem.N);
        if all(old_X.obj<=X.obj)
            dec    = old_X.dec;
            dec(i) = dec(i) - 0.5*SR(i)*(rand*2-1);
            X      = Problem.Evaluation(dec);
            [grade,improve,AppSet] = Grading(X,old_X,grade,improve,AppSet,Problem.N);
            if all(old_X.obj<=X.obj)
                X = old_X;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% MP_GEO_READ_NURBS: create a multipatch geometry structure and the interfaces structure from a file in a format similar to the one of the NURBS toolbox.
%
% [geometry, boundaries, interfaces, subdomains] = mp_geo_read_nurbs (filename)
%
% INPUT :
%
%  filename: the name of the file to be read
%
% OUTPUT:
%
%   geometry:   an array of structures that contain the following fields
%               map:     a function handle to evaluate the parametrization
%               map_der: a function handle to evaluate the derivatives of the parametrization
%               nurbs:   a structure compatible with the NURBS toolbox
%
%   boundaries: an array of structures that contain the following fields
%               nsides:  number of faces conforming each boundary
%               patches: patches to which these faces belong
%               faces:   local numbering of the faces in their patch
%
%   interfaces: an array of structures that contain the following fields
%      2D Case
%               ref:    reference number of each interface
%               patch1: the index of the first coinciding patch, and the index of its interface boundary edge
%               patch2: the index of the second coinciding patch, and the index of its interface boundary edge
%               ornt:   a flag telling if the direction on patch1 matches the direction on patch2
%      3D Case
%               ref:    reference number of each interface
%               patch1: the index of the first coinciding patch, and the index of its interface boundary face
%               patch2: the index of the second coinciding patch, and the index of its interface boundary face
%               flag:   a flag telling if the u- and v- parameter directions on the first face coincide with the u- and v- directions on the second face
%               ornt1:  a flag telling if the u- direction on the first face matches the direction on patch2
%               ornt2:  a flag telling if the v- direction on the first face matches the direction on patch2
%
%   subdomains: an array of structures that contains the following fields
%               name:    the name of the subdomain
%               patches: indices of the patches belonging to the subdomain
%
% The interface information is based on the paper
% [1] T. Dokken, E. Quak, V. Skytt, Requirements from Isogeometric Analysis for Changes in Product Design Ontologies. Proceedings of the Focus K3D Conference on Semantic 3D Media and Content, Sophia Antipolis (France), 2010.
%
% An explanation of the file format can be found in the file
%  geopdes_multipatch/doc/geo_specs_mp_v10.txt
%
% Copyright (C) 2009 Carlo de Falco
% Copyright (C) 2010, 2011, 2015 Rafael Vazquez
%
% 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [geom, boundaries, interfaces, subdomains] = mp_geo_read_nurbs (filename)

  fid = fopen (filename, 'r');
  if (fid < 0)
    error ('mp_geo_read_nurbs: cannot open file %s', filename);
  end

  line = fgetl (fid);
  while (line ~= -1)
    if (line(1) ~= '#')
      vec = str2num(line);
      ndim = vec(1);
      rdim = vec(2);
      npatches = vec(3);
      if (numel (vec) == 4)
        nintrfc = vec(4);
        nsubd = 0;
        if (nargout == 4)
          warning ('mp_geo_read_nurbs: no subdomain information is given in the file. The patch number will be used')
          for ii = 1:npatches
            subdomains(ii).name = sprintf ('Patch %i', ii);
            subdomains(ii).patches = ii;
          end
        end
      elseif (numel (vec) == 5)
        nintrfc = vec(4);
        nsubd = vec(5);
      else
        nintrfc = 0;
        interfaces = [];
        nsubd = 0;
        subdomains.name = 'SUBDOMAIN 1';
        subdomains.patches = 1;
      end
      break
    end
    line = fgetl (fid);
  end

  geom(1:npatches) = struct ('map', [], 'map_der', [], 'nurbs', []);
  for iptc = 1:npatches
    geom(iptc).nurbs.form = 'B-NURBS';
    geom(iptc).nurbs.dim = 4;

    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    if (isempty (str2num (line)))
      geom(iptc).name = line;
      line = fgetl (fid);
      while (line(1) == '#')
        line = fgetl (fid);
      end
    end
    vec = str2num (line);
    geom(iptc).nurbs.order = vec+1;

    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    vec = str2num (line);
    geom(iptc).nurbs.number = vec;

    for idim = 1:ndim
      line = fgetl (fid);
      while (line(1) == '#')
        line = fgetl (fid);
      end
      geom(iptc).nurbs.knots{idim} = str2num (line);
    end

    for idim = 1:rdim
      line = fgetl (fid);
      while (line(1) == '#')
        line = fgetl (fid);
      end
      cp_coord = reshape(str2num (line), geom(iptc).nurbs.number);
      geom(iptc).nurbs.coefs(idim,:,:,:) = cp_coord;
    end
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    geom(iptc).nurbs.coefs(4,:,:,:) = reshape(str2num (line), geom(iptc).nurbs.number);

    geom(iptc).map      = @(PTS) geo_nurbs (geom(iptc).nurbs, PTS, 0);
    geom(iptc).map_der  = @(PTS) geo_nurbs (geom(iptc).nurbs, PTS, 1);
    geom(iptc).map_der2 = @(PTS) geo_nurbs (geom(iptc).nurbs, PTS, 2);
  end

  interfaces = [];
  for intrfc = 1:nintrfc
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    interfaces(intrfc).ref = line;
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    vec = str2num (line);
    interfaces(intrfc).patch1 = vec(1);
    interfaces(intrfc).side1  = vec(2);
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    vec = str2num (line);
    interfaces(intrfc).patch2 = vec(1);
    interfaces(intrfc).side2  = vec(2);
    switch (ndim)
     case 2
       line = fgetl (fid);
       while (line(1) == '#')
         line = fgetl (fid);
       end
       interfaces(intrfc).ornt = str2num (line);
     case 3
       line = fgetl (fid);
       while (line(1) == '#')
         line = fgetl (fid);
       end
       vec = str2num (line);
       interfaces(intrfc).flag  = vec(1);
       interfaces(intrfc).ornt1 = vec(2);
       interfaces(intrfc).ornt2 = vec(3);
    end
  end

  subdomains = [];
  for isub = 1:nsubd
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    subdomains(isub).name = line;
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    subdomains(isub).patches = str2num (line);
  end

  boundaries = [];
  line = fgetl (fid); bnd = 0;
  while (line(1) == '#')
    line = fgetl (fid);
  end
  while (line ~= -1)
    bnd = bnd + 1;
    boundaries(bnd).name = line;
    line = fgetl (fid);
    while (line(1) == '#')
      line = fgetl (fid);
    end
    boundaries(bnd).nsides = str2num (line);
    boundaries(bnd).patches = zeros (boundaries(bnd).nsides, 1);
    boundaries(bnd).faces = zeros (boundaries(bnd).nsides, 1);
    for isides = 1:boundaries(bnd).nsides
      line = fgetl (fid);
      while (line(1) == '#')
        line = fgetl (fid);
      end
      vec = str2num (line);
      boundaries(bnd).patches(isides) = vec(1);
      boundaries(bnd).faces(isides) = vec(2);
    end
    line = fgetl (fid);
    while (isempty(line) || line(1) == '#')
      line = fgetl (fid);
    end
  end

  if (isempty(boundaries) && npatches == 1)
    for iface = 1:2*ndim
      boundaries(iface).name = ['BOUNDARY ' num2str(iface)];
      boundaries(iface).nsides = 1;
      boundaries(iface).patches = 1;
      boundaries(iface).faces = iface;
    end
  elseif (isempty (boundaries))
    warning ('The boundary information has not been assigned.\n')
  end

  if (isempty(subdomains))
    subdomains.name = 'WHOLE DOMAIN';
    subdomains.patches = 1:npatches;
  end

  fclose (fid);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%%*****************************************************************
%% sqlpmisc:
%% unscale and produce infeasibility certificates if appropriate
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [X,y,Z,termcode,resid,reldist,msg] = sqlpmisc(blk,At,C,b,X,y,Z,permZ,param);

   termcode    = param.termcode;
   iter        = param.iter;
   obj         = param.obj;
   relgap      = param.relgap;
   prim_infeas = param.prim_infeas;
   dual_infeas = param.dual_infeas;
   homRd       = param.homRd;
   homrp       = param.homrp;
   AX          = param.AX;
   ZpATynorm   = param.ZpATynorm;
   m0          = param.m0;
   indeprows   = param.indeprows;
   normX0      = param.normX0;
   normZ0      = param.normZ0;
   inftol      = param.inftol;
   maxit       = param.maxit;
   scale_data  = param.scale_data;
   printlevel  = param.printlevel;
%%
   resid = []; reldist = []; msg = [];
   if (scale_data)
      normA = param.normA; normC = param.normC; normb = param.normb;
   else
      normA = 1; normC = 1; normb = 1;
   end
   Anorm = ops(At,'norm'); xnorm = ops(X,'norm'); ynorm = norm(y);
   infeas = max(prim_infeas,dual_infeas);
%%
   if (iter >= maxit)
      termcode = -6;
      msg = 'sqlp stop: maximum number of iterations reached';
      if (printlevel); fprintf('\n  %s',msg); end
   end
   if (termcode <= 0)
      %%
      %% To detect near-infeasibility when the algorithm provides
      %% a "better" certificate of infeasibility than of optimality.
      %%
      err = max(infeas,relgap);
      iflag = 0;
      if (obj(2) > 0)
         if (homRd < 0.1*sqrt(err*inftol))
            iflag = 1;
            msg = sprintf('prim_inf,dual_inf,relgap = %3.2e, %3.2e, %3.2e',...
		  prim_infeas,dual_infeas,relgap);
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = 1;
         end
      end
      if (obj(1) < 0)
         if (homrp < 0.1*sqrt(err*inftol))
            iflag = 1;
            msg = sprintf('prim_inf,dual_inf,relgap = %3.2e, %3.2e, %3.2e',...
                  prim_infeas,dual_infeas,relgap);
            if (printlevel); fprintf('\n  %s',msg); end
            termcode = 2;
         end
      end
      if (iflag == 0)
         if (scale_data == 1)
            X = ops(ops(X,'./',normA),'*',normb);
            y = y*normC;
            Z = ops(ops(Z,'.*',normA),'*',normC);
         end
      end
   end
   if (termcode == 1) & (iter > 3)
      msg = 'sqlp stop: primal problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
      if (scale_data == 1)
         X = ops(X,'./',normA); b = b*normb;
      end
      rby = 1/(b'*y); y = rby*y; Z = ops(Z,'*',rby);
      resid = ZpATynorm * rby;
      reldist = ZpATynorm/(Anorm*ynorm);
   end
   if (termcode == 2) & (iter > 3)
      msg = 'sqlp stop: dual problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
      if (scale_data == 1)
         C = ops(C,'.*',normC);
         Z = ops(Z,'.*',normA);
      end
      tCX = blktrace(blk,C,X);
      X = ops(X,'*',1/(-tCX));
      resid = norm(AX)/(-tCX);
      reldist = norm(AX)/(Anorm*xnorm);
   end
   if (termcode == 3)
      maxblowup = max(ops(X,'norm')/normX0,ops(Z,'norm')/normZ0);
      msg = sprintf('sqlp stop: primal or dual is diverging, %3.1e',maxblowup);
      if (printlevel); fprintf('\n  %s',msg); end
   end
   [X,Z] = unperm(blk,permZ,X,Z);
   if ~isempty(indeprows)
      ytmp = zeros(m0,1);
      ytmp(indeprows) = y;
      y = ytmp;
   end
%%*****************************************************************************
%% unperm: undo the permutations applied in validate.
%%
%% [X,Z] = unperm(blk,permZ,X,Z);
%%
%% undoes the permutation introduced in validate.
%%*****************************************************************************

  function [X,Z] = unperm(blk,permZ,X,Z);
%%
  for p = 1:size(blk,1)
     if (strcmp(blk{p,1},'s') & ~isempty(permZ{p}))
        per = permZ{p};
        X{p} = X{p}(per,per);
        Z{p} = Z{p}(per,per);
     end
  end
%%*****************************************************************************

author: yu-jiang

repo: radpbook

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


text:
% Finding a point that satisfies many linear inequalities
% Section 11.4.1, Boyd & Vandenberghe "Convex Optimization"
% Written for CVX by Almir Mutapcic - 02/18/06
%
% We consider a set of linear inequalities A*x <= b which are
% infeasible. Here A is a matrix in R^(m-by-n) and b belongs
% to R^m. We apply a heuristic to find a point x that violates
% only a small number of inequalities.
%
% We use the sum of infeasibilities heuristic:
%
%   minimize   sum( max( Ax - b ) )
%
% which is equivalent to the following LP (book pg. 580):
%
%   minimize   sum( s )
%       s.t.   Ax <= b + s
%              s >= 0
%
% with variables x in R^n and s in R^m.

% problem dimensions (m inequalities in n-dimensional space)
m = 150;
n = 10;

% fix random number generator so we can repeat the experiment
seed = 0;
randn('state',seed);

% construct infeasible inequalities
A = randn(m,n);
b = randn(m,1);

fprintf(1, ['Starting with an infeasible set of %d inequalities ' ...
            'in %d variables.\n'],m,n);

% sum of infeasibilities heuristic
cvx_begin
   variable x(n)
   minimize( sum( max( A*x - b, 0 ) ) )
cvx_end

% full LP version of the sum of infeasibilities heuristic
% cvx_begin
%   variables x(n) s(m)
%   minimize( sum( s ) )
%   subject to
%     A*x <= b + s;
%     s >= 0;
% cvx_end

% number of satisfied inequalities
nv = length( find( A*x > b ) );
fprintf(1,'\nFound an x that violates %d out of %d inequalities.\n',nv,m);

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%  Copyright (C) 2004 Yokogawa Electric Corporation, All Rights Reserved.
%
%  usage:
%    patient_info = GetMeg160PatientInfoM( fid )
%
%  arguments:
%    fid             : file ID
%
%  return values:
%    patient_info        : structure of patient information.
%                        .patient_id         : patient id
%                        .patient_name       : patient name
%                        .patient_birthday   : patient birthday
%                        .patient_sex        : patient sex
%                        .patient_handed     : patient handed
%
%  confirmation of revision:
%   GetMeg160PatientInfoM( Inf ) will show and return revision of this function.
%


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (jtree)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector).
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% soft    - a cell array of soft/virtual evidence;
%           soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ]
%
% e.g., engine = enter_evidence(engine, ev, 'soft', soft_ev)

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes(:);
N = length(bnet.dag);

engine.evidence = evidence; % store this for marginal_nodes with add_ev option
engine.maximize = 0;

% set default params
exclude = [];
soft_evidence = cell(1,N);

% parse optional params
args = varargin;
nargs = length(args);
for i=1:2:nargs
  switch args{i},
   case 'soft',    soft_evidence = args{i+1};
   case 'maximize', engine.maximize = args{i+1};
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end

onodes = find(~isemptycell(evidence));
hnodes = find(isemptycell(evidence));
pot_type = determine_pot_type(bnet, onodes);
 if strcmp(pot_type, 'cg')
  check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag);
end

if is_mnet(bnet)
  pot = engine.user_pot;
  clqs = engine.nums_ass_to_user_clqs;
else
  % Evaluate CPDs with evidence, and convert to potentials
  pot = cell(1, N);
  for n=1:N
    fam = family(bnet.dag, n);
    e = bnet.equiv_class(n);
    if isempty(bnet.CPD{e})
      error(['must define CPD ' num2str(e)])
    else
      pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), ...
                              evidence);
    end
  end
  clqs = engine.clq_ass_to_node(1:N);
end

% soft evidence
soft_nodes = find(~isemptycell(soft_evidence));
S = length(soft_nodes);
if S > 0
  assert(pot_type == 'd');
  assert(mysubset(soft_nodes, bnet.dnodes));
end
for i=1:S
    n = soft_nodes(i);
    % Modif RD - 2006/12/22
    % Why end+1 it doesn't work for me...
    % replace with n and it is ok
    pot{ n } = dpot( n, ns( n ), soft_evidence{ n } );
end
% Modif RD - 2006/12/22
% But now we have to comment this line to ensure dimension matching
%clqs = [clqs engine.clq_ass_to_node(soft_nodes)];

[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes);
[clpot, seppot] = collect_evidence(engine, clpot, seppot);
[clpot, seppot] = distribute_evidence(engine, clpot, seppot);

C = length(clpot);
ll = zeros(1, C);
for i=1:C
  [clpot{i}, ll(i)] = normalize_pot(clpot{i});
end
loglik = ll(1); % we can extract the likelihood from any clique

engine.clpot = clpot;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt