text:
%kminvdiag 'Compute inverse of diagonal matrix '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros minvdiag.pane file
%
% Parameters:
% InputFile: i 'Input Matrix', required: 'Input matrix'
% OutputFile: o 'Output Matrix', required: 'Output matrix'
%
% Example: o = kminvdiag(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% minvdiag - Compute Inverse of Diagonal Matrix
%
%  DESCRIPTION
% .I minvdiag
% computes the inverse of a diagonal matrix. The matrix can be non-square.
% Output is of type KDCOMPLEX if the input matrix is a complex type, otherwise
% the output is of type KDOUBLE.
%
% No checking is done for a divide by zero condition.
%
% Note that only the diagonal elements of the input matrix are modified.
%
% "minvdiag\fR is a pane object pointing to \fIminvert\fR.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kminvdiag(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kminvdiag(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'minvert"  -diag'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


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

nsteps = 1;
Vfrac_final = 0.75;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
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 y = compute_orthogonal_projection( x, A, C )

% compute_orthogonal_projection - computer orthogonal projection
%
%   y = compute_orthogonal_projection( x, A, C );
%
%   y is the pojection of x onto the intersection of the sets
%       { y \ <x-cj,aj>=0 }
%   where the cj and the aj are the columns of C and A.
%
%   Copyright (c) 2007 Gabriel Peyre

% compute solution as
%   y = x + A*lambda

D = (A'*A);
if abs(det(D))<1e-9
    D = D+rand(size(D))*1e-9;
end
lambda = D \ ( -A'*x + diag(A'*C) );
y = x + A*lambda;

author: gpeyre

repo: matlab-toolboxes


text:
function result = matc(sf,tf)
% Compare two mat files
% matc matfile1.mat matfile2.mat

% Copyright 2004, Dr. Xianyao Chen, xianyaochen@yahoo.com.cn
% Last Modified by Xianyao Chen, 08-June-2004

if (nargin ~= 2)
  return
end
if exist(sf) ~=2
    msg = sprintf('%s%s','Not valid mat file: ',sf);
%     fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
   fprintf('%s',msg)
    return
end
if exist(tf) ~=2
    msg = sprintf('%s%s','Not valid mat file: ',tf);
%     fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
   fprintf('%s',msg)
    return
end
save #~tmp1.mat -mat
clear
load #~tmp1 sf

load(sf)
clear sf

ws = whos;
save #~tmp2.mat -mat
clear

load #~tmp1
load #~tmp2
delete('#~tmp2.mat');

m1 = size(ws,1);

warning off
% nPos = 0;
for k = 1:m1
    var = ws(k).name;
    as = load(sf,var);
    at = load(tf,var);
    atn = fieldnames(at);
    if k == 1
        msg = sprintf('%s%s','Source File: ',sf,' --- Target File: ',tf');
        fprintf('%s\n',msg)
        msg = sprintf('%50s','--------------------------------------------------');
        fprintf('%s\n',msg)
%         msg = sprintf('%s','----------------------------------------------------------------------------------------------------------------------------');
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
    end

    if isempty(atn)
        msg = sprintf('%12s%12s%6s%s','Variable :: ',var,'  ||  ',' not exist');
        fprintf('%s\n',msg)
%         msg = sprintf('%12s%12s%20s%20s%20s%20s%s','Variable :: ',var,'    |  Source File: ',sf,' ---  Target File: ',tf,' :: not exist');
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
        continue
    end
    aas = struct2cell(as);
    aat = struct2cell(at);
    if ~isequal(aas{1},aat{1})
        msg = sprintf('%12s%12s%6s%s','Variable :: ',var,'  ||  ',' not equal');
        fprintf('%s\n',msg)
%         msg = sprintf('%12s%12s%20s%20s%20s%20s%s','Variable :: ',var,'    |  Source File: ',sf,' ---  Target File: ',tf,' :: not equal');
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
        continue
    else
        msg = sprintf('%12s%12s%6s%s','Variable :: ',var,'  ||  ',' equal');
        fprintf('%s\n',msg)
%         msg = sprintf('%12s%12s%20s%20s%20s%20s%s','Variable :: ',var,'    |  Source File: ',sf,' ---  Target File: ',tf,' :: Exactly Equal');
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
        continue
    end
end
msg = sprintf('%50s','--------------------------------------------------');
fprintf('%s\n\n',msg)

% msg = sprintf('%s','----------------------------------------------------------------------------------------------------------------------------');
% fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
% fprintf('\n')

clear
load #~tmp1 tf

load(tf)
clear tf

ws = whos;
save #~tmp2.mat -mat
clear

load #~tmp1
load #~tmp2
delete('#~tmp2.mat');

m1 = size(ws,1);
warning off

for k = 1:m1
    var = ws(k).name;
    as = load(tf,var);
    at = load(sf,var);
    atn = fieldnames(at);
    if k == 1
        msg = sprintf('%s%s','Source File: ',tf,' --- Target File: ',sf');
        fprintf('%s\n',msg)
        msg = sprintf('%50s','--------------------------------------------------');
        fprintf('%s\n',msg)
        %fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
        %fprintf('\n')
    end

    if isempty(atn)
        msg = sprintf('%12s%12s%6s%s','Variable :: ',var,'  ||  ',' not exist');
        fprintf('%s\n',msg)
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
        continue
    end
    aas = struct2cell(as);
    aat = struct2cell(at);
    if ~isequal(aas{1},aat{1})
        msg = sprintf('%12s%12s%6s%s','Variable :: ',var,'  ||  ',' not equal');
        fprintf('%s\n',msg)
%         msg = sprintf('%12s%12s%20s%20s%20s%20s%s','Variable :: ',var,'    |  Source File: ',tf,' ---  Target File: ',sf,' :: not equal');
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
        continue
    else
        msg = sprintf('%12s%12s%6s%s','Variable :: ',var,'  ||  ',' equal');
        fprintf('%s\n',msg)
%         msg = sprintf('%12s%12s%20s%20s%20s%20s%s','Variable :: ',var,'    |  Source File: ',tf,' ---  Target File: ',sf,' :: Exactly Equal');
%         fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
%         fprintf('\n')
        continue
    end
end
warning on

msg = sprintf('%50s','--------------------------------------------------');
fprintf('%s\n',msg)
% fprintf( [ repmat( '\b', 1, nPos ), '%s' ], msg),
% fprintf('\n')

delete('#~tmp1.mat')

author: Sable

repo: mcbench-benchmarks


text:
function value = lrline ( xu, yu, xv1, yv1, xv2, yv2, dv )

%*****************************************************************************80
%
%% LRLINE determines if a point is left of, right or, or on a directed line.
%
%  Discussion:
%
%    The directed line is parallel to, and at a signed distance DV from
%    a directed base line from (XV1,YV1) to (XV2,YV2).
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, real XU, YU, the coordinates of the point whose
%    position relative to the directed line is to be determined.
%
%    Input, real XV1, YV1, XV2, YV2, the coordinates of two points
%    that determine the directed base line.
%
%    Input, real DV, the signed distance of the directed line
%    from the directed base line through the points (XV1,YV1) and (XV2,YV2).
%    DV is positive for a line to the left of the base line.
%
%    Output, integer VALUE, the result:
%    +1, the point is to the right of the directed line;
%     0, the point is on the directed line;
%    -1, the point is to the left of the directed line.
%
  tol = 100.0 * r8_epsilon ( );

  dx = xv2 - xv1;
  dy = yv2 - yv1;
  dxu = xu - xv1;
  dyu = yu - yv1;

  tolabs = tol * max ( abs ( dx ), ...
                 max ( abs ( dy ), ...
                 max ( abs ( dxu ), ...
                 max ( abs ( dyu ), abs ( dv ) ) ) ) );

  t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy );

  if ( tolabs < t )
    value = 1;
  elseif ( -tolabs <= t )
    value = 0;
  else
    value = -1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [p,f,p1,f1]=optimitaestep(n)
%OPTIMITAESTEP Calculate the optimal ITAE trasnfer function coefficient for a step input
%
% Usage:
% [p,itae,p0,itae0] = optimitaestep(n) returns the optimal coefficients of a
% nth-order transfer function based on the ITAE criterion for a step input.
% Input:  n - the order of transfer function
% Output: p - the optimal coefficients of the nth-order transfer function
%         itae - the optimal ITAE criterion.
%         p0 - original coefficients derivided by Graham and Lathrop.
%         itae0 - the ITAE citerion correcponding to p0.
%
% Background:
% The original ITAE (stands for Integral of Time multiplied by Absolute
% Error) coefficient table was drived by D. Graham and R.C. Lathrop through
% an analog computer. The table has been widely adopted as a standard in
% most undergraduate level control engineering text books. However, the
% optimality of these coefficients has been questioned by several researchers.
% This code provides a means to re-calculate these coefficient using more
% advanced numerical optimization techniques in digital computers.
%
% References:
% 1. D. Graham and R.C. Lathrop, "The Synthesis of Optimum Response: Criteria
%    ans Standard Forms, Part 2", Transactions of the AIEE 72, Nov. 1953, pp. 273-288
% 2. Y. Cao, "Correcting the minimum ITAE standard forms of
%    zero-displaceemnt-error systems", Journal of Zhejiang University
%    (Natural Science) Vol. 23, N0o.4, pp. 550-559, 1989.
%
% Example, the 7th-order optimal and original ITAE coefficients
%{
format compact
format short g
[p,f,p1,f1] = optimitaestep(7)
%}
% p =
%    1  2.2172  6.7447  9.3493  11.58   8.6799  4.3233  1
% f =
%    10.585
% p1 =
%    1  4.475   10.42   15.08   15.54   10.64   4.58    1
% f1 =
%    14.953
%

%Input check
error(nargchk(1,2,nargin));
if n<2 || n>8
    error('n must be within range of [2 8].')
end
itaeOrg={ 1.4, [1.75 2.15], [2.1 3.4 2.7],...
       [2.8 5.0 5.5 3.4], [3.25 6.60 8.60 7.45 3.95],...
       [4.475 10.42 15.08 15.54 10.64 4.58],...
       [5.2 12.8 21.6 25.75 22.2 13.3 5.15]};
p1=itaeOrg{n-1};
p0=fliplr(p1);
opt=optimset('display','off','TolX',1e-9,'TolFun',1e-9,'LargeScale','off');
dt=0.01;
tf=50;
flag=0;
while ~flag
    cost=@(x)itaecost(x,dt,tf);
    [p0,f,flag]=fminunc(cost,p0,opt);
    p=[1 fliplr(p0) 1];
    if any(real(roots(p))>0)
        tf=tf*2;
        flag=0;
        p0=fliplr(itaeOrg{n-1});
    end
end
if nargout>2
    f1=cost(fliplr(p1));
    p1=[1 p1 1];
end

function E=itaecost(p,dt,tf)
n=numel(p);
A=[zeros(n,1) eye(n);-1 -p];
B=[zeros(n,1);1];
A=expm([A B;zeros(1,n+2)]*dt);
x=[zeros(n+1,1);1];
E=0;
for t=0:dt:tf
    tdt=t*dt;
    x=A*x;
    e=1-x(1);
    E=E+abs(e)*tdt;
end


text:
function [g,f] = octspec(B,A,Fs,Fc,s,n);
% OCTSPEC Plots an octave filter characteristics.
%    OCTSPEC(B,A,Fs,Fc) plots the attenuation of the filter defined by
%    B and A at sampling frequency Fs. Fc is the center frequency of
%    the octave filter. The plot covers one decade on both sides of Fc.
%
%    OCTSPEC(B,A,Fs,Fc,'ANSI',N) superposes the ANSI Order-N analog
%    specification for comparison. Default is N = 3.
%
%    OCTSPEC(B,A,Fs,Fc,'IEC',N) superposes the characteristics of the
%    IEC 61260 class N specification for comparison. Default is N = 1.
%
%    [G,F] = OCTSPEC(B,A,Fs,Fc) returns two 512-point vectors with
%    the gain (in dB) in G and logarithmically spaced frequencies in F.
%    The plot can then be obtained by SEMILOGX(F,G)
%
%    See also OCTDSGN, OCT3SPEC, OCT3DSGN.

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
%         couvreur@thor.fpms.ac.be
% Last modification: Sept. 4, 1997, 10:30am.

% References:
%    [1] ANSI S1.1-1986 (ASA 65-1986): Specifications for
%        Octave-Band and Fractional-Octave-Band Analog and
%        Digital Filters, 1993.
%    [2] IEC 61260 (1995-08):  Electroacoustics -- Octave-Band and
%        Fractional-Octave-Band Filters, 1995.

if (nargin < 4) | (nargin > 6)
  error('Invalid number of input arguments.');
end

ansi = 0;
iec = 0;
if nargin > 4
  if strcmp(lower(s),'ansi')
    ansi = 1;
    if nargin == 5
      n = 3;
    end
  elseif strcmp(lower(s),'cei') | strcmp(lower(s),'iec')
    iec = 1;
    if nargin == 5
      n = 1
    end
    if (n < 0) | (n > 3)
      error('IEC class must be 0, 1, or 2');
    end
  end
end

N = 512;
pi = 3.14159265358979;
F = logspace(log10(Fc/10),log10(min(Fc*10,Fs/2)),N);
H = freqz(B,A,2*pi*F/Fs);
G = 20*log10(abs(H));

% Set output variables.
if nargout ~= 0
  g = G; f = F;
  return
end

% Generate the plot
if (ansi) 				% ANSI Order-n specification
  f = logspace(log10(Fc/10),log10(Fc*10),N);
  f1 = Fc/sqrt(2);
  f2 = Fc*sqrt(2);
  Qr = Fc/(f2-f1);
  Qd = (pi/2/n)/(sin(pi/2/n))*Qr;
  Af = 10*log10(1+Qd^(2*n)*((f/Fc)-(Fc./f)).^(2*n));
  semilogx(F,G,f,-Af,'--');
  legend('Filter',['ANSI order-' int2str(n)],0);
elseif (iec) 					% CEI specification
  semilogx(F,G);
  hold on
  if n == 0
    tolup =  [ .15 .15 .15 .15 .15 -2.3 -18.0 -42.5 -62 -75 -75 ];
    tollow = [ -.15 -.2 -.4 -1.1 -4.5 -realmax -inf -inf -inf -inf -inf ];
  elseif n == 1
    tolup =  [ .3 .3 .3 .3 .3 -2 -17.5 -42 -61 -70 -70 ];
    tollow = [ -.3 -.4 -.6 -1.3 -5 -realmax -inf -inf -inf -inf -inf ];
  elseif n == 2
    tolup =  [ .5 .5 .5 .5 .5 -1.6 -16.5 -41 -55 -60 -60  ];
    tollow = [ -.5 -.6 -.8 -1.6 -5.5 -realmax -inf -inf -inf -inf -inf ];
  end
  U = 2;
  f = Fc * U.^[ 0 1/8 1/4 3/8 1/2 1/2 1 2 3 4 NaN ];
  ff = Fc * U.^[ 0 -1/8 -1/4 -3/8 -1/2 -1/2 -1 -2 -3 -4 NaN ];
  f(length(f)) = realmax;
  ff(length(ff)) = realmin;
  semilogx(F,G,f,tolup,'--');
  semilogx(F,G,f,tollow,'--');
  semilogx(F,G,ff,tolup,'--');
  semilogx(F,G,ff,tollow,'--');
  hold off
  legend('Filter',['IEC class ' int2str(n)],0);
else
  semilogx(F,G);
end
xlabel('Frequency [Hz]'); ylabel('Gain [dB]');
title(['Octave filter: Fc =',int2str(Fc),' Hz, Fs = ',int2str(Fs),' Hz']);
axis([Fc/10 Fc*10 -80 5]);
grid on


text:
function cycle_brent_test05 ( )

%*****************************************************************************80
%
%% CYCLE_BRENT_TEST05 tests CYCLE_BRENT for F5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLE_BRENT_TEST05\n' );
  fprintf ( 1, '  Test CYCLE_BRENT for F5().\n' );
  fprintf ( 1, '  f5(i) = mod ( 16383 * i + 1, 65536 ).\n' );

  x0 = 1;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Starting argument X0 = %d\n', x0 );

  [ lam, mu ] = cycle_brent ( @f5, x0 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reported cycle length is %d\n', lam );
  fprintf ( 1, '  Expected value is 8\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reported distance to first cycle element is %d\n', mu );
  fprintf ( 1, '  Expected value is 0\n' );

  i = 0;
  x0 = 1;
  fprintf ( 1, '  %d  %d\n', i, x0 );
  for i = 1 : 10
    x0 = f5 ( x0 );
    fprintf ( 1, '  %d  %d\n', i, x0 );
  end

  return
end


text:
function W = todouble(W)

for j = 1:numel(W)
  if (iscell(W(j).W))
    for k=1:numel(W(j).W)
      W(j).W{k} = double(W(j).W{k});
    end
  else
    W(j).W = double(W(j).W);
  end
  W(j).params.l = double(W(j).params.l);
end

author: norouzi

repo: hdml

save_path: github-repos/MATLAB/norouzi-hdml


text:
function opt = mpoption_info_cplex(selector)
%MPOPTION_INFO_CPLEX  Returns MATPOWER option info for CPLEX.
%
%   DEFAULT_OPTS = MPOPTION_INFO_CPLEX('D')
%   VALID_OPTS   = MPOPTION_INFO_CPLEX('V')
%   EXCEPTIONS   = MPOPTION_INFO_CPLEX('E')
%
%   Returns a structure for CPLEX options for MATPOWER containing ...
%   (1) default options,
%   (2) valid options, or
%   (3) NESTED_STRUCT_COPY exceptions for setting options
%   ... depending on the value of the input argument.
%
%   This function is used by MPOPTION to set default options, check validity
%   of option names or modify option setting/copying behavior for this
%   subset of optional MATPOWER options.
%
%   See also MPOPTION.

%   MATPOWER
%   Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    selector = 'D';
end
if have_feature('cplex')
    switch upper(selector)
        case {'D', 'V'}     %% default and valid options
            opt = struct(...
                'cplex',        struct(...
                    'lpmethod',     0, ...
                    'qpmethod',     0, ...
                    'opts',         [], ...
                    'opt_fname',    '', ...
                    'opt',          0 ...
                ) ...
            );
        case 'E'            %% exceptions used by nested_struct_copy() for applying
            opt = struct(...
                'name',         { 'cplex.opts' }, ...
                'check',        0 ...
                );
%                 'copy_mode',    { @cplex_options } ...
        otherwise
            error('mpoption_info_cplex: ''%s'' is not a valid input argument', selector);
    end
else
    opt = struct([]);       %% CPLEX is not available
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function transpose(x)
%TRANSPOSE is not defined on tensors.
%
%   See also TENSOR, TENSOR/PERMUTE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

error('Transpose on tensor is not defined');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function bnet = mk_dbn(intra, inter, node_sizes, varargin)
% MK_DBN Make a Dynamic Bayesian Network.
%
% BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs
% from i in slice t to j in slice t iff intra(i,j) = 1, and
% from i in slice t to j in slice t+1 iff inter(i,j) = 1,
% for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1.
% node_sizes(i) is the number of values node i can take on.
% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
% See also mk_bnet.
%
% Optional arguments [default in brackets]
% 'discrete' - list of discrete nodes [1:n]
% 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ]
% 'eclass1' - equiv class for slice 1 [1:n]
% 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]
%    equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j},
%    i.e., nodes i and j have tied parameters.
% 'intra1' - topology of first slice, if different from others
% 'names' - a cell array of strings to be associated with nodes 1:n [{}]
%    This creates an associative array, so you write e.g.
%     'evidence(bnet.names{'bar'}) = 42' instead of  'evidence(2} = 42'
%     assuming names = { 'foo', 'bar', ...}.
%
% For backwards compatibility with BNT2, arguments can also be specified as follows
%   bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1)
%
% After calling this function, you must specify the parameters (conditional probability
% distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc.

n = length(intra);
ss = n;
bnet.nnodes_per_slice = ss;
bnet.intra = intra;
bnet.inter = inter;
bnet.intra1 = intra;
dag = zeros(2*n);
dag(1:n,1:n) = bnet.intra1;
dag(1:n,(1:n)+n) = bnet.inter;
dag((1:n)+n,(1:n)+n) = bnet.intra;
bnet.dag = dag;
bnet.names = {};

directed = 1;
if ~acyclic(dag,directed)
  error('graph must be acyclic')
end

bnet.eclass1 = 1:n;
%bnet.eclass2 = (1:n)+n;
bnet.eclass2 = bnet.eclass1;
for i=1:ss
  if isequal(parents(dag, i+ss), parents(dag, i)+ss)
    %fprintf('%d has isomorphic parents, eclass %d\n', i, bnet.eclass2(i))
  else
    bnet.eclass2(i) = max(bnet.eclass2) + 1;
    %fprintf('%d has non isomorphic parents, eclass %d\n', i, bnet.eclass2(i))
  end
end

dnodes = 1:n;
bnet.observed = [];

if nargin >= 4
  args = varargin;
  nargs = length(args);
  if ~isstr(args{1})
    if nargs >= 1, dnodes = args{1}; end
    if nargs >= 2, bnet.eclass1 = args{2}; end
    if nargs >= 3, bnet.eclass2 = args{3}; end
    if nargs >= 4, bnet.intra1 = args{4}; end
  else
    for i=1:2:nargs
      switch args{i},
       case 'discrete', dnodes = args{i+1};
       case 'observed', bnet.observed = args{i+1};
       case 'eclass1',  bnet.eclass1 = args{i+1};
       case 'eclass2',  bnet.eclass2 = args{i+1};
       case 'intra1',  bnet.intra1 = args{i+1};
       %case 'ar_hmm',  bnet.ar_hmm = args{i+1};  % should check topology
       case 'names',  bnet.names = assocarray(args{i+1}, num2cell(1:n));
       otherwise,
	error(['invalid argument name ' args{i}]);
      end
    end
  end
end

bnet.observed = sort(bnet.observed); % for comparing sets
ns = node_sizes;
bnet.node_sizes_slice = ns(:)';
bnet.node_sizes = [ns(:) ns(:)];

cnodes = mysetdiff(1:n, dnodes);
bnet.dnodes_slice = dnodes;
bnet.cnodes_slice = cnodes;
bnet.dnodes = [dnodes dnodes+n];
bnet.cnodes = [cnodes cnodes+n];

bnet.equiv_class = [bnet.eclass1(:) bnet.eclass2(:)];
bnet.CPD = cell(1,max(bnet.equiv_class(:)));
eclass = bnet.equiv_class(:);
E = max(eclass);
bnet.rep_of_eclass = zeros(1,E);
for e=1:E
  mems = find(eclass==e);
  bnet.rep_of_eclass(e) = mems(1);
end

ss = n;
onodes = bnet.observed;
hnodes = mysetdiff(1:ss, onodes);
bnet.hidden_bitv = zeros(1,2*ss);
bnet.hidden_bitv(hnodes) = 1;
bnet.hidden_bitv(hnodes+ss) = 1;

bnet.parents = cell(1, 2*ss);
for i=1:ss
  bnet.parents{i} = parents(bnet.dag, i);
  bnet.parents{i+ss} = parents(bnet.dag, i+ss);
end

bnet.auto_regressive = zeros(1,ss);
% ar(i)=1 means (observed) node i depends on i in the  previous slice
for o=bnet.observed(:)'
  if any(bnet.parents{o+ss} <= ss)
    bnet.auto_regressive(o) = 1;
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Local Regression and Likelihood, Figure 8.3.
%
% Discrimination/Classification, iris data for different
% smooting paramters.
%
% Note that the iris.mat file contains the full iris dataset;
% only Versicolor and Virginica are used in this example.
%
% The `Species' variable contains species names. `Specn' has
% them numerically, Setosa=1, Versicolor=2, Virginica=3.
%
% Author: Catherine Loader
%
% Need: get contour plot correct.
% Need: distinguish colors in plot.

load iris;
a = (2:9)/10;
z = zeros(size(a));

u = find(Specn >= 2);
pw = PetalWid(u);
pl = PetalLen(u);
y = (Specn(u)==3);

for i = 1:length(a)
  fit = locfit([pw pl],y,'deg',1,'alpha',a(i),'ev','cros','scale',0,'family','binomial');
  fv = fitted(fit);
  tb = tabulate(10*y+(fv>=0.5))
  z(i) = sum(y == (fv<=0.5));
end;

[a; z]

fit = locfit([pw pl],y,'deg',1,'scale',0,'family','binomial');
figure('Name','fig8_3: iris classification');
lfplot(fit,'contour');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function bouts = SegsToBouts(segs, moat, flength)
if ~exist('flength', 'var') || isempty(flength)
    flength = segs(end,2);
end
segsPlus = bsxfun(@plus, segs, [-moat moat]);
% remove gaps bigger than moat
candStarts = segsPlus(:,1);
candEnds = segsPlus(:,2);
dontStop = candStarts(2:end)<candEnds(1:end-1);
candEnds(find(dontStop)) = [];
candStarts(find(dontStop)+1) = [];
bouts = [candStarts candEnds];
% bouts(find((segsPlus(2:end,1) - segsPlus(1:end-1,2))<0),:) = []; % gap is too big
bouts(:,1) = max(1,bouts(:,1));
bouts(:,2) = min(flength,bouts(:,2));

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
% Let the user draw a square with the mouse,
% and then click on the corners to do a manual segmentation

ss = 6;
Q1 = 1; Q2 = 2; Q3 = 3; obsvel = 6;
CLOCKWISE = 1; ANTICLOCK = 2;
LR = 1; UD = 2; RL = 3; DU = 4;

% repeat this block manually incrementing the sequence number
% and setting ori.
% (since I don't know how to call getmouse as a call-return function).
seq = 4;
%ori = CLOCKWISE
ori = ANTICLOCK;
clear xpos ypos
getmouse
% end block

% manual segmentation with the mouse
startseg(1) = 1;
for i=2:4
  fprintf('click on start of segment %d\n', i);
  [x,y] = ginput(1);
  plot(x,y,'ro')
  d = dist2([xpos; ypos]', [x y]);
  startseg(i) = argmin(d);
end

% plot corners in green
%ti = first point in (i+1)st segment
t1 = startseg(1); t2 = startseg(2); t3 = startseg(3); t4 = startseg(4);
plot(xpos(t2), ypos(t2), 'g*')
plot(xpos(t3), ypos(t3), 'g*')
plot(xpos(t4), ypos(t4), 'g*')

xvel = xpos(2:end) - xpos(1:end-1);
yvel = ypos(2:end) - ypos(1:end-1);
speed = [xvel(:)'; yvel(:)'];
pos_data{seq} = [xpos(:)'; ypos(:)'];
vel_data{seq} = [xvel(:)'; yvel(:)'];
T = length(xvel);
Q1label{seq} = num2cell(repmat(ori, 1, T));
Q2label{seq} = zeros(1, T);
if ori == CLOCKWISE
  Q2label{seq}(t1:t2) = LR;
  Q2label{seq}(t2+1:t3) = UD;
  Q2label{seq}(t3+1:t4) = RL;
  Q2label{seq}(t4+1:T) = DU;
else
  Q2label{seq}(t1:t2) = RL;
  Q2label{seq}(t2+1:t3) = UD;
  Q2label{seq}(t3+1:t4) = LR;
  Q2label{seq}(t4+1:T) = DU;
end

% pos_data{seq}(:,t), vel_data{seq}(:,t) Q1label{seq}(t) Q2label{seq}(t)
save 'square4' pos_data vel_data Q1label Q2label

nseq = 4;
cases = cell(1,nseq);
for seq=1:nseq
  T = size(vel_data{seq},2);
  ev = cell(ss,T);
  ev(obsvel,:) = num2cell(vel_data{seq},1);
  ev(Q1,:) = Q1label{seq};
  ev(Q2,:) = num2cell(Q2label{seq});
  cases{seq} = ev;
end
save 'square4_cases' cases


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

function [lower, upper] = UpperBound3(S, g, df, B, f, Nr, NSim, NSSim, ...
    getpaths, payoff)
% method from Broadie for computing upper bounds, see Chapter 8

iVec = 1:NSSim;

v = g(:,end);   % start for backward induction
c = zeros(NSim,Nr-1);   % continuation value

% backward induction and regression from t_{Nr-1} up to t_1
for i = Nr-1:-1:1
        index = find(g(:,i) > 0); % all ITM paths
        s = S(index,i+1);         % values of S at given time point
        v = v * df(i+1);          % option value at t_i

        Acell = B(s);             % evaluate basis function in cell array B
        A = cell2mat(Acell{:,:}); % convert to matrix

        c(index,i) = A*f(:,i);                  % continuation value

        exercise = g(index,i) >= c(index,i);    % early exercise
        v(index(exercise)) = g(index(exercise),i);
end

lower = mean(v * df(1));    % final option value

% Computing the martingale numerically, martingale = pi

L = zeros(NSim,1); %erster Teil von delta_i vgl. Formel (6.6)
expectation = zeros(NSim,Nr); %zweiter Teil von delta_i vgl. Formel (6.6)
expectation(:,1) = lower * ones(NSim,1);

for i=1:1:Nr-1
    for j=1:NSim
        expectation(j,i+1) = subsimulation(S(j,i+1), df, B, NSSim, Nr-i,...
            f(:,i:end), iVec(1:NSSim), getpaths, payoff) * prod(df(1:i));
    end
end

pi = zeros(NSim,Nr+1); % stores the values of the constructed martingale

% first time step and last are different
i_exercise = g(:,1) >= c(:,1) & c(:,1) > 0; % exercise in this case

L(i_exercise) = g(i_exercise, 1) * prod(df(1:1));
L(~i_exercise) = expectation(~i_exercise,2);

pi(:,2) = pi(:,1) + L - expectation(:,1);

for i=2:1:Nr-1
    % exercise in this case if not already exercised
    i_exercise = g(:,i) >= c(:,i) & c(:,i) > 0 & ~i_exercise;

    L(i_exercise) = g(i_exercise, i) * prod(df(1:i));
    L(~i_exercise) = expectation(~i_exercise,i+1);

    pi(:,i+1) = pi(:,i) + L - expectation(:,i);
end
% finally exercise if in the money and not already exercised
i_exercise = g(:,Nr) > 0 & ~i_exercise;

pi(i_exercise,Nr+1) = pi(i_exercise,Nr) + g(i_exercise,Nr)* prod(df(1:Nr));
pi(~i_exercise,Nr+1) = pi(~i_exercise,Nr);

% upper bound using the martingale pi

maximum = zeros(NSim,1);

for j=1:1:NSim
    maximum(j) = max(g(j,:) - pi(j,2:end)); %vgl. Formel (6.5)
end
upper = mean(maximum);

end

function y = subsimulation(S0, df, B, NSim, Nr, beta, iVec, gp,payoff)
    S2 = gp(S0,NSim,Nr); S2 = S2(:,2:end);   % paths
    g2 = payoff(S2);                         % payoff

    exercise = Nr * ones(NSim,1);            % exercise per path
    for i=1:1:Nr-1
        i_nexercised = exercise == Nr;
        I_nexercise = iVec(i_nexercised);

        s = S2(i_nexercised,i);
        Acell = B(s);
        A = cell2mat(Acell{:,:});
        c = A * beta(:,i);

        i_exercise = g2(i_nexercised,i) >= c & g2(i_nexercised,i) > 0;
        exercise(I_nexercise(i_exercise)) = i;
    end
    summe=0;
    for j=1:1:NSim
        summe = summe + g2(j,exercise(j)) * prod(df(1:exercise(j)));
    end
    y = summe / NSim;                         % MC value from subsimulation
end


text:
function P = perms_(x)
%PERMS_       Functionality like MATLAB\PERMS, but fast
%
%   P = perms_(x)
%
% permuations of x in array of size n! x n for length(x)=n
%
%Sample timing comparing perms_(1:n) with built-in Matlab routine perms(1:n)
%in seconds on 120 MHz Pentium Laptop:
%
%
%      n     t(perms)  t(perms_)
% --------------------------------
%   2.0000    0.0052    0.0044
%   3.0000    0.0198    0.0091
%   4.0000    0.0875    0.0154
%   5.0000    0.4390    0.0253
%   6.0000    2.6650    0.0495
%   7.0000   18.8400    0.1870
%   8.0000  167.4100    1.3400
%
%
%Matlab description:
%   PERMS(1:N), or PERMS(V) where V is a vector of length N, creates a
%   matrix with N! rows and N columns containing all possible
%   permutations of the N elements.
%
%   This function is only practical for situations where N is less
%   than about 15.
%
%   See also NCHOOSEK, RANDPERM, PERMUTE.

% written  03/19/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 05/30/07     S.M. Rump  typo
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  n=length(x);
  if n==0
    P = [];
  else
    P = zeros(prod(1:n),n);
    P(1,1)=1;
    fac = 1;
    for k=2:n
      Pold = P(1:fac,1:k-1);
      w = 1:k;
      v = ((k-1)*fac+1):(k*fac);
      for i=k:-1:1
        P(v,w) = [ i*ones(fac,1) Pold+(Pold>=i) ];
        v = v-fac;
      end
      fac = k*fac;
    end
    P = x(P);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function acclaimPlayFile(fileNameAsf, fileNameAmc, frameLength)

% ACCLAIMPLAYFILE Play motion capture data from a asf and amc file.
% FORMAT
% DESC plays motion capture data by reading in a asf and amc files from disk.
% ARG fileNameAsf : the name of the ASF file to read in.
% ARG fileNameAmc : the name of the AMC file to read in.
% ARG frameLength : the length of the frames.
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% SEEALSO : acclaimPlayFile, bvhReadFile, skelPlayData

% MOCAP

skel = acclaimReadSkel(fileNameAsf);
[channels, skel] = acclaimLoadChannels(fileNameAmc, skel);
skelPlayData(skel, channels, frameLength);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function acq = acqf_vbmc(Xs,vp,gp,optimState,fmu,fs2,fbar,vtot)
%ACQF_VBMC Acquisition fcn. for prospective uncertainty search.

% Xs is in *transformed* coordinates

% Probability density of variational posterior at test points
p = max(vbmc_pdf(vp,Xs,0),realmin);

% Prospective uncertainty search
z = optimState.ymax;
acq = -vtot .* exp(fbar-z) .* p;

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [c,Ls,g,shift,M] = erblett(f,bins,fs,varargin)
%ERBLETT  ERBlet non-stationary Gabor filterbank
%   Usage: [c,Ls,g,shift,M] = erblett(f,bins,fs,varargin)
%          [c,Ls,g,shift] = erblett(...)
%          [c,Ls] = erblett(...)
%          c = erblett(...)
%
%   Input parameters:
%         f         : The signal to be analyzed (For multichannel
%                     signals, input should be a matrix which each
%                     column storing a channel of the signal)
%         bins      : Desired bins per ERB
%         fs        : Sampling rate of f (in Hz)
%         varargin  : Optional input pairs (see table below)
%   Output parameters:
%         c         : Transform coefficients (matrix or cell array)
%         Ls        : Original signal length (in samples)
%         g         : Cell array of Fourier transforms of the analysis
%                     windows
%         shift     : Vector of frequency shifts
%         M         : Number of time channels
%
%   This function computes an ERBlet constant-Q transform via non-stationary
%   Gabor filterbanks. Given the signal *f*, the ERBlet parameter *bins*,
%   as well as the sampling rate *fs* of *f*, the corresponding ERBlet
%   coefficients *c* are given as output. For reconstruction, the length of
%   *f* and the filterbank parameters can be returned also.
%
%   The transform produces phase-locked coefficients in the
%   sense that each filter is considered to be centered at
%   0 and the signal itself is modulated accordingly.
%
%   Optional input arguments arguments can be supplied like this::
%
%       erblett(f,bins,fs,'Qvar',Qvar)
%
%   The arguments must be character strings followed by an
%   argument:
%
%     'Qvar',Qvar              Bandwidth variation factor
%
%     'M_fac',M_fac            Number of time channels are rounded to
%                              multiples of this
%
%     'winfun',winfun          Filter prototype (see |firwin| for available
%                              filters)
%
%   Examples:
%   ---------
%
%   The following example shows analysis and synthesis with |erblett| and
%   |ierblett|:::
%
%       [f,fs] = gspi;
%       binsPerERB = 4;
%       [c,Ls,g,shift,M] = erblett(f,binsPerERB,fs);
%       fr = ierblett(c,g,shift,Ls);
%       rel_err = norm(f-fr)/norm(f)
%       plotfilterbank(c,Ls./M,[],fs,'dynrange',60);
%
%   See also:  ierblett, firwin
%
%   References:  ltfatnote027

% Authors: Thibaud Necciari, Nicki Holighaus
% Date: 10.04.13

warning(['LTFAT: ERBLETT has been deprecated and will be removed',...
         ' in the future releases, please use AUDFILTERS with ''erb'' and FILTERBANK instead.']);

%% Check input arguments
if nargin < 3
    error('Not enough input arguments');
end

[f,Ls,W]=comp_sigreshape_pre(f,upper(mfilename),0);

% Set defaults

definput.keyvals.usrM = [];
definput.keyvals.Qvar = 1;
definput.keyvals.M_fac = 1;
definput.keyvals.winfun = 'nuttall';

% Check input arguments

[flags,keyvals,usrM]=ltfatarghelper({'usrM'},definput,varargin);

%% Create the ERBlet dictionary

df = fs/Ls; % frequency resolution in the FFT

fmin = 0;
fmax = fs/2;

% Convert fmin and fmax into ERB
erblims = freqtoerb([fmin,fmax]);

% Determine number of freq. channels
Nf = bins*ceil(erblims(2)-erblims(1));

% Determine center frequencies
fc = erbspace(fmin,fmax,Nf)';

% Concatenate "virtual" frequency positions of negative-frequency windows
fc = [fc ; flipud(fc(1:end-1))];

gamma = audfiltbw(fc); % ERB scale

% Convert center frequencies in Hz into samples

posit = round(fc/df);% Positions of center frequencies in samples
posit(Nf+1:end) = Ls-posit(Nf+1:end);% Extension to negative freq.

% Compute desired essential (Gaussian) support for each filter
Lwin = 4*round(gamma/df);

% Nuttall windows are slightly broader than Gaussians, this is offset by
% the factor 1.1

M = round(keyvals.Qvar*Lwin/1.1);

% Compute cell array of analysis filters
g = arrayfun(@(x) firwin(keyvals.winfun,x)/sqrt(x),M,'UniformOutput',0);

g{1}=1/sqrt(2)*g{1};
g{end}=1/sqrt(2)*g{end};

M = keyvals.M_fac*ceil(M/keyvals.M_fac);
N = length(posit);  % The number of frequency channels

if ~isempty(usrM)
    if numel(usrM) == 1
        M = usrM*ones(N,1);
    else
        M = usrM;
    end
end

%% The ERBlet transform

% some preparation

f = fft(f);

c=cell(N,1); % Initialisation of the result

% The actual transform

for ii = 1:N
    Lg = length(g{ii});

    idx = [ceil(Lg/2)+1:Lg,1:ceil(Lg/2)];
    win_range = mod(posit(ii)+(-floor(Lg/2):ceil(Lg/2)-1),Ls)+1;

    if M(ii) < Lg % if the number of frequency channels is too small,
        % aliasing is introduced
        col = ceil(Lg/M(ii));
        temp = zeros(col*M(ii),W,assert_classname(f));

        temp([end-floor(Lg/2)+1:end,1:ceil(Lg/2)],:) = ...
            bsxfun(@times,f(win_range,:),g{ii}(idx));
        temp = reshape(temp,M(ii),col,W);

        c{ii} = squeeze(ifft(sum(temp,2)));

        % Using c = cellfun(@(x) squeeze(ifft(x)),c,'UniformOutput',0);
        % outside the loop instead does not provide speedup; instead it is
        % slower in most cases.
    else
        temp = zeros(M(ii),W,assert_classname(f));
        temp([end-floor(Lg/2)+1:end,1:ceil(Lg/2)],:) = ...
            bsxfun(@times,f(win_range,:),g{ii}(idx));

        c{ii} = ifft(temp);
    end
end

if max(M) == min(M)
    c = cell2mat(c);
    c = reshape(c,M(1),N,W);
end

if nargout > 3
    shift = [Ls-posit(end); diff(posit)];% Frequency hop sizes in samples
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function y = fun(x)

y=-20*exp(-0.2*sqrt((x(1)^2+x(2)^2)/2))-exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)+20+2.71289;

author: HuangCongQing


text:
function [cMap]=che(varargin)

% function [cMap]=che(n)
% ------------------------------------------------------------------------
% Creates the colormap data for n levels for the che colormap.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2016/04/14
%------------------------------------------------------------------------

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

%% Define colors
cMap=[122 141 155; 0 49 79; 214 20 28; 254 229 163; ]./255;

%% Resample colormap
[cMap]=resampleColormap(cMap,n);

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
heatmap = zeros(64, 64)*255;
for i = 1:64
    for j = 1:64
heatmap(i, j) = 1- heatmaps(1, i , j);

    end
end
hmo = HeatMap(heatmap);

colormap('hot')
imagesc(heatmap)
colorbar

author: Guanghan

repo: GNet-pose

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


text:
function [X,Y,polygon] = UISelect

%UISelect - Interactively select polygon zone in an existing plot.
%
%  USAGE
%
%    [X,Y,p] = UISelect
%
%    X,Y            polygon coordinates
%    p              handle for polygon object in figure
%
%  SEE
%
%    See also UIInPolygon
%

% Copyright (C) 2009-2011 by Michaël Zugaro
%
% 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.

hold on;
X = [];
Y = [];
polygon = [];
last = false;

while ~last,
	% Get next mouse click
	[x,y,button] = ginput(1);
	switch button,
		case 1,
			% Left button adds this point
			X(end+1,1) = x;
			Y(end+1,1) = y;
		case 2,
			% Middle button closes the selection
			last = true;
			if ~isempty(X),
				X(end+1,1) = X(1);
				Y(end+1,1) = Y(1);
			end
		case 3,
			% Right button cancels last point
			if ~isempty(X),
				X(end) = [];
				Y(end) = [];
			end
	end
	% Update existing polygon
	if isempty(polygon),
		polygon = plot(X,Y);
	else
		warning('off','MATLAB:hg:line:XDataAndYDataLengthsMustBeEqual');
		set(polygon,'XData',X);
		set(polygon,'YData',Y);
		warning('on','MATLAB:hg:line:XDataAndYDataLengthsMustBeEqual');
	end
end

hold off;


text:
function n = vectorNorm3d(v)
%VECTORNORM3D Norm of a 3D vector or of set of 3D vectors.
%
%   N = vectorNorm3d(V);
%   Returns the norm of vector V.
%
%   When V is a N-by-3 array, compute norm for each vector of the array.
%   Vectors are given as rows. Result is then a N-by-1 array.
%
%   NOTE: Computes only the Euclidean norm.
%
%   See also
%     vectors3d, normalizeVector3d, vectorAngle3d, hypot3
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-21
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

n = sqrt(sum(v.*v, ndims(v)));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [F,sA] = spm_MDP_log_evidence(qA,pA,rA)
% Bayesian model reduction for Dirichlet hyperparameters
% FORMAT [F,sA] = spm_MDP_log_evidence(qA,pA,rA)
%
% qA  - sufficient statistics of posterior of full model
% pA  - sufficient statistics of prior of full model
% rA  - sufficient statistics of prior of reduced model
%
% F   - free energy or (negative) log evidence of reduced model
% sA  - sufficient statistics of reduced posterior
%
% This routine computes the negative log evidence of a reduced model of a
% categorical distribution parameterised in terms of Dirichlet
% hyperparameters (i.e., concentration parameters encoding probabilities).
% It uses Bayesian model reduction to evaluate the evidence for models with
% and without a particular parameter.
%
% It is assumed that all the inputs are column vectors.
%
% A demonstration of the implicit pruning can be found at the end of this
% routine
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_log_evidence.m 7326 2018-06-06 12:16:40Z karl $

% change in free energy or log model evidence
%--------------------------------------------------------------------------
sA = qA + rA - pA;
F  = spm_betaln(qA) + spm_betaln(rA) - spm_betaln(pA) - spm_betaln(sA);

return

% notes: Illustration of synaptic homoeostasis in terms of Bayesian model
% reduction that considers a competition between two inputs:
%--------------------------------------------------------------------------
x     = linspace(1,32,128);
pA    = [1; 1];
rA    = pA;
rA(2) = 8;
for i = 1:numel(x)
    for j = 1:numel(x)
        qA = [x(i);x(j)];
        F(i,j) = spm_MDP_log_evidence(qA,pA,rA);
    end
end

subplot(2,2,1), imagesc(x,x,F + (F > 0)*4),
title('Free energy landscape','FontSize',16), axis square xy
xlabel('concentration parameter'), ylabel('concentration parameter')
subplot(2,2,2), plot(x,F'),  title('log evidence','FontSize',16)
xlabel('concentration parameter'), ylabel('Log-evidence'), axis square

author: spm


text:
function RegTranslation(handles, metric)
%
% 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/>.

global planC stateS;
indexS = planC{end};

    [originF, spacingF, centerF] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));

    [originM, spacingM, centerM] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationMovDataset));

    minstep = str2double(get(handles.para.minstep, 'string'));
    maxstep = str2double(get(handles.para.maxstep, 'string'));
    iternum = str2double(get(handles.para.iternum, 'string'));

    output = cell(1, 8);

    FImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
    MImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;

%     downsample the input datasets
    if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))

        tic;
        FImg = imdesample3d(FImg,2,2);
        MImg = imdesample3d(MImg,2,2);
        toc;
        DesampleTime = num2str(floor(toc/60));

    end

%     call registration method
    tic;
    switch (metric)
        case 'mean squares'
            [im, Rotation, Offset] = MeanSquare3D(int16(FImg), originF, spacingF, ... %MeanSquare3D,NormalizedCorrelation3D
                                                int16(MImg), originM, spacingM, ...
                                                minstep, maxstep, iternum);
        case 'normalized correlation'
            [im, Rotation, Offset] = NormalizedCorrelation3D(int16(FImg), originF, spacingF, ... %MeanSquare3D,NormalizedCorrelation3D
                                                int16(MImg), originM, spacingM, ...
                                                minstep, maxstep, iternum);
    end
    toc;
    RegTime = num2str(floor(toc/60));

    output{1} = ['Translation X = ' num2str(Offset(1))];
    output{2} = ['Translation Y = ' num2str(Offset(2))];
    output{3} = ['Translation Z = ' num2str(Offset(3))];
    output{4} = ['Number of Iterations = ' num2str(Offset(4))];
    output{5} = ['Best Value = ' num2str(Offset(5))];
    output{6} = ['Desample time = ' DesampleTime 'm'];
    output{7} = ['Register time = ' RegTime 'm'];

    set(handles.OutputList, 'string', output);

    %update the transM;
    dx = Offset(1); dy = Offset(2); dz = Offset(3);

    scanSetM = stateS.imageRegistrationMovDataset;
    oldTransM = getTransM(stateS.imageRegistrationMovDatasetType, scanSetM, planC);

    newTransform = eye(4);
    newTransform(:,4) = [dx dy dz 1];

