text:
function errorxy(Data,varargin);
%----------------------------------------------------------------------------
% errorxy function           Plot graphs with error bars on both axes.
%                          (This is an improved version of errorxy_o.m).
% Input  : - Data matrix.
%          * an arbitrary pairs of keyword and value parameters, where
%            keyword is one of the followings:
%            (In case that the second parameter is a number then call the
%             old version errorxy_o).
%            'ColX'     - Column number containing X values, default is 1.
%            'ColY'     - Column number containing Y values, default is 2.
%            'ColXe'    - Column number containing error X values.
%                         If one parameter is given then assume upper
%                         and lower errors are the same, if two values
%                         are given then the first is the left and the
%                         second is the right errors,
%                         default is NaN (no X errors).
%            'ColYe'    - Column number containing error Y values.
%                         If one parameter is given then assume right
%                         and left errors are the same, if two values
%                         are given then the first is the lower and the
%                         second is the upper errors,
%                         default is 3.
%            'Marker'   - Marker type, see plot for options, default is 'o'.
%            'MarkSize' - Marker size, default is 6.
%            'EdgeColor'- Marker edge color, default is 'b'.
%            'FaceColor'- Marker face color, default is NaN.
%            'ColorEB'  - Error bars color, default is the same as 'EdgeColor'.
%            'LineEB'   - Error bars line type, default is '-'.
%            'WidthEB'  - Error bars line width, default is 0.5.
%            'EdgeEB'   - Length [ticks] of error bar edge, default is 0.
%            'Con'      - Connect data points {'y' | 'n'}, default is 'n'.
%            'ConColor' - Color of Connecting line, default is the
%                         same as 'EdgeColor'
%            'ConLine'  - Type of connecting line, default is '-'.
%            'ConWidth' - Width of connecting line, default is 0.5.
%            'Hist'     - Plot histogram {'y' | 'n'}, default is 'n'.
%            'HistFaceColor'- Histogram face color, default is the same as 'EdgeColor'.
%            'HistEdgeColor'- Histogram edge color, default is the same as 'EdgeColor'.
%            'ScaleX'   - X scale {'linear' | 'log'}, default is 'linear'.
%            'ScaleY'   - Y scale {'linear' | 'log'}, default is 'linear'.
%            'DirX'     - X direction {'normal' | 'reverse'}, default is 'normal'.
%            'DirY'     - Y direction {'normal' | 'reverse'}, default is 'normal'.
% Plot   : errorxy plot
% Example: errorxy(Data,'ColYe',[5 6],'Marker','^');
% Tested : Matlab 5.3
%     By : Eran O. Ofek               March 2004
%    URL : http://wise-obs.tau.ac.il/~eran/matlab.html
%----------------------------------------------------------------------------
Ndots = 450;    % approximate number of dots in paper

if (nargin>1 & ischar(varargin{1})==0),
   %--- call old erroxy_o ---
   errorxy_o(Data,varargin{:});
else
   %--- new version of errorxy ---

   Nvarg = length(varargin);
   if (0.5.*Nvarg~=floor(0.5.*Nvarg)),
      error('Illegal number of input arguments');
   end

   %--- set default values ---
   ColX         = 1;
   ColY         = 2;
   ColXe        = [NaN, NaN];
   ColYe        = [3, 3];
   Marker       = 'o';
   MarkSize     = 6;
   EdgeColor    = 'b';
   FaceColor    = 'n';
   ColorEB      = NaN;
   LineEB       = '-';
   WidthEB      = 0.5;
   EdgeEB       = 0;
   Con          = 'n';
   ConColor     = NaN;
   ConLine      = '-';
   ConWidth     = 0.5;
   Hist         = 'n';
   HistFaceColor= NaN;
   HistEdgeColor= NaN;
   ScaleX       = 'linear';
   ScaleY       = 'linear';
   DirX         = 'normal';
   DirY         = 'normal';

   for I=1:2:Nvarg-1,
      switch varargin{I}
       case 'ColX'
          ColX        = varargin{I+1};
       case 'ColY'
          ColY        = varargin{I+1};
       case 'ColXe'
          ColXe       = varargin{I+1};
       case 'ColYe'
          ColYe       = varargin{I+1};
       case 'Marker'
          Marker      = varargin{I+1};
       case 'MarkSize'
          MarkSize    = varargin{I+1};
       case 'EdgeColor'
          EdgeColor   = varargin{I+1};
       case 'FaceColor'
          FaceColor   = varargin{I+1};
       case 'ColorEB'
          ColorEB     = varargin{I+1};
       case 'LineEB'
          LineEB      = varargin{I+1};
       case 'WidthEB'
          WidthEB     = varargin{I+1};
       case 'EdgeEB'
          EdgeEB      = varargin{I+1};
       case 'Con'
          Con         = varargin{I+1};
       case 'ConColor'
          ConColor    = varargin{I+1};
       case 'ConLine'
          ConLine    = varargin{I+1};
       case 'ConWidth'
          ConWidth   = varargin{I+1};
       case 'Hist'
          Hist       = varargin{I+1};
       case 'HistFaceColor'
          HistFaceColor  = varargin{I+1};
       case 'HistEdgeColor'
          HistEdgeColor  = varargin{I+1};
       case 'ScaleX'
          ScaleX     = varargin{I+1};
       case 'ScaleY'
          ScaleY     = varargin{I+1};
       case 'DirX'
          DirX       = varargin{I+1};
       case 'DirY'
          DirY       = varargin{I+1};
       otherwise
          error('Unknown keyword option');
      end
   end

   %--- set NaN colors ---
   if (isnan(ColorEB)==1),
      ColorEB = EdgeColor;
   end
   if (isnan(ConColor)==1),
      ConColor = EdgeColor;
   end
   if (isnan(HistFaceColor)==1),
      HistFaceColor = EdgeColor;
   end
   if (isnan(HistEdgeColor)==1),
      HistEdgeColor = EdgeColor;
   end

   if (length(ColXe)==1),
      ColXe = [ColXe, ColXe];
   end
   if (length(ColYe)==1),
      ColYe = [ColYe, ColYe];
   end

   %---------------------
   %--- Hold handling ---
   %---------------------
   NextPlot = get(gcf,'NextPlot');
   hold on;
   box on;

   %-----------------------------------
   %--- Data span for edge plotting ---
   %-----------------------------------
   SpanX = max(Data(:,ColX)) - min(Data(:,ColX));
   SpanY = max(Data(:,ColY)) - min(Data(:,ColY));
   EdgeX = SpanX.*EdgeEB./Ndots;
   EdgeY = SpanY.*EdgeEB./Ndots;

   %----------------------
   %--- plot histogram ---
   %----------------------
   switch Hist
    case 'n'
       % do nothing
    case 'y'
       Hhist = bar(Data(:,ColX),Data(:,ColY));
       set(Hhist,'FaceColor',HistFaceColor,...
                 'EdgeColor',HistEdgeColor);
    otherwise
       error('Unknown Hist option');
   end

   %------------------------
   %--- plot data points ---
   %------------------------
   N = size(Data,1);

   for I=1:1:N,
      H         = plot(Data(I,ColX),Data(I,ColY),'o');
      Hpoint(I) = H;
      set(Hpoint,'Marker',Marker,...
            'MarkerSize',MarkSize,...
            'MarkerEdgeColor',EdgeColor,...
            'MarkerFaceColor',FaceColor);
      %--- plot Y error bars ---
      if (isnan(ColYe(1))==0),
         H        = plot([Data(I,ColX); Data(I,ColX)],...
                         [Data(I,ColY)-Data(I,ColYe(1)); Data(I,ColY)+Data(I,ColYe(2))],...
                         '-');
         Hyerr(I) = H;
         set(Hyerr,'Color',ColorEB,...
     	        'LineStyle',LineEB,...
   	        'LineWidth',WidthEB);

         %--- Error bar edge ---
         if (EdgeEB>0),
            %--- Y error bar lower edge ---
    	 H        = plot([Data(I,ColX)-EdgeX; Data(I,ColX)+EdgeX],...
                            [Data(I,ColY)-Data(I,ColYe(1)); Data(I,ColY)-Data(I,ColYe(1))],...
   		         '-');
            Heyel(I) = H;
   	 set(Heyel,'Color',ColorEB,...
     	           'LineStyle',LineEB,...
   	           'LineWidth',WidthEB);
            %--- Y error bar upper edge ---
    	 H        = plot([Data(I,ColX)-EdgeX; Data(I,ColX)+EdgeX],...
                            [Data(I,ColY)+Data(I,ColYe(2)); Data(I,ColY)+Data(I,ColYe(2))],...
   		         '-');
            Heyeu(I) = H;
   	 set(Heyeu,'Color',ColorEB,...
     	           'LineStyle',LineEB,...
   	           'LineWidth',WidthEB);
         end
      end
      %--- plot X error bars ---
      if (isnan(ColXe(1))==0),
         H        = plot([Data(I,ColX)-Data(I,ColXe(1)); Data(I,ColX)+Data(I,ColXe(2))],...
                         [Data(I,ColY); Data(I,ColY)],...
                         '-');
         Hxerr(I) = H;
         set(Hxerr,'Color',ColorEB,...
     	        'LineStyle',LineEB,...
   	        'LineWidth',WidthEB);

         %--- Error bar edge ---
         if (EdgeEB>0),
            %--- X error bar lower edge ---
    	 H        = plot([Data(I,ColX)-Data(I,ColXe(1)); Data(I,ColX)-Data(I,ColXe(1))],...
                            [Data(I,ColY)-EdgeY; Data(I,ColY)+EdgeY],...
   		         '-');
            Hexel(I) = H;
   	 set(Hexel,'Color',ColorEB,...
     	           'LineStyle',LineEB,...
   	           'LineWidth',WidthEB);
            %--- X error bar upper edge ---
    	 H        = plot([Data(I,ColX)+Data(I,ColXe(2)); Data(I,ColX)+Data(I,ColXe(2))],...
                            [Data(I,ColY)-EdgeY; Data(I,ColY)+EdgeY],...
   		         '-');
            Hexer(I) = H;
   	 set(Hexer,'Color',ColorEB,...
     	           'LineStyle',LineEB,...
   	           'LineWidth',WidthEB);
         end
      end

   end

   %---------------------------
   %--- connect data points ---
   %---------------------------
   switch Con
    case 'n'
       % do nothing
    case 'y'
       Hcon = plot(Data(:,ColX),Data(:,ColY),'-');
       set(Hcon,'Color',ConColor,...
    	     'LineStyle',ConLine,...
      	     'LineWidth',ConWidth);
    otherwise
       error('Unknown Con option');
   end

   %---------------------------
   %--- set scale/direction ---
   %---------------------------
   set(gca,'XScale',ScaleX,...
           'YScale',ScaleY,...
           'XDir',DirX,...
           'YDir',DirY);

   %--- return hold to riginal state ---
   set(gcf,'NextPlot',NextPlot);
end


text:
classdef StressDescriptor < FieldNameDescriptor

    methods (Access = protected)

        function loadFieldVariable(obj)
            obj.fieldName = 'stress';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


% LEARNALGOKISSME Keep it simple and straightforward metric learning!
%
% LearnAlgoKISSME properties:
%   p  - Parameter struct.
%   p.roccolor - default 'g'.
%   p.pmetric  - Backproject on the cone of p.s.d. matrices? Default 1.
%   p.lambda   - Influence of the dissimilar pairs. Range 0 to 1, default
%   1.
%   s - struct to store the last result
%   available - if the LearnAlgoKISSME is ready
%
% LearnAlgoKISSME properties (Constant):
%   type - string identifier default 'kissme'
%
% LearnAlgoKISSME methods:
%   LearnAlgoKISSME  - ctor.
%   learnPairwise    - learn from pairwise equivalence labels
%   learn            - learn with fully supervised data
%   dist             - calculate distance between pairs of samples
%
% See also LearnAlgoLMNN,LearnAlgoITML,LearnAlgoLDML
%
% copyright by Martin Koestinger (2011)
% Graz University of Technology
% contact koestinger@icg.tugraz.at
%
% For more information, see the ICG Web site.
%
classdef LearnAlgoKISSME < LearnAlgo

    properties
       p %parameters
       s %struct
       available %are we ready?
    end

    properties (Constant)
        type = 'kissme'
    end

    methods
        function obj = LearnAlgoKISSME(p)
            % obj = LearnAlgoKISSME(p)
            % Constructor
            %
            % parameters:
            %   p  - Parameter struct.
            %   p.roccolor - default 'g'.
            %   p.pmetric - Backproject on the cone of p.s.d. matrices? Default 1.
            %   p.lambda - Influence of the dissimilar pairs. Range 0 to 1, default 1.
            %
            % return:
            %   obj       - object handle to instance of class LearnAlgoKISSME
            %
           if nargin < 1
               p = struct();
           end

           if ~isfield(p,'lambda')
             p.lambda =  1;
           end

           if ~isfield(p,'roccolor')
                p.roccolor = 'g';
           end

           if ~isfield(p,'pmetric')
               p.pmetric = 1;
           end

           obj.p  = p;
           check(obj);
        end

        function bool = check(obj)
           % bool = check(obj)
           % Checks if all dependencies are satisfied
           %
           bool = exist('SOPD') ~= 0;
           if ~bool
               className = class(obj);
               fprintf('Sorry %s not available\n',className);
           end
           obj.available = bool;
        end

        function s = learnPairwise(obj,X,idxa,idxb,matches)
            % s = learnPairwise(obj,X,idxa,idxb,matches)
            % Learn from pairwise equivalence labels
            %
            % parameters:
            %   obj       - instance of class LearnAlgoKISSME
            %   X         - input matrix, each column is an input vector
            %   [DxN*2]. N is the number of pairs. D is the feature
            %   dimensionality
            %   idxa      - index of image A in X [1xN]
            %   idxb      - index of image B in X [1xN]
            %   matches   - matches defines if a pair is similar (1) or
            %   dissimilar (0)
            %
            % return:
            %   s         - Result data struct
            %   s.M       - Trained quadratic distance metric
            %   s.t       - Training time in seconds
            %   s.p       - Used parameters, see LearnAlgoKISSME properties for details.s
            %   s.learnAlgo - class handle to obj
            %   s.roccolor  - line color for ROC curve, default 'g'
            %
            if ~obj.available
                s = struct();
                return;
            end

            idxa = double(idxa);
            idxb = double(idxb);

            %--------------------------------------------------------------
            %   KISS Metric Learning CORE ALGORITHM
            %

            tic;
            % Eqn. (12) - sum of outer products of pairwise differences (similar pairs)
            % normalized by the number of similar pairs.
            covMatches    = SOPD(X,idxa(matches),idxb(matches)) / sum(matches);
            % Eqn. (13) - sum of outer products of pairwise differences (dissimilar pairs)
            % normalized by the number of dissimilar pairs.
            covNonMatches = SOPD(X,idxa(~matches),idxb(~matches)) / sum(~matches);
            t = toc;

            tic;
            % Eqn. (15-16)
            s.M = inv(covMatches) - obj.p.lambda * inv(covNonMatches);
            if obj.p.pmetric
                % to induce a valid pseudo metric we enforce that  M is p.s.d.
                % by clipping the spectrum
                s.M = validateCovMatrix(s.M);
            end
            s.t = toc + t;

            %
            %   END KISS Metric Learning CORE ALGORITHM
            %--------------------------------------------------------------

            s.learnAlgo = obj;
            s.roccolor = obj.p.roccolor;
        end

        function s = learn(obj,X,y)
            % not implemented yet, sorry.
            if ~obj.available
                s = struct();
                return;
            end

            s.roccolor = obj.p.roccolor;
            error('not implemented yet!');
        end

        function d = dist(obj, s, X, idxa,idxb)
            % d = dist(obj, s, X, idxa,idxb)
            % Calculate the distance between pairs of samples
            %
            % parameters:
            %   obj       - instance of class LearnAlgoKISSME
            %   s         - struct with member M, quadratic form
            %   X         - Input matrix (each column is an input vector)
            %   idxa,idxb - idxa(c),idxb(c) index of images a,b, pair c in X
            %   matches   - matches(c) defines if pair c is similar (1) or dissimilar (0)
            %
            % return:
            %   d         - distance for each pair specified in idxa,idxb
            %
            d = cdistM(s.M,X,idxa,idxb);
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Eugenio Alcala Baselga
% Date: 02/06/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function errors_pos_ex = errors_Lyapunov(input)

x_d = input(1);
y_d = input(2);
phi_d = input(3);

x_r = input(4);
y_r = input(5);
phi_r = input(6);

% Lyapunov error computation
e_x = (x_d-x_r)*cos(phi_d)  + (y_d-y_r)*sin(phi_d);
e_y = -(x_d-x_r)*sin(phi_d) + (y_d-y_r)*cos(phi_d);
e_phi = limitare_unghi(phi_d - phi_r);

% SMC error computation
% Lh = -2.5;
% e_x = (x_d-x_r + Lh*cos(phi_r))*cos(phi_d) + (y_d-y_r + Lh*sin(phi_r))*sin(phi_d);
% e_y = -(x_d-x_r + Lh*cos(phi_r))*sin(phi_d) + (y_d-y_r + Lh*sin(phi_r))*cos(phi_d);
% e_phi = limitare_unghi(phi_d - phi_r);

errors_pos_ex = [e_x; e_y; e_phi];


text:
clear;
clc;

s = iio_sys_obj_matlab; % Constructor
s.ip_address = '10.66.99.200';
s.dev_name = 'ad9361';
s.in_ch_no = 2;
s.out_ch_no = 2;
s.in_ch_size = 8192;
s.out_ch_size = 8192;

s = s.setupImpl();

input = cell(1, s.in_ch_no + length(s.iio_dev_cfg.cfg_ch));
Fs = 30.72e6;
Fc = 1e6;
t = 1/Fs:1/Fs:s.in_ch_size/Fs;
for i=1:s.in_ch_no
    input{i} = sin(2*pi*Fc*t+(i-1)*pi/2)*1024;
end
input{s.getInChannel('RX_LO_FREQ')} = 2.4e9;
input{s.getInChannel('RX_SAMPLING_FREQ')} = 30.72e6;
input{s.getInChannel('RX_RF_BANDWIDTH')} = 18.0e6;
input{s.getInChannel('RX1_GAIN_MODE')} = 'slow_attack';
input{s.getInChannel('RX1_GAIN')} = 0;
input{s.getInChannel('RX2_GAIN_MODE')} = 'slow_attack';
input{s.getInChannel('RX2_GAIN')} = 0;
input{s.getInChannel('TX_LO_FREQ')} = 2.4e9;
input{s.getInChannel('TX_SAMPLING_FREQ')} = 30.72e6;
input{s.getInChannel('TX_RF_BANDWIDTH')} = 18.0e6;

output = cell(1, s.out_ch_no + length(s.iio_dev_cfg.mon_ch));

for i = 1:5
output = stepImpl(s, input);
rssi1 = output{s.out_ch_no+1};
rssi2 = output{s.out_ch_no+2};
end

s.releaseImpl();

figure % new figure
ax1 = subplot(2,1,1); % top subplot
ax2 = subplot(2,1,2); % bottom subplot

plot(ax1,output{1});
title(ax1,'I');
xlabel('Sample');
ylabel('Amplitude');

plot(ax2,output{2});
title(ax2,'Q');
xlabel('Sample');
ylabel('Amplitude');

author: analogdevicesinc


text:
function net = vgg16_intialization_down5(warpsize,cdim,ndim, RoIRoD)

net = load('data/pretrained_models/imagenet-vgg-verydeep-16.mat');
net = vl_simplenn_tidy(net);
net.layers = net.layers(1:31);

% Convert to DagNN.
net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;

% Dimension reduction
net.addLayer('dimred', convBlocknoPad(1,1,512,cdim), {'x31'}, {'xRed'}, {'dimred_filter','dimred_bias'});
net = initLayerParams(net);

% RoIAlign
if strcmp(RoIRoD,'RoIRoD') || strcmp(RoIRoD,'RoIOnly')
    roiSize = [warpsize,warpsize];
    net.addLayer('roialign', dagnn.RoiAlign('subdivisions',roiSize,'transformation',1/32), {'xRed','rois'}, 'xRoi');
    net.addLayer('roipool', dagnn.Pooling('method', 'avg', 'poolSize', [2 2], 'stride', 1), {'xRoi'}, {'pool_xRoi'}, {});
end

% RoDAlign
if strcmp(RoIRoD,'RoIRoD') || strcmp(RoIRoD,'RoDOnly')
    rodSize = [warpsize+1,warpsize+1];
    net.addLayer('rodalign', dagnn.RodAlign('subdivisions',rodSize,'transformation',1/32), {'xRed','rois'}, 'xRod');
    net.addLayer('rodpool', dagnn.Pooling('method', 'avg', 'poolSize', [2 2], 'stride', 1), {'xRod'}, {'pool_xRod'}, {});
end

switch RoIRoD
    case 'RoIRoD'
        net.addLayer('concat', dagnn.Concat('dim', 3), {'pool_xRoi', 'pool_xRod'}, {'xConc'}, {});
        net.addLayer('fc1', convBlocknoPad(warpsize,warpsize,cdim*2,ndim), {'xConc'}, {'fc1'}, {'fc1_filter', 'fc1_bias'});
    case 'RoIOnly'
        net.addLayer('fc1', convBlocknoPad(warpsize,warpsize,cdim,ndim), {'pool_xRoi'}, {'fc1'}, {'fc1_filter', 'fc1_bias'});
    case 'RoDOnly'
        net.addLayer('fc1', convBlocknoPad(warpsize,warpsize,cdim,ndim), {'pool_xRod'}, {'fc1'}, {'fc1_filter', 'fc1_bias'});
    otherwise
        error('Unknown RoIRoD choice.');
end
net = initLayerParams(net);
net.addLayer('relu_fc1', dagnn.ReLU(), {'fc1'}, {'relu_fc1'}, {});
net.addLayer('fc2', convBlocknoPad(1,1,ndim,ndim), {'relu_fc1'}, {'fc2'}, {'fc2_filter', 'fc2_bias'});
net = initLayerParams(net);
net.addLayer('relu_fc2', dagnn.ReLU(), {'fc2'}, {'relu_fc2'}, {});

net.addLayer('dropout',dagnn.DropOut('rate',0.5),'relu_fc2','dropout',{});

net.addLayer('predcls',convBlocknoPad(1,1,ndim,1), 'dropout','predcls',{'predcls_filter','predcls_bias'});
net = initLayerParams(net);

% net.print({'input', [256 256 3]}, 'all', true);
net.addLayer('losscls',dagnn.RegressionLoss('lossType','Huber'), {'predcls','label'}, 'losscls',{});

% Meta parameters
net.meta.trainOpts.learningRate = [logspace(-5,-4,5) logspace(-4,-4,35)];
net.meta.trainOpts.weightDecay = 0.0001;
net.meta.trainOpts.batchSize = 1 ;
net.meta.trainOpts.solver = @solver.adam ;
net.meta.trainOpts.numEpochs = numel(net.meta.trainOpts.learningRate) ;
net.meta.normalization.averageImage = reshape([122.7717 102.9801 115.9465],[1 1 3]);

end

function convObj = convBlocknoPad(fh,fw,fc,k)
    convObj = dagnn.Conv('size', [fh fw fc k], 'hasBias', true);
end

function net = initLayerParams(net)
    p = net.getParamIndex(net.layers(end).params) ;
    params = net.layers(end).block.initParams();
    [net.params(p).value] = deal(params{:}) ;
end

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
%% RATE OF CONVERGENCE OF CUBIC FINITE ELEMENT METHOD
%
% This example is to show the rate of convergence of cubic finite element
% approximation of the Poisson equation on the regular polygons:
%
% $$- \Delta u = f \; \hbox{in } Omega$$
%
% for the Dirichlet boundary condition.

clear variables
%% Set up problem
% PDE and Boundary condition.
pde = simpledata; % f = 1, g_D = 0
% FEM
option.elemType = 'P3';

%% Options
option.L0 = 0;
option.maxIt = 4;
option.printlevel = 1;
option.plotflag = 1;

%% Case 1: Triangle
[node,elem] = regpolygon(3,0.5);
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
showmesh(node,elem);
femPoisson(mesh,pde,option);

%% Case 2: Square
[node,elem] = squaremesh([0,1,0,1],0.25);
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
showmesh(node,elem);
femPoisson(mesh,pde,option);

%% Case 3: Pentagon
[node,elem] = regpolygon(5,0.5);
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
showmesh(node,elem);
femPoisson(mesh,pde,option);

%% Case 4: Hexagon
[node,elem] = regpolygon(6,0.5);
showmesh(node,elem);
bdFlag = setboundary(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
showmesh(node,elem);
femPoisson(mesh,pde,option);

repo: ifem


text:
% UNDER construction...
%
%
%
% Usage
% >> waverage = stdbin(ERP, bop)
%
% For Bin Operations GUI/scripting purposes
%
% Author: Javier
%

function  standard = stdbin(ERP, bop)

if nargin<1
        help stdbin
        return
end
if nargin<2
        error('ERPLAB says: You must specify bin indexes to be processed!')
end
if ~iserpstruct(ERP)
        error('ERPLAB says: Your data structure is not an ERP structure!')
end
if ischar(bop)
        bop = str2num(char(regexp(bop,'\d+','match')'))';
end
binarray = unique_bc2(bop);
if length(binarray)~=length(bop)
        fprintf('\n*** WARNING: Repeated bins were ignored.\n\n')
end
if length(binarray)<2
        error('ERPLAB says: You must specify 2 bin indexes at leat!')
end
if max(binarray)>ERP.nbin
        error('ERPLAB says: Some specified bins do not exist!')
end
datavg   = ERP.bindata(:,:,binarray);
standard = std(datavg, 0, 3);

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function checkCoordsNodes(flat)
%
%AUTHOR:  Wandell & Huk
%DATE:    12.01.00
%PURPOSE:
%   We verify that the gray matter nodes in the segmentation
% is consistent with the coords and nodes in the mrLoadRet
% directories.
%
% The nodes in Gray/coords are leftNodes and rightNodes.  These
% are supposed to be a subset of the nodes in the segmentation
% at leftPath and rightPath.  So, we load the two sets and we
% make sure that every one of the leftNodes is inside the set of
% nodes read from the segmentation.  Then we do it again for the
% right.
%
% bw, 01/12/01  Updated Flat coords reading so that it checks
%   the proper Flat sub-directory.
% djh, 2/5/2001
%    Use fullfile instead of hardcoded 'load Gray/coords'
%    Added check to compare mrLoadRet's gray nodes with the full set of gray nodes

mrGlobals

% open hiddenGray to get the relevant info
gray = getSelectedGray;
if isempty(gray)
    gray = initHiddenGray;
end
if isempty(flat)
    flat = initHiddenFlat;
end

% yank out relevant rows and reorder them as necessary
grayNodes = gray.nodes(1:3,:);
if ~isempty(gray.allLeftNodes)
    grayAllLeftNodes = gray.allLeftNodes(1:3,:);
else
    grayAllLeftNodes = [];
end
if ~isempty(gray.allRightNodes)
    grayAllRightNodes = gray.allRightNodes(1:3,:);
else
    grayAllRightNodes = [];
end
if ~isempty(flat.grayCoords{1})
    flatLeftGrayCoords = flat.grayCoords{1}([2 1 3],:);
else
    flatLeftGrayCoords = [];
end
if ~isempty(flat.grayCoords{2})
    flatRightGrayCoords = flat.grayCoords{2}([2 1 3],:);
else
    flatRightGrayCoords = [];
end

% check that gray.nodes are a subset of the full set of gray nodes
if isempty(gray.nodes)
    disp('Empty gray.nodes.');
else
    [cLeft,iaLeft,ibLeft] = intersectCols(grayNodes,grayAllLeftNodes);
    [cRight,iaRight,ibRight] = intersectCols(grayNodes,grayAllRightNodes);
    numDifference = size(grayNodes,2) - (size(cLeft,2)+size(cRight,2));
    if numDifference~=0
        disp('Gray.nodes are NOT all contained in the segmentation.');
        disp(['Number missing:  ',num2str(numDifference)]);
        mismatch(gray,flat);
    else
        disp('Gray.nodes are all contained in the segmentation.');
    end
end

% check that flat.grayCoords are contained in the gray.nodes
if isempty(flat.grayCoords)
    disp('Empty flat.grayCoords.');
else
    [cLeft,iaLeft,ibLeft] = intersectCols(flatLeftGrayCoords,grayNodes);
    [cRight,iaRight,ibRight] = intersectCols(flatRightGrayCoords,grayNodes);
    numDifference = (size(flatLeftGrayCoords,2) - size(cLeft,2)) + ...
                    (size(flatRightGrayCoords,2) - size(cRight,2));
    if numDifference~=0
        disp('Flat.grayCoords are NOT all contained in the gray.nodes');
        disp(['Number missing: ',num2str(numDifference)]);
        mismatch(gray,flat);
    else
        disp('Flat.grayCoords are all contained in the gray.nodes');
    end
end

% check that flat nodes are contained in segmentation
if isempty(flat.grayCoords)
    disp('Empty flat.grayCoords.');
else
    [cLeft,iaLeft,ibLeft] = intersectCols(flatLeftGrayCoords,grayAllLeftNodes);
    [cRight,iaRight,ibRight] = intersectCols(flatRightGrayCoords,grayAllRightNodes);
    numDifference = (size(flatLeftGrayCoords,2) - size(cLeft,2)) + ...
                    (size(flatRightGrayCoords,2) - size(cRight,2));
    if numDifference~=0
        disp('flat.grayCoords are NOT all contained in the segmentation.');
        disp(['Number missing:  ',num2str(numDifference)]);
        mismatch(gray,flat);
    else
        disp('flat.grayCoords are all contained in the segmentation.');
    end
end

msgbox([viewDir(flat),' and ',viewDir(gray),' are OK.'],'Verify Gray-Flat match');

return;

function mismatch(gray,flat)
myErrorDlg(['Mismatch detected between ',viewDir(flat),' and ',viewDir(gray),...
        '. Run Install Segmentation and/or Install Unfold from Xform menu.']);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [cepstra,aspectrum,pspectrum] = melfcc(samples, sr, varargin)
%[cepstra,aspectrum,pspectrum] = melfcc(samples, sr[, opts ...])
%  Calculate Mel-frequency cepstral coefficients by:
%   - take the absolute value of the STFT
%   - warp to a Mel frequency scale
%   - take the DCT of the log-Mel-spectrum
%   - return the first <ncep> components
%  This version allows a lot of options to be controlled, as optional
%  'name', value pairs from the 3rd argument on: (defaults in parens)
%    'wintime' (0.025): window length in sec
%    'hoptime' (0.010): step between successive windows in sec
%    'numcep'     (13): number of cepstra to return
%    'lifterexp' (0.6): exponent for liftering; 0 = none; < 0 = HTK sin lifter
%    'sumpower'    (1): 1 = sum abs(fft)^2; 0 = sum abs(fft)
%    'preemph'  (0.97): apply pre-emphasis filter [1 -preemph] (0 = none)
%    'dither'      (0): 1 = add offset to spectrum as if dither noise
%    'minfreq'     (0): lowest band edge of mel filters (Hz)
%    'maxfreq'  (4000): highest band edge of mel filters (Hz)
%    'nbands'     (40): number of warped spectral bands to use
%    'bwidth'    (1.0): width of aud spec filters relative to default
%    'dcttype'     (2): type of DCT used - 1 or 2 (or 3 for HTK or 4 for feac)
%    'fbtype'  ('mel'): frequency warp: 'mel','bark','htkmel','fcmel'
%    'usecmp'      (0): apply equal-loudness weighting and cube-root compr.
%    'modelorder'  (0): if > 0, fit a PLP model of this order
%    'broaden'     (0): flag to retain the (useless?) first and last bands
%    'useenergy'   (0): overwrite C0 with true log energy
% The following non-default values nearly duplicate Malcolm Slaney's mfcc
% (i.e. melfcc(d,16000,opts...) =~= log(10)*2*mfcc(d*(2^17),16000) )
%       'wintime': 0.016
%     'lifterexp': 0
%       'minfreq': 133.33
%       'maxfreq': 6855.6
%      'sumpower': 0
% The following non-default values nearly duplicate HTK's MFCC
% (i.e. melfcc(d,16000,opts...) =~= 2*htkmelfcc(:,[13,[1:12]])'
%  where HTK config has PREEMCOEF = 0.97, NUMCHANS = 20, CEPLIFTER = 22,
%  NUMCEPS = 12, WINDOWSIZE = 250000.0, USEHAMMING = T, TARGETKIND = MFCC_0)
%     'lifterexp': -22
%        'nbands': 20
%       'maxfreq': 8000
%      'sumpower': 0
%        'fbtype': 'htkmel'
%       'dcttype': 3
% For more detail on reproducing other programs' outputs, see
% http://www.ee.columbia.edu/~dpwe/resources/matlab/rastamat/mfccs.html
%
% 2005-04-19 dpwe@ee.columbia.edu after rastaplp.m.
% Uses Mark Paskin's process_options.m from KPMtools
% $Header: /Users/dpwe/matlab/rastamat/RCS/melfcc.m,v 1.3 2012/09/03 14:01:26 dpwe Exp dpwe $

if nargin < 2;   sr = 16000;    end

% Parse out the optional arguments
[wintime, hoptime, numcep, lifterexp, sumpower, preemph, dither, ...
 minfreq, maxfreq, nbands, bwidth, dcttype, fbtype, usecmp, modelorder, ...
 broaden, useenergy] = ...
    process_options(varargin, 'wintime', 0.025, 'hoptime', 0.010, ...
          'numcep', 13, 'lifterexp', 0.6, 'sumpower', 1, 'preemph', 0.97, ...
	  'dither', 0, 'minfreq', 0, 'maxfreq', 4000, ...
	  'nbands', 40, 'bwidth', 1.0, 'dcttype', 2, ...
	  'fbtype', 'mel', 'usecmp', 0, 'modelorder', 0, ...
          'broaden', 0, 'useenergy', 0);

if preemph ~= 0
  samples = filter([1 -preemph], 1, samples);
end

% Compute FFT power spectrum
[pspectrum,logE] = powspec(samples, sr, wintime, hoptime, dither);

aspectrum = audspec(pspectrum, sr, nbands, fbtype, minfreq, maxfreq, sumpower, bwidth);

if (usecmp)
  % PLP-like weighting/compression
  aspectrum = postaud(aspectrum, maxfreq, fbtype, broaden);
end

if modelorder > 0

  if (dcttype ~= 1)
    disp(['warning: plp cepstra are implicitly dcttype 1 (not ', num2str(dcttype), ')']);
  end

  % LPC analysis
  lpcas = dolpc(aspectrum, modelorder);

  % convert lpc to cepstra
  cepstra = lpc2cep(lpcas, numcep);

  % Return the auditory spectrum corresponding to the cepstra?
%  aspectrum = lpc2spec(lpcas, nbands);
  % else return the aspectrum that the cepstra are based on, prior to PLP

else

  % Convert to cepstra via DCT
  cepstra = spec2cep(aspectrum, numcep, dcttype);

end

cepstra = lifter(cepstra, lifterexp);

if useenergy
  cepstra(1,:) = logE;
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [Adraw]=sampleA_H(yData,Psi,B,Hvars,T,priorValues,dataValues)
% This function takes a draw from the conditional posterior distribution of
% A. The algorithm is based on Cogley and Sargent (2005).

%% Initialize
M=size(yData,2);
p = size(B,1)/M;
Hscaling=Hvars(p+1:T,:).^0.5; %sqrt of the diagonal elements of the H matrix

% priors
priorVarAscaling=priorValues.priorVarAscaling_H;  %prior variance for the below diagonal elements of A
priorMeanAscaling=priorValues.priorMeanAscaling_H; %prior mean for the below diagonal elements of A

%% Prepare data
Y_Psi=yData-Psi; %remove local mean
Y_Psi(1:p,:)=yData(1:p,:)-ones(p,1)*mean(yData(1:p,:));
X_Psi = lagmatrix(Y_Psi,1:p);
X_Psi = X_Psi(p+1:end,:);
Y_Psi=Y_Psi(p+1:end,:);

V=Y_Psi-X_Psi*B; %generate residuals

%% Sample row-by-row
Atemp=eye(M);

for i=2:M %the first row of A is by construction [1 0 ...M]

    Hscaling_i=Hscaling(:,i)*ones(1,M);
    Vscaled_i=V./Hscaling_i;

    % row data
    zi=Vscaled_i(:,i);  %rescaled residuals
    Zi=-Vscaled_i(:,1:i-1);

    % make row priors
    priorVar_ai=priorVarAscaling*eye(i-1);
    priorMean_ai=priorMeanAscaling*ones(i-1,1);

    % determine posterior
    postVar_ai=(Zi'*Zi+(priorVar_ai\eye(i-1)))\eye(i-1); %posterior variance implied by the normality assumption on prior and likelihood
    postMean_ai=postVar_ai*(Zi'*zi+(priorVar_ai\eye(i-1))*priorMean_ai); %posterior mean implied by the normality assumption

    % obtain cholvar
    [Cvar,testPD]=chol(postVar_ai);
    if testPD>0
        Cvar= cholred(postVar_ai);
        disp('NPD!')
    end
    Cvar=Cvar'; % transpose to obtain lower triangular matrix

    % sample from posterior
    aiDraw=postMean_ai+Cvar*randn(i-1,1);

    % save aiDraw
    Atemp(i,1:i-1)=aiDraw';

end

%% Save the draw of A
Adraw=Atemp;

end

author: european-central-bank

repo: BEAR-toolbox


text:
function v = cvx_value( x )
v = value( x );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
classdef TestSeamlessClone
    %TestSeamlessClone

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','balloon.jpg'));
            dst = img;
            mask = zeros(size(img,1), size(img,2), 'uint8');
            mask(20:120, 5:100) = 255;
            p = [400 200];

            out = cv.seamlessClone(img, dst, mask, p, 'Method','NormalClone');
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.seamlessClone();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


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

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

% Last Modified by GUIDE v2.5 11-Jan-2021 18:38:54

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

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

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

    % Choose default command line output for lvd_EditGeometricRefFramesGUI
    handles.output = hObject;

    lvdData = varargin{1};
    setappdata(hObject,'lvdData',lvdData);

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, lvdData)
    set(handles.refFrameListBox,'String',lvdData.geometry.refFrames.getListboxStr());

    numRefFrames = length(get(handles.refFrameListBox,'String'));

    if(handles.refFrameListBox.Value < 1)
        handles.refFrameListBox.Value = 1;
    elseif(handles.refFrameListBox.Value > numRefFrames)
        handles.refFrameListBox.Value = numRefFrames;
    end

    if(numRefFrames <= 0)
        handles.removeRefFrameButton.Enable = 'off';
    else
        handles.removeRefFrameButton.Enable = 'on';
    end

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

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        close(handles.lvd_EditGeometricRefFramesGUI);
    end

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

% --- Executes on selection change in refFrameListBox.
function refFrameListBox_Callback(hObject, eventdata, handles)
% hObject    handle to refFrameListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns refFrameListBox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from refFrameListBox
    lvdData = getappdata(handles.lvd_EditGeometricRefFramesGUI,'lvdData');
    if(strcmpi(get(handles.lvd_EditGeometricRefFramesGUI,'SelectionType'),'open') && ...
       lvdData.geometry.refFrames.getNumRefFrames() > 0)
        selRefFrameInd = hObject.Value;
        selRefFrame = lvdData.geometry.refFrames.getRefFrameAtInd(selRefFrameInd);

        selRefFrame.openEditDialog();

        populateGUI(handles, lvdData);
    end

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

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in addRefFrameButton.
function addRefFrameButton_Callback(hObject, eventdata, handles)
% hObject    handle to addRefFrameButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditGeometricRefFramesGUI,'lvdData');

    listboxStr = GeometricRefFrameEnum.getListBoxStr();
    [Selection,ok] = listdlgARH('ListString', listboxStr, ...
                                'SelectionMode', 'single', ...
                                'ListSize', [300 300], ...
                                'InitialValue',1, ...
                                'Name', 'Select New Reference Frame Type', ...
                                'PromptString', 'Select the type of the reference frame to be created:');

    if(ok)
        enum = GeometricRefFrameEnum.getEnumForListboxStr(listboxStr{Selection});

        switch enum
            case GeometricRefFrameEnum.CoordSystemOrigin
                numPoints = lvdData.geometry.points.getNumPoints();
                numCoordSystems = lvdData.geometry.coordSyses.getNumCoordSyses();
                if(numPoints >= 1 && numCoordSystems >= 1)
                    origin = lvdData.geometry.points.getPointAtInd(1);
                    coordSys = lvdData.geometry.coordSyses.getCoordSysAtInd(1);
                    newRefFrame = CoordSysPointRefFrame(coordSys, origin, 'New Reference Frame', lvdData);

                else
                    errordlg('Cannot create Coordinate System at Origin Reference Frame: There must be at least one geometric point and one coordinate system available.','Cannot Create Reference Frame');
                    return;
                end

            otherwise
                error('Unknown Reference Frame Type: %s', enum.name)
        end

        useTF = newRefFrame.openEditDialog();

        if(useTF)
            lvdData.geometry.refFrames.addRefFrame(newRefFrame);

            populateGUI(handles, lvdData);
        end
    end

% --- Executes on button press in removeRefFrameButton.
function removeRefFrameButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeRefFrameButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditGeometricRefFramesGUI,'lvdData');

    selRefFrameInd = handles.refFrameListBox.Value;
    selRefFrame = lvdData.geometry.refFrames.getRefFrameAtInd(selRefFrameInd);

	tf = selRefFrame.isInUse(lvdData);

    if(tf == false)
        lvdData.geometry.refFrames.removeRefFrame(selRefFrame);

        lvdData.viewSettings.removeGeoRefFrameFromViewProfiles(selRefFrame);

        populateGUI(handles, lvdData);

        numRefFrames = lvdData.geometry.refFrames.getNumRefFrames();
        if(handles.refFrameListBox.Value > numRefFrames)
            set(handles.refFrameListBox,'Value',numRefFrames);
        end
    else
        warndlg(sprintf('Could not delete the geometric reference frame "%s" because it is in use elsewhere.  Remove the reference frame dependencies before attempting to delete the reference frame.', selRefFrame.getName()),'Cannot Delete Reference Frame','modal');
    end

% --- Executes on key press with focus on lvd_EditGeometricRefFramesGUI or any of its controls.
function lvd_EditGeometricRefFramesGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditGeometricRefFramesGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            uiresume(handles.lvd_EditGeometricRefFramesGUI);
        case 'enter'
            uiresume(handles.lvd_EditGeometricRefFramesGUI);
        case 'escape'
            uiresume(handles.lvd_EditGeometricRefFramesGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function neqn = p03_neqn ( )

%*****************************************************************************80
%
%% P03_NEQN returns the number of equations for problem P03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function title = p07_title ( title )

%*****************************************************************************80
%
%% P07_TITLE returns the title of problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'F(X) = X**3, only linear Newton convergence.';

  return
end


text:
% op_takeextras.m
% Jamie Near, McGill University 2018.
%
% USAGE:
% out=op_takeextras(in,index);
%
% DESCRIPTION:
% Extract the extras with indices corresponding to the 'index' input
% array.
%
% INPUTS:
% in     = input data in matlab structure format.
% index  = vector indicating the indices of the extras you would like to extract.
%
% OUTPUTS:
% out    = Output dataset consisting of extras extracted from the input.

function out=op_takeextras(in,index);

if in.dims.extras==0
    %Can't take extras because there are none:
    error('ERROR:  There are no averages in this dataset!  Aborting!');
elseif in.dims.extras==1
    %SHOULD NEVER HAPPEN (Time dimension should always be dim=1)
    error('ERROR:  dims.subSpecs==1.  This should never happen!  Aborting!');
elseif in.dims.extras==2
    fids=in.fids(:,index,:,:,:);
elseif in.dims.extras==3;
    fids=in.fids(:,:,index,:,:);
elseif in.dims.extras==4;
    fids=in.fids(:,:,:,index,:);
elseif in.dims.extras==5
    fids=in.fids(:,:,:,:,index);
end

fids=squeeze(fids);

%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%change the dims variables
dims.t=in.dims.t;
if length(index)==1
    dims.extras=0;
    if in.dims.coils > in.dims.extras
        dims.coils=in.dims.coils-1;
    else
        dims.coils=in.dims.coils;
    end
    if in.dims.subSpecs > in.dims.extras
        dims.subSpecs=in.dims.subSpecs-1;
    else
        dims.subSpecs=in.dims.subSpecs;
    end
    if in.dims.averages > in.dims.extras
        dims.averages=in.dims.averages-1;
    else
        dims.averages=in.dims.averages;
    end
elseif length(index)>1
    dims.coils=in.dims.coils;
    dims.averages=in.dims.averages;
    dims.subSpecs=in.dims.subSpecs;
    dims.extras=in.dims.extras;
end

%re-calculate the sz variable
sz=size(fids);

%change the number of extras
if dims.extras==0
    extras=1;
else
    extras=sz(dims.extras);
end

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.dims=dims;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [ x, w ] = cce ( l )

%*****************************************************************************80
%
%% CCE computes a Clenshaw Curtis Exponential quadrature rule based on level.
%
%  Discussion:
%
%    Our convention is that the abscissas are numbered from left to right.
%
%    The rule is defined on [0,1].
%
%    The integral to approximate:
%
%      Integral ( 0 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%    The input value of L selects the size of the rule as follows:
%    L = 1, N = 1;
%    1 < L, N = 2^(L-1)+1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer L, the level of the rule.
%    1 <= L.
%
%    Output, real X(N,1), the abscissas.
%
%    Output, real W(N,1), the weights.
%
  if ( l < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CCE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of L = %d\n', l );
    error ( 'CCE - Fatal error!' );
  end
%
%  Find the value of N according to the level.
%
  n = cce_order ( l );
%
%  Compute the points and weights.
%
  [ x, w ] = cc ( n );

  return
end


text:
function [I J] = itriu(sz, k)
% function [I J] = itriu(sz) % OR
% I = itriu(sz) OR
%
% Return the subindices [I J] (or linear indices I if single output call)
% in the purpose of extracting an upper triangular part of the matrix of
% the size SZ. Input k is optional shifting. For k=0, extract from the main
% diagonal. For k>0 -> above the diagonal, k<0 -> below the diagonal
%
% This returnd same as [...] = find(triu(ones(sz),k))
% - Output is a column and sorted with respect to linear indice
% - No intermediate matrix is generated, that could be useful for large
%   size problem
% - Mathematically, A(itriu(size(A)) is called (upper) "half-vectorization"
%   of A
%
% Example:
%
% A = [ 7     5     4
%       4     2     3
%       9     1     9
%       3     5     7 ]
%
% I = itriu(size(A))  % gives [1 5 6 9 10 11]'
% A(I)                % gives [7 5 2 4  3  9]' OR A(triu(A)>0)
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Date: 21/March/2009

if isscalar(sz)
    sz = [sz sz];
end
m=sz(1);
n=sz(2);

% Main diagonal by default
if nargin<2
    k=0;
end

nc = n-max(k,0); % number of columns of the triangular part
lo = ones(nc,1); % lower row indice for each column
hi = min((1:nc).'-min(k,0),m); % upper row indice for each column

if isempty(lo)
    I = zeros(0,1);
    J = zeros(0,1);
else
    c=cumsum([0; hi-lo]+1); % cumsum of the length
    I = accumarray(c(1:end-1), (lo-[0; hi(1:end-1)]-1), ...
                   [c(end)-1 1]);
    I = cumsum(I+1); % row indice
    J = accumarray(c,1);
    J(1) = 1 + max(k,0); % The row indices starts from this value
    J = cumsum(J(1:end-1)); % column indice
end

if nargout<2
    % convert to linear indices
    I = sub2ind([m n], I, J);
end

end % itriu

author: Sable

repo: mcbench-benchmarks


text:
function C=HestonCall(St,K,r,T,vt,kap,th,sig,rho,lda)
%--------------------------------------------------------------------------
%PURPOSE: computes the option price using Heston's model.
%--------------------------------------------------------------------------
%USAGE: C=HestonCall(St,K,r,sig,T,vt,kap,th,lda,rho)
%--------------------------------------------------------------------------
%INPUT: St - scalar or vector, price of underlying at time t
%       K - scalar or vector, strike price
%       r - scalar or vector, continuously compound risk free rate expressed as a
%       positive decimal number.
%       sig- scalar or vector, volatility of the volatility of the
%       underlying(same time units as for r)
%       T - scalar or vector, time to maturity (same time units as for r)
%       vt - scalar or vector, instantaneous volatility
%       kap - scalar or vector, is the rate at which vt reverts to th
%       th - scalar or vector, is the long vol, or long run average price
%       volatility; as t tends to infinity, the expected value of ?t tends to ?
%       lda - scalar or vector, risk premium for volatility
%       rho - scalar or vector, correlation between underlying and
%       volatility (rho<0 generates the leverage effect)
%--------------------------------------------------------------------------
%OUTPUT: C - scalar or vector, Heston's model call option price
%--------------------------------------------------------------------------

dphi=0.01;
maxphi=50;
phi=(eps:dphi:maxphi)';

%f1 = CF_SVj(log(St),vt,T,r,kap*th,0.5,kap+lda-rho*sig,rho,sig,phi);
%P1 = 0.5+(1/pi)*sum(real(exp(-i*phi*log(K)).*f1./(i*phi))*dphi);
%f2 = CF_SVj(log(St),vt,T,r,kap*th,-0.5,kap+lda,rho,sig,phi);
%P2 = 0.5+(1/pi)*sum(real(exp(-i*phi*log(K)).*f2./(i*phi))*dphi);
%C = St*P1 -K*exp(-r*T)*P2;

f1 = CF_SVj(log(St),vt,T,0,kap*th,0.5,kap+lda-rho*sig,rho,sig,phi);
P1 = 0.5+(1/pi)*sum(real(exp(-i*phi*log(K)).*f1./(i*phi))*dphi);
f2 = CF_SVj(log(St),vt,T,0,kap*th,-0.5,kap+lda,rho,sig,phi);
P2 = 0.5+(1/pi)*sum(real(exp(-i*phi*log(K)).*f2./(i*phi))*dphi);
C = St*P1 -K*exp(-r*T)*P2;


text:
function [Y, R, E] = IsomapII(D, n_fcn, n_size, options)
% ISOMAPII   Computes Isomap embedding using an advanced version of
%             the algorithm in Tenenbaum, de Silva, and Langford (2000),
%             which can take advantage of sparsity in the graph and
%             redundancy in the distances.
%
% [Y, R, E] = isomapII(D, n_fcn, n_size, options);
%
% Input:
%    D = input-space distances between pairs of N points, which can
%     take 1 of 3 forms:
%      (1) a full N x N matrix (as in isomap.m)
%      (2) a sparse N x N matrix (missing entries are treated as INF)
%      (3) the name of a function (e.g. 'd_fun') that takes
%            one argument, i, and returns a row vector containng the
%            distances from all N points to point i.
%
%    n_fcn = neighborhood function ('epsilon' or 'k')
%    n_size = neighborhood size (value for epsilon or k)
%
%    options = optional structure of options:
%      options.dims = (row) vector of embedding dimensionalities to use
%                        (1:10 = default)
%      options.comp = which connected component to embed, if more than one.
%                        (1 = largest (default), 2 = second largest, ...)
%      options.display = plot residual variance and 2-D embedding?
%                        (1 = yes (default), 0 = no)
%      options.overlay = overlay graph on 2-D embedding?
%                        (1 = yes (default), 0 = no)
%      options.verbose = display progress reports?
%                        (1 = yes (default), 0 = no)
%      options.dijkstra = use dijkstra's algorithm for shortest paths with
%                         full N x N distance matrix.
%                         (1 = yes (default), 0 = use Floyd; Floyd should
%                          be used only if you are unable to MEX dijkstra.cpp)
%      options.Kmax = maximum number of neighbors (used for sparse versions
%                        of epsilon; by default, estimated by random sample)
%      options.landmarks = (row) vector of landmark points to use in MDS.
%                 (MDS finds the configuration that best approximates
%                  the distances from all points to the landmark points.
%                  The default landmark points are 1:N (i.e. all the points),
%                  which is equivalent to classical MDS.  Good
%                  results may often be obtained using a number of
%                  landmarks that is much smaller than N, but much
%                  larger than the data's intrinsic dimensionality.
%                  Note that this extension is experimental!  For
%                  discussion, see Steyvers, de Silva, and Tenenbaum
%                  (in preparation).)
%
% Output:
%    Y = Y.coords is a cell array, with coordinates for d-dimensional embeddings
%         in Y.coords{d}.  Y.index contains the indices of the points embedded.
%    R = residual variances for embeddings in Y
%    E = edge matrix for neighborhood graph
%

%    BEGIN COPYRIGHT NOTICE
%
%    Isomap II code -- (c) 1998-2000 Josh Tenenbaum
%
%    This code is provided as is, with no guarantees except that
%    bugs are almost surely present.  Published reports of research
%    using this code (or a modified version) should cite the
%    article that describes the algorithm:
%
%      J. B. Tenenbaum, V. de Silva, J. C. Langford (2000).  A global
%      geometric framework for nonlinear dimensionality reduction.
%      Science 290 (5500): 2319-2323, 22 December 2000.
%
%    Comments and bug reports are welcome.  Email to jbt@psych.stanford.edu.
%    I would also appreciate hearing about how you used this code,
%    improvements that you have made to it, or translations into other
%    languages.
%
%    You are free to modify, extend or distribute this code, as long
%    as this copyright notice is included whole and unchanged.
%
%    END COPYRIGHT NOTICE
%
% Modified by Ramon Casero <rcasero@gmail.com>, University of Oxford.
%
% Version: 0.1.1
%
% This file is distributed as a derivative work of a third-party function
% with project Gerardus.
%
% http://code.google.com/p/gerardus/
%
% Original code downloaded from the Isomap Homepage
%
% http://isomap.stanford.edu/

%%%%% Step 0: Initialization and Parameters %%%%%

if nargin < 3
     error('Too few input arguments');
elseif nargin < 4
     options = struct('dims',1:10,'overlay',1,'comp',1,'display',1,'dijkstra',1,'verbose',1);
end

if isa(D, 'function_handle')
     mode = 3;
     d_func = D;
     N = length(d_func(1));
elseif issparse(D)
     mode = 2;
     N = size(D,1);
     if ~(N==size(D,2))
          error('D must be a square matrix');
     end;
else
     mode = 1;
     N = size(D,1);
     if ~(N==size(D,2))
         error('D must be a square matrix');
     end;
end

if strcmp(n_fcn, 'k')
     K = n_size;
     if ~(K==round(K))
         error('Number of neighbors for k method must be an integer');
     end
     if ((mode==2) && ~(min(sum(D'>0))>=K))
         error('Sparse D matrix must contain at least K nonzero entries in each row');
     end
elseif strcmp(n_fcn, 'epsilon')
     epsilon = n_size;
     if isfield(options,'Kmax')
         K = options.Kmax;
     elseif (mode==3)    %% estimate maximum equivalent K %%
         tmp = zeros(10,N);
         for i=1:10
             tmp(i,:) = d_func(ceil(N*rand));
         end
         K = 2*max(sum(tmp'<epsilon));    % just to be safe
     end
else
     error('Neighborhood function must be either epsilon or k');
end

if (mode == 3)
     INF = inf;
else
     INF =  1000*max(max(D))*N;  %% effectively infinite distance
end

if ~isfield(options,'dims')
     options.dims = 1:10;
end
if ~isfield(options,'overlay')
     options.overlay = 1;
end
if ~isfield(options,'comp')
     options.comp = 1;
end
if ~isfield(options,'display')
     options.display = 1;
end
if ~isfield(options,'verbose')
     options.verbose = 1;
end
if ~isfield(options,'landmarks')
     options.landmarks = 1:N;
end
if ~isfield(options,'dijkstra')
     options.dijkstra = 1;
end
dims = options.dims;
comp = options.comp;
overlay = options.overlay;
displ = options.display;
verbose = options.verbose;
landmarks = options.landmarks;
use_dijk = options.dijkstra;

Y.coords = cell(length(dims),1);
R = zeros(1,length(dims));

%%%%% Step 1: Construct neighborhood graph %%%%%
disp('Constructing neighborhood graph...');

if ((mode == 1) && (use_dijk == 0))
     if strcmp(n_fcn, 'k')
         [tmp, ind] = sort(D);
         tic;
         for i=1:N
             D(i,ind((2+K):end,i)) = INF;
             if ((verbose == 1) && (rem(i,50) == 0))
                 disp([' Iteration: ' num2str(i) '     Estimated time to completion: ' num2str((N-i)*toc/60/50) ' minutes']); tic;
             end
         end
     elseif strcmp(n_fcn, 'epsilon')
         warning off    %% Next line causes an unnecessary warning, so turn it off
         D =  D./(D<=epsilon);
         D = min(D,INF);
         warning on
     end
     D = min(D,D');    %% Make sure distance matrix is symmetric
elseif ((mode == 1) && (use_dijk == 1))
     if n_fcn == 'k'
         [tmp, ind] = sort(D);
         tic;
         for i=1:N
             D(i,ind((2+K):end,i)) = 0;
             if ((verbose == 1) && (rem(i,50) == 0))
                 disp([' Iteration: ' num2str(i) '     Estimated time to completion: ' num2str((N-i)*toc/60/50) ' minutes']); tic;
             end
         end
     elseif strcmp(n_fcn, 'epsilon')
         D =  D.*(D<=epsilon);
     end
     D = sparse(D);
     D = max(D,D');    %% Make sure distance matrix is symmetric
elseif (mode == 2)
     if n_fcn == 'k'
         Di = zeros(N*K,1);      Dj = zeros(N*K,1);       Ds = zeros(N*K,1);
         counter = 0;
         [a,b,c] = find(D);
         tic;
         for i=1:N
             l = find(a==i);
             [g,f] = sort(c(l));
             Di(counter+(1:K)) = i;
             Dj(counter+(1:K)) = b(l(f(1:K)));
             Ds(counter+(1:K)) = g(1:K);
             counter = counter+K;
             if ((verbose == 1) && (rem(i,50) == 0))
                  disp([' Iteration: ' num2str(i) '     Estimated time to completion: ' num2str((N-i)*toc/60/50) ' minutes']); tic;
             end
         end
         D = sparse(Di(1:counter), Dj(1:counter), Ds(1:counter));
         clear Di Dj Ds counter;
     elseif strcmp(n_fcn, 'epsilon')
         D =  D.*(D<=epsilon);
     end
     D = max(D,D');    %% Make sure distance matrix is symmetric
elseif (mode == 3)
     Di = zeros(N*(K+1),1);      Dj = zeros(N*(K+1),1);       Ds = zeros(N*(K+1),1);
     counter = 0;
     tic;
     for i=1:N
         d = d_func(i);
         if n_fcn == 'k'
             [c,b] = sort(d);
             Di(counter+(1:(K+1))) = i;
             Dj(counter+(1:(K+1))) = b(1:(K+1));
             Ds(counter+(1:(K+1))) = c(1:(K+1));
             counter = counter+(K+1);
         elseif strcmp(n_fcn, 'epsilon')
             [a,b,c] = find(d.*(d<=epsilon));
             l = length(a);
             Di(counter+(1:l)) = i;
             Dj(counter+(1:l)) = b;
             Ds(counter+(1:l)) = c;
             counter = counter+l;
         end
         if ((verbose == 1) && (rem(i,50) == 0))
              disp([' Iteration: ' num2str(i) '     Estimated time to completion: ' num2str((N-i)*toc/60/50) ' minutes']); tic;
         end
     end
     D = sparse(Di(1:counter), Dj(1:counter), Ds(1:counter));
     clear Di Dj Ds counter;
     D = max(D,D');    %% Make sure distance matrix is symmetric
end

if (overlay == 1)
     if ((mode == 1) && (use_dijk == 0))
         E = int8(1-(D==INF));  %%  Edge information for subsequent graph overlay
     else
         [a,b,c] = find(D);
         E = sparse(a,b,ones(size(a)));
     end
end

%%%%% Step 2: Compute shortest paths %%%%%
disp('Computing shortest paths...');

if ((mode==1) && (use_dijk == 0))
     tic;
     for k=1:N
         D = min(D,repmat(D(:,k),[1 N])+repmat(D(k,:),[N 1]));
         if ((verbose == 1) && (rem(k,20) == 0))
              disp([' Iteration: ' num2str(k) '     Estimated time to completion: ' num2str((N-i)*toc/i/60) ' minutes']);
         end
     end
else
     D = dijkstra(D, landmarks);
end

%%%%% Step 3: Construct low-dimensional embeddings (Classical MDS) %%%%%
disp('Constructing low-dimensional embeddings (Classical MDS)...');

%%%%% Remove outliers from graph %%%%%
disp('  Checking for outliers...');

if ((mode == 1) && (use_dijk == 0))
     [tmp, firsts] = min(D==INF);     % first point each point connects to
else
     [tmp, firsts] = min(D==inf);     % first point each point connects to
end
[comps, I, J] = unique(firsts);    % first point in each connected component
n_comps = length(comps);           % number of connected components
size_comps = sum((repmat(firsts,n_comps,1)==((1:n_comps)'*ones(1,N)))');
                                   % size of each connected component
[tmp, comp_order] = sort(size_comps);  %% sort connected components by size
comps = comps(comp_order(end:-1:1));
size_comps = size_comps(comp_order(end:-1:1));
if (comp>n_comps)
     comp=1;                              %% default: use largest component
end
Y.index = find(firsts==comps(comp)); %% list of points in relevant component
Y.index = setdiff(Y.index,find(isinf(min(D)))); % prune points that don't connect
                                                % to any landmarks
N = length(Y.index);
[tmp, landmarks, land_ind] = intersect(landmarks,Y.index);
                                       % list of landmarks in component
nl = length(landmarks);
D = full(D(landmarks,Y.index))';
disp(['    Number of connected components in graph: ' num2str(n_comps)]);
disp(['    Embedding component ' num2str(comp) ' with ' num2str(length(Y.index)) ' points.']);

dims = unique(min(dims,nl-1));    % don't embed in more dimensions than landmarks-1
if (nl==N)
     opt.disp = 0;
     [vec, val] = eigs(-.5*(D.^2 - sum(D.^2)'*ones(1,N)/N - ones(N,1)*sum(D.^2)/N + sum(sum(D.^2))/(N^2)), max(dims), 'LR', opt);
else
     subB = -.5*(D.^2 - sum(D'.^2)'*ones(1,nl)/nl - ones(N,1)*sum(D.^2)/N+sum(sum(D.^2))/(N*nl));
     opt.disp = 0;
     [alpha,beta] = eigs(subB'*subB, max(dims), 'LR', opt);
     val = beta.^(1/2);
     vec = subB*alpha*inv(val);
     clear subB alpha beta;
end
h = real(diag(val));
[foo,sorth] = sort(h);  sorth = sorth(end:-1:1);
val = real(diag(val(sorth,sorth)));
vec = vec(:,sorth);

D = reshape(D,N*nl,1);
for di = 1:length(dims)
     Y.coords{di} = real(vec(:,1:dims(di)).*(ones(N,1)*sqrt(val(1:dims(di)))'))';
     r2 = 1-corrcoef(reshape(real(L2_distance(Y.coords{di}, Y.coords{di}(:,land_ind))),N*nl,1),D).^2;
     R(di) = r2(2,1);
     if (verbose == 1)
         disp(['  Isomap on ' num2str(N) ' points with dimensionality ' num2str(dims(di)) '  --> residual variance = ' num2str(R(di))]);
     end
end

clear D;

%%%%%%%%%%%%%%%%%% Graphics %%%%%%%%%%%%%%%%%%

if (displ==1)
     %%%%% Plot fall-off of residual variance with dimensionality %%%%%
     figure;
     hold on
     plot(dims, R, 'bo');
     plot(dims, R, 'b-');
     hold off
     ylabel('Residual variance');
     xlabel('Isomap dimensionality');

     %%%%% Plot two-dimensional configuration %%%%%
     twod = find(dims==2);
     if ~isempty(twod)
         figure;
         hold on;
         plot(Y.coords{twod}(1,:), Y.coords{twod}(2,:), 'ro');
         if (overlay == 1)
             gplot(E(Y.index, Y.index), [Y.coords{twod}(1,:); Y.coords{twod}(2,:)]');
             title('Two-dimensional Isomap embedding (with neighborhood graph).');
         else
             title('Two-dimensional Isomap.');
         end
         hold off;
     end
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% GEO_LOAD: create a geometry structure from a file or a cell-array of functions.
%
% geometry = geo_load (input)
%
% INPUT :
%
%   The input variable may be either
%   - a structure representing a NURBS surface or volume, as in the NURBS toolbox
%   - a string variable with the name of the file to be read (see doc/geo_specs_v21.txt)
%   - a cell-array of function handles, to evaluate the function, the first order
%       derivatives, and eventually the second order derivatives
%   - a 4x4 matrix representing an affine transformation
%
% OUTPUT:
%
%   geometry: a structure that contains, at least, the following fields
%             map:      a function handle to evaluate the parameterization
%             map_der:  a function handle to evaluate the derivatives of the parameterization
%             map_der2: a function handle to evaluate the second derivatives of the parameterization
%   The structure may contain further information. See the documentation.
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2013 Rafael Vazquez
% Copyright (C) 2014 Elena Bulgarello, Carlo de Falco, Sara Frizziero
% Copyright (C) 2015 Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function geometry = geo_load (in)

  if (isstruct (in) && isfield (in, 'form') && strcmpi (in.form, 'B-NURBS'))
%% geometry is given as a NURBS struct
    geometry.nurbs   = in;

  elseif (ischar (in))
%% load geometry from a file
    if (strcmpi (in(end-3:end), '.mat'))
      tmp  = load (in);
      geometry.nurbs   = tmp.geo;
    elseif (strcmpi (in(end-3:end), '.txt'))
      geometry = geo_read_nurbs (in);
%% load geometry from an xml file
    elseif (strcmpi (in(end-3:end), '.xml'))
      geometry.nurbs = xml_import (in);
    else
      error ('geo_load: unknown file extension');
    end

  elseif (iscell (in))
%% geometry is given as a cell array of functions
    geometry.map     =  in{1};
    geometry.map_der =  in{2};
    if (length (in) > 2)
      geometry.map_der2 =  in{3};
    end

  elseif (isnumeric (in) && all (size (in) == [4, 4]))
%% geometry is given as a 4x4 matrix representing an affine transformation
    geometry.map = @(ps) affine_map  (ps, in);
    geometry.map_der = @(ps) affine_map_der  (ps, in);
    geometry.map_der2 = @affine_map_der2;
  else
    error ('geo_load: wrong input type');
  end

% In case the geometry is a NURBS structure, we use the NURBS toolbox
% This allows to use a rectangular parametric domain, instead of [0,1]^n.
  if (isfield (geometry, 'nurbs'))
    if (any (abs(geometry.nurbs.coefs(3,:)) > 1e-12))
      rdim = 3;
    elseif (any (abs(geometry.nurbs.coefs(2,:)) > 1e-12) || (isfield(in, 'boundary_flag') && in.boundary_flag))
      rdim = 2;
    else
      rdim = 1;
    end
    geometry.rdim = rdim;

    wrn_struct = warning ('query', 'nrbderiv:SecondDerivative');
    warning ('off', 'nrbderiv:SecondDerivative')
    [deriv, deriv2] = nrbderiv (geometry.nurbs);
    geometry.dnurbs = deriv;
    geometry.dnurbs2 = deriv2;

    geometry.map      =  @(PTS) geo_nurbs (geometry.nurbs, deriv, deriv2, PTS, 0, rdim);
    geometry.map_der  =  @(PTS) geo_nurbs (geometry.nurbs, deriv, deriv2, PTS, 1, rdim);
    geometry.map_der2 =  @(PTS) geo_nurbs (geometry.nurbs, deriv, deriv2, PTS, 2, rdim);

    if (numel (geometry.nurbs.order) > 1)
      bnd = nrbextract (geometry.nurbs);
      for ibnd = 1:numel (bnd)
        [deriv, deriv2] = nrbderiv (bnd(ibnd));
        geometry.boundary(ibnd).nurbs    = bnd(ibnd);
        geometry.boundary(ibnd).dnurbs   = deriv;
        geometry.boundary(ibnd).dnurbs2  = deriv2;
        geometry.boundary(ibnd).rdim     = rdim;
        geometry.boundary(ibnd).map      = @(PTS) geo_nurbs (bnd(ibnd), deriv, deriv2, PTS, 0, rdim);
        geometry.boundary(ibnd).map_der  = @(PTS) geo_nurbs (bnd(ibnd), deriv, deriv2, PTS, 1, rdim);
        geometry.boundary(ibnd).map_der2 = @(PTS) geo_nurbs (bnd(ibnd), deriv, deriv2, PTS, 2, rdim);
      end
    end
    if (strcmpi (wrn_struct.state, 'on'))
      warning ('on', 'nrbderiv:SecondDerivative')
    end
  else
    for ibnd = 1:6 % This loop should be until 2*ndim, but ndim is not known
      geometry.boundary(ibnd).map     = @(PTS) boundary_map (geometry.map, ibnd, PTS);
      geometry.boundary(ibnd).map_der = @(PTS) boundary_map_der (geometry.map, geometry.map_der, ibnd, PTS);
    end
  end

end

function mps = affine_map  (ps, in)
  if (iscell (ps))
    ndim = numel (ps);
    npts = cellfun (@numel, ps);
    nps = prod (npts);
    if (ndim == 2)
      u = reshape (repmat (ps{1}(:), 1, npts(2)), 1, []);
      v = reshape (repmat (ps{2}(:)', npts(1), 1), 1, []);
      ps = [u; v];
    elseif (ndim == 3)
      u = reshape (ps{1}, npts(1), 1, 1);
      u = reshape (repmat (u, [1, npts(2), npts(3)]), 1, []);
      v = reshape (ps{2}, 1, npts(2), 1);
      v = reshape (repmat (v, [npts(1), 1, npts(3)]), 1, []);
      w = reshape (ps{3}, 1, 1, npts(3));
      w = reshape (repmat (w, [npts(1), npts(2), 1]), 1, []);
      ps = [u; v; w];
    end
  else
    ndim = size (ps, 1);
    nps  = size (ps, 2);
  end
  mps  = in([1:ndim, 4], [1:ndim, 4]) * [ps; ones(1, nps)];
  mps  = mps (1:ndim, :);
end

function mps = affine_map_der  (ps, in)
  if (iscell (ps))
    ndim = numel (ps);
    npts = cellfun (@numel, ps);
    nps = prod (npts);
    if (ndim == 2)
      u = reshape (repmat (ps{1}(:), 1, npts(2)), 1, []);
      v = reshape (repmat (ps{2}(:)', npts(1), 1), 1, []);
      ps = [u; v];
    elseif (ndim == 3)
      u = reshape (ps{1}, npts(1), 1, 1);
      u = reshape (repmat (u, [1, npts(2), npts(3)]), 1, []);
      v = reshape (ps{2}, 1, npts(2), 1);
      v = reshape (repmat (v, [npts(1), 1, npts(3)]), 1, []);
      w = reshape (ps{3}, 1, 1, npts(3));
      w = reshape (repmat (w, [npts(1), npts(2), 1]), 1, []);
      ps = [u; v; w];
    end
  else
    ndim = size (ps, 1);
    nps  = size (ps, 2);
  end
  mps  = repmat (in(1:ndim, 1:ndim), [1, 1, nps]);
end

function mps = affine_map_der2  (ps)
  if (iscell (ps))
    ndim = numel (ps);
    nps = prod (cellfun (@numel, ps));
  else
    ndim = size (ps, 1);
    nps  = size (ps, 2);
  end
  mps  = zeros (ndim, ndim, ndim, nps);
end

% These two functions are to compute boundary entities from the global ones
function F = boundary_map (map, iside, pts)

%%    ind  = [2 3; 2 3; 1 3; 1 3; 1 2; 1 2] in 3D, %ind  = [2 2 1 1] in 2D;
%%    ind2 = [1 1 2 2 3 3] in 3D,                  %ind2 = [1 1 2 2] in 2D
  ind2 = ceil (iside/2);

  if (iscell (pts))
    ndim = numel (pts) + 1;
    ind = setdiff (1:ndim, ind2);

    pts_aux(ind) = pts;
    if (mod (iside, 2) == 1)
      pts_aux{ind2} = 0;
    else
      pts_aux{ind2} = 1;
    end
  else
    error ('For the boundary, a cell array should be passed as the argument')
  end

  F = map (pts_aux);

end

function varargout = boundary_map_der (map, map_der, iside, pts)

%%    ind  = [2 3; 2 3; 1 3; 1 3; 1 2; 1 2] in 3D, %ind  = [2 2 1 1] in 2D;
%%    ind2 = [1 1 2 2 3 3] in 3D,                  %ind2 = [1 1 2 2] in 2D
  ind2 = ceil (iside/2);

  if (iscell (pts))
    ndim = numel (pts) + 1;
    ind = setdiff (1:ndim, ind2);

    pts_aux(ind) = pts;
    if (mod (iside, 2) == 1)
      pts_aux{ind2} = 0;
    else
      pts_aux{ind2} = 1;
    end
  else
    error ('For the boundary, a cell array should be passed as the argument')
  end

  DF = map_der (pts_aux);
  DF = DF(:,ind,:);
  if (nargout == 1)
    varargout{1} = DF;
  elseif (nargout == 2)
    F = map (pts_aux);
    varargout{1} = F;
    varargout{2} = DF;
  end

end

%!shared g1,g2,x
%!test
%! g1 = geo_load ('geo_ring.txt');
%! g2 = geo_load ('ring.mat');
%! for [v, k] = g1.nurbs
%!  assert (isequal (v, g2.nurbs.(k)))
%! endfor
%! x = rand (2, 10);
%!test
%! assert (g1.map (x), g2.map (x))
%!test
%! assert (g1.map_der (x), g2.map_der (x))

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ a, ipvt, info ] = chifa ( a, lda, n )

%*****************************************************************************80
%
%% CHIFA factors a complex hermitian matrix.
%
%  Discussion:
%
%    CHIFA performs the factoring by elimination with symmetric pivoting.
%
%    To solve A*X = B, follow CHIFA by CHISL.
%
%    To compute inverse(A)*C, follow CHIFA by CHISL.
%
%    To compute determinant(A), follow CHIFA by CHIDI.
%
%    To compute inertia(A), follow CHIFA by CHIDI.
%
%    To compute inverse(A), follow CHIFA by CHIDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%
%  Parameters:
%
%    Input, complex A(LDA,N); the hermitian matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); a block diagonal matrix and the multipliers which
%    were used to obtain it.  The factorization can be written
%    A = U*D*hermitian(U) where U is a product of permutation and unit upper
%    triangular matrices, hermitian(U) is the conjugate transpose of U, and
%    D is block diagonal with 1 by 1 and 2 by 2 blocks.  Only the diagonal
%    and upper triangle are used.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO.
%    0, normal value.
%    K, if the K-th pivot block is singular.  This is not an error condition
%    for this subroutine, but it does indicate that CHISL or CHIDI may
%    divide by zero if called.
%

%
%  Initialize.
%
%  ALPHA is used in choosing pivot block size.
%
  alpha = ( 1.0 + sqrt ( 17.0 ) ) / 8.0;

  info = 0;
%
%  Main loop on K, which goes from N to 1.
%
  k = n;

  while ( 1 )
%
%  Leave the loop if K = 0 or K = 1.
%
    if ( k == 0 )
      break
    end

    if ( k == 1 )
      ipvt(1) = 1;
      if ( cabs1 ( a(1,1) ) == 0.0 )
        info = 1;
      end
      break
    end
%
%  This section of code determines the kind of
%  elimination to be performed.  When it is completed,
%  KSTEP will be set to the size of the pivot block, and
%  SWAP will be set to TRUE if an interchange is
%  required.
%
    km1 = k - 1;
    absakk = cabs1 ( a(k,k) );
%
%  Determine the largest off-diagonal element in column K.
%
    imax = icamax ( k-1, a(1:k-1,k), 1 );
    colmax = cabs1 ( a(imax,k) );

    if ( alpha * colmax <= absakk )

      kstep = 1;
      swap = 0;

    else
%
%  Determine the largest off-diagonal element in row IMAX.
%
      rowmax = 0.0;

      for j = imax + 1 : k
        rowmax = max ( rowmax, cabs1 ( a(imax,j) ) );
      end

      if ( imax ~= 1 )
        jmax = icamax ( imax-1, a(1:imax-1,imax), 1 );
        rowmax = max ( rowmax, cabs1 ( a(jmax,imax) ) );
      end

      if ( alpha * rowmax <= cabs1 ( a(imax,imax) )  )
        kstep = 1;
        swap = 1;
      elseif ( alpha * colmax * ( colmax / rowmax ) <= absakk )
        kstep = 1;
        swap = 0;
      else
        kstep = 2;
        swap = ( imax ~= km1 );
      end

    end
%
%  Column K is zero.  Set INFO and iterate the loop.
%
    if ( max ( absakk, colmax ) == 0.0 )
      ipvt(k) = k;
      info = k;
      k = k - kstep;
      continue
    end

    if ( kstep ~= 2 )
%
%  1 x 1 pivot block.
%
      if ( swap )

        temp           = a(1:imax,imax);
        a(1:imax,imax) = a(1:imax,k);
        a(1:imax,k)    = temp;

        for jj = imax : k
          j = k + imax - jj;
          t         = conj ( a(j,k) );
          a(j,k)    = conj ( a(imax,j) );
          a(imax,j) = t;
        end

      end
%
%  Perform the elimination.
%
      for jj = 1 : km1
        j = k - jj;
        mulk = -a(j,k) / a(k,k);
        t = conj ( mulk );
        a(1:j,j) = a(1:j,j) + t * a(1:j,k);
        a(j,j) = real ( a(j,j) );
        a(j,k) = mulk;
      end
%
%  Set the pivot array.
%
      ipvt(k) = k;

      if ( swap )
        ipvt(k) = imax;
      end

    else
%
%  2 x 2 pivot block.
%
      if ( swap )

        temp           = a(1:imax,imax);
        a(1:imax,imax) = a(1:imax,k-1);
        a(1:imax,k-1)  = temp;

        for jj = imax : km1
          j = km1 + imax - jj;
          t         = conj ( a(j,k-1) );
          a(j,k-1)  = conj ( a(imax,j) );
          a(imax,j) = t;
        end

        t         = a(k-1,k);
        a(k-1,k)  = a(imax,k);
        a(imax,k) = t;

      end
%
%  Perform the elimination.
%
      km2 = k - 2;

      if ( 0 < k - 2 )

        ak = a(k,k) / a(k-1,k);
        akm1 = a(k-1,k-1) / conj ( a(k-1,k) );
        denom = 1.0 - ak * akm1;

        for jj = 1 : k - 2

          j = km1 - jj;
          bk = a(j,k) / a(k-1,k);
          bkm1 = a(j,k-1) / conj ( a(k-1,k) );
          mulk = ( akm1 * bk - bkm1 ) / denom;
          mulkm1 = ( ak * bkm1 - bk ) / denom;
          t = conj ( mulk );
          a(1:j,j) = a(1:j,j) + t * a(1:j,k);
          t = conj ( mulkm1 );
          a(1:j,j) = a(1:j,j) + t * a(1:j,k-1);
          a(j,k) = mulk;
          a(j,k-1) = mulkm1;
          a(j,j) = real ( a(j,j) );

        end

      end
%
%  Set the pivot array.
%
      if ( swap )
        ipvt(k) = -imax;
      else
        ipvt(k) = 1 - k;
      end

      ipvt(k-1) = ipvt(k);

    end

    k = k - kstep;

  end

  return
end


text:
report_this_filefun(mfilename('fullpath'));

%l = isnan(tmap);
%tmap(l) = 1;

%l = tmap < 1;
%tmap(l) = 1;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);

%  re4 = re3;

l = re4 < -2.0;
re4(l) = -2.0;
l = re4 > 2.0;
re4(l) = 2.0;
re4 = re4+1200;

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-20;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[33.0 36.2],'MapLonLimit',[-118.1 -115.3])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',1.5);
tightmap
%view([10 30])
hl = camlight ; lighting phong
set(gca,'projection','perspective');

aa = a;

l = aa(:,6) > 1.5 & aa(:,3) < 1992.2;
a = aa.subset(l);
for i = 1:length(a)

    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'or');
    hold on
    fac = 64/max(a.Depth);

    facm = 9/max(a.Magnitude);
    sm = aa(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm);
end

l = aa(:,6) > 1.5 & aa(:,3) > 1992.55;
a = aa.subset(l);
for i = 1:length(a)

    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'^b');
    hold on
    fac = 64/max(a.Depth);

    facm = 9/max(a.Magnitude);
    sm = aa(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm);
end
cd /alaskah2/stefan/ZMAP/agu99
load hec.mat

for i = 1:length(a)
    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'ok');
    hold on
    fac = 64/max(a.Depth);

    facm = 6/max(a.Magnitude);
    sm = a(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm,'markerfacecolor','y');
end

load land.mat

for i = 1:length(a)
    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'sk');
    ;
    hold on
    fac = 64/max(a.Depth);

    facm = 6/max(a.Magnitude);
    sm = a(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm,'markerfacecolor','y');
end
load bb.mat

for i = 1:length(a)
    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'sk');

    hold on
    fac = 64/max(a.Depth);

    facm = 6/max(a.Magnitude);
    sm = a(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm,'markerfacecolor','y');
end

for i = 1:length(maepi(:,3))
    pl =plotm(maepi(i,2),maepi(i,1),'hk');
    set(pl,'Markersize',18,'markerfacecolor','w');
end

plotm(faults(:,2), faults(:,1),'k','Linewidth',1.0);

zdatam(handlem('allline'),10000) % keep line on surface
j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j;  [ 0.85 0.9 0.9] ]
%caxis([ min(min(re4)) max(max(re4)) ]);

colormap(j); brighten(0.1);
caxis([-2.10 2.10]);
axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',3);

setm(gca,'mlabellocation',1)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',1)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',12,'Labelunits','dm')

set(gcf,'Inverthardcopy','off');

a = aa;

author: CelsoReyes


text:
function [] = analyzeParVariation(logprefix, signalname, signalname_ref)
%
% Author: Thomas Herrmann <thomas.herrmann@tum.de>     Date: 22-10-2020
%
% Description:
%   extracts log data from all Simulink simulation logs and displays in one
%   plot. Purpose of this script is the tuning of the TMPC algorithm.
%
% Input:
%   logprefix:      Prefix of all log files from simulation
%   signalname:     name of signal to be compared
%   signalname_ref: name of reference signal
%
%% Algorithm

% Get all log files with specific prefix
F = dir([logprefix, '*.*']);
% Allocate legend
lgd = struct();
F_data_struct = struct();

figure; hold on; grid on;
legend('-DynamicLegend');

for i = 1:length(F)
    % Get data from log file
    data_tmp = load(F(i).name);

    % get file identifier with specific parameter set
    fname = strrep(...
        strrep(F(i).name,[logprefix, '_'],''), ...
        '.mat', '');

    % append to struct for later use
    F_data_struct.(fname) = getfield(data_tmp.debug, signalname);

    % plot different paramter set signals over s coordinate
    plot(data_tmp.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_s_glob_m.Data, ...
        F_data_struct.(fname).Data, ...
        'DisplayName', fname);

    set(0, 'DefaultLegendInterpreter', 'none')

end

signal_ref = getfield(data_tmp.debug, signalname_ref);
plot(data_tmp.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_s_glob_m.Data, ...
    signal_ref.Data, ...
    'DisplayName', 'ref');
set(0, 'DefaultLegendInterpreter', 'none')
xlabel('s in m'); ylabel('v in mps');

end

author: TUMFTM


text:
function combined_waveforms = combine (waveformlist)
   %COMBINE merges waveforms based on start/end times and ChannelTag info.
   % combined_waveforms = combine (waveformlist) takes a vector of waveforms
   % and combines them based on SCNL information and start/endtimes.
   % DOES NO OTHER CHECKS

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$
   % Glenn Thompson 2018/01/09: replaced splice_waveform routine as it was
   % failing as part of drumplot.plot

   if numel(waveformlist) < 2  %nothing to do
      combined_waveforms = waveformlist;
      return
   end

   channelinfo = get(waveformlist,'channeltag');
   [uniquescnls, idx, scnlmembers] = unique(channelinfo);

   %preallocate
   combined_waveforms = repmat(waveform,size(uniquescnls));

   for i=1:numel(uniquescnls)
      w = waveformlist(scnlmembers == i);
      w = timesort(w);
      for j=(numel(w)-1):-1:1
         w(j) = piece_together(w(j:j+1));
         w(j+1) = waveform;
      end
      combined_waveforms(i) = w(1);
   end
end

function w = piece_together(w)
   if numel(w) > 2
      for i = numel(w)-1: -1 : 1
         w(i) = piece_together(w(i:i+1));
      end
      w = w(1);
      return;
   elseif numel(w) == 1
      return
   end
   if isempty(w(1))
      w = w(2);
      return;
   end;
   dt = dt_seconds(w(1),w(2));  %time overlap in seconds.
   sampleRates = round(get(w,'freq'));
   sampleInterval = 1 ./ sampleRates(1);

   if overlaps(dt, sampleInterval)
      w = spliceWaveform(w(1), w(2));
   else
      paddingAmount = round((dt * sampleRates(1))-1);
      w = spliceAndPad(w(1),w(2), paddingAmount);
   end
   w = w(1);
end

function w = spliceAndPad(w1, w2, paddingAmount)
   if paddingAmount > 0 && ~isinf(paddingAmount)
      toAdd = nan(paddingAmount,1);
   else
      toAdd = [];
   end

   w = set(w1,'data',[w1.data; toAdd; w2.data]);
end

% function w = spliceWaveform(w1, w2)
%
%    % NOTE * Function uses direct field access
%    timesToGrab = sum(get(w1,'timevector') < get(w2,'start'));
%
%    samplesRemoved = numel(w1.data) - timesToGrab;
%
%    w = set(w1,'data',[double(extract(w1,'index',1,timesToGrab)); w2.data]);
%
%    w= addhistory(w,'SPLICEPOINT: %s, removed %d points (overlap)',...
%       datestr(get(w2,'start')),samplesRemoved);
% end

function w = spliceWaveform(w1, w2)
%
%    % NOTE * Function uses direct field access
    snum1 = get(w1,'start');
    snum2 = get(w2,'start');
    enum1 = get(w1,'end');
    enum2 = get(w2,'end');
    [snummin, snuminindex] = min([snum1 snum2]);
    [enummax, enumaxindex] = max([enum1 enum2]);
    dnum1 = get(w1,'timevector');
    dnum2 = get(w2,'timevector');
    data1 = get(w1,'data');
    data2 = get(w2,'data');
    w = w1;
    dnum=dnum1;
    data=data1;
    ind = find(dnum2<snum1);
    if numel(ind)>0
        dnum = [dnum2(ind); dnum];
        data = [data2(ind); data];
    end
    ind = find(dnum2>enum1);
    if numel(ind)>0
        dnum = [dnum; dnum2(ind)];
        data = [data; data2(ind)];
    end
    snum = dnum(1);
    w = set(w, 'start', snum, 'data', data);

end

function result = overlaps(dt, sampleInterval)
   result = (dt- sampleInterval .* 1.25) < 0;
end

function t = dt_seconds(w1,w2)
   %  w1----] t [----w2
   firstsampleT = get(w2,'start');
   lastsampleT = get(w1,'timevector'); lastsampleT = lastsampleT(end);
   t = firstsampleT*86400 - lastsampleT * 86400;
end

function w = timesort(w)
   [~, I] = sort(get(w,'start'));
   w = w(I);
end

author: geoscience-community-codes

repo: GISMO