text:
function [cl_ext_spm, fwhm] = cl_ext_spm_grf(corrected_p, prim_p, residual_images, mask, varargin)
% This function is designed to estimate a cluster extent size for
% the correction for multiple comparisons based on a Gaussian Random Field
% Theory using SPM toolboxes.
%
% :Usage:
% ::
%
%     [cl_ext, fwhm] = cl_ext_spm_grf(corrected_p, prim_p, residual_images, mask, varargin)
%
% :Inputs:
%
%   **corrected_p:**
%        corrected p value
%
%        e.g.) cluster-extent corrected p < .05: corrected_p = .05
%
%   **prim_p:**
%        primary threshold for height (i.e., cluster-defining threshold)
%
%        e.g.) prim_p = [0.01 0.005 0.001];
%
%   **residual_images:**
%        residual image names; if you used
%
%        cl_ext_make_resid.m, this should be 'Res4d.nii'
%
%   **mask:**
%        mask image name
%
% :Optional Inputs:
%
%   **'doplot':**
%
%   **'twotail':**
%        default is one-tail - with this option, primary_p/2 will be used
%        for all clsuter extent estimations.
%
% Output:
%
%   **cl_ext_spm:**
%        cl_ext_spm is the cluster size that makes a corrected p value under
%       corrected_p (e.g., 0.05).
%
%   **fwhm (x, y, z in voxels):**
%        intrinsic smoothness level estimated by SPM (spm_est_smoothness.m)
%       If you want to convert this into mm, you need to multiply these
%       values by voxel sizes in mm.
%
% ..
%    Choong-Wan (Wani) Woo, 08/13/2012
%    modified by Wani, 05/18/2013
% ..

doplot = false;
isTwoTailed = false;

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'doplot', doplot = true;
            case 'twotail', isTwoTailed = true;
        end
    end
end

con_num = size(expand_4d_filenames(residual_images),1);
ndf = [1 con_num-1];

[fwhm, dummy, r] = spm_est_smoothness(residual_images, mask, [con_num con_num-1]);

V2R = 1/prod(fwhm);

cl_ext = zeros(length(prim_p),2);

if doplot
    scrsz = get(0, 'ScreenSize');
    create_figure('cluster_extent_spm');
    set(gcf, 'Position', [1 scrsz(4)/2 scrsz(3)/1.5 scrsz(4)/2]);
end

for i = 1:length(prim_p)
    if isTwoTailed
        u(i) = spm_u(prim_p(i)/2, ndf, 'T');
    else
        u(i) = spm_u(prim_p(i), ndf, 'T');
    end
end

for i = 1:length(prim_p) % cluster-extent threshold (k) based on RF in SPM

    P = [];

    for j = 1:100000
        k = j*V2R; % convert the number of voxels into the number of resels

        P(end+1,1) = spm_P_RF(1,k,u(i),ndf,'T',r,1);
        P(end,2) = j;

        if P(end,1) <= corrected_p
            cl_ext(i,1) = j;
            cl_ext(i,2) = P(end,1);
            break
        end
    end

    if doplot
        hh(i) = subplot(1,length(prim_p),i);
        plot(P(:,2), P(:,1), '-b', 'LineWidth', 1.5);
        xlabel('cluster extent size', 'FontSize', 16);
        ylabel('corrected P value', 'FontSize', 16);
        eval(['title(hh(i), ''Primary P:' num2str(prim_p(i)), ''', ''fontsize'', 16);']);
        set(gca, 'FontSize', 15);
        hold on;
        plot(cl_ext(i,1), cl_ext(i,2), 'r+', 'MarkerSize', 15);
        eval(['text(cl_ext(i,1)/2, cl_ext(i,2), ''+ cl size:' num2str(cl_ext(i,1)) ''', ''FontSize'', 14);']);
    end

end

cl_ext_spm = cl_ext(:,1);

return

author: canlab

repo: CanlabCore


text:
function k = ratquadKernDiagCompute(kern, x)

% RATQUADKERNDIAGCOMPUTE Compute diagonal of RATQUAD kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the rational quadratic kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : ratquadKernParamInit, kernDiagCompute, kernCreate, ratquadKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

k = repmat(kern.variance, size(x, 1), 1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Ft = CIELabFunction(t, inverse)
%
%       Ft = CIELabFunction(t, inverse)
%
%
%        Input:
%           -t: image
%           -inverse: takes as values 0 or 1. If it is set to 1 the
%                     transformation from XYZ to CIE Lab is applied, otherwise
%                     the transformation from CIE Lab to XYZ
%
%        Output:
%           -Ft: application of CIE Lab f or f^{-1} (if inverse = 1) to X
%
%     Copyright (C) 2013  Francesco Banterle
%
%     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/>.
%

if(inverse == 0) %forward function
    Ft = zeros(size(t));

    c1 = (6 / 29)^3;
    c2 = ((29 / 6)^2)/3;
    c3 = 4 / 29;

    Ft(t >  c1) = t(t >  c1).^(1 / 3);
    Ft(t <= c1) = t(t <= c1) * c2 + c3;
end

if(inverse == 1) %inverse function
    Ft = zeros(size(t));

    c1 = 6 / 29;
    c2 = ((6 / 29)^2) * 3;
    c3 = 4 / 29;

    Ft(t >  c1) =  t(t >  c1).^3;
    Ft(t <= c1) = (t(t <= c1) - c3) * c2;
end

end

author: banterle


text:
function linpack_d_test25 ( )

%*****************************************************************************80
%
%% TEST25 tests DSIFA and DSISL.
%
%  Discussion:
%
%    DSIFA and DSISL are for symmetric indefinite matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST25\n' );
  fprintf ( 1, '  For a symmetric indefinite matrix,\n' );
  fprintf ( 1, '  DSIFA factors the matrix,\n' );
  fprintf ( 1, '  DSISL solves a factored linear system,\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to the matrix A and the right hand side B.
%
  b(1:n-1) = 0.0;
  b(n)= n + 1;
%
%  Force B to be a column vector.
%
  b = b';

  a(1:n,1:n) = 0.0;

  for i = 1 : n
    a(i,i) = 2.0;
    if ( i < n )
      a(i,i+1) = -1.0;
    end
  end
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, ipvt, info ] = dsifa ( a, lda, n );

  if ( info ~= 0 )
    fprintf ( 1, '  Error!  DSIFA returns INFO = %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  b = dsisl ( a, lda, n, ipvt, b );
%
%  Print the result.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last five entries of the solution:\n' );
  fprintf ( 1, '  (Should be 1,2,3,4,5,...,n-1,n):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function kern = nddisimKernParamInit(kern)

% NDDISIMKERNPARAMINIT NDDISIM kernel parameter initialisation.
% The driven input single input motif (DISIM) kernel is specifically designed for
% working with gene networks where there is assumed to be a single
% transcription factor controlling several genes. This transcription
% factor, in turn, is driven by its own gene's RNA. The model takes the
% following form: each gene is
% related to the transcription factor through the following
% differential equation,
%
% dx(t)/dt = B + C f(t-delta) - D x(t),
%
% where D is a decay term, C is a response term, delta is a time delay
% and B is an initial level. Then if f(t) is assumed to be the result of
% a further differential equation,
%
% df(t)/dt = Sx'(t) - D' x'(t)
%
% where x'(t) is assumed to come from a Gaussian process with an RBF
% covariance function f(t) is a Gaussian process with a covariance function
% provided by the single input motif kernel (SIM) and x(t) is a Gaussian
% process with covariance function provided by this kernel, the DISIM kernel.
%
% The kernel is designed to interoperate with the multiple output
% block kernel so that f(t) can be inferred given several different
% instantiations of x(t) (associated with different genes).
%
% The parameters (B, C, delta, S, D and D') are constrained positive.
%
% FORMAT
% DESC initialises the single input motif
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit, simKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007, 2009
%
% COPYRIGHT: Jaakko Peltonen, 2011

% KERN

if kern.inputDimension > 1
  error('NDDISIM kernel is only valid for one-dimensional input.')
end

if isfield(kern, 'options') && isfield(kern.options, 'gaussianInitial') && ...
      kern.options.gaussianInitial,
  kern.gaussianInitial = 1;
  kern.initialVariance = 1;
else
  kern.gaussianInitial = 0;
end

kern.inverseWidth = 1;
kern.di_variance = 1;
kern.decay = 1;
kern.variance = 1;
kern.delay = 1e-7;

if kern.gaussianInitial,
  kern.nParams = 6;
else
  kern.nParams = 5;
end

if isfield(kern, 'options') && isfield(kern.options, 'paramTransform'),
  paramTransform = kern.options.paramTransform;
else
  paramTransform = 'sigmoidab';
end

switch paramTransform,
 case 'sigmoidab',
  for k=1:kern.nParams,
    kern.transforms(k).index = k;
    kern.transforms(k).type = 'sigmoidab';
    kern.transforms(k).transformsettings = [0 1e6];
  end;
 case 'bounded',
  for k=1:kern.nParams,
    kern.transforms(k).index = k;
    kern.transforms(k).type = optimiDefaultConstraint('bounded');
    kern.transforms(k).transformsettings = [0 1e6];
  end;
 case 'identity',
  for k=1:kern.nParams,
    kern.transforms(k).index = k;
    kern.transforms(k).type = 'identity';
    kern.transforms(k).transformsettings = [0 1e6];
  end;
 case 'positive',
  for k=1:kern.nParams,
    kern.transforms(k).index = k;
    kern.transforms(k).type = optimiDefaultConstraint('positive');
  end;
 case 'none',
 otherwise,
  error('Unknown paramTransform');
end

kern.isStationary = false;


text:
function [Series, Sigmas] = garchsim(parameters, model, distr, p, q, H, y, NumSamples, NumPaths)
%{
-----------------------------------------------------------------------
 PURPOSE:
 Simulation of GARCH models
-----------------------------------------------------------------------
  USAGE:
  [Series, Sigmas] = garchsim(parameters, model, distr, p, q, H, NumSamples, NumPaths)

 INPUTS:
 parameters:    a vector of parameters (i.e. constant, ARCH, GARCH, Gamma, Delta, DoF)
 model:         'GARCH', 'GJR', 'EGARCH', 'NARCH', 'NGARCH, 'AGARCH', 'APGARCH',
                'NAGARCH'
 distr:         'GAUSSIAN', 'T', 'GED', 'CAUCHY', 'HANSEN' and 'GC'
 p:             positive scalar integer representing the order of ARCH
 q:             positive scalar integer representing the order of GARCH
 H:             a vector of time series of positive pre-sample
                conditional standard deviations.
 y:             a vector of factors for the volatility process, must be positive!
 NumSamples:    The number of samples. If is left empty 100 samples will be generated.
 NumPaths:      The number of paths. If it is left empty 1 path will be estimated.

 OUTPUTS:
 Series:       (NumSamples x NumPaths) simulated series with GARCH variances
 Sigmas:       (NumSamples x NumbPaths) vector of conditional standard deviations.

-----------------------------------------------------------------------
 Author:
 Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
 Date:     08/2011
-----------------------------------------------------------------------
%}

if nargin == 0
    error('Parameters, GARCH Model, Distribution, ARCH, GARCH, Sigmas, Number of Samples,  Number of Paths')
end

if nargin < 8 | isempty(NumPaths)
    NumPaths=1;
end

% Verifying that the vector of parameters is a column vector
[r,c]=size(parameters);
if r<c
    parameters=parameters';
end

m  =  max([p,q]);

if strcmp(distr,'GAUSSIAN')
    RandomNumbers = randn(NumSamples+m, NumPaths);
elseif strcmp(distr,'T')
    RandomNumbers = trnd(parameters(end), NumSamples+m, NumPaths);
elseif strcmp(distr,'GED')
     error('Sorry: this distribution is not yet supported, it will be shortly!')
elseif strcmp(distr,'CAUCHY')
      error('Sorry: this distribution is not yet supported, it will be shortly!')
elseif strcmp(distr,'HANSEN')
      error('Sorry: this distribution is not yet supported, it will be shortly!')
elseif strcmp(distr,'GC')
      error('Sorry: this distribution is not yet supported, it will be shortly!')
elseif strcmp(distr,'Logistic')
      error('Sorry: this distribution is not yet supported, it will be shortly!')
elseif strcmp(distr,'Laplace')
      error('Sorry: this distribution is not yet supported, it will be shortly!')
elseif strcmp(distr,'Rayleigh')
      error('Sorry: this distribution is not yet supported, it will be shortly!')
end

resids=[];
h=[];
series=[];

if isscalar(y)
     factors = 0;
   else
     factors = sum(parameters(2+p+q:1+p+q+size(y,2)).*y(end,:));
end

% Set initial value to the unconditional
if isequal(strcat(model), 'GARCH')
   h(1:m,1:NumPaths) = (parameters(1)+factors)/(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q)));
elseif  isequal(strcat(model), 'GJR')
   h(1:m,1:NumPaths) = (parameters(1)+factors)/(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q)) - sum(parameters(4:1+p+q)));
elseif isequal(strcat(model), 'EGARCH')
   h(1:m,1:NumPaths) = exp((parameters(1)+factors)/(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q)) - sum(parameters(4:1+p+q))));
elseif isequal(strcat(model), 'NARCH')
   h(1:m,1:NumPaths) = (parameters(1)+factors)/(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q))); % This is only an approximation
elseif isequal(strcat(model), 'NGARCH')
   h(1:m,1:NumPaths) = nthroot((parameters(1)+factors)/(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q))),parameters(end));
elseif isequal(strcat(model), 'AGARCH')
   h(1:m,1:NumPaths) = (parameters(1)-sum(parameters(2:1+p).*parameters(4:1+p+q))+factors)./(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q)));
elseif isequal(strcat(model), 'APGARCH')
   h(1:m,1:NumPaths) = (parameters(1)+factors)/(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q))); % This is only an approximation
elseif isequal(strcat(model), 'NAGARCH')
    h(1:m,1:NumPaths) = (parameters(1)-sum(parameters(2:1+p).*parameters(4:1+p+q)) + factors)./(1-sum(parameters(2:1+p)) - sum(parameters(3:1+p+q)));
else error('Invalid GARCH Model');
end

if isempty(H)
    series=randn(1:m,NumPaths).*sqrt(h);
else
    series=sqrt(H(end-m+1:end))*ones(1,NumPaths).*RandomNumbers(1:m,:);
end

if strcmp(model,'GARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = parameters(1:1+p+q)'* [1 ; series(t-(1:p),i).^2;  h(t-(1:q),i) ];
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'GJR')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) =  parameters(1:1+2*p+q)'*[1; series(t-(1:p),i).^2; h(t-(1:q),i); (series(t-(1:p),i)<0).*(series(t-(1:p),i).^2)];
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'EGARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = exp(parameters(1:1+2*p+q)'*[1 ; (abs(series(t-(1:p),i))./sqrt(h(t-(1:p),i))-sqrt(2/pi)); log(h(t-(1:q),i)); series(t-(1:p),i)./sqrt(h(t-(1:p),i))]);
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'NARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = parameters(1:1+p+q)'*[1; (abs(series(t-(1:p),i))).^parameters(2+p+q); h(t-(1:q),i)];
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'NGARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = (parameters(1:1+p+q)'*[1; (series(t-(1:p),i)).^parameters(2+p+q); h(t-(1:q),i).^parameters(2+p+q)]).^(1/parameters(2+p+q));
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'AGARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = parameters(1:1+p+q)'*[1; (series(t-(1:p),i) - parameters(2+p+q:1+2*p+q)).^2; h(t-(1:q),i)];
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'APGARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = (parameters(1:1+p+q)'*[1; (abs(series(t-(1:p),i)) - parameters(2+p+q:1+2*p+q).*series(t-(1:p),i)).^parameters(2+2*p+q); h(t-(1:q),i).^parameters(2+2*p+q)]).^(1/parameters(2+2*p+q));
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
elseif strcmp(model,'NAGARCH')
    for i = 1:NumPaths
        for t = (m+1):(size(RandomNumbers,1)-m)
            h(t,i) = parameters(1:1+p+q)'*[1; (series(t-(1:p),i)./sqrt(h(t-(1:p),i)) + parameters(2+p+q:1+2*p+q)).^2; h(t-(1:q),i)];
            series(t,i)=RandomNumbers(t,i)*sqrt(h(t,i));
        end
    end
else error('Invalid GARCH Model');
    end

Series=series;
Sigmas=h;

end

author: Sable

repo: mcbench-benchmarks


text:
filename='Bridge_triangle_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function scan_data = rdprfile(ufilename,plotmode,verbose)
%RDPRFILE Reads data from a profilometer data file
%
% DATA = RDPRFILE(UFILENAME,VERBOSE)
%
% RDPRFILE examines the header in a profilometer data file and attempts to
%  determine the type of data file. If a supported type is found, it then
%  uses the appropriate function to open the file and return the data.
%  Returns SCAN_DATA=[] if unable to determine the file type.
%
% Supported profilometers:
%  Dektak IIA
%  AlphaStep 200
%  Alphastep 500
%  Dektak 3ST
%  Veeco Dektak 8
%  Bruker Dektak XT
%  Tencor P-10
%
% Options:
% If VERBOSE == 0, no messages will be displayed.
% If VERBOSE == 1, the data summary will be displayed. (default)
% If VERBOSE == 2, the file header and data summary will be displayed.
% If PLOTMODE == 0, the data will not be displayed.
% if PLOTMODE == 1, the data will be displayed in a new figure. (default)
%
% See also dk2rd.m, as2rd.m, as5rd.m, dk3rd.m, dk8rd.m, dkXTrd.m, p10rd.m
%
% M.A. Hopcroft
%  mhopeng at gmail dot com
%
% The license governing use of this code is in the accompanying file
%  "license.txt".
%
versionstr='rdprfile 1.1';

%
% MH Jul2013
% v1.1  add Dektak XT
% MH Sep2012
% v1.0  works

% handle input arguments - set defaults
% default to verbose mode 1 (summary only)
if nargin < 3
    verbose = 1;
end

% default to plot mode 1
if nargin < 2
    plotmode = 1;
end

% get interactive if not being called by another program
if (nargout == 0)
    if (verbose >= 1 && nargin == 0)
        clc;
        fprintf(1, '\n\n  %s\n  %s\n\n', versionstr, 'Load profilometer data files');
    end
end

%%%%
%% open user's file
if nargin >= 1 % if the user already specified a file, use it
    userfilename=ufilename;
    % check to see if file exists
    cfid=fopen(userfilename);
    if cfid== -1
        if cfid== -1
            disp(pwd);
            fprintf(1,'ERROR: data file ["%s"] not found in current directory',userfilename);
            scan_data=[];
            return
        end
        fclose(cfid);
    end

else % if user did not specify the file name, ask for one
    [ufilename, ufilepath]=uigetfile({'.txt','.csv'},'Select the Profilometer data file');
    if isequal(ufilename,0)
        scan_data=[];
        return
    end
    userfilename=fullfile(ufilepath,ufilename);
    %userfilename=input(' Enter the name of the profilometer data file: ','s');
end

% open the file
userfile=fopen(userfilename);

% % %
%% read lines from the header until we have an answer
if (nargout == 0), fprintf(1,'\n '); end % screen formatting for interactive mode

%% read the first line in the file
linein1=fgetl(userfile);
while isempty(linein1), linein1=fgetl(userfile); end % handle spurious carriage returns
if verbose >= 2, fprintf(1,'%s\n',linein1); end

% is this a Dektak II file?
if strfind(linein1,'SLOAN DEKTAK II')
    fprintf(1,'rdprfile: This appears to be a data file from a Dektak IIA\n');
    fclose(userfile);
    scan_data=dk2rd(userfilename,plotmode,verbose);
    return
end

% is this a Dektak XT file?
if strfind(linein1,'Meta Data')
    fprintf(1,'rdprfile: This appears to be a data file from a Dektak XT\n');
    fclose(userfile);
    scan_data=dkXTrd(userfilename,plotmode,verbose);
    return
end

%% read another line
linein2=fgetl(userfile);
while isempty(linein2), linein2=fgetl(userfile); end % handle spurious carriage returns
if verbose >= 2, fprintf(1,'%s\n',linein2); end

% is this a AlphaStep 200 file?
if length(linein2)>=4 && strcmp(linein2(1:4),'VERT')
    fprintf(1,'rdprfile: This appears to be a data file from an AlphaStep 200\n');
    fclose(userfile);
    scan_data=as2rd(userfilename,plotmode,verbose);
    return
end

% is this a AlphaStep 500 file?
dbl1=strtrim(linein1);
dbl2=strtrim(linein2);
if strcmp(dbl1(1),'A') && strcmp(dbl2(1),'B')
    fprintf(1,'rdprfile: This appears to be a data file from an AlphaStep 500\n');
    fclose(userfile);
    scan_data=as5rd(userfilename,plotmode,verbose);
    return
end

% is this a Tencor P-10 file?
if strcmp(linein1(1),'*') && strfind(linein2,'data points')
    fprintf(1,'rdprfile: This appears to be a data file from a Tencor P-10\n');
    fclose(userfile);
    scan_data=p10rd(userfilename,plotmode,verbose);
    return
end

% is this a Dektak 3ST file?
if length(linein1)>=9 && strcmp(linein1(1:9),'DATA FILE') && length(linein2)>=9 && strcmp(linein2(1:9),'AUTO PROG')
    fprintf(1,'rdprfile: This appears to be a data file from a Dektak 3ST\n');
    fclose(userfile);
    scan_data=dk3rd(userfilename,plotmode,verbose);
    return
end

%% read more lines
linein3=fgetl(userfile);
while isempty(linein3), linein3=fgetl(userfile); end
if verbose >= 2, fprintf(1,'%s\n',linein3); end
linein4=fgetl(userfile);
while isempty(linein4), linein4=fgetl(userfile); end
if verbose >= 2, fprintf(1,'%s\n',linein4); end
linein5=fgetl(userfile);
while isempty(linein5), linein5=fgetl(userfile); end
if verbose >= 2, fprintf(1,'%s\n',linein5); end

% is this a Dektak 8 file?
if length(linein4)>=9 && strcmp(linein4(1:9),'DATA FILE') && length(linein5)>=4 && strcmpi(linein5(1:4),'Date')
    fprintf(1,'rdprfile: This appears to be a data file from an Dektak 8\n');
    fclose(userfile);
    scan_data=dk8rd(userfilename,plotmode,verbose);
    return
end

% fail
fprintf(1,'rdprfile: Unable to determine file type!\n');
scan_data=[];


text:
function Err = errornormalize(Err)
%ERRORNORMALIZE  Normalization of error term
%
%For internal use only

%
%  Err = errornormalize(Err)
%
%with Err = Err.mant * beta^Err.exp and Err.mant normalized to
%approximately 1.
%

% written  12/30/98     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_LONG_BETA = getappdata(0,'INTLAB_LONG_BETA');

  % adapt error mantissa and exponent such that Err.mant around 1
  index = ( Err.mant~=0 );
  if any(index)
    q = floor( log( Err.mant(index) ) / log( INTLAB_LONG_BETA ) );
    Err.mant(index) = Err.mant(index) ./ INTLAB_LONG_BETA.^q;
    Err.exp(index) = Err.exp(index) + q;
  end
  Err.exp(~index) = 0;

repo: OpenMAS


text:
function res = parallelLine3d(line, point)
%PARALLELLINE3D  Create 3D line parallel to another one.
%
%   L2 = parallelLine3d(L, P)
%   Creates the 3D line L2, parallel to the line L, and containing the
%   point P.
%
%   Example
%
%   See also
%   geom3d, parallelLine, parallelPlane

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2012-08-23, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

res = [point line(:, 4:6)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function CreateViolaPanel(handles)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

    %Viola Panel
    handles.violaPanel = uipanel('Parent',handles.toolPanel,'Units','normalize','FontSize',10,'Title','Viola Optimizer Setup',...
                                    'Tag','violaPanel','Clipping','off',...
                                    'Position',[0 0 1 1],...
                                    'Visible','on');
    dx = 0; dy = 0;
%%
    h282 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[9.87096774193548+dx 305+dy 167.562724014337 25.0167597765363],...
    'String','Standard Deviation Fixed:',...
    'Style','text',...
    'Tag','text43');

    handles.viola_sdF = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 310+dy 70.573476702509 20],...
    'String','0.4',...
    'Style','edit',...
    'Tag','viola_sdF');

    h284 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6.91397849462366+dx 285+dy 170.519713261649 20],...
    'String','Standard Deviation Moving:',...
    'Style','text',...
    'Tag','text44');

    handles.viola_sdM = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 285+dy 70.573476702509 20],...
    'String','0.4',...
    'Style','edit',...
    'Tag','viola_sdM');

    h286 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[7.89964157706093+dx 260+dy 169.534050179211 20],...
    'String','Smoother Variance Fixed:',...
    'Style','text',...
    'Tag','text45');

    handles.viola_svF = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 260+dy 70.573476702509 20],...
    'String','2.0',...
    'Style','edit',...
    'Tag','viola_svF');

    h288 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6.91397849462366+dx 235.837988826816+dy 170.519713261649 20],...
    'String','Smoother Variance Moving:',...
    'Style','text',...
    'Tag','text46');

    handles.viola_svM = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 234.798882681564+dy 70.573476702509 20],...
    'String','2.0',...
    'Style','edit',...
    'Tag','viola_svM');
%%

    h290 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[9.87096774193548+dx 197.391061452514+dy 167.562724014337 20],...
    'String','Number of Samples(%):',...
    'Style','text',...
    'Tag','text47');

    h291 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8.88530465949821+dx 67.5027932960894+dy 168.548387096774 20],...
    'String','Maximum Step:',...
    'Style','text',...
    'Tag','text48');

    handles.viola_sn = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 200.508379888268+dy 70.573476702509 20],...
    'String','0.01',...
    'Style','edit',...
    'Tag','viola_sn');

    handles.viola_maxstep = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 69.5810055865922+dy 70.573476702509 20],...
    'String','16.0',...
    'Style','edit',...
    'Tag','viola_maxstep');

    h294 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6.91397849462366+dx 38.4078212290503+dy 170.519713261649 20],...
    'String','Iteration Number :',...
    'Style','text',...
    'Tag','text49');

    handles.viola_in = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 38.4078212290503+dy 70.573476702509 20],...
    'String','50',...
    'Style','edit',...
    'Tag','viola_in');

    h296 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6.91397849462366+dx 7.23463687150838+dy 170.519713261649 20],...
    'String','Registration Levels :',...
    'Style','text',...
    'Tag','text116');

    handles.viola_reglevels = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 7.23463687150838+dy 70.573476702509 20],...
    'String','3',...
    'Style','edit',...
    'Tag','viola_reglevels');

    h298 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8.88530465949821+dx 100.754189944134+dy 168.548387096774 20],...
    'String','Minimum Step:',...
    'Style','text',...
    'Tag','text117');

    handles.viola_minstep = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[187.290322580645+dx 102.832402234637+dy 70.573476702509 20],...
    'String','2.5',...
    'Style','edit',...
    'Tag','viola_minstep');

    h300 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[7.89964157706093+dx 132.966480446927+dy 168.548387096774 20],...
    'String','Translation Scale:',...
    'Style','text',...
    'Tag','text118');

    handles.viola_tscale = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[186.304659498208+dx 135.04469273743+dy 70.573476702509 20],...
    'String','0.001',...
    'Style','edit',...
    'Tag','viola_tscale');

    h302 = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[7.89964157706093+dx 166.217877094972+dy 168.548387096774 20],...
    'String','Rotation Scale:',...
    'Style','text',...
    'Tag','text119');

    handles.viola_rscale = uicontrol(...
    'Parent',handles.violaPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[186.304659498208+dx 168.296089385475+dy 70.573476702509 20],...
    'String','1.0',...
    'Style','edit',...
    'Tag','viola_rscale');

    guidata(handles.mainframe, handles);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function power_loss(y, ons, X)

%  'true' model fit

% assume 'true' is FIR estimate
DX = tor_make_deconv_mtx3(ons, 20, 1);

% get residual variance and rss from best model
[bdx, fdx, rdx, rss, r_var, N, kdx] = fit_model(DX, y);

% if you had specified the right canonical shape at the start
Xbest = conv(ons, bdx); Xbest = Xbest(1:N);
Xbest(:, end+1) = 1;

% get residual variance and rss from best model
[bdx, fdx, rdx, rss, r_var, N, kdx] = fit_model(DX, y);

% E(t), expected t-value if you had gotten the HRF exactly right
% --------------------------------------------------------
c = eye(2);
c(:, end) = []; % remove intercept; assumes intcpt is at end!

[bbest, fbest, rbest, rss, r_var, N, kbest] = fit_model(Xbest, y);

[tbest, pbest] = get_t_values(c, bbest, Xbest, r_var);

delta = tbest;
alpha = .001;
df = N - kbest;
t_thresh = tinv(1 - alpha ./ 2, df);

power_best = 1 - nctcdf(t_thresh, df, delta);

% E(t), expected t-value with canonical HRF
% --------------------------------------------------------
hrf = spm_hrf(2);
CX = conv(ons, hrf); CX = CX(1:N);
CX(:, end+1) = 1;

% get residual variance and rss from canonical model
[bc, fc, rc, rss_c, r_var_c, N, kc] = fit_model(CX, y);

[tc, pc] = get_t_values(c, bc, CX, r_var_c);

delta = tc;
df = N - kc;
power_c = 1 - nctcdf(t_thresh, df, delta);

% Difference: Power loss
power_loss_val = 100 * (power_best - power_c) ./ power_best;

% significance for FIR model: full vs. reduced
% --------------------------------------------------------
Xred = ones(N, 1);

[F, p, dfb, dfe] = full_vs_reduced(y, DX, Xred, rss, N, kdx);

f_thresh = finv(1 - alpha, dfb, dfe);
fdelta = F;
power_fir = ncfcdf(f_thresh, dfb, dfe, fdelta);

power_loss_fir = 100 * (power_best - power_fir) ./ power_best;

% Plot
% --------------------------------------------------------
create_figure('plot'); plot(y, 'k'); hold on;
plot(fbest, 'r');
plot(fc, 'b');

fprintf(1,'With best (ideal HRF) model: t = %3.2f, p = %3.4f, power = %3.4f\n', tbest, pbest, power_best);
fprintf(1,'With canonical HRF model: t = %3.2f, p = %3.4f, power = %3.4f\n', tc, pc, power_c);
fprintf(1,'With FIR model F-test: F = %3.2f, p = %3.4f, power = %3.4f\n', F, p, power_fir);

fprintf(1, 'Power loss: Canonical vs. ideal: %3.2f%%\n', power_loss_val);
fprintf(1, 'Power loss: FIR vs. ideal: %3.2f%%\n', power_loss_fir);

end

%
% full vs. reduced
function [F, p, dfb, dfe] = full_vs_reduced(y, Xfull, Xred, rss, N, kfull)

[bred, fred, rred, rss_red] = fit_model(Xred, y);

dfb = (kfull - 1);
dfe = (N - kfull);

F = ( (rss_red - rss) ./ dfb ) ./ ( rss ./ dfe );

p = 1 - fcdf(F, dfb, dfe);

end

function [b, f, r, rss, r_var, N, k] = fit_model(X, y)
N = size(y, 1);
k = size(X, 2);

b = pinv(X) * y;
f = X * b;
r = y - f;

rss = r' * r;
r_var = rss ./ (N - k);

end

function [t, p, cb, var_cbeta] = get_t_values(c, b, X, r_var)

var_cbeta = diag( sqrt(r_var) * c' * inv(X' * X) * c );

% effect
cb = c' * b;

% t-values
t = cb ./ var_cbeta;

p = 2 * ( 1 - tcdf(abs(t), size(X,1) - size(X,2)) );

end

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function L = ivmLikelihoods(model, x, y);

% IVMLOGLIKELIHOODS Return the likelihood for each point for the IVM.
% FORMAT
% DESC computes the likelihood of each given point under the given
% IVM model structure.
% ARG model : the IVM for which the likelihoods are to be computed.
% ARG x : the input points where the likelihoods are to be
% computed.
% ARG y : the target points where the likelihoods are to be
% computed.
%
% SEEALSO : noiseLikelihood, ivmPosteriorMeanVar
%
% COPYRIGHT : Neil D. Lawrence, 2005

if nargin < 3
  % This implies evaluate for the traing data.
  mu = model.mu;
  varsigma = model.varSigma;
  y = model.y;
else
  [mu, varsigma] = ivmPosteriorMeanVar(model, x);
end

L = noiseLikelihood(model.noise, mu, varsigma, y);


text:
clc;
clear;
close all;
warning off;
addpath(genpath(cd));

I=double(imread('E:\lichang\1.Image fusion total variation\img\tank\tank_IR.png'))/255;
%V=double(imread('E:\lichang\1.Image fusion total variation\img\tank\tank_VIS.png'))/255;
% I=histeq(I);
% V=histeq(V);
% figure; imshow(I);
% figure; imshow(V);
% imwrite(I,'new_4917_IR.png','png');
% imwrite(V,'new_4917_VIS.png','png');
% I=rgb2gray(I);
% V=rgb2gray(V);
%proposed
nmpdef;
pars_irn = irntvInputPars('l2tv');

pars_irn.adapt_epsR   = 1;
pars_irn.epsR_cutoff  = 0.01;   % This is the percentage cutoff
pars_irn.adapt_epsF   = 1;
pars_irn.epsF_cutoff  = 0.05;   % This is the percentage cutoff

pars_irn.pcgtol_ini = 1e-4;

pars_irn.loops      = 5;
pars_irn.U0         = I;

pars_irn.variant       = NMP_TV_SUBSTITUTION;
pars_irn.weight_scheme = NMP_WEIGHTS_THRESHOLD;
pars_irn.pcgtol_ini    = 1e-2;
pars_irn.adaptPCGtol   = 1;
Final_Metric=zeros(6,1);
tic;
for i=[0.1:0.1:0.9]%[0.01:0.01:0.09 0.1:0.1:0.9 1:10 20:10:100]

    X = irntv(I, {}, i, pars_irn);
    X=im2gray(X);
    imwrite(X,strcat('Proposed_tank_',strcat(num2str(i),'.png')),'png');
    %figure;
    %imshow(X);
%     Result = Metric(uint8(abs(I)*255),uint8(abs(V)*255),uint8(abs(X*255)));
%     temp=Result.Total;
%     Final_Metric=max(Final_Metric,temp);
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [data, Ht] = bekk_simulate(T,k,parameters,p,o,q,type)
% Simulation of symmetric and asymmetric BEKK(p,o,q) multivariate volatility models
%
% USAGE:
%  [DATA,HT] = bekk_simulate(T,K,PARAMETERS,P,O,Q,TYPE)
%
% INPUTS:
%   T          - Either a scalar containing the length of the series to simulate, or a T by K matrix
%                  of simulated random variables.  The default is to use standard normals.  Providing
%                  a T by K matrix allows other distributions to be used.
%   PARAMETERS - Vector of parameters.  The form of the parameters depends on the TYPE.
%                  'Scalar':
%                  [CC' a(1) ... a(p) g(1) ... g(o) b(1) ... b(q)]'  (all scalars)
%                  'Diagonal'
%                  [CC' diag(A(:,:,1))' ... diag(A(:,:,p))' diag(G(:,:,1))' ... diag(G(:,:,o))' diag(B(:,:,1))' ... diag(B(:,:,p))']'
%                  'Full'
%                  [CC' f(A(:,:,1)) ... f(A(:,:,p)) f(G(:,:,1)) ... f(G(:,:,o)) f(B(:,:,1)) ... f(B(:,:,q))]'
%                  where CC = chol2vec(C')' and f(M) = M(:)'
%   P          - Positive, scalar integer representing the number of symmetric innovations
%   O          - Non-negative, scalar integer representing the number of asymmetric innovations
%   Q          - Non-negative, scalar integer representing the number of conditional covariance lags
%   TYPE       - String, one of :
%                  'Scalar' (Default)
%                  'Diagonal'
%                  'Full'
%
% OUTPUTS:
%   DATA   - A T by K matrix of simulated data
%   HT     - A [K K T] dimension matrix of conditional covariances
%
% COMMENTS:
%   The dynamics of a BEKK are given by
%
%   H(:,:,t) = C*C' +
%       A(:,:,1)'*OP(:,:,t-1)*A(:,:,1) + ... + A(:,:,p)'*OP(:,:,t-1)*A(:,:,p) +
%       G(:,:,1)'*OPA(:,:,t-1)*G(:,:,1) + ... + G(:,:,o)'*OPA(:,:,t-1)*G(:,:,o) +
%       B(:,:,1)'*G(:,:,t-1)*B(:,:,1) + ... + B(:,:,q)'*OP(:,:,t-1)*B(:,:,q)
%
%   where in the scalar model A(:,:,i) = a(i)*eye(K) (similarly for G and B).
%
%  EXAMPLES:
%    % Scalar with A.^2=.05, G.^2=.1 and B.^2=.88
%    CCp = [1 .5;.5 4];
%    parameters = [chol2vec(chol(CCp)');sqrt([.05,.10,.88])']
%    [data,Ht] = bekk_simulate(1000,2,parameters,1,1,1,'Scalar')
%    % Diagonal
%    parameters = [chol2vec(chol(CCp)');sqrt([.05 .07 .93 .88])']
%    [data,Ht] = bekk_simulate(1000,2,parameters,1,0,1,'Diagonal')
%
% See also BEKK

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isscalar(T)
    e = randn(2*T,k);
else
    e = T;
    if size(e,2)~=k
        error('T must have K columns when providing simulated random numbers.')
    end
    T = size(e,1);
    e = [e(ceil(rand(T,1)*T),:);e];
end

if strcmpi(type,'Scalar')
    type = 1;
elseif strcmpi(type,'Diagonal')
    type = 2;
elseif strcmpi(type,'Full')
    type = 3;
else
    error('TYPE must be ''Scalar'', ''Diagonal'' or ''Full''.')
end

if p<=0 || floor(p)~=p
    error('P must be a positive scalar integer.')
end
if o<0 || floor(o)~=o
    error('O must be a positive scalar integer.')
end
if q<0 || floor(q)~=q
    error('Q must be a positive scalar integer.')
end

k2 = k*(k+1)/2;
switch type
    case 1
        count = p+o+q;
    case 2
        count = (p+o+q)*k;
    case 3
        count = (p+o+q)*k*k;
end
count = count + k2;
if length(parameters)~=count
    error('PARAMETERS does not have the expected number of elements.')
end
[C,A,G,B] = bekk_parameter_transform(parameters,p,o,q,k,type);
m = zeros(k*k);
for i=1:p
    m = m + kron(A(:,:,i),A(:,:,i));
end
for i=1:o
    m = m + 0.5*kron(G(:,:,i),G(:,:,i));
end
for i=1:q
    m = m + kron(B(:,:,i),B(:,:,i));
end

if max(eigs(m))>1
    backCast = C/.001;
    warning('MFE:nonstationary','The parameters do not correspond to the stationary region.')
else
    backCast = ((eye(k*k)-m)\eye(k*k))*C(:);
    backCast = reshape(backCast,k,k);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ht = repmat(eye(k),[1 1 2*T]);
eta = zeros(size(e));
for i=1:(2*T)
    Ht(:,:,i) = C;
    for j=1:p
        if (i-j)<=0
            Ht(:,:,i) = Ht(:,:,i) + A(:,:,j)'*backCast*A(:,:,j);
        else
            Ht(:,:,i) = Ht(:,:,i) + A(:,:,j)'*(e(i-j,:)'*e(i-j,:))*A(:,:,j);
        end
    end
    for j=1:o
        if (i-j)<=0
            Ht(:,:,i) = Ht(:,:,i) + G(:,:,j)'*backCast*G(:,:,j);
        else
            Ht(:,:,i) = Ht(:,:,i) + G(:,:,j)'*(eta(i-j,:)'*eta(i-j,:))*G(:,:,j);
        end
    end
    for j=1:q
        if (i-j)<=0
            Ht(:,:,i) = Ht(:,:,i) + B(:,:,j)'*backCast*B(:,:,j);
        else
            Ht(:,:,i) = Ht(:,:,i) + B(:,:,j)'*Ht(:,:,i-j)*B(:,:,j);
        end
    end
    Ht12 = Ht(:,:,i)^(0.5);
    e(i,:) = e(i,:)*Ht12;
    eta(i,:) = e(i,:).*(e(i,:)<0);
end

data = e(T+1:2*T,:);
Ht = Ht(:,:,T+1:2*T);

author: bashtage

repo: mfe-toolbox


text:
function LU = extract_bounds_from_norm_operator(LU,extstruct,extvariables,i);
arg = extstruct(i).arg{1};
if min(size(arg))==1 & length(extstruct(i).arg)>1
    p = extstruct(i).arg{2};
    % Maximize each element conservatively, and use the norm of
    % that worst-case vector
    xmax = abs(getbase(arg))*[1;max(abs(LU(getvariables(arg),:)),[],2)];
    LU(extvariables(i),2) = min([norm(xmax,p) LU(extvariables(i),2)]);
end
% norms are positive...
LU(extvariables(i),1) = max([0 LU(extvariables(i),1)]);

author: shamilmamedov


text:
function lagplot(c);

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% PREP PLOT
figure('Color','w','Position',[50 50 600 500]);
set(gcf,'DefaultAxesFontSize',14);
imagesc(c.L);
title('Lag time for maximum correlation (s)');

% ADD DATES TO AXES
n = length(c.trig);
set(gca,'XTick',[1:round(n/25):n]);
set(gca,'YTick',[1:round(n/25):n]);
yt = get(gca,'YTick');
set(gca,'YTickLabel',datestr(c.trig(yt),'yyyy-mm-dd HH:MM'),'FontSize',6);

% DRESS UP THE FIGURE
cmap = load('colormap_lag.txt');
colormap(cmap);
colorbar;
xlabel('Event number');
ylabel('Event date');

%PRINT OUT FIGURE
set(gcf, 'paperorientation', 'portrait');
set(gcf, 'paperposition', [1.25 2.5 6 6] );
%print(gcf, '-depsc2', 'FIG_tartan.ps')

author: geoscience-community-codes

repo: GISMO


text:
function pdvparamsgui(varargin)
%PDVPARAMSGUI Allows one to set PDV parameters through a GUI
%
% Written by: Tim Cox, NRL, and Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Parse input arguments
p = inputParser; % Extract parameter-value pairs
p.KeepUnmatched=true;
p.addParamValue('filename','');
p.addParamValue('aoi',[0 0 0 0]);
p.addParamValue('segment',1);
p.parse(varargin{:});

%% Draw figure
fig_hand = figure('MenuBar','none','Toolbar','none','NumberTitle','off',...
    'Name','PDV','Units','Pixels','Position',[0 0 300 250]);
movegui(fig_hand,'center');

% Pixel shift
uicontrol('Style','Text','String','Pixel Shift:',...
    'Units','normalized','Position',[.05 .7 .2 .1],'parent',fig_hand,...
    'BackgroundColor',get(fig_hand,'Color'));
edit_ps = uicontrol('Style','Edit','String','0.25',...
    'Units','normalized','Position',[.25 .7 .15 .11],'parent',fig_hand);

% Direction
bg_dir = uibuttongroup('Units','normalized','Position',[.05 .25 .35 .35],...
    'Title','Direction','BackgroundColor',get(fig_hand,'Color'));
radio_st = uicontrol('Style','Radio','String','Slow-Time',...
    'Units','normalized','Position',[.1 .65 .9 .2],'parent',bg_dir,...
    'BackgroundColor',get(fig_hand,'Color'));
radio_ft = uicontrol('Style','Radio','String','Fast-Time',...
    'Units','normalized','Position',[.1 .25 .9 .2],'parent',bg_dir,...
    'BackgroundColor',get(fig_hand,'Color'));

% Window size
uip_smooth = uipanel('Units','normalized','Position',[.45 .6 .5 .35],'parent',fig_hand,...
    'Title','Smoothing Window','BackgroundColor',get(fig_hand,'Color'));
uicontrol('Style','Text','String','Azimuth:',...
    'Units','normalized','Position',[.05 .6 .4 .3],'parent',uip_smooth,...
    'BackgroundColor',get(fig_hand,'Color'));
edit_az = uicontrol('Style','Edit','String','5',...
    'Units','normalized','Position',[.45 .6 .2 .35],'parent',uip_smooth);
uicontrol('Style','Text','String','Pixels',...
    'Units','normalized','Position',[.7 .6 .25 .3],'parent',uip_smooth,...
    'BackgroundColor',get(fig_hand,'Color'));
uicontrol('Style','Text','String','Range:',...
    'Units','normalized','Position',[.05 .1 .4 .3],'parent',uip_smooth,...
    'BackgroundColor',get(fig_hand,'Color'));
edit_rng = uicontrol('Style','Edit','String','5',...
    'Units','normalized','Position',[.45 .1 .2 .35],'parent',uip_smooth);
uicontrol('Style','Text','String','Pixels',...
    'Units','normalized','Position',[.7 .1 .25 .3],'parent',uip_smooth,...
    'BackgroundColor',get(fig_hand,'Color'));

% Fitler type
bg_ft = uibuttongroup('Units','normalized','Position',[.45 .25 .5 .35],...
    'Title','Filter Type','BackgroundColor',get(fig_hand,'Color'));
radio_mean = uicontrol('Style','Radio','String','Mean',...
    'Units','normalized','Position',[.1 .65 .9 .2],'parent',bg_ft,...
    'BackgroundColor',get(fig_hand,'Color'));
radio_med = uicontrol('Style','Radio','String','Median',...
    'Units','normalized','Position',[.1 .25 .9 .2],'parent',bg_ft,...
    'BackgroundColor',get(fig_hand,'Color'));
uicontrol('Style','pushbutton','String','Run','FontWeight','bold',...
    'Units','normalized','Position',[.35 .05 .3 .15],'parent',fig_hand,...
    'Callback',@Run_Callback);

    %% Callback for button press
    function Run_Callback(hObject, eventdata)
        if get(radio_st,'Val')
            Dim = 1;
        else
            Dim = 2;
        end

        if get(radio_mean,'Value')
            filtertype = 'mean';
        else
            filtertype = 'median';
        end

        if ~isempty(p.Results.filename)
            %now generate PDV for specified file and AOI
            if (p.Results.aoi(1) == 0)
                AzLimits = 'full';
                RgLimits = 'full';
            else
                AzLimits = round([p.Results.aoi(1)  p.Results.aoi(1)+p.Results.aoi(3)-1]);
                RgLimits = round([p.Results.aoi(2)  p.Results.aoi(2)+p.Results.aoi(4)-1]);
            end
            tmp_name=tempname;
            tmp_name2=tempname;

            pdvfile(p.Results.filename,tmp_name,...
                'azlimits',AzLimits,'rnglimits',RgLimits,...
                'framenumber',p.Results.segment,...
                'deltax',str2double(get(edit_ps,'String')),...
                'dim',Dim,...
                'filtersize',[str2double(get(edit_az,'String')) str2double(get(edit_rng,'String'))],...
                'filtertype',filtertype);
            chipfile(p.Results.filename,tmp_name2,...
                'azlimits',AzLimits,'rnglimits',RgLimits,...
                'framenumber',p.Results.segment);

            %launch MITM viewer to view PDV
            mitm_viewer({tmp_name tmp_name2},'initialRemap','linearremap');
            set(gcf,'DeleteFcn',@file_cleanup); % Delete temp files on closing
        else
            pdvdemo('deltax',str2double(get(edit_ps,'String')),...
                'dim',Dim,...
                'filtersize',[str2double(get(edit_az,'String')) str2double(get(edit_rng,'String'))],...
                'filtertype',filtertype);
        end
        close(fig_hand);

        function file_cleanup(obj, eventdata)
            delete(get(obj,'Children')); % Calls destructor for hg_mitm_viewer, which releases files
            % Clean up intermediate files
            delete(tmp_name);
            delete(tmp_name2);
        end
    end
end

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

author: ngageoint


text:
function plotb2east
% set the default plot direction of the x-axis

setMTEXpref('bAxisDirection','east');
setMTEXpref('aAxisDirection','');
setMTEXpref('aStarAxisDirection','');

author: mtex-toolbox

repo: mtex

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


text:
%computes the spectral flatness from the magnitude spectrum
%> called by ::ComputeFeature
%>
%> @param X: spectrogram (dimension FFTLength X Observations)
%> @param f_s: sample rate of audio data (unused)
%>
%> @retval vtf spectral flatness
% ======================================================================
function [vtf] = FeatureSpectralFlatness (X, f_s)

    % replace product sum with log/exp
    XLog = log(X + 1e-20);
    vtf = exp(mean(XLog, 1)) ./ (mean(X, 1));

    % avoid NaN for silence frames
    vtf (min(X, [], 1) == 0) = 0;

end

author: alexanderlerch

repo: ACA-Code


text:
function filterSteerable( theta )
% Steerable 2D Gaussian derivative filter (for visualization).
%
% This function is a demonstration of steerable filters.  The directional
% derivative of G in an arbitrary direction theta can be found by taking a
% linear combination of the directional derivatives dxG and dyG.
%
% USAGE
%  filterSteerable( theta )
%
% INPUTS
%  theta   - orientation in radians
%
% OUTPUTS
%
% EXAMPLE
%  filterSteerable( pi/4 );
%
% See also filterGauss
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% Get G
[x,y]=meshgrid(-1:.1:1, -1:.1:1 );
r = sqrt( x.^2 + y.^2 );
G = exp( -r .* r *2 );

% get first derivatives of G.  note: d/dx(G)=-2x*G
phi = atan2( y, x );
dxG = r .* cos(phi) .* G;
dyG = r .* sin(phi) .* G;

% get directional derivative by taking linear comb in theta
Gtheta = cos(theta)*dxG + sin(theta)*dyG;

% dislpay (scale for visualization purposes)
GS = cat(3,G,dxG*2,dyG*2,Gtheta*2);
figure(1); montage2( GS );

author: zouchuhang

repo: LayoutNet


text:
function item = subsasgn(item, subs, varargin)

% function item = subsasgn(item, subs, varargin)
% This function implements subsasgn for all classes derived from cfg_item.
% It relies on the capability of each class constructor to re-classify a
% struct object after a new value has been assigned to its underlying
% struct (This capability has to be implemented in the derived class).
% The structure of a configuration tree does not permit any arrays of
% cfg_item objects. Therefore, the only subscript reference and
% assignment within an cfg_item is a dot assignment to fields of this
% cfg_item.
% Subscript references we have to deal with are:
% one level
% item.(field)   - i.e. struct('type',{'.'} ,'subs',{field})
%
% to be dealt with elsewhere
% item.(field){fidx}
%
% In a future version, '()' and '{}' subscripts may be supported to
% access val fields of a cfg_item tree as if they were part of a
% harvested job. For cfg_branch objects (where dot assignments are used
% for val fields in their job tree) it is mandatory to index the job as a
% struct array to access harvested fields.
% This function is identical for all classes derived from cfg_item. A
% copy of it must be present in each derived class to be able to access
% derived fields.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: subsasgn.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

persistent local_mysubs_fields;
persistent par_class;
persistent par_fields;
if ~iscell(local_mysubs_fields)
    local_mysubs_fields = mysubs_fields;
    citem = class(item);
    switch citem
        case 'cfg_exbranch',
            par_class = 'cfg_branch';
            pf1 = subs_fields(item.cfg_branch);
            pf2 = subs_fields(cfg_item);
            par_fields = [pf1(:); pf2(:)]';
        case 'cfg_item',
            par_class = '';
            par_fields = {};
        otherwise
            par_class = 'cfg_item';
            par_fields = subs_fields(item.cfg_item);
    end;
end
if numel(item) ~= 1
    cfg_message('matlabbatch:subsasgn', ...
          'Arrays of cfg_item objects not supported.');
end;

switch subs(1).type,
    case {'.'},
        if numel(item) == 1 && nargin == 3
            if numel(subs) > 1
                % Only part of field value(s) assigned. Get old field value(s),
                % assign values to it.
                if any(strcmp(subs(1).subs,par_fields))
                    val0 = item.(par_class).(subs(1).subs);
                else
                    val0 = item.(subs(1).subs);
                end
                val1 = cfg_callbuiltin('subsasgn',val0,subs(2:end),varargin{1});
            else
                % Field values to assign
                val1 = varargin{1};
            end
            [ok, val] = subsasgn_check(item,subs,val1);
            if ok,
                switch subs(1).subs
                    case par_fields,
                        item.(par_class).(subs(1).subs) = val;
                    case local_mysubs_fields,
                        item.(subs(1).subs) = val;
                    otherwise
                        cfg_message('matlabbatch:subsasgn', ...
                                    ['Reference to unknown field ''%s''.\n' ...
                                     'To assign to a field in the job structure, use a reference like ' ...
                                     '''(x).%s''.'], subs(1).subs, subs(1).subs);
                end;
            end;
        else
            cfg_message('matlabbatch:subsasgn', ...
                        ['Array assignments not supported for fields of cfg_item objects.\n' ...
                         'To assign to a field in the job structure, use a reference like ' ...
                         '''(x).%s''.'], subs(1).subs);

        end;
    case {'()','{}'},
        cfg_message('matlabbatch:subsasgn', ...
                    'Subscript type ''%s'' reserved for future use.', subs(1).type);
    otherwise
        cfg_message('matlabbatch:subsasgn', ...
                    'Unknown subsref type: ''%s''. This should not happen.', subs(1).type);
end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [compiler,options] = mexcompiler
% mexcompiler returns the name of the compiler used by mex.

% Written by Tom Minka

compiler = '';
options = struct;
options.vcvarsopts = '';
if ispc
	mexopts = fullfile(prefdir,'mexopts.bat');
	if ~exist(mexopts,'file')
		return
	end
	fid = fopen(mexopts);
	while 1
		txt = fgetl(fid);
		if ~ischar(txt), break, end
		if isempty(compiler) && strncmp('rem ',txt,4)
			compiler = txt(5:end-8);
		elseif ~isempty(strmatch('set ',txt))
			txt = txt(5:end);
			pos = strfind(txt,'=');
			if ~isempty(pos)
				pos = pos(1);
				field = txt(1:(pos-1));
				txt = txt((pos+1):end);
				txt = strrep(txt,'%MATLAB%',matlabroot);
				options.(field) = txt;
				if strcmp(field,'VSINSTALLDIR')
					vsinstalldir = txt;
					if ~isempty(strmatch(':\Program Files (x86)\',vsinstalldir(2:end)))
						% http://msdn.microsoft.com/en-us/library/x4d2c09s(VS.80).aspx
						if exist(fullfile(vsinstalldir,'VC','bin','x86_amd64','cl.exe'))
							options.vcvarsopts = 'x86_amd64';
						else
							options.vcvarsopts = 'x86';
						end
					end
				end
			end
		end
	end
	fclose(fid);
else
	mexopts = fullfile(prefdir,'mexopts.sh');
	if ~exist(mexopts,'file')
		return
	end
	fid = fopen(mexopts);
	while 1
		txt = fgetl(fid);
		if ~ischar(txt), break, end
		if isempty(compiler) && strncmp('rem ',txt,4)
			compiler = txt(5:end-8);
		else
			pos = strfind(txt,'=');
			if ~isempty(pos)
				pos = pos(1);
				field = txt(1:(pos-1));
				txt = txt((pos+1):end);
				txt = strrep(txt,'%MATLAB%',matlabroot);
				options.(field) = txt;
			end
		end
	end
	fclose(fid);
end

author: andrewssobral

repo: lrslibrary


text:
function geometry_test0327 ( )

%*****************************************************************************80
%
%% TEST0327 tests LINE_EXP_NORMAL_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0327\n' );
  fprintf ( 1, '  LINE_EXP_NORMAL_2D determines a unit normal vector\n' );
  fprintf ( 1, '  to a given explicit line.\n' );

  p1(1:2,1) = [ 1.0; 3.0 ];
  p2(1:2,1) = [ 4.0; 0.0 ];

  r8vec_print ( 2, p1, '  Point 1: ' );
  r8vec_print ( 2, p2, '  Point 2: ' );

  normal = line_exp_normal_2d ( p1, p2 );

  r8vec_print ( 2, normal, '  Normal vector N:' );

  return
end


text:
function [rho, stationary]=inverse_ar_roots(phi)
% Computes the inverted roots of the characteristic equation of an AR(P)
%
% USAGE:
%  [RHO, STATIONARY] = inverse_ar_roots(PARAMETERS)
%
% INPUTS:
%   PARAMETERS - A column vector containing the AR parameters in the order
%                  y(t) = mu + phi(1) y(t-1) + phi(2)y(t-2) + ... + phi(P) y(t-P) + e(t)
%
% OUTPUTS:
%   RHO        - A max(P) by 1 vector containing the inverted roots of the characteristic equation
%                  corresponding to the AR model input
%   STATIONARY - Logical indicating whether the AR process is stationary
%
% COMMENTS:
% Process is stationary if min(abs(rho)) > 1
%
%  See also ARMAROOTS, ACF

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 1/1/2007

%Make sure phi is a row vector;

p=length(phi);
%Make sure phi is a row vector
if size(phi,1)>=size(phi,2) && min(size(phi))==1
    phi=phi';
else
    error('Phi should be a column vector.')
end
rho=roots([-fliplr(phi) 1]);
stationary=min(abs(rho))>1;


text:
classdef prtUiListDlg < prtUiManagerPanel

    % clear classes
    % close all
    % strs = {'asdf', 'qwer','poiu','lkhj'};
    % obj = prtUiListDlg(strs);
    %
    % if isempty(obj.enteredSelection)
    %     disp('Nothin');
    % else
    %     disp(strs{obj.enteredSelection});
    % end

    properties
        inputStruct
        tableFieldNames = {};
        dataCell = {};

        enteredSelection = [];

        windowSize = [500 300];
        textFontSize = 25;
        tableFontSize = 20;

        titleStr = '';
        messageStr = '
Make A Choice
';

        cancelStr = 'Cancel';
        okStr = 'OK';

        enableMultiSelect = false;

        handleStruct
    end

    properties (SetAccess = 'protected')
        madeThisWindow = false;
    end

    methods
        function self = prtUiListDlg(varargin)

            if nargin == 1
                self.inputStruct = varargin{1};
            elseif ~mod(nargin,2)
                self = prtUtilAssignStringValuePairs(self,varargin{:});
            else
                self.inputStruct = varargin{1};
                self = prtUtilAssignStringValuePairs(self,varargin{2:end});
            end

            if nargin~=0 && ~self.hgIsValid
               self.create()
            end

            init(self);

            waitfor(self.handleStruct.table);

        end
        function create(self)
            ss = get(0,'ScreenSize');
            screenCenter = ss(3:4)/2;

            self.handleStruct.figureHandle = figure('units','pixels',...
                'position',[screenCenter(1)-self.windowSize(1)/2 screenCenter(2)-self.windowSize(2)/2 self.windowSize],...
                'menubar','none',...
                'toolbar','none',...
                'numberTitle','off',...
                'Name',self.titleStr,...
                'Interruptible','off',...
                'BusyAction','cancel',...
                'DockControls','off');

            self.managedHandle = uipanel(self.handleStruct.figureHandle, ...
                'units','normalized',...
                'BorderType','none',...
                'Position',[0 0 1 1]);

            self.madeThisWindow = true;

        end

        function init(self)

            self.handleStruct.table = uitable('parent',self.managedHandle,...
                'units','pixels','position',[1 1 50 50]); % Dummy position
            if ~isfield(self.handleStruct,'figureHandle') || ~ishandle(self.handleStruct.figureHandle)
                self.handleStruct.figureHandle = gcf; %Is this ok? It depens how we use this.
            end
            drawnow;

            self.handleStruct.jScrollPane = findjobj(self.handleStruct.table);
            self.handleStruct.jTable = self.handleStruct.jScrollPane.getViewport.getView;

            % Set the actual data in the cell
            structArrayToCell(self);
            set(self.handleStruct.table,'data',self.dataCell,...
                'columnName',self.tableFieldNames,...
                'RearrangeableColumns','on',...
                'RowName',{},...
                'TooltipString','Available Options',...
                'ColumnEditable',false(1,length(self.tableFieldNames)),...
                'KeyPressFcn',@self.tableKeyPressFcn,...
                'CellSelectionCallback',@self.cellCelectionFcn);

            % Force selection of full rows only
            self.handleStruct.jTable.setNonContiguousCellSelection(false)

            % Force selection of a single row only
            if ~self.enableMultiSelect
                self.handleStruct.jTable.setSelectionMode(javax.swing.ListSelectionModel.SINGLE_SELECTION);
            end

            % Set the column widths to take up the whole area
            self.handleStruct.jTable.setAutoResizeMode(self.handleStruct.jTable.AUTO_RESIZE_SUBSEQUENT_COLUMNS);

            [self.handleStruct.textHandle, self.handleStruct.textHandleContainer]  = javacomponent('javax.swing.JLabel',[1 1 1 1], self.managedHandle);
            backgroundColor = get(self.managedHandle,'BackgroundColor');
            set(self.handleStruct.textHandle,'Text',self.messageStr,'Background',java.awt.Color(backgroundColor(1),backgroundColor(2),backgroundColor(3)),...
                'Font',java.awt.Font('sansserif',java.awt.Font.PLAIN,self.textFontSize));
            self.handleStruct.textHandle.setVerticalAlignment(javax.swing.JLabel.TOP);

            self.handleStruct.cancelButton = uicontrol('style','pushbutton',...
                'parent',self.managedHandle,...
                'units','pixels',...
                'position',[1 1 10 10],...
                'callback',@self.cancelCallback,...
                'string',self.cancelStr);

            self.handleStruct.okButton = uicontrol('style','pushbutton',...
                'parent',self.managedHandle,...
                'units','pixels',...
                'position',[1 1 10 10],...
                'callback',@self.okCallback,...
                'string',self.okStr);

            % Set and call the resize function to set the element positions
            set(self.managedHandle,'ResizeFcn',@self.resizeFunction);
            resizeFunction(self);

            drawnow;
            set(self.handleStruct.jTable,'Font',java.awt.Font('sansserif',java.awt.Font.PLAIN,self.tableFontSize));
        end

        function close(self)
            if self.madeThisWindow
                try %#ok<TRYNC>
                    close(self.handleStruct.figureHandle);
                end
            else
                try
                    delete(self.handleStruct.table);
                end
            end

        end

        function okCallback(self, varargin)
            self.enteredSelection = getSelection(self);
            close(self)
        end

        function cancelCallback(self, varargin)
            close(self);
        end

        function structArrayToCell(self)
            self.tableFieldNames = fieldnames(self.inputStruct)';
            self.dataCell = struct2cell(self.inputStruct)';
        end

        function resizeFunction(self, varargin)
            panelPosition = getpixelposition(self.managedHandle);

            border = 8;
            textHeight = 50;
            buttonHeight = 20;
            buttonWidth = 50;

            left = border+1;
            width = panelPosition(3);
            bottom = border+1;
            height= panelPosition(4);

            textBottom = max(height-textHeight-border,1);

            okButtonLeft = width/2-border/2-buttonWidth;
            cancelButtonLeft = width/2+border/2;

            tableBottom = bottom+border+buttonHeight;
            tableHeight = max(height-buttonHeight-border*4-textHeight,1);

            okButtonPos = [okButtonLeft bottom buttonWidth, buttonHeight];
            cancelButtonPos = [cancelButtonLeft bottom buttonWidth, buttonHeight];

            tablePos = [left tableBottom width-border*2-2 tableHeight];
            textPos = [left textBottom width-border*2-2 textHeight];

            set(self.handleStruct.okButton,'units','pixels');
            set(self.handleStruct.okButton,'position',okButtonPos);
            set(self.handleStruct.okButton,'units','normalized');

            set(self.handleStruct.cancelButton,'units','pixels');
            set(self.handleStruct.cancelButton,'position',cancelButtonPos);
            set(self.handleStruct.cancelButton,'units','normalized');

            set(self.handleStruct.table,'units','pixels');
            set(self.handleStruct.table,'position',tablePos);
            set(self.handleStruct.table,'units','normalized');

            set(self.handleStruct.textHandleContainer,'units','pixels')
            set(self.handleStruct.textHandleContainer,'position',textPos)
            set(self.handleStruct.textHandleContainer,'units','normalized')

            %drawnow;
        end

        function tableKeyPressFcn(self, varargin)
            event = varargin{2};

            if strcmpi(event.Key,'Return')

                % They hit return, so select this entry and exit
                % When you hit return MATLAB (I think), moves the selection
                % to the "next" cell. If you have the last one selected it
                % goes back around to zero.
                reportedSelection = getSelection(self);
                reportedSelection = reportedSelection - 1; % undo the wrap arround
                if reportedSelection==0 % We actually selected the end and they wrapped us around to 1
                    reportedSelection = length(self.inputStruct);
                end

                self.enteredSelection = reportedSelection;

                drawnow;
                self.handleStruct.jTable.setRowSelectionInterval(reportedSelection-1,reportedSelection-1);
                drawnow;

                close(self);

            elseif strcmpi(event.Key,'Escape')
                % Leave with an empty selection
                close(self)
            else
                % Nada, probably an arrow
            end
        end
        function reportedSelection = getSelection(self)
            if ~self.enableMultiSelect
                reportedSelection = self.handleStruct.jTable.getSelectedRow()+1;% Java is 0 referenced
            else
                reportedSelection = self.handleStruct.jTable.getSelectedRows()+1;% Java is 0 referenced
            end
        end
        function set.inputStruct(self, val)
            assert(isstruct(val) || iscellstr(val),'inputStruct must be a structure array or a cellstr');

            if iscellstr(val)
                cellOfStrs = val; cellOfStrs = cellOfStrs(:);
                self.inputStruct = struct('Choices',cellOfStrs);
            else
                self.inputStruct = val(:);
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
 function ab = BrokenStickRegression(xx, yy, nstick)
%BROKENSTICKREGRESSION  piecewise linear regression. Fits a line
% consisting of connected straight sections to a cloud of data points.
%
% AB = BrokenStickRegression(XX, YY, NSTICK); XX and YY are the data
% points; NSTICK is the number of connected straight sections. AB(:, 1)
% are the x-coordinates of endpoints and breakpoints in ascending order.
% AB(:, 2) are the corresponding y-coordinates. XX need not be in a
% monotonic order.
%
% AB = BrokenStickRegression(XX, YY, BREAKPOINTS); BREAKPOINTS is either
% a vector of at least two abscissa values chosen as starting
% breakpoints or one non-integer starting breakpoint on the abscissa.
% NUMEL(BREAKPOINTS) + 1 is the number of straight sections of the
% fitting curve. Choosing starting breakpoints helps in some cases to
% obtain better fits.
%
% Example 1:
% ---------
%    nstick = 4;
%    nn = 800;
%    xx = linspace(1.0, 11.5, nn)';
%    y0 = sin(xx);
%    yy = y0 + randn(nn, 1) * 0.4;
%    ab = BrokenStickRegression(xx, yy, nstick);
%    plot(xx, yy, 'b.', xx, y0, 'k', ab(:, 1), ab(:, 2), 'r-o')
%    title(['BrokenStickRegression(x, sin(x) + noise, ', ...
%          int2str(nstick), ')'])
%
% Example 2:
% ---------
%   xx = 0:100;
%   yy = [ones(1, 70), 1:10, 11:-1:0];
%   yy = [yy, zeros(1, 101 - length(yy))] + 0.8 * randn(1, 101);
%   bp = [65, 75, 90];                     % breakpoints.
%   plot(xx, yy, '.')
%   hold on
%   ab1 = BrokenStickRegression(xx, yy, numel(bp) + 1);
%   plot(ab1(:, 1), ab1(:, 2), 'k-o')
%   ab2 = BrokenStickRegression(xx, yy, bp);
%   plot(ab2(:, 1), ab2(:, 2), 'r-o')
%   legend('data points', 'NSTICK scalar', 'NSTICK vectorial', ...
%      'Location', 'NW')
%
% See also: POLYFIT.

% The algorithm uses POLYFIT and FMINSEARCH.
%
% pmwnave@yahoo.de
% 2010-11-05, started. Submitted to FEX on 2010-11-13, #29387.
% 2010-11-16, simplified. Update submitted to FEX on 2010-11-15.
% 2010-12-01, case recognized in which polyfit is supplied with only one
%             data point; case recognized in which the first breakpoint
%             slips to the left of the minimum data point. Both bugs
%             were discovered by Carlos Romero, EPFL, Lausanne. Thanks!
%             Update submitted to FEX on 2010-12-05.
% 2010-12-07, penalty scheme simplified.
% 2011-09-04, one breakpoint can now be specified, as suggested by Atul
%             Ingle on 2011-09-01.
%----------------------------------------------------------------------O
% Initialize the coordinates ab of the end- and breakpoints.
%
  xx = xx(:);
  yy = yy(:);
  a0 = [];

  if numel(xx) ~= numel(yy),
     error(' ### %s: XX and YY are not equally long.', mfilename)
  end
  if ~all(isfinite(xx)) || ~all(isfinite(yy)),
     error(' ### %s: XX or YY contain non-finite elements.', mfilename)
  end
  if (nargin == 2) || isempty(nstick),
     nstick = 1;
  end
  if numel(nstick) > 1 || round(nstick) ~= nstick,
     a0 = nstick;
     nstick = numel(a0) + 1;
  end
  if numel(xx) < nstick + 1,
     error(' ### %s: too few data points.', mfilename)
  end

  ab         = zeros(nstick + 1, 2);
  ab(1, 1)   = min(xx);
  ab(end, 1) = max(xx);

  if nstick == 1,
%
% Normal linear regression.
%
     pp       = polyfit(xx, yy, 1);
     ab(:, 2) = pp(1) * ab(:, 1) + pp(2);

  else
%
% Two or more sticks. Find the breakpoints that minimize the residuals.
% The appropriate function would be FMINCON which is part of the
% optimization toolbox, which not everybody owns. A work-around consists
% in applying penalties whenever constraints are violated. The penalties
% are defined in MinRes. A and B define constraints: A solution point X
% is admissible if A * X < B. a0 is(are) the start point(s) of the
% search process.
%
     A  = zeros(nstick, nstick - 1);
     A(1:nstick + 1:end) = -1;
     A(2:nstick + 1:end) =  1;
     B  = zeros(nstick, 1);
     B(1)      = -ab(1, 1);
     B(nstick) =  ab(end, 1);
     if isempty(a0),
        a0 = (ab(end, 1) - ab(1, 1)) / nstick * (1:nstick - 1)' + ...
           ab(1, 1);
     end
     aa = fminsearch(@MinRes, a0);
  end

function rr = MinRes(aa)                 % nested function.
%MinRes  calculates the residuals rr for given breakpoints aa.

%----------------------------------------------------------------------O
  ab(2:nstick, 1) = aa;
%
% Impose a penalty for violating constraining conditions. This approach
% can possibly be improved. Other constraints may be added as penalties.
%
  if nstick > 2 && any(A * ab(2:nstick, 1) >= B),
     rr = 1e20;
     return
  end
%
% Regression on the leftmost section.
%
  kk = find(xx <= ab(2, 1));
  kt = numel(kk);
  if kt == 1,
     tmp = 0;
  else
     [pp, ss, mu] = polyfit(xx(kk), yy(kk), 1);
     ab(1:2, 2)   = polyval(pp, ab(1:2, 1), [], mu);
     tmp          = yy(kk) - polyval(pp, xx(kk), [], mu);
  end
%
% All other sections.
%
  for ii = 2:nstick,
     kk    = find((xx > ab(ii, 1)) & (xx <= ab(ii + 1, 1)));
     lh    = xx(kk) - ab(ii, 1);
     rh    = (yy(kk) - ab(ii, 2)) * (ab(ii + 1, 1) - ab(ii, 1)) + ...
             ab(ii, 2) * (xx(kk) - ab(ii, 1));
             ab(ii + 1, 2) = lh' * rh / (lh' * lh);
     slope = (ab(ii + 1, 2) - ab(ii, 2)) / (ab(ii + 1, 1) - ab(ii, 1));
     tmp   = [tmp; yy(kk) - slope * (xx(kk) - ab(ii, 1)) - ab(ii, 2)];
  end
  rr = tmp' * tmp;

end % MinRes.
end % BrokenStickRegression.
% EOF BrokenStickRegression.


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

% Copyright 2002-2003 The MathWorks, Inc.

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

% Last Modified by GUIDE v2.5 24-Mar-2005 13:43:51

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @scoutdetails_OpeningFcn, ...
                   'gui_OutputFcn',  @scoutdetails_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 scoutdetails is made visible.
function scoutdetails_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to scoutdetails (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

%Plot helicopter pictures on startup
scout = imread('scout.bmp');
axes(handles.axes1);
image(scout);
axis off;

yaw_pedals = imread('Yaw Pedals.jpg');
axes(handles.axes3);
image(yaw_pedals);
axis off;

collective_lever = imread('Collective Stick.jpg');
axes(handles.axes4);
image(collective_lever);
axis off;

cyclic_stick = imread('Cockpit Instruments.jpg');
axes(handles.axes5);
image(cyclic_stick);
axis off;

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

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

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

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

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

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

close scoutdetails;

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

handles = guidata(scoutdetails);

set(handles.tab1,'Visible','on');
set(handles.tab2,'Visible','off');
set(handles.tab3,'Visible','off');

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

handles = guidata(scoutdetails);

set(handles.tab1,'Visible','off');
set(handles.tab2,'Visible','on');
set(handles.tab3,'Visible','off');

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

set(handles.tab1,'Visible','off');
set(handles.tab2,'Visible','off');
set(handles.tab3,'Visible','on');


text:
function [ha] = ft22ha(ft2)
% Convert area from square feet to hectares.
% Chad A. Greene 2012
ha = ft2*0.000009290304 ;


text:
function test_ft_plot_mesh

% MEM 2gb
% WALLTIME 00:15:00
% DEPENDENCY ft_plot_mesh

%% the first is a simple triangluar mesh

bnd.pos = randn(20,3);
bnd.tri = delaunay(bnd.pos(:,1:2));

figure
ft_plot_mesh(bnd);

%% the following was made by Lilla Magyari in 2013 and pertains to FEM meshes

% create an example segmentation
example.dim = [91 104 111]; % slightly different numbers
example.transform = eye(4);
example.coordsys = 'ctf';
example.unit = 'mm';
example.seg = zeros(example.dim);

% adjusting transformation matrix: center of the head-coordinates [0 0 0] should
% be the center of volume

% center of volume in voxel coordinates
x = round(example.dim(1)/2);
y = round(example.dim(2)/2);
z = round(example.dim(3)/2);

x = round(x);
y = round(y);
z = round(z);

origin = [x y z];

example.transform(1:4,4) = [-origin(:); 1]; % head-coordinate [0 0 0] is in the center of
% the volume (x y z in voxel-coordinates)

% compute position for each voxel in voxelspace and in headspace

[X, Y, Z] = ndgrid(1:example.dim(1), 1:example.dim(2), 1:example.dim(3));
voxelpos = [X(:) Y(:) Z(:)];
headpos = ft_warp_apply(example.transform, voxelpos);

% create 3 spheres

radius1 = 40;
radius2 = 30;
radius3 = 20;

for i = 1:size(headpos,1)
  % from small to large
  if norm(headpos(i,:))<radius3
    example.seg(i) = 3;
  elseif norm(headpos(i,:))<radius2
    example.seg(i) = 2;
  elseif norm(headpos(i,:))<radius1
    example.seg(i) = 1;
  end
end

% convert it to probabilistic
seg = ft_datatype_segmentation(example,'segmentationstyle','probabilistic');

% create smaller segmentation

example.dim = [10 11 9]; % slightly different numbers
example.transform = eye(4);
example.coordsys = 'ctf';
example.unit = 'mm';
example.seg = zeros(example.dim);

% adjusting transformation matrix: center of the head-coordinates [0 0 0] should
% be the center of volume

% center of volume in voxel coordinates
x = round(example.dim(1)/2);
y = round(example.dim(2)/2);
z = round(example.dim(3)/2);

x = round(x);
y = round(y);
z = round(z);

origin = [x y z];

example.transform(1:4,4) = [-origin(:); 1]; % head-coordinate [0 0 0] is in the center of
% the volume (x y z in voxel-coordinates)

% compute position for each voxel in voxelspace and in headspace

[X, Y, Z] = ndgrid(1:example.dim(1), 1:example.dim(2), 1:example.dim(3));
voxelpos = [X(:) Y(:) Z(:)];
headpos = ft_warp_apply(example.transform, voxelpos);

% create 3 spheres

radius1 = 4;
radius2 = 3;
radius3 = 2;

for i = 1:size(headpos,1)
  % from small to large
  if norm(headpos(i,:))<radius3
    example.seg(i) = 3;
  elseif norm(headpos(i,:))<radius2
    example.seg(i) = 2;
  elseif norm(headpos(i,:))<radius1
    example.seg(i) = 1;
  end
end

% convert it to probabilistic
seg_small = ft_datatype_segmentation(example,'segmentationstyle','probabilistic');

%% create mesh
cfg = [];
cfg.method = 'hexahedral';
%cfg.tissue =
mesh = ft_prepare_mesh(cfg,seg);

% 1 cube mesh
mesh0.pos = mesh.pos;
mesh0.hex = mesh.hex(1,:);
mesh0.tissue = mesh.tissue(1,:);
mesh0.tissuelabel = mesh.tissuelabel;
mesh0.unit = mesh.unit;

% few hundreds cube mesh

cfg = [];
cfg.method = 'hexahedral';
%cfg.tissue =
mesh2 = ft_prepare_mesh(cfg,seg_small);

%% plot mesh

% larger mesh
figure; ft_plot_mesh(mesh,'surfaceonly','yes');
% smaller mesh
figure; ft_plot_mesh(mesh2,'surfaceonly','yes');
% smallest mesh
figure; ft_plot_mesh(mesh0,'surfaceonly','yes');


text:
% Small technical example.
% Demonstrates how to show the details on the rotated quincunx grid.
%
disp('Small technical example.');
disp('Demonstrates how to show the details on the rotated quincunx grid.');
disp('FOR MORE INFORMATION:  help rota1001fill');
disp(' ');
disp('See also the report http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04178D.pdf');
disp('Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>');
disp(' (C) 1998-2006 Stichting CWI, Amsterdam, The Netherlands');
disp(' ');
%---INSERT YOUR IMAGE HERE-----------------------------------------------------
if exist('imread','file') == 2
  Orig = double(imread('zenithgray.TIF','tiff'));
else
  load zenithgray; Orig = zenithgray; clear zenithgray;
end
%
disp([' Dimensions of original      ' int2str( size(Orig) )]);
%---PARAMETERS-----------------------------------------------------------------
% How to execute, set parameters
%
 N=2;
%
 filtername='maxmin';
%
% Manage output, set preferences: see printshop
 dops   =3;
 dotiff =0;
%------------------------------------------------------------------------------
% Show original image
printshop('figOriginal', ' Original ', Orig, dops, dotiff, []);
%---DECOMPOSITION--------------------------------------------------------------
disp([' Filter type is ' filtername]);
disp([' Number of scales asked for is ' int2str(N)]);
[C,S] = QLiftDec2(Orig,N,filtername);
%
% Show the details (at level 1) on the rotated quincunx grid
[Detail10, Detail01] = retrieveQ1001(1, 'd', C, S);
background=max(max(max(Detail10)), max(max(Detail01)));
RotaDet=rota1001fill(Detail10, Detail01, background);
printshop('figDetail', 'Detail on the rotated quincunx grid', ...
          RotaDet, dops, dotiff, []);
%
Approx = retrieveR(N, 'a', C, S);
printshop('figApprox', 'Approximation', Approx, dops, dotiff, []);


text:
function y = cimage5(x,xsat,i0,i1)

% function y = cimage5(x,xsat,i0,i1)
% Plot a complex matrix x as a colour image, using Matlab5 RGB mode.
% xsat, if given, defines the amplitude of x which corresponds
% to saturated colours.  (If xsat<0, xmax is used.)
% i0,i1, if given, define the grey and sat. colour intensities.
% i0 is the minimum grey intensity at the centre of the circle
% and i1 is the intensity at the circumference.
% The colours are mapped around the circle as:
%           yellow
%       ?        orange
%   green            red
%      cyan     magenta
%            blue
%
% If an output is not requested, the image is plotted.
% Otherwise the image of map indices is returned as y.
%
% For a plot of the colour circle with a white centre, use:
% t=ones(33,1)*[-16:16]; cimage5(t-j*t.',16,0.9,1.0)
% For a dark centre use:
% t=ones(33,1)*[-16:16]; cimage5(t-j*t.',16,0.3,1.2)

% Nick Kingsbury, Cambridge University, July 1998.

if nargin < 4, i1 = 1.0; end
if nargin < 3, i0 = 0.9; end

xmax = max(max(abs(x)));
if nargin < 2, xsat = xmax; end
if xsat<0, xsat = xmax; end

fprintf(1,'Fig %.0f: xsat = %f, xmax = %f\n', gcf, xsat, xmax);

% Scale x and limit to a max amplitude of 1.
x = x ./ xsat;
x = x ./ max(abs(x),1);

% Calculate red, green and blue intensities in cx.
ax = (i0 - 0.5*i1) * (1 - abs(x));
cx(:,:,1) = (0.5 * i1) * (real(x) + 1) + ax;
cx(:,:,2) = (0.25 * i1) * (2 - real(x) + imag(x)) + ax;
cx(:,:,3) = (0.5 * i1) * (1 - imag(x)) + ax;
cx = max(min(cx,1),0);

if nargout == 0,
% Draw the image.
  image(cx)
  axis image
else
  y = cx;
end

return


text:
function [fidQRS, QRScorr, fidST, Tcorr] = Adjust_Fiducials(ecg, fidBase, q2f, f2s, qs_templ, st, st_templ, Param)
% OVERVIEW, The function adjusts the fiducial base marker for each
% beat so that the correlation between the QRS template and ST template is
% maximium.
%
% INPUTS        MANDATORY       DESCRIPTION
%               ecg             N by M array of ecg data. Contains M
%                               channels of ecg with N data points each.
%
%               fidBase         The location of R peaks for an upright QRS complex.
%
%               q2f             The number of samples between q point and fiducial base in the qrs complex template.
%
%               f2s             The number of samples between fiducial base
%                               and the s point in the qrs complex template.
%
%               qs_templ        qrs complex template.
%
%               st              s-t segment length estimate.
%
%               st_templ        s-t segment template.
%
%               Param           Structure containing parameters using for
%                               computing TWAs.
%
% OUTPUTS
%               fidQRS          The adjusted fiducial base so that the
%                               correlation between the QRS for each beat and template
%                               QRS is maximum.
%
%               QRScorr         The maximum correlation between QRS complex
%                               for each beat and template QRS.
%
%               fidST           The adjusted fiducial base so that the
%                               correlation between the S-T segment for each beat and template S-T segment is maximum.
%
%               Tcorr           The maximum correlation between S-T segment
%                               for each beat and the template S-T segment.
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Shamim Nemati,
%       editted by Ismail Sadiq on 10/26/2019.
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

NInapr = 0;

qoff = -q2f - Param.stAdjIntv - 1;
soff = f2s - Param.stAdjIntv - 1;
toff = f2s + st - Param.stAdjIntv - 1;
Tcorr=[];
warning off MATLAB:divideByZero     % to supress corrcoef warning when data is constant, in that case corrcoef is NaN and max index is any

for i = 1:length(fidBase)
    cc = zeros(1, 4.5 * Param.stAdjIntv + 1);
    ccc = zeros(1, 4.5 * Param.stAdjIntv + 1);
    for j = 1:(4.5 * Param.stAdjIntv + 1)
        if((fidBase(i) + j + qoff)<1), continue; end
        qs_ecg = ecg((fidBase(i) + j + qoff):(fidBase(i) + j + soff));
        if(length(qs_templ)~=length(qs_ecg)), continue;
        else
            a = corrcoef(qs_templ, qs_ecg);
            cc(j) = a(1, 2);
        end

        if strcmp(Param.Alignment, 'st')
            if((fidBase(i) + j + soff)<1), continue;end
            st_ecg = ecg((fidBase(i) + j + soff):(fidBase(i) + j + toff));
            if (length(st_templ)~=length(st_ecg)), continue;
            else
                b = corrcoef(st_templ, st_ecg);
                ccc(j) = b(1, 2);
            end
        end;
    end;
    [m ind] = max(cc);
    indQRS = ind;
    if (strcmp(Param.Alignment, 'st'))
        [mm ind] = max(ccc);
    end;

    if (m < Param.corrQRS || (strcmp(Param.Alignment, 'st') && mm < Param.corrT))
        NInapr = NInapr + 1;
    end;

    if (NInapr > 0.1 * length(fidBase))
        fidST = [];
        fidQRS = [];
        QRScorr = [];
        Tcorr = [];
        break;
    else
        fidST(i) = fidBase(i) + ind - Param.stAdjIntv - 1; % with respect to ST alignment
        fidQRS(i) = fidBase(i) + indQRS - Param.stAdjIntv - 1; % with respect to QRS alignment
        QRScorr(i) = m;
        if (strcmp(Param.Alignment, 'st'))
            Tcorr(i) = mm;
        end;

    end;
end;

return;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function opt=varargin2struct(varargin)
%
% opt=varargin2struct('param1',value1,'param2',value2,...)
%   or
% opt=varargin2struct(...,optstruct,...)
%
% convert a series of input parameters into a structure
%
% authors:Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2012/12/22
%
% input:
%      'param', value: the input parameters should be pairs of a string and a value
%       optstruct: if a parameter is a struct, the fields will be merged to the output struct
%
% output:
%      opt: a struct where opt.param1=value1, opt.param2=value2 ...
%
% 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)
%

len=length(varargin);
opt=struct;
if(len==0) return; end
i=1;
while(i<=len)
    if(isstruct(varargin{i}))
        opt=mergestruct(opt,varargin{i});
    elseif(ischar(varargin{i}) && i<len)
        opt=setfield(opt,varargin{i},varargin{i+1});
        i=i+1;
    else
        error('input must be in the form of ...,''name'',value,... pairs or structs');
    end
    i=i+1;
end

author: Guanghan

repo: GNet-pose

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


text:
classdef SO3KernelHandle < SO3Kernel
  % defines a kernel function as a function of the rotational angle

  properties
    fun = @(x) 1;
  end

  methods

    function psi = SO3KernelHandle(fun,varargin)

      % extract parameter and halfwidth
      if nargin == 0, return;end

      psi.fun = fun;
      psi.A = calcFourier(psi,getMTEXpref('maxS2Bandwidth'),varargin{:});

    end

    function c = char(psi)
      c = ['custom, halfwidth ' ...
        xnum2str(psi.halfwidth/degree) mtexdegchar];
    end

    function value = eval(psi,co2)
      % the kernel function on SO(3)
      co2 = cut2unitI(co2);
      value   =  psi.fun(co2);
    end

  end
end


text:
function [A,M] = Poisson2D();

% function [A,M] = Poisson2D()
% Purpose: Set up matrix for 2D Poisson equation based on stabilized
%          internal fluxes on symmetric form

Globals2D;
g = zeros(K*Np,1);
A = spalloc(K*Np, K*Np, 3*Np);  M = spalloc(K*Np, K*Np, 3*Np);

% Build matrix -- one column at a time
for i=1:K*Np
    g(i) = 1.0;
    gmat = reshape(g,Np,K);
    [Avec,Mvec] = PoissonRHS2D(gmat);

    ids = find(Avec); A(ids,i) = Avec(ids);
    ids = find(Mvec); M(ids,i) = Mvec(ids);
    g(i)=0.0;
end
return

repo: Aero-matlab


text:
function btv_test09 ( )

%*****************************************************************************80
%
%% BTV_TEST09 tests BURGERS_TIME_VISCOUS with the spike initial condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTV_TEST09\n' );
  fprintf ( 1, '  Test BURGERS_TIME_VISCOUS with the spike initial condition.\n' );

  nx = 81;
  nt = 400;
  t_max = 4.0;
  nu = 0.01;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial condition: spike\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Viscosity = %g\n', nu );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_viscous ( @ic_spike, nx, nt, t_max, nu, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 9 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers equation solutions over time, initial condition spike' )

  filename = 'btv_test09.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function metric = getStatisticDose(setDoseNumber1, setDoseNumber2, planC);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

dose1 = planC{planC{end}.dose}(setDoseNumber1).doseArray;

[xV yV zV] = getDoseXYZVals(planC{planC{end}.dose}(setDoseNumber1));

[x,y,z] = meshgrid(xV,yV,zV);

[dose2] = getDoseAt(setDoseNumber2, x, y, z);

maxDose1 = max(dose1(:));
maxDose2 = max(dose2(:));

dose1(dose1 == 0) = NaN;
dose2(dose2 == 0) = NaN;

dose1Norm = dose1/maxDose1;
dose2Norm = dose2/maxDose2;

dose2CorMax = (maxDose1/maxDose2)*dose2;

diffAbsNorm = abs(dose1Norm - dose2Norm);

metric.Max_Abs_Diff_Norm = max(diffAbsNorm(:));

diffAbsNormPr = 100*diffAbsNorm;

Frac = diffAbsNormPr(~isnan(diffAbsNormPr));

metric.Mean_Abs_Diff_Norm = mean(Frac/100);

Frac2 = find(Frac > 2);

Frac5 = find(Frac > 5);

Frac10 = find(Frac > 10);

metric.Fr2Pr = 100*length(Frac2)/length(Frac);

metric.Fr5Pr = 100*length(Frac5)/length(Frac);

metric.Fr10Pr = 100*length(Frac10)/length(Frac);

diffNormSq = (dose1Norm - dose2Norm).^2;

FracSq = diffNormSq(~isnan(diffNormSq));

metric.rmse = 100*sqrt(sum(FracSq(:))/length(FracSq));