text:
% moviethresh() - threshold cell arrays that will be sent to brainmovie
%
% Usage:
%   >> data = moviethresh( data, threshold, continous, dimcont);
%
% Inputs:
%   data       - input cell array containing data arrays
%   threshold  - absolute threshold (0=none)
%   continuous - impair number of data points that must be sequentially
%                organized (3,5,7,9 ...)
%   dimcont    - dimention (1=rows, 2=columns) of the data arrays to
%                consider for sequential continuity (default: 2)
%
% Outputs:
%   data       - modified cell array
%
% Example:
%   d = moviethresh( data, 0.1, 3, 2);
%   % remove (i.e) any (absolute) value in the array contained in data
%   % that are below 0.1 and remove any value that is flanked by
%
% See also: BRAINMOVIE, TIMECROSSF

% arno@salk.edu, Arnaud Delorme, CNL / Salk Institute, 2001

% This program is free software; you can redistribute it and/or
% modify it.
% 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.

function data = moviethresh( data, threshold, cnt, dimcontinous);

if nargin < 4
	help moviethresh;
end;

% scan datas
% ----------
nbremoved = 0;
total     = 0;
for index = 1:length(data(:))
	tmp = data{ index };

	% remove values below threshold
	% -----------------------------
	tmp( find(abs(tmp) < threshold) ) = 0;

	% %%%%%%%%%%%%%%%%% TRANSPOSE if dimention 2
	if dimcontinous == 2, tmp = tmp'; end;

	% remove values that are sparse along dimention 1
	% -----------------------------------------------
	s = size(tmp, 1);
	for dimcont = 1:s
		nbcontrigth = zeros(1, size(tmp,2));
		nbcontleft  = zeros(1, size(tmp,2));
		for dim2 = 1:size(tmp,2)
			for ind=1:cnt 									% cont elements on the left
				if dimcont-ind > 0
					if tmp(dimcont-ind , dim2) ~= 0, nbcontleft(dim2)  = nbcontleft(dim2)  + 1;
					else, break; end;
				else, break; end;
			end;
			for ind=1:cnt  									% count elements on the right
				if dimcont+ind <= s
					if tmp(dimcont+ind , dim2) ~= 0, nbcontrigth(dim2)  = nbcontrigth(dim2)  + 1;
					else, break; end;
				else, break; end;
			end;
		end;

		total     = total + sum(tmp(dimcont,:) ~= 0); % non nul values

		nbremoved = nbremoved + sum( (tmp(dimcont,:) .*  ((nbcontleft+nbcontrigth) < cnt)) > 0);
		tmp(dimcont,:) = tmp(dimcont,:) .* ( (nbcontleft+nbcontrigth) >= cnt);
		%fprintf( 'Pos %d: %s ** %s\n', dimcont, int2str( nbcontleft), int2str( nbcontrigth));
	end;

	% %%%%%%%%%%%%%%%%% TRANSPOSE if dimention 2
	if dimcontinous == 2, tmp = tmp'; end;

	data{ index } = tmp;

end;
fprintf('%d removed values out of %d\n', nbremoved, total)

%debug fprintf('\tval:%f\tdestval:%f\tdimcont:%d dim2:%d ind:%d left:%d\n', tmp(dimcont, dim2), tmp(dimcont-ind , dim2), dimcont, dim2, ind, nbcontleft(dim2));

author: PatternRecognition

repo: OpenBMI


text:
function [rVectECI, vVectECI] = getInertialVectFromFixedFrameVect(ut, rVectECEF, bodyInfo, vVectECEF)
%getInertialVectFromFixedFrameVect Summary of this function goes here
%   Detailed explanation goes here
    numElems = length(ut);

    spinAngle = getBodySpinAngle(bodyInfo, ut);

    cSA = reshape(cos(spinAngle),1,1,numElems);
    sSA = reshape(sin(spinAngle),1,1,numElems);
    zero = zeros(1,1,numElems);
    one = zero + 1;

%     R = [cos(spinAngle) -sin(spinAngle) 0;
%          sin(spinAngle) cos(spinAngle) 0;
%          0 0 1];

    R = [cSA, -sSA, zero;
         sSA, cSA, zero;
         zero, zero, one];

    rVectECEF = reshape(rVectECEF,3,1,numElems);
%     rVectECI = R * rVectECEF;
    rVectECI = mtimesx(R,rVectECEF);
    rVectECI = reshape(rVectECI,3,numElems);

    rotRateRadSec = 2*pi/bodyInfo.rotperiod;
    omegaRI = repmat([0;0;rotRateRadSec],1,1,numElems);
    vVectECEF = reshape(vVectECEF,3,1,numElems);
%     vVectECI = R*(vVectECEF + cross(omegaRI, rVectECEF));
    vVectECI = mtimesx(R,(vVectECEF + cross(omegaRI, rVectECEF)));
    vVectECI = reshape(vVectECI,3,numElems);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [y,c]=hp(z,lambda);
% PURPOSE: Hodrick-Prescott filtering
% ------------------------------------------------------------
% SYNTAX: [y,c]=hp(z,lambda);
% ------------------------------------------------------------
% OUTPUT: y: nx1 --> long-term trend
%         c: nx1 --> cycle as deviation from long-term trend (as %)
% ------------------------------------------------------------
% INPUT: z: nx1      --> short-term trend or s.a. time series
%        lambda: 1x1 --> balance between adjustment and smoothness
% ------------------------------------------------------------
% LIBRARY: dif
% ------------------------------------------------------------

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>

% Version 2.0 [May 2009]

% Size of the input data
[n,m] = size(z);

% Second difference matrix
D = dif(2,n);
D(1:2,:) = [];

% Matrix form of the filter (non inverted)
H = eye(n) + lambda*(D'*D);

% Computing long-term trend
y = H \ z;

% Computing cycle as deviation from long term trend
c = 100 * (z - y) ./ y;

author: Sable

repo: mcbench-benchmarks


text:
function rbfinfwhiteKernDisplay(kern, spacing)

% RBFINFWHITEKERNDISPLAY Display parameters of the RBF-WHITE kernel (with
% integration limits between minus infinity and infinity).
% FORMAT
% DESC displays the parameters of the RBF-WHITE kernel and the kernel type
% to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : rbfinfwhiteKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('RBF-INF-WHITE inverse width: %2.4f (length scale %2.4f)\n', ...
        kern.inverseWidth, 1/sqrt(kern.inverseWidth));
fprintf(spacing);
fprintf('RBF-INF-WHITE variance: %2.4f\n', kern.variance)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function graph_arc_print ( nedge, inode, jnode, title )

%*****************************************************************************80
%
%% GRAPH_ARC_PRINT prints out a graph from an edge list.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEDGE, the number of edges.
%
%    Input, integer INODE(NEDGE), JNODE(NEDGE), the beginning
%    and end nodes of the edges.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for i = 1 : nedge
    fprintf ( 1, '  %2d  %4d  %4d\n', i, inode(i), jnode(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function loglik = computeEMLoglik(xInit,PInit,xSmooth,y,F,H,Q,R,B,u)
% COMPUTEEMLOGLIK Compute EM loglikelihood p(y_{1:t},x_{1:t} | theta), where
% theta = {F,H,Q,R,B}.
%
% INPUTS:   xInit       The (xDim x 1) initial state mean for time 1
%           PInit       The (xDim x xDim) initial state covariance for time 1
%           xSmooth     The (xDim x N) matrix of smoothed estimates for time 1:N
%           y           The (yDim x N) matrix of measurements for time 1:N
%           F           The (xDim x xDim) state transition matrix
%           H           The (yDim x xDim) measurement matrix
%           Q           The (xDim x xDim) process noise covariance
%           R           The (yDim x yDim) measurement noise covariance
%           B           The (xDim x uDim) control input gain matrix
%           u           The (uDim x N) control input matrix for time 1:N
%
% OUTPUTS:  loglik   A boolean stating if EM has converged
%
% [1] A. W. Blocker, An EM algorithm for the estimation of affine state-space systems with or without known inputs, (2008)
%
% January 2018 Lyudmil Vladimirov, University of Liverpool.

    [xDim, N] = size(xSmooth);
    yDim = size(y,1);

    loglik = -sum((y-H*xSmooth).^2)/(2*R) - N/2*log(abs(R))...
        -1/2*sum((xSmooth(2:N)-F*xSmooth(1:N-1)-B*u(2:N)).^2)/(Q) - (N-1)/2*log(abs(Q))...
        -sum((xSmooth(1)-xInit).^2)/(2*PInit) - 0.5*log(abs(PInit)) - N*log(2*pi);

%     loglik = -.5*(y-H*xSmooth)/R*(y-H*xSmooth)' - N/2*log(abs(R))...
%         -1/2*(xSmooth(:,2:N)-F*xSmooth(1:N-1)-B*u(2:N))/Q - (N-1)/2*log(abs(Q))...
%         -sum((xSmooth(1)-xInit).^2)/(2*PInit) - 0.5*log(abs(PInit)) - N*log(2*pi);

%     loglik = cell(1,3);
%     loglik{1} = -.5*(xSmooth(:,1) - xInit)/PInit*(xSmooth(:,1) - xInit)' ...
%                   -.5*log(abs(PInit)) - .5*(N-1)*log(abs(Q)) - .5*N*log(abs(R)) ...
%                   -.5*N*(xDim+yDim)*log(2*pi);
%     loglik{2} = 0;
%     loglik{3} = 0;
%     for k = 1:N
%         if(k>1)
%             loglik{2} = loglik{2} + .5*(xSmooth(:,k) ...
%                           - F*xSmooth(:,k-1) - B*u(:,k))/Q*(xSmooth(:,k) - F*xSmooth(:,k-1) - B*u(:,k))';
%         end
%         loglik{3} = loglik{3} + .5*(y(:,k) - H*xSmooth(:,k))/R*(y(:,k) - H*xSmooth(:,k))';
%     end
%     loglik = loglik{1} - loglik{2} - loglik{3};
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function invA = cholinv(A)

% SYNTAX:
%   invA = cholinv(A);
%
% INPUT:
%   A = positive definite matrix to be inverted
%
% OUTPUT:
%   invA = inverse of matrix A
%
% DESCRIPTION:
%   Inverse of a positive definite matrix by using Cholesky decomposition.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Nardo
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

% compute cholesky decomposition
n = size(A,1);
U    = chol(A);
invU = U\speye(n);
%L    = inv(L);
invA = invU * invU';

author: goGPS-Project


text:
function outputDir = GetOutputDataDir()

file_masterdir = GetMasterFileDir;
outputDir = fullfile(file_masterdir,'\!Proj FishExplorer\output');

end

author: xiuyechen

repo: FishExplorer


text:
function derivativeVector = helmholtzHelper(inputVector, q, alpha)
%Generates the helmholtz derivative vector for forward/backward convolving.
%
%helmholtzHelper takes a row or column input vector and a float, q, that is
%the order of the derivative to be taken, and generates a shifted operator
%operator that is modified for the Helmholtz equation. The vector returned
%is the equivalent of [I - A] where A is the operator.

    %Generates the derivative operator vector.
    [step, lowerBound, upperBound] = defaultBoundaries(numel(inputVector));
    derivativeVector = -additiveConvolutionVector(step, lowerBound, ...
                                                 upperBound, q);

    %Undoes a derivative.
    derivativeVector = helmholtzModifier(derivativeVector, alpha);
    derivativeVector = shiftOperator(derivativeVector);

end


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function n_particles = unique_particles(binary_particles)
% UNIQUE_PARTICLES: Function determines the number of unique particles
% in the current system.

% Determine the number of total particles
[n_models, ~] = size(binary_particles);

% Find unique particles
unique_particles = unique(binary_particles,'rows');

% Find n_particles
n_particles = size(unique_particles,1)/n_models;

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function r = isfinite(a)
%ISFINITE     Array of 1's for finite components
%
%   r = isfinite(a)
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  r = reshape( isfinite(a.x(:)') & all(isfinite(a.dx),1) & all(isfinite(a.hx),1) , size(a.x) );

repo: OpenMAS


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MATLAB Code for                                                  %
%                                                                   %
%  Multi-Objective Particle Swarm Optimization (MOPSO)              %
%  Version 1.0 - Feb. 2011                                          %
%                                                                   %
%  According to:                                                    %
%  Carlos A. Coello Coello et al.,                                  %
%  "Handling Multiple Objectives with Particle Swarm Optimization," %
%  IEEE Transactions on Evolutionary Computation, Vol. 8, No. 3,    %
%  pp. 256-279, June 2004.                                          %
%                                                                   %
%  Developed Using MATLAB R2009b (Version 7.9)                      %
%                                                                   %
%  Programmed By: S. Mostapha Kalami Heris                          %
%                                                                   %
%         e-Mail: sm.kalami@gmail.com                               %
%                 kalami@ee.kntu.ac.ir                              %
%                                                                   %
%       Homepage: http://www.kalami.ir                              %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function i=RouletteWheelSelection(p)

    r=rand;
    c=cumsum(p);
    i=find(r<=c,1,'first');

end

author: sfvsfv

repo: Mathematical-modeling


text:
function newdata=struct2jdata(data,varargin)
%
% newdata=struct2jdata(data,opt,...)
%
% convert a JData object (in the form of a struct array) into an array
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%      data: a struct array. If data contains JData keywords in the first
%            level children, these fields are parsed and regrouped into a
%            data object (arrays, trees, graphs etc) based on JData
%            specification. The JData keywords are
%               "_ArrayType_", "_ArraySize_", "_ArrayData_"
%               "_ArrayIsSparse_", "_ArrayIsComplex_"
%      opt: (optional) a list of 'Param',value pairs for additional options
%           The supported options include
%               'Recursive', if set to 1, will apply the conversion to
%                            every child; 0 to disable
%
% output:
%      newdata: the covnerted data if the input data does contain a JData
%               structure; otherwise, the same as the input.
%
% examples:
%      obj=struct('_ArrayType_','double','_ArraySize_',[2 3],
%                 '_ArrayIsSparse_',1 ,'_ArrayData_',null);
%      ubjdata=struct2jdata(obj);
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

fn=fieldnames(data);
newdata=data;
len=length(data);
if(jsonopt('Recursive',0,varargin{:})==1)
  for i=1:length(fn) % depth-first
    for j=1:len
        if(isstruct(getfield(data(j),fn{i})))
            newdata(j)=setfield(newdata(j),fn{i},jstruct2array(getfield(data(j),fn{i})));
        end
    end
  end
end
if(~isempty(strmatch('x0x5F_ArrayType_',fn)) && ~isempty(strmatch('x0x5F_ArrayData_',fn)))
  newdata=cell(len,1);
  for j=1:len
    ndata=cast(data(j).x0x5F_ArrayData_,data(j).x0x5F_ArrayType_);
    iscpx=0;
    if(~isempty(strmatch('x0x5F_ArrayIsComplex_',fn)))
        if(data(j).x0x5F_ArrayIsComplex_)
           iscpx=1;
        end
    end
    if(~isempty(strmatch('x0x5F_ArrayIsSparse_',fn)))
        if(data(j).x0x5F_ArrayIsSparse_)
            if(~isempty(strmatch('x0x5F_ArraySize_',fn)))
                dim=double(data(j).x0x5F_ArraySize_);
                if(iscpx && size(ndata,2)==4-any(dim==1))
                    ndata(:,end-1)=complex(ndata(:,end-1),ndata(:,end));
                end
                if isempty(ndata)
                    % All-zeros sparse
                    ndata=sparse(dim(1),prod(dim(2:end)));
                elseif dim(1)==1
                    % Sparse row vector
                    ndata=sparse(1,ndata(:,1),ndata(:,2),dim(1),prod(dim(2:end)));
                elseif dim(2)==1
                    % Sparse column vector
                    ndata=sparse(ndata(:,1),1,ndata(:,2),dim(1),prod(dim(2:end)));
                else
                    % Generic sparse array.
                    ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3),dim(1),prod(dim(2:end)));
                end
            else
                if(iscpx && size(ndata,2)==4)
                    ndata(:,3)=complex(ndata(:,3),ndata(:,4));
                end
                ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3));
            end
        end
    elseif(~isempty(strmatch('x0x5F_ArraySize_',fn)))
        if(iscpx && size(ndata,2)==2)
             ndata=complex(ndata(:,1),ndata(:,2));
        end
        ndata=reshape(ndata(:),data(j).x0x5F_ArraySize_);
    end
    newdata{j}=ndata;
  end
  if(len==1)
      newdata=newdata{1};
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x,t]=meanLinTrajConstEnds(x0,x2,deltaT,q0,order,N,FFun,QFun)
%%MEANLINTRAJCONSTENDS Generate a the mean trajectory based on a linear
%               dynamic model with additive noise, as is typically used in
%               tracking systems, given known locations of the start and
%               end points. A polynomial dynamic model is used by default,
%               but one can specify other models.
%
%INPUTS: x0, x2 The xDimX1 target states at the beginning and end of the
%               region over which paths are to be simulated. These are
%               treated as known values.
%        deltaT The time between when state x0 occurred and when state x2
%               occurred. This is a positive value.
%            q0 The power spectral density of the process noise. This is
%               passed to QFun, as described below and by default would
%               corresponds to the same-named input in QPolyKal. With the
%               default polynomial dynamic model, FPolyKal and QPolyKal, q0
%               has no effect on the result and is just set to 1 if an
%               empty matrix is passed.
%         order The order >=0 of the filter. If order=1, then it is
%               constant velocity, 2 means constant acceleration, 3 means
%               constant jerk, etc. This is passed to FFun and QFun. If
%               custom FFun and QFun functions are used and this is not
%               needed, then an empty matrix can be passed. The default if
%               omitted or an empty matrix is passed is 1.
%             N The number of points in the random path to generate between
%               x0 and x2. The default if omitted or an empty matrix is
%               passed is 100; N>=0.
%    FFun, QFun Optional function handles to functions that respectively
%               provide the state transition matrix and the process noise
%               covariance matrix as a function of the time duration over
%               which the prediction is taken. The defaults if omitted or
%               empty matrices are passed are handles to FPolyKal and
%               QPolyKal. The functions are called as FFun(T,xDim,order)
%               and QFun(T,xDim,order,q0), where T is the duration over
%               which the prediction is taken.
%
%OUTPUTS: x The xDimX(N+2) set of mean trajectory according to the
%           specified model. x(:,1) is always x0 and x(:,N+2) is always x2.
%         t The (N+2)X1 set of time offsets (starting at 0) where the
%           values in x are generated.
%
%Consider a process at points 0, 1, and 2, where we are given the value at
%time x0. The transition models are:
%x1=F0*x0+w1
%x2=F1*x1+w2
%Where x0 is given (deterministic) and w1 and w2 are independent zero-mean
%Gaussian random variables with covariance matrices Q1 and Q1.
%x1 and x2 are jointly normal with mean mu and covariance matrix Sigma:
%mu=[x1Pred; = [F0*x0
%    x2Pred]    F1*F0*x0];
%Sigma=[   Q1,        Q1*F1'; =[Pxx, Pxz;
%       F1*Q1,  F1*Q1*F1'+Q2]   Pxz',Pzz];
%When given the value of x2, the conditional distribution of x1|x0,x2 is
%the same as the Kalman filter update taking the measurement covariance
%matrix to be Q2, the measurement matrix H to be F1, the prior state to be
%F0*x0, and the covariance of the prior state to be Q1. Thus, as derived,
%for example, in [Ch. 5,1], [2], the posterior distribution of x1 is
%Gaussian with mean and covariance matrix
%x1Bar=x1Pred+Pxz/Pzz*(x2-x2Pred)
%This function is similar to randLinTrajConstEnds.
%
%EXAMPLE:
%This plots a trajectory that eventually has to turn around, since the
%endpoint is not close to the direction of the initial velocity vector.
% x0=[1e3;0;30;20];
% x2=[709;
%     527;
%     -39;
%       4];
% deltaT=100;
% x=meanLinTrajConstEnds(x0,x2,deltaT);
% figure(1)
% clf
% plot(x(1,:),x(2,:))
% h1=xlabel('x');
% h2=ylabel('y');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(N))
    N=100;
end

if(nargin<5||isempty(order))
    order=1;
end

if(nargin<4||isempty(q0))
    q0=1;
end

xDim=size(x0,1);

if(nargin<8||isempty(QFun))
    Q=@(T)QPolyKal(T,xDim,order,q0);
else
    Q=@(T)QFun(T,xDim,order,q);
end

if(nargin<7||isempty(FFun))
    F=@(T)FPolyKal(T,xDim,order);
else
    F=@(T)FFun(T,xDim,order);
end

x=zeros(xDim,N+2);

%The increment between the sample points.
deltaTInc=deltaT/(N+1);

x(:,1,:)=x0;
x(:,N+2,:)=x2;

F0=F(deltaTInc);
Q1=Q(deltaTInc);
for curPoint=1:N
    deltaT0=deltaTInc*(curPoint-1);
    deltaT2=deltaT-deltaT0-deltaTInc;

    F1=F(deltaT2);
    Q2=Q(deltaT2);

    Pxz=Q1*F1';
    Pzz=F1*Q1*F1'+Q2;
    x1Pred=F0*x(:,curPoint);
    x2Pred=F1*F0*x(:,curPoint);
    x1Mean=x1Pred+Pxz*inv(Pzz)*(x2-x2Pred);

    x(:,curPoint+1)=x1Mean;
end

if(nargout>1)
    t=[0;deltaTInc*(1:(N+1)).'];
    t(end)=deltaT;%Make it exact despite finite precision limitations.
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the CHANNEL_CHANNEL-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Channel_Channel()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  512;       % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  512;       % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
ds= min(Lx/(2*Nx),Ly/(2*Ny));  % Lagrangian spacing
L = 0.9*Lx;                    % Length of Channel
w = 0.2*Ly;                    % Width of Channel
x0 = 0.3;                      % x-Center for Cylinder
y0 = 0.5;                      % y-Center for Cylinder
r = w/10;                      % Radii of Cylinder
struct_name = 'channel'; % Name for .vertex, .spring, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Channel_Immsersed_Boundary_Geometry(ds,L,w,Lx,Ly);
[xLag_C,yLag_C] = give_Me_Cylinder_Immsersed_Boundary_Geometry(ds,r,x0,y0);
[xLag_V,yLag_V] = give_Me_Piston_IB_Geometry(ds,w,Lx,Ly);

Npts = length(xLag_C);
NPistonStart = length(xLag);
fprintf('\n\n          POSSIBLE WARNING!!!      \n\n');
fprintf('# of pts in a cell: %d\n',Npts);
fprintf('If it is not EVEN...there will be a problem :)\n');
fprintf('Change radius, r, or ds slightly to try to make it even\n\n\n');

% Compute Curvatures
C =  compute_Curvatures(xLag_C,yLag_C);

% Plot Geometry to test
plot(xLag(1:end/2),yLag(1:end/2),'r-'); hold on;
plot(xLag(end/2+1:end),yLag(end/2+1:end),'r-'); hold on;
plot(xLag_C,yLag_C,'r-'); hold on;
plot(xLag_V,yLag_V,'k.-'); hold on;

plot(xLag,yLag,'*'); hold on;
plot(xLag_C,yLag_C,'g*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Ly]);

% Create different Bubbles
xLag_C1 = xLag_C; yLag_C1 = yLag_C;              % Bubble 1
xLag_C2 = xLag_C+33.5*r; yLag_C2 = yLag_C+2*r;      % Bubble 2
xLag_C3 = xLag_C+1.5*r; yLag_C3 = yLag_C-2.5*r;  % Bubble 3

% Combine Geometry Files into ONE Vector
%xLag = [xLag_C1 xLag_C2 xLag_C3 xLag xLag_V];
%yLag = [yLag_C1 yLag_C2 yLag_C3 yLag yLag_V];
xLag = [xLag_C1 xLag_C2 xLag_C3 xLag_V];
yLag = [yLag_C1 yLag_C2 yLag_C3 yLag_V];
yLag = yLag - 0.375;

first_Index = 3*length(xLag_C1)+1;
Ncell_Lag_Pts = [length(xLag_C1) length(xLag_C2) length(xLag_C3)];

% Combine CELL Geometry files into ONE vector for springs/beams (not necessary...just easier to count)
xLag_Cells = [xLag_C1 xLag_C2 xLag_C3];
yLag_Cells = [yLag_C1 yLag_C2 yLag_C3];

plot(xLag_Cells,yLag_Cells,'k*'); hold on;

%
% Prints .vertex file!
%
print_Lagrangian_Vertices(xLag,yLag,struct_name);

%
% Prints .spring file!
%
k_Spring = 2.5e7;
NPistonStart = first_Index; % Gives starting index of piston instead of num of channel pts
NPiston = length(xLag_V);
print_Lagrangian_Springs(xLag_C,yLag_C,xLag,yLag,k_Spring,ds,r,struct_name,NPistonStart,NPiston);
plot(xLag(NPistonStart),yLag(NPistonStart),'r*'); hold on;
plot(xLag(NPistonStart+NPiston-1),yLag(NPistonStart+NPiston-1),'g*'); hold on

%
% Prints .beam file!
%
k_Beam = 1e12;
print_Lagrangian_Beams(xLag_C,yLag_C,k_Beam,C,struct_name);

%
% Prints .coagulation file!
%
starting_index = 1;         % index of first cellular Lagrangian Pt.
threshold_radius = 0.075;   % how close do cells need to be to form a bond
bond_strength = 1e3;        % bond strength
fracture_force = 1e2;       % how strong of force to break bond
print_Lagrangian_Coagulation_Inputs(struct_name,Ncell_Lag_Pts,starting_index,threshold_radius,fracture_force,bond_strength);

%
% Prints .target file!
%
k_Target = 1e7;
fprintf('For update_target_pts: last target pt index before piston: %d\n\n\n',0);
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name,first_Index);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called struct.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints COAGULATION Parameters to the file struct.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Coagulation_Inputs(struct_name,Ncell_Lag_Pts,starting_index,threshold_radius,fracture_force,bond_strength)

    % Ncell_Lag_Pts:    -gives # of Lag. Pts in each cell (# of Lag. Pts that compose each cell)
    % starting_index:   -gives first index of first cell Lag. Pt. in vertex file
    %                   -NOTE: IB2d Assumes all cell Lag. Pts are next to each other in .vertex file
    % threshold_radius: -how close do cells need to be to form a bond
    % fracture_force:   -how strong of force to break bond between cells

    Ncells = length(Ncell_Lag_Pts);   % gives total number of cells

    coag_fid = fopen([struct_name '.coagulation'], 'w');

    fprintf(coag_fid, '%d\n', starting_index );   % Prints first index of a Lag. Pt. associated to a Cell
    fprintf(coag_fid, '%d\n', threshold_radius ); % Prints threshold radii
    fprintf(coag_fid, '%d\n', bond_strength );    % Prints bond strength
    fprintf(coag_fid, '%d\n', fracture_force );   % Prints fracture_force

    %Loops over all CELLs.
    for s = 1:Ncells
        fprintf(coag_fid, '%1.16e\n', Ncell_Lag_Pts(s) );
    end

    fclose(coag_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Target points to a file called struct.target
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name,first_index)

    N = length(xLag) - first_index + 1;  % # of target points (e.g., the walls of channel)

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 0:N-1                   % Loops over N target pts...just starts at zero bc of how first_index is defined
        fprintf(target_fid, '%d %1.16e\n', s+first_index, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called struct.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,xLagAll,yLagAll,k_Spring,ds,r,struct_name,NPistonStart,NPiston)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', 3*N + 3/2*N + (NPiston-1) ); %N MUST BE EVEN!

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:3*N
            if s < N
                x1 = xLag(s);    y1 = yLag(s);
                x2 = xLag(s+1);  y2 = yLag(s+1);
                ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif s==N
                %Case s=N
                x1 = xLag(s);    y1 = yLag(s);
                x2 = xLag(1);    y2 = yLag(1);
                ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            elseif ( ( s>N ) && ( s<2*N ) )
                x1 = xLagAll(s);    y1 = yLagAll(s);
                x2 = xLagAll(s+1);  y2 = yLagAll(s+1);
                ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif s==2*N
                x1 = xLagAll(s);    y1 = yLagAll(s);
                x2 = xLagAll(N+1);  y2 = yLagAll(N+1);
                ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, N+1,   k_Spring, ds_Rest);
            elseif ( ( s>2*N ) && ( s<3*N ) )
                x1 = xLagAll(s);    y1 = yLagAll(s);
                x2 = xLagAll(s+1);  y2 = yLagAll(s+1);
                ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif s==3*N
                x1 = xLagAll(s);      y1 = yLagAll(s);
                x2 = xLagAll(2*N+1);  y2 = yLagAll(2*N+1);
                ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 2*N+1,   k_Spring, ds_Rest);
            end
    end

    for s=1:N/2
        x1 = xLagAll(s);      y1 = yLagAll(s);
        x2 = xLagAll(s+N/2);  y2 = yLagAll(s+N/2);
        ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+N/2, k_Spring/5000, ds_Rest);
    end

    for s=1:N/2
        x1 = xLagAll(N+s);      y1 = yLagAll(N+s);
        x2 = xLagAll(N+s+N/2);  y2 = yLagAll(N+s+N/2);
        ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', N+s, N + s+N/2, k_Spring/5000, 2*r);
    end

    for s=1:N/2
        x1 = xLagAll(2*N+s);      y1 = yLagAll(2*N+s);
        x2 = xLagAll(2*N+s+N/2);  y2 = yLagAll(2*N+s+N/2);
        ds_Rest = sqrt( (x1-x2)^2 + (y1-y2)^2 );
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', 2*N+s, 2*N + s+N/2, k_Spring/5000, 2*r);
    end

    for s=NPistonStart:NPistonStart+NPiston-2
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called struct.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', 3*N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 1:N
            if s==1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',N, s, s+1,   k_Beam, C(s) );
            elseif  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s) );
            elseif s==N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C(s) );
            end
    end

    for s = N+1:2*N
            if s==N+1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',2*N, s, s+1,   k_Beam, C(s-N) );
            elseif  s <= 2*N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s-N) );
            elseif s==2*N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, N+1,   k_Beam, C(s-N) );
            end
    end

    for s = 2*N+1:3*N
            if s==2*N+1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',3*N, s, s+1,   k_Beam, C(s-2*N) );
            elseif  s <= 3*N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s-2*N) );
            elseif s==3*N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 2*N+1,   k_Beam, C(s-2*N) );
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes "curvature" of ellipse
%
% NOTE: not curvature in the traditional geometric sense, in the 'discrete'
% sense through cross product.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = compute_Curvatures(xLag,yLag)

%a-x component (rmin)
%b-y component (rmax)
%C = ab / ( sqrt( a^2*sin(t)^2 + b^2*cos(t)^2  )  )^3

N = length(xLag);
C = zeros( N );

%Note: needs to be done same order as you print .beam file!
for i=1:N

   % Pts Xp -> Xq -> Xr (same as beam force calc.)

   if ( (i > 1) && (i < N) )

        Xp = xLag(i-1); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(i-1); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==1)

        Xp = xLag(N); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(N); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==N)

        Xp = xLag(N-1); Xq = xLag(N); Xr = xLag(1);
        Yp = yLag(N-1); Yq = yLag(N); Yr = yLag(1);

   end

   C(i) = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp); %Cross product btwn vectors

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for channel
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Channel_Immsersed_Boundary_Geometry(ds,L,w,Lx,Ly)

% The immsersed structure is a channel %
x = (Lx-L)/2:ds:(L+(Lx-L)/2);  %xPts
yBot = (Ly-w)/2;               %yVal for bottom of Channel
yTop = Ly - (Ly-w)/2;          %yVal for top of Channel

xLag = [x x];
yLag = [yBot*ones(1,length(x)) yTop*ones(1,length(x))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives Lag. Pts for Piston
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Piston_IB_Geometry(ds,w,Lx,Ly)

yBot = (Ly-w)/2;               %yVal for bottom of Channel
yTop = Ly - (Ly-w)/2;          %yVal for top of Channel

yLag = Ly/2-0.05:ds:Ly/2+0.05;
xLag = 0.25*ones(1,length(yLag));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for cylinder
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Cylinder_Immsersed_Boundary_Geometry(ds,r,x0,y0)

% The immsersed structure is a cylinder %

dtheta = ds/ (2.08*r);
theta = 0; i=1;
while theta < 2*pi
   xLag(i) = x0 - r*cos(theta);
   yLag(i) = y0 - r*sin(theta);
   theta = theta + dtheta;
   i=i+1;
end

%ds_Rest = sqrt( ( xLag(2) - xLag(3) )^2 + ( xLag(2) - xLag(3) )^2 );

author: nickabattista


text:
% GP_COV_PP - Piecewise polynomial covariance function with compact support
% for Gaussian processes.
%
% [K, DK_LOGTHETA, DK_X2] = GP_COV_CS(X1, X2, LOGTHETA)

% Last modified 2010-10-28
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function covfunc = gp_cov_pp_proto(x1,x2)

if nargin < 2
  x2 = x1;
end

% Dimensions
d = rows(x2); % dimensionality of inputs
m = cols(x1); % number of other inputs
n = cols(x2); % number of inputs

q = 2; % PP parameter, fix to 2 for now..

% Distance matrix
D = sqrt(sq_dist(x1,x2));

covfunc = @get_covariance;

  function [K,L,M] = get_covariance(theta)

  if nargin == 0
    % Return the number of parameters and the size of the covariance matrix
    K = 1;
    L = m;
    M = n;
    return
  end

  % Threshold
  thres = theta(1); % length scale, threshold
  R = D./thres;
  I = sparse(R<1);
  NNZ = nnz(I); % number of nonzero elements

  j = floor(rows(x1)) + q + 1;
  switch q
    %%%%%%%%%%
   case 0
    error('not yet implemented for q=0')
    K(I) = (1-R(I)) .^ j;
    %%%%%%%%%%
   case 1
    error('not yet implemented for q=1')
    %%%%%%%%%%

   case 2
    a = j^2 + 4*j + 3;
    b = 3*j + 6;
    c = 3;
    z = 1/3;

    % Helper
    R2 = spalloc(m,n,NNZ);
    R2(I) = D(I).^2 ./ (thres^2);

    % Covariance matrix
    K = spalloc(m,n,NNZ);
    K(I) = z * (1-R(I)).^(j+2) .* (a*R2(I) + b*R(I) + c);

    dK_dR = spalloc(m,n,NNZ); % copy sparseness
    dK_dR(I) = -z * (j+2) * (1-R(I)).^(j+1) .* (a*R2(I) + b*R(I) + c) ...
        + z * (1-R(I)).^(j+2) .* (2*a*R(I) + b);

    % Gradient for hyperparameters
    if nargout >= 2
      % TODO 3D SPARSE MATRICES DO NOT WORK. USE CELL ARRAYS?
      dK_dlogtheta = spalloc(m,n,NNZ);
      dK_dlogtheta(I) = dK_dR(I) .* (-D(I) ./ thres);
      dK_dlogtheta = full(dK_dlogtheta); % blaaah.. :(
    end

    % Gradients for inputs x2
    if nargout >= 3
      if isempty(x2)
        error('Can''t calculate gradient: x2 not given');
      end
      % TODO 3D SPARSE MATRICES DO NOT WORK. USE CELL ARRAYS?
      % TODO: you have dK_dR and dD_dx but NOT dR_dD!!
      dK_dx2 = bsxfun(@times, reshape(full(dK_dR./thres),[1,m,n]), dD_dx2);
      % blaah the need for fullness.. :(
    end

    %%%%%%%%%%

   case 3
    error('not yet implemented for q=3')
    %%%%%%%%%%
   otherwise
    error('q not valid')
  end

  if ~issparse(K)
    error('K not sparse, it should, wtf?!');
  end

  if nargout >= 2
    L = ldlchol2chol( ldlchol(K) );
  end
  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [filename,xform] = ctrSave(ctr,filename,xform)
% Save Contrack (ctr) structure to file.
%
%  [filename,xform] = ctrSave(ctr,filename)
%
% Saves the ConTrack parameters (in the ctr structure) to a file in ASCII
% format. This file is read by ConTrack on execution.  This text file is
% normally saved relative to the dt6 directory in the parallel directory
% of fibers, ..\fibers\contrack\XXX.
%
% contrack itself is called from the dt6 file directory (p.dt6Dir)
%
% See also ctrCreate, ctrGet, ctrSet
%

% Check input parameters
if notDefined('ctr'), error(' ConTrack structure required.'); end
if (notDefined('filename'))
    filename = mrvSelectFile('w','txt',[],'ConTrack parameters filename');
    if isempty(filename), disp('User canceled'); return; end
end
if notDefined('xform'), xform = eye(4); end

% Open up the file, checking that the user sent in a good file path
fileDir = fileparts(filename);
if isempty(fileDir), fileDir=pwd; end

if ~isdir(fileDir), error('No directory %s',fileDir);
else fid = fopen(filename,'wt');
    if(fid == -1)
        error('Can not write ConTrack parameters file.');
    end
end

% Start the writing
fprintf(fid,'Params: %.0f\n',ctrGet(ctr,'version'));  % Should this be a version parameter
fprintf(fid,'Image Directory: %s/ \n',ctrGet(ctr,'image_directory'));
fprintf(fid,'WM/GM Mask Filename: %s\n',ctrGet(ctr,'fa_filename'));
fprintf(fid,'PDF Filename: %s\n',ctrGet(ctr,'pdf_filename'));
fprintf(fid,'ROI MASK Filename: %s\n',ctrGet(ctr,'mask_filename'));
fprintf(fid,'Desired Samples: %g\n',ctrGet(ctr,'desired_samples'));
fprintf(fid,'Max Pathway Nodes: %g\n',ctrGet(ctr,'max_nodes'));
fprintf(fid,'Min Pathway Nodes: %g\n',ctrGet(ctr,'min_nodes'));
fprintf(fid,'Step Size (mm): %g\n',ctrGet(ctr,'step_size'));
fprintf(fid,'Start Is Seed VOI: %s\n',ctrGet(ctr,'roi',1,'seed_region'));
fprintf(fid,'End Is Seed VOI: %s\n',ctrGet(ctr,'roi',2,'seed_region'));
fprintf(fid,'Save Out Spacing: %g\n',ctrGet(ctr,'save_out_spacing'));
fprintf(fid,'Threshold for WM/GM specification: %g\n',ctrGet(ctr,'fa_absorb'));
fprintf(fid,'Absorption Rate WM: %g\n',ctrGet(ctr,'abs_normal'));
fprintf(fid,'Absorption Rate NotWM: %g\n',ctrGet(ctr,'abs_penalty'));
fprintf(fid,'Local Path Segment Smoothness Standard Deviation: %g\n',ctrGet(ctr,'smooth_std'));
fprintf(fid,'Local Path Segment Angle Cutoff: %g\n',ctrGet(ctr,'angle_cutoff'));

% Can't we use num2str here?  Is it an issue for Tony's read?
fprintf(fid,'ShapeFunc Params (LinMidCl,LinWidthCl,UniformS): [ %s ]\n',getStrFromVec(ctrGet(ctr,'shape_params_vec')));

% Close up
fclose(fid);

return;

%------------------------------------------------
function str = getStrFromVec(vec)
% Shouldn't this just be num2str?
str = [];
for ii = 1:length(vec)
    str = [str num2str(vec(ii)) ', ']; %#ok<AGROW>
end
str = str(1:end-2);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function writeWIRE(filename, V,E)
  % WRITEWIRE writes an .wire file with vertex/edge information (see PyMesh
  % documentation).
  %
  % writeWIRE(filename,V,F,UV,N)
  %
  % Input:
  %  filename  path to .wire file
  %  V  #V by 3 list of vertices
  %  E  #E by 2 list of edge indices into E
  %

  f = fopen( filename, 'w' );
  if size(V,2) == 2
    warning('Appending 0s as z-coordinate');
    V(:,end+1:3) = 0;
  else
    assert(size(V,2) == 3);
  end

  fprintf( f, 'v %0.17g %0.17g %0.17g\n', V');
  assert(size(E,2) == 2,'Edge list should have two columns');
  fprintf( f, 'l %d %d\n', E');

  fclose(f);
end

author: alecjacobson

repo: gptoolbox


text:
function [l,m] =  findnei(k)
    % This script finds overlapping alarms in space-time
    % and groups them togetehr
    %
    % Stefan Wiemer    4/95
    global abo iala

    report_this_filefun(mfilename('fullpath'));

    d = sqrt(((abo(k,1) - abo(:,1))).^2 + ((abo(k,2) - abo(:,2))).^2);
    m = d < abo(:,3)+abo(k,3) &  abs(abo(:,5)-abo(k,5)) < iala;
    l = find(m == 1);

author: CelsoReyes


text:
function header = dtiReadPDBHeader (fid)

if (~exist('fid','var'))
    header = struct();
    header.xformToAcPc = repmat (0, 4);
    header.numPaths = 0;
else
    offset = fread (fid, 1, 'uint');
    mx = fread (fid, 16, 'double');
    header.xformToAcPc = reshape (mx, 4,4);
    fseek (fid, offset-4, -1);
    header.numPaths = fread (fid, 1, 'uint');
end

% xxx should contain information about statistics too!


text:
function [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
  = lsqr( m, n, A, b, damp, atol, btol, conlim, itnlim, show )
%
%        [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
% = lsqr( m, n, A, b, damp, atol, btol, conlim, itnlim, show );
%
% LSQR solves  Ax = b  or  min ||b - Ax||_2  if damp = 0,
% or   min || (b)  -  (  A   )x ||   otherwise.
%          || (0)     (damp I)  ||2
% A  is an m by n matrix defined or a function handle of aprod( mode,x ),
% that performs the matrix-vector operations.
% If mode = 1,   aprod  must return  y = Ax   without altering x.
% If mode = 2,   aprod  must return  y = A'x  without altering x.

%-----------------------------------------------------------------------
% LSQR uses an iterative (conjugate-gradient-like) method.
% For further information, see
% 1. C. C. Paige and M. A. Saunders (1982a).
%    LSQR: An algorithm for sparse linear equations and sparse least squares,
%    ACM TOMS 8(1), 43-71.
% 2. C. C. Paige and M. A. Saunders (1982b).
%    Algorithm 583.  LSQR: Sparse linear equations and least squares problems,
%    ACM TOMS 8(2), 195-209.
% 3. M. A. Saunders (1995).  Solution of sparse rectangular systems using
%    LSQR and CRAIG, BIT 35, 588-604.
%
% Input parameters:
% atol, btol  are stopping tolerances.  If both are 1.0e-9 (say),
%             the final residual norm should be accurate to about 9 digits.
%             (The final x will usually have fewer correct digits,
%             depending on cond(A) and the size of damp.)
% conlim      is also a stopping tolerance.  lsqr terminates if an estimate
%             of cond(A) exceeds conlim.  For compatible systems Ax = b,
%             conlim could be as large as 1.0e+12 (say).  For least-squares
%             problems, conlim should be less than 1.0e+8.
%             Maximum precision can be obtained by setting
%             atol = btol = conlim = zero, but the number of iterations
%             may then be excessive.
% itnlim      is an explicit limit on iterations (for safety).
% show = 1    gives an iteration log,
% show = 0    suppresses output.
%
% Output parameters:
% x           is the final solution.
% istop       gives the reason for termination.
% istop       = 1 means x is an approximate solution to Ax = b.
%             = 2 means x approximately solves the least-squares problem.
% r1norm      = norm(r), where r = b - Ax.
% r2norm      = sqrt( norm(r)^2  +  damp^2 * norm(x)^2 )
%             = r1norm if damp = 0.
% anorm       = estimate of Frobenius norm of Abar = [  A   ].
%                                                    [damp*I]
% acond       = estimate of cond(Abar).
% arnorm      = estimate of norm(A'*r - damp^2*x).
% xnorm       = norm(x).
% var         (if present) estimates all diagonals of (A'A)^{-1} (if damp=0)
%             or more generally (A'A + damp^2*I)^{-1}.
%             This is well defined if A has full column rank or damp > 0.
%             (Not sure what var means if rank(A) < n and damp = 0.)
%
%
%        1990: Derived from Fortran 77 version of LSQR.
% 22 May 1992: bbnorm was used incorrectly.  Replaced by anorm.
% 26 Oct 1992: More input and output parameters added.
% 01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.
%              Print log reformatted.
% 14 Jun 1997: show  added to allow printing or not.
% 30 Jun 1997: var   added as an optional output parameter.
% 07 Aug 2002: Output parameter rnorm replaced by r1norm and r2norm.
%              Michael Saunders, Systems Optimization Laboratory,
%              Dept of MS&E, Stanford University.
% 03 Jul 2007: Modified 'aprodname' to A, which can either be an m by n
%              matrix, or a function handle.
%              Ewout van den Berg, University of British Columbia
% 03 Jul 2007: Modified 'test2' condition, omitted 'test1'.
%              Ewout van den Berg, University of British Columbia
%-----------------------------------------------------------------------

%     Initialize.

msg=['The exact solution is  x = 0                              '
     'Ax - b is small enough, given atol, btol                  '
     'The least-squares solution is good enough, given atol     '
     'The estimate of cond(Abar) has exceeded conlim            '
     'Ax - b is small enough for this machine                   '
     'The least-squares solution is good enough for this machine'
     'Cond(Abar) seems to be too large for this machine         '
     'The iteration limit has been reached                      '];

wantvar= nargout >= 6;
if wantvar, var = zeros(n,1); end

if show
   disp(' ')
   disp('LSQR            Least-squares solution of  Ax = b')
   str1 = sprintf('The matrix A has %8g rows  and %8g cols', m, n);
   str2 = sprintf('damp = %20.14e    wantvar = %8g', damp,wantvar);
   str3 = sprintf('atol = %8.2e                 conlim = %8.2e', atol, conlim);
   str4 = sprintf('btol = %8.2e                 itnlim = %8g'  , btol, itnlim);
   disp(str1);   disp(str2);   disp(str3);   disp(str4);
end

itn    = 0;		 istop  = 0;        nstop  = 0;
ctol   = 0;		 if conlim > 0, ctol = 1/conlim; end;
anorm  = 0;		 acond  = 0;
dampsq = damp^2; ddnorm = 0;        res2   = 0;
xnorm  = 0;	     xxnorm = 0;        z      = 0;
cs2    = -1;     sn2    = 0;

% Set up the first vectors u and v for the bidiagonalization.

% These satisfy  beta*u = b,  alfa*v = A'u.

u      = b(1:m);	x    = zeros(n,1);
alfa   = 0;		beta = norm( u );
if beta > 0
   u = (1/beta) * u;	v = Aprod(u,2);
   alfa = norm( v );
end
if alfa > 0
   v = (1/alfa) * v;    w = v;
end

arnorm = alfa * beta;
if arnorm == 0
   if show, disp(msg(1,:)); end
   return
end
arnorm0= arnorm;

rhobar = alfa;		phibar = beta;		bnorm  = beta;
rnorm  = beta;
r1norm = rnorm;
r2norm = rnorm;
head1  = '   Itn      x(1)       r1norm     r2norm ';
head2  = ' Compatible   LS      Norm A   Cond A';

if show
   disp(' ')
   disp([head1 head2])
   test1  = 1;		test2  = alfa / beta;
   str1   = sprintf( '%6g %12.5e',        itn,   x(1) );
   str2   = sprintf( ' %10.3e %10.3e', r1norm, r2norm );
   str3   = sprintf( '  %8.1e %8.1e',   test1,  test2 );
   disp([str1 str2 str3])
end

%------------------------------------------------------------------
%     Main iteration loop.
%------------------------------------------------------------------
while itn < itnlim
      itn = itn + 1;
%     Perform the next step of the bidiagonalization to obtain the
%     next  beta, u, alfa, v.  These satisfy the relations
%                beta*u  =  a*v   -  alfa*u,
%                alfa*v  =  A'*u  -  beta*v.

      u    = Aprod(v,1)  -  alfa*u;
      beta = norm( u );
      if beta > 0
         u     = (1/beta) * u;
         anorm = norm([anorm alfa beta damp]);
         v     = Aprod(u, 2)  -  beta*v;
         alfa  = norm( v );
         if alfa > 0,  v = (1/alfa) * v; end
      end

%     Use a plane rotation to eliminate the damping parameter.
%     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.

      rhobar1 = norm([rhobar damp]);
      cs1     = rhobar / rhobar1;
      sn1     = damp   / rhobar1;
      psi     = sn1 * phibar;
      phibar  = cs1 * phibar;

%     Use a plane rotation to eliminate the subdiagonal element (beta)
%     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.

      rho     = norm([rhobar1 beta]);
      cs      =   rhobar1/ rho;
      sn      =   beta   / rho;
      theta   =   sn * alfa;
      rhobar  = - cs * alfa;
      phi     =   cs * phibar;
      phibar  =   sn * phibar;
      tau     =   sn * phi;

%     Update x and w.

      t1      =   phi  /rho;
      t2      = - theta/rho;
      dk      =   (1/rho)*w;

      x       = x      +  t1*w;
      w       = v      +  t2*w;
      ddnorm  = ddnorm +  norm(dk)^2;
      if wantvar, var = var  +  dk.*dk; end

%     Use a plane rotation on the right to eliminate the
%     super-diagonal element (theta) of the upper-bidiagonal matrix.
%     Then use the result to estimate  norm(x).

      delta   =   sn2 * rho;
      gambar  = - cs2 * rho;
      rhs     =   phi  -  delta * z;
      zbar    =   rhs / gambar;
      xnorm   =   sqrt(xxnorm + zbar^2);
      gamma   =   norm([gambar theta]);
      cs2     =   gambar / gamma;
      sn2     =   theta  / gamma;
      z       =   rhs    / gamma;
      xxnorm  =   xxnorm  +  z^2;

%     Test for convergence.
%     First, estimate the condition of the matrix  Abar,
%     and the norms of  rbar  and  Abar'rbar.

      acond   =   anorm * sqrt( ddnorm );
      res1    =   phibar^2;
      res2    =   res2  +  psi^2;
      rnorm   =   sqrt( res1 + res2 );
      arnorm  =   alfa * abs( tau );

%     07 Aug 2002:
%     Distinguish between
%        r1norm = ||b - Ax|| and
%        r2norm = rnorm in current code
%               = sqrt(r1norm^2 + damp^2*||x||^2).
%        Estimate r1norm from
%        r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
%     Although there is cancellation, it might be accurate enough.

      r1sq    =   rnorm^2  -  dampsq * xxnorm;
      r1norm  =   sqrt( abs(r1sq) );   if r1sq < 0, r1norm = - r1norm; end
      r2norm  =   rnorm;

%     Now use these norms to estimate certain other quantities,
%     some of which will be small near a solution.

      test1   =   rnorm / bnorm;
      test2   =   arnorm / arnorm0;
%     test2   =   arnorm/( anorm * rnorm );
      test3   =       1 / acond;
      t1      =   test1 / (1    +  anorm * xnorm / bnorm);
      rtol    =   btol  +  atol *  anorm * xnorm / bnorm;

%     The following tests guard against extremely small values of
%     atol, btol  or  ctol.  (The user may have set any or all of
%     the parameters  atol, btol, conlim  to 0.)
%     The effect is equivalent to the normal tests using
%     atol = eps,  btol = eps,  conlim = 1/eps.

      if itn >= itnlim,   istop = 7; end
      if 1 + test3  <= 1, istop = 6; end
      if 1 + test2  <= 1, istop = 5; end
      if 1 + t1     <= 1, istop = 4; end

%     Allow for tolerances set by the user.

      if  test3 <= ctol,  istop = 3; end
      if  test2 <= atol,  istop = 2; end
%     if  test1 <= rtol,  istop = 1; end

%     See if it is time to print something.

      prnt = 0;
      if n     <= 40       , prnt = 1; end
      if itn   <= 10       , prnt = 1; end
      if itn   >= itnlim-10, prnt = 1; end
      if rem(itn,10) == 0  , prnt = 1; end
      if test3 <=  2*ctol  , prnt = 1; end
      if test2 <= 10*atol  , prnt = 1; end
%     if test1 <= 10*rtol  , prnt = 1; end
      if istop ~=  0       , prnt = 1; end

      if prnt == 1
         if show
            str1 = sprintf( '%6g %12.5e',        itn,   x(1) );
            str2 = sprintf( ' %10.3e %10.3e', r1norm, r2norm );
            str3 = sprintf( '  %8.1e %8.1e',   test1,  test2 );
            str4 = sprintf( ' %8.1e %8.1e',    anorm,  acond );
            disp([str1 str2 str3 str4])
         end
      end
      if istop > 0, break, end
end

%     End of iteration loop.
%     Print the stopping condition.

if show
   disp(' ')
   disp('LSQR finished')
   disp(msg(istop+1,:))
   disp(' ')
   str1 = sprintf( 'istop =%8g   r1norm =%8.1e',   istop, r1norm );
   str2 = sprintf( 'anorm =%8.1e   arnorm =%8.1e', anorm, arnorm );
   str3 = sprintf( 'itn   =%8g   r2norm =%8.1e',     itn, r2norm );
   str4 = sprintf( 'acond =%8.1e   xnorm  =%8.1e', acond, xnorm  );
   disp([str1 '   ' str2])
   disp([str3 '   ' str4])
   disp(' ')
end

%-----------------------------------------------------------------------
% End of lsqr.m
%-----------------------------------------------------------------------

function z = Aprod(x,mode)
   if mode == 1
      if isnumeric(A), z = A*x;
      else             z = A(x,1);
      end
   else
      if isnumeric(A), z = (x'*A)';
      else             z = A(x,2);
      end
   end
end % function Aprod

end

author: mpf


text:
classdef EMObject < handle
	% EMObject is the combination of Shape and Material.

	properties (SetAccess = immutable)
		shape
		material
	end

	methods
		function this = EMObject(shape, material)
			if nargin ~= 0
				chkarg(istypesizeof(shape, 'Shape', [1 0]), '"shape" should be row vector with Shape as elements.');
				chkarg(istypesizeof(material, 'Material'), ...
					'"material" should be instance of Material.');

				n = length(shape);
				this(n) = EMObject;
				for i = 1:n
					this(i).material = material;
					this(i).shape = shape(i);
				end
			end
		end

		% A function assigning eps and mu could have been a member function of
		% EMObject, but there are too many subpixel smoothing schemes for eps and
		% mu.  Suppoting all those schemes would make EMObject very heavyweight. A
		% decision of taking out those schemes from EMObject was made to leave
		% EMObject lightweight.
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%% <../index.html MATLAB xUnit Test Framework>: How to Run Tests Silently and Query the Results
% When you run a test suite using |runtests|, the results are
% summarized in the Command Window.  This example shows you how to
% run a test suite so that nothing prints to the Command Window, and
% it shows you how to write a program to automatically determine the
% results of running the test suite.
%
% There are four steps to follow.
%
% 1. Construct a |TestSuite| object.  In this example we'll use the |fromPwd|
% method of the |TestSuite| class to construct a test suite using all the test
% cases found in the |examples_general| directory.

cd examples_general
suite = TestSuite.fromPwd();

%%
% You can look up information about the individual test cases.

suite.TestComponents{1}

%%
% You can see above that the first test component in the test suite is itself
% another test suite, which contains the test cases defined by the M-file named
% TestUsingTestCase. Here's what one of these individual test cases looks like:

suite.TestComponents{1}.TestComponents{1}

%%
% 2. Construct a TestLogger object.  This object can receive
% notifications about what happens when a test suite is executed.

logger = TestRunLogger;

%%
% 3. Call the |run| method of the |TestSuite| object, passing it the
% logger.

suite.run(logger);

%%
% The |TestLogger| object can now be queried to determine what
% happened during the test.

logger

%%
% There were eight test cases run (logger.NumTestCases), resulting in
% one test failure and one test error.  Detailed information about
% what went wrong can be found in |logger.Faults|.

logger.Faults(1)

%%

logger.Faults(2)

%%
% You can drill further to determine the names of the failing tests,
% as well as the complete stack trace associated with each failure.

logger.Faults(1).TestCase

%%

logger.Faults(1).Exception.stack(1)

%%

logger.Faults(1).Exception.stack(2)

%%
% <../index.html Back to MATLAB xUnit Test Framework>

%%
% Copyright 2008-2010 The MathWorks, Inc.

author: ricedsp


text:
function nav=adjnav(nav,opt)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%adjust the wavelength of priority frequencies for BDS2 and BDS3
%eph and geph struct to eph and geph matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global glc

%adjust eph
if nav.n>0
    eph0=zeros(nav.n,40);
    for i=1:nav.n
        eph=nav.eph(i);
        eph0(i,:)=[eph.sat, eph.iode, eph.iodc, eph.sva, eph.svh, eph.week, eph.code, eph.flag, eph.toc.time, eph.toc.sec,...
            eph.toe.time, eph.toe.sec, eph.ttr.time, eph.ttr.sec, eph.A, eph.e, eph.i0, eph.OMG0, eph.omg, eph.M0,...
            eph.deln, eph.OMGd, eph.idot, eph.crc, eph.crs, eph.cuc, eph.cus, eph.cic, eph.cis, eph.toes,...
            eph.fit, eph.f0, eph.f1, eph.f2, eph.tgd, eph.Adot, eph.ndot];
    end
    nav=rmfield(nav,'eph');
    nav.eph=eph0;
end

%adjust geph
if nav.ng>0
    geph0=zeros(nav.ng,22);
    for i=1:nav.ng
        geph=nav.geph(i);
        geph0(i,:)=[geph.sat, geph.iode, geph.frq, geph.svh, geph.sva, geph.age, geph.toe.time, geph.toe.sec, geph.tof.time, geph.tof.sec,...
            geph.pos', geph.vel', geph.acc', geph.taun,...
            geph.gamn, geph.dtaun];
    end
    nav=rmfield(nav,'geph');
    nav.geph=geph0;
end

% adjust wavelength
bds_frq_flag=1;
lam=nav.lam; nav=rmfield(nav,'lam'); nav.lam=zeros(glc.MAXSAT,glc.NFREQ);
if glc.NFREQ>3&&(size(opt.bd2frq,2)<=3||size(opt.bd3frq,2)<=3)&&~isempty(strfind(opt.navsys,'C'))
    bds_frq_flag=0;
    fprintf('Warning:Specified frequency of BDS less than used number of frequency!\n');
end
for i=1:glc.MAXSAT
    [sys,prn]=satsys(i);
    if sys==glc.SYS_BDS
        if prn<19 %BD2
            if ~bds_frq_flag,continue;end
            frq=opt.bd2frq;
            for j=1:glc.NFREQ
               nav.lam(i,j)=lam(i,frq(j));
            end
        else %BD3
            if ~bds_frq_flag,continue;end
            frq=opt.bd3frq;
            for j=1:glc.NFREQ
                nav.lam(i,j)=lam(i,frq(j));
            end
        end
    else %GPS GLO GAL QZS
        nav.lam(i,:)=lam(i,1:3);
    end
end

return

repo: GINav


text:
% *********************************************************************
% Identify and plot block boundaries
% (requires job object with calculated variant Ids)
% *********************************************************************
home; close all; clearvars -except job
%% Application - Determine Variant pairing
%Derive variant types 1 to 6 from variants 1 to 24
laths = job.transformedGrains;
laths.prop.variantType = laths.variantId - (laths.packetId-1) * max(job.variantId)/max(job.packetId);
%Get the boundary segments and neighbouring laths
lathboundaries = laths.boundary(job.csChild,job.csChild);
boundaryIds = lathboundaries.grainId;

[c,ia] = ismember(boundaryIds,laths.id);
boundaryIds(any(~c,2),:) = [];
lathboundaries(any(~c,2)) = [];
%Identify special boundaries, here V1-V2, V1-V3(V5), V1-V6 and V1-V4
%[S. Morito, A.H. Pham, T. Hayashi, T. Ohba, Mater. Today Proc. 2 (2015)
%S913–S916.]
varTypes = laths(laths.id2ind(boundaryIds)).variantType;
cond(1,:) = any(ismember(varTypes,1),2) & any(ismember(varTypes,2),2) | ...
            any(ismember(varTypes,3),2) & any(ismember(varTypes,4),2) | ...
            any(ismember(varTypes,5),2) & any(ismember(varTypes,6),2);     %V1-V2
cond(2,:) = any(ismember(varTypes,1),2) & any(ismember(varTypes,3),2) | ...
            any(ismember(varTypes,1),2) & any(ismember(varTypes,5),2) | ...
            any(ismember(varTypes,2),2) & any(ismember(varTypes,4),2) | ...
            any(ismember(varTypes,2),2) & any(ismember(varTypes,6),2) | ...
            any(ismember(varTypes,3),2) & any(ismember(varTypes,5),2) | ...
            any(ismember(varTypes,4),2) & any(ismember(varTypes,6),2);     %V1-V3(V5)
cond(3,:) = any(ismember(varTypes,1),2) & any(ismember(varTypes,6),2) | ...
            any(ismember(varTypes,2),2) & any(ismember(varTypes,3),2) | ...
            any(ismember(varTypes,4),2) & any(ismember(varTypes,5),2);     %V1-V6
cond(4,:) = any(ismember(varTypes,1),2) & any(ismember(varTypes,4),2) | ...
            any(ismember(varTypes,2),2) & any(ismember(varTypes,5),2) | ...
            any(ismember(varTypes,3),2) & any(ismember(varTypes,6),2);     %V1-V4

%% Application - Plot a map of variant pairing
%Plot the boundaries
colors = {'r','g','b','k'};
ebsdC = job.ebsdPrior;
try data = ebsdC.bc; catch data = ebsdC.imagequality; end
figure; plot(ebsdC,data);
setColorRange([mean(data)-2*std(data),mean(data)+2*std(data)])
mtexColorMap black2white
hold on
labelstr = {'V1-V2','V1-V3(V5)','V1-V6','V1-V4'};
for i = 1:size(cond,1)
    plot(lathboundaries(cond(i,:)),'linewidth',2,'linecolor',colors{i},'DisplayName',labelstr{i});
end
legend

%% Application - Plot the fractions of variant pairing
%Plot bar diagram to show the special boundary fractions
for i = 1:size(cond,1)
    frac(i) = sum(lathboundaries(cond(i,:)).segLength)/sum(lathboundaries(any(cond)).segLength);
end
figure;
bar(categorical(labelstr),frac,'r');
ylabel('Fraction of block boundaries')

author: mtex-toolbox

repo: mtex

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


text:
function varargout = crosslegtable(X)
% Cross-leg table function
%
%   CROSSLEGTABLE([x1, x2]) returns the value of the Cross-legged
%   table function at the specified points. [x1] and [x2] may be vectors.
%   The search domain is
%
%               -10 < x_i < 10
%
%   The global minimum is found on the planes x = 0 and y = 0, with
%
%                   fmin = -1.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;  % # dims
        varargout{2} = [-10, -10]; % LB
        varargout{3} = [+10, +10]; % UB
        varargout{4} = NaN; % solution (too complicated)
        varargout{5} = -1; % function value at solution

    % otherwise, output function value
    else

        % keep values in the serach interval
        X(X < -10) = inf;     X(X > 10) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = -(abs(sin(x1).*sin(x2).*exp(abs(100 - sqrt(x1.^2 + x2.^2)/pi))) + 1).^(-0.1);

    end

end

author: Sable


text:
function test_headmodel_interpolate

% WALLTIME 00:45:00
% MEM 6gb
% DEPENDENCY mesh_sphere ft_headmodeltype ft_headmodel_interpolate ft_prepare_vol_sens ft_compute_leadfield leadfield_interpolate ft_apply_transform

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% create a set of electrodes that nicely covers the upper half of a sphere
[pnt, tri] = mesh_sphere(162);
pnt = pnt .* 10; % convert to cm
sel = find(pnt(:,3)>0);
elec1 = [];
elec1.elecpos = pnt(sel,:);
for i=1:length(sel)
  elec1.label{i} = sprintf('elec%d', i);
end
elec1.unit = 'cm';

% create another set of electrodes covering the upper half of a sphere
pnt = elec1.elecpos(1:42,:) + randn(42,3);
elec2 = [];
for i=1:size(pnt,1)
  elec2.elecpos(i,:) = pnt(i,:)./norm(pnt(i,:)) * 10;
  elec2.label{i} = sprintf('rand%d', i);
end
elec2.unit = 'cm';

% create another set of electrodes representing a bipolar montage
bipolar = [];
bipolar.labelold = elec1.label;
bipolar.tra = zeros(length(elec1.label)-1, length(elec1.label));
for i=1:(length(bipolar.labelold)-1)
  bipolar.labelnew{i} = sprintf('%s-%s', bipolar.labelold{i}, bipolar.labelold{i+1});
  bipolar.tra(i,i  ) = +1;
  bipolar.tra(i,i+1) = -1;
end
elec3 = ft_apply_montage(elec1, bipolar);

% create an identical set of electrodes but with other channel names
elec4 = elec1;
for i=1:length(elec4.label)
  elec4.label{i} = sprintf('chan%d', i);
end

% update the electrode sets to the latest standards
elec1 = ft_datatype_sens(elec1);
elec2 = ft_datatype_sens(elec2);
elec3 = ft_datatype_sens(elec3);
elec4 = ft_datatype_sens(elec4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% construct a singlesphere volume conduction model
volA   = [];
volA.c = 1;
volA.r = 10;
volA.o = [0 0 0];

% compute the precomputed leadfield
cfg = [];
cfg.headmodel = volA;
cfg.elec = elec1;
cfg.resolution = 1;
leadfield = ft_prepare_leadfield(cfg);

% remember one position
firstindx = find(leadfield.inside, 1, 'first');
pos1 = leadfield.pos(firstindx,:);
pos2 = [0 0 7];

% remember the precomputed leadfield for the first dipole
lf1 = leadfield.leadfield{firstindx};

% the following call generates the volume conduction model and writes all files to disk
filename = fullfile(tempname, 'leadfield');
ft_headmodel_interpolate(filename, elec1, leadfield, 'smooth', false);

% the next day you would start by reading it from disk
volB = ft_read_headmodel([filename '.mat']); % this is a mat file containing the "vol" structure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use the same electrodes
[volAA, elecAA] = ft_prepare_vol_sens(volA, elec1);
[volBB, elecBB] = ft_prepare_vol_sens(volB, elec1);

% compare the leadfields
lfa = ft_compute_leadfield(pos1, elecAA, volAA); % the original
lfb = ft_compute_leadfield(pos1, elecBB, volBB); % the interpolation

assert(all(mean(lfa, 1)<eps(double(1))), 'the leadfield is not average referenced'); % this is in double precision
assert(all(mean(lfb, 1)<eps(single(1))), 'the leadfield is not average referenced'); % this is in single precision
assert(isalmostequal(lf1, lfb, 'reltol', 1e-4), 'the leadfields are different');
assert(isalmostequal(lfa, lfb, 'reltol', 1e-4), 'the leadfields are different');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use a subset of electrodes
[volAA, elecAA] = ft_prepare_vol_sens(volA, elec1, 'channel', elec1.label(1:50));
[volBB, elecBB] = ft_prepare_vol_sens(volB, elec1, 'channel', elec1.label(1:50));

% compare the leadfields
lfa = ft_compute_leadfield(pos1, elecAA, volAA); % the original
lfb = ft_compute_leadfield(pos1, elecBB, volBB); % the interpolation

assert(all(mean(lfa, 1)<eps), 'the leadfield is not average referenced');
assert(all(mean(lfb, 1)<eps), 'the leadfield is not average referenced');
assert(isalmostequal(lfa, lfb, 'reltol', 1e-3), 'the leadfields are different');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use the same electrodes with with different names
[volAA, elecAA] = ft_prepare_vol_sens(volA, elec4);
[volBB, elecBB] = ft_prepare_vol_sens(volB, elec4);

% compare the leadfields
lfa = ft_compute_leadfield(pos1, elecAA, volAA); % the original
lfb = ft_compute_leadfield(pos1, elecBB, volBB); % the interpolation

assert(all(mean(lfa, 1)<eps), 'the leadfield is not average referenced');
assert(all(mean(lfb, 1)<eps), 'the leadfield is not average referenced');
assert(isalmostequal(lfa, lfb, 'reltol', 1e-4), 'the leadfields are different');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use electrodes with a different placement
[volAA, elecAA] = ft_prepare_vol_sens(volA, elec2);
[volBB, elecBB] = ft_prepare_vol_sens(volB, elec2);

% compare the leadfields
lfa = ft_compute_leadfield(pos2, elecAA, volAA); % the original
lfb = ft_compute_leadfield(pos2, elecBB, volBB); % the interpolation

assert(all(mean(lfa, 1)<eps), 'the leadfield is not average referenced');
assert(all(mean(lfb, 1)<eps), 'the leadfield is not average referenced');
c = corrcoef(lfa(:), lfb(:)); c = c(1,2);
assert(c>0.8, 'the leadfields are different');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use a bipolar electrode montage
[volAA, elecAA] = ft_prepare_vol_sens(volA, elec3);
[volBB, elecBB] = ft_prepare_vol_sens(volB, elec3);

% compare the leadfields
lfa = ft_compute_leadfield(pos2, elecAA, volAA); % the original
lfb = ft_compute_leadfield(pos2, elecBB, volBB); % the interpolation

assert(~all(mean(lfa, 1)<eps), 'the leadfield is average referenced'); % it should NOT be average referenced
assert(~all(mean(lfb, 1)<eps), 'the leadfield is average referenced'); % it should NOT be average referenced
assert(isalmostequal(lfa, lfb, 'abstol', inf, 'reltol', inf, 'relnormtol', 1e-3), 'the leadfields are different');


text:
function b = isPointOnRay(point, ray, varargin)
%ISPOINTONRAY Test if a point belongs to a ray
%
%   B = isPointOnRay(PT, RAY);
%   Returns 1 if point PT belongs to the ray RAY.
%   PT is given by [x y] and RAY by [x0 y0 dx dy].
%
%   If PT is a N-by-2 array, and RAY is a M-by-4 array, then the result is
%   a N-by-M array containing the result of each pair-wise test.
%
%   B = isPointOnRay(PT, RAY, TOL);
%   Specifies the tolerance to use for testing if point is on the ray.
%
%   Example
%     ray = [10 20 3 4];
%     % test for a point on the ray
%     p1 = [16 28];
%     isPointOnRay(p1, ray)
%     ans =
%       logical
%        0
%     % test for a point on the supporting line but "before" the origin
%     p2 = [7 16];
%     isPointOnRay(p1, ray)
%     ans =
%       logical
%        0
%
%   See also:
%   rays2d, points2d, isPointOnLine, isPointOnEdge
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   07/07/2005 normalize condition to test if on the line and add support
%       of multiple rays or points
%   22/05/2009 rename to isPointOnRay, add psb to specify tolerance
%   26/01/2010 was drawing a line before making test

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% number of rays and points
Nr = size(ray, 1);
Np = size(point, 1);

% if several rays or several points, adapt sizes of arrays
x0 = repmat(ray(:,1)', Np, 1);
y0 = repmat(ray(:,2)', Np, 1);
dx = repmat(ray(:,3)', Np, 1);
dy = repmat(ray(:,4)', Np, 1);
xp = repmat(point(:,1), 1, Nr);
yp = repmat(point(:,2), 1, Nr);

% test if points belongs to the supporting line
b1 = abs((xp-x0).*dy - (yp-y0).*dx) ./ (dx.*dx + dy.*dy) < tol;

% check if points lie the good direction on the rays
ind     = abs(dx) > abs(dy);
t       = zeros(size(b1));
t(ind)  = (xp(ind) - x0(ind)) ./ dx(ind);
t(~ind) = (yp(~ind) - y0(~ind)) ./ dy(~ind);

% combine the two tests
b = b1 & (t >= 0);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function d = double(npd)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% double(kde) -- evaluates true (1) if the kde is non-empty
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

  if (npd.N > 0) d = 1;            % return 1 if the density exists
  else d = 0; end;

author: ShapeNet

repo: RenderForCNN


text:
% Lawn sprinker example from Russell and Norvig p454
% See www.cs.berkeley.edu/~murphyk/Bayes/usage.html for details.

rand('state', 0);
randn('state', 0);

N = 4;
dag = zeros(N,N);
C = 1; S = 2; R = 3; W = 4;
dag(C,[R S]) = 1;
dag(R,W) = 1;
dag(S,W)=1;

false = 1; true = 2;
ns = 2*ones(1,N); % binary nodes

bnet = mk_bnet(dag, ns);
bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);

nsamples = 500;
samplesM = cell(N, nsamples);
for i=1:nsamples
  samplesM(:,i) = sample_bnet(bnet);
end

hide = rand(N, nsamples) > 0.9;
[I,J]=find(hide);
for k=1:length(I)
  samplesM{I(k), J(k)} = [];
end

% Make a initial chain like dag
G0 = zeros(N,N);
for i=1:N-1
   G0(i, i+1) = 1;
end

figure;
draw_graph(G0);

B0 = mk_bnet(G0, ns);
% use random params
for i=1:N
  B0.CPD{i} = tabular_CPD(B0, i, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
end

max_loop = 30;
%profile on -detail mmex
[B0, order, best_score] = learn_struct_EM(B0, samplesM, max_loop);
%profile report

dag1 = B0.dag;
dag1 = dag1(order,order);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [this, saveFileName] = save(this, varargin)
% Saves 4D MrImage in single file.
%
%   Y = MrImageSpm4D()
%   Y.save('fileName', fileName)
%
% This is a method of class MrImageSpm4D.
%
% IN
%   fileName    string: default via get_filename
% OUT
%
% EXAMPLE
%   save
%
%   See also MrImageSpm4D MrDataNd.save MrDataNd.get_filename

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2018-05-23
% Copyright (C) 2018 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

defaults.fileName = this.get_filename();
args = tapas_uniqc_propval(varargin, defaults);

tapas_uniqc_strip_fields(args);

saveFileName = this.write_single_file(fileName);
end

author: translationalneuromodeling

repo: tapas


text:
function [graph,nodeNames,nodeInfo] = prtUtilReadGml(file)
%[graph,names,etc] = prtUtilReadGml(file)
%  Read graph-spec files with .gml formating.

etc = [];

fid = fopen(file,'r');
directed = false;

nodeList = [];
edgeList = [];
while ~feof(fid)
    str = fgetl(fid);
    str = lower(str);

    if ~isempty(strfind(str,'directed'))
        temp = regexp(str,'directed (?<directed>[0-9]*)','names');
        directed = logical(str2double(temp.directed));
    end

    if ~isempty(strfind(str,']')) || ~isempty(strfind(str,'['))
        continue;
    else
        if ~isempty(strfind(str,'node'))
            newNode = readNode(fid);
            if isempty(nodeList)
                nodeList = newNode;
            else
                nodeList = cat(1,nodeList,newNode);
            end
        end
        if ~isempty(strfind(str,'edge'));
            newEdge = readEdge(fid);
            if isempty(edgeList)
                edgeList = newEdge;
            else
                edgeList = cat(1,edgeList,newEdge);
            end
        end
    end
end
fclose all;

nodeNames = {nodeList.id}';
nodeVals = cat(1,nodeList.val);
edgeVals = cat(1,edgeList.connection);

edgeIndices = zeros(size(edgeVals));
for i = 1:size(edgeVals,1)
    edgeIndices(i,:) = [find(edgeVals(i,1) == nodeVals),find(edgeVals(i,2) == nodeVals)];
end

graph = sparse(edgeIndices(:,1),edgeIndices(:,2),1,length(nodeVals),length(nodeVals));
if ~directed
    graph = graph + sparse(edgeIndices(:,2),edgeIndices(:,1),1,length(nodeVals),length(nodeVals));
end
nodeInfo = nodeList;

%
function newEdge = readEdge(fid)

newEdge = struct('connection',[]);
while true
    str = fgetl(fid);
    if ~isempty(strfind(str,'['))

    elseif ~isempty(strfind(str,']'))
        return;
    elseif ~isempty(strfind(str,'source'))
        temp1 = regexp(str,'source (?<nodeNumber>[0-9].*)','names');
        int1 = str2double(temp1.nodeNumber);

        str = fgetl(fid);
        temp2 = regexp(str,'target (?<nodeNumber>[0-9].*)','names');
        int2 = str2double(temp2.nodeNumber);
        newEdge.connection = [int1,int2];
    end
end

function newNode = readNode(fid)

newNode = struct('val',[],'id','');
while true
    str = fgetl(fid);
    str = strtrim(str);
    if ~isempty(strfind(str,'['))

    elseif ~isempty(strfind(str,']'))
        return;
    elseif ~isempty(regexp(str,'^id','once'))
        temp = regexp(str,'id (?<nodeNumber>[0-9].*)','names');
        int = str2double(temp.nodeNumber);
        newNode.val = int;
        newNode.id = sprintf('Node %d',int);
    elseif ~isempty(regexp(str,'^label','once'))
        temp = regexp(str,'label (?<nodeLabel>["A-Z].*)','names');
        temp.nodeLabel = strrep(temp.nodeLabel,'"','');
        newNode.id = temp.nodeLabel;
    elseif ~isempty(regexp(str,'^value','once'))
        temp = regexp(str,'value (?<value>.*)','names');
        newNode.value = str2double(temp.value);
    else
        warning('I didnt know what to do');
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function c=vertcat(varargin)
%
%           vertical concatination of the psfMatrix objects
%
%  Input:  any number of psfMatrix objects, or a multiPsfMatrix object
%
%  Output:  c is a multiPsfMatrix object of size (number of argin) by 1
%

%  K. Lee  2/28/02

if isa( varargin{1}, 'cell')
  A=varargin{1};  			   % A is a cell
  A=A{1};         			   % A is a psfMatrix object
  bdry=A.boundary;
  Apsf=A.psf;        			   % Apsf is a cell
  [r,d]=imageSize(A);
else
  A=varargin{1};
  bdry=A.boundary;
  Apsf=A.psf;
  [r,d]=imageSize(A);
end
c=multiPsfMatrix(varargin{1});
for i=2:nargin
  A=varargin{i};
  [r1,c1]=imageSize(A);
  bdry2=A.boundary;
    if ~isa(varargin{i},'psfMatrix'),
      error('psfMatrix objects can only be concatenated with other psfMatrix objects')
  elseif  (strcmp(bdry,bdry2)==1)
        c=multiPsfMatrix(c,varargin{i});
    else
      error('psfMatrix objects must have the same boundary')

    end
end

repo: IRtools


text:
function [out] = rainfall_2(In,T,p1,p2)
%rainfall_2

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Rainfall based on a temperature threshold interval
% Constraints:  -
% @(Inputs):    In   - incoming precipitation flux [mm/d]
%               T    - current temperature [oC]
%               p1   - midpoint of the combined rain/snow interval [oC]
%               p2   - length of the mixed snow/rain interval [oC]

out = min(In,max(0,In.*(T-(p1-0.5*p2))/p2));

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function varargout = Read_MapGUI(varargin)
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Read_MapGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @Read_MapGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before Read_MapGUI is made visible.
function Read_MapGUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);
% --- Outputs from this function are returned to the command line.
function varargout = Read_MapGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% --- Executes on button press in browser_button.
function browser_button_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile('*.mat', 'Pick a .mat file');      %Browse for file
handles.fullpathname=strcat(pathname, filename);                    %Connect path name and file name to one string
handles.filename=filename;
set(handles.continue_button, 'Enable', 'Inactive');                 %Disable Continue button
data_type=lower(filename(end-2:end));                               %Extract data type and convert to string
if data_type =='mat'                                                 %#ok<BDSCA>
    set(handles.text4, 'String', handles.fullpathname);         %Set GUI.full_path to string fullpathname
    if filename(1:5)=="Track"
        set(handles.text4, 'String', 'ERROR');
    else
        set(handles.continue_button, 'Enable', 'On');                   %Enable Continue button
    end
else
    set(handles.text4, 'String', 'ERROR');                      %Display "ERROR"
end
guidata(hObject, handles);

% --- Executes on slider movement.
function yellow_slider_Callback(hObject, eventdata, handles)
a=get(handles.yellow_slider,'Value');                                       %Get slider value of "global slider"
a=round(a*100)/100;                                                               %Round value of "global slider" to 1/50
set(handles.yellow_slider,'Value', a);                                          %Set slider value of "global slider" to rounded value
set(handles.text_yellow,'String',num2str(a));                                   %Set text above "global slider" to slider value
% --- Executes during object creation, after setting all properties.
function yellow_slider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function green_slider_Callback(~, eventdata, handles)
a=get(handles.green_slider,'Value');                                       %Get slider value of "global slider"
a=round(a*100)/100;                                                               %Round value of "global slider" to 1/50
set(handles.green_slider,'Value', a);                                          %Set slider value of "global slider" to rounded value
set(handles.text_green,'String',num2str(a));                                   %Set text above "global slider" to slider value
% --- Executes during object creation, after setting all properties.
function green_slider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function blue_slider_Callback(hObject, eventdata, handles)
a=get(handles.blue_slider,'Value');                                       %Get slider value of "global slider"
a=round(a*100)/100;                                                               %Round value of "global slider" to 1/50
set(handles.blue_slider,'Value', a);                                          %Set slider value of "global slider" to rounded value
set(handles.text_blue,'String',num2str(a));                                   %Set text above "global slider" to slider value
% --- Executes during object creation, after setting all properties.
function blue_slider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function red_slider_Callback(hObject, eventdata, handles)
a=get(handles.red_slider,'Value');                                       %Get slider value of "global slider"
a=round(a*100)/100;                                                               %Round value of "global slider" to 1/50
set(handles.red_slider,'Value', a);                                          %Set slider value of "global slider" to rounded value
set(handles.text_red,'String',num2str(a));                                   %Set text above "global slider" to slider value
% --- Executes during object creation, after setting all properties.
function red_slider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function white_slider_Callback(hObject, eventdata, handles)
a=get(handles.white_slider,'Value');                                       %Get slider value of "global slider"
a=round(a*100)/100;                                                               %Round value of "global slider" to 1/50
set(handles.white_slider,'Value', a);                                          %Set slider value of "global slider" to rounded value
set(handles.text_white,'String',num2str(a));                                   %Set text above "global slider" to slider value
% --- Executes during object creation, after setting all properties.
function white_slider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on button press in continue_button.
function continue_button_Callback(~, eventdata, handles) %#ok<*INUSL>
assignin ('base', 'path_map', handles.fullpathname);         %Assignin referenceline path to workspace
assignin ('base', 'filename_map', handles.filename);
yellow=get(handles.yellow_slider,'Value');
green=get(handles.green_slider,'Value');
blue=get(handles.blue_slider,'Value');
red=get(handles.red_slider,'Value');
white=get(handles.white_slider,'Value');
assignin ('base', 'yellow', yellow);
assignin ('base', 'green', green);
assignin ('base', 'blue', blue);
assignin ('base', 'red', red);
assignin ('base', 'white', white);
set(handles.continue_button, 'Enable', 'Off');
set(handles.browser_button, 'Enable', 'Off');
set(findall(handles.uipanel5, '-property', 'enable'),'enable', 'off');
GIMP_Correction;

author: TUMFTM


text:
function sout = cvx_solver( flag )

%CVX_SOLVER    CVX solver selection.
%   CVX_SOLVER SDPT3  selects SDPT3  as the solver CVX uses to solve models.
%   CVX_SOLVER SeDuMi selects SeDuMi as the solver CVX uses to solve models.
%   The solver string is case-insensitive; e.g., 'sdpt3' is fine.
%
%   The default solver is SeDuMi. For most models, this will be the best
%   choice. Nevertheless, for your problems you may find that SDPT3 is
%   more accurate or reliable, hence the need for this command.
%
%   If CVX_SOLVER is called within a model---that is, between the statements
%   CVX_BEGIN and CVX_END---then the new solver selection applies only to that
%   particular model. In this mode, the change takes place for *both*
%   successive approximation models and normal models.
%
%   If called outside of a model, then the change applies to
%   all subsequent models; that is, it modifies the default solver.
%
%   On exit, CVX_SOLVER returns the *previous* solver. This can be used to
%   restore that solver selection at a later time, in a manner similar to
%   CVX_PRECISION.
%
%   CVX_SOLVER, with no arguments, returns the current solver.

global cvx___
cvx_global
if nargin > 0,
    if isempty( flag ),
        flag = 'sedumi';
    elseif ~ischar( flag ) || size( flag, 1 ) ~= 1,
        error( 'Argument must be a string.' );
    elseif ~isfield( cvx___.path.solvers, lower( flag )  ),
        error( 'Unknown, unusable, or missing solver: %s', flag );
    end
end
if isempty( cvx___.problems ),
    s = cvx___.solver;
    if nargin > 0,
        cvx___.solver = flag;
    end
else
    s = cvx___.problems(end).solver;
    if nargin > 0,
        if ~strcmpi( s, flag ) && ~isa( evalin( 'caller', 'cvx_problem', '[]' ), 'cvxprob' ),
            warning( 'CVX:Solver', 'The global CVX solver selection cannot be changed while a model is being constructed.' );
        else
            cvx___.problems(end).solver = flag;
        end
    end
end
if nargin > 0 && cvx___.path.hold,
    cvx_setspath( flag );
end
if nargin == 0 || nargout > 0,
    sout = s;
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST