text:
function [ Model2 ] = NNTest()

    % Extract Features From the Images
    [Features] = LoadImages();
    % Load The Image Classes
    [yTrain] = LoadLabels();

    % Tweak these parameters to get optimum accuracy
    varianceThreshold = 0.90;
    thresholdDifference = 0.0004;
    regularizationRate = 6.5;
    learningRate = 1;
    hiddenNodes = 175;

    [xTrain, projection] = BestFeats(Features, varianceThreshold);

    % Train the Neural Net.
    [ weights1, weights2, ~ ] = NNTrain(xTrain, yTrain, hiddenNodes, ...
                                    learningRate, regularizationRate, ...
                                    thresholdDifference);
    Model2 = struct('weights1', weights1, 'weights2', weights2, ...
             'projection', projection);
    save('Model2.mat', 'Model2');

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function gen_cropped_MPII_validate_images()
% Author: Guanghan Ning, Feb 2017
clear all;
clc;

target_dist = 1.45;
boxsize = 368;

flag_showScaledImg = 1;
% Get [input images paths] and [output images path]
mpii_img_folder_path = '../dataset/MPI/images/';
img_output_folder = 'dataset_mpii/images_validate_cropped/qualitative/';
if exist(img_output_folder) ~= 7  % if not a folder path
    mkdir(img_output_folder);
end

% Load annotations
load('eval_MPII/mpii_Tompson_validation', 'annolist');
load('eval_MPII/mpii_Tompson_validation', 'keypointsAll');
load('eval_MPII/mpii_Tompson_validation', 'RELEASE_img_index');
load('eval_MPII/mpii_Tompson_validation', 'RELEASE_person_index');

annolist_test = annolist;
num_test_imgs = size(annolist_test, 2);

load('eval_MPII/mpii_human_pose_v1_u12_1','RELEASE');
annolist_all = RELEASE.annolist;

%% Start from here if everything already read from mat file.
for img_id = 1:num_test_imgs

    img_index = RELEASE_img_index(img_id);

    % (1) Deal with rect ids

        % Get rect id. Note the difference between rect_ct and rect_id
        rect_id = RELEASE_person_index(img_id);

        % Find scale and mid-points in order to convert prediction
        x_mid = annolist_test(img_id).annorect.objpos.x;
        y_mid = annolist_test(img_id).annorect.objpos.y;
        objpos_org = [x_mid, y_mid];
        scale_provided = annolist_test(img_id).annorect.scale;

        % (2) Deal with images
        img_name = annolist_all(img_index).image.name;
        img_path = strcat(mpii_img_folder_path, img_name);
        if ~exist(img_path, 'file')
            continue
        end

        img = imread(img_path);

        scale_in_cpp = target_dist/scale_provided;
        center_s = objpos_org * scale_in_cpp;
        fprintf('scale: %f\n', scale_in_cpp);

        img = imresize(img, scale_in_cpp);
        disp(size(img));
        [img, pad] = padAround(img, boxsize, center_s);
        disp(size(img));

        [pathstr, name, ext] = fileparts(img_path);
        img_output_path = strcat(img_output_folder, name, '_', num2str(rect_id), ext);
        imwrite(img, img_output_path);

        output_scale(img_id) = scale_in_cpp;
end
save('scale_in_cpp_MPII_validate.mat', 'output_scale');

function [img_padded, pad] = padAround(img, boxsize, center)
    center = round(center);
    h = size(img, 1);
    w = size(img, 2);
    pad(1) = boxsize/2 - center(2); % up
    pad(3) = boxsize/2 - (h-center(2)); % down
    pad(2) = boxsize/2 - center(1); % left
    pad(4) = boxsize/2 - (w-center(1)); % right

    pad_pixel_value = 128;

    pad_up = repmat(img(1,:,:)*0, [pad(1) 1 1]) + pad_pixel_value;
    img_padded = [pad_up; img];
    pad_left = repmat(img_padded(:,1,:)*0, [1 pad(2) 1]) + pad_pixel_value;
    img_padded = [pad_left img_padded];
    pad_down = repmat(img_padded(end,:,:)*0, [pad(3) 1 1]) + pad_pixel_value;
    img_padded = [img_padded; pad_down];
    pad_right = repmat(img_padded(:,end,:)*0, [1 pad(4) 1]) + pad_pixel_value;
    img_padded = [img_padded pad_right];

    center = center + [max(0,pad(2)) max(0,pad(1))];
    img_padded = img_padded(center(2)-(boxsize/2-1):center(2)+boxsize/2, center(1)-(boxsize/2-1):center(1)+boxsize/2, :); %cropping if needed

function [x,y] = findMaximum(map)
    [~,i] = max(map(:));
    [x,y] = ind2sub(size(map), i);

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function varargout = squeeze(varargin)
%SQUEEZE   Squeeze a CHEBFUN2 to one variable, if possible.
%   G = squeeze(F) returns a CHEBFUN2 if F depends on x and y. If F depends only
%   on the x-variable a row CHEBFUN is returned and if it depends on just the
%   y-variable a column CHEBFUN is returned.

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

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

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p=baselineParameters()

p=struct();

p.h = 0.8;
p.gz = 0.001;
p.pi_ss = 1.005;
p.beta = 0.9975;
p.kappa1 = 5;
p.kappa2 = 0.5;
p.epsilon = 6;
p.phi = 10;
p.eta = 2;
p.xi = 0.5;
p.rho_r = 0.8;
p.rho_a = 0.9;
p.sig_r = 0.01;
p.sig_a = 0.01;
p.sig_pi = 0.01;
p.sig_z = 0.01;
p.sigma = 2;
p.d = 1;

end

author: jmaih


text:
function [s,fs]=bat()
%BAT  Load the 'bat' test signal
%   Usage:  s=bat;
%
%   `bat` loads the 'bat' signal. It is a 400 samples long recording
%   of a bat chirp sampled with a sampling period of 7 microseconds.
%   This gives a sampling rate of 143 kHz.
%
%   `[sig,fs]=bat` additionally returns the sampling frequency *fs*.
%
%   The signal can be obtained from
%   `<http://dsp.rice.edu/software/bat-echolocation-chirp>`_
%
%   Please acknowledge use of this data in publications as follows:
%
%     The author wishes to thank Curtis Condon, Ken White, and Al Feng of
%     the Beckman Institute of the University of Illinois for the bat data
%     and for permission to use it in this paper.
%
%   Examples:
%   ---------
%
%   Plot of 'bat' in the time-domain:::
%
%     plot((1:400)/143000,bat);
%     xlabel('Time (seconds)');
%     ylabel('Amplitude');
%
%   Plot of 'bat' in the frequency-domain:::
%
%     plotfftreal(fftreal(bat),143000,90);
%
%   Plot of 'bat' in the time-frequency-domain:::
%
%     sgram(bat,143000,90);
%
%   See also:  batmask

%   AUTHOR : Peter L. Søndergaard
%   TESTING: TEST_SIGNALS
%   REFERENCE: OK

if nargin>0
  error('This function does not take input arguments.')
end;

f=mfilename('fullpath');

s=load('-ascii',[f,'.asc']);
fs=143000;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function varargout = process_opposite( varargin )
% PROCESS_OPPOSITE: Opposite values.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Opposite values';
    sProcess.FileTag     = 'opposite';
    sProcess.Category    = 'Filter';
    sProcess.SubGroup    = 'Pre-process';
    sProcess.Index       = 74;
    sProcess.Description = '';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Default values for some options
    sProcess.processDim  = 2;    % Process time by time
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data', 'raw'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function sInput = Run(sProcess, sInput) %#ok<DEFNU>
    % Opposite values
    sInput.A = -sInput.A;
end

author: brainstorm-tools


text:
%SerialLink.perturb Perturb robot parameters
%
% RP = R.perturb(P) is a new robot object in which the dynamic parameters (link
% mass and inertia) have been perturbed.  The perturbation is multiplicative so
% that values are multiplied by random numbers in the interval (1-P) to (1+P).
% The name string of the perturbed robot is prefixed by 'P/'.
%
% Useful for investigating the robustness of various model-based control
% schemes. For example to vary parameters in the range +/- 10 percent is:
%    r2 = p560.perturb(0.1);
%
% See also SerialLink.rne.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function  r2 = perturb(r, p)

	if nargin == 1
		p = 0.1;	% 10 percent disturb by default
	end

    r2 = SerialLink(r);

    links = r2.links;
	for i=1:r.n
		s = (2*rand-1)*p + 1;
		links(i).m = links(i).m * s;

		s = (2*rand-1)*p + 1;
		links(i).I = links(i).I * s;
	end

	r2.name = ['P/' r.name];

author: petercorke

repo: robotics-toolbox-matlab


text:
%% Dimensional Splitting: Front Tracking
% In this example, we consider the 2D scalar conservation law
%
% $$ u_t + (u^2/2)_x + (u^3/3)_y = 0, \qquad u(x,y,0)=u_0(x,y)$$
%
% which we solve using dimensional splitting. That is, we introduce two
% one-dimensional operators
%
% $$ S_x(t): \quad v_t + (v^2/2)_x = 0, \qquad v(x,0)=v_0(x)$$
%
% $$ S_y(t): \quad w_t + (w^3/3)_y = 0, \qquad w(y,0)=w_0(y)$$
%
% and construct the approximate solution from the formula
%
% $$u(x,t)\approx [S_x(0.5\Delta t)\circ S_y(\Delta t)\circ S_x(0.5\Delta t) ]^n u_0(x)$$
%
% As a one-dimensional solver, we will use front tracking. This algorithm
% makes a piecewise constant approximation to the flux function
% and then solves the 1-D hyperbolic problem *exactly*, given piecewise
% constant initial data.
%
% *Disclaimer*: The front-tracking code is built using arrays and is
% admittedly very slow, e.g., compared with the <testNTds.html NT
% central-difference code>. It is only meant to illustrate the concepts.
% For serious computing, one should use e.g., a C/C++ code built using
% pointers.

%% Initial setup
xmin=-1.0; xmax=1.0;
T = 2.0;
N = 50;
h = (xmax-xmin)/N;
x = xmin:h:xmax; y = 0.5*(x(1:end-1)+x(2:end));
[X,Y] = meshgrid(y,y);

%%
% The initial function u0 is equal -1 inside a circle of radius 0.4
% centered at (-1/2,-1/2), equal 1 inside a circle of radius 0.4 centered
% at (1/2,1/2), and zero otherwise. For later use, we define an anonymous
% function to do the computation of u0.
initData = (@(x,y) 1.0*((x+0.5).^2 + (y+0.5).^2 < 0.16) ...
   - 1.0*((x-0.5).^2 + (y-0.5).^2 < 0.16));
u0 = initData(X,Y);
contourf(X,Y,u0,-1:0.25:1); axis equal; colorbar, title('Initial data')

%% Number of time steps
% In the first example, we fix the grid resolution and compare the
% approximate solutions generated with four different splitting steps (n=1,
% 4, 16, 64). Generally, it makes sense to match the resolution of the
% front-tracking algorithm with the grid resolution. The resolution of the
% front-tracking method is determined by segments in the piecewise linear
% approximation of the flux functions. Here, we simply scale these segments
% to be inversely proportional with sqrt(N).
for i=1:4,
   tic;
   u=DimSplit(u0,x,x,'burgerRsol','cubRsol',0.1/sqrt(N),4^(i-1), T,...
      'wbar','periodic',[xmin xmax],[xmin xmax]);
   t=toc; st=sprintf('%4.2f',t);

   subplot(2,2,i);
   contourf(y,y,u,30), axis equal image; colorbar
   title( sprintf('Grid: %dx%d  Steps: %d', N, N, 4^(i-1)));
   xlabel(['Time used: ', st,' sec.']);
   p=get(gca,'position'); p([3 4])=p([3 4])+0.03;
   set(gca,'position',p,'XTickLabel',[]);
end
%%
% Looking at the plots, it is amazing to observe how well the operator
% splitting method resolves the dynamics of the problem using only one
% splitting step. The qualitative behavior of the red and blue "plumes" is
% captured almost correctly, even though the shape and position of the
% curved  shock is not correct. Increasing the number of splitting steps to
% four (CFL number of 12.5) ensures that all qualitative features of the
% solution are captured correctly. For 64 time steps, the effective CFL
% number is less than one and the front-tracking method is equivalent to a
% standard Godunov method. By considering accuracy versus runtime, the
% optimal number of splitting steps is most likely somewhere between four
% and sixteen.

