text:
function CreateEulerPanel(handles)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

    %Euler Panel
    handles.eulerPanel = uipanel('Parent',handles.toolPanel,'Units','normalize','FontSize',10,'Title','Registration Parameters',...
                                    'Tag','eulerPanel','Clipping','off',...
                                    'Position',[0 0 1 1],...
                                    'Visible','on');
    dx = 0; dy = 80;
    h380 = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[19.5308641975309+dx 214.470588235294+dy 115.086419753087 19],...
    'String','Rotation Scale:',...
    'Style','text',...
    'Tag','text97');

    handles.Euler_RS = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[143.395061728395+dx 210+dy 77.0493827160494 24.5882352941177],...
    'String','1.0',...
    'Style','edit',...
    'Tag','Euler_RS');

    h382 = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[19.5308641975309+dx 179.823529411765+dy 115.086419753087 19],...
    'String','Translation Scale:',...
    'Style','text',...
    'Tag','text98');

    handles.Euler_TS = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[143.395061728395+dx 175.352941176471+dy 77.0493827160494 23.4705882352941],...
    'String','0.001',...
    'Style','edit',...
    'Tag','Euler_TS');

    h384 = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[10.7530864197531+dx 146.294117647059+dy 123.864197530864 20.1176470588235],...
    'String','Initial Angle(Degree):',...
    'Style','text',...
    'Tag','text99');

    handles.Euler_InitAngle = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[143.395061728395+dx 142.941176470588+dy 77.0493827160494 23.4705882352941],...
    'String','0',...
    'Style','edit',...
    'Tag','Euler_InitAngle');

    handles.Euler_minstep = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[143.395061728395+dx 99.3529411764706+dy 76.0740740740741 26.8235294117647],...
    'String','0.01',...
    'Style','edit',...
    'Tag','Euler_minstep');

    handles.Euler_maxstep = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[143.395061728395+dx 62.4705882352942+dy 76.0740740740741 25.7058823529412],...
    'String','1.0',...
    'Style','edit',...
    'Tag','Euler_maxstep');

    h388 = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[19.5308641975309+dx 106.058823529412+dy 115.086419753087 19],...
    'String','Minimum Step:',...
    'Style','text',...
    'Tag','text100');

    h389 = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[19.5308641975309+dx 69.1764705882353+dy 115.086419753087 19],...
    'String','Maximum Step:',...
    'Style','text',...
    'Tag','text101');

    h390 = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[19.5308641975309+dx 31.1764705882353+dy 115.086419753087 19],...
    'String','Iteration Number :',...
    'Style','text',...
    'Tag','text102');

    handles.Euler_IterNum = uicontrol(...
    'Parent',handles.eulerPanel,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[143.395061728395+dx 25.5882352941177+dy 76.0740740740741 25.7058823529412],...
    'String','200',...
    'Style','edit',...
    'Tag','Euler_IterNum');

    guidata(handles.mainframe, handles);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ fea, out ] = ex_navierstokes3( varargin )
%EX_NAVIERSTOKES3 2D Example for incompressible stationary flow around a cylinder.
%
%   [ FEA, OUT ] = EX_NAVIERSTOKES3( VARARGIN ) Stationary flow around a cylinder. References:
%
%   [1] John V, Matthies G. Higher-order finite element discretizations in a
%       benchmark problem for incompressible flows. International Journal for
%       Numerical Methods in Fluids 2001.
%
%   [2] Nabh G. On higher order methods for the stationary incompressible
%       Navier-Stokes equations. PhD Thesis, Universitaet Heidelberg, 1998.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       rho         scalar {1}             Density
%       miu         scalar {0.001}         Molecular/dynamic viscosity
%       umax        scalar {0.3}           Maximum magnitude of inlet velocity
%       igrid       scalar {2}             Grid type: >0 regular (igrid refinements)
%                                                     <0 unstruc. grid (with hmax=|igrid|)
%       sf_u        string {sflag1}        Shape function for velocity
%       sf_p        string {sflag1}        Shape function for pressure
%       iphys       scalar 0/{1}           Use physics mode to define problem (=1)
%       solver      string {}              Solver selection default, openfoam, su2, fenics
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct
%
%   See also EX_NAVIERSTOKES3B

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
            'rho',      1;
            'miu',      0.001;
            'umax',     0.3;
            'igrid',    2;
            'sf_u',     'sflag1';
            'sf_p',     'sflag1';
            'iphys',    1;
            'solver',   '';
            'iplot',    1;
            'tol',      [0.05 0.35 0.1];
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Model parameters.
rho       = opt.rho;     % Density.
miu       = opt.miu;     % Molecular/dynamic viscosity.
umax      = opt.umax;    % Maximum magnitude of inlet velocity.
umean     = 2/3*umax;    % Mean inlet velocity.
% Geometry and grid parameters.
h         = 0.41;        % Height of rectangular domain.
l         = 2.2;         % Length of rectangular domain.
xc        = 0.2;         % x-coordinate of cylinder center.
yc        = 0.2;         % y-coordinate of cylinder center.
diam      = 0.1;         % Diameter of cylinder.
% Discretization parameters.
sf_u      = opt.sf_u;    % FEM shape function type for velocity.
sf_p      = opt.sf_p;    % FEM shape function type for pressure.

% Grid generation.
fea.sdim = { 'x' 'y' };
if( opt.igrid>=1 )
  fea.grid = cylbenchgrid( opt.igrid );
  fea.grid.s(:) = 1;
else
  gobj1 = gobj_rectangle( 0, 2.2, 0, 0.41, 'R1' );
  gobj2 = gobj_circle( [0.2 0.2], 0.05, 'C1' );
  fea.geom.objects = { gobj1 gobj2 };
  fea = geom_apply_formula( fea, 'R1-C1' );
  fea.grid = gridgen( fea, 'hmax', abs(opt.igrid), 'fid', fid );
end
if( opt.iphys==1 && strcmp(opt.solver,'fenics') &&  size(fea.grid.c,1)==4 )
  warning( 'Converting quadrilateral mesh to triangular to support FEniCS.' )
  fea.grid = quad2tri( fea.grid );
end
n_bdr = max(fea.grid.b(3,:));   % Number of boundaries.

% Boundary conditions.
dtol      = sqrt(eps)*1e3;
i_inflow  = findbdr( fea, ['x<=',num2str(dtol)] );     % Inflow boundary number.
i_outflow = findbdr( fea, ['x>=',num2str(l-dtol)] );   % Outflow boundary number.
s_inflow  = ['4*',num2str(umax),'*(y*(',num2str(h),'-y))/',num2str(h),'^2'];   % Definition of inflow profile.
i_cyl     = findbdr( fea, ['sqrt((x-',num2str(xc),').^2+(y-',num2str(yc),').^2)<=(',num2str(diam/2+dtol),')'] );    % Cylinder boundary number.

% Problem definition.
if ( opt.iphys==1 )

  fea = addphys(fea,@navierstokes);     % Add Navier-Stokes equations physics mode.
  fea.phys.ns.eqn.coef{1,end} = { rho };
  fea.phys.ns.eqn.coef{2,end} = { miu };
  fea.phys.ns.sfun            = { sf_u sf_u sf_p };     % Set shape functions.

  fea.phys.ns.bdr.sel(i_inflow)  = 2;
  fea.phys.ns.bdr.sel(i_outflow) = 4;
  fea.phys.ns.bdr.coef{2,end}{1,i_inflow} = s_inflow;   % Set inflow profile.
  fea = parsephys(fea);                 % Check and parse physics modes.

else

  fea.dvar  = { 'u'  'v'  'p'  };       % Dependent variable name.
  fea.sfun  = { sf_u sf_u sf_p };       % Shape function.

  % Define equation system.
  cvelx = [num2str(rho),'*',fea.dvar{1}];   % Convection velocity in x-direction.
  cvely = [num2str(rho),'*',fea.dvar{2}];   % Convection velocity in y-direction.
  fea.eqn.a.form = { [2 3 2 3;2 3 1 1]       [2;3]                   [1;2];
                     [3;2]                   [2 3 2 3;2 3 1 1]       [1;3];
                     [2;1]                   [3;1]                   []   };
  fea.eqn.a.coef = { {2*miu miu cvelx cvely}  miu                    -1;
                     miu                     {miu 2*miu cvelx cvely} -1;
                     1                       1                      [] };
  fea.eqn.f.form = { 1 1 1 };
  fea.eqn.f.coef = { 0 0 0 };

  % Define boundary conditions.
  fea.bdr.d = cell(3,n_bdr);
  [fea.bdr.d{1:2,:}]         = deal( 0);

  fea.bdr.d{1,i_inflow}      = s_inflow;

  [fea.bdr.d{:,i_outflow  }] = deal([]);

  fea.bdr.n = cell(3,n_bdr);
end

% Parse and solve problem.
fea = parseprob(fea);             % Check and parse problem struct.
if( opt.iphys==1 && strcmp(opt.solver,'openfoam') )
  logfid = fid; if( ~got.fid ), fid = []; end
  fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid );
  fid = logfid;
elseif( opt.iphys==1 && strcmp(opt.solver,'su2') )
  logfid = fid; if( ~got.fid ), fid = []; end
  fea.sol.u = su2( fea, 'fid', fid, 'logfid', logfid );
  fid = logfid;
elseif( opt.iphys==1 && strcmp(opt.solver,'fenics') )
  fea = fenics( fea, 'fid', fid );
else
  jac.form  = {[1;1] [1;1] [];[1;1] [1;1] []; [] [] []};
  jac.coef  = {[num2str(rho),'*ux'] [num2str(rho),'*uy'] []; [num2str(rho),'*vx'] [num2str(rho),'*vy'] []; [] [] []};
  fea.sol.u = solvestat( fea, 'fid', fid, 'nsolve', 2, 'jac', jac, 'nlrlx', '(1+(it>2))/2' );   % Call to stationary solver.
end

% Postprocessing.
s_velm = 'sqrt(u^2+v^2)';
if ( opt.iplot>0 )
  figure
  subplot(2,1,1)
  postplot( fea, 'surfexpr', s_velm )
  title( 'Velocity field' )
  subplot(2,1,2)
  postplot( fea, 'surfexpr', 'p' )
  title( 'Pressure' )
end