%     planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = (newTransform * oldTransM);
%     planC{indexS.dose}(scanSetM).transM = (newTransform * oldTransM);
    planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = newTransform;
    sliceCallBack('refresh');

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function op = linop_handles( sz, Af, At, cmode )
%LINOP_HANDLES Linear operator from user-supplied function handles.
%OP = LINOP_HANDLES( SZ, AF, AT, CMODE )
%    Constructs a TFOCS-compatible linear operator from separate function
%    handles that compute the forward and adjoint operations. The first
%    argument, SZ, gives the size of the linear operator; and the forward
%    and adjoint handles are AF and AT, respectively.
%
%    If the inputs and outputs are simple vectors, then SZ can take the
%    standard Matlab form [M,N], where N is the input length and M is the
%    output length. If the input or output is a matrix or array, then SZ
%    should take the form { S_in, S_out }, where S_in is the size of the
%    input and S_out is the size of the output.
%
%    If the input or output space of the operator is complex, then the
%    CMODE string must be supplied, which describes the forward operation:
%       'R2R': real input and output
%       'R2C': real input, complex output
%       'C2R': complex input, real output
%       'C2C': complex input, complex output
%    If CMODE is not supplied, then 'R2R' is assumed. If the operator
%    detects a complex input or output when it is not expected, then an
%    error results. Therefore, you must make sure that your operators
%    return real values when they are expected to do so.

error(nargchk(3,4,nargin));
if numel( sz ) ~= 2,
    error( 'Size must have two elements.' );
elseif isnumeric( sz ),
    sz = { [sz(2),1], [sz(1),1] };
elseif ~isa( sz, 'cell' ),
    error( 'Invalid operator size specification.' );
end
if ~isa( Af, 'function_handle' ) || ~isa( At, 'function_handle' ),
    error( 'Second and third arguments must be function handles.' );
end
if nargin < 4 || isempty( cmode ),
    cmode = 'R2R';
elseif ~ischar( cmode ) || size( cmode, 1 ) ~= 1,
    error( 'Fourth argument must be a string.' );
else
    cmode = upper( cmode );
end

switch cmode,
    case 'R2R', op = @(x,mode)linop_handles_r2r( sz, Af, At, x, mode );
    case {'R2C','R2CC'}, op = @(x,mode)linop_handles_r2c( sz, Af, At, x, mode );
    case {'C2R','CC2R'}, op = @(x,mode)linop_handles_c2r( sz, Af, At, x, mode );
    case 'C2C', op = @(x,mode)linop_handles_c2c( sz, Af, At, x, mode );
    otherwise,
        error( 'Invalid complex mode: %s', cmode );
end

function y = linop_handles_r2r(sz, Af, At, x, mode )
switch mode,
    case 0, y = sz;
    case 1, y = realcheck( Af( realcheck( x ) ) );
    case 2, y = realcheck( At( realcheck( x ) ) );
end

function y = linop_handles_r2c(sz, Af, At, x, mode )
switch mode,
    case 0, y = sz;
    case 1, y = Af( realcheck( x ) );
    case 2, y = realcheck( At( x ) );
end

function y = linop_handles_c2r(sz, Af, At, x, mode )
switch mode,
    case 0, y = sz;
    case 1, y = realcheck( Af( x ) );
    case 2, y = At( realcheck( x ) );
end

function y = linop_handles_c2c(sz, Af, At, x, mode )
switch mode,
    case 0, y = sz;
    case 1, y = Af( x );
    case 2, y = At( x );
end

function y = realcheck( y )
if ~isreal( y ),
    error( 'Unexpected complex value in linear operation.' );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function th = threeshold_estimation(distance,method)
    if nargin < 2
        method = 'mean';
    end

    switch method
        case 'mean'
            th = mean(distance)+2.5*std(distance);
        case 'median'
            th = median(distance)+2.5*1.4826*mad(distance,1);
    end

author: alexandrebarachant

repo: covariancetoolbox


text:
classdef Tensor2VoigtPSrepresentation < handle

    properties (Access = protected)
        indexTransformer
    end

    methods (Access = protected)

        function createIndexTransformer(obj)
            obj.indexTransformer = TensorVoigtIndexTransformerPS();
        end

    end