%% Increasing grid resolution
% In the second example, we fix the CFL number to 20 and consider four
% different grid resolutions.
CFL=20;
for i=1:4,
   N= 32*(2^(i-1));
   h = (xmax-xmin)/N; delta=0.1/sqrt(N);
   x = xmin:h:xmax;
   y = 0.5*(x(1:end-1)+x(2:end));
   [X,Y] = meshgrid(y,y);
   u0=initData(X,Y);

   tic;
   Nstep=ceil(T/(CFL*h));
   u=DimSplit(u0,x,x,'burgerRsol','cubRsol',delta,Nstep,T,...
      'wbar','periodic',[xmin xmax],[xmin xmax]);
   t=toc; st=sprintf('%4.2f',t);

   subplot(2,2,i);
   pcolor(y,y,u), axis equal image; shading flat; colorbar
   title( sprintf('Grid: %dx%d, Steps: %d', N, N, Nstep));
   xlabel(['Time used: ', st,' sec.']);
   p=get(gca,'position'); p([3 4])=p([3 4])+0.03;
   set(gca,'position',p,'XTickLabel',[]);
end

repo: Aero-matlab


text:
function fresnel_cos_values_test ( )

%*****************************************************************************80
%
%% FRESNEL_COS_VALUES_TEST demonstrates the use of FRESNEL_COS_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FRESNEL_COS_VALUES_TEST:\n' );
  fprintf ( 1, '  FRESNEL_COS_VALUES stores values of\n' );
  fprintf ( 1, '  the Fresnel cosine integral C(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           C(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = fresnel_cos_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [charData] = srcCharLayerForward(W_rnn, W_emb, input, mask, charMap, vocabSize, params, isTest)
% Running char layer forward to compute word representations
% Input:
%   W_rnn: recurrent connections of multiple layers, e.g., W_rnn{ll}.
%   input: indices for the current batch
%   isTest: 1 -- don't store intermediate results in each state
%   isDecoder: 0 -- encoder, 1 -- decoder
% Output:
%   charData
%
% Thang Luong @ 2015, <lmthang@stanford.edu>

  charData.params = params;
  charData.params.numLayers = params.charNumLayers;

  % find rare words
  charData.rareFlags = input > params.srcCharShortList;
  charData.rareWords = unique(input(charData.rareFlags));

  % sample from frequent words
  if isTest == 0 && params.charSrcSample > 0
    % select by types: not masked and not unk
    freqWords = unique(input(~charData.rareFlags & mask & (input ~= params.srcUnk)));
    numSelect = floor(length(freqWords)*params.charSrcSample);
    perm = randperm(length(freqWords));
    selectFreqWords = freqWords(perm(1:numSelect));

    % assert
    if params.assert
      assert(isempty(intersect(charData.rareWords, selectFreqWords)) == 1);
      assert(ismember(params.srcSos, selectFreqWords) == 0);
    end

    % update rare words and flags
    charData.rareWords = union(charData.rareWords, selectFreqWords);
    charData.rareFlags = ismember(input, charData.rareWords);

    if params.assert
      assert(isempty(find(charData.rareWords == params.srcUnk, 1)));
    end
  end

  charData.numRareWords = length(charData.rareWords);
  charSeqs = charMap(charData.rareWords);
  seqLens = cellfun(@(x) length(x), charSeqs);

  if charData.numRareWords > 0
    charData.params.curBatchSize = charData.numRareWords;
    [charData.batch, charData.mask, charData.maxLen, charData.numSeqs] = leftPad(charMap(charData.rareWords), seqLens, params.srcCharSos, params.srcCharEos);

    charData.rnnFlags = struct('decode', 0, 'test', isTest, 'attn', 0, 'feedInput', 0, 'charSrcRep', 0, 'charTgtGen', 0);
    zeroState = createZeroState(charData.params);
    [charData.states, ~, ~] = rnnLayerForward(W_rnn, W_emb, zeroState, charData.batch, charData.mask, charData.params, charData.rnnFlags, [], [], []);
    charData.rareWordReps = charData.states{end}{end}.h_t;
    charData.rareWordMap = zeros(vocabSize, 1);
    charData.rareWordMap(charData.rareWords) = 1:charData.numRareWords;
  end
end

author: lmthang


text:
function neq = sys_init(params)

global N problem

% Here is a list of the different system information paramters.
% neq = 1  : number of state variables.
% neq = 2  : number of inputs.
% neq = 3  : number of parameters.
% neq = 4  : reserved.
% neq = 5  : reserved.
% neq = 6  : number of objective functions.
% neq = 7  : number of nonlinear trajectory constraints.
% neq = 8  : number of linear trajectory constraints.
% neq = 9  : number of nonlinear endpoint inequality constraints.
% neq = 10 : number of linear endpoint inequality constraints.
% neq = 11 : number of nonlinear endpoint equality constraints.
% neq = 12 : number of linear endpoint equality constraints.
% neq = 13 : 0 => nonlinear, 1 => linear, 2 => LTI, 3 => LQR, 4 => LQR and LTI.

% The default value is 0 for all except neq = 6 which defaults to 1.

% if params == [] then setup neq.  Otherwise the system parameters are
% getting passed.

if isempty(params),
    if problem == 1
        %LTI
        neq = [1, 2*N+1; 2, 1; 3, 4];
    elseif problem == 2
        %LTV
        neq = [1, 2*N+1; 2, 1; 3, 4];
    elseif problem == 3
        %Bang-bang
        neq = [1, 2*N+3; 2, 1; 3, 4; 12 2];
    end
else
    global sys_params
    sys_params = params;
end

author: Sable

repo: mcbench-benchmarks


text:
function [K, Knovar, argExp] = linardVardistPsi1Compute(linardkern, vardist, Z)

% LINARDVARDISTPSI1COMPUTE description.

% VARGPLVM

% variational means
N  = size(vardist.means,1);
%  inducing variables
M = size(Z,1);

A = rbfardKern.inverseWidth;

scales = sparse(diag(sqrt(kern.inputScales)));
x = x*scales;

if nargin < 3
  sk = x*x';
else
  x2 = x2*scales;
  sk = x*x2';
end
k = sk*kern.variance;

author: SheffieldML

repo: vargplvm


text:
function [u] = spm_u(a,df,STAT)
% uncorrected critical height threshold at a specified significance level
% FORMAT [u] = spm_u(a,df,STAT)
% a     - critical probability - {alpha}
% df    - [df{interest} df{error}]
% STAT  - Statistical field
%               'Z' - Gaussian field
%               'T' - T - field
%               'X' - Chi squared field
%               'F' - F - field
%               'P' - P - value
%
% u     - critical height {uncorrected}
%__________________________________________________________________________
%
% spm_u returns the uncorrected critical threshold at a specified
% significance
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_u.m 2690 2009-02-04 21:44:28Z guillaume $

if     STAT == 'Z'

    u   = spm_invNcdf(1 - a      );

elseif STAT == 'T'

    u   = spm_invTcdf(1 - a,df(2));

elseif STAT == 'X'

    u   = spm_invXcdf(1 - a,df(2));

elseif STAT == 'F'

    u   = spm_invFcdf(1 - a,df   );

elseif STAT == 'P'

    u   = a;

end

author: spm


text:
function [ft3] = oz2ft3(oz)
% Convert volume from US liquid ounces to cubic feet.
% Chad Greene 2012
ft3 = oz*0.0010443793403;


text:
function [datout, tim, Fnew] = ft_preproc_resample(dat, Fold, Fnew, method)

% FT_PREPROC_RESAMPLE resamples all channels in the data matrix
%
% Use as
%   dat = ft_preproc_resample(dat, Fold, Fnew, method)
% where
%   dat    = matrix with the input data (Nchans X Nsamples)
%   Fold   = scalar, original sampling frequency in Hz
%   Fnew   = scalar, desired sampling frequency in Hz
%   method = string, can be 'resample', 'decimate', 'downsample', 'fft'
%
% The resample method applies an anti-aliasing (lowpass) FIR filter to
% the data during the resampling process, and compensates for the filter's
% delay. For the other two methods you should apply an anti-aliassing
% filter prior to calling this function.
%
% If the data contains NaNs, these are ignored for the computation, but
% retained in the output.
%
% See also PREPROC, FT_PREPROC_LOWPASSFILTER

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

[nchans, nsamples] = size(dat);

if nargout>1
  tim = 1:nsamples;
  tim = ft_preproc_resample(tim, Fold, Fnew, method);
end

if Fold==Fnew
  ft_warning('requested sampling rate after resampling is the same as the sampling rate of the data, no resampling is performed');
  datout = dat;
  return
end

% resample and decimate require double formatted input
if ~strcmp(method, 'downsample')
  typ = class(dat);
  dat = cast(dat, 'double');
end

% preprocessing fails on channels that contain NaN
nanchan = any(isnan(dat),2);
if any(nanchan)
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
  if ft_platform_supports('matlabversion', '2020a', inf)
    % temporarily replace with zero, this is not needed for older MATLAB versions
    dat(nanchan,:) = 0;
  end
end

switch method
  case 'resample'
    [fold, fnew] = rat(Fold./Fnew); %account for non-integer fs
    Fnew         = Fold.*(fnew./fold); %get new fs exact

    % the actual implementation resamples along columns
    datout = resample(dat', fnew, fold)';

  case 'decimate'
    fac         = round(Fold/Fnew);
    % this only works one channel at the time
    nresampled  = ceil(nsamples/fac);
    datout      = zeros(nchans, nresampled);
    for i=1:nchans
      datout(i,:) = decimate(dat(i,:), fac);
    end

  case 'downsample'
    fac = round(Fold/Fnew);
    % the actual implementation resamples along columns
    datout = downsample(dat', fac)';

  case 'fft'
    % Code written for SPM by Jean Daunizeau
    fac         = Fnew/Fold;
    nresampled  = floor(nsamples*fac);
    fac         = nresampled/nsamples;
    datfft      = fftshift(fft(dat,[],2),2);
    middle = floor(size(datfft,2)./2)+1;
    if fac>1 % upsample
      npad = floor((nresampled-nsamples)./2);

      if nsamples/2 == floor(nsamples/2)
        datfft(:,1) = []; % throw away non symmetric DFT coef
      end

      datfft  = [zeros(size(datfft,1),npad), datfft,zeros(size(datfft,1),npad)];
    else % downsample
      ncut    = floor(nresampled./2);
      datfft  = datfft(:,middle-ncut:middle+ncut);
    end
    datout      = fac*ifft(ifftshift(datfft,2),[],2);
  otherwise
    ft_error('unsupported resampling method');
end

if any(nanchan) && ft_platform_supports('matlabversion', '2020a', inf)
  % replace the zeros back to nan, this is not needed for older MATLAB versions
  datout(nanchan,:) = nan;
end

if ~strcmp(method, 'downsample')
  % convert back into the original input format
  datout = cast(datout, typ);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
% [fh, varfh] = gppred(X, f, Covf, Xh, logtheta, covfunc, varargin)
%
function [fh, varfh] = gppred(X, f, Covf, Xh, logtheta, covfunc, ...
                              varargin)

warning('Deprecated. Use gp_predict or gp_predict_pseudo.');

% If you need more than just variances for PCA, you could give some sparse
% matrix to indicate which covariances should be evaluated. But again, this
% covariances should be zero for the latent function values fh a
% priori.. Things would be much easier if I just factored the principal
% components too.. :)

if ~iscell(covfunc)
  covfunc = {covfunc};
end

opts = struct('cholkx', [], 'khx', [], 'kh', []);

[ opts, errmsg, wrnmsg ] = argschk( opts, varargin{:} );
if ~isempty(errmsg), error( errmsg ), end
if ~isempty(wrnmsg), warning( wrnmsg ), end

nh = cols(Xh);

I = diag(sparse(ones(cols(X),1)));

if isempty(opts.cholkx)
  Kpp = feval(covfunc{:}, logtheta, X, X);

  % DEBUG: REGULARIZE
  %Kpp = regularize(Kpp);

  Kpp = regularize(Kpp);
  Lp = chol(Kpp, 'lower');
%  [Lp,Kpp] = safechol(Kpp, 1e-10, 'lower');
  % DEBUG:
  %Kpp = Lp*Lp';
else
  Lp = opts.cholkx;
end

if isempty(opts.khx)
  Kxp = feval(covfunc{:}, logtheta, Xh, X);
else
  Kxp = opts.khx;
end

if isempty(opts.kh)
  Kx = feval(covfunc{:}, logtheta, Xh, []);
else
  Kx = opts.kh;
end

fh = Kxp * solve_triu(Lp', solve_tril(Lp, f));

if nargout >= 2
  varfh = zeros(nh,1);
  for i=1:nh
    r = solve_tril(Lp, Kxp(i,:)');
    s = solve_triu(Lp', r);
    varfh(i) = Kx(i) - r'*r + s'*Covf*s;
  end
end

% $$$ if isequal(covfunc{1}, @gpcovPP)
% $$$ %if strcmpi(covfunc{1}, @gpcovPP)
% $$$   if ~issparse(Lp)
% $$$     error('Lp should be sparse');
% $$$   else
% $$$     density_in_gppred = nnz(Lp) / prod(size(Lp))
% $$$   end
% $$$ end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function fem2d_project ( sample_prefix, fem_prefix )

%*****************************************************************************80
%
%% MAIN is the main program for FEM2D_PROJECT.
%
%  Discussion:
%
%    FEM2D_PROJECT reads files defining a sampling of a (scalar or vector)
%    function of 2 arguments, and a list of nodes and triangular elements
%    to use for a finite element representation of the data.
%
%    It computes a set of finite element coefficients to be associated with
%    the given finite element mesh, and writes that information to a file
%    so that an FEM representation is formed by the node, element and value
%    files.
%
%  Usage:
%
%    fem2d_project ( 'sample_prefix', 'fem_prefix' )
%
%    where 'sample_prefix' is the common prefix for the SAMPLE files:
%
%    * sample_prefix_nodes.txt,     the node coordinates where samples were taken;
%    * sample_prefix_elements.txt,  the nodes that make up each element;
%    * sample_prefix_values.txt,    the sample values.
%
%    and 'fem_prefix' is the common prefix for the FEM files:
%
%    * fem_prefix_nodes.txt,    the node coordinates;
%    * fem_prefix_elements.txt, the nodes that make up each element;
%    * fem_prefix_values.txt,   the values defined at each node,
%                               (computed by this program).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string SAMPLE_PREFIX, the common prefix for sample files.
%
%    Input, string FEM_PREFIX, the common prefix for FEM files.
%
  timestamp ( );

  fprintf ( '\n' );
  fprintf ( 'FEM2D_PROJECT\n' );
  fprintf ( '  MATLAB version.\n' );
  fprintf ( '\n' );
  fprintf ( '  Read files defining a sampling of a function of 2 arguments.\n' );
  fprintf ( '  Read files defining a finite element mesh.\n' );
  fprintf ( '  Project the sample data onto the mesh, and\n' );
  fprintf ( '  write a file of FEM coefficient values.\n' );
%
%  Get the number of command line arguments.
%
  if ( nargin < 1 )

    fprintf ( '\n' );
    sample_prefix = input ( 'Enter the sample file prefix:  ' );

  else

  end

  if ( nargin < 2 )

    fprintf ( '\n' );
    fem_prefix = input ( 'Enter the FEM file prefix:  ' );

  else

  end
%
%  Create the filenames.
%
  sample_node_filename    = strcat ( sample_prefix, '_nodes.txt' );
  sample_element_filename = strcat ( sample_prefix, '_elements.txt' );
  sample_value_filename   = strcat ( sample_prefix, '_values.txt' );

  fem_node_filename       = strcat ( fem_prefix, '_nodes.txt' );
  fem_element_filename    = strcat ( fem_prefix, '_elements.txt' );
  fem_value_filename      = strcat ( fem_prefix, '_values.txt' );
%
%  Read the SAMPLE NODE, ELEMENT and VALUE data.
%
  sample_node_xy = load ( sample_node_filename );
  sample_node_xy = sample_node_xy';
  [sample_node_dim, sample_node_num ] = size ( sample_node_xy );

  fprintf ( '\n' );
  fprintf ( '  Sample node spatial dimension is %d\n', sample_node_dim );
  fprintf ( '  Sample node number is            %d\n', sample_node_num );

  if ( sample_node_dim ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'FEM2D_PROJECT - Fatal error!\n' );
    fprintf ( '  Spatial dimension of the sample nodes is not 2.\n' );
    return
  end

  sample_element_node = load ( sample_element_filename );
  sample_element_node = sample_element_node';
  [ sample_element_order, sample_element_num ] = size ( sample_element_node );

  fprintf ( '\n' );
  fprintf ( '  Sample element order is  %d\n', sample_element_order );
  fprintf ( '  Sample element number is %d\n', sample_element_num );

  if ( sample_element_order ~= 3 )
    fprintf ( '\n' );
    fprintf ( 'FEM2D_PROJECT - Fatal error!\n' );
    fprintf ( '  The sample elements must be of order 3.\n' );
    return
  end

  sample_value = load ( sample_value_filename );
  sample_value = sample_value';
  [ sample_value_dim, sample_value_num ] = size ( sample_value );

  fprintf ( '\n' );
  fprintf ( '  Sample value dimension is        %d\n', sample_value_dim );
  fprintf ( '  Sample value number is           %d\n', sample_value_num );

  if ( sample_value_num ~= sample_node_num )
    fprintf ( '\n' );
    fprintf ( 'FEM2D_PROJECT - Fatal error!\n' );
    fprintf ( '  Number of sample nodes and values are not equal.\n' );
    return
  end

  sample_value_min = min ( sample_value, [], 2 );
  sample_value_max = max ( sample_value, [], 2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Index      SAMPLE Min      SAMPLE Max\n' );
  fprintf ( 1, '\n' );
  for i = 1 : sample_value_dim
    fprintf ( 1, '  %8d  %14f  %14f\n', i, sample_value_min(i), sample_value_max(i) );
  end
%
%  Create the sample element neighbor array.
%
  sample_element_neighbor = triangulation_order3_neighbor_triangles ( ...
    sample_element_num, sample_element_node );

  fprintf ( '\n' );
  fprintf ( '  The element neighbor array has been computed.\n' );
%
%  Read the FEM NODE and ELEMENT data.
%
  fem_node_xy = load ( fem_node_filename );
  fem_node_xy = fem_node_xy';
  [ fem_node_dim, fem_node_num ] = size ( fem_node_xy );

  fprintf ( '\n' );
  fprintf ( '  The FEM node dimension is        %d\n', fem_node_dim );
  fprintf ( '  The FEM node number is           %d\n', fem_node_num );

  if ( fem_node_dim ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'FEM2D_PROJECT - Fatal error!\n' );
    fprintf ( '  Spatial dimension of the nodes is not 2.\n' );
    return
  end

  fem_element_node = load ( fem_element_filename );
  fem_element_node = fem_element_node';
  [ fem_element_order, fem_element_num ] = size ( fem_element_node );

  fprintf ( '  The FEM element order is         %d\n', fem_element_order );
  fprintf ( '  The FEM element number is        %d\n', fem_element_num );

  if ( fem_element_order ~= 3 )
    fprintf ( '\n' );
    fprintf ( 'FEM2D_PROJECT - Fatal error!\n' );
    fprintf ( '  The FEM elements must be of order 3.\n' );
    return
  end
%
%  Compute the FEM values.
%
  fem_value_dim = sample_value_dim;
  fem_value_num = fem_node_num;

  fem_value = fem2d_transfer ( sample_element_order, ...
    sample_element_num, sample_value_dim, ...
    sample_node_xy, sample_element_node, sample_element_neighbor, sample_value, ...
    fem_node_num, fem_element_num, fem_value_dim, fem_node_xy, fem_element_node );
%
%  Print Min, Max.
%
  fem_value_min = min ( fem_value, [], 2 );
  fem_value_max = max ( fem_value, [], 2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Index         FEM Min         FEM Max\n' );
  fprintf ( 1, '\n' );
  for i = 1 : fem_value_dim
    fprintf ( 1, '  %8d  %14f  %14f\n', i, fem_value_min(i), fem_value_max(i) );
  end
%
%  Write the FEM values.
%
  r8mat_write ( fem_value_filename, fem_value_dim, fem_value_num, fem_value );

  fprintf ( '\n' );
  fprintf ( '  FEM value data written to "%s".\n',  fem_value_filename );
%
%  Terminate.
%
  fprintf ( '\n' );
  fprintf ( 'FEM2D_PROJECT\n' );
  fprintf ( '  Normal end of execution.\n' );

  fprintf ( '\n' );
  timestamp ( );

  return
end
function [ phi, dphidx, dphidy ] = basis_mn_t3 ( t, n, p )

%*****************************************************************************80
%
%% BASIS_MN_T3: all bases functions at N points for a T3 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices of a triangle.
%    It works directly with these coordinates, and does not refer to a
%    reference element.
%
%    The sides of the triangle DO NOT have to lie along a coordinate
%    axis.
%
%    The routine evaluates the basis functions associated with each vertex,
%    and their derivatives with respect to X and Y.
%
%  Physical Element T3:
%
%            3
%           / \
%          /   \
%         /     \
%        /       \
%       1---------2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the vertices of the triangle.  It is common to list
%    these points in counter clockwise order.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real P(2,N), the coordinates of the evaluation points.
%
%    Output, real PHI(3,N), the basis functions at the evaluation points.
%
%    Output, real DPHIDX(3,N), DPHIDY(3,N), the basis derivatives
%    at the evaluation points.
%
%  Local parameters:
%
%    Local, real AREA, is (twice) the area of the triangle.
%
  area = t(1,1) * ( t(2,2) - t(2,3) ) ...
       + t(1,2) * ( t(2,3) - t(2,1) ) ...
       + t(1,3) * ( t(2,1) - t(2,2) );

  phi(1,1:n) =     (  ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) )     ...
                    - ( t(2,3) - t(2,2) ) * ( p(1,1:n) - t(1,2) ) );
  dphidx(1,1:n) =   - ( t(2,3) - t(2,2) );
  dphidy(1,1:n) =     ( t(1,3) - t(1,2) );

  phi(2,1:n) =     (  ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) )     ...
                    - ( t(2,1) - t(2,3) ) * ( p(1,1:n) - t(1,3) ) );
  dphidx(2,1:n) =   - ( t(2,1) - t(2,3) );
  dphidy(2,1:n) =     ( t(1,1) - t(1,3) );

  phi(3,1:n) =     (  ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) )     ...
                    - ( t(2,2) - t(2,1) ) * ( p(1,1:n) - t(1,1) ) );
  dphidx(3,1:n) =   - ( t(2,2) - t(2,1) );
  dphidy(3,1:n) =     ( t(1,2) - t(1,1) );
%
%  Normalize.
%
  phi(1:3,1:n)    = phi(1:3,1:n) / area;
  dphidx(1:3,1:n) = dphidx(1:3,1:n) / area;
  dphidy(1:3,1:n) = dphidy(1:3,1:n) / area;

  return
end
function [ phi, dphidx, dphidy ] = basis_mn_t6 ( t, n, p )

%*****************************************************************************80
%
%% BASIS_MN_T6: all bases for N points in a T6 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices and midside
%    nodes of a triangle.  It works directly with these coordinates, and does
%    not refer to a reference element.
%
%    This routine requires that the midside nodes be "in line"
%    with the vertices, that is, that the sides of the triangle be
%    straight.  However, the midside nodes do not actually have to
%    be halfway along the side of the triangle.
%
%  Physical element T6:
%
%    This picture indicates the assumed ordering of the six nodes
%    of the triangle.
%
%             3
%            / \
%           /   \
%          6     5
%         /       \
%        /         \
%       1-----4-----2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,6), the nodal oordinates of the element.
%    It is common to list these points in counter clockwise order.
%
%    Input, real P(2,N), the evaluation points.
%
%    Output, real PHI(6,N), the basis functions at the evaluation points.
%
%    Output, real DPHIDX(6,N), DPHIDY(6,N), the basis derivatives at the
%    evaluation points.
%
%  Local Parameters:
%
%    Local, real AREA, is (twice) the area of the triangle.
%

%
%  Basis function 1: PHI(X,Y) = G(3,2) * H(6,4) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3)   - t(2,2) ) ...
          - ( t(1,3)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  gn(1:n) = ( t(1,1) - t(1,2) ) * ( t(2,3) - t(2,2) ) ...
          - ( t(1,3) - t(1,2) ) * ( t(2,1) - t(2,2) );

  hx(1:n) = ( p(1,1:n) - t(1,4) ) * ( t(2,6)   - t(2,4) ) ...
          - ( t(1,6)   - t(1,4) ) * ( p(2,1:n) - t(2,4) );

  hn(1:n) = ( t(1,1) - t(1,4) ) * ( t(2,6) - t(2,4) ) ...
          - ( t(1,6) - t(1,4) ) * ( t(2,1) - t(2,4) );

  phi(1,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(1,1:n) =  ( ( t(2,3) - t(2,2) ) * hx(1:n) ...
                    + gx(1:n) * ( t(2,6) - t(2,4) ) ) ...
                  ./ ( gn(1:n) .* hn(1:n) );

  dphidy(1,1:n) = -( ( t(1,3) - t(1,2) ) * hx(1:n) ...
                   + gx * ( t(1,6) - t(1,4) ) ) ...
                  ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 2: PHI(X,Y) = G(3,1) * H(4,5) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3)   - t(2,1) ) ...
          - ( t(1,3)   - t(1,1) ) * ( p(2,1:n) - t(2,1) );

  gn(1:n) = ( t(1,2) - t(1,1) ) * ( t(2,3) - t(2,1) ) ...
          - ( t(1,3) - t(1,1) ) * ( t(2,2) - t(2,1) );

  hx(1:n) = ( p(1,1:n) - t(1,5) ) * ( t(2,4)   - t(2,5) ) ...
          - ( t(1,4)   - t(1,5) ) * ( p(2,1:n) - t(2,5) );

  hn(1:n) = ( t(1,2) - t(1,5) ) * ( t(2,4) - t(2,5) ) ...
          - ( t(1,4) - t(1,5) ) * ( t(2,2) - t(2,5) );

  phi(2,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(2,1:n) =  ( ( t(2,3) - t(2,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,4) - t(2,5) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(2,1:n) = -( ( t(1,3) - t(1,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,4) - t(1,5) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 3: PHI(X,Y) = G(1,2) * H(5,6) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1)   - t(2,2) ) ...
          - ( t(1,1)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  gn(1:n) = ( t(1,3) - t(1,2) ) * ( t(2,1) - t(2,2) ) ...
          - ( t(1,1) - t(1,2) ) * ( t(2,3) - t(2,2) );

  hx(1:n) = ( p(1,1:n) - t(1,6) ) * ( t(2,5)   - t(2,6) ) ...
          - ( t(1,5)   - t(1,6) ) * ( p(2,1:n) - t(2,6) );

  hn(1:n) = ( t(1,3) - t(1,6) ) * ( t(2,5) - t(2,6) ) ...
          - ( t(1,5) - t(1,6) ) * ( t(2,3) - t(2,6) );

  phi(3,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(3,1:n) =  ( ( t(2,1) - t(2,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,5) - t(2,6) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(3,1:n) = -( ( t(1,1) - t(1,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,5) - t(1,6) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 4: PHI(X,Y) = G(1,3) * H(2,3) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,1)   - t(2,3) ) ...
          - ( t(1,1)   - t(1,3) ) * ( p(2,1:n) - t(2,3) );

  gn(1:n) = ( t(1,4) - t(1,3) ) * ( t(2,1) - t(2,3) ) ...
          - ( t(1,1) - t(1,3) ) * ( t(2,4) - t(2,3) );

  hx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,2)   - t(2,3) ) ...
          - ( t(1,2)   - t(1,3) ) * ( p(2,1:n) - t(2,3) );

  hn(1:n) = ( t(1,4) - t(1,3) ) * ( t(2,2) - t(2,3) ) ...
          - ( t(1,2) - t(1,3) ) * ( t(2,4) - t(2,3) );

  phi(4,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(4,1:n) =  ( ( t(2,1) - t(2,3) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,2) - t(2,3) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(4,1:n) = -( ( t(1,1) - t(1,3) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,2) - t(1,3) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 5: PHI(X,Y) = G(2,1) * H(3,1) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,2)   - t(2,1) ) ...
          - ( t(1,2)   - t(1,1) ) * ( p(2,1:n) - t(2,1) );

  gn(1:n) = ( t(1,5) - t(1,1) ) * ( t(2,2) - t(2,1) ) ...
          - ( t(1,2) - t(1,1) ) * ( t(2,5) - t(2,1) );

  hx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3)   - t(2,1) ) ...
          - ( t(1,3)   - t(1,1) ) * ( p(2,1:n) - t(2,1) );

  hn(1:n) = ( t(1,5) - t(1,1) ) * ( t(2,3) - t(2,1) ) ...
          - ( t(1,3) - t(1,1) ) * ( t(2,5) - t(2,1) );

  phi(5,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(5,1:n) =  ( ( t(2,2) - t(2,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,3) - t(2,1) ) ) ...
               ./ ( gn(1:n) .* hn(1:n) );

  dphidy(5,1:n) = -( ( t(1,2) - t(1,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,3) - t(1,1) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 6: PHI(X,Y) = G(1,2) * H(3,2) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1)   - t(2,2) ) ...
          - ( t(1,1)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  gn(1:n) = ( t(1,6) - t(1,2) ) * ( t(2,1) - t(2,2) ) ...
          - ( t(1,1) - t(1,2) ) * ( t(2,6) - t(2,2) );

  hx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3)   - t(2,2) ) ...
          - ( t(1,3)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  hn(1:n) = ( t(1,6) - t(1,2) ) * ( t(2,3) - t(2,2) ) ...
          - ( t(1,3) - t(1,2) ) * ( t(2,6) - t(2,2) );

  phi(6,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(6,1:n) =  ( ( t(2,1) - t(2,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,3) - t(2,2) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(6,1:n) = -( ( t(1,1) - t(1,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,3) - t(1,2) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  return
end
function fem_value = fem2d_transfer ( sample_element_order, ...
  sample_element_num, sample_value_dim, sample_node_xy, ...
  sample_element_node, sample_element_neighbor, sample_value, fem_node_num, ...
  fem_element_num, fem_value_dim, fem_node_xy, fem_element_node )

%*****************************************************************************80
%
%% FEM2D_TRANSFER "transfers" from one finite element mesh to another.
%
%  Discussion:
%
%    1) the linear system A*X=B is defined with A being a full storage matrix.
%       This can easily be fixed using MATLAB's SPARSE facility.
%
%    2) the quadrature rule used is low order.
%
%    3) the triangular elements are assumed to be linear.
%
%
%    We are also given a set of "sample" finite element function defined
%    by SAMPLE_NODE_XY, SAMPLE_ELEMENT, and SAMPLE_VALUE.
%
%    We are given a second finite element mesh, FEM_NODE_XY and
%    FEM_ELEMENT_NODE.
%
%    Our aim is to "project" the sample data values into the finite element
%    space, that is, to come up with a finite element function FEM_VALUE which
%    well approximates the sample data.
%
%    Now let W(x,y) represent a function interpolating the sample data, and
%    let Vk(x,y) represent the finite element basis function associated with
%    node K.
%
%    Then we seek the coefficient vector U corresponding to a finite element
%    function U(x,y) of the form:
%
%      U(x,y) = sum ( 1 <= K <= N ) Uk * Vk(x,y)
%
%    To determine the coefficent vector entries U, we form a set of
%    projection equations.  For node K at grid point (I,J), the associated
%    basis function Vk(x,y) is used to pose the equation:
%
%      Integral U(x,y) Vk(x,y) dx dy = Integral W(x,y) Vk(x,y) dx dy
%
%    The left hand side is the usual stiffness matrix times the desired
%    coefficient vector U.  To complete the system, we simply need to
%    determine the right hand side, that is, the integral of the data function
%    W against the basis function Vk.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SAMPLE_NODE_NUM, the number of nodes.
%
%    Input, integer SAMPLE_ELEMENT_ORDER, the element order.
%
%    Input, integer SAMPLE_ELEMENT_NUM, the number of elements.
%
%    Input, integer SAMPLE_VALUE_DIM, the value dimension.
%
%    Input, integer SAMPLE_VALUE_NUM, the number of values.
%
%    Input, real SAMPLE_NODE_XY(2,SAMPLE_NODE_NUM), the nodes.
%
%    Input, integer SAMPLE_ELEMENT_NODE(SAMPLE_ELEMENT_ORDER,SAMPLE_ELEMENT_NUM),
%    the nodes that make up each element.
%
%    Input, integer SAMPLE_ELEMENT_NEIGHBOR(3,SAMPLE_ELEMENT_NUM),
%    the neighbor triangles.
%
%    Input, real SAMPLE_VALUE(SAMPLE_VALUE_DIM,SAMPLE_NODE_NUM),
%    the values.
%
%    Input, integer FEM_NODE_NUM, the number of nodes.
%
%    Input, integer FEM_ELEMENT_ORDER, the element order.
%
%    Input, integer FEM_ELEMENT_NUM, the number of elements.
%
%    Input, integer FEM_VALUE_DIM, the value dimension.
%
%    Input, integer FEM_VALUE_NUM, the number of values.
%
%    Input, real FEM_NODE_XY(2,FEM_NODE_NUM), the nodes.
%
%    Input, integer FEM_ELEMENT_NODE(FEM_ELEMENT_ORDER,FEM_ELEMENT_NUM),
%    the nodes that make up each element.
%
%    Output, real FEM_VALUE(FEM_VALUE_DIM,FEM_VALUE_NUM),
%    the values.
%
  project_node_num = 1;
  quad_num = 3;
%
%  Assemble the coefficient matrix A and the right-hand side B.
%
  b = zeros (fem_node_num,fem_value_dim);
%
%  Define A as a sparse matrix.
%
  a = sparse ( [], [], [], fem_node_num, fem_node_num );

% a = zeros (fem_node_num,fem_node_num);

  for element = 1 : fem_element_num

    i1 = fem_element_node(1,element);
    i2 = fem_element_node(2,element);
    i3 = fem_element_node(3,element);

    area = 0.5 * ...
      ( fem_node_xy(1,i1) * ( fem_node_xy(2,i2) - fem_node_xy(2,i3) ) ...
      + fem_node_xy(1,i2) * ( fem_node_xy(2,i3) - fem_node_xy(2,i1) ) ...
      + fem_node_xy(1,i3) * ( fem_node_xy(2,i1) - fem_node_xy(2,i2) ) );
%
%  Consider each quadrature point.
%  Here, we use the midside nodes as quadrature points.
%
    for quad = 1 : quad_num

      q1 =       quad;
      q2 = mod ( quad, quad_num ) + 1;

      nq1 = fem_element_node(q1,element);
      nq2 = fem_element_node(q2,element);

      xq = 0.5 * ( fem_node_xy(1,nq1) + fem_node_xy(1,nq2) );
      yq = 0.5 * ( fem_node_xy(2,nq1) + fem_node_xy(2,nq2) );
      wq = 1.0 / 3.0;
%
%  Consider each test function in the element.
%
      for ti1 = 1 : 3

        ti2 = mod ( ti1,     3 ) + 1;
        ti3 = mod ( ti1 + 1, 3 ) + 1;

        nti1 = fem_element_node(ti1,element);
        nti2 = fem_element_node(ti2,element);
        nti3 = fem_element_node(ti3,element);

        qi = 0.5 * ( ...
            ( fem_node_xy(1,nti3) - fem_node_xy(1,nti2) ) ...
          * ( yq - fem_node_xy(2,nti2) ) ...
          - ( fem_node_xy(2,nti3) - fem_node_xy(2,nti2) ) ...
          * ( xq - fem_node_xy(1,nti2) ) ) / area;
%
%  The projection takes place here.  The finite element code needs the value
%  of the sample function at the point (XQ,YQ).  The call to PROJECTION
%  locates (XQ,YQ) in the triangulated mesh of sample data, and returns a
%  value produced by piecewise linear interpolation.
%
        project_node_xy(1,1) = xq;
        project_node_xy(2,1) = yq;

        project_value = projection ( sample_node_xy, ...
          sample_element_order, sample_element_num, sample_element_node, ...
          sample_element_neighbor, sample_value_dim, sample_value, ...
          project_node_num, project_node_xy );

        b(nti1,1:fem_value_dim) = b(nti1,1:fem_value_dim) ...
          + area * wq * ( project_value(1:fem_value_dim,1)' * qi );
%
%  Consider each basis function in the element.
%
        for tj1 = 1 : 3

          tj2 = mod ( tj1,     3 ) + 1;
          tj3 = mod ( tj1 + 1, 3 ) + 1;

          ntj1 = fem_element_node(tj1,element);
          ntj2 = fem_element_node(tj2,element);
          ntj3 = fem_element_node(tj3,element);

          qj = 0.5 * ( ...
              ( fem_node_xy(1,ntj3) - fem_node_xy(1,ntj2) ) ...
            * ( yq - fem_node_xy(2,ntj2) ) ...
            - ( fem_node_xy(2,ntj3) - fem_node_xy(2,ntj2) ) ...
            * ( xq - fem_node_xy(1,ntj2) ) ) / area;

          a(nti1,ntj1) = a(nti1,ntj1) + area * wq * ( qi * qj );

        end

      end

    end

  end
%
%  Solve the linear system A * X = B.
%
  fem_value = ( a \ b )';

  return
end
function isgn = i4col_compare ( m, n, a, i, j )

%*****************************************************************************80
%
%% I4COL_COMPARE compares columns I and J of a integer array.
%
%  Example:
%
%    Input:
%
%      M = 3, N = 4, I = 2, J = 4
%
%      A = (
%        1  2  3  4
%        5  6  7  8
%        9 10 11 12 )
%
%    Output:
%
%      ISGN = -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an array of N columns of vectors of length M.
%
%    Input, integer I, J, the columns to be compared.
%    I and J must be between 1 and N.
%
%    Output, integer ISGN, the results of the comparison:
%    -1, column I < column J,
%     0, column I = column J,
%    +1, column J < column I.
%

%
%  Check.
%
  if ( i < 1)
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Column index I = %d < 1.\n', i );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  if ( n < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  N = %d < column index I = %d.\n', n, i );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  if ( j < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Column index J = %d < 1.\n', j );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  if ( n < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  N = %d < column index J = %d.\n', n, j );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  isgn = 0;

  if ( i == j )
    return
  end

  k = 1;

  while ( k <= m )

    if ( a(k,i) < a(k,j) )
      isgn = -1;
      return
    elseif ( a(k,j) < a(k,i) )
      isgn = +1;
      return
    end

    k = k + 1;

  end

  return
end
function a = i4col_sort_a ( m, n, a )

%*****************************************************************************80
%
%% I4COL_SORT_A ascending sorts an I4COL.
%
%  Discussion:
%
%    In lexicographic order, the statement "X < Y", applied to two real
%    vectors X and Y of length M, means that there is some index I, with
%    1 <= I <= M, with the property that
%
%      X(J) = Y(J) for J < I,
%    and
%      X(I) < Y(I).
%
%    In other words, the first time they differ, X is smaller.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of A, and the length of
%    a vector of data.
%
%    Input, integer N, the number of columns of A.
%
%    Input, integer A(M,N), the array of N columns of M-vectors.
%
%    Output, integer A(M,N), the columns of A have been sorted in ascending
%    lexicographic order.
%
  if ( m <= 0 )
    return
  end

  if ( n <= 1 )
    return
  end
%
%  Initialize.
%
  indx = 0;
  isgn = 0;
%
%  Call the external heap sorter.
%
  while ( 1 )

    [ indx, i, j ] = sort_heap_external ( n, indx, isgn );
%
%  Interchange the I and J objects.
%
    if ( 0 < indx )

      a = i4col_swap ( m, n, a, i, j );
%
%  Compare the I and J objects.
%
    elseif ( indx < 0 )

      isgn = i4col_compare ( m, n, a, i, j );

    elseif ( indx == 0 )

      break

    end

  end

  return
end
function a = i4col_swap ( m, n, a, i, j )

%*****************************************************************************80
%
%% I4COL_SWAP swaps columns I and J of a integer array of column data.
%
%  Example:
%
%    Input:
%
%      M = 3, N = 4, I = 2, J = 4
%
%      A = (
%        1  2  3  4
%        5  6  7  8
%        9 10 11 12 )
%
%    Output:
%
%      A = (
%        1  4  3  2
%        5  8  7  6
%        9 12 11 10 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer A(M,N), an array of N columns of length M.
%
%    Input, integer I, J, the columns to be swapped.
%
%    Output, integer A(M,N), the array, with columns I and J swapped.
%
  if ( i < 1 || n < i || j < 1 || n < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_SWAP - Fatal error!\n' );
    fprintf ( 1, '  I or J is out of bounds.\n' );
    fprintf ( 1, '  I =    %d\n', i );
    fprintf ( 1, '  J =    %d\n', j );
    fprintf ( 1, '  N =    %d\n', n );
    error ( 'I4COL_SWAP - Fatal error!' );
  end

  if ( i == j )
    return
  end

  col(1:m) = a(1:m,i)';
  a(1:m,i) = a(1:m,j);
  a(1:m,j) = col(1:m)';

  return
end
function sample_value = projection ( fem_node_xy, ...
  fem_element_order, fem_element_num, fem_element_node, ...
  fem_element_neighbor, fem_value_dim, fem_value, sample_node_num, ...
  sample_node_xy )

%*****************************************************************************80
%
%% PROJECTION evaluates an FEM function on a T3 or T6 triangulation.
%
%  Discussion:
%
%    Note that the sample values returned are true values of the underlying
%    finite element function.  They are NOT produced by constructing some
%    other function that interpolates the data at the finite element nodes
%    (something which MATLAB's griddata function can easily do.)  Instead,
%    each sampling node is located within one of the associated finite
%    element triangles, and the finite element function is developed and
%    evaluated there.
%
%    MATLAB's scattered data interpolation is wonderful, but it cannot
%    be guaranteed to reproduce the finite element function corresponding
%    to nodal data.  This routine can (or at least tries to%).
%
%    So if you are using finite elements, then using THIS routine
%    (but not MATLAB's griddata function), what you see is what you have%
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FEM_NODE_NUM, the number of nodes.
%
%    Input, real FEM_NODE_XY(2,FEM_NODE_NUM), the coordinates
%    of the nodes.
%
%    Input, integer FEM_ELEMENT_ORDER, the order of the elements,
%    either 3 or 6.
%
%    Input, integer FEM_ELEMENT_NUM, the number of triangles.
%
%    Input, integer FEM_ELEMENT_NODE(FEM_ELEMENT_ORDER,FEM_ELEMENT_NUM), the
%    nodes that make up each triangle.
%
%    Input, integer FEM_ELEMENT_NEIGHBOR(3,FEM_ELEMENT_NUM), the
%    index of the neighboring triangle on each side, or -1 if no neighbor there.
%
%    Input, integer FEM_VALUE_DIM, the "dimension" of the values.
%
%    Input, real FEM_VALUE(FEM_VALUE_DIM,FEM_NODE_NUM), the
%    finite element coefficient values at each node.
%
%    Input, integer SAMPLE_NODE_NUM, the number of sample nodes.
%
%    Input, real SAMPLE_NODE_XY(2,SAMPLE_NODE_NUM), the sample nodes.
%
%    Output, real SAMPLE_VALUE(FEM_VALUE_DIM,SAMPLE_NODE_NUM),
%    the sampled values.
%
  sample_value = zeros(fem_value_dim,sample_node_num);
%
%  For each sample point: find the triangle T that contains it,
%  and evaluate the finite element function there.
%
  for j = 1 : sample_node_num

    p_xy(1:2,1) = sample_node_xy(1:2,j);
%
%  Find the triangle T that contains the point.
%
    t = triangulation_search_delaunay ( ...
      fem_node_xy, fem_element_num, fem_element_node, ...
      fem_element_neighbor, p_xy );

    if ( t == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PROJECTION - Fatal error!\n' );
      fprintf ( 1, '  Triangulation search failed!\n' );
      error ( 'PROJECTION - Fatal error!' );
    end
%
%  Evaluate the finite element basis functions at the point in T.
%
    t_node(1:fem_element_order) = fem_element_node(1:fem_element_order,t);

    t_xy(1:2,1:fem_element_order) = fem_node_xy(1:2,t_node);

    if ( fem_element_order == 3 )
      b = basis_mn_t3 ( t_xy, 1, p_xy );
    elseif ( fem_element_order == 6 )
      b = basis_mn_t6 ( t_xy, 1, p_xy );
    end
%
%  Multiply by the finite element values to get the sample values.
%
    for i = 1 : fem_value_dim
      sample_value(i,j) = fem_value(i,t_node(1:fem_element_order)) ...
        * b(1:fem_element_order);
    end

  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For greater precision, try:
%
%     fprintf ( output_unit, '  %24,16f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %14f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ indx, i, j ] = sort_heap_external ( n, indx, isgn )

%*****************************************************************************80
%
%% SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.
%
%  Discussion:
%
%    The actual list of data is not passed to the routine.  Hence this
%    routine may be used to sort integers, reals, numbers, names,
%    dates, shoe sizes, and so on.  After each call, the routine asks
%    the user to compare or interchange two items, until a special
%    return value signals that the sorting is completed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2004
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf.
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of items to be sorted.
%
%    Input, integer INDX, the main communication signal.
%    The user must set INDX to 0 before the first call.
%    Thereafter, the user should set the input value of INDX
%    to the output value from the previous call.
%
%    Input, integer ISGN, results of comparison of elements I and J.
%    (Used only when the previous call returned INDX less than 0).
%    ISGN <= 0 means I is less than or equal to J;
%    0 <= ISGN means I is greater than or equal to J.
%
%    Output, integer INDX, the main communication signal.
%    If INDX is
%
%      greater than 0, the user should:
%      * interchange items I and J;
%      * call again.
%
%      less than 0, the user should:
%      * compare items I and J;
%      * set ISGN = -1 if I < J, ISGN = +1 if J < I;
%      * call again.
%
%      equal to 0, the sorting is done.
%
%    Output, integer I, J, the indices of two items.
%    On return with INDX positive, elements I and J should be interchanged.
%    On return with INDX negative, elements I and J should be compared, and
%    the result reported in ISGN on the next call.
%
  persistent i_save;
  persistent j_save;
  persistent k;
  persistent k1;
  persistent n1;

  if ( isempty ( i_save ) )
    i_save = -1;
  end

  if ( isempty ( j_save ) )
    j_save = -1;
  end
%
%  INDX = 0: This is the first call.
%
  if ( indx == 0 )

    k = floor ( n / 2 );
    k1 = k;
    n1 = n;
%
%  INDX < 0: The user is returning the results of a comparison.
%
  elseif ( indx < 0 )

    if ( indx == -2 )

      if ( isgn < 0 )
        i_save = i_save + 1;
      end

      j_save = k1;
      k1 = i_save;
      indx = -1;
      i = i_save;
      j = j_save;
      return;
    end

    if ( 0 < isgn )
      indx = 2;
      i = i_save;
      j = j_save;
      return;
    end

    if ( k <= 1 )

      if ( n1 == 1 )
        i_save = 0;
        j_save = 0;
        indx = 0;
      else
        i_save = n1;
        n1 = n1 - 1;
        j_save = 1;
        indx = 1;
      end

      i = i_save;
      j = j_save;
      return;

    end

    k = k - 1;
    k1 = k;
%
%  0 < INDX, the user was asked to make an interchange.
%
  elseif ( indx == 1 )

    k1 = k;

  end

  while ( 1 )

    i_save = 2 * k1;

    if ( i_save == n1 )
      j_save = k1;
      k1 = i_save;
      indx = -1;
      i = i_save;
      j = j_save;
      return;
    elseif ( i_save <= n1 )
      j_save = i_save + 1;
      indx = -2;
      i = i_save;
      j = j_save;
      return;
    end

    if ( k <= 1 )
      break;
    end

    k = k - 1;
    k1 = k;

  end

  if ( n1 == 1 )
    i_save = 0;
    j_save = 0;
    indx = 0;
    i = i_save;
    j = j_save;
  else
    i_save = n1;
    n1 = n1 - 1;
    j_save = 1;
    indx = 1;
    i = i_save;
    j = j_save;
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function triangle_neighbor = triangulation_order3_neighbor_triangles ( ...
  triangle_num, triangle_node )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER3_NEIGHBOR_TRIANGLES determines triangle neighbors.
%
%  Discussion:
%
%    A triangulation of a set of nodes can be completely described by
%    the coordinates of the nodes, and the list of nodes that make up
%    each triangle.  However, in some cases, it is necessary to know
%    triangle adjacency information, that is, which triangle, if any,
%    is adjacent to a given triangle on a particular side.
%
%    This routine creates a data structure recording this information.
%
%    The primary amount of work occurs in sorting a list of 3 * TRIANGLE_NUM
%    data items.
%
%    This routine was modified to use columns instead of rows.
%
%  Example:
%
%    The input information from TRIANGLE_NODE:
%
%    Triangle   Nodes
%    --------   ---------------
%     1         3      4      1
%     2         3      1      2
%     3         3      2      8
%     4         2      1      5
%     5         8      2     13
%     6         8     13      9
%     7         3      8      9
%     8        13      2      5
%     9         9     13      7
%    10         7     13      5
%    11         6      7      5
%    12         9      7      6
%    13        10      9      6
%    14         6      5     12
%    15        11      6     12
%    16        10      6     11
%
%    The output information in TRIANGLE_NEIGHBOR:
%
%    Triangle  Neighboring Triangles
%    --------  ---------------------
%
%     1        -1     -1      2
%     2         1      4      3
%     3         2      5      7
%     4         2     -1      8
%     5         3      8      6
%     6         5      9      7
%     7         3      6     -1
%     8         5      4     10
%     9         6     10     12
%    10         9      8     11
%    11        12     10     14
%    12         9     11     13
%    13        -1     12     16
%    14        11     -1     15
%    15        16     14     -1
%    16        13     15     -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes that make up each triangle.
%
%    Output, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the three triangles that are direct
%    neighbors of a given triangle.  TRIANGLE_NEIGHBOR(1,I) is the index of the triangle
%    which touches side 1, defined by nodes 2 and 3, and so on.  TRIANGLE_NEIGHBOR(1,I)
%    is negative if there is no neighbor on that side.  In this case, that
%    side of the triangle lies on the boundary of the triangulation.
%

%
%  Step 1.
%  From the list of nodes for triangle T, of the form: (I,J,K)
%  construct the three neighbor relations:
%
%    (I,J,1,T) or (J,I,1,T),
%    (J,K,2,T) or (K,J,2,T),
%    (K,I,3,T) or (I,K,3,T)
%
%  where we choose (I,J,1,T) if I < J, or else (J,I,1,T)
%
  col = zeros ( 4, 3 * triangle_num );

  for tri = 1 : triangle_num

    i = triangle_node(1,tri);
    j = triangle_node(2,tri);
    k = triangle_node(3,tri);

    if ( i < j )
      col(1:4,1+3*(tri-1)) = [ i, j, 1, tri ]';
    else
      col(1:4,1+3*(tri-1)) = [ j, i, 1, tri ]';
    end

    if ( j < k )
      col(1:4,2+3*(tri-1)) = [ j, k, 2, tri ]';
    else
      col(1:4,2+3*(tri-1)) = [ k, j, 2, tri ]';
    end

    if ( k < i )
      col(1:4,3+3*(tri-1)) = [ k, i, 3, tri ]';
    else
      col(1:4,3+3*(tri-1)) = [ i, k, 3, tri ]';
    end

  end
%
%  Step 2. Perform an ascending dictionary sort on the neighbor relations.
%  We only intend to sort on rows 1 and 2; the routine we call here
%  sorts on rows 1 through 4 but that won't hurt us.
%
%  What we need is to find cases where two triangles share an edge.
%  Say they share an edge defined by the nodes I and J.  Then there are
%  two columns of COL that start out ( I, J, ?, ? ).  By sorting COL,
%  we make sure that these two columns occur consecutively.  That will
%  make it easy to notice that the triangles are neighbors.
%
  col = i4col_sort_a ( 4, 3*triangle_num, col );
%
%  Step 3. Neighboring triangles show up as consecutive columns with
%  identical first two entries.  Whenever you spot this happening,
%  make the appropriate entries in TRIANGLE_NEIGHBOR.
%
  triangle_neighbor(1:3,1:triangle_num) = -1;

  icol = 1;

  while ( 1 )

    if ( 3 * triangle_num <= icol )
      break
    end

    if ( col(1,icol) ~= col(1,icol+1) || col(2,icol) ~= col(2,icol+1) )
      icol = icol + 1;
      continue
    end

    side1 = col(3,icol);
    tri1 =  col(4,icol);
    side2 = col(3,icol+1);
    tri2 =  col(4,icol+1);

    triangle_neighbor(side1,tri1) = tri2;
    triangle_neighbor(side2,tri2) = tri1;

    icol = icol + 2;

  end

  return
end
function [ triangle_index, edge ] = ...
  triangulation_search_delaunay ( node_xy, ...
  triangle_num, triangle_node, triangle_neighbor, p )

%*****************************************************************************80
%
%% TRIANGULATION_SEARCH_DELAUNAY searches a Delaunay triangulation for a point.
%
%  Purpose:
%
%    The algorithm "walks" from one triangle to its neighboring triangle,
%    and so on, until a triangle is found containing point P, or P is found
%    to be outside the convex hull.
%
%    The algorithm computes the barycentric coordinates of the point with
%    respect to the current triangle.  If all three quantities are positive,
%    the point is contained in the triangle.  If the I-th coordinate is
%    negative, then P lies on the far side of edge I, which is opposite
%    from vertex I.  This gives a hint as to where to search next.
%
%    For a Delaunay triangulation, the search is guaranteed to terminate.
%    For other triangulations, a cycle may occur.
%
%    Note the surprising fact that, even for a Delaunay triangulation of
%    a set of points, the nearest point to P need not be one of the
%    vertices of the triangle containing P.
%
%    The code can be called for triangulations of any order, but only
%    the first three nodes in each triangle are considered.  Thus, if
%    higher order triangles are used, and the extra nodes are intended
%    to give the triangle a polygonal shape, these will have no effect,
%    and the results obtained here might be misleading.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the vertices.
%
%    Input, integer TRIANGLE_ORDER, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles in the triangulation.
%
%    Input, integer TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM),
%    the nodes that make up each triangle.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle
%    neighbor list.
%
%    Input, real P(2), the coordinates of a point.
%
%    Output, integer TRIANGLE_INDEX, the index of the triangle where the
%    search ended.  If a cycle occurred, then TRIANGLE_INDEX = -1.
%
%    Output, integer EDGE, indicates the position of the point P in
%    triangle TRIANGLE_INDEX:
%    0, the interior or boundary of the triangle;
%    -1, outside the convex hull of the triangulation, past edge 1;
%    -2, outside the convex hull of the triangulation, past edge 2;
%    -3, outside the convex hull of the triangulation, past edge 3.
%
  persistent triangle_index_save;

  if ( isempty ( triangle_index_save ) )
    triangle_index_save = -1;
  end

  count = 0;
  edge = 0;

  if ( triangle_index_save < 1 || triangle_num < triangle_index_save )
    triangle_index = floor ( ( triangle_num + 1 ) / 2 );
  else
    triangle_index = triangle_index_save;
  end

  while ( 1 )

    count = count + 1;

    if ( triangle_num < count )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TRIANGULATION_SEARCH_DELAUNAY - Fatal error!\n' );
      fprintf ( 1, '  The algorithm seems to be cycling.\n' );
      triangle_index = -1;
      edge = -1;
      return
    end
%
%  Get the vertices of triangle TRIANGLE_INDEX.
%
    a = triangle_node(1,triangle_index);
    b = triangle_node(2,triangle_index);
    c = triangle_node(3,triangle_index);
%
%  Using vertex C as a base, compute the distances to vertices A and B,
%  and the point P.
%
    dxa = node_xy(1,a) - node_xy(1,c);
    dya = node_xy(2,a) - node_xy(2,c);

    dxb = node_xy(1,b) - node_xy(1,c);
    dyb = node_xy(2,b) - node_xy(2,c);

    dxp = p(1)         - node_xy(1,c);
    dyp = p(2)         - node_xy(2,c);

    det = dxa * dyb - dya * dxb;
%
%  Compute the barycentric coordinates of the point P with respect
%  to this triangle.
%
    alpha = ( dxp * dyb - dyp * dxb ) / det;
    beta =  ( dxa * dyp - dya * dxp ) / det;
    gamma = 1.0 - alpha - beta;
%
%  If the barycentric coordinates are all positive, then the point
%  is inside the triangle and we're done.
%
    if ( 0.0 <= alpha && 0.0 <= beta && 0.0 <= gamma )
      break
    end
%
%  At least one barycentric coordinate is negative.
%
%  If there is a negative barycentric coordinate for which there exists
%  an opposing triangle neighbor closer to the point, move to that triangle.
%
%  (Two coordinates could be negative, in which case we could go for the
%  most negative one, or the most negative one normalized by the actual
%  distance it represents).
%
    if ( alpha < 0.0 && 0 < triangle_neighbor(2,triangle_index) )
      triangle_index = triangle_neighbor(2,triangle_index);
      continue;
    elseif ( beta < 0.0 && 0 < triangle_neighbor(3,triangle_index) )
      triangle_index = triangle_neighbor(3,triangle_index);
      continue;
    elseif ( gamma < 0.0 && 0 < triangle_neighbor(1,triangle_index) )
      triangle_index = triangle_neighbor(1,triangle_index);
      continue;
    end
%
%  All negative barycentric coordinates correspond to vertices opposite
%  sides on the convex hull.
%
%  Note the edge and exit.
%
    if ( alpha < 0.0 )
      edge = -2;
      break
    elseif ( beta < 0.0 )
      edge = -3;
      break
    elseif ( gamma < 0.0 )
      edge = -1;
      break
    end

  end

  triangle_index_save = triangle_index;

  return
end


text:
% test for denoising using Non Local Means vs Wavelets
%
%   Copyright (c) 2007 Gabriel Peyre

path(path, 'toolbox/');

if not(exist('name'))
    name = 'peppers';
    name = 'polygons_blurred';
    name = 'boat';
    name = 'mandrill';
    name = 'lenacoul';
    name = 'lena';
    name = 'images/corral';
    name = 'barb';
end

%% load image
n = 64*2;
n0 = [];
M = load_image(name, n0);
c = size(M)/2;
if strcmp(name, 'lena') && n==128
    c = [120 200]; % lena hat
end
if strcmp(name, 'mandrill') && n==128
    c = [350 350];
end
eta = 0.12;
M0 = rescale( crop(M,n,c), eta, 1-eta );
s = size(M0,3); % number of colors

%% add some Gaussian noise
sigma = 0.06 * max( M0(:) );
randn('state',1234);    % to have reproductible results
M = M0 + randn(n,n,s)*sigma;
% approx wavelet threshold
T = 3*sigma;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% TI Wavelets
clear options;
options.wavelet_type = 'biorthogonal';
options.wavelet_vm = 3;
Jmin = 4;
options.decomp_type = 'quad';
MW = perform_atrou_transform(M,Jmin,options);
disp('--> Computing best threshold (wavelets).');
Twav = compute_best_threshold('wavelet',M,M0,sigma, options);
MWT = perform_thresholding(MW, T);
MW1 = perform_atrou_transform(MWT,Jmin,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Portilla and Simoncelli method
options.sigma = sigma;
options.repres2 = 'daub3';
options.parent = 1;
disp('--> BLS-GSM denoising.');
ML1 = perform_blsgsm_denoising(M, options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Non-local means
% options of NL means  (15,30)
options.max_dist = 5;  % search width, the smaller the faster the algorithm will be
options.ndims = 40;     % number of dimension used for distance computation (PCA dim.reduc. to speed up)
options.do_patchwise = 0;
options.mask = 'linear';
options.mask = 'cst';
disp('--> NLMeans denoising.');
options.Tlist = linspace(0.02,0.07,12);
options.klist = [4 5];
[tmp,options] = compute_best_threshold('nlmeans',M,M0,sigma, options);
[MN1,Wx,Wy] = perform_nl_means(M, options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Non-local patchwise
do_patchwise = 1;
if do_patchwise
options.do_patchwise = 1;
disp('--> NLMeans denoising.');
[tmp,options] = compute_best_threshold('nlmeans',M,M0,sigma, options);
[MNpwise1,Wx,Wy] = perform_nl_means(M, options);
end

pnoisy = psnr(M,M0);
pwav = psnr(M0,MW1);
pbls= psnr(M0,ML1);
pnlm= psnr(M0,MN1);
if do_patchwise
    pnlpwise = psnr(M0,MNpwise1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Display
imgs = {M0 M MW1 ML1 MN1};
opts =  {'Original' sprintf('Noisy,psnr=%.2f', pnoisy) sprintf('Wav,psnr=%.2f', pwav) ...
        sprintf('BLS,psnr=%.2f(+%.2f)', pbls, pbls-pwav) ...
        sprintf('NLMeans,psnr=%.2f(+%.2f)', pnlm, pnlm-pwav) };
if do_patchwise
    imgs{end+1} = MNpwise1;
    opts{end+1} = sprintf('NLMedian,psnr=%.2f(+%.2f)', pnlpwise, pnlpwise-pwav);
end
display_image_layout(imgs, opts );

% save
repimg = ['results/denoising/'];
if not(exist(repimg))
    mkdir(repimg);
end
saveas(gcf, [repimg name '-denoising.png'], 'png');

repimg = ['results/denoising/' name '/'];
if not(exist(repimg))
    mkdir(repimg);
end

%% write results
fid = fopen([repimg name '-results.txt'], 'a');
fprintf(fid, '--> %s: n=%d, max_dist=%d, ndims=%d, sigma=%.2f, k=%d, T=%.2f, mask=%s\n', name, n, ...
        options.max_dist, options.ndims, sigma, options.k, options.T, options.mask);
fprintf(fid, 'Noisy    %.4f\n', pnoisy);
fprintf(fid, 'Wavelets %.4f\n', pwav);
fprintf(fid, 'BLS-GSM  %.4f (+%.4f)\n', pbls, pbls-pwav);
fprintf(fid, 'NL-Means %.4f (+%.4f)\n', pnlm, pnlm-pwav);
if do_patchwise
    fprintf(fid, 'NL-Median %.4f (+%.4f)\n', pnlpwise, pnlpwise-pwav);
end
fclose(fid);

% scaling for method noise
s = std(M0(:)-MW1(:))*6;

%% save images
warning off;
imwrite(clamp(M0), [repimg name '-original.png'], 'png');
imwrite(clamp(M), [repimg name '-noisy.png'], 'png');
imwrite(clamp(MW1), [repimg name '-wavelets.png'], 'png');
imwrite(clamp(ML1), [repimg name '-blsgsm.png'], 'png');
imwrite(clamp(MN1), [repimg name '-nlmeans.png'], 'png');
imwrite(clamp( 0.5+ (M-MW1)/s ) , [repimg name '-wavelets-methnoise.png'], 'png');
imwrite(clamp( 0.5+ (M-ML1)/s ) , [repimg name '-blsgsm-methnoise.png'], 'png');
imwrite(clamp( 0.5+ (M-MN1)/s ) , [repimg name '-nlmeans-methnoise.png'], 'png');
if do_patchwise
imwrite(clamp(MNpwise1), [repimg name '-nlpatchwise.png'], 'png');
imwrite(clamp( 0.5+ (M-MNpwise1)/s ) , [repimg name '-nlpatchwise-methnoise.png'], 'png');
end
warning off;

author: gpeyre

repo: matlab-toolboxes


text:
function rgb = magenta

% This returns a predefined color as [red green blue] values
%   red               = [255   0   0]/255;
%   green             = [  0 192   0]/255;
%   blue              = [  0   0 255]/255;
%   magenta           = [255 255   0]/255;
%   cyan              = [  0 255 255]/255;
%   yellow            = [255 255   0]/255;
%   white             = [255 255 255]/255;
%   black             = [  0   0   0]/255;
%
%   skull             = [140  85  85]/255
%   cortex            = [255 213 119]/255;
%   cortex_light      = [199 194 169]/255;
%   cortex_dark       = [100  97  85]/255;
%   skin              = [249 223 192]/255;
%   skin_light        = [249 223 192]/255;
%   skin_medium_light = [225 194 158]/255;
%   skin_medium       = [188 142 106]/255;
%   skin_medium_dark  = [155 102	65]/255;
%   skin_dark         = [ 91  71  61]/255;

rgb = [255 0 255]/255;


text:
classdef IpOptRecalcYEnum < matlab.mixin.SetGet
    %IpOptRecalcYEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        No('No','no')
        Yes('Yes','yes')
    end

    properties
        name(1,:) char
        optStr(1,:) char
    end

    methods
        function obj = IpOptRecalcYEnum(name, optStr)
            obj.name = name;
            obj.optStr = optStr;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('IpOptRecalcYEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('IpOptRecalcYEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('IpOptRecalcYEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function grad = B_beamforming(future_layers, input_layers, curr_layer)
weight = input_layers{1}.a;
nBin = length(curr_layer.freqBin);
input = input_layers{2}.a;

[Di,Ti,Ni] = size(input);
[Dw,Tw,Nw] = size(weight);

future_grad = GetFutureGrad(future_layers, curr_layer);
future_grad = permute(future_grad, [1 2 4 3]);
input = reshape(input, nBin,Di/nBin,Ti,Ni);

input2 = permute(input, [1 3 2 4]);
grad{1} = bsxfun(@times, conj(input2), future_grad);
grad{1} = permute(grad{1}, [1 3 2 4]);
grad{1} = reshape(grad{1}, Di,Ti,Ni);

if Tw==1
    grad{1} = sum(grad{1},2);
end

grad{2} = [];

end

author: singaxiong

repo: SignalGraph


text:
function [H,R] = myqr (A)
%MYQR QR factorization using Householder reflections
% uses function [v,beta,xnorm] = hmake1 (x)
% and function hx = happly (v, beta, x)
%
% Example
%   [H,R] = myqr (A)
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

[m n] = size (A) ;

H = zeros (m,n) ;
R = zeros (m,n) ;

for k = 1:n

    % apply prior H's
    % fprintf ('\n-----------------init %d\n', k) ;
    x = A (:,k) ;
    for i = 1:k-1
        v = H(((i+1):m),i) ;
        v = [1 ; v] ;                                                       %#ok
        beta = H (i,i) ;
        % n1 = norm (x (i:m)) ;
        x (i:m) = happly (v, beta, x (i:m)) ;
        % n2 = norm (x (i:m)) ;
        % fprintf ('=============== i %d %g %g\n', i, n1, n2) ;
        % beta
        % v'
        % X = x'
        % pause
        % i
        % x
    end
    % k
    % x

    % make Hk
    % fprintf ('x(k:m) = ') ; x (k:m)
    [v,beta,xnorm] = hmake1 (x (k:m)) ;

    H (k,k) = beta ;
    H (k+1:m, k) = v (2:end) ;

    R (1:(k-1),k) = x (1:(k-1)) ;
    R (k,k) = xnorm ;
    % full (R)
    % pause
end

% s2 = svd (full (R)) ;
% [s1 s2 s1-s2]


text:
function f=dft(f,N,dim);
%DFT   Normalized Discrete Fourier Transform
%   Usage: f=dft(f);
%          f=dft(f,N,dim);
%
%   `dft` computes a normalized or unitary discrete Fourier transform. The
%   unitary discrete Fourier transform is computed by
%
%   ..                     L-1
%     c(k+1) = 1/sqrt(L) * sum f(l+1)*exp(-2*pi*i*k*l/L)
%                          l=0
%
%   .. math:: c\left(k+1\right)=\frac{1}{\sqrt{L}}\sum_{l=0}^{L-1}f\left(l+1\right)e^{-2\pi ikl/L}
%
%   for $k=0,\ldots,L-1$.
%
%   The output of `dft` is a scaled version of the output from `fft`. The
%   function takes exactly the same arguments as `fft`. See the help on `fft`
%   for a thorough description.
%
%   See also:  idft

%   AUTHOR: Peter L. Søndergaard, Jordy van Velthoven
%   TESTING: TEST_DFT
%   REFERENCE: REF_DFT

complainif_argnonotinrange(nargin,1,3,mfilename);

if nargin<3
  dim=[];
end;

if nargin<2
  N=[];
end;

[f,N,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,N,dim,'DFT');

% Force FFT along dimension 1, since we have permuted the dimensions
% manually
f=fft(f,N,1)/sqrt(N);

f=assert_sigreshape_post(f,dim,permutedsize,order);


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reduce an image applying Gaussian Pyramid.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IResult = GPReduce(I,displayflag)

if ~exist('displayflag','var')
	displayflag = 1;
end

Wt = [0.0500    0.2500    0.4000    0.2500    0.0500];

dim = size(I);
newdim = ceil(dim*0.5);
IResult = zeros(newdim,class(I)); % Initialize the array in the beginning ..
I = single(I);
m = [-2:2];n=m;

switch length(dim)
	case 1
		%% Pad the boundaries.
		I = [ I(1) ; I(1) ;  I ; I(dim(1));  I(dim(1)) ];  % Add two rows towards the beginning and the end.
		for i = 0 : newdim(1) -1
			A = I(2*i+m+3).*Wt;;
			IResult(i + 1)= sum(A(:));
		end
	case 2
		%% Pad the boundaries.
		I = [ I(1,:) ; I(1,:) ;  I ; I(dim(1),:);  I(dim(1),:) ];  % Add two rows towards the beginning and the end.
		I = [ I(:,1)  I(:,1)     I   I(:,dim(2))  I(:,dim(2)) ];  % Add two columns towards the beginning and the end.

		Wt2 = Wt'*Wt;

		for i = 0 : newdim(1) -1
			for j = 0 : newdim(2) -1
				A = I(2*i+m+3,2*j+m+3).*Wt2;
				IResult(i + 1, j + 1) = sum(A(:));
			end
		end

	case 3
		Wt3 = ones(5,5,5);
		for i = 1:5
			Wt3(i,:,:) = Wt3(i,:,:) * Wt(i);
			Wt3(:,i,:) = Wt3(:,i,:) * Wt(i);
			Wt3(:,:,i) = Wt3(:,:,i) * Wt(i);
		end

		%% Pad the boundaries.
		I2 = zeros(dim+4,class(I));
		I2(3:2+dim(1),3:2+dim(2),3:2+dim(3)) = I;
		I2(1,:,:)=I2(3,:,:);I2(1,:,:)=I2(3,:,:);I2(end,:,:)=I2(end-2,:,:);I2(end-1,:,:)=I2(end-2,:,:);
		I2(:,1,:)=I2(:,3,:);I2(:,2,:)=I2(:,3,:);I2(:,end,:)=I2(:,end-2,:);I2(:,end-1,:)=I2(:,end-2,:);
		I2(:,:,1)=I2(:,:,3);I2(:,:,2)=I2(:,:,3);I2(:,:,end)=I2(:,:,end-2);I2(:,:,end-1)=I2(:,:,end-2);
		I=I2; clear I2;

		if( displayflag==1)
			H = waitbar(0,'Progress');
			set(H,'Name','GPReduce ...');
		end
		for k = 0 : newdim(3) - 1
			if( displayflag==1)
				waitbar(k/newdim(3),H,sprintf('(%d%%) %d out of %d',round(k/newdim(3)*100),k,newdim(3)));
			else
				fprintf('.');
			end
			for j = 0 : newdim(2) -1
				for i = 0 : newdim(1) -1
					A = I(2*i+m+3,2*j+m+3,2*k+m+3).*Wt3;
					IResult(i+1,j+1,k+1) = sum(A(:));
				end
			end
		end
		if( displayflag==1)
			close(H);
		else
			fprintf('\n');
		end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [w,b,sv,obj] = linear_primal_svm(lambda,wInit,bInit,D,noneg, maxIteration,opt)
% Solves the following SVM optimization problem in the primal (with quatratic
%   penalization of the training errors). Default solved by Newton.
%
% min_{w,e}  lambda/2 * w'w  + 1/2 sum_i out_i^2
%      s.t.  Y_i * (w .* X_i + b) >= D_i - out_i
%            w(nonneg)>=0
%
% A global variable X containing the training inputs
%   should be defined. X is an n x d matrix (n = number of points).
%   X can be either normal matrix or sparse matrix.
% A global variable Y is the target vector of size nx1. Normal SVM will
%   have +1 and -1 value, but it can be actually aribitury value
% A global variable n is the number of elements that you want to use for training.
% LAMBDA is the regularization parameter ( = 1/C)
% wInit is an optional input for the initial value of [w;b]
% dvec is an optional input, usually it is 1 for standard SVM
% maxIteration is the number of iterations allowd
%
% W is the hyperplane w (vector of length d).
% B is the bias
% The outputs on the training points are either X*W+B
% SV is the support vector index number
% OBJ is the objective function value
% OPT is a structure containing the options (in brackets default values):
%   cg: Do not use Newton, but nonlinear conjugate gradients [0]
%   lin_cg: Compute the Newton step with linear CG
%           [0 unless solving sparse linear SVM]
%   iter_max_Newton: Maximum number of Newton steps [20]
%   prec: Stopping criterion
%   cg_prec and cg_it: stopping criteria for the linear CG.

% Original written by Olivier Chapelle @ http://olivier.chapelle.cc/primal/
% Modified by Jianxiong Xiao to have several advance features @ http://mit.edu/jxiao/

if ~exist('maxIteration','var') || maxIteration==Inf     % Assign the options to their default values
    maxIteration = 10000000;
end

if ~exist('opt','var')       % Assign the options to their default values
    opt = [];
end
if ~isfield(opt,'cg'),                opt.cg = 0;                        end;
if ~isfield(opt,'lin_cg'),            opt.lin_cg = 0;                    end;
if ~isfield(opt,'iter_max_Newton'),   opt.iter_max_Newton = 20;          end; % used to be 20
if ~isfield(opt,'prec'),              opt.prec = 1e-6;                   end;
if ~isfield(opt,'cg_prec'),           opt.cg_prec = 1e-4;                end;
if ~isfield(opt,'cg_it'),             opt.cg_it = 20;                    end;

global X;
global Y;

if ~exist('noneg','var')
    noneg = [];
end

if ~exist('dvec','var') || isempty(D)
    D = ones(numel(Y),1);
end

if isempty(X), error('Global variable X undefined'); end;

if ~exist('bInit','var')
    bInit=0;
end
if ~exist('wInit','var')
    d = size(X,2);
    wInit = zeros(d,1);
end
if issparse(X)
    opt.lin_cg = 1;
end;
if ~opt.cg
    [sol,obj, sv] = primal_svm_linear   (lambda,maxIteration,wInit,bInit,D,noneg,opt);
else
    [sol,obj, sv] = primal_svm_linear_cg(lambda,maxIteration,wInit,bInit,D,noneg,opt);
end;

% The last component of the solution is the bias b.
b = sol(end);
w = sol(1:end-1);
fprintf('\n');

% -------------------------------
% Train a linear SVM using Newton
% -------------------------------
function  [w,obj,sv] = primal_svm_linear(lambda,maxIteration,wInit,bInit,D,noneg,opt)

global X;
global Y;
global n;
d = size(X,2);

w = [wInit; bInit]; % The last component of w is b.
w(noneg) = max(w(noneg),0);
%out = ones(n,1); % Vector containing 1-Y.*(X*w)
out = D(1:n) - Y(1:n).*(X(1:n,:)*w(1:end-1)+w(end));

for iter=1:maxIteration
    if iter > opt.iter_max_Newton;
        warning('PrimalSVM:MaxNumNewton','Maximum number of Newton steps reached. Try larger lambda');
        break;
    end;

    [obj, grad, sv] = obj_fun_linear(w,lambda,out);

    % Compute the Newton direction either exactly or by linear CG
    if opt.lin_cg
        % Advantage of linear CG when using sparse input: the Hessian is never computed explicitly.
        [step, foo, relres] = minres(@hess_vect_mult, -grad, opt.cg_prec,opt.cg_it,[],[],[],sv,lambda);
    else
        Xsv = X(sv,:);
        hess = lambda*diag([ones(d,1); 0]) + [[Xsv'*Xsv sum(Xsv,1)']; [sum(Xsv) length(sv)]];   % Hessian
        step  = - hess \ grad;   % Newton direction
    end;

    % Do an exact line search
    [t,out, sv] = line_search_linear(w,step,out, lambda);

    w = w + t*step;
    w(noneg) = max(w(noneg),0);
    fprintf('Iter = %d, Obj = %f, Nb of sv = %d, Newton decr = %.3f,  Line search = %.3f',iter,obj,length(sv),-step'*grad/2,t);
    if opt.lin_cg
        fprintf(', Lin CG acc = %.4f     \n',relres);
    else
        fprintf('      \n');
    end;

    if -step'*grad < opt.prec * obj
        % Stop when the Newton decrement is small enough
        break;
    end;
end;

% -----------------------------------------------------
% Train a linear SVM using nonlinear conjugate gradient
% -----------------------------------------------------
function  [w, obj, sv] = primal_svm_linear_cg(lambda,maxIteration,wInit,bInit, D,noneg,opt)
global X;
global Y;
global n;
d = size(X,2);

w = [wInit; bInit]; % The last component of w is b.
w(noneg) = max(w(noneg),0);
%out = ones(n,1); % Vector containing 1-Y.*(X*w)
out = D(1:n) - Y(1:n).*(X(1:n,:)*w(1:end-1)+w(end));

%go = [X(1:n,:)'*Y(1:n); sum(Y(1:n))];  % -gradient at w=0, need to be change for w!=0 initialization
[~, grad] = obj_fun_linear(w,lambda,out); go = -grad; % -gradient

s = go; % The first search direction is given by the gradient
for iter=1:maxIteration
    if iter > opt.cg_it * min(n,d)
        warning('PrimalSVM:MaxNumCG','Maximum number of CG iterations reached. Try larger lambda');
        break;
    end;

    % Do an exact line search
    [t,out,sv] = line_search_linear(w,s,out,lambda);
    w = w + t*s;
    w(noneg) = max(w(noneg),0);

    % Compute the new gradient
    [obj, gn, sv] = obj_fun_linear(w,lambda,out); gn=-gn;
    fprintf('Iter = %d, Obj = %f, Norm of grad = %.3f     \n',iter,obj,norm(gn));

    % Stop when the relative decrease in the objective function is small
    if t*s'*go < opt.prec*obj, break; end;

    % Flecher-Reeves update. Change 0 in 1 for Polack-Ribiere
    be = (gn'*gn - 0*gn'*go) / (go'*go);
    s = be*s+gn;
    go = gn;
end;

function [obj, grad, sv] = obj_fun_linear(w,lambda,out)
% Compute the objective function, its gradient and the set of support vectors
% Out is supposed to contain 1-Y.*(X*w)
global X;
global Y;
global n;
out = max(0,out);
wb0 = w; wb0(end) = 0;  % Do not penalize b <= Very important for object detection
obj = sum(out.^2)/2 + lambda*(wb0')*wb0/2; % L2 penalization of the errors
grad = lambda*wb0 - [((out.*Y(1:n))'*X(1:n,:))'; sum(out.*Y(1:n))]; % Gradient
sv = find(out>0);

function [t,out,sv] = line_search_linear(w,d,out,lambda)
% From the current solution w, do a line search in the direction d by
% 1D Newton minimization
global X;
global Y;
global n;
t = 0;
% Precompute some dots products
Xd = X(1:n,:)*d(1:end-1)+d(end);
wd = lambda * w(1:end-1)'*d(1:end-1);
dd = lambda * d(1:end-1)'*d(1:end-1);
while 1
    out2 = out - t*(Y(1:n).*Xd); % The new outputs after a step of length t
    sv = find(out2>0);
    g = wd + t*dd - (out2(sv).*Y(sv))'*Xd(sv); % The gradient (along the line)
    h = dd + Xd(sv)'*Xd(sv); % The second derivative (along the line)
    t = t - g/h; % Take the 1D Newton step. Note that if d was an exact Newton
    % direction, t is 1 after the first iteration.
    if g^2/h < 1e-10, break; end;
    %    fprintf('%f %f\n',t,g^2/h)
end;
out = out2;

function y = hess_vect_mult(w,sv,lambda)
% Compute the Hessian times a given vector x.
% hess = lambda*diag([ones(d-1,1); 0]) + (X(sv,:)'*X(sv,:));
global X;
global n;
y = lambda*w;
y(end) = 0;
z = (X(1:n,:)*w(1:end-1)+w(end));  % Computing X(sv,:)*x takes more time in Matlab :-(
zz = zeros(length(z),1);
zz(sv)=z(sv);
y = y + [(zz'*X(1:n,:))'; sum(zz)];

author: jianxiongxiao

repo: ProfXkit


text:
classdef Results
% A class that produces various measures describing system accuracy.

properties (SetAccess = private)
  tar
  non
  numtar
  numnon
  Pmiss
  Pfa
  eer
  pb
end

methods
  % res = Results(scr,key);
  % res = Results(tar,non);
  function res = Results(param1,param2)
    if isnumeric(param1)
        assert(isnumeric(param2))
        res.tar = param1;
        res.non = param2;
    else
        scr = param1;
        key = param2;
        assert(isa(scr,'Scores'))
        assert(isa(key,'Key'))
        [res.tar, res.non] = scr.get_tar_non(key);
    end
    [res.Pmiss,res.Pfa] = rocch(res.tar,res.non);
    res.eer = rocch2eer(res.Pmiss,res.Pfa);
    res.numtar = length(res.tar);
    res.numnon = length(res.non);
    Nmiss = res.Pmiss * res.numtar;
    Nfa = res.Pfa * res.numnon;
    res.pb = rocch2eer(Nmiss,Nfa);
  end
  function actdcf = get_act_dcf(res,prior)
    actdcf = fast_actDCF(res.tar,res.non,logit(prior),false);
  end
  function mindcf = get_min_dcf(res,prior)
    Ptar = prior;
    Pnon = 1-prior;
    cdet = [Ptar,Pnon]*[res.Pmiss(:)';res.Pfa(:)'];
    mindcf = min(cdet,[],2);
  end
  function actdcf = get_norm_act_dcf(res,prior)
    actdcf = fast_actDCF(res.tar,res.non,logit(prior),true);
  end
  function mindcf = get_norm_min_dcf(res,prior)
    Ptar = prior;
    Pnon = 1-prior;
    cdet = [Ptar,Pnon]*[res.Pmiss(:)';res.Pfa(:)'];
    mindcf = min(cdet,[],2) ./ min(Ptar,Pnon);
  end
  function pb = get_prbep(res)
    pb = res.pb;
  end
  function eer_val = get_eer(res)
    eer_val = res.eer;
  end
  function ntar = num_tar(res)
    ntar = res.numtar;
  end
  function nnon = num_non(res)
    nnon = res.numnon;
  end
end

end

author: nesl


text:
function updateTaxonomyInfoAGORA
% This function retrieves the newest taxonomy information for each AGORA
% strain from NCBI Taxonomy to facilitate mapping taxonomic assignments
% from metagenomic sequencing data to AGORA. An updated version of the
% AGORA information table is saved in spreadsheet format.
%
% .. Author: Almut Heinken, 04/2022

% AGORA information table
infoFile = readInputTableForPipeline('AGORA_infoFile.xlsx');

ncbiCol = find(strcmp(infoFile(1,:),'NCBI Taxonomy ID'));
strainCol = find(strcmp(infoFile(1,:),'Strain'));
speciesCol = find(strcmp(infoFile(1,:),'Species'));
genusCol = find(strcmp(infoFile(1,:),'Genus'));
familyCol = find(strcmp(infoFile(1,:),'Family'));
orderCol = find(strcmp(infoFile(1,:),'Order'));
classCol = find(strcmp(infoFile(1,:),'Class'));
phylumCol = find(strcmp(infoFile(1,:),'Phylum'));
kingdomCol = find(strcmp(infoFile(1,:),'Kingdom'));

taxa = {
    strainCol 'strain'
    speciesCol 'species'
    genusCol 'genus'
    familyCol 'family'
    orderCol 'order'
    classCol 'class'
    phylumCol 'phylum'
    kingdomCol 'superkingdom'
    };

for i=2:size(infoFile,1)
    i
    taxonomy = parseNCBItaxonomy(infoFile{i,ncbiCol});
    for j=1:size(taxa,1)
        if isfield(taxonomy,taxa{j,2})
            if j==1
                if ~strcmp(taxonomy.('strain'),'root')
                    if isfield(taxonomy,'species')
                        if ~strcmp(taxonomy.('strain'),taxonomy.('species'))
                            if ~isempty(taxonomy.(taxa{j,2}))
                                infoFile{i,taxa{j,1}} = taxonomy.(taxa{j,2});
                            end
                        end
                    end
                end
            else
                if ~isempty(taxonomy.(taxa{j,2}))
                    infoFile{i,taxa{j,1}} = taxonomy.(taxa{j,2});
                end
            end
        end
    end
end
writetable(cell2table(infoFile),'AGORA_infoFile.xlsx','FileType','spreadsheet','writeVariableNames',false);

% AGORA2 information table
infoFile = readInputTableForPipeline('AGORA2_infoFile.xlsx');

ncbiCol = find(strcmp(infoFile(1,:),'NCBI Taxonomy ID'));
strainCol = find(strcmp(infoFile(1,:),'Strain'));
speciesCol = find(strcmp(infoFile(1,:),'Species'));
genusCol = find(strcmp(infoFile(1,:),'Genus'));
familyCol = find(strcmp(infoFile(1,:),'Family'));
orderCol = find(strcmp(infoFile(1,:),'Order'));
classCol = find(strcmp(infoFile(1,:),'Class'));
phylumCol = find(strcmp(infoFile(1,:),'Phylum'));
kingdomCol = find(strcmp(infoFile(1,:),'Kingdom'));

taxa = {
    strainCol 'strain'
    speciesCol 'species'
    genusCol 'genus'
    familyCol 'family'
    orderCol 'order'
    classCol 'class'
    phylumCol 'phylum'
    kingdomCol 'superkingdom'
    };

for i=2:size(infoFile,1)
    i
    taxonomy = parseNCBItaxonomy(infoFile{i,ncbiCol});
    for j=1:size(taxa,1)
        if isfield(taxonomy,taxa{j,2})
            if j==1
                if ~strcmp(taxonomy.('strain'),'root')
                    if isfield(taxonomy,'species')
                        if ~strcmp(taxonomy.('strain'),taxonomy.('species'))
                            if ~isempty(taxonomy.(taxa{j,2}))
                                infoFile{i,taxa{j,1}} = taxonomy.(taxa{j,2});
                            end
                        end
                    end
                end
            else
                if ~isempty(taxonomy.(taxa{j,2}))
                    infoFile{i,taxa{j,1}} = taxonomy.(taxa{j,2});
                end
            end
        end
    end
end
writetable(cell2table(infoFile),'AGORA2_infoFile.xlsx','FileType','spreadsheet','writeVariableNames',false );

end

author: opencobra

repo: cobratoolbox