text:
function coeffVar = compute_coeff_var(this, varargin)
% computes standard deviation image (coeff_var) over 4th dimension of MrImage
% NOTE: short-cut for compute_stat_image('coeff_var')
%
%   Y = MrImage()
%   coeffVar = Y.compute_coeff_var('PropertyName', PropertyValue)
%
% This is a method of class MrImage.
%
% IN
%   'PropertyName'
%               'selectedVolumes'       [1,nVols] vector of selected
%                                       volumes for statistical calculation
% OUT
%   coeffVar         MrImage holding voxel-wise coefficient of variation
%                     image (coeff_var), i.e. 1./snr
%                     (with thresholding to avoid Inf-values)
%
% EXAMPLE
%   Y = MrImage()
%   coeffVar = Y.compute_coeff_var('selectedVolumes', [6:100])
%
%   See also MrImage compute_stat_image

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

coeffVar = this.compute_stat_image('coeff_var', varargin{:});

author: translationalneuromodeling

repo: tapas


text:
%       Created by Dwight Nwaigwe  April 2011
%	This program uses the ensemble kalman filter to estimate a system's state.
%	The state is x_new=f(x,u)+w, where u some input, w the
%	Gaussian distributed process noise, and f is a nonlinear function. The measurement
%	is y_new=h(x)+v where h is a nonlinear function and v Gaussian distributed measurement noise.

%       The algorithm used in this code is referenced from the following:
%       S Gillijns et. al., "What Is the Ensemble Kalman Filter and How Well Does it Work?"
%       Proceedings of the 2006 American Control Conference,
%       Minneapolis, Minnesota, USA, June 14-16, 2006, pp 4448-4453.

function [x_tr,x_estbar,ybar]= ensemblekfilter(f,h,x_tr,x_ini,w,z,num_iterations)

%      Example

%      The state and measurement in this example  is taken from Dan Simon, "Kalman Filtering",
%      Embedded Systems Programming,2001.
%
%
%	syms  x1 x2;   %variables must be named x1...xn
%	f=[x1+.1*x2+.005;x2+.1];
%	h=[x1];
%	x_tr=[1;1]; %initial value of state
%	x_ini=ones(2,20); %ensemble of initial estimate of the state
%	w=[10^-3; .02];  %process noise standard deviation
%	z=[10];  %measurement noise standard deviation
%	num_iterations=600;
% 	num_members=20;
% 	[a,b,c]=ensemblekfilter(f,h,x_tr,x_ini,w,z,num_iterations);

[dummy,num_members]=size(x_ini);
p1=length(f);
m1=length(h);
var_vector=[];
xvec=[];
x_estvec=[];
yvec=[];
x_est=x_ini;

for j=1:p1  %create vector containing variables x1 to xn
  eval(sprintf(' syms x%d', j));
  var_vector=[var_vector sprintf('x%d ',j)];
  end
var_vector=strcat('[',var_vector);
var_vector=strcat(var_vector,']');

Zcov=eye(m1); %create measurement noise covariance matrix
for j=1:m1
  Zcov(j,j)=z(j)^2;
end

for i=1:num_iterations

   x_tr=subs(f,var_vector,x_tr)+w.*randn(p1,1); %compute true value of state at next time step

   for j=1:num_members
     W(:,j)=w.*randn(p1,1);                          %create process noise
     Z(:,j)=z.*randn(m1,1);                          %create measurement noise
     x_est(:,j)=subs(f,var_vector,x_est(:,j))+W(:,j);      %forecast state
     y(:,j)=subs(h,var_vector,x_tr)+Z(:,j);                 %make measurement
     y_for(:,j)=subs(h,var_vector,x_est(:,j));              %forecast measurement
   end

   x_estbar=mean(x_est,2);
   ybar=mean(y,2);
   y_forbar=mean(y_for,2);

   for j=1:p1
     Ex(j,:)=[x_est(j,:)-x_estbar(j)];
   end

   for j=1:m1
     Ey(j,:)=[y_for(j,:)-y_forbar(j)];
   end

   Pxy=Ex*Ey'/(num_members-1);
   Pyy=Ey*Ey'/(num_members-1)+Zcov;                     %The addition of Zcov to Pyy is not done in Gillijns et. al but I use it here in case num_members=2 or Pyy is nearly singular
   K=Pxy*inv(Pyy);
   x_est=x_est+K*(y-y_for);
   xvec=[xvec x_tr];
   x_estvec=[x_estvec x_estbar];
   yvec=[yvec ybar];
   if i==num_iterations
   x_estbar=mean(x_est,2);
   end
end

author: Sable

repo: mcbench-benchmarks


text:
%% (Internal) Calculate the median of absolute deviations from the median (MEDA)
%
%   meda = nanmeda(data)
%
% Arguments:
%
%      + data: data
%
% Output:
%
%      + meda:  median of absolute deviations from the median.
%
% Example:
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function meda = nanmeda(data)

median_data = nanmedian(data);
meda = nanmedian(abs(bsxfun(@minus, data, median_data)));

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function fd1d_advection_lax ( )

%*****************************************************************************80
%
%% FD1D_ADVECTION_LAX solves the advection equation using the LAX method.
%
%  Discussion:
%
%    The Lax method is stable for the advection problem, if the time step
%    satisifies the Courant-Friedrichs-Levy (CFL) condition:
%
%      dt <= dx / c
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_ADVECTION_LAX:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the constant-velocity advection equation in 1D,\n' );
  fprintf ( 1, '    du/dt = - c du/dx\n' );
  fprintf ( 1, '  over the interval:\n' );
  fprintf ( 1, '    0.0 <= x <= 1.0\n' );
  fprintf ( 1, '  with periodic boundary conditions, and\n' );
  fprintf ( 1, '  with a given initial condition\n' );
  fprintf ( 1, '    u(0,x) = (10x-4)^2 (6-10x)^2 for 0.4 <= x <= 0.6\n' );
  fprintf ( 1, '           = 0 elsewhere.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We modify the FTCS approach using the Lax method:\n' );
  fprintf ( 1, '    du/dt = (u(t+dt,x)-0.5*u(t,x-dx)-0.5*u(t,x+dx))/dt\n' );
  fprintf ( 1, '    du/dx = (u(t,x+dx)-u(t,x-dx))/2/dx\n' );

  nx = 101;
  dx = 1.0 / ( nx - 1 );
  x = linspace ( 0.0, 1.0, nx );
  nt = 1000;
  dt = 1.0 / nt;
  c = 1.0;

  u = zeros ( 1, nx );
  i = find ( 0.4 <= x & x <= 0.6 );
  u(i) = ( 10.0 * x(i) - 4.0 ).^2 .* ( 6.0 - 10.0 * x(i) ).^2;

  iplot = 1;
  uplot(iplot,:) = u(:)';
  tplot(iplot) = 0.0;
  plotstep = ceil ( nt / 50 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes NX = %d\n', nx );
  fprintf ( 1, '  Number of time steps NT = %d\n', nt );
  fprintf ( 1, '  Constant velocity C = %g\n', c );
  fprintf ( 1, '  CFL condition: dt (%g) <= dx / c (%g)\n', dt, dx / c );

  im1 = [ nx, 1:nx-2, nx-1 ];
  i   = [ 1,  2:nx-1, nx   ];
  ip1 = [ 2,  3:nx,   1    ];

  for j = 1 : nt

    unew(i) = 0.5 * ( u(im1) + u(ip1) ) - c * dt / dx / 2.0 * ( u(ip1) - u(im1) );
    u(i) = unew(i);

    if ( rem ( j, plotstep ) < 1 )
      iplot = iplot + 1;
      uplot(iplot,:) = u(:)';
      tplot(iplot) = j * dt;
    end

  end
%
%  Plot.
%
  mesh ( x, tplot, uplot );
  xlabel ( '<--X-->' );
  ylabel ( '<--T-->');
  title ( 'U(X,T)');

  filename = 'fd1d_advection_lax.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot as "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_ADVECTION_LAX\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [m,n,newE] = grValidation(E);
% The validation of array E - auxiliary function for GrTheory Toolbox.
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

if ~isnumeric(E),
  error('The array E must be numeric!')
end
if ~isreal(E),
  error('The array E must be real!')
end
se=size(E); % size of array E
if length(se)~=2,
  error('The array E must be 2D!')
end
if (se(2)<2),
  error('The array E must have 2 or 3 columns!')
end
if ~all(all(E(:,1:2)>0))
  error('1st and 2nd columns of the array E must be positive!')
end
if ~all(all((E(:,1:2)==round(E(:,1:2))))),
  error('1st and 2nd columns of the array E must be integer!')
end
m=se(1);
if se(2)<3, % not set the weight
  E(:,3)=1; % all weights =1
end
newE=E(:,1:3);
n=max(max(newE(:,1:2))); % number of vertexes
return

author: HuangCongQing


text:
function vtuk_puvw_write ( output_unit, title, node_num, element_num, ...
  element_order, xyz, element_node, p, uvw )

%*****************************************************************************80
%
%% VTU_PUVW_WRITE writes pressure and velocity data to a VTU file.
%
%  Discussion:
%
%    The data is assumed to have been computed by a finite element program
%    for a 3D geometry which has been meshed using tetrahedral elements
%    of 4 or 10 nodes.
%
%    The solution data includes the pressure and velocity vector at each node.
%
%    The VTU format used here is the modern XML-based format used by the
%    Visual Toolkit for unstructured grid data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OUTPUT_UNIT, the output unit.
%
%    Input, string TITLE, a title for the data.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, real XYZ(3,NODE_NUM), the node coordinates.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the
%    nodes that make up each element.  The node indices are zero based.
%
%    Input, real P(1,NODE_NUM), the pressure at each node.
%
%    Input, real UVW(3,NODE_NUM), the velocity at each node.
%
  if ( element_order == 10 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VTU_PUVW_WRITE - Note:\n' );
    fprintf ( 1, '  As a temporary measure, we are handling quadratic tets\n' );
    fprintf ( 1, '  as though they were linear.  That is, the element information\n' );
    fprintf ( 1, '  will only use the vertices, not the midside nodes.\n' );
    element_order = 4;
  end

  fprintf ( output_unit, '<VTKFile type="unstructuredGrid" version="0.1" byte_order="BigEndian">\n' );
  fprintf ( output_unit, '  <UnstructuredGrid>\n' );
  fprintf ( output_unit, '    <Piece NumberOfPoints="%d" NumberOfCells="%d">\n', ...
    node_num, element_num );
  fprintf ( output_unit, '      <PointData Scalars="scalars">\n' );
  fprintf ( output_unit, '        <DataArray type="Float64" Name="pressure" Format="ascii">\n' );
  for node = 1 : node_num
    fprintf ( output_unit, '%f\n', p(node ) );
  end
  fprintf ( output_unit, '        </DataArray>\n' );
  fprintf ( output_unit, '        <DataArray type="Float64" NumberOfComponents="3" Name="velocity" Format="ascii">\n' );
  for node = 1 : node_num
    fprintf ( output_unit, '%f  %f  %f\n', uvw(1:3,node) );
  end
  fprintf ( output_unit, '        </DataArray>\n' );
  fprintf ( output_unit, '      </PointData>\n' );
  fprintf ( output_unit, '      <Points>\n' );
  fprintf ( output_unit, '        <DataArray type="Float64" NumberOfComponents="3" Format="ascii">\n' );
  for node = 1 : node_num
    fprintf ( output_unit, '%f  %f  %f\n', xyz(1:3,node) );
  end
  fprintf ( output_unit, '        </DataArray>\n' );
  fprintf ( output_unit, '      </Points>\n' );
  fprintf ( output_unit, '      <Cells>\n' );

  fprintf ( output_unit, '        <DataArray type="Int32" Name="connectivity" Format="ascii">\n' );
  for element = 1 : element_num
    for order = 1 : element_order
      fprintf ( output_unit, '  %d', element_node(order,element) );
    end
    fprintf ( output_unit, '\n' );
  end
  fprintf ( output_unit, '        </DataArray>\n' );

  fprintf ( output_unit, '        <DataArray type="Int32" Name="offsets" Format="ascii">\n' );
  offset = 0;
  for element = 1 : element_num
    offset = offset + element_order;
    fprintf ( output_unit, '%d\n', offset );
  end
  fprintf ( output_unit, '        </DataArray>\n' );
  fprintf ( output_unit, '        <DataArray type="Int32" Name="types" Format="ascii">\n' );
  if ( element_order == 4 )
    for element = 1 : element_num
      fprintf ( output_unit, '10\n' );
    end
  elseif ( element_order == 10 )
    for element = 1 : element_num
      fprintf ( output_unit, '24\n' );
    end
  end
  fprintf ( output_unit, '        </DataArray>\n' );

  fprintf ( output_unit, '      </Cells>\n' );
  fprintf ( output_unit, '    </Piece>\n' );
  fprintf ( output_unit, '  </UnstructuredGrid>\n' );
  fprintf ( output_unit, '</VTKFILE>\n' );

  return
end


text:
%  computing the Householder transformation  w = H*v  with
%            H = I - 2*u*u'
%
%  syntax  >> w = HouseholderTransform(u,v)

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = liop(varargin)
% VL_LIOP Local Intensity Order Pattern descriptor
%   D = VL_LIOP(I) computes the LIOP descriptor of an image I, as
%   described by [1]. I is a gray-scale square image with odd side
%   length of class SINGLE. D is a column vector containing the LIOP
%   descriptor of I. Note that LIOP is also integrated in the VL_COVDET()
%   function for feature extraction.
%
%   VL_LIOP() accepts the following options:
%
%   NumNeighbours:: 4
%     Set the number of neighbours sampled to construct the order
%     pattern of each image pixel.
%
%   Radius:: 5
%     Set the radius of the circular neighbourhood used to sampled
%     the local order pattern of each pixel.
%
%   NumSpatialBins:: 6
%     Set the number of spatial pooling regions. The LIOP descriptor
%     has dimension factorial(NumNeighbours) * NumSpatialBins.
%
%   IntensityThreshold:: -0.02
%     Set the intensity threshold used to weight order patterns as they
%     are pooled into a histogram. A negative value is interpreted
%     as a fraction of the difference between the maximum and minimum
%     intensity in each local patch.
%
%   Verbose::
%     If specified, be verbose
%
%   REFERENCES::
%   [1] Z. Wang, B. Fan, F. Wu. Local Intensity Order Pattern for feature
%   description. In ICCV, 2011
%
%   See: LIOP, VL_COVDET(),
%   VL_HELP().
[varargout{1:nargout}] = vl_liop(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function r = isfinite(p)
%ISFINITE     Logical result:  polynomial p contains only finite components, y/n
%
%   r = isfinite(p)
%
%Result 1 iff all coefficients of p are finite
%

% written  08/07/02     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
%

  r = all(isfinite(p.c));

repo: OpenMAS


text:
clear;
t0=0;
tf=10;
b0=[0.2 0.2];
[t,b]=ode45('dfun3',[t0, tf],b0);
plot(b(:,1),b(:,2));
hold on
b0=[0.1 0.1];
[t,b]=ode45('dfun3',[t0, tf],b0);
plot(b(:,1),b(:,2));
hold on
b0=[0.2 0.15];
[t,b]=ode45('dfun3',[t0, tf],b0);
plot(b(:,1),b(:,2));
hold on
xlabel('x1');
ylabel('x2');

author: sfvsfv

repo: Mathematical-modeling