text:
%% FUNCTION Least_Dirty
%   Dirty Multi-Task Learning with Least Squares Loss.
%
%% OBJECTIVE
%   argmin_W ||X(P+Q) - Y||_F^2 + lambda1*||P||_{1,inf} + lambda2*||Q||_{1,1}
%    s.t. W = P + Q
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   rho1: group sparsity regularization parameter
%   rho2: elementwise sparsity regularization parameter
%
%% OUTPUT
%   W: model: d * t
%   funcVal: function value vector.
%   P: group sparsity structure (joint feature slection)
%   Q: elementwise sparsity component
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Pinghua Gong and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Jalali, A. and Ravikumar, P. and Sanghavi, S. and Ruan, C. A dirty
%       model for multi-task learning, NIPS 2010.
%
%% RELATED FUNCTIONS
%   init_opts, combine_input (utils), prf_lbm (c_files)

function [W, funcVal, P, Q] = Least_Dirty(X, Y, lambda1, lambda2, opts)

if nargin <4
    error('\n Inputs: X, Y, and lambda1, and lambda2 should be specified!\n');
end
if nargin <5
    opts = [];
end

% initialize options.
opts=init_opts(opts);

% initial Lipschiz constant.
if isfield(opts, 'lFlag')
    lFlag = opts.lFlag;
else
    lFlag = false;
end

task_num = length(X);
[X, y, ~, samplesize] = combine_input(X, Y);
dimension = size(X, 2);

% initialize a starting point
if opts.init==2
    P0 = zeros(dimension, task_num);
    Q0 = zeros(dimension, task_num);
elseif opts.init == 0
    P0 = randn(dimension, task_num);
    Q0 = randn(dimension, task_num);
else
    if isfield(opts,'P0')
        P0=opts.P0;
        if (nnz(size(P0)-[dimension, task_num]))
            error('\n Check the input .P0');
        end
    else
        P0=zeros(dimension, task_num);
    end

    if isfield(opts,'Q0')
        Q0=opts.Q0;
        if (nnz(size(Q0)-[dimension, task_num]))
            error('\n Check the input .Q0');
        end
    else
        Q0=zeros(dimension, task_num);
    end
end

% Set an array to save the objective value
funcVal = [];

P = P0; Q = Q0;

%[d,m] = size(P); % d: dimension, m: the number of tasks
X = diagonalize(X,samplesize);
XtX = X'*X; Xty = X'*y;

Pn = P;  Qn = Q;
t_new = 1;
L1norm = max(sum(abs(X),1)); Linfnorm = max(sum(abs(X),2));

if lFlag
    % Upper bound for largest eigenvalue of Hessian matrix
    L = 2*min([L1norm*Linfnorm; size(X,1)*Linfnorm*Linfnorm; size(X,2)*L1norm*L1norm; size(X,1)*size(X,2)*max(abs(X(:)))]);
else
    % Lower bound for largest eigenvalue of Hessian matrix
    L = 2*max(L1norm*L1norm/size(X,1),Linfnorm*Linfnorm/size(X,2));
end
% Initial function value
funcVal = cat(1, funcVal, norm(X*(P(:) + Q(:)) - y)^2 + lambda1*L1infnorm(P) + lambda2*L11norm(Q));
%fun(1) = norm(X*(P(:) + Q(:)) - y)^2 + lambda1*L1infnorm(P) + lambda2*L11norm(Q);

