text:
function outStruct = mergeStructures(s1, s2)
% outStruct = mergeStructures(firstStruct, secondStruct)
%
% returns a structure with fields from both input structures.
% If both input structures have the same fields, the s2
% fields overwrite the s1 fields.
%
% ras, 09/06/2005, imported into mrVista 2.0
% ras, 07/21/2007, deals with struct arrays and empty structs.

% if s1 is empty, our job is easy: we can overwrite it with s2:
if isempty(s1)
	outStruct = s2;
	return
end

% deal with different sizes of structures:
% we can do this if s2 is a struct array, but not if
% s1 is (for now)
if length(s1) > 1
	error('firstStruct needs to be length 1 or empty.')
end

if length(s2) > 1
	% we allow only for up to 2 dimensions here
	for i = 1:size(s2, 1)
		for j = 1:size(s2, 2)
			outStruct(i,j) = mergeStructures(s1, s2(i,j));
		end
	end
	return
end

%% core part of merge:
% make sure all fields which belong to either s1 or s2 are assigned to
% outStruct:
outStruct = s1;

if ~isempty(s2)
	fieldNames = fieldnames(s2);
	for i = 1:length(fieldNames)
		f = fieldNames{i};
		outStruct.(f) = s2.(f);
	end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [hpb] = TW2hpb(TW)
% Convert power from terawatts to boiler horsepower.
% Chad A. Greene 2012
hpb = TW*101941995.00484;

author: Sable

repo: mcbench-benchmarks


text:
function [view, mapvol, covol] = polarAngleMap(view, dt, scans, params, legend, W);
%
% [view, map, co] = polarAngleMap(view, <dt, scans, params>, <legend>, <W>);
%
% AUTHOR: rory
% PURPOSE:
% Given corAnal data for a polar angle ("meridian")-mapping experiment
% (single scan or set of scans), produce a parameter map of preferred
% polar angle in units of degrees of visual field.
%
% If a single scan is provided as input, this function saves the
% parameter map in the scan's data type, assigning it only for that
% scan.
%
% However, if multiple input scans are provided (see below),
% the code saves the results in a new data type 'Meta_Analysis'.
% corAnal data from each scan are first converted into real-world
% units, then overlapping data are averaged together in a weighted
% average, based on each scan's coherence. (I.e., if one scan's
% co values for a given voxel are much higher than other scans,
% it will dominate the determination of what angle is represented.)
%
% I wrote this code to use in conjunction with my across-session
% tools (createCombinedSession, importTSeries) to run meta-analyses
% on retinotopy data.
%
% ARGUMENTS:
%   INPUT:
%   view: mrVista view. <Defaults to selected gray view>
%   dt: for a single scan, name or number of the data type
%       from which the input data come. If analyzing multiple
%       scans, a cell of length nScans of data type names/numbers.
%       <default: cur data type>
%
%   scans: scan or scans to use as input. <default: cur scan>
%
%   params: struct (or nScans-long struct array) specifying how
%       the stimulus mapped polar angle during each scan. Needs
%       the following fields:
%       params.startAngle: angle of center of wedge stimulus, measured
%           in degrees clockwise from 12-o-clock, at the start of each
%           cycle;
%       params.width: width of wedge stimulus in degrees.
%       params.direction: 'cw' or 'ccw', direction in which the stimulus
%       proceeded. (cw=clockwise or ccw=counterclockwise).
%       params.visualField: number of degrees the stimulus traversed
%       each cycle (e.g., 360 if it went all the way around).
%       <default: get these params using retinoCheckParams>
%
%   legend: optional flag which, if set to 1, provides for a separate
%       figure with a legend image to go with the polar angle map.
%       <default 0, don't show this>
%
%   W: optional vector of weights for each input scan, for use when
%       doing a meta-analysis across scans. The vector should be
%       the same length as the input scans, and should specify the
%       overall weight, on top of the coherence, that that scan's voxels
%       get. Useful for me, when high-res data produces lower co values,
%       but is actually more reliable at identifying meridian
%       representations. <default: all ones.>
%
%   OUTPUT:
%   view: mrVista view, set to the relevant data type / scan and with
%         the map loaded and set to map mode.
%
%   map: the map volume produced (but not the cell-of-scans set in the
%        view, the numeric matrix).
%
%   co: the maximum coherence at each voxel, across all the input scans.
%       (It may be more sensible to make this the mean, but I'm trying
%       max for now.) Same format as map.
%
%
% ras, 01/10/06.
if notDefined('view'),  view = getSelectedGray;                 end
if notDefined('dt'),    dt = viewGet(view, 'curDataType');      end
if notDefined('scans'), scans = viewGet(view, 'curScan');       end
if notDefined('legend'), legend = 0;                            end
if notDefined('W'),      W = ones(size(scans));                 end
if notDefined('params')
    params = retinoCheckParams(view, dt, scans);
end

mapName = 'Polar Angle (clock position)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Deal with single input scan instances separately %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if length(scans)==1
    % this should be easy: just convert the corAnal values into
    % degrees of polar angle
    view = selectDataType(view, dt);
    corAnalPath = fullfile(dataDir(view), 'corAnal.mat');
    if ~exist(corAnalPath, 'file')
        error('corAnal not found. Run computeCorAnal first.');
    end
    load(corAnalPath, 'ph', 'co');
    srcPh = ph{scans}; srcCo = co{scans}; clear ph co;

    % map from corAnal ph to polar angle
    mapvol = polarAngle(srcPh, params) ./ 30;

    % make and set the parameter map
    mapPath = fullfile(dataDir(view), 'Polar_Angle_Map.mat');
    if exist(mapPath, 'file')
        load(mapPath, 'map', 'co');
    else
        map = cell(1, numScans(view));
    end

    map{scans} = mapvol;

    % let's set the map colormap to be the same as the phase mode
    % colormap, and save this with the map
    if checkfields(view, 'ui', 'phMode')
        view.ui.mapMode.cmap = view.ui.phMode.cmap;
        view.ui.mapMode.clipMode = [0 12];
    end

    view = setParameterMap(view, map, mapName);
    saveParameterMap(view, mapPath, 1, 1);

    if legend, polarAngleMapLegend(view); end

    % that should be it!
    if nargout>=3, covol = co{scans}; end
    view = refreshScreen(view);
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if we get here, we have multiple scans: parse the arguments to    %
% be cell arrays, and get set up:                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nScans = length(scans);

if ~iscell(dt)
    mrGlobals;
    % assume a single dt is specified, and all the scans are
    % coming from that data type
    for i = 1:nScans
        if isnumeric(dt), tmp{i} = dataTYPES(dt).name;
        else,             tmp{i} = dt;
        end
    end
    dt = tmp; clear tmp;
else
    for i = 1:nScans
        if isnumeric(dt{i}), dt{i} = dataTYPES(dt{i}).name; end
    end
end

%%%%%Get corAnal volumes for each input scan
srcCo = cell(1, nScans); srcPh = cell(1, nScans);
uniqueDts = unique(dt);
for i = 1:length(uniqueDts)
    corAnalPath = fullfile(viewDir(view), uniqueDts{i}, 'corAnal.mat');
    if ~exist(corAnalPath, 'file')
        error('corAnal not found. Run computeCorAnal first.');
    end
    load(corAnalPath, 'co', 'ph', 'amp')

    I = cellfind(dt, uniqueDts{i});
    srcCo(I) = co(scans(I));
    srcPh(I) = ph(scans(I));
    srcAmp(I) = amp(scans(I));
end

%%%%%Set up the target data type for the multi-scan meta-analysis
view = initScan(view, 'Meta_Analysis', [], {dt{1} scans(1)});
view = selectDataType(view, 'Meta_Analysis');
view = setCurScan(view, numScans(view));
view = setAnnotation(view, sprintf('Meta Analysis for %s scans %s', ...
                                dt{1}, num2str(scans)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Polar Angle Map  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%Set NaNs to zero -- will be ignored
for i = 1:nScans
    srcCo{i}(isnan(srcCo{i})) = 0;
    srcAmp{i}(isnan(srcAmp{i})) = 0;
    srcPh{i}(isnan(srcPh{i})) = 0;
end

%%%%%Convert each phase map into real-word units
for i = 1:nScans
    srcPh{i} = polarAngle(srcPh{i}, params(i)) ./ 30;
end

%%%%%because we allow an additional level of user-defined weights
%%%%%(orig. b/c I wanted to weigh high-res scans higher than low-res),
%%%%%adjust the coherence weights accordingly for each scan
for i = 1:nScans, srcCo{i} = srcCo{i} .* W(i); end

%%%%%Set up the weighted average
% we'll need a volume representing the sum of the coherence
% for each voxel, across input scans. This will serve as the
% denominator of the weight formula for each input scan.
coSum = zeros(size(srcCo{1})); coMax = zeros(size(srcCo{1}));
for i = 1:nScans,
    coSum = coSum + srcCo{i};
    coMax = max(coMax, srcCo{i});
end

%%%%%initialize the map and co volumes
mapvol = zeros(size(srcCo{1}));
covol = zeros(size(srcCo{1}));

%%%%%compute co volume as the mean across all input co volumes
for i = 1:nScans, covol = covol + srcCo{i}; end
covol = covol ./ length(srcCo);

%%%%%Compute the weighted average, iteratively across input scans
for i = 1:nScans
    mapvol = mapvol + (srcPh{i} .* srcCo{i} ./ coSum);

%     % alternate attempt: use winner-take-all: scan with the
%     % highest co value for a given voxel determines the map
%     % value at that voxel.
%     Imax = find(srcCo{i}==coMax);
%     mapvol(Imax) = srcPh{i}(Imax);
end

covol = coMax;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output the parameter map %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize the output map and co data fields, loading it if it
% already exists:
mapPath = fullfile(dataDir(view), 'Polar_Angle_Map.mat');
if exist(mapPath, 'file')
    load(mapPath, 'map', 'mapName', 'co');
else
    map = cell(1, numScans(view));
    co = cell(1, numScans(view));
end

% append the map volume for the new scans
map{numScans(view)} = mapvol;
co{numScans(view)} = covol;

% before saving the map, copy the view's phase mode color map
% to the map mode, which will also be saved:
if checkfields(view, 'ui', 'phMode')
    view.ui.mapMode.cmap = view.ui.phMode.cmap;
    view.ui.mapMode.clipMode = [0 12];
end

% set the map in the view, and save it
view = setParameterMap(view, map, mapName);
saveParameterMap(view, mapPath, 1, 1);
save(mapPath, 'co', '-append'); % also add the co field
view.co = co;
refreshScreen(view);

% also set corAnal amp and ph fields, and save a corAnal, so
% we can view in that mode as well (nice colorbar, can ph-restrict)
% we map the ph back from degrees to radians:
phvol = deg2rad(mapvol.*30);
ampvol = srcAmp{1};
if exist(fullfile(dataDir(view), 'corAnal.mat'), 'file')
    view = loadCorAnal(view);
end
view.co{numScans(view)} = covol;
view.ph{numScans(view)} = phvol;
view.amp{numScans(view)} = ampvol;
view = saveCorAnal(view, 1);

newParams.type = 'polar_angle';     % set retino params such
newParams.startAngle = 0;           % that the default HSV color map
newParams.direction = 'clockwise';  % produces a nice wedge color bar
newParams.visualField = 360;
newParams.width = 0;
retinoSetParams(view, 'Meta_Analysis', numScans(view), newParams);

% show a legend if requested
if legend, polarAngleMapLegend(view); end

% ok, think that's it!

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function A = polarAngleMapLegend(view);
% img = polarAngleMapLegend(view);
% Using the current map mode settings, produce a legend
% for a polar angle parameter map and plot in a separate figure.
% Returns a truecolor image if requested.
mode = view.ui.mapMode;

% generate an angle map A and a radius map R
% A will start at the 12-o-clock and run clockwise back to 12,
% ranging from 1 to the number of colors in the cmap.
[X Y] = meshgrid(1:256, 1:256);
X = X-128; Y = Y-128;
A = atan2(X, Y);
A = fliplr(mod(A-pi, 2*pi));
A = rescale(A, [], [1 mode.numColors]);
R = sqrt(X.^2 + Y.^2);

% take cmap from color part of map mode
cmap = mode.cmap(mode.numGrays+1:end,:);

% convert A to truecolor
A = ind2rgb(A, cmap);

% for each color channel, mask out region outside radius (128 pixels)
[I J] = find(R>128);
for ch = 1:3
    ind = sub2ind(size(A), I, J, repmat(ch, size(I)));
    A(ind) = 1;
end

% put up the image
figure('Color', 'w', 'Name', 'Polar Angle Map Legend');
imshow(A);

return


text:
% function thetaOpt = StochasticGradientDescent (gradFunc, theta0, maxiter)
% runs gradient descent until convergence, returning the optimal parameters thetaOpt.
%
% Inputs:
% gradFunc      function handle to a function [cost, grad] = gradFunc(theta, i)
%               that computes the LR cost / objective function and the gradient
%               of the negative log likelihood of the ith data instance, given
%               parameters theta.
% theta0        initial value of theta to start gradient descent from.
% maxIter       number of iterations to run SGD for.
%
% Output:
% thetaOpt      optimal value of theta.
%
%
% Note - function handles may be new to some of you. Briefly, function handles
% are a way of passing a function to another function as an argument. The
% syntax for calling a function handle is exactly the same as for calling any
% other function.
%
% For more information, refer to the official documentation:
% Matlab - http://www.mathworks.com/help/techdoc/matlab_prog/f2-38133.html
% Octave - http://www.gnu.org/software/octave/doc/interpreter/Function-Handles.html
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

function thetaOpt = StochasticGradientDescent (gradFunc, theta0, maxIter)

    % The grader will accept all answers that are near enough
    % to the optimal value, so don't worry about being off by one
    % iteration etc.

    thetaOpt = zeros(size(theta0));

    %%%%%%%%%%%%%%
    %%% Student code
    for i = 1:maxIter
	    [cost, grad] = gradFunc(theta0,i);
	    theta0 = theta0 - (0.1)*(grad)/(1+sqrt(i));
    end
    thetaOpt = theta0;
    %%%%%%%%%%%

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function p = power(a,b)
	if isscalar(a)
		a = a*ones(size(b),'uint64');
	elseif isscalar(b)
		b = a*ones(size(b),'uint64');
	end
	if ~isequal(size(a),size(b))
		error('Size mismatch');
	end

	p = ones(size(a),'uint64');
	for k = 1:length(a)
		if b(k)>63
			error('Power>63');
		end
		if a(k)==0 && b(k)==0
			error('0^0');
		end
		for pow = 1:double(b(k))
			p(k) = p(k)*a(k);
		end
	end
end


text:
function rtk=udstate(rtk,obsr,obsb,nav,ind)

global glc
tt=rtk.tt;

% update the position/velocity/acceleration
rtk=udpos(rtk,tt);

% update ionosphereic parameter
if rtk.opt.ionoopt==glc.IONOOPT_EST
    bl=norm(rtk.x(1:3)-rtk.basepos);
    rtk=udion(rtk,tt,bl,ind);
end

% update tropspheric parameter
if rtk.opt.tropopt>=glc.TROPOPT_EST
    rtk=udtrop(rtk,tt,bl);
end

% update reciever inter-frequency bias for glonass
if rtk.opt.glomodear==2&&rtk.mask(2)==1
    rtk=udrcvbias(rtk,tt);
end

% update the ambiguity
if rtk.opt.mode>glc.PMODE_DGNSS
    rtk=udbias(rtk,obsr,obsb,nav,ind);
end

return

repo: GINav


text:
function [Rob,Sen,Lmk,Obs,Frm,Fac] = solveGraphCholesky(Rob,Sen,Lmk,Obs,Frm,Fac,options)

% SOLVEGRAPHCHOLESKY Solves the SLAM graph using Cholesky decomposition.
%   [Rob,Sen,Lmk,Obs,Frm,Fac] = solveGraphCholesky(Rob,Sen,Lmk,Obs,Frm,Fac)
%   solves the graph-SLAM problem using Cholesky decomposition of the
%   Hessian matrix.
%
%   IMPORTANT NOTE: This method is illustrative and constitutes the
%   motivation for this toolbox. One can achieve better performances, both
%   in computing time and possibly in robustness and accuracy, by using
%   Matlab's built-in nonlinear optimization tools, such as LSQNONLIN.
%
%   See courseSLAM.pdf in the documentation for details about the Cholesky
%   decomposition for solving the graph-SLAM problem.
%
%   See also ERRORSTATEJACOBIANS, UPDATESTATES, COMPUTEERROR,
%   COMPUTERESIDUAL, COLAMD, '\', MLDIVIDE, LSQNONLIN.

% Copyright 2015-    Joan Sola @ IRI-UPC-CSIC.

global Map

% Control of iterations and exit conditions
n_iter      = options.niterations; % exit criterion of number of iterations
target_dres = options.target_dres; % exit criterion for error variation
target_res  = options.target_res;  % exit criterion for current residual
res_old     = 1e10;                   % last iteration's error

% Map states range
Map.mr = find(Map.used);

for it = 1:n_iter

%     fprintf('----------------\nIteration: %d; \n',it)

    % Compute Jacobians for projection onto the manifold
    [Frm,Lmk] = errorStateJacobians(Frm,Lmk);

    % Build Hessian H and rhs vector b, in global Map
    Fac       = buildProblem(Rob,Sen,Lmk,Obs,Frm,Fac);

    if it == 1 % do this only once:

        % Column permutation
        p  = colamd(Map.H(Map.mr,Map.mr))';

        % Permutated map range
        pr = Map.mr(p);
    end

    % Decomposition
    [Map.R, ill] = chol(Map.H(pr,pr));

    if ill
        error('Ill-conditioned Hessian')
    end

    % Solve for dx and reorder:
    %   - dx is Map.x(mr)
    %   - reordered dx is Map.x(pr)
    y         = -Map.R'\Map.b(pr); % solve for y
    Map.x(pr) = Map.R\y;

    % NOTE: Matlab is able to do all the reordering and Cholesky
    % factorization for you. If you just use the operator '\', as in
    % 'dx = -H\b', Matlab will reorder H, then factor it as R'R, then solve
    % the two subproblems, then reorder back the result into dx. Use the
    % following line to accomplish this, and comment out the code from line
    % 'if it == 1' until here:
    %
    %     Map.x(Map.mr) = -Map.H(Map.mr,Map.mr)\Map.b(Map.mr);

    % Update nominal states
    [Rob,Lmk,Frm] = updateStates(Rob,Lmk,Frm);

    % Check resulting errors
    [res, err_max] = computeResidual(Rob,Sen,Lmk,Obs,Frm,Fac);
    dres           = res - res_old;
    res_old        = res;

    %         fprintf('Residual: %.2e; variation: %.2e \n', res, dres)

    if ( ( -dres <= target_dres ) || (err_max <= target_res) ) %&& ( abs(derr) < target_derr) )
        break;
    end

end

end

function Fac = buildProblem(Rob,Sen,Lmk,Obs,Frm,Fac)

% BUILDPROBLEM Build least squares problem's matrix H and vector b
%   Fac = BUILDPROBLEM(Rob,Sen,Lmk,Obs,Frm,Fac) Builds the least squares
%   problem's matrix H and vector b for a solution using sparse Cholesky
%   factorization of H.

global Map

% Reset Hessian and rhs vector
Map.H(Map.mr,Map.mr) = 0;
Map.b(Map.mr)    = 0;

% Iterate all factors
for fac = find([Fac.used])

    % Extract some pointers
    rob    = Fac(fac).rob;
    sen    = Fac(fac).sen;
    lmk    = Fac(fac).lmk;
    frames = Fac(fac).frames;

    % Compute factor error, info mat, and Jacobians
    [Fac(fac), e, W, ~, J1, J2, r1, r2] = computeError(...
        Rob(rob),       ...
        Sen(sen),       ...
        Lmk(lmk),       ...
        Obs(sen,lmk),   ...
        Frm(frames),    ...
        Fac(fac));

    % Compute sparse Hessian blocks
    H_11 = J1' * W * J1;
    H_12 = J1' * W * J2;
    H_22 = J2' * W * J2;

    % Compute rhs vector blocks
    b1   = J1' * W * e;
    b2   = J2' * W * e;

    % Update H and b
    Map.H(r1,r1) = Map.H(r1,r1) + H_11;
    Map.H(r1,r2) = Map.H(r1,r2) + H_12;
    Map.H(r2,r1) = Map.H(r2,r1) + H_12';
    Map.H(r2,r2) = Map.H(r2,r2) + H_22;

    Map.b(r1,1)  = Map.b(r1,1)  + b1;
    Map.b(r2,1)  = Map.b(r2,1)  + b2;

end

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [ y, m, d ] = deflate_common ( y, m, d )

%*****************************************************************************80
%
%% DEFLATE_COMMON "deflates" dates in the Common Calendar transition month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, the YMD date.
%
%    Output, integer Y, integer M, integer D, the YMD date.
%

  if ( y == 1582 )
    if ( m == 10 )
      if ( 15 <= d )
        d = d - 10;
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function view=initLevelSlider(view,curLevel)
%
%       view=initLevelSlider(view,[curLevel])
%
% Installs the slider and text box that control the gray level.
% The slider is placed on the top right, near the slice control slider.
%
% The first time through, the callback and related textbox are installed
%
% It is also called repeatedly as dataTypes are switched to re-display
%   the slider that represents level information
%
% If you change this function make parallel changes in:
%    initSliceSlider
%
% bw,  01/16/01  Adjusted sliderstep
% bw,  12/25/00  Wrote it.
% ras, 09/30/04  adapted to initLevelSlider
if ieNotDefined('curLevel')
    curLevel = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add an edit field for the # of gray levels to show %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
color = get(gcf,'Color');

% text label
htmp = uicontrol('Style','text',...
                 'Units','Normalized',...
                 'Position',[0 0.75 0.12 0.03],...
                 'BackgroundColor',color,...
                 'HorizontalAlignment','left',...
                 'String','And next:');
view.ui.level.numLevelLabel = htmp;

% edit field
callBackStr = sprintf('%s = refreshScreen(%s);',view.name,view.name);
htmp = uicontrol('Style','edit',...
                 'Units','Normalized',...
                 'Position',[0.10 0.75 0.03 0.03],...
                 'BackgroundColor',color,...
                 'Callback',callBackStr,...
                 'String','0');
view.ui.level.numLevelEdit = htmp;

% also need to get text label for level slider
% (to turn off if looking across levels)
htmp = findobj('Parent',gcf,'Style','text','String','level:');
view.ui.level.levelLabel = htmp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize slider stuff                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The sliderHandle always exists because it is created in openXXWindow
% Set the step size for clicking on the arrows and the trough
sliderHandle = view.ui.level.sliderHandle;

nScans=numScans(view);
if nScans < 2
    %   Can't have a slider value of 1.  Sorry.
    nScans=2;
    set(sliderHandle,'Enable','off');
else
    set(sliderHandle,'Enable','on');
end

hemi = findSelectedButton(view.ui.sliceButtons);
nLevels = view.numLevels(hemi); % update for diff. # levels

% By setting the max to a little more than nScans, we should
% never quite reach it by clicking the arrow or the trough
set(sliderHandle,'Min',1,'Max',nLevels);

% Clicking an arrow moves by an amount (max - min)*sliderstep(1)
% Cllicking in the trough moves by an amount (max - min)*sliderstep(2)
% setSlider traps any values == to the max value and queries the user
%    this is a cheap way of setting a value by hand.
sliderStep(1) = 1/(nLevels - 1);
sliderStep(2) = 2*sliderStep(1);
set(sliderHandle,'sliderStep',sliderStep);

% update the callback
callBackStr = ...
    ['val = get(',view.name,'.ui.level.sliderHandle,''value'');'...
        'val = round(val);'...
        'setSlider(',view.name,',',view.name,'.ui.level,val);'...
        view.name,'=refreshScreen(',view.name,');'];
set(sliderHandle,'CallBack',callBackStr);

% If we are first opening the window, create the text box.
if ~isfield(view.ui.level,'textHandle')

   % Position the text box to the right of the slider
   pos = get(sliderHandle,'Position');
   l2 = pos(1) + pos(3) + .01;
   w2 = pos(3)/4; b = pos(2); h = pos(4);
   position = [l2 b w2 h];

   view.ui.level.textHandle = ...
      uicontrol('Style','text',...
      'Units','normalized',...
      'Position',position);
end

% This text tells the user how many levels there are
% for the given gray mesh.
%
str = sprintf('%.0f',nLevels);
set(view.ui.level.textHandle,'String',str);

setSlider(view,view.ui.level,curLevel);

return


text:
function stats = extractStats(net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
stats = struct() ;
for i = 1:numel(sel)
    stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

author: aimerykong


text:
function D = driving_function_mono_nfchoa_ps(x0,xs,f,N,conf)
%DRIVING_FUNCTION_MONO_NFCHOA_PS driving signal for a point source in NFC-HOA
%
%   Usage: D = driving_function_mono_nfchoa_ps(x0,xs,f,N,conf)
%
%   Input parameters:
%       x0          - position of the secondary sources / m [nx3]
%       xs          - position of virtual point source / m [nx3]
%       f           - frequency of the monochromatic source / Hz
%       N           - maximum order of spherical harmonics
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   See also: driving_function_mono_nfchoa, driving_function_imp_nfchoa_ps

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargmatrix(x0,xs);
isargpositivescalar(f,N);
isargstruct(conf);

%% ===== Configuration ==================================================
xref = conf.xref;
c = conf.c;
dimension = conf.dimension;
driving_functions = conf.driving_functions;
X0 = conf.secondary_sources.center;

%% ===== Computation ====================================================

% Secondary source positions
x00 = bsxfun(@minus,x0,X0);
[phi0,theta0,r0] = cart2sph(x00(:,1),x00(:,2),x00(:,3));

% Point source
xs0 = bsxfun(@minus,xs,X0);
[phi,theta,r] = cart2sph(xs0(:,1),xs0(:,2),xs0(:,3));

% Wavenumber
omega = 2*pi*f;

% modal window
win = modal_weighting(N,conf);

% Initialize empty driving signal
D = zeros(size(x0,1),1);

if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2D point source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    % Reference point
    xref = repmat(xref,[size(x0,1) 1]);

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % 2.5D point source
        %
        %                     _N_    (2)
        %               1     \     h|m|(w/c r)
        % D(phi0,w) = ------  /__  ------------- e^(i m (phi0-phi))
        %             2pi r0  m=-N   (2)
        %                           h|m|(w/c r0)
        %
        % https://sfs.rtfd.io/en/3.2/d_nfchoa/#equation-fd-nfchoa-point-25d
        for m=-N:N
            D = D + 1 ./ (2.*pi.*r0) ...
                .* win(abs(m)+1) ...
                .* sphbesselh(abs(m),2,omega./c.*r) ...
                ./ sphbesselh(abs(m),2,omega./c.*r0) ...
                .* exp(1i.*m.*(phi0-phi));
        end
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2.5D point source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('3D',dimension)

    % === 3-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % 3D point source
        %
        %                              _N_  _n_   (2)
        %                       1      \    \    hn(w/c r)   -m
        % D(theta0,phi0,w) = -------   /__  /__ ----------- Yn(theta,phi) ...
        %                    2pi r0^2  n=0 m=-n   (2)
        %                                        hn(w/c r0)
        %                       m
        %                    x Yn(theta0,phi0)
        %
        % https://sfs.rtfd.io/en/3.2/d_nfchoa/#equation-fd-nfchoa-point-3d
        %
        for n=0:N
            for m=-n:n
                D = D + 1 ./ (2.*pi.*r0.^2) ...
                    .* win(n+1) ...
                    .* sphbesselh(n,2,omega./c.*r) ...
                    ./ sphbesselh(n,2,omega./c.*r0) ...
                    .* sphharmonics(n,-m,theta,phi) ...
                    .* sphharmonics(n,m,theta0,phi0);
            end
        end
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 3D point source.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function out=acos(x)

precision=x(1).precision;
mpPi=mppi(precision);
out=mpPi/2+i.*log(i*x+sqrt(1-x.^2));

%%%precision=x(1).precision;
%%%out=mp(zeros(size(x)));
%%%
%%%mpPi=mppi(precision);
%%%for ii=1:numel(x)
%%% out(ii)=mpPi/2+i*log(i*x(ii)+sqrt(1-x(ii)^2));
%%%end % for ii=1:max(ex,

author: opencobra

repo: cobratoolbox


text:
function alpha = normalizeAngle(alpha, varargin)
%NORMALIZEANGLE  Normalize an angle value within a 2*PI interval.
%
%   ALPHA2 = normalizeAngle(ALPHA);
%   ALPHA2 is the same as ALPHA modulo 2*PI and is positive.
%
%   ALPHA2 = normalizeAngle(ALPHA, CENTER);
%   Specifies the center of the angle interval.
%   If CENTER==0, the interval is [-pi ; +pi]
%   If CENTER==PI, the interval is [0 ; 2*pi] (default).
%
%   Example:
%   % normalization between 0 and 2*pi (default)
%   normalizeAngle(5*pi)
%   ans =
%       3.1416
%
%   % normalization between -pi and +pi
%   normalizeAngle(7*pi/2, 0)
%   ans =
%       -1.5708
%
%   See also
%   vectorAngle, lineAngle
%

% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2008-03-10,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% HISTORY
% 2010-03-31 rename as normalizeAngle, and add psb to specify interval
%   center

center = pi;
if ~isempty(varargin)
    center = varargin{1};
end

alpha = mod(alpha-center+pi, 2*pi) + center-pi;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function exact = p00_exact ( prob )

%*****************************************************************************80
%
%% P00_EXACT returns the exact integral for any problem.
%
%  Discussion:
%
%    This routine provides a "generic" interface to the exact integral
%    routines for the various problems, and allows a problem to be called
%    by number (PROB) rather than by name.
%
%    In some cases, the "exact" value of the integral is in fact
%    merely a respectable approximation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the desired test problem.
%
%    Output, real EXACT, the exact value of the integral.
%
  if ( prob == 1 )
    exact = p01_exact ( );
  elseif ( prob == 2 )
    exact = p02_exact ( );
  elseif ( prob == 3 )
    exact = p03_exact ( );
  elseif ( prob == 4 )
    exact = p04_exact ( );
  elseif ( prob == 5 )
    exact = p05_exact ( );
  elseif ( prob == 6 )
    exact = p06_exact ( );
  elseif ( prob == 7 )
    exact = p07_exact ( );
  elseif ( prob == 8 )
    exact = p08_exact ( );
  elseif ( prob == 9 )
    exact = p09_exact ( );
  elseif ( prob == 10 )
    exact = p10_exact ( );
  elseif ( prob == 11 )
    exact = p11_exact ( );
  elseif ( prob == 12 )
    exact = p12_exact ( );
  elseif ( prob == 13 )
    exact = p13_exact ( );
  elseif ( prob == 14 )
    exact = p14_exact ( );
  elseif ( prob == 15 )
    exact = p15_exact ( );
  elseif ( prob == 16 )
    exact = p16_exact ( );
  elseif ( prob == 17 )
    exact = p17_exact ( );
  elseif ( prob == 18 )
    exact = p18_exact ( );
  elseif ( prob == 19 )
    exact = p19_exact ( );
  elseif ( prob == 20 )
    exact = p20_exact ( );
  elseif ( prob == 21 )
    exact = p21_exact ( );
  elseif ( prob == 22 )
    exact = p22_exact ( );
  elseif ( prob == 23 )
    exact = p23_exact ( );
  elseif ( prob == 24 )
    exact = p24_exact ( );
  elseif ( prob == 25 )
    exact = p25_exact ( );
  elseif ( prob == 26 )
    exact = p26_exact ( );
  elseif ( prob == 27 )
    exact = p27_exact ( );
  elseif ( prob == 28 )
    exact = p28_exact ( );
  elseif ( prob == 29 )
    exact = p29_exact ( );
  elseif ( prob == 30 )
    exact = p30_exact ( );
  elseif ( prob == 31 )
    exact = p31_exact ( );
  elseif ( prob == 32 )
    exact = p32_exact ( );
  elseif ( prob == 33 )
    exact = p33_exact ( );
  elseif ( prob == 34 )
    exact = p34_exact ( );
  elseif ( prob == 35 )
    exact = p35_exact ( );
  elseif ( prob == 36 )
    exact = p36_exact ( );
  elseif ( prob == 37 )
    exact = p37_exact ( );
  elseif ( prob == 38 )
    exact = p38_exact ( );
  elseif ( prob == 39 )
    exact = p39_exact ( );
  elseif ( prob == 40 )
    exact = p40_exact ( );
  elseif ( prob == 41 )
    exact = p41_exact ( );
  elseif ( prob == 42 )
    exact = p42_exact ( );
  elseif ( prob == 43 )
    exact = p43_exact ( );
  elseif ( prob == 44 )
    exact = p44_exact ( );
  elseif ( prob == 45 )
    exact = p45_exact ( );
  elseif ( prob == 46 )
    exact = p46_exact ( );
  elseif ( prob == 47 )
    exact = p47_exact ( );
  elseif ( prob == 48 )
    exact = p48_exact ( );
  elseif ( prob == 49 )
    exact = p49_exact ( );
  elseif ( prob == 50 )
    exact = p50_exact ( );
  elseif ( prob == 51 )
    exact = p51_exact ( );
  elseif ( prob == 52 )
    exact = p52_exact ( );
  elseif ( prob == 53 )
    exact = p53_exact ( );
  elseif ( prob == 54 )
    exact = p54_exact ( );
  elseif ( prob == 55 )
    exact = p55_exact ( );
  elseif ( prob == 56 )
    exact = p56_exact ( );
  elseif ( prob == 57 )
    exact = p57_exact ( );
  elseif ( prob == 58 )
    exact = p58_exact ( );
  elseif ( prob == 59 )
    exact = p59_exact ( );
  elseif ( prob == 60 )
    exact = p60_exact ( );
  elseif ( prob == 61 )
    exact = p61_exact ( );
  elseif ( prob == 62 )
    exact = p62_exact ( );
  elseif ( prob == 63 )
    exact = p63_exact ( );
  elseif ( prob == 64 )
    exact = p64_exact ( );
  elseif ( prob == 65 )
    exact = p65_exact ( );
  elseif ( prob == 66 )
    exact = p66_exact ( );
  elseif ( prob == 67 )
    exact = p67_exact ( );
  elseif ( prob == 68 )
    exact = p68_exact ( );
  elseif ( prob == 69 )
    exact = p69_exact ( );
  elseif ( prob == 70 )
    exact = p70_exact ( );
  elseif ( prob == 71 )
    exact = p71_exact ( );
  elseif ( prob == 72 )
    exact = p72_exact ( );
  elseif ( prob == 73 )
    exact = p73_exact ( );
  elseif ( prob == 74 )
    exact = p74_exact ( );
  elseif ( prob == 75 )
    exact = p75_exact ( );
  elseif ( prob == 76 )
    exact = p76_exact ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_EXACT - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_EXACT - Fatal error!' );
  end

  return
end


text:
classdef SparseEA2 < ALGORITHM
% <multi> <real/integer/binary> <large/none> <constrained/none> <sparse>
% Improved SparseEA

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

    methods
        function main(Algorithm,Problem)
            %% Population initialization
            % Calculate the fitness of each decision variable
            TDec    = [];
            TMask   = [];
            TempPop = [];
            Fitness = zeros(1,Problem.D);
            for i = 1 : 1+4*any(Problem.encoding~=4)
                Dec = unifrnd(repmat(Problem.lower,Problem.D,1),repmat(Problem.upper,Problem.D,1));
                Dec(:,Problem.encoding==4) = 1;
                Mask       = eye(Problem.D);
                Population = Problem.Evaluation(Dec.*Mask);
                TDec       = [TDec;Dec];
                TMask      = [TMask;Mask];
                TempPop    = [TempPop,Population];
                Fitness    = Fitness + NDSort([Population.objs,Population.cons],inf);
            end
            % Generate initial population
            Dec = unifrnd(repmat(Problem.lower,Problem.N,1),repmat(Problem.upper,Problem.N,1));
            Dec(:,Problem.encoding==4) = 1;
            Mask = false(Problem.N,Problem.D);
            for i = 1 : Problem.N
                Mask(i,TournamentSelection(2,ceil(rand*Problem.D),Fitness)) = 1;
            end
            Population = Problem.Evaluation(Dec.*Mask);
            [Population,Dec,Mask,FrontNo,CrowdDis] = EnvironmentalSelection([Population,TempPop],[Dec;TDec],[Mask;TMask],Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool       = TournamentSelection(2,2*Problem.N,FrontNo,-CrowdDis);
                [OffDec,OffMask] = Operator(Problem,Dec(MatingPool,:),Mask(MatingPool,:),Fitness);
                Offspring        = Problem.Evaluation(OffDec.*OffMask);
                [Population,Dec,Mask,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],[Dec;OffDec],[Mask;OffMask],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ gpos ] = point_image_to_ground( im_points, sicd_meta, varargin )
%POINT_IMAGE_TO_GROUND Transforms pixel row, col to ground ECF coordinate
% This function implements the SICD Image Projections Description Document:
% http://www.gwg.nga.mil/ntb/baseline/docs/SICD/index.html
%
% gpos = point_image_to_ground(im_points, sicd_meta, 'PropertyName', PropertyValue, ...)
%
% Inputs:
%    im_points   - [2xN] (row; column) coordinates of N points in image (or
%                  subimage if FirstRow/FirstCol are nonzero).  Zero-based,
%                  following SICD convention (rather than MATLAB
%                  convention, which is one-based); that is, upper-left
%                  pixel is [0;0].
%    sicd_meta   - SICD meta data structure
%
%       Property name     Description
%       projection_type   'plane', 'hae', 'dem', 'best' (try DEM first
%                         then HAE).  Default plane.
%       gref              Ground plane reference point ECF coordinates (m).
%                         Default is SCP.  Only valid if projection_type is
%                         'plane'.
%       ugpn              Ground plane unit normal vector.  Default is
%                         tangent to the surface of constant geodetic
%                         height above the WGS-84 reference ellipsoid
%                         passing through GREF.  Only valid if
%                         projection_type is 'plane'.
%       hae0              Surface height (m) above the WGS-84 reference
%                         ellipsoid for projection point.  Only valid if
%                         projection_type is 'hae'.
%       delta_hae_max     Height threshold for convergence of iterative
%                         constant HAE computation (m).  Default 1 meter.
%                         Only valid if projection_type is 'hae'.
%       hae_nlim          Maximum number of iterations allowed for constant
%                         hae computation.  Default 3.  Only valid if
%                         projection_type is 'hae'.
%       dem               SRTM pathname or structure with
%                         lats/lons/elevations fields where are all are
%                         arrays of same size and elevation is height above
%                         WGS-84 ellipsoid.
%                         Only valid if projection_type is 'dem'.
%       del_DISTrrc       Maximum distance between adjacent points along
%                         the R/Rdot contour. Recommended value: 10.0 m.
%                         Only valid if projection_type is 'dem'.
%       del_HDlim         Height difference threshold for determining if a
%                         point on the R/Rdot contour is on the DEM
%                         surface (m).  Recommended value: 0.001 m.  Only
%                         valid if projection_type is 'dem'.
%       delta_arp         ARP position adjustable parameter (ECF, m).  Default 0.
%       delta_varp        VARP position adjustable parameter (ECF, m/s).  Default 0.
%       range_bias        Range bias adjustable parameter (m).  Default 0.
%       adj_params_frame  Coordinate frame used for expressing delta_arp
%                         and delta_varp adjustable parameters.  Allowed
%                         values: 'ECF', 'RIC_ECF', 'RIC_ECI'. Default ECF.
%
% Outputs:
%    gpos        - [3xN] ECF Ground Points along the R/Rdot contour
%
% Contributors: Thomas McDowall, Harris Corporation
%               Lisa Talbot, NGA/IB
%               Rocco Corsetti, NGA/IB
%               Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Parse input parameters
if isvector(im_points), im_points = im_points(:); end % Assure orientation
scp = [sicd_meta.GeoData.SCP.ECF.X; ...
    sicd_meta.GeoData.SCP.ECF.Y; ...
    sicd_meta.GeoData.SCP.ECF.Z]; % Default values for some input arguments require this
scp_llh = ecf_to_geodetic(scp);  % Also should be available in sicd_meta.GeoData.SCP.LLH
scp_hae = scp_llh(3);
p = inputParser;
p.addParamValue('projection_type','plane', @(x) any(strcmp(x,{'plane','hae','hae_newton','dem','best'})));
p.addParamValue('gref',scp, @(x) (size(x,1)==3)||(numel(x)==3)); % ECF (meters)
if isfield(sicd_meta,'PFA')&&isfield(sicd_meta.PFA,'FPN')
    default_ugpn=[sicd_meta.PFA.FPN.X; sicd_meta.PFA.FPN.Y; sicd_meta.PFA.FPN.Z];
else
    default_ugpn=wgs_84_norm(scp);
end
p.addParamValue('ugpn',default_ugpn, @(x) (size(x,1)==3)||(numel(x)==3)); % ECF
p.addParamValue('hae0',scp_hae, @isscalar); % meters
p.addParamValue('delta_hae_max', 1, @(x) isscalar(x) && (x>0)); % meters
p.addParamValue('hae_nlim', 3, @isscalar);
p.addParamValue('dem',[]);
p.addParamValue('del_DISTrrc', 10, @isscalar);
p.addParamValue('del_HDlim', 0.001, @isscalar);
% Adjustable parameters
p.addParamValue('delta_arp',[0 0 0], @(x) numel(x)==3); % ECF? (meters)
p.addParamValue('delta_varp',[0 0 0], @(x) numel(x)==3); % ECF? (meters/s)
p.addParamValue('range_bias',0, @isscalar); % meters
p.addParamValue('adj_params_frame','ECF', @(x) any(strcmpi(x,{'ECF','RIC_ECF','RIC_ECI'})));
p.FunctionName = mfilename;
p.parse(varargin{:});

%% Chapter 4: Compute projection parameters
try
    % Computation of r/rdot is specific to image formation type
    [r, rdot, arp_coa, varp_coa] = coa_projection_set(sicd_meta, im_points);
catch % If full metadata not available, try to make the approximation of uniformly sampled grid
    sicd_meta.Grid.Type = 'PLANE';
    if ~isfield(sicd_meta.Grid, 'TimeCOAPoly') % Another approximation that may have to be made
        sicd_meta.Grid.TimeCOAPoly = sicd_meta.Timeline.CollectDuration/2;
    end
    [r, rdot, arp_coa, varp_coa] = coa_projection_set(sicd_meta, im_points);
    warning('point_image_to_ground:IncompleteMetadata',...
        'Unable to compute precise position due to incomplete metadata.  Resorting to approximation.');
end
% After r/rdot is computed, the rest is generic to all SAR.

%% Apply adjustable parameters
if strcmpi(p.Results.adj_params_frame,'ECF')
    % No transformation necessary
    delta_arp = p.Results.delta_arp(:);
    delta_varp = p.Results.delta_varp(:);
else % Translate from RIC frame to ECF frame
    % Use the RIC frame at SCP COA time, not at COA time for im_points
    ARP_SCP_COA = [sicd_meta.SCPCOA.ARPPos.X; ...
        sicd_meta.SCPCOA.ARPPos.Y; ...
        sicd_meta.SCPCOA.ARPPos.Z];
    VARP_SCP_COA = [sicd_meta.SCPCOA.ARPVel.X; ...
        sicd_meta.SCPCOA.ARPVel.Y; ...
        sicd_meta.SCPCOA.ARPVel.Z];
    if strcmpi(p.Results.adj_params_frame,'RIC_ECI')
        T_ECEF_RIC = ric_ecf_mat(ARP_SCP_COA, VARP_SCP_COA, 'eci');
    else % RIC_ECF
        T_ECEF_RIC = ric_ecf_mat(ARP_SCP_COA, VARP_SCP_COA, 'ecf');
    end
    delta_arp = T_ECEF_RIC * p.Results.delta_arp(:);
    delta_varp = T_ECEF_RIC * p.Results.delta_varp(:);
end
arp_coa = arp_coa + repmat(delta_arp,1,size(arp_coa,2));
varp_coa = varp_coa + repmat(delta_varp,1,size(arp_coa,2));
r = r + p.Results.range_bias;

%% Perform actual projection
switch p.Results.projection_type
    case 'plane' % Chapter 5 of SICD document
        gpos = point_to_ground_plane(r, rdot, arp_coa, varp_coa, ...
            p.Results.gref, p.Results.ugpn);
    case 'hae' % Chapter 9 of SICD projection document (draft)
        gpos = point_to_hae(r, rdot, arp_coa, varp_coa, scp, ...
            p.Results.hae0, p.Results.delta_hae_max, p.Results.hae_nlim);
    case 'hae_newton' % Spotlight Synthetic Aperture Radar (SAR) Sensor Model, NGA.SIG.0005_1.0, 2010-03-30
        if ~strcmp(sicd_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
            error('point_image_to_ground:invalid_projection_mode','HAE_NEWTON projection method only valid for spotlight data.');
        end
        gpos = point_to_hae_newton(r, rdot, arp_coa, varp_coa, scp, ...
            p.Results.hae0);
    case 'dem' % Chapter 10 of SICD document (draft)
        gpos = point_to_DEM(r, rdot, arp_coa, varp_coa, scp, ...
            p.Results.delta_hae_max, p.Results.hae_nlim, p.Results.dem, ...
            p.Results.del_DISTrrc, p.Results.del_HDlim);
    case 'best' % Best available.  Use DEM if valid, otherwise use HAE
        try
            gpos = point_to_DEM(r, rdot, arp_coa, varp_coa, scp, ...
                p.Results.delta_hae_max, p.Results.hae_nlim, p.Results.dem, ...
                p.Results.del_DISTrrc, p.Results.del_HDlim);
        catch
            warning('point_image_to_ground:invalid_dem','Best projection reverted to HAE.');
            gpos = point_to_hae(r, rdot, arp_coa, varp_coa, scp, ...
                p.Results.hae0, p.Results.delta_hae_max, p.Results.hae_nlim);
        end
    otherwise
        % Unrecognized projection type
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
%KNAPSACK Solves the 0-1 knapsack problem for positive integer weights
%
%  [BEST AMOUNT] = KNAPSACK(WEIGHTS, VALUES, CONSTRAINT)
%
%       WEIGHTS    : The weight of every item (1-by-N)
%       VALUES     : The value of every item (1-by-N)
%       CONSTRAINT : The weight constraint of the knapsack (scalar)
%
%       BEST       : Value of best possible knapsack (scalar)
%       AMOUNT     : 1-by-N vector specifying the amount to use of each item (0 or 1)
%
%
%	EXAMPLE :
%
%		weights = [1 1 1 1 2 2 3];
%		values  = [1 1 2 3 1 3 5];
%		[best amount] = KNAPSACK(weights, values, 7)
%
%		best =
%
%	    	13
%
%
%		amount =
%
%	    	 0     0     1     1     0     1     1
%
%
%   See Knapsack problem on Wikipedia.
%
%   Copyright 2009 Petter Strandmark
%   petter.strandmark@gmail.com
function [best amount] = knapsack(weights, values, W)
    if ~all(is_positive_integer(weights)) || ...
       ~is_positive_integer(W)
        error('Weights must be positive integers');
    end
    %We work in one dimension
    [M N] = size(weights);
    weights = weights(:);
    values = values(:);
    if numel(weights) ~= numel(values)
        error('The size of weights must match the size of values');
    end
    if numel(W) > 1
        error('Only one constraint allowed');
    end

    % Solve the problem

    % Note that A would ideally be indexed from A(0..N,0..W) but MATLAB
    % does not allow this.
    A = zeros(length(weights)+1,W+1);
    % A(j+1,Y+1) means the value of the best knapsack with capacity Y using
    % the first j items.
    for j = 1:length(weights)
        for Y = 1:W
            if weights(j) > Y
                A(j+1,Y+1) = A(j,Y+1);
            else
                A(j+1,Y+1) = ...
                    max( A(j,Y+1), values(j) + A(j,Y-weights(j)+1));
            end
        end
    end

   best = A(end,end);

   %Now backtrack
   amount = zeros(length(weights),1);
   a = best;
   j = length(weights);
   Y = W;
   while a > 0
       while A(j+1,Y+1) == a
           j = j - 1;
       end
       j = j + 1; %This item has to be in the knapsack
       amount(j) = 1;
       Y = Y - weights(j);
       j = j - 1;
       a = A(j+1,Y+1);
   end

    amount = reshape(amount,M,N);
end

function yn = is_positive_integer(X)
    yn = X>0 & floor(X)==X;
end

author: kezhang-cs

repo: Video-Summarization-with-LSTM


text:
function ui = mrViewGrayscalePanel(ui, dockFlag, vis);
%
% ui = mrViewGrayscalePanel(ui, [dockFlag], [vis=0]);
%
% Add a panel for setting the grayscale display
% values of the mr image (underlay). Also initializes
% the panel to be hidden, since it may not be needed
% right away.
%
% dockFlag: if 1, will attach panel to the
% mrViewer figure. Otherwise, will make a
% separate figure for the panel.
%
% vis: if 0, will initialize the panel to be hidden (default when
% starting up mrViewer); otherwise, leave it visible (as when
% docking/undocking the panel).
%
% ras, 07/05.
% ras, 07/15/05 -- now places the panel in a separate figure
% by default.
% ras, 06/07 -- now a horizontal panel; added vis toggle.
if ~exist('ui','var') | isempty(ui), ui = mrViewGet; end
if ~exist('dockFlag','var') | isempty(dockFlag), dockFlag = 0; end
if ~exist('vis', 'var') | isempty(vis),  vis = 0; end

% add the panel
if dockFlag==1
    ui.panels.grayscale = mrvPanel('above', .2, ui.fig, 'normalized');
else
    hfig = figure('Name', sprintf('Grayscale [%s]',ui.tag),...
        'Units', 'normalized', 'Position', [.2 .4 .12 .3],...
        'MenuBar', 'none', 'NumberTitle', 'off', ...
        'UserData', ui.tags);
    ui.panels.grayscale = uipanel('Parent',hfig,'Units','normalized',...
        'Position',[0 0 1 1]);

	% close request function for this figure: just toggle visibility
    crf = 'tmp = findobj(''Parent'', gcf, ''Type'', ''uipanel''); ';
    crf = [crf 'mrvPanelToggle(tmp); clear tmp;'];
    set(hfig, 'CloseRequestFcn', crf);
end

set(ui.panels.grayscale, 'BackgroundColor', [.9 .9 .9],...
    'ShadowColor', [1 1 1], 'Title', 'Grayscale');

% add a slider for brightness
ui.controls.brightness = mrvSlider([.02 .2 .2 .5], 'Brightness', ...
    'Callback', 'mrViewSetGrayscale;', 'Range', [-1 1], 'Val', 0,...
    'Parent', ui.panels.grayscale);

% add a toggle button for manually setting clip values
cb = 'mrViewSetGrayscale(gcf, ''manualclip'', gcbo); ';
ui.controls.manualClip = mrvButtons('Manual Clip Values',...
    [.25 .55 .15 .4], ui.panels.grayscale, cb, 0);

% add a button to guess clip val
cb = 'mrViewSetGrayscale([],''clip'',''guess'');';
ui.controls.guessClip = uicontrol('Parent', ui.panels.grayscale, ...
    'Style', 'pushbutton', ...
    'String', 'Guess Clip Vals', ...
    'Units', 'normalized', ...
    'Position', [.25 .05 .15 .4], ...
    'Callback', cb);

% add sliders for clip min / max
rng = [min(ui.mr.data(:)) max(ui.mr.data(:))];
ui.controls.clipMin = mrvSlider([.45 .55 .2 .4], 'Clip Min Val',...
    'Callback', 'mrViewSetGrayscale;', ...
    'Range', rng, 'Val', rng(1),...
    'Visible', 'off', ...
    'Parent', ui.panels.grayscale);

ui.controls.clipMax = mrvSlider([.45 .05 .2 .4], 'Clip Max Val',...
    'Callback', 'mrViewSetGrayscale;',...
    'Range', rng, 'Val', rng(2),...
    'Visible', 'off',...
    'Parent', ui.panels.grayscale);

% hide the panel
mrvPanelToggle(ui.panels.grayscale,'off');

%%%%%hide the panel
if vis==0
	mrvPanelToggle(ui.panels.grayscale, 'off');
end

return


text:
% EEG_INTERP - interpolate data channels
%
% Usage: EEGOUT = eeg_interp(EEG, badchans, method, t_range);
%
% Inputs:
%     EEG      - EEGLAB dataset
%     badchans - [integer array] indices of channels to interpolate.
%                For instance, these channels might be bad.
%                [chanlocs structure] channel location structure containing
%                either locations of channels to interpolate or a full
%                channel structure (missing channels in the current
%                dataset are interpolated).
%     method   - [string] method used for interpolation (default is 'spherical').
%                'invdist'/'v4' uses inverse distance on the scalp
%                'spherical' uses superfast spherical interpolation.
%                'spacetime' uses griddata3 to interpolate both in space
%                and time (very slow and cannot be interrupted).
%     t_range  - [integer array with just two elements] time interval of the
%                badchans which should be interpolated. First element is
%                the start time and the second element is the end time.
% Output:
%     EEGOUT   - data set with bad electrode data replaced by
%                interpolated data
%
% Author: Arnaud Delorme, CERCO, CNRS, Mai 2006-

% Copyright (C) Arnaud Delorme, CERCO, 2006, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function EEG = eeg_interp(ORIEEG, bad_elec, method, t_range)

    if nargin < 2
        help eeg_interp;
        return;
    end
    EEG = ORIEEG;

    if nargin < 3
        disp('Using spherical interpolation');
        method = 'spherical';
    end

    if nargin < 4
        t_range = [ORIEEG.xmin ORIEEG.xmax];
    end

    % check channel structure
    tmplocs = ORIEEG.chanlocs;
    if isempty(tmplocs) || isempty([tmplocs.X])
        error('Interpolation require channel location');
    end

    if isstruct(bad_elec)

        % remove duplicate channel labels
        % -------------------------------
        allLabels = { bad_elec.labels };
        if length(unique(allLabels)) ~= length(allLabels)
            [~,uniqueInd] = unique(allLabels);
            bad_elec = bad_elec(uniqueInd);
        end

        % add missing channels in interpolation structure
        % -----------------------------------------------
        lab1 = { bad_elec.labels };
        tmpchanlocs = EEG.chanlocs;
        lab2 = { tmpchanlocs.labels };
        [~, tmpchan] = setdiff_bc( lab2, lab1);
        tmpchan = sort(tmpchan);

        % From 'bad_elec' using only fields present on EEG.chanlocs
        fields = fieldnames(bad_elec);
        [~, indx1] = setxor(fields,fieldnames(EEG.chanlocs)); clear tmp;
        if ~isempty(indx1)
            bad_elec = rmfield(bad_elec,fields(indx1));
            fields = fieldnames(bad_elec);
        end

        if ~isempty(tmpchan)
            newchanlocs = [];
            for index = 1:length(fields)
                if isfield(bad_elec, fields{index})
                    for cind = 1:length(tmpchan)
                        fieldval = getfield( EEG.chanlocs, { tmpchan(cind) },  fields{index});
                        newchanlocs = setfield(newchanlocs, { cind }, fields{index}, fieldval);
                    end
                end
            end
            newchanlocs(end+1:end+length(bad_elec)) = bad_elec;
            bad_elec = newchanlocs;
        end
        if length(EEG.chanlocs) == length(bad_elec), return; end

        lab1 = { bad_elec.labels };
        tmpchanlocs = EEG.chanlocs;
        lab2 = { tmpchanlocs.labels };
        [~, badchans] = setdiff_bc( lab1, lab2);
        fprintf('Interpolating %d channels...\n', length(badchans));
        if isempty(badchans), return; end
        goodchans      = sort(setdiff(1:length(bad_elec), badchans));

        % re-order good channels
        % ----------------------
        [~, tmp2, neworder] = intersect_bc( lab1, lab2 );
        [~, ordertmp2] = sort(tmp2);
        neworder = neworder(ordertmp2);
        EEG.data = EEG.data(neworder, :, :);

        % looking at channels for ICA
        % ---------------------------
        %[tmp sorti] = sort(neworder);
        %{ EEG.chanlocs(EEG.icachansind).labels; bad_elec(goodchans(sorti(EEG.icachansind))).labels }

        % update EEG dataset (add blank channels)
        % ---------------------------------------
        if ~isempty(EEG.icasphere)

            [~, sorti] = sort(neworder);
            EEG.icachansind = sorti(EEG.icachansind);
            EEG.icachansind = goodchans(EEG.icachansind);
            EEG.chaninfo.icachansind = EEG.icachansind;

            % TESTING SORTING
            %icachansind = [ 3 4 5 7 8]
            %data = round(rand(8,10)*10)
            %neworder = shuffle(1:8)
            %data2 = data(neworder,:)
            %icachansind2 = sorti(icachansind)
            %data(icachansind,:)
            %data2(icachansind2,:)
        end
        % { EEG.chanlocs(neworder).labels; bad_elec(sort(goodchans)).labels }
        %tmpdata                  = zeros(length(bad_elec), size(EEG.data,2), size(EEG.data,3));
        %tmpdata(goodchans, :, :) = EEG.data;

        % looking at the data
        % -------------------
        %tmp1 = mattocell(EEG.data(sorti,1));
        %tmp2 = mattocell(tmpdata(goodchans,1));
        %{ EEG.chanlocs.labels; bad_elec(goodchans).labels; tmp1{:}; tmp2{:} }
        %EEG.data      = tmpdata;

        EEG.chanlocs  = bad_elec;

    else
        badchans  = bad_elec;
        goodchans = setdiff_bc(1:EEG.nbchan, badchans);
        if strcmpi(method, 'sphericalfast')
            EEG.data(badchans,:) = [];
            EEG.nbchan = length(goodchans);
        else
            oldelocs  = EEG.chanlocs;
            EEG       = pop_select(EEG, 'nochannel', badchans);
            EEG.chanlocs = oldelocs;
            disp('Interpolating missing channels...');
        end
    end

    % find non-empty good channels
    % ----------------------------
    origoodchans = goodchans;
    chanlocs     = EEG.chanlocs;
    nonemptychans = find(~cellfun('isempty', { chanlocs.theta }));
    [tmp indgood ] = intersect_bc(goodchans, nonemptychans);
    goodchans = goodchans( sort(indgood) );
    datachans = getdatachans(goodchans,badchans);
    badchans  = intersect_bc(badchans, nonemptychans);
    if isempty(badchans), return; end

    % scan data points
    % ----------------
    if strcmpi(method, 'spherical') || strcmpi(method, 'sphericalfast')
        % get theta, rad of electrodes
        % ----------------------------
        tmpgoodlocs = EEG.chanlocs(goodchans);
        xelec = [ tmpgoodlocs.X ];
        yelec = [ tmpgoodlocs.Y ];
        zelec = [ tmpgoodlocs.Z ];
        rad = sqrt(xelec.^2+yelec.^2+zelec.^2);
        xelec = xelec./rad;
        yelec = yelec./rad;
        zelec = zelec./rad;
        tmpbadlocs = EEG.chanlocs(badchans);
        xbad = [ tmpbadlocs.X ];
        ybad = [ tmpbadlocs.Y ];
        zbad = [ tmpbadlocs.Z ];
        rad = sqrt(xbad.^2+ybad.^2+zbad.^2);
        xbad = xbad./rad;
        ybad = ybad./rad;
        zbad = zbad./rad;
        if isempty(xbad)
            error('Trying to interpolate a channel without coordinates assigned to it');
        end

        EEG.data = reshape(EEG.data, EEG.nbchan, EEG.pnts*EEG.trials);
        %[tmp1 tmp2 tmp3 tmpchans] = spheric_spline_old( xelec, yelec, zelec, EEG.data(goodchans,1));
        %max(tmpchans(:,1)), std(tmpchans(:,1)),
        %[tmp1 tmp2 tmp3 EEG.data(badchans,:)] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, EEG.data(goodchans,:));
        [tmp1 tmp2 tmp3 badchansdata] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, EEG.data(datachans,:));
        %max(EEG.data(goodchans,1)), std(EEG.data(goodchans,1))
        %max(EEG.data(badchans,1)), std(EEG.data(badchans,1))
        EEG.data = reshape(EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
    elseif strcmpi(method, 'spacetime') % 3D interpolation, works but x10 times slower
        disp('Warning: if processing epoch data, epoch boundary are ignored...');
        disp('3-D interpolation, this can take a long (long) time...');
        tmpgoodlocs = EEG.chanlocs(goodchans);
        tmpbadlocs = EEG.chanlocs(badchans);
        [xbad ,ybad]  = pol2cart([tmpbadlocs.theta],[tmpbadlocs.radius]);
        [xgood,ygood] = pol2cart([tmpgoodlocs.theta],[tmpgoodlocs.radius]);
        pnts = size(EEG.data,2)*size(EEG.data,3);
        zgood = [1:pnts];
        zgood = repmat(zgood, [length(xgood) 1]);
        zgood = reshape(zgood,prod(size(zgood)),1);
        xgood = repmat(xgood, [1 pnts]); xgood = reshape(xgood,prod(size(xgood)),1);
        ygood = repmat(ygood, [1 pnts]); ygood = reshape(ygood,prod(size(ygood)),1);
        tmpdata = reshape(EEG.data, prod(size(EEG.data)),1);
        zbad = 1:pnts;
        zbad = repmat(zbad, [length(xbad) 1]);
        zbad = reshape(zbad,prod(size(zbad)),1);
        xbad = repmat(xbad, [1 pnts]); xbad = reshape(xbad,prod(size(xbad)),1);
        ybad = repmat(ybad, [1 pnts]); ybad = reshape(ybad,prod(size(ybad)),1);
        badchansdata = griddata3(ygood, xgood, zgood, tmpdata,...
                                              ybad, xbad, zbad, 'nearest'); % interpolate data
    else
        % get theta, rad of electrodes
        % ----------------------------
        tmpchanlocs = EEG.chanlocs;
        [xbad ,ybad]  = pol2cart([tmpchanlocs( badchans).theta],[tmpchanlocs( badchans).radius]);
        [xgood,ygood] = pol2cart([tmpchanlocs(goodchans).theta],[tmpchanlocs(goodchans).radius]);

        fprintf('Points (/%d):', size(EEG.data,2)*size(EEG.data,3));
        badchansdata = zeros(length(badchans), size(EEG.data,2)*size(EEG.data,3));

        for t=1:(size(EEG.data,2)*size(EEG.data,3)) % scan data points
            if mod(t,100) == 0, fprintf('%d ', t); end
            if mod(t,1000) == 0, fprintf('\n'); end;

            %for c = 1:length(badchans)
            %   [h EEG.data(badchans(c),t)]= topoplot(EEG.data(goodchans,t),EEG.chanlocs(goodchans),'noplot', ...
            %        [EEG.chanlocs( badchans(c)).radius EEG.chanlocs( badchans(c)).theta]);
            %end
            tmpdata = reshape(EEG.data, size(EEG.data,1), size(EEG.data,2)*size(EEG.data,3) );
            if strcmpi(method, 'invdist'), method = 'v4'; end
            [Xi,Yi,badchansdata(:,t)] = griddata(ygood, xgood , double(tmpdata(datachans,t)'),...
                                                    ybad, xbad, method); % interpolate data
        end
        fprintf('\n');
    end

    tmpdata               = zeros(length(bad_elec), EEG.pnts, EEG.trials);
    tmpdata(origoodchans, :,:) = EEG.data;
    %if input data are epoched reshape badchansdata for Octave compatibility...
    if length(size(tmpdata))==3
        badchansdata = reshape(badchansdata,length(badchans),size(tmpdata,2),size(tmpdata,3));
    end

    if ~isempty(t_range)
        t_start = t_range(1);
        t_end   = t_range(2);
        if t_start~=ORIEEG.xmin || t_end~=ORIEEG.xmax
            if length(size(tmpdata))==2
                times_2b_ignored = [1:floor(t_start*ORIEEG.srate), floor(t_end*ORIEEG.srate):floor(ORIEEG.xmax*ORIEEG.srate)];
                badchansdata(:, times_2b_ignored) = ORIEEG.data(badchans, times_2b_ignored);
            end
        end
    end

    tmpdata(badchans,:,:) = badchansdata;
    EEG.data = tmpdata;
    EEG.nbchan = size(EEG.data,1);
    if ~strcmpi(method, 'sphericalfast')
        EEG = eeg_checkset(EEG);
    end

% get data channels
% -----------------
function datachans = getdatachans(goodchans, badchans);
      datachans = goodchans;
      badchans  = sort(badchans);
      for index = length(badchans):-1:1
          datachans(find(datachans > badchans(index))) = datachans(find(datachans > badchans(index)))-1;
      end

% -----------------
% spherical splines
% -----------------
function [x, y, z, Res] = spheric_spline_old( xelec, yelec, zelec, values);

SPHERERES = 20;
[x,y,z] = sphere(SPHERERES);
x(1:(length(x)-1)/2,:) = []; x = [ x(:)' ];
y(1:(length(y)-1)/2,:) = []; y = [ y(:)' ];
z(1:(length(z)-1)/2,:) = []; z = [ z(:)' ];

Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec);
Gsph  = computeg(x,y,z,xelec,yelec,zelec);

% equations are
% Gelec*C + C0  = Potential (C unknown)
% Sum(c_i) = 0
% so
%             [c_1]
%      *      [c_2]
%             [c_ ]
%    xelec    [c_n]
% [x x x x x]         [potential_1]
% [x x x x x]         [potential_ ]
% [x x x x x]       = [potential_ ]
% [x x x x x]         [potential_4]
% [1 1 1 1 1]         [0]

% compute solution for parameters C
% ---------------------------------
meanvalues = mean(values);
values = values - meanvalues; % make mean zero
C = pinv([Gelec;ones(1,length(Gelec))]) * [values(:);0];

% apply results
% -------------
Res = zeros(1,size(Gsph,1));
for j = 1:size(Gsph,1)
    Res(j) = sum(C .* Gsph(j,:)');
end
Res = Res + meanvalues;
Res = reshape(Res, length(x(:)),1);

function [xbad, ybad, zbad, allres] = spheric_spline( xelec, yelec, zelec, xbad, ybad, zbad, values);

newchans = length(xbad);
numpoints = size(values,2);

%SPHERERES = 20;
%[x,y,z] = sphere(SPHERERES);
%x(1:(length(x)-1)/2,:) = []; xbad = [ x(:)'];
%y(1:(length(x)-1)/2,:) = []; ybad = [ y(:)'];
%z(1:(length(x)-1)/2,:) = []; zbad = [ z(:)'];

Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec);
Gsph  = computeg(xbad,ybad,zbad,xelec,yelec,zelec);

% compute solution for parameters C
% ---------------------------------
meanvalues = mean(values);
values = values - repmat(meanvalues, [size(values,1) 1]); % make mean zero

values = [values;zeros(1,numpoints)];
C = pinv([Gelec;ones(1,length(Gelec))]) * values;
clear values;
allres = zeros(newchans, numpoints);

% apply results
% -------------
for j = 1:size(Gsph,1)
    allres(j,:) = sum(C .* repmat(Gsph(j,:)', [1 size(C,2)]));
end
allres = allres + repmat(meanvalues, [size(allres,1) 1]);

% compute G function
% ------------------
function g = computeg(x,y,z,xelec,yelec,zelec)

unitmat = ones(length(x(:)),length(xelec));
EI = unitmat - sqrt((repmat(x(:),1,length(xelec)) - repmat(xelec,length(x(:)),1)).^2 +...
                (repmat(y(:),1,length(xelec)) - repmat(yelec,length(x(:)),1)).^2 +...
                (repmat(z(:),1,length(xelec)) - repmat(zelec,length(x(:)),1)).^2);

g = zeros(length(x(:)),length(xelec));
%dsafds
m = 4; % 3 is linear, 4 is best according to Perrin's curve
for n = 1:7
    if ismatlab
        L = legendre(n,EI);
    else % Octave legendre function cannot process 2-D matrices
        for icol = 1:size(EI,2)
            tmpL = legendre(n,EI(:,icol));
            if icol == 1, L = zeros([ size(tmpL) size(EI,2)]); end
            L(:,:,icol) = tmpL;
        end
    end
    g = g + ((2*n+1)/(n^m*(n+1)^m))*squeeze(L(1,:,:));
end
g = g/(4*pi);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function subpak_test06 ( )

%*****************************************************************************80
%
%% TEST06 tests BAR_CHECK, BAR_CODE, BAR_DIGIT_CODE_LEFT and BAR_DIGIT_CODE_RIGHT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  BAR_CHECK checks digits for a barcode;\n' );
  fprintf ( 1, '  BAR_CODE computes the barcode for a string of\n' );
  fprintf ( 1, '  11 digits;\n' );
  fprintf ( 1, '  BAR_DIGIT_CODE_LEFT returns the left digit code.\n' );
  fprintf ( 1, '  BAR_DIGIT_CODE_RIGHT returns the right digit code.\n' );

  for i = 1 : 11
    digit(i) = mod ( i-1, 10 );
  end

  check = bar_check ( digit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The check digit is %d\n', check );

  digit(12) = check;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The left and right digit codes:\n' );
  fprintf ( 1, '\n' );
  for i = 0 : 9
    codel = bar_digit_code_left ( i );
    coder = bar_digit_code_right ( i );
    fprintf ( 1, '  %2d  %s  %s\n', i, codel, coder );
  end

  bar = bar_code ( digit );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Bar code:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    %s\n', bar(1:9) );
  fprintf ( 1, '    %s\n', bar(10:12) );
  fprintf ( 1, '    %s\n', bar(13:19) );
  fprintf ( 1, '    %s\n', bar(20:26) );
  fprintf ( 1, '    %s\n', bar(27:33) );
  fprintf ( 1, '    %s\n', bar(34:40) );
  fprintf ( 1, '    %s\n', bar(41:47) );
  fprintf ( 1, '    %s\n', bar(48:54) );
  fprintf ( 1, '    %s\n', bar(55:59) );
  fprintf ( 1, '    %s\n', bar(60:66) );
  fprintf ( 1, '    %s\n', bar(67:73) );
  fprintf ( 1, '    %s\n', bar(74:80) );
  fprintf ( 1, '    %s\n', bar(81:87) );
  fprintf ( 1, '    %s\n', bar(88:94) );
  fprintf ( 1, '    %s\n', bar(95:101) );
  fprintf ( 1, '    %s\n', bar(102:104) );
  fprintf ( 1, '    %s\n', bar(105:113) );

  return
end


text:
function Y = multiplyOneRegion( psfMatData, X );
%
%           Y = multiplyOneRegion( psfMatData, X );
%
%  This function computes the multiplication of a point spread
%  function matrix times an image vector; that is,
%
%       Y = psfMatrix * X
%
%  Here we assume the dimension of the PSF is essentially the same
%  as the dimension of the image.
%
%  Input:
%   psfMatData  -  complex array containing the matrix data, usually
%                  computed from onePsfMatrix.m
%            X  -  array containing the image to which the psf matrix
%                  is to be multiplied
%
%  Output:
%          Y  -  contains the result after multiplication.
%

%  J. Nagy  1/6/02

%
%  First we determine if any padding is needed.
%
Xpad = padarray(X, size(psfMatData) - size(X), 'post');

%
%  Now we perform the multiplications.
%  Note that this should work for 2D and 3D images.
%

Y = ifftn( psfMatData .* fftn( Xpad ) );

[nx, ny, nz] = size( X );
Y = real( Y(1:nx, 1:ny, 1:nz) );

repo: IRtools


text:
function legendre_exactness ( n, x, w, p_max )

%*****************************************************************************80
%
%% LEGENDRE_EXACTNESS: quadrature rule exactness for Legendre integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points in the rule.
%
%    Input, real X(N), the quadrature points.
%
%    Input, real W(N), the quadrature weights.
%
%    Input, integer P_MAX, the maximum exponent.
%    0 <= P_MAX.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature rule of order %d for the Legendre integral.\n', n );
  fprintf ( 1, '  Degree          Relative Error\n' );
  fprintf ( 1, '\n' );

  for p = 0 : p_max

    e = legendre_monomial_quadrature ( n, x, w, p );

    fprintf ( 1, '  %6d  %24.16f\n', p, e );

  end

  return
end


text:
function sparse_unique_index = sgmga_unique_index ( dim_num, level_weight, ...
  level_max, rule, growth, np, p, tol, point_num, point_total_num )

%*****************************************************************************80
%
%% SGMGA_UNIQUE_INDEX maps nonunique to unique points of an SGMGA grid.
%
%  Discussion:
%
%    The sparse grid usually contains many points that occur in more
%    than one product grid.
%
%    When generating the point locations, it is easy to realize that a point
%    has already been generated.
%
%    But when it's time to compute the weights of the sparse grids, it is
%    necessary to handle situations in which weights corresponding to
%    the same point generated in multiple grids must be collected together.
%
%    This routine generates ALL the points, including their multiplicities,
%    and figures out a mapping from them to the collapsed set of unique points.
%
%    This mapping can then be used during the weight calculation so that
%    a contribution to the weight gets to the right place.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    An Anisotropic Sparse Grid Stochastic Collocation Method for Partial
%    Differential Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2411-2442.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real LEVEL_WEIGHT(DIM_NUM), the anisotropic weights.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Input, integer RULE(DIM_NUM), the rule in each dimension.
%     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
%     2, "F2",  Fejer Type 2, Open Fully Nested.
%     3, "GP",  Gauss Patterson, Open Fully Nested.
%     4, "GL",  Gauss Legendre, Open Weakly Nested.
%     5, "GH",  Gauss Hermite, Open Weakly Nested.
%     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
%     7, "LG",  Gauss Laguerre, Open Non Nested.
%     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
%     9, "GJ",  Gauss Jacobi, Open Non Nested.
%    10, "HGK", Hermite Genz-Keister, Open Fully Nested.
%    11, "UO",  User supplied Open, presumably Non Nested.
%    12, "UC",  User supplied Closed, presumably Non Nested.
%
%    Input, integer GROWTH(DIM_NUM), the growth in each dimension.
%    0, "DF", default growth associated with this quadrature rule;
%    1, "SL", slow linear, L+1;
%    2  "SO", slow linear odd, O=1+2((L+1)/2)
%    3, "ML", moderate linear, 2L+1;
%    4, "SE", slow exponential;
%    5, "ME", moderate exponential;
%    6, "FE", full exponential.
%
%    Input, integer NP(DIM_NUM), the number of parameters used by each rule.
%
%    Input, real P(*), the parameters needed by each rule.
%
%    Input, real TOL, the tolerance for point equality.
%
%    Input, integer POINT_NUM, the number of unique points in
%    the grid.
%
%    Input, integer POINT_TOTAL_NUM, the total number of points
%    in the grid.
%
%    Output, integer SPARSE_UNIQUE_INDEX(POINT_TOTAL_NUM), lists,
%    for each (nonunique) point, the corresponding index of the same point in
%    the unique listing.
%

%
%  Special cases.
%
  if ( level_max < 0 )
    sparse_unique_index = [];
    return
  end

  if ( level_max == 0 )
    sparse_unique_index(1) = 1;
    return
  end
%
%  Generate SPARSE_TOTAL_ORDER and SPARSE_TOTAL_INDEX arrays
%  for the TOTAL set of points.
%
  sparse_total_order = zeros(dim_num,point_total_num);
  sparse_total_index = zeros(dim_num,point_total_num);

  point_total_num2 = 0;
%
%  Initialization for SGMGA_VCN_ORDERED.
%
  level_weight_min_pos = r8vec_min_pos ( dim_num, level_weight );
  q_min = level_max * level_weight_min_pos - sum ( level_weight(1:dim_num) );
  q_max = level_max * level_weight_min_pos;
  level_1d_max = zeros(dim_num,1);
  for dim = 1 : dim_num
    if ( 0.0 < level_weight(dim) )
      level_1d_max(dim) = floor ( q_max / level_weight(dim) ) + 1;
      if ( q_max <= ( level_1d_max(dim) - 1 ) * level_weight(dim) )
        level_1d_max(dim) = level_1d_max(dim) - 1;
      end
    else
      level_1d_max(dim) = 0;
    end
  end
  more_grids = 0;
  level_1d = [];
%
%  Seek all vectors LEVEL_1D which satisfy the constraint:
%
%    LEVEL_MAX * LEVEL_WEIGHT_MIN_POS - sum ( LEVEL_WEIGHT )
%      < sum ( 1 <= I <= DIM_NUM ) LEVEL_WEIGHT(I) * LEVEL_1D(I)
%      <= LEVEL_MAX * LEVEL_WEIGHT_MIN_POS.
%
  while ( 1 )

    [ level_1d, more_grids ] = sgmga_vcn_ordered ( dim_num, level_weight, ...
      level_1d_max, level_1d, q_min, q_max, more_grids );

    if ( ~more_grids )
      break
    end
%
%  Compute the combinatorial coefficient.
%
    coef = sgmga_vcn_coef ( dim_num, level_weight, level_1d, q_max );

    if ( coef == 0.0 )
      continue
    end
%
%  Transform each 1D level to a corresponding 1D order.
%
    order_1d = level_growth_to_order ( dim_num, level_1d, rule, growth );
%
%  The inner loop generates a POINT of the GRID of the LEVEL.
%
    point_index = [];
    more_points = 0;

    while ( 1 )

      [ point_index, more_points ] = vec_colex_next3 ( dim_num, order_1d, ...
        point_index, more_points );

      if ( ~more_points )
        break;
      end

      point_total_num2 = point_total_num2 + 1;
      sparse_total_order(1:dim_num,point_total_num2) = order_1d(1:dim_num);
      sparse_total_index(1:dim_num,point_total_num2) = point_index(1:dim_num);

    end

  end
%
%  Now compute the coordinates of the TOTAL set of points.
%
  sparse_total_point = zeros(dim_num,point_total_num);
  sparse_total_point(1:dim_num,1:point_total_num) = r8_huge ( );

  level_weight_min_pos = r8vec_min_pos ( dim_num, level_weight );
  q_max = level_max * level_weight_min_pos;

  p_index = 1;

  for dim = 1 : dim_num

    if ( 0 < level_weight(dim) )
      level_1d_max(dim) = floor ( q_max / level_weight(dim) ) + 1;
      if ( q_max <= ( level_1d_max(dim) - 1 ) * level_weight(dim) )
        level_1d_max(dim) = level_1d_max(dim) - 1;
      end
    else
      level_1d_max(dim) = 0;
    end

    for level = 0 : level_1d_max(dim)

      order = level_growth_to_order  ( 1, level, rule(dim), growth(dim) );

      if ( rule(dim) == 1 )
        points = clenshaw_curtis_compute_points ( order );
      elseif ( rule(dim) == 2 )
        points = fejer2_compute_points ( order );
      elseif ( rule(dim) == 3 )
        points = patterson_lookup_points ( order );
      elseif ( rule(dim) == 4 )
        points = legendre_compute_points ( order );
      elseif ( rule(dim) == 5 )
        points = hermite_compute_points ( order );
      elseif ( rule(dim) == 6 )
        alpha = p(p_index);
        points = gen_hermite_compute_points ( order, alpha );
      elseif ( rule(dim) == 7 )
        points = laguerre_compute_points ( order );
      elseif ( rule(dim) == 8 )
        alpha = p(p_index);
        points = gen_laguerre_compute_points ( order, alpha );
      elseif ( rule(dim) == 9 )
        alpha = p(p_index);
        beta = p(p_index+1);
        points = jacobi_compute_points ( order, alpha, beta );
      elseif ( rule(dim) == 10 )
        points = hermite_genz_keister_lookup_points ( order );
      elseif ( rule(dim) == 11 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SGMGA_UNIQUE_INDEX - Fatal error!\n' );
        fprintf ( 1, '  Do not know how to assign points for rule 11.\n' );
        error ( 'SGMGA_UNIQUE_INDEX - Fatal error!' );
      elseif ( rule(dim) == 12 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SGMGA_UNIQUE_INDEX - Fatal error!\n' );
        fprintf ( 1, '  Do not know how to assign points for rule 12.\n' );
        error ( 'SGMGA_UNIQUE_INDEX - Fatal error!' );
      else
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SGMGA_UNIQUE_INDEX - Fatal error!\n' );
        fprintf ( 1,'  Unexpected value of RULE = %d\n', rule(dim) );
        error ( 'SGMGA_UNIQUE_INDEX - Fatal error!' );
      end

      index = find ( sparse_total_order(dim,1:point_total_num) == order );

      sparse_total_point(dim,index) = points ( sparse_total_index(dim,index) );

    end

    p_index = p_index + np(dim);

  end
%
%  Merge points that are too close.
%
  seed = 123456789;

  [ point_num, undx, sparse_unique_index, seed ] = ...
    point_radial_tol_unique_index ( dim_num, point_total_num, ...
    sparse_total_point, tol, seed );

  for point = 1 : point_total_num
    rep = undx(sparse_unique_index(point));
    if ( point ~= rep )
      sparse_total_point(1:dim_num,point) = sparse_total_point(1:dim_num,rep);
    end
  end
%
%  Construct an index that indicates the "rank" of the unique points.
%
  [ undx, sparse_unique_index ] = point_unique_index ( dim_num, ...
    point_total_num, sparse_total_point, point_num );

  return
end


text:
function [ ] = tapas_huge_compile( )
% Compile integrator for DCM fMRI forward model. This function assumes that
% source code file 'tapas_huge_int_euler.c' is in the same folder as this
% function and places the compiled mex-file in the same folder. Running
% this function requires a C compiler to be installed. For more information
% see:
%
% https://www.mathworks.com/support/requirements/supported-compilers.html
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This 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 non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

%%
% check if mex file already exists
if exist('tapas_huge_int_euler', 'file') == 3
    return
end

% get path to toolbox
tbPath = mfilename('fullpath');
tbPath = fileparts(tbPath);

% build full path to source file
scPath = fullfile(tbPath, 'tapas_huge_int_euler.c');

% check if source file present
if exist(scPath,'file')~=2
    error('TAPAS:HUGE:missingSource',...
        ['Could not locate source file. Make sure ' ...
         'tapas_huge_int_euler.c is in the huge toolbox folder.']);
end

% compile
try
    mex(scPath, '-outdir', tbPath);
catch err
    disp('tapas_huge_compile: Failed to compile mex function.');
    disp('Make sure you have selected a C language compiler for mex.');
    disp('For more information, enter mex -help on the command line.');
    rethrow(err)
end

author: translationalneuromodeling

repo: tapas


text:
function [ dat ] = visual_spectrum( data , varargin )
% Description:
%   It is the feature plots of data. This function can be useful to see an
%   overview of data according to time, frequency, and channel.
%
% Example Code:
%    visuspect = visual_spectrum(SMT , {'Xaxis' , 'Frequency'; 'Yaxis' , 'Channel'});
%    visuspect = visual_spectrum(SMT , {'Xaxis' , 'Time'; 'Yaxis' , 'Channel'; 'Band' ,[8 10]});
%    visuspect = visual_spectrum(SMT , {'Xaxis' , 'Time'; 'Yaxis' , 'Frequency'; 'Channel' ,{'C4'}});
%
% Input:
%   [dat] = visual_spectrum(data , <OPT>)
%   data: Data structrue (ex) Epoched data structure
%
% Option:
%      .Xaxis - selecting the domain what you interested in x axis
%                 (e.g. {'Xaxis' , 'Frequency'},{'Xaxis' , 'time'})
%      .Yaxis - selecting the domain what you interested in y axis
%                 (e.g. {'Yaxis' , 'Channel'},{'Yaxis' , 'Frequency'})
%      .Band -  Selecting the interested frequency band in Time-Channel
%               domian
%                 (e.g. {'Band', [8 10]})
%      .Channel - Selecting the interested channel in Time-Frequency domain
%                 (e.g. {'Channel', {'C4'}})
%
% Return:
%    data:  Epoched data structure
%
% See also:
%    opt_cellToStruct
%
% Reference:
%           G. Schalk, D.J. McFarland, T. Hinterberger, N. Birbaumer, and
%           J. R. Wolpaw,"BCI2000: A General-Purpose Brain-Computer
%           Interface (BCI) System, IEEE Transactions on Biomedical
%           Engineering, Vol. 51, No. 6, 2004, pp.1034-1043.
%
% Ji Hoon, Jeong
% jh_jeong@korea.ac.kr
%

%% Load the data
dat = data;
opt = opt_cellToStruct(varargin{:});
sampleFrequency=dat.fs;

%% Setting the parameter
% We refer to visualization matlab code in BCI 2000 open source toolbox
%
% modelOrder = 18+round(sampleFrequency/100);
modelOrder = 20;
opt.lowPassCutoff=(sampleFrequency/2);
opt.freqBinWidth= 2;
opt.highPassCutoff=-1;
opt.trend=1;
% Setting the parameter about spectral moving (using autoregression)
opt.spectralSize= 500;
opt.spectralStep = 100;
opt.spectralSize = round(opt.spectralSize/1000 * sampleFrequency);
opt.spectralStep = round(opt.spectralStep/1000 * sampleFrequency);
params = [modelOrder, opt.highPassCutoff+opt.freqBinWidth/2, opt.lowPassCutoff, opt.freqBinWidth,round(opt.freqBinWidth/.2), ...
    opt.trend, sampleFrequency];
params(8) = opt.spectralStep;
params(9) = opt.spectralSize/opt.spectralStep;
%%
C1_idx=find(dat.y_dec==1);
C2_idx = find(dat.y_dec==2);
C1=prep_selectTrials(dat,{'Index', C1_idx});
C2=prep_selectTrials(dat,{'Index',C2_idx});
for trial=1:size(C1.x,2)
    plotData1=C1.x(:,trial,:);
    plotData1=squeeze(plotData1);
    %trialspectrum = [frequency by channel by time by trial]
    [trialspectrum(:,:,:,trial), C1_freqBins] = mem( plotData1, params );
    plotData2=C2.x(:,trial,:);
    plotData2=squeeze(plotData2);
    [trialspectrum2(:,:,:,trial), C2_freqBins] = mem( plotData2, params );
end
%%
if strcmp(opt.Xaxis, 'Frequency') == 1 && strcmp(opt.Yaxis, 'Channel') == 1
    for tri = 1: size(C1_idx,2)
        tmp=mean(trialspectrum(:,:,:,tri),3);
        C1Data(:,:,tri)=tmp;
        tmp2=mean(trialspectrum2(:,:,:,tri),3);
        C2Data(:,:,tri)=tmp2;
    end
    C1Data = double(C1Data); C2Data = double(C2Data);
    for ch=1:size(C1Data, 2)
        for samp=1:size(C1Data, 1)
            ressq(samp, ch)=rsqu(C1Data(samp, ch, :), C2Data(samp, ch, :));
        end
    end
    plotData = {mean(C1Data,3),mean(C2Data,3), ressq};
    freqBins = C1_freqBins - opt.freqBinWidth/2;
    freqBins(end+1) = freqBins(end) + diff(freqBins(end-1:end));
    for i = 1:size(plotData,2)
        plotData{i}=cat(2, plotData{i}, zeros(size(plotData{i}, 1), 1));
        plotData{i}=cat(1, plotData{i}, zeros(1, size(plotData{i}, 2)));
        figure(i);
        surf(freqBins, [1:size(plotData{i}, 2)], plotData{i}(1:end,:)'); hold on;
        axis tight; set(gcf,'Renderer','Zbuffer');
        xlim([freqBins(2) freqBins(end)]);
        set(gca,'XTick',[freqBins(2):2: freqBins(end)]);
        set(gca,'YTick',[1:1: size(C1Data,2)]);
        set(gca,'Yticklabel',char(dat.chan));
        view(2);colormap jet;colorbar;
        xlabel('Frequency [Hz]'); ylabel('Channels');
        switch(i)
            case 1
                caxis([0 300]);
                title('Frequency-channel power spectrum per channel in Class 1');
            case 2
                caxis([0 300]);
                title('Frequency-channel power spectrum per channel in Class 2');
            case 3
                title('Frequency-channel power spectrum per channels between classes using r^2 values');
        end
    end
%%
elseif strcmp(opt.Xaxis, 'Time') == 1 && strcmp(opt.Yaxis, 'Channel') == 1
    fre_inx = find(C1_freqBins == opt.Band(2));
    for tri = 1: size(C1_idx,2)
        tmp = trialspectrum(fre_inx,:,:,tri);
        tmp = squeeze(tmp);
        C1Data(:,:,tri)=tmp;
        tmp2 = trialspectrum2(fre_inx,:,:,tri);
        tmp2 = squeeze(tmp2);
        C2Data(:,:,tri)=tmp2;
    end
    res1Data = double(permute(C1Data,[2 1 3]));
    res2Data = double(permute(C2Data,[2 1 3]));
    for ch=1:size(res1Data, 2)
        for samp=1:size(res1Data, 1)
            ressq(samp, ch)=rsqu(res1Data(samp, ch, :), res2Data(samp, ch, :));
        end
    end
    plotData = {mean(C1Data,3),mean(C2Data,3) , ressq'};
    a = dat.ival(1)/100;b = size(tmp,2)-abs(dat.ival(1)/100)-1;
    a = a*100;b = b*100;
    for i = 1:size(plotData,2)
        figure(i);
        surf([a:100:b], [1:size(plotData{i}, 1)], plotData{i}(1:end,:)); hold on;
        axis tight; set(gcf,'Renderer','Zbuffer');
        set(gca,'XTick',[a:500:b]);
        set(gca,'YTick',[1:1: size(C1Data,2)]);
        set(gca,'Yticklabel',char(dat.chan));
        view(2);colormap jet;colorbar;
        xlabel('Time'); ylabel('Channels');
        switch(i)
            case 1
                caxis([0 1000]);
                title('Time-channel power spectrum per channel in Class 1');
            case 2
                caxis([0 1000]);
                title('Time-channel power spectrum per channel in Class 2');
            case 3
                title('Time-channel power spectrum per channels between classes using r^2 values');
        end
    end
%%
elseif strcmp(opt.Xaxis, 'Time') == 1 && strcmp(opt.Yaxis, 'Frequency') == 1
    Chanidx = find(strcmp(dat.chan, opt.Channel)== 1 );
    for tri = 1: size(C1_idx,2)
        tmp = trialspectrum(:,Chanidx,:,tri);
        tmp = squeeze(tmp);
        C1Data(:,:,tri)=tmp;
        tmp2 = trialspectrum2(:,Chanidx,:,tri);
        tmp2 = squeeze(tmp2);
        C2Data(:,:,tri)=tmp2;
    end
    plotData = {mean(C1Data,3),mean(C2Data,3)};
    a = dat.ival(1)/100;b = size(tmp,2)-abs(dat.ival(1)/100)-1;
    a = a*100;b = b*100;
    freqBins = C1_freqBins - opt.freqBinWidth/2;
    freqBins(end+1) = freqBins(end) + diff(freqBins(end-1:end));
    for i = 1:size(plotData,2)
        plotData{i}=cat(1, plotData{i}, zeros(1, size(plotData{i}, 2)));
        figure(i);
        surf([a:100:b]', freqBins', plotData{i}(1:end,:)); hold on;
        axis tight; set(gcf,'Renderer','Zbuffer');
        set(gca,'XTick',[a:500:b]);
        ylim([freqBins(2) freqBins(end)]);
        set(gca,'YTick',[freqBins(2):2: freqBins(end)]);
        view(2);colormap jet;colorbar;caxis([0 1000]);
        xlabel('Time'); ylabel('Frequency');
        switch(i)
            case 1
                title(['Time-frequency power spectrum in channel ', char(opt.Channel),  ' and class 1']);
            case 2
                title(['Time-frequency power spectrum in channel ', char(opt.Channel),  ' and class 2']);
        end
    end
end
end

author: PatternRecognition

repo: OpenBMI


text:
function A = slkfd(K, nums, varargin)
%SLKFD Perform Kernelized Fisher Discriminant Analysis
%
% $ Syntax $
%   - A = slkfd(K, nums, ...)
%
% $ Arguments $
%   - K:        the kernel gram matrix of the training samples
%   - nums:     the numbers of samples in all classes
%   - A:        the projection coefficient matrix
%
% $ Description $
%   - A = slkfd(K, nums, ...) performs Kernerlized Fisher discriminant
%     analysis on the samples X according to the specified properties.
%     \*
%     \t   Table 1.  The properties of Fisher Discriminant Analysis   \\
%     \h     name     &     description                                \\
%           'sol'     &  The cell containing the arguments for solving
%                        the generalized eigen-problem by slsymgeig.
%                        default = {}.                                 \\
%           'dimset'  &  The cell containing the arguments for determining
%                        the output feature dimension. default = {}.
%                        (refer to sldim_by_eigval).                   \\
%           'Sb'      &  The pre-computed kernelized between-class
%                        scattering matrix or the cell containing
%                        the arguments for computing the kernelized
%                        scatter matrix in the form {type, ...},
%                        which is input to slkernelscatter.     \\
%           'Sw'      &  The pre-computed kernelized within-class
%                        scattering matrix or the cell containing
%                        the arguments for computing the scatter
%                        matrix in the form {type, ...}, which is
%                        input to slkernelscatter.     \\
%         'weights'   &  The sample weights. default = [].   \\
%     \*
%
% $ History $
%   - Created by Dahua Lin on May 03, 2006
%

%% parse and verify input arguments

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

% for K
n = size(K, 1);
if ~isequal(size(K), [n, n])
    error('sltoolbox:invaliddims', ...
        'K should be a square matrix');
end

% for nums
nc = length(nums);
if ~isequal(size(nums), [1 nc])
    error('sltoolbox:sizmismatch', ...
        'nums should be a 1 x nc row vector');
end

% for options

opts.sol = {};
opts.dimset = {};
opts.Sb = {'Sb'};
opts.Sw = {'Sw'};
opts.weights = [];
opts = slparseprops(opts, varargin{:});

has_Sb = ~isempty(opts.Sb) && isnumeric(opts.Sb);
has_Sw = ~isempty(opts.Sw) && isnumeric(opts.Sw);
if has_Sb
    Sb = opts.Sb;
    if ~isequal(size(Sb), [n n])
        error('sltoolbox:sizmismatch', ...
            'Sb should be a n x n matrix');
    end
end
if has_Sw
    Sw = opts.Sw;
    if ~isequal(size(Sw), [n n])
        error('sltoolbox:sizmismatch', ...
            'Sw should be a n x n matrix');
    end
end

if ~isempty(opts.weights)
    w = opts.weights;
    if ~isequal(size(w), [1 n])
        error('sltoolbox:sizmismatch', ...
            'The weights should be a 1 x n row vector');
    end
else
    w = [];
end

%% Compute

%% Step 1: Construct the eigen-problem

if ~has_Sb
    Sb = slscatter(K, opts.Sb{:}, 'sweights', w, 'nums', nums);
end

if ~has_Sw
    Sw = slscatter(K, opts.Sw{:}, 'sweights', w, 'nums', nums);
end

%% Step 2: Resolve the eigen-problem

[evs, A] = slsymgeig(Sb, Sw, opts.sol{:});

rk = sldim_by_eigval(evs, opts.dimset{:});
A = A(:, 1:rk);

author: lmthang


text:
function [MFCC] = VAD_MFCC(wave,Fs)

%% Extracts MFCC features from the speech signal

frameL=30;
frameshift=10;

Nfft=1024;

Start=1;
Stop=round(frameL/1000*Fs);
Shift=round(frameshift/1000*Fs);
Ind=1;
Win=blackman(Stop);
Es=[];

wave=filter([1 -0.97],1,wave);

while Stop<length(wave)
    Seg=wave(Start:Stop);
    Seg=Seg.*Win;

    Spec=fft(Seg,Nfft);

    mfcc = hspec2fwcep(abs(Spec(1:end/2+1)), Fs, 12);

    MFCC(Ind,:)=mfcc;

    Start=Start+Shift;
    Stop=Stop+Shift;
    Ind=Ind+1;
end

for i = 1 : 13
    MFCC( :, i ) = MFCC( : , i ) - mean(MFCC( : , i));
end

MFCC=MFCC';

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep