text:
function p03_default ( dim_num )

%*****************************************************************************80
%
%% P03_DEFAULT sets default values for problem 03.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p03_i4 ( 'D', '*', i4 );

  return
end


text:
function tests = test_imEuler1d(varargin)
%TEST_IMEULER1D  One-line description here, please.
%
%   output = test_imEuler1d(input)
%
%   Example
%   test_imEuler1d
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-10-08,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function testBasic(testCase)

img = [0 0 1 1 1 0 1 0];

len = imEuler1d(img);

exp = 2;
assertEqual(testCase, exp, len);

function testAddBorder(testCase)

img = ones(1, 5);
imgb = padarray(img, [0 1]);

len = imEuler1d(img);
len2 = imEuler1d(imgb);

assertEqual(testCase, len, len2);

function testCalib(testCase)

img = [0 0 1 1 1 0 1 0];

len = imEuler1d(img, .4);

exp = 2;
assertEqual(testCase, exp, len);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ap, ph2] = phmap(score,label,labelvalue)
hd2_ind=find(score<=2);
if isempty(hd2_ind)
    ph2=0;
else
    ph2=sum(label(hd2_ind)==labelvalue)/length(hd2_ind);
end
ap = apcal(score,label,labelvalue);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function out = riverplot_line(x, y, color, thickness, steepness, varargin)
% out = riverplot_line(x, y, color, thickness, steepness, varargin)
%
% Draw sigmoidal line from 2-coord point x to 2-coord point y
% x = [xcoord ycoord] for left-hand point
% y = [xcoord ycoord] for right-hand point
%
% x = rect1.topright
% y = rect2.topleft

% Programmers' notes:
% 8/21/2017 Stephan Geuter
% changed steepnes treatment computation. Steepness coefficient is now used
% in sigmoid function to control steepnes and curvature at tangets (box
% edges). shift in xstartpoints is removed. Makes ribbons smoother.
%
%

% x, y inputs can be [x,y,z] coordinate triplets or clusters.
if isstruct(x), x = x.mm_center(1:2); end
if isstruct(y), y = y.mm_center(1:2); end

% make middle points with sigmoid function

xdiff = y(1) - x(1);  % difference in x position
ydiff = y(2) - x(2);

% xstartpoints = [x(1) + xdiff * steepness(1) y(1) - xdiff * steepness(1)]';
xstartpoints = [x(1) y(1)]';

% sigmoid reference curve
% sigmoid = inline('p(1) .* ( 1 ./ (1 + p(2)*exp(-p(3)*x)) )','p','x');
sigmoid = @(p,x) (p(1) .* ( 1 ./ (1 + p(2)*exp(-p(3)*x)) ));
xx = linspace(-5, 5, 50); % generate 50 points with standard sigmoid
%                         % Range determines steepness bounds, larger range
%                         = more steep
yy = sigmoid([1 1 1+steepness(1)], xx)'; % add steepness coefficient to have ribbons meet the box edges

xmidpoints = linspace(xstartpoints(1), xstartpoints(end), 50)';
ymidpoints = x(2) + yy * ydiff;

xcoords = [x(1); xstartpoints(1); xmidpoints; xstartpoints(end); y(1)];
ycoords = [x(2); x(2); ymidpoints; y(2); y(2)];

%     % make x, y bend percents
%     if length(steepness) == 1
%         steepness = repmat(steepness, 1, 2);
%     end
%
%     % make 4 coords, adding 2 middle points, so we can bend
%     xdiff = y(1) - x(1);  % difference in x position
%     xmidpoints = [x(1) + xdiff * steepness(1) x(1) + xdiff * .5  y(1) - xdiff * steepness(1)]';
%
%     xcoords = [x(1); xmidpoints; y(1)];
%
%     ydiff = (y(2) - x(2)) ./ 2;
%     ycoords = [x(2); x(2); x(2) + ydiff; y(2); y(2)];

%     plot(xcoords, ycoords, 'go');
%
%      if any(steepness)
%
%     %[p,sse,fit] = nonlin_fit(ycoords, xcoords - mean(xcoords), 'start',[1 1 1]);
%
% sigmoid = inline('p(1) .* ( 1 ./ (1 + p(2)*exp(-p(3)*x)) )','p','x');
% xcoords2 = linspace(xcoords(1), xcoords(end), 100);
%
% ycoords2 = sigmoid([1 1 1], xcoords2 - mean(xcoords2));
% ycoords2 = ycoords(1) + ydiff * ycoords2;
%
% plot(xcoords2, ycoords2, 'rx');

%         % bow out: curved line
%         n = length(xcoords);
%
%         nsamples = [];
%         if length(varargin) > 0
%             nsamples = varargin{1};
%         end
%
%         if isempty(nsamples)
%             nsamples = 10 * n;
%         end
%
%         t = 1:n;
%         ts = 1:((n-1)/(nsamples-1)):n;          % spline grid
%
%         xcoords = spline(t, xcoords, ts);
%         ycoords = spline(t, ycoords, ts);
%
%     end

    h = plot(xcoords, ycoords,'Color',color,'LineWidth',thickness);

    out.h = h;
    out.xcoords = xcoords;
    out.ycoords = ycoords;

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function line = UIAddLine(dir,varargin)

%UIAddLine - Interactively add horizontal/vertical line to existing plot.
%
%  USAGE
%
%    line = UIAddLine(dir,<options>)
%
%    dir            'v' for vertical, 'h' for horizontal
%    <options>      options for function plot
%
%  NOTE
%
%    Click where the line should appear, then hit 'ENTER'.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help UIAddLine'' for details).');
end

hold on;
x = ginput;
if strcmp(lower(dir),'v'),
	line = plot([x(end,1) x(end,1)],ylim,varargin{:});
else
	line = plot(xlim,[x(end,2) x(end,2)],varargin{:});
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [xhat,e,k,theta0,P,b] = celp9600(x,N,L,M,c,cb,Pidx)
%  celp --> 9600 bps CELP analyzer and synthesizer.
%
%    [xhat,e,k,theta0,P,b] = celp9600(x,N,L,M,c,cb,Pidx)
%
%    The function implements a 9600 bps CELP analyzer and synthesizer,
%    if the speech signal is sampled at 8 kHz, the frame size N is 20 ms
%    (160 samples), and the block duration L for the excitation sequence
%    selection is 5 ms (40 samples). Furthermore, assume that the
%    codebook has 1024 sequences which require 10 bit to send the index k,
%    and that the lag of the pitch filter, P, is searched in the range 16
%    to 160 which require 8 bit to represent exactly. Thus, the quantization
%    procedure only affects the M = 10 LPC parameters (inverse sine) coded
%    by 6 bits, the gain Theta0 coded by 7 bit, and the pitch filter coeff.
%    b coded by 8 bit.

Nx = length(x);                         % Signal length.
F  = fix(Nx/N);                         % No. of frames.
J  = N/L;                               % No. blocks per frame.

% Initialize output signals.
xhat   = zeros(Nx,1);                   % Synthesized signal.
e      = zeros(Nx,1);                   % Excitation signal.
k      = zeros(J,F);                    % Columns are excitation
theta0 = zeros(J,F);                    % parameters per frame.
P      = zeros(J,F);
b      = zeros(J,F);

ebuf  = zeros(Pidx(2),1);               % Vectors with previous excitation
ebuf2 = ebuf; bbuf = 0;                 % samples.
Zf = []; Zw = []; Zi = [];              % Memory hangover in filters.

for (f=1:F)
  %fprintf(1,'... Frame no. %g out of %g.\n',f,F);
  n = (f-1)*N+1:f*N;                    % Time index of current speech frame.

  [kappa,kf,theta0f,Pf,bf,ebuf,Zf,Zw] = celpana(x(n),L,M,c,cb,Pidx,bbuf,...
                                                                ebuf,Zf,Zw);

  sigma  = 2/pi*asin(kappa);
  sigma  = udecode(uencode(sigma,6),6);
  kappa  = sin(pi/2*sigma);
  theta0 = udecode(uencode(theta0,7,0.2),7,0.2);
  b      = udecode(uencode(b,8,1.4),8,1.4);

  [xhat(n),ebuf2,Zi] = celpsyn(cb,kappa,kf,theta0f,Pf,bf,ebuf2,Zi);

  % Output excitation signal and parameters for current frame.
  e(n)        = ebuf(Pidx(2)-N+1:Pidx(2));
  k(:,f)      = kf;
  theta0(:,f) = theta0f;
  P(:,f)      = Pf;
  b(:,f)      = bf; bbuf = bf(J);       % Last estimated b used in next frame.
end

author: Sable

repo: mcbench-benchmarks


text:
function  [n,dspchans] = plx_chanmap(filename)
% plx_chanmap(filename): return map of raw DSP channel numbers for each channel
%
% [n,dspchans] = plx_chanmap(filename)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%
% OUTPUT:
%   n - number of spike channels
%   dspchans - 1 x n array of DSP channel numbers
%
% Normally, there is one channel entry in the .plx for for each raw DSP channel,
% so the mapping is trivial dspchans[i] = i.
% However, for certain .plx files saved in some ways from OFS (notably after
% loading data files from other vendors), the mapping can be more complex.
% E.g. there may be only 2 non-empty channels in a .plx file, but those channels
% correspond to raw DSP channel numbers 7 and 34. So in this case NChans = 2,
% and dspchans[1] = 7, dspchans[2] = 34.
% The plx_ routines that return arrays always return arrays of size NChans. However,
% routines that take channels numbers as arguments always expect the raw DSP
% channel number.  So in the above example, to get the timestamps from unit 4 on
% the second channel, use
%   [n,ts] = plx_ts(filename, dspchans[2], 4 );

if nargin ~= 1
    error 'expected 1 input argument';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    pl2 = PL2GetFileIndex(filename);
    n = numel(pl2.SpikeChannels);
    dspchans = 1:n;
    return;
end

[n,dspchans] = mexPlex(26, filename);

end


text:
function e = boundedges ( p, t )

%*****************************************************************************80
%
%% BOUNDEDGES finds the boundary edges in a triangular mesh.
%
%  Discussion:
%
%    You may need this routine if you need to enforce boundary
%    conditions in a PDE.
%
%    The 3D version of this code is called SURFTRI.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,2), the coordinates of a set of nodes.
%
%    Input, integer T(NT,1:3), a list of the nodes which make up each triangle
%    of a triangulation of the nodes in P.
%
%    Output, integer E(*,*), ?
%

%
%  Form all edges, non-duplicates are boundary edges
%
  edges = [t(:,[1,2]);
           t(:,[1,3]);
           t(:,[2,3])];

  node3 = [t(:,3);t(:,2);t(:,1)];
  edges = sort(edges,2);
  [foo,ix,jx] = unique(edges,'rows');
  vec = histc(jx,1:max(jx));
  qx = find(vec==1);
  e = edges(ix(qx),:);
  node3 = node3(ix(qx));
%
%  Orientation
%
  v1 = p(e(:,2),:)-p(e(:,1),:);
  v2 = p(node3,:)-p(e(:,1),:);
  ix = find(v1(:,1).*v2(:,2)-v1(:,2).*v2(:,1)>0);
  e(ix,[1,2]) = e(ix,[2,1]);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [surface_handles, colors] = isosurface(r, varargin)
% Create a series of surfaces in different colors, one for each region
% - Options for single color
%
% [surface_handles, colors] = isosurface(r, [optional arguments])
%
% optional arguments:
% - Any optional inputs to imageCluster, e.g., 'alpha'
% - 'colors', followed by single color in { } or cell array of multiple colors
% - 'nomatchleftright' or 'nosymmetric', do not match colors across hemispheres (left/right)
%               Note: The default matches, and may override your colors.
%
% Examples:
% atlasfile = which('Morel_thalamus_atlas_object.mat');
% load(atlasfile)
%
% surface_handles = isosurface(r);
% surface_handles = isosurface(r, 'alpha', .5);
% surface_handles = isosurface(r, 'alpha', .5, 'nomatchleftright');
%
% view(135, 30);
% lightRestoreSingle;
% lightFollowView;
%
% load(which('CIT168_MNI_subcortical_atlas_object.mat'));
% r = atlas2region(atlas_obj);
%
% surface_handles = isosurface(r, 'alpha', .5, 'color', {[.3 .6 .4] [.5 .4 .2]});
% surface_handles = isosurface(r, 'alpha', .5, 'color', {[.3 .6 .4] [.5 .4 .2]}, 'nomatchleftright');
% p = addbrain('hires right');
% lightFollowView;

k = length(r);

% ..
%    DEFAULTS AND INPUTS
% ..

% input_args = varargin;           % save for later
colors = scn_standard_colors(k); % generate colors
matchcolorsleftright = true;

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            % do nothing for inputs passed on to imageCluster
            case {'alpha'}

            case {'color', 'colors'}
                colors = varargin{i+1}; varargin{i+1} = []; varargin{i} = [];
                if length(colors) == 1 % single color; matchcolorsleftright will overwrite this
                    matchcolorsleftright = false;
                end

            case {'nomatchleftright', 'nosymmetric'}, matchcolorsleftright = false; varargin{i} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

if iscolumn(colors), colors = colors'; end

% Handle single-color input, and other cases where too few colors:
while length(colors) < k, colors = [colors colors]; end

if matchcolorsleftright
    colors = match_colors_left_right(r);
end

colors = colors(1:length(r));

% -------------------------------------------------------------------------
% Make surfaces
% -------------------------------------------------------------------------

cl = region2struct(r);

surface_handles = [];

for i = 1:k

    try
        out = imageCluster('cluster', cl(i), 'color', colors{i}, varargin{:});

        surface_handles(i) = out;
    catch
        disp('Error imaging isosurface; too few voxels?')
        %surface_handles(i) = [];
    end

%     % Isocaps, if needed
%
%     % pp = [];
%     isocap = isocaps(mesh_struct.X, mesh_struct.Y, mesh_struct.Z, V, mythresh);
%     p(end + 1) = patch(isocap, 'FaceColor', 'interp','EdgeColor', 'none', 'FaceAlpha',1);

end

% -------------------------------------------------------------------------
% Lighting, etc.
% -------------------------------------------------------------------------

lightRestoreSingle;
lighting gouraud;
axis vis3d image
material dull

end % function


text:
function varargout = legcoeffs(varargin) %#ok<STOUT>
%LEGCOEFFS   LEGCOEFFS does not support TRIGTECH objects.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% TODO: Add support for this.

error('CHEBFUN:TRIGTECH:legcoeffs:nosupport', ...
      'LEGCOEFFS does not support TRIGTECH objects.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [SACdata,SeisData,sacfiles]=sacpc2mat(varargin)
% [SACdata,SeisData,filenames] = SACPCMAT('file1','file2',..., 'filen' )
%
% reads n SAC files file1, file2, filen (SAC files are assumed to have
% PC byte order) and converts them to matlab
% format. The filenames can contain globbing characters (e.g. * and ?).
% These are expanded and all matching files loaded.
%
% SACPCMAT( cellarray ) where cellarray={'file1','file2',...,'filen'}
% is equivalent to the standard form.
%
% SACdata is an n x 1 struct array containing the header variables
%         in the same format as is obtained by using MAT function
%         of SAC2000.
%         SACdata(i).trcLen contains the number of samples.
%
% SeisData is an m x n array (where m=max(npts1, npts2, ...) )
%         containing the actual data.
%
% filenames is a n x 1 string cell array with the filenames actually read.
%
% Note that writing
%
%  [SACdata,SeisData] = sacsun2mat('file1','file2',..., 'filen' )
%
% is equivalent to the following sequence
%
% sac2000
% READ file1 file2 .. filen
% MAT
%
% (in fact the failure of above sequence to work properly on my
% system motivated this script).
%
%
% SACPC2MAT was written by F Tilmann (tilmann@esc.cam.ac.uk)
% based on sac_sun2pc_mat  by C. D. Saragiotis (I copied the
% routines doing the actual work from this code but
% used a different header structure and made the routine
% flexible).
% It was tested on MATLAB5 on a PC but
% should work on newer versions, too.
%
% (C) 2004
%

F = 4-1; % float byte-size minus 1;
K = 8-1; % alphanumeric byte-size minus 1
L = 4-1; % long integer byte-size minus 1;

fnames={};
for i=1:nargin
  if ischar(varargin{i})
    fnames=cell([fnames; cellstr(varargin{i})]);
  elseif iscellstr(varargin{i}) & size(varargin{i},1)==1
    fnames=cell([fnames; varargin{i}']);
  elseif iscellstr(varargin{i}) & size(varargin{i},2)==1
    fnames=cell([fnames; varargin{i}]);
  end
end
% expand globs
sacfiles={};k=1;
for i=1:length(fnames)
  dirlist=dir(fnames{i});
  for j=1:length(dirlist)
    if ~dirlist(j).isdir
      sacfiles{k,1}=dirlist(j).name;
      k=k+1;
    end
  end
end

maxnpts=0;
for i=1:length(sacfiles)
  fid=fopen(sacfiles{i},'rb');
  if fid==-1
    error(sprintf('Could not open SAC file %s',fnames{i}))
  end
  SACdata(i,1)=readSacHeader(fid,F,K,L);
  npts=SACdata(i).trcLen;
  if npts>maxnpts
    maxnpts=npts;
  end
  fprintf('Processing file %d: %s\n',i,sacfiles{i});
  SeisData(npts,i)=0;   % Magnify seis matrix if necessary
  SeisData(:,i)=[ readSacData(fid,npts,F+1); zeros(maxnpts-npts,1)];    % Have to pad with zeros if new data have less data points than some previously read file
end

function hdr = readSacHeader(fileId,F,K,L)
% hdr = readSacHeader(FID)
% sacReadAlphaNum reads the SAC-format header fields and returns most of them.
%
%    The output variable, 'hdr' is a structure, whose fields are
%    the fields as in the SACdata structure generated by SAC's
%    matlab command MAT
%
%    Created by C. D. Saragiotis, August 5th, 2003, modified by F Tilmann
headerBytes = 632;
chrctr = fread(fileId,headerBytes,'uchar');
chrctr = chrctr(:)';
% Read floats
hdr.times.delta  = sacReadFloat(chrctr(1:1+F)); % increment between evenly spaced samples
%hdr.DEPMIN = sacReadFloat(chrctr(5:5+F)); % MINimum value of DEPendent variable
%hdr.DEPMAX = sacReadFloat(chrctr(9:9+F)); % MAXimum value of DEPendent variable
% (not currently used) SCALE  = sacReadFloat(chrctr(13:13+F)); % Mult SCALE factor for dependent variable
%hdr.ODELTA = sacReadFloat(chrctr(17:17+F)); % Observd increment if different than DELTA
hdr.times.b      = sacReadFloat(chrctr(21:21+F)); % Begining value of the independent variable
hdr.times.e      = sacReadFloat(chrctr(25:25+F)); % Ending value of the independent variable
hdr.times.o      = sacReadFloat(chrctr(29:29+F)); % event Origin time
hdr.times.a      = sacReadFloat(chrctr(33:33+F)); % first Arrival time
hdr.times.t0 = sacReadFloat(chrctr(41:41+F));
hdr.times.t1 = sacReadFloat(chrctr(45:45+F));
hdr.times.t2 = sacReadFloat(chrctr(49:49+F));
hdr.times.t3 = sacReadFloat(chrctr(53:53+F));
hdr.times.t4 = sacReadFloat(chrctr(57:57+F));
hdr.times.t5 = sacReadFloat(chrctr(61:61+F));
hdr.times.t6 = sacReadFloat(chrctr(65:65+F));
hdr.times.t7 = sacReadFloat(chrctr(69:69+F));
hdr.times.t8 = sacReadFloat(chrctr(73:73+F));
hdr.times.t9 = sacReadFloat(chrctr(77:77+F));
hdr.times.f      = sacReadFloat(chrctr(81:81+F)); % Fini of event time

hdr.response = sacReadFloat(reshape(chrctr(85:85+10*(F+1)-1),F+1,10)');

hdr.station.stla   = sacReadFloat(chrctr(125:125+F)); % STation LAttitude
hdr.station.stlo   = sacReadFloat(chrctr(129:129+F)); % STation LOngtitude
hdr.station.stel   = sacReadFloat(chrctr(133:133+F)); % STation ELevation
hdr.station.stdp   = sacReadFloat(chrctr(137:137+F)); % STation DePth below surface

hdr.event.ecla   = sacReadFloat(chrctr(141:141+F)); % EVent LAttitude
hdr.event.evlo   = sacReadFloat(chrctr(145:145+F)); % EVent LOngtitude
hdr.event.evel   = sacReadFloat(chrctr(149:149+F)); % EVent ELevation
hdr.event.evdp   = sacReadFloat(chrctr(153:153+F)); % EVent DePth below surface
hdr.event.mag   = sacReadFloat(chrctr(157:157+F)); % EVent DePth below surface

userdata=sacReadFloat(reshape(chrctr(161:161+10*(F+1)-1),F+1,10)');

hdr.evsta.dist   = sacReadFloat(chrctr(201:201+F)); % station to event DISTance (km)
hdr.evsta.az     = sacReadFloat(chrctr(205:205+F)); % event to station AZimuth (degrees)
hdr.evsta.baz    = sacReadFloat(chrctr(209:209+F)); % station to event AZimuth (degrees)
hdr.evsta.gcarc  = sacReadFloat(chrctr(213:213+F)); % station to event Great Circle ARC length (degrees)

%hdr.DEPMEN = sacReadFloat(chrctr(225:225+F)); % MEaN value of DEPendent variable

hdr.station.cmpaz  = sacReadFloat(chrctr(229:229+F)); % CoMPonent AZimuth (degrees clockwise from north)
hdr.station.cmpinc = sacReadFloat(chrctr(233:233+F)); % CoMPonent INCident angle (degrees from vertical)

hdr.llnl.xminimum = sacReadFloat(chrctr(237:237+L));
hdr.llnl.xmaximum = sacReadFloat(chrctr(241:241+L));
hdr.llnl.yminimum = sacReadFloat(chrctr(245:245+L));
hdr.llnl.ymaximum = sacReadFloat(chrctr(249:249+L));

% Read long integers
hdr.event.nzyear = sacReadLong(chrctr(281:281+L)); % GMT YEAR
hdr.event.nzjday = sacReadLong(chrctr(285:285+L)); % GMT julian DAY
hdr.event.nzhour = sacReadLong(chrctr(289:289+L)); % GMT HOUR
hdr.event.nzmin  = sacReadLong(chrctr(293:293+L)); % GMT MINute
hdr.event.nzsec  = sacReadLong(chrctr(297:297+L)); % GMT SECond
hdr.event.nzmsec = sacReadLong(chrctr(301:301+L)); % GMT MilliSECond

hdr.llnl.norid = sacReadLong(chrctr(309:309+L));
hdr.llnl.nevid = sacReadLong(chrctr(313:313+L));

hdr.trcLen   = sacReadLong(chrctr(317:317+L)); % Number of PoinTS per data component

hdr.llnl.nwfid = sacReadLong(chrctr(325:325+L));
hdr.llnl.nxsize = sacReadLong(chrctr(329:329+L));
hdr.llnl.nysize = sacReadLong(chrctr(333:333+L));

hdr.descrip.iftype = sacReadLong(chrctr(341:341+L)); % File TYPE
hdr.descrip.idep   = sacReadLong(chrctr(345:345+L)); % type of DEPendent variable
hdr.descrip.iztype  = sacReadLong(chrctr(349:349+L)); % reference time equivalence

hdr.descrip.iinst  = sacReadLong(chrctr(357:357+L)); % type of recording INSTrument

% Before there were floats read here?
hdr.descrip.istreg = sacReadLong(chrctr(361:361+F)); % STation geographic REGion
hdr.descrip.ievreg = sacReadLong(chrctr(365:365+F)); % EVent geographic REGion
hdr.descrip.ievtyp = sacReadLong(chrctr(369:369+F)); % EVent geographic REGion
hdr.descrip.iqual  = sacReadLong(chrctr(373:373+F)); % QUALity of data
hdr.descrip.isynth = sacReadLong(chrctr(377:377+F)); % SYNTHetic data flag

hdr.event.imagtyp  = sacReadLong(chrctr(381:381+F));
hdr.event.imagsrc  = sacReadLong(chrctr(385:385+F));

% no logical SAC header variables in matlab SACdata structure!
% previous version set these defaults
% $$$ % SAC defaults
% $$$ hdr.IEVTYPE = 'IUNKN';
% $$$ %IEVTYPE= sacReadFloat(chrctr(369:369+F)); % EVent TYPE
% $$$ hdr.LEVEN = 'TRUE';  % true, if data are EVENly spaced (required)
% $$$ %LEVEN= sacReadFloat(chrctr(421:521+F)); % true, if data are EVENly spaced (required)
% $$$ hdr.LPSPOL = 'FALSE'; % true, if station components have a PoSitive POLarity
% $$$ %LPSPOL= sacReadFloat(chrctr(425:425+F)); % true, if station components have a PoSitive POLarity
% $$$ hdr.LOVROK = 'FALSE'; % true, if it is OK to OVeRwrite this file in disk
% $$$ %LOVROK= sacReadFloat(chrctr(429:429+F)); % true, if it is OK to OVeRwrite this file in disk
% $$$ hdr.LCALDA = 'TRUE'; % true, if DIST, AZ, BAZ and GCARC are to be calculated from station and event coordinates
% $$$ %LCALDA= sacReadFloat(chrctr(433:433+F)); % true, if DIST, AZ, BAZ and GCARC are to be calculated from station and event coordinates

% Read alphanumeric data
hdr.station.kstnm  = sacReadAlphaNum(chrctr(441:441+K)); % STation NaMe
hdr.event.kevnm  = sacReadAlphaNum(chrctr(449:449+2*(K+1)-1)); % EVent NaMe
%hdr.KHOLE  = sacReadAlphaNum(chrctr(465:465+K)); % HOLE identification, if nuclear event
hdr.times.ko     = sacReadAlphaNum(chrctr(473:473+K)); % event Origin time identification
hdr.times.ka     = sacReadAlphaNum(chrctr(481:481+K)); % first Arrival time identification

hdr.times.kt0 = sacReadAlphaNum(chrctr(489:489+K));
hdr.times.kt1 = sacReadAlphaNum(chrctr(497:497+K));
hdr.times.kt2 = sacReadAlphaNum(chrctr(505:505+K));
hdr.times.kt3 = sacReadAlphaNum(chrctr(513:513+K));
hdr.times.kt4 = sacReadAlphaNum(chrctr(521:521+K));
hdr.times.kt5 = sacReadAlphaNum(chrctr(529:529+K));
hdr.times.kt6 = sacReadAlphaNum(chrctr(537:537+K));
hdr.times.kt7 = sacReadAlphaNum(chrctr(545:545+K));
hdr.times.kt8 = sacReadAlphaNum(chrctr(553:553+K));
hdr.times.kt9 = sacReadAlphaNum(chrctr(561:561+K));

hdr.times.kf     = sacReadAlphaNum(chrctr(569:569+K)); % Fini identification
kuser0 = sacReadAlphaNum(chrctr(577:577+K)); % USER-defined variable storage area
kuser1 = sacReadAlphaNum(chrctr(585:585+K)); % USER-defined variable storage area
kuser2 = sacReadAlphaNum(chrctr(593:593+K)); % USER-defined variable storage area
hdr.station.kcmpnm = sacReadAlphaNum(chrctr(601:601+K)); % CoMPonent NaMe
hdr.station.knetwk = sacReadAlphaNum(chrctr(609:609+K)); % name of seismic NETWorK
%hdr.KDATRD = sacReadAlphaNum(chrctr(617:617+K)); % DATa Recording Date onto the computer
%hdr.KINST  = sacReadAlphaNum(chrctr(625:625+K)); % generic name of recording INSTrument

usercell=num2cell(userdata);
[usercell{find(userdata==-12345)}]=deal([]);
[hdr.user(1:10).data]=deal(usercell{:});
[hdr.user(1:10).label]=deal(kuser0, kuser1,kuser2,[], [], [], [], [], [], []);

function X = readSacData(fid,N,F)
% function data = readSacData('filename',NPTS,floatByteSize)
chrctr = fread(fid,N*F,'uchar');
x=reshape(chrctr,F,N)';
%x
X = sacReadFloat(x);

function lNumber = sacReadLong(cb)
% reads long integers (4 bytes long)
% cb is the character buffer
cb = cb(:);
% SUN sac format use: lNumber = (256.^(3:-1:0))*cb;
lNumber = (256.^(0:3))*cb;        % changed line
if lNumber == -12345, lNumber = []; end

function fNumber = sacReadFloat(cb)
% reads floats (4 bytes long)
% cb is the character buffer
C = size(cb,1);
% SUN sac format use: stringOfBitSequence = [dec2bin(cb(:,1),8) dec2bin(cb(:,2),8) dec2bin(cb(:,3),8) dec2bin(cb(:,4),8)];
stringOfBitSequence = [dec2bin(cb(:,4),8) dec2bin(cb(:,3),8) ...
            dec2bin(cb(:,2),8) dec2bin(cb(:,1),8)];  %changed line
bitSequence = stringOfBitSequence=='1';
fSign = -2*bitSequence(:,1)+1;
fExp = bitSequence(:,2:9)*(2.^(7:-1:0)') - 127;
fMantissa = [ones(C,1) bitSequence(:,10:32)]*(2.^(0:-1:-23)');
fNumber = fSign.*fMantissa.*(2.^fExp);
isZeroCheck = sum(bitSequence')'==0;
fNumber = fNumber.*(~isZeroCheck);
if C==1 & fNumber == -12345, fNumber = []; end

function alphaNum = sacReadAlphaNum(cb)
% reads alphanumeric data (8 or 16 bytes long). If it cb is empty, it returns a ' '
% cb is the character buffer
K = max(size(cb));
alphaNumTemp = char(cb);
if K == 8
    if alphaNumTemp == '-12345  '
        alphaNum = [];
    else
        alphaNum = alphaNumTemp;
    end
else
    if K == 16
        if alphaNumTemp == '-12345   -12345 ' | alphaNumTemp == '-12345          '
            alphaNum = [];
        else
            alphaNum = alphaNumTemp;
        end
    end
end

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function TC = arduino_lab1(heater)

persistent icount a

if (isempty(icount))
   try
      a = arduino;
      disp(a)
   catch
      warning('Unable to connect, user input required')
      disp('For Windows:')
      disp('  Open device manager, select "Ports (COM & LPT)"')
      disp('  Look for COM port of Arduino such as COM4')
      disp('For MacOS:')
      disp('  Open terminal and type: ls /dev/*.')
      disp('  Search for /dev/tty.usbmodem* or /dev/tty.usbserial*. The port number is *.')
      disp('For Linux')
      disp('  Open terminal and type: ls /dev/tty*')
      disp('  Search for /dev/ttyUSB* or /dev/ttyACM*. The port number is *.')
      disp('')
      com_port = input('Specify port (e.g. COM4 for Windows or /dev/ttyUSB0 for Linux): ','s');
      a = arduino(com_port,'Uno');
      disp(a)
   end
   icount = 0;
end

% voltage read functions
v1 = @() readVoltage(a, 'A0');
v2 = @() readVoltage(a, 'A2');

% temperature calculations as a function of voltage for TMP36
TC = @(V) (V - 0.5)*100.0;          % Celsius
TK = @(V) TC(V) + 273.15;           % Kelvin
TF = @(V) TK(V) * 9.0/5.0 - 459.67; % Fahrenheit

% temperature read functions
T1C = @() TC(v1());
T2C = @() TC(v2());

% LED function (0 <= level <= 1)
led = @(level) writePWMDutyCycle(a,'D9',max(0,min(1,level)));  % ON

% heater output (0 <= heater <= 100)
% limit to 0-0.9 (0-100%)
h1 = @(level) writePWMDutyCycle(a,'D3',max(0,min(100,level))*0.9/100);
% limit to 0-0.5 (0-100%)
h2 = @(level) writePWMDutyCycle(a,'D5',max(0,min(100,level))*0.5/100);

% increment counter
icount = icount + 1;

% read temperature
TC = T1C();
% write heater
h1(heater);

% indicate high temperature with LED
if TC > 40
    led(1)
else
    led(0)
end

end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function s = i4_to_binary ( i )

%*****************************************************************************80
%
%% I4_TO_BINARY produces the binary representation of an integer.
%
%  Example:
%
%     I       S
%    --  ------
%     1      '1'
%     2     '10'
%     3     '11'
%     4    '100'
%    -9  '-1001'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, an integer to be represented.
%
%    Output, string S, the binary representation.
%
  s = [];

  if ( i == 0 )
    s = strcat ( '0', s );
    return
  end
%
%  Write the string backwards.
%
  i_copy = abs ( i );

  while ( 0 < i_copy )

    if ( mod ( i_copy, 2 ) == 1 )
      s = strcat ( '1', s );
    else
      s = strcat ( '0', s );
    end

    i_copy = floor ( i_copy / 2 );

  end
%
%  Attach a minus sign, if needed.
%
  if ( i < 0 )
    s = strcat ( '-', s );
  end

  return
end


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%

sim = exp(-sumsq(x1 - x2) / (2 * sigma^2));

% =============================================================

end

author: zsiciarz

repo: ml-coursera


text:
function [ a, b ] = p47_lim ( )

%*****************************************************************************80
%
%% P47_LIM returns the integration limits for problem 47.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end


text:
%% OPT05_RUN
%
%  Modified:
%
%    09 January 2008
%
   %---------------------------------------------------------------------
   %  Running testcase, from D+S, pp. 157-8; x_* = [.1,10].
   %  This tests the scaling implementation.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_5 as scalar optimization problem:\n' );
   fprintf('Exact solution (0.1, 10)\n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt05_fgh';
   options = [];
   options.verbose            = 0;
   options.method             = 'newton';
   options.globalization      = 'trust_region';
   options.step_tolerance     = 1.e-14;
   options.gradient_tolerance = 1.e-14;
   options.max_fevals         = 151;
   options.max_iterations     = 151;

   alpha = 10;
   fprintf ( 'Program scale ALPHA = %f\n', alpha );
   x0 = [-1.2/alpha, alpha];
   fprintf('Unscaled:\n')
%
% should require ~94 iterations
%
   x = entrust(fname, x0, options, alpha);

   fprintf('Newton produced  (%10.7e, %10.7e)\n\n',x(1),x(2))
   f = opt05_fgh ( x, 'f', alpha );
   fprintf('Value of F(X) = %f\n', f );

   options.scale_x            = [ 1/alpha; alpha ];
   fprintf('Scaled:\n')

   % should require ~24 iterations
   x = entrust(fname, x0, options, alpha);

   fprintf('Newton produced  (%10.7e, %10.7e)\n\n',x(1),x(2))
   f = opt05_fgh ( x, 'f', alpha );
   fprintf('Value of F(X) = %f\n', f );

   %---------------------------------------------------------------------
   %  Test Gauss-Newton strategies.
   %---------------------------------------------------------------------
   fprintf('---------------------------------------------------------\n')
   fprintf('Running testcase_5 as least squares problem: \n')
   fprintf('Exact solution (0.1, 10)\n')
   fprintf('---------------------------------------------------------\n')
   fname = 'opt05_rj';
   options = [];
   options.verbose            = 0;
   options.method             = 'gauss_newton';
   options.step_tolerance     = 1.e-15;
   options.globalization      = 'none';
   options.gradient_tolerance = 1.e-10;
   options.max_iterations     = 40;

   alpha = 10;
   fprintf ( 'Program scale ALPHA = %f\n', alpha );
   x0 = [-1.2/alpha, alpha];
   fprintf('Unscaled:\n')

   x = entrust(fname, x0, options, alpha);

   fprintf('Gauss-Newton produced  (%10.7e, %10.7e)\n\n',x(1),x(2))
   [ res, jac ] = opt05_rj ( x, 'f', alpha );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );

   options.scale_x            = [ 1/alpha; alpha ];
   fprintf('Scaled:\n')

   x = entrust(fname, x0, options, alpha);

   fprintf('Gauss-Newton produced  (%10.7e, %10.7e)\n\n',x(1),x(2))
   [ res, jac ] = opt05_rj ( x, 'f', alpha );
   fprintf('Norm of RES(X) = %f\n', norm ( res ) );


text:
% [INPUT]
% handle = A handle to the plotting function to be executed.

function safe_plot(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('handle',@(x)validateattributes(x,{'function_handle'},{'scalar'}));
    end

    ip.parse(varargin{:});
    ipr = ip.Results;

    nargoutchk(0,0);

    safe_plot_internal(ipr.handle);

end

function safe_plot_internal(handle)

    persistent ids;

    name = func2str(handle);
    name = regexprep(name,'^@\([^)]*\)','');
    name = regexprep(name,'\([^)]*\)$','');

    try
        id = [upper(name) '-' upper(char(java.util.UUID.randomUUID()))];
    catch
        id = randi([0 10000000]);

        while (ismember(id,ids))
            id = randi([0 100000]);
        end

        ids = [ids; id];
        id = [upper(name) '-' sprintf('%08s',num2str(id))];
    end

    try
        handle(id);
    catch e
        delete(findobj('Type','Figure','Tag',id));

        try
            r = getReport(e,'Extended','Hyperlinks','off');
            r = strsplit(r,new_line());
            r = cellfun(@(x)['  ' x],r,'UniformOutput',false);
            r = strjoin(r,new_line());
            r = clear_text(r);

            warning('MATLAB:SystemicRisk',['The following exception occurred in the plotting function ''' name ''':' new_line() r]);
        catch
            warning('MATLAB:SystemicRisk',['The following exception occurred in the plotting function ''' name ''':' new_line() e.message]);
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [structdis,iMinusMu,mu,sigma]= divisiveNormalization(imdist)

    window = fspecial('gaussian',7,7/6);
    window = window/sum(sum(window));

    mu = filter2(window, imdist, 'same');
    mu_sq = mu.*mu;

    sigma = sqrt(abs(filter2(window, imdist.*imdist, 'same') - mu_sq));
    iMinusMu = (imdist-mu);
    structdis =iMinusMu./(sigma +1);
end

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function c = long2intval(C)
%LONG2INTVAL  Conversion long to intval (with correct rounding)
%
%  c = long2intval(C)
%

% written  12/30/98     S.M. Rump
% modfied  02/09/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 07/31/05     S.M. Rump  header corrected
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  INTLAB_LONG_BETA = getappdata(0,'INTLAB_LONG_BETA');
  INTLAB_LONG_LOGBETA = getappdata(0,'INTLAB_LONG_LOGBETA');
  INTLAB_LONG_ERROR = getappdata(0,'INTLAB_LONG_ERROR');

  % extra treatment of zero component
  indexzero = ( C.exponent==-inf );

  n = size(C.mantissa,1);
  % take maximum first 80 bits
  precC = size(C.mantissa,2);
  p = min( precC , ceil(80/INTLAB_LONG_LOGBETA)+1 );
  if INTLAB_LONG_ERROR
    if p<precC
      C.error = errorupdate( 1 , C.error , 0 , ...
                  1 , any(C.mantissa(:,p+1:precC)~=0,2) , C.exponent-precC );
    end
    Cerr = C.error.mant .* INTLAB_LONG_BETA.^(C.error.exp-C.exponent+p);
  else
    Cerr = 0;
  end

  factor = ( INTLAB_LONG_BETA.^(-p:-1) )';
  exppos = ( C.exponent>=0 );
  expneg = ~exppos;
  E = floor( 600/INTLAB_LONG_LOGBETA );
  F = INTLAB_LONG_BETA ^ E;

  Cmantp = C.mantissa(:,p);
  setround(-1)
  C.mantissa(:,p) = Cmantp - Cerr;
  cinf = C.mantissa(:,p:-1:1) * factor ;
  if any(exppos)
    cinf(exppos) = ( cinf(exppos)*F ) .* ...
                      ( INTLAB_LONG_BETA.^(C.exponent(exppos)-E) );
  end
  if any(expneg)
    cinf(expneg) = ( cinf(expneg)/F ) .* ...
                      ( INTLAB_LONG_BETA.^(C.exponent(expneg)+E) );
  end

  setround(1)
  C.mantissa(:,p) = Cmantp + Cerr;
  csup = C.mantissa(:,p:-1:1) * factor ;
  if any(exppos)
    csup(exppos) = ( csup(exppos)*F ) .* ...
                      ( INTLAB_LONG_BETA.^(C.exponent(exppos)-E) );
  end
  if any(expneg)
    csup(expneg) = ( csup(expneg)/F ) .* ...
                      ( INTLAB_LONG_BETA.^(C.exponent(expneg)+E) );
  end
  if INTLAB_LONG_ERROR & any(indexzero)
    csup(indexzero) = C.error.mant(indexzero) .* ...
                        INTLAB_LONG_BETA.^C.error.exp(indexzero);
    cinf(indexzero) = - csup(indexzero);
  end

  c = hull( C.sign.*cinf , C.sign.*csup );

  setround(rndold)

repo: OpenMAS


text:
function [vk, f] = ADMM3D_solver(psf,b,solverSettings)
% ADMM solver to compute 3D diffusercam images
% H: Impulse stack 3D array
% b: measurement image from camera
% solverSettings: user defined params. See DiffuserCam_settings.m for
% details

assert(size(psf,1) == size(b,1) || size(psf,2) == size(b,2),'image and impulse have different dimensions');

if ~isfield(solverSettings,'print_interval')
    solverSettings.print_interval = 1;
end

if ~isfield(solverSettings,'disp_figs')
    solverSettings.disp_figs = 1;
end

if ~isfield(solverSettings,'autotune')
    solverSettings.autotune = 1;
end

if ~isfield(solverSettings,'maxIter')
    solverSettings.maxIter = 200;
end

if ~isfield(solverSettings,'regularizer')
    solverSettings.regularizer = 'TV';
end

if ~isfield(solverSettings,'display_func')
    solverSettings.display_func = @(x)x;
end

if ~isfield(solverSettings,'cmap')
    solverSettings.cmap = 'gray';
end

if isfield(solverSettings,'save_vars')
    assert(iscell(solverSettings.save_vars),'solverSettings.save_vars must be a cell array of strings')
else
    solverSettings.save_vars = {'vk'};
end

if strcmpi(solverSettings.save_vars,'')
    solverSettings.save_vars = {'vk'};
end

mu1 = solverSettings.mu1;   %Set initial ADMM parameters
mu2 = solverSettings.mu2;
mu3 = solverSettings.mu3;

[Ny, Nx, Nz] = size(psf);   %Get problem size

% Setup convolutional forward op
p1 = floor(Ny/2);
p2 = floor(Nx/2);
%h(p1,p2,Nz/2) = 1;
pad2d = @(x)padarray(x,[p1,p2],'both');  %2D padding
crop2d = @(x)x(p1+1:end-p1,p2+1:end-p2,:); %2D cropping
crop3d = @(x)crop2d(x(:,:,1));   %3D cropping. This is D
vec = @(X)reshape(X,numel(X),1);
pad3d = @(x)padarray(pad2d(x),[0 0 Nz-1],'post');
psf = circshift(flip(psf,3),ceil(Nz/2)+1,3)/norm(psf(:));  %Shift impulse stack and normalize
Hs = fftn(ifftshift(pad2d(psf)));  %Compute 3D spectrum
Hs_conj = conj(Hs);
clear psf
Hfor = @(x)real((ifftn(Hs.*fftn((x)))));
Hadj = @(x)real((ifftn(Hs_conj.*fftn((x)))));
HtH = abs(Hs.*Hs_conj);

vk = 0*real(Hs);   %Initialize variables. vk is the primal (this is the image you want to find)
xi = vk;  % Dual associated with Mv = nu (boundary condition variables)
rho = vk;  % Dual associated with v = w   (nonnegativity)
Dtb = pad3d(b);

switch lower(solverSettings.regularizer)
    case('tv')
        PsiTPsi = generate_laplacian(vk);
        eta_1 = vk(1:end-1,:,:);  %Duals associatd with Psi v = u (TV sparsity)
        eta_2 = vk(:,1:end-1,:);   %zeros(2*Ny,2*Nx-1,Nz);
        eta_3 = vk(:,:,1:end-1);   %zeros(2*Ny,2*Nx,Nz-1);
        PsiT = @(P1,P2,P3)cat(1,P1(1,:,:),diff(P1,1,1),-P1(end,:,:)) + ...
            cat(2,P2(:,1,:),diff(P2,1,2),-P2(:,end,:)) + ...
            cat(3,P3(:,:,1),diff(P3,1,3),-P3(:,:,end));

        % Sparsifying map
        Psi = @(x)deal(-diff(x,1,1),-diff(x,1,2),-diff(x,1,3));
        [uk1, uk2, uk3] = Psi(vk);
        Lvk1 = uk1;
        Lvk2 = uk2;
        Lvk3 = uk3;
    case('tv_native')
        PsiTPsi = generate_laplacian(vk);
        PsiT = @(P1,P2,P3,P4)cat(1,P1(1,:,:),diff(P1,1,1),-P1(end,:,:)) + ...
            cat(2,P2(:,1,:),diff(P2,1,2),-P2(:,end,:)) + ...
            cat(3,P3(:,:,1),diff(P3,1,3),-P3(:,:,end)) + ...
            solverSettings.tau_n/solverSettings.tau*P4;

        %Sparsifying with gradient and l1
        Psi = @(x)deal(-diff(x,1,1),-diff(x,1,2),-diff(x,1,3),...
            solverSettings.tau_n/solverSettings.tau*x);
        [uk1, uk2, uk3, uk4] = Psi(vk);
        Lvk1 = uk1;
        Lvk2 = uk2;
        Lvk3 = uk3;
        Lvk4 = uk4;
        eta_1 = vk(1:end-1,:,:);  %Duals associatd with Psi v = u (TV sparsity)
        eta_2 = vk(:,1:end-1,:);   %zeros(2*Ny,2*Nx-1,Nz);
        eta_3 = vk(:,:,1:end-1);   %zeros(2*Ny,2*Nx,Nz-1);
        eta_4 = vk;
        PsiTPsi = PsiTPsi + solverSettings.tau_n^2/solverSettings.tau^2;
    case('native')

        PsiTPsi = 1;
        PsiT = @(x)x;
        Psi = @(x)x;   %Identity operator for native sparsity
        uk = vk;
        Lvk = uk;
        eta = uk;
end

v_mult = 1./(mu1*HtH + mu2*PsiTPsi + mu3);  %Denominator of v update (in 3D frequency space)

DtD = pad3d(ones(Ny, Nx, 'like', b)); % Initialize DtD with same datatype as input
nu_mult = 1./(DtD + mu1);   %denominator of nu update

n = 0;  %Initialize number of steps to 0

% Store solver parameters in structure, f
% Initialize residuals with NaNs
f.dual_resid_s = zeros(1,solverSettings.maxIter)./0;
f.primal_resid_s = zeros(1,solverSettings.maxIter)./0;
f.dual_resid_u = f.dual_resid_s;
f.primal_resid_u = f.dual_resid_u;
f.dual_resid_w = f.dual_resid_s;
f.primal_resid_w = f.dual_resid_s;
f.objective = f.primal_resid_u;
f.data_fidelity = f.primal_resid_u;
f.regularizer_penalty = f.primal_resid_u;
Hvkp = vk;
tic
while n<solverSettings.maxIter
    n = n+1;
    Hvk = Hvkp;
    nukp = nu_mult.*(mu1*(xi/mu1 + Hvk) + Dtb);
    wkp = max(rho/mu3 + vk,0);
    switch lower(solverSettings.regularizer)
        case('tv')
            [uk1, uk2, uk3] = DiffuserCam_soft_3d(Lvk1+eta_1/mu2, Lvk2+eta_2/mu2, Lvk3+eta_3/mu2,solverSettings.tau/mu2);
            vkp_numerator = mu3*(wkp-rho/mu3) + ...
                mu2*PsiT(uk1 - eta_1/mu2,uk2 - eta_2/mu2, uk3 - eta_3/mu2) + ...
                mu1*Hadj(nukp - xi/mu1);
        case('tv_native')
            [uk1, uk2, uk3, uk4] = DiffuserCam_soft_3d(Lvk1 + eta_1/mu2, Lvk2 + eta_2/mu2, ...
                Lvk3 + eta_3/mu2, solverSettings.tau/mu2, Lvk4 + eta_4/mu2);
            vkp_numerator = mu3*(wkp-rho/mu3) + ...
                mu2*PsiT(uk1 - eta_1/mu2,uk2 - eta_2/mu2, uk3 - eta_3/mu2, uk4 - eta_4/mu2) + ...
                mu1*Hadj(nukp - xi/mu1);
        case('native')
            uk = DiffuserCam_soft_3d([],[],[],solverSettings.tau_n/mu2,Lvk + eta/mu2);
            vkp_numerator = mu3*(wkp-rho/mu3) + mu2*PsiT(uk - eta/mu2) + mu1*Hadj(nukp - xi/mu1);
    end

    vkp = real(ifftn(v_mult .* fftn(vkp_numerator)));

    %Update dual and parameter for Hs=v constraint
    Hvkp = Hfor(vkp);
    r_sv = Hvkp-nukp;
    xi = xi + mu1*r_sv;
    f.dual_resid_s(n) = gather(mu1*norm(vec(Hvk - Hvkp)));
    f.primal_resid_s(n) = gather(norm(vec(r_sv)));
    [mu1, mu1_update] = ADMM3D_update_param(mu1,solverSettings.resid_tol,solverSettings.mu_inc,solverSettings.mu_dec,f.primal_resid_s(n),f.dual_resid_s(n));

    % Update dual and parameter for Ls=v
    f.data_fidelity(n) = gather(.5*norm(crop3d(Hvkp)-b,'fro')^2);
    switch lower(solverSettings.regularizer)
        case('tv')
            Lvk1_ = Lvk1;
            Lvk2_ = Lvk2;
            Lvk3_ = Lvk3;
            [Lvk1, Lvk2, Lvk3] = Psi(vkp);
            r_su_1 = Lvk1 - uk1;
            r_su_2 = Lvk2 - uk2;
            r_su_3 = Lvk3 - uk3;
            eta_1 = eta_1 + mu2*r_su_1;
            eta_2 = eta_2 + mu2*r_su_2;
            eta_3 = eta_3 + mu2*r_su_3;
            f.dual_resid_u(n) = gather(mu2*sqrt(norm(vec(Lvk1_ - Lvk1))^2 + norm(vec(Lvk2_ - Lvk2))^2 + norm(vec(Lvk3_ - Lvk3))^2));
            f.primal_resid_u(n) = gather(sqrt(norm(vec(r_su_1))^2 + norm(vec(r_su_2))^2 + norm(vec(r_su_3))^2));
            f.regularizer_penalty(n) = gather(solverSettings.tau*(sum(vec(abs(Lvk1))) + sum(vec(abs(Lvk2))) + sum(vec(abs(Lvk3)))));

        case('tv_native')
            Lvk1_ = Lvk1;
            Lvk2_ = Lvk2;
            Lvk3_ = Lvk3;
            Lvk4_ = Lvk4;
            [Lvk1, Lvk2, Lvk3, Lvk4] = Psi(vkp);
            r_su_1 = Lvk1 - uk1;
            r_su_2 = Lvk2 - uk2;
            r_su_3 = Lvk3 - uk3;
            r_su_4 = Lvk4 - uk4;
            eta_1 = eta_1 + mu2*r_su_1;
            eta_2 = eta_2 + mu2*r_su_2;
            eta_3 = eta_3 + mu2*r_su_3;
            eta_4 = eta_4 + mu2*r_su_4;
            f.dual_resid_u(n) = gather(mu2*sqrt(norm(vec(Lvk1_ - Lvk1))^2 + norm(vec(Lvk2_ - Lvk2))^2 + ...
                norm(vec(Lvk3_ - Lvk3))^2 + norm(vec(Lvk4_ - Lvk4))^2));
            f.primal_resid_u(n) = gather(sqrt(norm(vec(r_su_1))^2 + norm(vec(r_su_2))^2 + ...
                norm(vec(r_su_3))^2 + norm(vec(r_su_4))^2));
           f.regularizer_penalty(n) = gather(solverSettings.tau*(sum(vec(abs(Lvk1))) +...
               sum(vec(abs(Lvk2))) + sum(vec(abs(Lvk3)))) + ...
               solverSettings.tau_n*sum(vec(abs(Lvk4))));
        case('native')
            Lvk_ = Lvk;
            Lvk = Psi(vkp);
            r_su = Lvk - uk;
            eta = eta + mu2*r_su;
            f.dual_resid_u(n) = gather(mu2*norm(vec(Lvk_ - Lvk)));
            f.primal_resid_u(n) = gather(norm(vec(r_su)));
            f.regularizer_penalty(n) = gather(solverSettings.tau_n*(sum(vec(abs(Lvk)))));
    end
    f.objective(n) = f.data_fidelity(n) + f.regularizer_penalty(n);

    [mu2, mu2_update] = ADMM3D_update_param(mu2,solverSettings.resid_tol,...
        solverSettings.mu_inc,solverSettings.mu_dec,...
        f.primal_resid_u(n),f.dual_resid_u(n));

    % Update nonnegativity dual and parameter (s=w)
    r_sw = vkp-wkp;
    rho = rho + mu3*r_sw;
    f.dual_resid_w(n) = gather(mu3*norm(vec(vk - vkp)));
    f.primal_resid_w(n) = gather(norm(vec(r_sw)));
    [mu3, mu3_update] = ADMM3D_update_param(mu3,solverSettings.resid_tol,solverSettings.mu_inc,solverSettings.mu_dec,f.primal_resid_w(n),f.dual_resid_w(n));

    %Update filters
    if mu1_update || mu2_update || mu3_update
        %fprintf('Mu updates: %i \t %i \t %i\n',mu1_update, mu2_update, mu3_update);
        mu_update = 1;
    else
        mu_update = 0;
    end
    if mu_update
        v_mult = 1./(mu1*HtH + mu2*PsiTPsi + mu3);  %This is the frequency space division fo S update
        nu_mult = 1./(DtD + mu1);
    end

    vk = vkp;

    if mod(n,solverSettings.save_every) == 0
        fprintf('saving state %i...\n',n)
        out_file = save_state(solverSettings,n);
        save(out_file,solverSettings.save_vars{:});
        fprintf('done saving\n')
    end

    if mod(n,solverSettings.print_interval) == 0
        t_iter = toc/solverSettings.print_interval;
         fprintf('iter: %i \t t: %.2g \t cost: %.2g \t data_fidelity: %.2g \t norm: %.2g \t Primal v: %.2g \t Dual v: %.2g \t Primal u: %.2g \t Dual u: %.2g \t Primal w: %.2g \t Dual w: %.2g \t mu1: %.2g \t mu2: %.2g \t mu3: %.2g \n',...
            n,t_iter,f.objective(n),f.data_fidelity(n),f.regularizer_penalty(n),f.primal_resid_s(n), f.dual_resid_s(n),f.primal_resid_u(n), f.dual_resid_u(n),f.primal_resid_w(n), f.dual_resid_w(n),mu1,mu2,mu3)
            %disp([n,f.objective(n),f.data_fidelity(n),f.regularizer_penalty(n),f.primal_resid_s(n), f.dual_resid_s(n),f.primal_resid_u(n), f.dual_resid_u(n),f.primal_resid_w(n), f.dual_resid_w(n),mu1,mu2,mu3])
        tic;
    end
    if mod(n,solverSettings.disp_figs) == 0
        draw_figures(vk,solverSettings)
    end
end
end

% Private function to display figures
function draw_figures(xk, solverSettings)
set(0,'CurrentFigure',solverSettings.fighandle)
if numel(size(xk))==2
    imagesc(solverSettings.display_func(xk))
    axis image
    colorbar
    colormap(solverSettings.color_map);

elseif numel(size(xk))==3
    xk = solverSettings.disp_crop(xk);
    subplot(1,3,1)

    im1 = squeeze(sum(xk,3));
    imagesc(solverSettings.display_func(im1));
    hold on
    axis image
    colormap (solverSettings.cmap)
    %colorbar
    caxis([0 prctile(im1(:),solverSettings.disp_percentile)])
    set(gca,'fontSize',6)
    axis off
    title('XY')
    hold off

    subplot(1,3,2)
    im2 = squeeze(max(xk,[],1));
    imagesc(im2);
    hold on
    %axis image
    colormap (solverSettings.cmap)
    %colorbar
    set(gca,'fontSize',8)
    caxis([0 prctile(im2(:),solverSettings.disp_percentile)])
    title('XZ')
    axis off
    hold off

    subplot(1,3,3)
    im3 = squeeze(max(xk,[],2));
    imagesc(solverSettings.disp_func(im3));
    hold on
    %axis image
    colormap (solverSettings.cmap)
    title('YZ')
    colorbar
    set(gca,'fontSize',8)
    caxis([0 prctile(im3(:),solverSettings.disp_percentile)]);
    axis off
    hold off

end

drawnow
end

function PsiTPsi = generate_laplacian(lapl)  %Takes in an array and makes laplacian on same grid (fft shifted)
    %lapl = zeros(2*Ny,2*Nx,Nz);    %Compute laplacian in closed form. This is the kernal to compute Psi'Psi
    lapl(1) = 6;
    lapl(1,2,1) = -1;
    lapl(2,1,1) = -1;
    lapl(1,1,2) = -1;
    lapl(1,end,1) = -1;
    lapl(end,1,1) = -1;
    lapl(1,1,end) = -1;
    PsiTPsi = abs(fftn(lapl));   %Compute power spectrum of laplacian
end

function [mu_out, mu_update] = ADMM3D_update_param(mu,resid_tol,mu_inc,mu_dec,r,s)
    if r > resid_tol*s
        mu_out = mu*mu_inc;
        mu_update = 1;
    elseif r*resid_tol < s
        mu_out = mu/mu_dec;
        mu_update = -1;
    else
        mu_out = mu;
        mu_update = 0;
    end
end

author: Waller-Lab

repo: DiffuserCam


text:
classdef HypergeometricD
%%HYPERGEOMETRICD Functions to handle the scalar hypergeometric
%                 distribution. Given an urn containg N1+N2 objects, N2 of
%                 type 1 and N2 of type 2, choose N balls. The number of
%                 balls of type 1 that is chosen is given by a
%                 hypergeometric distribution with parameters n, N1 and N2.
%Implemented methods are: mean, var, PMF, CDF, rand
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)
function val=mean(n,N1,N2)
%%MEAN Obtain the mean of the hypergeometric distribution.
%
%INPUTS: n The total number of objects drawn. n<=N1+N2.
%       N1 The number of objects of type 1.
%       N2 The number of objectes of type 2.
%
%OUTPUTS: val The mean of the hypergeometric distribution under
%             consideration.
%
%The hypergeometric distribution is given in Chapter 2.8 of [1].
%
%EXAMPLE:
%We verify that the computed variance by comparing it to the sample
%variance.
% n=20;
% N1=63;
% N2=28;
% varVal=HypergeometricD.mean(n,N1,N2)
% numSamp=1e3;
% sampVarVal=mean(HypergeometricD.rand([1,numSamp],n,N1,N2))
%One will find both values are about 13.836
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=n*N1/(N1+N2);
end

function val=var(n,N1,N2)
%%VAR Obtain the variance of the hypergeometric distribution.
%
%INPUTS: n The total number of objects drawn. n<=N1+N2.
%       N1 The number of objects of type 1.
%       N2 The number of objectes of type 2.
%
%OUTPUTS: val The variance of the hypergeometric distribution under
%             consideration.
%
%The hypergeometric distribution is given in Chapter 2.8 of [1].
%
%EXAMPLE:
%We verify that the computed variance by comparing it to the sample
%variance.
% n=20;
% N1=63;
% N2=28;
% varVal=HypergeometricD.var(n,N1,N2)
% numSamp=1e3;
% sampVarVal=var(HypergeometricD.rand([1,numSamp],n,N1,N2))
%One will find both values are about 3.36.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=n*N1*N2/(N1+N2)^2*(1-(n-1)/(N1+N2-1));
end

function val=PMF(x,n,N1,N2)
%%PMF Evaluate the hypergeometric probability mass function (PMF) at given
%     points.
%
%INPUTS: x The point(s) at which the hypergeometric PMF is evaluated.
%        n The total number of objects drawn. n<=N1+N2.
%        N1 The number of objects of type 1.
%        N2 The number of objectes of type 2.
%
%OUTPUTS: val The PMF of the hypergeometric distribution evaluated at the
%             desired point(s).
%
%The hypergeometric distribution is given in Chapter 2.8 of [1].
%
%EXAMPLE:
%Here, we validate the PMF by generating random samples and comparing the
%PMF plot with a histogram of the random samples.
% n=20;
% N1=25;
% N2=15;
% numSamples=1e5;
% figure(1)
% clf
% histogram(HypergeometricD.rand([numSamples,1],n,N1,N2),'BinWidth',1,'Normalization','pdf')
% hold on
% x=0:20;
% vals=HypergeometricD.PMF(x,n,N1,N2);
% stem(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(x));
    numEls=numel(x);

    denom=binomial(N1+N2,n);
    for curEl=1:numEls
        i=x(curEl);
        if(i>0&&i<=n)
            val(curEl)=binomial(N1,i).*binomial(N2,n-i)/denom;
        end
    end
end

function val=CDF(x,n,N1,N2)
%%CDF Evaluate the cumulative distribution function of the hypergeometric
%     distribution at desired points.
%
%INPUTS: x The point(s) at which the hypergeometric CDF is evaluated.
%        n The total number of objects drawn. n<=N1+N2.
%        N1 The number of objects of type 1.
%        N2 The number of objectes of type 2.
%
%OUTPUTS: val The CDF of the hypergeometric distribution evaluated at the
%             desired point(s).
%
%The hypergeometric distribution is given in Chapter 2.8 of [1].
%
%EXAMPLE:
%We validate the CDF value by comparing it to a value computed from random
%samples.
% n=20;
% N1=25;
% N2=15;
% x=12;
% numSamples=1e5;
% prob=HypergeometricD.CDF(x,n,N1,N2)
% probSamp=mean(HypergeometricD.rand([numSamples,1],n,N1,N2)<=x)
%One will find the values of both be about 0.5.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(x));
    numEls=numel(x);

    denom=binomial(N1+N2,n);
    for curEl=1:numEls
        i=floor(x(curEl));
        if(i>0&&i<=n)
            %Just sum the numerators in the PDF and normalize. This is
            %slow if N1+N2 is large.

            sumVal=0;
            for k=0:i
                sumVal=sumVal+binomial(N1,k).*binomial(N2,n-k);
            end
            val(curEl)=sumVal/denom;
        elseif(i>n)
           val(curEl)=1;
        end
    end
end

function x=rand(N,n,N1,N2)
%%RAND Generate hypergeometric random variables.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%        n The total number of objects drawn. n<=N1+N2.
%        N1 The number of objects of type 1.
%        N2 The number of objectes of type 2.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated hypergeometric random variables.
%
%The inverse transform method fo Chapter 4.1 of [1] is used.
%
%REFERENCES:
%[1] S. M. Ross, Simulation, 4th ed. Amsterdam: Academic Press, 2006.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N N];
    else
        dims=N;
    end

    denom=binomial(N1+N2,n);

    x=zeros(dims);
    numEls=prod(dims);

    for curEl=1:numEls
        U=rand();
        cumProb=0;
        for k=0:n
            cumProb=cumProb+binomial(N1,k).*binomial(N2,n-k)/denom;

            if(U<=cumProb)
                break;
            end
        end
        x(curEl)=k;
    end
end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [M_mean,n]=meanfilt3(M,IND,k,v)

%Creating spherical mask environment with radius k
[kI,kJ,kK]=cart2im(k,k,k,v);
kI=round(kI); kJ=round(kJ); kK=round(kK);
[MASK_J,MASK_I,MASK_K]=meshgrid(-kJ:kJ,-kI:kI,-kK:kK);
% [MASK_X,MASK_Y,MASK_Z]=im2cart(MASK_I,MASK_J,MASK_K,v);
MASK_X=MASK_J.*v(2); MASK_Y=MASK_I.*v(1); MASK_Z=MASK_K.*v(3);

R=sqrt(MASK_X.^2 + MASK_Y.^2 + MASK_Z.^2);

Lv=R<=k;
MASK_I=MASK_I(Lv); MASK_J=MASK_J(Lv); MASK_K=MASK_K(Lv);

%Getting mask indices
[IND_mask]=maskfind(M,IND(:),MASK_I(:),MASK_J(:),MASK_K(:));
L_valid=IND_mask>0;
INT_valid=M(IND_mask(L_valid));
INT_mask=nan(size(IND_mask));
INT_mask(L_valid)=INT_valid;

%Calculating median, ignoring NaN's
M_mean=gnanmean(INT_mask,2);

%Calculating number of elements used in median calculation
n=sum(~isnan(INT_mask),2);

%         %Creating spherical mask environment
%         k=k+iseven(k);
%         k_offset=round(k/2)-1;
%         [MASK_J,MASK_I,MASK_K]=meshgrid(-k_offset:k_offset);
%         R=sqrt(MASK_J.^2 + MASK_I.^2 + MASK_K.^2);
%         Lv=R<=k_offset;
%         MASK_I=MASK_I(Lv); MASK_J=MASK_J(Lv); MASK_K=MASK_K(Lv);
%
%         %Getting mask indices
%         [IND_mask]=maskfind(M,IND(:),MASK_I(:),MASK_J(:),MASK_K(:));
%         L_valid=IND_mask>0;
%         INT_valid=M(IND_mask(L_valid));
%         INT_mask=nan(size(IND_mask));
%         INT_mask(L_valid)=INT_valid;
%
%         %Calculating median, ignoring NaN's
%         M_median=nanmedian(INT_mask,2);
%
%         %Calculating number of elements used in median calculation
%         n=sum(~isnan(INT_mask),2);
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = prtUiPlot(varargin)
% obj = prtUiPlot(input, arguments, to, plot, ...)
%   obj is a prtUiPlotManager

obj = prtUiManagerPlot;
obj.plot(varargin{:});

if nargout
    varargout = {obj};
else
    varargout = {};
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% MAIN.m
%
% Solve the cart-pole swing-up problem  --  minimum time
%
% Note:  This problem is much more difficult to solve than the
% minimum-force version. This is because most of the control trajectory is
% sitting on a constraint: the maximum or minimum control force. This is
% generally true of minimum-time trajectories: they have bang-bang
% solutions. To get the exact solution, you would need to do many steps of
% mesh refinement. Here I only do two iterations, to keep total time
% reasonable. Another problem with minimum-time objective functions is that
% they sometimes have singular arcs: solutions where there is no single
% best control trajectory. This will manifest itself as "chattering" in the
% control trajectory and slow convergence. One solution is to include a
% regularization term, such as force squared with a very small coefficient,
% which forces a unique solution along the singular arc.
%

clc; clear;
addpath ../../

p.m1 = 2.0;  % (kg) Cart mass
p.m2 = 0.5;  % (kg) pole mass
p.g = 9.81;  % (m/s^2) gravity
p.l = 0.5;   % (m) pendulum (pole) length

dist = 1.0;  %How far must the cart translate during its swing-up
maxForce = 50;  %Maximum actuator forces

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Set up function handles                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.func.dynamics = @(t,x,u)( cartPoleDynamics(x,u,p) );
problem.func.pathObj = @(t,x,u)( ones(size(t)) );

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Set up problem bounds                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.bounds.initialTime.low = 0;
problem.bounds.initialTime.upp = 0;
problem.bounds.finalTime.low = 0.01;
problem.bounds.finalTime.upp = inf;

problem.bounds.initialState.low = zeros(4,1);
problem.bounds.initialState.upp = zeros(4,1);
problem.bounds.finalState.low = [dist;pi;0;0];
problem.bounds.finalState.upp = [dist;pi;0;0];

problem.bounds.state.low = [-2*dist;-2*pi;-inf;-inf];
problem.bounds.state.upp = [2*dist;2*pi;inf;inf];

problem.bounds.control.low = -maxForce;
problem.bounds.control.upp = maxForce;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Initial guess at trajectory                          %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.guess.time = [0,2];
problem.guess.state = [problem.bounds.initialState.low, problem.bounds.finalState.low];
problem.guess.control = [0,0];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Solver options                                  %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.options(1).nlpOpt = optimset(...
    'Display','iter',...
    'TolFun',1e-3,...
    'MaxFunEvals',1e5);
problem.options(1).method = 'trapezoid';
problem.options(1).trapezoid.nGrid = 10;

problem.options(2).nlpOpt = optimset(...
    'Display','iter',...
    'TolFun',1e-6,...
    'MaxFunEvals',1e5);
problem.options(2).method = 'trapezoid';
problem.options(2).trapezoid.nGrid = 30;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                            Solve!                                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

soln = optimTraj(problem);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Display Solution                                 %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%%%% Unpack the simulation
t = linspace(soln(end).grid.time(1), soln(end).grid.time(end), 150);
z = soln(end).interp.state(t);
u = soln(end).interp.control(t);

%%%% Plots:

%%%% Draw Trajectory:
[p1,p2] = cartPoleKinematics(z,p);

figure(2); clf;
nFrame = 9;  %Number of frames to draw
drawCartPoleTraj(t,p1,p2,nFrame);

%%%% Show the error in the collocation constraint between grid points:
%
if strcmp(soln(end).problem.options.method,'trapezoid') || strcmp(soln(end).problem.options.method,'hermiteSimpson')
    % Then we can plot an estimate of the error along the trajectory
    figure(5); clf;

    % NOTE: the following commands have only been implemented for the direct
    % collocation(trapezoid, hermiteSimpson) methods, and will not work for
    % chebyshev or rungeKutta methods.
    cc = soln(end).interp.collCst(t);

    subplot(2,2,1);
    plot(t,cc(1,:))
    title('Collocation Error:   dx/dt - f(t,x,u)')
    ylabel('d/dt cart position')

    subplot(2,2,3);
    plot(t,cc(2,:))
    xlabel('time')
    ylabel('d/dt pole angle')

    idx = 1:length(soln(end).info.error);
    subplot(2,2,2); hold on;
    plot(idx,soln(end).info.error(1,:),'ko');
    title('State Error')
    ylabel('cart position')

    subplot(2,2,4); hold on;
    plot(idx,soln(end).info.error(2,:),'ko');
    xlabel('segment index')
    ylabel('pole angle');
end

%%%% Plot the state and control against time
figure(1); clf;
plotPendulumCart(t,z,u,p);

author: MatthewPeterKelly

repo: OptimTraj


text:
function [x, fval, exitflag, output, lambda] = cplexlp (f, Aineq, bineq, Aeq, beq, lb, ub, x0, options)
%%
% Purpose
% Solve linear programming problems.
%
% Syntax
%    x = cplexlp(f, Aineq, bineq)
%    x = cplexlp(f, Aineq, bineq, Aeq, beq)
%    x = cplexlp(f, Aineq, bineq, Aeq, beq, lb, ub)
%    x = cplexlp(f, Aineq, bineq, Aeq, beq, lb, ub, x0)
%    x = cplexlp(f, Aineq, bineq, Aeq, beq, lb, ub, x0, options)
%    x = cplexlp(problem)
%    [x, fval] = cplexlp(...)
%    [x, fval, exitflag] = cplexlp(...)
%    [x, fval, exitflag, output] = cplexlp(...)
%    [x, fval, exitflag, output, lambda] = cplexlp(...)
%
% Description
% Finds the minimum of a problem specified by
%    min      f*x
%    st.      Aineq*x <= bineq
%             Aeq*x    = beq
%             lb <= x <= ub
%
% f, bineq, beq, lb, and ub are column vectors.
% Aineq and Aeq are matrices.
%
% x = cplexlp(f,Aineq,bineq) solves the linear programming problem min f*x
% such that Aineq*x <= bineq.
%
% x = cplexlp(f,Aineq,bineq,Aeq,beq) solves the preceding problem with the
% additional equality constraints Aeq*x = beq. If no inequalities exist,
% set Aineq=[] and bineq=[].
%
% x = cplexlp(f,Aineq,bineq,Aeq,beq,lb,ub) defines a set of lower and upper
% bounds on the design variables, x, so that the solution is always in the
% range lb <= x <= ub. If no equalities exist, set Aeq=[] and beq=[].
%
% x = cplexlp(f,Aineq,bineq,Aeq,beq,lb,ub,x0) sets the starting point for
% the algorithm to x0. If no bounds exist, set lb=[] and ub=[].
%
% x = cplexlp(f,Aineq,bineq,Aeq,beq,lb,ub,x0,options) minimizes with the
% default optimization options replaced by values in the structure options,
% which can be created using the function cplexoptimset. If you do not wish
% to give an initial point, set x0=[].
%
% x = cplexlp(problem) where problem is a structure.
%
% [x,fval] = cplexlp(...) returns the value of the objective function at
% the solution x: fval = f*x.
%
% [x,fval,exitflag] = cplexlp(...) returns a value exitflag that describes
% the exit condition of cplexlp.
%
% [x,fval,exitflag,output] = cplexlp(...) returns a structure output that
% contains information about the optimization.
%
% [x,fval,exitflag,output,lambda] = cplexlp(...) returns a structure lambda
% whose fields contain the Lagrange multipliers at the solution x.
%
% Input Arguments
% problem   Structure containing the following fields:
%           f         Double column vector for objective function
%           Aineq     Double matrix for linear inequality constraints
%           bineq     Double column vector for linear inequality
%                     constraints
%           Aeq       Double matrix for linear equality constraints
%           beq       Double column vector for linear equality constraints
%           lb        Double column vector of lower bounds
%           ub        Double column vector of upper bounds
%           x0        Double column vector of initial point of x
%           options   Options structure created with cplexoptimset
%
% Output Arguments
% x         Solution found by the optimization function. If exitflag > 0,
%           then x is a solution; otherwise, x is the value of the
%           optimization routine when it terminated prematurely.
% fval      Value of the objective function at the solution x
% exitflag  Integer identifying the reason the optimization algorithm
%           terminated
% output    Structure containing information about the optimization. The
%           fields of the structure are:
%           iterations         Number of iterations
%           algorithm          Optimization algorithm used
%           message            Exit message
%           time               Execution time of the algorithm
%           cplexstatus        Status code of the solution
%           cplexstatusstring  Status string of the solution
% lambda    Structure containing the Lagrange multipliers at the solution x
%           (separated by constraint type). The fields of the structure
%           are:
%           lower              Lower bounds lb
%           upper              Upper bounds ub
%           ineqlin            Linear inequalities
%           eqlin              Linear equalities
%
%
%  See also cplexoptimset
%

% ---------------------------------------------------------------------------
% File: cplexlp.m
% ---------------------------------------------------------------------------
% Licensed Materials - Property of IBM
% 5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55
% Copyright IBM Corporation 2008, 2010. All Rights Reserved.
%
% US Government Users Restricted Rights - Use, duplication or
% disclosure restricted by GSA ADP Schedule Contract with
% IBM Corp.
% ---------------------------------------------------------------------------

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function d = ddiff ( d1, d2 )

%*****************************************************************************80
%
%% DDIFF returns the signed distance to a region that is the difference of two regions.
%
%  Discussion:
%
%    The author comments that this is not the true signed distance function for
%    the difference set, in particular, around corners.
%
%  Copyright:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real D1, D2, the signed distances to region 1 and 2.
%
%    Output, real D, the signed distance to the region formed by
%    removing from region 1 its intersection with region 2.
%
  d = max ( d1, -d2 );

  return
end


text:
%% Stream programming in Matlab
%  The stream programming is a computer programming methodology that is
% born in the field of mathematical languages using lambda calculus toole
% and exploiting the lazy evaluation methodology. Here a unefficient
% library, really useful to explain lazy evaluation and stream programming,
% is given. The set of tools is really wide to give examples of pure lazy
% evaluation, simple stream generation, the composition pattern, the
% accumulation pattern, the filtering pattern and a wide set of examples
% that allow the user to learn what stream programming is. This library
% isn't useful to create systems that must work with large ammounts of
% data, but can used where the number of data elements is small.. for
% example to generate streams of images (i.e. a stream of the
% approximations of a given image using a set of fourier armonics, this is
% a classical accumulation example).

%% What lazy evaluation is
%  The lazy evaluation is the technique of managing entities that aren't
% computed yet. For example we can think to the gradient of an image that
% can be taken as a variable and passed in functions and in workspaces, but
% that isn't computed yet, when we acces, for example, to a tail of the
% gradient that isn't computed automatically the lazy evaluation system
% computes that part of the matrix. This is usefull to optimize the code
% for execution time (parts of the gradient aren't computed) and for memory
% usage (because the parts of the gradient that aren't used arent allocated
% in the gradient matrix).
%
%  For example we can delay the evaluation of a function with a set of
% parameters:

% Function evaluation delayed:
dsin = delayEval(@sin,{2}),

%% Forcing the evaluation
%  After the delayed evaluation composition, the entity "dsin" can be thing
% to the promise of evaluation of the function "sin" with parameter "2".
% This isn't done yet, we can force the evaluation when we wish to get it.
% If we have tools to guess if an entity must be evaluated or not we can
% automatize the evaluation so that it is done only when an entity is
% accessed the first time, like the example of the gradient.
%
%  Here we force the evaluation exlpicitly:

% Forcing the evaluation:
res = forceEval(dsin),

%% The Stream Programming (SP)
%  The stream programming is a technique, based on lazy evaluation, to
% manage infinitely long sequence of elements, for examples the sequences
% of numbers, or the sequences of functions, series or other things. To
% manage an infinitely long sequence, that cannot be stored in a finite
% computer memory, we think to the sequence as a first element  followed
% with the promise to compute the other part of the sequence: this is a
% stream:
%
%     % A sample stream generating function generating the integers sequence:
%     function s = integersFrom(n)
%     s = {n, delayEval(@integersFrom,{n+1})};

% Generating an integers sequence:
integers = integersFrom(0),

% Getting another integer:
otherIntegers = forceEval(integers{2}),

% Computing the first ten integers:
firstTenInts = streamHead(integers,10),

%% Fibonacci numbers
%  The Fibonacci numbers can be constructed as a sequence where each number
% can be constructed as the sum of the previous two numbers. The first two
% numbers ar 0 and 1. With the following stream constructor the Fibonacci
% sequence can be constructed:
%
%     % The constructor of the Fibonacci sequence:
%     function s = fiboFrom(a1,a2)
%     s = {a1,delayEval(@fiboFrom,{a2,a1+a2})};
%
%  using this constructor we can obtain:

% Creating the sequence:
fib = fiboFrom(0,1);

% Getting the first ten values:
firstTenInts = streamHead(fib,10),

%% Functional list operators (from Schema language)
%  In the lisp and Schema languages the car and cdr operators are the bases
% for the management of lists. Here we have used "streamHead" to get the
% first elements of a stream, that is an infinite list; functions "head"
% and "tail" are defined as aliases of "streamCar" and "streamCdr", this
% two functions can be easily defined as follows:
%
%     % Head of a stream:
%     function v = streamCar(s)
%     v = s{1};
%
%     % Tail of a stream:
%     function so = streamCdr(si)
%     so = forceEval(si{2});
%
%  the streamHead function simply iterates with this two functions:
%
%     % Iteratively obtain a sequence of values:
%     function vals = streamHead(s,n)
%     s1 = s;
%     vals = cell(1,n);
%     for i=1:n
%         vals{i} = streamCar(s1);
%         s1 = streamCdr(s1);
%     end
%
%  other helper functions can be written with the same philosophy, like
% "streamRange" or "streamCons" to concatenate a value to a stream:
%
%     % Cut a range of a stream:
%     function vals = streamRange(s,start,nElem)
%     s1 = s;
%     for i=1:start-1
%         s1 = streamCdr(s1);
%     end
%     vals = cell(1,nElem);
%     for i=1:nElem
%         vals{i} = streamCar(s1);
%         s1 = streamCdr(s1);
%     end
%
%     % Concatenate a value with a stream:
%     function so = streamCons(v,si)
%     so = {v, delayEval(@fIdentity,{si})};

%% SP patterns: accumulate
%  In the software engeneering the design patterns are fundamental bricks
% for a correct and robust design of a complex system; generally desogn
% patterns are defined for OOP, but the same can be done for functional
% programming, procedural programming... and SP.
%  The accumulate pattern allows to define elements that operates over
% streams and that can be used as accumulators of data. Accumulators can be
% used to do a wide number of works like integration, sum of series...
% and for example the generation of streams containing successive
% approximations of series of what we wish to manage, for example the
% approximations of an image using its Fourier coefficients. Here the
% accumulator tool is defined as follows:
%
%     % Given a stream, a starting value and an accumulator function an
%     % accumulated stream is generated with the starting value as first
%     % element and the accumulator function that computes the next element.
%     function so = streamAccumulate(si,start,accumulator)
%     so = {start,delayEval(@streamAccumulate, ...
%             {streamCdr(si),
%              feval(accumulator,start,streamCar(si)),accumulator})};
%
%  Using this simple function (but take a deep look to it) we can define
% easily the integration:
%
%     % Integration of a sequence of numbers:
%     function so = integrateStream(si,dt,c)
%     so = scaleStream(streamAccumulate(si,c/dt,@plus),dt);
%
%  Let's try it:

% A discretization step:
dt = 0.1;

% A sequence of values from 0:
X = scaleStream(integers,dt);

% A sine function:
Y = filterStream(X,@sin);

% Getting sample values:
Xs = streamHead(X,10),
Ys = streamHead(Y,10),

% Computing the integral stream:
Yint = integrateStream(Y,dt,-1);
Yints = streamHead(Yint,10),
figure; hold on;
plotStream(X,Y,100,0),
plotStream(X,Yint,100,0),

%% SP patterns: accumulate
%  The composition of streams can be used to do different computations in
% different streams and then mixing them int a single stream. A simple
% example is the stream sum of other streams, we can try to sum the sin
% stream and it's integral.
%  A composition can be done given a sequence of streams and a composing
% function as follows:
%
%     % Composing n streams:
%     function so = streamCompose(streams,func)
%     argsv = {func};
%     argsf = {};
%     for i=1:size(streams,2);
%         argsv{i+1} = streamCar(streams{i});
%         argsf{i} = streamCdr(streams{i});
%     end
%     val = builtin('feval',argsv{:});
%     so = {val,delayEval(@streamCompose,{argsf,func})};
%
%  This tool allows to generate easily new composition operators:
%
%     % Adding the streams values:
%     function so = addStreams(varargin)
%     so = streamCompose(varargin,@sumVals);
%     function ris = sumVals(varargin)
%     ris = sum(cell2mat(varargin));
%
%     % Subtracting the stream values:
%     function so = subStreams(s1,s2)
%     so = streamCompose({s1,s2},@subVals);
%     function ris = subVals(v1,v2)
%     ris = v1-v2;
%
%     % Multiply streams:
%     function so = mulStreams(varargin)
%     so = streamCompose(varargin,@mulVals);
%     function ris = mulVals(varargin)
%     ris = prod(cell2mat(varargin));
%
%     % Divide streams:
%     function so = divStreams(s1,s2)
%     so = streamCompose({s1,s2},@divVals);
%     function ris = divVals(v1,v2)
%     ris = v1/v2;
%
%  Whit this tool all the composition of streams that we wish to do can be
% easily added to our operators set.

% Generating the sum of the two streams:
Ysum = addStreams(Y,Yint);

% Plotting:
figure; plotStream(X,Ysum,100,0),

%% SP patterns: filter
%  A stream can be filtered in a sense different from the previous one
% (that is the mapping of functional programming), a stream can be filtered
% with a filter that decides if an element must be discarded or not. For
% example immagine that you wish to get a sequence of prime numbers, you
% can get them filtering the stream of integers from 2 removing the
% integers that can be divided by the previous primes. We can design a
% stream filter function that recives another boolean function that can
% decide if an element can be mantained or must be rejected:
%
%     % A filtering element for streams:
%     function so = streamFilter(si,filter,varargin)
%     while ~builtin('feval',filter,head(si),varargin{:})
%         si = tail(si);
%     end
%     so = {head(si),delayEval(@streamFilter,{tail(si),filter,varargin{:}})};
%
%  This can filter element by element a stream (and can generate infinite
% loops). This can be usefull to generate streams from other discarding
% elements; the prime numbers stream is an extreme case because the filters
% chain grows of one for each element (prime) generated:
%
%     % Generation of all the primes:
%     function so = primes
%     so = PrimesFrom(integersFrom(2));
%
%     % The filtering function add a new filter to the stream:
%     function so = PrimesFrom(si)
%     % The first number is ok:
%     val=head(si); so=tail(si);
%     % The others are filtered:
%     so = streamFilter(so,@(v,p)(mod(v,p)~=0),val);
%     % Composing:
%     so = {val,delayEval(@PrimesFrom,{so})};
%
%  This is really inefficient, also because the use of the recursion is
% heavy, but it's an approach that generates the required results:

% Creating the primes stream:
pr = primes;

% Extracting values:
streamHead(pr,20),

%% Conclusions
%  The stream programming is a powerful technique to manage infinite
% sequences of elements, the lazy evaluation allows to manage this
% particular entities easily, this implementation of streams in Matlab does
% not look at computational complexity but an efficient version can be
% implemented. Tools like Simulink, Labview, HP-VEE gives the stream
% programming as a base programming tool and a graphical environment to
% generate programs connecting stream generators, operators and consumers.
% The stream programming can be an useful tool also in programs written in
% other procedural languages like Python or Java, the power is big.. but
% the diffusion of this technique is not :(


text:
%% Physics-informed dynamic mode decompositions
%
% Computes a dynamic mode decomposition when the solution matrix is
% constrained to lie in a matrix manifold.
% The options available for the "method" so far are
%
% - "exact", "exactSVDS"
% - "orthogonal"
% - "uppertriangular, "lowertriangular"
% - "diagonal", "diagonalpinv", "diagonaltls", "symtridiagonal"
% - "circulant", "circulantTLS", "circulantunitary", "circulantsymmetric",
% "circulantskewsymmetric"
% - "BCCB", "BCCBtls", "BCCBskewsymmetric",
% "BCCBunitary", "hankel", "toeplitz"
% - "symmetric", "skewsymmetric"
%
function [A, varargout] = piDMD(X,Y,method,varargin)

[nx, nt] = size(X);

if strcmp(method,'exact') || strcmp(method,'exactSVDS')

    if nargin>3
        r = varargin{1};
    else
        r = min(nx,nt);
    end

    if strcmp(method,'exact')
        [Ux,Sx,Vx] = svd(X,0);
        Ux = Ux(:,1:r); Sx = Sx(1:r,1:r); Vx = Vx(:,1:r);
    elseif strcmp(method,'exactSVDS')
        [Ux,Sx,Vx] = svds(X,r);
    end

    Atilde = (Ux'*Y)*Vx*pinv(Sx);
    A = @(v) Ux*(Atilde*(Ux'*v));

    if nargout==2
    varargout{1} = eig(Atilde);
    elseif nargout>2
    [eVecs,eVals] = eig(Atilde);
    eVals = diag(eVals);
    eVecs = Y*Vx*pinv(Sx)*eVecs./eVals.';
    varargout{1} = eVals;
    varargout{2} = eVecs;
    end

elseif strcmp(method,'orthogonal')

    if nargin>3
        r = varargin{1};
    else
        r = min(nx,nt);
    end
    [Ux,~,~] = svd(X,0); Ux = Ux(:,1:r);
    Yproj = Ux'*Y; Xproj = Ux'*X; % Project X and Y onto principal components
    [Uyx, ~, Vyx] = svd(Yproj*Xproj',0);
    Aproj = Uyx*Vyx';
    A = @(x) Ux*(Aproj*(Ux'*x));

    if nargout==2
        eVals = eig(Aproj);
        varargout{1} = eVals;
    elseif nargout>2
        [eVecs,eVals] = eig(Aproj);
        varargout{1} = diag(eVals);
        varargout{2} = Ux*eVecs;
        if nargout > 3; varargout{3} = Aproj; end
    end

elseif strcmp(method,'uppertriangular')

    [R,Q] = rq(X); % Q*Q' = I
    Ut = triu(Y*Q');
    A = Ut/R;

elseif strcmp(method,'lowertriangular')

    A = rot90(piDMD(flipud(X),flipud(Y),'uppertriangular'),2);

% The codes allows for matrices of variable banded width. The fourth input,
% a 2xn matrix called d, specifies the upper and lower bounds of the
% indices of the non-zero elements. The first column corresponds to the width of
% the band below the diagonal and the second column is the width of the
% band above. For example, a diagonal matrix would have d = ones(nx,2) and
% a tridiagonal matrix would have d = [2 2]+zeros(nx,2). If you only specify
% d as a scalar then the algorithm converts the input to obtain a banded
% diagonal matrix of width d.
elseif startsWith(method,'diagonal')
    if nargin>3
        d = varargin{1}; % arrange d into an nx-by-2 matrix
        if numel(d) == 1
            d = d*ones(nx,2);
        elseif numel(d) == nx
             d = repmat(d,[1,2]);
        elseif any(size(d)~=[nx,2])
            error('Diagonal number is not in an allowable format.')
        end
    else
        d = ones(nx,2); % default is for a diagonal matrix
    end
    % Allocate cells to build sparse matrix
    Icell = cell(1,nx); Jcell = cell(1,nx); Rcell = cell(1,nx);
    for j = 1:nx
    l1 = max(j-(d(j,1)-1),1); l2 = min(j+(d(j,2)-1),nx);
    C = X(l1:l2,:); b = Y(j,:); % preparing to solve min||Cx-b|| along each row
    if strcmp(method,'diagonal')
            sol = b/C;
    elseif strcmp(method,'diagonalpinv')
            sol = b*pinv(C);
    elseif strcmp(method,'diagonaltls')
            sol = tls(C.',b.').';
    end
    Icell{j} = j*ones(1,1+l2-l1); Jcell{j} = l1:l2; Rcell{j} = sol;
    end
    Imat = cell2mat(Icell); Jmat = cell2mat(Jcell); Rmat = cell2mat(Rcell);
    Asparse = sparse(Imat,Jmat,Rmat,nx,nx);
    A = @(v) Asparse*v;

    if nargout==2
        eVals = eigs(Asparse,nx);
        varargout{1} = eVals;
    elseif nargout>2
        [eVecs, eVals] = eigs(Asparse,nx);
        varargout{1} = diag(eVals); varargout{2} = eVecs;
    end

elseif strcmp(method,'symmetric') || strcmp(method,'skewsymmetric')

[Ux,S,V] = svd(X,0);
C = Ux'*Y*V;
C1 = C;
if nargin>3; r = varargin{1}; else; r = rank(X); end
Ux = Ux(:,1:r);
Yf = zeros(r);
    if strcmp(method,'symmetric')
    for i = 1:r
        Yf(i,i) = real(C1(i,i))/S(i,i);
        for j = i+1:r
            Yf(i,j) = (S(i,i)*conj(C1(j,i)) + S(j,j)*C1(i,j)) / (S(i,i)^2 + S(j,j)^2);
        end
    end
    Yf = Yf + Yf' - diag(diag(real(Yf)));
    elseif strcmp(method,'skewsymmetric')
    for i = 1:r
        Yf(i,i) = 1i*imag(C1(i,i))/S(i,i);
        for j = i+1:nx
            Yf(i,j) = (-S(i,i)*conj(C1(j,i)) + S(j,j)*(C1(i,j))) / (S(i,i)^2 + S(j,j)^2);
        end
    end
    Yf = Yf - Yf' - 1i*diag(diag(imag(Yf)));
    end

    A = @(v) Ux*Yf*(Ux'*v);

    if nargout==2
        varargout{1} = eig(Yf);
    elseif nargout>2
        [eVecs,eVals] = eig(Yf);
        eVals = diag(eVals);
        eVecs = Ux*eVecs;
        varargout{1} = eVals;
        varargout{2} = eVecs;
    end

elseif strcmp(method,'toeplitz') || strcmp(method,'hankel')
   if  strcmp(method,'toeplitz'); J = eye(nx);
   elseif strcmp(method,'hankel'); J = fliplr(eye(nx)); end
    Am = fft([eye(nx) zeros(nx)].',[],1)'/sqrt(2*nx); % Define the left matrix
    B = fft([(J*X)' zeros(nt,nx)].',[],1)'/sqrt(2*nx); % Define the right matrix
    BtB = B'*B;
    AAt = ifft(fft([eye(nx) zeros(nx); zeros(nx,2*nx)]).').'; % Fast computation of A*A'
    y = diag(Am'*conj(Y)*B)'; % Construct the RHS of the linear system
    L = (AAt.*BtB.')'; % Construct the matrix for the linear system
    d = [y(1:end-1)/L(1:end-1,1:end-1) 0]; % Solve the linear system
    newA = ifft(fft(diag(d)).').'; % Convert the eigenvalues into the circulant matrix
    A = newA(1:nx,1:nx)*J; % Extract the Toeplitz matrix from the circulant matrix

elseif startsWith(method,'circulant')

 fX = fft(X); fY = fft(conj(Y));

 d = zeros(nx,1);
    if endsWith(method,'TLS') % Solve in the total least squares sense
        for j = 1:nx
        d(j) = tls(fX(j,:)',fY(j,:)');
        end
    elseif ~endsWith(method,'TLS') % Solve the other cases
    d = diag(fX*fY')./vecnorm(fX,2,2).^2;
        if endsWith(method,'unitary'); d = exp(1i*angle(d));
        elseif endsWith(method,'symmetric'); d = real(d);
        elseif endsWith(method,'skewsymmetric'); d = 1i*imag(d);
        end
    end
  eVals = d; % These are the eigenvalues
  eVecs = fft(eye(nx)); % These are the eigenvectors
 if nargin>3
    r = varargin{1}; % Rank constraint
    res = diag(abs(fX*fY'))./vecnorm(fX')'; % Identify least important eigenvalues
    [~,idx] = mink(res,nx-r); % Remove least important eigenvalues
    d(idx) = 0; eVals(idx) = []; eVecs(:,idx) = [];
 end

 if nargout>1; varargout{1} = eVals; end
 if nargout>2; varargout{2} = eVecs; end

 A = @(v) fft(d.*ifft(v)); % Reconstruct the operator in terms of FFTs

elseif strcmp(method,'BCCB') || strcmp(method,'BCCBtls') || strcmp(method,'BCCBskewsymmetric') || strcmp(method,'BCCBunitary')

    if isempty(varargin); error('Need to specify size of blocks.'); end
    s = varargin{1}; p = prod(s);
    % Equivalent to applying the block-DFT matrix F
    % defined by F = kron(dftmtx(M),dftmtx(N)) to the
    % matrix X
    aF =  @(x) reshape(     fft2(reshape(x ,[s,size(x,2)])) ,[p,size(x,2)])/sqrt(p);
    aFt = @(x) conj(aF(conj(x)));
    fX = aF(conj(X)); fY = aF(conj(Y));
    d = zeros(p,1);

    if strcmp(method,'BCCB')
    for j = 1:p; d(j) = conj(fX(j,:)*fY(j,:)')/norm(fX(j,:)').^2; end
    elseif strcmp(method,'BCCBtls')
    for j = 1:p; d(j) = tls(fX(j,:)',fY(j,:)')'; end
    elseif strcmp(method,'BCCBskewsymmetric')
    for j = 1:p; d(j) = 1i*imag(fY(j,:)/fX(j,:)); end
    elseif strcmp(method,'BCCBsymmetric')
    for j = 1:p; d(j) = real(fY(j,:)/fX(j,:)); end
    elseif strcmp(method,'BCCBunitary')
    for j = 1:p; d(j) = exp(1i*angle(fY(j,:)/fX(j,:))); end
    end

    % Returns a function handle that applies A
     if nargin>4
        r = varargin{2};
        res = diag(abs(fX*fY'))./vecnorm(fX')';
        [~,idx] = mink(res,nx-r);
        d(idx) = 0;
    end
    A = @(x) aF((conj(d).*aFt(x)));
    varargout{1} = d;
    % Eigenvalues are given by d

elseif strcmp(method,'BC') || strcmp(method,'BCtri') || strcmp(method,'BCtls')

    s = varargin{1}; p = prod(s);
        M = s(2); N = s(1);
    if isempty(s); error('Need to specify size of blocks.'); end
    % Equivalent to applying the block-DFT matrix F
    % defined by F = kron(dftmtx(M),eye(N)) to the
    % matrix X
aF  =  @(x) reshape(fft(reshape(x,[s,size(x,2)]),[],2) ,[p,size(x,2)])/sqrt(M);
aFt =  @(x) conj(aF(conj(x)));

fX = aF(X); fY = aF(Y);
    d = cell(M,1);

for j = 1:M
    ls = (j-1)*N + (1:N);
    if strcmp(method,'BC')
        d{j} = fY(ls,:)/fX(ls,:);
    elseif strcmp(method,'BCtri')
        d{j} = piDMD(fX(ls,:),fY(ls,:),'diagonal',2);
    elseif strcmp(method,'BCtls')
        d{j} = tls(fX(ls,:)',fY(ls,:)')';
    end
end

    BD = blkdiag(d{:});
    A = @(v) aFt(BD*aF(v));

elseif strcmp(method,'symtridiagonal')

    T1e = vecnorm(X,2,2).^2; % Compute the entries of the first block
    T1 = spdiags(T1e,0,nx,nx); % Form the leading block
    T2e = dot(X(2:end,:),X(1:end-1,:),2); % Compute the entries of the second block
    T2 = spdiags([T2e T2e],-1:0,nx,nx-1); % Form the second and third blocks
    T3e = [0; dot(X(3:end,:),X(1:end-2,:),2)]; % Compute the entries of the final block
    T3 = spdiags(T1e(1:end-1) + T1e(2:end),0,nx-1,nx-1) ...
         + spdiags(T3e,1,nx-1,nx-1) + spdiags(T3e,1,nx-1,nx-1)'; % Form the final block
    T = [T1 T2; T2' T3]; % Form the block tridiagonal matrix
    d = [dot(X,Y,2); dot(X(1:end-1,:),Y(2:end,:),2) + dot(X(2:end,:),Y(1:end-1,:),2)]; % Compute the RHS vector
    c = real(T)\real(d); % Take real parts then solve linear system
    % Form the solution matrix
    A = spdiags(c(1:nx),0,nx,nx) + spdiags([0;c(nx+1:end)],1,nx,nx) + spdiags([c(nx+1:end); 0],-1,nx,nx);
else
    error('The selected method doesn''t exist.');

end

author: baddoo

repo: piDMD

save_path: github-repos/MATLAB/baddoo-piDMD


text:
fileII = 'test2d_triangle';
b.fileName = fileII;
d = FemDataContainer(b);
femII = FEM.create(d);
femII.solve();
femII.print(fileII);

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [neighborIds neighborDistances] = kNearestNeighbors(dataMatrix, queryMatrix, k)
%--------------------------------------------------------------------------
% Program to find the k - nearest neighbors (kNN) within a set of points.
% Distance metric used: Euclidean distance
%
% Usage:
% [neighbors distances] = kNearestNeighbors(dataMatrix, queryMatrix, k);
% dataMatrix  (N x D) - N vectors with dimensionality D (within which we search for the nearest neighbors)
% queryMatrix (M x D) - M query vectors with dimensionality D
% k           (1 x 1) - Number of nearest neighbors desired
%
% Example:
% a = [1 1; 2 2; 3 2; 4 4; 5 6];
% b = [1 1; 2 1; 6 2];
% [neighbors distances] = kNearestNeighbors(a,b,2);
%
% Output:
% neighbors =
%      1     2
%      1     2
%      4     3
%
% distances =
%          0    1.4142
%     1.0000    1.0000
%     2.8284    3.0000
%--------------------------------------------------------------------------

neighborIds = zeros(size(queryMatrix,1),k);
neighborDistances = neighborIds;

numDataVectors = size(dataMatrix,1);
numQueryVectors = size(queryMatrix,1);
for i=1:numQueryVectors,
    dist = sum((repmat(queryMatrix(i,:),numDataVectors,1)-dataMatrix).^2,2);
    [sortval sortpos] = sort(dist,'ascend');
    neighborIds(i,:) = sortpos(1:k);
    neighborDistances(i,:) = sqrt(sortval(1:k));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB