text:
function opts = updateMovie(S, dt, p, opts, t, v, compGrid, plotGrid)
%UPDATEMOVIE   Update the movie when solving a PDE specified by a SPINOP.

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

% Set-up:
nVars = S.numVars;
Ylim = opts{1};
dataToPlot = opts{2};
xx = compGrid{1};
xxx = plotGrid{1};
N = size(xx, 1) - 1;
Nplot = size(xxx, 1) - 1;

for k = 1:nVars

    % Extract each variable:
    idx = (k-1)*N + 1;
    vv = dataToPlot(v(idx:idx+N-1));
    vv = [vv; vv(1)]; %#ok<*AGROW> add repeated values (periodic endpoints)

    % Change axes if necessary:
    if ( nargout == 1 )
        minvnew = min(dataToPlot(vv));
        maxvnew = max(dataToPlot(vv));
        if ( maxvnew > Ylim(2*(k-1) + 2) )
            vscalenew = max(abs(minvnew), maxvnew);
            Ylim(2*(k-1) + 2) = maxvnew + .1*vscalenew;
        end
        if ( minvnew < Ylim(2*(k-1) + 1) )
            vscalenew = max(abs(minvnew), maxvnew);
            Ylim(2*(k-1) + 1) = minvnew - .1*vscalenew;
        end
    end

    % Interpolate each variable on a finer grid:
    if ( Nplot > N )
        vvv = interp1(xx, vv, xxx, 'spline');
    else
        vvv = vv;
    end

    % Update each variable:
    set(p{1,k}, 'ydata', vvv)
    set(p{1,k}.Parent, 'ylim', [Ylim(2*(k-1) + 1), Ylim(2*(k-1) + 2)])

    % Update each title:
    titleString = sprintf('N = %i (DoFs = %i), dt = %1.1e, t = %.4f', N, ...
        nVars*N, dt, t);
    set(p{2,k}, 'String', titleString)
    drawnow

end

% Update outputs:
opts{1} = Ylim;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [vtoe,nne]=VertToEle(t)
np = max(t(:)); ne = size(t,1);
nne = zeros(np,1);
for ie = +1 : ne                                                            %--Go through once and find max nnz.
    %for iv = +1 : +3
    %kjr loop unrolled.
    nm1=t(ie,1);
    nm2=t(ie,2);
    nm3=t(ie,3);
    %vrtx = t(ie,iv);
    nne(nm1,1) = nne(nm1,1) + 1;
    nne(nm2,1) = nne(nm2,1) + 1;
    nne(nm3,1) = nne(nm3,1) + 1;
    %end
end
mnz  = max(nne);                                                           %--max number of non-zeros
vtoe = zeros(mnz,np);                                                      %--vertex to element connectivity
nne = zeros(np,1);                                                         %--number of neighboring elements
for ie = +1 : ne
    nm1=t(ie,1);
    nm2=t(ie,2);
    nm3=t(ie,3);

    %kjr loop unrolled.
    nne(nm1,1) = nne(nm1,1) +1;
    nne(nm2,1) = nne(nm2,1) +1;
    nne(nm3,1) = nne(nm3,1) +1;

    vtoe(nne(nm1,1),nm1) = ie;
    vtoe(nne(nm2,1),nm2) = ie;
    vtoe(nne(nm3,1),nm3) = ie;

    %vtoe(nm1,nne(nm1,1)) = ie;
    %vtoe(nm2,nne(nm2,1)) = ie;
    %vtoe(nm3,nne(nm3,1)) = ie;
end
nne  = nne';
%vtoe = vtoe';

end

author: CHLNDDEV


text:
function [norot, rot] = nc2010_plot_rppca_toyexperiment(n,m,d,ncomps)
%
% nc2010_plot_toyexperiment(model,n,m,d,ncomps,flatW)
%

model = 'rppca';

filename = sprintf(['/home/jluttine/matlab/neurocomputing2010/' ...
                    'nc2010_%s_toyexperiment_n=%d_m=%d_d=%d_ncomps=%d'], ...
                   model,n,m,d,ncomps);

norot = load([filename,'_rotate=0.mat'], 'duration', 'cost', 'rmse', 'rmse_test');
rot = load([filename,'_rotate=1.mat'], 'duration', 'cost', 'rmse', 'rmse_test');
plot_results(norot, rot, n,m,d,ncomps);

if nargout < 1
  clear norot;
  clear rot;
end
return
% $$$
% $$$ norot = load('rotationpaper_experiment_n=100_m=10_d=3_ncomps=9_rotate=0.mat');
% $$$ rot = load('rotationpaper_experiment_n=100_m=10_d=3_ncomps=9_rotate=1.mat');
% $$$ plot_results(norot, rot, 100,10,3,9);
% $$$
% $$$ norot = load('rotationpaper_experiment_n=1000_m=50_d=10_ncomps=49_rotate=0.mat');
% $$$ rot = load('rotationpaper_experiment_n=1000_m=50_d=10_ncomps=49_rotate=1.mat');
% $$$ plot_results(norot, rot, 1000,50,10,49);
% $$$
% $$$ norot = load('rotationpaper_experiment_n=1000_m=100_d=10_ncomps=30_rotate=0.mat');
% $$$ rot = load('rotationpaper_experiment_n=1000_m=100_d=10_ncomps=30_rotate=1.mat');
% $$$ plot_results(norot, rot, 1000,100,10,30);
% $$$
% $$$ return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_results(norot, rot, n,m,d,ncomps)

% $$$ norot.cost(norot.cost==-inf) = -1e10;
% $$$ rot.cost(rot.cost==-inf) = -1e10;
suffix = sprintf(' (n=%d m=%d d=%d ncomps=%d)',n,m,d,ncomps);

figure
semilogx([norot.duration(:), rot.duration(:)], -[norot.cost(:), rot.cost(:)])
%title(['Negative loglikelihood', suffix]);
filename = sprintf('fig_cost_n=%d_m=%d_d=%d_ncomps=%d', n,m,d,ncomps);
set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1:2), 8,8]);
pos = get(gcf, 'paperposition');
set(gcf, 'paperposition', [pos(1:2),8,8])

figure
semilogx([norot.duration(:), rot.duration(:)], [norot.rmse(:), rot.rmse(:)])
%title(['Training set RMSE', suffix]);
filename = sprintf('fig_rmse_n=%d_m=%d_d=%d_ncomps=%d', n,m,d,ncomps);
set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1:2), 8,8]);
pos = get(gcf, 'paperposition');
set(gcf, 'paperposition', [pos(1:2),8,8])

figure
semilogx([norot.duration(:), rot.duration(:)], [norot.rmse_test(:), rot.rmse_test(:)])
%title(['Test set RMSE', suffix]);
filename = sprintf('fig_rmsetest_n=%d_m=%d_d=%d_ncomps=%d', n,m,d,ncomps);
set(gcf, 'units', 'centimeters', 'paperunits', 'centimeters');
pos = get(gcf, 'position');
set(gcf, 'position', [pos(1:2), 8,8]);
pos = get(gcf, 'paperposition');
set(gcf, 'paperposition', [pos(1:2),8,8])

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [cfg] = ft_topoplotCC(cfg, freq)

% FT_TOPOPLOTCC plots the coherence between channel pairs
%
% Use as
%  ft_topoplotCC(cfg, freq)
%
% The configuration should contain:
%   cfg.feedback    = string (default = 'textbar')
%   cfg.layout      = specification of the layout, see FT_PREPARE_LAYOUT
%   cfg.foi         = the frequency of interest which is to be plotted (default is the first frequency bin)
%   cfg.widthparam  = string, parameter to be used to control the line width (see below)
%   cfg.alphaparam  = string, parameter to be used to control the opacity (see below)
%   cfg.colorparam  = string, parameter to be used to control the line color
%
% The widthparam should be indicated in pixels, e.g. usefull numbers are 1
% and larger.
%
% The alphaparam should be indicated as opacity between 0 (fully transparent)
% and 1 (fully opaque).
%
% The default is to plot the connections as lines, but you can also use
% bidirectional arrows:
%    cfg.arrowhead    = string, 'none', 'stop', 'start', 'both' (default = 'none')
%    cfg.arrowsize    = scalar, size of the arrow head in figure units,
%                       i.e. the same units as the layout (default is automatically determined)
%    cfg.arrowoffset  = scalar, amount that the arrow is shifted to the side in figure units,
%                       i.e. the same units as the layout (default is automatically determined)
%    cfg.arrowlength  = scalar, amount by which the length is reduced relative to the complete line (default = 0.8)
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure. For this particular function, the input should be
% structured as a cell-array.
%
% See also FT_PREPARE_LAYOUT, FT_MULTIPLOTCC, FT_CONNECTIVITYPLOT

% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar freq
ft_preamble provenance freq
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
freq = ft_checkdata(freq, 'cmbrepresentation', 'sparse');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'required', {'foi', 'layout'});

% set the defaults
cfg.feedback    = ft_getopt(cfg, 'feedback',    'text');
cfg.alphaparam  = ft_getopt(cfg, 'alphaparam',  []);
cfg.widthparam  = ft_getopt(cfg, 'widthparam',  []);
cfg.colorparam  = ft_getopt(cfg, 'colorparam',  'cohspctrm');
cfg.newfigure   = ft_getopt(cfg, 'newfigure',   'yes');
cfg.arrowhead   = ft_getopt(cfg, 'arrowhead',   'none');  % none, stop, start, both
cfg.arrowsize   = ft_getopt(cfg, 'arrowsize',   nan);     % length of the arrow head, should be in in figure units, i.e. the same units as the layout
cfg.arrowoffset = ft_getopt(cfg, 'arrowoffset', nan);     % absolute, should be in figure units, i.e. the same units as the layout
cfg.arrowlength = ft_getopt(cfg, 'arrowlength', 0.8);     % relative to the complete line
cfg.linestyle   = ft_getopt(cfg, 'linestyle',   []);
cfg.colormap    = ft_getopt(cfg, 'colormap',    colormap);
cfg.renderer    = ft_getopt(cfg, 'renderer',    []);      % let MATLAB decide on the default

tmpcfg = keepfields(cfg, {'layout', 'rows', 'columns', 'commentpos', 'scalepos', 'elec', 'grad', 'opto', 'showcallinfo'});
lay = ft_prepare_layout(tmpcfg, freq);

beglabel = freq.labelcmb(:,1);
endlabel = freq.labelcmb(:,2);
ncmb     = size(freq.labelcmb,1);

% select the data to be used in the figure
fbin = nearest(freq.freq, cfg.foi);

if isfield(freq, cfg.widthparam)
  widthparam = freq.(cfg.widthparam)(:,fbin);
else
  widthparam = ones(ncmb,1);
end

if isfield(freq, cfg.alphaparam)
  alphaparam = freq.(cfg.alphaparam)(:,fbin);
else
  alphaparam = [];
end

if isfield(freq, cfg.colorparam)
  colorparam = freq.(cfg.colorparam)(:,fbin);
else
  colorparam = [];
end

if strcmp(cfg.newfigure, 'yes')
  figure
end

hold on
axis equal

if isnan(cfg.arrowsize)
  % use the size of the figure to estimate a decent number
  siz = axis;
  cfg.arrowsize = (siz(2) - siz(1))/50;
  ft_warning('using an arrowsize of %f', cfg.arrowsize);
end

if isnan(cfg.arrowoffset)
  % use the size of the figure to estimate a decent number
  siz = axis;
  cfg.arrowoffset = (siz(2) - siz(1))/100;
  ft_warning('using an arrowoffset of %f', cfg.arrowoffset);
end

rgb = cfg.colormap;
if ~isempty(colorparam)
  cmin = min(colorparam(:));
  cmax = max(colorparam(:));

  % this line creates a sorting vector that cause the most extreme valued arrows to be plotted last
  [srt, srtidx] = sort(abs(colorparam));

  colorparam = (colorparam - cmin)./(cmax-cmin);
  colorparam = round(colorparam * (size(rgb,1)-1) + 1);
end

if strcmp(cfg.newfigure, 'yes')
  ft_plot_layout(lay, 'label', 'no', 'box', 'off');
end

% fix the limits for the axis
axis(axis);

if ~exist('srtidx', 'var')
  srtidx = 1:ncmb;
end

ft_progress('init', cfg.feedback, 'plotting connections...');

for i=srtidx(:)'

  if strcmp(beglabel{i}, endlabel{i})
    % skip autocombinations
    continue
  end

  if widthparam(i)>0 && (isempty(alphaparam)||alphaparam(i)>0)
    ft_progress(i/ncmb, 'plotting connection %d from %d (%s -> %s)\n', i, ncmb, beglabel{i}, endlabel{i});

    begindx = strcmp(beglabel{i}, lay.label);
    endindx = strcmp(endlabel{i}, lay.label);
    xbeg = lay.pos(begindx,1);
    ybeg = lay.pos(begindx,2);
    xend = lay.pos(endindx,1);
    yend = lay.pos(endindx,2);

    if isempty(cfg.linestyle)
      if strcmp(cfg.arrowhead, 'none')
        x = [xbeg xend]';
        y = [ybeg yend]';
        % h = line(x, y);
        h = patch(x, y, 1);
      else
        arrowbeg  = [xbeg ybeg];
        arrowend  = [xend yend];
        center    = (arrowbeg+arrowend)/2;
        direction = (arrowend - arrowbeg);
        direction = direction/norm(direction);
        offset    = [direction(2) -direction(1)];
        arrowbeg  = cfg.arrowlength * (arrowbeg-center) + center + cfg.arrowoffset * offset;
        arrowend  = cfg.arrowlength * (arrowend-center) + center + cfg.arrowoffset * offset;

        h = arrow(arrowbeg, arrowend, 'Ends', cfg.arrowhead, 'length', 0.05);

      end % if arrow
      if ~isempty(widthparam)
        set(h, 'LineWidth', widthparam(i));
      end

      if ~isempty(alphaparam)
        set(h, 'EdgeAlpha', alphaparam(i));
        set(h, 'FaceAlpha', alphaparam(i)); % for arrowheads
      end

      if ~isempty(colorparam)
        set(h, 'EdgeColor', rgb(colorparam(i),:));
        set(h, 'FaceColor', rgb(colorparam(i),:)); % for arrowheads
      end
    elseif ~isempty(cfg.linestyle)

      % new style of plotting, using curved lines, this does not allow for
      % alpha mapping on the line segment
      switch cfg.linestyle
        case 'curve'
          tmp = cscvn([xbeg mean([xbeg,xend])*0.8 xend;ybeg mean([ybeg,yend])*0.8 yend]);
          pnt = fnplt(tmp);
          h   = line(pnt(1,:)', pnt(2,:)');

          if ~isempty(widthparam)
            set(h, 'LineWidth', widthparam(i));
          end
          if ~isempty(colorparam)
            set(h, 'Color', rgb(colorparam(i),:));
          end

          % deal with the arrow
          if ~strcmp(cfg.arrowhead, 'none')
            arrowbeg  = pnt(:,1)';
            arrowend  = pnt(:,end)';
            directionbeg = (arrowbeg - pnt(:,2)');
            directionend = (arrowend - pnt(:,end-1)');

            directionbeg = directionbeg/norm(directionbeg);
            directionend = directionend/norm(directionend);

            switch cfg.arrowhead
              case 'stop'
                pnt1 = arrowend - 0.05*directionend + 0.02*[directionend(2) -directionend(1)];
                pnt2 = arrowend;
                pnt3 = arrowend - 0.05*directionend - 0.02*[directionend(2) -directionend(1)];
                h_arrow = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', [0 0 0]);

              case 'start'
                pnt1 = arrowbeg - 0.05*directionbeg + 0.02*[directionbeg(2) -directionbeg(1)];
                pnt2 = arrowbeg;
                pnt3 = arrowbeg - 0.05*directionbeg - 0.02*[directionbeg(2) -directionbeg(1)];
                h_arrow = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', [0 0 0]);

              case 'both'
                pnt1 = arrowbeg - 0.05*directionbeg + 0.02*[directionbeg(2) -directionbeg(1)];
                pnt2 = arrowbeg;
                pnt3 = arrowbeg - 0.05*directionbeg - 0.02*[directionbeg(2) -directionbeg(1)];
                h_arrow(1) = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', [0 0 0]);

                pnt1 = arrowend - 0.05*directionend + 0.02*[directionend(2) -directionend(1)];
                pnt2 = arrowend;
                pnt3 = arrowend - 0.05*directionend - 0.02*[directionend(2) -directionend(1)];
                h_arrow(2) = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', [0 0 0]);

            end
          else
            h_arrow = [];
          end

          if ~isempty(widthparam)
            set(h, 'LineWidth', widthparam(i));
            if ~isempty(h_arrow)
              set(h_arrow, 'LineWidth', widthparam(i));
            end
          end
          if ~isempty(colorparam)
            set(h, 'Color', rgb(colorparam(i),:));
            if ~isempty(h_arrow)
              set(h_arrow, 'Edgecolor', rgb(colorparam(i),:));
              set(h_arrow, 'Facecolor', rgb(colorparam(i),:));
            end
          end

          if ~isempty(alphaparam)
            if ~isempty(h_arrow)
              set(h_arrow, 'EdgeAlpha', alphaparam(i));
              set(h_arrow, 'FaceAlpha', alphaparam(i)); % for arrowheads
            end
          end

        otherwise
          ft_error('unsupported linestyle specified');
      end

    end % if empty linestyle
  end
end % for srtindx

ft_progress('close');

% this is needed for the figure title
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

% set the figure window title
if ~isempty(dataname)
  set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), mfilename, join_str(', ', dataname)));
else
  set(gcf, 'Name', sprintf('%d: %s', double(gcf), mfilename));
end
set(gcf, 'NumberTitle', 'off');

% set renderer if specified
if ~isempty(cfg.renderer)
  set(gcf, 'renderer', cfg.renderer)
end

% improve the fit in the axis
axis tight

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous freq
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for plotting arrows, see also fieldtrip/private/arrow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = arrow(arrowbeg, arrowend, varargin)
ends   = ft_getopt(varargin, 'ends');
length = ft_getopt(varargin, 'length'); % the length of the arrow head, in figure units
color  = [0 0 0]; % in RGB

direction = (arrowend - arrowbeg);
direction = direction/norm(direction);
offset    = [direction(2) -direction(1)];

pnt1 = arrowbeg;
pnt2 = arrowend;
h = patch([pnt1(1) pnt2(1)], [pnt1(2) pnt2(2)], color);

switch ends
  case 'stop'
    pnt1 = arrowend - length*direction + 0.4*length*offset;
    pnt2 = arrowend;
    pnt3 = arrowend - length*direction - 0.4*length*offset;
    h(end+1) = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', color);

  case 'start'
    pnt1 = arrowbeg + length*direction + 0.4*length*offset;
    pnt2 = arrowbeg;
    pnt3 = arrowbeg + length*direction - 0.4*length*offset;
    h(end+1) = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', color);

  case 'both'
    pnt1 = arrowend - length*direction + 0.4*length*offset;
    pnt2 = arrowend;
    pnt3 = arrowend - length*direction - 0.4*length*offset;
    h(end+1) = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', color);

    pnt1 = arrowbeg + length*direction + 0.4*length*offset;
    pnt2 = arrowbeg;
    pnt3 = arrowbeg + length*direction - 0.4*length*offset;
    h(end+1) = patch([pnt1(1) pnt2(1) pnt3(1)]', [pnt1(2) pnt2(2) pnt3(2)]', color);

  case 'none'
    % don't draw arrow heads
end

author: spm


text:
%%
% Test for power diagrams

addpath('../toolbox/');
addpath('power_bounded/');
addpath('power_diagrams/');

rep = 'results/semi-discrete/';
[~,~] = mkdir(rep);

% number of sites
N = 200;
% x and y coordinate of the sites
xy = rand(N, 2);
sigma = .1;
xy = sigma*randn(N, 2)+1/2;

% weights
w = ones(N,1);

% the bounding box in clockwise order
bb = [0,0; 0,1; 1,1; 1,0];
% get the power diagram
[V,C] = power_bounded(xy(:,1),xy(:,2), w, bb);
% [PD, PDinf] = powerDiagramWrapper(xy, w);
% draw the resulted power diagram
plot_power(xy,V,C,bb);

q = .01;
axis([-q 1+q -q 1+q]);

% target histogram
Atgt = ones(N,1)/N;

% optimization run
tau = .0004 * N; % gradient step size
niter = 200;
iter_disp = round( [1:8 niter/16 niter/8 niter/4 niter/2 niter] );
iter_disp = 1:2:niter;
err = []; kdisp = 0;
for i=1:niter
    progressbar(i,niter);
    [V,C] = power_bounded(xy(:,1),xy(:,2), w, bb);
    A = area_power(xy,V,C,bb);
    if intersect(i,iter_disp)
        kdisp = kdisp+1;
        clf; plot_power(xy,V,C,bb); axis([-q 1+q -q 1+q]);
        saveas(gcf, [rep 'iteration-' num2str(kdisp), '.eps'], 'epsc');
        drawnow;
    end
    A = A/sum(A); % be sure to be normalized ...
    wnew = w - tau * (A(:)-Atgt(:));
    err(i) = norm(w-wnew);
    w = wnew;
end

% error decay
plot(err); axis tight;

% plot final OT
clf; plot_power(xy,V,C,bb, 2);
saveas(gcf, [rep 'matching.eps'], 'epsc');
axis([-q 1+q -q 1+q]);

author: optimaltransport


text:
function [chd,Ichd] = ffmSubdivide(X,prt,edg,fig)
%+========================================================================+
%|                                                                        |
%|         OPENFFM - LIBRARY FOR FAST AND FREE MEMORY CONVOLUTION         |
%|           openFfm is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : ffmSubdivide.m                                |
%|    #    |   VERSION    : 0.6                                           |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Recursive subdivision using an octree         |
%|  `---'  |                                                              |
%+========================================================================+

% Initialisation des enfants par subdivision type Octree
Nprt    = length(prt.ind);
chd.ctr = zeros(8*Nprt,3);
chd.ind = cell(8*Nprt,1);
chd.nbr = zeros(8*Nprt,1);
Ichd    = zeros(Nprt,8);

% Boucle sur les parents
l = 0;
for i = 1:Nprt
    % Subdivision centrale
    ind  = prt.ind{i};
    cutx = X(ind,1)<=prt.ctr(i,1);
    cuty = X(ind,2)<=prt.ctr(i,2);
    cutz = X(ind,3)<=prt.ctr(i,3);

    % Repartition des indices
    cut = logical([ ...
        cutx  .*  cuty  .*  cutz , ...
        ~cutx  .*  cuty  .*  cutz , ...
        cutx  .* ~cuty  .*  cutz , ...
        ~cutx  .* ~cuty  .*  cutz , ...
        cutx  .*  cuty  .* ~cutz , ...
        ~cutx  .*  cuty  .* ~cutz , ...
        cutx  .* ~cuty  .* ~cutz , ...
        ~cutx  .* ~cuty  .* ~cutz ] );
    Ncut = sum(cut,1);

    % Centres des enfants
    ctr = ones(8,1)*prt.ctr(i,1:3) + 0.25*edg .* [...
        -1 -1 -1 ; 1 -1 -1; -1 1 -1 ; 1 1 -1 ;...
        -1 -1  1 ; 1 -1  1; -1 1  1 ; 1 1  1 ];

    % Ajout d'une boite pour chaque enfant non vide
    for j = find(Ncut)
        chd.ctr(l+1,:) = ctr(j,:);
        chd.ind{l+1}   = ind(cut(:,j));
        chd.nbr(l+1)   = Ncut(j);
        Ichd(i,j)      = l+1;
        l = l + 1;
    end
end

% Extraction des boites vides
chd.ctr = chd.ctr(1:l,:);
chd.ind = chd.ind(1:l);
chd.nbr = chd.nbr(1:l);

% Verification que chaque pot a son couvercle
if abs( sum(chd.nbr) - sum(prt.nbr) ) > 1e-12
    error('ffmSubDivide.m - error 1')
end
if abs( sum(cell2mat(chd.ind)) - sum(cell2mat(prt.ind)) ) > 1e-12
    error('ffmSubDivide.m - error 2')
end

% Representation graphique
if fig
    % Configuration
    figure(fig); clf
    hold on
    title('Box plot')
    grid on ; axis equal

    % Representation des centres
    plot3(chd.ctr(:,1),chd.ctr(:,2),chd.ctr(:,3),'ok','MarkerSize',8)

    % Representation des points par boites
    for i = 1:length(chd.ind)
        plot3(X(chd.ind{i},1),X(chd.ind{i},2),X(chd.ind{i},3),'.','Color',rand(1,3))
    end
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function domisfit(catalog,sig,plu,az,phi,R)
    % domisfit calculates the misfit for each EQ to a given stress tensor orientation.
    % The actual calculation is done using a call to a fortran program.
    %
    % Stefan Wiemer 07/95

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    global mi mif1 mif2 newcat2 mi2
    global cumu2
    report_this_filefun();

    hodis = fullfile(ZG.hodi, 'external');
    cd(hodis);

    % prepare the focal; mechnism in Gephard format ...
    tmp = [catalog(:,10:12) ];
    tmp(tmp(:,2) >89.999,2) = 89.;

    try
        save data.inp tmp -ascii
    catch ME
        error_handler(ME, ['Error - could not save file ' ZmapGlobal.Data.Directories.output 'data.inp - permission?']);
    end

    infi =  'data.inp';
    outfi = 'tmpin.dat';
    fid = fopen('inmifi.dat','w');
    fprintf(fid,'%s\n',infi);
    fprintf(fid,'%s\n',outfi);
    fclose(fid);

    delete(outfi)

    [status, result] = system('datasetupDD < inmifi.dat');

    fid = ('tmpin.dat');
    format = '%f%f%f%f%f';
    %[d1, d2, d3, d4, d5] = textread(fid,format,'headerlines',1);
    C = textscan(fid,format,'HeaderLines',1); %Problem: "Errorlines" cause crashes.
    dall=[C{:}];

    %dall = [d1, d2, d3, d4, d5];
    save tmpin.dat dall -ascii

    infi = 'tmpin.dat';
    outfi = 'tmpout.dat';

    fid = fopen('inmifi.dat','w');

    fprintf(fid,'%s\n',infi);
    fprintf(fid,'%s\n',outfi);
    fprintf(fid,'%2.0f\n',sig);
    fprintf(fid,'%6.2f\n',plu);
    fprintf(fid,'%6.2f\n',az);
    fprintf(fid,'%6.2f\n',phi);
    fprintf(fid,'%3.2f\n',R);
    le = catalog.Count;
    fprintf(fid,'%6i\n',le);

    fclose(fid);
    try
        delete outfi
    catch ME
        warning(ME.message);
    end

    comm = 'testfm < inmifi.dat'
    [status,result]=system(comm)
    try
        % CGR: it looks like the format will be:
        % line 1: ndata kdata  [where ndata is # of data, kdata is # of fault planes]
        % lines 2-end : az1 dip1 az2 dip2 wt
        % load('tmpout.dat')
        s=importdata('../external/tmpout.dat', ' ', 1);

        headernumbers=str2num(s.textdata{1});
        nData=headernumbers(1);
        kData=headernumbers(2);

        mi = s.data; % probably as [az1, dip1, az2, dip2, wt ; ...]
    catch ME
        warning(ME.message)
    end

    % mi = tmpout; % mi gets results from the fortran progarm

    misfitAngle=mi(:,2);

    mif1=findobj('Type','Figure','-and','Name','Misfit Map');

    if isempty(mif1)
        mif1 = figure_w_normalized_uicontrolunits( ...
            'Name','Misfit Map',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'NextPlot','add', ...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        create_my_menu_1();

        labelList=['Size | Size + Thickness | Size +Thickness +color  '];
        labelPos = [0.2 0.93 0.35 0.05];
        hndl2=uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'String',labelList,...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@callbackfun_005);

        decimationList=['1 | 1/2 | 1/3 | 1/4 | 1/5 | 1/6| 1/7| 1/8 | 1/9 | 1/10'];
        labelPos = [0.9 0.93 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',4,...
            'String',decimationList,...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@store_decimation_cb);

        uicontrol(...
            'Style','pushbutton',...
            'Units','normalized',...
            'Position',[0.9 0.6 0.08 0.08],...
            'String','X-sec',...
            'callback',@callbackfun_007);
        set(gca,'NextPlot','add')
        %end killed
        uicontrol(...
            'Style','pushbutton',...
            'Units','normalized',...
            'Position',[0.9 0.7 0.08 0.08],...
            'String','Map',...
            'callback',@callbackfun_008);
        set(gca,'NextPlot','add')
    end

    figure(mif1)

    plotmima(4, mi)

    mif2=findobj('Type','Figure','-and','Name','Misfit ');

    if isempty(mif2)
        mif2 = figure_w_normalized_uicontrolunits( ...
            'Name','Misfit ',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'NextPlot','add', ...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        create_my_menu_2();
        listFields={'Longitude','Latitude','Time','Magnitude','Depth','Strike','Default'};
        labelPos = [0.7 0.9 0.25 0.08];
        hFieldChoice=uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'String',listFields,...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@callbackfun_012);
        set(gca,'NextPlot','add')
    end

    figure(mif2)
    delete(findobj(mif2,'Type','axes'));

    plotmi(listFields{1}, catalog, mi)

    %% ui functions
    function create_my_menu_1() %TODO rename to something more intelligent
        add_menu_divider();    %
        omp2= uimenu('Label','Tools');
        uimenu(omp2,'label','Misfit-Magnitude',...
            'MenuSelectedFcn',@cb_misfitmag);
        uimenu(omp2,'label','Misfit-Depth',...
            'MenuSelectedFcn',@cb_misfitdep);
        uimenu(omp2,'label','Earthquake-Depth',...
            'MenuSelectedFcn',@cb_eqdep);
        uimenu(omp2,'label','Earthquake-Strike',...
            'MenuSelectedFcn',@cb_eqstrike);
        %
    end

    function create_my_menu_2() %TODO rename to something more intelligent
        add_menu_divider();
        omp1= uimenu('Label','Tools');
        uimenu(omp1,'label','Save sorted catalog',...
            'MenuSelectedFcn',@callbackfun_009);
        uimenu(omp1,'label','AS Function',...
            'MenuSelectedFcn',@cb_astmisfit);
        uimenu(omp1,'label','Compare',...
            'MenuSelectedFcn',@cb_comparemisfit);
    end

    %% callback functions
    function cb_misfitmag(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        mi_ma(misfitAngle);
    end

    function cb_misfitdep(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        mi_dep(misfitAngle);
    end

    function cb_eqdep(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        eq_dep(misfitAngle);
    end

    function cb_eqstrike(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        earthquake_strike();
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        in2=hndl2.Value;
        plotmima(in2, mi);
    end

    function store_decimation_cb(mysrc,myevt)
        global oneOfHowManyPopupIdx
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        in3=mysrc.Value;
        oneOfHowManyPopupIdx=in3;
        in2=hndl2.Value;
        plotmima(in2, mi) ;
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        var1 = 3;
        plotmimac(mi, inde); % No idea what inde is or where it comes from
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        var1 = 1;
        mifigrid(var1,mi);
    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        save_sortpere;
    end

    function cb_astmisfit(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ast_misfit();
    end

    function cb_comparemisfit(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        compare_misfit();
    end

    function callbackfun_012(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        in2=hFieldChoice.Value;
        plotmi(listfields{in2}, catalog, mi);
    end

end

function earthquake_strike()
    % plot the earthquake number along the strike on the map view
    %	August 1995 by Zhong Lu

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();
    myFigName='Earthquake Number Map';
    mif55=findobj('Type','Figure','-and','Name',myFigName);

    if isempty(mif55)
        mif55 = figure_w_normalized_uicontrolunits( ...
            'Name',myFigName,...
            'NumberTitle','off', ...
            'backingstore','on',...
            'NextPlot','add', ...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));
    end
    figure(mif55)

    set(gca,'NextPlot','add')

    tt = newcat2;
    [ts,ti] = sort(tt(:,15));
    tt = tt(ti(:,1),:);

    for i = 1:length(tt)
        pt = plot(tt(i,1),tt(i,2),'o');
        set(gca,'NextPlot','add')
    end

end

function eq_dep(misfitAngle)
    %  earthquake_depth
    % August 95 by Zhong Lu

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();
    myFigName='Depth vs Earthquake Number';

    mif66=findobj('Type','Figure','-and','Name',myFigName);

    if isempty(mif66)
        mif66 = figure_w_normalized_uicontrolunits( ...
            'Name',myFigName,...
            'NumberTitle','off', ...
            'backingstore','on',...
            'NextPlot','add', ...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));
    end
    figure(mif66)
    set(gca,'NextPlot','add')

    x = [1:length(mmi)]';
    [ss,ssi]=sort(catalog.Depth);
    plot(x,ss,'go');

    grid on

    ylabel('Depth of Earthquake','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    xlabel('Earthquake Number','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    set(gca,'NextPlot','replace');
end

function mi_dep(misfitAngle)
    %  misfit_magnitude
    % August 95 by Zhong Lu

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    mif77=findobj('Type','Figure','-and','Name','Misfit as a Function of Depth');

    if isempty(mif77)
        mif77 = figure_w_normalized_uicontrolunits( ...
            'Name','Misfit as a Function of Depth',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'NextPlot','add', ...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        set(gca,'NextPlot','add')

    end
    figure_w_normalized_uicontrolunits(mif77)
    set(gca,'NextPlot','add')

    plot(catalog.Depth,misfitAngle,'go');

    grid
    %set(gca,'box','on',...
    %        'SortMethod','childorder','TickDir','out','FontWeight',...
    %        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2);

    xlabel('Depth of Earthquake','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    ylabel('Misfit Angle ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    set(gca,'NextPlot','replace');
end

function mi_ma(misfitAngle)
    %  misfit_magnitude
    % August 95 by Zhong Lu

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    mif88=findobj('Type','Figure','-and','Name','Misfit as a Function of Magnitude');

    if isempty(mif88)
        mif88 = figure_w_normalized_uicontrolunits( ...
            'Name','Misfit as a Function of Magnitude',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'NextPlot','add', ...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        set(gca,'NextPlot','add')

    end
    figure(mif88)
    set(gca,'NextPlot','add')

    plot(catalog.Magnitude,misfitAngle,'go');

    grid
    %set(gca,'box','on',...
    %        'SortMethod','childorder','TickDir','out','FontWeight',...
    %        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2);

    xlabel('Magnitude of Earthquake','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    ylabel('Misfit Angle ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    set(gca,'NextPlot','replace');
end

function ast_misfit()
    %  ast_misfit calculates A as(t) value for a cumulative number curve and displayed in the plot.
    %

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    dat(:,2)=mi2(:,2);
    dat(:,1)=[1:length(mi2(:,1))]';
    cumu = dat(:,2);
    xt = dat(:,1);
    cumu2 = cumsum(cumu);

    %  winlen_days is the cutoff at the beginning and end of the analyses
    %  to avoid spikes at the end
    winlen_days = 5;

    %
    % calculate mean and z value
    ncu = length(xt);
    as = zeros(1,ncu);

    t0b = dat(1,1);
    n = length(dat(:,1));
    teb = dat(n,1);
    tdiff = ncu;

    for i = winlen_days+1:tdiff-winlen_days
        mean1 = mean(cumu(1:i));
        mean2 = mean(cumu(i+1:ncu));
        var1 = cov(cumu(1:i));
        var2 = cov(cumu(i+1:ncu));
        as(i) = (mean1 - mean2)/(sqrt(var1/i+var2/(tdiff-i)));
    end     % for i

    %  Plot the as(t)
    %clf
    figure;
    orient landscape
    % orient tall
    rect = [0.1,  0.10, 0.8, 0.7];
    axes('position',rect);
    yyaxis('left')
    plot(xt,as);
    yyaxis('right')
    plot(xt,cumu2);
    xlabel('Event');
    ylabel('z-value');
    grid

    set(gca,'NextPlot','add');

    %  show option from here
    %
    uicontrol('Units','normal','Position',[.9 .86 .10 .05],'String','Close', 'callback',@(~,~)close())

    str2 = 'AS of Earthquake Number';
    title(str2);
end

function compare_misfit()
    % Compare is used to compare the significance of two segments
    % in the plot of cumulative misfit as a function of earthquake number.
    %  --- Zhong Lu, June 1994.
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    dat(:,2)=mi2(:,2);
    dat(:,1)=[1:length(mi2(:,1))]';
    cumu = dat(:,2);
    xt = dat(:,1);
    cumu2 = cumsum(cumu);

    set(gca,'NextPlot','add');
    ZG.bin_dur = days(0.5);
    par2 = 1.0;
    choice = input('type 1 to select range with cursor, 2 to input event numbers  ');
    if choice == 1
        t1 = [];
        t1 = ginput(1);
        t1(1)=round(t1(1));
        t1p = [  t1 ; t1(1) t1(2)-ZG.bin_dur];
        plot(t1p(:,1),t1p(:,2),'r');
        text( t1(1),t1(2)-par2,['t1: ', num2str(t1p(1))] );

        t2 = [];
        t2 = ginput(1);
        t2(1)=round(t2(1));
        t2p = [  t2 ; t2(1) t2(2)-ZG.bin_dur];
        plot(t2p(:,1),t2p(:,2),'r');
        text( t2(1),t2(2)-par2,['t2: ', num2str(t2p(1))] );

        t3 = [];
        t3 = ginput(1);
        t3(1)=round(t3(1));
        t3p = [  t3 ; t3(1) t3(2)+ZG.bin_dur];
        plot(t3p(:,1),t3p(:,2),'r');
        text( t3(1),t3(2)+par2,['t3: ', num2str(t3p(1))] );

        t4 = [];
        t4 = ginput(1);
        t4(1)=round(t4(1));
        t4p = [  t4 ; t4(1) t4(2)+ZG.bin_dur];
        plot(t4p(:,1),t4p(:,2),'r');
        text( t4(1),t4(2)+par2,['t4: ', num2str(t4p(1))] );
    else
        %tmp = 't1(1),t2(1),t3(1),t4(1)';
        t1(1) = str2double(input('type the 1st event number, then return    ','s'));
        t2(1) = str2double(input('type the 2nd event number, then return    ','s'));
        t3(1) = str2double(input('type the 3rd event number, then return    ','s'));
        t4(1) = str2double(input('type the last event number, then return    ','s'));
    end  % if
    set(gca,'NextPlot','add');

    mean1 = mean(cumu(t1(1):t2(1)));
    mean2 = mean(cumu(t3(1):t4(1)));
    var1  = cov(cumu(t1(1):t2(1)));
    var2  = cov(cumu(t3(1):t4(1)));
    zvalue = (mean1 - mean2)/(sqrt(var1/(t2(1)-t1(1)+1)+var2/(t4(1)-t3(1)+1)))

    if abs(zvalue) >= 2.58 %99%
        S = sprintf('Significant at 99%% ');
        disp(S);
    elseif abs(zvalue) >= 1.96 %95%
        S = sprintf('Significant at 95%% ');
        disp(S);
    elseif abs(zvalue) >= 1.64 %90%
        S = sprintf('Significant at 90%% ');
        disp(S);
    elseif abs(zvalue) >= 1.44 %85%
        S = sprintf('Significant at 85%% ');
        disp(S);
    else
        S = sprintf('May Significant below 85%% ');
        disp(S);
    end % if

    % use the t-test
    tvalue=(mean1 - mean2) * sqrt(t2(1)-t1(1)+t4(1)-t3(1)) / sqrt((t2(1)-t1(1)) * var1+(t4(1)-t3(1))*var2) / sqrt(1.0/(t2(1)-t1(1)+1)+1.0/(t4(1)-t3(1)+1))

    N=t2(1)-t1(1)+t4(1)-t3(1)
    disp('N=n1+n2-2');

end

function plotmima(var1, mi)
    report_this_filefun();

    ZG=ZmapGlobal.Data;
    global mif1
    global oneOfHowManyPopupIdx

    sc = oneOfHowManyPopupIdx;
    angMisfit = mi(:,2)+1; % added 1 because it's used as sizes
    figure(mif1) %TODO figure out where mif1 comes from
    delete(findobj(mif1,'Type','axes'));
    rect = [0.15,  0.20, 0.75, 0.65];
    axes('position',rect)
    watchon

    if var1 == 1

        for i = 1:catalog.Count
            pl =  plot(catalog.Longitude(i),catalog.Latitude(i),'ro');
            set(gca,'NextPlot','add')
            set(pl,'MarkerSize',angMisfit(i)/sc)
        end

    elseif var1 == 2

        for i = 1:catalog.Count
            pl =  plot(catalog.Longitude(i),catalog.Latitude(i),'bx');
            set(gca,'NextPlot','add')
            set(pl,'MarkerSize',angMisfit(i)/sc,'LineWidth',angMisfit(i)/sc)
        end

    elseif var1 == 3

        for i = 1:catalog.Count
            pl =  plot(catalog.Longitude(i),catalog.Latitude(i),'bx');
            set(gca,'NextPlot','add')
            c = angMisfit(i)/max(angMisfit);
            set(pl,'MarkerSize',angMisfit(i)/sc,'LineWidth',angMisfit(i)/sc,'Color',[ c c c ] )
        end

    elseif var1 == 4
        pl =  plot(catalog.Longitude,catalog.Latitude,'bx');
    end

    set(gca,'NextPlot','add')
    %axis([ s2_west s1_east s4_south s3_north])
    %zmap_update_displays();

    xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    strib = 'Misfit Map ';
    title(strib,'FontWeight','bold',...
        'FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')

    set(gca,'Color',color_bg);
    set(gca,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    watchoff
end

function plotmimac(mi,inde)

    % TODO maybe move into domisfit

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    global mif2 mif1
    global oneOfHowManyPopupIdx

    %var1 = 4;
    sc = oneOfHowManyPopupIdx;
    figure(UNK) % FIXME: really? this figure? unsure
    delete(findobj(UNK,'Type','axes'));
    rect = [0.15,  0.20, 0.75, 0.65];
    axes('position',rect)
    watchon

    % check if cross-section exists
    figNumber=findobj('Type','Figure','-and','Name','Cross -Section');

    if isempty(figNumber)
        errordlg('Please create a cross-section first, then rerun the last selection');
        nlammap
        return
    end

    % check if cross-section is still current
    if max(mi(:,1)) > length(mi(:,1))
        errordlg('Please rerun the cross-section first, then rerun the last selection');
        nlammap
        return
    end

    mic = mi(inde,:);
    le = size(newa,2); %FIXME where does newa come from? ZG.newa? input parameter?  Needs to be treated like a ZmapCatalog

    if var1 == 1
        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'ro');
            set(gca,'NextPlot','add')
            set(pl,'MarkerSize',mic(i,2)/sc)
        end

    elseif var1 == 2

        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'bx');
            set(gca,'NextPlot','add')
            set(pl,'MarkerSize',mic(i,2)/sc,'LineWidth',mic(i,2)/sc)
        end

    elseif var1 == 3

        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'bx');
            set(gca,'NextPlot','add')
            c = mic(i,2)/max(mic(:,2));
            %c = newa(i,15)*10;
            set(pl,'MarkerSize',mic(i,2)/sc+3,'LineWidth',mic(i,2)/sc+0.5,'Color',[ c c c ] )
        end

    elseif var1 == 4

        g = jet;
        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'bx');
            set(gca,'NextPlot','add')
            c = floor(mic(i,2)/max(mic(:,2))*63+1);
            set(pl,'MarkerSize',4,'LineWidth',2,'Color',[ g(c,:) ] )
        end
        colorbar
        colormap(jet)
    end

    if exist('maex', 'var')
        set(gca,'NextPlot','add')
        pl = plot(maex,-maey,'*m');
        set(pl,'MarkerSize',8,'LineWidth',2)
    end

    if exist('maex', 'var')
        set(gca,'NextPlot','add')
        pl = plot(maex,-maey,'*m');
        set(pl,'MarkerSize',8,'LineWidth',2)
    end

    xlabel('Distance [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    ylabel('Depth [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    strib = [  'Misfit '];
    title(strib,'FontWeight','bold',...
        'FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')

    set(gca,'Color',color_bg);
    set(gca,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position',[0.9 0.7 0.08 0.08],...
        'String','Grid',...
        'callback',@(~,~)mificrgr(mi,inde));

    uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position',[0.9 0.6 0.08 0.08],...
        'String','Sel EQ',...
        'callback',@cb_pickinv);

    watchoff

    function cb_pickinv(~,~)
        newa2=crosssel(newa);
        ZG.newt2=newa2;
        ZG.newcat=newa2;
        ctp=CumTimePlot(ZG.newt2);
        ctp.plot();
    end
end

function newcat2=plotmi(var1, newcat2, mi)
    %plot misfit (?)
    % TODO make this work with the new catalogs
    report_this_filefun();

    global  mif2 mif1
    global tmp % REALLY? global tmp?  "tmp" is 1:nEvents
    % cumu2 mi2
    figNumber=findobj('Type','Figure','-and','Name','Misfit ');
    figure(figNumber);
    delete(findobj(figNumber,'Type','axes'));
    rect = [0.15,  0.15, 0.75, 0.65];
    axes('position',rect)
    ax=gca;
    nEvents=newcat2.Count;
    tmp=1:nEvents;
    sixSlices=round(0 : nEvents/5 : nEvents);
    sixSlices(1)=1;

    var2=var1;

    misfitAngle = mi(:,2);
    X = 1:nEvents;
    xtitle=sprintf('Number of Eqs (sorted by %s)',lower(var1));
    switch (var1)
        case {'Longitude','Latitude','Magnitude','Depth'}
            % plot_by_lon(); %by lon
            plot_by_field(var1);
        case 'Time'
            plot_by_time(); % by date
        case 'Strike'
            plot_by_strike(); % along strike
        case 'Default'
            option_7(); %unsorted
        otherwise
            error('unknown choice for plotmi');
    end

    grid('on')
    set(ax,'Color',color_bg);
    set(ax,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    ylabel('Cumulative Misfit ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    xlabel(xtitle,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

    if var1 ~= "Default"
        set(gca,'NextPlot','add')
        for i=1:6
            plot(ax,tmp(sixSlices(i)),cumu2(sixSlices(i)),'xr');
            str=['  ',num2str(newcat2(sixSlices(i),var2))];
            te=text(tmp(sixSlices(i)),cumu2(sixSlices(i)),str);
            set(te,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
        end
    else

    end

    function plot_by_field(name)
        % assumes that misfit matrix (mi) has same number of rows as
        % number of earthquakes in catalog
        [~,is] = sort(newcat2.(name));
        newcat2.sort(name); % sort the catalog itself by this field

        cumu2=cumsum(misfitAngle(is));
        plot(1:nEvents , cumu2 , 'o');
        xtitle=sprintf('Number of Eqs (sorted by %s)',lower(name));
    end

    function plot_by_time()
        [~,is] = sort(newcat2.Date);
        newcat2.sort(Date);
        cumu2=cumsum(misfitAngle(is));
        pl = plot(tmp,cumu2,'o');
        xtitle='Number of Eqs (sorted by time)';
    end

    function plot_by_strike()
        % [~,is] = sort(newcat2(:,15));
        [~,is] = sort(newcat2(:,end));
        newa2 = newcat2.subset(is) ;
        cumu2=cumsum(misfitAngle(is));
        pl = plot(newa2(:,16)-18.6,cumu2,'o');
        xtitle='Number of Eqs (sorted along strike)';
        var2=15;
    end

    function option_7()
        mi2 = mi ;
        cumu2=cumsum(mi2(:,2));
        pl = plot(tmp,cumu2,'o');
        xtitle='Number of Eqs ';
    end
end

function mifigrid(var1,mi)

    % This function creates a grid with spacing dx, dy (in degrees)
    % The size is selected interactively in an input window.
    % The relative quiescence will be calculated for every grid point
    % for a specific time and plotted in a Seismolap-Quiescence map
    %
    % turned into function by Celso G Reyes 2017
    %
    %MIFIGRID(var1,mi)
    %
    % Alexander Allmann

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    global  ni me1 va1
    global a h1 map dx dy Mmin lap1 seismap
    global normlap1 normlap2 mif1 mifmap

    if var1==1

        %input window
        %
        %default parameters
        dx= .5;                      %grid spacing east-west
        dy= .5;                      %grid spacing north-south
        % ldx=100;                     %side length of interaction zone in km (for seislap)
        % tlap=300;                    %interaction time in days (for seislap)
        Mmin=3;                      %minimum magnitude

        %create a input window
        figure_w_normalized_uicontrolunits(...
            'Name','Grid Input Parameter',...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'units','points',...
            'Visible','off', ...
            'Position',[ ZG.welcome_pos + [200, -200], 450, 250]);
        axis off

        %create a dialog box for the input
        freq_field1=uicontrol('Style','edit',...
            'Position',[.60 .36 .15 .08],...
            'Units','normalized','String',num2str(dx),...
            'callback',@callbackfun_001);

        freq_field2=uicontrol('Style','edit',...
            'Position',[.60 .27 .15 .08],...
            'Units','normalized','String',num2str(dy),...
            'callback',@callbackfun_002);

        freq_field3=uicontrol('Style','edit',...
            'Position',[.60 .48 .15 .08],...
            'Units','normalized','String',num2str(ni),...
            'callback',@callbackfun_003);

        close_button=uicontrol('Style','Pushbutton',...
            'Position',[.70 .05 .15 .12 ],...
            'Units','normalized', 'Callback', @callbackfun_004,'String','Cancel');

        go_button1=uicontrol('Style','Pushbutton',...
            'Position',[.20 .05 .15 .12 ],...
            'Units','normalized',...
            'callback',@callbackfun_005,...
            'String','Go');

        txt4 = text(...
            'Position',[0.50 0.74 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.l ,...
            'FontWeight','bold',...
            'String',' Grid Parameter');
        txt5 = text(...
            'Position',[0. 0.35 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','Spacing in x (dx) in deg:');

        txt6 = text(...
            'Position',[0. 0.25 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String','Spacing in y (dy) in deg:');

        txt2 = text(...
            'Position',[0. 0.5 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            'FontWeight','bold',...
            'String',' # of EQ Ni:');

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

    elseif var1==2           %area selection

        figure(map);
        set(gca,'NextPlot','add')
        ax=findobj(gcf,'Tag','mainmap_ax');
        [x,y, mouse_points_overlay] = select_polygon(ax);

        %figure_w_normalized_uicontrolunits(mif1)

        plos2 = plot(x,y,'b-');        % plot outline
        sum3 = 0.;
        pause(0.3)

        %create a rectangular grid
        xvect=[min(x):dx:max(x)];
        yvect=[min(y):dy:max(y)];
        tmpgri=zeros((length(xvect)*length(yvect)),2);
        n=0;
        for i=1:length(xvect)
            for j=1:length(yvect)
                n=n+1;
                tmpgri(n,:)=[xvect(i) yvect(j)];
            end
        end
        %extract all gridpoints in chosen polygon
        XI=tmpgri(:,1);
        YI=tmpgri(:,2);
        ll = polygon_filter(x,y, XI, YI, 'inside');
        %grid points in polygon
        newgri=tmpgri(ll,:);

        % Plot all grid points
        gcf
        plot(newgri(:,1),newgri(:,2),'+k')
        drawnow

        if length(xvect) < 2  ||  length(yvect) < 2
            errordlg('Selection too small! (not a matrix)');
            return
        end

        %calculate lap1(relative quiescence) at every grid point
        %
        ZG.newcat=a;                   %ZG.newcat is only a local variable
        bcat=ZG.newcat;

        me1=zeros(length(newgri(:,1)),1);
        va1=zeros(length(newgri(:,1)),1);

        wai = waitbar(0,' Please Wait ...  ');
        set(wai,'NumberTitle','off','Name','Makegrid - Percent completed');
        drawnow

        for i= 1:length(me1)   %all eqs which are in spacewindow in east-west direction
            l = sqrt(((ZG.newcat.Longitude-newgri(i,1))*cosd(newgri(i,2))*111).^2 +...
                ((ZG.newcat.Latitude-newgri(i,2))*111).^2) ;
            [s,is] = sort(l);
            b = ZG.newcat.subset(is) ;       % re-orders matrix to agree row-wise
            mi2 = mi(is(:,1),2);    % take first ni points
            mi2 = mi2(1:ni);
            me1(i) = mean(mi2);
            va1(i) = std(mi2);
            if rem(i,20)==0;  waitbar(i/length(me1));end

        end

        close(wai)
        %make a color map
        % Find out if figure already exists
        %
        mifmap=findobj('Type','Figure','-and','Name','Misfit-Map 2');

        % Set up the Seismicity Map window Enviroment
        %
        if isempty(mifmap)
            mifmap = figure_w_normalized_uicontrolunits( ...
                'Name','Misfit-Map 2',...
                'NumberTitle','off', ...
                'NextPlot','replace', ...
                'backingstore','on',...
                'Visible','off', ...
                'Position',[ 600 400 500 650]);
            % make menu bar

            set(gca,'NextPlot','add')
        end
        figure(mifmap)
        delete(findobj(mifmap,'Type','axes'));

        set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on','SortMethod','childorder')

        %minimum and maximum of normlap2 for automatic scaling
        ZG.maxc = max(normlap2);
        ZG.minc = min(normlap2);

        %construct a matrix for the color plot
        normlap1=ones(length(tmpgri(:,1)),1);
        normlap2=nan(length(tmpgri(:,1)),1)
        normlap3=nan(length(tmpgri(:,1)),1)
        normlap1(ll)=me1;
        normlap2(ll)=normlap1(ll);
        normlap1(ll)=va1;
        normlap3(ll)=normlap1(ll);

        normlap2=reshape(normlap2,length(yvect),length(xvect));
        normlap3=reshape(normlap3,length(yvect),length(xvect));

        %plot color image
        orient tall
        memifig2

        return

        rect = [0.25,  0.60, 0.7, 0.35];
        axes('position',rect)
        set(gca,'NextPlot','add')
        pco1 = pcolor(xvect,yvect,normlap2);
        shading interp
        colormap(jet)
        axis([ s2_west s1_east s4_south s3_north])
        set(gca,'NextPlot','add')
        colorbar
        overlay
        title('Mean of the Misfit','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

        set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on','TickDir','out')

        rect = [0.25,  0.10, 0.7, 0.35];
        axes('position',rect)
        set(gca,'NextPlot','add')
        pco1 = pcolor(xvect,yvect,normlap3);
        axis([ s2_west s1_east s4_south s3_north])
        set(gca,'NextPlot','add')
        shading interp
        colormap(jet)
        set(gca,'NextPlot','add')
        colorbar
        title(' Variance of the Misfit','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

        set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on','TickDir','out')

        overlay
        memifig2
    end

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dx=str2double(freq_field1.String);
        freq_field1.String=num2str(dx);
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dy=str2double(freq_field2.String);
        freq_field2.String=num2str(dy);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ni=str2double(freq_field3.String);
        freq_field3.String=num2str(ni);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;
        var1 = 2;
        mifigrid(var1, mi);
    end

end

function memifig2()
    % Misfitmap 2 (?)

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    %input window
    %
    %default parameters

    %make a color map
    % Find out if figure already exists
    %
    mifmap=findobj('Type','Figure','-and','Name','Misfit-Map 2');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(mifmap)
        mifmap = figure_w_normalized_uicontrolunits( ...
            'Name','Misfit-Map 2',...
            'NumberTitle','off', ...
            'NextPlot','replace', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',[ 600 400 500 350]);
        % make menu bar
        set(gca,'NextPlot','add')
    end

    figure(mifmap);
    delete(findobj(mifmap,'Type','axes'));

    set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','SortMethod','childorder')

    %minimum and maximum of normlap2 for automatic scaling
    ZG.maxc = max(normlap2);
    ZG.minc = min(normlap2);

    %construct a matrix for the color plot
    normlap1=ones(length(tmpgri(:,1)),1);
    normlap2=nan(length(tmpgri(:,1)),1)
    normlap3=nan(length(tmpgri(:,1)),1)
    normlap1(ll)=me1;
    normlap2(ll)=normlap1(ll);
    normlap1(ll)=va1;
    normlap3(ll)=normlap1(ll);

    normlap2=reshape(normlap2,length(yvect),length(xvect));
    normlap3=reshape(normlap3,length(yvect),length(xvect));

    %plot color image
    orient tall
    gx = xvect; gy = yvect;

    set(gca,'NextPlot','add')
    pco1 = pcolor(xvect,yvect,normlap2);
    shading interp
    colormap(flipud(jet(10)));
    axis([ min(gx) max(gx) min(gy) max(gy)])
    axis image

    set(gca,'NextPlot','add')
    h5 = colorbar('vert');
    set(h5,'Pos',[0.82 0.46 0.02 0.20],...
        'FontSize',12)

    if exist('maex', 'var')
        set(gca,'NextPlot','add')
        pl = plot(maex,-maey,'*k');
        set(pl,'MarkerSize',6,'LineWidth',2)
    end

    overlay
    title(['Mean of the Misfit (' num2str(sig) '/' num2str(az)  '/' num2str(plu) '/' num2str(phi) '/' num2str(R) ')'] ,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    xlabel('Longitude in [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('latitude in [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

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

function mificrgr( mi, inde)
    % This function creates a grid with spacing dx, dy (in degrees)
    % The size is selected interactively in an input window.
    % The relative quiescence will be calculated for every grid point
    % for a specific time and plotted in a Seismolap-Quiescence map
    %  Alexander Allmann
    %
    % turned into function by Celso G Reyes 2017

    %global freq_field1 freq_field2 freq_field3 freq_field4 freq_field5
    %global freq_field6 ni me1 va1
    %global h1 map dx dy Mmin lap1 seismap
    %global normlap1 normlap2 mif1 mifmap

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    %input window
    %
    %default parameters
    dx= .5;                      %grid spacing east-west
    dy= .5;                      %grid spacing north-south
    %ldx=100;                     %side length of interaction zone in km (for seislap)
    %tlap=300;                    %interaction time in days (for seislap)
    Mmin=3;                      %minimum magnitude
    ni=100; % number of events

    %create a input window
    fig=figure_w_normalized_uicontrolunits(...
        'Name','Grid Input Parameter',...
        'NumberTitle','off', ...
        'NextPlot','new', ...
        'units','points',...
        'Visible','off', ...
        'MenuBar','none',...
        'Position',[ ZG.welcome_pos + [200, -200], 700, 250]);
    axis off

    gridOpts = ZG.gridopt;
    esp = EventSelectionParameters('NumClosestEvents',ni);
    selOpts = EventSelectionChoice(fig,'evsel',esp);

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.70 .05 .15 .12 ],...
        'Units','normalized', 'Callback', @callbackfun_cancel,'String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'callback',@callbackfun_go,...
        'String','Go');

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

    function my_calculate()

        xsec_h=xsec_fig();

        if isempty(xsec_h) % best guess at figure handle variable
            nlammap
        else
            figure(xsec_h);
        end

        set(gca,'NextPlot','add')
        ax=findobj(gcf,'Tag','mainmap_ax');
        [x,y, mouse_points_overlay] = select_polygon(ax);

        figure(xsec_h)

        plos2 = plot(x,y,'b-');        % plot outline
        sum3 = 0.;
        pause(0.3)

        %create a rectangular grid
        xvect=[min(x):dx:max(x)];
        yvect=[min(y):dy:max(y)];
        tmpgri=zeros((length(xvect)*length(yvect)),2);
        n=0;
        for i=1:length(xvect)
            for j=1:length(yvect)
                n=n+1;
                tmpgri(n,:)=[xvect(i) yvect(j)];
            end
        end
        %extract all gridpoints in chosen polygon
        XI=tmpgri(:,1);
        YI=tmpgri(:,2);

        ll = polygon_filter(x,y, XI, YI, 'inside');
        %grid points in polygon
        newgri=tmpgri(ll,:);

        % Plot all grid points
        gcf
        plot(newgri(:,1),newgri(:,2),'+k')
        drawnow

        if length(xvect) < 2  ||  length(yvect) < 2
            errordlg('Selection too small! (not a matrix)');
            return
        end

        %
        ZG.newcat=a;                   %ZG.newcat is only a local variable
        bcat=ZG.newcat;

        me1=zeros(length(newgri(:,1)),1);
        va1=zeros(length(newgri(:,1)),1);
        mic = mi(inde,:);

        wai = waitbar(0,' Please Wait ...  ');
        set(wai,'NumberTitle','off','Name','Makegrid - Percent completed');
        drawnow

        for i= 1:length(me1)   %all eqs which are in spacewindow in east-west direction
            x = newgri(i,1);y = newgri(i,2);

            l = sqrt(((xsecx' - x)).^2 + ((xsecy + y)).^2) ;
            [s,is] = sort(l);
            b = ZG.newcat.subset(is) ;       % re-orders matrix to agree row-wise
            mi2 = mic(is(:,1),2);    % take first ni points
            mi2 = mi2(1:ni);
            me1(i) = mean(mi2);
            va1(i) = std(mi2);
            if rem(i,20)==0;  waitbar(i/length(me1));end
        end

        close(wai)
        %make a color map
        % Find out if figure already exists
        %
        mifmap=findobj('Type','Figure','-and','Name','Misfit-Map 2');

        % Set up the Seismicity Map window Enviroment
        %
        if isempty(mifmap)
            mifmap = figure_w_normalized_uicontrolunits( ...
                'Name','Misfit-Map 2',...
                'NumberTitle','off', ...
                'NextPlot','replace', ...
                'backingstore','on',...
                'Visible','off', ...
                'Position',[ 600 400 500 650]);
            % make menu bar

            set(gca,'NextPlot','add')
        end

        figure(mifmap);
        delete(findobj(mifmap,'Type','axes'));

        set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on','SortMethod','childorder')

        %minimum and maximum of normlap2 for automatic scaling
        ZG.maxc = max(normlap2);
        ZG.minc = min(normlap2);

        %construct a matrix for the color plot
        normlap1=ones(length(tmpgri(:,1)),1);
        normlap2=nan(length(tmpgri(:,1)),1)
        normlap3=nan(length(tmpgri(:,1)),1)
        normlap1(ll)=me1;
        normlap2(ll)=normlap1(ll);
        normlap1(ll)=va1;
        normlap3(ll)=normlap1(ll);

        normlap2=reshape(normlap2,length(yvect),length(xvect));
        normlap3=reshape(normlap3,length(yvect),length(xvect));

        %plot color image
        orient tall
        gx = xvect;
        gy = yvect;
        memifig

        return

        rect = [0.25,  0.60, 0.7, 0.35];
        axes('position',rect)
        set(gca,'NextPlot','add')
        pco1 = pcolor(xvect,yvect,normlap2);
        shading interp
        colormap(jet)
        %axis([ s2_west s1_east s4_south s3_north])
        axis([ min(gx) max(gx) min(gy) max(gy)])
        axis image

        set(gca,'NextPlot','add')
        colorbar
        if exist('maex', 'var')
            set(gca,'NextPlot','add')
            pl = plot(maex,-maey,'*m');
            set(pl,'MarkerSize',8,'LineWidth',2)
        end

        %overlay
        title('Mean of the Misfit','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        xlabel('Distance in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

        set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on','TickDir','out')

        rect = [0.25,  0.10, 0.7, 0.35];
        axes('position',rect)
        set(gca,'NextPlot','add')
        pco1 = pcolor(xvect,yvect,normlap3);
        axis([ min(gx) max(gx) min(gy) max(gy)])
        axis image

        if exist('maex', 'var')
            set(gca,'NextPlot','add')
            pl = plot(maex,-maey,'*w');
            set(pl,'MarkerSize',8,'LineWidth',2)
        end

        set(gca,'NextPlot','add')
        shading interp
        colormap(jet)
        set(gca,'NextPlot','add')
        colorbar
        title(' Variance of the Misfit','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        xlabel('Distance in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

        set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
            'FontWeight','bold','LineWidth',1.5,...
            'Box','on','TickDir','out')

    end

    function callbackfun_cancel(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;

    end

    function callbackfun_go(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        dx=gridOpts.dx;
        dy=gridOpts.dy;
        ni=selOpts.ni;
        delete(selOpts);
        close;
        my_calculate();
    end

end

function memifig()
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    %input window
    %
    %default parameters

    %make a color map
    % Find out if figure already exists
    %
    mifmap=findobj('Type','Figure','-and','Name','Misfit-Map 2');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(mifmap)
        mifmap = figure_w_normalized_uicontrolunits( ...
            'Name','Misfit-Map 2',...
            'NumberTitle','off', ...
            'NextPlot','replace', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',[ 600 400 500 350]);
        % make menu bar

        set(gca,'NextPlot','add')
    end
    figure(mifmap);
    delete(findobj(mifmap,'Type','axes'));

    set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on','SortMethod','childorder')

    %minimum and maximum of normlap2 for automatic scaling
    ZG.maxc = max(normlap2);
    ZG.minc = min(normlap2);

    %construct a matrix for the color plot
    normlap1=ones(length(tmpgri(:,1)),1);
    normlap2=nan(length(tmpgri(:,1)),1)
    normlap3=nan(length(tmpgri(:,1)),1)
    normlap1(ll)=me1;
    normlap2(ll)=normlap1(ll);
    normlap1(ll)=va1;
    normlap3(ll)=normlap1(ll);

    normlap2=reshape(normlap2,length(yvect),length(xvect));
    normlap3=reshape(normlap3,length(yvect),length(xvect));

    %plot color image
    orient tall
    gx = xvect; gy = yvect;

    set(gca,'NextPlot','add')
    pco1 = pcolor(xvect,yvect,normlap2);
    shading interp
    colormap( flipud(jet(64)) );
    axis([ min(gx) max(gx) min(gy) max(gy)])
    axis image

    set(gca,'NextPlot','add')
    h5 = colorbar('vert');
    set(h5,'Pos',[0.82 0.46 0.03 0.10],...
        'FontSize',2)

    if exist('maex', 'var')
        set(gca,'NextPlot','add')
        pl = plot(maex,-maey,'*w');
        set(pl,'MarkerSize',6,'LineWidth',1)
    end

    %overlay
    title('Mean of the Misfit','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    xlabel('Distance in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)

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

author: CelsoReyes


text:
function msg = comparisonMessage(user_message, assertion_message, A, B)
%comparisonMessage Generate assertion message when comparing two arrays.
%   msg = comparisonMessage(user_message, assertion_message, A, B) returns a
%   string appropriate to use in a call to throw inside an assertion function
%   that compares two arrays A and B.
%
%   The string returned has the following form:
%
%       <user_message>
%       <assertion_message>
%
%       First input:
%       <string representation of value of A>
%
%       Second input:
%       <string representation of value of B>
%
%   user_message can be the empty string, '', in which case user_message is
%   skipped.

%   Steven L. Eddins
%   Copyright 2009 The MathWorks, Inc.

msg = sprintf('%s\n\n%s\n%s\n\n%s\n%s', ...
    assertion_message, ...
    'First input:', ...
    xunit.utils.arrayToString(A), ...
    'Second input:', ...
    xunit.utils.arrayToString(B));

if ~isempty(user_message)
    msg = sprintf('%s\n%s', user_message, msg);
end

author: ricedsp


text:
% polar angle of closest point on the wall
function [data,units] = compute_arena_angle(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = atan2(trx(fly).y_mm-trx.landmark_params{n}.arena_center_mm_y(n),...
    trx(fly).x_mm-trx.landmark_params{n}.arena_center_mm_x(n));

end
units = parseunits('rad');

author: kristinbranson

repo: JAABA


text:
function imdb = os_get_database(osDir)

load(fullfile(osDir, 'imdb/imdb.mat'), 'imdb') ;
imdb.imageDir512 = fullfile(osDir, imdb.imageDir512) ;
imdb.imageDir1024 = fullfile(osDir, imdb.imageDir1024) ;
imdb.maskDir512 = fullfile(osDir, imdb.maskDir512) ;
imdb.maskDir1024 = fullfile(osDir, imdb.maskDir1024) ;

% use these by default
imdb.imageDir = imdb.imageDir512 ;
imdb.maskDir = imdb.maskDir512 ;
imdb.segmDir = fullfile(osDir, 'segm/512') ;

% split images in train, val, test
n = numel(imdb.images.id) ;
m = round(n/3) ;
sets = [1 * ones(1,m), 2 * ones(1,m), 3 * ones(1,n-2*m)] ;
rng(0) ;
imdb.images.set = sets(randperm(n)) ;
[~,i] = ismember(imdb.segments.imageId, imdb.images.id) ;
imdb.segments.set = imdb.images.set(i) ;
imdb.segments.label = imdb.segments.materialClass ;

% now remove all the segments that belong to classes that are not in use
ok = logical(imdb.meta.inUse(imdb.segments.label)) ;
imdb.segments = soaSubsRef(imdb.segments, ok) ;

% finally, merge the background classes
% bkg = [18 25]
imdb.segments.label(imdb.segments.label == 25) = 18 ;
imdb.meta.inUse(25) = false ;
imdb.meta.classes{18} = 'other' ;
imdb.meta.classes{25} = 'other' ;

% no difficult regions by default
imdb.segments.difficult = false(1, numel(imdb.segments.id)) ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function result = evaluate(output_dir, target_dir, task)
  ts = get_task(output_dir, target_dir, task);
  [p, s, a] = run_eval_template(ts{1,2},ts{1,3},ts{1,4},get_path(ts{1,5}),ts{1,6},ts{1,7});
  result = [ts{1,1} ' psnr,ssim,abs= ' num2str(p) ', ' num2str(s) ', ' num2str(a)];
  result
end

repo: toflow


text:
% Wrapper Feature Selection Toolbox by Jingwei Too - 9/12/2020

function model = jfs(type,feat,label,opts)
switch type
  % 2020
  case 'mpa'  ; fun = @jMarinePredatorsAlgorithm;
  case 'gndo' ; fun = @jGeneralizedNormalDistributionOptimization;
  case 'sma'  ; fun = @jSlimeMouldAlgorithm;
  case 'eo'   ; fun = @jEquilibriumOptimizer;
  case 'mrfo' ; fun = @jMantaRayForagingOptimization;
  % 2019
  case 'aso'  ; fun = @jAtomSearchOptimization;
  case 'hho'  ; fun = @jHarrisHawksOptimization;
  case 'hgso' ; fun = @jHenryGasSolubilityOptimization;
  case 'pfa'  ; fun = @jPathFinderAlgorithm;
  case 'pro'  ; fun = @jPoorAndRichOptimization;
  % 2018
  case 'boa'  ; fun = @jButterflyOptimizationAlgorithm;
  case 'epo'  ; fun = @jEmperorPenguinOptimizer;
  case 'tga'  ; fun = @jTreeGrowthAlgorithm;
  % 2017
  case 'abo'  ; fun = @jArtificialButterflyOptimization;
  case 'ssa'  ; fun = @jSalpSwarmAlgorithm;
  case 'sbo'  ; fun = @jSatinBowerBirdOptimization;
  case 'wsa'  ; fun = @jWeightedSuperpositionAttraction;
  % 2016
  case 'ja'   ; fun = @jJayaAlgorithm;
  case 'csa'  ; fun = @jCrowSearchAlgorithm;
  case 'sca'  ; fun = @jSineCosineAlgorithm;
  case 'woa'  ; fun = @jWhaleOptimizationAlgorithm;
  % 2015
  case 'alo'  ; fun = @jAntLionOptimizer;
  case 'hlo'  ; fun = @jHumanLearningOptimization;
  case 'mbo'  ; fun = @jMonarchButterflyOptimization;
  case 'mfo'  ; fun = @jMothFlameOptimization;
  case 'mvo'  ; fun = @jMultiVerseOptimizer;
  case 'tsa'  ; fun = @jTreeSeedAlgorithm;
  % 2014
  case 'gwo'  ; fun = @jGreyWolfOptimizer;
  case 'sos'  ; fun = @jSymbioticOrganismsSearch;
  % 2012
  case 'fpa'  ; fun = @jFlowerPollinationAlgorithm;
  case 'foa'  ; fun = @jFruitFlyOptimizationAlgorithm;
  % 2009 - 2010
  case 'ba'   ; fun = @jBatAlgorithm;
  case 'fa'   ; fun = @jFireflyAlgorithm;
  case 'cs'   ; fun = @jCuckooSearchAlgorithm;
  case 'gsa'  ; fun = @jGravitationalSearchAlgorithm;
  % Traditional
  case 'abc'  ; fun = @jArtificialBeeColony;
  case 'hs'   ; fun = @jHarmonySearch;
  case 'de'   ; fun = @jDifferentialEvolution;
  case 'aco'  ; fun = @jAntColonyOptimization;
  case 'acs'  ; fun = @jAntColonySystem;
  case 'pso'  ; fun = @jParticleSwarmOptimization;
  case 'gat'  ; fun = @jGeneticAlgorithmTour;
  case 'ga'   ; fun = @jGeneticAlgorithm;
  case 'sa'   ; fun = @jSimulatedAnnealing;
end
tic;
model = fun(feat,label,opts);
% Computational time
t = toc;

model.t = t;
fprintf('\n Processing Time (s): %f % \n',t); fprintf('\n');
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
function L = pureShear(exp,comp,e)
% defines velocityGradientTensor representing pure shear
%
% Syntax
%   L = velocityGradientTensor.uniaxialCompression(exp,comp)
%
% Input
%  exp  - @vector3d expansion direction
%  comp - @vector3d compression direction
%  e - strain rate
%
% Output
%  L - @velocityGradientTensor
%

if nargin == 2, e = 1; end

v1 = normalize(exp - comp);
v2 = normalize(exp + comp);

L = 2*velocityGradientTensor(e .* (dyad(v1,v2) + dyad(v1,v2)'));

end

author: mtex-toolbox

repo: mtex

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


text:
function opts = getGlobalDefaults()
% Returns a structure with default option values for Manopt.
%
% function opts = getGlobalDefaults()
%
% Returns a structure opts containing the global default options such as
% verbosity level etc. Typically, global defaults are overwritten by solver
% defaults, which are in turn overwritten by user-specified options.
% See the online Manopt documentation for details on options.
%
% See also: mergeOptions

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   Aug. 2, 2018 (NB):
%       Changed default storedepth to 2 from 20, since solvers should now
%       use storedb.erase() to keep the cache lean.

    % There should be no reason to modify this file.
    % For better compatibility with future Manopt versions,
    % use the options structure of solvers.
    %
    % Really: don't modify it.

    % Verbosity level: 0 is no output at all. The higher the verbosity, the
    % more info is printed / displayed during solver execution.
    opts.verbosity = 3;

    % If debug is set to true, additional computations may be performed and
    % debugging information is outputed during solver execution.
    opts.debug = false;

    % Maximum number of store structures to store. If set to 0, caching
    % capabilities are not disabled, but the cache will be emptied at each
    % iteration of iterative solvers (more specifically: every time the
    % solver calls to purge the storedb).
    opts.storedepth = 2;

    % Maximum amount of time a solver may execute, in seconds.
    opts.maxtime = inf;

end

author: NicolasBoumal

repo: manopt


text:
function plotFFTSeriesAngle(view,scan)
%
% plotFFTSeriesAngle(view,[scan])
%
% Plots the FFT of the tSeries for the current scan, averaging across
% all pixels (in all slices) in the current ROI.

if ~exist('scan','var')
    scan = getCurScan(view);
end
scan = getCurScan(view);
nCycles = numCycles(view,scan);
maxCycles = round(numFrames(view,scan)/3); % number of frequencies to plot

% Compute the mean tSeries
if view.selectedROI
  ROIcoords = getCurROIcoords(view);
else
  myErrorDlg('No current ROI');
end
tSeries = meanTSeries(view,scan,ROIcoords);

% Calulate the FFT;
absFFT=2*abs(fft(tSeries)) / length(tSeries);
angleFFT=angle(fft(tSeries));

% Header
ROIname = view.ROIs(view.selectedROI).name;
headerStr = ['Mean tSeries, ROI ',ROIname,', scan ',num2str(scan)];
set(gcf,'Name',headerStr);

% plot it
x= [1:maxCycles];
y1 =[absFFT(2:maxCycles+1)];
y2 =[angleFFT(2:maxCycles+1)];

subplot(2,1,1);
plot(x(1:nCycles-1),y1(1:nCycles-1),'b','LineWidth',2)
hold on
plot(x(nCycles-1:nCycles+1),y1(nCycles-1:nCycles+1),'r','LineWidth',2)
plot(x(nCycles+1:maxCycles),y1(nCycles+1:maxCycles),'b','LineWidth',2)
plot(x,y1,'bo','LineWidth',2);
hold off

% Ticks
fontSize = 14;
xtick=nCycles:nCycles:(maxCycles+1);
set(gca,'xtick',xtick);
set(gca,'FontSize',fontSize)
xlabel('Cycles per scan','FontSize',fontSize)
ylabel('Percent modulation','FontSize',fontSize)
grid on

% plot the phase as well
subplot(2,1,2);

plot(x(1:nCycles-1),y2(1:nCycles-1),'b','LineWidth',2)
hold on
plot(x(nCycles-1:nCycles+1),y2(nCycles-1:nCycles+1),'r','LineWidth',2)
plot(x(nCycles+1:maxCycles),y2(nCycles+1:maxCycles),'b','LineWidth',2)
plot(x,y2,'bo','LineWidth',2);
hold off
% Ticks
fontSize = 14;
xtick=nCycles:nCycles:(maxCycles+1);
set(gca,'xtick',xtick);
set(gca,'FontSize',fontSize)
xlabel('Cycles per scan','FontSize',fontSize)
ylabel('Phase','FontSize',fontSize)
grid on

% Save the data in gca('UserData')
data.x = x(1:maxCycles);
data.y1  =  y1(1:maxCycles);
data.y2 = 	y2(1:maxCycles);

set(gca,'UserData',data);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Example_3_ECGC: Mesh the greater US East Coast and Gulf of Mexico region
% with a high resolution inset around New York

clearvars; clc;

addpath('..')
addpath(genpath('../utilities/'));
addpath(genpath('../datasets/'));
addpath(genpath('../m_map/'));

% WJP: 08/02/2019: Updated to demonstrate using non-box
% (arbitrary polygon's) bbox's in both outer and inner meshes

%% STEP 1: set mesh extents and set parameters for mesh.
%% The greater US East Coast and Gulf of Mexico region

bbox = [-71.6 42.7; -64 30; -80 24; -85 38; -71.6 42.7]; %polygon boubox
min_el    = 1e3;  		        % minimum resolution in meters.
max_el    = 50e3; 		        % maximum resolution in meters.
wl        = 30;                 % 60 elements resolve M2 wavelength.
dt        = 0;                  % Automatically set timestep based on nearshore res
grade     = 0.35;               % mesh grade in decimal percent.
R         = 3; 			        % Number of elements to resolve feature.

%% STEP 2: specify geographical datasets and process the geographical data
%% to be used later with other OceanMesh classes...
dem       = 'SRTM15+.nc';
coastline = 'GSHHS_f_L1';
gdat1 = geodata('shp',coastline,'dem',dem,'h0',min_el,...
                'bbox',bbox);

%% STEP 3: create an edge function class
fh1 = edgefx('geodata',gdat1,...
             'fs',R,'wl',wl,'max_el',max_el,...
             'dt',dt,'g',grade);

%% Repeat STEPS 1-3 for a high resolution domain for High Res New York Part
min_el    = 30;  		% minimum resolution in meters.
max_el    = 1e3; 		% maximum resolution in meters.
max_el_ns = 240;  		% maximum resolution nearshore.

coastline = 'PostSandyNCEI';
dem       = 'PostSandyNCEI.nc';

%polygon boubox
bbox2 = [-74.25 40.5; -73.75 40.55; -73.75 41; -74 41; -74.25 40.5];
gdat2 = geodata('shp',coastline,'dem',dem,'h0',min_el,'bbox',bbox2);

fh2 = edgefx('geodata',gdat2,'fs',R,'wl',wl,...
             'max_el',max_el,'max_el_ns',max_el_ns,...
             'dt',dt,'g',grade);

%% STEP 4: Pass your edgefx class object along with some meshing options
%% and build the mesh...
mshopts = meshgen('ef',{fh1 fh2},'bou',{gdat1 gdat2},...
                  'plot_on',1,'proj','lam');
mshopts = mshopts.build;

%% Plot and save the msh class object/write to fort.14
m = mshopts.grd; % get out the msh object
m = interp(m,{gdat1 gdat2},'mindepth',1); % interpolate bathy to the mesh with minimum depth of 1 m
m = make_bc(m,'auto',gdat1);               % make the nodestring boundary conditions
plot(m,'type','bd');  % plot mesh on native projection with boundary conditions
plot(m,'type','b');   % plot bathy on native projection
save('ECGC_w_NYHR.mat','m'); write(m,'ECGC_w_NYHR');


text:
function obj = t2_pnp_func_new(in1,in2,in3)
coefs_tq1_1 = in2(1);
coefs_tq1_2 = in2(4);
coefs_tq2_1 = in2(2);
coefs_tq1_3 = in2(7);
coefs_tq2_2 = in2(5);
coefs_tq3_1 = in2(3);
coefs_tq1_4 = in2(10);
coefs_tq2_3 = in2(8);
coefs_tq3_2 = in2(6);
coefs_tq1_5 = in2(13);
coefs_tq2_4 = in2(11);
coefs_tq3_3 = in2(9);
coefs_tq1_6 = in2(16);
coefs_tq2_5 = in2(14);
coefs_tq3_4 = in2(12);
coefs_tq1_7 = in2(19);
coefs_tq2_6 = in2(17);
coefs_tq3_5 = in2(15);
coefs_tq1_8 = in2(22);
coefs_tq2_7 = in2(20);
coefs_tq3_6 = in2(18);
coefs_tq1_9 = in2(25);
coefs_tq2_8 = in2(23);
coefs_tq3_7 = in2(21);
coefs_tq2_9 = in2(26);
coefs_tq3_8 = in2(24);
coefs_tq3_9 = in2(27);
coefs_tq1_10 = in2(28);
coefs_tq2_10 = in2(29);
coefs_tq3_10 = in2(30);
pinvG2_1 = in1(2);
pinvG2_2 = in1(5);
pinvG2_3 = in1(8);
q0 = in3(1,:);
q1 = in3(2,:);
q2 = in3(3,:);
q3 = in3(4,:);
obj = q0.^2.*(coefs_tq1_1.*pinvG2_1+coefs_tq2_1.*pinvG2_2+coefs_tq3_1.*pinvG2_3)+q1.^2.*(coefs_tq1_5.*pinvG2_1+coefs_tq2_5.*pinvG2_2+coefs_tq3_5.*pinvG2_3)+q2.^2.*(coefs_tq1_8.*pinvG2_1+coefs_tq2_8.*pinvG2_2+coefs_tq3_8.*pinvG2_3)+q3.^2.*(coefs_tq1_10.*pinvG2_1+coefs_tq2_10.*pinvG2_2+coefs_tq3_10.*pinvG2_3)+q0.*q1.*(coefs_tq1_2.*pinvG2_1+coefs_tq2_2.*pinvG2_2+coefs_tq3_2.*pinvG2_3)+q0.*q2.*(coefs_tq1_3.*pinvG2_1+coefs_tq2_3.*pinvG2_2+coefs_tq3_3.*pinvG2_3)+q0.*q3.*(coefs_tq1_4.*pinvG2_1+coefs_tq2_4.*pinvG2_2+coefs_tq3_4.*pinvG2_3)+q1.*q2.*(coefs_tq1_6.*pinvG2_1+coefs_tq2_6.*pinvG2_2+coefs_tq3_6.*pinvG2_3)+q1.*q3.*(coefs_tq1_7.*pinvG2_1+coefs_tq2_7.*pinvG2_2+coefs_tq3_7.*pinvG2_3)+q2.*q3.*(coefs_tq1_9.*pinvG2_1+coefs_tq2_9.*pinvG2_2+coefs_tq3_9.*pinvG2_3);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function varargout = harris(varargin)
% VL_HARRIS  Harris corner strength
%  H = VL_HARRIS(I,SI) computes the Harris corner strength of the image I
%  at ``integration'' scale SI.
%
%  The Harris strength [1] of a pixel is a statistic of the gradient
%  of the image integrated in a neighborhood of that pixel. This
%  neighborhood is a Gaussian window of variance SI.
%
%  In computing the Harris corner strength, there is a second scale
%  parameter, the ``derivation'' scale SD, which is the variance of
%  the Gaussian kernel used to pre-smooth the image I before computing
%  its gradient. SI and SD are independent parameters and VL_HARRIS(I,SI)
%  assumes that I is already smoothed at level SD.
%
%  VL_HARRIS(I,SI) uses Noble's variation [2] of the Harris score. If
%  SIGMAP and SIGMAM are respectively the biggest and smallest
%  eigenvalue of the structure tensor at a pixel, the score is given
%  by (SIGMAP*SIGMAM) / (SIGMAP+SIGMAM/2).  Let GAMMA = SIGMAM/SIGMAP
%  the ratio between the eigenvalues, which measures the degree of
%  anisotropy of the tensor and is always comprised in the range
%  [0,1].  Noble's score can be decomposed in two factors: the biggest
%  eigenvalue SIGMAP and the number
%
%    RHO = (2 GAMMA) / (GAMMA + 1).
%
%  RHO is another measure of isotropy that has value one for a
%  symmetric tensor and and zero for maximally anisotropic tensor.
%  [H,DETAILS] = VL_HARRIS(I,SIGMA) returns the additional structure
%  DETAILS with the following fields:
%
%    DETAILS.SIGMAP
%    DETAILS.RHO
%
%  VL_HARRIS(I,SI,ALPHA) uses Harris' original score [1], defined to be
%  SIGMAP*SIGMAM - ALPHA*(SIGMAP+SIGMAM)^2. This can be decomposed in
%  the factors SIGMAP^2 (note the square) and
%
%    RHO = GAMMA - ALPHA (1+GAMMA)^2.
%
%  Note that RHO is equal to -ALPHA for a maximally anisotropic
%  tensor. Typically ALPHA=0.04 and this is what is used by
%  VL_HARRIS(I,SI,[]).
%
%  REMARK. The gradient of the image I, used to compute the structure
%  tensor, is computed using central differencies. This means that a
%  function line [+1,-1,+1,...] has null Harris' score. This is
%  generally assumed to be a sampling artifact, and might be
%  avoided by oversampling the image.
%
%  EXAMPLE::
%    To extacts Harris points from image I:
%      idx = vl_localmax( vl_harris( vl_imsmooth( I, sd ), si ) ) ;
%      [i,j] = ind2sub( size(I), idx )
%
%  REFERENCES::
%    [1] C. Harris and M. Stephens, "A combined corner and edge detector,"
%    in Proceedings of The Fourth Alvey Vision Conference, pp. 147-151,
%    1988.
%
%    [2] J. A. Noble, "Finding corners, "Image Vision Computing, vol. 6,
%    no. 2, pp. 121-128, 1988.
%
%  See also: VL_HELP().
[varargout{1:nargout}] = vl_harris(varargin{:});

author: yihui-he

repo: panorama

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


text:
%%*****************************************************************
%% HKMrhsfun: compute the right-hand side vector of the
%%            Schur complement equation for the  HKM direction.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [rhs,EinvRc,hRd] = HKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ)

    m = length(rp);
    if (nargin > 8)
       corrector = 1;
    else
       corrector = 0;
       hRd = zeros(m,1);
    end
    hEinvRc = zeros(m,1);
    EinvRc = cell(size(blk,1),1);
    rhsfree = [];
%%
    for p = 1:size(blk,1)
       pblk = blk(p,:);
       n = sum(pblk{2});
       if strcmp(pblk{1},'l')
	  if iscell(sigmu)
             EinvRc{p} = sigmu{p}./Z{p} -X{p};
	  else
             EinvRc{p} = sigmu./Z{p} -X{p};
          end
          Rq = sparse(n,1);
          if (corrector) & (norm(par.parbarrier{p})==0)
             Rq = dX{p}.*dZ{p}./Z{p};
          else
             tmp  = par.dd{p}.*Rd{p};
             tmp2 = mexMatvec(At{p},tmp,1);
             hRd = hRd + tmp2;
          end
	  EinvRc{p} = EinvRc{p} - Rq;
          tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
          hEinvRc = hEinvRc + tmp2;
       elseif strcmp(pblk{1},'q')
	  if iscell(sigmu)
             EinvRc{p} = qops(pblk,sigmu{p},par.Zinv{p},3) -X{p};
	  else
             EinvRc{p} = sigmu*par.Zinv{p} -X{p};
          end
          Rq = sparse(n,1);
          if (corrector) & (norm(par.parbarrier{p})==0)
  	     ff{p} = qops(pblk,1./par.gamz{p},Z{p},3);
             hdx = qops(pblk,par.gamz{p},ff{p},5,dX{p});
             hdz = qops(pblk,par.gamz{p},ff{p},6,dZ{p});
             hdxdz = Arrow(pblk,hdx,hdz);
             Rq = qops(pblk,par.gamz{p},ff{p},6,hdxdz);
          else
             tmp = par.dd{p}.*Rd{p} ...
                   + qops(pblk,qops(pblk,Rd{p},par.Zinv{p},1),X{p},3) ...
                   + qops(pblk,qops(pblk,Rd{p},X{p},1),par.Zinv{p},3);
             tmp2 = mexMatvec(At{p,1},tmp,1);
             hRd = hRd + tmp2;
          end
	  EinvRc{p} = EinvRc{p} - Rq;
          tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
          hEinvRc = hEinvRc + tmp2;
       elseif strcmp(pblk{1},'s')
	  if iscell(sigmu)
             %%ss = [0,cumsum(pblk{2})];
             %%sigmuvec = zeros(n,1);
             %%for k = 1:length(pblk{2});
             %%   sigmuvec(ss(k)+1:ss(k+1)) = sigmu{p}(k)*ones(pblk{2}(k),1);
             %%end
             sigmuvec = mexexpand(pblk{2},sigmu{p});
             EinvRc{p} = par.Zinv{p}*spdiags(sigmuvec,0,n,n) -X{p};
          else
             EinvRc{p} = sigmu*par.Zinv{p} -X{p};
          end
          Rq = sparse(n,n);
          if (corrector) & (norm(par.parbarrier{p})==0)
             Rq = Prod3(pblk,dX{p},dZ{p},par.Zinv{p},0);
	     Rq = 0.5*(Rq+Rq');
          else
             tmp = Prod3(pblk,X{p},Rd{p},par.Zinv{p},0,par.nzlistAy{p});
             tmp = 0.5*(tmp+tmp');
             tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),{tmp});
             hRd = hRd + tmp2;
          end
	  EinvRc{p} = EinvRc{p} - Rq;
          tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
          hEinvRc = hEinvRc + tmp2;
       elseif strcmp(pblk{1},'u')
          rhsfree = [rhsfree; Rd{p}];
       end
    end
%%
    rhs = rp + hRd - hEinvRc;
    rhs = full([rhs; rhsfree]);
%%*******************************************************************

author: yu-jiang

repo: radpbook

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


text:
function [u,p,info] =  tripreMaxwellsaddle(A,G,f,g,node,elem,bdFlag,M,grad,option)
%% tripreMaxwellsaddle solve the Maxwell system with divgence free condition,
%
%         [A  G] [u]  = [f]
%         [G' O] [p]  = [g]
%
% where  G = M_e*grad with the mass matrix for the edge element Me.
%
% This system can be rewritten as
%         [A+G*DMinv*G'       G] [u]  = f +G*DMinv*g0
%         [G'                 O] [p]  = g0
%
% where DMinv is the inverse of the diagonal matrix of the mass. Then
%
%  [A+G*DMinv*G'  G] [I    grad]                = [Abar   O  ]
%  [G'            O] [0  -Dminv*grad'*Me*grad ] = [G'     Ap]
%
%  We solve the system by GMRES with  the preconditioner
%
%       [I   grad                ]  [Abar O   ]^{-1}
%       [O  -Dminv*grad'*M_e*grad]  [G'   A_p ]
%
% where we compute the inverse Abar by mgHodgeLapE and Ap by mg.
%
% Created by Long chen and Jie Zhou on Aug,2015.

Nf = length(f);
Ng = length(g);
t = cputime;

%% Parameters
if ~exist('option','var'), option = []; end
Ndof = Nf + Ng;
option = mgoptions(option,Ndof);    % parameters
x0 = option.x0;
tol = option.tol;
maxIt = option.solvermaxit;
printlevel = option.printlevel;
d = size(node,2);

%% Set up auxiliary matrices
if d == 2
    area = simplexvolume(node,elem);
    Mvlump = accumarray([elem(:,1);elem(:,2);elem(:,3)],[area;area;area]/3,...
                        [max(elem(:)),1]);
elseif d == 3
    volume = abs(simplexvolume(node,elem)); % uniform refine in 3D is not orientation presereved
    Mvlump = accumarray([elem(:,1);elem(:,2);elem(:,3);elem(:,4)],...
                        [volume;volume;volume;volume]/4,[max(elem(:)),1]);
end
DMinv = spdiags(1./Mvlump(option.isFreeNode),0,Ng,Ng);
f = f + G*(DMinv*g);  % add second equation to the first one
Abar = A + G*DMinv*G'; % Hodge Laplacian
Ap = grad'*M*grad;    % scalar Laplacian

%% Set up matrices for multigrid
setupOption.solver = 'NO';
setupOption.isFreeEdge = option.isFreeEdge;
setupOption.freeDof = option.isFreeNode;
[~,~,Ai_N,Bi_N,BBi_N,Res_N,Pro_N] = mg(Ap,g,elem,setupOption);
[x,info,Ai,Bi,BBi,Res,Pro] = mgHodgeLapE(Abar,f,node,elem,bdFlag,setupOption); %#ok<*ASGLU>

%% Form a big matrix equation
bigA = [Abar G; G' sparse(Ng,Ng)];
bigF = [f; g];

%% Preconditioned GMRES
% options for V-cycle of Schur complement
if strcmp(option.solver,'CG') % change default set up in mgoption
   option.solver = 'Vcycle';
end
if isfield(option,'Vit')
   option.solvermaxIt = option.Vit;
else
   option.solvermaxIt = 1;
end
option.setupflag = false;
option.x0 = zeros(Nf,1);
option.printlevel = 0;
Apmgoption.printlevel = 0;
Apmgoption.printlevel = 0;
Apmgoption.setupflag = false;
Apmgoption.maxIt     = 1;
% minres for the saddle point system
[x,flag,stopErr,itStep,err] = gmres(bigA,bigF,20,tol,maxIt,@tripre,[],x0);
itStep = (itStep(1)-1)*20 + itStep(2); % total iteration
u = x(1:Nf);
p = x(Nf+1:end);

%% Output
time = cputime - t;
if printlevel >= 1
    fprintf('Triangular Preconditioned GMRES \n');
    fprintf('#dof: %8.0u,  #nnz: %8.0u, V-cycle: %2.0u, iter: %2.0u,   err = %8.2e,   time = %4.2g s\n',...
                 Ndof, nnz(bigA), option.solvermaxIt, itStep, stopErr, time)
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',stopErr);

%% Preconditioner
    function e  = tripre(r)
        r1 = r(1:Nf);
        r2 = r(Nf+1:end);
        e1 = mgHodgeLapE(Abar,r1,node,elem,bdFlag,option,Ai,Bi,BBi,Res,Pro);
        e2 = mg(Ap,r2-G'*e1,elem,Apmgoption,Ai_N,Bi_N,BBi_N,Res_N,Pro_N);
        e  = [e1 + grad*e2; -DMinv*(Ap*e2)];
    end
end

repo: ifem


text:
function [flatArray, index] = flattenCellMatrix(nestedArray)
% flatten nested (TRAFO) cell array
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

    if(nnz(size(nestedArray) ~= 1) > 1 || any(any(any(any(cellfun(@iscell,nestedArray))))))
        [flatArray,maxdim] = flattenCell(nestedArray);

        index = flatArray(2,:);
        index = cellfun(@(x) x(:,1:maxdim), index, 'UniformOutput', false);
        flatArray = flatArray(1,:);
    else
        index = [];
        flatArray = nestedArray(:).';
    end

end

function [flatArray,maxdim] = flattenCell(nestedArray,cnt)

% maximal 4D cell array as input

narginchk(1,2);
if(~exist('cnt','var'))
    cnt = 1;
end
if ~iscell(nestedArray),
    error('Must be a cell array.');
end
flatArray{2,1} = [];
for i=1:numel(nestedArray)
    if iscell(nestedArray{i})
        [y, maxdim] = flattenCell(nestedArray{i},cnt+1);
        [flatArray{:,end+1:end+size(y,2)}] = deal(y{:});
        for j=size(flatArray,2):-1:size(flatArray,2)-size(y,2)+1
            [row, col, slice, page] = ind2sub(size(nestedArray),i);
            flatArray{2,j}(cnt,:) = [row, col, slice, page];
            maxdim = max([maxdim, ndims(nestedArray)]);
        end
    else
        flatArray{1,end+1} = nestedArray{i};
        flatArray{2,end} = zeros(cnt,4);
        maxdim = ndims(nestedArray);
        [row, col, slice, page] = ind2sub(size(nestedArray),i);
        flatArray{2,end}(cnt,:) = [row, col, slice, page];
    end
end
flatArray(:,1) = [];

end

author: thomaskuestner