text:
function [obj_subset, to_extract] = select_atlas_subset(obj, varargin)
% Select a subset of regions in an atlas by name or integer code, with or without collapsing regions together
%
% Select a subset of regions in an atlas using:
% - a cell array of one or more strings to search for in labels
% - a vector of one or more integers indexing regions
%
% Output:
%       a new object, obj_subset
%       to_extract, a logical vector indicating which regions were selected
%
% obj_subset = select_atlas_subset(obj, varargin)
%
% Options:
% 'flatten' : Flatten integer vector in .dat to a single 1/0 mask. Good for
% creating masks from combinations of regions, or adding a set to another
% atlas as a single atlas region. .probability_maps reduced to single max
% map
%
% 'labels_2' : If you enter any field name in the object, e.g., labels_2,
% the function will search for keyword matches here instead of in obj.labels.
%
% Examples:
%
% atlasfile = which('Morel_thalamus_atlas_object.mat');
% load(atlasfile)
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, [1 3]);
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, [1 3], {'VPL'});
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'VPL'});             % Sensory VPL thalamus, when using morel atlas
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'VPL' 'VPM' 'VPI'}); % Sensory thalamus, both lat and medial
%
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'Hb'});              % Habenula
%
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'CL' 'CeM' 'CM' 'Pf'});  % Intralaminar
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'Pv' 'SPf'});           % Midline group
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'CL' 'CeM' 'CM' 'Pf' 'Pv' 'SPf'}); % Intralaminar and Midline group
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'MD'});              % Mediodorsal nuc.
%
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'Pu'});              % Pulvinar
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'LGN'});             % LGN
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'MGN'});             % MGN
%
% Create a mask from a set of regions:
% obj_subset = select_atlas_subset(atlas_obj, {'CL' 'CeM' 'CM' 'Pf'}, 'flatten');  % Intralaminar
% r = atlas2region(obj_subset);
% orthviews(r)
%
% Load a canonical atlas and select only the default mode, limbic, and brainstem regions
% based on the .labels_2 property
% -------------------------------------------------------
% atlas_obj = load_atlas('canlab2018_2mm');
% [obj_subset, to_extract] = select_atlas_subset(atlas_obj, {'Cerebellum', 'Def'}, 'labels_2');
% montage(obj_subset);
%
% see also: image_vector.select_voxels_by_value

% Programmers' notes:
% Created by tor wager
% Edited by tor, 12/2019, to use any field to select labels; and update
% auxiliary label fields too.
%

% -------------------------------------------------------------------------
% DEFAULTS AND INPUTS
% -------------------------------------------------------------------------

strings_to_find = [];
integers_to_find = [];
doflatten = false;

% for entry of optional field to search for keywords
myfields = fieldnames(obj);
mylabelsfield = 'labels';

% optional inputs with default values
for i = 1:length(varargin)

    if iscell(varargin{i})
        strings_to_find = varargin{i};

    elseif isnumeric(varargin{i})
        integers_to_find = varargin{i};

    elseif ischar(varargin{i})
        switch varargin{i}

             case 'flatten', doflatten = true;

                 %             case 'xxx', xxx = varargin{i+1}; varargin{i+1} = [];

            otherwise

                if any(strcmp(varargin{i}, myfields))
                    mylabelsfield = varargin{i};

                else

                    warning(['Unknown input string option:' varargin{i} '. Assuming it might be an atlas label. Place atlas labels in a cell array']);
                    strings_to_find{1} = varargin{i};

                end
        end
    end
end

% -------------------------------------------------------------------------
% INIT
% -------------------------------------------------------------------------

k = num_regions(obj);
to_extract = false(1, k);

% -------------------------------------------------------------------------
% FIND BY STRING
% -------------------------------------------------------------------------

for i = 1:length(strings_to_find)

    % Find which names match
    wh = ~cellfun(@isempty, strfind(obj.(mylabelsfield), strings_to_find{i}));
    if strmatch(mylabelsfield, 'label_descriptions')
        wh=wh';
    end

    to_extract = to_extract | wh;

end

% -------------------------------------------------------------------------
% FIND BY NUMBERS
% -------------------------------------------------------------------------

to_extract(integers_to_find) = true;

if ~any(to_extract)
    error('No regions identified to extract.');
end

% -------------------------------------------------------------------------
% Extract
% -------------------------------------------------------------------------

obj_subset = obj;

% labels and names

obj_subset.labels = obj_subset.labels(to_extract);

% obj_subset.label_descriptions = obj_subset.label_descriptions(to_extract);

% if ~isempty(obj_subset.image_names) && size(obj_subset.image_names, 1) == k
%     obj_subset.image_names = obj_subset.image_names(to_extract, :);
% end

my_strings = {'image_names' 'fullpath'};

for i = 1:length(my_strings)

    if  ~isempty(obj_subset.(my_strings{i})) && size(obj_subset.(my_strings{i}), 1) == k
        obj_subset.(my_strings{i}) = obj_subset.(my_strings{i})(to_extract, :);
    end

end

my_strings = {'label_descriptions' 'image_names' 'fullpath' 'labels_2' 'labels_3' 'labels_4' 'labels_5'};

for i = 1:length(my_strings)
    if  ~isempty(obj_subset.(my_strings{i})) && length(obj_subset.(my_strings{i})) == k
        obj_subset.(my_strings{i}) = obj_subset.(my_strings{i})(to_extract);
    end

end

if ~isempty(obj.probability_maps) && size(obj.probability_maps, 2) == k  % valid p maps

    obj_subset.probability_maps(:, ~to_extract) = [];

    obj_subset = probability_maps_to_region_index(obj_subset);

else % must use .dat vector with integers

    % rebuild integer index

    dat = zeros(size(obj_subset.dat));
    wh_cols = find(to_extract);

    for i = 1:length(wh_cols)

        wh_vals = obj_subset.dat == wh_cols(i);

        dat(wh_vals) = i;

    end

    obj_subset.dat = dat;

end

if doflatten
    % Flatten into mask
    % -------------------------------------------------------------------
    obj_subset.dat = single(obj_subset.dat > 0);

    obj_subset.probability_maps = max(obj_subset.probability_maps, [], 2);

    % Add labels for combined mask

    mylabels = [];
    if iscell(strings_to_find)
        mylabels = cat(2, strings_to_find{:});
    end

    if ~isempty(integers_to_find)
        morelabels = cat(2, obj.labels{integers_to_find});
        mylabels = [mylabels morelabels];
    end

    obj_subset.labels = {mylabels};

    obj_subset.label_descriptions = {sprintf('Combined subset: %s', mylabels)};

        % Make sure they are length of prob maps, to conform to standard atlas
    % object
%     np = size(obj_subset.probability_maps, 2);
%     if np
%         obj_subset.labels = repmat(obj_subset.labels, 1, np);
%         obj_subset.label_descriptions = repmat(obj_subset.label_descriptions, np, 1);
%     end

end % flatten

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [pnt, tri] = triangulate_seg(seg, npnt, origin)

% TRIANGULATE_SEG constructs a triangulation of the outer surface of a
% segmented volume. It starts at the center of the volume and projects the
% vertices of an evenly triangulated sphere onto the outer surface. The
% resulting surface is star-shaped from the origin of the sphere.
%
% Use as
%   [pnt, tri] = triangulate_seg(seg, npnt, origin)
%
% Input arguments:
%  seg    = 3D-matrix (boolean) containing segmented volume. If not boolean
%           seg = seg(~=0);
%  npnt   = requested number of vertices
%  origin = 1x3 vector specifying the location of the origin of the sphere
%           in voxel indices. This argument is optional. If undefined, the
%           origin of the sphere will be in the centre of the volume.
%
% Output arguments:
%  pnt = Nx3 matrix of vertex locations
%  tri = Mx3 matrix of triangles
%
% Seg will be checked for holes, and filled if necessary. Also, seg will be
% checked to consist of a single boolean blob. If not, only the outer surface
% of the largest will be triangulated. SPM is used for both the filling and
% checking for multiple blobs.
%
% See also KSPHERE

% Copyright (C) 2005-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

% impose it to be boolean
seg = (seg~=0);
dim = size(seg);
len = ceil(sqrt(sum(dim.^2))/2);

if ~any(seg(:))
  ft_error('the segmentation is empty')
end

% define the origin if it is not provided in the input arguments
if nargin<3
  origin(1) = dim(1)/2;
  origin(2) = dim(2)/2;
  origin(3) = dim(3)/2;
end

% ensure that the seg consists of only one filled blob.
% if not filled: throw a warning and fill
% if more than one blob: throw a warning and use the biggest

% look for holes
seg = volumefillholes(seg);

% ensure that SPM is available, needed for spm_bwlabel
ft_hastoolbox('spm8up', 3) || ft_hastoolbox('spm2', 1);

% look for >1 blob
[lab, num] = spm_bwlabel(double(seg), 26);
if num>1,
  ft_warning('the segmented volume consists of more than one compartment, using only the biggest one for the segmentation');

  for k = 1:num
    n(k) = sum(lab(:)==k);
  end
  [m,ix] = max(n);
  seg(lab~=ix) = false;
end

% start with a unit sphere with evenly distributed vertices
[pnt, tri] = mesh_sphere(npnt, 'ksphere');

ishollow = false;

for i=1:npnt
  % construct a sampled line from the center of the volume outward into the direction of the vertex
  lin = (0:0.5:len)' * pnt(i,:);
  lin(:,1) = lin(:,1) + origin(1);
  lin(:,2) = lin(:,2) + origin(2);
  lin(:,3) = lin(:,3) + origin(3);
  % round the sampled line towards the nearest voxel indices, which allows
  % a quick nearest-neighbour interpolation/lookup
  lin = round(lin);
  % exclude indices that do not lie within the volume
  sel = lin(:,1)<1 | lin(:,1)>dim(1) | ...
        lin(:,2)<1 | lin(:,2)>dim(2) | ...
        lin(:,3)<1 | lin(:,3)>dim(3);
  lin = lin(~sel,:);
  sel = sub2ind(dim, lin(:,1), lin(:,2), lin(:,3));

  % interpolate the segmented volume along the sampled line
  int = seg(sel);

  % the value along the line is expected to be 1 at first and then drop to 0
  % anything else suggests that the segmentation is hollow
  ishollow = any(diff(int)==1);

  % find the last sample along the line that is inside the segmentation
  sel = find(int, 1, 'last');
  % this is a problem if sel is empty. If so, use the edge of the volume
  if ~isempty(sel)
    % take the last point inside and average with the first point outside
    pnt(i,:) = lin(sel,:);
  else
    % take the edge
    pnt(i,:) = lin(end,:);
  end
end

if ishollow
  % this should not have hapened, especially not after filling the holes
  ft_warning('the segmentation is not star-shaped, please check the surface mesh');
end

% undo the shift of the origin from where the projection is done
% pnt(:,1) = pnt(:,1) - origin(1);
% pnt(:,2) = pnt(:,2) - origin(2);
% pnt(:,3) = pnt(:,3) - origin(3);

% fast unconditional re-implementation of the standard MATLAB function
function [s] = sub2ind(dim, i, j, k)
s = i + (j-1)*dim(1) + (k-1)*dim(1)*dim(2);

author: spm


text:
function a = real(a)
%REAL         Real part of Taylor
%
%   c = real(a)
%
%

% written  05/21/09     S.M. Rump
%

  a.t = real(a.t);

repo: OpenMAS


text:
function [sbc, fpe, logdp, np] = ARFIT_arord(R, m, mcor, ne, pmin, pmax)
%ARORD	Evaluates criteria for selecting the order of an AR model.
%
%  [SBC,FPE]=ARORD(R,m,mcor,ne,pmin,pmax) returns approximate values
%  of the order selection criteria SBC and FPE for models of order
%  pmin:pmax. The input matrix R is the upper triangular factor in the
%  QR factorization of the AR model; m is the dimension of the state
%  vectors; the flag mcor indicates whether or not an intercept vector
%  is being fitted; and ne is the number of block equations of size m
%  used in the estimation. The returned values of the order selection
%  criteria are approximate in that in evaluating a selection
%  criterion for an AR model of order p < pmax, pmax-p initial values
%  of the given time series are ignored.
%
%  ARORD is called by ARFIT.
%
%  See also ARFIT, ARQR.

%  For testing purposes, ARORD also returns the vectors logdp and np,
%  containing the logarithms of the determinants of the (scaled)
%  covariance matrix estimates and the number of parameter vectors at
%  each order pmin:pmax.

%  Modified 17-Dec-99
%  Author: Tapio Schneider
%          tapio@gps.caltech.edu

  imax 	  = pmax-pmin+1;        % maximum index of output vectors

  % initialize output vectors
  sbc     = zeros(1, imax);     % Schwarz's Bayesian Criterion
  fpe     = zeros(1, imax);     % log of Akaike's Final Prediction Error
  logdp   = zeros(1, imax);     % determinant of (scaled) covariance matrix
  np      = zeros(1, imax);     % number of parameter vectors of length m
  np(imax)= m*pmax+mcor;

  % Get lower right triangle R22 of R:
  %
  %   | R11  R12 |
  % R=|          |
  %   | 0    R22 |
  %
  R22     = R(np(imax)+1 : np(imax)+m, np(imax)+1 : np(imax)+m);

  % From R22, get inverse of residual cross-product matrix for model
  % of order pmax
  invR22  = inv(R22);
  Mp      = invR22*invR22';

  % For order selection, get determinant of residual cross-product matrix
  %       logdp = log det(residual cross-product matrix)
  logdp(imax) = 2.*log(abs(prod(diag(R22))));

  % Compute approximate order selection criteria for models of
  % order pmin:pmax
  i = imax;
  for p = pmax:-1:pmin
    np(i)      = m*p + mcor;	% number of parameter vectors of length m
   if p < pmax
      % Downdate determinant of residual cross-product matrix
      % Rp: Part of R to be added to Cholesky factor of covariance matrix
      Rp       = R(np(i)+1:np(i)+m, np(imax)+1:np(imax)+m);

      % Get Mp, the downdated inverse of the residual cross-product
      % matrix, using the Woodbury formula
      L        = chol(eye(m) + Rp*Mp*Rp')';
      N        = L \ Rp*Mp;
      Mp       = Mp - N'*N;

      % Get downdated logarithm of determinant
      logdp(i) = logdp(i+1) + 2.* log(abs(prod(diag(L))));
   end

   % Schwarz's Bayesian Criterion
   sbc(i) = logdp(i)/m - log(ne) * (ne-np(i))/ne;

   % logarithm of Akaike's Final Prediction Error
   fpe(i) = logdp(i)/m - log(ne*(ne-np(i))/(ne+np(i)));

   % Modified Schwarz criterion (MSC):
   % msc(i) = logdp(i)/m - (log(ne) - 2.5) * (1 - 2.5*np(i)/(ne-np(i)));

   i      = i-1;                % go to next lower order
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [X,Y,vals,labI]=mp_azim(optn,varargin)
% MP_AZIM Azimuthal projections
%           This function should not be used directly; instead it is
%           is accessed by various high-level functions named M_*.

% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997
%
%         13/5/97 - Added satellite perspective
%          1/6/97 - Another stab at removing some /0 errors.
%         10/8/00 - Rotation for projections?
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
% Mathematical formulas for the projections and their inverses are taken from
%
%      Snyder, John P., Map Projections used by the US Geological Survey,
%      Geol. Surv. Bull. 1532, 2nd Edition, USGPO, Washington D.C., 1983.
%
% These are azimuthal projections, best suited for circular areas. The
% stereographic is commonly used for polar regions.
%   Stereographic  - conformal
%   Orthographic   - neither conformal nor equal-area, but looks like globe
%                    with viewpoint at infinity.
%   Azimuthal Equal-area  - equal area, but not conformal (by Lambert)
%   Azimuthal Equidistant - distance and direction from center are true
%   Gnomonic       - all great circles are straight lines.
%   Satellite      - a perspective view from a finite distance

global MAP_PROJECTION MAP_VAR_LIST

name={'Stereographic','Orthographic','Azimuthal Equal-area','Azimuthal Equidistant','Gnomonic','Satellite'};

pi180=pi/180;

switch optn

  case 'name'
     X=name;

  case {'usage','set'}
     X=char({['     ''' varargin{1} ''''],...
              '     <,''lon<gitude>'',center_long>',...
              '     <,''lat<itude>'', center_lat>',...
              '     <,''rad<ius>'', ( degrees | [longitude latitude] ) | ''alt<itude>'', alt_frac >',...
              '     <,''rec<tbox>'', ( ''on'' | ''off'' | ''circle'' )>',...
	      '     <,''rot<angle>'', degrees CCW>'});

  case 'get'

     X=char([' Projection: ' MAP_PROJECTION.name '  (function: ' MAP_PROJECTION.routine ')'],...
            [' center longitude: ' num2str(MAP_VAR_LIST.ulong) ],...
            [' center latitude: ' num2str(MAP_VAR_LIST.ulat) ],...
            [' radius/altitude : ' num2str(MAP_VAR_LIST.uradius) ],...
            [' Rectangular border: ' MAP_VAR_LIST.rectbox ],...
	    [' Rotation angle: ' num2str(MAP_VAR_LIST.rotang) ]);

  case 'initialize'

    MAP_VAR_LIST=[];
    MAP_PROJECTION.name=varargin{1};
    MAP_VAR_LIST.ulong=0;
    MAP_VAR_LIST.ulat=60;
    MAP_VAR_LIST.rectbox='circle';
    MAP_VAR_LIST.uradius=90;
    MAP_VAR_LIST.rotang=0;
    MAP_VAR_LIST.ellipsoid='normal';
    k=2;
    while k<length(varargin)
      switch varargin{k}(1:3)
         case 'lon'
           MAP_VAR_LIST.ulong=varargin{k+1}(:)';
         case 'lat'
           MAP_VAR_LIST.ulat=varargin{k+1}(:)';
         case {'rad','alt'}
           MAP_VAR_LIST.uradius=varargin{k+1};
         case 'rec'
           MAP_VAR_LIST.rectbox=varargin{k+1};
	 case 'rot'
	   MAP_VAR_LIST.rotang=varargin{k+1};
         otherwise
           disp(['Unknown option: ' varargin{k}]);
      end
      k=k+2;
    end
    if strcmp(MAP_VAR_LIST.rectbox,'off'), MAP_VAR_LIST.rectbox='circle'; end

    MAP_VAR_LIST.rlong=MAP_VAR_LIST.ulong*pi180;
    MAP_VAR_LIST.rlat=MAP_VAR_LIST.ulat*pi180;

    % Compute the various limits
    % For the perspective viewpoints, we can't go *quite* to the horizon because this causes
    % problems in the clipping later on - but we can get pretty close (within .995) of it.
    % This is a fudge factor that can probably be changed if we increase the number of points
    % in grid lines...

    if length(MAP_VAR_LIST.uradius)==1
      if strcmp(MAP_PROJECTION.name,name{2}) && abs(MAP_VAR_LIST.uradius-90)<.2
        MAP_VAR_LIST.radius=89.8;
      elseif strcmp(MAP_PROJECTION.name,name{5})
        MAP_VAR_LIST.radius=min(80,MAP_VAR_LIST.uradius);
      elseif strcmp(MAP_PROJECTION.name,name{6})
        MAP_VAR_LIST.radius=acos(1/(1+MAP_VAR_LIST.uradius))/pi180*.98; % uradius is the height fraction here
      else
        MAP_VAR_LIST.radius=MAP_VAR_LIST.uradius;
      end
      rradius=MAP_VAR_LIST.radius*pi180;
    else
      % do some sperical trig
      edge=MAP_VAR_LIST.uradius*pi180 - [MAP_VAR_LIST.rlong 0];
      cosc=sin(MAP_VAR_LIST.rlat)*sin(edge(2))+cos(MAP_VAR_LIST.rlat)*cos(edge(2))*cos(edge(1));
      sinc=sqrt( ( cos(edge(2))*sin(edge(1)))^2 + ...
                 (cos(MAP_VAR_LIST.rlat)*sin(edge(2))-sin(MAP_VAR_LIST.rlat)*cos(edge(2))*cos(edge(1)))^2);
      rradius=atan2(sinc,cosc);
      MAP_VAR_LIST.radius=rradius/pi180;
    end
    MAP_VAR_LIST.cosradius=cos(rradius);

    switch MAP_PROJECTION.name
      case name(1)
        MAP_VAR_LIST.rhomax=2*tan(rradius/2);
      case name(2)
        MAP_VAR_LIST.rhomax=sin(rradius);
      case name(3)
        MAP_VAR_LIST.rhomax=2*sin(rradius/2);
      case name(4)
        MAP_VAR_LIST.rhomax=rradius;
      case name(5)
        MAP_VAR_LIST.rhomax=tan(rradius);
      case name(6)
        MAP_VAR_LIST.rhomax=sin(rradius)/(1+(1-cos(rradius))/MAP_VAR_LIST.uradius);
    end

    if strcmp(MAP_VAR_LIST.rectbox,'on')
      if length(MAP_VAR_LIST.uradius)==1
        MAP_VAR_LIST.xlims=[-1 1]/sqrt(2)*MAP_VAR_LIST.rhomax;
        MAP_VAR_LIST.ylims=[-1 1]/sqrt(2)*MAP_VAR_LIST.rhomax;
      else
        [X,Y]=mp_azim('ll2xy',MAP_VAR_LIST.uradius(1),MAP_VAR_LIST.uradius(2),'clip','off');
        MAP_VAR_LIST.xlims=[-abs(X) abs(X)];
        MAP_VAR_LIST.ylims=[-abs(Y) abs(Y)];
      end
    else
      MAP_VAR_LIST.xlims=[-MAP_VAR_LIST.rhomax MAP_VAR_LIST.rhomax];
      MAP_VAR_LIST.ylims=[-MAP_VAR_LIST.rhomax MAP_VAR_LIST.rhomax];
    end

    mu_util('lllimits');

  case 'll2xy'

    long=varargin{1}*pi180-MAP_VAR_LIST.rlong;
    lat=varargin{2}*pi180;
    vals=zeros(size(long));

    pi180=pi/180;
    cosc     =sin(MAP_VAR_LIST.rlat)*sin(lat)+cos(MAP_VAR_LIST.rlat)*(cos(lat).*cos(long));
    sinAzsinc=sin(long).*cos(lat);
    cosAzsinc=cos(MAP_VAR_LIST.rlat)*sin(lat)-sin(MAP_VAR_LIST.rlat)*(cos(lat).*cos(long));
    sinc=sqrt(sinAzsinc.^2+cosAzsinc.^2);

    switch MAP_PROJECTION.name
      case name(1)
        cosc(cosc==-1)=-1+eps;
        rho=2*sinc./(1+cosc);  % = 2*tan(c/2)
      case name(2)
        rho=sinc;   % = sinc
      case name(3)
        cosc(cosc==-1)=-1+eps;
        rho=sqrt(2)*sinc./sqrt(1+cosc);   % = 2*sin(c/2)
      case name(4)
        rho=atan2(sinc,cosc); % = c
      case name(5)
        rho=sinc./cosc; % = tan(c)
      case name(6)
        rho=sinc./(1+(1-cosc)/MAP_VAR_LIST.uradius); %
    end

    sinc(sinc==0)=eps;
    Az=(sinAzsinc+sqrt(-1)*cosAzsinc)./sinc;
    Az(abs(Az)==0)=-1;

    % Clip out-of-range values. We test against cos(c) (where c is the angular
    % distance from map center) rather than directly against rhomax, because
    % in the orthographic map rho->0 for points on the other side of the
    % globe whereas c does not!

    % Also, we clip on rho even if we later clip on X/Y because in some projections (e.g. the
    % orthographic) the X/Y locations wrap back.
    if ~strcmp(varargin{4},'off')
        vals = vals | cosc<=MAP_VAR_LIST.cosradius+eps*10;
        [rho,Az]=mu_util('clip',varargin{4},rho,MAP_VAR_LIST.rhomax,cosc<MAP_VAR_LIST.cosradius,Az);
        Az=Az./abs(Az);
    end

     X=rho.*real(Az*exp(i*pi180*MAP_VAR_LIST.rotang));
     Y=rho.*imag(Az*exp(i*pi180*MAP_VAR_LIST.rotang));

    if strcmp(MAP_VAR_LIST.rectbox,'on')  && ~strcmp(varargin{4},'off')
        vals= vals | X<=MAP_VAR_LIST.xlims(1)+eps*10 | X>=MAP_VAR_LIST.xlims(2)-eps*10 | ...
                     Y<=MAP_VAR_LIST.ylims(1)+eps*10 | Y>=MAP_VAR_LIST.ylims(2)-eps*10;
        [X,Y]=mu_util('clip',varargin{4},X,MAP_VAR_LIST.xlims(1),X<MAP_VAR_LIST.xlims(1) | isnan(X),Y);
        [X,Y]=mu_util('clip',varargin{4},X,MAP_VAR_LIST.xlims(2),X>MAP_VAR_LIST.xlims(2) | isnan(X),Y);
        [Y,X]=mu_util('clip',varargin{4},Y,MAP_VAR_LIST.ylims(1),Y<MAP_VAR_LIST.ylims(1) | isnan(Y),X);
        [Y,X]=mu_util('clip',varargin{4},Y,MAP_VAR_LIST.ylims(2),Y>MAP_VAR_LIST.ylims(2) | isnan(Y),X);
    end

  case 'xy2ll'

    rho=sqrt(varargin{1}.^2+varargin{2}.^2);
    Z=exp(i*(atan2(varargin{2},varargin{1})-MAP_VAR_LIST.rotang*pi180));
    V1=rho.*real(Z);
    V2=rho.*imag(Z);

    ir=rho==0;  % To prevent /0 warnings when rho is 0
    rho(ir)=eps;

    switch MAP_PROJECTION.name
      case name(1)
        c=2*atan(rho/2);
      case name(2)
        c=asin(rho);
        c(abs(rho)>1.0)=NaN;  % points outside the map
      case name(3)
        c=2*asin(rho/2);
         c(abs(rho)>2.0)=NaN;  % points outside the map
     case name(4)
        c=rho;
      case name(5)
        c=atan(rho);
      case name(6)
        arg1=(MAP_VAR_LIST.uradius+1)./sqrt(1+(MAP_VAR_LIST.uradius./rho).^2);
        c=asin(arg1) - atan(rho/MAP_VAR_LIST.uradius);
        c(arg1>1.0)=NaN;
    end
    c(ir)=eps; % we offset this slightly so that the correct limit is achieved in the
               % division below:

%    Y=(asin(cos(c)*sin(MAP_VAR_LIST.rlat) + ...
%            cos(MAP_VAR_LIST.rlat)*sin(c).*varargin{2}./rho))/pi180;
%
%    switch MAP_VAR_LIST.ulat,
%      case 90,
%        X=(MAP_VAR_LIST.rlong+atan2(varargin{1},-varargin{2}))/pi180;
%      case -90,
%        X=(MAP_VAR_LIST.rlong+atan2(varargin{1},varargin{2}))/pi180;
%      otherwise
%        X=(MAP_VAR_LIST.rlong+atan2( varargin{1}.*sin(c), ...
%          cos(MAP_VAR_LIST.rlat)*cos(c).*rho - sin(MAP_VAR_LIST.rlat)*varargin{2}.*sin(c) ) )/pi180;
%     end;

   % Can be problem if the argument is slightly larger than 1 - then the asin
   % returns a complex number.
   arg=cos(c)*sin(MAP_VAR_LIST.rlat) + ...
            cos(MAP_VAR_LIST.rlat)*sin(c).*V2./rho;
    arg=min(max(arg,-1),1);

    Y=(asin(arg))/pi180;

    switch MAP_VAR_LIST.ulat
      case 90
        X=(MAP_VAR_LIST.rlong+atan2(V1,-V2))/pi180;
      case -90
        X=(MAP_VAR_LIST.rlong+atan2(V1,V2))/pi180;
      otherwise
        X=(MAP_VAR_LIST.rlong+atan2( V1.*sin(c), ...
          cos(MAP_VAR_LIST.rlat)*cos(c).*rho - sin(MAP_VAR_LIST.rlat)*V2.*sin(c) ) )/pi180;
    end

  case 'xgrid'

    [X,Y,vals,labI]=mu_util('xgrid',MAP_VAR_LIST.longs,MAP_VAR_LIST.lats,varargin{1},31,varargin{2:3});

  case 'ygrid'

    [X,Y,vals,labI]=mu_util('ygrid',MAP_VAR_LIST.lats,MAP_VAR_LIST.longs,varargin{1},91,varargin{2:3});

  case 'box'

     [X,Y]=mu_util('box',31);

end

author: goGPS-Project


text:
function [TAR,FAR,threshold] = rocplot(match, nonmatch)
%TAR=Ver

%match=-match;
%nonmatch=-nonmatch;

 R1=min([match; nonmatch]);
 R2=max([nonmatch;match]);

a=100/(R2-R1);
b=-a*R1;

matchN=a*match+b;
matchN(matchN<0)=0;
nonmatch_score=a*nonmatch+b;
TAR=[1:length(matchN)]'/length(matchN);
FRR=1-TAR;

nonmatchL=zeros(size(matchN));
smatch=sort(matchN);
for i=1:length(nonmatch_score)
    low=1;
    high=length(smatch);
	while low <=high
		mid=ceil((low+high)*.5);
%[nonmatch_score(i),smatch(mid)]
		if nonmatch_score(i)>smatch(mid) low=mid+1;
		else high=mid-1;
		end
	end

 %   [val,idx]=min((nonmatch_score(i)>smatch));
 %    nonmatchL(idx)=nonmatchL(idx)+1;
	if low > length(smatch)
	low= length(smatch);
	end
    nonmatchL(low)=nonmatchL(low)+1;
end
cnonmatchL=cumsum(nonmatchL);
FAR=cnonmatchL/cnonmatchL(end);
iind=find(FAR>=0.1);
VER1= TAR(iind(1));
iind=find(FAR>=0.01);
VER01= TAR(iind(1));
iind=find(FAR>=0.001);
VER001= TAR(iind(1));

threshold=(smatch-b)/a;

author: happynear

repo: FaceVerification


text:
function A = allcomb(varargin)

% ALLCOMB - All combinations
%    B = ALLCOMB(A1,A2,A3,...,AN) returns all combinations of the elements
%    in the arrays A1, A2, ..., and AN. B is P-by-N matrix is which P is the product
%    of the number of elements of the N inputs. This functionality is also
%    known as the Cartesian Product. The arguments can be numerical and/or
%    characters, or they can be cell arrays.
%
%    Examples:
%       allcomb([1 3 5],[-3 8],[0 1]) % numerical input:
%       % -> [ 1  -3   0
%       %      1  -3   1
%       %      1   8   0
%       %        ...
%       %      5  -3   1
%       %      5   8   1 ] ; % a 12-by-3 array
%
%       allcomb('abc','XY') % character arrays
%       % -> [ aX ; aY ; bX ; bY ; cX ; cY] % a 6-by-2 character array
%
%       allcomb('xy',[65 66]) % a combination
%       % -> ['xA' ; 'xB' ; 'yA' ; 'yB'] % a 4-by-2 character array
%
%       allcomb({'hello','Bye'},{'Joe', 10:12},{99999 []}) % all cell arrays
%       % -> {  'hello'  'Joe'        [99999]
%       %       'hello'  'Joe'             []
%       %       'hello'  [1x3 double] [99999]
%       %       'hello'  [1x3 double]      []
%       %       'Bye'    'Joe'        [99999]
%       %       'Bye'    'Joe'             []
%       %       'Bye'    [1x3 double] [99999]
%       %       'Bye'    [1x3 double]      [] } ; % a 8-by-3 cell array
%
%    ALLCOMB(..., 'matlab') causes the first column to change fastest which
%    is consistent with matlab indexing. Example:
%      allcomb(1:2,3:4,5:6,'matlab')
%      % -> [ 1 3 5 ; 1 4 5 ; 1 3 6 ; ... ; 2 4 6 ]
%
%    If one of the arguments is empty, ALLCOMB returns a 0-by-N empty array.
%
%    See also NCHOOSEK, PERMS, NDGRID
%         and NCHOOSE, COMBN, KTHCOMBN (Matlab Central FEX)

% Tested in Matlab R2015a
% version 4.1 (feb 2016)
% (c) Jos van der Geest
% email: samelinoa@gmail.com

% History
% 1.1 (feb 2006), removed minor bug when entering empty cell arrays;
%     added option to let the first input run fastest (suggestion by JD)
% 1.2 (jan 2010), using ii as an index on the left-hand for the multiple
%     output by NDGRID. Thanks to Jan Simon, for showing this little trick
% 2.0 (dec 2010). Bruno Luong convinced me that an empty input should
% return an empty output.
% 2.1 (feb 2011). A cell as input argument caused the check on the last
%      argument (specifying the order) to crash.
% 2.2 (jan 2012). removed a superfluous line of code (ischar(..))
% 3.0 (may 2012) removed check for doubles so character arrays are accepted
% 4.0 (feb 2014) added support for cell arrays
% 4.1 (feb 2016) fixed error for cell array input with last argument being
%     'matlab'. Thanks to Richard for pointing this out.

narginchk(1,Inf) ;

NC = nargin ;

% check if we should flip the order
if ischar(varargin{end}) && (strcmpi(varargin{end},'matlab') || strcmpi(varargin{end},'john')),
    % based on a suggestion by JD on the FEX
    NC = NC-1 ;
    ii = 1:NC ; % now first argument will change fastest
else
    % default: enter arguments backwards, so last one (AN) is changing fastest
    ii = NC:-1:1 ;
end

args = varargin(1:NC) ;
% check for empty inputs
if any(cellfun('isempty',args)),
    warning('ALLCOMB:EmptyInput','One of more empty inputs result in an empty output.') ;
    A = zeros(0,NC) ;
elseif NC > 1
    isCellInput = cellfun(@iscell,args) ;
    if any(isCellInput)
        if ~all(isCellInput)
            error('ALLCOMB:InvalidCellInput', ...
                'For cell input, all arguments should be cell arrays.') ;
        end
        % for cell input, we use to indices to get all combinations
        ix = cellfun(@(c) 1:numel(c), args,'un',0) ;

        % flip using ii if last column is changing fastest
        [ix{ii}] = ndgrid(ix{ii}) ;

        A = cell(numel(ix{1}),NC) ; % pre-allocate the output
        for k=1:NC,
            % combine
            A(:,k) = reshape(args{k}(ix{k}),[],1) ;
        end
    else
        % non-cell input, assuming all numerical values or strings
        % flip using ii if last column is changing fastest
        [A{ii}] = ndgrid(args{ii}) ;
        % concatenate
        A = reshape(cat(NC+1,A{:}),[],NC) ;
    end
elseif NC==1,
    A = args{1}(:) ; % nothing to combine

else % NC==0, there was only the 'matlab' flag argument
    A = zeros(0,0) ; % nothing
end

author: epfl-lasa


text:
function featureVector = gaborFeatures(img,gaborArray,d1,d2)

% GABORFEATURES extracts the Gabor features of an input image.
% It creates a column vector, consisting of the Gabor features of the input
% image. The feature vectors are normalized to zero mean and unit variance.
%
%
% Inputs:
%       img         :	Matrix of the input image
%       gaborArray	:	Gabor filters bank created by the function gaborFilterBank
%       d1          :	The factor of downsampling along rows.
%       d2          :	The factor of downsampling along columns.
%
% Output:
%       featureVector	:   A column vector with length (m*n*u*v)/(d1*d2).
%                           This vector is the Gabor feature vector of an
%                           m by n image. u is the number of scales and
%                           v is the number of orientations in 'gaborArray'.
%
%
% Sample use:
%
% img = imread('cameraman.tif');
% gaborArray = gaborFilterBank(5,8,39,39);  % Generates the Gabor filter bank
% featureVector = gaborFeatures(img,gaborArray,4,4);   % Extracts Gabor feature vector, 'featureVector', from the image, 'img'.
%
%
%
%   Details can be found in:
%
%   M. Haghighat, S. Zonouz, M. Abdel-Mottaleb, "CloudID: Trustworthy
%   cloud-based and cross-enterprise biometric identification,"
%   Expert Systems with Applications, vol. 42, no. 21, pp. 7905-7916, 2015.
%
%
%
% (C)	Mohammad Haghighat, University of Miami
%       haghighat@ieee.org
%       PLEASE CITE THE ABOVE PAPER IF YOU USE THIS CODE.

if (nargin ~= 4)        % Check correct number of arguments
    error('Please use the correct number of input arguments!')
end

if size(img,3) == 3     % Check if the input image is grayscale
    warning('The input RGB image is converted to grayscale!')
    img = rgb2gray(img);
end

img = double(img);

%% Filter the image using the Gabor filter bank

% Filter input image by each Gabor filter
[u,v] = size(gaborArray);
gaborResult = cell(u,v);
for i = 1:u
    for j = 1:v
        gaborResult{i,j} = imfilter(img, gaborArray{i,j});
    end
end

%% Create feature vector

% Extract feature vector from input image
featureVector = [];
for i = 1:u
    for j = 1:v

        gaborAbs = abs(gaborResult{i,j});
        gaborAbs = downsample(gaborAbs,d1);
        gaborAbs = downsample(gaborAbs.',d2);
        gaborAbs = gaborAbs(:);

        % Normalized to zero mean and unit variance. (if not applicable, please comment this line)
        gaborAbs = (gaborAbs-mean(gaborAbs))/std(gaborAbs,1);

        featureVector =  [featureVector; gaborAbs];

    end
end

%% Show filtered images (Please comment this section if not needed!)

% % Show real parts of Gabor-filtered images
% figure('NumberTitle','Off','Name','Real parts of Gabor filters');
% for i = 1:u
%     for j = 1:v
%         subplot(u,v,(i-1)*v+j)
%         imshow(real(gaborResult{i,j}),[]);
%     end
% end
%
% % Show magnitudes of Gabor-filtered images
% figure('NumberTitle','Off','Name','Magnitudes of Gabor filters');
% for i = 1:u
%     for j = 1:v
%         subplot(u,v,(i-1)*v+j)
%         imshow(abs(gaborResult{i,j}),[]);
%     end
% end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [A,B, S, err] = wls(x,y)
    %WEIGHTED LINEAR LEAST SQUARES REGRESSION
    %	[A, B, S, err] = WLS(x,y) finds the A and B coefficients of A log cumalative frequency
    %curve and the error.
    %      A, B: a- and b- values of a weighted linear regression fit.
    %    err: estimate of the std deviation of the error in predicting a future observation at X by A and B
    %      S: contains fields for triangular factor(R) from QR decomp... see polyfit (used in polyval)
    %
    %report_this_filefun();
    %partially vectorized version

%    global S % output of POLYFIT used for error estimates
    %mima = min(x);

    S=[];
    err=inf;

    if any(size(x) ~= size(y))
        error('X and Y vectors must be the same size.')
    end
    x = x(:);
    y = y(:);
    l = isinf(y);
    y(l) = [];
    x(l) = [];
    % weight the values
    teny= 10.^y;
    wx = ones(1,ceil(sum(teny)));
    wy = wx;
    fteny=floor(teny);
    ks = cumsum([1; fteny]);
    for i = 1:length(x)
        wx(ks(i):ks(i+1)-1) = wx(ks(i):ks(i+1)-1) * x(i);
        wy(ks(i):ks(i+1)-1) = wy(ks(i):ks(i+1)-1) * teny(i);
    end
    %x = wx;
    %y = log10(wy);

    l = wx  > min(x);%mima;

    %[B, A,err ] = ma(x',y');
    %b2 = -abs(B);

    if sum(l) <= 5
        p = [NaN NaN] ;
    elseif nargout > 2
        [p,S] = polyfit(wx(l),log10(wy(l)),1);
        if nargout==4 && S.df > 0
            [~,err] = polyval(p,wx,S);
            err = mean(err);
        end
    else
        p = polyfit(wx(l),log10(wy(l)),1);
    end
    A = p(2);
    B = p(1) ;
end

author: CelsoReyes


text:
function obj = setIndices(obj, index)
     % Sets indinces of the Nlp Variables
    %
    % Parameters:
    %  index: the indices of the variables @type integer

    assert(length(index) == obj.Dimension, ...
        'The length of the variable indices must be equal to the dimension of the NlpVariable object.');

    obj.Indices = index;
end

author: ayonga

repo: frost-dev

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


text:
function res = TTeMPS_op_laplace_to_TTeMPS_op( A )
    %TTeMPS_op_laplace_to_TTeMPS Convert to TT_op_laplace to TTeMPS_op format.
    %
    %   See also TTeMPS_to_TT, TTeMPS_op_to_TT, TTeMPS_op_laplace_to_TTeMPS_op.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    C = cell(1, A.order);
    for i = 1:A.order
        % make a 4D tensor out of it again
        tmp = reshape( full(A.U{i}), [A.rank(i), A.rank(i+1), A.size_col(i), A.size_row(i)] );
        % inverse permute the indices
        C{i} = ipermute( tmp, [1 4 2 3] );
    end
    res = TTeMPS_op(C);
end

author: NicolasBoumal

repo: manopt


text:
function [curImg] = read_current_image(curFolderPath, imgFilesPathList, iCurImg)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    % curFolderPath is the path of current sub folder
    % imgFilesPathList is the path set of all image files in current sub
    % folder
    % iCurImg is the index of current image in the path list

    % get the current image name
    curImgName = imgFilesPathList(iCurImg).name;

    % read current image
    curImg =  imread(strcat(curFolderPath,curImgName));

end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function [ x, obj_val] = QCQP_LB1( H_wave,y_wave,N,l,u)
%UNTITLED9 Summary of this function goes here
%   Detailed explanation goes here
Q = 2*H_wave*H_wave.';
f = -2*H_wave*y_wave;
c = y_wave.'*y_wave;
A = zeros(N,2*N);
for ii = 1:N
    A(ii,ii) = cos((l(ii)+u(ii))/2);
    A(ii,ii+N) = sin((l(ii)+u(ii))/2);         %Linear Constraint
end

function [y,grady] = quadobj(x,Q,f,c)         %Objective Function
    y = 1/2*x.'*Q*x+f.'*x+c;
    if nargout>1
        grady = Q*x + f;
    end
end

function [y_con,yeq,grady,gradyeq] = quadconstr(x)      %Nonlinear Constraint
    for nn = 1:N
        y_con(nn) = x(nn)^2+x(nn+N)^2 - 1;
    end
    yeq = [];
    if nargout>2
       gradyeq = zeros(length(x),N);
       for nn = 1:N
           grady(nn,nn) = 2*x(nn);
           grady(nn+N,nn) = 2*x(nn+N);
       end
    end
    gradyeq = [];
end

Acon = -A;

for ii = 1:N
    b(ii) = -cos((u(ii)-l(ii))/2);
end

% options = optimoptions(@fmincon,'Algorithm','sqp',...
%     'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
%     'Display','off','OptimalityTolerance',1e-6,'MaxIterations',400);

options = optimoptions('fmincon','Algorithm','sqp', 'GradObj','on');

fun = @(x)quadobj(x,Q,f,c);
nonlconstr = @(x)quadconstr(x);
x1 = ones(2*N,1);
x = fmincon(fun,x1,Acon,b,[],[],[],[],nonlconstr,options);
obj_val = objval_func(x,H_wave,y_wave);
end

author: yuanhao-cui

repo: Must-Reading-on-ISAC

path: github-repos/MATLAB/yuanhao-cui-Must-Reading-on-ISAC/Must-Reading-on-ISAC-34cd6615c52ebca121428a979e756c608b195040/Codes/Fan2018TSP/Codes for DFRC Waveform Design/Constant Modulus/QCQP_LB1.m


text:
%
% rho=findrhoc(n,p)
%
% Finds an intensity rho such that C(n,rho)=p.
%
% Note: Must have 0<p<1.  Returns NaN if p is not in this range.
%
function rho=findrhoc(n,p)
%
% Sanity check- make sure that n is a positive integer.
%
  if ((floor(n) ~= n) || (n < 1))
    warning('n is not a positive integer');
    rho=NaN;
    return;
  end;
%
% Sanity check- make sure that p is a probability with 0<p<1.
%
  if ((p<0.0) || (p>1.0))
     warning('Invalid p value!');
     rho=NaN;
     return;
  end;
%
% We know that at rho=0, p=0, and at rho=+Inf, p=1.   We start by finding
% an interval [0,a] containing the root.
%
a=1.0;
testp=erlangc(n,a);
while (testp < p),
  a=a*2.0;
  testp=erlangc(n,a);
end;
%
% Now, the root is somewhere between 0 and a.  Use bisection to find it.
%
%
left=0.0;
right=a;
mid=(left+right)/2;
midp=erlangc(n,mid);
while ((right-left) > 0.0001*max([1 left])),
  if (midp < p),
    left=mid;
    mid=(left+right)/2;
    midp=erlangc(n,mid);
  else
    right=mid;
    mid=(left+right)/2;
    midp=erlangc(n,mid);
  end;
end;
%
% Return the left end point of the current interval, which has prob < p.
%
rho=left;

author: Sable

repo: mcbench-benchmarks


text:
function [ o, x, w ] = en_r2_09_1 ( n, option )

%*****************************************************************************80
%
%% EN_R2_09_1 implements the Stroud rule 9.1 for region EN_R2.
%
%  Discussion:
%
%    The rule has order O = ( 2 * N^4 - 4 * N^3 + 22 * N^2 - 8 * N + 3 ) / 3.
%
%    The rule has precision P = 9.
%
%    EN_R2 is the entire N-dimensional space with weight function
%
%      w(x) = exp ( - x1^2 - x2^2 ... - xn^2 )
%
%    There are two versions of each rule, chosen by setting the
%    OPTION variable to 1 or 2.
%
%    The rule as tabulated by Stenger is available for N = 2 through 20.
%    This function accepts N = 3 through 6.
%
%     N    O
%    __  ___
%     3   77
%     4  193
%     5  421
%     6  825
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%    3 <= N <= 6.
%
%    Input, integer OPTION, chooses rule option 1 or 2.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( n < 3 | 6 < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EN_R2_09_1 - Fatal error!\n' );
    fprintf ( 1, '  3 <= N <= 6 required.\n' );
    error ( 'EN_R2_09_1 - Fatal error!' )
  end

  if ( nargin < 2 )
    option = 1;
  end

  if ( option < 1 | 2 < option )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'EN_R2_09_1 - Fatal error!\n' );
    fprintf ( 1, '  1 <= OPTION <= 2 required.\n' );
    error ( 'EN_R2_09_1 - Fatal error!' )
  end

  o = ( 2 * n^4 - 4 * n^3 + 22 * n^2 - 8 * n + 3 ) / 3;
  volume = sqrt ( pi^n );

  if ( n == 3 )
    u =    0.202018287045609E+01;
    v =    0.958572464613819E+00;
    b0 =   0.676448734429924E+00;
    b1 =   0.511989106291551E-02;
    b2 =   0.448595723493744E+00;
    b3 =   0.235223454595606E-03;
    b4 =   0.915390713080005E-01;
    b5 =   0.139208199920793E-01;
    b6 =   0.235223454595606E-03;
    b7 =   0.915390713080008E-01;
    b8 =   0.000000000000000E+00;
  elseif ( n == 4 & option == 1 )
    u =    0.202018287045609E+01;
    v =    0.958572464613819E+00;
    b0 = - 0.860452945007048E+00;
    b1 = - 0.405511998533795E-01;
    b2 =   0.107026475449715E+01;
    b3 =   0.138974239307092E-03;
    b4 = - 0.162248779448181E+00;
    b5 =   0.246740110027234E-01;
    b6 =   0.138974239307094E-03;
    b7 =   0.162248779448181E+00;
    b8 =   0.138974239307094E-03;
  elseif ( n == 4 & option == 2 )
    u =    0.958572464613819E+00;
    v =    0.202018287045609E+01;
    b0 =   0.265029088766810E-02;
    b1 =   0.637601342635332E+00;
    b2 = - 0.394394059389228E-01;
    b3 =   0.540829264827264E-01;
    b4 = - 0.416922717921281E-03;
    b5 =   0.246740110027234E-01;
    b6 =   0.540829264827270E-01;
    b7 =   0.416922717921281E-03;
    b8 =   0.540829264827269E-01;
  elseif ( n == 5 & option == 1 )
    u =    0.202018287045609E+01;
    v =    0.958572464613819E+00;
    b0 = - 0.827347006200826E+01;
    b1 = - 0.160820174530905E+00;
    b2 =   0.353499863758467E+01;
    b3 =   0.738976276909564E-03;
    b4 = - 0.862735421812943E+00;
    b5 =   0.437335458190621E-01;
    b6 = - 0.246325425636523E-03;
    b7 =   0.287578473937648E+00;
    b8 =   0.246325425636523E-03;
  elseif ( n == 5 & option == 2 )
    u =    0.958572464613819E+00;
    v =    0.202018287045609E+01;
    b0 = - 0.624416791055272E+00;
    b1 =   0.467494915583104E+00;
    b2 = - 0.152937760910536E+00;
    b3 =   0.287578473937646E+00;
    b4 = - 0.221692883072871E-02;
    b5 =   0.437335458190621E-01;
    b6 = - 0.958594913125490E-01;
    b7 =   0.738976276909568E-03;
    b8 =   0.958594913125492E-01;
  elseif ( n == 6 & option == 1 )
    u =    0.202018287045609E+01;
    v =    0.958572464613819E+00;
    b0 = - 0.361840434143098E+02;
    b1 = - 0.447936529138517E+00;
    b2 =   0.112077863004144E+02;
    b3 =   0.392940404320855E-02;
    b4 = - 0.254859786784158E+01;
    b5 =   0.775156917007496E-01;
    b6 = - 0.130980134773619E-02;
    b7 =   0.509719573568315E+00;
    b8 =   0.436600449245395E-03;
  elseif ( n == 6 & option == 2 )
    u =    0.958572464613819E+00;
    v =    0.202018287045609E+01;
    b0 =   0.448873836333650E+01;
    b1 = - 0.238473566140736E+01;
    b2 = - 0.413008493198885E+00;
    b3 =   0.152915872070494E+01;
    b4 = - 0.654900673868093E-02;
    b5 =   0.775156917007496E-01;
    b6 = - 0.509719573568314E+00;
    b7 =   0.130980134773618E-02;
    b8 =   0.169906524522772E+00;
  end

  x = zeros(n,o);
  w = zeros(o,1);

  k = 0;
%
%  1 point.
%
  k = k + 1;
% x(1:n,k) = 0.0;
  w(k) = b0;
%
%  2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - u;
    w(k) = b1;
    k = k + 1;
    x(i,k) = + u;
    w(k) = b1;
  end
%
%  2 * N points.
%
  for i = 1 : n
    k = k + 1;
    x(i,k) = - v;
    w(k) = b2;
    k = k + 1;
    x(i,k) = + v;
    w(k) = b2;
  end
%
%  4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = - u;
      w(k) = b3;
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = + u;
      w(k) = b3;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = - u;
      w(k) = b3;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = + u;
      w(k) = b3;
    end
  end
%
%  4 * ( N * ( N - 1 ) / 2 ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = - v;
      w(k) = b4;
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = + v;
      w(k) = b4;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = - v;
      w(k) = b4;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = + v;
      w(k) = b4;
    end
  end
%
%  4 * ( N * ( N - 1 ) ) points.
%
  for i = 1 : n - 1
    for j = i + 1 : n
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = - v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = - u;
      x(j,k) = + v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = - v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = + u;
      x(j,k) = + v;
      w(k) = b5;
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = - u;
      w(k) = b5;
      k = k + 1;
      x(i,k) = - v;
      x(j,k) = + u;
      w(k) = b5;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = - u;
      w(k) = b5;
      k = k + 1;
      x(i,k) = + v;
      x(j,k) = + u;
      w(k) = b5;
    end
  end
%
%  8 * ( N * ( N - 1 ) * ( N - 2 ) / 6 ) points.
%
  for i = 1 : n - 2
    for j = i + 1 : n - 1
      for l = j + 1 : n
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - u;
        x(l,k) = - u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = - u;
        x(l,k) = + u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + u;
        x(l,k) = - u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = - u;
        x(j,k) = + u;
        x(l,k) = + u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - u;
        x(l,k) = - u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = - u;
        x(l,k) = + u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + u;
        x(l,k) = - u;
        w(k) = b6;
        k = k + 1;
        x(i,k) = + u;
        x(j,k) = + u;
        x(l,k) = + u;
        w(k) = b6;
      end
    end
  end
%
%  8 * ( N * ( N - 1 ) * ( N - 2 ) / 6 ) points.
%
  for i = 1 : n - 2
    for j = i + 1 : n - 1
      for l = j + 1 : n
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = - v;
        x(l,k) = - v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = - v;
        x(l,k) = + v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = + v;
        x(l,k) = - v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = - v;
        x(j,k) = + v;
        x(l,k) = + v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = - v;
        x(l,k) = - v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = - v;
        x(l,k) = + v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = + v;
        x(l,k) = - v;
        w(k) = b7;
        k = k + 1;
        x(i,k) = + v;
        x(j,k) = + v;
        x(l,k) = + v;
        w(k) = b7;
      end
    end
  end
%
%  16 * ( N * ( N - 1 ) * ( N - 2 ) * ( N - 3 ) / 24 ) points.
%
  for i = 1 : n - 3
    for j = i + 1 : n - 2
      for l = j + 1 : n - 1
        for m = l + 1 : n
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = - u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = - u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = - u;
          x(m,k) = + u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = - u;
          w(k) = b8;
          k = k + 1;
          x(i,k) = + u;
          x(j,k) = + u;
          x(l,k) = + u;
          x(m,k) = + u;
          w(k) = b8;
        end
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Globally Optimal Solutions,
%           Solvability/Observability Analysis and Uncertainty Description.
%           IEEE Transactions on Robotics.
%           https://doi.org/10.1109/TRO.2022.3155880
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com
%
%
% test_rel_att.m: The QPEP illustration of range-based
%                 relative attitude estimtion

clear all
close all
clc

if(verLessThan('matlab', '8.0.0'))
   error('The MATLAB version is too old to be supported.');
end

addpath('func_files');
addpath('solvers');
addpath('utils');
addpath('homotopy');

p00 = randn(3, 1);
d00 = sqrt(p00.' * p00);
C0 = orthonormalize(randn(3, 3));

len = 10;
p1 = zeros(3, len);
p2 = zeros(3, len);
epsilon = zeros(len, 1);
noise = 1e-2; % Noise level
for i = 1 : len
    p1(:, i) = randn(3, 1);
    p2(:, i) = randn(3, 1);
    epsilon(i) = (p1(:, i) - p00).' * C0 * p2(:, i) + p1(:, i).' * p00 + noise * randn(1, 1);
end

syms q0 q1 q2 q3
q = [q0; q1; q2; q3];
syms t1 t2 t3;
t = [t1; t2; t3];
R = q2R(q);
rr = R.' * t;
syms r1 r2 r3
r = [r1; r2; r3];
eqs = sym(zeros(len + 2, 1));
for i = 1 : len
    eqs(i) = p1(:, i).' * R * p2(:, i) + p1(:, i).' * t - r.' * p2(:, i) - epsilon(i);
end
eqs(len + 1) = q.' * q - 1;
eqs(len + 2) = r.' * r - t.' * t;
eqs = expand(eqs);
x = [q; t; r];
H = expand(jacobian(eqs, x).' * eqs);
assumeAlso(q.' * q == 1);
assumeAlso(r.' * r == t.' * t);
eq = vpa(expand(simplify(H)), 32);
ss = vpasolve(eq(5 : 10), [t; r]);
eq_ = eq(1 : 4);
eq_ = subs(eq_, t1, ss.t1);
eq_ = subs(eq_, t2, ss.t2);
eq_ = subs(eq_, t3, ss.t3);
t_func = matlabFunction([ss.t1; ss.t2; ss.t3], 'Vars', {q});
eq_ = subs(eq_, r1, ss.r1);
eq_ = subs(eq_, r2, ss.r2);
eq_ = subs(eq_, r3, ss.r3);
eq_ = vpa(expand(eval(eq_)), 32);
syms lambda
eqs = [
    eq_;
    q.' * q - 1;
    ]

str = '';
for i = 1 : length(eqs)
    str = strcat(str, sprintf(' PP{%d} = char(vpa(%%s, 32));', i));
end

str_ = sprintf(str, char(eqs(1)), ...
                    char(eqs(2)), ...
                    char(eqs(3)), ...
                    char(eqs(4)), ...
                    char(eqs(5)));
eval(str_);
[S, vars] = psolve(PP);
S = S.';
SS = S;
for i = 1 : length(vars)
    if(strcmp(vars{i}, 'q0'))
        SS(:, 1) = S(:, i);
    elseif(strcmp(vars{i}, 'q1'))
        SS(:, 2) = S(:, i);
    elseif(strcmp(vars{i}, 'q2'))
        SS(:, 3) = S(:, i);
    elseif(strcmp(vars{i}, 'q3'))
        SS(:, 4) = S(:, i);
    elseif(strcmp(vars{i}, 'lambda'))
        SS(:, 5) = S(:, i);
    end
end
S = real(SS);
xs_ = S;
sols = SS.';

num = size(sols, 2);
sol = zeros(4, num);
ts = zeros(3, num);
Ls = zeros(len, 1);
for i = 1 : num
    sol(:, i) = real(sols(1 : 4, i));
    sol(:, i) = sol(:, i) ./ norm(sol(:, i));
    C = q2R(sol(:, i));
    t = t_func(sol(:, i));
    ts(:, i) = t;
    loss = 0;
    for j = 1 : len
        loss = loss + (epsilon(j) - (p1(:, j) - t).' * C * p2(:, j) - p1(:, j).' * t)^2;
    end
    Ls(i) = loss;
end
[~, idx] = sort(Ls);

q_ = positive_quat(sol(:, idx(1)).')
q_true = positive_quat(dcm2quat(C0))

t_ = ts(:, idx(1)).'
t_true = p00.'

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function G = slnngraph(X, X2, nnparams, varargin)
%SLNNGRAPH Constructs a nearest neighborhood based graph
%
% $ Syntax $
%   - G = slnngraph(X, [], nnparams, ...)
%   - G = slnngraph(X, X2, nnparams, ...)
%
% $ Arguments $
%   - X:        The sample matrix of the (referenced) nodes
%   - X2:       The sample matrix of the (query) nodes
%   - nnparams: The cell array of parameters to slfindnn for determining
%               neighborhoods.
%   - G:        The adjacency matrix of the constructed graph
%
% $ Description $
%   - G = slnngraph(X, [], nnparams, ...) constructs a graph with adjacency
%     matrix representation using specified nearest neighbor strategy.
%     You can further specify the following properties to control the
%     process of adjacency matrix construction.
%       \*
%       \t   The Properties of Graph Matrix construction           \\
%       \h     name         &      description                      \\
%             'valtype'     & The type of the matrix values          \\
%                             - 'logical':  using logical value
%                             - 'numeric':  using numeric value (default)
%             'sparse'      & Whether to construct sparse matrix
%                             (default = true)                      \\
%             'tfunctor'    & The function to transform the distance
%                             values to edge values. (default = [])  \\
%             'sym'         & whether to symmetrizes the graph
%                             (default = false)                       \\
%             'symmethod'   & The method used to symmetrization
%                             (default = [])                          \\
%       \*
%
%   - G = slnngraph(X, X2, nnparams, ...) constructs a bigraph with the
%     source and target set respectively specified.
%
%   - There are three configurations on X and X2:
%       - construct graph with the same set X with all edges connecting
%         the same node excluded. You can use the following syntax:
%           slnngraph(X, [], ...)
%       - construct graph with the same set X with the edges connecting
%         the same node preserved. You can use the following syntax:
%           slnngraph(X, X, ...)
%       - construct graph with the different source and target sets X and
%         X2, You can use the following syntax:
%           slnngraph(X, X2, ...)
%
% $ Remarks $
%   - tfunctor only takes effect for numeric value type.
%
%   - The option sym and symmethod only take effect when X and X2 has the
%     same number of elements.
%
% $ History $
%   - Created by Dahua Lin, on Sep 8th, 2006
%   - Modified by Dahua Lin, on Sep 11st, 2006
%       - revise to conform to the defined neighborhood system
%         representation
%       - fix small bugs
%

%% parse input

if nargin < 3
    raise_lackinput('slnngraph', 3);
end

if ~isnumeric(X) || ndims(X) ~= 2
    error('sltoolbox:invalidarg', ...
        'X should be a 2D numeric matrix');
end

if isempty(X2)
    X2 = [];
else
    if ~isnumeric(X2) || ndims(X2) ~= 2
        error('sltoolbox:invalidarg', ...
            'A non-empty X2 should be a 2D numeric matrix');
    end
end

if ~iscell(nnparams)
    error('stoolbox:invalidarg', ...
        'The parameters for slfindnn should be groupped in a cell array');
end

opts.valtype = 'numeric';
opts.sparse = true;
opts.tfunctor = [];
opts.sym = false;
opts.symmethod = [];
opts = slparseprops(opts, varargin{:});

switch opts.valtype
    case 'logical'
        islogic = true;
    case 'numeric';
        islogic = false;
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid value type for graph: %s', opts.valtype);
end

%% Main

n0 = size(X, 2);
if isempty(X2)
    nq = n0;
else
    nq = size(X2, 2);
end
can_sym = (n0 == nq);

% find nearest neighbor
if ~islogic
    [nnidx, dists] = slfindnn(X, X2, nnparams{:});
else
    nnidx = slfindnn(X, X2, nnparams{:});
end

% group edges and vectorizes distances
edges = sladjlist2edgeset(nnidx, 0);
edges = edges(:, [2,1]);  % flip in order to place source to left, target to right
clear nnidx;
if ~islogic
    dists = vertcat(dists{:});
else
    dists = [];
end

% compute edge values
if ~isempty(dists)
    if isempty(opts.tfunctor)
        vals = dists;
    else
        vals = slevalfunctor(opts.tfunctor, dists);
    end
else
    vals = [];
end
clear dists;

% make graph matrix
G = slmakeadjmat(n0, nq, edges, vals, islogic, opts.sparse);

% symmetrize the graph
if opts.sym && can_sym
    G = slsymgraph(G, opts.symmethod);
end

author: lmthang


text:
function f = previewfig(h,varargin)
    %PREVIEWFIG  Preview a figure to be exported using EXPORTFIG.
    %   F = PREVIEWFIG(H) creates a preview of H with the default
    %   EXPORTFIG options and returns the preview's figure handle in F.
    %   F = PREVIEWFIG(H,OPTIONS) creates a preview with OPTIONS as
    %   described in EXPORTFIG.
    %   PREVIEWFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) creates a preview
    %   with the specified parameter-value pairs to H as described in
    %   EXPORTFIG.
    %
    %   See also EXPORTFIG, APPLYTOFIG, RESTOREFIG.

    %  Copyright 2000 Ben Hinkle
    %  Email bug reports and comments to bhinkle@mathworks.com

    filename = [tempname, '.png'];
    args = {'resolution',0,'format','png'};
    if nargin > 1
        exportfig(h, filename, varargin{:}, args{:});
    else
        exportfig(h, filename, args{:});
    end

    X = imread(filename,'png');
    height = size(X,1);
    width = size(X,2);
    delete(filename);
    f = figure_w_normalized_uicontrolunits( 'Name', 'Preview', ...
        'Menubar', 'none', ...
        'NumberTitle', 'off', ...
        'Visible', 'off');
    image(X);
    axis image;
    ax = findobj(f, 'type', 'axes');
    axesPos = [0 0 width height];
    set(ax, 'Units', 'pixels', ...
        'Position', axesPos, ...
        'Visible', 'off');
    figPos = get(f,'Position');
    rootSize = get(groot,'ScreenSize');
    figPos(3:4) = axesPos(3:4);
    if figPos(1) + figPos(3) > rootSize(3)
        figPos(1) = rootSize(3) - figPos(3) - 50;
    end
    if figPos(2) + figPos(4) > rootSize(4)
        figPos(2) = rootSize(4) - figPos(4) - 50;
    end
    set(f, 'Position',figPos, ...
        'Visible', 'on');


text:
%% THE FACTORIZE OBJECT for solving linear systems
%
% Copyright 2009, Timothy A. Davis, University of Florida.  May 19, 2009.
% davis@cise.ufl.edu
% http://www.cise.ufl.edu/~davis
%
% This is a demonstration of the FACTORIZE object for solving linear
% systems and least-squares problems, and for computations with the
% matrix inverse and pseudo-inverse.

%% Rule Number One: never multiply by the inverse, inv(A)
%
% Use backslash or a matrix factorization instead (LU, CHOL, or QR).

%% Rule Number Two:  never break Rule Number One
%
% However, the problem with Rule Number One is that it can be hard to
% figure out which matrix factorization to use and how to use it.
% BACKSLASH (MLDIVIDE) is great, but it can't be reused when solving
% multiple systems (x=A\b and y=A\c).  Its syntax doesn't match the use of
% the inverse in mathematical expressions, either.
%
% The goal of the FACTORIZE object is to solve this problem ...
%
% "Don't let that INV go past your eyes; to solve that system, FACTORIZE!"

%% How to use BACKSLASH solve A*x=b
%
% First, let's create a square matrix A and a right-hand-side b for a
% linear system A*x=b.  There are many ways to solve this system.  The best
% way is to use x=A\b.  The residual r is a vector of what's left over in
% each equation, and its norm tells you how accurately the system was
% solved.

format compact ;
A = rand (3)
b = rand (3,1)
x = A\b
r = b-A*x ;
norm (r)

%% BACKSLASH versus INV ... let the battle begin
%
% The backslash operation x=A\b is mathematically the same as x=inv(A)*b.
% However, backslash is faster and more accurate since it uses a matrix
% factorization instead of multiplying by the inverse.  Even though your
% linear algebra textbook might write x=A^(-1)*b as the solution to the
% system A*x=b, your textbook author never means for you to compute the
% inverse.
%
% These next statements give the same answer, so what's the big deal?

S = inv(A) ;
x = S*b
x = A\b

%%
% The big deal is that you should care about speed and you should care even
% more about accuracy.  BACKSLASH relies on matrix factorization (LU, CHOL,
% QR, or other specialized methods).  It's faster and more reliable than
% multiplying by the inverse, particularly for large matrices and sparse
% matrices.  Here's an illustration of how pathetic inv(A)*b can be.

A = gallery ('frank',16) ; xtrue = ones (16,1) ; b = A*xtrue ;

x = inv(A)*b ; norm (A*x-b)
x = A\b      ; norm (A*x-b)

%%
% The performance difference between BACKSLASH and INV for even small
% sparse matrices is striking.

load west0479 ;
A = west0479 ;
n = size (A,1)
b = rand (n,1) ;
tic ; x = A\b ; toc
norm (b-A*x)
tic ; x = inv(A)*b ; toc
norm (b-A*x)

%%
% What if you want to solve multiple systems?  Use a matrix factorization.
% But which one?  And how do you use it?  Here are some alternatives using
% LU for the sparse west0479 matrix, but some are faster than others.

tic ; [L,U]     = lu(A) ; x1 = U \ (L \ b)         ; t1=toc ; nz1=nnz(L+U);
tic ; [L,U,P]   = lu(A) ; x2 = U \ (L \ P*b)       ; t2=toc ; nz2=nnz(L+U);
tic ; [L,U,P,Q] = lu(A) ; x3 = Q * (U \ (L \ P*b)) ; t3=toc ; nz3=nnz(L+U);

fprintf ('1: nnz(L+U): %5d time: %8.4f resid: %e\n', nz1,t1, norm(b-A*x1));
fprintf ('2: nnz(L+U): %5d time: %8.4f resid: %e\n', nz2,t2, norm(b-A*x2));
fprintf ('3: nnz(L+U): %5d time: %8.4f resid: %e\n', nz3,t3, norm(b-A*x3));

%% LU and LINSOLVE are fast and accurate but complicated to use
%
% A quick look at ``help lu'' will scroll off your screen.  For full
% matrices, [L,U,p] = lu (A,'vector') is fastest.  Then for the
% forward/backsolves, use LINSOLVE instead of BACKSLASH for even faster
% performance.  But for sparse matrices, use the optional 'Q' output of LU
% so you get a good fill-reducing ordering.  But you can't use 'Q' if the
% matrix is full.  But LINSOLVE doesn't work on sparse matrices.
%
% But ... Ack!  That's getting complicated ...
%
% Here's the best way to solve A*x=b and A*y=c when A is full and
% unsymmetric:

n = 1000 ;
A = rand (n) ;
b = rand (n,1) ;
c = rand (n,1) ;
tic ; [L,U,p] = lu (A, 'vector') ; LUtime = toc

tic ; x = U \ (L \ b (p,:)) ;
      y = U \ (L \ c (p,:)) ; toc

tic ; opL = struct ('LT', true) ;
      opU = struct ('UT', true) ;
      x = linsolve (U, linsolve (L, b(p,:), opL), opU) ;
      y = linsolve (U, linsolve (L, c(p,:), opL), opU) ; toc

%% INV is easy to use, but slow and inaccurate
%
% Oh bother!  Using LU and LINSOLVE is too complicated.  You just want to
% solve your system.  Let's just compute inv(A) and use it twice.  Easy to
% write, but slower and less accurate ...

S = inv (A) ;
x = S*b ; norm (b-A*x)
y = S*c ; norm (c-A*y)

%%
% Sometimes using the inverse seems inevitable.  For example, your textbook
% might show the Schur complement formula as S = A-B*inv(D)*C.  This can be
% done without inv(D) in one of two ways: SLASH or BACKSLASH (MRDIVIDE or
% MLDIVIDE to be precise).
%
% inv(A)*B and A\B are mathematically equivalent, as are B*inv(A) and B/A,
% so these three methods give the same results (ignoring computational
% errors, which are worse for inv(D)).  Only the first equation looks like
% the equation in your textbook, however.

A = rand (200) ; B = rand (200) ; C = rand (200) ; D = rand (200) ;

tic ; S1 = A - B*inv(D)*C ; toc ;
tic ; S2 = A - B*(D\C) ;    toc ;
tic ; S3 = A - (B/D)*C ;    toc ;

%% So the winner is ... nobody
%
% BACKSLASH: mostly simple to use (except remember that Schur complement
%       formuala?).  Fast and accurate ... but slow if you want to solve
%       two linear systems with the same matrix A.
%
% LU, QR, CHOL: fast and accurate.  Awful syntax to use.  Drag out your
%       linear algebra textbook if you want to use these in MATLAB.
%       Whenever I use them I have to derive them from scratch, even
%       though I *wrote* most of the sparse factorizations used in MATLAB!
%
% INV: slow and inaccurate.  Wins big on ease-of-use, though, since it's a
%       direct plug-in for all your nice mathematical formulas.
%
% No method is best on all three criterion: speed, accuracy, and ease of
% use.
%
% Is there a solution?  Yes ... keeping reading ...

%% The FACTORIZE object to the rescue
%
% The FACTORIZE method is just as easy to use as INV, but just as fast and
% accurate as BACKSLASH, LU, QR, CHOL, and LINSOLVE.
%
% F = factorize(A) computes the factorization of A and returns it as an
% object that you can reuse to solve a linear system with x=F\b.  It picks
% LU, QR, or Cholesky for you, just like BACKSLASH.
%
% S = inverse(A) is simpler yet.  It does NOT compute inv(A), but
% factorizes A.  When multiplying S*b, it doesn't mulitply by the inverse,
% but uses the correct forward/backsolve equations to solve the linear
% system.

n = 1000 ;
A = rand (n) ;
b = rand (n,1) ;
c = rand (n,1) ;

tic ;                       x = A\b ; y = A\c ; toc
tic ; S = inv(A) ;          x = S*b ; y = S*c ; toc
tic ; F = factorize(A) ;    x = F\b ; y = F\c ; toc
tic ; S = inverse(A) ;      x = S*b ; y = S*c ; toc

%% Least-squares problems
%
% Here are some different methods for solving a least-squares problem when
% your system is over-determined.  The last two methods are the same.

A = rand (1000,200) ;
b = rand (1000,1) ;

tic ; x = A\b            ; toc, norm (A'*A*x-A'*b)
tic ; x = pinv(A)*b      ; toc, norm (A'*A*x-A'*b)
tic ; x = inverse(A)*b   ; toc, norm (A'*A*x-A'*b)
tic ; x = factorize(A)\b ; toc, norm (A'*A*x-A'*b)

%%
% FACTORIZE is better than BACKSLASH because you can reuse the
% factorization for different right-hand-sides.  For full-rank matrices,
% it's better than PINV because it's faster (and PINV fails for sparse
% matrices).

A = rand (1000,200) ;
b = rand (1000,1) ;
c = rand (1000,1) ;

tic ;                  ; x = A\b ; y = A\c ; toc
tic ; S = pinv(A)      ; x = S*b ; y = S*c ; toc
tic ; S = inverse(A)   ; x = S*b ; y = S*c ; toc
tic ; F = factorize(A) ; x = F\b ; y = F\c ; toc

%% Underdetermined systems
%
% The under-determined system A*x=b where A has more columns than rows has
% many solutions.  x=A\b finds a basic solution (some of the entries in x
% are zero).  pinv(A)*b finds a minimum 2-norm solution, but it's slow.  QR
% factorization will do the same if A has full rank.  That's what the
% factorize(A) and inverse(A) methods do.

A = rand (200,1000) ;
b = rand (200,1) ;

tic ; x = A\b            ; toc, norm (x)
tic ; x = pinv(A)*b      ; toc, norm (x)
tic ; x = inverse(A)*b   ; toc, norm (x)
tic ; x = factorize(A)\b ; toc, norm (x)

%% Computing selected entries in the inverse or pseudo-invers
%
% If you want just a few entries from the inverse, it's still better to
% formulate the problem as a system of linear equations and to use a
% matrix factorization instead of computing inv(A).  The FACTORIZE object
% does this for you, by overloading the subsref operator.

A = rand (1000) ;

tic ; S = inv (A)     ; S (2:3,4), toc
tic ; S = inverse (A) ; S (2:3,4), toc

%% Computing the entire inverse or pseudo-inverse
%
% Rarely, and I mean RARELY, you really do need the inverse.  More
% frequently what you want is the pseudo-inverse.  You can force a
% factorization to become a plain matrix by converting it to double.  Note
% that inverse(A) only handles full-rank matrices (either dense or sparse),
% whereas pinv(A) works for all dense matrices (not sparse).

A = rand (500) ;
tic ; S1 = inv (A) ;            ; toc
tic ; S2 = double (inverse (A)) ; toc
norm (S1-S2)

A = rand (500,400) ;
tic ; S1 = pinv (A)             ; toc
tic ; S2 = double (inverse (A)) ; toc
norm (S1-S2)

%% Update/downdate of a dense Cholesky factorization
%
% Wilkinson considered the update/downdate of a matrix factorization to be
% a key problem in computational linear algebra.  The idea is that you
% first factorize a matrix.  Next, make a low-rank change to A, and patch
% up (or down...) the factorization so that it becomes the factorization of
% the new matrix.  In MATLAB, this only works for dense symmetric positive
% definite matrices, via cholupdate.  This is much faster than computing
% the new factorization from scratch.

n = 1000 ;
A = rand (n) ;
A = A*A' + n*eye (n) ;
w = rand (n,1) ; t = rand (n,1) ; b = rand (n,1) ;
F = factorize (A) ;

tic ; F = F + w ; x = F\b ; toc
tic ; y = (A+w*w')\b ;      toc
norm (x-y)

tic ; F = F - t ; x = F\b ; toc
tic ; y = (A+w*w'-t*t')\b ; toc
norm (x-y)

%% Caveat Executor
%
% One caveat:  If you have a large number of very small systems to solve,
% the object-oriented overhead of creating and using an object can dominate
% the run time, at least in MATLAB R2009a.  For this case, if you want the
% best performance, stick with BACKSLASH, or LU and LINSOLVE (just extract
% the appropriate formulas from factorize.m and mldivide.m in the
% @factorize directory).
%
% Hopefully the object-oriented overhead will drop in future versions of
% MATLAB, and you can ignore this caveat.

A = rand (10) ; b = rand (10,1) ; F = factorize (A) ;

tic ; for k = 1:10000, x = F\b ; end ; toc

tic ; for k = 1:10000, x = A\b ; end ; toc

[L,U,p] = lu (A, 'vector') ;
opL = struct ('LT', true) ;
opU = struct ('UT', true) ;
tic ;
for k = 1:10000
    x = linsolve (U, linsolve (L, b(p,:), opL), opU) ;
end
toc

%% Summary
%
% So ... don't use INV, and don't worry about how to use LU, CHOL, or QR
% factorization.  Just install the FACTORIZE package, and you're on your
% way.  Assuming you are now in the Factorize/ directory, cut-and-paste
% these commands into your command window:
%
%   addpath (pwd)
%   savepath
%
% And remember ...
%
% "Don't let that INV go past your eyes; to solve that system, FACTORIZE!"

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% problem 7 - function that computes the convolution of x[n] and h[n]

function [y,n] = convd(x,n1,h,n2)
a = n1(1)+n2(1);
b = n1(end) + n2(end);
n =a:b;
y = conv(x,h);
stem(n,y);
legend('y[n]')


text:
function y = trace_inv( varargin )

% TRACE_INV   Trace of the inverse of a PSD matrix.
%     For square matrix X, TRACE_INV(X) is TRACE(INV(X)) if X is Hermitian
%     or symmetric and positive definite; and +Inf otherwise.
%
%     An error results if X is not a square matrix.
%
%     Disciplined convex programming information:
%         TRACE_INV is convex and nonmonotonic (at least with respect to
%         elementwise comparison), so its argument must be affine.

persistent P
if isempty( P ),
    P.nargs     = 1;
    P.args      = [];
    P.empty     = 0;
    P.constant  = @trace_inv_diag;
    P.diagonal  = @trace_inv_diag;
    P.affine    = @trace_inv_aff;
    P.structure = 'psdeig';
end
y = cvx_matrix_op( P, varargin );

function z = trace_inv_diag( D )
z = sum( 1.0 ./ D );

function z = trace_inv_aff( X ) %#ok
cvx_begin sdp
    epigraph variable z nonnegative_
    variable Y(size(X)) hermitian_if(X)
    real(trace(Y)) <= z;
    [Y,eye(size(X));eye(size(X)),X] >= 0; %#ok
cvx_end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% This subroutine assigns creates a grid with
% spacing dx,dy (in degreees). The size will
% be selected interactiVELY. The bvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/95

display ('This is /src/bdepth_ratio_auto.m');

global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dx = 0.1;
    dy = 0.1 ;
    ni = 1000;
    Nmin = 50;
    stan2 = nan;
    stan = nan;
    prf = nan;
    av = nan;
    mid_point = 5;
    top_zonet = 0;
    top_zoneb = 5;
    bot_zonet = 7;
    bot_zoneb = 15;
    topstan2 = nan;
    botstan2 = nan;
    topstan = nan;
    botstan = nan;
    topav = nan;
    botav = nan;

    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Depth Ratio Grid Input Parameter',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ wex+200 wey-200 650 300]);
    axis off
    labelList2=[' Automatic Mcomp (max curvature) | Fixed Mc (Mc = Mmin) | Automatic Mcomp (90% probability) | Automatic Mcomp (95% probability) | Best (?) combination (Mc95 - Mc90 - max curvature)'];

    labelPos = [0.2 0.60 0.6  0.08];
    hndl2=uicontrol(...
        'Style','popup',...
        'Position',labelPos,...
        'Units','normalized',...
        'String',labelList2,...
        'Callback','inb2 =get(hndl2,''Value''); ');

    % set(hndl2,'value',5);

    % creates a dialog box to input grid parameters
    %
    %

    % mid_point_field=uicontrol('Style','edit',...
    %     'Position',[.47 .80 .12 .08],...
    %     'Units','normalized','String',num2str(mid_point),...
    %     'Callback','mid_point=str2double(get(mid_point_field,''String'')); set(mid_point_field,''String'',num2str(mid_point));');

    top_zonet_field=uicontrol('Style','edit',...
        'Position',[.36 .80 .06 .06],...
        'Units','normalized','String',num2str(top_zonet),...
        'Callback','top_zonet=str2double(get(top_zonet_field,''String'')); set(top_zonet_field,''String'',num2str(top_zonet));');
    top_zoneb_field=uicontrol('Style','edit',...
        'Position',[.36 .74 .06 .06],...
        'Units','normalized','String',num2str(top_zoneb),...
        'Callback','top_zoneb=str2double(get(top_zoneb_field,''String'')); set(top_zoneb_field,''String'',num2str(top_zoneb));');

    bot_zonet_field=uicontrol('Style','edit',...
        'Position',[.78 .80 .06 .06],...
        'Units','normalized','String',num2str(bot_zonet),...
        'Callback','bot_zonet=str2double(get(bot_zonet_field,''String'')); set(bot_zonet_field,''String'',num2str(bot_zonet));');
    bot_zoneb_field=uicontrol('Style','edit',...
        'Position',[.78 .74 .06 .06],...
        'Units','normalized','String',num2str(bot_zoneb),...
        'Callback','bot_zoneb=str2double(get(bot_zoneb_field,''String'')); set(bot_zoneb_field,''String'',num2str(bot_zoneb));');

    freq_field=uicontrol('Style','edit',...
        'Position',[.30 .50 .12 .08],...
        'Units','normalized','String',num2str(ni),...
        'Callback','ni=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ni));set(tgl2,''value'',0); set(tgl1,''value'',1)');

    freq_field0=uicontrol('Style','edit',...
        'Position',[.70 .50 .12 .08],...
        'Units','normalized','String',num2str(ra),...
        'Callback','ra=str2double(get(freq_field0,''String'')); set(freq_field0,''String'',num2str(ra)) ; set(tgl2,''value'',1); set(tgl1,''value'',0)');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.30 .40 .12 .08],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.30 .30 .12 .08],...
        'Units','normalized','String',num2str(dy),...
        'Callback','dy=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dy));');

    tgl1 = uicontrol('Style','checkbox',...
        'string','Number of Events:',...
        'Position',[.09 .50 .2 .08], 'Callback','set(tgl2,''value'',0)',...
        'Units','normalized');

    set(tgl1,'value',1);

    tgl2 =  uicontrol('Style','checkbox',...
        'string','OR: Constant Radius',...
        'Position',[.47 .50 .2 .08], 'Callback','set(tgl1,''value'',0)',...
        'Units','normalized');

    freq_field4=uicontrol('Style','edit',...
        'Position',[.30 .20 .12 .08],...
        'Units','normalized','String',num2str(Nmin),...
        'Callback','Nmin=str2double(get(freq_field4,''String'')); % set(freq_field4,''String'',num2str(Nmin));');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .15 .12 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'Callback',' inb1 =get(hndl2,''Value'');tgl1 =get(tgl1,''Value'');tgl2 =get(tgl2,''Value'');close,sel =''ca'', bdepth_ratio',...
        'String','Go');
    Nmin;
    text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.20 .75 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Please choose an Mc estimation option   ');

    mid_txt=text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.24 1.0 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Depth limits for depth ratio calculation  ');

    top_txt=text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.10 .85 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Top and bottom for TOP zone(km):');

    bot_txt=text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.40 .85 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Top and bottom for BOTTOM zone(km):');

    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.30 0.64 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String',' Grid Parameter');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.06 0.40 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in x (dx) in deg:');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.06 0.29 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in y (dy) in deg:');

    txt7 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.06 0.17 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Min. No. of events > Mc:');

    set(gcf,'visible','on');
    watchoff

end   % if nargin ==0

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'
    selgp
    itotal = length(newgri(:,1));
    if length(gx) < 4  ||  length(gy) < 4
        errordlg('Selection too small! (Dx and Dy are in degreees! ');
        return
    end

    zmap_message_center.set_info(' ','Running bdepth_ratio... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    bvg = [];
    allcount = 0.;
    nobv = 0;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','b-value grid - percent done');;
    drawnow

    % sort by depth

    % [s,is] = sort(a.Depth);
    % adepth = a(is(:,1),:);

    % find row index of ratio midpoint
    l = a.Depth >= top_zonet & a.Depth <  top_zoneb;
    top_zone = a.subset(l);

    l = a.Depth >= bot_zonet & a.Depth <  bot_zoneb;
    bot_zone = a.subset(l);

    %
    % overall b-value
    [bv magco stan av me mer me2,  pr] =  bvalca3(top_zone,inb1,inb2);
    tbo1 = bv; tno1 = length(top_zone(:,1));

    [bv magco stan av me mer me2,  pr] =  bvalca3(bot_zone,inb1,inb2);
    bbo1 = bv; bno1 = length(bot_zone(:,1));

    depth_ratio = tbo1/bbo1;

    disp(depth_ratio);

    % loop over all points
    for i= 1:length(newgri(:,1))
        x = newgri(i,1);y = newgri(i,2);
        allcount = allcount + 1.;
        i2 = i2+1;

        % calculate distance from center point and sort wrt distance
        l = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2) ;
        [s,is] = sort(l);
        b = a(is(:,1),:) ;       % re-orders matrix to agree row-wise

        if tgl1 == 0   % take point within r
            l3 = l <= ra;
            b = a.subset(l3);      % new data per grid point (b) is sorted in distanc  (from center point)
            rd = ra;
        else
            % take first ni points
            b = b(1:ni,:);      % new data per grid point (b) is sorted in distance
            l2 = sort(l); rd = l2(ni);

        end

        %estimate the completeness and b-value
        newt2 = b;

        % sort by depth

        l = b(:,7) >= top_zonet & b(:,7) <  top_zoneb;
        topb = b(l,:);

        l = b(:,7) >= bot_zonet & b(:,7) <  bot_zoneb;
        botb = b(l,:);

       if length(topb) >= Nmin  && length(botb) >= Nmin

            if inb1 == 3
                newt2 = topb;
                mcperc_ca3;
                l = topb(:,6) >= Mc90-0.05; magco = Mc90;
                n1 = length(topb(l,:));
                if length(topb(l,:)) >= Nmin
                    [topbv magco stan av me mer me2,  pr] =  bvalca3(topb(l,:),2,2);
                    [topav2 topbv2 topstan2 ] =  bmemag(topb(l,:));
                else  topbv = nan; topbv2 = nan, magco = nan; av = nan; topav2 = nan;
                    nobv = nobv + 1;
                end
                newt2 = botb;
                mcperc_ca3;
                l = botb(:,6) >= Mc90-0.05; magco = Mc90;
                n2 = length(botb(l,:));
                if length(botb(l,:)) >= Nmin
                    [botbv magco stan av me mer me2,  pr] =  bvalca3(botb(l,:),2,2);
                    [botav2 botbv2 botstan2 ] =  bmemag(botb(l,:));

                else  botbv = nan; botbv2 = nan; magco = nan; av = nan; botav2 = nan;
                    nobv = nobv + 1;
                end

            elseif inb1 == 4
                newt2 = topb;
                mcperc_ca3;
                l = topb(:,6) >= Mc95-0.05; magco = Mc95;
                n1 = length(topb(l,:));

                if length(topb(l,:)) >= Nmin
                    [topbv magco topstan topav topme topmer topme2,  toppr] =  bvalca3(topb(l,:),2,2);
                    [topav2 topbv2 topstan2 ] =  bmemag(topb(l,:));
                else
                    topbv = nan; topbv2 = nan, magco = nan; topav = nan; topav2 = nan;
                    nobv = nobv + 1;
                end
                newt2 = botb;
                mcperc_ca3;
                l = botb(:,6) >= Mc90-0.05; magco = Mc95;
                n2 = length(botb(l,:));

                if length(botb(l,:)) >= Nmin
                    [botbv magco botstan botav botme botmer botme2,  botpr] =  bvalca3(botb(l,:),2,2);
                    [botav2 botbv2 botstan2 ] =  bmemag(botb(l,:));
                else
                    botbv = nan; botbv2 = nan, magco = nan; botav = nan; botav2 = nan;
                    nobv = nobv + 1;
                end

            elseif inb1 == 5
                newt2 = topb;
                mcperc_ca3;
                if isnan(Mc95) == 0
                    magco = Mc95;
                elseif isnan(Mc90) == 0
                    magco = Mc90;
                else
                    [bv magco stan av me mer me2,  pr] =  bvalca3(b,1,1);
                end

                l = topb(:,6) >= magco-0.05;
                n1 = length(topb(l,:));
                if length(topb(l,:)) >= Nmin
                    [topbv magco topstan topav topme topmer topme2,  toppr] =  bvalca3(topb(l,:),2,2);
                    [topav2 topbv2 topstan2 ] =  bmemag(topb(l,:));
                else
                    topbv = nan; topbv2 = nan, magco = nan; topav = nan; topav2 = nan;
                    nobv = nobv + 1;
                end

                newt2 = botb;
                mcperc_ca3;
                l = botb(:,6) >= magco-0.05;
                n2 = length(botb(l,:));

                if length(botb(l,:)) >= Nmin
                    [botbv magco botstan botav botme botmer botme2,  botpr] =  bvalca3(botb(l,:),2,2);
                    [botav2 botbv2 botstan2 ] =  bmemag(botb(l,:));
                else
                    botbv = nan; botbv2 = nan, magco = nan; botav = nan; bottopav2 = nan;
                    nobv = nobv + 1;
                end

            elseif inb1 == 1
                % newt2 = topb;
                % mcperc_ca3;
                [topbv magco topstan topav topme topmer topme2,  toppr] =  bvalca3(topb,1,1);
                l = topb(:,6) >= magco-0.05;
                if length(topb(l,:)) >= Nmin
                    [topav2 topbv2 topstan2 ] =  bmemag(topb(l,:));
                    n1 = length(topb(l,:));
                else
                    topbv = nan; topbv2 = nan; magco = nan; topav = nan; topav2 = nan;
                    n1 = nan;

                    nobv = nobv + 1;
                end
                %    newt2 = botb;
                %    mcperc_ca3;
                [botbv magco botstan botav botme botmer botme2,  botpr] =  bvalca3(botb,1,1);
                l = botb(:,6) >= magco-0.05;
                if length(botb(l,:)) >= Nmin
                    [botav2 botbv2 botstan2 ] =  bmemag(botb(l,:));
                    n2 = length(botb(l,:));

                else
                    botbv = nan; botbv2 = nan; magco = nan; botav = nan; botav2 = nan;
                    n2 = nan;
                    nobv = nobv + 1;
                end

            elseif inb1 == 2
                [topbv magco topstan topav topme topmer topme2,  toppr] =  bvalca3(topb,2,2);
                n1 = length(topb);
                [topav2 topbv2 topstan2 ] =  bmemag(topb);

                [botbv magco botstan botav botme botmer botme2,  botpr] =  bvalca3(botb,2,2);

                n2 = length(botb);
                [botav2 botbv2 botstan2 ] =  bmemag(botb);

            end

        else
            topbv = nan; topbv2 = nan; magco = nan; topav = nan; topav2 = nan;
            botbv = nan; botbv2 = nan; magco = nan; botav = nan; botav2 = nan;
            nobv = nobv + 1;

        end

        bv = topbv/botbv; bv2 = topbv2/botbv2; stan2 = topstan2/botstan2; av = topav/botav;
        stan = topstan/botstan;
        n = n1+n2; bo1 = topbv; b2 = botbv;
        da = -2*n*log(n) + 2*n1*log(n1+n2*bo1/b2) + 2*n2*log(n1*b2/bo1+n2) - 2;
        pr = (1  -  exp(-da/2-2))*100;

        bvg = [bvg ; bv magco x y rd bv2 av pr topbv botbv];
        waitbar(allcount/itotal)

    end

    snobv = num2str(nobv);
    spnobv = num2str((nobv/allcount)*100.0);
    disp(['Not enough EQs to calculate b values for ', snobv, ' gridnodes. This is ',spnobv,'% of the total'])

    save tmpncal.mat
    quit

    % for newgr
    %save cnssgrid.mat
    %quit
    % save data
    %
    catSave3 =...
        [ 'zmap_message_center.set_info(''Save Grid'',''  '');think;',...
        '[file1,path1] = uiputfile(fullfile(hodi, ''eq_data'', ''*.mat''), ''Grid Datafile Name?'') ;',...
        ' sapa2 = [''save '' path1 file1 '' bvg gx gy dx dy par1 tdiff t0b teb a main faults mainfault coastline yvect xvect tmpgri ll depth_ratio top_zonet top_zoneb bot_zoneb bot_zonet''];',...
        ' if length(file1) > 1, eval(sapa2),end , done']; eval(catSave3)

    close(wai)
    watchoff

    % plot the results
    % old and re3 (initially ) is the b-value matrix
    %
    normlap2=ones(length(tmpgri(:,1)),1)*nan;
    normlap2(ll)= bvg(:,1);
    re3=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,5);
    r=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,6);
    meg=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,2);
    old1=reshape(normlap2,length(yvect),length(xvect));

    %  normlap2(ll)= bvg(:,7);
    %  pro=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,7);
    avm=reshape(normlap2,length(yvect),length(xvect));

    %  normlap2(ll)= bvg(:,9);
    %  stanm=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,8);
    Prmap=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,9);
    top_b=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,10);
    bottom_b=reshape(normlap2,length(yvect),length(xvect));

    old = re3;

    % View the b-value map
    view_bdepth

end   % if sel = na

% Load exist b-grid
if sel == 'lo'
    [file1,path1] = uigetfile(['*.mat'],'b-value gridfile');
    if length(path1) > 1
        think
        load([path1 file1])
        normlap2=ones(length(tmpgri(:,1)),1)*nan;

        normlap2(ll)= bvg(:,1);
        re3=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,5);
        r=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,6);
        meg=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,2);
        old1=reshape(normlap2,length(yvect),length(xvect));

        %  normlap2(ll)= bvg(:,7);
        %  pro=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,7);
        avm=reshape(normlap2,length(yvect),length(xvect));

        %  normlap2(ll)= bvg(:,9);
        % stanm=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,8);
        Prmap=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,9);
        top_b=reshape(normlap2,length(yvect),length(xvect));

        normlap2(ll)= bvg(:,10);
        bottom_b=reshape(normlap2,length(yvect),length(xvect));

        old = re3;

        view_bdepth
    else
        return
    end
end


text:
function [ know, x ] = p22_sol ( n )

%*****************************************************************************80
%
%% P22_SOL returns the solution for problem 22.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 December 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = zeros ( n, 1 );

  return
end


text:
function [y,deriv] = bsxtimes(w,m,n)
% This is an MV2DF
%
% w = [vec(A); vec(b) ] --> vec(bsxfun(@times,A,b)),
%
%      where A is an m-by-n matrix and
%            b is a 1-by-n row.
%

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w) bsxtimes(w,m,n);
    return;
end

if isa(w,'function_handle')
    f = bsxtimes([],m,n);
    y = compose_mv(f,w,[]);
    return;
end

[A,b] = extract(w,m,n);
y = bsxfun(@times,A,b);
y = y(:);

deriv = @(Dy) deriv_this(Dy,A,b);

function  [g,hess,linear] = deriv_this(Dy,A,b)
g = gradient(Dy,A,b);
linear = false;

hess = @(v) hess_this(v,Dy,A,b);

function [h,Jv] = hess_this(v,Dy,A,b)
[m,n] = size(A);
[vA,vb] = extract(v,m,n);
h = gradient(Dy,vA,vb);
if nargout>1
  Jv = bsxfun(@times,vA,b);
  Jv = Jv + bsxfun(@times,A,vb);
  Jv = Jv(:);
end

function [A,b] = extract(w,m,n)
A = reshape(w(1:m*n),m,n);
b = w(m*n+1:end).';

function g = gradient(Dy,A,b)
Dy = reshape(Dy,size(A));
gA = bsxfun(@times,Dy,b);
gb = sum(Dy.*A,1);
g = [gA(:);gb(:)];

function test_this()
m = 5;
n = 10;
A = randn(m,n);
b = randn(1,n);
w = [A(:);b(:)];

f = bsxtimes([],m,n);
test_MV2DF(f,w);

author: nesl