text:
clear; clc; close all; warning off;
addpath(genpath('./.'));
addpath(genpath('/home/./caffe/')) ;
caffe.set_mode_gpu();
caffe.set_device(0);

scale      = 4;
ext        =  {'*.JPG','*.png','*.bmp', '*.tif'};
CameraType = {'Canon', 'Nikon'};
folder     = './.';
filepaths  = dir(fullfile(folder, '*.caffemodel'));
weights = fullfile(folder,filepaths.name);
model = '*.prototxt';
net = caffe.Net(model, weights,'test');

count = 1;
for q = 1:2
    folder2   =  ['Test/', CameraType{q}, '/', num2str(scale)];
    filepaths2 =  [];
    for p = 1 : length(ext)
        filepaths2 = cat(1,filepaths2, dir(fullfile(folder2, ext{p})));
    end
    for i = 1 : 2: length(filepaths2)
        disp(i)
        HR = imread(fullfile(folder2,filepaths2(i).name));
        HR = modcrop(HR, 4);
        HR_Ycbcr = im2double(rgb2ycbcr(HR));
        HR_Y = im2single(HR_Ycbcr(:, :, 1));

        LR = imread(fullfile(folder2,filepaths2(i+1).name));
        LR = modcrop(LR, 4);
        LR_Ycbcr = im2double(rgb2ycbcr(LR));
        LR_Y = im2single(LR_Ycbcr(:, :, 1));

        % The testing image is too large and will cause out-of-memory issue.
        size_img = size(LR_Y);
        if ((size_img(1) > 1200) && (size_img(2) > 1200))
            size_patch = [1200 1200];
            size_skip = [800 800];
            [patch] = im2patch(LR_Y, size_patch, size_skip);
            [patchY] = im2patch(HR_Y, size_patch, size_skip);
        else
            patch = LR_Y;
            patchY = HR_Y;
        end

        for k = 1:size(patch,3)
            net.blobs('data').reshape([size(patch(:,:,k),1) size(patch(:,:,k),2), 1, 1]);
            net.reshape();
            res = net.forward({patch(:,:,k)});
            result = res{1};

            [PSNRCur, SSIMCur] = Cal_PSNRSSIM(im2uint8(patchY(:,:,k)),im2uint8(result),0,0);
            PSNRs(count) = PSNRCur;
            SSIMs(count) = SSIMCur;
            count = count + 1;
        end
    end
end
mean(PSNRs)
mean(SSIMs)
caffe.reset_all()

author: csjcai

repo: RealSR

save_path: github-repos/MATLAB/csjcai-RealSR


text:
function [transform] = pos2transform(pos, dim)

% POS2TRANSFORM reconstructs a transformation matrix from an ordered list
% of positions.
%
% Use as
%   [transform] = pos2transform(pos, dim)
% where pos is an ordered list of positions that should specify a full 3D volume.
%
% The output transform is a 4x4 homogenous transformation matrix which transforms
% from 'voxelspace' into the positions provided in the input
%
% See also POS2DIM

% Copyright (C) 2009, Jan-Mathijs Schoffelen

if nargin>1
  % do nothing
else
  dim = pos2dim(pos);
end
x   = 1:dim(1);
y   = 1:dim(2);
z   = 1:dim(3);
[X,Y,Z] = ndgrid(x, y, z);
ind = [X(:) Y(:) Z(:)];
ind = ind'; ind(4,:) = 1;
pos = pos'; pos(4,:) = 1;

% build in some robustness against nans
sel = sum(isfinite(pos))==4;

transform = pos(:,sel)/ind(:,sel);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef diskfun < separableApprox
%DISKFUN   DISKFUN class for representing functions on the unit disk.
%
%   Class for approximating smooth functions defined on the unit disk. The
%   functions should be smooth. Diskfun objects can be constructed in a
%   variety of ways:
%
%    1. Y = diskfun(F), where F is a function handle in Cartesian
%       coordinates (x,y), e.g., @(x,y) x.*y + cos(x).
%    2. Y = diskfun(F, 'polar'), where F is a function handle in polar
%       coordinates (theta,r). Here, theta is the angular variable satisfying
%       -pi <= theta <= pi, and r is the radial variable satisfying 0 <= r < 1,
%       e.g., @(theta,r) cos(r.*sin(theta))
%    3. Y = diskfun(F), where F is a matrix of numbers.
%
% If F is a function handle then it should allow for vectorized evaluations.
%
% If F is a matrix, F = (f_ij), the numbers fij are used as function values
% at the tensor Fourier-Chebyshev points on the rectangle [-pi,pi]x[0,1].
%
% DISKFUN(F, k) returns a rank k approximation to F.
%
% DISKFUN(F, [m n]) returns a representation of F using a degree m
% Chebyshev approximation of F in the radial direction and a degree n
% trigonometric approximation in the angular direction. The result is
% compressed in low rank form and the rank k is still determined adaptively
% (satisfying k<=min(m,n)+1).
%
% The DISKFUN software system is based on:
%
% H. Wilber, A. Townsend, and G. Wright, Computing with functions in
% spherical and polar geometries II: The disk, SIAM. J. Sci. Comput., 39-4
% (2017), C238-C262.
%
% See also CHEBFUN2, SPHEREFUN, DISKFUNV.

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

% TODO: Include documentation of fixed eps construction

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function f = diskfun(varargin)
            % The main diskfun constructor!

            % Return an empty DISKFUN:
            if ( (nargin == 0) || isempty(varargin{1}) )
                f.domain = [-pi pi 0 1];
                return
            end
            % Call the constructor, all the work is done here:
            f = constructor(f, varargin{:});
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% HIDDEN METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false, Hidden = true )
        % Make f have BMC-II symmetry.
        f = projectOntoBMCII(f);

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% PUBLIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        % The main bulk of the DISKFUN constructor:
        g = constructor(g, op, dom, varargin);

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Fast Poisson & Helmholtz solvers:
        u = poisson( f, bc, m, n );
        u = helmholtz( f, k, bc, m, n);

        % Converts a function in polar coordinates to one in Cartesian
        % coordinates on the disk
        fdf = pol2cartf(f, r, th);

        % Disk harmonics
        Y = harmonic(l,m,type);

        % Convert coeffs to a diskfun
        f = coeffs2diskfun(CFS);

        varargout = coeffs2vals(U, varargin);
        varargout = vals2coeffs(U, varargin);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Private Static methods implemented by DISKFUN class.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = private, Static = true )

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties (Access = public)
        % DOMAIN: default is [-pi,pi] x [0,1] which corresponds to using
        % polar coordinates.
        % Doubled-up disk will have a domain of [-pi,pi] x [-1,1].
        idxPlus
        idxMinus
        nonZeroPoles = 0;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Private constant properties
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Constant )
        % TODO: Add support for this constant here.
        %alpha = 50;  % Growth factor control.
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function M = moveUDB(F, d)
%------------------------------------------------------------------------------
% Moves gridfunction F in vertical direction.
% Excess area is filled by extension of boundaryvalues.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: July 7, 2001.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[n, m] = size(F);
if d > 0
  if (d+1) > n
    error(' moveUDB    d >= n ')
  else
%   M = [flipud(F(2:(d+1),:)); F(1:(n-d),:)];
    B = F(1,:);
    M = [B(ones(1,d),:); F(1:(n-d),:)];
  end
elseif d < 0
  if (-d+1) > n
    error(' moveUDB   -d >= n ')
  else
%   M = [F((-d+1):n,:); flipud(F((n+d):(n-1),:))];
    B = F(n,:);
    M = [F((-d+1):n,:); B(ones(1,-d),:)];
  end
else
  M = F;
end
%------------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function p = chebyshev(n,var)
%CHEBYSHEV    Chebyshev polynomial (exactly representable up to degree n<=81)
%
%   p = chebyshev(n,var)
%
%Computed by recursion
%  p(0,x) = 1
%  p(1,x) = x
%  p(n,x) = 2*x*p(n-1,x) - p(n-2,x)  for n>1
%
%Parameter var for dependent variable optional (default 'x')
%

% written  07/30/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

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

  if nargin==1
    var = 'x';
  else
    if ~ischar(var)
      error('variable must be string')
    end
  end

  if n==0
    p = polynom(1,var);
  elseif n==1
    p = polynom([1 0],var);
  else
    t1 = [1 zeros(1,n)];
    t2 = [0 1 zeros(1,n-1)];
    for i=3:n+1
      t3 = [0 2*t2(1:n)] - t1;
      t1 = t2;
      t2 = t3;
    end
    p = polynom(fliplr(t3),var);
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
close all; clear;
d = 2;
k = 3;
n = 5000;
%% Generate data
[X,label] = kmeansRnd(d,k,n);
plotClass(X,label);
%% kmeans init with kmeans++ seeding (recomended)
y = kmeans(X,kseeds(X,k));
figure;
plotClass(X,y);
%% kmeans with random initialization
y = kmeans(X,k);
figure;
plotClass(X,y);
%% kmeans init with labels
y = kmeans(X,label);
figure;
plotClass(X,y);
%% kmeans init with centers
mu = rand(d,k);
y = kmeans(X,mu);
figure;
plotClass(X,y);
%% kmeans++ seeding
mu = kseeds(X,k);
[~,y] = min(dot(mu,mu,1)'/2-mu'*X,[],1); % assign sample labels
figure;
plotClass(X,y);

author: dipum

repo: dipum-toolbox

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


text:
function [E,ER] = arap_energy(V,F,U,varargin)
  % ARAP_ENERGY This function is meant as a very human readable groundtruth for
  % evaluating ARAP energies.
  % For faster energy computation see `arap_gradient.m`
  %
  % Inputs:
  %    V  #V by dim list of mesh vertex rest positions
  %    F  #F by simplex list of element indices into V
  %    U  #V by dim list of deformed vertex positions
  %   Optional:
  %     'Energy' followed by 'spokes','spokes-and-rims', or 'elements'
  %       {'elements' or 'spokes-and-rims'} for tets or triangles respectively.
  % Outputs:
  %    E  scalar energy
  %    ER  #R list of "per-rotation" energy contributions
  %
  %
  % Known issues:
  % If E_G is the energy returned by arap_gradient, then
  % 'spokes-and-rims' energy E = 3 E_G
  % 'spokes' energy E = 2 E_G
  % 'elements' energy E = E_G
  %

  switch size(F,2)
  case 4
    energy = 'elements';
  case 3
    energy = 'spokes-and-rims';
  end
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'Energy'}, ...
    {'energy'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  dim = size(V,2);
  C = cotangent(V,F);

  n = size(V,1);
  m = size(F,1);
  switch energy
  case 'spokes'
    nr = size(V,1);
    edge_set = cell(nr,1);
    switch size(F,2)
    case 4
      AC = sparse( ...
        F(:,[2 3 1 4 4 4 3 1 2 1 2 3]),F(:,[3 1 2 1 2 3 2 3 1 4 4 4]), ...
        [C C],n,n);
    case 3
      AC = sparse(F(:,[2 3 1 3 1 2]),F(:,[3 1 2 2 3 1]),[C C],n,n);
    end
    % loop over vertices
    for r = 1:nr
      [Jr,~,Cr] = find(AC(:,r));
      Ir = repmat(r,numel(Jr),1);
      edge_set{r}.E = [Ir Jr];
      edge_set{r}.C = Cr;
    end
  case {'elements','spokes-and-rims'}
    nr = size(F,1);
    switch size(F,2)
    case 3
      I = F(:,[2 3 1]);
      J = F(:,[3 1 2]);
    case 4
      I = F(:,[2 3 1 4 4 4]);
      J = F(:,[3 1 2 1 2 3]);
    end
    edge_set = cell(nr,1);
    % loop over elements
    for r = 1:nr
      edge_set{r}.E = [I(r,:)' J(r,:)'];
      edge_set{r}.C = C(r,:);
    end
    if strcmp(energy,'spokes-and-rims')
      Fedge_set = edge_set;
      nr = size(V,1);
      edge_set = cell(nr,1);
      V2F = sparse( ...
        F,repmat(1:size(F,1),size(F,2),1)',repmat(1:size(F,2),size(F,1),1),n,m);
      % loop over vertices
      for r = 1:nr
        [~,Fr] = find(V2F(r,:));
        edge_set{r}.E = [];
        edge_set{r}.C = [];
        for f = Fr
          edge_set{r}.E = [edge_set{r}.E;Fedge_set{f}.E];
          edge_set{r}.C = [edge_set{r}.C Fedge_set{f}.C];
        end
      end

    end

  end

  ER = zeros(nr,1);
  % loop over rotations
  for r = 1:nr
    % get edges for this rotation
    Er = edge_set{r}.E;
    Cr = edge_set{r}.C;
    S = zeros(dim,dim);
    % loop over edges to build covariance matrix for rotation fitting
    for e = 1:size(Er,1)
      i = Er(e,1);
      j = Er(e,2);
      ev = V(j,:)-V(i,:);
      eu = U(j,:)-U(i,:);
      S = S + Cr(e) * eu'*ev;
    end
    R = fit_rotation(S);
    % Loop over edges to compute energy contribution
    for e = 1:size(Er,1)
      i = Er(e,1);
      j = Er(e,2);
      ev = V(j,:)-V(i,:);
      eu = U(j,:)-U(i,:);
      %ER(r) = ER(r) + 0.5 * Cr(e) * (eu - ev*R)*(eu - ev*R)';
      %ER(r) = ER(r) + 0.5 * Cr(e) * (eu*eu' - 2*ev*R*eu' + ev*R*R'*ev');
      ER(r) = ER(r) + 0.5 * Cr(e) * (eu*eu' - 2*ev*R*eu' + ev*ev');
    end
  end
  E = sum(ER);
end

author: alecjacobson

repo: gptoolbox


text:
function [ xd, yd, xdp, ydp ] = hermite_interpolant ( n, x, y, yp )

%*****************************************************************************80
%
%% HERMITE_INTERPOLANT sets up a divided difference table from Hermite data.
%
%  Discussion:
%
%    The polynomial represented by the divided difference table can be
%    evaluated by calling HERMITE_INTERPOLANT_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer N, the number of items of data
%    ( X(I), Y(I), YP(I) ).
%
%    Input, real X(N), the abscissas.
%    These values must be distinct.
%
%    Input, real Y(N), YP(N), the function and
%    derivative values.
%
%    Output, real XD(2*N), YD(2*N), the divided difference table
%    for the interpolant value.
%
%    Output, real XDP(2*N-1), YDP(2*N-1), the divided difference
%    table for the interpolant derivative.
%
  x = x ( : );
  y = y ( : );
  yp = yp ( : );
%
%  Copy the data:
%
  nd = 2 * n;
  xd = zeros ( nd, 1 );
  xd(1:2:nd-1) = x(1:n);
  xd(2:2:nd  ) = x(1:n);
%
%  Carry out the first step of differencing.
%
  yd = zeros ( nd, 1 );
  yd(1) = y(1);
  yd(3:2:nd-1) = ( y(2:n) - y(1:n-1) ) ./ ( x(2:n) - x(1:n-1) );
  yd(2:2:nd  ) = yp(1:n);
%
%  Carry out the remaining steps in the usual way.
%
  for i = 3 : nd
    for j = nd : -1 : i
      yd(j) = ( yd(j) - yd(j-1) ) / ( xd(j) - xd(j+1-i) );
    end
  end
%
%  Compute the difference table for the derivative.
%
  [ ndp, xdp, ydp ] = dif_deriv ( nd, xd, yd );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function SlabAnalysis(mCatalog,vCoastlines,vFaults,vVolcano)
% % Example: SlabAnalysis(a,coastline,faults,vo)
% -------------------------------------------------------------------------------------------------------------
% Creates a cell-array with subcatalogs for every grid node defined by mPolygon. These subcatalogs
%   contain only indices to the earthquake "rows" in mCatalog.
%
% Input parameters:
%   params.mCatalog             Earthquake catalog
%   params.vCoastlines          Coastlines
%   params.vFaults              Faults
%   params.vVolcano             Volcano
%
% Output parameters:
%   none
%
% Thomas van Stiphout
% February 28, 2006
% -------------------------------------------------------------------------------------------------------------
%

if ~strcmp(get(gcf,'name'),'3DResult plot')
% Launch GUI
hMenuFig = gui_slab;
handles = guidata(hMenuFig);

% Analyze Output
if ~ishandle(hMenuFig)
    %     handles.answer = 0;
elseif (handles.answer==0)
    disp('Canceled SlabAnalysis - Goodbye');
elseif (handles.answer==1)
    % define equidistance btw gridpoints
    params.nDx=str2double(get(handles.nDx,'String'));
    params.nDy=params.nDx;
    % Define the percentile of the depth values
    params.vPercentiles=str2double(get(handles.vPerc,'String'));
    % Radius for Map-Gridding for Columns
    params.fColRadius=str2double(get(handles.fRadius,'String'));
    % Minimum number of events that have to be in the column
    params.nColMinEvents=str2double(get(handles.nMinColEvents,'String'));
    % Volume sampling is cylindrical
    params.bCylSmp=get(handles.chkCylSmp,'Value');
    params.bCylSmpModeN=get(handles.radCylSmpModeN,'Value');
    params.bCylSmpModeR=get(handles.radCylSmpModeR,'Value');
    params.fCylSmpValue=str2double(get(handles.fCylSmpValue,'String'));
    params.fCylSmpBnd=str2double(get(handles.fCylSmpBnd,'String'));
    if ((params.bCylSmpModeN==1) && (params.bCylSmpModeR==0))
        params.nGriddingMode=1; % constant numbers
        params.nNumberEvents=params.fCylSmpValue;
    elseif ((params.bCylSmpModeN==0) && (params.bCylSmpModeR==1))
        params.nGriddingMode=0; % constant radius
    else
        disp('Something is going wrong');
    end
    % Volume sampling is 3D-spherical
    params.b3DSmp=get(handles.chk3DSmp,'Value');
    params.b3DSmpModeN=get(handles.radCylSmpModeN,'Value');
    params.b3DSmpModeR=get(handles.radCylSmpModeR,'Value');
    params.f3DSmpValue=str2double(get(handles.fCylSmpValue,'String'));
    params.f3DSmpBnd=str2double(get(handles.fCylSmpBnd,'String'));
    if ((params.bCylSmpModeN==1) && (params.bCylSmpModeR==0))
        params.n3DGriddingMode=3; % constant numbers
        params.nNumberEvents=params.fCylSmpValue;
    elseif ((params.bCylSmpModeN==0) && (params.bCylSmpModeR==1))
        params.n3DGriddingMode=4; % constant radius
    else
        disp('Something is going wrong');
    end
    %     close(get(handles.figure1,'Name'));

    if (params.bCylSmp==1)
        if ~exist('Name')
            Name=cellstr('Cylindrical Volume Sampling');
        else
            Name(size(Name,1)+1,1)=cellstr('Cylindrical Volume Sampling');
        end
    end
    if params.b3DSmp
        if ~exist('Name')
            Name=cellstr('3D-Spherical Volume Sampling');
        else
            Name(size(Name,1)+1,1)=cellstr('3D Spherical Volume Sampling');
        end
    end

    params.Name=Name;
    params.mCatalog=mCatalog;
    params.vCoastline=vCoastlines;
    params.vFaults=vFaults;
    params.vVolcanoes=vVolcano;
    params.chkVolcanoes=1;

    clear mCatalog vCoastlines vFaults vVolcano
    % Create grid
    params.bMap=1;
    params.nGriddingMode=1;
    params.fSizeRectHorizontal=nan;
    params.fSizeRectDepth=nan;
    params.fSizeRectHorizontal=nan;
    params.fSizeRectDepth=nan;

    % Create map-gridding
    % by defining equidistance btw gridpoints

    % Short distance conversion
    params.fDLon=km2deg(params.nDx);
    params.fDLat=km2deg(params.nDy);
    %
    hFigure=gcf;

    % Select grid from seismicity map
    [params.mPolygon, params.vX, params.vY, params.vUsedNodes] = ex_selectgrid(hFigure, params.fDLon, params.fDLat, 0);

    % Sampling Radius is set equally to the grid raster; indices for each grid point
    [params.caNodeIndices params.vResolution] = ex_CreateIndexCatalog2(params.mCatalog, params.mPolygon, 1, 1, ...
        [], params.fColRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);

    % preparing the matrix containing the depth-levels according to the
    % percentiles defined
    params.mDepths=ones(length(params.mPolygon),length(params.vPercentiles))*nan;
    params.mValueGrid=ones(length(params.mPolygon),1)*nan;

    % calculation the depths-levels according to the percentiles for each gridpoint
    for i=1:length(params.mPolygon)
        params.mDepths(i,:)=prctile(params.mCatalog(params.caNodeIndices{i},7),params.vPercentiles) ;
        %     params.mValueGrid(i)=length(params.mCatalog(params.caNodeIndices{i}));
    end

    % Take only depth values that are well constrained,
    % i.e. with at least 50 earthquakes per map-node; set others to NaN
    params.mDepths((params.vResolution <= params.nColMinEvents),:)=nan;
    % Create polygon with with slab-surface values x,y,z
end
delete(hMenuFig);

elseif strcmp(get(gcf,'name'),'3DResult plot')
    disp('here we go for a cross section');
    view(0,90);
    vSecLim1 = ginput(1)
    vSecLim2 = ginput(1)
    dist_=deg2km(distance(vSecLim1(2),vSecLim1(1),vSecLim2(2),vSecLim2(1)));
    azimuth_=azimuth(vSecLim1(2),vSecLim1(1),vSecLim2(2),vSecLim2(1));

    view(3);
end
for nLayer=1:length(params.vPercentiles);
    params.mPolygon(:,3)=params.mDepths(:,nLayer);
    % overgive depth values to the result-matrix
    params.mValueGrid(:,1)=params.mDepths(:,nLayer);

    for i=1:length(params.vY)
        params.mX(i,:)=params.vX;
    end
    for i=1:length(params.vX)
        params.mY(:,i)=params.vY;
    end

    params.mValueGrid(isnan(params.mPolygon(:,3)),:)=nan;

    % prepare 3D figure by reshaping
    params.mX=reshape(params.mX,length(params.vUsedNodes),1);
    params.mY=reshape(params.mY,length(params.vUsedNodes),1);
    params.mZ=ones(length(params.vUsedNodes),1)*nan;
    params.mZ(params.vUsedNodes)=params.mValueGrid(:,1);
    params.mX=reshape(params.mX,length(params.vY),length(params.vX));
    params.mY=reshape(params.mY,length(params.vY),length(params.vX));
    params.mZ=reshape(params.mZ,length(params.vY),length(params.vX));

    % transform spherical coordinates to xyz-coordinates
    % by applying distances-command

    % center of data as reference/origin point
    fOriginX_=(max(max(params.mX))+min(min(params.mX)))./2;
    fOriginY_=(max(max(params.mY))+min(min(params.mY)))./2;
    mOriginX_=ones(size(params.mX))*fOriginX_;
    mOriginY_=ones(size(params.mX))*fOriginY_;

    mYkm=deg2km(distance(params.mY,mOriginX_,mOriginY_,mOriginX_)).*sign(params.mY-mOriginY_);
    mXkm=deg2km(distance(mOriginY_,params.mX,mOriginY_,mOriginX_))*sin(deg2rad(fOriginY_)).*sign(params.mX-mOriginX_);

    % mZkm=ones(size(params.mZ))*50
    mZkm=params.mZ;

    % mZkm=ones(size(params.mZ))*50
    % mZkm(isnan(params.mZ))=nan
    % surfnorm(mXkm,mYkm,-mZkm);

    [mNxkm,mNykm,mNzkm]=surfnorm(mXkm,mYkm,-mZkm);

    % make figure
    % figure;
    % surf(params.mX,params.mY,-params.mZ);
    % surf(params.mX,params.mY,-params.mZ);
    % colorbar;
    % set(gca,'CLim',[0.6 1.4]);
    % box on
    % hold on; plot3(params.mPolygon(:,1),params.mPolygon(:,2),-params.mPolygon(:,3),'s','MarkerSize',5);
    [mNx,mNy,mNz]=surfnorm(params.mX,params.mY,-params.mZ);
    [mNxkm,mNykm,mNzkm]=surfnorm(mXkm,mYkm,-mZkm);
    % determin the gradient of the slab at each grid point
    [mGx,mGy]=gradient(params.mZ);
    [mGxkm,mGykm]=gradient(params.mZ);
    mGz=zeros(size(mGx,1),size(mGx,2));
    mGzkm=zeros(size(mGxkm,1),size(mGxkm,2));
    % take only valid gridpoints
    % mGx(isnan(mNx))=nan; mGy(isnan(mNx))=nan; mGz(isnan(mNx))=nan;

    vGx=reshape(mGx,size(mGx,1)*size(mGx,2),1);vGy=reshape(mGy,size(mGy,1)*size(mGy,2),1);vGz=reshape(mGz,size(mGz,1)*size(mGz,2),1);
    vNx=reshape(mNx,size(mNx,1)*size(mNx,2),1);vNy=reshape(mNy,size(mNy,1)*size(mNy,2),1);vNz=reshape(mNz,size(mNz,1)*size(mNz,2),1);
    % same for kmgrid
    vGxkm=reshape(mGxkm,size(mGxkm,1)*size(mGxkm,2),1);vGykm=reshape(mGykm,size(mGykm,1)*size(mGykm,2),1);vGzkm=reshape(mGzkm,size(mGzkm,1)*size(mGzkm,2),1);
    vNxkm=reshape(mNxkm,size(mNxkm,1)*size(mNxkm,2),1);vNykm=reshape(mNykm,size(mNykm,1)*size(mNykm,2),1);vNzkm=reshape(mNzkm,size(mNzkm,1)*size(mNzkm,2),1);

    mask=(~isnan(vGx) .* ~isnan(vGy) .* ~isnan(vGz) .* ~isnan(vNx) .* ~isnan(vNy) .* ~isnan(vNy));
    vGx(~mask)=nan;vGy(~mask)=nan;vGz(~mask)=nan;vNx(~mask)=nan;vNy(~mask)=nan;vNz(~mask)=nan;
    % same for km-grid
    mask=(~isnan(vGxkm) .* ~isnan(vGykm) .* ~isnan(vGzkm) .* ~isnan(vNxkm) .* ~isnan(vNykm) .* ~isnan(vNykm));
    vGxkm(~mask)=nan;vGykm(~mask)=nan;vGzkm(~mask)=nan;vNxkm(~mask)=nan;vNykm(~mask)=nan;vNzkm(~mask)=nan;
    % create a vector perpendicular to the surface define by the normal and the
    % gradient, mT is tangent and represents the axes of the sampling cylinder.
    params.mT=cross([vGx vGy vGz],[vNx vNy vNz]);
    % same for km-grid
    mTkm=cross([vGxkm vGykm vGzkm],[vNxkm vNykm vNzkm]);
%     hold on;quiver3(mXkm(:),mYkm(:),-mZkm(:),mTkm(:,1),mTkm(:,2),mTkm(:,3),2,'r');
    % hold on;quiver3(params.mX(:),params.mY(:),-params.mZ(:),params.mT(:,1),params.mT(:,2),params.mT(:,3),2,'r');
    % params.mT=params.mT(params.vUsedNodes,:)

    % we only want to have horizontal cylinder axis
    params.mT(~isnan(params.mT(:,3)),3)=0;

    for ii=1:length(Name)
        if strcmp(Name(ii),'Cylindrical Volume Sampling')
            disp(Name(ii));
            params.sComment=strcat(Name(ii),' Depth=',num2str(params.vPercentiles(nLayer)));
            [params.caNodeIndices2 params.vResolution2]=ex_CreateIndexCatalogCylinder(params.mCatalog,params.mPolygon,params.mT,...
                params.vUsedNodes,params.bCylSmpModeN,params.fCylSmpValue,params.fCylSmpBnd);
            if ((params.bCylSmpModeN==1) && (params.bCylSmpModeR==0))
                params.vcsGridNames(6) = cellstr(char('Resolution [R]'));
            else
                params.vcsGridNames(6) = cellstr(char('Resolution [N]'));
            end
        elseif strcmp(Name(ii),'3D Spherical Volume Sampling')
            disp(Name(ii));
            params.sComment=strcat(Name(ii),' Depth=',num2str(params.vPercentiles(nLayer)));
            % nGriddingMode = 3 means 3D sphere around polygonpoint
            params.nGriddingMode=3;
            % select
            % Create Indices to catalog and select quakes in time period
            [params.caNodeIndices2 params.vResolution2] = ex_CreateIndexCatalog3D(params.mCatalog, params.mPolygon, '1', params.n3DGriddingMode, ...
                params.f3DSmpValue, params.f3DSmpBnd, params.fSizeRectHorizontal, params.fSizeRectDepth);
            if ((params.b3DSmpModeN==1) && (params.b3DSmpModeR==0))
                params.vcsGridNames(6) = cellstr(char('Resolution [R]'));
            else
                params.vcsGridNames(6) = cellstr(char('Resolution [N]'));
            end
            % clear indices where not enough quakes for depth estimation
            for i=1:length(params.mPolygon)
                if isnan(params.mPolygon(i,3))
                    params.caNodeIndices2{i}=[];
                end
            end
        else
            disp('Something is going wrong');
        end

        % Calculation of bValue
        for i=1:length(params.mPolygon)
            if ~isempty(params.mCatalog(params.caNodeIndices2{i},:))
                % function [fBValue, fStdDev, fAValue, fMeanMag] =  calc_bvalue(mCatalog, fBinning)
                [params.mValueGrid(i,2),params.mValueGrid(i,3),params.mValueGrid(i,4),params.mValueGrid(i,5),] =  calc_bvalue(params.mCatalog(params.caNodeIndices2{i},:));
                params.mValueGrid(i,6)=max(params.vResolution2{i});
            else
                % create NaN's in mValueGrid, where strike is not defined
                params.mValueGrid(i,:)=nan;
            end

        end

        params.vcsGridNames(1:5) = cellstr(char(...
            'Depth Level',...               %   1
            'b-Value',...                   %   2
            'Standard Deviation',...        %   3
            'A-Value',...                   %   4
            'Mean Magnitude'));             %   5

        params.bMap=2;
        params.mValueGrid(isnan(params.mValueGrid(:,1)),:)=nan;
        vResults((nLayer-1)*length(Name)+ii)=params;

    end
end
save vResults
gui_result2(vResults);
disp('surf done');

% % plot the gradient on the slab
% hold on;quiver3(params.mX(:),params.mY(:),-params.mZ(:),mGx(:),mGy(:),0,2);
% % plot strike
% hold on;quiver3(params.mX(:),params.mY(:),-params.mZ(:),params.mT(:,1),params.mT(:,2),params.mT(:,3),3,'r');

% prepare Value that will be illustrated on the 3D surface defined above
% params.mC=ones(length(params.vUsedNodes),1)*nan;
% params.mC(params.vUsedNodes)=params.mValueGrid(:,2);
% params.mC(params.mC == 0)=nan;
% params.mC=reshape(params.mC,length(params.vY),length(params.vX));
% params.mZ(isnan(params.mC))=nan;

% make figure
% figure;
% surf(params.mX,params.mY,-params.mZ);
% surf(params.mX,params.mY,-params.mZ,params.mC);
% colorbar;
% set(gca,'CLim',[0.6 1.4]);
% box on
% shading interp

author: CelsoReyes


text:
% Saving Input Files In Parts Example
%
% This example demonstrates how to save the HDF5 input files required by
% the C++ code in parts. It builds on the Running C++ Simulations Example.
%
% author: Bradley Treeby
% date: 3rd December 2013
% last update: 13th February 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% modify this parameter to run the different examples
example_number = 1;
% 1: Save the input data to disk in parts
% 2: Reload the output data from disk

% input and output filenames (these must have the .h5 extension)
input_filename  = 'example_input.h5';
output_filename = 'example_output.h5';

% pathname for the input and output files
pathname = tempdir;

% remove input file if it already exists
if example_number == 1 && exist([pathname input_filename], 'file')
    delete([pathname input_filename]);
end

% load HDF5 constants
run([getkWavePath 'private/getH5Literals']);

% =========================================================================
% SIMULATION SETTINGS
% =========================================================================

% set the properties of the computational grid
Nx = 256;                   % number of grid points in the x direction
Ny = 128;                   % number of grid points in the y direction
Nz = 64;                    % number of grid points in the z direction
dx = 0.1e-3;                % grid point spacing in the x direction [m]
dy = 0.1e-3;                % grid point spacing in the y direction [m]
dz = 0.1e-3;                % grid point spacing in the z direction [m]
Nt = 1200;                  % number of time steps
dt = 15e-9;                 % time step [s]

% set the properties of the perfectly matched layer
pml_x_size  = 10;           % [grid points]
pml_y_size  = 10;           % [grid points]
pml_z_size  = 10;           % [grid points]
pml_x_alpha = 2;            % [Nepers/grid point]
pml_y_alpha = 2;            % [Nepers/grid point]
pml_z_alpha = 2;            % [Nepers/grid point]

% define a scattering ball
ball_radius = 20;           % [grid points]
ball_x      = Nx/2 + 40;    % [grid points]
ball_y      = Ny/2;         % [grid points]
ball_z      = Nz/2;         % [grid points]

% define the properties of the medium
c0_background   = 1500;     % [kg/m^3]
c0_ball         = 1800;     % [kg/m^3]
rho0_background = 1000;     % [kg/m^3]
rho0_ball       = 1200;     % [kg/m^3]
alpha_coeff     = 0.75;     % [dB/(MHz^y cm)]
alpha_power     = 1.5;

% define a the properties of a single square source element facing in the
% x-direction
source_y_size   = 60;       % [grid points]
source_z_size   = 30;       % [grid points]
source_freq     = 2e6;      % [Hz]
source_strength = 0.5e6;    % [Pa]

% =========================================================================
% WRITE THE INPUT FILE
% =========================================================================

if example_number == 1

    % ---------------------------------------------------------------------
    % WRITE THE MEDIUM PARAMETERS
    % ---------------------------------------------------------------------

    % update command line status
    tic; fprintf('Writing medium parameters... ');

    % :::---:::---:::---:::---:::---:::---:::---:::---:::---:::---:::---:::

    % create the scattering ball and density matrix
    ball       = makeBall(Nx, Ny, Nz, ball_x, ball_y, ball_z, ball_radius, [], true);
    rho0       = rho0_background*ones(Nx, Ny, Nz, 'single');
    rho0(ball) = rho0_ball;

    % make sure the input is in the correct data format
    eval(['rho0 = ' MATRIX_DATA_TYPE_MATLAB '(rho0);']);

    % save the density matrix to both regular and staggered grid parameters
    writeMatrix([pathname input_filename], rho0, 'rho0');
    writeMatrix([pathname input_filename], rho0, 'rho0_sgx');
    writeMatrix([pathname input_filename], rho0, 'rho0_sgy');
    writeMatrix([pathname input_filename], rho0, 'rho0_sgz');

    % clear variable to free memory
    clear('rho0');

    % :::---:::---:::---:::---:::---:::---:::---:::---:::---:::---:::---:::

    % create the sound speed matrix
    c0       = c0_background*ones(Nx, Ny, Nz, 'single');
    c0(ball) = c0_ball;

    % set the reference sound speed to the maximum in the medium
    c_ref = max(c0(:));

    % get the sound speed at the location of the source
    c_source = min(c0(:));

    % make sure the input is in the correct data format
    eval(['c0 = ' MATRIX_DATA_TYPE_MATLAB '(c0);']);

    % save the sound speed matrix
    writeMatrix([pathname input_filename], c0, 'c0');

    % clear variable to free memory
    clear('c0', 'ball');

    % :::---:::---:::---:::---:::---:::---:::---:::---:::---:::---:::---:::

    % make sure the inputs are in the correct data format
    eval(['alpha_coeff = ' MATRIX_DATA_TYPE_MATLAB '(alpha_coeff);']);
    eval(['alpha_power = ' MATRIX_DATA_TYPE_MATLAB '(alpha_power);']);

    % save the absorption variables
    writeMatrix([pathname input_filename], alpha_coeff, 'alpha_coeff');
    writeMatrix([pathname input_filename], alpha_power, 'alpha_power');

    % clear variables to free memory
    clear('alpha_coeff', 'alpha_power');

    % ---------------------------------------------------------------------
    % WRITE THE SOURCE PARAMETERS
    % ---------------------------------------------------------------------

    % update command line status
    toc; tic; fprintf('Writing source parameters... ');

    % define a square source mask facing in the x-direction using the
    % normal k-Wave syntax
    p_mask = false(Nx, Ny, Nz);
    p_mask(1 + pml_x_size, Ny/2 - source_y_size/2:Ny/2 + source_y_size/2, Nz/2 - source_z_size/2:Nz/2 + source_z_size/2) = 1;

    % find linear source indices
    p_source_index = find(p_mask == 1);
    p_source_index = reshape(p_source_index, [], 1);

    % make sure the input is in the correct data format
    eval(['p_source_index = ' INTEGER_DATA_TYPE_MATLAB '(p_source_index);']);

    % save the source index matrix
    writeMatrix([pathname input_filename], p_source_index, 'p_source_index');

    % clear variables to free memory
    clear p_mask p_source_index;

    % define a time varying sinusoidal source
    p_source_input  = source_strength.*sin(2*pi*source_freq*(0:(Nt-1))*dt);

    % apply an cosine ramp to the beginning to avoid startup transients
    ramp_length = round((2*pi/source_freq)/dt);
    p_source_input(1:ramp_length) = p_source_input(1:ramp_length).*(-cos( (0:(ramp_length-1))*pi/ramp_length ) + 1)/2;

    % scale the source magnitude to be in the correct units for the code
    p_source_input = p_source_input .* (2*dt./(3*c_source*dx));

    % cast matrix to single precision
    eval(['p_source_input = ' MATRIX_DATA_TYPE_MATLAB '(p_source_input);']);

    % save the input signal
    writeMatrix([pathname input_filename], p_source_input, 'p_source_input');

    % clear variables to free memory
    clear('p_source_input');

    % ---------------------------------------------------------------------
    % WRITE THE SENSOR PARAMETERS
    % ---------------------------------------------------------------------

    % update command line status
    toc; tic; fprintf('Writing sensor parameters... ');

    % define a sensor mask through the central plane
    sensor_mask = false(Nx, Ny, Nz);
    sensor_mask(:, :, Nz/2) = 1;

    % extract the indices of the active sensor mask elements
    sensor_mask_index = find(sensor_mask);
    sensor_mask_index = reshape(sensor_mask_index, [], 1);

    % make sure the input is in the correct data format
    eval(['sensor_mask_index = ' INTEGER_DATA_TYPE_MATLAB '(sensor_mask_index);']);

    % save the sensor mask
    writeMatrix([pathname input_filename], sensor_mask_index, 'sensor_mask_index');

    % clear variables to free memory
    clear('sensor_mask', 'sensor_mask_index');

    % ---------------------------------------------------------------------
    % WRITE THE GRID PARAMETERS AND FILE ATTRIBUTES
    % ---------------------------------------------------------------------

    % update command line status
    toc; tic; fprintf('Writing grid parameters and attributes... ');

    % write grid parameters
    writeGrid([pathname input_filename], [Nx, Ny, Nz], [dx, dy, dz], ...
        [pml_x_size, pml_y_size, pml_z_size], [pml_x_alpha, pml_y_alpha, pml_z_alpha], ...
        Nt, dt, c_ref);

    % write flags
    writeFlags([pathname input_filename]);

    % set additional file attributes
    writeAttributes([pathname input_filename]);

    toc;

    % display the required syntax to run the C++ simulation
    disp(['Using a terminal window, navigate to the ' filesep 'binaries folder of the k-Wave Toolbox']);
    disp('Then, use the syntax shown below to run the simulation:');
    if isunix
        disp(['./kspaceFirstOrder3D-OMP -i ' pathname input_filename ' -o ' pathname output_filename ' --p_final --p_max']);
    else
        disp(['kspaceFirstOrder3D-OMP.exe -i ' pathname input_filename ' -o ' pathname output_filename ' --p_final --p_max']);
    end

    return

% =========================================================================
% READ THE OUTPUT FILE AND PLOT VISUALISATION
% =========================================================================

else

    % load output data from the C++ simulation
    sensor_data.p_final = h5read([pathname output_filename], '/p_final');
    sensor_data.p_max   = h5read([pathname output_filename], '/p_max');

    % take an x-y slice through the final pressure output (this is recorded
    % over the entire grid)
    sensor_data.p_final = squeeze(sensor_data.p_final(:, :, Nz/2));

    % reshape the maximum pressure output (this is recorded at the grid points
    % specified by the sensor mask)
    sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Ny]);

    % add a display mask
    ball_outline = makeCircle(Nx, Ny, ball_x, ball_y, ball_radius);
    sensor_data.p_max  (ball_outline == 1) = max(sensor_data.p_max(:));
    sensor_data.p_final(ball_outline == 1) = max(sensor_data.p_final(:));

    % remove the pml
    sensor_data.p_max   = sensor_data.p_max  (1 + pml_x_size:end - pml_x_size, 1 + pml_y_size:end - pml_y_size);
    sensor_data.p_final = sensor_data.p_final(1 + pml_x_size:end - pml_x_size, 1 + pml_y_size:end - pml_y_size);

    % get a suitable plot scale
    x_vec = (0:(Nx-2*pml_x_size-1))*dx;
    y_vec = (0:(Ny-2*pml_y_size-1))*dy;
    [x_sc, scale, prefix] = scaleSI(max([x_vec, y_vec]));

    % plot the final pressure field in the x-y plane
    figure;
    subplot(1, 2, 1);
    imagesc(y_vec*scale, x_vec*scale, sensor_data.p_final, [-1, 1]*source_strength);
    colormap(getColorMap);
    xlabel(['y [' prefix 'm]']);
    ylabel(['x [' prefix 'm]']);
    axis image;
    title('Final Pressure Field');

    % plot the maximum pressure field in the x-y plane
    subplot(1, 2, 2);
    imagesc(y_vec*scale, x_vec*scale, sensor_data.p_max, [-2, 2]*source_strength);
    colormap(getColorMap);
    xlabel(['y [' prefix 'm]']);
    ylabel(['x [' prefix 'm]']);
    axis image;
    title('Maximum Pressure');

end

repo: Aero-matlab


text:
 function sn = nufft_scale(Nd, Kd, alpha, beta, Nmid)
%function sn = nufft_scale(Nd, Kd, alpha, beta, Nmid)
%|
%| Compute scaling factors for NUFFT
%|
%| in
%|	Nd,Kd
%|	alpha	{d}
%|	beta	{d}
%|
%| option
%|	Nmid	[d]		midpoint: floor(Nd/2) or default (Nd-1)/2
%|
%| out
%|	sn	[[Nd]]		scaling factors
%|
%| Copyright 2004-7-8, Jeff Fessler, University of Michigan

if nargin == 1 && streq(Nd, 'test'), nufft_scale_test, return, end
if nargin < 4, ir_usage, end

if nargin < 5, Nmid = (Nd-1)/2; end

dd = length(Nd);
if dd == 1 && ~iscell(alpha) % 1D case
	sn = nufft_scale1(Nd(1), Kd(1), alpha, beta, Nmid(1));
return
end

% scaling factors: "outer product" of 1D vectors
sn = 1;
for id=1:dd
	tmp = nufft_scale1(Nd(id), Kd(id), alpha{id}, beta{id}, Nmid(id));
	sn = sn(:) * tmp';
end
if length(Nd) > 1
	sn = reshape(sn, Nd);	% [(Nd)]
else
	sn = sn(:);	% [(Nd)]
end

% Compute scaling factors for 1D NUFFT (from Fourier series coefficients)
% in:
%	N,K
%	alpha
%	beta
% out:
%	sn	[N]		scaling factors
%
% Copyright 2001-10-4, Jeff Fessler, The University of Michigan

function sn = nufft_scale1(N, K, alpha, beta, Nmid)

if ~isreal(alpha(1)), error 'need real alpha_0', end
L = length(alpha) - 1;

%
% compute scaling factors from Fourier coefficients
%
if L > 0
	sn = zeros(N,1);
	n = [0:(N-1)]';
	i_gam_n_n0 = 1i * (2*pi/K) * (n - Nmid) * beta;

	for l1=-L:L
		alf = alpha(abs(l1)+1);
		if l1 < 0, alf = conj(alf); end
		sn = sn + alf * exp(i_gam_n_n0 * l1);
	end

else
	sn = alpha * ones(N,1);
end

% self test
function nufft_scale_test

N = 100;
K = 2*N;
alpha = [1.0 -0.0 -0.2];
sn = nufft_scale(N, K, alpha, 1);
if im
	n = [0:(N-1)]';
	clf, plot(n, real(sn), 'r-', n, imag(sn), 'b-')
	ir_legend({'$s[n]$ real', '$s[n]$ imag'})
	ylabelf('$s[n]$')
	xlabelf('$n$')
end

pr minmax(real(sn))
pr minmax(imag(sn))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ labels ] = SegmentToLabels( segment )
%SEGMENTTOLABELS Convert output so labels start from 1 and are consecutive

values = unique(segment);
labels = zeros(size(segment));
for i = 1:length(values)
    labels(segment==values(i)) = i;
end

author: Cloud-CV

repo: object-proposals


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (gaussian_inf_engine)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes;
O = find(~isemptycell(evidence));
H = find(isemptycell(evidence));
vals = cat(1, evidence{O});

% Compute Pr(H|o)
[Hmu, HSigma, loglik] = condition_gaussian(engine.mu, engine.Sigma, H, O, vals(:), ns);

engine.Hmu = Hmu;
engine.HSigma = HSigma;
engine.hnodes = H;

%%%%%%%%

function [mu2, Sigma2, loglik] = condition_gaussian(mu, Sigma, X, Y, y, ns)
% CONDITION_GAUSSIAN Compute Pr(X|Y=y) where X and Y are jointly Gaussian.
% [mu2, Sigma2, ll] = condition_gaussian(mu, Sigma, X, Y, y, ns)

if isempty(y)
  mu2 = mu;
  Sigma2 = Sigma;
  loglik = 0;
  return;
end

use_log = 1;

if length(Y)==length(mu) % instantiating every variable
  mu2 = y;
  Sigma2 = zeros(length(y));
  loglik = gaussian_prob(y, mu, Sigma, use_log);
  return;
end

[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns);
K = SXY*inv(SYY);
mu2 = muX + K*(y-muY);
Sigma2 = SXX - K*SYX;
loglik = gaussian_prob(y, muY, SYY, use_log);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clc;
clearvars;
n = 8;
A = gallery('lotkin',n)
H0 = triu(A, -1)
[H1, Z] = spx.la.schur.francis_step(H0)
% for i=1:4
%     H = spx.la.schur.francis_step(H)
%     H(n, n-1)
% end

HD = H0 - Z * H1 * Z'
max(max(abs(HD)))

author: indigits

repo: sparse-plex


text:
function varargout = waterfall(varargin)
%WATERFALL   Waterfall plot of a CHEBFUN2.
%   WATERFALL(F) displays the waterfall plot of F.
%
%   WATERFALL(F, S) displays the column and row chebfuns of F that are used for
%   its approximation.  This is a 3D version of plot(f,S), where S is a string
%   (see PLOT).
%
%   WATERFALL(F, S, 'nslices', N) displays the first min(N,length(f)) columns
%   and rows.
%
%   WATERFALL supports passing options to the plot, similar to standard Matlab
%   plot commands. The options supported are:
%       'color':      Color of lines and markers plotted.
%       'marker':     Marker for pivot points.
%       'markersize': Size of markers plotted at pivot points.
%
%   H = WATERFALL(...) returns a handle to a waterfall plot object.
%
% See also PLOT.

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

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

end


text:
function [features,F,M] = lpc2spec(lpcas, nout)
% [features,F,M] = lpc2spec(lpcas,nout)
%    Convert LPC coeffs back into spectra
%    nout is number of freq channels, default 17 (i.e. for 8 kHz)
% 2003-04-11 dpwe@ee.columbia.edu  part of rastamat

if nargin < 2
  nout = 17;
end

[rows, cols] = size(lpcas);
order = rows - 1;

gg = lpcas(1,:);
aa = lpcas./repmat(gg,rows,1);

% Calculate the actual z-plane polyvals: nout points around unit circle
zz = exp((-j*[0:(nout-1)]'*pi/(nout-1))*[0:order]);

% Actual polyvals, in power (mag^2)
features =  ((1./abs(zz*aa)).^2)./repmat(gg,nout,1);

F = zeros(cols, floor(rows/2));
M = F;

for c = 1:cols;
  aaa = aa(:,c);
  rr = roots(aaa');
  ff = angle(rr');
%  size(ff)
%  size(aaa)
  zz = exp(j*ff'*[0:(length(aaa)-1)]);
  mags = sqrt(((1./abs(zz*aaa)).^2)/gg(c))';

  [dummy, ix] = sort(ff);
  keep = ff(ix) > 0;
  ix = ix(keep);
  F(c,1:length(ix)) = ff(ix);
  M(c,1:length(ix)) = mags(ix);
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function settings = meshMatchSettings(src, tgt, varargin);
% Adjust the view settings on one or more target meshes to match that on
% the source mesh.
%
%   settings = meshMatchSettings([sourceMesh=current mesh], [targetMeshes], [gray/volume view]);
%
% INPUTS:
%	sourceMesh: mesh whose view settings you want to replicate on other
%	meshes. You can provide a mesh structure, the number of a mesh in the
%	current gray view, or the number of the mesh window as a character (e.g.,
%	use '1' for a mesh window named 'mrMesh 1', but use 1 for the
%	first mesh in the view). [Default: selected mesh in gray view].
%
%	targetMeshes: specification of one or more meshes whose view settings
%	will match that of the source mesh. Can be specified in the same way as
%	sourceMesh. If you want to specify multiple target meshes, you can
%	specify them as a cell array (or, in the case of numeric mesh indexes,
%	a numeric array, or a struct array in the case of full mesh strucures).
%	[Default: prompt user to select meshes].
%
%	[gray/volume view]: a gray or volume view structure which has the
%	relevant meshes loaded. This isn't needed if you're providing the full
%	mesh structures as inputs. [Default: use selected gray view]
%
% RETURNS:
%	settings: settings structure from the source mesh.
%
% SEE ALSO: meshSettings, meshStoreSettings, meshRetrieveSettings.
%
% ras, 08/2009.

% I allow several target meshes to be input as varargin (allowing me to type the
% IDs as strings): the view G is always the last of these, if it's specified.
if ~isempty(varargin)
	if isstruct(varargin{end})
		G = varargin{end};
	else
		G = [];
	end
	tgt = [{tgt} varargin];
end

if notDefined('G')
	% don't call getSelectedGray if both src and target are provided as
	% structures -- we don't want to force this code to depend on the
	% mrVista data structure if we don't need to.
	if ( isstruct(src) | (iscell(src) & isstruct(src{1})) ) & ...
		( isstruct(tgt) | (iscell(tgt) & isstruct(tgt{1})) )
		G = [];
	else
		G = getSelectedGray;
	end
end

if notDefined('src'),
	src = G.mesh{G.meshNum3d};
end

% parse the source mesh structure
if ~isstruct(src)
	src = parseMeshStructure(src, G);
end

% get the source mesh settings
settings = meshSettings(src);

if notDefined('tgt'),
	% user dialog
	dlg.fieldName = 'whichMeshes';
	dlg.style = 'listbox';
	dlg.string = sprintf(['Apply Mesh Settings from Mesh %i to which ' ...
						  'other meshes?'], src.id);
	for n = 1:length(G.mesh)
		dlg.list{n} = [num2str(G.mesh{n}.id) '. ' G.mesh{n}.name];
	end
	dlg.value = 1;

	[resp ok] = generalDialog(dlg, mfilename);
	if ~ok, fprintf('[%s]: User aborted.\n', mfilename); return; end

	for n = 1:length(resp.whichMeshes)
		tgt(n) = cellfind(dlg.list, resp.whichMeshes{n});
	end
end

% parse the target mesh specification
if ~ischar(tgt) & length(tgt) > 1
	% iteratively apply to several meshes
	for n = 1:length(tgt)
		if iscell(tgt)
			currTarget = parseMeshStructure(tgt{n}, G);
		else
			currTarget = parseMeshStructure(tgt(n), G);
		end
		meshApplySettings(currTarget, settings);
	end

else
	tgt = parseMeshStructure(tgt, G);
	meshApplySettings(tgt, settings);

end

return
% /--------------------------------------------------------------------/ %

% /--------------------------------------------------------------------/ %
function msh = parseMeshStructure(msh, G);
% given a mesh specification (which could be one of many things) and a view
% with meshes, return a mesh structure.
if isstruct(msh),
	% we're good...
	return;
end

if iscell(msh)
	for n = 1:length(msh)
		msh(n) = parseMeshStructure(msh{n}, G);
	end
	return
end

if ischar(msh)
	% id of mesh window ... find it in this view:
	targetID = str2num(msh);
	for n = 1:length(G.mesh)
		ids(n) = G.mesh{n}.id;
	end
	msh = find(ids==targetID);
	if isempty(msh)
		error('Couldn''t find mesh id %i in view.', targetID);
	end

	% now it's numeric, so it should go through the indexing code below...
end

if isnumeric(msh)
	% numeric index into G's meshes
	whichMeshes = msh;  clear msh
	for n = 1:length(whichMeshes)
		msh(n) = G.mesh{whichMeshes(n)};
	end
end

if isempty(msh)
	error('Empty mesh specification.')
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function sS=transpose(sS)
% transpose list of slipSystem

sS.b = sS.b.';
sS.l = sS.l.';
sS.u = sS.u.';

author: mtex-toolbox

repo: mtex

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


text:
function ker0_values_test ( )

%*****************************************************************************80
%
%% KER0_VALUES_TEST demonstrates the use of KER0_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'KER0_VALUES_TEST:\n' );
  fprintf ( 1, '  KER0_VALUES stores values of \n' );
  fprintf ( 1, '  the Kelvin function KER of order 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = ker0_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the string to be incremented.
%
%    Output, string FILE_NAME, the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end


text:
function h = fh03 ( p, varargin )

%*****************************************************************************80
%
%% FH03 returns a uniform mesh size function.
%
%  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 P(NP,ND), the point coordinates.
%
%    Input, VARARGIN, room for extra arguments.
%
%    Output, real H(NP,1), the mesh size function.
%
  np = size ( p, 1 );
  h = ones ( np, 1 );

  return
end


text:
clear all, close all, clc
%% Create a simple signal with two frequencies
dt = .001;
t = 0:dt:1;
f = sin(2*pi*50*t) + sin(2*pi*120*t); % Sum of 2 frequencies
f = f + 2.5*randn(size(t));  %  Add some noise

%% Compute the Fast Fourier Transform FFT
n = length(t);
fhat = fft(f,n);       % Compute the fast Fourier transform
PSD = fhat.*conj(fhat)/n; % Power spectrum (power per freq)
freq = 1/(dt*n)*(0:n); % Create x-axis of frequencies in Hz
L = 1:floor(n/2);   % Only plot the first half of freqs

%% Use the PSD to filter out noise
indices = PSD>100;  % Find all freqs with large power
PSDclean = PSD.*indices;  % Zero out all others
fhat = indices.*fhat;  % Zero out small Fourier coeffs. in Y
ffilt = ifft(fhat); % Inverse FFT for filtered time signal

%% PLOTS
subplot(3,1,1)
plot(t,f,'r','LineWidth',1.2), hold on
plot(t,f,'k','LineWidth',1.5)
legend('Noisy','Clean')

subplot(3,1,2)
plot(t,f,'k','LineWidth',1.5), hold on
plot(t,ffilt,'b','LineWidth',1.2)
legend('Clean','Filtered')

subplot(3,1,3)
plot(freq(L),PSD(L),'r','LineWidth',1.5), hold on
plot(freq(L),PSDclean(L),'-b','LineWidth',1.2)
legend('Noisy','Filtered')

%
% %% PLOTS
% subplot(3,1,1)
% plot(t,y,'r','LineWidth',1.2)
% hold on
% plot(t,x,'k','LineWidth',1.5)
% axis([0 .25 -5 5])
% legend('Noisy','Clean')
%
% subplot(3,1,2)
% plot(t,x,'k','LineWidth',1.5)
% hold on
% plot(t,yfilt,'b','LineWidth',1.2)
% axis([0 .25 -5 5])
% legend('Clean','Filtered')
%
% subplot(3,1,3)
% plot(freq(L),PSD(L),'r','LineWidth',1.5)
% hold on
% plot(freq(L),PSDclean(L),'-b','LineWidth',1.2)
% legend('Noisy','Filtered')

author: dynamicslab


text:
function T=dewPointTemp4Pres(p,algChoice)
%%DEWPOINTTEMP4PRES For a given partial vapor pressure of water, find the
%                   dew point temperature. that is the temperature at which
%                   the air would be saturated with water, in equilibrium.
%
%INPUTS: p The (partial) pressure of the water vapor for which the dew
%          point temperature is desired in units of Newtons per square
%          meter (Pascals).
% algChoice An optional parameter specifying the algorithm to use. The
%          choices are
%          0 The corrected version of the Clausius-Clapeyron equation for
%            use over land or in the upper air.
%          1 The empirical Magnus-type equation for use over water or in
%            the upper air (-40C to 50C temperature).
%          2 The empirical Magnus-type equation for use over ice (-80C to
%            0C temperature)
%          If algChoice is omitted, then the default value of 0, the
%          corrected Clausius-Clapeyron equation is used. Choice 0 is taken
%          from [1]. Choices 1 and 2, are taken from [2].
%
%As is the case in the function dewPointPres4Temp, which is the inverse of
%this function, formulae 0 is from [1], where the numerical inversion
%method given in equations 44 and 45 in Section 5 of the paper are used.
%
%Formulas 1 and 2 are obtained by solving for the inverse of the Magnus
%approximations given in [2]. Simple, explicit solutions for the function
%inverses can be found.
%
%Note that the function dewPointPres4Temp is the inverse of this function.
%
%REFERENCES:
%[1] D. Koutsoyiannis, "Clausius-Clapeyron equation and saturation vapour
%    pressure: simple theory reconcided with practice," European Journal of
%    Physics, vol. 33, no. 2, pp. 295-305, Mar. 2012.
%[2] O. A. Alduchov and R. E. Eskridge, "Improved Magnus Form Approximation
%    of Saturation Vapor Pressure," Journal of Applied Meteorology, vol.
%    35, no. 4, pp. 601-609, Apr. 1996.
%
%February 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(algChoice))
    algChoice=0;
end

switch(algChoice)
    case 1
        %This evaluates the inverse of Equation 21 in the Alduchov and
        %Eskridge paper.

        p=p/100;%Convert to hectoPascals.
        logRat=log(p/6.1094);

        T=-243.04*logRat./(logRat-17.625);
        %Convert the temperature to degrees Kelvin.
        T=T-Constants.absoluteZero;
    case 2
        %This evaluates the inverse of Equation 23 in the Alduchov and
        %Eskridge paper for the saturation temperature over ice.

        p=p/100;%Convert to hectoPascals.
        logRat=log(p/6.1121);

        T=-273.86*logRat./(logRat-22.587);
        %Convert the temperature to degrees Kelvin.
        T=T-Constants.absoluteZero;
    otherwise
        %The inverse of Equation 23 in the Koutsoyiannis paper is not
        %simple to find. However, the paper provides iterative solutions in
        %the form of Equations 44 and 45 for specific values related to
        %water. For general values, the initial estimate in Equation 44 is
        %T0/T=1+1/(alpha/(R*T0)-(cL-cP)/R)*ln(p0/p);
        %and the iteration to refine the ratio is
        %T0/T_{new}=1+(R*T0/alpha)*ln(p0/p)+((cL-cP)/R)/(alpha/(R*T0))*ln(T0/T_{old})
        %where R is the specific gas constant of water vapor with units of
        %J/(kg*K), T0 is the temperature at the triple point of water, with
        %units of Kelvin, p0 is the pressure at the triple point of water
        %in hecoPascals, cP and cL are the specific heat of water vapor and
        %liquid water at constant pressure with units of J/(kg*K), and
        %alpha=L0+(cL-cP)*T0, where L0 is the latent heat of water for a
        %constant pressure at the triple point with units of J/kg.
        %
        %In the implementation here, the specific numerical from the paper
        %are used. This means that alpha has been tweaked a little bit to
        %better match their real data.

        numIter=27;

        p0=6.11657*100;%The pressure at the triple point of water in Pascals.
        T0=273.16;%The temperature at the triple point of water in Kelvin.

        LPRat=log(p0./p);

        TRat=1+1/(24.921-5.06)*LPRat;
        for curIter=1:numIter
            TRat=1+(1/(24.921))*LPRat+(5.06/24.921).*log(TRat);
        end
        T=T0./TRat;
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 [ mD ] = CalcDistanceMatrixBRows( mX )
% ----------------------------------------------------------------------------------------------- %
% [ mD ] = CalcDistanceMatrixBRows( mX )
%   Calculates the distance matrix for the input data. The distance matrix
%   is a symmetric matrix where 'mD(ii, jj) = dist(mX(ii, :), mX(jj, :));'.
%   This function uses the squared Euclidean Distance for the distance
%   metric.
% Input:
%   - mX            -   Data Matrix.
%                       Each data sample is a row of the matrix.
%                       Structure: Matrix (numVars x varDim).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - mD            -   Distance Matrix.
%                       A symmetric matrix where 'mD(ii, jj) = dist(mX(ii,
%                       :), mX(jj, :));'.
%                       Structure: Matrix (numVars x numVars).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, inf).
% References
%   1.  A
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     01/01/2021  Royi Avital	RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

numVar = size(mX, 1);
varDim = size(mX, 2);

mD = sum((reshape(mX, numVar, 1, varDim) - reshape(mX, 1, numVar, varDim)) .^ 2, 3);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [nodes, edges, faces] = boundedVoronoi2d(box, germs)
%BOUNDEDVORONOI2D Computes a bounded voronoi diagram as a graph structure.
%
%   [NODES, EDGES, FACES] = boundedVoronoi2d(BOX, GERMS)
%   GERMS an array of points with dimension 2
%   NODES, EDGES, FACES: usual graph representation, FACES as cell array
%
%   Example
%     % clip a graph with
%     box = [0 100 0 100];
%     [n, e, f] = boundedVoronoi2d(box, rand(100, 2)*100);
%     [n, e, f] = clipGraph(n, e, f, box);
%     drawGraph(n, e);
%
%   See also
%     graphs, boundedCentroidalVoronoi2d, clipGraph, clipGraphPolygon

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2007-01-12
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

% uniformize input for box.
box = box';
box = box(:);

% add points far enough
boxSizeX  = box(2) - box(1);
boxSizeY  = box(4) - box(3);
farPoints = [...
    box(2)+2*boxSizeX  box(4)+3*boxSizeY;...
    box(1)-3*boxSizeX  box(4)+2*boxSizeY;...
    box(1)-2*boxSizeX  box(3)-3*boxSizeY;...
    box(2)+3*boxSizeX  box(3)-2*boxSizeY;...
    ];

% extract voronoi vertices and face structure
[V, C] = voronoin([germs ; farPoints]);

% initialize graph structure, without edges and without faces
nodes = V(2:end, :);
edges = zeros(0, 2);
faces = cell(1, size(germs, 1));

% for each cell associated to a germ, we retrieve the nodes of the Voronoi
% diagram, remove 1 because first vertex at infinity is removed, keep the
% list of vertices that constitute the face, and generate the set of edges
% around face.
for i = 1:size(germs, 1)
    face = C{i};
    face = face-1;
    faces{i} = face;
    edges = [edges; sort([face' face([2:end 1])'], 2)]; %#ok<AGROW>
end

% remove duplicate edges as they were created twice (once for each face)
edges = unique(edges, 'rows');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function cvx_pushobj( dir, x )

global cvx___
try
    np = length(cvx___.problems);
    pstr = cvx___.problems(np);
catch
    cvx_throw( 'No CVX model is present.' );
end

%
% Check problem
%

if ~isinf( pstr.direction )
    if isnan( pstr.direction )
        cvx_throw( 'Objective functions cannot be added to sets.' );
    else
        cvx_throw( 'An objective has already been supplied for this problem.' );
    end
end

%
% Check direction
%

if ~ischar( dir ) || size( dir, 1 ) ~= 1,
    cvx_throw( 'The second argument must be a string.' );
end

%
% Check objective expression
%

persistent remap_min remap_max
if isempty( remap_max ),
    remap_min = cvx_remap( { { 'constant' } ; { 'l_convex' } ; { 'convex' } } );
    remap_max = cvx_remap( { { 'constant' } ; { 'l_concave' } ; { 'concave' } } );
end
if ~isa( x, 'cvx' ) && ~isnumeric( x ),
    cvx_throw( 'Cannot accept an objective of type ''%s''.', class( arg ) );
elseif ~isreal( x ),
    cvx_throw( 'Expressions in objective functions must be real.' );
elseif isempty( x ),
    warning( 'CVX:EmptyObjective', 'Empty objective.' );
end
cx = cvx_classify( x );
switch dir,
    case { 'minimize', 'minimise' },
        vx = remap_min( :, cx );
        dir = +1;
    case { 'maximize', 'maximise' },
        vx = remap_max( :, cx );
        dir = -1;
    otherwise,
        cvx_throw( 'Invalid objective type: %s', dir );
end
avx = all( vx, 2 );
if ~any( avx ),
    vx = any( vx, 1 );
    if all( vx ),
        cvx_throw( 'Cannot combine logarithmic and linear objectives.' );
    else
        cvx_dcp_error( type, 'unary', cvx_fastref( x, ~vx ) );
    end
end

%
% Store the objective
%

isgeo = avx(2) & ( ~avx(1) | pstr.gp );
if isgeo,
    x = log( x );
    dir = 2 * dir;
end
if np > 1,
    bx = cvx_sparsify( cvx_basis( x ), [], 'magnitude' );
    x = cvx( size( x ), bx );
end
if isgeo,
    x = exp( x );
end
pstr.objective = x;
pstr.direction = dir;
cp = pstr.checkpoint(4);
if cp > 1,
    x = find( any( cvx_basis( x ), 2 ), 2, 'first' );
    if any( x > 1 & x < cp ),
        pstr.checkpoint(4) = x(1+(x(1)==1)) - 1;
    end
end
cvx___.problems( end ) = pstr;

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

author: yu-jiang

repo: radpbook

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


text:
function output = callipopt(model)

options = [];
try
    options.ipopt = optiRemoveDefaults(model.options.ipopt,ipoptset());
catch
    options.ipopt = model.options.ipopt;
end
options.ipopt.print_level = 2+model.options.verbose;

% Standard NONLINEAR setup
tempF = model.F_struc;
tempK = model.K;
tempx0 = model.x0;
model = yalmip2nonlinearsolver(model);
% % Try to propagate initial values now that nonlinear evaluation logic is
% % set up. Do this if some of the intials are nan
if any(isnan(model.x0)) & ~all(isnan(model.x0))
    [~,~,xevaledout] = fmincon_fun(model.x0,model);
    startNan = nnz(isnan(xevaledout));
    goon = 1;
    while goon
        temp = propagatex0(xevaledout,tempF,tempK);
        model.x0 = temp(model.linearindicies);
        [~,~,xevaledout] = fmincon_fun(model.x0,model);
        goon =  nnz(isnan(xevaledout)) < startNan;
        startNan = nnz(isnan(xevaledout));
    end
end
model.x0(isnan(model.x0))=0;

if ~model.derivative_available
    disp('Derivate-free call to ipopt not yet implemented')
    error('Derivate-free call to ipopt not yet implemented')
end
if model.options.savedebug
    save ipoptdebug model
end
showprogress('Calling IPOPT',model.options.showprogress);

% Figure out which variables are artificially introduced to normalize
% arguments in callback operators to simplify chain rules etc. We can do
% our function evaluations and gradient computations in our lifted world,
% but only expose the model in the original variables to the nonlinear
% solver.
% model = compressLifted(model);

Fupp = [ repmat(0,length(model.bnonlinineq)+length(model.K.q)*(model.K.q(1)>0),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(0,length(model.b),1);
    repmat(0,length(model.beq),1)];

Flow = [ repmat(-inf,length(model.bnonlinineq)+length(model.K.q)*(model.K.q(1)>0),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(-inf,length(model.b),1);
    repmat(0,length(model.beq),1)];

if isempty(Flow)
    Flow = [];
    Fupp = [];
end

% Since ipopt react strangely on lb>ub, we should bail if that is detected
% (ipopt creates an exception)
if ~isempty(model.lb)
    if any(model.lb>model.ub)
        problem = 1;
        solverinput = [];
        solveroutput = [];
        output = createoutput(model.c*0,[],[],problem,'IPOPT',solverinput,solveroutput,0);
        return
    end
end

% These are needed to avoid recomputation due to ipopts double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G = [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

funcs.objective = @(x)ipopt_callback_f(x,model);
funcs.gradient = @(x)ipopt_callback_df(x,model);
if ~isempty(Fupp)
    funcs.constraints = @(x)ipopt_callback_g(x,model);
    funcs.jacobian  = @(x)ipopt_callback_dg(x,model);
end

options.lb = model.lb(:)';
options.ub = model.ub(:)';
if ~isempty(Fupp)
    options.cl = Flow;
    options.cu = Fupp;
end

if ~isempty(Fupp)
    Z = jacobiansparsityfromnonlinear(model);
    funcs.jacobianstructure = @() Z;
end

if ~model.options.usex0
    model.x0 = (options.lb+options.ub)/2;
    model.x0(isinf(options.ub)) = options.lb(isinf(options.ub))+1;
    model.x0(isinf(options.lb)) = options.ub(isinf(options.lb))-1;
    model.x0(isinf(model.x0)) = 0;
    if any(model.variabletype == 4)
        problematic = find(any(model.monomtable(:,model.linearindicies) < 0 ,1));
        if ~isempty(problematic)
            problematic = problematic(find(model.x0(problematic)==0));
            Oneisfeas = problematic(find(model.ub(problematic) > 1));
            model.x0(Oneisfeas) = 1;
        end
    end
    model.x0(find(model.lb==model.ub)) = model.lb(find(model.lb==model.ub));
end

% If quadratic objective and no nonlinear constraints, we can supply an
% Hessian of the Lagrangian
usedinObjective = find(model.c | any(model.Q,2));
if ~any(model.variabletype(usedinObjective)) & any(model.Q)
    if  length(model.bnonlinineq)==0 & length(model.bnonlineq)==0
        H = model.Q(:,model.linearindicies);
        H = H(model.linearindicies,:);
        funcs.hessian = @(x,s,l) tril(2*H);
        funcs.hessianstructure = @()tril(sparse(double(H | H)));
        options.ipopt.hessian_approximation = 'exact';
    end
end

solvertime = tic;
[xout,info] = ipopt(model.x0,funcs,options);
solvertime = toc(solvertime);

% Duals currently not supported
lambda = [];

if ~isempty(xout) && ~isempty(model.lift);
    x = zeros(length(model.linearindicies),1);
    x(model.lift.linearIndex) = xout(:);
    x(model.lift.liftedIndex) = model.lift.T*xout(:) + model.lift.d;
    x = RecoverNonlinearSolverSolution(model,x);
else
    x = RecoverNonlinearSolverSolution(model,xout);
end

switch info.status
    case {0,1}
        problem = 0;
    case {2}
        problem = 1;
    case {-1}
        problem = 3;
    case {3,4,-2,-3}
        problem = 4;
    case {-11,-12,-13}
        problem = 7;
    case {-10,-100,-101,-102,-199}
        problem = 11;
    otherwise
        problem = -1;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(x,D_struc,[],problem,'IPOPT',solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

  predictions =  X * theta;
    updates = X' * (predictions - y);
    theta = theta - alpha * (1/m) * updates;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

author: vugsus

repo: coursera-machine-learning


text:
function [ nb, ni, np, ns, pr_bud, pr_bud_angular, r, seed ] = ...
  parameter_set ( nb, ni, np, ns, pr_bud, pr_bud_angular, r, seed )

%*****************************************************************************80
%
%% PARAMETER_SET allows the user to reset parameter values;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer NB, the number of boundary points.
%
%    Input/output, integer NI, the number of interior points.
%
%    Input/output, integer NP, the total number of points.
%
%    Input/output, integer NS, the number of sample points.
%
%    Input/output, real PR_BUD, the probability that a boundary point will bud.
%
%    Input/output, real PR_BUD_ANGULAR, the probability that a boundary point bud
%    will be constrained to the boundary.
%
%    Input/output, real R, the radius of the circle.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use the "who" command to see current variables\n' );
  fprintf ( 1, '  Reset any value using the usual MATLAB assignment statements.\n' );
  fprintf ( 1, '  To terminate assignments and resume the program, type "return".\n' );
  fprintf ( 1, '  To terminate the program, type "dbquit".\n' );

  keyboard

  return
end


text:
function handles = two_link_workcell_init(vrep,id)

handles = struct('id',id);

armJoints = -1*ones(2,1);
links = -1*ones(2,1);
forcesensor = -1;
for i=1:2
[res,armJoints(i)] = vrep.simxGetObjectHandle(id,['joint',num2str(i)],vrep.simx_opmode_oneshot_wait);vrchk(vrep, res);
[res,links(i)] = vrep.simxGetObjectHandle(id,['link',num2str(i)],vrep.simx_opmode_oneshot_wait);vrchk(vrep, res);
end
[res,forcesensor] = vrep.simxGetObjectHandle(id,'Force_sensor',vrep.simx_opmode_oneshot_wait);vrchk(vrep, res);
handles.armJoints = armJoints;
handles.links = links;
handles.forcesensor = forcesensor;

%% streaming
for i=1:2
vrep.simxGetJointPosition(id,armJoints(i),vrep.simx_opmode_streaming);vrchk(vrep, res, true);
vrep.simxGetObjectFloatParameter(id,armJoints(i),2012,vrep.simx_opmode_streaming);vrchk(vrep, res, true);
vrep.simxGetJointForce(id,armJoints(i),vrep.simx_opmode_streaming);vrchk(vrep, res, true);
end
vrep.simxReadForceSensor(id,forcesensor,vrep.simx_opmode_streaming);vrchk(vrep, res, true);
vrep.simxGetObjectOrientation(id,forcesensor,-1,vrep.simx_opmode_streaming);vrchk(vrep, res, true);
% Make sure that all streaming data has reached the client at least once
vrep.simxGetPingTime(id);
end

author: xuhuairuogu

repo: V-REP-Simulation-Projects


text:
function [centres, options, post, errlog] = kmeans(centres, data, options)
%KMEANS	Trains a k means cluster model.
%
%	Description
%	 CENTRES = KMEANS(CENTRES, DATA, OPTIONS) uses the batch K-means
%	algorithm to set the centres of a cluster model. The matrix DATA
%	represents the data which is being clustered, with each row
%	corresponding to a vector. The sum of squares error function is used.
%	The point at which a local minimum is achieved is returned as
%	CENTRES.  The error value at that point is returned in OPTIONS(8).
%
%	[CENTRES, OPTIONS, POST, ERRLOG] = KMEANS(CENTRES, DATA, OPTIONS)
%	also returns the cluster number (in a one-of-N encoding) for each
%	data point in POST and a log of the error values after each cycle in
%	ERRLOG.    The optional parameters have the following
%	interpretations.
%
%	OPTIONS(1) is set to 1 to display error values; also logs error
%	values in the return argument ERRLOG. If OPTIONS(1) is set to 0, then
%	only warning messages are displayed.  If OPTIONS(1) is -1, then
%	nothing is displayed.
%
%	OPTIONS(2) is a measure of the absolute precision required for the
%	value of CENTRES at the solution.  If the absolute difference between
%	the values of CENTRES between two successive steps is less than
%	OPTIONS(2), then this condition is satisfied.
%
%	OPTIONS(3) is a measure of the precision required of the error
%	function at the solution.  If the absolute difference between the
%	error functions between two successive steps is less than OPTIONS(3),
%	then this condition is satisfied. Both this and the previous
%	condition must be satisfied for termination.
%
%	OPTIONS(14) is the maximum number of iterations; default 100.
%
%	See also
%	GMMINIT, GMMEM
%

%	Copyright (c) Ian T Nabney (1996-2001)

[ndata, data_dim] = size(data);
[ncentres, dim] = size(centres);

if dim ~= data_dim
  error('Data dimension does not match dimension of centres')
end

if (ncentres > ndata)
  error('More centres than data')
end

% Sort out the options
if (options(14))
  niters = options(14);
else
  niters = 100;
end

store = 0;
if (nargout > 3)
  store = 1;
  errlog = zeros(1, niters);
end

% Check if centres and posteriors need to be initialised from data
if (options(5) == 1)
  % Do the initialisation
  perm = randperm(ndata);
  perm = perm(1:ncentres);

  % Assign first ncentres (permuted) data points as centres
  centres = data(perm, :);
end
% Matrix to make unit vectors easy to construct
id = eye(ncentres);

% Main loop of algorithm
for n = 1:niters

  % Save old centres to check for termination
  old_centres = centres;

  % Calculate posteriors based on existing centres
  d2 = dist2(data, centres);
  % Assign each point to nearest centre
  [minvals, index] = min(d2', [], 1);
  post = id(index,:);

  num_points = sum(post, 1);
  % Adjust the centres based on new posteriors
  for j = 1:ncentres
    if (num_points(j) > 0)
      centres(j,:) = sum(data(find(post(:,j)),:), 1)/num_points(j);
    end
  end

  % Error value is total squared distance from cluster centres
  e = sum(minvals);
  if store
    errlog(n) = e;
  end
  if options(1) > 0
    fprintf(1, 'Cycle %4d  Error %11.6f\n', n, e);
  end

  if n > 1
    % Test for termination
    if max(max(abs(centres - old_centres))) < options(2) & ...
        abs(old_e - e) < options(3)
      options(8) = e;
      return;
    end
  end
  old_e = e;
end

% If we get here, then we haven't terminated in the given number of
% iterations.
options(8) = e;
if (options(1) >= 0)
  disp(maxitmess);
end


text:
%% FUNCTION Logistic_Lasso
%   Sparse Structure-Regularized Learning with Logistic Loss.
%
%% OBJECTIVE
%   argmin_{W,C} { sum_i^t (- sum(log (1./ (1+ exp(-X{i}*W(:, i) - Y{i} .* C(i)))))/length(Y{i}))
%            + rho1 * \|W\|_1 + opts.rho_L2 * \|W\|_F^2}
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   rho1: sprasity controlling parameter
%   opts.rho_L2: L2-norm regularization parameter
%
%% OUTPUT
%   W: model: d * t
%   C: model: 1 * t
%   funcVal: function value vector.
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% Related papers
%
%   [1] Tibshirani, J. Regression shrinkage and selection via
%   the Lasso, Journal of the Royal Statistical Society. Series B 1996
%
%% Related functions
%   Least_Lasso, init_opts

%% Code starts here
function [W, C, funcVal] = Logistic_Lasso(X, Y, rho1, opts)

if nargin <3
    error('\n Inputs: X, Y, and rho1 should be specified!\n');
end
X = multi_transpose(X);

if nargin <4
    opts = [];
end

% initialize options.
opts=init_opts(opts);

if isfield(opts, 'rho_L2')
    rho_L2 = opts.rho_L2;
else
    rho_L2 = 0;
end

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];

%initialize a starting point
C0_prep = zeros(1, task_num);
for t_idx = 1: task_num
    m1 = nnz(Y{t_idx} == 1);
    m2 = nnz(Y{t_idx} == -1);
    if ( m1==0 || m2==0 )
        C0_prep(t_idx) = 0;
    else
        C0_prep(t_idx) = log(m1/m2);
    end
end

if opts.init==2
    W0 = zeros(dimension, task_num);
    %C0 = zeros(1, task_num);
    C0 = C0_prep;
elseif opts.init== 0
    W0 = randn(dimension, task_num);
    C0 = C0_prep;
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .W0');
        end
    else
        W0 = zeros(dimension, task_num);
    end
    if isfield(opts,'C0')
        C0=opts.C0;
    else
        C0=C0_prep;
    end
end

bFlag=0; % this flag tests whether the gradient step only changes a little

Wz= W0;
Cz= C0;
Wz_old = W0;
Cz_old = C0;

t = 1;
t_old = 0;
iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;
    Cs = (1 + alpha) * Cz - alpha * Cz_old;

    %Ws = Wz;
    %Cs = Cz;

    % compute function value and gradients of the search point
    [gWs, gCs, Fs ]  = gradVal_eval(Ws, Cs);

    % the Armijo Goldstein line search scheme
    while true
        [Wzp l1c_wzp] = l1_projection(Ws - gWs/gamma, 2 * rho1 / gamma);
        Czp = Cs - gCs/gamma;
        Fzp = funVal_eval  (Wzp, Czp);

        delta_Wzp = Wzp - Ws;
        delta_Czp = Czp - Cs;
        nrm_delta_Wzp = norm(delta_Wzp, 'fro')^2;
        nrm_delta_Czp = norm(delta_Czp, 'fro')^2;
        r_sum = (nrm_delta_Wzp+nrm_delta_Czp)/2;

        %         Fzp_gamma = Fs + trace(delta_Wzp' * gWs)...
        %             + trace(delta_Czp' * gCs)...
        %             + gamma/2 * nrm_delta_Wzp ...
        %             + gamma/2 * nrm_delta_Czp;

        Fzp_gamma = Fs + sum(sum(delta_Wzp .* gWs))...
            + sum(sum(delta_Czp .* gCs))...
            + gamma/2 * nrm_delta_Wzp ...
            + gamma/2 * nrm_delta_Czp;

        if (r_sum <=1e-20)
            bFlag=1; % this shows that, the gradient step makes little improvement
            break;
        end

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Cz_old = Cz;
    Wz = Wzp;
    Cz = Czp;

    funcVal = cat(1, funcVal, Fzp + rho1 * l1c_wzp);

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;
C = Czp;

% private functions

    function [z, l1_comp_val] = l1_projection (v, beta)
        % this projection calculates
        % argmin_z = \|z-v\|_2^2 + beta \|z\|_1
        % z: solution
        % l1_comp_val: value of l1 component (\|z\|_1)

        z = sign(v).*max(0,abs(v)- beta/2);

        l1_comp_val = sum(sum(abs(z)));
    end

    function [grad_W, grad_C, funcVal] = gradVal_eval(W, C)
        grad_W = zeros(dimension, task_num);
        grad_C = zeros(1, task_num);
        lossValVect = zeros (1 , task_num);

        for i = 1:task_num
            [ grad_W(:, i), grad_C(:, i), lossValVect(:, i)] = unit_grad_eval( W(:, i), C(i), X{i}, Y{i});
        end

        grad_W = grad_W + rho_L2 * 2 * W;
        % here when computing function value we do not include
        % l1 norm.
        funcVal = sum(lossValVect) + rho_L2 * norm(W, 'fro')^2;
    end

    function [funcVal] = funVal_eval (W, C)
        funcVal = 0;
        for i = 1: task_num
            funcVal = funcVal + unit_funcVal_eval( W(:, i), C(i), X{i}, Y{i});
        end

        % here when computing function value we do not include
        % l1 norm.
        funcVal = funcVal + rho_L2 * norm(W, 'fro')^2;
    end

end

function [ grad_w, grad_c, funcVal ] = unit_grad_eval( w, c, x, y)
%gradient and logistic evaluation for each task
m = length(y);
weight = ones(m, 1)/m;
weighty = weight.* y;
aa = -y.*(x'*w + c);
bb = max( aa, 0);
%funcVal = weight' * log(exp(-bb) + exp(aa-bb) + bb); %bug!!!
funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
pp = 1./ (1+exp(aa));
b = -weighty.*(1-pp);
grad_c = sum(b);
grad_w = x * b;
end

function [ funcVal ] = unit_funcVal_eval( w, c, x, y)
%function value evaluation for each task
m = length(y);
weight = ones(m, 1)/m;
aa = -y.*(x'*w + c);
bb = max( aa, 0);
%funcVal = weight' * log(exp(-bb) + exp(aa-bb) + bb); % BUG!!!!
funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function [H,g,Aineq,bineq,Aeq,beq] = fcn_get_QP_form_eta(Xt,Ut,Xd,Ud,p)
% min. 0.5 * x' * H *x + g' * x
% s.t. Aineq *x <= bineq
%      Aeq * x <= beq
% X = [pc dpc vR wb pf]': [30,1]
% q = [pc dpc eta wb]: [12 1]
% lb/ub - [4,n_hor]

%% parameters
mu = p.mu;
n_hor = p.predHorizon;
Umax = p.Umax;
decayRate = p.decayRate;

R = p.R;
Q = p.Q;
Qf = p.Qf;
[Qx,Qv,Qeta,Qw] = deal(Q(1:3,1:3),Q(4:6,4:6),Q(7:9,7:9),Q(10:12,10:12));
[Qxf,Qvf,Qetaf,Qwf] = deal(Qf(1:3,1:3),Qf(4:6,4:6),Qf(7:9,7:9),Qf(10:12,10:12));

nX = 12;
nU = 12;

%%%%%%% A,B,d matrices for linear dynamics %%%%%%%%%%%
[A,B,d] = fcn_get_ABD_eta(Xt,Ut,p);

%% Decompose
Rt = reshape(Xt(7:15,1),[3,3]);
qt = [Xt(1:6);[0;0;0];Xt(16:18)];

% lb <= Fz <= ub
Fzd = Ud([3 6 9 12],:);
lb = -1 * Fzd;
ub = 2 * Fzd;

%% Matrices for QP
H = zeros((nX + nU) * n_hor);
g = zeros(size(H,1),1);
Aeq = zeros(nX * n_hor,(nX+nU) * n_hor);
beq = zeros(size(Aeq,1),1);
if p.gait == -2
    Aineq_unit = [1 0 0;-1 0 0;0 1 0;0 -1 0;0 0 1;0 0 -1];
else
    Aineq_unit = [1 0 -mu;-1 0 -mu;0 1 -mu;0 -1 -mu;0 0 1; 0 0 -1];
end
nAineq_unit = size(Aineq_unit,1);
Aineq = zeros(4*nAineq_unit*n_hor,(nX+nU)*n_hor);
bineq = zeros(size(Aineq,1),1);
for i_hor = 1:n_hor
    xd = Xd(1:3,i_hor);
    vd = Xd(4:6,i_hor);
    Rd = reshape(Xd(7:15,i_hor),[3,3]);
    wd = Xd(16:18,i_hor);

    %% Objective function
    idx_u = (i_hor-1) * (nX + nU) + (1:nU);
    idx_x = (i_hor-1) * (nX + nU) + nU + (1:nX);
    if i_hor == n_hor
        H(idx_x,idx_x) = Qf * decayRate^(i_hor-1);
        g(idx_x) = [-Qxf * xd;
                    -Qvf * vd;
                     Qetaf * veeMap(logm(Rd' * Rt));
                    -Qwf * wd] * decayRate^(i_hor-1);
    else
        H(idx_x,idx_x) = Q * decayRate^(i_hor-1);
        g(idx_x) = [-Qx * xd;
                    -Qv * vd;
                     Qeta * veeMap(logm(Rd' * Rt));
                    -Qw * wd] * decayRate^(i_hor-1);
    end
    H(idx_u,idx_u) = R * decayRate^(i_hor-1);
    g(idx_u) = R' * (Ut - Ud(:,i_hor)) * decayRate^(i_hor-1);

    %% Equality constraints
    if i_hor == 1
        Aeq(1:nX,1:(nU+nX)) = [-B,eye(nX)];
        beq(1:nX) = A * qt + d;
    else
        Aeq((i_hor-1)*nX+(1:nX),(i_hor-2)*(nX+nU)+nU+(1:(2*nX+nU)))= [-A -B eye(nX)];
        beq((i_hor-1)*nX+(1:nX)) = d;
    end

    %% Inequality constraints
    Fi = zeros(4*nAineq_unit,12);
    hi = zeros(size(Fi,1),1);
    for i_leg = 1:4
        idx_F = (i_leg-1)*nAineq_unit + (1:nAineq_unit);
        idx_u = (i_leg-1)*3 + (1:3);
        Fi(idx_F,idx_u) = Aineq_unit;
        if p.gait == -2
            hi(idx_F) = [Umax-Ut(idx_u(1));Umax+Ut(idx_u(1));...
                         Umax-Ut(idx_u(2));Umax+Ut(idx_u(2));...
                         Umax-Ut(idx_u(3));Umax+Ut(idx_u(3))];
        else
            hi(idx_F) = [mu*Ut(idx_u(3))-Ut(idx_u(1));
                         mu*Ut(idx_u(3))+Ut(idx_u(1));
                         mu*Ut(idx_u(3))-Ut(idx_u(2));
                         mu*Ut(idx_u(3))+Ut(idx_u(2));
                         ub(i_leg,i_hor)-Ut(idx_u(3))+Ud(idx_u(3),i_hor);
                        -lb(i_leg,i_hor)+Ut(idx_u(3))-Ud(idx_u(3),i_hor)];
        end
    end
    idx_A = (i_hor-1) * 4*nAineq_unit + (1:4*nAineq_unit);
    idx_z = (i_hor-1) * (nX+nU) + (1:nU);
    Aineq(idx_A,idx_z) = Fi;
    bineq(idx_A) = hi;
end

end

author: YanranDing

repo: RF-MPC

save_path: github-repos/MATLAB/YanranDing-RF-MPC


text:
function p = mtexODFPath
% returns the default path to ODF-files
%
% Syntax
%   mtexODFPath %
%   getMTEXpref('ODFPath') % long form of mtexODFPath
%   setMTEXpref('ODFPath',path) % change the path to ODF files
%
% See also
% mtexDataPath mtexCifPath mtexEBSDPath mtexPoleFigurePath mtexODFPath
% mtexTensorPath

DataPath = getMTEXpref('ODFPath');

if exist(DataPath,'dir')
  p = DataPath;
else
  error('Data package not installed!');
end


text:
% This plots the convergence of the human observer study in terms of the
% Kendall coefficients of agreement versus different numbers of workers and
% error thresholds for the sanity checks.
function plotObserverStudyConvergence(resultDir, outputDir)

    % Get precomputed coefficients.
    load([resultDir, '/observerStudy/observerStudyConvergence.mat']);

    % Mean and standard deviation of the Kendall coefficient of agreement
    % over the samples of the Monte Carlo simulation.
    kendallMean = mean(kendall, 3);
    kendallStd = std(kendall, [], 3);

    % Plot Kendall coefficient of agreement vs. number of workers.
    lineStyle = {'-', '--', '-.', ':'};
    lineColor = {[0 0 1], [1 0 0], [1 0.65 0], [0 1 0]};
    figure;
    for maxSanityCheckFailsIdx = 1:length(maxSanityCheckFails)
        errorbar(numWorkers, kendallMean(maxSanityCheckFailsIdx,:), kendallStd(maxSanityCheckFailsIdx,:), ...
            lineStyle{maxSanityCheckFailsIdx}, 'Color', lineColor{maxSanityCheckFailsIdx}, 'LineWidth', 1.5);
        legendStr{maxSanityCheckFailsIdx} = sprintf('$n_f$ = %s', num2str(maxSanityCheckFails(maxSanityCheckFailsIdx)));
        hold on;
    end
    hold off;
    grid on;
    xlabel('Number of user sessions');
    ylabel('Coefficient of agreement');
    xlim([numWorkers(1) numWorkers(end)]);

    % Add legend
    legend(legendStr, 'Location', 'NorthEast', 'Orientation', 'horizontal', 'Interpreter', 'latex');

    % Save to TikZ.
    matlab2tikz([outputDir, '/', 'userStudyConvergence.tikz'], ...
        'height', '\figureheight', 'width', '\figurewidth', 'showInfo', false, ...
        'extraaxisoptions', ['xlabel near ticks,', ...
        'ylabel near ticks,', ...
        'scaled y ticks=false,', ...
        'yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=2},']);

author: thomas-koehler

repo: SupER


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% init0.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% self-defined initialization list
L = 3*ones(n,1);
l = 2*ones(n,1);
x0 = [-10*ones(n,1) zeros(n,1) 10*ones(n,1)];

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function mtr = mtrSet(mtr,param,val,varargin)
% Set data for MetroTrac (mtr) structure
%
%  val = mtrSet(mtr,param,val,varargin)
%
% Set the parameters of a mtr structure.
% Access this structure should go through this routine and through
%    mtrGet.
%
%

if ieNotDefined('mtr'), error('No MetroTrac structure defined.'); end
if ieNotDefined('param'), error('No parameter defined'); end
if ieNotDefined('val'), val = []; end

switch lower(param)
    case {'name'}
        mtr.name = val;
    case {'type'}
        mtr.type = val;
    case {'tensors_filename'}
        mtr.tensors_filename = val;
    case {'fa_filename'}
        mtr.fa_filename = val;
    case {'pdf_filename'}
        mtr.pdf_filename = val;
    case {'mask_filename'}
        mtr.mask_filename = val;
    case {'xmask_filename'}
        mtr.Xmask_filename = val;
    case {'compute_fa'}
        mtr.compute_fa = val;
    case {'require_way'}
        mtr.require_way = val;
    case {'sampler_type'}
        mtr.sampler_type = val;
    case {'desired_samples'}
        mtr.desired_samples = val;
    case {'burnin'}
        mtr.burnin = val;
    case {'max_nodes'}
        mtr.max_nodes = val;
    case {'min_nodes'}
        mtr.min_nodes = val;
    case {'step_size'}
        mtr.step_size = val;
    case {'skip_samples'}
        mtr.skip_samples = val;
    case {'roi'}
        if length(varargin) ~= 2, error('Incomplete mtrSet'); end
        arg1 = varargin{1};
        arg2 = varargin{2};
        switch lower(arg2)
            case {'coords'}
                mtr.roi{arg1}.coords = val;
            case {'valid_cortex'}
                mtr.roi{arg1}.valid_cortex = val;
            case {'seed_region'}
                mtr.roi{arg1}.seed_region = val;
        end
    case {'translate'}
        mtr.translate = val;
    case {'isr'}
        mtr.isr = val;
    case {'esr'}
        mtr.esr = val;
    case {'temp_swap'}
        mtr.temp_swap = val;
    case {'inv_temp'}
        mtr.inv_temp = val;
    case {'start_path_tries'}
        mtr.start_path_tries = val;
    case {'save_out_spacing'}
        mtr.save_out_spacing = val;
    case {'fa_absorb'}
        mtr.fa_absorb = val;
    case {'abs_normal'}
        mtr.abs_normal = val;
    case {'abs_penalty'}
        mtr.abs_penalty = val;
    case {'smooth_std'}
        mtr.smooth_std = val;
    case {'angel_cutoff'}
        mtr.angle_cutoff = val;
    case {'shape_params_vec'}
        mtr.shape_params_vec = val;
end

return;


text:
function rad=sec2rad(sec)
% SEC2RAD  Converts seconds of arc to radians. Vectorized.
% Version: 2 Feb 98
% Useage:  rad=sec2rad(sec)
% Input:   sec - vector of angles in seconds of arc
% Output:  rad - vector of angles in radians

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

rad=sec.*pi./180./3600;

author: Sable

repo: mcbench-benchmarks


text:
function [TS_DataMat,Operations] = FilterFeatures(TS_DataMat,Operations,cfnParams)
% FilterFeatures filters down to a reduced feature set
%-------------------------------------------------------------------------------
if isempty(fieldnames(cfnParams)) || isempty(cfnParams.reducedFeatureSet) || strcmp(cfnParams.reducedFeatureSet,'all')
    return
else
    opIDs = TS_GiveMeFeatureSet(cfnParams.reducedFeatureSet,Operations);
    keepMe = ismember(Operations.ID,opIDs);
    Operations = Operations(keepMe,:);
    TS_DataMat = TS_DataMat(:,keepMe);
    fprintf(1,'Only using %u features from the %s set...\n',sum(keepMe),cfnParams.reducedFeatureSet);
end

% Check there are some left:
numFeatures = height(Operations);
if numFeatures==0
    error('No features to do classification with...');
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef (Abstract) ZmapResultsPlugin
    % base class for displaying results
    properties
       Parent
    end

    properties(Dependent)
        Data
        PlotTag
    end

    properties(Abstract)
    end

    methods(Abstract)
    end

    methods
        function r = get.Data(obj)
            r = obj.Parent.Result.values;
        end
        function r = get.PlotTag(obj)
            r = obj.Parent.PlotTag;
        end
    end
end

author: CelsoReyes