end

author: SwanLab

repo: Swan


text:
function [Y,W,SetupStruc] = Process_LCMV(s,Transfer,SetupStruc)
K = SetupStruc.LCMV.K;
hop = SetupStruc.LCMV.hop;
win = hanning(K,'periodic');
win = win/sqrt(sum(win(1:hop:K).^2));
SetupStruc.LCMV.win = win;  % Preserve 'win' in 'SetupStruc'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(s,2);
for i = 1:N
    X(:,:,i) = fft(enframe(s(:,i),win,hop)');
end
frame_N = size(X,2);
K_m = K/2+1;
Num = size(Transfer,3);
Y = zeros((frame_N-1)*hop+K,Num);
Y_f = zeros(size(X,1),size(X,2),Num);
WNG =zeros(K_m,Num);
%%%%%%%%%%%%%%%%%%%%%%%%%% Obtain processing matrix 'W'
theta = 10^-4;
W = zeros(Num,N,K_m);
for i = 2:K_m
    X_f = permute(X(i,:,:),[3 2 1]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% LCMV processing, imposing the linear constraints to the angle of sources
    Steer = permute(Transfer(i,:,:),[2 3 1]);
    R = X_f*X_f'/frame_N;
    if rcond(R)<theta
        R = R+eye(N)*min(diag(R))*theta;
    end
    R_denominator = Steer'/R*Steer;
    if rcond(R_denominator)<theta
        R_denominator = R_denominator+eye(Num)*min(diag(R_denominator))*theta;
    end
    W_f = R\Steer/R_denominator;
    W_f = W_f';
    W(:,:,i) = W_f;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    WNG(i,:) = -10*log10(diag(W_f*W_f')');
    Y_ = W_f*X_f;
    Y_f(i,:,:) = Y_.';
    if(i~=K_m)
        Y_f(K+2-i,:,:) = Y_';
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Recover signals
if(K/hop==2)
    win = ones(K,1);
end
for i = 1:Num
    Y(:,i) = overlapadd(real(ifft(Y_f(:,:,i)))',win,hop);
end
% autoPlot(WNG,'LCMV');
return;

repo: Voice-Separation-and-Enhancement


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

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

% Last Modified by GUIDE v2.5 06-Jun-2012 18:22:54

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Leveling_OpeningFcn, ...
                   'gui_OutputFcn',  @Leveling_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 Leveling is made visible.
function Leveling_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 Leveling (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = Leveling_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;

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

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

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

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

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

[FileName,PathName] = uigetfile({'*.dat';'*.txt';'*.*'},'Select the data file');
if FileName~=0
  directory=strcat(PathName ,FileName);
  setappdata(Leveling,'myfile',directory);
   setappdata(Leveling,'file',FileName);
  [ben,V1]=benmark(directory);%ben=bnckmark V=benmark value
[sta1,sta2,V,std]=station(directory);%sta1 and sta1 -stn1 and 2 V=value std=standard deviation
n=length(V1);
fid=fopen('Info','w');
fprintf(fid,'Benchmark\n');
fprintf(fid,'station    Elevation value\n');
num=ones(n,1);
for i=1 : n
    fprintf(fid,'%s           %.4f\n',ben{1,i},V1(i,1));
end
fprintf(fid,'********************************\n');
fprintf(fid,'observations\n');
fprintf(fid,'From   To     Elevation value    (+/-)std\n');
n1=length(V);
num1=ones(n1,1);
num2=ones(n1,1);
for i=1 : n1
    C1=sta1{1,i};
      num1(i,1)=unicode2native(C1);
C2=sta2{1,i};
 num2(i,1)=unicode2native(C2);
end
 for i=1 : n1
    fprintf(fid,'%s      %s       %.5f          %.4f \n',sta1{1,i},sta2{1,i},V(i,1),std(i,1));
 end
  to_tum=[num1
       num2];
 to=unique(to_tum);
 %number of station
 numT=length(to);
 %fprintf('number of station:%d\n',numT);
 %number of parameter

numA=numT-n;
TT=to;
  %fprintf('number of parameter:%d\n',numA);
  for i=1 :  n
      for j=1 : numT
           if num(i,1)==TT(j,1)
              TT(j,1)=0;
           end

          %fprintf('%d : %d \n',i,j)
      end
  end
  K=0;
  para=ones(numA,1);
  for j=1 : numT
  if TT(j,1) ~=0
      K=K+1;
      para(K)=TT(j,1);
  end
  end

  %A matrrik

  A=zeros(n1,numA);
  for i=1 : numA
      for j=1 : n1
           if para(i,1)==num1(j,1)
               A(j,i)=1;

           end
          %fprintf('%d : %d \n',i,j)
      end
  end
  for i=1 : numA
      for j=1 : n1
           if para(i,1)==num2(j,1)
               A(j,i)=-1;
           end
          %fprintf('%d : %d \n',i,j)
      end
  end
 %B matric
B=zeros(n1,1);
for i=1 : n
    for j=1 : n1
        if num(i,1)==num1(j,1)
            B(j,1)=V1(i,1);
        end
    end
end
 for i=1 : n
    for j=1 : n1
        if num(i,1)==num2(j,1)
            B(j,1)=-V1(i,1);
        end
    end
end
%fprintf('%.4f\n',B)
fprintf(fid,'number of benchmark:%d\n',n);
fprintf(fid,'number of station:%d\n',numT);
fprintf(fid,'number of observation:%d\n',n1);
fprintf(fid,'number of parameter:%d\n',numA);
fclose(fid);
fid1=fopen('Info');
indic = 1;
    while 1
         tline = fgetl(fid1);
         if ~ischar(tline),
             break
         end
         strings{indic}=tline;
         indic = indic + 1;
    end
    fclose(fid1);
set(handles.edit1,'String', strings)
end

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
fc=getappdata(Leveling,'myfile');
fN=getappdata(Leveling,'file');
[a1,a2]=size(fN);
%disp(a2)
if a2 > 0
ind=findstr(fN,'.');
%disp(fN(1:ind-1))
output=strcat(fN(1:ind-1) ,'.out') ;
 adj(fc,output)
 fid1=fopen(output);
indic = 1;
    while 1
         tline = fgetl(fid1);
         if ~ischar(tline),
             break
         end
         strings{indic}=tline;
         indic = indic + 1;
    end
    fclose(fid1);
set(handles.edit2,'String', strings)
set(handles.pushbutton3,'Enable','on');
end

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

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

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

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

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.edit1,'String', '')
set(handles.edit2,'String', '')

clear fN
set(handles.pushbutton3,'Enable','off');

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

delete(Leveling)

author: Sable

repo: mcbench-benchmarks


text:
%%*************************************************************
%% Asum: compute the matrix
%%
%%  Ay = Asum(blk,A,y)
%%
%%  input:  A = a CELL ARRAY with m columns.
%%          y = mx1 vector.
%%          permAy = a permutation of [1:m] coding the order to
%%                   sum the matrices y(k)*Ak, k = 1,...m.
%%          iscmp = 1,  if Ay is complex
%%                  0,  otherwise.
%%
%%  output: Ay = sum_{k=1}^m y(k)*Ak, a column CELL ARRAY
%%               with the same structure as A{:,1}.
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 2 Feb 01
%%*************************************************************

   function Ay = Asum(blk,A,y)

   Ay = cell(size(blk,1),1);
   m = length(y);
   for p = 1:size(blk,1)
      n = sum(blk{p,2});
      blktmp = blk{p,2};
      if strcmp(blk{p,1},'s');
         tmp = sparse(n,n);
         for k = 1:m; tmp = tmp + A{p,k}*y(k); end;
         if (length(blktmp) == 1)
            if (nnz(tmp) > 0.15*n*n);
               if issparse(tmp); tmp = full(tmp); end;
            else;
               if ~issparse(tmp); tmp = sparse(tmp); end;
            end;
         elseif (length(blktmp) > 1);
            if ~issparse(tmp); tmp = sparse(tmp); end;
         end;
      elseif strcmp(blk{p,1},'l');
         tmp = zeros(n,1);
         for k = 1:m; tmp = tmp + A{p,k}*y(k); end;
         if (nnz(tmp) > 0.15*n);
            if issparse(tmp); tmp = full(tmp); end;
         else;
            if ~issparse(tmp); tmp = sparse(tmp); end;
         end;
      end;
      Ay{p} = tmp;
   end;
%%-------------------------------------------------------------

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function obj=MOcovMFileCollection(root_dir, method, monitor, exclude_pat)
% instantiate MOcovMFileCollection
%
% obj=MOcovMFileCollection(root_dir, method, monitor)
%
% Inputs:
%   root_dir                root directory containing m-files to be
%                           covered.
%   method                  Coverage method, one of:
%                           - 'file'    rewrite m-files after adding
%                                       statements that record which lines
%                                       are covered
%                           - 'profile' use Matlab profiler
%                           default: 'file'
%   monitor                 optional MOcovProgressMonitor instance
%   exclude_pat             Optional cell array of patterns to exclude.
%
% See also: mocov

    if nargin<4 || isempty(exclude_pat)
        exclude_pat={};
    end

    if nargin<3 || isempty(monitor)
        monitor=MOcovProgressMonitor();
    end

    if nargin<2 || isempty(method);
        method='file';
    end

    props=struct();
    props.root_dir=root_dir;
    props.monitor=monitor;
    props.exclude_pat=exclude_pat;
    props.mfiles=[];
    props.orig_path=[];
    props.temp_dir=[];
    props.method=method;
    obj=class(props,'MOcovMFileCollection');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [ FV ] = mesh_refine(FV,Nface)

% mesh_refine - creates smaller triangles from a triangle mesh
%
% [ FV ] = mesh_refine( FV, Nface )
%
% FV.vertices   - vertex matrix (Nx3)
% FV.faces      - face matrix (Mx3), indices into vertex matrix rows
%
% Nface         - subdivide faces into 4 or 6 faces,
%                 the default 4 provides an even subdivision
%
% This function calls mesh_refine_tri4 or mesh_refine_tri6.  See
% these for more details.
%

% This can be done until some minimal distance (D) of the mean
% distance between vertices of all triangles is achieved.  If
% no D argument is given, the function refines the mesh once.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  08/2002, Darren.Weber_at_radiology.ucsf.edu, created
%                    adapted this function as a wrapper to
%                    mesh_refine_tri4 & mesh_refine_tri6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('FV','var'),
    error('MESH_REFINE: NO input FV struct');
elseif isempty(FV),
    error('MESH_REFINE: NO input FV struct');
end

if ~exist('Nface','var'),
    Nface = 4;
elseif isempty(Nface),
    Nface = 4;
end

switch Nface,

case 4,
    FV = mesh_refine_tri4(FV);
case 6,
    FV = mesh_refine_tri6(FV);
otherwise
    FV = mesh_refine_tri4(FV);
end

return

author: PatternRecognition

repo: OpenBMI


text:
function J = GenToneMap(im)
% ==============================================
%   Compute the tone map 'T'
%
%   Paras:
%   @im        : input image ranging value from 0 to 1.
%

    %% Parameters
    Ub = 225;
    Ua = 105;

    Mud = 90;

    DeltaB = 9;
    DeltaD = 11;

    % groups from dark to light
    % 1st group
    Omega1 = 42;
    Omega2 = 29;
    Omega3 = 29;
    % 2nd group
    Omega1 = 52;
    Omega2 = 37;
    Omega3 = 11;
    % 3rd group
    Omega1 = 76;
    Omega2 = 22;
    Omega3 = 2;

    %% Compute the target histgram
    histgramTarget = zeros(256, 1);
    total = 0;
    for ii = 0 : 255
        if ii < Ua || ii > Ub
            p = 0;
        else
            p = 1 / (Ub - Ua);
        end

        histgramTarget(ii+1, 1) = (...
            Omega1 * 1/DeltaB * exp(-(255-ii)/DeltaB) + ...
            Omega2 * p + ...
            Omega3 * 1/sqrt(2 * pi * DeltaD) * exp(-(ii-Mud)^2/(2*DeltaD^2))) * 0.01;

        total = total + histgramTarget(ii+1, 1);
    end
    histgramTarget(:, 1) = histgramTarget(:, 1)/total;

%     %% Smoothing
%     im = medfilt2(im, [5 5]);

    %% Histgram matching
    J = histeq(im, histgramTarget);

    %% Smoothing
    G = fspecial('average', 10);
    J = imfilter(J, G,'same');
end

repo: PencilDrawing


text:
function y = fcube ( x )

%*****************************************************************************80
%
%% FCUBE sets the cubic data values.
%
%  Discussion:
%
%    Y(X) = ( ( X + 2 ) * X + 3 ) * X + 4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real Y, the value of the function.
%
  y = ( ( (       1.0E+00 ) ...
                * x + 2.0E+00 ) ...
                * x + 3.0E+00 ) ...
                * x + 4.0E+00;

  return
end

author: johannesgerer

repo: jburkardt-m


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function  [Jc,para,dJ,H] = LDDMMobjFctn(T,Rc,omega,m,beta,M,wRef,xc,vc)
%
% Objective Function for LDDMM using Lagrangian PDE solver
%
% computes J(vc) = D(T(yc(vc)),Rc) + S(vc - vRef), where
%
% vc       - is a velocity field (stationary or instationary)
% yc(vc)   = trafo(wc,xc) is obtained by tracing characteristics using RK4
%            scheme
% Tc       = T(yc) = imgModel(T,omega,yc),
% D(Tc,Rc) = distance(Tc,Rc,omega,m)
% S        = regularizer, e.g., S(uc) = 0.5*uc'*B*uc
% vRef     = reference guess for velocities
%
% For more details see the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
%
% Input:
%   T      - data for template image, Tc = imgModel(T,omega,yc)
%   Rc     - reference image on grid, Rc = imgModel(R,omega,xc)
%   omega  - representation of computational domain
%   m      - discretization size
%   vRef   - reference velocity
%   xc     - discretization of Omega
%   omegaV - computational domain for velocity field (can be larger)
%   mV     - discretization size for velocities
%   N      - number of time steps for RK4 scheme
%   vc     -  current velocities
%
% Output:
%  Jc      - current function value J(vc)
%  para    - struct {Tc=T(y(vc)), Rc, omega, m, yc=y(vc,xc), Jc}, for plots
%  dJ      - gradient of J
%  H       - approximation to Hessian of J
%
% see also ELDDMM_Hand2D
%==============================================================================
function [Jc,para,dJ,H] = LDDMMobjFctn(T,Rc,omega,m,vRef,xc,omegaV,mV,N,vc)

para = struct([]);
if nargin == 0,
  help(mfilename);
  runMinimalExample;
  return;
elseif ~exist('vc','var') || isempty(vc),
  % if wc is not an input argument, reports status
  if nargout == 1, Jc = 'LDDMM';  return; end;
  % report current settings
  dimstr  = @(m) sprintf('[%s]',sprintf(' %d',m));
  vc      = trafo('w0');
  nt      = regularizer('get','nt');
  fprintf('Large Deformation Diffeomorphic Mapping (LDDMM):\n');
  fprintf('   J(vc) = D(T(yc(vc)),Rc) + S(vc-vRef) != min\n');
  fprintf('  %20s : %s\n','m',dimstr(m));
  fprintf('  %20s : %s\n','omega',dimstr(omega));
  fprintf('  %20s : %s\n','IMAGE MODEL',imgModel);
  fprintf('  %20s : %s\n','DISTANCE',distance);
  fprintf('  %20s : %s\n','TRAFO',trafo);
  fprintf('  %20s : %s\n','#timeSteps',num2str(N));
  fprintf('  %20s : %s\n','nt',num2str(nt));
  fprintf('  %20s : %s\n','mV',dimstr(mV));
  fprintf('  %20s : %s\n','omegaV',dimstr(omegaV));
  fprintf('  %20s : %s\n','length(vc)',num2str(length(vc)));
  Jc = vc; % return starting guess
  return;
end;
tspan = [1 0];
dim   = numel(omega)/2;
nt    = round(numel(vc)/(dim*prod(m)))-1;

% do the work ------------------------------------------------------------
matrixFree   = regularizer('get','matrixFree');
doDerivative = (nargout>2);            % flag for necessity of derivatives

% compute transformation, distance, and regularization and combine these
if nt<1
    [yc,dy,pTrafo] = getTrafoFromVelocityRK4(vc,xc,'omega',omegaV,'m',mV,'N',N,'tspan',tspan,'doDerivative',doDerivative);
else
    [yc,dy,pTrafo] = getTrafoFromInstationaryVelocityRK4(vc,xc,omegaV,[mV,nt],'N',N,'tspan',tspan,'doDerivative',doDerivative);
end
[Tc,dT] = imgModel(T,omega,center(yc,m),'doDerivative',doDerivative);

% compute distance
[Dc,~,dD,dres,d2psi] = distance(Tc,Rc,omega,m,'doDerivative',doDerivative);

% compute regularizer
[Sc,dS,d2S] = regularizer(vc-vRef,omegaV,mV,'doDerivative',doDerivative,'tspan',tspan);

Jc = Dc + Sc;

% collect variables for plots
para = struct('Tc',Tc,'Rc',Rc,'omega',omega,'m',m,'yc',yc,'Jc',Jc,'Sc',Sc,'Dc',Dc,...
    'N',pTrafo.N);

if ~doDerivative, return; end;
dD = dD*dT*dy;
dJ = dD + dS;
if nargout<4, return; end;

% approximation to Hessian
dres = dres*dT*dy;

% multiply outer and inner derivatives, note: dy might be sparse
if not(matrixFree),
  H  = dres'*d2psi*dres + d2S;
else
    % approximation to d2D in matrix free mode
    % d2D   = dr'*d2psi*dr
    % P and P' are operators matrix free
    H.omega     = omegaV;
    H.m         = mV;
    H.nt        = nt;
    H.tspan     = tspan;
    H.d2D.how   = '*dr''*d2psi*dr';
    H.d2D.P     = @(x) x;
    H.d2D.dr    = dres;
    H.d2D.d2psi = d2psi;

    H.d2S = d2S;

end;

function runMinimalExample
setup2DGaussianData;
lvl = 5;
regularizer('reset','regularizer','mbElastic','alpha',1)
v0 = getVelocityStartingGuess(omega,m);
xc = getCellCenteredGrid(omega,ML{lvl}.m);
T  = reshape(imgModel(ML{lvl}.T,omega,center(xc,ML{lvl}.m)),ML{lvl}.m);
Rc = imgModel(ML{lvl}.R,omega,center(xc,ML{lvl}.m));

fctn = @(vc) LDDMMobjFctn(T,Rc,omega,ML{lvl}.m,v0,xc,omega,m,10,vc);
checkDerivative(fctn,v0)


text:
% DEMUNLABELLEDONEIVM2 Test IVM code on a toy crescent data.
%
% Recreates the toy crescent data example shown in the NIPS paper.

% GPMAT

randn('seed', 1e6)
rand('seed', 1e6)

% Generate a toy data-set
dataSetName = 'unlabelledOne';
experimentNo = 2;

labelledProb = 0.1;
[X, y] = mapLoadData(['semi:' dataSetName ':' num2str(labelledProb)], 1e6);

ind = find(isnan(y));
X(ind, :) = [];
y(ind, :) = [];

options = ivmOptions;
options.noise = 'probit';
options.kern = {'rbf', 'white'};
options.display = 2;
options.numActive = 100;
prior = 0;

% Initialise the model.
model = ivmCreate(size(X, 1), size(y, 2), X, y, options);

prior.type = 'gamma';
prior = priorParamInit(prior);
prior.a = 1;
prior.b = 1;
prior.index = 2;
model.kern.comp{1}.priors(1) = prior;
prior.index = 1;
model.kern.comp{2}.priors(1) = prior;
if options.display > 1
  ivm3dPlot(model, 'ncnmContour', i); %incnmTwoDPlot(model, i);
end
for i = 1:15

  % Plot the data.
  % Select the active set.
  model = ivmOptimiseIvm(model, options.display);
  if options.display > 1
    ivm3dPlot(model, 'ncnmContour', i); %incnmTwoDPlot(model, i);
  end
  % Optimise the kernel parameters.
  model = ivmOptimiseKernel(model, options.display, options.kernIters);
  ivmDisplay(model);

end
model = ivmOptimiseIvm(model, options.display);
if options.display > 1
  ivm3dPlot(model, 'ncnmContour', i); %incnmTwoDPlot(model, i);
end
model = ivmOptimiseIvm(model, options.display);
% Display the final model.
ivmDisplay(model);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
[kern, noise, ivmInfo] = ivmDeconstruct(model);
save(['dem' capName num2str(experimentNo) '.mat'], ...
     'kern', ...
     'noise', ...
     'ivmInfo');

if exist('printDiagram') & printDiagram
  ivmPrintPlot(model, 'ivmContour', [], [], [], capName, experimentNo);
end

author: SheffieldML

repo: GPmat


text:
function mastercorr_plot_stats(W)

%MASTERCORR_PLOT_STATS plot statistics from master correlation algorithm
% W = MASTERCORR_PLOT_STATS(W) plots summary statistics following an
% application of MASTERCORR_SCAN. This summary is useful for assessing the
% quality and the nature of the detected events. The input waveform W may
% be any dimension. However it must contain the fields produced by
% MASTERCORR_SCAN including: MASTERCORR_TRIG, MASTERCORR_CORR,
% MASTERCORR_ADJACENT_CORR
%
% *** NOTE ABOUT MULTIPLE WAVEFORMS ***
% This function is designed to accept NxM waveform matrices as input. In
% this case the plot will contain relevent data from all element waveforms
% of W. This is useful, for example, when W is a 24x1
% matrix of hourly waveforms. However, unexpected (or clever!) results may
% be produced when W is complicated by elements with different channels or
% master waveform snippets. For some uses it may prove wise to pass only
% selected elements of W to this function. For example:
% C = MASTERCORR_EXTRACT(W(1:5))
%
% See also mastercorr_scan, mastercorr_extract

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

% CHECK INPUTS
if nargin ~= 1
    error('Incorrect number of inputs');
end
if ~isa(W,'waveform')
    error('First argument must be waveform objects');
end

% READ MASTERCORR FIELDS
trig = [];
corr = [];
corrAdj = [];
if numel(W)>1
    for n=1:numel(W)
        trig = [trig ; get(W(n),'MASTERCORR_TRIG')];
        corr = [corr ; get(W(n),'MASTERCORR_CORR')];
        corrAdj = [corrAdj ; get(W(n),'MASTERCORR_ADJACENT_CORR')];

    end
else
    trig = get(W,'MASTERCORR_TRIG');
    corr = get(W,'MASTERCORR_CORR');
    corrAdj = get(W,'MASTERCORR_ADJACENT_CORR');
end

if numel(trig)==0
   disp('No MASTERCORR_TRIG field in this data');
   return;
end

% SORT BY TIME
[tmp,index] = sort(trig);
trig = trig(index);
corr = corr(index);
corrAdj = corrAdj(index);

% MAKE PLOT
figure('Color','w','Position',[50 50 850 1100]);
box on; hold on;
set(gcf,'DefaultLineLineWidth',1);
set(gcf,'DefaultAxesFontSize',12);
set(gcf,'DefaultLineMarkerSize',5);

% TIME HISTORY
subplot(3,1,1);
plot([trig' ; trig'],[corr' ; corrAdj'],'-','Color',[0.7 0.7 0.7])
hold on;
plot(trig,corr,'ok','MarkerFaceColor','y')
xlim([min(trig) max(trig)]);
ylim([0.5 1.0]);
datetick('x','KeepLimits');
set(gca,'YGrid','on')
set(gca,'YTick',[0:.1:1]);
ylabel('Correlation coefficient')
legend('adjacent peak','location','SouthWest');

subplot(3,1,2);
plot(trig(2:end),86400*(trig(2:end)-trig(1:end-1)),'ko','MarkerFaceColor','y')
set(gca,'YScale','log');
xlim([min(trig) max(trig)]);
datetick('x','KeepLimits');
set(gca,'YGrid','on')
%set(gca,'YTick',10.^[1:10]);
ylabel('Time (s)')
legend('Event spacing','location','SouthWest');

subplot(3,1,3)
hist(corr,[0.5:0.01:1]);
xlim([0.5 1]);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','y','EdgeColor','k')
legend('Correlation','location','SouthWest');

% PLOT INTEREVENT TIME AGAINST CORRELATION VALUE
%figure
%tDiff = 86400*(trig(2:end)-trig(1:end-1));
%plot([tDiff' ; tDiff'],[corr(2:end)' ; corr(1:end-1)'],'k-')
%hold on;
%plot(tDiff,max([corr(2:end)' ; corr(1:end-1)']),'ko');
%set(gca,'XScale','log');

%PRINT OUT FIGURE
set(gcf, 'paperorientation', 'portrait');
set(gcf, 'paperposition', [.25 .25 8 10.5] );

author: geoscience-community-codes

repo: GISMO


text:
function[varargout]=cum2mom(varargin)
%CUM2MOM  Convert cumulants to moments.
%
%   [M0,M1,...MN]=CUM2MOM(K0,K1,...KN) converts the first N cumulants
%   K0,K1,...KN into the first N moments M0,M1,...MN.
%
%   The KN and MN are all scalars or arrays of the same size.
%
%   Note for a probability density function, M0=1 and K0=0.
%
%   MCELL=CUM2MOM(KCELL), where KCELL is a cell array whose (N+1)th
%   is the Nth cumulant, returns a similar cell array of moments.
%
%   CUM2MOM is inverted by MOM2CUM.
%
%   See also MOM2CUM, BELLPOLY.
%
%   Usage: [m0,m1,m2]=cum2mom(k0,k1,k2);
%          mcell=cum2mom(kcell);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008 J.M. Lilly --- type 'help jlab_license' for details

%Tests for CUM2MOM are run by MOM2CUM.
if strcmpi(varargin{1}, '--t'),return,end

if nargin==1&&iscell(varargin{1})
    cum=varargin{1};
else
    cum=varargin;
end
mom=cell(size(cum));
mom{1}=exp(cum{1});

cum=cum(2:end);
mom(2:end)=bellpoly(cum);

for i=2:length(mom)
    mom{i}=mom{i}.*mom{1};
end

if nargin==1&&iscell(varargin{1})
    varargout{1}=mom;
else
    varargout=mom;
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab