text:
function [p] = gui_mesh_save(p,mesh_save_command,parent)
% GUI_MESH_OPEN: Load an emse mesh (wfr) file into matlab workspace
%
% Usage: [p] = gui_mesh_save(p,[command],[parent])
%
% p is a structure, generated by 'eeg_toolbox_defaults'
% command is either 'init' or 'save'
% parent is a handle to the gui that calls this gui, useful
% for updating the UserData field of the parent from this gui.
% The p structure is returned to the parent when the parent
% handle is given.
%

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

% Licence:  GNU GPL, no express or implied warranties
% History:  02/02 Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p','var'),
   [p] = eeg_toolbox_defaults;
elseif isempty(p),
   [p] = eeg_toolbox_defaults;
end

if ~exist('mesh_save_command','var'), mesh_save_command = 'init'; end

switch mesh_save_command,

case 'init',

    if exist('parent','var'),
        MESHSave = INIT(p,parent);
    else
        MESHSave = INIT(p);
    end

case 'save',

    MESHSave = get(gcbf,'Userdata');
    gui_mesh_save(MESHSave.p,'init',MESHSave.gui);

    if isequal(get(MESHSave.handles.Bhold,'Value'),1),
       [p] = gui_updateparent(MESHSave,0);
    else
       [p] = gui_updateparent(MESHSave);
        close gcbf;
    end

    fprintf('GUI_MESH_OPEN: Save As not yet implemented\n');

otherwise,

    close gcbf;

end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [MESHSave] = INIT(p,parent),

    % GUI General Parameters

    GUIwidth  = 500;
    GUIheight = 120;

    GUI = figure('Name','MESH File Load','Tag','MESH_OPEN',...
                 'NumberTitle','off',...
                 'MenuBar','none','Position',[1 1 GUIwidth GUIheight]);
    movegui(GUI,'center');

    MESHSave.gui = GUI;

    Font.FontName   = 'Helvetica';
    Font.FontUnits  = 'Pixels';
    Font.FontSize   = 12;
    Font.FontWeight = 'normal';
    Font.FontAngle  = 'normal';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Mesh Data Selection and Parameters

    G.Title_data = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
        'Position',[.01 .75 .17 .2],...
        'TooltipString','Mesh/Tesselation Formats (see, ''help mesh_open'')',...
        'String','Mesh Type:','HorizontalAlignment','left');

    switch p.mesh.type
        case 'BrainStorm',       meshType = 1;
        otherwise,               meshType = 1;
    end
    G.PmeshType = uicontrol('Tag','PmeshType','Parent',GUI,'Style','popupmenu',...
        'Units','Normalized',Font,  ...
        'Position',[.20 .75 .28 .2],...
        'String',{'BrainStorm'},'Value',meshType,...
        'TooltipString','Mesh/Tesselation Formats (see, ''help mesh_save'')',...
        'Callback',strcat('MESHSave = get(gcbf,''Userdata'');',...
                          'MESHSave.p.mesh.type = popupstr(MESHSave.handles.PmeshType);',...
                          'set(gcbf,''Userdata'',MESHSave); clear MESHSave;'));

    G.Title_path = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
        'Position',[.01 .50 .17 .2],...
        'String','Path','HorizontalAlignment','left');
    G.EmeshPath = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
    	'Position',[.20 .50 .58 .2], 'String',p.mesh.path,...
        'Callback',strcat('MESHSave = get(gcbf,''Userdata'');',...
                          'MESHSave.p.mesh.path = get(MESHSave.handles.EmeshPath,''String'');',...
                          'set(gcbf,''Userdata'',MESHSave); clear MESHSave;'));

    G.Title_file = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
        'Position',[.01 .25 .17 .2],...
        'String','File','HorizontalAlignment','left');
    G.EmeshFile = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
    	'Position',[.20 .25 .58 .2], 'String',p.mesh.file,...
        'Callback',strcat('MESHSave = get(gcbf,''Userdata'');',...
                          'MESHSave.p.mesh.file = get(MESHSave.handles.EmeshFile,''String'');',...
                          'set(gcbf,''Userdata'',MESHSave); clear MESHSave;'));

    Font.FontWeight = 'bold';

    % BROWSE: Look for the data
    browsecommand = strcat('MESHSave = get(gcbf,''Userdata'');',...
        'cd(MESHSave.p.mesh.path);',...
        '[file, path] = uigetfile(',...
        '{''*.dat;*.txt'', ''ASCII Tesselations (*.dat,*.txt)'';', ...
        ' ''*.asc;*.tri'', ''FreeSurfer Tesselations (*.asc,*.tri)'';', ...
        ' ''*.wfr'', ''EMSE Tesselations (*.wfr)'';', ...
        ' ''*.mat'', ''BrainStorm Tesselations (*.mat)'';', ...
        ' ''*.*'',   ''All Files (*.*)''},', ...
        '''Select Mesh/Tesselation File'');',...
        'if ~isequal(path,0), MESHSave.p.mesh.path = path; end;',...
        'if ~isequal(file,0), MESHSave.p.mesh.file = file; end;',...
        'set(MESHSave.handles.EmeshPath,''String'',MESHSave.p.mesh.path);',...
        'set(MESHSave.handles.EmeshFile,''String'',MESHSave.p.mesh.file);',...
        'if ~isempty(findstr(file,''.dat'')) | ~isempty(findstr(file,''.txt'')), ',...
            'set(MESHSave.handles.PmeshType,''Value'',1); ',...
            'MESHSave.p.mesh.type = popupstr(MESHSave.handles.PmeshType);',...
        'elseif ~isempty(findstr(file,''.asc'')) | ~isempty(findstr(file,''.tri'')), ',...
            'set(MESHSave.handles.PmeshType,''Value'',2); ',...
            'MESHSave.p.mesh.type = popupstr(MESHSave.handles.PmeshType);',...
        'elseif findstr(file,''.wfr''), ',...
            'set(MESHSave.handles.PmeshType,''Value'',3); ',...
            'MESHSave.p.mesh.type = popupstr(MESHSave.handles.PmeshType);',...
        'else, ',...
            'set(MESHSave.handles.PmeshType,''Value'',4); ',...
            'MESHSave.p.mesh.type = popupstr(MESHSave.handles.PmeshType);',...
        'end;',...
        'set(gcbf,''Userdata'',MESHSave); clear MESHSave;');
    G.BmeshFile = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized',Font, ...
        'Position',[.01 .01 .17 .2], 'String','BROWSE',...
        'BackgroundColor',[0.8 0.8 0.0],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback', browsecommand );

    % SAVE AS: Load & convert the data!
    G.Bsave = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.40 .01 .18 .2],...
        'String','SAVE AS','BusyAction','queue',...
        'TooltipString','Load and convert the mesh (not implemented yet).',...
        'BackgroundColor',[0.0 0.0 0.5],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback',strcat('MESHSave = get(gcbf,''Userdata'');',...
                          'p = mesh_write(MESHSave.p,''save'');',...
                          'clear MESHSave;'));

    % Cancel
    G.Bcancel = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.80 .01 .18 .2],...
        'String','CANCEL','BusyAction','queue',...
        'TooltipString','Close, do not return the p struct.',...
        'BackgroundColor',[0.75 0.0 0.0],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback','gui_mesh_open([],''cancel'');');

    % Store userdata
    if exist('parent','var'), MESHSave.parent.gui = parent; end
    MESHSave.handles = G;
    MESHSave.p = p;
    set(GUI,'HandleVisibility','callback');
    set(GUI,'Userdata',MESHSave);

return

author: PatternRecognition

repo: OpenBMI


text:
function h = p37_h ( n, x )

%*****************************************************************************80
%
%% P37_H evaluates the Hessian for problem 37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the N by N Hessian matrix.
%
  h = zeros ( n, n );

  arg = - ( x(1) - pi )^2 - ( x(2) - pi )^2;
  dargdx1 = - 2.0 * ( x(1) - pi );
  dargdx2 = - 2.0 * ( x(2) - pi );

  factor = cos ( x(2) ) * ( sin ( x(1) ) - cos ( x(1) ) * dargdx1 );
  dfdx1 = cos ( x(2) ) * ( cos ( x(1) ) + sin ( x(1) ) * dargdx1 + 2.0 * cos ( x(1) ) );
  dfdx2 = - sin ( x(2) ) * ( sin ( x(1) ) - cos ( x(1) ) * dargdx1 );

  h(1,1) = ( dfdx1 + factor * dargdx1 ) * exp ( arg );
  h(1,2) = ( dfdx2 + factor * dargdx2 ) * exp ( arg );

  factor = cos ( x(1) ) * ( sin ( x(2) ) - cos ( x(2) ) * dargdx2 );
  dfdx1 = - sin ( x(1) ) * ( sin ( x(2) ) - cos ( x(2) ) * dargdx2 );
  dfdx2 = cos ( x(1) ) * ( cos ( x(2) ) + sin ( x(2) ) * dargdx2 ...
    + 2.0 * cos ( x(2) ) );

  h(2,1) = ( dfdx1 + factor * dargdx1 ) * exp ( arg );
  h(2,2) = ( dfdx2 + factor * dargdx2 ) * exp ( arg );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F = curl(F)
%CURL  curl of a CHEBFUN2V
%   S = CURL(F) returns the CHEBFUN2 of the curl of F. If F is a CHEBFUN2V with
%   two components then it returns the CHEBFUN2 representing
%         CURL(F) = F(2)_x - F(1)_y,
%   where F = (F(1),F(2)).  If F is a CHEBFUN2V with three components then it
%   returns the CHEBFUN2V representing the 3D curl operation.

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

% Empty check:
if ( isempty( F ) )
    F = chebfun2v;
    return
end

Fc = F.components;

if ( F.nComponents == 2 )  % 2D curl
    F = diff(Fc{2}, 1, 2) - diff(Fc{1}, 1, 1);
else                       % 3D curl
    F = [ diff(Fc{3},1,1) ; -diff(Fc{3},1,2) ;...
          diff(Fc{2},1,2) - diff(Fc{1},1,1) ];
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%krf_polar 'polar transform'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros rf_polar.pane file
%
% Parameters:
% InputFile: i 'Input (rect)', required: 'input image, rectangular coordinates'
% OutputFile: o 'Output (polar)', required: 'output image, log-polar coordinates'
% Double: rmin 'min', default: 1: 'minimum radius'
% Double: rmax 'max', default: 64: 'maximum radius'
% Double: radinc 'sampling increment', default: 1: 'radius sampling increment'
% Double: amin 'min', default: -90: 'minimum angle'
% Double: amax 'max', default: 90: 'maximum angle'
% Double: anginc 'sampling increment', default: 1: 'angle sampling increment'
%
% Example: o = krf_polar(i, {'i','';'o','';'rmin',1;'rmax',64;'radinc',1;'amin',-90;'amax',90;'anginc',1})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% rf_polar - polar transform
%
%  DESCRIPTION
% This is a spacial transform (warping) of an image.  Specifically, we
% map the specified section of an annulus in the input image to a
% rectangle in the output image.  This is done by mapping annulus pixels
% (m,n) to (r, theta) pixels, where r is the radius and theta is the
% angle from m-axis.
%
% The center of the annulus is also specified as the upper left corner
% of the image (0,0), or the center (N/2, N/2).  In both cases, angles
% are measured from the m (horizontal) axis, and positive in the
% clockwise direction (to remain consistant with the left-hand
% coordinate system typically used in image processing).
%
% The dimensions of the output image are set to either the minimum size
% that will retain the resolution of the input image, or that size
% rounded up to the nearest power of 2 depending on the -logradsz and
% -anglesz flags.
%
%
%
%  EXAMPLES
% See the html tutorial in $REGISTER/examples/html/README_FIRST.html
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996 - 1997, University of New Mexico.  All rights reserved.
%

function varargout = krf_polar(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = krf_polar(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'rmin', 1;'rmax', 64;'radinc', 1;'amin', -90;'amax', 90;'anginc', 1};
maxval={0,0,2,2,2,360,360,2};
minval={0,0,2,2,2,-360,-360,2};
istoggle=[0,0,0,0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Double','Double','Double','Double','Double','Double'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'rf_polar"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function demoStarTrails()
%%DEMOSTARTRAILS This demonstrates how to determine the apparent locations
%                of stars in the sky at a particular location and time
%                using the Hipparcos catalog and a low-fidelity atmospheric
%                refraction model. The stars are plotted over a number of
%                times on the same plot such that their motion in the sky
%                becomes apparent. The planets, Sun and Moon are not
%                plotted.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The location of the observer on the surface of the Earth in ellipsoidal
%coordinates.
obsLoc=[-pi/8;pi-pi/8+pi/8;0];

%Properties of the imaging sensor:
pixX=3264;%Pixels wide
pixY=2448;%Pixels high.
pixWidth=1.5e-6;%Width of a pixel in meters
f=4.12e-3;%Focal length, meters.
%The size of the sensor.
sensX=pixX*pixWidth;
sensY=pixY*pixWidth;
%The angular width of the field of view when viewing something at infinity.
radAz=2*atan(sensX/(2*f));%Width in azimuth (radians)
radEl=2*atan(sensY/(2*f));%Width in elevation (radians)

%The half-azimuth and half-elevation of the observer's field of view from
%the center of the field.
locFOVHalfWidth=[radAz/2;radEl/2];

%It is assumed that the observer's true orientation is a rotation 60
%degrees about the local ellipsoidal ENU x-axis-up followed by a rotation
%of 15 degrees about the rotated z-axis. So, we will now build the
%ENU2Local rotation matrix.
R1=axisAng2RotMat([1;0;0],60*pi/180);
R2=axisAng2RotMat([0;0;1],15*pi/180);
ENU2LocRotMat=R2*R1;

disp('Loading the star catalog.')
%Load stars with a magnitude of at most 6 (the visible region). The stars
%are formatted such that the data file can be directly passed to the
%changeEpoch function, which allows the data to be used at different times.
%The catalog is provided at epoch 1991.25.
dataFile=getHipparcosCat([],6,false,true);

%Change the epoch from 1995.25 to J2000.0
dataFile=changeEpoch(dataFile);

%Set the initial date of display. Minutes will be advanced to demonstrate
%the motion of the stars.
year=2014;
month=1;
day=1;
hour=0;
second=0;

disp('Plotting the stars at different times in the local coordinate')
disp('system and field of view of the observer')
figure(1)
hold on
for minute=0:5:20
    [Jul1,Jul2]=Cal2UTC(year,month,day,hour,minute,second);
    %This returns the apparent star locations in the local ENU coordinate
    %system, where zSpherENU is the location in spherical azimuth and
    %elevation in ENU coordinate and uObs is unit vectors in local ENU
    %coordinates.
    [zSpherENU,uStarsENU]=starCat2Obs(dataFile,Jul1,Jul2,obsLoc);

    %We have to determine which stars are in the observer's field of view.
    %First, we need to rotate the apparent locations into the observer's
    %local coordinate system
    uStarLocal=ENU2LocRotMat*uStarsENU;

    %Next, obtain the azimuth and elevation of the stars in the local
    %coordinate system.
    zStarLocal=Cart2Sphere(uStarLocal);
    %Discard the useless range component.
    zStarLocal=zStarLocal(2:3,:);

    %The thing in the middle of the field of view is at local azimuth and
    %elevation (0,0). We need to figure out which stars and celestial
    %bodies are in the true field of view and get rid of the others.
    sel=zStarLocal(1,:)<locFOVHalfWidth(1);
    sel=sel&(zStarLocal(1,:)>-locFOVHalfWidth(1));
    sel=sel&(zStarLocal(2,:)<locFOVHalfWidth(2));
    sel=sel&(zStarLocal(2,:)>-locFOVHalfWidth(2));

    uStarLocal=uStarLocal(:,sel);
    zStarLocal=zStarLocal(:,sel);

    %Plot the stars in terms of azimuth and elevation in degrees.
    scatter(zStarLocal(1,:)*(180/pi),zStarLocal(2,:)*(180/pi),'filled')
end
axis tight
h0=title('Stars Observed Every 5 Minutes for 20 Minutes');
h1=xlabel('Local Azimuth^\circ');
h2=ylabel('Local Elevation^\circ');
set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h0,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% FUNCTION Construct_iMSF
%   Prepares input for iMSF optimization. This function is not intended to
%   directly called by users. See the Logistic_iMSF function.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Lei Yuan, Jiayu Zhou, and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 12, 2012.
%
%% RELATED FUNCTIONS
%   init_opts, Construct_iMSF, Logistic_iMSF, MultiSource_LogisticR

function [A_Set, Y_Set, W, G, ind, PS, A_Set_Complete] = Construct_iMSF(X_Set, Y)
%

X = [];
temp = [];
IndMissing = [];
nF = [];
FeatureTypes = cell(length(X_Set), 1);

for i = 1:length(X_Set)
    temp = X_Set{i};
    FeatureTypes{i} = sprintf('%d', i);

    temp(~isnan(temp(:, 1)), :) = zscore(temp(~isnan(temp(:, 1)), :));

    IndMissing = [IndMissing ~isnan(temp(:, 1))];
    X = [X temp];
    nF = [nF size(temp, 2)];
end

IndAllZero = sum(IndMissing, 2) == 0;
IndMissing = IndMissing(~IndAllZero, :);
numIndMissing = bin2dec(num2str(IndMissing));
% IndMissing denotes what type of feature is missing for each particular
% sample.

X = X(~IndAllZero, :);
Y = Y(~IndAllZero, :);

PS = nchoose(FeatureTypes);
% PS is the collection of all combinations.
ind_valid = true(length(PS), 1);

% The starting index of each feature type.
nF_Start = [0 cumsum(nF)];
nF_Start = nF_Start(1:end - 1) + 1;

% We shall first remove the combinations that have too few samples.
for i = 1:length(PS)
    curTypes = PS{i};
    temp = zeros(1, length(FeatureTypes));

    for j = 1:length(curTypes)
        type = curTypes{j};
        temp(strcmp(type, FeatureTypes)) = 1;
    end

    temp = bin2dec(num2str(temp));

    ind_temp = numIndMissing == temp;
    Y_Cur = Y(ind_temp);

    if length(unique(Y_Cur)) < 2
        % If there is not enough training label (at least 1 pos 1 neg), we
        % remove this combination.
        ind_valid(i) = false;
    end
end

PS = PS(ind_valid);

IndMatrix = zeros(length(PS), sum(nF));

IndStart = 1;

for i = 1:length(PS)
    curTypes = PS{i};

    for j = 1:length(curTypes)
        type = curTypes{j};
        loc_type = find(strcmp(type, FeatureTypes));
        IndCur = nF_Start(loc_type):nF_Start(loc_type) + nF(loc_type) - 1;
        IndMatrix(i, IndCur) = IndStart:IndStart + nF(loc_type) - 1;
        IndStart = IndStart + nF(loc_type);
    end
end

W = [];
G = [];
IndStart = 1;

for i = 1:size(IndMatrix, 2)
    curGroup = IndMatrix(:, i)';
    curGroup(curGroup == 0) = [];
    G = [G curGroup];
    W = [W [IndStart; IndStart + length(curGroup) - 1]];
    IndStart = IndStart + length(curGroup);
end

A_Set = cell(length(PS), 1);
A_Set_Complete = cell(length(PS), 1);
Y_Set = cell(length(PS), 1);
ind = zeros(length(PS) + 1, 1);

for i = 1:length(PS)
    curTypes = PS{i};
    temp = zeros(1, length(FeatureTypes));
    ind_col = false(1, sum(nF));

    for j = 1:length(curTypes)
        type = curTypes{j};
        loc_type = find(strcmp(type, FeatureTypes));
        temp(loc_type) = 1;
        IndCur = nF_Start(loc_type):nF_Start(loc_type) + nF(loc_type) - 1;
        ind_col(IndCur) = true;
    end

    ind(i + 1) = ind(i) + nnz(ind_col);

    temp = bin2dec(num2str(temp));

    ind_temp = numIndMissing == temp;

    A_Set{i} = X(ind_temp, ind_col);
    A_Set_Complete{i} = X(ind_temp, :);
    Y_Set{i} = Y(ind_temp);
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function vis_sticks_waf(img,pred_stickmen,gt_stickmen,evalidx,scoredetail,fname)

% plot predicted sticks, upper body detections and GT sticks
% predicted and GT sticks are only plotted for upper body detections
% that match the GT upper body bounding boxes

clf; imagesc(img); axis equal; axis  off; hold on;
for ridx = 1:length(gt_stickmen)
    if (evalidx(ridx) ~= 0) % upper body detection matches GT bounding box
        % plot detection
        sc = scoredetail;
%         if (sc{ridx}(6) == 0) % head prediction is wrong
        if (sc{ridx}(1) == 0) % torso prediction is wrong
            c = 'y';
        else
            c = 'r';
        end
        b = pred_stickmen(evalidx(ridx)).det;
        rectangle('Pos',[b(1) b(2) b(3)-b(1) b(4)-b(2)],'lineWidth',2,'edgeColor',c);

        % plot number of parts estimated correctly / total number of parts
        text(b(1),b(2),sprintf('%d/%d',sum(sc{ridx}),numel(sc{ridx})),'FontSize',10,'color','k','BackgroundColor','g',...
            'verticalalignment','top','horizontalalignment','left');
        DrawStickman(pred_stickmen(evalidx(ridx)).coor,[]);
    end
end
print(gcf,'-dpng',[fname '.png']);
% draw a subset of GT poses for which evaluation is performed
clf; imagesc(img); axis equal; axis  off; hold on;
for ridx = 1:length(gt_stickmen)
    if (evalidx(ridx) ~= 0) % upper body detection matches GT bounding box
        DrawStickman(gt_stickmen(ridx).coor,[]);
    end
end
print(gcf,'-dpng',[fname '_gt.png']);
unix(['montage ' fname '.png' ' ' fname '_gt.png' ' -geometry +0+0 ' fname '_combo.png']);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
  function ob = Gbspline1(varargin)
%|function ob = Gbspline1(options)
%|
%| Construct Gbspline1 object for 1D B-spline interpolation
%|
%| See Gbspline1_test.m for example usage.
%|
%| in
%|
%| options
%|	'type'	char		'val2coef' convert signal values to coefficients
%|				'
%|				'coef2val' | 'synth' conv
%|	'mask'	size(image)	logical array of object support.
%|	'chat'		verbose printing of debug messages
%|	'psf'		point spread function (aka impulse response)
%|	'type'		type of bspline1:
%|				'conv,same'	usual case (default)
%|				'fft,same'	(periodic end conditions)
%|				todo: allow replicated end conditions!
%|				'imfilter,same'	(requires image toolbox)
%|				'imfilter,circ'	(requires image toolbox)
%|	'imfilter_options'	options to 'imfilter' (if used)
%|
%| out
%|	ob [nd np]	np = sum(mask(:)), so it is already "masked"
%|			nd = prod(size(mask)) for 'conv,same' type
%|
%| Copyright 2011-11-13, Jeff Fessler, University of Michigan

fail 'not done'

if nargin == 1 && streq(mask, 'test'), Gbspline1_test, return, end
if nargin < 1, ir_usage, end

% option defaults
arg.type = 'todo';
arg.mask = mask;
arg.chat = 0;

% options specified by name/value pairs
arg = vararg_pair(arg, varargin);

does_many = false;
switch arg.type
case 'todo'
	arg.odim =
	handle_forw = @(arg,x)
	handle_back = @(arg,y)
	does_many = true;

otherwise
	error 'unknown bspline1 type'
end

ob = fatrix2('idim', arg.odim, 'arg', arg, 'does_many', does_many, ...
	'forw', handle_forw, 'back', handle_back);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ g ] = gsp_design_held(G, param)
%GSP_DESIGN_HELD Create a Held filterbank
%   Usage: g = gsp_design_held( G );
%          g = gsp_design_held( G, param );
%
%   Inputs parameters:
%       G       : Graph structure or lmax
%       param   : Structure of optional parameters
%
%   Outputs parameters:
%       g       : filterbank
%
%   This function create a parseval filterbank of $2$ filters. The low-pass
%   filter is defined by a function $f_l(x)$:
%
%   ..                  /  1                                if x <= a
%   ..        f_l(x) = |   sin(2*pi*mu(val(r2ind)/(8*a)))   if a < x <= 2a
%   ..                  \  0                                if x > 2a
%
%   .. math:: f_{l}=\begin{cases} 1 & \mbox{if }x\leq a\\ \sin\left(2\pi\mu\left(\frac{x}{8a}\right)\right) & \mbox{if }a<x\leq2a\\ 0 & \mbox{if }x>2a \end{cases}
%
%   with
%
%   ..        mu(x) = -1 + 24*x - 144*x^2 + 256*x^3
%
%   .. math:: \mu(x) = -1+24x-144*x^2+256*x^3
%
%   The high pass filter is adaptated to obtain a tight frame.
%
%   This function will compute the maximum eigenvalue of the laplacian. To
%   be more efficient, you can precompute it using::
%
%       G = gsp_estimate_lmax(G);
%
%   Example:::
%
%         G = gsp_sensor(100);
%         G = gsp_estimate_lmax(G);
%         g = gsp_design_held(G);
%         gsp_plot_filter(G,g);
%         [A,B] = gsp_filterbank_bounds(G,g)
%
%   *param* is an optional structure containing the following fields
%
%   * *param.verbose*: verbosity level. 0 no log - 1 display warnings.
%     (default 1)
%   * *param.a*: see equations above for this parameter. Note that the
%     spectrum is scaled between 0 and 2 (default 2/3).
%

% Author: Nathanael Perraudin, David Shuman
% Date  : 21 June 2014
% Testing: test_filter

if nargin < 2
    param = struct;
end

if ~isfield(param,'verbose'), param.verbose = 1; end
if ~isfield(param,'a'), param.a = 2/3; end

if isstruct(G)
    if ~isfield(G,'lmax')
        if param.verbose
            fprintf('GSP_DESIGN_HELD has to compute lmax \n')
        end
        G = gsp_estimate_lmax(G);
    end
   lmax = G.lmax;
else
   lmax = G;
end

a = param.a;

g = cell(2,1);
g{1} = @(x) held(x*(2/lmax),a);
g{2} = @(x) real(sqrt(1-(held(x*(2/lmax),a)).^2));

end

function y = held(val,a)

y = zeros(size(val));

l1 = a;
l2 = 2*a;
mu = @(x) -1+24*x-144*x.^2+256*x.^3;

r1ind = val >= 0     &    val < l1;
r2ind = val >= l1    &    val < l2;
r3ind = val >= l2;

y(r1ind) = 1;
y(r2ind) = sin(2*pi*mu(val(r2ind)/(8*a)));
y(r3ind) = 0;

end

repo: gspbox


text:
function r = model_has_colors(model)
r = isfield(model,'r') && length(model.r) > 1;

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function drawscene(X,C,R,fig,ctypehint,scenetitle,camsId)

% drawscene ... plots calibration points, cameras and their viewing axes
%
% drawscene(X,C,R,fig,ctypehint,scenetitle)
%
% X ............ 4xPOINTS matrix containg POINTS object points
% C ............ 3xCAMS matrix containing the camera centers (in world coord.)
% R ............ 3*CAMSx3 matrix containing camera rotation matrices
%                (needed for drawing the viewing axes)
% fig .......... figure handle (defaults to 1)
% ctypehint .... calibration object type of X (defaults to 'cloud')
% scenetitle ... title of the plot (defaults to '')
% camsIs ....... 1xCAMS vector with cameas Id (default is 1:CAMS

% $Author: svoboda $
% $Revision: 2.0 $
% $Id: drawscene.m,v 2.0 2003/06/19 12:07:03 svoboda Exp $
% $State: Exp $

POINTS = size(X,2);
CAMS   = size(C,2);

if nargin < 7
  camsId = [1:CAMS];
end

if (nargin < 3)
  error('not enough input arguments');
end
if (nargin < 5)
  scenetitle = '';
end
if (nargin < 4)
  ctypehint = 'cloud';
end

figure(fig); clf
title(scenetitle)
grid on
axis equal

% plot camera positions (blue)
drawcloud(C,fig,'b');

% plot calibration object (red)
drawobject(X,ctypehint,fig,'r');

% Mean of all points
centroid = mean(X(1:3,:)');

% plot viewing axes
for i=1:CAMS
  axis_dir = -R(3*i,:); % 3rd row of i-th rotation matrix
  axis_len = 0.6*norm(C(1:3,i)-centroid');
  endpoint = C(1:3,i)+axis_len*axis_dir';
  line([C(1,i),endpoint(1)],[C(2,i),endpoint(2)],[C(3,i),endpoint(3)]);
  text(C(1,i),C(2,i),C(3,i),sprintf('%4d',camsId(i)),'Color','k');
end

author: strawlab

repo: MultiCamSelfCal


text:
clc;
clear;
%path='C:\Users\Administrator\Desktop\mat-cvpr17_v1\ECSSD\';
path='./PR_results/';
dirpath=dir([path '*.mat']);
method1 = cell(length(dirpath),3);

str=['r','r','b','b','c','c','m','m','k','k','y','y','g','g','r','r','m','m','k','k','r','r','r','r','r'];

rr=[];
Fmeasure=[];
aAUC=[];
method2 = cell(length(dirpath),1);
for i=1:length(dirpath)
  load([path dirpath(i).name]);
  Fmeasure=[Fmeasure,mFmeasure];
  aAUC=[aAUC,AUC];                     %AUC 图像修改部分
%  method2(i)=alg_dir{i}(1);
  method2{i}=dirpath(i).name(1:end-4);

end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        figure(5);
        barMsra = [Fmeasure' aAUC'];
        bar(barMsra );
        set( gca, 'xtick', 1:1:length(dirpath) ),

        set( gca ,'xticklabels',  method2 , 'fontsize', 8 );
        legend('Precision','Recall','Fmeasure','AUC');      %AUC 图像修改部分
        %legend('Precision','Recall','Fmeasure');
              grid on;
   %     saveas(  figure(5), [basedir,'SED1_bar.fig']);

% xlabel('Recall');
% ylabel('Precision');
%legend(dirpath(1).name(5:end),dirpath(2).name(5:end),dirpath(3).name(5:end),dirpath(4).name(5:end),dirpath(5).name(5:end),dirpath(6).name(5:end),dirpath(7).name(5:end),dirpath(8).name(5:end),dirpath(9).name(5:end),dirpath(10).name(5:end),dirpath(11).name(5:end),dirpath(12).name(5:end),dirpath(13).name(5:end),dirpath(14).name(5:end),dirpath(15).name(5:end),dirpath(16).name(5:end),dirpath(17).name(5:end),dirpath(18).name(5:end),dirpath(19).name(5:end),dirpath(20).name(5:end),dirpath(21).name(5:end));
%legend(dirpath(1).name(8:end-4),dirpath(2).name(8:end-4),dirpath(3).name(8:end-4),dirpath(4).name(8:end-4),dirpath(5).name(8:end-4),dirpath(6).name(8:end-4),dirpath(7).name(8:end-4),dirpath(8).name(8:end-4),dirpath(9).name(8:end-4),dirpath(10).name(8:end-4),dirpath(11).name(8:end-4),dirpath(12).name(8:end-4),dirpath(13).name(8:end-4),dirpath(14).name(8:end-4));
%legend(dirpath(1).name(1:end-4),dirpath(2).name(1:end-4),dirpath(3).name(1:end-4),dirpath(4).name(1:end-4),dirpath(5).name(1:end-4),dirpath(6).name(1:end-4),dirpath(7).name(1:end-4),dirpath(8).name(1:end-4))%,dirpath(9).name(1:end-4),dirpath(10).name(1:end-4),dirpath(11).name(1:end-4));%,dirpath(12).name(1:end-4),dirpath(13).name(1:end-4));
%legend(str(1,:),str(2,:),str(3,:),str(4,:),str(5,:),str(6,:),str(7,:),str(8,:),str(9,:),str(10,:),str(11,:),str(12,:),str(13,:),str(14,:),str(15,:),str(16,:),str(17,:));
%legend(dirpath(1).name(2:end-4),dirpath(2).name(2:end-4),dirpath(3).name(2
%legend(dirpath(1).name(1:end-4),dirpath(2).name(1:end-4),dirpath(3).name(1:end-4),dirpath(4).name(1:end-4),...
%dirpath(5).name(1:end-4),dirpath(6).name(1:end-4))%,dirpath(7).name(1:end-4),dirpath(8).name(1:end-4),dirpath(9).name(1:end-4),...
% dirpath(10).name(1:end-4),dirpath(11).name(1:end-4),dirpath(12).name(1:end-4),dirpath(13).name(1:end-4),...
% dirpath(14).name(1:end-4),dirpath(15).name(1:end-4),dirpath(16).name(1:end-4),dirpath(17).name(1:end-4),...
% dirpath(18).name(1:end-4),dirpath(19).name(1:end-4),dirpath(20).name(1:end-4),dirpath(21).name(1:end-4))%,...
%dirpath(22).name(1:end-4),dirpath(23).name(1:end-4));
grid on;
% save 300 AUC;

repo: Saliency-Evaluation-Toolbox


text:
function [optExists, value, vargs] = peel_option(vargs, searchValue, searchClass, minPos)
   optExists = false; value = [];
   if minPos > numel(vargs)
      return
   end
   switch searchClass
      case 'char'
         for n = minPos:numel(vargs)
            if ischar(vargs{n})
               if isempty(searchValue) || strcmp(searchValue,vargs{n})
                  optExists = true; value = vargs{n}; vargs(n) = [];
                  return
               end
            end
         end
      otherwise
         for n = minPos: numel(vargs)
            if isa(vargs{n},searchClass)
               if isempty(searchValue) || vargs{n} == searchValue
                  optExists = true; value = vargs{n}; vargs(n) = [];
                  return
               end
            end
         end
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function [rhol,rhoul,El] = ModEuler(h,h_next,q,fr,v)
%% Global Variables
global CFL r_time theta dt dtdx nx %#ok<*NUSED>
global w k nv
global gamma etpfix

%% Compute primitive Variables @ cell center variables
% Build vector U and dU, and F
    % Define:
    L = 1:nx-1; R = 2:nx; % inner nodes / middle points
    % here U = [U1 U2 U3]' is defined for every cell center
    U = q;
        % U1 = Density, U2 = Momentum, U3 = Total Energy

   % we need to review this initial formulation!
    r = q(1,:);
    u = q(2,:)./q(1,:);
    E = q(3,:);
    p = (gamma-1).*(E-r.*u.^2/2);
    H = (E + p)./r;

%% Main Loop
% Compute Roe Averages
    % Velovity 'u'
    u_bar = (sqrt (r(L)).*u(L) + sqrt (r(R)).*u(R)) ...
        ./ (sqrt(r(L))+sqrt(r(R)));
    % Total Entalpy 'H'
    H_bar = (sqrt (r(L)).*H(L) + sqrt (r(R)).*H(R)) ...
        ./ (sqrt(r(L))+sqrt(r(R)));
    % Sound Speed 'a'
    a_bar = sqrt((gamma-1)*(H_bar-0.5*u_bar.^2));

% Compute Delta U's
% dU = U(R)-U(L) at the cell boundaries { x_{i},x_{i+1}, ... }
dU = U(:,R)-U(:,L);

% Compute Fluxes at the cell centers
% F = [F1 F2 F3]
F = [r.*u; r.*u.^2 + p; r.*u.*H];
% Fluxes to the left and right of 'F_ {i+1/2}'
FL = F(:,L); FR = F(:,R);

% Scaled Right Eigenvectors are given by:
k1_bar = [1*ones(1,nx-1); u_bar-a_bar; H_bar-u_bar.*a_bar];
k2_bar = [1*ones(1,nx-1); u_bar      ; 1/2*u_bar.^2      ];
k3_bar = [1*ones(1,nx-1); u_bar+a_bar; H_bar+u_bar.*a_bar];

% compute Roe waves strength alpha_bar
alpha2_bar = (gamma-1)./(a_bar.^2).*(dU (1,:).*(H_bar-u_bar.^2) ...
    + u_bar.*dU(2,:) - dU(3,:));
alpha1_bar = 1./(2*a_bar).*(dU (1,:).*(u_bar+a_bar) ...
    - dU(2,:)-a_bar.*alpha2_bar);
alpha3_bar = dU(1,:)-(alpha1_bar + alpha2_bar);

% Eigenvalues of A_bar are (same as the original A mat)
lambda_bar(1,:) = abs(u_bar - a_bar);
lambda_bar(2,:) = abs(u_bar);
lambda_bar(3,:) = abs(u_bar + a_bar);

% Entropy Fix
boolv1 = lambda_bar(1,:) < etpfix;
lambda_bar(1,:) = 0.5*(etpfix + lambda_bar (1,:).^2/etpfix).*boolv1 ...
    + lambda_bar(1,:).*(1-boolv1);
boolv2 = lambda_bar(3,:) < etpfix;
lambda_bar(3,:) = 0.5*(etpfix + lambda_bar (3,:).^2/etpfix).*boolv2 ...
    + lambda_bar(3,:).*(1-boolv2);

% Conditioning data
alpha1_bar  = repmat(alpha1_bar,3,1);
alpha2_bar  = repmat(alpha2_bar,3,1);
alpha3_bar  = repmat(alpha3_bar,3,1);
lambda1_bar = repmat(lambda_bar(1,:),3,1);
lambda2_bar = repmat(lambda_bar(2,:),3,1);
lambda3_bar = repmat(lambda_bar(3,:),3,1);

% Roe Fluxes
Flux = 0.5*(FL+FR)-0.5*(alpha1_bar.*lambda1_bar.*k1_bar + ...
    alpha2_bar.*lambda2_bar.*k2_bar + ...
    alpha3_bar.*lambda3_bar.*k3_bar );

%% Compute next time step
% compute difference of the integral(v*fl)
[M1,M2,M3] =integrate_vf(k,w,fr,v); M = [M1;M2;M3]; dM = M(:,R) - M(:,L);

% Remap from h(nx,nx) to h(3,nx)
hh = h(1:3,:);
hh_next = h_next(1:3,:);
dhh = hh_next - hh;

% Compute next time step
U_next = zeros(3,nx);
for i = 2 : nx-1
    U_next(:,i) = U(:,i) - dtdx.*(1-hh(:,i)).*(Flux(:,i) - Flux(:,i-1)) - ...
                           dtdx.*(1-hh(:,i)).* dM(i)- ...
                           q(:,i).*dhh(:,i);
end

% BCs
 U_next(:,1)  = U(:,2); % Neumann condition to the left
 U_next(:,nx) = U(:,nx-1); % Neumann condition to the right

% Compute variables of the new time step
rhol = U_next(1,:);     % Density
rhoul = U_next (2,:);   % Momentum
El = U_next(3,:);       % Total Energy

repo: Aero-matlab


text:
function data2num_finviz(finviz_dir)
clearvars -except finviz_dir
addpath('functions');

% Finviz data is a cell matrix (n_tickers x n_vars)
finviz_filename=strcat(finviz_dir,'/finviz.csv');
[finviz_data,stop]=grab_all_finviz_data(finviz_filename);

FINHEAD = set_finviz_labels(); % global index for each column in finviz.csv
%FINHEAD_names = fieldnames(FINHEAD); % cell of field name strings

% Numbers file
% FHEAD is the local index we decide to output (FHEAD_names is its
% fieldnames)
ntickers=size(finviz_data,1);
if(1)

  data = [];
  i=1;

  MarketCap_vec = convert2numvec( finviz_data(:,FINHEAD.MarketCap) ); data = [data; MarketCap_vec']; FHEAD=struct('MarketCap',{i}); i=i+1;
  PE_vec = convert2numvec( finviz_data(:,FINHEAD.PE) ); data = [data; PE_vec']; [FHEAD(:).PE]=i; i=i+1;
  PEForward_vec = convert2numvec( finviz_data(:,FINHEAD.PEForward) ); data = [data; PEForward_vec']; [FHEAD(:).PEForward]=i; i=i+1;
  PEG_vec = convert2numvec( finviz_data(:,FINHEAD.PEG) ); data = [data; PEG_vec']; [FHEAD(:).PEG]=i; i=i+1;

  PS_vec = convert2numvec( finviz_data(:,FINHEAD.PS) ); data = [data; PS_vec']; [FHEAD(:).PS]=i; i=i+1;
  PB_vec = convert2numvec( finviz_data(:,FINHEAD.PB) ); data = [data; PB_vec']; [FHEAD(:).PB]=i; i=i+1;
  PCash_vec = convert2numvec( finviz_data(:,FINHEAD.PCash) ); data = [data; PCash_vec']; [FHEAD(:).PCash]=i; i=i+1;
  PFreeCash_vec = convert2numvec( finviz_data(:,FINHEAD.PFreeCash) ); data = [data; PFreeCash_vec']; [FHEAD(:).PFreeCash]=i; i=i+1;
  DividendYield_vec = convert2numvec( finviz_data(:,FINHEAD.DividendYield) ); data = [data; DividendYield_vec']; [FHEAD(:).DividendYield]=i; i=i+1;
  PayoutRatio_vec = convert2numvec( finviz_data(:,FINHEAD.PayoutRatio) ); data = [data; PayoutRatio_vec']; [FHEAD(:).PayoutRatio]=i; i=i+1;

  EPS_vec = convert2numvec( finviz_data(:,FINHEAD.EPS) ); data = [data; EPS_vec']; [FHEAD(:).EPS]=i; i=i+1;
  EPSGrowthThisYear_vec = convert2numvec( finviz_data(:,FINHEAD.EPSGrowthThisYear) ); data = [data; EPSGrowthThisYear_vec']; [FHEAD(:).EPSGrowthThisYear]=i; i=i+1;
  EPSGrowthNextYear_vec = convert2numvec( finviz_data(:,FINHEAD.EPSGrowthNextYear) ); data = [data; EPSGrowthNextYear_vec']; [FHEAD(:).EPSGrowthNextYear]=i; i=i+1;
  EPSGrowthPast5Years_vec = convert2numvec( finviz_data(:,FINHEAD.EPSGrowthPast5Years) ); data = [data; EPSGrowthPast5Years_vec']; [FHEAD(:).EPSGrowthPast5Years]=i; i=i+1;
  EPSGrowthNext5Years_vec = convert2numvec( finviz_data(:,FINHEAD.EPSGrowthNext5Years) ); data = [data; EPSGrowthNext5Years_vec']; [FHEAD(:).EPSGrowthNext5Years]=i; i=i+1;

  SalesGrowth5Years_vec = convert2numvec( finviz_data(:,FINHEAD.SalesGrowth5Years) ); data = [data; SalesGrowth5Years_vec']; [FHEAD(:).SalesGrowth5Years]=i; i=i+1;
  EPSGrowthLastQuarter_vec = convert2numvec( finviz_data(:,FINHEAD.EPSGrowthLastQuarter) ); data = [data; EPSGrowthLastQuarter_vec']; [FHEAD(:).EPSGrowthLastQuarter]=i; i=i+1;
  SalesGrowthLastQuarter_vec = convert2numvec( finviz_data(:,FINHEAD.SalesGrowthLastQuarter) ); data = [data; SalesGrowthLastQuarter_vec']; [FHEAD(:).SalesGrowthLastQuarter]=i; i=i+1;
  SharesOutstanding_vec = convert2numvec( finviz_data(:,FINHEAD.SharesOutstanding) ); data = [data; SharesOutstanding_vec']; [FHEAD(:).SharesOutstanding]=i; i=i+1;

  Float_vec = convert2numvec( finviz_data(:,FINHEAD.Float) ); data = [data; Float_vec']; [FHEAD(:).Float]=i; i=i+1;
  InsiderOwnership_vec  = convert2numvec( finviz_data(:,FINHEAD.InsiderOwnership) ); data = [data; InsiderOwnership_vec']; [FHEAD(:).InsiderOwnership]=i; i=i+1;
  InsiderTransactions_vec  = convert2numvec( finviz_data(:,FINHEAD.InsiderTransactions) ); data = [data; InsiderTransactions_vec']; [FHEAD(:).InsiderTransactions]=i; i=i+1;
  InstitutionalOwnership_vec = convert2numvec( finviz_data(:,FINHEAD.InstitutionalOwnership) ); data = [data; InstitutionalOwnership_vec']; [FHEAD(:).InstitutionalOwnership]=i; i=i+1;
  InstitutionalTransactions_vec = convert2numvec( finviz_data(:,FINHEAD.InstitutionalTransactions) ); data = [data; InstitutionalTransactions_vec']; [FHEAD(:).InstitutionalTransactions]=i; i=i+1;

  FloatShort_vec  = convert2numvec( finviz_data(:,FINHEAD.FloatShort) ); data = [data; FloatShort_vec']; [FHEAD(:).FloatShort]=i; i=i+1;
  ShortRatio_vec = convert2numvec( finviz_data(:,FINHEAD.ShortRatio) ); data = [data; ShortRatio_vec']; [FHEAD(:).ShortRatio]=i; i=i+1;

  ReturnOnAssets_vec = convert2numvec( finviz_data(:,FINHEAD.ReturnOnAssets) ); data = [data; ReturnOnAssets_vec']; [FHEAD(:).ReturnOnAssets]=i; i=i+1;
  ReturnOnEquity_vec = convert2numvec( finviz_data(:,FINHEAD.ReturnOnEquity) ); data = [data; ReturnOnEquity_vec']; [FHEAD(:).ReturnOnEquity]=i; i=i+1;
  ReturnOnInvestment_vec = convert2numvec( finviz_data(:,FINHEAD.ReturnOnInvestment) ); data = [data; ReturnOnInvestment_vec']; [FHEAD(:).ReturnOnInvestment]=i; i=i+1;
  CurrentRatio_vec = convert2numvec( finviz_data(:,FINHEAD.CurrentRatio) ); data = [data; CurrentRatio_vec']; [FHEAD(:).CurrentRatio]=i; i=i+1;
  QuickRatio_vec = convert2numvec( finviz_data(:,FINHEAD.QuickRatio) ); data = [data; QuickRatio_vec']; [FHEAD(:).QuickRatio]=i; i=i+1;
  LTDebtEquity_vec = convert2numvec( finviz_data(:,FINHEAD.LTDebtEquity) ); data = [data; LTDebtEquity_vec']; [FHEAD(:).LTDebtEquity]=i; i=i+1;
  TotalDebtEquity_vec = convert2numvec( finviz_data(:,FINHEAD.TotalDebtEquity) ); data = [data; TotalDebtEquity_vec']; [FHEAD(:).TotalDebtEquity]=i; i=i+1;

  GrossMargin_vec  = convert2numvec( finviz_data(:,FINHEAD.GrossMargin) ); data = [data; GrossMargin_vec']; [FHEAD(:).GrossMargin]=i; i=i+1;
  OperatingMargin_vec  = convert2numvec( finviz_data(:,FINHEAD.OperatingMargin) ); data = [data; OperatingMargin_vec']; [FHEAD(:).OperatingMargin]=i; i=i+1;
  ProfitMargin_vec  = convert2numvec( finviz_data(:,FINHEAD.ProfitMargin) ); data = [data; ProfitMargin_vec']; [FHEAD(:).ProfitMargin]=i; i=i+1;

  PerformanceWeek_vec  = convert2numvec( finviz_data(:,FINHEAD.PerformanceWeek) ); data = [data; PerformanceWeek_vec']; [FHEAD(:).PerformanceWeek]=i; i=i+1;
  PerformanceMonth_vec  = convert2numvec( finviz_data(:,FINHEAD.PerformanceMonth) ); data = [data; PerformanceMonth_vec']; [FHEAD(:).PerformanceMonth]=i; i=i+1;
  PerformanceQuarter_vec  = convert2numvec( finviz_data(:,FINHEAD.PerformanceQuarter) ); data = [data; PerformanceQuarter_vec']; [FHEAD(:).PerformanceQuarter]=i; i=i+1;
  PerformanceHalfYear_vec  = convert2numvec( finviz_data(:,FINHEAD.PerformanceHalfYear) ); data = [data; PerformanceHalfYear_vec']; [FHEAD(:).PerformanceHalfYear]=i; i=i+1;
  PerformanceYear_vec  = convert2numvec( finviz_data(:,FINHEAD.PerformanceYear) ); data = [data; PerformanceYear_vec']; [FHEAD(:).PerformanceYear]=i; i=i+1;
  PerformanceYTD_vec  = convert2numvec( finviz_data(:,FINHEAD.PerformanceYTD) ); data = [data; PerformanceYTD_vec']; [FHEAD(:).PerformanceYTD]=i; i=i+1;

  Beta_vec = convert2numvec( finviz_data(:,FINHEAD.Beta) ); data = [data; Beta_vec']; [FHEAD(:).Beta]=i; i=i+1;
  AverageTrueRange_vec = convert2numvec( finviz_data(:,FINHEAD.AverageTrueRange) ); data = [data; AverageTrueRange_vec']; [FHEAD(:).AverageTrueRange]=i; i=i+1;
  VolatilityWeek_vec = convert2numvec( finviz_data(:,FINHEAD.VolatilityWeek) ); data = [data; VolatilityWeek_vec']; [FHEAD(:).VolatilityWeek]=i; i=i+1;
  VolatilityMonth_vec = convert2numvec( finviz_data(:,FINHEAD.VolatilityMonth) ); data = [data; VolatilityMonth_vec']; [FHEAD(:).VolatilityMonth]=i; i=i+1;

  SMA20day_vec  = convert2numvec( finviz_data(:,FINHEAD.SMA20day) ); data = [data; SMA20day_vec']; [FHEAD(:).SMA20day]=i; i=i+1;
  SMA50day_vec  = convert2numvec( finviz_data(:,FINHEAD.SMA50day) ); data = [data; SMA50day_vec']; [FHEAD(:).SMA50day]=i; i=i+1;
  SMA200day_vec  = convert2numvec( finviz_data(:,FINHEAD.SMA200day) ); data = [data; SMA200day_vec']; [FHEAD(:).SMA200day]=i; i=i+1;

  High50day_vec  = convert2numvec( finviz_data(:,FINHEAD.High50day) ); data = [data; High50day_vec']; [FHEAD(:).High50day]=i; i=i+1;
  Low50day_vec  = convert2numvec( finviz_data(:,FINHEAD.Low50day) ); data = [data; Low50day_vec']; [FHEAD(:).Low50day]=i; i=i+1;
  High52week_vec  = convert2numvec( finviz_data(:,FINHEAD.High52week) ); data = [data; High52week_vec']; [FHEAD(:).High52week]=i; i=i+1;
  Low52week_vec  = convert2numvec( finviz_data(:,FINHEAD.Low52week) ); data = [data; Low52week_vec']; [FHEAD(:).Low52week]=i; i=i+1;

  RSI14_vec = convert2numvec( finviz_data(:,FINHEAD.RSI14) ); data = [data; RSI14_vec']; [FHEAD(:).RSI14]=i; i=i+1;
  ChangeFromOpen_vec = convert2numvec( finviz_data(:,FINHEAD.ChangeFromOpen) ); data = [data; ChangeFromOpen_vec']; [FHEAD(:).ChangeFromOpen]=i; i=i+1;
  Gap_vec = convert2numvec( finviz_data(:,FINHEAD.Gap) ); data = [data; Gap_vec']; [FHEAD(:).Gap]=i; i=i+1;
  AnalystRecommendation_vec = convert2numvec( finviz_data(:,FINHEAD.AnalystRecommendation) ); data = [data; AnalystRecommendation_vec']; [FHEAD(:).AnalystRecommendation]=i; i=i+1;
  AverageVolume_vec = convert2numvec( finviz_data(:,FINHEAD.AverageVolume) ); data = [data; AverageVolume_vec']; [FHEAD(:).AverageVolume]=i; i=i+1;
  RelativeVolume_vec = convert2numvec( finviz_data(:,FINHEAD.RelativeVolume) ); data = [data; RelativeVolume_vec']; [FHEAD(:).RelativeVolume]=i; i=i+1;
  Price_vec = convert2numvec( finviz_data(:,FINHEAD.Price) ); data = [data; Price_vec']; [FHEAD(:).Price]=i; i=i+1;
  ChangePercentage_vec = convert2numvec( finviz_data(:,FINHEAD.ChangePercentage) ); data = [data; ChangePercentage_vec']; [FHEAD(:).ChangePercentage]=i; i=i+1;
  Volume_vec = convert2numvec( finviz_data(:,FINHEAD.Volume) ); data = [data; Volume_vec']; [FHEAD(:).Volume]=i; i=i+1;
  %TargetPrice_vec = convert2numvec( finviz_data(:,FINHEAD.TargetPrice) ); data = [data; TargetPrice_vec']; [FHEAD(:).TargetPrice]=i; i=i+1;

  FHEAD_names = fieldnames(FHEAD); % cell of field name strings

  % Save name cells
  name_cell = cell(ntickers,5);
  for i=1:ntickers
    cell_print = finviz_data(i,FINHEAD.Ticker);
    name_cell{i,1} = cell_print{1};
    cell_print = finviz_data(i,FINHEAD.Company);
    name_cell{i,2} = cell_print{1};
    cell_print = finviz_data(i,FINHEAD.Sector);
    name_cell{i,3} = cell_print{1};
    cell_print = finviz_data(i,FINHEAD.Industry);
    name_cell{i,4} = cell_print{1};
    cell_print = finviz_data(i,FINHEAD.Country);
    name_cell{i,5} = cell_print{1};
  end

  % Save all data into data2.mat
  output_matfile = strcat(finviz_dir,'/finviz.mat');
  save(output_matfile,'FHEAD','FHEAD_names','name_cell','data');

end

rmpath('functions');
end

function [numvec]=convert2numvec(cellvec)
fprintf('Converting 2 numvec...\n')
n = length(cellvec);
numvec=zeros(n,1);
for i=1:n
  cellvec_i = cellvec(i);
  numstring_i = cellvec_i{1};
  if(strcmp(numstring_i,'')==1)
    numvec(i) = NaN;
  elseif(strcmp(numstring_i(end),'%')==1)
    numstring_i = numstring_i(1:end-1);
    numvec(i) = str2num(numstring_i);
  else
    numvec(i) = str2num(numstring_i);
  end
end
end

function [FINHEAD]=set_finviz_labels()
FINHEAD = struct('Ticker',{1});
[FINHEAD(:).Company]=2;
[FINHEAD(:).Sector]=3;
[FINHEAD(:).Industry]=4;
[FINHEAD(:).Country]=5;

[FINHEAD(:).MarketCap]=6;
[FINHEAD(:).PE]=7;
[FINHEAD(:).PEForward]=8;
[FINHEAD(:).PEG]=9;

[FINHEAD(:).PS]=10;
[FINHEAD(:).PB]=11;
[FINHEAD(:).PCash]=12;
[FINHEAD(:).PFreeCash]=13;
[FINHEAD(:).DividendYield]=14;
[FINHEAD(:).PayoutRatio]=15;

[FINHEAD(:).EPS]=16;%
[FINHEAD(:).EPSGrowthThisYear]=17;%
[FINHEAD(:).EPSGrowthNextYear]=18;%
[FINHEAD(:).EPSGrowthPast5Years]=19;%
[FINHEAD(:).EPSGrowthNext5Years]=20;%

[FINHEAD(:).SalesGrowth5Years]=21; %
[FINHEAD(:).EPSGrowthLastQuarter]=22;%
[FINHEAD(:).SalesGrowthLastQuarter]=23;%
[FINHEAD(:).SharesOutstanding]=24;%

[FINHEAD(:).Float]=25;%
[FINHEAD(:).InsiderOwnership]=26; %
[FINHEAD(:).InsiderTransactions]=27; %
[FINHEAD(:).InstitutionalOwnership]=28; %
[FINHEAD(:).InstitutionalTransactions]=29; %

[FINHEAD(:).FloatShort]=30; %
[FINHEAD(:).ShortRatio]=31; %

[FINHEAD(:).ReturnOnAssets]=32;
[FINHEAD(:).ReturnOnEquity]=33;
[FINHEAD(:).ReturnOnInvestment]=34;
[FINHEAD(:).CurrentRatio]=35;
[FINHEAD(:).QuickRatio]=36;
[FINHEAD(:).LTDebtEquity]=37;
[FINHEAD(:).TotalDebtEquity]=38;

[FINHEAD(:).GrossMargin]=39; %
[FINHEAD(:).OperatingMargin]=40; %
[FINHEAD(:).ProfitMargin]=41; %

[FINHEAD(:).PerformanceWeek]=42; %
[FINHEAD(:).PerformanceMonth]=43; %
[FINHEAD(:).PerformanceQuarter]=44;%
[FINHEAD(:).PerformanceHalfYear]=45; %
[FINHEAD(:).PerformanceYear]=46; %
[FINHEAD(:).PerformanceYTD]=47; %

[FINHEAD(:).Beta]=48;
[FINHEAD(:).AverageTrueRange]=49;
[FINHEAD(:).VolatilityWeek]=50;
[FINHEAD(:).VolatilityMonth]=51;

[FINHEAD(:).SMA20day]=52; %
[FINHEAD(:).SMA50day]=53; %
[FINHEAD(:).SMA200day]=54; %

[FINHEAD(:).High50day]=55; %
[FINHEAD(:).Low50day]=56; %
[FINHEAD(:).High52week]=57; %
[FINHEAD(:).Low52week]=58;%

[FINHEAD(:).RSI14]=59; %
[FINHEAD(:).ChangeFromOpen]=60; %
[FINHEAD(:).Gap]=61; %
[FINHEAD(:).AnalystRecommendation]=62; %
[FINHEAD(:).AverageVolume]=63; %
[FINHEAD(:).RelativeVolume]=64; %
[FINHEAD(:).Price]=65; %
[FINHEAD(:).ChangePercentage]=66; %
[FINHEAD(:).Volume]=67; %
[FINHEAD(:).EarningsDate]=68;
[FINHEAD(:).TargetPrice]=69;
[FINHEAD(:).IPODate]=70;
end

author: gudbrandtandberg


text:
function fig2simulinkmaskicon(hFig,colors)
% Creates a mask for simulink subsystem from a given figure. This subsystem is
% masked with a ICON correspoding to lines on the figure. This is useful to illustrate
% the behaviour of the subsystem by a icon.
%
% You can use the simplot command to create it from the actual output. See the
% documentation of simplot (doc simplot) for further info.
%
% INPUTS:   hFig      - handle to a figure to be used (e.g. gcf)
%           colors    - (OPTIONAL) colors for different lines.
%                       Supported (default) colors: blue, red, green, magenta, yellow, cyan, black
%
% OUTPUTS:  Creates a simulink subsystem with a given mask
%           Prints out the Mask/Icon string
%
% Author: Paul Proteus (e-mail: proteus.paul (at) yahoo (dot) com)
% Version: 1.0
% Changes tracker:  13.07.2010  - First version
%
% Example:  h = figure; hold on;
%           plot([0:0.1:10],sin([0:0.1:10]));
%           plot([0:0.5:12],cos([0:0.5:12]));
%           fig2simulinkmaskicon(gcf);
%           fig2simulinkmaskicon(gcf,'yellow');
%           fig2simulinkmaskicon(gcf,{'red','magenta'});

%% Check required inputs
if ~strcmp(get(hFig,'Type'),'figure'),
    % Checking the validity of the figure handle
    error('hFig argument is not a valid figure handle');
end

% Find the lines
hLine = findall(hFig,'type','line');

% Check if there is at least one line
if isempty(hLine),
    % There is no line to be used. Figure has to have at least one.
    error('There are no lines in this figure. Cannot create subsystem mask.');
end

% If the colors are not provided, use the default
if ~exist('colors','var'),
    % define the default colors
    colors = { ...
        'blue',    ...
        'red',     ...
        'green',   ...
        'magenta', ...
        'yellow',  ...
        'cyan',    ...
        'black',   ...
        };
end

% Check if the colors where not provided as a string
if ischar(colors),
    % Make it cell-array for further use
    colors = {colors};
end

%% Initialization
XData = get(hLine,'XData');
YData = get(hLine,'YData');

if ~iscell(XData),
    % In case there is only one line specified in this hLine
    XData  = {XData};
    YData  = {YData};
end

nGraphs = length(XData);

%% Prepare and print out the MASK/ICON string
disp('Copy following lines to the Mask/Icon of the selected subsystem:');
disp('------- START OF COPY SECTION');

allPreparedString = [];
% Loop through all lines in the figure
for idx=1:nGraphs,
    % Color selection
    colorIdx = min(idx,length(colors)); % if too many lines - used the last color
    preparedString = ['color(''' colors{colorIdx} ''');'];
    disp(preparedString);

    allPreparedString = [allPreparedString preparedString];

    % Plot the line
    preparedString = 'plot(';
    preparedString = [preparedString ...
                      '[' num2str(XData{idx}) '],' ...
                      '[' num2str(YData{idx}) '],'];
    preparedString = [preparedString(1:end-1) ');'];
    % get rid of the too many spaces
    preparedString = regexprep(preparedString, '\s*', ' ');

    disp(preparedString);
    allPreparedString = [allPreparedString preparedString];
end

disp('------- END');

%% Create a new simulink subsystem block and set it up
systemName = ['model_fig2simulinkmaskicon_' num2str(round(rand*100000))];
new_system(systemName);
open_system(systemName);

add_block('built-in/SubSystem',[systemName '/SS_Masked']);

% Enlarge it and mask it
hBlock = get_param([systemName '/SS_Masked'],'handle');
set(hBlock,'Position',[45 35 181 113]);
set(hBlock,'Mask','on');
set(hBlock,'MaskDisplay',allPreparedString);

author: Sable

repo: mcbench-benchmarks


text:
%% Creating test problems and initial guesses
% We demonstrate how to use Tensor Toolbox |create_problem| and
% |create_guess| functions to create test problems for fitting algorithms.

%% Creating a CP test problem
% The |create_problem| function allows a user to generate a test problem
% with a known solution having a pre-specified solution. The
% |create_problem| function generates both the solution (as a |ktensor| for
% CP) and the test data (as a |tensor|). We later show that a
% pre-specificed solution can be used as well.

% Create a problem
info = create_problem('Size', [5 4 3], 'Num_Factors', 3, 'Noise', 0.10);

%%

% Display the solution created by create_problem
soln = info.Soln

%%

% Display the data created by create_problem
data = info.Data

%%

% The difference between true solution and measured data should match the
% specified 10% noise.
diff = norm(full(info.Soln) - info.Data)/norm(full(info.Soln))

%% Creating a Tucker test problem
% The |create_problem| function can also be used to create Tucker problems
% by specifying the |'Type'| as |'Tucker'|. In this case, the
% |create_problem| function generates both the solution (as a |ttensor| for
% Tucker) and the test data (as a |tensor|).

% Create a problem
info = create_problem('Type', 'Tucker', 'Size', [5 4 3], 'Num_Factors', [3 3 2]);

%%

% Display the Tucker-type solution created by create_problem
soln = info.Soln

%%

% Difference between true solution and measured data (default noise is 10%)
diff = norm(full(info.Soln) - info.Data)/norm(full(info.Soln))

%% Recreating the same test problem
% We can recreate exactly the same test problem when we use the same random
% seed and other parameters.

% Set-up, including specifying random state
sz = [5 4 3]; %<- Size
nf = 2; %<- Number of components
state = RandStream.getDefaultStream.State; %<- Random state

%%

% Generate first test problem
info1 = create_problem('Size', sz, 'Num_Factors', nf, 'State', state);

%%

% Generate second identical test problem
info2 = create_problem('Size', sz, 'Num_Factors', nf, 'State', state);

%%

% Check that the solutions are identical
tf = isequal(info1.Soln, info2.Soln)

%%

% Check that the data are identical
diff = norm(info1.Data - info2.Data)

%% Checking default parameters and recreating the same test problem
% The |create_problem| function returns the parameters that were used to
% generate it. These can be used to see the defaults. Additionally, if
% these are saved, they can be used to recreate the same test problems for
% future experiments.

% Generate test problem and use second output argument for parameters.
[info1,params] = create_problem('Size', [5 4 3], 'Num_Factors', 2);

%%

% Here are the parameters
params

%%

% Recreate an identical test problem
info2 = create_problem(params);

%%

% Check that the solutions are identical
tf = isequal(info1.Soln, info2.Soln)

%%

% Check that the data are identical
diff = norm(info1.Data - info2.Data)

%% Options for creating factor matrices, core tensors, and lambdas
% Any function with two arguments specifying the size can be used to
% generate the factor matrices. This is specified by the
% |'Factor_Generator'| option to |create_problem|.
%
% Pre-defined options for |'Factor_Generator'| for creating factor matrices
% (for CP or Tucker) include:
%
% * |'rand'| - Uniform on [0,1]
% * |'randn'| - Gaussian with mean 0 and std 1
% * |'orthogonal'| - Generates a random orthogonal matrix. This option only
% works when the number of factors is less than or equal to the smallest
% dimension.
% * |'stochastic'| - Generates nonnegative factor matrices so that each
% column sums to one.
%
% Pre-defined options for |'Lambda_Generator'| for creating lambda vector
% (for CP) include:
%
% * |'rand'| - Uniform on [0,1]
% * |'randn'| - Gaussian with mean 0 and std 1
% * |'orthogonal'| - Creates a random vector with norm one.
% * |'stochastic'| - Creates a random nonnegative vector whose entries sum
% to one.
%
% Pre-defined options for |'Core_Generator'| for creating core tensors (for
% Tucker) include:
%
% * |'rand'| - Uniform on [0,1]
% * |'randn'| - Gaussian with mean 0 and std 1

% Here is ane example of a custom factor generator
factor_generator = @(m,n) 100*rand(m,n);
info = create_problem('Size', [5 4 3], 'Num_Factors', 2, ...
    'Factor_Generator', factor_generator, 'Lambda_Generator', @ones);
first_factor_matrix = info.Soln.U{1}

%%

% Here is an example of a custom core generator for Tucker:
info = create_problem('Type', 'Tucker', 'Size', [5 4 3], ...
    'Num_Factors', [2 2 2], 'Core_Generator', @tenones);
core = info.Soln.core

%%

% Here's another example for CP, this time using a function to create
% factor matrices such that the inner products of the columns are
% prespecified.
info = create_problem('Size', [5 4 3], 'Num_Factors', 3, ...
    'Factor_Generator', @(m,n) tt_ccong(m,n,.9));
U = info.Soln.U{1};
congruences = U'*U

%% Generating data from an existing solution
% It's possible to skip the solution generation altogether and instead just
% generate appropriate test data.

% Manually generate a test problem (or it comes from some
% previous call to |create_problem|.
soln = ktensor({rand(50,3), rand(40,3), rand(30,3)});

% Use that soln to create new test problem.
info = create_problem('Soln', soln);

% Check whether solutions is equivalent to the input
iseq = isequal(soln,info.Soln)

%% Creating dense missing data problems
% It's possible to create problems that have a percentage of missing data.
% The problem generator randomly creates the pattern of missing data.

% Specify 25% missing data as follows:
[info,params] = create_problem('Size', [5 4 3], 'M', 0.25);

%%

% Here is the pattern of known data (1 = known, 0 = unknown)
info.Pattern

%%

% Here is the data (incl. noise) with missing entries zeroed out
info.Data

%% Creating sparse missing data problems.
% If |Sparse_M| is set to true, then the data returned
% is sparse. Moreover, the dense versions are never explicitly created.
% This option only works when M >= 0.8.

% Specify 80% missing data and sparse
info = create_problem('Size', [5 4 3], 'M', 0.80, 'Sparse_M', true);

%%

% Here is the pattern of known data
info.Pattern

%%

% Here is the data (incl. noise) with missing entries zeroed out
info.Data

%% Create missing data problems with a pre-specified pattern
% It's also possible to provide a specific pattern (dense or sparse) to be
% used to specify where data should be missing.

% Create pattern
P = tenrand([5 4 3]) > 0.5;
% Create test problem with that pattern
info = create_problem('Size', size(P), 'M', P);
% Show the data
info.Data

%% Creating sparse problems (CP only)
% If we assume each model parameter is the input to a Poisson process, then
% we can generate a sparse test problems. This requires that all the factor
% matrices and lambda be nonnegative. The default factor generator
% ('randn') won't work since it produces both positive and negative values.

% Generate factor matrices with a few large entries in each column; this
% will be the basis of our soln.
sz = [20 15 10];
nf = 4;
A = cell(3,1);
for n = 1:length(sz)
    A{n} = rand(sz(n), nf);
    for r = 1:nf
        p = randperm(sz(n));
        idx = p(1:round(.2*sz(n)));
        A{n}(idx,r) = 10 * A{n}(idx,r);
    end
end
S = ktensor(A);
S = normalize(S,'sort',1);
%%

% Create sparse test problem based on provided solution. The
% 'Sparse_Generation' says how many insertions to make based on the
% provided solution S. The lambda vector of the solution is automatically
% rescaled to match the number of insertions.
info = create_problem('Soln', S, 'Sparse_Generation', 500);
num_nonzeros = nnz(info.Data)
total_insertions = sum(info.Data.vals)
orig_lambda_vs_rescaled = S.lambda ./ info.Soln.lambda

%% Generating an initial guess
% The |create_guess| function creates a random initial guess as a cell
% array of matrices. Its behavior is very similar to |create_problem|. A
% nice option is that you can generate an initial guess that is a
% pertubation of the solution.

info = create_problem;

% Create an initial guess to go with the problem that is just a 5%
% pertubation of the correct solution.
U = create_guess('Soln', info.Soln, 'Factor_Generator', 'pertubation', ...
    'Pertubation', 0.05);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function exact = p06_exact ( )

%*****************************************************************************80
%
%% P06_EXACT returns the exact integral for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 1.460447131787105;

  return
end


text:
% KALMAN_FILTER - Under construction..

function [X, CovX] = kalman_filter_old(measurement, ...
                                   transition_matrix, ...
                                   process_noise, ...
                                   measurement_model_matrix, ...
                                   measurement_noise, ...
                                   x, ...
                                   C)

warning('This function is deprecated')

error('Not ready yet..')
return

N = 0;
D = 1000;

X = zeros(D,N);
if nargout >= 2
  CovX = zeros(D,D,10);
end

for n=1:N
  % Get the parameters
  A = transition_matrix(n);
  Q = process_noise(n);
  y = measurement(n);
  H = measurement_model_matrix(n);
  R = measurement_noise(n);

  % Prediction step
  x = A*x;
  C = A*C*A' + Q;

  % Update step
  if length(y) < D
    v = y - H*x;
    S = H*C*H' + R;
    K = C*H'/S;
    x = m + K*v;
    C = C - K*S*K';
  else
    K = C*H';
    L = chol(C +
  end

  % Store the distribution
  X(:,n) = x;
  if nargout >= 2
    CovX(:,:,n) = C;
  end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% DEMROBOTWIRELESS1 Wireless Robot data from University of Washington, without dynamics and without back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('ftc');
latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [signal, state] = flt_dynamicloreta(varargin)
% Return the current source density for a given head model and data using
% the cortically-constrained LORETA (low resolution electrical
% tomographic analysis) with a Bayesian update scheme for hyperparameters.
% The reconstructed CSD time-series (or source potential maps) will be
% stored in signal.srcpot. This matrix has dimension [num_voxels x num_samples].
%
% Author: Tim Mullen, Jan 2013, SCCN/INC/UCSD
%         Alejandro Ojeda, Jan 2013, SCCN/INC/UCSD
%         Christian Kothe, Jan 2013, SCCN/INC/UCSD

if ~exp_beginfun('filter'), return; end

declare_properties('name','Dynamic LORETA', 'experimental',true, 'independent_channels',false, 'independent_trials',false);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg_nogui({'K','ForwardModel'},[],[],'Forward model (matrix)','shape','matrix'), ...
    arg_nogui({'L','LaplacianOperator'},[],[],'Laplacian operator. Sparse matrix of N sources x N sources, this is matrix is used as the square root of the precision matrix of the sources.'), ...
    arg_sub({'options','LoretaOptions'},{},...
        { ...
        arg({'maxTol','MaxTolerance'},1e-12,[0 Inf],'Tolerance for hyperparameter update loop','cat','Loreta Options'), ...
        arg({'maxIter','MaxIterations'},100,[1 Inf],'Maximum iterations for hyperparameter update loop','cat','Loreta Options'), ...
        arg({'gridSize','GridSize'},100,[1 Inf],'Lambda grid size.'), ...
        arg({'history','TrackHistory'},false,[],'Track history for hyperparameters'), ...
        arg({'verbose','VerboseOutput'},false,[],'Verbosity','cat','Loreta Options'), ...
        arg({'initNoiseFactor','InitialNoiseFactor'},0.001,[0 Inf],'Fraction of noise level. Used for initializing alpha parameter','cat','Loreta Options') ...
        arg({'block_size','BlockSize'},5, [], 'Block granularity for processing. The inverse operator will be updated using blocks of this many samples. This assumes that the inverse solution is spatially stationary over this many samples.'), ...
        arg({'skipFactor','SkipFactor'},0,[0 Inf],'Number of blocks to skip'), ...
        arg({'maxblocks','MaxBlocks'},Inf,[0 Inf],'Maximum number of blocks'), ...
        arg({'standardize','Standardize'},'all',{'none','channels','all'},'Rescale data to unit variance. If ''channels'', standardization is carried out across channels for each time point. If ''all'' each data sample is normalized by the standard deviation taken over all data.'), ...
        arg({'useGPU','UseGPU'},false,[],'Use GPU to accelerate computation.'), ...
        },'Additional options for Loreta function'), ...
    arg({'verb','Verbosity'},false,[],'Verbose output'), ...
    arg_nogui({'state','State'},[],[],'State object. When provided, hyperparameters will be estimated adaptively from prior state'));
if verb
    fprintf('Estimating current source density using cLORETA (%s)\n',mfilename);
end

[nchs, npnts, ntrs] = size(signal.data);
if isempty(block_size) || block_size > npnts
    block_size = npnts;
end
numsplits    = floor(npnts/block_size);

% if necessary, cast to double-precision
if ~strcmpi(class(signal.data),'double')
    signal.data = double(signal.data);
end

% normData the data
if ~strcmpi(normData,'none')
    switch normData
        case 'channels'
            scale = std(signal.data,[],1);
        case 'time'
            scale = std(signal.data,[],2);
        case 'all'
            scale = std(signal.data(:));
    end
    signal.data = bsxfun(@rdivide,signal.data,scale);
%     scale = std(signal.data(:));
%     signal.data = signal.data./scale;
end

if isempty(state) || ~isfield(state,'iLV') || isempty(state.iLV)
    if verb
        fprintf('...computing SVD of LFM.\n');
    end
    % mode is offline or we are initializing online filter
    % perform one-time SVD for faster computation.
    [U,S,V]      = svd(K/L,'econ');
    state.iLV    = L\V;
    state.s2     = diag(S).^2; %s^2
    state.Ut     = U';
    state.sigma2 = repmat({options.sigma2},1,ntrs);
    state.tau2   = repmat({options.tau2},1,ntrs);
end

if npnts == 0
    % no data
    signal.srcpot    = [];
    state.srcweights = [];
    exp_endfun; return;
end

signal.srcpot    = zeros([size(K,2), npnts, ntrs]);
state.srcweights = zeros(size(L,1),nchs);
sum_srcweights   = zeros(size(L,1),nchs);
signal.loretaHistory = struct([]);

if verb
    fprintf('...assuming %d stationary blocks of length %d\n',numsplits,block_size);
end

% loop over all trials
for tr=1:ntrs
    if verb
        fprintf('\nTrial (%d\%d).',tr,ntrs);
    end
    k = 0;
    % loop over sub-blocks and estimate CSD for each block
    for i=0:skipFactor+1:numsplits-1
        if verb
            if i+1 >= floor(numsplits*(k+1)/10)
                k = k + 1;
                fprintf('%0.3g%%...',round((i/numsplits)*100));
            end
        end
        range = 1+floor(i*npnts/numsplits) : min(npnts,floor((i+1)*npnts/numsplits));
        % call (dynamic bayesian) loreta estimator
        [signal.srcpot(:,range,tr), state.sigma2{tr}, state.tau2{tr}, state.srcweights, tmpHist] ...
            = dynamicLoreta( signal.data(:,range,tr), state.Ut, state.s2, state.iLV,...
                             state.sigma2{tr}, state.tau2{tr}, options);
        if ~isempty(tmpHist)
            signal.loretaHistory{tr} = [signal.loretaHistory{tr},tmpHist];
        end

        if skipFactor > 0
            % estimate CSD for samples between blocks using current inverse operator
            range = 1+floor((i+1)*npnts/numsplits) : min(npnts,floor((i+skipFactor+1)*npnts/numsplits));
            signal.srcpot(:,range,tr) = state.srcweights*signal.data(:,range,tr);
        end

        % running sum
        sum_srcweights = sum_srcweights + state.srcweights;
    end
end

if numsplits > 1
    % store the mean inverse operator over all splits
    state.srcweights = sum_srcweights/(numsplits*ntrs);
end

if ~strcmpi(normData,'none')
    % recale data to original units
%     signal.srcpot     = signal.srcpot*scale;
%     state.srcweights  = state.srcweights/scale;
    signal.srcpot = bsxfun(@times,signal.srcpot,scale);
%     signal.srcpot = bsxfun(@rdivide,signal.srcpot,std(signal.srcpot,[],1));
%     state.srcweights  = bsxfun(@times,state.srcweights,scale'); %state.srcweights/mean(scale);
end

if verb
    fprintf('done.\n');
end

exp_endfun;

repo: REST


text:
function [varargout] = calc_stiff_matrix_val(varargin)

% CALC_STIFF_MATRIX is implemented as mex file
%
% $Id$

error('The mex file %s is missing', [mfilename '.' mexext]);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ThermalComputer < handle

    properties (Access = public)
        computation
    end

    properties (Access = private)
        testName
    end

    methods (Access = public)
        function obj = ThermalComputer(cParams)
            obj.testName = cParams.testName;
        end

        function compute(obj)
            s = obj.createFEMparameters();
            obj.computation = FEM.create(s);
            obj.computation.solve();
        end
    end

    methods (Access = private)

        function s = createFEMparameters(obj)
            gidParams = obj.createGiDparameters();
            s.dim       = gidParams.pdim;
            s.type      = gidParams.ptype;
            s.scale     = gidParams.scale;
            s.mesh      = gidParams.mesh;
            s.dirichlet = gidParams.dirichlet;
            s.pointload = gidParams.pointload;
        end

        function gidParams = createGiDparameters(obj)
            file = obj.testName;
            gidReader = FemInputReader_GiD();
            gidParams = gidReader.read(file);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function f = round(f)
%ROUND   Pointwise round function of a CHEBTECH.
%   G = FLOOR(F) returns the CHEBTECH G such that G(X) = FLOOR(F(x)) for each x
%   in [-1, 1].
%
%   If F is complex, then the G = ROUND(REAL(F)) + 1i*ROUND(IMAG(F)).
%
%   Note that ROUND() assumes the output G(X) is a constant. If it is not, then
%   garbage is returned with no warning.
%
% See also CEIL, FLOOR, FIX.

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

% Evaluate at the two end points, and an arbitrary interior point:
arbitraryPoint = 0.1273881594;
fx = feval(f, [-1 ; arbitraryPoint ; 1]);
% Take the mean:
meanfx = mean(fx, 1);
% Compute the round:
f.coeffs = round(meanfx);

end

repo: chebfun


text:
function Y = sne(X, d, perplexity)
%SNE Implementation of Stochastic Neighbor Embedding
%
%   mappedX = sne(X, no_dims, perplexity)
%
% Runs the Stochastic Neighbor Embedding algorithm using a Gaussian kernel
% with fixed perplexity. The high-dimensional datapoints are specified by X.
% The target dimensionality is specified in no_dims, and the variance of
% the Gaussian kernel may be specified through perplexity (default = 30).
% The function returns the embedded points in Y.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    if ~exist('d', 'var') || isempty(d)
        d = 2;
    end
    if ~exist('perplexity', 'var') || isempty(perplexity)
        perplexity = 30;
    end

    % Initialize some variables
    n = size(X, 1);                 % number of instances
    eta = .05;                      % learning rate
    max_iter = 2000;                % maximum number of iterations
    jitter = 0.3;                   % initial jitter
    jitter_decay = 0.99;            % jitter decay
    momentum = 0.5;                 % initial momentum
    final_momentum = 0.8;           % final momentum
    mom_switch_iter = 750;          % iteration where momentum changes

    % Initialize embedding coordinates randomly (close to origin)
    Y = 0.0001 * rand(n, d);
    dC = zeros(n, d);
    y_incs = zeros(n, d);

    % Compute Gaussian kernel for high-dimensional data representation
    P = x2p(X, perplexity, 1e-5);                                        % use fixed perplexity
    P = max(P, eps);

    % Iterating loop
    for iter=1:max_iter

        % Compute Gaussian kernel for low-dimensional data representation
        sum_Y = sum(Y .^ 2, 2);                                                     % precomputation for pairwise distances
        Q = exp(-bsxfun(@plus, sum_Y, bsxfun(@plus, sum_Y', -2 * Y * Y')) ./ 2);    % Gaussian probabilities
        Q(1:n+1:end) = 0;
        Q = Q ./ repmat(sum(Q, 2), [1 n]);
        Q = max(Q, eps);

        % Compute cost function between P and Q
        if ~rem(iter, 20)
            costs = sum(P .* log((P + eps) ./ (Q + eps)), 2) ./ n;              % division by n corrects for # of datapoints
            cost = sum(costs);
            disp(['Iteration ' num2str(iter) ': error is ' num2str(cost)]);
        end

        % Compute gradient
        PQ = P - Q + P' - Q';
        for i=1:n
            dC(i,:) = sum(bsxfun(@times, bsxfun(@minus, Y(i,:), Y), PQ(i,:)'), 1);
        end

        % Perform the gradient search
        y_incs = momentum * y_incs - eta * dC;
        Y = Y + y_incs;
		Y = Y + jitter * randn(size(Y));
        Y = bsxfun(@minus, Y, mean(Y, 1));

        % Reduce jitter over time and change momentum
        jitter = jitter * jitter_decay;
        if iter == mom_switch_iter
            momentum = final_momentum;
        end
    end

repo: cluster


text:
function files = scan_get_files(n, filt, mesg, wd)
% Drop-in replacement for spm_get.
% Main advantage is that it can use spm_select for SPM5+ to mimic older spm_get behavior.
%
% :Usage:
% ::
%
%     files = scan_get_files(n, filt, mesg, wd)

    if(~exist('wd', 'var') || isempty(wd))
        wd = pwd();
    end

    switch(spm('Ver', [], 1))
        case 'SPM2'
            files = spm_get(n, filt, mesg, wd);
        case {'SPM5', 'SPM8'}
            if(any(n < 0))
                typ = 'dir';
                n = abs(n);
                select_filt = strrep(filt, '*', '.*');
            else
                switch(deblank(filt))
                    case {'*.img' '*img' '*.nii' '*nii' '*.hdr' '*hdr'}
                        typ = 'image';
                        select_filt = '.*';
                    case {'*.mat' '*mat'}
                        typ = 'mat';
                        select_filt = '.*';
                    otherwise
                        typ = 'any';
                        select_filt = strrep(filt, '*', '.*');
                end
            end
            files = spm_select(n, typ, mesg, [], wd, select_filt);
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef synthetic

methods(Static)

function signal = picket_fence(N)
    import spx.discrete.number;
    if ~number.is_perfect_square(N)
        error('N must be perfect square');
    end
    n = sqrt(N);
    signal = zeros(1, N);
    signal(1:n:end) = 1;
end

function X = uniform(N, S)
    % Generates uniformly distributed signals

    % Generate Gaussian vectors
    X = randn(N, S);
    % Make them unit length
    X = spx.norm.normalize_l2(X);
end

end

end

author: indigits

repo: sparse-plex


text:
function varargout = min(varargin)
%MIN   Minimum value of a SPHEREFUN in one direction.
%   MIN(f) returns a chebfun representing the minimum of the SPHEREFUN along the
%   latitude direction, i.e, MIN(f) = @(lambda) max( f ( lambda, : ) )
%
%   MIN(f, [], dim) returns a CHEBFUN representing the minimum of f along the
%   DIM direction. If DIM = 1 is along the latitude-direction and DIM = 2 is along the
%   longitude-direction.
%
%   WARNING: This function is not always accurate to the expected precision.
%
%   For the global minimum use MIN2.
%
% See also SPHEREFUN/MAX, SPHEREFUN/MAX2, SPHEREFUN/MIN2, SPHEREFUN/MINANDMAX2.

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

[varargout{1:nargout}] = min@separableApprox(varargin{:});

end

author: chebfun


text:
function L=Cell_watershed(I0, t4)
hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I0), hy, 'replicate');
Ix = imfilter(double(I0), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
%figure, imshow(gradmag,[]), title('gradmag')

% first parameter
se = strel('disk', t4);
Io = imopen(I0, se);
%figure, imshow(Io), title('Io')

Ie = imerode(I0, se);
Iobr = imreconstruct(Ie, I0);
% figure, imshow(Iobr), title('Iobr')

Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
% figure, imshow(Iobrcbr), title('Iobrcbr')

fgm = imregionalmax(Iobrcbr,8);
%figure, imshow(fgm), title('fgm')

I2 = I0;
I2(fgm) = 255;
% figure, imshow(I2), title('fgm superimposed on original image')

% Second parameter
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

% Third parameter
fgm4 = bwareaopen(fgm3, 5);
I3 = I0;
I3(fgm4) = 255;
% figure, imshow(I3)
% title('fgm4 superimposed on original image')

bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
%figure, imshow(bw), title('bw')
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
% figure, imshow(bgm), title('bgm')

gradmag2 = imimposemin(gradmag, bgm | fgm4);

L = watershed(gradmag2);
I4 = I0;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;
% figure, imshow(I4)
% title('Markers and object boundaries superimposed on original image')

Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
% figure, imshow(Lrgb)
% title('Lrgb')
figure, imshow(I0), hold on
himage = imshow(Lrgb);
set(himage, 'AlphaData', 0.3);
title('nurse cell blobs')

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function r8vec_index_delete_one_test ( )

%*****************************************************************************80
%
%% R8VEC_INDEX_DELETE_ONE_TEST tests R8VEC_INDEX_DELETE_ONE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 25;
  n = 0;
  x = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_INDEX_DELETE_ONE_TEST\n' );
  fprintf ( 1, '  R8VEC_INDEX_DELETE_ONE deletes one copy of a\n' );
  fprintf ( 1, '  particular value.\n' );
  fprintf ( 1, '\n' );

  xval = 8.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  xval = 7.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  seed = 123456789;

  for i = 1 : 20
    [ xval, seed ] = r8_uniform_ab ( 0.0, 20.0, seed );
    xval = round ( xval );
    fprintf ( 1, '  %f\n', xval );
    [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );
  end

  xval = 7.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  xval = 8.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %3d  %3d  %9f  %9f\n', i, indx(i), x(i), x(indx(i)) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call R8VEC_INDEX_DELETE_ONE to delete one value of 8:\n' );

  xval = 8.0;
  [ n, x, indx ] = r8vec_index_delete_one ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %3d  %3d  %9f  %9f\n', i, indx(i), x(i), x(indx(i)) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ma_updateDispAxis(handles, stateLog, orbitNumToPlot)
%ma_updateDispAxis Summary of this function goes here
%   Detailed explanation goes here
    hDispAxes = handles.dispAxes;
    hShowSoICheckBox = handles.showSoICheckBox;
    hShowChildrenCheckBox = handles.showChildrenCheckBox;
    hDispAxisTitleLabel = handles.dispAxisTitleLabel;
    hShowOtherSpacecraftCheckBox = handles.showOtherSpacecraftCheckBox;
    hShowChildBodyMarker = handles.showChildBodyMarker;

    if(isempty(orbitNumToPlot))
        orbitNumToPlot = 1;
        set(hDispAxes,'UserData',orbitNumToPlot);
    end

    showSoICheckBoxINT = get(hShowSoICheckBox,'value');
    if(showSoICheckBoxINT == 1)
        showSoI = true;
    else
        showSoI = false;
    end

    hShowChildrenCheckBoxINT = get(hShowChildrenCheckBox,'value');
    if(hShowChildrenCheckBoxINT == 1)
        showChildBodies = true;
    else
        showChildBodies = false;
    end

    hShowOtherSpacecraftCheckBoxINT = get(hShowOtherSpacecraftCheckBox,'value');
    if(hShowOtherSpacecraftCheckBoxINT == 1)
        showOtherSC = true;
    else
        showOtherSC = false;
    end

    showChildBodyMarkerINT = get(hShowChildBodyMarker,'value');
    if(showChildBodyMarkerINT == 1)
        showChildMarker = true;
    else
        showChildMarker = false;
    end

    axes(hDispAxes);
    cla(gca);
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    celBodyData = getappdata(handles.ma_MainGUI,'celBodyData');
    plotStateLog(stateLog, handles, showSoI, showChildBodies, showChildMarker, showOtherSC, orbitNumToPlot, hDispAxisTitleLabel, maData, celBodyData);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = c8_cube_root ( x )

%*****************************************************************************80
%
%% C8_CUBE_ROOT returns the principal cube root of a C8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex X, the number whose cube root is desired.
%
%    Output, complex VALUE, the cube root of X.
%
  a = real ( x );
  b = imag ( x );
  mag = sqrt ( a * a + b * b );

  if ( mag == 0.0 )

    value = 0.0;

  else

    theta = atan2 ( b, a );

    value = mag.^( 1 / 3 ) * ( cos ( theta / 3 ) + i * sin ( theta / 3 ) );

  end

  return
end


text:
function K=compute_permutation_constraint4(V)

% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, September 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

%[B11,B12,...,B33]=lambda1*v1+lambda2*v2+lambda3*v3

N=size(V,2); %dimension of the kernel
n=4; %dimension of Bij
idx=[1 2 3 4; 2 5 6 7; 3 6 8 9; 4 7 9 10];

%1.-Generation of the first set of equations Bii.Bjj=Bij.Bii  (n(n-1)/2 eqs).
nrowsK=n*(n-1)/2+n*(n-1)*n/2;
ncolsK=N*(N+1)/2;
K=zeros(nrowsK,ncolsK);

t=1;
for i=1:n
    for j=i+1:n
        offset=1;
        for a=1:N
            for b=a:N
                if a==b
                    K(t,offset)=V(idx(i,i),a)*V(idx(j,j),a)-V(idx(i,j),a)*V(idx(i,j),a);
                else
                    K(t,offset)=V(idx(i,i),a)*V(idx(j,j),b)-V(idx(i,j),a)*V(idx(i,j),b)+...
                                V(idx(i,i),b)*V(idx(j,j),a)-V(idx(i,j),b)*V(idx(i,j),a);
                end
                offset=offset+1;
            end

        end
        t=t+1;
        %fprintf('t:%d\t offset:%d\n',t,offset);
    end
end

for k=1:n
    for j=k:n
        for i=1:n
            if (i~=j & i~=k)
                offset=1;
                for a=1:N
                    for b=a:N
                        if a==b
                            K(t,offset)=V(idx(i,j),a)*V(idx(i,k),a)-V(idx(i,i),a)*V(idx(j,k),a);
                        else
                            K(t,offset)=V(idx(i,j),a)*V(idx(i,k),b)-V(idx(i,i),a)*V(idx(j,k),b)+...
                                        V(idx(i,j),b)*V(idx(i,k),a)-V(idx(i,i),b)*V(idx(j,k),a);
                        end
                        offset=offset+1;
                    end

                end
                t=t+1;
                %fprintf('t:%d\t offset:%d\n',t,offset);
            end
        end
    end
end

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
function KLDemo
% Illustration of information gains with Bayesian fusion
% FORMAT KLDemo)
%
%--------------------------------------------------------------------------
% This routine  illustrates the benefit of multimodal or Bayesian fusion in
% terms of conditional dependencies among parameters. In other words, it
% shows that even if one data modality contains no information about a
% particular set of parameters, it can help resolve uncertainty about
% another set and thereby disclose information contained in the other
% modality. This is illustrated here using a simple linear model with
% neuronal and haemodynamic parameters to show that EEG can provide some
% information gain, in relation to haemodynamic parameters.
%
% comment the orthogonalisation of the fMRI design matrix below to see the
% effect of conditional dependencies on the haemodynamic information gain
% afforded by EEG data
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Setup and preliminaries
%--------------------------------------------------------------------------
m    = 64;                         % number of observations
n    = 8;                          % number of parameters
j{1} = 1:n;                        % indices of both parameters
j{2} = 1:(n/2);                    % indices of neuronal parameters
j{3} = (1 + n/2):n;                % indices of haemodynamic parameters

% design matrices
%--------------------------------------------------------------------------
XEEG      = randn(m,n);
XMRI      = randn(m,n);

% make EEG design uninformative about haemodynamic parameters
%--------------------------------------------------------------------------
XEEG(:,j{3}) = 0;

% and orthogonalised fMRI design with respect to the EEG design
%--------------------------------------------------------------------------
% XMRI      = spm_orth(XMRI);

B         = randn(n,1);            % parameters
YEEG      = XEEG*B + randn(m,1)/8; % EEG data
YMRI      = XMRI*B + randn(m,1)/4; % MRI data

% model inversion using parametric empirical Bayes
%==========================================================================
PEEG{1}.X = XEEG;
PEEG{1}.C = {eye(m,m)};
PEEG{2}.X = zeros(n,1);
PEEG{2}.C = eye(n,n);

CEEG = spm_PEB(YEEG,PEEG,1);       % inversion using EEG data

PMRI{1}.X = XMRI;
PMRI{1}.C = {eye(m,m)};
PMRI{2}.X = zeros(n,1);
PMRI{2}.C = eye(n,n);

CMRI = spm_PEB(YMRI,PMRI,1);       % inversion using MRI data

PMRI{1}.X = XMRI;
PMRI{1}.C = {eye(m,m)};
PMRI{2}.X = CEEG{2}.E;             % Bayesian belief updating
PMRI{2}.C = CEEG{2}.C;             % using posteriors from EEG inversion

CMRE = spm_PEB(YMRI,PMRI,1);       % inversion using EEG and MRI data

% evaluate fMRI posteriors using Bayesian model reduction
%==========================================================================
[F,sE,sC] = spm_log_evidence(CMRE{2}.E,CMRE{2}.C,CEEG{2}.E,CEEG{2}.C,PEEG{2}.X,PEEG{2}.C);
CMRR.E    = sE;
CMRR.C    = sC;

%  evaluate information gain is entailed divergence
%==========================================================================
%  crucially, we will use the Bayesian model reduction estimate which means
%  we never have to actually invert the fMRI data
%--------------------------------------------------------------------------
for i = 1:length(j)          %  loop over different subsets of parameters

    % subsets
    %----------------------------------------------------------------------
    k      = j{i};

    % information gains
    %----------------------------------------------------------------------
    D(1,i) = spm_kl_normal(CEEG{2}.E(k),CEEG{2}.C(k,k),PEEG{2}.X(k),PEEG{2}.C(k,k));
    D(2,i) = spm_kl_normal(   CMRR.E(k),   CMRR.C(k,k),PEEG{2}.X(k),PEEG{2}.C(k,k));
    D(3,i) = spm_kl_normal(CMRE{2}.E(k),CMRE{2}.C(k,k),CEEG{2}.E(k),CEEG{2}.C(k,k));
    D(4,i) = spm_kl_normal(CMRE{2}.E(k),CMRE{2}.C(k,k),CMRI{2}.E(k),CMRI{2}.C(k,k));

    %  cow divergence between reduced and direct MRI posterior
    %----------------------------------------------------------------------
    D(5,i) = spm_kl_normal(CMRR.E(k),CMRR.C(k,k),CMRI{2}.E(k),CMRI{2}.C(k,k));

end

% illustrate results
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1'); clf;
str{1} =  '0 -> EEG';
str{2} =  '0 -> MRI';
str{3} =  'EEG -> MRI & EEG';
str{4} =  'MRI -> MRI & EEG';

%  show for information gains
%--------------------------------------------------------------------------
for i = 1:4
    subplot(3,2,i), bar(D(i,:))
    title(str{i},'fontsize',16)
    set(gca,'XTickLabel',{'both','neuronal','haemodynamic'})
    axis square
end
a = axis;

%  show  divergence between direct and reduced posterior
%--------------------------------------------------------------------------
subplot(3,1,3),bar(D(5,:))
title('KL between estimated and direct fMRI','fontsize',12)
set(gca,'XTickLabel',{'both','neuronal','haemodynamic'})
axis (a), axis square

author: spm


text:
function A = slgda(K, nums, sol)
%SLGDA Performs Baudat's Generalized Discriminant Analysis
%
% $ Syntax $
%   - A = slgda(K, nums)
%   - A = slgda(K, nums, sol)
%
% $ Arguments $
%   - K:        the kernel gram matrix
%   - nums:     the numbers of samples in each classes
%   - sol:      the cell containing the parameter for generalized eigen
%               decomposition
%   - A:        the resulting projection coefficient matrix
%
% $ Description $
%   - A = slgda(K, nums) performs Generalized Discriminant Analysis(GDA),
%     an representative work in using kernel method to extend LDA,
%     proposed by Baudat et al. The generalized eigen-problem is
%     solved in a default way by slsymgeig.
%
%   - A = slgda(K, nums, sol) in the function, slsymgeig will be invoked
%     to solve the generalized eigen-decomposition problem. sol is
%     a cell containing the parameters for slsymgeig.
%
% $ Remarks $
%   - The function follows the instructions given in the original paper
%     on GDA.
%
%   - The projection is learned after the kernel gram matrix is
%     centralized.
%
%   - The aim of the function is to give an exact implementation of
%     of representative work GDA, so it does not offer other facilities
%     such as weighting and other ways of scatter computation. For higher
%     flexibility, please use the function slkfd.
%
% $ History $
%   - Created by Dahua Lin on May 3rd, 2005
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('slkernelscatter', 2);
end

if ndims(K) ~= 2 || size(K, 1) ~= size(K, 2)
    error('sltoolbox:invaliddims', ...
        'The gram matrix K should be a square matrix');
end

M = size(K, 1);         % number of samples
N = length(nums);       % number of classes
if ~isequal(size(nums), [1, N])
    error('sltoolbox:invaliddims', ...
        'The nums should be a 1 x N row vector');
end
if sum(nums) ~= M
    error('sltoolbox:sizmismatch', ...
        'The total number in nums is inconsistent with that in K');
end

if nargin < 3
    sol = {};
end

%% Compute

%% Centralize

K = slcenkernel(K);

%% Construct the eigen-problem

[sp, ep] = slnums2bounds(nums);
W = zeros(M, M);
for i = 1 : N
    ni = nums(i);
    spi = sp(i); epi = ep(i);
    W(spi:epi, spi:epi) = 1 / ni;
end
clear sp ep;

B = K * W * K;
clear W;

V = K * K;

%% Resolve the eigen-problem

[evs, A] = slsymgeig(B, V, sol{:});
d = sldim_by_eigval(evs);
A = A(:, 1:d);

author: lmthang