% Calculate benchmark quantities (line integration method).
s_tfx = ['nx*p+',num2str(miu),'*(-2*nx*ux-ny*(uy+vx))'];
s_tfy = ['ny*p+',num2str(miu),'*(-nx*(vx+uy)-2*ny*vy)'];
s_cd  = ['2*(',s_tfx,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
s_cl  = ['2*(',s_tfy,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
i_cub = 10;
c_d1  = intbdr(s_cd,fea,i_cyl,i_cub);
c_l1  = intbdr(s_cl,fea,i_cyl,i_cub);
dp    = evalexpr('p',[0.15 0.25;0.2 0.2],fea);

% Calculate benchmark quantities (volume integration method).
bdrm   = fea.bdr.bdrm{1};
ind_b  = [];
ind_bm = [];
for ii=i_cyl
  ind_b  = [ind_b  find(fea.grid.b(3,:)==ii)];
  ind_bm = [ind_bm find(bdrm(3,:)==ii)];
end
ind_c    = fea.grid.b(1,ind_b);
ind_gdof = bdrm(4,ind_bm);

% Create field 'a' with values one on the cylinder and zero everywhere else.
fea.dvar = [ fea.dvar, {'a'}       ];
fea.sfun = [ fea.sfun, fea.sfun(1) ];
fea      = parseprob(fea);
n_dof    = max(fea.eqn.dofm{1}(:));
u_a      = zeros(n_dof,1);
u_a(ind_gdof) = 1;
fea.sol.u= [fea.sol.u;u_a];
fea.eqn  = struct;
fea.bdr  = struct;
fea      = parseprob(fea);

s_tfx    = ['ax*p+',num2str(miu),'*(-2*ax*ux-ay*(uy+vx))-(u*ux+v*uy)*a'];
s_tfy    = ['ay*p+',num2str(miu),'*(-ax*(vx+uy)-2*ay*vy)-(u*vx+v*vy)*a'];
s_cd     = ['2*(',s_tfx,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
s_cl     = ['2*(',s_tfy,')/(',num2str(rho),'*',num2str(umean),'^2*',num2str(diam),')'];
c_d2     = intsubd(s_cd,fea,[],[],3);
c_l2     = intsubd(s_cl,fea,[],[],3);

if( ~isempty(fid) )
  fprintf(fid,'\n\nBenchmark quantities:\n\n')

  fprintf(fid,'Drag coefficient,    cd = %6f (l), %6f (v) (Ref: 5.579535)\n',c_d1,c_d2)
  fprintf(fid,'Lift coefficient,    cl = %6f (l), %6f (v) (Ref: 0.010619)\n',c_l1,c_l2)
  fprintf(fid,'Pressure,            dp = %6f (Ref: 0.117520)\n',dp(1)-dp(2))
end

% Error checking.
out.cd   = [c_d1 c_d2];
out.cl   = [c_l1 c_l2];
out.dp   = dp(1)-dp(2);
out.err  = [abs(out.cd-5.579535)/5.579535;
            abs(out.cl-0.010619)/0.010619;
            abs(dp(1)-dp(2)-0.117520)/0.117520 0];
out.pass = (out.err(1,2)<opt.tol(1))&&(out.err(2,2)<opt.tol(2))&&(out.err(3,1)<opt.tol(3));
if ( nargout==0 )
  clear fea out
end

%------------------------------------------------------------------------------%
function [ grid ] = cylbenchgrid( nlev )
% CYLBENCHGRID Generate 2d quadrilateral grid for the DFG cylinder benchmark.

  ns = 8*2^(nlev-1);
  r  = [0.05 0.06 0.08 0.11 0.15];
  x  = [0.41 0.5 0.7 1 1.4 1.8 2.2];
  for ilev=2:nlev
    r = sort( [ r (r(1:end-1)+r(2:end))/2 ] );
    x = sort( [ x (x(1:end-1)+x(2:end))/2 ] );
  end

  grid1 = ringgrid( r, 4*ns, [], [], [0.2;0.2] );
  grid2 = holegrid( ns, 2^(nlev-1), [0 0.41;0 0.41], 0.15, [0.2;0.2] );
  grid2 = gridmerge( grid1, 5:8, grid2, 1:4 );
  grid3 = rectgrid( x, ns, [0.41 2.2;0 0.41] );
  grid  = gridmerge( grid3, 4, grid2, 6 );

author: precise-simulation

repo: featool-multiphysics


text:
function [ n_data, nu, x, fx ] = bessel_in_values ( n_data )

%*****************************************************************************80
%
%% BESSEL_IN_VALUES returns some values of the In Bessel function.
%
%  Discussion:
%
%    The modified Bessel functions In(Z) and Kn(Z) are solutions of
%    the differential equation
%
%      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
%
%    In Mathematica, the function can be evaluated by:
%
%      BesselI[n,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer NU, the order of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 28;

  fx_vec = [ ...
     0.5016687513894678E-02, ...
     0.1357476697670383E+00, ...
     0.6889484476987382E+00, ...
     0.1276466147819164E+01, ...
     0.2245212440929951E+01, ...
     0.1750561496662424E+02, ...
     0.2281518967726004E+04, ...
     0.3931278522104076E+08, ...
     0.2216842492433190E-01, ...
     0.2127399592398527E+00, ...
     0.1033115016915114E+02, ...
     0.1758380716610853E+04, ...
     0.2677764138883941E+21, ...
     0.2714631559569719E-03, ...
     0.9825679323131702E-02, ...
     0.2157974547322546E+01, ...
     0.7771882864032600E+03, ...
     0.2278548307911282E+21, ...
     0.2752948039836874E-09, ...
     0.3016963879350684E-06, ...
     0.4580044419176051E-02, ...
     0.2189170616372337E+02, ...
     0.1071597159477637E+21, ...
     0.3966835985819020E-24, ...
     0.4310560576109548E-18, ...
     0.5024239357971806E-10, ...
     0.1250799735644948E-03, ...
     0.5442008402752998E+19 ];

  nu_vec = [ ...
     2,  2,  2,  2, ...
     2,  2,  2,  2, ...
     3,  3,  3,  3, ...
     3,  5,  5,  5, ...
     5,  5, 10, 10, ...
    10, 10, 10, 20, ...
    20, 20, 20, 20 ];

  x_vec = [ ...
      0.2E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      2.5E+00, ...
      3.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     20.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    nu = 0;
    x = 0.0;
    fx = 0.0;
  else
    nu = nu_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function kern = ggwhiteKernExpandParam(kern, params)

% GGWHITEKERNEXPANDPARAM Create kernel structure from GG white kernel's parameters.
% FORMAT
% DESC returns a gaussian gaussian white
%	kernel structure filled with the parameters in the given vector.
%	This is used as a helper function to enable parameters to be
%	optimised in, for example, the NETLAB optimisation functions.
% RETURN kern : kernel structure with the given parameters in the relevant
%	   locations.
% ARG kern : the kernel structure in which the parameters are to be
%	   placed.
% ARG param : vector of parameters which are to be placed in the kernel
%	   structure.
%
% SEEALSO : ggwhiteKernParamINit, ggwhiteKernExtractParam, kernExpandParam
%
% COPYRIGHT : Mauricio A. Alvarez and Neil D. Lawrence, 2008
%
% MODIFICATIONS : Mauricio A Alvarez, 2009

% KERN

kern.precisionG = params(1:end-2)';
kern.sigma2Noise = params(end-1);
kern.variance = params(end);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [output1,output2] = function_name(input1,input2,input3)
%FUNCTION_NAME - One line description of what the function or script performs (H1 line)
%Optional file header info (to give more details about the function than in the H1 line)
%Optional file header info (to give more details about the function than in the H1 line)
%Optional file header info (to give more details about the function than in the H1 line)
%
% Syntax:  [output1,output2] = function_name(input1,input2,input3)
%
% Inputs:
%    input1 - Description
%    input2 - Description
%    input3 - Description
%
% Outputs:
%    output1 - Description
%    output2 - Description
%
% Example:
%    Line 1 of example
%    Line 2 of example
%    Line 3 of example
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: OTHER_FUNCTION_NAME1,  OTHER_FUNCTION_NAME2

% Author: Denis Gilbert, Ph.D., physical oceanography
% Maurice Lamontagne Institute, Dept. of Fisheries and Oceans Canada
% email address: gilbertd@dfo-mpo.gc.ca
% Website: http://www.qc.dfo-mpo.gc.ca/iml/
% December 1999; Last revision: 12-May-2004

%------------- BEGIN CODE --------------

Enter your executable matlab commands here

%------------- END OF CODE --------------

author: Sable

repo: mcbench-benchmarks


text:
% Test whether stable conditional Gaussian inference works
% Make a linear dynamical system
%   X1 -> X2
%   |     |
%   v     v
%   Y1    Y2

intra = zeros(2);
intra(1,2) = 1;
inter = zeros(2);
inter(1,1) = 1;
n = 2;

X = 2; % size of hidden state
Y = 2; % size of observable state

ns = [X Y];
bnet = mk_dbn(intra, inter, ns, 'discrete', [], 'observed', 2);

x0 = rand(X,1);
V0 = eye(X);
C0 = rand(Y,X);
R0 = eye(Y);
A0 = rand(X,X);
Q0 = eye(X);

bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', x0, 'cov', V0);
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(Y,1), 'cov', R0, 'weights', C0);
bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', zeros(X,1), 'cov', Q0, 'weights', A0);

T = 5; % fixed length sequences

engine = {};
engine{end+1} = kalman_inf_engine(bnet);
engine{end+1} = scg_unrolled_dbn_inf_engine(bnet, T);
engine{end+1} = jtree_unrolled_dbn_inf_engine(bnet, T);

inf_time = cmp_inference_dbn(bnet, engine, T, 'check_ll', 0);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function tests = test_ft_preproc_detrend

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preproc_detrend

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testOptions(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nchan   = 8;
nsample = 1000;
dat     = randn(nchan, nsample) + 1;

result = [];
result{end+1} = ft_preproc_detrend(dat, 1, 1000);
result{end+1} = ft_preproc_detrend(dat, 1, 900);
result{end+1} = ft_preproc_detrend(dat, 101, 1000);
result{end+1} = ft_preproc_detrend(dat, 101, 900);

% all iterations were done with (slightly) different options, hence the results should not be equal
for i=1:numel(result)
  for j=(i+1):numel(result)
    assert(~isequal(result{i}, result{j}), 'the results %d and %d should not be equal', i, j);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function sumsq = sumsqcomp(dpixc, res, peval)
ROIx=peval.ROIx(1):peval.ROIx(2);
ROIy=peval.ROIy(1):peval.ROIy(2);
ROIz=peval.ROIz(1):peval.ROIz(2);
dpixcs = dpixc(ROIx,ROIy,ROIz);
dveccs = reshape(dpixcs, peval.numpix, peval.nt);
htr =res.htrace;
for ii=1:size(htr,1)
    htr(ii,3,:)=res.h(3,:);
    itmp=res.w*squeeze(htr(ii,:,:));
    isq=(dveccs-itmp).^2;
    sumsq(ii)=sum(isq(:));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function G = gaborPyramid(s, theta, sz);
% Create a gabor wavelet pyramid.
%
% G = gaborPyramid(<s=0:2>, <theta=22.5:22.5:180>, <sz=[256 256]>);
%
% BACKGROUND: this function creates a wavelet pyramid analogous to that
% used by Kay et al (Nature 2000) to represent images for application to
% fMRI activity patterns. This pyramid contains a series of Gabor wavelets,
% each similar to the one described by Tai Sing Lee, "Image Representation
% Using 2D Gabor Wavelets".
% (IEEE Transactions 1996: online at www.cnbc.cmu.edu/~tai/papers/pami.pdf)
%
% "Gabor" wavelets (or "Weyl-Heisenberg groups" in physics) are a combined
% harmonic and Gaussian function which the physicist Dennis Gabor showed
% convey an optimal combination of spatial (or temporal) position localization
% and frequency-domain localization. These functions serve as models of
% receptive fields for early visual neurons, esp. simple cells. See the Lee
% paper for more background.
%
% "Pyramid" refers to the fact that, for different spatial frequencies,
% there will be a different number of channels, spanning the different
% positions that can be localized. This code follows the Kay et al
% convention of having spatial frequencies be powers of 2 in cycles per
% image: for 2 cycles per image, there will be 2^2=4 possible positions;
% for 4 cycles per image, there will be 4^2 possible positions.
%
%
% INPUTS:
% s: vector of spatial frequencies to use. s should be an integer, and specify
% the log base 2 of the spatial frequency, in units of cycles / image. That is,
% 2^s is the number of cycles in the image for a given channel.
% E.g.: s=0 indicates a single cycle per image; s=1 indicates two cycles
% per image. A DC offset is automatically added to the list of channels as well,
% which has undefined spatial frequency (2^s = 0).
%
% theta: vector of orientations to use, in degrees counterclockwise from
% the x0 axis. That is, 0=horizontal, 45 = increasing in x0/y0 direction, 90 =
% vertical, 135 = increasing in x0 while decreasing in y0 direction. Note
% that theta is cyclic with modulo 180 (theta=200 is the same as theta=20).
%
% sz: size [x y] of the resulting wavelets.
%
% OUTPUTS:
% Return results in a structure G, with the following fields:
%	odd: 3-D matrix of size [imgSize(1) imgSize(2) nChannels].
%		odd-symmetric (sine) component channels.
%		The # of channels depends on omega and theta (see above), and
%		represents the even channels only in the complex function.
%	even: same size as odd, containing even-symmetric (cosine) component
%	channels.
%	cyclesPerImage: vector of wavelengths corresponding to each channel.
%	orientation: vector of orientations corresponding to each channel.
%	x0: vector of x0 position centers for each channel
%	y0: vector of y0 position centers for each channel.
%
% REFERENCES: We implement a version of the wavelet decomposition used in
% Kay et al., Nature 2008. See also www.cnbc.cmu.edu/~tai/papers/pami.pdf.
%
% ras, 03/2008.
if notDefined('s'),			s = 0:2;					end
if notDefined('theta'),		theta = 22.5:22.5:180;		end
if notDefined('sz'),		sz = [128 128];				end

%%%%% params / setup
%% constants
kappa = 2.5; % pi;  % constrains the spat. freq. bandwith to be 1 octave

%% parse inputs
if length(sz)==1  % square image, only 1 size specified
	sz = [sz sz];
end

%% init empty struct
G.odd = [];
G.even = [];
G.layer = [];
G.cyclesPerImage = [];
G.orientation = [];
G.x0 = [];
G.y0 = [];

%%%%% (1) determine the number of channels
% first, convert the spatial frequency param s into units of cycles per
% image (cpi):
cpi = 2 .^ s;

% we add one channel for the DC component:
cpi = [0 cpi];

% for each # of cycles/image, we have a square grid of positions for that
% frequency (see Kay et al, Sup Fig 2).
tmp = sum( cpi .^ 2 ); % # of joint spatial-frequency/position combinations

% for each layer except the first one, we have one channel for each
% orientation. So, figure out the total # of channels:
nChannels = tmp * length(theta) + 1;
G.nChannels = nChannels;

%%%%% (2) determine parameters for each channel
% now, assign to each channel a unique combination of spatial frequency
% (cyclesPerImage), orientation, and position (x0 and y0). We'll call the
% set of channels corresponding to this spat. freq a "layer" in the pyramid:
for layer = 1:length(cpi)
	% the 0 cycles/image layer is a special case: only 1 DC channel
	if cpi(layer)==0
		G.layer = 1;
		G.cyclesPerImage = 0;
		G.orientation = 0; % could be NaN; this will help w/ plotting
		G.x0 = 0;
		G.y0 = 0;

	else
		% # positions along x/y axis for this layer
		nGrid = cpi(layer);

		% # channels in this layer
		n = length(theta) * nGrid ^ 2;

		% append cpi, orientation values
		G.layer = [G.layer repmat(layer, [1 n])];
		G.cyclesPerImage = [ G.cyclesPerImage, repmat(cpi(layer), [1 n]) ];
		G.orientation = [ G.orientation, repmat(theta, [1 nGrid^2]) ];

		% get a grid of (x0, y0) centers for each position for this layer
		% these positions are in pixels from the uppper left-hand corner of
		% the image:
		yRange = linspace(0, sz(1), nGrid+2);  yRange = yRange(2:end-1);
		xRange = linspace(0, sz(2), nGrid+2);  xRange = xRange(2:end-1);
		[x0 y0] = meshgrid( xRange, yRange );

		% we tile the position space separately for each orientation
		x0 = repmat(x0(:)', [length(theta) 1]);
		y0 = repmat(y0(:)', [length(theta) 1]);
		G.x0 = [G.x0, x0(:)'];
		G.y0 = [G.y0, y0(:)'];
	end
end

%%%%% (3) create the channels
% initialize the empty channels
G.odd	= zeros(sz(1), sz(2), nChannels);
G.even	= zeros(sz(1), sz(2), nChannels);

% main loop
for n = 1:nChannels
	% get a map of (x,y) pixel positions relative to the center of the
	% Gabor function for this channel (x0, y0) across the image
	[x y] = meshgrid( [1:sz(2)] - G.x0(n), [1:sz(1)] - G.y0(n) );

	% express the orientation in radians
	theta = G.orientation(n) * pi/180;

	% convert spatial frequency from units of cycles per image to units of
	% radians per unit distance: we take as the image size the larger of
	% the two edge lengths for the image:
	omega = 2*pi * G.cyclesPerImage(n)*2 / max(sz);

	% compute the wavelet
	phi = gaborWavelet(x, y, theta, omega, kappa);

	% grab the odd and even components
	G.even(:,:,n)	= real(phi);
	G.odd(:,:,n)	= imag(phi);
end

% to save memory, let's make everything single-precision
for f = fieldnames(G)'
	G.(f{1}) = single( G.(f{1}) );
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function extras = write_extras(fname,extras)
% Write extra bits of information
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

[pth,nam,ext] = fileparts(fname);
switch ext
case {'.hdr','.img','.nii'}
    mname = fullfile(pth,[nam '.mat']);
case {'.HDR','.IMG','.NII'}
    mname = fullfile(pth,[nam '.MAT']);
otherwise
    mname = fullfile(pth,[nam '.mat']);
end
if isstruct(extras) && ~isempty(fieldnames(extras)),
    savefields(mname,extras);
end;

function savefields(fnam,p)
if length(p)>1, error('Can''t save fields.'); end;
fn = fieldnames(p);
for i_=1:length(fn),
    eval([fn{i_} '= p.' fn{i_} ';']);
end;
if str2num(version('-release'))>=14,
    fn = {'-V6',fn{:}};
end;
if numel(fn)>0,
    save(fnam,fn{:});
end;
return;


text:
function fig = draw_bndry(fig, bndryy, bndryx, clr)
% fig = draw_bndry(fig, bndryy, bndryx, clr)
% draws specified points (boundaries) into picture matrix (figure)

    bndryy = round(bndryy);
    bndryx = round(bndryx);
    c = size(fig, 1);
    d = size(fig, 2);
    e = size(fig, 3);
    % draw only pixel that are in the image region
    draw = bndryy >= 1 & bndryy <= c & bndryx >= 1 & bndryx <= d;

    % Manipulate Figure
    ind = sub2ind([c,d], bndryy(draw), bndryx(draw));
    for j = 1:e
        fig(ind) = clr(j);
        ind = ind + c*d;
    end
end


text:
function [cum] = lgpdens_cum(bb,x1,x2)
% [CUM] = LGPDENS_CUM(BB)
%
% Description
%   Given Bayesian Bootstrap estimated density, integrates it from point x1
%   to point x2

% Copyright (c) 2012 Ernesto Ulloa
% Copyright (c) 2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.addRequired('bb',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('x1',@(x) ~isempty(x) && isreal(x))
ip.addRequired('x2', @(x) ~isempty(x) && isreal(x))
ip.parse(bb,x1,x2)

 [p,pq,xt]=lgpdens(bb);
 I1=min(find(xt>x1));
 I2=max(find(xt<x2));
 sd=xt(2)-xt(1);
 cum=sd*trapz(p(I1:I2));

end


text:
function s = char(ebsd)
% ebsd -> char

s = [];
for p = 1:numel(ebsd.phaseMap)

  ss = [];
  ss = [ss, 'phase ',num2str(ebsd.phaseMap(p)) ' ']; %#ok<AGROW>

  CS = ebsd.CSList{p};
  if ischar(CS)
    ss = [ss, '(not Indexed): ' CS ', ']; %#ok<AGROW>
  else

    if ~isempty(CS.mineral)
      ss = [ss, '(' CS.mineral '): ']; %#ok<AGROW>
    else
      ss = [ss, ': ']; %#ok<AGROW>
    end

    ss = [ss 'symmetry ' CS.pointGroup ', ']; %#ok<AGROW>
  end

  ss = [ss, num2str(nnz(ebsd.phaseId == p)),' orientations '];    %#ok<AGROW>

  s = strvcat(s,ss);
end

author: mtex-toolbox

repo: mtex

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


text:
function nlp = truck_opt_cost(nlp, bounds, varargin)
    % add custom cost function
    sys = nlp.Plant;
    x = sys.States.x;
    dx = sys.States.dx;

    % 1. target states: 200*(x_f-x_des)'*P1*(x_f-x_des)
    rd = bounds.rd;
    x_des_val = sys.x_des(rd);
    x_f   = sys.States.x;
    x_d = SymVariable('xd',[8,1]);

    Q = diag([6 1 10 1 10 1]);
    R = 100;
    P1 = sys.pcare(Q,R);
    tar_pos = 100*(x_f-x_d).'*P1*(x_f-x_d);

    tar_pos_fn = SymFunction('target_pos',tar_pos,{x_f},{x_d});
    nlp = addNodeCost(nlp, tar_pos_fn, 'x', 'last', x_des_val');

    % 2. square intgeral of ddy and y
    ya = sys.VirtualConstraints.y.ActualFuncs;
    ya2 = ya{1};
    Dya2 = ya{3};
    ddya2 = Dya2*dx;

    yd2_square = SymFunction('yd2_square',tovector(500*ya2.^2),{x});
    nlp = addRunningCost(nlp, yd2_square, {'x'});

    ddyd2_square = SymFunction('ddyd2_square',tovector(0.1*ddya2.^2),{x,dx});
    nlp = addRunningCost(nlp, ddyd2_square, {'x','dx'});

    alpha_last = SymFunction('alpha_last',tovector(100*ya2.^2),{x});
    nlp = addNodeCost(nlp, alpha_last, {'x'},'last');

    % 6. norm(y) -> y[1]^2 + y[2]^2 + ... + y[n]^2 (remove the sqrt operator)
    y_squre = SymFunction('y_square',tovector(0.1*x(3).^2),{x});
    nlp = addNodeCost(nlp, y_squre, 'x', 'all');

    % 7. var(u) @todo how to compute the mean value

    % 8. max(abs(y))
    y_max = sys.Params.ymax;
    ymax_cost_fn = SymFunction('ymax_cost',tovector(20*y_max),{y_max});
    nlp = addNodeCost(nlp, ymax_cost_fn, 'ymax', 'first');
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function EEG = highPassAndICA(EEG, varargin)
% Perform a high-pass filter and ICA on data that has been Prepped
try
    params = vargin2struct(varargin);
    if isfield(EEG.etc, 'noiseDetection') && ...
            isfield(EEG.etc.noiseDetection, 'detrend') && ...
            ~isfield(params, 'detrendChannels')
        params.detrendChannels = EEG.etc.noiseDetection.detrend.detrendChannels;
    end
    [EEG, detrend]  = removeTrend(EEG, params);
    EEG.etc.highPassAndICA.detrend = detrend;

    %% Perform ICA on the data
    if isfield(EEG.etc, 'noiseDetection') && ...
            isfield(EEG.etc.noiseDetection, 'reference')
        referenceChannels = ...
            EEG.etc.noiseDetection.reference.referenceChannels;
        interpolatedChannels = ...
            EEG.etc.noiseDetection.reference.interpolatedChannels.all;
        channelsLeft = setdiff(referenceChannels, interpolatedChannels);
        pcaDim = length(channelsLeft) - 1;
        fprintf('%s: pcaDim = %d\n', EEG.setname, pcaDim);
        EEG = pop_runica(EEG, 'icatype', 'cudaica', 'extended', 1, ...
            'chanind', referenceChannels, 'pca', pcaDim);
        EEG.etc.highPassAndICA.ICA = ...
            ['extended infomax with pca Dim ' num2str(pcaDim)];
        EEG.etc.noiseDetection.reference = ...
                   cleanupReference(EEG.etc.noiseDetection.reference);
%     else
%         EEG = pop_runica(EEG, 'icatype', 'runica', 'extended', 1);
%         EEG.etc.highPassAndICA.ICA = 'runica extended infomax';
    end

catch mex
    errorMessages.highPassAndICA = ['failed highPassAndICA: ' getReport(mex)];
    errorMessages.status = 'unprocessed';
    EEG.etc.highPassAndICA.errors = errorMessages;
    fprintf(2, '%s\n', errorMessages.highPassAndICA);
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [X,Y,Width,Height,Lbl] = read_lay(layoutname);

% READ_LAY reads an electrode or gradiometer layout file
% Layout files are used for topoplotting and multiplotting.
%
% Use as
%   [X, Y, Width, Height, Lbl] = read_lay(layoutname)

% Copyright (C) 2003, Ole Jensen
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_lay.m 945 2010-04-21 17:41:20Z roboos $

if ~exist(layoutname, 'file')
  error(sprintf('could not open layout file: %s', layoutname));
end

% discard the channel number
[chNum,X,Y,Width,Height,Lbl] = textread(layoutname,'%f %f %f %f %f %q');

author: PatternRecognition

repo: OpenBMI


text:
%% (Internal) Reshape a matrix into a row vector
% Return x as a row vector. This function is useful when a function returns a
% column vector or matrix and you want to immediately reshape it in a functional
% way. Suppose f(a,b,c) returns a column vector, matlab will not let you write
% f(a,b,c)(:)' - you would have to first store the result. With this function you
% can write rowvec(f(a,b,c)) and be assured that the result is a row vector.
%
%    x = colvec(x)
%
% Arguments:
%
%      + x: data
%
% Output:
%
%      + x: columnized data
%
% Example:
%
% See also colvec
%
% This file is from pmtk3.googlecode.com
% Copyright 2008-2015
%
function x = rowvec(x)
    x = x(:)';
end

author: marianux

repo: ecg-kit

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


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  qdd = call_direct_dynamics_fast(input)
%  Auxiliar function for the simulink model SIMULATE_ROBOT_AND_CONTROLLER.
%  Rearranges the inputs coming from the simulink model and calls the
%  function accel.
%  As a result the instantaneous acceleration at each joint is returned.
%
%  See also ACCEL.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function qdd = call_direct_dynamics_2dofplanar_s3b(input)

global robot

% torque = input(1:2);   % Input torque at each joint
% q   = input(3:4);	   % Joint positions
% qd  = input(5:6);	   % Joint speeds
% fe=[0 0 0 0 0 0]'; %external forces applied
% torque
%
% % Compute acceleration
% qdd = forwarddynamics_2dofplanar(robot, q, qd, torque, 9.81, fe)

% global robot

%set friction to zero
%robot.friction = 0;

torque = input(1:2);   % Input torque at each joint
q   = input(3:4);	   % Joint positions
qd  = input(5:6);	   % Joint speeds

g = [0 -9.81 0]';
% Compute acceleration
qdd = accel(robot, q, qd, torque, g);

repo: ARTE


text:
classdef TestVigdergauzMicroStructure < handle

    properties (Access = private)
        designVariable
        unfittedMesh
        volume
        numericalVolume
        computation
        settings
        testName;
    end

    methods (Access = public)

        function obj = TestVigdergauzMicroStructure(cParams)
            obj.init(cParams);
            obj.createUnfittedMesh();
            obj.computeFractionVolume();
        end

        function error = computeError(obj)
            v = obj.volume;
            nv = obj.numericalVolume;
            error = abs(v - nv)/v;
        end

    end

    methods (Access = private)

        function init(obj, cParams)
            obj.volume   = cParams.volume;
            obj.testName = cParams.testName;
            s.testName = obj.testName;
            computer = TopOptComputer(s);
            computer.compute();
            obj.computation = computer.computation;
            obj.settings = computer.settings;
            obj.designVariable = obj.computation.designVariable;
        end

        function createUnfittedMesh(obj)
            meshBackground = obj.computation.designVariable.mesh;
            s.backgroundMesh  = meshBackground.innerMeshOLD;
            s.boundaryMesh    = meshBackground.boxFaceMeshes;
            cParams = SettingsMeshUnfitted(s);
            obj.unfittedMesh = UnfittedMesh(cParams);
            obj.unfittedMesh.compute(obj.designVariable.value); % peta aqui
        end

        function computeFractionVolume(obj)
            v = obj.unfittedMesh.computeMass();
            obj.numericalVolume = v;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function dy_dx = parabolderiv( x, y, k, graphics_flag )

% function dy_dx = parabolderiv( x, y, k, graphics_flag )
%
% default: example from Ref. [1], p. 322
%
% or try: parabolderiv( ( -2 : 0.01 : 2 )', exp( -2 : 0.01 : 2 )', 2 )
%         parabolderiv( ( -2 : 0.01 : 2 )', ( exp( -2 : 0.01 : 2 ) + 0.01 * rand( 1, 401 ) )', 25 )
%         parabolderiv( ( ( 0 : 200 * pi ) / 100 )', ( sin( 0 : 0.01 : 2 * pi ) + 0.01 * rand( 1, 629 ) )', 40 )
%           for demonstration of noisy input and edge deviations
%
% graphics_flag = 1: show result in graph (default)
%                    (includes comparison to 2 other methods)
%                 0: show no graph
%
% This program differentiates an empirical function (an array of evenly
% spaced values), for instance periodically sampled data from an
% experiment.
%
% Method used: basically, a parabola is fit from k points to the left to k
% points to the right of the point where the derivative is required, this
% is then analytically differentiated. These calculations are not actually
% performed, but the method as given by Lanczos in Ref. [1] is used: only
% 1 fit parameter of the parabola is needed, which is calculated directly
% from the data; at the edges a parabola is fit through the first (last)
% 2k points, from which the derivative is calculated directly. Additional
% information (on accuracy, &c.) in Ref. [2]. Noise is handled well. An
% example from Ref. [1] is included.
%
% The datapoints need to be equidistant. If the graphics_flag is set, the
% result of this method is compared to applying the standard Matlab diff()
% function on the raw data as well as on adjacent averaged filtered data
% (zero phase delay). The example commands given above illustrate this.
%
% References:
%
% [1] Title                    : Applied analysis / by Cornelius Lanczos
%     Author                   : Cornelius Lanczos
%     Edition                  : 3rd print.
%     Publisher                : London : Pitman, 1964
%     Pages                    : 539 p.
%     Bibliographic annotation : 1st print: 1957
%     see pp. 321 - 324
%
% [2] Title                    : Digital filters / by Richard W. Hamming
%     Author                   : Richard W. Hamming
%     Edition                  : 3rd ed.
%     Publisher                : Englewoood Cliffs : Prentice-Hall, 1989
%     Pages                    : XIV, 284 p.
%     Bibliographic annotation : 1st print: 1977
%     ISBN                     : 0-13-212812-8
%     see p. 137
%
% Last update 24-03-2004 by Robert Klein-Douwel
%
% mail: robertkdkd@yahoo.co.uk, R.J.H.Klein-Douwel@tue.nl
% web:  http://www.sci.kun.nl/mlf/robertkd/

linewidth = 1; % can be changed for plotting purposes

% fill in some default values

if nargin < 4
    graphics_flag = 1;
end

if nargin < 3 % take example from Ref. [1], p. 322
    y = [ 0
          4
         25
         50
         67.4
        124.9
        172.0
        201.4
        288.1
        321.3
        387.1 ];

    x = ( 0 : 1 : length( y ) - 1 )';

    k = 2;
end

% initialise stuff

if length( x ) ~= length( y )
    error( 'x and y vectors have different lengths (RKD)' );
end

n = length( y ); %number of points

if 2 * k + 1 > n
    error( 'k too large or too few datapoints (RKD)' );
end

dy = zeros( n, 1 );

% check equidistancy of x
dx = diff( x );
d2x = diff( dx );

rel_error_d2x = max( abs( d2x ) ) / min( dx );
if rel_error_d2x > 401 * eps
    % vectors like ( -2 : 0.01 : 2 ) are not equidistant, but have a small variation
    % in dx and d2x; if this is not more than 401 * eps, then it is neglected
    disp( [ 'relative error in d2x = ' num2str( rel_error_d2x ) ] );
    disp( [ 'relative error / eps = ' num2str( rel_error_d2x / eps ) ] );
    disp( [ 'min( d2x ) = ' num2str( min( d2x ) ) ] );
    disp( [ 'max( d2x ) = ' num2str( max( d2x ) ) ] );
    error( 'non-equidistant data points (RKD)' );
end
h = dx( 1 ); % stepsize

% start calculations
dy_denominator = 2 * h * sum( ( 1 : k ).^2 );

a = - k : + k;
for i = k + 1 : n - k
    dy( i ) = sum( a .* y( i + a )' );
end

dy = dy / dy_denominator;

% first and last k points:
% fit parabola over first and last 2k points
kk = 2 * k;

% first kk points
p1 = polyfit( x( 1 : kk ), y( 1 : kk ), 2 );
q1 = polyder( p1 ); % take derivative
dy( 1 : k ) = polyval( q1, x( 1 : k ) );

% last kk points
p2 = polyfit( x( n - kk + 1 : n ), y( n - kk + 1 : n ), 2 );
q2 = polyder( p2 ); % take derivative
dy( n - k + 1 : n ) = polyval( q2, x( n - k + 1 : n ) );

if nargout ~= 0 % output
    dy_dx = dy;
end

% end of calculations (everything below is just for making comparisons and nice output)

if graphics_flag == 1 % make some comparisons
    xdiff_range = x( 1 : n - 1 ) + h / 2; % x coordinates for results of diff() function

    % compare with diff()
    dydx_matlab = diff( y ) / h;

    if nargin >=3 % not enough data points in example to apply this filter

        % another test: try first adjacent averaging and then diff()
        % use k points to left and to right
        kkk = 2 * k + 1;
        b = ones( 1, kkk ) / kkk; % k point averaging filter
        y_filt = filtfilt( b, 1, y ); % noncausal filtering, zero phase distortion

        % calculate dy_filt/dx
        dy_filtdx = diff( y_filt ) / h;
    end

    if nargin < 3 % show results of example
        disp( '      x         y         dy' );
        disp( [ x, y, dy ] );
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% show results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if graphics_flag == 1
    if nargout == 0
        close all;
    end

    dx_idx_range = k + 1 : n - k;

    % parabolic fits to first (last) kk points
    x1_fit = min( x( 1 : kk ) ) : h / 100 : max( x( 1 : kk ) );
    y1_fit = polyval( p1, x1_fit );

    x2_fit = min( x( n - kk + 1 : n ) ) : h / 100 : max( x( n - kk + 1 : n ) );
    y2_fit = polyval( p2, x2_fit );

    figure;
    hold on;
    plot( x, y, '-ob', 'LineWidth', linewidth );
    plot( x( dx_idx_range ), dy( dx_idx_range ), '-or', 'LineWidth', linewidth );
    plot( x( 1 : k ), dy( 1 : k ), '-sr', 'LineWidth', linewidth );
    plot( x1_fit, y1_fit, '--r', 'LineWidth', linewidth );
    plot( xdiff_range, dydx_matlab, ':sg', 'LineWidth', linewidth );

    if nargin >= 3
        plot( xdiff_range, dy_filtdx, ':dm', 'LineWidth', linewidth );
    end

    legend( 'y', [ 'dy/dx_{Lanczos: {\pm}' num2str( k ) '}' ], ...
        [ 'first/last ' num2str( k ) ' pts' ], ...
        [ 'y_{fit, parabola, ' num2str( kk ) ' pts}' ], 'diff( y )/ dx', ...
        [ 'diff( y_{adj avg ({\pm}' num2str( k ) ')} )/ dx' ], 0 );

    plot( x( n - k + 1 : n ), dy( n - k + 1 : n ), '-sr', 'LineWidth', linewidth );
    plot( x2_fit, y2_fit, '--r', 'LineWidth', linewidth );

    plot( x, dy, '-r', 'LineWidth', linewidth );

    % replot to put these on top again
    plot( x( dx_idx_range ), dy( dx_idx_range ), '-or', 'LineWidth', linewidth );
    plot( x( 1 : k ), dy( 1 : k ), '-sr', 'LineWidth', linewidth );
    hold off;

    xlabel( 'x' );
    ylabel( 'y' );

    title( 'parabolic derivative (Lanczos)' );
end


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

function particle_system = move_particles(particle_system)
% MOVE_PARTICLES: Function moves the current particles according to an MCMC
% kernel with an adaptive logistic regression proposal.

% Extract function inputs
n_vars     = particle_system.n_vars;
model_orig = particle_system.models;
model_val  = particle_system.model_val;
rho        = particle_system.rho;
logistic_A = particle_system.logistic_A;

% Determine the number of current models, and dimensions
[n_models, ~] = size(model_orig);

% Compute the old particle diversity ratio
diversity_old = unique_particles(model_orig);

% Declare vector for accept_ratios
mh_accept_vect = [];

total_fn_count = particle_system.total_fn_count(end);

iter_temp = 0;
while(1)

    % Generate matrix to store new samples and data
    new_model = zeros(n_models, n_vars);
    new_mval  = zeros(n_models, 1);

    % Declare vector to store acceptance ratio
    accept_ratio = 0;

    for i=1:n_models

        % Propose new model
        [prop_model, ~] = logistic_sample(n_vars, logistic_A);

        % Evaluate proposal density under both models
        q_new = evaluate_logit(particle_system, prop_model);
        q_old = evaluate_logit(particle_system, model_orig(i,:));

        % Evaluate probability function for both models
        total_fn_count = total_fn_count + 1;
        model_obj_new = particle_system.objective(prop_model);
        log_post_new = rho*model_obj_new;
        log_post_old = rho*model_val(i);

        % Compute acceptance probability
        accept_prob = exp(log_post_new - log_post_old)*(q_old/q_new);

        % Accept proposed sample/data or assign old values
        if rand < min(1,accept_prob)
            new_model(i,:) = prop_model;
            new_mval(i)    = model_obj_new;
            accept_ratio   = accept_ratio + 1;
        else
            new_model(i,:) = model_orig(i,:);
            new_mval(i)    = model_val(i);
        end

    end

    %% Update Sample Information

    % Update variables
    model_orig     = new_model;
    model_val      = new_mval;
    mh_accept_vect = [mh_accept_vect, accept_ratio/n_models];

    % Compute the new particle diversity
    diversity_new = unique_particles(new_model);

    % Compare the diversity
    if (abs(diversity_new - diversity_old) < 0.02 || diversity_new > 0.95)
        break
    end

    % Update diversity_old
    diversity_old = diversity_new;

    iter_temp = iter_temp + 1;
    if iter_temp > 100
        break
    end

end

% Update particle_system
particle_system.models     = model_orig;
particle_system.model_val  = model_val;
particle_system.mh_accept  = min(mh_accept_vect);
particle_system.total_fn_count = [particle_system.total_fn_count, total_fn_count];

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
% This .m file, "view_mcgrid.m", plots ratechanges calculated with mcgrid.m
% needs re3, gx, gy
%
% define size of the plot etc.
%
if isempty(name) >  0
    name = '  '
end
think
report_this_filefun(mfilename('fullpath'));
%co = 'w';

% This is the info window text
%
ttlStr='The b and p -Value Map Window                 ';
hlpStr1zmap= ...
    ['                                                '
    ' This window displays b-values and p-values     '
    ' using a color code.                            '
    ' Some of the menu-bar options are               '
    ' described below:                               '
    '                                                '
    ' Threshold: You can set the maximum size that   '
    '   a volume is allowed to have in order to be   '
    '   displayed in the map. Therefore, areas with  '
    '   a low seismicity rate are not displayed.     '
    '   edit the size (in km) and click the mouse    '
    '   outside the edit window.                     '
    'FixAx: You can chose the minimum and maximum    '
    '        values of the color-legend used.        '
    'Polygon: You can select earthquakes in a        '
    ' polygon either by entering the coordinates or  '
    ' defining the corners with the mouse            '];
hlpStr2zmap= ...
    ['                                                '
    'Circle: Select earthquakes in a circular volume:'
    '      Ni, the number of selected earthquakes can'
    '      be edited in the upper right corner of the'
    '      window.                                   '
    ' Refresh Window: Redraws the figure, erases     '
    '       selected events.                         '

    ' zoom: Selecting Axis -> zoom on allows you to  '
    '       zoom into a region. Click and drag with  '
    '       the left mouse button. type <help zoom>  '
    '       for details.                             '
    ' Aspect: select one of the aspect ratio options '
    ' Text: You can select text items by clicking.The'
    '       selected text can be rotated, moved, you '
    '       can change the font size etc.            '
    '       Double click on text allows editing it.  '
    '                                                '
    '                                                '];

% Set up the Seismicity Map window Enviroment
%

% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('rc-value-map',1);
newrcmapWindowFlag=~existFlag;

if newrcmapWindowFlag
    oldfig_button = 0
end

if oldfig_button == 0
    rcmap = figure_w_normalized_uicontrolunits( ...
        'Name','rc-value-map',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'backingstore','on',...
        'Visible','off', ...
        'Position',[ (fipo(3:4) - [600 400]) ZmapGlobal.Data.map_len]);
    % make menu bar
    matdraw

    %lab1 = 'p-value:';

    add_symbol_menu('eq_plot');

    %    uicontrol('Units','normal',...
    %       'Position',[.0 .93 .08 .06],'String','Info ',...
    %        'Callback',' web([''file:'' hodi ''/zmapwww/chp11.htm#996756'']) ');

    options = uimenu('Label',' Analyze ');
    uimenu(options,'Label','Refresh ', 'Callback','view_tcgrid')
    %    uimenu(options,'Label','Select EQ in Circle',...
    %        'Callback','h1 = gca;met = ''ni''; ho=false;cirpva;watchoff(rcmap)')
    uimenu(options,'Label','Select EQ in Circle - Constant R',...
         'Callback','h1 = gca;met = ''ra''; ho=false;plot_circ_FMD2Periods;watchoff(rcmap)')
    uimenu(options,'Label','Select EQ with const. number',...
         'Callback','h1 = gca;ho2=true;ho=true;plot_constnrbootfit_a2;watchoff(rcmap)')

    op1 = uimenu('Label',' Maps ');

    %Menu for adjusting several parameters.
    adjmenu =  uimenu(op1,'Label','Adjust Map Display Parameters'),...
        uimenu(adjmenu,'Label','Adjust Mmin cut',...
         'Callback','asel = ''mag''; adju2; view_mcgrid ')
    uimenu(adjmenu,'Label','Adjust Rmax cut',...
         'Callback','asel = ''rmax''; adju2; view_mcgrid')
    uimenu(adjmenu,'Label','Adjust goodness of fit cut',...
         'Callback','asel = ''gofi''; adju2; view_mcgrid ')
    uimenu(adjmenu,'Label','Adjust p-value st. dev. cut',...
         'Callback','asel = ''pstdc''; adju2; view_mcgrid ')
    % Menu MAPS
    uimenu(op1,'Label','Magnitude shift',...
         'Callback',' lab1=''dM'';re3 = mMagShift; view_mcgrid')
    uimenu(op1,'Label','KS-Test result',...
         'Callback','lab1=''H'';re3 = mHkstest; view_mcgrid')
    uimenu(op1,'Label','Magnitude shift validated',...
         'Callback',' lab1=''dM'';re3 = mMagShift_valid; view_mcgrid')
    uimenu(op1,'Label','Normalized \Delta_{FMD}',...
         'Callback','lab1=''\Delta_{FMD}'';re3 = mChFMD; view_mcgrid')
    uimenu(op1,'Label','Resolution Map (Number of events period 1)',...
         'Callback','lab1=''Number of events'';re3 = mNumevents1; view_mcgrid')
    uimenu(op1,'Label','Resolution Map (Number of events period 2)',...
         'Callback','lab1=''Number of events'';re3 = mNumevents2; view_mcgrid')
    uimenu(op1,'Label','Mc period 1',...
         'Callback','lab1=''Mc'';re3 = mMc1; view_mcgrid')
    uimenu(op1,'Label','Mc period 2',...
         'Callback','lab1=''Mc'';re3 = mMc2; view_mcgrid')
    uimenu(op1,'Label','Mc change',...
         'Callback','lab1=''dMc'';re3 = mdMc; view_mcgrid')
    uimenu(op1,'Label','Mc change validated',...
         'Callback','lab1=''dMc'';re3 = mdMc_val; view_mcgrid')
    uimenu(op1,'Label','AIC difference Utsu',...
         'Callback','lab1=''dAIC'';re3 = mdAIC_Utsu; view_mcgrid')
    uimenu(op1,'Label','Probability of stationarity (favoring stationarity)',...
         'Callback','lab1=''Mc'';re3 = mStationary1_Utsu; view_mcgrid')
    uimenu(op1,'Label','Probability of stationarity (favoring non-stationarity)',...
         'Callback','lab1=''Prob'';re3 = mStationary2_Utsu; view_mcgrid')

    % Menu DISPLAY
    add_display_menu(1);

    %re3 = pvalg;
    tresh = nan; re4 = re3;
    oldfig_button = 1;

    colormap(jet)
    tresh = nan; minpe = nan; Mmin = nan; minsd = nan;

end   % This is the end of the figure setup.

% Now lets plot the color-map!
%
figure_w_normalized_uicontrolunits(rcmap)
delete(gca)
delete(gca)
delete(gca)
dele = 'delete(sizmap)';er = 'disp('' '')'; eval(dele,er);
reset(gca)
cla
hold off
watchon;
set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on','SortMethod','childorder')

rect = [0.18,  0.10, 0.7, 0.75];
rect1 = rect;

% find max and min of data for automatic scaling
maxc = max(max(re3));
maxc = fix(maxc)+1;
minc = min(min(re3));
minc = fix(minc)-1;

% set values greater tresh = nan
%
re4 = re3;

% plot image
%
orient landscape
%set(gcf,'PaperPosition', [0.5 1 9.0 4.0])

%Plots re4, which contains the filtered values.
axes('position',rect)
hold on
pco1 = pcolor(gx,gy,re4);

axis([ min(gx) max(gx) min(gy) max(gy)])
axis image
hold on
if sha == 'fl'
    shading flat
else
    shading interp
end
% make the scaling for the recurrence time map reasonable
if lab1(1) =='T'
    l = isnan(re3);
    re = re3;
    re(l) = [];
    caxis([min(re) 5*min(re)]);
end

%If the colorbar is freezed.
if fre == 1
    caxis([fix1 fix2])
end

title2([name ';  '   num2str(t0b) ' to ' num2str(teb) ],'FontSize',ZmapGlobal.Data.fontsz.s,...
    'Color','r','FontWeight','bold')

xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

% plot overlay
%
hold on
overlay_
ploeq = plot(a.Longitude,a.Latitude,'k.');
set(ploeq,'Tag','eq_plot','MarkerSize',ms6,'Marker',ty,'Color',co,'Visible',vi)

set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on','TickDir','out')
h1 = gca;
hzma = gca;

% Create a colorbar
%
h5 = colorbar('horiz');
set(h5,'Pos',[0.35 0.05 0.4 0.02],...
    'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s,'TickDir','out')

rect = [0.00,  0.0, 1 1];
axes('position',rect)
axis('off')
%  Text Object Creation
txt1 = text(...
    'Color',[ 0 0 0 ],...
    'EraseMode','normal',...
    'Units','normalized',...
    'Position',[ 0.33 0.06 0 ],...
    'HorizontalAlignment','right',...
    'Rotation',[ 0 ],...
    'FontSize',ZmapGlobal.Data.fontsz.s,....
    'FontWeight','bold',...
    'String',lab1);

% Make the figure visible
%
set(gca,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on','TickDir','out')
figure_w_normalized_uicontrolunits(rcmap);
%sizmap = signatur('ZMAP','',[0.01 0.04]);
%set(sizmap,'Color','k')
axes(h1)
watchoff(rcmap)
%whitebg(gcf,[ 0 0 0 ])
done

author: CelsoReyes


text:
function linplus_test572 ( )

%*****************************************************************************80
%
%% TEST572 tests R8SP_MXV, R8SP_VXM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 7;
  n = 5;
  nz_num = 10;
  col = [ 2, 5, 1, 5, 1, 2, 3, 4, 4, 1 ];
  row = [ 1, 1, 2, 2, 4, 4, 4, 5, 6, 7 ];
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST572\n' );
  fprintf ( 1, '  R8SP_MXV multiplies a R8SP matrix by a vector;\n' );
  fprintf ( 1, '  R8SP_VXM multiplies a vector by a R8SP matrix;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );
  fprintf ( 1, '  Matrix nonzeros =  %d\n', nz_num );
%
%  Set the matrix.
%
  [ a, seed ] = r8sp_random ( m, n, nz_num, row, col, seed );
%
%  Make a R8GE copy.
%
  c = r8sp_to_r8ge ( m, n, nz_num, row, col, a );
%
%  Print the R8GE copy.
%
  r8ge_print ( m, n, c, '  The R8SP matrix, in R8GE form:' );

  x(1) = 1.0E+00;
  x(2:n-1) = 0.0E+00;
  x(n) = -1.0E+00;

  r8vec_print ( n, x, '  The vector x:' );

  b = r8sp_mxv ( m, n, nz_num, row, col, a, x );

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

  x(1) = 1.0E+00;
  x(2:m-1) = 0.0E+00;
  x(m) = -1.0E+00;

  r8vec_print ( m, x, '  The vector x:' );

  b = r8sp_vxm ( m, n, nz_num, row, col, a, x );

  r8vec_print ( n, b, '  The product A'' * x:' );

  return
end


text:
function [els_sun,els_ode,els_spm] = mci_compare_forward (model)
% Compare integration methods
% FORMAT [els_sun,els_ode,els_spm] = mci_compare_forward (model)
%
% model     'phase', 'nmm-r2p2'
%
% Run integration 9 times - compare speed and accuracy
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_compare_forward.m 6548 2015-09-11 12:39:47Z will $

[P,M,U,Y] = mci_compare_setup (model);

M.reltol=1e-2;
M.abstol=1e-4;

M = spm_mci_minit (M);
if isstruct(M.pC)
    pC=full(diag(spm_vec(M.pC)));
else
    pC = M.pC;
end

UI.u=U';
UI.dt=M.T/M.N;

Nsamp=9;
for s=1:Nsamp,
    % New point
    P=spm_normrnd(M.vpE,pC,1);

    % Parameters in reduced space
    Pr = M.V'*(P-M.vpE);

    M.int='sundials';
    tic;
    y = spm_mci_fwd (P,M,U);
    ysun(:,s)=y(:,1);
    els_sun(s)=toc;

    M.int='ode15';
    tic;
    y = spm_mci_fwd (P,M,U);
    yode(:,s)=y(:,1);
    els_ode(s)=toc;

    % DCM for fMRI uses spm_int
    % DCM for ERP uses spm_int_L
    tic;
    y=spm_int_L(P,M,UI);
    yspm(:,s)=y(:,1);
    els_spm(s)=toc;

end

hs=figure;
set(hs,'Name','First Time Series');
lw=2;
rN=ceil(sqrt(Nsamp));
for s=1:Nsamp,
    subplot(rN,rN,s);
    plot(ysun(:,s),'k');
    hold on
    grid on
    plot(yode(:,s),'b');
    plot(yspm(:,s),'r');
end
legend('Sundials','ODE15','spm-int-L');

hs=figure;
set(hs,'Name','Integration Speed');
k=1; lw=2;
plot(els_ode,els_sun,'kx','MarkerSize',10);
hold on
mo=min(els_ode);
ma=max(els_ode);
plot([mo ma],[mo ma],'k','LineWidth',lw);
set(gca,'FontSize',18);
grid on
xlabel('ODE15');
ylabel('Sundials');

hs=figure;
set(hs,'Name','Integration Speed');
k=1; lw=2;
plot(els_spm,els_sun,'kx','MarkerSize',10);
hold on
mo=min(els_spm);
ma=max(els_spm);
plot([mo ma],[mo ma],'k','LineWidth',lw);
set(gca,'FontSize',18);
grid on
xlabel('spm-int-L');
ylabel('Sundials');

author: spm


text:
function q = subSet(q,ind)
% indexing of quaternions
%
% Syntax
%   subSet(q,ind) %
%

q.a = q.a(ind);
q.b = q.b(ind);
q.c = q.c(ind);
q.d = q.d(ind);


text:
function [h,h_normals] = visualize_model_and_normals(model,vector_scale,NS,point_size,color)

if(~exist('vector_scale','var'))
    vector_scale = 0;
end

if(~exist('NS','var'))
    NS = 1000;
end

if(~exist('point_size','var'))
    point_size = 20;
end

if(~exist('color','var'))
    color = 'black';
end

% Downsample the model to have a lighter version of it
N = length(model.x);

if(NS ~= 0)
    st = ceil(N/NS);
    model = downsample_model(model,st);
end

%--- display the model and the vectors
if(isfield(model,'r') && ~isempty(model.r))
    h = scatter3(model.x,model.y,model.z,point_size,[model.r ,model.g ,model.b ],'fill');

else
    h = scatter3(model.x,model.y,model.z,point_size,color,'fill');
end

if(vector_scale > 0)
    hold on;
    h_normals = quiver3(model.x,model.y,model.z,model.nx,model.ny,model.nz,vector_scale);
end

xlabel('x');
ylabel('y');
zlabel('z');

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function varargout = cosh(varargin)
%COSH   Hyperbolic cosine of a DISKFUN.
%   COSH(F) returns the hyperbolic cosine of F.
%
% See also DISKFUN/SINH, DISKFUN/COS.

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

[varargout{1:nargout}] = cosh@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function handle = skelVargplvmVisualise(yVal, chan, skel, padding)

% SKELVISUALISE For drawing a skel representation of 3-D data.
% FORMAT
% DESC draws a skeleton representation in a 3-D plot.
% ARG channels : the channels to update the skeleton with.
% ARG skel : the skeleton structure.
% RETURN handle : a vector of handles to the plotted structure.
%
% SEEALSO : skelVargplvmModify
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
% MODIFICATIONS: Andreas C. Damianou, 2014

% SHEFFIELDML

if nargin<4
  padding = 0;
end
if isempty(yVal)
    channels = chan(1,:);
else
    channels = demCmu35VargplvmLoadChannels(yVal,skel);
end
channels = [channels zeros(1, padding)];
vals = skel2xyz(skel, channels);
connect = skelConnectionMatrix(skel);

indices = find(connect);
[I, J] = ind2sub(size(connect), indices);
handle(1) = plot3(vals(:, 1), vals(:, 3), vals(:, 2), '.');
axis ij % make sure the left is on the left.
set(handle(1), 'markersize', 20);
%/~
%set(handle(1), 'visible', 'off')
%~/
hold on
grid on
for i = 1:length(indices)
  handle(i+1) = line([vals(I(i), 1) vals(J(i), 1)], ...
              [vals(I(i), 3) vals(J(i), 3)], ...
              [vals(I(i), 2) vals(J(i), 2)]);
  set(handle(i+1), 'linewidth', 2);
end
axis equal
xlabel('x')
ylabel('z')
zlabel('y')
axis on

author: SheffieldML

repo: vargplvm


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% 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.

function CstFR4lossy(mws)

%'@ define material: FR-4 (lossy)

material = invoke(mws,'material');
invoke(material,'Reset');
invoke(material,'Name','FR-4 (lossy)');
invoke(material,'FrqType','all');
invoke(material,'Type','Normal');
invoke(material,'SetMaterialUnit','GHz','mm');
invoke(material,'Epsilon','4.3');
invoke(material,'Mue','1.0');
invoke(material,'Kappa','0.0');
invoke(material,'TanD','0.025');
invoke(material,'TanDFreq','10.0');
invoke(material,'TanDGiven','True');
invoke(material,'TanDModel','ConstTanD');
invoke(material,'KappaM','0.0');
invoke(material,'TanDM','0.0');
invoke(material,'TanDMFreq','0.0');
invoke(material,'TanDMGiven','False');
invoke(material,'TanDMModel','ConstKappa');
invoke(material,'DispModelEps','None');
invoke(material,'DispModelMue','None');
invoke(material,'DispersiveFittingSchemeEps','General 1st');
invoke(material,'DispersiveFittingSchemeMue','General 1st');
invoke(material,'UseGeneralDispersionEps','False');
invoke(material,'UseGeneralDispersionMue','False');
invoke(material,'Rho','0.0');
invoke(material,'ThermalType','Normal');
invoke(material,'ThermalConductivity','0.3');
invoke(material,'SetActiveMaterial','all');
invoke(material,'Colour','0.94','0.82','0.76');
invoke(material,'Wireframe','False');
invoke(material,'Transparency','0');
invoke(material,'Create');
release(material);
end

repo: CST-MATLAB-API


text:
function [LIFT] = Lift2D()

% function [LIFT] = Lift2D()
% Purpose  : Compute surface to volume lift term for DG formulation

Globals2D;
Emat = zeros(Np, Nfaces*Nfp);

% face 1
faceR = r(Fmask(:,1));
V1D = Vandermonde1D(N, faceR);
massEdge1 = inv(V1D*V1D');
Emat(Fmask(:,1),1:Nfp) = massEdge1;

% face 2
faceR = r(Fmask(:,2));
V1D = Vandermonde1D(N, faceR);
massEdge2 = inv(V1D*V1D');
Emat(Fmask(:,2),Nfp+1:2*Nfp) = massEdge2;

% face 3
faceS = s(Fmask(:,3));
V1D = Vandermonde1D(N, faceS);
massEdge3 = inv(V1D*V1D');
Emat(Fmask(:,3),2*Nfp+1:3*Nfp) = massEdge3;

% inv(mass matrix)*\I_n (L_i,L_j)_{edge_n}
LIFT = V*(V'*Emat);
return

repo: Aero-matlab


text:
function [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,y,k,numCentralKnots,useSparse)
%%BSPLINEPOLYFITMULTIDIM Given a set of real vector (or scalar) points as
%         well as real or complex scalar function values, this function
%         returns the knots and coefficients for a set of b-splines that
%         can be used to interpolate to a desired polynomial degree within
%         the sampled region of the function. B-splines can have continuous
%         derivatives across the entire sampled region and can be easily
%         differentiated. This function determines the knots to use by
%         taking the tensor product of the approximate knots from
%         Equation 10 of Chapter 13 of [1]. The knots are not optimal. The
%         first and last knots in each dimensions are repeated. This
%         function can provide interpolation coefficients for multiple sets
%         of values at once.
%
%INPUTS: tau A numPointsXnumDims matrix of the coordinates of the points in
%            each dimension in ascending order. The full grid of points
%            in all dimensions could be found using
%            [tau1,tau2,etc.]=ndgrid(tau(1:tauLengths(1),1),tau(1:tauLengths(2),2),...)
%            and the dimensions can be put together into the full set of
%            numDimsXtotalNumPoints points as
%            tauTotal=[tau1(:).';tau2(:).';...]. Note that using meshgrid
%            in 2D/3D will put the points in the wrong order; one has to
%            use the ordering of ndgrid. These values must be real.
% tauLengths A numDimsX1 or a 1XnumDims vector this lists the number of
%            points in each dimension of tau. The total number of sample
%            points is totalNumPoints=prod(tauLengths).
%          y This is the totalNumPointsXnumSets full sets of function
%            values taken at all of the points implied by the vectors in
%            tau for numSets different interpolation problems. The ordering
%            of the points is the same as the ordering of tauTotal, which
%            is described above. These values can be real or complex.
%          k A numDimsX1 or 1XnumDims set of the order of the b-splines in
%            each dimension. The value k-1 is the polynomial order of the
%            approximation being performed. If the order is the same in all
%            dimensions, then a single scalar can be passed.
% numCentralKnots This optional parameter is used if one passes a higher
%            density of points in tau than is needed, because one wants to
%            perform least squares interpolation. This is the numDimsX1 or
%            1XnumDims list of the number of knots in the central region of
%            knots in each dimension. It varies from 0 to n-k. The higher
%            the number, the more knots are used total. A total of
%            numCentralKnots+2*k will be returned by this function, though,
%            as explained below, the first and last knots in each dimension
%            are not used.
%  useSparse The linear system being solved for the coefficients can be
%            quite large, but is not necessarily dense. If useSparse is
%            true, then the matrix for the system will be allocated as a
%            sparse matrix. This can slow down the algorithm, but it can
%            also keep Matlab from running out of memory on large systems.
%            The default if this parameter is omitted or an empty matrix is
%            passed is true.
%
%OUTPUTS: a A hypermatrix with numDim+1 indices containing the set of
%           coefficients for the b-splines covering all of the
%           interpolation intervals with the final index selecting the set
%           if numSets>1. If numSets=1, then the final index is unitary,
%           meaning that there are effectively only numDim dimensions. This
%           and the following outputs can be passed to
%           BSplineInterpValMultiDim to perform interpolation over the
%           region spanned by the tau values.
%         t The maxNumKnotsXnumDims set of knots for the interpolation
%           function. The first and last k(curDim)-1 knots in each
%           dimension are outside of the ends of the region with data or
%           mark the ends of the region with data. The number of rows is
%           the maximum number needed for all of the dimensions. tLength
%           says how many items in each column are actually used. The knots
%           are the same for each set.
%   tLength A numDimX1 vector where tLength(i) says the number of elements
%           in t(:,i). The lengths are the same for each set.
%
%The comments ot the function BSplinePolyFit describe how Equation 7 in
%Chapter IX is used in one dimension. The issue of multiple dimensions is
%addressed in Chapter 17 of [1] with a focus on solving the problem in
%2D. The basic equation to be solved in 2D is given by Equation 11 in
%Chapter 17 of [1] and it can be easily extended to an arbitrary number
%of dimensions. Chapter 17 gives a solution that is more efficient in 2D
%than solving Equation 11. However, we choose to solve the equation to make
%the problem valid for an arbitrary number of dimensions.
%
%The ability to do a least-squares fit involves making fewer inner knots
%than suggested in the text. Rather than making n-k knots, we make
%numCentralKnots. However, this invalidates the approximation in Equation
%10 of Chapter 13 of [1] for selecting knots. The correction applied is to
%make the average taken in the equation to involve more points.
%Additionally, the knots in Chapter 13 are univariate. We make them
%multivariate by taking a tensor product.
%
%EXAMPLE 1:
%Here, we fit a curve to 2D data on an irregular grid. We choose a
%sufficiently high order that it recreates the polynomial exactly and we
%verify this by looking at the maximum error on a fine grid and plotting
%the results.
% f=@(x,y)(x.^4-2*x.^2+x).*(y.^4-2*y.^2+y);
% numDims=2;
% numPointsX=10;
% numPointsY=11;
% tauLengths=[numPointsX;numPointsY];
% tau=zeros(max(tauLengths),numDims);
% tau(1:tauLengths(1),1)=linspace(-1.5,1.5,numPointsX);
% tau(1:tauLengths(2),2)=linspace(-1.5,1.5,numPointsY);
% %Note that using meshgrid would have put the elements in the wong order.
% %ndgrid must be used.
% [tau1,tau2]=ndgrid(tau(1:tauLengths(1),1),tau(1:tauLengths(2),2));
% y=f(tau1,tau2);
%
% k=[5;5];
% [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,y(:),k);
%
% numPointsX=100;
% numPointsY=101;
% pointsX=linspace(-1.5,1.5,numPointsX);
% pointsY=linspace(-1.5,1.5,numPointsY);
% [X,Y]=ndgrid(pointsX,pointsY);
% x=[X(:).';Y(:).'];
%
% zTrue=f(X,Y);
%
% z=BSplineInterpValMultiDim(x,t,tLength,a,k);
% z=reshape(z,[numPointsX,numPointsY]);
%
% figure(1)
% clf
% surface(X,Y,z,'EdgeColor','none')
% title('Interpolated Surface')
%
% figure(2)
% clf
% surface(X,Y,zTrue,'EdgeColor','none')
% title('True Surface')
% max(abs(zTrue(:)-z(:)))
%One sees that the curves are the same and that the difference at any point
%is on the order of finite precision error. This is consistent with the
%fact that the order of the spline fit was chosen to be the same as the
%order of the true polynomial.
%
%EXAMPLE 2:
%Here, we fit a curve to 3D data on an irregular grid. We choose a
%sufficiently high order that it recreates the polynomial exactly and we
%verify this by looking at the maximum error on a fine grid.
% f=@(x,y,z)(x.^4-2*x.^2+x).*(y.^4-2*y.^2+y).*(z.^4-3*z.^2+z);
% numDims=3;
% numPointsX=10;
% numPointsY=11;
% numPointsZ=9;
% tauLengths=[numPointsX;numPointsY;numPointsZ];
% tau=NaN(max(tauLengths),numDims);
% tau(1:tauLengths(1),1)=linspace(-1.5,1.5,numPointsX);
% tau(1:tauLengths(2),2)=linspace(-1.5,1.5,numPointsY);
% tau(1:tauLengths(3),3)=linspace(-1.5,1.5,numPointsZ);
% %Note that using meshgrid would have put the elements in the wong order.
% %ndgrid must be used.
% [tau1,tau2,tau3]=ndgrid(tau(1:tauLengths(1),1),tau(1:tauLengths(2),2),tau(1:tauLengths(3),3));
% y=f(tau1,tau2,tau3);
%
% k=[5;5;5];
% [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,y(:),k);
%
% numPointsX=50;
% numPointsY=51;
% numPointsZ=49;
% pointsX=linspace(-1.5,1.5,numPointsX);
% pointsY=linspace(-1.5,1.5,numPointsY);
% pointsZ=linspace(-1.5,1.5,numPointsZ);
% [X,Y,Z]=ndgrid(pointsX,pointsY,pointsZ);
% x=[X(:).';Y(:).';Z(:).'];
%
% wTrue=f(X,Y,Z);
% w=BSplineInterpValMultiDim(x,t,tLength,a,k);
% max(abs(wTrue(:)-w(:)))
%As was the case with Example 1, the maximum error is on the order of what
%one would expect due to finite precision limitations.
%
%EXAMPLE 3:
%In this case, a complex 2D function is fit and interpolated.
% numPoints=50;
% xMin=0;
% xMax=3;
% yMin=0;
% yMax=3;
% x=linspace(xMin,xMax,numPoints);
% y=linspace(yMin,yMax,numPoints);
% [X,Y]=ndgrid(x,y);
% f=@(x,y)((sin(5*y)+1j*cos(10*x)).*exp(1j*(x.^2-2*y.^2)));
%
% fxy=f(X(:),Y(:));
% k=[5;5];
% tau=[x(:),y(:)];
% tauLengths=[numPoints,numPoints];
% [a,t,tLength]=BSplinePolyFitMultiDim(tau,tauLengths,fxy,k);
%
% %Interpolate on a finer grid.
% numPoints=100;
% x=linspace(xMin,xMax,numPoints);
% y=linspace(xMin,xMax,numPoints);
% [X,Y]=ndgrid(x,y);
%
% fxy=f(X(:),Y(:));
% pts=[X(:).';Y(:).'];
% fxyInterp=reshape(BSplineInterpValMultiDim(pts,t,tLength,a,k),[numPoints,numPoints]);
%
% fxy=reshape(fxy,[numPoints,numPoints]);
%
% figure(1)
% clf
% hold on
% surface(X,Y,real(fxy),'EdgeColor','none')
% title('Real Function Values')
% colorbar()
%
% figure(2)
% clf
% hold on
% surface(X,Y,imag(fxy),'EdgeColor','none')
% title('Imaginary Function Values')
% colorbar()
%
% figure(3)
% clf
% hold on
% surface(X,Y,abs(real(fxyInterp)-real(fxy)),'EdgeColor','none')
% title('Real Interpolation Error')
% colorbar()
%
% figure(4)
% clf
% hold on
% surface(X,Y,abs(imag(fxyInterp)-imag(fxy)),'EdgeColor','none')
% title('Imaginary Interpolation Error')
% colorbar()
%
%REFERENCES:
%[1] C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag,
%    1978.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDims=size(tau,2);
n=tauLengths(:);
numSets=size(y,2);

if(nargin<6)
    useSparse=false;
end

%The total number of points.
NTotal=prod(n);
if(NTotal~=size(y,1))
    error('The dimensionalities of the inputs are not consistent.')
end

if(isscalar(k))
    k=repmat(k,[numDims,1]);
end
k=k(:);

if(nargin<5||isempty(numCentralKnots))
    numCentralKnots=n-k;
elseif(any(numCentralKnots<0)||any(numCentralKnots(:)>n-k))
    error('The value of numCentralKnots is invalid.')
end
numCentralKnots=numCentralKnots(:);

deltaVal=n-k-numCentralKnots;

tLength=numCentralKnots+2*k;
maxTLength=max(tLength);

t=zeros(maxTLength,numDims);
for curDim=1:numDims
    taLengthCur=tauLengths(curDim);
    tauCur=tau(1:taLengthCur,curDim);
    kCur=k(curDim);
    deltaValCur=deltaVal(curDim);
    numCentKnotsCur=numCentralKnots(curDim);

    t(1:kCur,curDim)=tauCur(1);
    t((numCentKnotsCur+kCur+1):(numCentKnotsCur+2*kCur),curDim)=tauCur(n(curDim));

    %The ad-hoc point choice from 10 of Chapter 13, modified to take more
    %points if a least-squares solution is being performed.
    for i=1:numCentKnotsCur
        t(k+i,curDim)=sum(tauCur((i+1):(i+kCur-1+deltaValCur)))/(kCur-1+deltaValCur);
    end
end

%Given the locations of the knots, we build the equations. The a values are
%of the form a(d1,d2,..dnumDims). To access them, we stack them one after
%another in the same order as a(:) would in Matlab.
totalNumA=prod(numCentralKnots+k);

if(useSparse)
    %The amount of memory used will grow as B is filled in.
    B=sparse(NTotal,totalNumA);
else
    B=zeros(NTotal,totalNumA);
end

spanMin=ones(numDims,1);
tIdxVals=k;

%tauIdxVals holds the index of the value of tau used in each dimension
tauIdxVals=ones(numDims,1);
tauVal=zeros(numDims,1);
%Fill in tauVal with the current value.
for curDim=1:numDims
    tauVal(curDim)=tau(1,curDim);
end

%The number of B values produced per call to evalBSplinePolysMultiDim.
numBCurEls=prod(k);

%We have to go through all of the indices.
for curEq=1:NTotal
    BCur=evalBSplinePolysMultiDim(t,tLength,tauVal,k,tIdxVals);
    %We have to put the elements in BCur into the proper spots in B. These
    %spots start at spanMin in each dimension.
    for curEl=1:numBCurEls
        dimIdx=index2NDim(k,curEl);
        %The multidimensional index of the a value by which the current B
        %value is multiplied.
        aDimIdx=spanMin+dimIdx-1;
        aIdx=nDim2Index(n-deltaVal,aDimIdx);

        B(curEq,aIdx)=BCur(curEl);
    end

    if(curEq==NTotal)
        break;
    end

    %Now, we increment tauIdxVals to move onto a new value of tau in each
    %dimension. the values are kept in the tauVal array. Each time a value
    %is incremented, we have to check whether the corresponding entry in
    %tIdxVals and spanMin should be incremented.
    curIdx=1;
    tauIdxVals(curIdx)=tauIdxVals(curIdx)+1;
    while(tauIdxVals(curIdx)>tauLengths(curIdx))
        %If it exceeds k, then go back to 1.
        tauIdxVals(curIdx)=1;
        tauVal(curIdx)=tau(1,curIdx);

        tIdxVals(curIdx)=k(curIdx);
        spanMin(curIdx)=1;

        %Increment tIdxVals and spanMin as necessary to reach the valid
        %region.
        while(tauVal(curIdx)>t(tIdxVals(curIdx)+1,curIdx))
            tIdxVals(curIdx)=tIdxVals(curIdx)+1;
            spanMin(curIdx)=spanMin(curIdx)+1;
        end

        curIdx=curIdx+1;
        tauIdxVals(curIdx)=tauIdxVals(curIdx)+1;
    end

    %Increment tIdxVals and spanMin as necessary to reach the valid
    %region.
    tauVal(curIdx)=tau(tauIdxVals(curIdx),curIdx);
    while(tauVal(curIdx)>t(tIdxVals(curIdx)+1,curIdx))
        tIdxVals(curIdx)=tIdxVals(curIdx)+1;
        spanMin(curIdx)=spanMin(curIdx)+1;
    end
end

a=B\y;
if(numDims>1)
    a=reshape(a,[(numCentralKnots+k)',numSets]);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function[]=makefigs_closedcurves
%MAKEFIGS_CLOSEDCURVES  Makes a sample figure for CLOSEDCURVES.

load qgsnapshot, use qgsnapshot
dx=qgsnapshot.x(2)-qgsnapshot.x(1);
[cv,zeta,N,S,P]=psi2fields(dx,qgsnapshot.psi);
P=frac(P,std(P(:)));

[xc,yc]=closedcurves(x,y,P,-2);
[xp,yp]=closedcurves(x,y,P,-2,'periodic',100,100);
[xpax,ypax,fp]=periodize(100,100,x,y,P);

figure,jpcolor(xpax,ypax,fp),axis equal,axis tight
colormap gray, flipmap, caxis([-5 5]),
vlines(xpax([100 end-100]),'w'),hlines(ypax([100 end-100]),'w')
xtick([-5:1:5]*1000),ytick([-5:1:5]*1000),noxlabels,noylabels
cellplot(xp,yp,'2b'),hold on,cellplot(xc,yc,'r'),fontsize 11
title('Periodized (blue) non-periodized (red) curves for a QG model')
h=colorbar;
h.Label.String='Okubo-Weiss Parameter normalized by its own standard deviation';

%To print
if 0
    currentdir=pwd;
    cd([whichdir('jlab_license') '/figures'])
    print -dpng closedcurves
    crop closedcurves.png
    cd(currentdir)
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function set1001a(m)
% sets parameters of DT3152 according to
% http://cmp.felk.cvut.cz/cmp/hardware/pulnix.html#TM-1001analog

% m=openfg('DT3152');
% setfg(m,'digioconfig',255);
%   setfg(m,'PCMode',1)
%   setfg(m,'PCValue',1)
%   setfg(m,'exposure',0)
   setfg(m,'InputChannel',0)
   setfg(m,'FrameType','NonInterlaced')
   setfg(m,'VideoType','Composite')
   setfg(m,'ClockSource','internal')
   setfg(m,'ClockFreq',20000000)
   setfg(m,'TotalPixPerLine',1272)

   setfg(m,'backporchstart',90)
   setfg(m,'clampstart',165)
   setfg(m,'clampend',170)

   setfg(m,'FirstActivePixel',235)
   setfg(m,'ActivePixelCount',988)

   setfg(m,'TotalLinesperFld',1050)
   setfg(m,'FirstActiveLine',30)
   setfg(m,'ActiveLineCount',1015)

   setfg(m,'FrameLeft',0)
   setfg(m,'HorFrameInc',1)
   setfg(m,'FrameWidth',988)
   setfg(m,'FrameTop',0)
   setfg(m,'VerFrameInc',1)
   setfg(m,'FrameHeight',1015)

%   setfg(m,'ExtOnLoToHi',0)
%   expozice=128+64+32;		% nejdelsi
%   mod16 = 5;
%   mod81 = 0;
%   mod82 = 4;
%   mod83 = 2;
%   mod84 = 6;
%   mod8  = 1;
%   integracestart=0;
%   integracestop=8;
%   integracepovolena =0;
%   integracenepovolena=16;
%   integracezakazana=integracestop+integracenepovolena;
%   setfg(m,'digitalio',mod16+integracezakazana+expozice)

setfg(m, 'SyncMaster', 'Off')
setfg(m, 'SyncValue','VerFreq', 30)		%TM9701
%setfg(m, 'SyncValue','VerFreq', 15)		%TM1001
setfg(m, 'SyncMaster', 'On')

author: Sable

repo: mcbench-benchmarks


text:
function linpack_c_test12 ( )

%*****************************************************************************80
%
%% TEST12 tests CHIFA and CHISL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST12\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  Hermitian matrix (HI):\n' );
  fprintf ( 1, '  CHIFA factors the matrix.\n' );
  fprintf ( 1, '  CHISL solves a linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  seed = 123456789;

  for i = 1 : n
    [ a(i,i), seed ] = r4_uniform_01 ( seed );
    for j = i+1 : n
      [ a(i,j), seed ] = c4_uniform_01 ( seed );
      a(j,i) = conj ( a(i,j) );
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  (%8f  %8f)', real ( a(i,j) ), imag ( a(i,j) ) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Set the values of the right hand side vector B.
%
  [ x, seed ] = c4vec_uniform_01 ( n, seed );

  b(1:n) = a(1:n,1:n) * transpose ( x(1:n) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The right hand side B:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  (%8f  %8f)\n', real ( b(i) ), imag ( b(i) ) );
  end
%
%  Factor the matrix A.
%
  [ a, ipvt, info ] = chifa ( a, lda, n );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  CHIFA returned an error flag INFO = %d', info );
    return
  end
%
%  Solve the system.
%
  b = chisl ( a, lda, n, ipvt, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed                     Exact\n' );
  fprintf ( 1, '  Solution                     Solution\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  (%12f  %12f)  (%12f  %12f)\n', ...
      real ( b(i) ), imag ( b(i) ), real ( x(i) ), imag ( x(i) ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function vrchk(vrep, res, buffer)
% Checks VREP return code. Set buffer to 1 if you are reading from a buffered
% call.

% (C) Copyright Renaud Detry 2013.
% Distributed under the GNU General Public License.
% (See http://www.gnu.org/copyleft/gpl.html)

  if (nargin < 3),
    buffer = false;
  end
  if (nargin < 2),
    error('Missing arguments.');
  end

  expl = 'Undefined error';

  if res == vrep.simx_error_noerror,
    % Nothing to say
    return;
  elseif res == vrep.simx_error_novalue_flag,
    if buffer,
      % No problem to report
      return;
    else
      expl = 'There is no command reply in the input buffer. This should not always be considered as an error, depending on the selected operation mode';
    end
  elseif res == vrep.simx_error_timeout_flag,
    expl = 'The function timed out (probably the network is down or too slow)';
  elseif res == vrep.simx_error_illegal_opmode_flag,
    expl = 'The specified operation mode is not supported for the given function';
  elseif res == vrep.simx_error_remote_error_flag,
    expl = 'The function caused an error on the server side (e.g. an invalid handle was specified)';
  elseif res == vrep.simx_error_split_progress_flag,
    expl = 'The communication thread is still processing previous split command of the same type';
  elseif res == vrep.simx_error_local_error_flag,
    expl = 'The function caused an error on the client side';
  elseif res == vrep.simx_error_initialize_error_flag,
    expl = 'simxStart was not yet called';
  end
  error(sprintf('Remote API function call returned with error code: %d. Explanation: %s.\n', res, expl));

end

author: xuhuairuogu

repo: V-REP-Simulation-Projects

path: github-repos/MATLAB/xuhuairuogu-V-REP-Simulation-Projects/V-REP-Simulation-Projects-841b944af4ea3a8fb250578d36434515f577f411/admittance control(adapted from an example in the book--A Systematic Approach to Learning Robot Programming with ROS)/vrep_matlab_communication/vrchk.m


text:
function [IR, VAR] = VARir(VAR,VARopt)
% =========================================================================
% Compute impulse responses (IRs) for a VAR model estimated with the
% VARmodel.m function. Four identification schemes can be specified:
% zero contemporaneous restrictions, zero long-run restrictions, sign
% restrictions, and external instrumenmts.
% =========================================================================
% [IRF, VAR] = VARir(VAR,VARopt)
% -------------------------------------------------------------------------
% INPUT
%   - VAR: structure, result of VARmodel.m
%   - VARopt: options of the VAR (result of VARmodel.m)
% -------------------------------------------------------------------------
% OUTPUT
%   - IR(:,:,:) : matrix with IRF (H horizons, N variables, N shocks)
%   - VAR: structure including VAR estimation results. Note here that the
%       structure VAR is an output of VARmodel, too. This fucntion adds to
%       VAR some additional results, e.g. VAR.B is the structural impact
%       matrix
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Check inputs
%==========================================================================
if ~exist('VAR','var')
    error('You need to provide VAR structure, result of VARmodel');
end
IV = VAR.IV;
if strcmp(VARopt.ident,'iv')
    if isempty(IV)
        error('You need to provide the data for the instrument in VAR (IV)');
    end
end

%% Retrieve and initialize variables
%==========================================================================
nsteps = VARopt.nsteps;
impact = VARopt.impact;
shut   = VARopt.shut;
recurs = VARopt.recurs;
Fcomp  = VAR.Fcomp;
nvar   = VAR.nvar;
nlag   = VAR.nlag;
sigma  = VAR.sigma;
IR     = nan(nsteps,nvar,nvar);

%% Compute Wold representation
%==========================================================================
% Initialize Wold multipliers
PSI = zeros(nvar,nvar,nsteps);
% Re-write F matrix to compute multipliers
VAR.Fp = zeros(nvar,nvar,nlag);
I = VAR.const+1;
for ii=1:nsteps
    if ii<=nlag
        VAR.Fp(:,:,ii) = VAR.F(:,I:I+nvar-1);
    else
        VAR.Fp(:,:,ii) = zeros(nvar,nvar);
    end
    I = I + nvar;
end
% Compute multipliers
PSI(:,:,1) = eye(nvar);
for ii=2:nsteps
    jj=1;
    aux = 0;
    while jj<ii
        aux = aux + PSI(:,:,ii-jj)*VAR.Fp(:,:,jj);
        jj=jj+1;
    end
    PSI(:,:,ii) = aux;
end
% Update VAR with Wold multipliers
VAR.PSI = PSI;

%% Identification: Recover B matrix
%==========================================================================
% B matrix is recovered with Cholesky decomposition
if strcmp(VARopt.ident,'short')
    [out, chol_flag] = chol(sigma);
    if chol_flag~=0; error('VCV is not positive definite'); end
    B = out';
% B matrix is recovered with Cholesky on cumulative IR to infinity
elseif strcmp(VARopt.ident,'long')
    Finf_big = inv(eye(length(Fcomp))-Fcomp);
    Finf = Finf_big(1:nvar,1:nvar);
    D  = chol(Finf*sigma*Finf')';
    B = Finf\D;
% B matrix is recovered with SR.m
elseif strcmp(VARopt.ident,'sign')
    if isempty(VAR.B)
        error('You need to provide the B matrix with SR.m and/or SignRestrictions.m')
    else
        B = VAR.B;
    end
% B matrix is recovered with external instrument IV
elseif strcmp(VARopt.ident,'iv')
    % Recover residuals (first variable is the one to be instrumented - order matters!)
    up = VAR.resid(:,1);     % residuals to be instrumented
    uq = VAR.resid(:,2:end); % residulas for second stage

    % Make sample of IV comparable with up and uq
    [aux, fo, lo] = CommonSample([up IV(VAR.nlag+1:end,:)]);
    p = aux(:,1);
    q = uq(end-length(p)+1:end,:); pq = [p q];
    Z = aux(:,2:end);

    % Run first stage regression and fitted
    FirstStage = OLSmodel(p,Z);
    p_hat = FirstStage.yhat;

    % Recover first column of B matrix with second stage regressions
    Biv(1,1) = 1;  % Start with impact IR normalized to 1
    sqsp = zeros(size(q,2),1);
    for ii=2:nvar
        SecondStage = OLSmodel(q(:,ii-1),p_hat);
        Biv(ii,1) = SecondStage.beta(2);
        sqsp(ii-1) = SecondStage.beta(2);
    end
    % Update size of the shock (ftn 4 of Gertler and Karadi (2015))
    sigma_b = (1/(length(pq)-VAR.ntotcoeff))*...
        (pq-repmat(mean(pq),size(pq,1),1))'*...
        (pq-repmat(mean(pq),size(pq,1),1));
    s21s11 = sqsp;
    S11 = sigma_b(1,1);
    S21 = sigma_b(2:end,1);
    S22 = sigma_b(2:end,2:end);
    Q = s21s11*S11*s21s11'-(S21*s21s11'+s21s11*S21')+S22;
    sp = sqrt(S11-(S21-s21s11*S11)'*(Q\(S21-s21s11*S11)));
    % Rescale Biv vector
    Biv = Biv*sp;
    B = zeros(nvar,nvar);
    B(:,1) = Biv;
% If none of the above, you've done somerthing wrong :)
else
    disp('---------------------------------------------')
    disp('Identification incorrectly specified.')
    disp('Choose one of the following options:');
    disp('- short: zero contemporaneous restrictions');
    disp('- long:  zero long-run restrictions');
    disp('- sign:  sign restrictions');
    disp('- iv:  external instrument');
    disp('---------------------------------------------')
    error('ERROR. See details above');
end

%% Compute the impulse response
%==========================================================================
for mm=1:nvar
    % Set to zero a row of the companion matrix if "shut" is selected
    if shut~=0
        Fcomp(shut,:) = 0;
    end
    % Initialize the impulse response vector
    response = zeros(nvar, nsteps);
    % Create the impulse vector
    impulse = zeros(nvar,1);
    % Set the size of the shock
    if impact==0
        impulse(mm,1) = 1; % one stdev shock
    elseif impact==1
        impulse(mm,1) = 1/B(mm,mm); % unitary shock
    else
        error('Impact must be either 0 or 1');
    end
    % First period impulse response (=impulse vector)
    response(:,1) = B*impulse;
    % Shut down the response if "shut" is selected
    if shut~=0
        response(shut,1) = 0;
    end
    % Recursive computation of impulse response
    if strcmp(recurs,'wold')
        for kk = 2:nsteps
            response(:,kk) = PSI(:,:,kk)*B*impulse;
        end
    elseif strcmp(recurs,'comp')
        for kk = 2:nsteps
            FcompN = Fcomp^(kk-1);
            response(:,kk) = FcompN(1:nvar,1:nvar)*B*impulse;
        end
    end
    IR(:,:,mm) = response';
end
% Update VAR with structural impact matrix
VAR.B = B;
if strcmp(VARopt.ident,'iv')
    VAR.FirstStage = FirstStage;
    VAR.sigma_b = sigma_b;
    VAR.Biv = Biv;
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function[]=flipx
% FLIPX    Flips the direction of the x-axis
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2002, 2004 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(get(gca,'xdir'),'normal')
  set(gca,'xdir','reverse')
  breverse=1;
elseif strcmpi(get(gca,'xdir'),'reverse')
  set(gca,'xdir','normal')
  breverse=0;
end

%Update for compatibility with ZOOL
if isappdata(gcf,'xlimitslock')
   if getappdata(gcf,'xlimitslock')
      h=gca;
      hall=axeshandles;
      for i=1:length(hall)
         if breverse
             set(hall(i),'xdir','reverse')
         else
             set(hall(i),'xdir','normal')
         end
      end
   end
end


text:
function uiimport(varargin)
% overload standard MATLAB uiimport

if nargin && ischar(varargin{1})
  [pathstr, name, ext] = fileparts(varargin{1});
  if any(strcmpi(getMTEXpref('poleFigureExtensions'),ext))
    import_wizard('PoleFigure',varargin{:});
  elseif any(strcmpi(getMTEXpref('EBSDExtensions'),ext))
    import_wizard('EBSD',varargin{:});
  else
    old_uiimport(varargin{:});
  end
else
  old_uiimport(varargin{:});
end

function old_uiimport(varargin)

s = which('uiimport','-all');
pathstr = fileparts(s{end});
opwd = pwd;
cd(pathstr);
uiimport(varargin{:});
cd(opwd);

author: mtex-toolbox

repo: mtex