text:
function som_trajectory(bmus,varargin)

%SOM_TRAJECTORY Launch a "comet" trajectory visualization GUI.
%
%  som_show(sM,'umat','all')
%  bmus = som_bmus(sM,sD);
%  som_trajectory(bmus)
%  som_trajectory(bmus, 'data1', sD, 'trajsize', [12 6 3 1]')
%  som_trajectory(bmus, 'data1', sD.data(:,[1 2 3]), 'name1', {'fii' 'faa' 'foo'})
%
% Input arguments ([]'s are optional):
%   bmus      (matrix) size Nx1, vector of BMUS
%   ['argID', (string) Other arguments can be given as 'argID', value
%    value]   (varies) pairs. See list below for valid values.
%
% NOTE: the GUI only works on a figure which has been made with SOM_SHOW.
%
% Here are the valid argument IDs (case insensitive) and associated values:
%  'color'      string 'xor' or ColorSpec, default: 'xor'.
%               (default: according to lattice as in som_cplane)
%  'TrajSize'   vector of size Nx1 to define the length of comet
%               (N) and size of the comet dots in points.
%               default: [16 12 10 8 6 4]'
%  'Data1'      SOM Toolbox data struct or matrix. The size of
%               data matrix (in data struct the field .data) is
%               Nxd, where N must be the same as the amount of
%               BMUS given in the first input argument 'bmus'
%               This data is shown in a new window in d subplots.
%               Default: BMU indices (first input argument)
%  'Name1'      cell array of d strings which contains names
%               for the components in 'Data1'. If 'Data1' is a SOM
%               Toolbox data struct, the existing component names
%               are overdone.
%  'Figure'     scalar that must be a handle to an existing figure
%               which has been made using SOM_SHOW function.
%               Default: current active figure (gcf).
%
% The following tools can be found in the 'Tools' -menu.
%
%  Remove Trajectory: removes trajectory from the map.
%  Dye Nodes        : opens GUI for selecting color for the nodes
%                     and points selected.
%  Clear Markers    : removes markers from map and data figure.
%  Save             : saves the current situation as a struct.
%  Load             : loads the struct from workspace and draws markers.
%
% Mouse operation
%
%  In data window: Left button is used to drag the operation point ruler
%                  if left button is on blank area, it starts
%  In map window : Left button starts a polygon; right button
%                  finishes; middle button toggles a unit.
%
% SOM_TRAJECTORY is an application for observing trajectory behavior.
%
% Using mouse the line in data figure can be dragged and the
% trajectory moves in the SOM SHOW figure. It is also possible to move
% trajectory by pressing keys '>' and '<' when mouse pointer is above
% data figure.
%
% Regions can be chosen from the data and the points in that region
% are mapped to the component planes. Regions can be chosen also in
% the map.  In this situation data points and map nodes are also
% marked (Left mouse button adds point to the polygon indicating the
% region and right button finals the polygon). By clicking a node (the
% middle button) that node is either added or removed from selection.
%
% It should be noticed that choosing intervals from data may cause
% situations that seem to be bugs. If there exisist marks of different
% color, removing them by clicking the map may left some marks in the
% data, because more than one point in the data is mapped to the same
% node in the map and the removing operation depends on the color of
% the marks. However, all the marks can be removed by using the 'Clear
% Markers' -operation.
%
% FEATURES
%
% The first input argument 'bmus' may also be a munits x N matrix
% In this case each column defines a "fuzzy response". That is,
% each column defines a hit histogram function). The element
% bmus(i,t) sets the size of marker on unit i at time t.
% NOTE: - in this case no regions can be selcted on the map!
%       - only > and < keys can be used to move the operation point
%         line: it can't be dragged
%       - "fuzzy response is always black (hope this will be fixed)
%
% It is possible to open a second data window showing different data:
% use indetifiers 'Data2' (and 'Name2'). The argument syntax is
% identical to 'Data1' (and 'Name1').
%
% See also SOM_SHOW, SOM_SHOW_ADD, SOM_BMUS.

% Contributed to SOM Toolbox 2.0, February 11th, 2000 by Johan
% Himberg and Juha Parhankangas
% Copyright (c) 2000 by the Johan Himberg and Juha Parhankangas
% http://www.cis.hut.fi/projects/somtoolbox/

% Check arguments

error(nargchk(1,Inf,nargin)); % Check no. of input arguments

%% Init input argument struct (see subfunction)
Traj=iniTraj(bmus);

% Check tentative BMU input validity

if ~vis_valuetype(bmus,{'nxm'}),
  error(['First input should be a vector of BMU indices or' ...
	 ' a "response matrix"']);
end

%% Check optional arguments
for i=1:2:length(varargin)
  identifier=lower(varargin{i});
  value=varargin{i+1};

  % Trajectory color
  switch identifier
   case 'color'
    if isempty(value)
      value='xor';
    end
    if vis_valuetype(value,{'colorstyle','xor'})
      Traj.color=value;
    else
      error('''Color'' has to be ColorSpec or string ''xor''.');
    end

   % 1st data
   case 'data1'
    if isempty(value),
      value=[];
    elseif vis_valuetype(value,{'nxm'})
      Traj.primary_data=value;
    elseif isstruct(value) && isfield(value,'type') && ...
	  ischar(value.type) && strcmp(value.type,'som_data'),
      Traj.primary_data=value.data;
      if isempty(Traj.primary_names),
	Traj.primary_names=value.comp_names;
      end
    end

    % 2nd data
   case 'data2'
    if isempty(value),
      value=[];
    elseif vis_valuetype(value,{'nxm'})
      Traj.secondary_data=value;
    elseif isstruct(value) && isfield(value,'type') && ...
	  ischar(value.type) && strcmp(value.type,'som_data'),
      Traj.secondary_data=value.data;
      if isempty(Traj.secondary_names),
	Traj.secondary_names=value.comp_names;
      end
     end

   % Trajectory length & size
   case 'trajsize'
    if isempty(value),
      Traj.size=[16 12 10 8 6 4]';
    end
    if vis_valuetype(value,{'nx1'})
      Traj.size=value
    else
      error('''TrajSize'' has to be a nx1 vector.');
    end

   % Names for first data variables
   case 'name1'
    if isempty(value),
      Traj.primary_names=[];
    elseif ~vis_valuetype(value,{'cellcolumn_of_char'}),
      error('''Name1'': variable names must be in a cell column array.')
    else
      Traj.primary_names = value;
    end
   % Names for 2nd data variables
   case 'name2'
    if isempty(value),
      Traj.secondary_names=[];
    elseif ~vis_valuetype(value,{'cellcolumn_of_char'}),
      error('''Name2'': variable names must be in a cell column array.')
    else
      Traj.secondary_names = value;
    end

   % Figure number
   case 'figure'
    if isempty(value)
      Traj.figure='gcf';
    end
    if vis_valuetype(value,{'1x1'})
      Traj.figure=value;
    else
      error('''Figure'' should be number of an existing figure.')
    end
  end
end

%% Get SOM data from figure
[h,msg,lattice,msize,dim]=vis_som_show_data('all',Traj.figure);

%% Not a SOM_SHOW figure?
if ~isempty(msg);
  error('Figure is invalid: use SOM_SHOW to draw the figure.');
end

% Get map size from figure data
Traj.lattice=lattice;
Traj.msize=msize;
if length(msize)>2,
  error(['This function works only with 2D maps: figure contains' ...
	 ' something else.']);
end
munits=prod(msize);

% Check BMU (or response) and map match

if vis_valuetype(bmus,{'nx1'});
  if max(bmus)>prod(msize) || min(bmus) <1
    error('BMU indexes out of range.')
  elseif any(round(bmus)~=bmus)
    error('BMU indexes must be integer.');
  elseif isempty(Traj.primary_data),
    Traj.primary_data=bmus;
  end
elseif size(bmus,1) ~= munits
  error(['Response matrix column number must match with the number of' ...
	 ' map units.']);
else
  bmus=bmus';
  if isempty(Traj.primary_data),
    Traj.primary_data=[1:size(bmus,1)]';
    Traj.primary_names={'BMU Index'};
  end
end

size1=size(Traj.primary_data);
size2=size(Traj.secondary_data);

% Data2 must not be defined alone

if isempty(Traj.primary_data)&&~isempty(Traj.secondary_data),
  error('If ''Data2'' is specified ''Data1'' must be specified, too.');
elseif ~isempty(Traj.secondary_data) ...
      && size1~= size2
  % If data1 and data2 exist both, check data1 and data2 match
  error('''Data1'' and ''Data2'' have different amount of data vectors.')
end

% Check BMU and data1 match (data2 matches with 1 anyway)

if ~isempty(Traj.primary_data) && size(bmus,1) ~= size1,
  error(['The number of data vectors in ''data1'' must match with' ...
	 ' the number of rows in the first input argument (bmus).']);
end

% Check that number of names and data dimension is consistent

if ~isempty(Traj.primary_names) && (size1(2)~=length(Traj.primary_names)),
  error('Number of component names and ''Data1'' dimension mismatch.');
end
if ~isempty(Traj.secondary_names) && ...
      (size2(2)~=length(Traj.secondary_names)),
  error('Number of component names and ''Data2'' dimension mismatch.');
end

%% Call the function that does the job
vis_trajgui(Traj);

%%% Subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Traj=iniTraj(bmus)

Traj.figure=gcf;
Traj.primary_data=[];
Traj.secondary_data=[];
Traj.primary_names = [];
Traj.secondary_names = [];
Traj.size=[16 12 10 8 6 4]';
Traj.bmus=bmus;
Traj.color='xor';
Traj.msize=[];
Traj.lattice=[];

author: ilarinieminen

repo: SOM-Toolbox


text:
function calpak_test689 ( )

%*****************************************************************************80
%
%% CALPAK_TEST689 tests YMD_TO_DECIMAL
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2013
%
%  Author:
%
%    John Burkardt
%
  dhi = 1;
  dlo = 1;
  fhi = 0.0;
  flo = 0.0;
  mhi = 1;
  mlo = 1;
  seed = 123456789;
  yhi = 1970;
  ylo = 1960;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST689\n' );
  fprintf ( 1, '  YMD_TO_DECIMAL converts a date to a year and decimal.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMD                            Y.F\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ y, m, d, f, seed ] = ymdf_uniform_common ( ylo, mlo, dlo, flo, yhi, ...
      mhi, dhi, fhi, seed );

    s = ymd_to_s_common ( y, m, d );

    yf = ymd_to_decimal ( y, m, d );

    fprintf ( 1, '  %10s  %14.4f\n', s, yf );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% output = bilateralFilter( data, edge, ...
%                          edgeMin, edgeMax, ...
%                          sigmaSpatial, sigmaRange, ...
%                          samplingSpatial, samplingRange )
%
% Bilateral and Cross-Bilateral Filter using the Bilateral Grid.
%
% Bilaterally filters the image 'data' using the edges in the image 'edge'.
% If 'data' == 'edge', then it the standard bilateral filter.
% Otherwise, it is the 'cross' or 'joint' bilateral filter.
% For convenience, you can also pass in [] for 'edge' for the normal
% bilateral filter.
%
% Note that for the cross bilateral filter, data does not need to be
% defined everywhere.  Undefined values can be set to 'NaN'.  However, edge
% *does* need to be defined everywhere.
%
% data and edge should be of the greyscale, double-precision floating point
% matrices of the same size (i.e. they should be [ height x width ])
%
% data is the only required argument
%
% edgeMin and edgeMax specifies the min and max values of 'edge' (or 'data'
% for the normal bilateral filter) and is useful when the input is in a
% range that's not between 0 and 1.  For instance, if you are filtering the
% L channel of an image that ranges between 0 and 100, set edgeMin to 0 and
% edgeMax to 100.
%
% edgeMin defaults to min( edge( : ) ) and edgeMax defaults to max( edge( : ) ).
% This is probably *not* what you want, since the input may not span the
% entire range.
%
% sigmaSpatial and sigmaRange specifies the standard deviation of the space
% and range gaussians, respectively.
% sigmaSpatial defaults to min( width, height ) / 16
% sigmaRange defaults to ( edgeMax - edgeMin ) / 10.
%
% samplingSpatial and samplingRange specifies the amount of downsampling
% used for the approximation.  Higher values use less memory but are also
% less accurate.  The default and recommended values are:
%
% samplingSpatial = sigmaSpatial
% samplingRange = sigmaRange
%
%

% Copyright (c) <2007> <Jiawen Chen, Sylvain Paris, and Fredo Durand>
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
%

function output = bilateralFilter( data, edge, edgeMin, edgeMax, sigmaSpatial, sigmaRange, samplingSpatial, samplingRange )

if( ndims( data ) > 2 ),
    error( 'data must be a greyscale image with size [ height, width ]' );
end

if( ~isa( data, 'double' ) ),
    error( 'data must be of class "double"' );
end

if ~exist( 'edge', 'var' ),
    edge = data;
elseif isempty( edge ),
    edge = data;
end

if( ndims( edge ) > 2 ),
    error( 'edge must be a greyscale image with size [ height, width ]' );
end

if( ~isa( edge, 'double' ) ),
    error( 'edge must be of class "double"' );
end

inputHeight = size( data, 1 );
inputWidth = size( data, 2 );

if ~exist( 'edgeMin', 'var' ),
    edgeMin = min( edge( : ) );
    %warning( 'edgeMin not set!  Defaulting to: %f\n', edgeMin );
end

if ~exist( 'edgeMax', 'var' ),
    edgeMax = max( edge( : ) );
    %warning( 'edgeMax not set!  Defaulting to: %f\n', edgeMax );
end

edgeDelta = edgeMax - edgeMin;

if ~exist( 'sigmaSpatial', 'var' ),
    %sigmaSpatial = min( inputWidth, inputHeight ) / 16;
    sigmaSpatial = min( inputWidth, inputHeight ) / 64;
    fprintf( 'Using default sigmaSpatial of: %f\n', sigmaSpatial );
end

if ~exist( 'sigmaRange', 'var' ),
    %sigmaRange = 0.1 * edgeDelta;
    sigmaRange = 0.025 * edgeDelta;
    fprintf( 'Using default sigmaRange of: %f\n', sigmaRange );
end

if ~exist( 'samplingSpatial', 'var' ),
    samplingSpatial = sigmaSpatial;
end

if ~exist( 'samplingRange', 'var' ),
    samplingRange = sigmaRange;
end

if size( data ) ~= size( edge ),
    error( 'data and edge must be of the same size' );
end

% parameters
derivedSigmaSpatial = sigmaSpatial / samplingSpatial;
derivedSigmaRange = sigmaRange / samplingRange;

paddingXY = floor( 2 * derivedSigmaSpatial ) + 1;
paddingZ = floor( 2 * derivedSigmaRange ) + 1;

% allocate 3D grid
downsampledWidth = floor( ( inputWidth - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
downsampledHeight = floor( ( inputHeight - 1 ) / samplingSpatial ) + 1 + 2 * paddingXY;
downsampledDepth = floor( edgeDelta / samplingRange ) + 1 + 2 * paddingZ;

gridData = zeros( downsampledHeight, downsampledWidth, downsampledDepth );
gridWeights = zeros( downsampledHeight, downsampledWidth, downsampledDepth );

% compute downsampled indices
[ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 );

% ii =
% 0 0 0 0 0
% 1 1 1 1 1
% 2 2 2 2 2

% jj =
% 0 1 2 3 4
% 0 1 2 3 4
% 0 1 2 3 4

% so when iterating over ii( k ), jj( k )
% get: ( 0, 0 ), ( 1, 0 ), ( 2, 0 ), ... (down columns first)

di = round( ii / samplingSpatial ) + paddingXY + 1;
dj = round( jj / samplingSpatial ) + paddingXY + 1;
dz = round( ( edge - edgeMin ) / samplingRange ) + paddingZ + 1;

% perform scatter (there's probably a faster way than this)
% normally would do downsampledWeights( di, dj, dk ) = 1, but we have to
% perform a summation to do box downsampling
for k = 1 : numel( dz ),

    dataZ = data( k ); % traverses the image column wise, same as di( k )
    if ~isnan( dataZ  ),

        dik = di( k );
        djk = dj( k );
        dzk = dz( k );

        gridData( dik, djk, dzk ) = gridData( dik, djk, dzk ) + dataZ;
        gridWeights( dik, djk, dzk ) = gridWeights( dik, djk, dzk ) + 1;

    end
end

% make gaussian kernel
kernelWidth = 2 * derivedSigmaSpatial + 1;
kernelHeight = kernelWidth;
kernelDepth = 2 * derivedSigmaRange + 1;

halfKernelWidth = floor( kernelWidth / 2 );
halfKernelHeight = floor( kernelHeight / 2 );
halfKernelDepth = floor( kernelDepth / 2 );

[gridX, gridY, gridZ] = meshgrid( 0 : kernelWidth - 1, 0 : kernelHeight - 1, 0 : kernelDepth - 1 );
gridX = gridX - halfKernelWidth;
gridY = gridY - halfKernelHeight;
gridZ = gridZ - halfKernelDepth;
gridRSquared = ( gridX .* gridX + gridY .* gridY ) / ( derivedSigmaSpatial * derivedSigmaSpatial ) + ( gridZ .* gridZ ) / ( derivedSigmaRange * derivedSigmaRange );
kernel = exp( -0.5 * gridRSquared );

% convolve
blurredGridData = convn( gridData, kernel, 'same' );
blurredGridWeights = convn( gridWeights, kernel, 'same' );

% divide
blurredGridWeights( blurredGridWeights == 0 ) = -2; % avoid divide by 0, won't read there anyway
normalizedBlurredGrid = blurredGridData ./ blurredGridWeights;
normalizedBlurredGrid( blurredGridWeights < -1 ) = 0; % put 0s where it's undefined

% for debugging
% blurredGridWeights( blurredGridWeights < -1 ) = 0; % put zeros back

% upsample
[ jj, ii ] = meshgrid( 0 : inputWidth - 1, 0 : inputHeight - 1 ); % meshgrid does x, then y, so output arguments need to be reversed
% no rounding
di = ( ii / samplingSpatial ) + paddingXY + 1;
dj = ( jj / samplingSpatial ) + paddingXY + 1;
dz = ( edge - edgeMin ) / samplingRange + paddingZ + 1;

% interpn takes rows, then cols, etc
% i.e. size(v,1), then size(v,2), ...
output = interpn( normalizedBlurredGrid, di, dj, dz );

author: jianxiongxiao

repo: ProfXkit


text:
function [data_matrix12, xvals2, yvals2, data_matrix22, bestx, besty, bestz] = surf_plot_tor(data_matrix1, xvals, yvals, xname, yname, zname, varargin)
% Stylized surface plot of one or two surfaces
%
% :Usage:
% ::
%
%    [data_matrix12, xvals2, yvals2, data_matrix22] = surf_plot_tor(data_matrix1, xvals, yvals, xname, yname, zname, [data_matrix2])
%
% :Examples: See classify_search_script3...in meta-analysis classification
% ::
%
%    surf_plot_tor(corrc_mean, mya, mys, 'Activation feature cutoff', 'Sensitivity feature cutoff', 'Classification accuracy', worstcat)
%
% xvals are columns, yvals are rows!
%
% ..
%    Tor Wager, Sept. 07
% ..

    data_matrix22 = [];

    ystep = 1/5 .* range(yvals) ./ length(yvals);
    xstep = 1/5 .* range(xvals) ./ length(xvals);

    [X,Y] = meshgrid(xvals,yvals);

%     if size(X, 1) == size(data_matrix1, 2)
%         error('Wrong sizes...are the x and y inputs flipped?');
%     end

    yvals2 = 0:ystep:max(yvals);
    xvals2 = 0:xstep:max(xvals);

    [X2,Y2] = meshgrid(xvals2,yvals2);
    data_matrix12 = interp2(X,Y,data_matrix1,X2,Y2,'cubic');

    create_figure('Surface plot');

    han = surf(X2,Y2,data_matrix12);

    xlabel('Activation feature cutoff')
    ylabel('Sensitivity feature cutoff')
    zlabel('Classification accuracy');

    xlabel(xname);
    ylabel(yname);
    zlabel(zname);

    set(han,'EdgeColor','none')
    grid on

    % lines
    plot3(X2(end,:),Y2(end,:),data_matrix12(end,:),'r','LineWidth',2);
    plot3(X2(end,:),Y2(end,:),data_matrix12(end,:),'r','LineWidth',3);

    for i = 10:10:size(X2,1)
        plot3(X2(i,:),Y2(i,:),data_matrix12(i,:),'Color',[.5 .5 .5],'LineWidth',1);
    end

    view(135, 30)
    drawnow

    scn_export_papersetup(600);

    if length(varargin) > 0
        data_matrix2 = varargin{1};

        data_matrix22 = interp2(X,Y,data_matrix2,X2,Y2,'cubic');
        hold on
        han2 = surf(X2,Y2,data_matrix22);
        set(han2,'EdgeColor','none')
        set(han,'FaceAlpha',.7)

        % lines
        plot3(X2(end,:),Y2(end,:),data_matrix22(end,:),'b','LineWidth',3);
        plot3(X2(1,:),Y2(1,:),data_matrix22(1,:),'k','LineWidth',1);

        for i = 10:10:size(X2,1)
            plot3(X2(i,:),Y2(i,:),data_matrix22(i,:),'Color',[.5 .5 .5],'LineWidth',1);
        end

    end

    % calc max and put max point on map
    % --------------------------------------------------
    if length(varargin) > 0
        % average two maps
        mymap = .5 * (data_matrix12 + data_matrix22);
        fprintf('Using average of two maps to calculate max.');
    else
        mymap = data_matrix12;
    end

    [bestsum] = max(mymap(:));
    [row, col] =  find(mymap == bestsum);
bestx = xvals2(col);
besty = yvals2(row);
bestz = mymap(row, col);

fprintf('Maximum: x = %3.3f, y = %3.3f, z = %3.3f\n', bestx, besty, bestz);

plotz = bestz;
if length(varargin) > 0
    max1 = data_matrix12(row, col);
    max2 = data_matrix22(row, col);
    plotz = max([max1 max2]);

    bestz = [bestz max1 max2];
    fprintf('Height at overall max: Map 1: %3.3f, Map 2: %3.3f\n', max1, max2);
end

% plot
plot3(bestx, besty, plotz, 'ko', 'MarkerSize', 10, 'MarkerFaceColor', 'k')

mylowerbound = get(gca,'ZLim'); mylowerbound = mylowerbound(1);
plot3([bestx bestx], [besty besty], [plotz mylowerbound], 'k-', 'LineWidth', 2);

myxlowerbound = get(gca,'XLim'); %myxlowerbound = myxlowerbound(1);
plot3([myxlowerbound], [besty besty], [mylowerbound mylowerbound], 'k-', 'LineWidth', 2);

myylowerbound = get(gca,'YLim'); %myylowerbound = myylowerbound(1);
plot3([bestx bestx], [myylowerbound], [mylowerbound mylowerbound], 'k-', 'LineWidth', 2);

set(gca, 'ZLim', get(gca, 'ZLim'));
axis vis3d

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function Y=repmat(varargin)
%REPMAT (overloaded)

try
  X = varargin{1};
  Y = X;
  Y.basis = [];
  n = Y.dim(1);
  m = Y.dim(2);
  for i = 1:length(Y.lmi_variables)+1
    temp = repmatfixed(reshape(X.basis(:,i),n,m),varargin{2:end});
    Y.basis(:,i) = temp(:);
  end
  Y.dim(1) = size(temp,1);
  Y.dim(2) = size(temp,2);
  % Reset info about conic terms
  Y.conicinfo = [0 0];
catch
  error(lasterr)
end

function B = repmatfixed(A,M,N)

if nargin < 2
    error('MATLAB:repmat:NotEnoughInputs', 'Requires at least 2 inputs.')
end

if nargin == 2
    if isscalar(M)
        siz = [M M];
    else
        siz = M;
    end
else
    siz = [M N];
end

if isscalar(A)
    nelems = prod(siz);
    if nelems>0
        % Since B doesn't exist, the first statement creates a B with
        % the right size and type.  Then use scalar expansion to
        % fill the array. Finally reshape to the specified size.
        B = spalloc(nelems,1,nnz(A));
        B(nelems) = A;
        if ~isequal(B(1), B(nelems)) | ~(isnumeric(A) | islogical(A))
            % if B(1) is the same as B(nelems), then the default value filled in for
            % B(1:end-1) is already A, so we don't need to waste time redoing
            % this operation. (This optimizes the case that A is a scalar zero of
            % some class.)
            B(:) = A;
        end
        B = reshape(B,siz);
    else
        B = A(ones(siz));
    end
elseif ndims(A) == 2 & numel(siz) == 2
    [m,n] = size(A);

    if (m == 1 & siz(2) == 1)
        B = A(ones(siz(1), 1), :);
    elseif (n == 1 & siz(1) == 1)
        B = A(:, ones(siz(2), 1));
    else
        mind = (1:m)';
        nind = (1:n)';
        mind = mind(:,ones(1,siz(1)));
        nind = nind(:,ones(1,siz(2)));
        B = A(mind,nind);
    end
else
    Asiz = size(A);
    Asiz = [Asiz ones(1,length(siz)-length(Asiz))];
    siz = [siz ones(1,length(Asiz)-length(siz))];
    for i=length(Asiz):-1:1
        ind = (1:Asiz(i))';
        subs{i} = ind(:,ones(1,siz(i)));
    end
    B = A(subs{:});
end

function a = isscalar(b)
[n,m] = size(b);
a = (n*m == 1);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%CUSTOMEQN Custom equation physics mode definition.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym rhs (@var{f})
%% Right-hand side of symbolic expression.
%%
%% Example:
%% @example
%% @group
%% syms x
%% eqn = 5*x <= 3*x + 6
%%   @result{} eqn = (sym) 5⋅x ≤ 3⋅x + 6
%% rhs(eqn)
%%   @result{} ans = (sym) 3⋅x + 6
%% @end group
%% @end example
%%
%% Gives an error if any of the symbolic objects have no right-hand side.
%%
%% @seealso{@@sym/lhs, @@sym/children, @@sym/formula, @@sym/argnames}
%% @end defmethod

function R = rhs(f)

  R = elementwise_op ('lambda a: a.rhs', f);

end

%% most tests are in lhs
%!test
%! syms x
%! f = x + 1 == 2*x;
%! assert (isequal (rhs(f), 2*x))

%!error <AttributeError>
%! syms x
%! rhs(x)

repo: octsympy


text:
function [xUpdate,PUpdate,innov,Pzz,W]=EKFUpdate(xPred,PPred,z,R,h,HJacob,numIter,HHessian,innovTrans,measPredTrans,stateDiffTrans,stateTrans)
%%EKFUPDATE Perform the measurement update step in the first- or second-
%           order Extended Kalman Filter (EKF) or iterated Extended Kalman
%           filter (IEKF), which is what one gets with numIter>0.
%
%INPUTS: xPred The xDimX1 predicted target state.
%        PPred The xDimXxDim predicted state covariance matrix.
%            z The zDimX1 vector measurement.
%            R The zDimXzDim measurement covariance matrix.
%            h A function handle for the measurement function that takes
%              the state as its argument.
%       HJacob A function handle for the measurement Jacobian matrix that
%              takes the target state as a parameter or the zDimXxDim
%              measurement Jacobian matrix itself. Each column is the
%              derivative vector of h with respect to the corresponding
%              element of x. If not supplied or an empty matrix is passed,
%              then HJacob will be found using numerical differentiation
%              via the numDiff function with default parameters.
%      numIter The number of iterations to perform if an iterated EKF is
%              desired. If this parameter is omitted or an empty matrix is
%              passed, the default value of zero is used. That is, just
%              use the standard update without any additional iterations.
%     HHessian This parameter is only provided if a second-order EKF is
%              desired. This is either a function handle for the
%              measurement Hessian hypermatrix, or it is the measurement
%              Hessian hypermatrix itself. The matrix is
%              xDimXxDimXzDim. The matrix HHess=HHessian(x) is such that
%              HHess(i,j,k) is the second derivative of the kth element
%              of the vector returned by h with respect to the ith and jth
%              components of x. For each k, the Hessian matrix is
%              symmetric. If this parameter is omitted or an empty matrix
%              is passed, a first-order EKF update is used.
%   innovTrans An optional function handle that computes and optionally
%              transforms the value of the difference between the
%              observation and any predicted points. This is called as
%              innovTrans(a,b) and the default if omitted or an empty
%              matrix is passed is @(a,b)bsxfun(@minus,a,b). This must be
%              able to handle sets of values. For a zDimX1 measurement,
%              either of the inputs could be zDimXN in size while one of
%              the inputs could be zDimX1 in size.  This only needs to be
%              supplied when a measurement difference must be restricted
%              to a certain range. For example, the innovation between two
%              angles will be 2*pi if one angle is zero and the other
%              2*pi, even though they are the same direction. In such an
%              instance, a function handle to the
%              wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%              appropriate parameters should be passed for innovTrans.
% measPredTrans An optional function handle that transforms the predicted
%              measurement into a particular range. The second order EKF
%              has a linear correction to the prediction of the
%              measurement, which might not be appropriate for all types of
%              measurements.
% stateDiffTrans An optional function handle that, like innovTrans does for
%              the measurements, a difference between states and transforms
%              it however might be necessary. For example, a state
%              containing angular components will generally need to be
%              transformed so that the difference between the angles is
%              wrapped to -pi/pi.
%   stateTrans An optional function handle that takes a state estimate and
%              transforms it. This is useful if one wishes the elements of
%              the state to be bound to a certain domain. For example, if
%              an element of the state is an angle, one should generally
%              want to bind it to the region +/-pi.
%
%OUTPUTS: xUpdate The xDim X 1 updated state vector.
%         PUpdate The updated xDim X xDim state covariance matrix.
%      innov, Pzz The zDimX1 innovation and the zDimXzDim innovation
%                 covariance matrix are returned in case one wishes to
%                 analyze the consistency of the estimator or use those
%                 values in gating or likelihood evaluation.
%               W The xDimXzDim gain used in the update. This can be
%                 useful when gating and using the function
%                 calcMissedGateCov.
%
%The first-order EKF is summarized in Figure 10.3.3-1 in Chapter 10.3.3 of
%[1]. The Joseph-form covariance update given in Chapter 5 of the same
%book is used for improved numerical stability. The expressions for the
%second-order EKF come from Chapter 10.3.2 of [1].
%
%The iteration of the measurement equation for the iterated EKF is given in
%Chapter 10.5.2 of the same book. In [2], it is noted that the iteration
%need only be done over each measurement update, not over an entire batch
%of measurements.
%
%The optional parameter innovTrans is not described in the above reference,
%but allow for possible modifications to the filter as described in [3].
%The parameters have been added to allow the filter to be used with
%angular quantities. For example, if the measurement consisted of
%range and angle, z=[r;theta], then
%innovTrans=@(a,b)[bsxfun(@minus,a(1,:),b(1,:));
%                  wrapRange(bsxfun(@minus,a(2,:),b(2,:)),-pi,pi)];
%should be used to approximately deal with the circular nature of the
%measurements.
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%[2] T. H. Kerr, "Streamlining Measurement Iteration for EKF Target
%    Tracking," IEEE Transactions on Aerospace and Electronic Systems,
%    vol. 27, no. 2, pp. 408-421, Mar. 1991.
%[3] D. F. Crouse, "Cubature/ unscented/ sigma point Kalman filtering with
%    angular measurement models," in Proceedings of the 18th International
%    Conference on Information Fusion, Washington, D.C., 6-9 Jul. 2015.
%
%April 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

zDim=size(z,1);
if(nargin<6||isempty(HJacob))
    HJacob=@(x)numDiff(x,h,zDim);
end

if(nargin<7||isempty(numIter))
    numIter=0;
end

if(nargin<8)
    HHessian=[];
end

if(nargin<9||isempty(innovTrans))
    %The function just returns the input.
    innovTrans=@(a,b)bsxfun(@minus,a,b);
end

if(nargin<10||isempty(measPredTrans))
    measPredTrans=@(x)x;
end

if(nargin<11||isempty(stateDiffTrans))
    stateDiffTrans=@(x)x;
end

if(nargin<12||isempty(stateTrans))
    stateTrans=@(x)x;
end

if(isa(HJacob,'function_handle'))
    H=HJacob(xPred);
else
    H=HJacob;%If the Jacobian matrix was directly given.
end

if(~isempty(HHessian))
    if(isa(HHessian,'function_handle'))
        HHess=HHessian(xPred);
    else
        HHess=HHessian;%If the Hessian matrix was directly given.
    end

    [zPredHessTerm,PzzHessTerm]=getHessianTerms(HHess,PPred,zDim);
else
    zPredHessTerm=0;
    PzzHessTerm=0;
end

zPred=measPredTrans(h(xPred)+zPredHessTerm);
innov=innovTrans(z,zPred);

Pzz=R+H*PPred*H'+PzzHessTerm;
%Ensure symmetry
Pzz=(Pzz+Pzz')/2;
W=PPred*H'/Pzz;

xUpdate=stateTrans(xPred+W*innov);

temp=W*H;
temp=eye(size(temp))-temp;
PUpdate=temp*PPred*temp'+W*R*W';
%Ensure symmetry
PUpdate=(PUpdate+PUpdate')/2;

for curIter=1:numIter
    if(isa(HJacob,'function_handle'))
        H=HJacob(xUpdate);
    else
        H=HJacob;%If the Jacobian matrix was directly given.
    end

    if(~isempty(HHessian))
        if(isa(HHessian,'function_handle'))
            HHess=HHessian(xUpdate);
        else
            HHess=HHessian;%If the hessian matrix was directly given.
        end

        [zPredHessTerm,PzzHessTerm]=getHessianTerms(HHess,PPred,zDim);
    else
        zPredHessTerm=0;
        PzzHessTerm=0;
    end

    %Update the PEstimate
    Pzz=R+H*PPred*H'+PzzHessTerm;
    %Ensure symmetry
    Pzz=(Pzz+Pzz')/2;
    W=PPred*H'/Pzz;
    temp=W*H;
    temp=eye(size(temp))-temp;
    PUpdate=temp*PPred*temp'+W*R*W';
    %Ensure symmetry
    PUpdate=(PUpdate+PUpdate')/2;
    zPred=measPredTrans(h(xUpdate)+zPredHessTerm);

    %Update the x estimate
    xUpdate=stateTrans(xUpdate+PUpdate*H'*lsqminnorm(R,innovTrans(z,zPred))-PUpdate*lsqminnorm(PPred,stateDiffTrans(xUpdate-xPred)));
end

end

function [zPredHessTerm,PzzHessTerm]=getHessianTerms(HHess,PPred,zDim)
%The function returns the terms associated with the second derivative in
%the second-order EKF.

PzzHessTerm=zeros(zDim,zDim);
zPredHessTerm=zeros(zDim,1);
for n=1:zDim
    en=zeros(zDim,1);
    en(n)=1;

    HPProdn=HHess(:,:,n)*PPred;
    zPredHessTerm=zPredHessTerm+en*trace(HPProdn);

    for m=1:zDim
        em=zeros(zDim,1);
        em(m)=1;
        PzzHessTerm=PzzHessTerm+en*em'*trace(HPProdn*HHess(:,:,m)'*PPred);
    end
end
PzzHessTerm=PzzHessTerm/2;
zPredHessTerm=zPredHessTerm/2;
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:
% remesh    Surface simplification
%
% Usage:
%
%  [surface]   = remesh (surface, opt)
%  [surface]   = remesh (TRI, X, Y, Z, opt)
%  [TRI,X,Y,Z] = remesh (surface, opt)
%  [TRI,X,Y,Z] = remesh (TRI, X, Y, Z, opt)
%
% Description:
%
%  Reduces the vertex and face count of a surface using the QSlim
%  algorithm.
%
% Input:
%
%  TRIV    - ntx3 triangulation matrix with 1-based indices (as the one
%            returned by the MATLAB function delaunay).
%  X,Y,Z   - vectors with nv vertex coordinates.
%  surface - alternative way to specify the mesh as a struct, having .TRIV,
%            .X, .Y, and .Z as its fields.
%  opt       - (optional) settings
%              .placement  [-O] Optimal placement policy (default: 3)
%                                0 - end points,
%                                1 - end or mid points,
%                                2 - line,
%                                3 - optimal
%              .boundary   [-B] Use boundary preservation planes with given weight
%                               (default: 1000)
%              .weight     [-W] Quadric weighting policy (default: 1)
%                                0 - uniform,
%                                1 - area,
%                                2 - angle
%              .contract   [-F] Contraction of faces ('faces') or edges ('edges')
%                               (default: 'edges')
%              .penalty    [-m] Penalty for bad meshes (default: 1)
%              .compact    [-c] Compactness ratio (default: 0)
%              .join       [-j] Join only without removing faces (default: 'false')
%              .faces      [-t] Target number of faces
%              .vertices        Target number of vertices (default: 2000)
%              .verbose         Verbosity level (default: 1)
%                                0 - no display,
%                                1 - display results
%
% Output:
%
%  TRIV    - ntx3 triangulation matrix of the simplified surface.
%  X,Y,Z   - list of vertices of the simplified surface.
%  surface - alternative way to specify the simplified surface.
%
% References:
%
% [1] http://www.cs.cmu.edu/~garland/quadrics/
%
% TOSCA = Toolbox for Surface Comparison and Analysis
% Web: http://tosca.cs.technion.ac.il
% Version: 0.9
%
% (C) QSlim 2.0 copyright Michael Garland, 1998.
% (C) Mex iterface copyright Alex Bronstein, 2007.

function [TRIV_, X_, Y_, Z_] = remesh(surface, options)

if isstruct(surface),
    TRI = surface.TRIV;
    X = surface.VERT(:,1);
    Y = surface.VERT(:,2);
    Z = surface.VERT(:,3);
else
    TRI = surface;
    X = varargin{1};
    Y = varargin{2};
    Z = varargin{3};
    varargin = varargin(4:end);
end

opt = [3 1000 1 0 1 0 0];
% verts = min(length(X), 2000);
verts = min(length(X), 1e4);
faces = 0;
verbose = 0;

if nargin > 1,
    if isfield(options, 'placement'),
        opt(1) = options.placement;
    end
    if isfield(options, 'boundary'),
        opt(2) = options.boundary;
    end
    if isfield(options, 'weight'),
        opt(3) = options.weight;
    end
    if isfield(options, 'contract'),
        if strcmpi(value,'face') | strcmpi(value,'faces'),
            opt(4) = 1;
        elseif strcmpi(value,'edge') | strcmpi(value,'edges'),
            opt(4) = 0;
        else
            error('Invalid value setting for .contract. Must be "faces" or "edges".');
        end
    end
    if isfield(options, 'penalty'),
        opt(5) = options.penalty;
    end
    if isfield(options, 'compact'),
        opt(6) = options.compact;
    end
    if isfield(options, 'join'),
        if strcmpi(options.join,'true') | value == 1,
            opt(7) = 1;
        elseif strcmpi(options.join,'false') | value == 0,
            opt(7) = 0;
        else
            error('Invalid value setting for .join. Must be "true" or "false".');
        end
    end
    if isfield(options, 'vertices'),
        verts = options.vertices;
    end
    if isfield(options, 'faces'),
        faces = options.faces;
    end
    if isfield(options, 'verbose'),
        verbose = options.verbose;
    end

end

if verts < 1 & faces < 1,
    error('Either face or vertex positive count targets have to be specified.');
end

switch opt(1),
    case 0, placement_policy = 'end point';
    case 1, placement_policy = 'end or mid point';
    case 2, placement_policy = 'line';
    case 3, placement_policy = 'optimal';
    otherwise, error('Invalid placement policy. Use .placement = 0-3.');
end

switch opt(3),
    case 0, weighting_policy = 'uniform';
    case 1, weighting_policy = 'area';
    case 2, weighting_policy = 'angle';
    otherwise, error('Invalid weighting policy. Use .weigh = 0-2.');
end

if opt(4), contract = 'face', else contract = 'edge'; end
if opt(7), join = 'on', else join = 'off'; end

if verbose > 0,
    fprintf (1, 'Surface simplification\n');
    fprintf (1, 'Placement policy:         %s\n', placement_policy);
    fprintf (1, 'Weighting policy:         %s\n', weighting_policy);
    fprintf (1, 'Boundary preserv. weight: %-8.6g\n', opt(2));
    fprintf (1, 'Contraction:              %s\n', contract);
    fprintf (1, 'Bad mesh penalty:         %-8.6g\n', opt(5));
    fprintf (1, 'Compactness ratio:        %-8.6g\n', opt(6));
    fprintf (1, 'Join only:                %s\n', join);
    fprintf (1, 'Original faces:           %-d\n', size(TRI,1));
    fprintf (1, 'Original vertices:        %-d\n', length(X));
    fprintf (1, 'Target faces:             %-d\n', faces);
    fprintf (1, 'Target vertices:          %-d\n', verts);
end

[TRIV_,X_,Y_,Z_] = qslim(TRI, X, Y, Z, faces, verts, opt);
if verbose > 0,
    fprintf (1, 'Done. faces = %d vertices = %d\n\n', size(TRIV_,1), length(X_));
end

idx = find( X_.^2 + Y_.^2 + Z_.^2 > 0 );
N = max(idx);
X_ = X_(1:N);
Y_ = Y_(1:N);
Z_ = Z_(1:N);

if nargout <= 1,
    surface_ = [];
    surface_.TRIV = TRIV_;
    surface_.VERT    = [X_ Y_ Z_];
    surface_.n = size(surface_.VERT,1);
    surface_.m = size(surface_.TRIV,1);
    TRIV_ = surface_;
    if isstruct(surface) & isfield(surface,'D'),
        warning('The returned result is an uninitialized surface. Use init_surface to initialize.');
    end
end

author: OshriHalimi


text:
function [ abd, info ] = zpbfa ( abd, lda, n, m )

%*****************************************************************************80
%
%% ZPBFA factors a complex hermitian positive definite band matrix.
%
%  Discussion:
%
%    ZPBFA is usually called by ZPBCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%
%  Band storage:
%
%    If A is a hermitian positive definite band matrix,
%    the following program segment will set up the input.
%
%      m = (band width above diagonal)
%      do j = 1, n
%        i1 = max ( 1, j-m )
%        do i = i1, j
%          k = i-j+m+1
%          abd(k,j) = a(i,j)
%        end do
%      end do
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex ABD(LDA,N); the matrix to be factored.
%    The columns of the upper triangle are stored in the columns of ABD
%    and the diagonals of the upper triangle are stored in the rows of ABD.
%
%    Input, integer LDA, the leading dimension of ABD.
%    LDA must be at least M+1.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer M, the number of diagonals above the main diagonal.
%    0 <= M < N.
%
%    Output, integer INFO.
%    0, for normal return.
%    K, if the leading minor of order K is not positive definite.
%
%    Output, complex ABD(LDA,N); an upper triangular matrix R, stored in
%    band form, so that A = hermitian(R)*R.
%
  info = 0;

  for j = 1 : n

    s = 0.0;
    ik = m + 1;
    jk = max ( j - m, 1 );
    mu = max ( m + 2 - j, 1 );

    for k = mu : m
      t = abd(k,j) - abd(ik:ik+k-mu-1,jk)' * abd(mu:mu+k-mu-1,j);
      t = t / abd(m+1,jk);
      abd(k,j) = t;
      s = s + real ( t * conj ( t ) );
      ik = ik - 1;
      jk = jk + 1;
    end

    s = real ( abd(m+1,j) ) - s;

    if ( s <= 0.0 | imag ( abd(m+1,j) ) ~= 0.0 )
      info = j;
      break
    end

    abd(m+1,j) = sqrt ( s );

  end

  return
end


text:
function Problem2_26
%Problem 2.26 Ternary condensation inside a vertical tube.
% Uses the function BVP4C to numerically solve the boundary-value problem.
% BVPINIT is used to form an initial guess for the solution on
% a mesh of ten equally spaced points. A guess for unknown parameters
% (the two fluxes) is the last argument of BVPINIT.
% BVP4C returns the solution as the structure 'sol'. The computed fluxes are
% returned in the field sol.parameters.
% The calculated concentration profiles are plotted.

% Calculate the mass-transfer coefficients
Re = 9574;
density = 0.882;
viscosity = 0.00001602;
diameter = 0.0254;
% MS diffusion coefficients
D = [1.6*10^-5 1.444*10^-5 3.873*10^-5];
% Calculate Schmidt numbers
Sc(1) = viscosity/(density*D(1));
Sc(2) = viscosity/(density*D(2));
Sc(3) = viscosity/(density*D(3));
% Calculate Sherwood numbers
Sh(1) = 0.023*Re^0.83*Sc(1)^0.44;
Sh(2) = 0.023*Re^0.83*Sc(2)^0.44;
Sh(3) = 0.023*Re^0.83*Sc(3)^0.44;
% Molecular weight of pure components
M = [60.1 18 28];
y1 = [0.1123
      0.4246
      0.4631];
Mav1 = M*y1;
y2 = [0.1457
      0.1640
      0.6903];
Mav2 = M*y2;
Mav = (Mav1+Mav2)/2;
% Calculate total molar density
c = density/Mav;
F12 = Sh(1)*c*D(1)*1000/diameter;
F13 = Sh(2)*c*D(2)*1000/diameter;
F23 = Sh(3)*c*D(3)*1000/diameter;
solinit = bvpinit(linspace(0,1,10),[1, 1], [1, 1]);
options = bvpset('stats','on');
sol = bvp4c(@prob226ode,@prob226bc,solinit,options,F12,F13,F23);
x = sol.x;
y = sol.y;
y(3,:) = 1-y(1,:) -y(2,:);
fprintf('\n');
fprintf('Computed Mass-Transfer Coefficients:');
fprintf('\n');
fprintf('F12, mole/m2-s = %7.4f.\n',F12);
fprintf('\n');
fprintf('F13, mole/m2-s = %7.4f.\n',F13);
fprintf('\n');
fprintf('F23, mole/m2-s = %7.4f.\n',F23);
fprintf('\n');
fprintf(' Computed N1 in mole/m2-s =%7.4f.\n',sol.parameters(1))
fprintf('\n')
fprintf(' Computed N2 in mole/m2-s=%7.4f.\n',sol.parameters(2))
fprintf('\n')
clf reset
plot(x,y(1,:),'k-*',x,y(2,:),'k:+',x,y(3,:),'k-.o')
axis([0 1 0 0.7])
title('Concentration profiles in Problem 2.26')
xlabel('Distance along the diffusion path, cm.')
ylabel('Mole fraction')
legend('Isopropanol', 'Water', 'Nitrogen')
shg

% --------------------------------------------------------------------------

function dydx = prob226ode(x,y,N,F12,F13,F23)
N3 = 0;

%F12 = 0.933;
%F13 = 0.880;
%F23 = 1.530;
dydx = [ (y(1)*N(2)-y(2)*N(1))/F12+(y(1)*N3-(1-y(1)-y(2))*N(1))/F13
         (y(2)*N(1)-y(1)*N(2))/F12+(y(2)*N3-(1-y(1)-y(2))*N(2))/F23 ];

% --------------------------------------------------------------------------

function res = prob226bc(ya,yb,N,F12,F13,F23)
res = [ya(2)-0.4246
       yb(2) - 0.1640
       ya(1) - 0.1123
       yb(1) - 0.1457 ];

% --------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function imageHandle = prtPlotUtilPlotGriddedEvaledClassifier(DS, linGrid, gridSize, cMap, prtActor)
% Internal function
% xxx Need Help xxx

% Copyright (c) 2013 New Folder Consulting
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.

% Check to see if we are plotting a decision
if ~isempty(prtActor.internalDecider)
    % When using a decider DS contains integers and the linear spacing
    % applied using the colormap may not be correct (depending on the
    % specific integer labels). Therefore we have to change DS to be
    % consecutive integers corresponding to the integer labels that were
    % used.

    [dontNeed, DS] = ismember(DS,prtActor.internalDecider.classList); %#ok<ASGLU>

%     if ~isempty(prtActor.internalDecider.dataSetSummary.uniqueClasses)
%         [dontNeed, DS] = ismember(DS,prtActor.internalDecider.dataSetSummary.uniqueClasses); %#ok<ASGLU>
%     else
%         % Have something unsupervised like a clusterer
%         % Just roll with it.
%     end
end
nDims = size(linGrid,2);

switch nDims
    case 1
        DS = repmat(DS,1,max(prtActor.dataSetSummary.uniqueClasses));
        imageHandle = imagesc(linGrid,[min(prtActor.dataSetSummary.uniqueClasses),max(prtActor.dataSetSummary.uniqueClasses)]',DS');
        set(gca,'YTickLabel',[])
        colormap(cMap)
    case 2
        % find decision boundary, center colormap around it; THIS CODE IS
        % BROKEN FOR MARY CLASSIFIERS
        % 		results = prtActor.run(prtActor.dataSet);
        % 		decider = prtDecisionBinaryMinPe; % can't use prtActor.internalDecider
        % 		decider = decider.train(results);
        % 		threshold = decider.threshold;
        % 		crange = threshold + min(max(results.X)-threshold,threshold-min(results.X))*[-1,1];

        xx = reshape(linGrid(:,1),gridSize);
        yy = reshape(linGrid(:,2),gridSize);
        % 		if all(isfinite(crange))
        imageHandle = imagesc(xx(1,:),yy(:,1),reshape(DS,gridSize));
        % 		else
        % 			imageHandle = imagesc(xx(1,:),yy(:,1),reshape(DS,gridSize));
        % 		end
        colormap(cMap)
    case 3
        xx = reshape(linGrid(:,1),gridSize);
        yy = reshape(linGrid(:,2),gridSize);
        zz = reshape(linGrid(:,3),gridSize);
        imageHandle = slice(xx,yy,zz,reshape(DS,gridSize),max(xx(:)),max(yy(:)),[min(zz(:)),mean(zz(:))]);
        view(3)
        colormap(cMap)
        imageHandle = imageHandle(1); % Sorry we need to throw the others away
end
axis tight;
axis xy;

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function bsf_test ( )

%*****************************************************************************80
%
%% BSF_TEST tests BSF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BSF_TEST:\n' );
  fprintf ( 1, '  A demonstration of the Black-Scholes formula\n' );
  fprintf ( 1, '  for option valuation.\n' );

  s0 = 2.0;
  t0 = 0.0;
  e = 1.0;
  r = 0.05;
  sigma = 0.25;
  t1 = 3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The asset price at time T0, S0    = %f\n', s0 );
  fprintf ( 1, '  The time                    T0    = %f\n', t0 );
  fprintf ( 1, '  The exercise price          E     = %f\n', e );
  fprintf ( 1, '  The interest rate           R     = %f\n', r );
  fprintf ( 1, '  The asset volatility        SIGMA = %f\n', sigma );
  fprintf ( 1, '  The expiry date             T1    = %f\n', t1 );

  c = bsf ( s0, t0, e, r, sigma, t1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The option value C = %f\n', c );

  return
end


text:
function y = cone(x)
%CONE Defines a low-level exponential cone constraint x(2)*exp(x(1)/x(2)) <= x(3)
%
% Input
%    x       : Linear 3x1 SDPVAR object
%
% Example
%
% Standard  exponential cone constraint x(2)*exp(x(1)/x(2)) <= x(3)
%    F = expcone(x)
%
% To quickly define several cones, the argument can be a matrix, and the
% command is then short-hand for
% for i = 1:size(x,2);F = [F,cone(x(:,i))];end
%
% See also  @SDPVAR/CONE

[n,m] = size(x);
if min([n m])==1
    x = reshape(x,3,1);
end
[n,m] = size(x);
if n ~=3
    error('x must be a vector or matrix of height 3')
end
y = x;
if min([n m])>1
	y.typeflag = 22;
else
	y.typeflag = 21;
end
y = lmi(y);

author: shamilmamedov


text:
%
% rf(x, y, z, errtol)
%
% Inputs:
%
%   x       Input vector size 1xN.
%   y       Input vector size 1xN.
%   z       Input vector size 1xN.
%   errtol  Error tolerance.
%
% Matlab function to compute Carlson's symmetric elliptic integral Rf.
% Implementation of Carlson's Duplication Algorithm 1 in "Computing
% Elliptic Integrals by Duplication," by B. C. Carlson, Numer. Math.
% 33, 1-16 (1979).
%
% Returns NaN's for any argument values outside input range.
%
% Algorithm is also from Carlson's ACM TOMS Algorithm 577.
%
% This code is a complete rewrite of the algorithm in vectorized form.
% It was not produced by running a FORTRAN to Matlab converter.
%
% The following text is copied from ACM TOMS Algorithm 577 FORTRAN code:
%
%   X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
%
%   ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
%   RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
%   16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
%
%   SAMPLE CHOICES:  ERRTOL     RELATIVE TRUNCATION
%                               ERROR LESS THAN
%                    1.D-3      3.D-19
%                    3.D-3      2.D-16
%                    1.D-2      3.D-13
%                    3.D-2      2.D-10
%                    1.D-1      3.D-7
%
% Note by TRH:
%
%   Absolute truncation error when the integrals are order 1 quantities
%   is closer to errtol, so be careful if you want high absolute precision.
%
% Thomas R. Hoffend Jr., Ph.D.
% 3M Company
% 3M Center Bldg. 236-GC-26
% St. Paul, MN 55144
% trhoffendjr@mmm.com
%

function f = rf(x, y, z, errtol)

% Argument limits as set by Carlson:
LoLim = 5.0 * realmin;
UpLim = 5.0 * realmax;

% Check input arguments for acceptability:
mask = (min([x; y; z]) >= 0) & ...
       (min([(x + y); (x + z); (y + z)]) >= LoLim) & ...
       (max([x; y; z]) < UpLim);

% Define internally acceptable variable ranges for iterations:
Xi = x(mask);
Yi = y(mask);
Zi = z(mask);

% Carlson's duplication algorithm for Rf:
Xn = Xi;
Yn = Yi;
Zn = Zi;
Mu = (Xn + Yn + Zn) / 3.0d+0;
Xndev = 2.0 - (Mu + Xn) ./ Mu;
Yndev = 2.0 - (Mu + Yn) ./ Mu;
Zndev = 2.0 - (Mu + Zn) ./ Mu;
epslon = max( abs([Xndev Yndev Zndev]) );
while (epslon >= errtol)
    Xnroot = sqrt(Xn);
    Ynroot = sqrt(Yn);
    Znroot = sqrt(Zn);
    lambda = Xnroot .* (Ynroot + Znroot) + Ynroot .* Znroot;
    Xn = 0.25 * (Xn + lambda);
    Yn = 0.25 * (Yn + lambda);
    Zn = 0.25 * (Zn + lambda);
    Mu = (Xn + Yn + Zn) / 3.0d+0;
    Xndev = 2.0 - (Mu + Xn) ./ Mu;
    Yndev = 2.0 - (Mu + Yn) ./ Mu;
    Zndev = 2.0 - (Mu + Zn) ./ Mu;
    epslon = max( abs([Xndev Yndev Zndev]) );
end
C1 = 1.0 / 24.0;
C2 = 3.0 / 44.0;
C3 = 1.0 / 14.0;
E2 = Xndev .* Yndev - Zndev .* Zndev;
E3 = Xndev .* Yndev .* Zndev;
S = 1.0 + (C1 * E2 - 0.1D0 - C2 * E3) .* E2 + C3 * E3;
f(mask) = S ./ sqrt(Mu);

% Return NaN's where input argument was out of range:
f(~mask) = NaN;


text:
function area = sphere_triangle_angles_to_area ( r, a, b, c )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_ANGLES_TO_AREA computes the area of a spherical triangle.
%
%  Discussion:
%
%    A sphere centered at 0 in 3D satisfies the equation:
%
%      X*X + Y*Y + Z*Z = R*R
%
%    A spherical triangle is specified by three points on the surface
%    of the sphere.
%
%    The area formula is known as Girard's formula.
%
%    The area of a spherical triangle is:
%
%      AREA = ( A + B + C - PI ) * R*R
%
%    where A, B and C are the (surface) angles of the triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real A, B, C, the angles of the triangle.
%
%    Output, real AREA, the area of the spherical triangle.
%

%
%  Apply Girard's formula.
%
  area = r * r * ( a + b + c - pi );

  return
end


text:
load hall1-200;
matlabpool open;
X = normalize(XO);
[P Q] = RPMF(X, 2, 1, 1, 1e-2);
show(X, P * Q, abs(X - P * Q), [144 176]);
matlabpool close;
%%
L = P * Q;
%S = abs(X - P * Q);
S = X - P * Q;
isize = [144 176];
for i=1:min([300,size(X,2)])
    subplot(1,3,1);imshow(reshape(X(:,i),isize),[]), title('X(Sample)');
    subplot(1,3,2);imshow(reshape(L(:,i),isize),[]), title('L(Low-rank)');
    subplot(1,3,3);imshow(reshape(S(:,i),isize),[]), title('S(Sparse)');
    pause(0.01);
end

author: andrewssobral

repo: lrslibrary


text:
% Compute the per-pixel perceptual energy over each window

function F = WinAvgPE(img, R)

    ff = fft2(img);
    ffc = fftshift(ff);

    [m, n]=size(img);
    [u,v]=freqspace([m,n],'meshgrid');
%     u = u/2*m; v = v/2*n;
%     u = u/m*R; v = v/n*R;
    u = u/2*R; v = v/2*R;
    r=sqrt(u.^2+v.^2);
    fw = MannosSkarision(r);

    affc = abs(ffc(:));
    energ = sqrt(sum(affc.*affc));
    F = abs(ffc.*fw);

    F = sqrt(sum(F(:).*F(:)))/(m*n);
end

%%
function fc = MannosSkarision(r)
    fc = 2.6*(0.0192+0.114*r).*exp(-(0.114*r).^1.1);
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function b = r85_vxm ( n, a, x )

%*****************************************************************************80
%
%% R85_VXM multiplies a vector by a R85 matrix.
%
%  Discussion:
%
%    The R85 storage format represents a pentadiagonal matrix as a 5
%    by N array, in which each row corresponds to a diagonal, and
%    column locations are preserved.  Thus, the original matrix is
%    "collapsed" vertically into the array.
%
%  Example:
%
%    Here is how a R85 matrix of order 6 would be stored:
%
%       *   *  A13 A24 A35 A46
%       *  A12 A23 A34 A45 A56
%      A11 A22 A33 A44 A55 A66
%      A21 A32 A43 A54 A65  *
%      A31 A42 A53 A64  *   *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the linear system.
%
%    Input, real A(5,N), the R85 matrix.
%
%    Input, real X(N), the vector to be multiplied by A'.
%
%    Output, real B(N), the product A' * x.
%
  b(1:n)   =            a(3,1:n)   * x(1:n);
  b(2:n)   = b(2:n)   + a(4,1:n-1) * x(1:n-1);
  b(3:n)   = b(3:n)   + a(5,1:n-2) * x(1:n-2);
  b(1:n-1) = b(1:n-1) + a(2,2:n)   * x(2:n);
  b(1:n-2) = b(1:n-2) + a(1,3:n)   * x(3:n);

  return
end


text:
% Get vertex coordinates for the patch type 'P'
% [C, N, S, W, E] = HealpixGetPatchVertexCoordsP(n, i, j, INFO)
%
% Parameters
% n : grid resolution
% i : ring index
% j : intra-ring index
% INFO : intermediate information (output of HealpixSelectPatchClass())
% C : intra-patch coordinates
% N : coordinates for north vertex
% S : coordinates for south vertex
% W : coordinates for west vertex
% E : coordinates for east vertex

function [C, N, S, W, E] = HealpixGetPatchVertexCoordsP(n, i, j, INFO)

% gradient of the border
grad = INFO.polar_part;

decimal_i = INFO.decimal_i_n;
int_i = INFO.int_i_n;

int_j = fix(j - grad * decimal_i);
decimal_j = j - int_j;
if int_j == 0
    int_j = 4 * int_i;
end

if decimal_j > (grad + 1) * decimal_i
    north_i = int_i - 1;
    south_i = int_i + 1;
    north_int_j = mod(int_j - grad - 1, 4 * north_i) + 1;
    south_int_j = mod(int_j + (grad + 1) - 1, 4 * south_i) + 1;
    east_int_j = mod(int_j + 1 - 1, 4 * int_i) + 1;

    N = [north_i, north_int_j];
    S = [south_i, south_int_j];
    W = [int_i  , int_j];
    E = [int_i  , east_int_j];

    offset_j = (grad + 0.5) * decimal_i;
    C = [decimal_i, decimal_j - offset_j - 0.5];
else
    curnt_i = int_i + 1;
    south_i = int_i + 2;
    if south_i > n
        south_i = n;
    end
    west_int_j = mod(int_j + grad - 1, 4 * curnt_i) + 1;
    east_int_j = mod(int_j + grad + 1 - 1, 4 * curnt_i) + 1;
    south_int_j = mod(west_int_j + (grad + 1) - 1, 4 * south_i) + 1;

    N = [int_i    , int_j];
    S = [int_i + 2, south_int_j];
    W = [int_i + 1, west_int_j];
    E = [int_i + 1, east_int_j];

    offset_j = (grad + 0.5) * (1 - decimal_i);
    C = [decimal_i - 1, decimal_j + offset_j - (grad + 0.5)];
end

if INFO.is_south_pole
    TMP = N;
    N = S;
    S = TMP;
    N(1) = 4 * n - N(1);
    S(1) = 4 * n - S(1);
    W(1) = 4 * n - W(1);
    E(1) = 4 * n - E(1);
    C(1) = -C(1);
end

author: mtex-toolbox

repo: mtex

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


text:
function [sp_fit_group,sp_profiles,sp_fit] = spectdecompose(sp_fit,options,chan)
%
% From a multitaper or MAR estimation of the states spectra (per subject),
% it factorises the spectral information (which is given by frequency bin)
% into a set of spectral components (which are much fewer than the
% frequency bins), in order to facilitate interpretation and visualisation.
% The decomposition can be done either using PCA or non-negative matrix
% factorisation or NNMF (the default). Note that if NNMF is used, the
% solution might somewhat vary every time the function is called.
%
% INPUTS:
%
% sp_fit                The output of hmmspectramar or hmmspectramt. It
%                       should be a cell, where each element corresponds
%                       to one subject
% options               Struct indicating the decomposing options
%        .Ncomp         Number of components (default 4)
%        .Method        'NNMF' or 'PCA' (default 'NNMF')
%        .Base          What to base the factorisation on: PSD if 'psd' or on
%                       coherence if 'coh' (default: 'coh')
%        .sp_profiles   Spectral profiles; if supplied, these will be used
%                       instead of computing them (default: empty)
%        .plot          If the spectral profiles are to be plotted (default: 0)
% chan                  If supplied, a vector of integers defining which
%                       channels will be used in the decomposion;
%                       if not specified, all will beused
%
% OUTPUT:
%
% sp_fit_group          A single struct containing the mean across subjects
%                       of sp_fit. If there is only one subject, then
%                       sp_fit_group equals sp_fit
% sp_profiles           The (frequency bins by spectral components) mixing
%                       matrix used to project from (no. of frequency bins
%                       by regions) to (no. of components by regions)
% sp_fit                A cell with the spectral components, with the same fields
%                       and dimensions that the input 'sp_fit', but with Ncomp
%                       components instead of No. frequency bins. If there
%                       is only one subject, then sp_fit is a single struct
%
% Authors: Mark Woolrich, OHBA, University of Oxford (2017)
%          Diego Vidaurre, OHBA, University of Oxford (2017)
%          Andrew Quinn, OHBA, University of Oxford (2018)

if ~iscell(sp_fit)
    warning('Only one subject was provided; recommended to run at the group level')
    aux = sp_fit;
    sp_fit = {};
    sp_fit{1} = aux;
    clear aux
end

if nargin < 2, options = struct(); end
if ~isfield(options,'Niterations'), options.Niterations = 10; end
if ~isfield(options,'Ncomp'), options.Ncomp = 4; end
if ~isfield(options,'Method'), options.Method = 'NNMF'; end
if ~isfield(options,'Base'), options.Base = 'coh'; end
if ~isfield(options,'plot'), options.plot = 0; end

ndim = size(sp_fit{1}.state(1).psd,2); % no. channels
if nargin < 3, chan = 1:ndim;
else, ndim = length(chan);
end
ndim2 = ndim*(ndim-1)/2;
Nf = size(sp_fit{1}.state(1).psd,1); % no. frequencies
N = length(sp_fit); % no. subjects
K = length(sp_fit{1}.state); % no. states
ind_offdiag = triu(true(ndim),1)==1;

% check that all estimations have the same number of freq bins
for n = 1:N
    if isempty(sp_fit{n})
        error('One or more elements of sp_fit are empty - remove them first')
    end
    if size(sp_fit{n}.state(1).psd,1) ~= Nf
        error(['It is necessary for the spectral estimation of all subjects ' ...
            'to have the same number of frequency bins. ' ...
            'In the case of the multitaper, this is can be done by setting ' ...
            'the options tapers and win to a fixed value.'])
    end
end

% put coh and psd in temporary arrays
psd_comps = zeros(N,K,Nf,ndim);
coh_comps = zeros(N,K,Nf,ndim,ndim);
for n = 1:N
    for k = 1:K
        psd = sp_fit{n}.state(k).psd(:,chan,chan);
        coh = sp_fit{n}.state(k).coh(:,chan,chan);
        for j = chan
            psd_comps(n,k,:,j) = psd(:,j,j);
            for l = chan
                coh_comps(n,k,:,j,l) = coh(:,j,l);
            end
        end
    end
end

% Build the matrix that is to be factorised
Xpsd = zeros(Nf,ndim*K);
keep_psd = true(1,ndim*K);
for k = 1:K
    ind = (1:ndim) + (k-1)*ndim;
    notnan = ~isnan(psd_comps(:,k,1,1)) & ~isinf(psd_comps(:,k,1,1));
    if any(notnan)
        Xpsd(:,ind)= squeeze(mean(abs(psd_comps(notnan,k,:,:)),1)); % mean across subjects
    else
        keep_psd(ind) = false;
    end
end
Xcoh = zeros(Nf,K*ndim2);
keep_coh = true(1,ndim2*K);
for k = 1:K
    ind = (1:ndim2) + (k-1)*ndim2;
    notnan = ~isnan(coh_comps(:,k,1,1,2)) & ~isinf(coh_comps(:,k,1,1,2));
    if any(notnan)
        ck = squeeze(mean(abs(coh_comps(notnan,k,:,:,:)),1));
        Xcoh(:,ind) = ck(:,ind_offdiag);
    else
        Xcoh(:,ind) = false;
    end

end
if strcmpi(options.Base,'psd')
    X = Xpsd(:,keep_psd);
else
    X = Xcoh(:,keep_coh);
end

% Specify fit function, a unimodal gaussian
gauss_func = @(x,f) f.a1.*exp(-((x-f.b1)/f.c1).^2);
% Default fit options
options_fit = fitoptions('gauss1');
% constrain lower and upper bounds
options_fit.Lower = [0,1,0];
options_fit.Upper = [Inf,size(X,1),size(X,1)];

if ~isfield(options,'sp_profiles') || isempty(options.sp_profiles)
    % Doing the decomposition
    if strcmpi(options.Method,'NNMF')
        bestval = Inf; fitval = zeros(options.Niterations,1);
        for ii = 1:options.Niterations
            try
                [A,~] = nnmf(X,options.Ncomp,'replicates',500,'algorithm','als');
            catch
                error('nnmf not found - perhaps the Matlab version is too old')
            end
            for jj = 1:size(A,2)
                f = fit( linspace(1,size(A,1),size(A,1))',A(:,jj), 'gauss1',options_fit);
                residuals = A(:,jj) - gauss_func(1:size(A,1),f)';
                fitval(ii) = fitval(ii) + sum( residuals.^2 ) / size(A,2);
            end
            if fitval(ii) < bestval
                bestval = fitval(ii);
                bestfit = A;
            end
        end
        %sp_profiles = pinv(X') * b'; % you lose non-negativity by doing this
        sp_profiles = bestfit;
        [~,ind] = max(sp_profiles,[],1);
        [~,neworder_auto] = sort(ind);
        sp_profiles = sp_profiles(:,neworder_auto);
    else
        try
            m = mean(X);
            X = X - repmat(m,Nf,1);
            [~,A] = pca(X,'NumComponents',options.Ncomp);
        catch
            error('Error running pca - maybe not matlab''s own?')
        end
        sp_profiles = A;
    end
else
    sp_profiles = options.sp_profiles;
end

% plot if required
if options.plot
    figure;
    if options.Ncomp > 4
        j1 = ceil(options.Ncomp/2); j2 = 2;
    else
        j1 = options.Ncomp; j2 = 1;
    end
    for j = 1:options.Ncomp
        subplot(j1,j2,j)
        plot(sp_profiles(:,j),'LineWidth',2.5)
    end
end

% group level
sp_fit_group = struct();
sp_fit_group.state = struct();

% NNMF / PCA
psd = zeros(ndim*K,options.Ncomp);
coh = zeros(ndim2*K,options.Ncomp);
if strcmpi(options.Method,'NNMF')
   psd(keep_psd,:) = (Xpsd(:,keep_psd)' * sp_profiles); % ndim by components
   coh(keep_coh,:) = (Xcoh(:,keep_coh)' * sp_profiles);
%     opt = statset('maxiter',0);
%     [~,b] = nnmf(Xpsd,options.Ncomp,'algorithm','als',...
%         'w0',sp_profiles,'Options',opt);
%     psd = b'; % regions by components
%     [~,b] = nnmf(Xcoh,options.Ncomp,'algorithm','als',...
%         'w0',sp_profiles,'Options',opt);
%     coh = b'; % pairs of regions by components
else
    Xpsd(:,keep_psd) = Xpsd(:,keep_psd) - repmat(mean(Xpsd(:,keep_psd)),Nf,1);
    Xcoh(:,keep_coh) = Xcoh(:,keep_coh) - repmat(mean(Xcoh(:,keep_coh)),Nf,1);
    psd(keep_psd,:) = (Xpsd(:,keep_psd)' * sp_profiles); % ndim by components
    coh(keep_coh,:) = (Xcoh(:,keep_coh)' * sp_profiles);
end
if any(isnan(psd(:))) || any(isnan(coh(:)))
   warning('There are NaNs in the estimations for those states that were not used')
end
for k = 1:K
    sp_fit_group.state(k).psd = zeros(options.Ncomp,ndim,ndim);
    sp_fit_group.state(k).coh = ones(options.Ncomp,ndim,ndim);
    ind = (1:ndim) + (k-1)*ndim;
    for i = 1:options.Ncomp
        sp_fit_group.state(k).psd(i,:,:) = diag(psd(ind,i));
    end
    ind = (1:ndim2) + (k-1)*ndim2;
    for i = 1:options.Ncomp
        graphmat = zeros(ndim);
        graphmat(ind_offdiag) = coh(ind,i);
        graphmat=(graphmat+graphmat') + eye(ndim);
        sp_fit_group.state(k).coh(i,:,:) = graphmat;
    end
end

% Subject level
if N > 1 && nargout == 3
    for n = 1:N
        sp_fit{n} = struct();
        sp_fit{n}.state = struct();
        % prepare matrix
        Xpsd = zeros(Nf,ndim*K);
        keep_psd = true(1,ndim*K);
        for k = 1:K
            ind = (1:ndim) + (k-1)*ndim;
            Xpsd(:,ind)= squeeze(abs(psd_comps(n,k,:,:)));
            if any(isnan(var(Xpsd(:,ind)))) || any(isinf(var(Xpsd(:,ind))))
                keep_psd(ind) = false;
                warning(['Session ' num2str(n) ' did not use state ' num2str(k) '; PSD set to NaN'])
            end
        end
        Xcoh = zeros(Nf,K*ndim2);
        keep_coh = true(1,K*ndim2);
        for k = 1:K
            ind = (1:ndim2) + (k-1)*ndim2;
            ck = squeeze(abs(coh_comps(n,k,:,:,:)));
            Xcoh(:,ind) = ck(:,ind_offdiag);
            if any(isnan(var(Xcoh(:,ind)))) || any(isinf(var(Xcoh(:,ind))))
                keep_coh(ind) = false;
                warning(['Session ' num2str(n) ' did not use state ' num2str(k) '; Coh set to NaN'])
            end
        end
        % NNMF / PCA
        psd = zeros(ndim*K,options.Ncomp);
        coh = zeros(ndim2*K,options.Ncomp);
        if strcmpi(options.Method,'NNMF')
            opt = statset('maxiter',1);
            [~,b] = nnmf(Xpsd(:,keep_psd),options.Ncomp,'algorithm','als',...
                'w0',sp_profiles,'Options',opt);
            psd(keep_psd,:) = b'; % regions by components
            [~,b] = nnmf(Xcoh(:,keep_coh),options.Ncomp,'algorithm','als',...
                'w0',sp_profiles,'Options',opt);
            coh(keep_coh,:) = b';
        else
            Xpsd(:,keep_psd) = Xpsd(:,keep_psd) - repmat(mean(Xpsd(:,keep_psd)),Nf,1);
            Xcoh(:,keep_coh) = Xcoh(:,keep_coh) - repmat(mean(Xcoh(:,keep_coh)),Nf,1);
            psd(keep_psd,:) = (Xpsd(:,keep_psd)' * sp_profiles);
            coh(keep_coh,:) = (Xcoh(:,keep_coh)' * sp_profiles);
        end
        % Reshape stuff
        for k = 1:K
            sp_fit{n}.state(k).psd = zeros(options.Ncomp,ndim,ndim);
            sp_fit{n}.state(k).coh = ones(options.Ncomp,ndim,ndim);
            ind = (1:ndim) + (k-1)*ndim;
            for i = 1:options.Ncomp
                sp_fit{n}.state(k).psd(i,:,:) = diag(psd(ind,i));
            end
            ind = (1:ndim2) + (k-1)*ndim2;
            for i = 1:options.Ncomp
                graphmat = zeros(ndim);
                graphmat(ind_offdiag) = coh(ind,i);
                graphmat=(graphmat+graphmat') + eye(ndim);
                sp_fit{n}.state(k).coh(i,:,:) = graphmat;
            end
        end
    end
else
    sp_fit = sp_fit_group;
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function varargout = piano2(varargin)
% PIANO2 M-file for piano2.fig
%      PIANO2, by itself, creates a new PIANO2 or raises the existing
%      singleton*.
%
%      H = PIANO2 returns the handle to a new PIANO2 or the handle to
%      the existing singleton*.
%
%      PIANO2('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in PIANO2.M with the given input arguments.
%
%      PIANO2('Property','Value',...) creates a new PIANO2 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before piano2_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to piano2_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help piano2

% Last Modified by GUIDE v2.5 29-Aug-2008 11:01:33

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @piano2_OpeningFcn, ...
                   'gui_OutputFcn',  @piano2_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% AUTHOR: MARK FRANKOSKY (MFRANKOSKY@GMAIL.COM)

% --- Executes just before piano2 is made visible.
function piano2_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to piano2 (see VARARGIN)

% Choose default command line output for piano2
handles.output      = hObject;
handles.SampleRate  = 1/20000;
handles.SoundVector = 0;
handles.TimeValue   = 0.3488;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes piano2 wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = piano2_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in GShparp3.
function GShparp3_Callback(hObject, eventdata, handles)
% hObject    handle to GShparp3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 207.65;
PlayNote(handles,freq);

% --- Executes on button press in ASharp3.
function ASharp3_Callback(hObject, eventdata, handles)
% hObject    handle to ASharp3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 233.08;
PlayNote(handles,freq);

% --- Executes on button press in FSharp5.
function FSharp5_Callback(hObject, eventdata, handles)
% hObject    handle to FSharp5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 739.99;
PlayNote(handles,freq);

% --- Executes on button press in FSharp3.
function FSharp3_Callback(hObject, eventdata, handles)
% hObject    handle to FSharp3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 185.00;
PlayNote(handles,freq);

% --- Executes on button press in ASharp4.
function ASharp4_Callback(hObject, eventdata, handles)
% hObject    handle to ASharp4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 466.16;
PlayNote(handles,freq);

% --- Executes on button press in DSharp5.
function DSharp5_Callback(hObject, eventdata, handles)
% hObject    handle to DSharp5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 622.25;
PlayNote(handles,freq);

% --- Executes on button press in CSharp5.
function CSharp5_Callback(hObject, eventdata, handles)
% hObject    handle to CSharp5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 554.37;
PlayNote(handles,freq);

% --- Executes on button press in GSharp4.
function GSharp4_Callback(hObject, eventdata, handles)
% hObject    handle to GSharp4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 415.3;
PlayNote(handles,freq);

% --- Executes on button press in CSharp4.
function CSharp4_Callback(hObject, eventdata, handles)
% hObject    handle to CSharp4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 277.18;
PlayNote(handles,freq);

% --- Executes on button press in G3.
function G3_Callback(hObject, eventdata, handles)
% hObject    handle to G3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 196;
PlayNote(handles,freq);

% --- Executes on button press in A3.
function A3_Callback(hObject, eventdata, handles)
% hObject    handle to A3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 220;
PlayNote(handles,freq);

% --- Executes on button press in B3.
function B3_Callback(hObject, eventdata, handles)
% hObject    handle to B3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 246.94;
PlayNote(handles,freq);

% --- Executes on button press in C4.
function C4_Callback(hObject, eventdata, handles)
% hObject    handle to C4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 261.625;
PlayNote(handles,freq);

% --- Executes on button press in D4.
function D4_Callback(hObject, eventdata, handles)
% hObject    handle to D4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 293.66;
PlayNote(handles,freq);

% --- Executes on button press in E4.
function E4_Callback(hObject, eventdata, handles)
% hObject    handle to E4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 329.63;
PlayNote(handles,freq);

% --- Executes on button press in F4.
function F4_Callback(hObject, eventdata, handles)
% hObject    handle to F4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 349.23;
PlayNote(handles,freq);

% --- Executes on button press in G4.
function G4_Callback(hObject, eventdata, handles)
% hObject    handle to G4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 392;
PlayNote(handles,freq);

% --- Executes on button press in A4.
function A4_Callback(hObject, eventdata, handles)
% hObject    handle to A4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 440;
PlayNote(handles,freq);

% --- Executes on button press in B4.
function B4_Callback(hObject, eventdata, handles)
% hObject    handle to B4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 493.88;
PlayNote(handles,freq);

% --- Executes on button press in C5.
function C5_Callback(hObject, eventdata, handles)
% hObject    handle to C5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 523.25;
PlayNote(handles,freq);

% --- Executes on button press in D5.
function D5_Callback(hObject, eventdata, handles)
% hObject    handle to D5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 587.33;
PlayNote(handles,freq);

% --- Executes on button press in E5.
function E5_Callback(hObject, eventdata, handles)
% hObject    handle to E5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 659.26;
PlayNote(handles,freq);

% --- Executes on button press in F5.
function F5_Callback(hObject, eventdata, handles)
% hObject    handle to F5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 698.46;
PlayNote(handles,freq);

% --- Executes on button press in FSharp4.
function FSharp4_Callback(hObject, eventdata, handles)
% hObject    handle to FSharp4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 369.99;
PlayNote(handles,freq);

% --- Executes on button press in DSharp4.
function DSharp4_Callback(hObject, eventdata, handles)
% hObject    handle to DSharp4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 311.13;
PlayNote(handles,freq);

% --- Executes on button press in G5.
function G5_Callback(hObject, eventdata, handles)
% hObject    handle to G5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 783.99;
PlayNote(handles,freq);

% --- Executes on button press in F3.
function F3_Callback(hObject, eventdata, handles)
% hObject    handle to F3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 174.61;
PlayNote(handles,freq);

% --- Executes on button press in E3.
function E3_Callback(hObject, eventdata, handles)
% hObject    handle to E3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 164.81;
PlayNote(handles,freq);

% --- Executes on button press in D3.
function D3_Callback(hObject, eventdata, handles)
% hObject    handle to D3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 146.83;
PlayNote(handles,freq);

% --- Executes on button press in A5.
function A5_Callback(hObject, eventdata, handles)
% hObject    handle to A5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 880;
PlayNote(handles,freq);

% --- Executes on button press in C3.
function C3_Callback(hObject, eventdata, handles)
% hObject    handle to C3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 130.81;
PlayNote(handles,freq);

% --- Executes on button press in B5.
function B5_Callback(hObject, eventdata, handles)
% hObject    handle to B5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 987.77;
PlayNote(handles,freq);

% --- Executes on button press in GSharp5.
function GSharp5_Callback(hObject, eventdata, handles)
% hObject    handle to GSharp5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 830.61;
PlayNote(handles,freq);

% --- Executes on button press in ASharp5.
function ASharp5_Callback(hObject, eventdata, handles)
% hObject    handle to ASharp5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 932.33;
PlayNote(handles,freq);

% --- Executes on button press in DSharp3.
function DSharp3_Callback(hObject, eventdata, handles)
% hObject    handle to DSharp3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 155.56;
PlayNote(handles,freq);

% --- Executes on button press in CSharp3.
function CSharp3_Callback(hObject, eventdata, handles)
% hObject    handle to CSharp3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

freq = 138.59;
PlayNote(handles,freq);

function PlayNote(handles,freq)
% plays the note that was pressed for 1 second
% also sets the frequency value on GUI for last played note.

strFreq = num2str(freq);
set(handles.STFreqValue, 'String', strFreq);

SampleRate  = handles.SampleRate;
TimeValue   = handles.TimeValue;
Samples     = 0:SampleRate:TimeValue;
SinOn       = get(handles.SinWave, 'Value');
SquareOn    = get(handles.SquareWave, 'Value');
SawtoothOn  = get(handles.SawtoothWave, 'Value');

if SinOn == 1
  soundVector = sin(2*pi*freq*Samples);
elseif SquareOn == 1
  soundVector = square(2*pi*freq*Samples);
elseif SawtoothOn == 1
  soundVector = sawtooth(2*pi*freq*Samples);
end

sound(soundVector, 1/SampleRate)

RecordOn = get(handles.Record, 'Value');

if RecordOn == 1
  SoundVectorLong = handles.SoundVector;
  if SoundVectorLong == 0
    SoundVectorLong = soundVector;
  else
    SoundVectorLong = cat(2, SoundVectorLong, soundVector);
  end
  handles.SoundVector = SoundVectorLong;
  guidata(handles.figure1, handles);
end

% --- Executes on button press in Play.
function Play_Callback(hObject, eventdata, handles)
% hObject    handle to Play (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

SoundVector   = handles.SoundVector;
SampleRate    = handles.SampleRate;

sound(SoundVector, 1/SampleRate)

% --- Executes on button press in Stop.
function Stop_Callback(hObject, eventdata, handles)
% hObject    handle to Stop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.SoundVector = 0;
guidata(handles.figure1, handles);

function Time_Callback(hObject, eventdata, handles)
% hObject    handle to Time (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Time as text
%        str2double(get(hObject,'String')) returns contents of Time as a double

Value     = str2double(get(hObject, 'String'));
TimeValue = 60/Value;

half    = get(handles.HalfNote, 'Value');
Whole   = get(handles.WholeNote, 'Value');
Eighth  = get(handles.EighthNote, 'Value');

if half == 1
    TimeValue = TimeValue * 2;
elseif Whole == 1
    TimeValue = TimeValue * 4;
elseif Eighth == 1
    TimeValue = TimeValue / 2;
end

handles.TimeValue = TimeValue;
guidata(handles.figure1, handles);

% --- Executes on button press in QuarterNote.
function QuarterNote_Callback(hObject, eventdata, handles)
% hObject    handle to QuarterNote (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.QuarterNote, 'Value', 1);
set(handles.HalfNote, 'Value', 0);
set(handles.EighthNote, 'Value', 0);
set(handles.WholeNote, 'Value', 0);

Value     = str2double(get(handles.Time, 'String'));
TimeValue = 60/Value;

handles.TimeValue = TimeValue;
guidata(handles.figure1, handles);

% --- Executes on button press in HalfNote.
function HalfNote_Callback(hObject, eventdata, handles)
% hObject    handle to HalfNote (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.QuarterNote, 'Value', 0);
set(handles.HalfNote, 'Value', 1);
set(handles.EighthNote, 'Value', 0);
set(handles.WholeNote, 'Value', 0);

Value     = str2double(get(handles.Time, 'String'));
TimeValue = 60/Value * 2;

handles.TimeValue = TimeValue;
guidata(handles.figure1, handles);

% --- Executes on button press in EighthNote.
function EighthNote_Callback(hObject, eventdata, handles)
% hObject    handle to EighthNote (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.QuarterNote, 'Value', 0);
set(handles.HalfNote, 'Value', 0);
set(handles.EighthNote, 'Value', 1);
set(handles.WholeNote, 'Value', 0);

Value     = str2double(get(handles.Time, 'String'));
TimeValue = (60/Value) / 2;

handles.TimeValue = TimeValue;
guidata(handles.figure1, handles);

% --- Executes on button press in WholeNote.
function WholeNote_Callback(hObject, eventdata, handles)
% hObject    handle to WholeNote (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(handles.QuarterNote, 'Value', 0);
set(handles.HalfNote, 'Value', 0);
set(handles.EighthNote, 'Value', 0);
set(handles.WholeNote, 'Value', 1);

Value     = str2double(get(handles.Time, 'String'));
TimeValue = 60/Value * 4;

handles.TimeValue = TimeValue;
guidata(handles.figure1, handles);

% --- Executes during object creation, after setting all properties.
function QuarterNote_CreateFcn(hObject, eventdata, handles)
% hObject    handle to QuarterNote (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

set(hObject, 'Value', 1);

% --- Executes on button press in OutToWav.
function OutToWav_Callback(hObject, eventdata, handles)
% hObject    handle to OutToWav (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

SoundVector     = handles.SoundVector;
SampleRate      = handles.SampleRate;

if SoundVector == 0
    return;
else
  Filename = get(handles.Filename, 'String');
  wavwrite(SoundVector, 1/SampleRate, 32, Filename);
end

% --- Executes on button press in SinWave.
function SinWave_Callback(hObject, eventdata, handles)
% hObject    handle to SinWave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of SinWave

set(handles.SinWave, 'Value', 1);
set(handles.SquareWave, 'Value', 0);
set(handles.SawtoothWave, 'Value', 0);

% --- Executes during object creation, after setting all properties.
function SinWave_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SinWave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

set(hObject, 'Value', 1);

% --- Executes on button press in SquareWave.
function SquareWave_Callback(hObject, eventdata, handles)
% hObject    handle to SquareWave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of SquareWave

set(handles.SinWave, 'Value', 0);
set(handles.SquareWave, 'Value', 1);
set(handles.SawtoothWave, 'Value', 0);

% --- Executes on button press in SawtoothWave.
function SawtoothWave_Callback(hObject, eventdata, handles)
% hObject    handle to SawtoothWave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of SawtoothWave

set(handles.SinWave, 'Value', 0);
set(handles.SquareWave, 'Value', 0);
set(handles.SawtoothWave, 'Value', 1);


text:
function dataout = ft_annotate(cfg, datain)

% FT_ANNOTATE returns the same output data as the user has provided as input, but allows
% to add comments to that data structure. These comments are stored along with the other
% provenance information and can be displayed with FT_ANALYSISPIPELINE. Adding comments
% is especially useful if you have manually (i.e. in plain MATLAB) modified the data
% structure, whereby some provenance information is missing.
%
% Use as
%   outdata = ft_examplefunction(cfg, indata)
% where the input data structure can be any of the FieldTrip data structures and
% the configuration structure should contain
%   cfg.comment    = string
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_ANALYSISPIPELINE, FT_MATH

% Copyright (C) 2013, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar datain
ft_preamble provenance datain
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% ensure that the required options are present
cfg = ft_checkconfig(cfg, 'required', 'comment');

% simply copy the input data to the output data
dataout = datain;

% this line is meant to provide some information
% but also to ensure that trackconfig does not remove the cfg.comment field
fprintf('adding the comment: %s\n', cfg.comment);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous datain
ft_postamble provenance dataout
ft_postamble history dataout
ft_postamble savevar dataout

author: spm


text:
function warpedScanNumV = getAssocWarpedScanNum(scanNumV,planC)
% Return warped scan created from input (moving) scanNum.
% ------------------------------------------------------------------------
% scanNumV   : Vector of scan nos.
% planC
% ------------------------------------------------------------------------
% AI 12/09/21

indexS = planC{end};

assocMovScanUIDc = {planC{indexS.scan}.assocMovingScanUID};

warpedScanNumV = nan(1,length(scanNumV)); % Initialize
for nScan = 1:length(scanNumV)
    movScanUID = {planC{indexS.scan}(scanNumV(nScan)).scanUID};
    assocIdxV = strcmpi(movScanUID,assocMovScanUIDc);
    if any(assocIdxV)
        warpedScanNumV(nScan) = find(assocIdxV);
    end
end
warpedScanNumV = warpedScanNumV(~isnan(warpedScanNumV));

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% Copyright (C) 2016-2017 Lagu
%% Copyright (C) 2017, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym [@var{K}, @var{E}] = ellipke (@var{m})
%% Complete elliptic integrals of the first and second kinds.
%%
%% Example:
%% @example
%% @group
%% syms m
%% [K, E] = ellipke (m)
%%   @result{} K = (sym) K(m)
%%   @result{} E = (sym) E(m)
%% @end group
%% @end example
%%
%% @seealso{ellipke, @@sym/ellipticK, @@sym/ellipticE}
%% @end defmethod

function varargout = ellipke(m)

  if (nargin ~= 1 || nargout > 2)
    print_usage ();
  end

  if (nargout == 0 || nargout == 1)
    varargout = {ellipticK(m)};
  else
    varargout = {ellipticK(m) ellipticE(m)};
  end

end

%!error ellipke (sym(1), 2)

%!test
%! for i = 2:10
%!   [K E] = ellipke (sym (1)/i);
%!   [k e] = ellipke (1/i);
%!   assert (double ([K E]), [k e], 2*eps)
%! end


text:
function element_num = grid_ql_element_num ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_QL_ELEMENT_NUM counts the elements in a grid of QL quadrilaterals.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NUM = NELEMX * NELEMY = 6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NUM, the number of elements in the grid.
%
  element_num = nelemx * nelemy;

  return
end