%count = 0;
for iter = 1 : opts.maxIter
    P_old = P; Q_old = Q;
    t_old = t_new;
    gradvec = 2*(XtX*(Pn(:)+Qn(:)) - Xty);
    gradmat = reshape(gradvec,dimension,task_num);
    % If we estimate the upper bound of Lipschitz constant, no line search
    % is needed.
    if lFlag
        P = proximalL1infnorm(Pn - gradmat/L, lambda1/L);
        Q = proximalL11norm(Qn - gradmat/L, lambda2/L);
    else
        % line search
        for inneriter = 1:20
            P = proximalL1infnorm(Pn - gradmat/L, lambda1/L);
            Q = proximalL11norm(Qn - gradmat/L, lambda2/L);
            dP = P - Pn;  dQ = Q - Qn;
            if 2*((dP(:) + dQ(:))'*XtX*(dP(:) + dQ(:))) <= L*sum(sum((dP.*dP + dQ.*dQ)))
                break;
            else
                L = L*2;
            end
        end
    end
    funcVal = cat(1, funcVal, norm(X*(P(:) + Q(:)) - y)^2 + lambda1*L1infnorm(P) + lambda2*L11norm(Q));
    %fun(iter+1) = norm(X*(P(:) + Q(:)) - y)^2 + lambda1*L1infnorm(P) + lambda2*L11norm(Q);

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    %     % stopping condition
    %     if abs(fun(iter) - fun(iter+1))/fun(iter+1) < tol
    %         count = count + 1;
    %     else
    %         count = 0;
    %     end
    %     if count >= 3
    %         break;
    %     end

    % Update the coefficient
    t_new = (1+sqrt(1+4*t_old^2))/2;
    Pn = P + (t_old-1)/t_new*(P - P_old);
    Qn = Q + (t_old-1)/t_new*(Q - Q_old);

end

W = P + Q;

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
load('../../word2vector_matlab/GoogleNews_words.mat');
load('../../word2vector_matlab/GoogleNews_vectors.mat');
load('caption_train.mat');
w_sum = cellfun(@(x) sum(x),w_names);
subset.is_appear = zeros(1,3000000,'single');

for i = 1:numel(caption_dic)
    disp(i);
    split_tline = strsplit(caption_dic{i},{'-',' ','.',',','(',')'},'CollapseDelimiters',true);
    for j=1:numel(split_tline)  % use sum to do hash
        word = split_tline{j};
        sub_index = find(w_sum == sum(word));
        if(isempty(sub_index))
            continue;
        end
        ind = cellfun(@(x) strcmp(x,word),w_names(sub_index));
        if(sum(ind(:))==0)
            continue;
        end
        index = sub_index(ind==1);
        %tmp = w_features(:,index);
        %sentence = sentence+tmp;
        subset.is_appear(index) = subset.is_appear(index) +1;
    end
end

sub = find(subset.is_appear>0);
subset.names = {w_names{sub}};
subset.features = w_features(:,sub);
save('CUHK-PEDES_dictionary.mat','subset');

author: layumi

repo: Image-Text-Embedding


text:
%%%  Replication files for:
%%%  ""Nowcasting", 2010, (by Marta Banbura, Domenico Giannone and Lucrezia Reichlin),
%%% in Michael P. Clements and David F. Hendry, editors, Oxford Handbook on Economic Forecasting.
%%%
%%% The software can be freely used in applications.
%%% Users are kindly requested to add acknowledgements to published work and
%%% to cite the above reference in any resulting publications
%
%Description:
%
%remNaNs    Treats NaNs in dataset for use in DFM.
%
%  Syntax:
%    [X,indNaN] = remNaNs(X,options)
%
%  Description:
%    remNaNs() processes NaNs in a data matrix X according to 5 cases (see
%    below for details). These are useful for running functions in the
%    'DFM.m' file that do not take missing value inputs.
%
%  Input parameters:
%    X (T x n): Input data where T gives time and n gives the series.
%    options: A structure with two elements:
%      options.method (numeric):
%      - 1: Replaces all missing values using filter().
%      - 2: Replaces missing values after removing trailing and leading
%           zeros (a row is 'missing' if >80% is NaN)
%      - 3: Only removes rows with leading and closing zeros
%      - 4: Replaces missing values after removing trailing and leading
%           zeros (a row is 'missing' if all are NaN)
%      - 5: Replaces missing values with spline() then runs filter().
%
%      options.k (numeric): remNaNs() relies on MATLAB's filter function
%      for the 1-D filter. k controls the rational transfer function
%      argument's numerator (the denominator is set to 1). More
%      specifically, the numerator takes the form 'ones(2*k+1,1)/(2*k+1)'
%      For additional help, see MATLAB's documentation for filter().
%
%  Output parameters:
%    X: Outputted data.
%    indNaN: A matrix indicating the location for missing values (1 for NaN).

function [X,indNaN] = remNaNs_spline(X,options)
[T,N]=size(X);    % Gives dimensions for data input
k=options.k;      % Inputted options
indNaN=isnan(X);  % Returns location of NaNs

switch options.method
    case 1  % replace all the missing values
        for i = 1:N  % loop through columns
            x = X(:,i);
            x(indNaN(:,i))= nan_median(x);  % Replace missing values series median
            x_MA =filter (ones(2*k+1,1)/(2*k+1), 1, [x(1)*ones(k,1);x;x(end)*ones(k,1)]);
            x_MA=x_MA(2*k+1:end);  % Match dimensions
            % Replace missing observations with filtered values
            x(indNaN(:,i))=x_MA(indNaN(:,i));
            X(:,i)=x;  % Replace vector
        end
    case 2 %replace missing values after removing leading and closing zeros

        % Returns row sum for NaN values. Marks true for rows with more
        % than 80% NaN
        rem1=(sum(indNaN,2)>N*0.8);
        nanLead =(cumsum(rem1)==(1:T)');
        nanEnd =(cumsum(rem1(end:-1:1))==(1:T)');
        nanEnd = nanEnd(end:-1:1);  % Reverses nanEnd
        nanLE = (nanLead | nanEnd);

        % Subsets X for for
        X(nanLE,:) = [];
        indNaN = isnan(X);  % Index for missing values
        % Loop for each series
        for i = 1:N
            x = X(:,i);
            isnanx = isnan(x);
            t1 = min(find(~isnanx));  % First non-NaN entry
            t2 = max(find(~isnanx));  % Last non-NaN entry
            % Interpolates without NaN entries in beginning and end
            x(t1:t2) = spline(find(~isnanx),x(~isnanx),(t1:t2)');
            isnanx = isnan(x);
            % replace NaN observations with median
            x(isnanx) = median(x,'omitnan');
            % Apply filter
            x_MA = filter (ones(2*k+1,1)/(2*k+1),1,[x(1)*ones(k,1);x;x(end)*ones(k,1)]);
            x_MA = x_MA(2*k+1:end);
            % Replace nanx with filtered observations
            x(isnanx) = x_MA(isnanx);
            X(:,i) = x;
        end
    case 3 %only remove rows with leading and closing zeros
        rem1=(sum(indNaN,2)==N);
        nanLead=(cumsum(rem1)==(1:T)');
        nanEnd=(cumsum(rem1(end:-1:1))==(1:T)');
        nanEnd=nanEnd(end:-1:1);
        nanLE=(nanLead | nanEnd);
        X(nanLE,:)=[];
        indNaN = isnan(X);
    case 4  %remove rows with leading and closing zeros & replace missing values
        rem1=(sum(indNaN,2)==N);
        nanLead=(cumsum(rem1)==(1:T)');
        nanEnd=(cumsum(rem1(end:-1:1))==(1:T)');
        nanEnd=nanEnd(end:-1:1);
        nanLE=(nanLead | nanEnd);
        X(nanLE,:)=[];
        indNaN=isnan(X);
        for i = 1:N
            x = X(:,i);
            isnanx = isnan(x);
            t1 = min(find(~isnanx));
            t2 = max(find(~isnanx));
            x(t1:t2) = spline(find(~isnanx),x(~isnanx),(t1:t2)');
            isnanx = isnan(x);
            x(isnanx)=nan_median(x);
            x_MA =filter (ones(2*k+1,1)/(2*k+1),1,[x(1)*ones(k,1);x;x(end)*ones(k,1)]);
            x_MA=x_MA(2*k+1:end);
            x(isnanx)=x_MA(isnanx);
            X(:,i)=x;
        end
    case 5 %replace missing values
        indNaN=isnan(X);
        for i = 1:N
            x = X(:,i);
            isnanx = isnan(x);
            t1 = min(find(~isnanx));
            t2 = max(find(~isnanx));
            x(t1:t2) = spline(find(~isnanx),x(~isnanx),(t1:t2)');
            isnanx = isnan(x);
            x(isnanx)=nan_median(x);
            x_MA =filter (ones(2*k+1,1)/(2*k+1),1,[x(1)*ones(k,1);x;x(end)*ones(k,1)]);
            x_MA=x_MA(2*k+1:end);
            x(isnanx)=x_MA(isnanx);
            X(:,i)=x;
        end
end

author: FRBNY-TimeSeriesAnalysis

repo: Nowcasting


text:
function [positions, time] = tracker_ransac(video_path, img_files, pos, target_sz, ...
	padding, kernel, lambda, output_sigma_factor, interp_factor, cell_size, ...
	features, show_visualization)
%TRACKER Kernelized/Dual Correlation Filter (KCF/DCF) tracking.
%   This function implements the pipeline for tracking with the KCF (by
%   choosing a non-linear kernel) and DCF (by choosing a linear kernel).
%
%   It is meant to be called by the interface function RUN_TRACKER, which
%   sets up the parameters and loads the video information.
%
%   Parameters:
%     VIDEO_PATH is the location of the image files (must end with a slash
%      '/' or '\').
%     IMG_FILES is a cell array of image file names.
%     POS and TARGET_SZ are the initial position and size of the target
%      (both in format [rows, columns]).
%     PADDING is the additional tracked region, for context, relative to
%      the target size.
%     KERNEL is a struct describing the kernel. The field TYPE must be one
%      of 'gaussian', 'polynomial' or 'linear'. The optional fields SIGMA,
%      POLY_A and POLY_B are the parameters for the Gaussian and Polynomial
%      kernels.
%     OUTPUT_SIGMA_FACTOR is the spatial bandwidth of the regression
%      target, relative to the target size.
%     INTERP_FACTOR is the adaptation rate of the tracker.
%     CELL_SIZE is the number of pixels per cell (must be 1 if using raw
%      pixels).
%     FEATURES is a struct describing the used features (see GET_FEATURES).
%     SHOW_VISUALIZATION will show an interactive video if set to true.
%
%   Outputs:
%    POSITIONS is an Nx2 matrix of target positions over time (in the
%     format [rows, columns]).
%    TIME is the tracker execution time, without video loading/rendering.
%
%   Joao F. Henriques, 2014

	%if the target is large, lower the resolution, we don't need that much
	%detail
	resize_image = (sqrt(prod(target_sz)) >= 100);  %diagonal size >= threshold
	if resize_image,
		pos = floor(pos / 2);
		target_sz = floor(target_sz / 2);
	end

	%window size, taking padding into account
	window_sz = floor(target_sz * (1 + padding));

% 	%we could choose a size that is a power of two, for better FFT
% 	%performance. in practice it is slower, due to the larger window size.
% 	window_sz = 2 .^ nextpow2(window_sz);

	%create regression labels, gaussian shaped, with a bandwidth
	%proportional to target size
	output_sigma = sqrt(prod(target_sz)) * output_sigma_factor / cell_size;
	yf = fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

	%store pre-computed cosine window
	cos_window = hann(size(yf,1)) * hann(size(yf,2))';

	if show_visualization,  %create video interface
		update_visualization = show_video(img_files, video_path, resize_image);
	end

	%note: variables ending with 'f' are in the Fourier domain.

	time = 0;  %to calculate FPS
	positions = zeros(numel(img_files), 2);  %to calculate precision

	for frame = 1:numel(img_files),
		%load image
		im = imread([video_path img_files{frame}]);
		if size(im,3) > 1,
			im = rgb2gray(im);
		end
		if resize_image,
			im = imresize(im, 0.5);
		end

		tic()

		if frame == 1
			camera_feat.detector = cv.FeatureDetector('SURF');
			camera_feat.extractor = cv.DescriptorExtractor('SURF');
			camera_feat.matcher = cv.DescriptorMatcher('FlannBased');
			camera_feat.last_keypoints = camera_feat.detector.detect(im);
			camera_feat.last_descriptors = camera_feat.extractor.compute(im, camera_feat.last_keypoints);
		end

		if frame > 1,
			[camera_feat, camera_H] = find_homography(camera_feat, im, pos, target_sz);
			pos = project_t(pos([2, 1]), camera_H);
            pos = pos([2, 1]);
            pos = min([size(im,1),size(im,2)], pos);
            pos = max([0,0], pos);
			%obtain a subwindow for detection at the position from last
			%frame, and convert to Fourier domain (its size is unchanged)
			patch = get_subwindow(im, pos, window_sz);
			zf = fft2(get_features(patch, features, cell_size, cos_window));

			%calculate response of the classifier at all shifts
			switch kernel.type
			case 'gaussian',
				kzf = gaussian_correlation(zf, model_xf, kernel.sigma);
			case 'polynomial',
				kzf = polynomial_correlation(zf, model_xf, kernel.poly_a, kernel.poly_b);
			case 'linear',
				kzf = linear_correlation(zf, model_xf);
			end
			response = real(ifft2(model_alphaf .* kzf));  %equation for fast detection

			%target location is at the maximum response. we must take into
			%account the fact that, if the target doesn't move, the peak
			%will appear at the top-left corner, not at the center (this is
			%discussed in the paper). the responses wrap around cyclically.
			[vert_delta, horiz_delta] = find(response == max(response(:)), 1);
			if vert_delta > size(zf,1) / 2,  %wrap around to negative half-space of vertical axis
				vert_delta = vert_delta - size(zf,1);
			end
			if horiz_delta > size(zf,2) / 2,  %same for horizontal axis
				horiz_delta = horiz_delta - size(zf,2);
			end
			pos = pos + cell_size * [vert_delta - 1, horiz_delta - 1];
		end

		%obtain a subwindow for training at newly estimated target position
		patch = get_subwindow(im, pos, window_sz);
		xf = fft2(get_features(patch, features, cell_size, cos_window));

		%Kernel Ridge Regression, calculate alphas (in Fourier domain)
		switch kernel.type
		case 'gaussian',
			kf = gaussian_correlation(xf, xf, kernel.sigma);
		case 'polynomial',
			kf = polynomial_correlation(xf, xf, kernel.poly_a, kernel.poly_b);
		case 'linear',
			kf = linear_correlation(xf, xf);
		end
		alphaf = yf ./ (kf + lambda);   %equation for fast training

		if frame == 1,  %first frame, train with a single image
			model_alphaf = alphaf;
			model_xf = xf;
		else
			%subsequent frames, interpolate model
			model_alphaf = (1 - interp_factor) * model_alphaf + interp_factor * alphaf;
			model_xf = (1 - interp_factor) * model_xf + interp_factor * xf;
		end

		%save position and timing
		positions(frame,:) = pos;
		time = time + toc();

		%visualization
		if show_visualization,
			box = [pos([2,1]) - target_sz([2,1])/2, target_sz([2,1])];
			stop = update_visualization(frame, box);
			if stop, break, end  %user pressed Esc, stop early

			drawnow
% 			pause(0.05)  %uncomment to run slower
		end

	end

	if resize_image,
		positions = positions * 2;
	end
end

author: flyers

repo: drone-tracking


text:
function [U2] = mle_lam0(phi2,PP,PnP)
%
%
% this function calculates the maximum likelihood function U2
% (lambda) with 0 Fourier harmonics
%
%
% zero harmonic
%
Np=sum(sum(PnP));
sumPP=sum(sum(PP));
U2=Np*log(phi2(1))-sumPP*phi2(1);

author: Sable

repo: mcbench-benchmarks


text:
function jed = nyt_to_jed ( volume, issue )

%*****************************************************************************80
%
%% NYT_TO_JED converts an NYT date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer VOLUME, ISSUE, the New York Times
%    volume and issue.
%
%    Output, real JED, the Julian Ephemeris Date.
%
  jed_epoch_50000 = 2449790.5;

  if ( 149 < volume )

    jed = jed_epoch_50000 + issue - 50000 + 500;
%
%  Take care of the bizarre case of the second half of Volume 149,
%  Jan 1 2000 to Sep 17 2000, issues 51254 through ?, which were also
%  lowered by 500.
%
  elseif ( volume == 149 && issue < 51600 )
    jed = jed_epoch_50000 + issue - 50000 + 500;
  elseif ( 44028 <= issue )
    jed = jed_epoch_50000 + issue - 50000;
%
%  Factor in the strike of 1978.
%
  else
    jed = jed_epoch_50000 + issue - 50000 - 88;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Test file for @chebfun/minandmax.m.

function pass = test_minandmax(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% Check empty case.
[y, x] = minandmax(chebfun());
pass(1) = isempty(y) && isempty(x);

% Check operation without breakpoints.
f = chebfun(@(x) ((x - 0.2).^3 - (x - 0.2) + 1).*sec(x - 0.2), pref);
[y, x] = minandmax(f);
y_exact = [0.710869767377087 ; 1.884217141925336];
pass(2) = all(abs(y - y_exact) <= 10*vscale(f)*eps) && ...
    all(abs(feval(f, x) - y_exact) <= 10*vscale(f)*eps);

% Check operation with breakpoints.
f = chebfun(@(x) ((x - 0.2).^3 - (x - 0.2) + 1).*sec(x - 0.2), ...
    linspace(-1, 1, 10), pref);
[y, x] = minandmax(f);
y_exact = [0.710869767377087 ; 1.884217141925336];
pass(3) = all(abs(y - y_exact) <= 10*vscale(f)*eps) && ...
    all(abs(feval(f, x) - y_exact) <= 10*vscale(f)*eps);

% Check operation for complex-valued chebfuns.
f = chebfun({@(x) exp((1 + 1i)*x), @(x) 1-x/10}, [-1 0 1],  pref);
[y, x] = minandmax(f);
y_exact = [exp(-1 - 1i) ; 1];
err1 = abs(y - y_exact);
err2 = abs(feval(f, x) - y_exact);
pass(4) = all(err1 <= 10*vscale(f)*eps) && ...
    all(err2 <= 10*vscale(f)*eps);

% Check operation for pointValues.
f = chebfun({-1, 1, 2}, [-1, 0, 1, 2]);
[y, ignored] = minandmax(f);
pass(5) = all(y == [-1 ; 2]);

f.pointValues(1,1) = 10;
f.pointValues(3,1) = -10;
[y, x] = minandmax(f);
pass(6) = all(y == [-10 ; 10]) && all(x == [1 ; -1]);

% Check computation of local extrema.
f = chebfun(@(x) sin(x).^2 + sin(x.^2), [0, 4]);
[y, x] = minandmax(f, 'local');
y_exact = [        0
   1.923771282655145
  -0.342247088203205
   1.117294907913736
  -0.971179645473729
   1.343997479566445
   0.284846700239241];
x_exact = [        0
   1.323339426259694
   2.220599667639221
   2.781195946808315
   3.308480466603983
   3.776766383330969
   4.000000000000000];
pass(7) = numel(y == 7) && norm(y - y_exact, inf) < 10*vscale(f)*eps;

% Check operation for array-valued chebfuns.
f = chebfun(@(x) [sin(10*x) cos(10*x) exp(x)], [-1 -0.5 0.5 1]);
[y, x] = minandmax(f);
y_exact = [-1 -1 exp(-1) ; 1 1 exp(1)];
fx = feval(f, x(:));
fx = [fx(1:2,1) fx(3:4,2) fx(5:6,3)];
pass(8) = all(abs(y(:) - y_exact(:)) <= 10*vscale(f)*eps) && ...
    all(abs(fx(:) - y_exact(:)) <= 10*vscale(f)*eps);

f = chebfun(@(x) [exp((1 + 1i)*x) sec(1i*(x - 0.5))], [-1 0 1], ...
    pref);
[y, x] = minandmax(f);
y_exact = [exp(-1 - 1i) sec(-1.5i) ; exp(1 + 1i) 1];
fx = feval(f, x(:));
fx = [fx(1:2,1) fx(3:4,2)];
pass(9) = all(abs(y(:) - y_exact(:)) <= 10*vscale(f)*eps) && ...
    all(abs(fx(:) - y_exact(:)) <= 10*vscale(f)*eps);

op = @(x) sin(x).^2 + sin(x.^2);
f = chebfun(@(x) [op(x) op(x/2)], [0, 4]);
[y, x] = minandmax(f, 'local');
y_exact = [        0                  0
   1.923771282655145  1.923771282655145
  -0.342247088203205  0.070019315123878
   1.117294907913736  NaN
  -0.971179645473729  NaN
   1.343997479566445  NaN
   0.284846700239241  NaN];
x_exact = [        0                  0
   1.323339426259694  2.646678852519388
   2.220599667639221  4.000000000000000
   2.781195946808315  NaN
   3.308480466603983  NaN
   3.776766383330969  NaN
   4.000000000000000  NaN];
fx1 = feval(f, x_exact(:,1));
fx2 = feval(f, x_exact(1:3,2));
pass(10) = isequal(size(y), [7 2]) && ...
    all(isnan(y(4:end,2))) && all(isnan(x(4:end,2)));
pass(11) = norm(y(:,1) - y_exact(:,1), inf) < 10*vscale(f)*eps && ...
    norm(y(1:3,2) - y_exact(1:3,2), inf) < 10*vscale(f)*eps && ...
    norm(fx1(:,1) - y_exact(:,1), inf) < 10*vscale(f)*eps && ...
    norm(fx2(:,2) - y_exact(1:3,2), inf) < 10*vscale(f)*eps;

%% Test on singular function: piecewise smooth chebfun - splitting on.

dom = [-2 7];
pow = -0.5;
op = @(x) (x - dom(1)).^pow.*(sin(300*x).^2);
f = chebfun(op, dom, 'exps', [pow 0], 'splitting', 'on');
[y, x] = minandmax(f);
y_exact = [0 ; Inf];
fx = op(x);
pass(12) = ((max(abs(y - y_exact)) < 1e4*eps) && ...
          (max(abs(fx - y_exact)) < 1e4*eps));

%% Tests on function defined on unbounded domain:

% Doubly-infinite domain:

% Set the domain:
dom = [-Inf Inf];

op = @(x) (1-exp(-x.^2))./x;
f = chebfun(op, dom);
[vals, pos] = minandmax(f);
% These exact solutions are obtained using Mathematica:
vExact = [-0.6381726863389515 ; 0.6381726863389515];
pExact = [-1.120906422778534 ; 1.120906422778534];
errV = vals - vExact;
errP = pos - pExact;
pass(13) = ( norm(errV, inf) < 1e2*eps*vscale(f) ) && ...
    ( norm(errP, inf) < 1e2*eps*vscale(f) );

%% from #824
gam = chebfun('gamma(x)',[-4 4],'blowup','on','splitting','on');
mm = minandmax(1./gam);
pass(14) = norm(mm - [-1.125953228398760;4.079508980001102]) < 1e4*eps;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function imu=readimu(opt,fname)

global glc gls
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading imu file %s...',fname0);

ts=str2time(opt.ts); te=str2time(opt.te);
imu=gls.imu; ins=opt.ins;
NMAX=200000; imu.data=repmat(gls.imud,NMAX,1);
imudata=csvread(fname,1,0); n=size(imudata,1);

if ins.data_format==1
    dt=1/ins.sample_rate;
    for i=1:n
        time=gpst2time(imudata(i,1),imudata(i,2));
        if ts.time~=0&&timediff(time,ts)<0,continue;end
        if te.time~=0&&timediff(time,te)>0,continue;end

        if imu.n+1>size(imu.data,1)
            imu.data(imu.n+1:imu.n+NMAX)=repmat(gls.imud,NMAX,1);
        end
        imu.data(imu.n+1).time=time;
        imu.data(imu.n+1).dw(1:3)=imudata(i,3:5)*dt;
        imu.data(imu.n+1).dv(1:3)=imudata(i,6:8)*dt;
        imu.n=imu.n+1;
    end
elseif ins.data_format==2
    for i=1:n
        time=gpst2time(imudata(i,1),imudata(i,2));
        if ts.time~=0&&timediff(time,ts)<0,continue;end
        if te.time~=0&&timediff(time,te)>0,continue;end

        if imu.n+1>size(imu.data,1)
            imu.data(imu.n+1:imu.n+NMAX)=repmat(gls.imud,NMAX,1);
        end
        imu.data(imu.n+1).time=time;
        imu.data(imu.n+1).dw(1:3)=imudata(i,3:5);
        imu.data(imu.n+1).dv(1:3)=imudata(i,6:8);
        imu.n=imu.n+1;
    end
end

if imu.n<size(imu.data,1)
    imu.data(imu.n+1:end,:)=[];
end

fprintf('over\n');

return

repo: GINav


text:
% --- Executes on selection change in popupmenuColorMap.
function popupmenuColorMap_Callback(hObject, eventdata, handles)
    handles.data.cmapName = get(handles.popupmenuColorMap, 'String');
    handles.data.cmapName = handles.data.cmapName(get(handles.popupmenuColorMap, 'Value'));
    handles.data.cmap = feval(handles.data.cmapName{1, 1}, 256);
    colormap(handles.focusWindow, handles.data.cmap);
    colormap(handles.spectogramWindow, handles.data.cmap);
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function v = sum2(f)
%SUM2   Double integral of a SPHEREFUN over its domain.
%   I = SUM2(F) returns the double definite integral of a SPHEREFUN.

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

% Split f into its plus/minus terms.  The minus terms have integral zero
% on the sphere since the rows in this case are anti-periodic with period
% pi. Thus, we only need to integrate the plus terms.

[cols, d, rows] = cdr(f);

% If there are no plus terms then the integral is zero
if ( isempty(f.idxPlus) )
    v = 0;
    return
end

cols = cols(:, f.idxPlus); % Just keep the plus terms
rows = rows(:, f.idxPlus);
d = diag(d(f.idxPlus, f.idxPlus)).';

% Integrate the rows over their domain.
intRows = sum(rows);

% One could use the following code to do the integrals in latitude (theta),
% but this can be slow due to the fact that the integrand (which are
% trigfuns) may first be converted to chebfuns when sum (with an interval
% other than the period of the integrand) is called.  Mathematically, this
% is, of course, unnecessary and I tried arguing that this should not be
% the case.  However, I was not successful in convincing everyone and so we
% are stuck with a potentially slow method to do definite integrals rather
% than a fast one. See ticket numbers #1004 and #1034 for more details.
%
% We are going to do the integral the fast way.

% Slow code: Left here in case someone ever makes sum(f,[a,b]) fast for
% trigfuns.
%
% % Create a trigfun of the measure for the sphere.
% measure = chebfun(@(x) sin(x), f.domain(3:4), 'trig');
%
% % Multiply the columns by the measure
% cols = cols.*(measure*ones(1, size(cols, 2)));
%
% % Integrate each column over the non-doubled up latitude coordinate.
% intCols = sum(cols, [0 pi]);

%
% Fast code: We know the columns are even functions, which means they have
% cosine series expansions:
% col(:, j) = sum_{k=0}^{m} a_k cos(k*t)
% So, in the case of the elevation angle being measured as co-latitude, we
% want to compute the integral
% int_{0}^{pi}col(:,j).*sin(t)dt = sum_{k=0}^{m} a_k int_{0}^{pi}cos(k*t).*sin(t)dt
% This simplifies down to
% int_{0}^{pi}col(:,j).*sin(t)dt = sum_{k=0}^{m} a_k (1+(-1)^k)/(1-k^2)

[a, ignore] = trigcoeffs(cols);

k = (0:size(a, 1)-1).';
intFactor = 2./(1-k(1:2:end).^2);
intCols = sum(bsxfun(@times, a(1:2:end, :), intFactor));

% Put the integrals together to get the final result.
v = sum(d.*intRows.*intCols);

% TODO: Add support for different domains.

end


text:
%% Copyright 2012 The MathWorks, Inc.
%% Face Tracking
% Open video file
videoFileReader = vision.VideoFileReader('MWface.avi');

% Read video frame
videoFrame = step(videoFileReader);

% Detect the face
faceDetector = vision.CascadeObjectDetector('MaxSize',[150 150]);
faceBBox     = step(faceDetector,videoFrame);

% Display results
figure; subplot(1,2,1);imshow(videoFrame); hold on;
rectangle('Position',faceBBox(1,:),'LineWidth',2,'EdgeColor',[1 1 0])
[hueChannel,~,~] = rgb2hsv(videoFrame);
subplot(1,2,2);imshow(hueChannel);hold on;
rectangle('Position',faceBBox(1,:),'LineWidth',2,'EdgeColor',[0 1 0])

%% Find a better object for tracking
% Detect the nose
noseDetector = vision.CascadeObjectDetector('Nose');
noseBBox     = step(noseDetector,videoFrame);

% Display results
figure; subplot(1,2,1);imshow(videoFrame); hold on;
rectangle('Position',noseBBox(1,:),'LineWidth',2,'EdgeColor',[1 1 0])
subplot(1,2,2); imshow(hueChannel); hold on;
rectangle('Position',noseBBox(1,:),'LineWidth',2,'EdgeColor',[0 1 0])

%% Set up tracking loop
% Create a tracker object.
tracker = vision.HistogramBasedTracker;

% Initialize the tracker histogram
initializeObject(tracker, hueChannel, noseBBox(1,:));
% figure;bar(tracker.ObjectHistogram);

%%
% Create a video player object for displaying video frames.
videoInfo    = info(videoFileReader);
videoPlayer  = vision.VideoPlayer('Position',[300 300 videoInfo.VideoSize+30]);
boxInserter  = vision.ShapeInserter('BorderColor','Custom',...
    'CustomBorderColor',[255 255 0]);

% Track the face over successive video frames until the video is finished.
while ~isDone(videoFileReader)
    videoFrame = step(videoFileReader);        % Extract the next video frame
    [hueChannel,~,~] = rgb2hsv(videoFrame);    % RGB -> HSV

    % Track using the Hue channel data
    [bbox, ~, score] = step(tracker, hueChannel);

    % Insert a bounding box around the object being tracked
    videoOut = step(boxInserter, videoFrame, bbox);
    % Display the annotated video frame using the video player object
    step(videoPlayer, videoOut);

end

%% Now try to reacquire the face when it reappears
release(tracker);
reset(videoPlayer);
initializeObject(tracker, hueChannel, noseBBox(1,:));

% Track the face over successive video frames until the video is finished.
reset(videoFileReader);
while ~isDone(videoFileReader)
    videoFrame = step(videoFileReader);        % Extract the next video frame
    [hueChannel,~,~] = rgb2hsv(videoFrame);    % RGB -> HSV

    % Track using the Hue channel data
    [bbox, ~, score] = step(tracker, hueChannel);

    % Determine if face is still in the scene using score value
    if score > 0.4
        videoOut = step(boxInserter, videoFrame, bbox);
    else
        videoOut = videoFrame;
        faceBBox = step(faceDetector,videoFrame);     % Find face again
        if ~isempty(faceBBox)
            release(noseDetector);
            noseBBox = step(noseDetector,imcrop(videoFrame,faceBBox(1,:)));
            noseBBox(1:2) = noseBBox(1:2) + faceBBox(1:2);
            release(tracker);
            initializeObject(tracker, hueChannel, noseBBox(1,:));
        end
    end
    % Display the annotated video frame using the video player object
    step(videoPlayer, videoOut);

end
% Release resources
release(videoFileReader);
release(videoPlayer);


text:
function [u, s, U_p, U_k, U_d] = pinHoleHmg(p,k,d)

% PINHOLEHMG Pin-hole camera model for HMG points, with optional radial distortion.
%   U = PINHOLEHMG(P) gives the projected pixel U of a homogeneous point P
%   in a canonical pin-hole camera, that is, with calibration parameters
%     u0 = 0
%     v0 = 0
%     au = 1
%     av = 1
%   It uses reference frames {RDF,RD} (right-down-front for the 3D world
%   points and right-down for the pixel), according to this scheme:
%
%         / z (forward)
%        /
%       +------- x                 +------- u
%       |                          |
%       |      3D : P=[x;y;z]      |     image : U=[u;v]
%       | y                        | v
%
%   U = PINHOLE(P,K) allows the introduction of the camera's calibration
%   parameters:
%     K = [u0 v0 au av]'.
%
%   U = PINHOLE(P,K,D) allows the introduction of the camera's radial
%   distortion parameters:
%     D = [K2 K4 K6 ...]'
%   so that the new pixel is distorted following the distortion equation:
%     U_D = U * (1 + K2*R^2 + K4*R^4 + ...)
%   with R^2 = sum(U.^2), being U the projected point in the image plane
%   for a camera with unit focal length.
%
%   [U,S] = PINHOLE(...) returns the depth S from the camera center.
%
%   If P is a points matrix, PINHOLE(P,...) returns a pixel matrix U and a
%   depths row-vector S. P, U and S are defined as
%     P = [P1 ... Pn];   Pi = [xi;yi;zi]
%     U = [U1 ... Un];   Ui = [ui;vi]
%     S = [S1 ... Sn]
%
%   [U,S,U_p,U_k,U_d] returns the Jacobians of U wrt P, K and D. It only
%   works for single points P=[x;y;z].
%
%   See also PINHOLE, INVPINHOLEHMG, PINHOLEIDP.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if nargout <= 2 % only pixel

    peuc = hmg2euc(p);

    switch nargin
        case 1
            [u, s] = pinHole(peuc);
        case 2
            [u, s] = pinHole(peuc,k);
        case 3
            [u, s] = pinHole(peuc,k,d);
    end

else % Jacobians

    if size(p,2) > 1
        error('Jacobians not available for multiple points')
    else

        [peuc, PEUC_p] = hmg2euc(p);

        switch nargin
            case 1
                [u, s, U_peuc] = pinHole(peuc);
                U_p = U_peuc*PEUC_p;

            case 2
                [u, s, U_peuc, U_k]  = pinHole(peuc,k);
                U_p = U_peuc*PEUC_p;

            case 3
                [u, s, U_peuc, U_k, U_d] = pinHole(peuc,k,d);
                U_p = U_peuc*PEUC_p;
        end

    end

end

return

%% jacobians
syms x y z t u0 v0 au av d2 d4 d6 real
p = [x;y;z;t];
k = [u0;v0;au;av];
d = [d2;d4;d6];

[u, s, U_p, U_k, U_d] = pinHoleHmg(p,k,d);

simplify(U_p - jacobian(u,p))
simplify(U_k - jacobian(u,k))
simplify(U_d - jacobian(u,d))

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [f,g] = spsurfun(y, z)
% SPSURFUN Evaluate the sparse grid interpolant at single point
%    IP = SPSURFUN(Y,Z)  Computes the interpolated value IP at
%    the single point [Y1, ..., YD] for the sparse grid
%    interpolant Z.
%
%    [IP,IPGRAD] = SPSURFUN(Y,Z)  Computes the interpolated
%    value IP and the gradient vector IPGRAD.
%
%    Example:
%       f = inline('x.^2 + y.^2 - 2.*z');
%       g1 = inline('2*x + 0*y + 0*z');
%       g2 = inline('2*y + 0*x + 0*z');
%       g3 = inline('-2  + 0*x + 0*y + 0*z');
%       z = spvals(f,3,[],spset('GridType','Chebyshev'));
%       [ip,ipgrad] = spsurfun([0.5, 0.2, 0.2], z)
%       f_exact = f(0.5, 0.2, 0.2)
%       g_exact = [g1(0.5, 0.2, 0.2); ...
%                  g2(0.5, 0.2, 0.2); ...
%                  g3(0.5, 0.2, 0.2)]
%
%    See also SPINTERP, SPVALS.
%
%    Note:
%       SPSURFUN is provided for conveniece to be used as an
%       alternative to SPINTERP, where the point Y to be
%       evaluated is given as a row or column vector. This
%       functional form is often adopted by multivariate
%       optimization algorithms (such as fminsearch) in
%       Matlab.
%       Note that this form allows the evaluation of the sparse
%       grid interpolant at a single point only. Therefore, It
%       is recommended to use SPINTERP instead if multiple
%       evaluations of the interpolant can be performed
%       simultaneously.

% Author : Andreas Klimke
% Version: 1.0
% Date   : September 8, 2006

% Change log:
% V1.0   : September 8, 2006
%          Initial version

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

if nargout <= 1
  f = spinterp(z, y);
elseif nargout == 2
 [f,g] = spinterp(z, y);
 g = g{1,1};
end


text:
function [ t1, t2, pi ] = lines_par_int_2d ( f1, g1, x1, y1, f2, g2, ...
  x2, y2 )

%*****************************************************************************80
%
%% LINES_PAR_INT_2D determines where two parametric lines intersect in 2D.
%
%  Discussion:
%
%    The parametric form of a line in 2D is:
%
%      X = X0 + F * T
%      Y = Y0 + G * T
%
%    We normalize by choosing F*F+G*G=1 and 0 <= F.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real F1, G1, X1, Y1, define the first parametric line.
%
%    Input, real F2, G2, X2, Y2, define the second parametric line.
%
%    Output, real T1, T2, the T parameters on the first and second
%    lines of the intersection point.
%
%    Output, real PI(2,1), the intersection point.
%
  det = f2 * g1 - f1 * g2;

  if ( det == 0.0 )
    t1 = 0.0;
    t2 = 0.0;
    pi = [];
  else
    t1 = ( f2 * ( y2 - y1 ) - g2 * ( x2 - x1 ) ) / det;
    t2 = ( f1 * ( y2 - y1 ) - g1 * ( x2 - x1 ) ) / det;
    pi(1,1) = x1 + f1 * t1;
    pi(2,1) = y1 + g1 * t1;
  end

  return
end


text:
function cookbook()
   %cookbook   Cookbook for TraceFilter
   %   Demonstrates the usage of TraceFilter

   disp('TraceFilter Cookbook')
   disp(' - not created, yet -');
end

author: geoscience-community-codes

repo: GISMO


text:
% AS-RPCA: Active Subspace: Towards Scalable Low-Rank Learning (Liu and Yan, 2012)
% process_video('RPCA', 'AS-RPCA', 'dataset/demo.avi', 'output/demo_AS-RPCA.avi');
lambda = 1/sqrt(min(size(M)));
[L,S] = as_rpca(M,lambda);

author: andrewssobral

repo: lrslibrary


text:
function p=internalpoint(v,aloop)
%
% p=internalpoint(v,aloop)
%
% imperical function to find an internal point
% of a planar polygon
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2008/04/08
%
% input:
%    v:     x,y,z coordinates of each node of the mesh
%    aloop:  input, a single vector separated by NaN, each segment
%             is a close-polygon consisted by node IDs
% output:
%    p:   output, [x y z] of an internal point of aloop
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

p=[];
nd=v(aloop,:);
boxfacet=find(sum(abs(diff(nd)))<1e-2); % find a flat loop
if(length(boxfacet))   % if the loop is flat along x/y/z dir
     bf=boxfacet(1);    % no degeneracy allowed
     idx=setdiff([1 2 3],bf);

     p0=(nd(1,:)+nd(2,:))/2;
     pvec=complex(p0(idx(1)),p0(idx(2)));
     vec=nd(2,:)-nd(1,:);
     vec=complex(vec(idx(1)),vec(idx(2)))*exp(i*pi/2)*(1e-5)/sqrt(sum(vec.*vec));
     testpt=[real(pvec+vec) imag(pvec+vec);real(pvec-vec) imag(pvec-vec)];
     in=inpolygon(testpt(:,1),testpt(:,2), nd(:,idx(1)),nd(:,idx(2)));
     p=testpt(find(in>0),:);
     p([bf,idx(1),idx(2)])=[nd(1,bf),p];
end

if(length(p)==0|length(p)==2)
    error('fail to find an internal point of curve');
end

repo: REST


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

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

% Last Modified by GUIDE v2.5 13-Jul-2011 13:55:22

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

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

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

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

% Update handles structure
guidata(hObject, handles);

set(hObject,'Position',[0.8706    0.8308    0.1263    0.1250]);

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

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

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

% --------------------------------------------------------------------
function menuFile_Callback(hObject, eventdata, handles)
return;

% --------------------------------------------------------------------
function menuLoad_Callback(hObject, eventdata, handles)
%
% Load a mesh file saved by mrVista/mrTools
%

uData = get(hObject,'userdata');

% It would be best to be able to read the list of currently open windows
% and assign a number based on that.
if ~checkfields(uData,'windowID'),
    nextWindow = 500;
    uData.windowID = nextWindow;
else
    nextWindow = max(uData.windowID(:))+1;
    uData.windowID(end+1) = nextWindow;
end

% Read a mesh starting in the Matlab working directory.
msh = mrmReadMeshFile(pwd);

if  isempty(msh), return; end

msh = meshSet(msh,'windowid',nextWindow);

% Display the mesh.  Turn off the origin lines.
msh = mrmInitMesh(msh);

mrmSet(msh,'hidecursor');
name = meshGet(msh,'name');
if isempty(name), [p,name] = fileparts(filename); end

mrmSet(msh,'title',name);
% Save the file information in the window object's user data.  We should be
% doing this for mrDiffusion, too.  And we should use this saved
% information when we ask about Editing windows.
set(hObject,'userdata',uData);

% Added in ability to open associated movie maker - RFB 06/2010
if (get(handles.OpenMovieMakerCheckbox, 'Value'))
    mrmMakeMovieGUI(nextWindow);
end

return;

% --------------------------------------------------------------------
function menuFileLoadMRD_Callback(hObject, eventdata, handles)
% Read a Matlab file containing a mrDiffusion mrMesh data set.  Then
% display the data.

% This should contain information about previously loaded data.  Here we
% store the new file information in uData.  --- Or we should ... not yet
% implemented.
uData = get(hObject,'userdata');

% Locate the file
persistent mrdPath;
persistent lastMshID;
curPath = pwd;

if(isempty(mrdPath) || isequal(mrdPath,0)), mrdPath = curPath; end
chdir(mrdPath);
[f, mrdPath] = uigetfile({'*.mat'}, 'Load MRD mesh...');
chdir(curPath);

if(isnumeric(f)), disp('Load Mesh (MRD) cancelled.'); return; end
fname = fullfile(mrdPath, f);

% Load the Mesh File.
d = load(fname);
msh = d.handles.mrMesh;
if isempty(lastMshID), lastMshID = msh.id;
else
    if msh.id == lastMshID
        lastMshID = lastMshID + 1;
        msh.id = lastMshID+1;
    end
end

% Make sure the mrMesh server is running
if ~mrmCheckServer('localhost'), mrmStart(msh.id,msh.host); end
mrmSet(msh,'refresh');

% Create the window and insert the data
msh = dtiInitMrMeshWindow(msh);
dtiMrMeshAddROIs(d.handles,msh);
dtiMrMeshAddFGs(d.handles,msh);
dtiMrMeshAddImages(d.handles,msh,d.origin,d.xIm,d.yIm,d.zIm);

% Make sure the window number is part of the title.  This way we can change the
% title later if we by an Edt | Set Window Title pull down.  Probably we
% should keep track of all the loaded meshes and the window numbers.
[p,fname,e] = fileparts(f);
str = mrmSet(msh,'windowtitle',sprintf('%.0f %s',meshGet(msh,'id'),fname));

% Added in ability to open associated movie maker - RFB 06/2010
if (get(handles.OpenMovieMakerCheckbox, 'Value'))
    mrmMakeMovieGUI(msh.id);
end

return;

% --------------------------------------------------------------------
function menuClose_Callback(hObject, eventdata, handles)
closereq;
return;

% --------------------------------------------------------------------
function menuQuit_Callback(hObject, eventdata, handles)
closereq;
return;

% --------------------------------------------------------------------
function menuEdit_Callback(hObject, eventdata, handles)

return;

% --------------------------------------------------------------------
function menuDelete_Callback(hObject, eventdata, handles)
return;

% --------------------------------------------------------------------
function menuHelp_Callback(hObject, eventdata, handles)
return;

% --------------------------------------------------------------------
function menuHelpGeneral_Callback(hObject, eventdata, handles)
% hObject    handle to menuHelpGeneral (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
helpMessage = 'Use mrmViewer to view Mesh Files.  Ordinarily these have been saved by mrVista''s 3D window viewer or mrDiffusion.';
hdl = mrMessage(helpMessage,'left',[0.7,0.85,0.15, 0.1]);

return

% --------------------------------------------------------------------
function menuHelpMovie_Callback(hObject, eventdata, handles)
% hObject    handle to menuHelpMovie (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
doc('mrmMakeMovieGUI');
return

% --- Executes on button press in btnLoad.
function btnLoad_Callback(hObject, eventdata, handles)
menuLoad_Callback(hObject, eventdata, handles)
return;

% --- Executes on button press in btnLoadMRD.
function btnLoadMRD_Callback(hObject, eventdata, handles)
menuFileLoadMRD_Callback(hObject, eventdata, handles)
return;

% --------------------------------------------------------------------
function menuOriginOff_Callback(hObject, eventdata, handles)
msh.host = 'localhost';
msh.id = whichMeshWindow;
msh.actor = 32;
mrmSet(msh,'hidecursor')
return;

% --------------------------------------------------------------------
function id = whichMeshWindow
%
%  This should be a check on uData values or something.

prompt={'Enter mesh window'};
def={'1'};
dlgTitle='Select mesh window';
lineNo=1;
answer=inputdlg(prompt,dlgTitle,lineNo,def);

if isempty(answer), id = [];
else id = str2num(answer{1});
end

return;

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

msh.host = 'localhost';
msh.id = whichMeshWindow;
dlgTitle='mrmViewer Edit|Title';
prompt={'Enter window title'};
lineNo=1;
def = {'Title'};
answer=inputdlg(prompt,dlgTitle,lineNo,def);

if isempty(answer), disp('User canceled.'); return;
else  mrmSet(msh,'title',answer{1}); end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = callsabr(k,f,limit)
% pricing a call using a density f and integration

    y = zeros(1,length(k));
    for j = 1:length(k)
        F = @(x) max(x-k(j),0) .* f(x);
        y(j) = quad(F,0,limit);
    end
end


text:
classdef SimpAllInterpolationImplicit2D < SimpAllInterpolationImplicit

   methods  (Access = public)

        function obj = SimpAllInterpolationImplicit2D(cParams)
            obj.init(cParams);
            obj.computeNstre();
            obj.computeSymbolicInterpolationFunctions();
            obj.dmu0 = obj.computeDmu0();
            obj.dmu1 = obj.computeDmu1();
            obj.dk0  = obj.computeDKappa0();
            obj.dk1  = obj.computeDKappa1();
        end

   end

    methods  (Access = protected)

        function [pMu,pKappa] = computePolarizationTensorAsMuKappa(eMatrix,eInclusion,nuMatrix,nuInclusion)
            coef = obj.compute2Dcoefficients(eMatrix,eInclusion,nuMatrix,nuInclusion);
            [p1,p2] = obj.computeP1P2(coef);
            [pMu,pKappa] = obj.computePkappaPmu(p1, p2);
        end

        function dmu0 = computeDmu0(obj)
            E1  = obj.matProp.E1;
            E0  = obj.matProp.E0;
            nu1 = obj.matProp.nu1;
            nu0 = obj.matProp.nu0;
            mu0 = obj.matProp.mu0;
            [pMu0,~] = obj.computePolarizationTensorAsMuKappa(E0,E1,nu0,nu1);
            dmu0     = -mu0*pMu0;
        end

        function dmu1 = computeDmu1(obj)
            E1  = obj.matProp.E1;
            E0  = obj.matProp.E0;
            nu1 = obj.matProp.nu1;
            nu0 = obj.matProp.nu0;
            mu1 = obj.matProp.mu1;
            [pMu1,~] = obj.computePolarizationTensorAsMuKappa(E1,E0,nu1,nu0);
            dmu1     = mu1*pMu1;
        end

        function dkappa0 = computeDKappa0(obj)
            E1     = obj.matProp.E1;
            E0     = obj.matProp.E0;
            nu1    = obj.matProp.nu1;
            nu0    = obj.matProp.nu0;
            kappa0 = obj.matProp.kappa0;
            [~,pKappa0] = obj.computePolarizationTensorAsMuKappa(E0,E1,nu0,nu1);
            dkappa0     = -kappa0*pKappa0;
        end

        function dkappa1 = computeDKappa1(obj)
            E1     = obj.matProp.E1;
            E0     = obj.matProp.E0;
            nu1    = obj.matProp.nu1;
            nu0    = obj.matProp.nu0;
            kappa1 = obj.matProp.kappa1;
            [~,pKappa1] = obj.computePolarizationTensorAsMuKappa(E1,E0,nu1,nu0);
            dkappa1     = kappa1*pKappa1;
        end

    end

    methods  (Access = protected, Static)

        function coef = compute2Dcoefficients(eMatrix,eInclusion,nuMatrix,nuInclusion)
            coef.a    = (1 + nuMatrix)/(1 - nuMatrix);
            coef.b    = (3 - nuMatrix)/(1 + nuMatrix);
            coef.gam  = eInclusion/eMatrix;
            coef.tau1 = (1 + nuInclusion)/(1 + nuMatrix);
            coef.tau2 = (1 - nuInclusion)/(1 - nuMatrix);
            coef.tau3 = (nuInclusion*(3*nuMatrix - 4) + 1)/(nuMatrix*(3*nuMatrix - 4) + 1);
        end

        function [p1,p2] = computeP(s)
            a    = s.a;
            b    = s.b;
            gam  = s.gam;
            tau1 = s.tau1;
            tau2 = s.tau2;
            tau3 = s.tau3;
            p1   = 1/(b*gam+tau1)*(1+b)*(tau1-gam);
            p2   = 0.5*(a-b)/(b*gam+tau1)*(gam*(gam-2*tau3)+tau1*tau2)/(a*gam+tau2);
        end

        function [pMu,pKappa] = computePKappaMu(p1, p2)
            pMu    = p1;
            pKappa = 2*p2 + p1;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function obj = GridData(geodata,obj,varargin)
% obj = GridData(geodata,obj,varargin);
% GridData: Uses the cell-averaged approach to interpolate the nodes
%           on the unstructured grid to the DEM.
%       Input : geodata - either a geodata class or a filename of the netcdf DEM
%                   obj - A msh class object containing:
%                     p - A vector of the locations of all the vertices
%                         in the mesh, length in number of nodes, nn x 2
%                     t - Matrix size ne x 3 of the element table of the mesh
%
%          K (optional) - vector of relevant nodes to search. This can
%                         significantly speed up the calculation by either
%                         only interpolating part of the mesh or
%                         intepolating whole mesh but in segments. Function
%                         automatically removes uncessary portions of DEM
%                         Example to call:
%                         K = find( obj.p(:,1) >= lon_min & ...
%                                   obj.p(:,2) <= lon_max);
%                         obj = GridData(fname,obj,'K',K);
%       type (optional) - type is either 'depth', 'slope' or 'all'
%                         'all' is the default (both slope and depth).
%                         'slope' to gets the gradients of DEM
%
%     interp (optional) - interp is either the normal griddedInterpolant
%                         options in MATLAB or is 'CA' (default). Note: 'CA'
%                         applies linear griddedInterpolant when the DEM
%                         and grid sizes are similar.
%
%          N (optional) - enlarge cell-averaging stencil by factor N (only
%                         relevant for CA interpolation method).
%                         default value N=1.
%
%        nan (optional) - 'fill' to fill in any NaNs everywhere
%                       - 'fillinside' to fill NaNs only in DEM extents
%
%   mindepth (optional) - ensure the minimum depth is bounded in the
%                         interpolated region
%
%   maxdepth (optional) - ensure the maximum depth is bounded in the
%                         interpolated region
%
%   ignoreOL (optional) - NaN overland data for more accurate seabed interpolation
%
%   lut (optional)      - A look up table (lut). See nlcd and ccap in
%                         datasets/ for examples
%
% slope_calc (optional) - 'rms' [default]: Compute cell-averaged slope based on root-mean-square
%                       - 'abs'          : Compute cell-averaged slope based on mean of absolute value
%
%      Output : obj     - A mesh class object with the following updated:
%               b       - The depth or vertical difference from the datum
%                         in the mesh if type is 'depth' or 'all'
%               bx      - The x direction gradient if type is 'slope' or
%                         'all'
%               by      - The y direction gradient if type is 'slope' or
%                         'all'
%
%	Author: William Pringle, and Keith Roberts CHL, Notre Dame University
%	Created: 2018-03-12
%
%   Edits by Keith Roberts, 2018-02-22 to avoid FillValue contaimation and
%   to only interpolate data inside extent of DEM.
%
%   Edits by Keith Roberts, 2018-10-18 to bound depth in interpolated
%   regions
%
%   Edits by Keith Roberts, 2019-4-4 to interpolate the floodplain
%
%% Test optional arguments
% default
nn = length(obj.p);
K  = (1:nn)';        % K is all of the grid
type = 'all';
interp = 'CA';
NaNs = 'ignore'; nanfill = false;
ignoreOL = 0 ;
N    = 1 ;
mindepth = -inf ;
maxdepth = +inf ;
rms_slope_calc = true;
slope_calc = 'rms';
lut = [] ;
if ~isempty(varargin)
    varargin=varargin{1} ;
    names = {'K','type','interp','nan','N','mindepth','maxdepth','ignoreOL','slope_calc','lut'};
    for ii = 1:length(names)
        ind = find(~cellfun(@isempty,strfind(varargin(1:2:end),names{ii})));
        if ~isempty(ind)
            if ii == 1
                K = varargin{ind*2};
                nn = length(K);
                if nn == 0
                   error('K index entry is empty!')
                end
            elseif ii == 2
                type = varargin{ind*2};
            elseif ii == 3
                interp = varargin{ind*2};
            elseif ii == 4
                NaNs = varargin{ind*2};
            elseif ii == 5
                N = varargin{ind*2} ;
            elseif ii ==6
                mindepth = varargin{ind*2} ;
            elseif ii ==7
                maxdepth = varargin{ind*2}  ;
            elseif ii ==8
                ignoreOL = varargin{ind*2} ;
            elseif ii ==9
                slope_calc = varargin{ind*2};
                if strcmp(slope_calc,'rms')
                   rms_slope_calc = true;
                elseif strcmp(slope_calc,'abs')
                   rms_slope_calc = false;
                else
                   error(['Invalid entry, ' temp ', for slope_calc'])
                end
            elseif ii == 10
                lut = varargin{ind*2};
            end
        end
    end
end

if ignoreOL
    disp('NaNing overland data before interpolating')
end
if N > 1
    disp(['Enlarging CA stencil by factor ',num2str(N)]) ;
end
if mindepth > -inf
    disp(['Bounding minimum depth to ',num2str(mindepth), ' meters.']) ;
end
if maxdepth < inf
    disp(['Bounding maximum depth to ',num2str(maxdepth), ' meters.']) ;
end

if strcmp(type,'slope')
    disp(['Computing cell-averaged slope using the ' slope_calc ' method'])
    if isempty(obj.b) || all(obj.b == 0)
        error(['You must have the bathymetry on the mesh before ' ...
               'calculating the slopes'])
    elseif sum(isnan(obj.b)) > 0
        warning(['There are some NaNs in the bathymetry that could ' ...
                 'affect the computation of the slopes'])
    end
end
if length(K) == 1
    if strcmp(type,'slope') || strcmp(type,'all')
        error('Cannot compute slopes on bathymetry with K of length 1')
    end
end

if strcmp(NaNs,'fill') || strcmp(NaNs,'fillinside')
    disp('Fill in NaNs using nearest neighbour interpolation.')
    nanfill = true;
end
if strcmp(NaNs,'fill')
   warning(['Note that nan "fill" option will try and put values everywhere ' ...
            'on mesh even outside of gdat/dem extents unless K logical is ' ...
            'set. Change nan to "fillinside" to avoid this.'])
end

if ~isempty(lut)
    disp('Using look up table...');
end

%% Let's read the LON LAT of DEM if not already geodata
flipUD = 0;
if ~isa(geodata,'geodata')
    [xvn, yvn, zvn] = getdemvarnames(geodata);
    DEM_XA = double(ncread(geodata,xvn));
    DEM_YA = double(ncread(geodata,yvn));
    DELTA_X = mean(diff(DEM_XA));
    DELTA_Y = mean(diff(DEM_YA));
    if DELTA_Y < 0
       flipUD = 1;
       DELTA_Y = -DELTA_Y;
    end
else
    DEM_XA = geodata.Fb.GridVectors{1};
    DEM_YA = geodata.Fb.GridVectors{2};
    DEM_Z  = -geodata.Fb.Values;
    if ignoreOL
       DEM_Z(DEM_Z <= 0) = NaN ;
    end
    DELTA_X = mean(diff(DEM_XA));
    DELTA_Y = mean(diff(DEM_YA));
    [DEM_X,DEM_Y] = ndgrid(DEM_XA,DEM_YA);
end

if max(DEM_XA) > 180
   lon_change = obj.p(:,1) < 0; lon_dir = 1;
elseif max(obj.p(:,1)) > 180
   lon_change = obj.p(:,1) > 180; lon_dir = -1;
else
   lon_change = false(length(obj.p),1); lon_dir = 0;
end
obj.p(lon_change,1) = obj.p(lon_change,1) + lon_dir*360;

% kjr edit 20180320
if length(K) == length(obj.p)
    % msh class may not have populated b, if this case populate with NaN to
    % detect issues
    if isempty(obj.b)
        obj.b = obj.p(:,1)*NaN;
    end
    if ~strcmp(NaNs,'fill')
        outside = obj.p(:,1) < min(DEM_XA)-DELTA_X | ...
                  obj.p(:,1) > max(DEM_XA)+DELTA_X | ...
                  obj.p(:,2) < min(DEM_YA)-DELTA_Y | ...
                  obj.p(:,2) > max(DEM_YA)+DELTA_Y;
        K(outside) = [];
        if isempty(K)
            warning('no mesh vertices contained within DEM bounds'); return;
        end
    end
end

% If the length of DEM_X and DEM_Y is too large then let's break it up by
% using K
% test
lon_min = min(obj.p(K,1));
lon_max = max(obj.p(K,1));
lat_min = min(obj.p(K,2));
lat_max = max(obj.p(K,2));

I = find(DEM_XA >= lon_min & DEM_XA <= lon_max);
J = find(DEM_YA >= lat_min & DEM_YA <= lat_max);
% single or integer array
L = length(I)*length(J)*4e-9;
% larger than 2 GB
if L > 2
    % divide by latitude
    times = ceil(L/2); lat_l = lat_min; dl = (lat_max - lat_min)/times;
    for t = 1:times
        lat_r = lat_l + dl;
        if t == times
            lat_r = lat_max;
        end
        KK{t} = K(obj.p(K,2) >= lat_l & obj.p(K,2) <= lat_r);
        lat_l = lat_r;
    end
else
    times = 1;
    KK{1} = K;
end
% This deletes any elements straddling the -180/180 boundary
xt = [obj.p(obj.t(:,1),1) obj.p(obj.t(:,2),1) ...
      obj.p(obj.t(:,3),1) obj.p(obj.t(:,1),1)];
dxt = diff(xt,[],2);
tt = obj.t;
tt(abs(dxt(:,1)) > 180 | abs(dxt(:,2)) > 180 |  abs(dxt(:,3)) > 180,:) = [];

% Do this once;
vtoe_o = VertToEle(tt); %find connecting elements to each node
pmid = squeeze(mean(reshape(obj.p(tt,:),[],3,2),2)); % get mid points of elements
pmid(end+1,:) = NaN;

K_o = K; % save the original K
disp(['Looping over the DEM ' num2str(times) ' time(s)'])
for t = 1:times
K = KK{t};
%% Get the grid size for each node
% Get the subset
vtoe = vtoe_o(:,K);
% make sure when vtoe is zero we just return a NaN
vtoe(vtoe == 0) = length(tt) + 1;
% the connecting element centers
pcon = reshape(pmid(vtoe,:),size(vtoe,1),[],2);

% delta is max and min bounds of the connecting element centers (or if only
% one connecting element then do difference with the vertex itself
if length(K) == 1
    pcon_max = max(max(squeeze(pcon)),2*obj.p(K,:)-max(squeeze(pcon)));
    pcon_min = min(min(squeeze(pcon)),2*obj.p(K,:)-min(squeeze(pcon)));
else
    pcon_max = max(squeeze(max(pcon,[],1)),2*obj.p(K,:)-squeeze(min(pcon,[],1)));
    pcon_min = min(squeeze(min(pcon,[],1)),2*obj.p(K,:)-squeeze(max(pcon,[],1)));
end
delta = pcon_max - pcon_min;

% kjr 10/17/2018 enlarge the cell-averaging stencil by a factor N
pcon_max = N*pcon_max+(1-N)*obj.p(K,:) ;
pcon_min = N*pcon_min+(1-N)*obj.p(K,:) ;

%% Now read the bathy (only what we need)
if ~isa(geodata,'geodata')
    BufferL = max(delta,[],1);
    lon_min = min(obj.p(K,1)) - BufferL(1);
    lon_max = max(obj.p(K,1)) + BufferL(1);
    lat_min = min(obj.p(K,2)) - BufferL(2);
    lat_max = max(obj.p(K,2)) + BufferL(2);
    I = find(DEM_XA > lon_min & DEM_XA < lon_max);
    J = find(DEM_YA > lat_min & DEM_YA < lat_max);
    if exist('DEM_X','var')
        clear DEM_X DEM_Y DEM_Z
    end
    [DEM_X,DEM_Y] = ndgrid(DEM_XA(I),DEM_YA(J));
    DEM_Z = single(ncread(geodata,zvn,...
                   [I(1) J(1)],[length(I) length(J)]));
    if flipUD
       % handle DEMS from packed starting from the bottom left
       DEM_YA = flipud(DEM_YA) ;
       DEM_Y = fliplr(DEM_Y);
       DEM_Z = fliplr(DEM_Z);
    end

    % make into depths (ADCIRC compliant)
    DEM_Z = -DEM_Z;

    if ignoreOL
       DEM_Z(DEM_Z <= 0) = NaN ;
    end
end

% bound all depths below mindepth
DEM_Z(DEM_Z < mindepth) = mindepth ;

% bound all depths above maxdepth
DEM_Z(DEM_Z > maxdepth) = maxdepth ;

% if look up table is passed
if ~isempty(lut)
    DEM_Z(isnan(DEM_Z) | DEM_Z == 0)=length(lut); % <--default value to NaN
    % Convert to Mannings
    DEM_Z = lut(abs(round(DEM_Z)));
end
%% Make the new bx, by, b and calculate gradients if necessary
if strcmp(type,'slope') || strcmp(type,'all')
    bx = NaN(length(K),1);
    by = NaN(length(K),1);
    % Get the dx and dy of the dem in meters
    DELTA_X1 = DELTA_X*111e3;
    DELTA_X1 = DELTA_X1*cosd(DEM_Y(1,:)); % for gradient function
    DELTA_Y1 = DELTA_Y*111e3;
    if exist('DEM_ZY','var')
        clear DEM_ZY DEM_ZX
    end
    % Calculate the gradient of the distance function.
    [DEM_ZY,DEM_ZX] = EarthGradient(DEM_Z,DELTA_Y1,DELTA_X1);
    % New method of averaging the absolute values
    DEM_ZY = abs(DEM_ZY); DEM_ZX = abs(DEM_ZX);
end
if strcmp(type,'depth') || strcmp(type,'all')
    b = NaN(length(K),1);
end

%% Calculate N - number of DEM grids to average - for each node
if strcmp(interp,'CA')
    [~,IDXX,IDXY] = FindLinearIdx(obj.p(K,1),obj.p(K,2),DEM_X,DEM_Y);
    [~,IDXR,IDXT] = FindLinearIdx(pcon_max(:,1),pcon_max(:,2),DEM_X,DEM_Y);
    [~,IDXL,IDXB] = FindLinearIdx(pcon_min(:,1),pcon_min(:,2),DEM_X,DEM_Y);
    % make sure inside box
    % y
    high = DEM_Y(IDXR(1),IDXT)' > pcon_max(:,2);
    IDXT(high) = max(1,IDXT(high) - 1);
    low = DEM_Y(IDXL(1),IDXB)' < pcon_min(:,2);
    IDXB(low) = min(size(DEM_Y,2),IDXB(low) + 1);
    % x
    high = DEM_X(IDXR,IDXT(1)) > pcon_max(:,1);
    IDXR(high) = max(1,IDXR(high) - 1);
    low = DEM_X(IDXL,IDXB(1)) < pcon_min(:,1);
    IDXL(low) = min(size(DEM_X,1),IDXL(low) + 1);
    % Make sure no negative or positive Nx, Ny
    I = IDXL > IDXR;
    IDXL(I) = IDXX(I); IDXR(I) = IDXX(I);
    I = IDXB > IDXT;
    IDXB(I) = IDXY(I); IDXT(I) = IDXY(I);
    % The span of x and y
    Nx = IDXR - IDXL;
    Ny = IDXT - IDXB;

    disp(['Performing cell-averaging with Ny(max) ' num2str(max(Ny)) ...
          ' and Nx(max)' num2str(max(Nx))])

    % Average for the depths
    if strcmp(type,'depth') || strcmp(type,'all')
        for ii = 1:length(K)
            pts = reshape(DEM_Z(IDXL(ii):IDXR(ii),...
                                IDXB(ii):IDXT(ii)),[],1);
            b(ii) = mean(pts,'omitnan');
        end
        if sum(~isnan(b)) == 0
            warning('All depths were NaNs. Doing nothing and returning');
            return
        end
        % Try and fill in the NaNs
        if nanfill
            if sum(isnan(b)) > 0
                localcoord = obj.p(K,:);
                [KI, ~] = ourKNNsearch(localcoord(~isnan(b),:),',localcoord(isnan(b)',1);
                %KI = knnsearch(localcoord(~isnan(b),:),localcoord(isnan(b),:));
                bb = b(~isnan(b));
                b(isnan(b)) = bb(KI); clear bb localcoord
            end
        end
    end

    % RMS for the slopes
    if strcmp(type,'slope') || strcmp(type,'all')
        for ii = 1:length(K)
            pts = reshape(DEM_ZX(IDXL(ii):IDXR(ii),...
                                 IDXB(ii):IDXT(ii)),[],1);
            if rms_slope_calc
               bx(ii) = sqrt(mean(pts.^2,'omitnan'));
            else
               bx(ii) = mean(abs(pts),'omitnan');
            end
            pts = reshape(DEM_ZY(IDXL(ii):IDXR(ii),...
                                 IDXB(ii):IDXT(ii)),[],1);
            if rms_slope_calc
               by(ii) = sqrt(mean(pts.^2,'omitnan'));
            else
               by(ii) = mean(abs(pts),'omitnan');
            end
        end
        if nanfill
            % Try and fill in the NaNs
            if ~isempty(find(isnan(bx),1))
                localcoord = obj.p(K,:);
                [KI, ~] = ourKNNsearch(localcoord(~isnan(bx),:)',localcoord(isnan(bx),:)',1);
                %KI = knnsearch(localcoord(~isnan(bx),:),localcoord(isnan(bx),:));
                bb = bx(~isnan(bx));
                bx(isnan(bx)) = bb(KI); clear bb localcoord
            end
            if ~isempty(find(isnan(by),1))
                localcoord = obj.p(K,:);
                [KI, ~] = ourKNNsearch(localcoord(~isnan(by),:)',localcoord(isnan(by),:)',1);
                %KI = knnsearch(localcoord(~isnan(by),:),localcoord(isnan(by),:));
                bb = by(~isnan(by));
                by(isnan(by)) = bb(KI); clear bb localcoord
            end
        end
    end
else
    %% Get gridded interpolant for non-CA interp
    disp(['Using griddedInterpolant with ' interp 'interp option'])
    if strcmp(type,'slope') || strcmp(type,'all')
        Fx = griddedInterpolant(DEM_X,DEM_Y,DEM_ZX,interp);
        Fy = griddedInterpolant(DEM_X,DEM_Y,DEM_ZY,interp);
        bx = Fx(obj.p(K,1),obj.p(K,2));
        by = Fy(obj.p(K,1),obj.p(K,2));
    end
    if strcmp(type,'depth') || strcmp(type,'all')
        F = griddedInterpolant(DEM_X,DEM_Y,DEM_Z,interp,'none');
        b = F(obj.p(K,1),obj.p(K,2));
    end
end
%% Put in the global msh class
if strcmp(type,'depth') || strcmp(type,'all')
    if isempty(obj.b)
        obj.b = zeros(length(obj.p),1);
    end
    % Only overwrite non-nan values
    obj.b(K(~isnan(b))) = b(~isnan(b));
end
if strcmp(type,'slope') || strcmp(type,'all')
    if isempty(obj.bx)
        obj.bx = zeros(length(obj.p),1);
    end
    obj.bx(K(~isnan(bx))) = bx(~isnan(bx));
    % Put in the global array
    if isempty(obj.by)
        obj.by = zeros(length(obj.p),1);
    end
    obj.by(K(~isnan(by))) = by(~isnan(by));
end
end
% New method of averaging asbolute values of slope before multiplying
% by the sign of the slope on unstructured grid
if strcmp(type,'slope') || strcmp(type,'all')
    [Hx,Hy] = Unstruc_Bath_Slope( tt,obj.p(:,1),obj.p(:,2),obj.b);
    obj.bx(K_o) = sign(Hx(K_o)).*obj.bx(K_o);
    obj.by(K_o) = sign(Hy(K_o)).*obj.by(K_o);
end
obj.p(lon_change,1) = obj.p(lon_change,1) - lon_dir*360;
%EOF
end

function [xvn, yvn, zvn] = getdemvarnames(fname)
    % Define well-known variables for longitude and latitude
    % coordinates in Digital Elevation Model NetCDF file (CF
    % compliant).
    xvn = []; yvn = []; zvn = [];
    wkv_x = {'x','Longitude','longitude','lon','lon_z'} ;
    wkv_y = {'y','Latitude', 'latitude','lat','lat_z'} ;
    finfo = ncinfo(fname);
    for ii = 1:length(finfo.Variables)
        if ~isempty(xvn) && ~isempty(yvn) && ~isempty(zvn); break; end
        if length(finfo.Variables(ii).Size) == 1
            if isempty(xvn) && ...
                    any(strcmp(finfo.Variables(ii).Name,wkv_x))
                xvn = finfo.Variables(ii).Name;
            end
            if isempty(yvn) && ...
                    any(strcmp(finfo.Variables(ii).Name,wkv_y))
                yvn = finfo.Variables(ii).Name;
            end
        elseif length(finfo.Variables(ii).Size) == 2
            if isempty(zvn)
                zvn = finfo.Variables(ii).Name;
            end
        end
    end
    if isempty(xvn)
        error('Could not locate x coordinate in DEM') ;
    end
    if isempty(yvn)
        error('Could not locate y coordinate in DEM') ;
    end
    if isempty(zvn)
        error('Could not locate z coordinate in DEM') ;
    end
end

author: CHLNDDEV


text:
function [ml] = cm32ml(cm3)
% Convert volume from cubic centimeters to milliliters.
% Chad Greene 2012
ml = cm3;


text:
function [params] = sv_calcMcNodes(params,nNodeStart, nNodeEnd)
% function [params] = sv_calcMcNodes(params,nNodeStart, nNodeEnd);
% ----------------------------------------------------------------
% Calculation of magnitude of completeness specifying the nodes to be calculated for distributing on different CPUs
%
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bRandom            Perform random simulation (=1) or real calculation (=0)
%   params.nCalculation       Number of random simulations
%   params.bMap               Calculate a map (=1) or a cross-section (=0)
%   params.bNumber            Use constant number (=1) or constant radius (=0)
%   params.nNumberEvents      Number of earthquakes if bNumber == 1
%   params.fMaxRadius         Maximum Radius using a constant number of events; works only with bNumber == 1
%   params.fRadius            Radius of gridnode if bNumber == 0
%   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
%   params.fMinMag            Lower limit of magnitude range for testing
%   params.fMaxMag            Upper limit of magnitude range for testing
%   params.bTimePeriod        Calculate seismicity difference for 2 periods (0) until start and end of catalog or
%                             a specific time period before and after fSplitTime (1)
%   params.fTimePeriod        Length of time periods
%   params.bTstart            Check for starting time of temporal mapping
%   params.fTstart            Starting time for temporal mapping
%   params.bBstnum            Check for boostrap sampling
%   params.fBstnum            Number of bootstrap samples
%   params.fBinning           Bin size for magnitude binning
%   params.sComment           Comment on calculation

% Output parameters:
%   Same as input parameters including
%   params.mValueGrid         Matrix of calculated values
%   params.vcsGridNames       Names of parameters calculated
%   Check sv_NodeCalcMc.m for a list of variables!!
%
% J. Woessner; woessner@seismo.ifg.ethz.ch
% last update: 27.03.03

global bDebug;
if bDebug
    report_this_filefun(mfilename('fullpath'));
end

% Initialize
vResults = [];
params.sComment = [];
if isempty(params.fBinning)
    params.fBinning = 0.1;
end

% Determine time period of catalog
params.fTminCat = min(params.mCatalog(:,3));
params.fTmaxCat = max(params.mCatalog(:,3));
% Adjust to decimal years
fTimePeriod =params.fTimePeriod/365;

% Init result matrix
mValueGrid_ = [];

% Temporary saving the original catalog
mCatalog = params.mCatalog;

% Check for bootstrapping or not
% ------------------------------
% Case of calculations with bootstrapping
if (params.bBstnum == 1)
    % Loop over time
    fTstart = params.fTstart;
    while fTstart < params.fTmaxCat
        mValueGrid_ = [];
        params.mCatalog = mCatalog;
        % Create Indices to catalog and select quakes in time period
        vSel = (fTstart <= params.mCatalog(:,3) & params.mCatalog(:,3) < fTstart+fTimePeriod);
        params.mCatalog = params.mCatalog(vSel,:);
        [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
         params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
        % Loop over all grid nodes
        hWaitbar1 = waitbar(0,'Calculating nodes...');
        set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
        for nNode_ = nNodeStart:nNodeEnd
            % Create node catalog
            mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
            % Check for constant number of events calculations
            if (params.nGriddingMode == 0)
                [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
            end
            [nX,nY] = size(mNodeCatalog_);
            if (nX < params.nMinimumNumber)
                mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN nX NaN NaN NaN NaN];
            else
                [rCalcNodeResult_] = sv_NodeCalcMc(params,mNodeCatalog_);
                % Store the results
                mValueGrid_= [mValueGrid_; rCalcNodeResult_.fMc_max rCalcNodeResult_.fMc_90 rCalcNodeResult_.fMc_95 rCalcNodeResult_.fMc_com...
                        rCalcNodeResult_.fMc_EMR rCalcNodeResult_.fMc_shi rCalcNodeResult_.fMc_Bst...
                        rCalcNodeResult_.fStd_Mc rCalcNodeResult_.fBvalue_Bst rCalcNodeResult_.fStd_B...
                        rCalcNodeResult_.fAvalue_Bst rCalcNodeResult_.fStd_A nX...
                        rCalcNodeResult_.bH rCalcNodeResult_.fPval rCalcNodeResult_.bH_Bst rCalcNodeResult_.fPval_Bst];
            end; % End of if on nNode_
            if rem(nNode_,floor(length(params.mPolygon(:,1))/10)) == 0
                waitbar(nNode_/length(params.mPolygon(:,1)))
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Temporary saving
                params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
                    'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)', 'Mc(Bst-mean)', 'Mc(Bst-2nd-moment)',...
                    'Mc(Bst-b)', 'Mc(b_2nd-moment)','Mc(Bst-a)', 'Mc(a_2nd-moment)','Number of events',...
                    'H(KST)','P(KST)','H(KST_Bst)','P(KST_Bst)'));
                params.mValueGrid = mValueGrid_;
                % Add parameter to params.sComment
                if  params.nGriddingMode == 0;   % Constant number
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                            num2str(params.fMaxRadius) ' km'];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                            '_Nmin_' num2str(params.nMinimumNumber) '_Node' num2str(nNode_) '.mat'], 'vResults');
                elseif params.nGriddingMode == 1;   % Constant radius
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                else  % Rectangle mode
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                            ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                            '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                end; % END of params.nGriddingmode
                vResults =[];
            end; % End updating waitbar
        end; % for nNode
        close(hWaitbar1);
        % Parameter description
        params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
            'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)', 'Mc(Bst-mean)', 'Mc(Bst-2nd-moment)',...
            'Mc(Bst-b)', 'Mc(b_2nd-moment)','Mc(Bst-a)', 'Mc(a_2nd-moment)','Number of events',...
                    'H(KST)','P(KST)','H(KST_Bst)','P(KST_Bst)'));
        params.mValueGrid = mValueGrid_;
        % Add parameter to params.sComment
        if  params.nGriddingMode == 0;   % Constant number
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        elseif params.nGriddingMode == 1;   % Constant radius
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        else  % Rectangle mode
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        end
        vResults =[];
        fTstart = fTstart+fTimePeriod;
    end; % End of while fTstart

    % Case of no bootstrapping
else
    % Loop over time
    fTstart = params.fTstart;
    while fTstart < params.fTmaxCat
        mValueGrid_ = [];
        params.mCatalog = mCatalog;
        % Create Indices to catalog and select quakes in time period
        vSel = (fTstart <= params.mCatalog(:,3) & params.mCatalog(:,3) < fTstart+fTimePeriod);
        params.mCatalog = params.mCatalog(vSel,:);
        [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
         params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
        % Loop over all grid nodes
        hWaitbar1 = waitbar(0,'Calculating nodes...');
        set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
        for nNode_ = nNodeStart:nNodeEnd
            % Create node catalog
            mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
            % Check for constant number of events calculations
            if (params.nGriddingMode == 0)
                [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
            end
            [nX,nY] = size(mNodeCatalog_);
            if (nX < params.nMinimumNumber)
                mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN nX];
            else
                [rCalcNodeResult_] = sv_NodeCalcMc(params,mNodeCatalog_);
                mValueGrid_= [mValueGrid_; rCalcNodeResult_.fMc_max rCalcNodeResult_.fMc_90 rCalcNodeResult_.fMc_95 rCalcNodeResult_.fMc_com...
                        rCalcNodeResult_.fMc_EMR rCalcNodeResult_.fMc_shi nX];
            end; % End of if on nX
            if rem(nNode_,500) == 0
                waitbar(nNode_/length(params.mPolygon(:,1)))
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Temporary saving
                params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
                    'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)','Number of events'));
                params.mValueGrid = mValueGrid_;
                % Add parameter to params.sComment
                if  params.nGriddingMode == 0;   % Constant number
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                            num2str(params.fMaxRadius) ' km'];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                            '_Nmin_' num2str(params.nMinimumNumber) '_Node' num2str(nNode_) '.mat'], 'vResults');
                elseif params.nGriddingMode == 1;   % Constant radius
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                else  % Rectangle mode
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                            ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                            '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                end; % END of params.nGriddingmode
                vResults =[];
            end; % End updating waitbar
        end; % for nNode
        close(hWaitbar1);
        % Parameter description
        params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
            'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)'));
        params.mValueGrid = mValueGrid_;
        % Add parameter to params.sComment
        if  params.nGriddingMode == 0;   % Constant number
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        elseif params.nGriddingMode == 1;   % Constant radius
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        else  % Rectangle mode
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period ' num2str(params.fTimePeriod) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        end
        vResults =[];
        fTstart = fTstart+fTimePeriod;
    end; % End of while fTstart
end; % END of if params.bBst

author: CelsoReyes


text:
% This function executes IPOPT to find the maximum likelihood solution to
% least squares regression with L1 regularization or the "Lasso". The inputs
% are the data matrix A (in which each row is an example vector), the vector
% of regression outputs y, and the penalty parameter lambda, a number
% greater than zero. The output is the estimated vector of regression
% coefficients.
%
% Copyright (C) 2008 Peter Carbonetto. All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% Author: Peter Carbonetto
%         Dept. of Computer Science
%         University of British Columbia
%         September 18, 2008
function w = lasso (A, y, lambda)

  % Get the number of examples (n) and the number of regression
  % coefficients (m).
  [n m] = size(A);

  % The starting point.
  x0 = { zeros(m,1)
         ones(m,1) };

  % The constraint functions are bounded from below by zero.
  options.cl = zeros(2*m,1);
  options.cu = repmat(Inf,2*m,1);

  % Set up the auxiliary data.
  options.auxdata = { n m A y lambda };

  % Set the IPOPT options.
  options.ipopt.jac_d_constant   = 'yes';
  options.ipopt.hessian_constant = 'yes';
  options.ipopt.mu_strategy      = 'adaptive';
  options.ipopt.max_iter         = 100;
  options.ipopt.tol              = 1e-8;

  % The callback functions.
  funcs.objective         = @objective;
  funcs.constraints       = @constraints;
  funcs.gradient          = @gradient;
  funcs.jacobian          = @jacobian;
  funcs.jacobianstructure = @jacobianstructure;
  funcs.hessian           = @hessian;
  funcs.hessianstructure  = @hessianstructure;

  % Run IPOPT.
  [x info] = ipopt(x0,funcs,options);
  w        = x{1};

% ------------------------------------------------------------------
function f = objective (x, auxdata)
  [n m A y lambda] = deal(auxdata{:});
  [w u] = deal(x{:});
  f     = norm(y - A*w)^2/2 + lambda*sum(u);

% ------------------------------------------------------------------
function c = constraints (x, auxdata)
  [w u] = deal(x{:});
  c     = [ w + u; u - w ];

% ------------------------------------------------------------------
function g = gradient (x, auxdata)
  [n m A y lambda] = deal(auxdata{:});
  w = x{1};
  g = { -A'*(y - A*w)
        repmat(lambda,m,1) };

% ------------------------------------------------------------------
function J = jacobianstructure (auxdata)
  m = auxdata{2};
  I = speye(m);
  J = [ I I
        I I ];

% ------------------------------------------------------------------
function J = jacobian (x, auxdata)
  m = auxdata{2};
  I = speye(m);
  J = [  I  I
        -I  I ];

% ------------------------------------------------------------------
function H = hessianstructure (auxdata)
  m = auxdata{2};
  H = [ tril(ones(m))  zeros(m)
          zeros(m)     zeros(m) ];
  H = sparse(H);

% ------------------------------------------------------------------
function H = hessian (x, sigma, lambda, auxdata)
  [n m A y lambda] = deal(auxdata{:});
  H = [ tril(A'*A)  zeros(m)
         zeros(m)   zeros(m) ];
  H = sparse(sigma * H);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%compute orientation of tailsmcentralvector

function [data,units]=compute_tailsmcentralang(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
tailsmcentralang=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    tailsmcentralang{1,i}=bsxfun(@atan2,trx(larva).ycentral_mm-trx(larva).ytailsm_mm,trx(larva).xcentral_mm-trx(larva).xtailsm_mm);
end
units=parseunits('rad');
data=tailsmcentralang;

author: kristinbranson

repo: JAABA


text:
function [ R ] = V2R( V )
%V2R converts a 1x3 angle-axis vector into a 3x3 rotation matrix
%   Inputs -
%   V - 1x3 vector of form [rx,ry,rz] where rx,ry,rz is an angle-axis
%   representation of the angle where the unit vector representing the axis
%   has been multipled by the angle of rotation about it
%
%   Outputs -
%   R - a standard 3x3 transformation matrix

validateattributes(V, {'numeric'},{'size',[1,3]});

V = double(V(:));

s = norm(V);
if(s == 0)
    R = eye(3);
    return;
end

V = [V/s; s];
R = vrrotvec2mat(V);

end

author: ZacharyTaylor

repo: Camera-to-Arm-Calibration


text:
function title = p34_title ( )

%*****************************************************************************80
%
%% P34_TITLE returns a title for problem 34.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Six-Hump Camel-Back Polynomial';

  return
end


text:
function linplus_test155 ( )

%*****************************************************************************80
%
%% TEST155 tests R8BTO_MXV, R8BTO_VXM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  l = 3;
  m = 2;
  n = m * l;
%
%  I don't know how MATLAB wants me to enter a triply-dimensioned array.
%  I've tried several different approaches with no success, so
%  here's the clodhopper's method:
%
  a = zeros ( m, m, 2*l-1 );

  a(1,1,1) = 1.0E+00;
  a(2,1,1) = 5.0E+00;
  a(1,2,1) = 2.0E+00;
  a(2,2,1) = 5.0E+00;
  a(1,1,2) = 3.0E+00;
  a(2,1,2) = 6.0E+00;
  a(1,2,2) = 4.0E+00;
  a(2,2,2) = 6.0E+00;
  a(1,1,3) = 5.0E+00;
  a(2,1,3) = 7.0E+00;
  a(1,2,3) = 6.0E+00;
  a(2,2,3) = 7.0E+00;
  a(1,1,4) = 7.0E+00;
  a(2,1,4) = 8.0E+00;
  a(1,2,4) = 8.0E+00;
  a(2,2,4) = 8.0E+00;
  a(1,1,5) = 9.0E+00;
  a(2,1,5) = 9.0E+00;
  a(1,2,5) = 0.0E+00;
  a(2,2,5) = 9.0E+00;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST155\n' );
  fprintf ( 1, '  For a real block Toeplitz matrix,\n' );
  fprintf ( 1, '  R8BTO_MXV computes A * x.\n' );
  fprintf ( 1, '  R8BTO_VXM computes A''* x.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Block order M =  %d\n', m );
  fprintf ( 1, '  Block number L = %d\n', l );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  r8bto_print ( m, l, a, '  The block Toeplitz matrix:' );

  x = r8ge_indicator ( m, l );

  r8ge_print ( m, l, x, '  The matrix x:' );

  b = r8bto_mxv ( m, l, a, x );

  r8ge_print ( m, l, b, '  The product A*x:' );

  b = r8bto_vxm ( m, l, a, x );

  r8ge_print ( m, l, b, '  The product A''*x:' );

  return
end