text:
function [Bt]=GetCoriolisMatrix1(T,Pcii,Icii,mcii,dq)
%% About the function: this is a function that is used to calculate
% Coriolis matrix of the serially linked manipulator, the return value of
% the function is (nxn) Coriolis matrix.

%% Arreguemnts:
% T is (4x4xn) transformation matrix of the serially linked robot, each
% (4x4) matrix represents the transform for each link in the base frame.
% Pcii is 3Xn matrix while each column represents the local coordinates
% of the center of mass of each link.
% Icii is (3x3xn) matrix, each 3x3 matrix of which represnets the
% associated link inertial tensor represented in its local inertial frame
% mcii is (1xn) vector, each element of which specifies a mass of one of
% the links

% Copyright Mohammad SAFEEA

n=max(size(mcii));
%% Initialization of Ai and Bi.
Bi=zeros(3,n,n);
Di=zeros(3,n,n);
L=zeros(3,3);
Lj_1=zeros(3,1);
wj=zeros(3,n);
half_wj=zeros(3,n);
%% Calculate === some auxuliary variables
Pcii_A=zeros(3,n);
w=zeros(3,n);
vci=zeros(3,n);
vi=zeros(3,n);
mcii_Pcii_A=zeros(3,n);
Pcii_A(:,1)=T(1:3,1:3,1)*Pcii(:,1);
w(:,1)=T(1:3,3,1)*dq(1);
wj(:,1)=w(:,1);
half_wj(:,1)=0.5*wj(:,1);
mcii_Pcii_A(:,1)=mcii(1)*Pcii_A(:,1);
double_kj=zeros(3,n);
double_kj(:,1)=2*T(1:3,3,1);
for i=2:n
        Pcii_A(:,i)=T(1:3,1:3,i)*Pcii(:,i);
        wj(:,i)=T(1:3,3,i)*dq(i);
        half_wj(:,i)=0.5*wj(:,i);
        w(:,i)=w(:,i-1)+wj(:,i);
        double_kj(:,i)=2*T(1:3,3,i);
        mcii_Pcii_A(:,i)=mcii(i)*Pcii_A(:,i);
end
%% calculating the links model, Mci and ddPci
for i=1:n
    %% calculating the Mci term
    Pci=Pcii_A(:,i)+T(1:3,4,i);
    L=T(1:3,1:3,i)*(trace(Icii(:,:,i))*eye(3)-2*Icii(:,:,i))*T(1:3,1:3,i)';
    Lj_1=L*T(1:3,3,i);
    for j=i:-1:2
    Bi(:,j,i)=Bi(:,j,i)+cross(Lj_1,half_wj(:,j));
    Lj_1=L*T(1:3,3,j-1);
    Bi(:,j-1,i)=Bi(:,j-1,i)+cross(Lj_1,w(:,i)-w(:,j-1));
    end
    j=1;
    %Bi(:,j,i)=Bi(:,j,i)+0.5*cross(Lj_1,wj(:,j));
    Bi(:,j,i)=Bi(:,j,i)+cross(Lj_1,half_wj(:,j));
    %% calculating the ddPci term
    vr=zeros(3,1);
    for j=i:-1:2
        Pcij=Pci-T(1:3,4,j);
        %Di(:,j,i)=Di(:,j,i)+(dq(j))*cross(T(1:3,3,j),cross(T(1:3,3,j),Pcij
        %));
        Di(:,j,i)=Di(:,j,i)+(dq(j))*(T(1:3,3,j)*(T(1:3,3,j)'*Pcij)-Pcij);
        vr=vr+cross(wj(:,j),Pcij);
        Di(:,j-1,i)=Di(:,j-1,i)+cross(double_kj(:,j-1),vr);
    end
    j=1;
    Pcij=Pci-T(1:3,4,j);
    Di(:,j,i)=Di(:,j,i)+(dq(j))*(T(1:3,3,j)*(T(1:3,3,j)'*Pcij)-Pcij);
end
Bt=zeros(n,n);
Fac_D=zeros(3,n);
Mac_B=zeros(3,n);
Pjp1_j=zeros(3,1);
%% calculating Mac for all of the links, then calculating two by filling At
%% and Bt

start=n-1;
j=n;
%% recursive rprocedure on moments and forces
        for k=1:n %% iterate through the matrix
            Mac_B(:,k)=Bi(:,k,j)+cross(mcii_Pcii_A(:,j),Di(:,k,j));
        end
        %% on forces
        Fac_D=mcii(j)*Di(:,:,j);
        Bt(j,:)=T(1:3,3,j)'*Mac_B;

for j=start:-1:1 %% iterate through the joints
    Pjp1_j=T(1:3,4,j+1)-T(1:3,4,j);
    %% recursive rprocedure on moments and forces
        for k=1:j %% iterate through the matrix
            Mac_B(:,k)=Mac_B(:,k)+Bi(:,k,j)+cross(Pjp1_j,Fac_D(:,k))+cross(mcii_Pcii_A(:,j),Di(:,k,j));
        end
        for k=j+1:n
            Mac_B(:,k)=Mac_B(:,k)+cross(Pjp1_j,Fac_D(:,k));
        end
    %% on forces
    Fac_D(:,1:j)=Fac_D(:,1:j)+mcii(j)*Di(:,1:j,j);
    Bt(j,:)=T(1:3,3,j)'*Mac_B;

end
%% close the function
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [ yval, ypval, yppval, left_out ] = spline_cubic_val2 ( ...
  n, t, y, ypp, left, tval )

%*****************************************************************************80
%
%% SPLINE_CUBIC_VAL2 evaluates a piecewise cubic spline at a point.
%
%  Discussion:
%
%    This routine is a modification of SPLINE_CUBIC_VAL; it allows the
%    user to speed up the code by suggesting the appropriate T interval
%    to search first.
%
%    SPLINE_CUBIC_SET must have already been called to define the
%    values of YPP.
%
%    In the LEFT interval, let RIGHT = LEFT+1.  The form of the spline is
%
%      SPL(T) =
%          A
%        + B * ( T - T(LEFT) )
%        + C * ( T - T(LEFT) )**2
%        + D * ( T - T(LEFT) )**3
%
%    Here:
%      A = Y(LEFT)
%      B = ( Y(RIGHT) - Y(LEFT) ) / ( T(RIGHT) - T(LEFT) )
%        - ( YPP(RIGHT) + 2 * YPP(LEFT) ) * ( T(RIGHT) - T(LEFT) ) / 6
%      C = YPP(LEFT) / 2
%      D = ( YPP(RIGHT) - YPP(LEFT) ) / ( 6 * ( T(RIGHT) - T(LEFT) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl de Boor,
%    A Practical Guide to Splines,
%    Springer Verlag, 1978.
%
%  Parameters:
%
%    Input, integer N, the number of knots.
%
%    Input, real T(N), the knot values.
%
%    Input, real Y(N), the data values at the knots.
%
%    Input, real YPP(N), the second derivatives of the spline at
%    the knots.
%
%    Input, integer LEFT, the suggested T interval to search.
%    LEFT should be between 1 and N-1.  If LEFT is not in this range,
%    then its value will be ignored.
%
%    Input, real TVAL, a point, typically between T(1) and T(N), at
%    which the spline is to be evalulated.  If TVAL lies outside
%    this range, extrapolation is used.
%
%    Output, real YVAL, YPVAL, YPPVAL, the value of the spline, and
%    its first two derivatives at TVAL.
%
%    Output, integer LEFT_OUT, the actual interval in which TVAL lies.
%
  left_out = left;
%
%  Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
%
%  What you want from R8VEC_BRACKET3 is that TVAL is to be computed
%  by the data in interval {T(LEFT), T(RIGHT)].
%
  left_out = r8vec_bracket3 ( n, t, tval, left_out );
  right = left_out + 1;
%
%  In the interval LEFT_OUT, the polynomial is in terms of a normalized
%  coordinate  ( DT / H ) between 0 and 1.
%
  dt = tval - t(left_out);
  h = t(right) - t(left_out);

  yval = y(left_out) + dt * ( ( y(right) - y(left_out) ) / h ...
              - ( ypp(right) / 6.0 + ypp(left_out) / 3.0 ) * h ...
       + dt * ( 0.5 * ypp(left_out) ...
       + dt * ( ( ypp(right) - ypp(left_out) ) / ( 6.0 * h ) ) ) );

  ypval = ( y(right) - y(left_out) ) / h ...
      - ( ypp(right) / 6.0 + ypp(left_out) / 3.0 ) * h ...
      + dt * ( ypp(left_out) ...
      + dt * ( 0.5 * ( ypp(right) - ypp(left_out) ) / h ) );

  yppval = ypp(left_out) + dt * ( ypp(right) - ypp(left_out) ) / h;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function corrderived_plot_swarm_history(c,varargin)

%CORRDERIVED_PLOT_SWARM_HISTORY plots relevent details of seismic swarms
% CORRDERIVED_PLOT_SWARM_HISTORY(C)reads correlation object C and plots various
% parameters to characterize the swarm. Typical input will be a correlation
% object containing waveforms from a single station for numerous events.
% This function requires that the CLUST and WAVEFORM fields of the
% correlation object be filled. The heavily lifting will have already been
% done by the LINKAGE and CLUSTER functions. This is predominantly a
% plotting function. Plots include:
%    - A histogram of all events and clustered events
%    - An amplitude measure of all events
%    - A plot showing the "lifespan" of each cluster
%
% CORRDERIVED_PLOT_SWARM_HISTORY(C,CLUSTERSIZE) specifies the minimum number of
% events to be counted as a cluster. The default is 5. Note that the
% clusters defined in the correaltion object may have a few as 1 event. For
% analysis purposes however, the user will often want to set a minimum
% size for significant clusters. CLUSTERSIZE is this value.
%
% The trace amplitude plot is based on the maximum value of the hilbert
% transform of the input data. The user may wish to narrow the window of
% data to a portion of the input trace by using the cop function. Example
%     c1 = crop(c,-1,3)
%     corrderived_plot_swarm_history(c1,10)
%
% See also correlation/linkage correlation/cluster

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% CHECK INPUTS
if nargin>2
    error('Incorrect number of inputs');
end
if ~isa(c,'correlation')
    error('First argument must be a correlation object');
end

if isempty(get(c,'CLUST'))
    error('CLUSTER field must be filled in input argument. See HELP CLUSTER');
end

if numel(varargin)==1
    clusterSize = varargin{1};
else
    clusterSize = 5;
end

% FIND EVENTS IN/OUT OF CLUSTERS
inMult = find(c,'BIG',clusterSize);
% find non-multiplets
for n = 1:get(c,'TRACES')
    if numel(find(inMult==n))==0
        notInMult(n) = 1;
    else
        notInMult(n) = 0;
    end
end
notInMult = find(notInMult)';
inMult = sort(inMult);

% GET STA_CHAN
if ~check(c,'STA') | ~check(c,'CHAN')
   warning('This function was written assuming data from a single sta_chan');
end
wtmp = waveform(c);
wtmp = wtmp(1);
nscl = [ get(wtmp,'NETWORK') '_' get(wtmp,'STATION') '_' get(wtmp,'CHANNEL') '_' get(wtmp,'LOCATION') ];

% PLOT IT
figure('Color','w','Position',[0 0 1100 850]);
box on; hold on;
set(gcf,'DefaultLineLineWidth',0.1);
set(gcf,'DefaultAxesFontSize',12);

% PLOT EVENT RATES
subplot(3,1,1)
disp('Preparing event rate histograms ...');
c1 = subset(c,inMult);
c1 = sort(c1);
trigMult = get(c1,'TRIG');
trigAll = get(c,'TRIG');
edges = [floor(min(trigAll)):1/24:ceil(max(trigAll))];
nMult = histc(trigMult,edges);
nAll = histc(trigAll,edges);
%
bar(edges,nAll,'y');
%h = findobj(gca,'Type','patch');
%set(h,'FaceColor',[.7 .7 .7])
hold on;
bar(edges,nMult,'r');
h = findobj(gca,'Type','patch');
set(h,'EdgeColor',[0 0 0],'LineWidth',0.1)
xlim([[min(trigAll)-1/24 max(trigAll)+1/24]]);
datetick('x','KeepLimits');
legend('All events','Clustered events');

title('Event rates');
ylabel('Events per hour');

% AMPLITUDE PLOT
disp('Preparing amplitude measures ...');
subplot(3,1,2)
c1 = crop(c,-1,4);
w = waveform(c1);
w = hilbert(w);
amp = max(double(w));
trig = get(c1,'TRIG');
plot(trig(notInMult),amp(notInMult),'ko','MarkerFaceColor','y','MarkerSize',4);
hold on;
plot(trig(inMult),amp(inMult),'ko','MarkerFaceColor','r','MarkerSize',4);
set(gca,'YScale','log');
xlim([[min(trigAll)-1/24 max(trigAll)+1/24]]);
ylim([min(amp) max(amp)]);
datetick('x','KeepLimits');
%set(gca,'YGrid','on')
ylabel('amplitude (nm/s)')
legend('All events','Clustered events');
title(['Event amplitudes    (derived from  ' nscl ')'],'Interpreter','none');

% PLOT EVENT RATES
subplot(3,1,3)
disp('Preparing cluster lifespan plot ...');
family = getclusterstat(c);
index = find(family.numel>=clusterSize);
[tmp,index] = sort(family.begin(index));

box on; hold on;
for n = 1:numel(index)
    plot([ family.begin(index(n)) family.finish(index(n)) ],[n n],'-','Color',[0.7 0.7 0.7],'LineWidth',1);
    plot(family.trig{index(n)},repmat(n,family.numel(index(n)),1),'o','Color',[0 0 0],'MarkerFaceColor','r','MarkerSize',4);

end
xlim([[min(trigAll)-1/24 max(trigAll)+1/24]]);
ylim([0 numel(index)+1]);
datetick('x','KeepLimits');
set(gca,'YTick',[1:numel(index)]);
set(gca,'YTickLabel',index);
ylabel('cluster rank')
title(['Cluster lifespan  (clusters of ' num2str(clusterSize) ' or more events)']);

%PRINT OUT FIGURE
set(gcf, 'paperorientation', 'portrait');
set(gcf, 'paperposition', [.25 .25 8 10.5] );
print(gcf,'-dpsc2','FIG_SWARM_HISTORY.ps');

author: geoscience-community-codes

repo: GISMO


text:
function test_suite = test_drawRect
%TESTDRAWRECT  One-line description here, please.
%   output = testDrawRect(input)
%
%   Example
%   testDrawRect
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-05-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSimple(testCase) %#ok<*DEFNU>

x0 = 20;
y0 = 30;
w = 40;
h = 50;

hf = 101;
figure(hf); clf;
hr = drawRect([x0 y0 w h]);
xth = [x0 x0+w x0+w x0 x0];
yth = [y0 y0 y0+h y0+h y0];
testCase.assertEqual(xth, get(hr, 'XData'), 'AbsTol', .01);
testCase.assertEqual(yth, get(hr, 'YData'), 'AbsTol', .01);

close(hf);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


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

%initialize q,
q = [0.1 0.1 0.1 0.1 0.1]
T = directkinematic(robot, q)
drawrobot3d(robot, q)

repo: ARTE


text:
function x = nonnegative( varargin ) %#ok

%NONNEGATIVE   The nonnegative orthant.
%   NONNEGATIVE(SX), where SX is a valid size vector, creates an array
%   of size SX and constrains each element to be nonnegative. Therefore,
%   given the declaration
%      variable x(sx)
%   the constraint
%      x == nonnegative(sx);
%   is equivalent to
%      x >= 0;
%   Obviously, the inequality form is simpler and is preferred in most
%   circumstances.
%
%   Disciplined convex programming information:
%       NONNEGATIVE is a cvx set specification. See the user guide for
%       details on how to use sets.

sx = cvx_get_dimlist( varargin ); %#ok
cvx_begin set
    variable x( sx ) nonnegative
cvx_end

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

author: yu-jiang

repo: radpbook

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


text:
function volume = checkSelectedVolume
% volume = checkSelectedVolume;
%
% checks for a selected volume view. If one exists, returns it;
% otherwise, initializes an empty volume view.
mrGlobals;
volume = getSelectedVolume;
if isempty(volume)
    % open a hidden view, but warn that you're doing this
    warning('No  VOLUME view found. Initializing a hidden volume view.')
    inplane = initHiddenVolume;
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function optionValueReturned = optionvalue(SpotPrice, StrikePrice, RiskFreeRate,TimeExpiry, Volatility,theOptionType,ButterflyRange)
% mcc -d compiled -B 'ccom:BSOptionModel,BSOptionModelClass,1' optionvalue.m webvizroutine.m
%CALCROUTINE Calculate the value of the option

OptionType = lower(theOptionType);

if (strcmpi(OptionType,'call') == 1) %Call option

     %Convert months to expiry to years to expiry
     TimeExpiry = TimeExpiry / 12;

     %Calculate the value of the option
     optionValueReturned = blsprice(SpotPrice, StrikePrice, RiskFreeRate, ...
          TimeExpiry, Volatility);

elseif (strcmpi(OptionType,'put') == 1) %Put option

     %Convert months to expiry to years to expiry
     TimeExpiry = TimeExpiry / 12;

     %Calculate the value of the put option
     optionValueReturned = blsprice(SpotPrice, StrikePrice, RiskFreeRate, ...
          TimeExpiry, Volatility);

elseif (strcmpi(OptionType,'straddle') == 1) %Straddle option

     %Convert months to expiry to years to expiry
     TimeExpiry = TimeExpiry / 12;

     %Compute the value of the straddle

     optionValueReturned = blsstrval(SpotPrice, StrikePrice, RiskFreeRate, ...
          TimeExpiry, Volatility);

elseif (strcmpi(OptionType,'butterfly') == 1) %Butterfly option

     %Convert months to expiry to years to expiry
     TimeExpiry = TimeExpiry / 12;

     %Compute the value of the butterfly
     optionValueReturned = blsbtyval(SpotPrice, StrikePrice, RiskFreeRate, ...
          TimeExpiry, Volatility, ButterflyRange);

end

%--------------------------------------------------------------------------

function StraddleValue = blsstrval(SpotPrice, StrikePrice, RiskFreeRate, ...
          TimeExpiry, Volatility)
%BLSSTRVAL Black Scholes value of a straddle option

%Calculate the value of both the call and put option
[CallValue, PutValue] = blsprice(SpotPrice, StrikePrice, RiskFreeRate, ...
     TimeExpiry, Volatility);

%Compute the value of the straddle
StraddleValue = CallValue + PutValue;

%end of BLSSTRVAL subroutine

function ButterflyValue = blsbtyval(SpotPrice, StrikePrice, RiskFreeRate, ...
          TimeExpiry, Volatility, ButterflyRange)
%BLSBTYVAL Black Scholes value of a butterfly option

%Set the different strike prices
LowStrike = StrikePrice .* (1 - ButterflyRange);
HighStrike = StrikePrice .* (1 + ButterflyRange);

%Value the long positions in the low and high struck calls
LowValue = blsprice(SpotPrice, LowStrike, RiskFreeRate, ...
     TimeExpiry, Volatility);
HighValue = blsprice(SpotPrice, HighStrike, RiskFreeRate, ...
     TimeExpiry, Volatility);

%Value the short position in the calls struck at the initial strike
%price
ShortValue = 2 .* -(blsprice(SpotPrice, StrikePrice, RiskFreeRate, ...
     TimeExpiry, Volatility));

%Calculate the total value of the butterfly
ButterflyValue = LowValue + HighValue + ShortValue;

%end of BLSBTYVAL subroutine

author: Sable

repo: mcbench-benchmarks


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [pathS1, pathS2] = MC_B_path(S0,r,d,sigmaB, sigmaBS,T,Z)
%
% Simulate a path within the Bachelier model model
%
%   S0 spot price
%   K the strike price
%   r the riskless rate
%   d the dividend yield
%   sigma the volatility
%   T the maturity
%   NTime number of timesteps
%   NSim number of simulations

NSim = size(Z,1);
NTime = size(Z,2);
Delta = T/NTime;                          % The discretisation step

lnS1 = zeros(NSim,NTime+1);               % init the logspot price path
lnS1(:,1)=log(S0*exp(-d*T));              % adjust due to dividend yield
S1 = zeros(NSim,NTime+1);
S1(:,1) = S0;

%dW = randn(NSim,NTime);                   % precompute all randoms

for i=1:NTime
    lnS1(:,i+1) = lnS1(:,i) + (r-d) * Delta + sigmaBS* sqrt(Delta)*Z(:,i);
    S1(:,i+1) = S1(:,i) + (r-d) * Delta + sigmaB * sqrt(Delta) * Z(:,i);
end

pathS1 = exp(lnS1);
pathS2 = S1;
clear dW;
end


text:
function pstruct = tapas_unitsq_sgm_mu3_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2018 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

return;

author: translationalneuromodeling

repo: tapas


text:
function [ipinfo,tt]=hostipinfo(m,t)
%HOSTIPINFO get host name and internet connection information
%  Usage:    hostname   = hostipinfo;        % computer name
%            ip_address = hostipinfo('i');   % IP address
%
%  Inputs:   m  mode input selects the information to output (see below)
%            t  optional output from an "ipconfig -all" command
%
%  Outputs:  ipinfo   cell array with one row per IP connection. The number
%                     of columns is detemined by the m input (see below).
%                     If only a single output value is selected by the m input, then
%                     ipinfo will be a string rather than a cell array.
%            tt       output from the "ipconfig -all" command
%
%  The "m" input selects whch IP connections to supply information about
%  (rows of IPINFO) and what information to supply (columns of IPINFO).
%
%  Row selection: A   include all connections
%                 L   Local Internet Connection [default if m non-empty]
%                 B   Bluetooth Connection
%                 P   PPP connection (e.g. a VPn connection)
%
%  Col selection: h   hostname (same for all connections) [default if m is empty]
%                 c   Connection type code: L=local, B=bluetooth, P=PPP or ?=other
%                 t   Connection type (e.g. "Ethernet" or "PPP")
%                 n   Connection name (e.g. "Local Area Connection")
%                 m   MAC address (e.g. "00-53-45-00-00-00"
%                 i   IP address (e.g. "192.168.1.1") [default if m non-empty]

%      Copyright (C) Mike Brookes 2012
%      Version: $Id: hostipinfo.m 1443 2012-02-22 09:10:03Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2
    t=evalc('!ipconfig -all');
    if nargin<1
        m='';
    end
end
tt=t;
if ~numel(m)
    m='Lh';
end
host=regexp(t,'Host Name[ .]*:\s*([\w-]*)','tokens','once'); % get host name
if ~numel(m)
    ipinfo=host{1};
else
    [iad,ad]=regexp(t,'\n(\w*) adapter ([ \w]*):','start','tokens'); % get adapters
    [iadx,adx]=regexp(t,'\n(\w* adapter [ \w]*):','start','tokens'); % get entire adapter names
    [imac,mac]=regexp(t,'Physical Address[ .]*:\s*(\S*)','start','tokens'); % get MAC addresses
    [iip,ip]=regexp(t,'IP Address[ .]*:\s*(\S*)','start','tokens'); % get IP addresses
    nad=numel(iad);  % number of adapters listed in ipconfig
    iad(nad+1)=length(t); % add an extra fake adapter at the end of the list
    fty='hctnmi'; % possible output fields
    nfty=numel(fty);    % number of fields
    ipinfo=cell(nad,nfty); % cell array to store the fields for each adapter
    nr=nad;
    aty='LBP';   % types of adapter
    naty=numel(aty);
    assx={'Ethernet adapter Local Area Connection'; 'Ethernet adapter Bluetooth Network Connection' ; 'PPP adapter '};
    for i=1:nad
        ipinfo{i,1}=host{1};
        ipinfo(i,3:4)=ad{i};
        ix=find((imac>iad(i)) & (imac<iad(i+1)),1);
        if numel(ix)
            ipinfo{i,5}=mac{ix}{1};
        else
            ipinfo{i,5}='';
        end
        ix=find((iip>iad(i)) & (iip<iad(i+1)),1);
        if numel(ix)
            ipinfo{i,6}=ip{ix}{1};
        else
            ipinfo{i,6}='';
        end
        jx=cell(length(assx),1);
        for j=1:length(assx)
            jx{j}=strfind(adx{i}{1},assx{j});    % compare adapter strings with adapter types
        end
        jx=find(~cellfun('isempty',jx));
        if numel(jx)
            ipinfo{i,2}=aty(jx);
        else
            ipinfo{i,2}='?';
        end
    end
    % Select the requested adapters from the rows of ipinfo{}
    if ~any(m=='A')
        ix=(1:naty)*(aty(ones(numel(m),1),:)'==m(ones(naty,1),:));
        ix=ix(ix>0);
        nr=numel(ix);
        if ~nr
            ix=1;  % default to 'L' = Local Internet
        end
        jx=ix;      % jx will hold first match
        for i=1:numel(ix)
            a=strmatch(aty(ix(i)),ipinfo(:,2));
            if numel(a)
                jx(i)=a(1);
            else
                jx(i)=0;
            end
        end
        jx=jx(jx>0);
        ipinfo=ipinfo(jx,:);
    end
    % Select the requested fields from the columns of ipinfo{}
    ix=(1:nfty)*(fty(ones(numel(m),1),:)'==m(ones(nfty,1),:));
    ix=ix(ix>0);
    nc=numel(ix);
    if ~nc
        ix=6;  % default to 'i' = IP address
    end
    ipinfo=ipinfo(:,ix);
    if nr<2 && nc<2
        ipinfo=ipinfo{1};  % just output a string if only one element requested
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function test_pull1248

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_read_header ft_read_data ft_read_event dataset2files

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/pull1248/1106miedo.vhdr');

%%

hdr = ft_read_header(filename);
dat = ft_read_data(filename);
evt = ft_read_event(filename); % this used to crash due to there not being a vmrk file

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Q_gppca = mohsst5_experiment_vbtdsepinit(data, D, maxiter, densityW, ...
                                             densityX, updatepseudow, ...
                                             updatepseudox)

if nargin < 3
  maxiter = 10;
end
if nargin < 6
  updatepseudow = true;
end
if nargin < 7
  updatepseudox = true;
end

[M,N] = size(data.observations);

% Divide data into test and train sets
testsize = 20; % testset size in percents
Itest = load(sprintf('/share/bayes/data/jaakko/mohsst5/ind%dtest.mat',testsize));
Itrain = load(sprintf('/share/bayes/data/jaakko/mohsst5/ind%dtrain.mat',testsize));
Ytest = data.observations;
Ytest(Itrain.Itrain) = nan;
Yobs = data.observations;
Yobs(Itest.Itest) = nan;

% Run VB PCA
disp('-- Run VB PCA --')
Q_vbpca = pca_full(Yobs,D,'maxiters',min(15,maxiter),'rotate2pca',true, ...
                   'algorithm','vb');

%tsplot(Q_vbpca.S)

% Run TDSEP
disp('-- Run TDSEP --')
sel = [0 1 2 3 5 7 10 15 20 25 30 40 50 70 100];
%sel = ceil(0:50);
%sel = ceil(exp(0:5));
R = tdsep2(Q_vbpca.S, sel);
R
[Q_vbpca.A, Q_vbpca.Av, Q_vbpca.S, Q_vbpca.Sv] = vbpca_rotate(R, Q_vbpca.A, ...
                                                  Q_vbpca.Av, Q_vbpca.S, ...
                                                  Q_vbpca.Sv);

%tsplot(Q_vbpca.S)

%return

% $$$ % Load VB PCA results with 80 components
% $$$ Q_vbpca = load('/home/alexilin/matlab/kaplan/vbresults80.mat');

disp('Initialize with VB PCA results');
Qinit.W = Q_vbpca.A(:,1:D);
Qinit.X = Q_vbpca.S(1:D,:);
Qinit.tau = 1/Q_vbpca.V;
Qinit.varW = zeros(size(Qinit.W));
Qinit.varX = zeros(size(Qinit.X));
for m=1:M
  v = diag(Q_vbpca.Av{m});
  Qinit.varW(m,:) = v(1:D);
end
for n=1:N
  v = diag(Q_vbpca.Sv{n});
  Qinit.varX(:,n) = v(1:D);
end
% Put the mean as the last component.. :|
Qinit.X(end,:) = 1;
Qinit.varX(end,:) = 0;
Qinit.W(:,end) = Q_vbpca.Mu;
Qinit.varW(:,end) = Q_vbpca.Muv;

% Covariance functions
% W distances in kilometers
switch 1
 case 1
  gpcovW = {@gpcovScale, @(logtheta,x1,x2) gpcov(logtheta,x1,x2,@sqdistEarth)};
  logthetaW = [log(2); log(1e3)];
  Wcov = 'se';
end
% X distances in days
switch 2
  case 2
   gpcovX = @gpcovPP;
   logthetaX = log(80);
   Xcov = 'cs';
end

for d=1:D
  covfunW{d} = gpcovW;
  covfunX{d} = gpcovX;
  initthetaW{d} = logthetaW;
  initthetaX{d} = logthetaX;
end

datestring = datestr(now, 'yyyymmdd');
filename = sprintf(['/share/bayes/data/jaakko/mohsst5/' ...
                    'mohsst5_vbtdsepinit_testsize=%d_D=%d_Wcov=%s_Wden=' ...
                    '%d_Xcov=%s_Xden=%d_iter=%d_date=%s'], testsize, D, ...
                   Wcov, ceil(100*densityW), Xcov, ceil(100*densityX), ...
                   maxiter, datestring)

disp('-- Run GP PCA --')
Q_gppca = vbgppcamv(Yobs, D,...
                    data.coordinates, data.time,...
                    covfunW, initthetaW,...
                    covfunX,initthetaX, ...
                    'init', Qinit, ...
                    'maxiter',maxiter, ...
                    'pseudodensityx', densityX, ...
                    'pseudodensityw', densityW, ...
                    'loglikelihood', true, ...
                    'updatehyper', [1 1 10:10:(maxiter-1)], ...
                    'updatepseudox', updatepseudox, ...
                    'updatepseudow', updatepseudow, ...
                    'maxsearchx', 3, ...
                    'maxsearchw', 3, ...
                    'checkgradx', false, ...
                    'autosavetime', 3600, ...
                    'autosavefile', filename);

Yh_gppca = Q_gppca.W * Q_gppca.X;
%Yh_vbpca = Q_vbpca.A * Q_vbpca.S + repmat(Q_vbpca.Mu,1,N);

trainrmse_gppca = rmse(Yobs, Yh_gppca)
%trainrmse_vbpca = rmse(Yobs, Yh_vbpca)
testrmse_gppca = rmse(Ytest, Yh_gppca)
%testrmse_vbpca = rmse(Ytest, Yh_vbpca)

% Huuuuuge matrices... Don't save them..
Q_gppca.CovXp = [];
Q_gppca.CovWp = [];

Q_gppca.filename = filename;

%save(filename, '-struct', 'Q_gppca');
save(filename, '-struct', 'Q_gppca');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [a] = tcoeffs(X,P,window,weight,nModes)
%TCOEFFS Continuously-discrete temporal expansion coefficients of SPOD
%modes
%   [A] = TCOEFFS(X,P,WINDOW,WEIGHT,NMODES) returns the
%   continuously-discrete temporal SPOD mode expansion coefficients of the
%   leading NMODES modes. P is the data matrix of SPOD modes returned by
%   SPOD. X, WINDOW and WEIGHT are the same variables as for SPOD. If no
%   windowing function is specified, a Hamming window of length WINDOW will
%   be used. If WEIGHT is empty, a uniform weighting of 1 is used.
%
%   Reference:
%     [1] A. Nekkanti, O. T. Schmidt, Frequency–time analysis, low-rank
%         reconstruction and denoising of turbulent flows using SPOD,
%         Journal of Fluid Mechanics 926, A26, 2021
%
% A. Nekkanti (aknekkan@eng.ucsd.edu), O. T. Schmidt (oschmidt@ucsd.edu)
% Last revision: 7-Oct-2022 Brandon Yeung <byeung@ucsd.edu>

dims        = size(X);
nt          = dims(1);
nGrid       = prod(dims(2:end));

window = window(:); weight = weight(:);

% default window size and type
if length(window)==1
    window  = hammwin(window);
end

nDFT        = length(window);
winWeight   = 1/mean(window);

% inner product weight
if isempty(weight)
    weight  = ones(nGrid,1);
end

X           = reshape(X,nt,nGrid);
X           = X-mean(X,1);              % subtract mean

ndims       = size(P);
P           = permute(P,[1 length(ndims) 2:length(ndims)-1]);
if isreal(X)
    nFreq   = ceil(nDFT/2)+1;
    P       = reshape(P,ceil(nDFT/2)+1,ndims(end),nGrid);
else
    nFreq   = nDFT;
    P       = reshape(P,nDFT,ndims(end),nGrid);
end

weight     = reshape(weight,1,nGrid);
% zero-padding
X          = [zeros(ceil(nDFT/2),nGrid); X; zeros(ceil(nDFT/2),nGrid);];
a          = zeros(nFreq,nModes,nt);
winCorr_fac= winWeight/nDFT;
disp(' ')
disp('Calculating expansion coefficients')
disp('------------------------------------')
for i=1:nt
    X_blk               = fft(X(i:i+nDFT-1,:).*window);
    X_blk               = X_blk(1:nFreq,:);
    % correction for windowing and zero-padding
    if (i<ceil(nDFT/2)+1)
        corr 	= 1/(winCorr_fac*sum(window(ceil(nDFT/2)-i+1:nDFT)));
    elseif (i>nt-ceil(nDFT/2)+1)
        corr	= 1/(winCorr_fac*sum(window(1:nt+ceil(nDFT/2)-i)));
    else
        corr 	= 1;
    end
    for l=1:nModes
        a(:,l,i) = corr*winCorr_fac*dot(squeeze(P(:,l,:)),weight.*X_blk,2);
    end
    disp(['time ' num2str(i) '/' num2str(nt)])
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [window] = hammwin(N)
%HAMMWIN Standard Hamming window of lenght N
    window = 0.54-0.46*cos(2*pi*(0:N-1)/(N-1))';
end

author: SpectralPOD


text:
function w = quadwts(n)
%QUADWTS   Quadrature weights for Chebyshev points of 2nd kind.
%   QUADWTS(N) returns the N weights for Clenshaw-Curtis quadrature on 2nd-kind
%   Chebyshev points.
%
% See also CHEBPTS, BARYWTS.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPER NOTE:
% We use a variant of Waldvogel's algorithm [1], due to Nick Hale. (See below)
% We note this is similar to Greg Von Winkel's approach, which can be found on
% the MathWorks File Exchange.
%
% Let $f(x) = \sum_{k=0}^nc_kT_k(x)$, then\vspace*{-3pt} }
%   I(f) = v.'*c
% where
%   v = \int_{-1}^1T_k(x)dx = { 2/(1-k^2) : k even
%                             { 0         : k odd
%     = v'*inv(TT)*f(x) where TT_{j,k} = T_k(x_j)
%     = (inv(TT)'*v)'*f(x)
% Therefore
%   I(f) = w.'f(x) => w = inv(TT).'*v;
% Here inv(TT).' = inv(TT) is an inverse discrete cosine transform of Type I.
%
% Furthermore, since odd entries in v are zero, can compute via FFT without
% doubling up from N to 2N (though we still need to double up from N/2 to N to
% facilitate the use of ifft).
%
% References:
%   [1] Joerg Waldvogel, "Fast construction of the Fejer and Clenshaw-Curtis
%       quadrature rules", BIT Numerical Mathematics 46 (2006), pp 195-202.
%   [2] Greg von Winckel, "Fast Clenshaw-Curtis Quadrature",
%       http://www.mathworks.com/matlabcentral/fileexchange/6911, (2005)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ( n == 0 )                      % Special case (no points!)
    w = [];
elseif ( n == 1 )                  % Special case (single point)
    w = 2;
else                               % General case
    c = 2./[1, 1-(2:2:(n-1)).^2];  % Exact integrals of T_k (even)
    c = [c, c(floor(n/2):-1:2)];   % Mirror for DCT via FFT
    w = ifft(c);                   % Interior weights
    w([1,n]) = w(1)/2;             % Boundary weights
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function c = mtimes(a,b)
% implements e * c and c * e
%
% Description
% If |e| is a matrix of embeddings and |c| is a matrix of coefficients
% then |e * c| is again a matrix of embeddings defined by
%
% $$ [\mathrm{e * c}]_{j\ell} = \sum_{k} \mathrm{e}_{jk} \mathrm{c}_{k \ell}$$
%
% Syntax
%   out = e * c
%   out = c * e
%
% Input
%  e - @embedding
%  c - double
%
% Output
%  out- @embedding
%

if isa(a,'embedding')
  for i = 1:length(a.u), a.u{i} = a.u{i} * b; end
  c = a;
else
  for i = 1:length(b.u), b.u{i} = a * b.u{i}; end
  c = b;
end

end

author: mtex-toolbox

repo: mtex

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


text:
function gray = grayFilter(img, varargin)
%GRAYFILTER Compute configuration map of a binary image
%
%   GRAY = grayFilter(IMG);
%   Returns a gray-scale image, with size dim(IMG)-1, containing values of
%   the 2x2 configuration of the original binary image.
%
%   Example:
%   img = [0 0 0 0 0;0 1 1 0 0;0 1 1 1 0;0 0 0 0 0];
%   grayFilter(img)
%   ans =
%      8    12     4     0
%     10    15    13     4
%      2     3     3     1
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 22/10/2004.
%

%   HISTORY
%   03/11/2004 : add 3x3 2D case.

% pre-processing
img = img~=0;
dim = size(img);
nd = length(dim);

% size of neighborhood to consider
nu=1;

% extract input parameters
if ~isempty(varargin)
    var = varargin{1};
    if length(var)==1
        nu = var;
    else
        coef = var;
        % find size of filter from length of coefficients.
        nu = power(log(length(coef))/log(2), 1/nd);
    end
end

if length(dim)==2
    % 2 dimensions
    NY = dim(1)-nu;
    NX = dim(2)-nu;

    if nu==1
        gray=1*img(1:NY,    1:NX) + ...
             2*img(1:NY,    2:NX+1) + ...
             4*img(2:NY+1,  1:NX) + ...
             8*img(2:NY+1,  2:NX+1) ;
    elseif nu==2
        gray =  1*img(1:NY,     1:NX) + ...
                2*img(1:NY,     2:NX+1) + ...
                4*img(1:NY,     3:NX+2) + ...
                8*img(2:NY+1,   1:NX) + ...
               16*img(2:NY+1,   2:NX+1) + ...
               32*img(2:NY+1,   3:NX+2) + ...
               64*img(3:NY+2,   1:NX) + ...
              128*img(3:NY+2,   2:NX+1) + ...
              256*img(3:NY+2,   3:NX+2) ;
    end

else
    % 3 dimensions
    NY = dim(1)-nu;
    NX = dim(2)-nu;
    NZ = dim(3)-nu;

    gray=   img(1:NY,    1:NX,   1:NZ) + ...
          2*img(1:NY,    2:NX+1, 1:NZ) + ...
          4*img(2:NY+1,  1:NX,   1:NZ) + ...
          8*img(2:NY+1,  2:NX+1, 1:NZ) + ...
         16*img(1:NY,    1:NX,   2:NZ+1) + ...
         32*img(1:NY,    2:NX+1, 2:NZ+1) + ...
         64*img(2:NY+1,  1:NX,   2:NZ+1) + ...
        128*img(2:NY+1,  2:NX+1, 2:NZ+1);
end

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [nll,g,H,T] = LogisticLoss(w,X,y)
% w(feature,1)
% X(instance,feature)
% y(instance,1)

[n,p] = size(X);

Xw = X*w;
yXw = y.*Xw;

nll = sum(	([zeros(n,1) -yXw]));

if nargout > 1
    if nargout > 2
        sig = 1./(1+exp(-yXw));
        g = -X.'*(y.*(1-sig));
    else
        g = -X.'*(y./(1+exp(yXw)));
    end
end

if nargout > 2
    H = X.'*diag(sparse(sig.*(1-sig)))*X;
end

if nargout > 3
    T = zeros(p,p,p);
    for j1 = 1:p
        for j2 = 1:p
            for j3 = 1:p
                T(j1,j2,j3) = sum(y(:).^3.*X(:,j1).*X(:,j2).*X(:,j3).*sig.*(1-sig).*(1-2*sig));
            end
        end
    end
end

author: singaxiong

repo: SignalGraph


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

function result = areFramesWithinRange(robot_data, frame_a, frame_b, range)

pose_a = ...
    robot_data.Sim_O_C{frame_a};
pos_a = pose_a(1:3, 4);
pose_b = ...
    robot_data.Sim_O_C{frame_b};
pos_b = pose_b(1:3, 4);

result = norm(pos_a - pos_b) < range;

end

author: uzh-rpg


text:
% Compute the pdf, cdf and rank distributions for a sequence of numbers.
%
% INPUTS: sequence of values: x, (1xn), 'plt' - 'on' or 'off'
%         bin - bin size, default is [], then it gets selected
%         automatically
% OUTPUTS: pdf, cdf and rank distribution values, plot is optional
%
% Note: pdf = frequency distribution, cdf = cumulative frequency,
%       rank = log-log scale of the sorted sequence
% GB: last updated, November 24 2012

function [xpdf,ypdf,xcdf,ycdf,logk,logx]=pdfCdfRank(x,plt='off',bin=[])

xx=unique(x);

if length(bin)==0; bin = length(xx)/5; end

for ii=1:numel(xx)
  xcdf(ii) = xx(ii);
  ycdf(ii) = length(find(x<=xx(ii)))/numel(x);

  % how many x's fall in the interval [xx(ii)-0.5*numel(xx)/bin,xx(ii)+0.5*numel(xx)/bin]
  xpdf(ii) = xx(ii);
  ypdf(ii) = length(find(abs(xx(ii)-x)<=0.5*numel(xx)/bin))/numel(x);
end

x=-sort(-x);
logk=log(1:length(x));
logx=log(x);

if strcmp(plt,'on')
     set(gcf,'color',[1,1,1])

     subplot(2,2,1)
     hist(x,bin,'facecolor',[0.5 0.5 0.5])
     title('histogram')

     subplot(2,2,3)
     plot(xpdf,ypdf,'k.')
     title('pdf')
     axis('tight')

     subplot(2,2,2)
     plot(xcdf,ycdf,'k.')
     title('cdf')
     axis('tight')

     subplot(2,2,4)
     plot(logk,logx,'k.')
     title('rank')
     axis('tight')
end

author: aeolianine

repo: octave-networks-toolbox


text:
function chebvar(varargin)
%CHEBVAR   Short-cut for constructing CHEBFUN variables.
%   CHEBVAR arg1 arg2 ...
%   is short-hand notation for creating symbolic variables
%          arg1 = chebfun('arg1');
%          arg2 = chebfun('arg2'); ...
%   The outputs are created in the current workspace.
%
%   CHEBVAR arg1 arg2 ... DOM constructs the CHEBFUN objects on the domain DOM,
%   i.e.,  arg1 = chebfun('arg1', DOM);
%          arg2 = chebfun('arg2', DOM); ...
%
%   In both cases, the CHEBFUN is created according to the currently stored
%   default preferences.
%
%   Example:
%     chebvar x
%     f = sin(x)
%
% See also CHEBFUN, CHEBFUNPREF.

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

% Trivial case:
if ( nargin == 0 )
    return
end
dom = [];

% Locate valid variable names:
isVar = cellfun(@isvarname, varargin);

% The last entry may be a domain:
if ( ~all(isVar) )
    dom = str2num(varargin{end}); %#ok<ST2NM>
    if ( ~isempty(dom) )
        varargin(end) = [];
        isVar(end) = [];
    end
end

% Check validity of variable name:
if ( ~all(isVar) )
    error('CHEBFUN:chebvar:badName', 'Not a valid variable name.');
end

% Acquire some preferences:
pref = chebfunpref();
if ( isempty(dom) )
    dom = pref.domain;
end

% Loop over each of the inputs:
for k = numel(varargin):-1:1
    op = str2op(varargin{k});
    f = chebfun(op, dom, pref);
    assignin('caller', varargin{k}, f);
end

end

function op = str2op(op)
% This is here as it's a clean function with no other variables hanging around
% in the scope.
depVar = symvar(op);
if ( numel(depVar) ~= 1 )
    error('CHEBFUN:chebvar:indepVars', ...
        'Incorrect number of independent variables in string input.');
end
op = eval(['@(' depVar{:} ')', op]);
end


text:
function [ str2 ] = newstr( str, pos )
% Function to a new string by flipping a single bit
% By Kyriakos Tsourapas
% You may contact me through the Mathworks site
% University of Essex 2002

str2 = str;

% FLIP THE POSITION
if str2(pos) == 0
    str2(pos) = 1;
else
    str2(pos) = 0;
end


text:
function [policy, iter, cpu_time] = mdp_value_iteration(P, R, discount, epsilon, max_iter, V0)

% mdp_value_iteration   Resolution of discounted MDP with value iteration algorithm
% Arguments --------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA)  = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a matrix (SxS) possibly sparse
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount  = discount rate in ]0; 1]
%               beware to check conditions of convergence for discount = 1.
%   epsilon   = epsilon-optimal policy search, upper than 0,
%               optional (default : 0.01)
%   max_iter  = maximum number of iteration to be done, upper than 0,
%               optional (default : computed)
%   V0(S)     = starting value function, optional (default : zeros(S,1))
% Evaluation -------------------------------------------------------------
%   policy(S) = epsilon-optimal policy
%   iter      = number of done iterations
%   cpu_time  = used CPU time
%--------------------------------------------------------------------------
% In verbose mode, at each iteration, displays the variation of V
% and the condition which stopped iterations: epsilon-optimum policy found
% or maximum number of iterations reached.

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

cpu_time = cputime;

global mdp_VERBOSE;

% check of arguments
if iscell(P); S = size(P{1},1); else S = size(P,1); end;
if discount <= 0 || discount > 1
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0; 1]')
    disp('--------------------------------------------------------')
elseif nargin > 3 && (epsilon <= 0)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: epsilon must be upper than 0')
    disp('--------------------------------------------------------')
elseif nargin > 4 && max_iter <= 0
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: The maximum number of iteration must be upper than 0')
    disp('--------------------------------------------------------')
elseif nargin > 5 && size(V0,1) ~= S
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: V0 must have the same dimension as P')
    disp('--------------------------------------------------------')
else

    if discount == 1
        disp('--------------------------------------------------------')
        disp('MDP Toolbox WARNING: check conditions of convergence.')
        disp('With no discount, convergence is not always assumed.')
        disp('--------------------------------------------------------')
    end;

    PR = mdp_computePR(P,R);

    % initialization of optional arguments
    if nargin < 6; V0 = zeros(S,1); end;
    if nargin < 4; epsilon = 0.01; end;
    % compute a bound for the number of iterations
    if discount ~= 1
       computed_max_iter = mdp_value_iteration_bound_iter(P, R, discount, epsilon, V0);
    end;
    if nargin < 5
        if discount ~= 1
            max_iter = computed_max_iter;
        else
            max_iter = 1000;
        end;
    else
        if discount ~= 1 && max_iter > computed_max_iter
            disp(['MDP Toolbox WARNING: max_iter is bounded by ' num2str(computed_max_iter,'%12.1f') ])
            max_iter = computed_max_iter;
        end;
    end;

    % computation of threshold of variation for V for an epsilon-optimal policy
    if discount ~= 1
        thresh = epsilon * (1-discount)/discount;
    else
        thresh = epsilon;
    end;

    if mdp_VERBOSE; disp('  Iteration    V_variation'); end;

    iter = 0;
    V = V0;
    is_done = false;
    while ~is_done
        iter = iter + 1;
        Vprev = V;

        [V, policy] = mdp_bellman_operator(P,PR,discount,V);

        variation = mdp_span(V - Vprev);
        if mdp_VERBOSE;
            disp(['      ' num2str(iter,'%5i') '         ' num2str(variation)]);
        end;
        if variation < thresh
            is_done = true;
            if mdp_VERBOSE
                disp('MDP Toolbox: iterations stopped, epsilon-optimal policy found')
            end;
        elseif iter == max_iter
            is_done = true;
            if mdp_VERBOSE
                disp('MDP Toolbox: iterations stopped by maximum number of iteration condition')
            end;
        end;
    end;

end;

cpu_time = cputime - cpu_time;


text:
function w = wave2gray(c, s, scale, border)
%WAVE2GRAY Display wavelet decomposition coefficients.
%   W = WAVE2GRAY(C, S, SCALE, BORDER) displays and returns a
%   wavelet coefficient image.
%
%   EXAMPLES:
%     wave2gray(c, s);                      Display w/defaults.
%     foo = wave2gray(c, s);                Display and return.
%     foo = wave2gray(c, s, 4);             Magnify the details.
%     foo = wave2gray(c, s, -4);            Magnify absolute values.
%     foo = wave2gray(c, s, 1, 'append');   Keep border values.
%
%   INPUTS/OUTPUTS:
%     [C, S] is a wavelet decomposition vector and bookkeeping
%     matrix.
%
%     SCALE       Detail coefficient scaling
%     ----------------------------------------------------------
%     0 or 1      Maximum range (default)
%     2,3...      Magnify default by the scale factor
%     -1, -2...   Magnify absolute values by abs(scale)
%
%     BORDER      Border between wavelet decompositions
%     ----------------------------------------------------------
%     'absorb'    Border replaces image (default)
%     'append'    Border increases width of image
%
%     Image W:   ------- ------ ------------ -------------------
%                |      |      |            |
%                | a(n) | h(n) |            |
%                |      |      |            |
%                ------- ------     h(n-1)  |
%                |      |      |            |
%                | v(n) | d(n) |            |        h(n-2)
%                |      |      |            |
%                ------- ------ ------------
%                |             |            |
%                |    v(n-1)   |   d(n-1)   |
%                |             |            |
%                -------------- ------------ -------------------
%                |                          |
%                |          v(n-2)          |        d(n-2)
%                |                          |
%
%     Here, n denotes the decomposition step scale and a, h, v, d are
%     approximation, horizontal, vertical, and diagonal detail
%     coefficients, respectively.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/11/21 15:01:44 $

% Check input arguments for reasonableness.
error(nargchk(2, 4, nargin));

if (ndims(c) ~= 2) | (size(c, 1) ~= 1)
  error('C must be a row vector.');   end

if (ndims(s) ~= 2) | ~isreal(s) | ~isnumeric(s) | (size(s,2) ~= 2)
  error('S must be a real, numeric two-column array.');   end

elements = prod(s, 2);
if (length(c) < elements(end)) | ...
      ~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
   error(['[C S] must be a standard wavelet ' ...
          'decomposition structure.']);
end

if (nargin > 2) & (~isreal(scale) | ~isnumeric(scale))
   error('SCALE must be a real, numeric scalar.');
end

if (nargin > 3) & (~ischar(border))
  error('BORDER must be character string.');
end

if nargin == 2
   scale = 1;  % Default scale.
end

if nargin < 4
   border = 'absorb';  % Default border.
end

% Scale coefficients and determine pad fill.
absflag = scale < 0;
scale = abs(scale);
if scale == 0
   scale = 1;
end

[cd, w] = wavecut('a', c, s);   w = mat2gray(w);
cdx = max(abs(cd(:))) / scale;
if absflag
   cd = mat2gray(abs(cd), [0, cdx]);   fill = 0;
else
   cd = mat2gray(cd, [-cdx, cdx]);   fill = 0.5;
end

% Build gray image one decomposition at a time.
for i = size(s, 1) - 2:-1:1
   ws = size(w);

   h = wavecopy('h', cd, s, i);
   pad = ws - size(h);     frontporch = round(pad / 2);
   h = padarray(h, frontporch, fill, 'pre');
   h = padarray(h, pad - frontporch, fill, 'post');

   v = wavecopy('v', cd, s, i);
   pad = ws - size(v);     frontporch = round(pad / 2);
   v = padarray(v, frontporch, fill, 'pre');
   v = padarray(v, pad - frontporch, fill, 'post');

   d = wavecopy('d', cd, s, i);
   pad = ws - size(d);     frontporch = round(pad / 2);
   d = padarray(d, frontporch, fill, 'pre');
   d = padarray(d, pad - frontporch, fill, 'post');

   % Add 1 pixel white border.
   switch lower(border)
   case 'append'
      w = padarray(w, [1 1], 1, 'post');
      h = padarray(h, [1 0], 1, 'post');
      v = padarray(v, [0 1], 1, 'post');
   case 'absorb'
      w(:, end) = 1;   w(end, :) = 1;
      h(end, :) = 1;   v(:, end) = 1;
   otherwise
      error('Unrecognized BORDER parameter.');
   end

   w = [w h; v d];                 % Concatenate coefs.
end

if nargout == 0
   imshow(w);                      % Display result.
end

repo: weiyanmin


text:
function designs = get_design_set(name, d_0)
%GET_DESIGN_SET Gets a set of arrays designs for simulations.
%Inputs:
%   name - 'SameAperture': four designs with the same aperture.
%          'SameNumSensor': four designs with the same number of sensors.
%          'Nested10': seven nested arrays with the same number of sensors.
%   d_0 - Minimal inter-element spacing.
%Output:
%   designs - A cell array of array designs.
switch lower(name)
    case 'sameaperture'
        % same aperture
        designs = { ...
            design_array_1d('coprime', [2 3], d_0, '2M', 'Co-prime (2,3)') ...
            design_array_1d('custom', [0 1 2 6 9]*d_0, d_0, 'MRA 5') ...
            design_array_1d('nested', [1 5], d_0, 'Nested (1,5)') ...
            design_array_1d('nested', [4 2], d_0, 'Nested (4,2)') ...
        };
    case 'samenumsensor'
        designs = { ...
            design_array_1d('coprime', [3 5], d_0, '2M', 'Co-prime (3,5)') ...
            design_array_1d('custom', [0 1 4 10 16 22 28 30 33 35]*d_0, d_0, 'MRA 10') ...
            design_array_1d('nested', [4 6], d_0, 'Nested (4,6)') ...
            design_array_1d('nested', [3 7], d_0, 'Nested (3,7)') ...
        };
    case 'nested10'
        designs = { ...
            design_array_1d('nested', [2 8], d_0, 'Nested (2,8)') ...
            design_array_1d('nested', [3 7], d_0, 'Nested (3,7)') ...
            design_array_1d('nested', [4 6], d_0, 'Nested (4,6)') ...
            design_array_1d('nested', [5 5], d_0, 'Nested (5,5)') ...
            design_array_1d('nested', [6 4], d_0, 'Nested (6,4)') ...
            design_array_1d('nested', [7 3], d_0, 'Nested (7,3)') ...
            design_array_1d('nested', [8 2], d_0, 'Nested (8,2)') ...
        };
    otherwise
        error('Unknown design set ''%s''.', name);
end
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function layer = genNetworkLSTMEncoderDecoderScale(para)

layer = genNetworkLSTMEncoderDecoder(para);

layerScale.name = 'Input';
layerScale.inputIdx = 3;
layerScale.dim = [1 1];
layerScale.next = 1;

layer = [layer(1:end-1) layerScale layer(end)];

layer{end-3}.next = 3;
layer{end-2}.next = 2;
layer{end}.prev = [-3 -2 -1];
% layer{end} = rmfield(layer{end}, 'costFrameSelection');

end


text:
function S=streakbar(X,Y,U,V,unit)

% H=streakbar(X,Y,U,V,unit) creates a colorbar for (but not exclusively)
% the function streakarrow.
% The arrays X and Y defines the coordinates for U and V.
% U and V are the same arrays used for streakarrow.
% The string variable unit is the unit of the vector magnitude
% Example:
%   streakbar(X,Y,U,V,'m/s')

Vmag=sqrt(U.^2+V.^2);
Vmin=min(Vmag(:)); Vmax=max(Vmag(:));

 P=get(gca,'position');
 %axes('position',[P(1)+P(3)+.02  P(2)+0.01  .01  P(4)-0.02]')
 axes('position',[P(1)+P(3)+.02  P(2)  .01  P(4)]')
 [X,Y]=meshgrid( [0 1], linspace(Vmin,Vmax,64));
 Q= [1:64; 1:64];
 S=pcolor(X', Y',Q); shading flat; set(gca,'XTickLabel',[],  'Yaxislocation', 'right')
 title(unit)


text:
function [success] = downsampleHdf5Movie(inputFilePath, varargin)
	% Downsamples the Hdf5 movie in inputFilePath piece by piece and appends it to already created hdf5 file.
	% Biafra Ahanonu
	% started: 2013.12.19
	% base append code based on work by Dinesh Iyer (http://www.mathworks.com/matlabcentral/newsreader/author/130530)
	% note: checked between this implementation and imageJ's scale, they are nearly identical (subtracted histogram mean is 1 vs each image mean of ~1860, so assuming precision error).
	% inputs
		% inputFilePath
	% options
		% datasetName = hierarchy where data is stored in HDF5 file

	[success] = ciapkg.hdf5.downsampleHdf5Movie(inputFilePath,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [c, ceq ] = multiFlybyNonlcon(x, fitnessfcn, minRadiiSingle, maxRadiiSingle, minXferRad, maxDepartVInf, maxArriveVInf, maxDeltaVSingle)
%multiFlybyNonlcon Summary of this function goes here
%   Detailed explanation goes here
    ceq = [];

    numPop = size(x,1);

    minRadii = [];
    maxRadii = [];
    maxDeltaV = [];
%     minRadii = NaN(1, numPop*numel(minRadiiSingle));
%     maxRadii = NaN(1, numPop*numel(maxRadiiSingle));
    for(i=1:length(minRadiiSingle))
        minRadii = [minRadii, minRadiiSingle(i) * ones(1, numPop)]; %#ok<AGROW>
        maxRadii = [maxRadii, maxRadiiSingle(i) * ones(1, numPop)]; %#ok<AGROW>
        maxDeltaV = [maxDeltaV, maxDeltaVSingle(i)*ones(1, numPop)]; %#ok<AGROW>
    end

    [~, rp, ~, ~, deltaVVect, vInfDNorm, ~, vInfArrive, ~, ~, ~, ~, ~, xferRp] = fitnessfcn(x);

    deltaV = sqrt(sum(abs(deltaVVect).^2,1));

    c = (minRadii(1:length(rp)) - rp)';
    c = reshape(c,numPop,size(c,1)/numPop);

    c2 = (rp - maxRadii(1:length(rp)))';
    c2 = reshape(c2,numPop,size(c2,1)/numPop);

    c3 = -(xferRp - minXferRad)'; %need the minus sign to make greater than values negative for the constraint
    c3 = reshape(c3,numPop,size(c3,1)/numPop);

%     c = horzcat(c,c2,c3);

    c4 = vInfDNorm - maxDepartVInf;
    c5 = (sqrt(sum(abs(vInfArrive).^2,1)))' - maxArriveVInf;

    c6 = (deltaV - maxDeltaV)';
    c6 = reshape(c6,numPop,size(c6,1)/numPop);

    c = horzcat(c,c2,c3,c4,c5,c6);

    c(isnan(c)) = 1;
    c(c == -Inf) = -realmax;
    c(c == Inf) = realmax;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function prior = truncatedGammaPriorExpandParam(prior, params)

% TRUNCATEDGAMMAPRIOREXPANDPARAM Expand truncated gamma prior structure from params.

% PRIOR

prior.a = params(1);
prior.b = params(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function binom_test ( )

%*****************************************************************************80
%
%% BINOM_TEST tests R4_BINOM and R8_BINOM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BINOM_TEST:\n' );
  fprintf ( 1, '  Test BINOM_VALUES, R4_BINOM, R8_BINOM.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             A               B     BINOM(A,B)\n' );
  fprintf ( 1, '                                R4_BINOM(A,B)         Diff\n' );
  fprintf ( 1, '                                R8_BINOM(A,B)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, fx1 ] = binomial_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_binom ( a, b );
    fx3 = r8_binom ( a, b );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14d  %14d  %14d\n', a, b, fx1 );
    fprintf ( 1, '                                  %14d  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                                  %14d  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%MDL_IRB140 Create model of ABB IRB 140 manipulator
%
% MDL_IRB140 is a script that creates the workspace variable robot which
% describes the kinematic characteristics of an ABB IRB 140 manipulator
% using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%   qd         lower arm horizontal as per data sheet
%
% Reference::
% - "IRB 140 data sheet", ABB Robotics.
% - "Utilizing the Functional Work Space Evaluation Tool for Assessing a
%   System Design and Reconfiguration Alternatives"
%   A. Djuric and R. J. Urbanic
%
% Notes::
% - SI units of metres are used.
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also SerialLink, mdl_fanuc10l, mdl_m16, mdl_motormanHP6, mdl_S4ABB2p8, mdl_puma560.

% MODEL: ABB, IRB140, 6DOF, standard_DH

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function r = mdl_irb140()

    deg = pi/180;

    % robot length values (metres)
    d1 = 0.352;
    a1 = 0.070;
    a2 = 0.360;
    d4 = 0.380;
    d6 = 0.065;

    % DH parameter table
    %     theta d a alpha
    dh = [0 d1 a1  -pi/2
          0 0  a2  0
          0 0  0   pi/2
          0 d4 0   -pi/2
          0 0  0   pi/2
          0 d6 0   pi/2];

    % and build a serial link manipulator

    robot = SerialLink(dh, 'name', 'IRB 140', ...
        'manufacturer', 'ABB', 'ikine', 'nooffset');

    % place the variables into the global workspace
    if nargin == 1
        r = robot;
    elseif nargin == 0
        assignin('base', 'irb140', robot);
        assignin('base', 'qz', [0 0 0 0 0 0]); % zero angles
        assignin('base', 'qd', [0 -90 180 0 0 -90]*deg); % data sheet pose, horizontal
        assignin('base', 'qr', [0 -90 90 0 90 -90]*deg); % ready pose, arm up
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function k = diagKernDiagCompute(kern, x)

% DIAGKERNDIAGCOMPUTE Compute diagonal of DIAG kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the diagonal noise covariance function kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : diagKernParamInit, kernDiagCompute, kernCreate, diagKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

  if size(x, 2)>1
    error('Diag kernel requires 1-dimensional input.')
  end
  trans = str2func([kern.trans, 'Transform']);
  k = kern.variance*trans(x, 'atox');
end


text:
function tests = test_spm_dcm_peb_bmc_fam
% Unit Tests for test_spm_dcm_peb_bmc
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_peb_bmc_fam.m 6946 2016-11-23 15:26:29Z peter $

tests = functiontests(localfunctions);

% -------------------------------------------------------------------------
function test_peb_bmc(testCase)

data_path = get_data_path();

% Load
GCM = load(fullfile(data_path,'models','GCM_simulated.mat'));
PEB = load(fullfile(data_path,'PEB_test.mat'));
GCM = GCM.GCM;
PEB = PEB.PEB;

% Model comparison
[BMA,BMR] = spm_dcm_peb_bmc(PEB(1),GCM(1,:));

% Family 1: model 1 (full) and model 2 (fwd modulation)
% Family 2: model 3 (bwd) and model 4 (bwd modulation)
families = [1 1 2 2];

% Model comparison under families
[BMA,fam] = spm_dcm_peb_bmc_fam(BMA,BMR,families);

% -------------------------------------------------------------------------
% Should be equal priors as families were balanced
Nm = size(GCM,2);
testCase.assertEqual([0.5 0.5], fam.family.prior)
testCase.assertTrue(all(fam.model.prior(:) == 1/(Nm*Nm)));

% Family 1 and model 2 should win
testCase.assertTrue(fam.family.post(1,1) > 0.9);
testCase.assertTrue(fam.model.Pw(2) > 0.9);
testCase.assertTrue(fam.model.Px(2) > 0.9);

% -------------------------------------------------------------------------
% Run again with unbalanced families
families = [1 1 1 2];
[BMA,fam] = spm_dcm_peb_bmc_fam(BMA,BMR,families);

% Families should have equal priors
testCase.assertEqual([0.5 0.5], fam.family.prior);

% Model joint prior matrix should sum to 1
testCase.assertEqual(1, sum(fam.model.prior(:)));

% Models within each family should have same total prior probability
family_total_prior = [];
for f = 1:2
    p = (fam.model.prior(families == f,families == f));
    family_total_prior(f) = sum(p(:));
end
testCase.assertTrue(...
    all(family_total_prior(:) - family_total_prior(1) < 1e-10));

% -------------------------------------------------------------------------
function data_path = get_data_path()

data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'fMRI', 'simulated_2region');

author: spm


text:
function [net, info] = cnn_plate()
global datadir;
startup;
if exist(opts.imdbPath,'file')
    imdb=load(opts.imdbPath);
else
    imdb=cnn_plate_setup_data(datadir);
    mkdir(opts.expDir) ;
    save(opts.imdbPath, '-struct', 'imdb') ;
end
net=cnn_plate_init();
net.meta.normalization.averageImage =imdb.images.data_mean ;
opts.train.gpus=1;
[net, info] = cnn_train(net, imdb, getBatch(opts), ...
  'expDir', opts.expDir, ...
  net.meta.trainOpts, ...
  opts.train, ...
  'val', find(imdb.images.set == 3)) ;

function fn = getBatch(opts)
% --------------------------------------------------------------------
    fn = @(x,y) getSimpleNNBatch(x,y) ;
end
function [images, labels]  = getSimpleNNBatch(imdb, batch)
    images = imdb.images.data(:,:,:,batch) ;
    labels = imdb.images.labels(1,batch) ;
    if opts.train.gpus > 0
        images = gpuArray(images) ;
    end
end
end

author: imistyrain

repo: MatConvNet-oneclick


text:
function var = getvariables(F)
var = getvariables(lmi(F));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function domain = CreateDomain(Lx, Ly, dx, dy)

% initialize dependant variables
x_intervals = Lx/dx + 1;
y_intervals = Ly/dy + 1;
%x_Step_Interval = (stepX/dx) + 1;
%y_Step_Interval = (stepY/dy) + 1;

% Create Domain Matrix
domain = zeros(x_intervals,y_intervals);

end

author: Sable

repo: mcbench-benchmarks


text:
function SD = SD_evaluation(F)
[m,n]=size(F);
u=mean(mean(F));
SD=sqrt(sum(sum((F-u).^2))/(m*n));
end

author: Linfeng-Tang

repo: Image-Fusion


text:
% demo for linear decomposition

clear all
close all

k = 5; % number of simulated components
n = 15; % 1st dim of data matrix
p = 7; % 2nd dim of data matrix

A0 = randn(n,k);
B0 = randn(k,p);
Y = A0*B0;
y = Y+randn(size(Y));

[A,B,out,posterior] = VBA_LinDecomp(y,k);

hf = figure('color',[1 1 1]);
ha = subplot(2,2,1,'parent',hf);
plot(ha,A,A0,'.')
ha = subplot(2,2,2,'parent',hf);
plot(ha,B',B0','.')
ha = subplot(2,2,3,'parent',hf);
bar(corr(A,A0),'parent',ha)
ha = subplot(2,2,4,'parent',hf);
bar(corr(B',B0'),'parent',ha)

author: MBB-team

repo: VBA-toolbox


text:
function bool = gsp_check_fourier(G)
%GSP_CHECK_FOURIER Check is the Fourier basis is computed
%   Usage:  bool = gsp_check_fourier(G):
%
%   Input parameters:
%       G           : Graph structure
%   Output parameters:
%       bool        : boolean
%
%   This function check is the Laplacian eigenvalues and eigenfunctions are
%   computed.
%

bool = ( isfield(G,'e') && isfield(G,'U') );

end

repo: gspbox


text:
function [SPM] = spm_fMRI_design(SPM,save_SPM)
% Assembles a design for fMRI studies
% FORMAT [SPM] = spm_fMRI_design(SPM)
%
% 1st level
%--------------------------------------------------------------------------
% SPM.
%       xY: [1x1 struct] - data structure
%    nscan: [1xs double] - nscan(s) = number of scans in session s
%      xBF: [1x1 struct] - Basis function structure
%     Sess: [1xs struct] - Session structure array
%       xX: [1x1 struct] - Design matrix structure
%
%
%    2nd level
%    ----------------------------------------------------------------------
%    SPM.xY
%           RT: - repetition time {seconds)
%
%    SPM.xBF
%            T: - microtime resolution (number of time bins per scan)
%           T0: - microtime onset (reference time bin, see slice timing)
%        UNITS: - 'scans'|'secs' - units in which onsets are specified
%     Volterra: - 1|2 - order of [Volterra] convolution
%           dt: - length of time bin {seconds}
%         name: - name of basis set
%       length: - support of basis set {seconds}
%        order: - order of basis set
%           bf: - basis set matrix
%
%    SPM.Sess(s)
%            U: - Input structure array
%            C: - User specified covariate structure
%          row: - scan   indices for session s
%          col: - effect indices for session s
%           Fc: - F Contrast information for input-specific effects
%
%    SPM.xX
%            X: - design matrix
%           iH: - vector of H partition (indicator variables) indices
%           iC: - vector of C partition (covariates)          indices
%           iB: - vector of B partition (block effects)       indices
%           iG: - vector of G partition (nuisance variables)  indices
%         name: - cellstr of names for design matrix columns
%
%
%        3rd level
%        ------------------------------------------------------------------
%        SPM.Sess(s).U
%               dt: - time bin length {seconds}
%             name: - {1 x j} cell of names for each input or cause
%              ons: - (q x 1) onsets for q  trials {in UNITS}
%              dur: - (q x 1) durations for trials {in UNITS}
%                P: - Parameter stucture
%                u: - (t x j) inputs or stimulus function matrix
%              pst: - (1 x k) peristimulus times (seconds)
%             orth: - boolean: orthogonalise inputs?
%
%
%        SPM.Sess(s).C
%
%                C: - [kx1 double] of user specified regressors
%             name: - {1xk} cellstr of regressor names
%
%
%        SPM.Sess(s).Fc
%
%                i: - indices pertaining to each input
%             name: - names pertaining to each input
%                p: - grouping of regressors per parameter
%
%
%            4th level
%            --------------------------------------------------------------
%            SPM.Sess(s).U(i).P(p)
%
%                 name: - parameter name
%                    P: - (q x 1) parameter matrix
%                    h: - order of polynomial expansion (0 = none)
%                    i: - sub-indices of U(i).u for plotting
%
%
% saves SPM.mat if save_SPM==1 (this is the default)
%__________________________________________________________________________
%
% spm_fMRI_design allows you to build design matrices with separable
% session-specific partitions.  Each partition may be the same (in which
% case it is only necessary to specify it once) or different.  Responses
% can be either event- or epoch related, where the latter model prolonged
% and possibly time-varying responses to state-related changes in
% experimental conditions.  Event-related response are modelled in terms
% of responses to instantaneous events.  Mathematically they are both
% modelled by convolving a series of delta (stick) or box-car functions,
% encoding the input or stimulus function. with a set of hemodynamic
% basis functions.
%
% spm_fMRI_design allows you to combine both event- and epoch-related
% responses in the same model and/or regressor. You specify the number
% of trial (event or epoch) types.  Epoch and event-related
% responses are modeled in exactly the same way by first specifying their
% onsets [in terms of onset times] and then their durations.  Events are
% specified with a duration of 0.  If you enter a single number for the
% durations it will be assumed that all trials conform to this duration.
%
% Interactions or response modulations can enter at two levels.  Firstly
% the stick function itself can be modulated by some parametric variate
% (this can be time or some trial-specific variate like reaction time)
% modeling the interaction between the trial and the variate or, secondly
% interactions among the trials themselves can be modeled using a Volterra
% series formulation that accommodates interactions over time (and therefore
% within and between trial types).  The first sort of interaction is
% specified by extra (modulated) stick functions in Sess(s).u.  If
% a polynomial expansion of the specified variate is requested there will
% be more than one column.  The corresponding name of the explanatory
% variables in X.name is Sn(s) trial(u)xparam(p)^q*bf(i) for the qth
% order expansion of the parameter convolved with the ith basis function
% for the uth trial in the sth session.  If no parametric variate is
% specified the name is simply Sn(s) trial(u)*bf(i).  Interactions among
% and within trials enter as new trial types but do not have .pst or .ons
% fields.  These interactions can be characterized later, in results, in
% terms of the corresponding second order Volterra Kernels.
%
% The design matrix is assembled on a much finer time scale (xBF.dt) than the
% TR and is then sub-sampled at the acquisition times.  After down-sampling
% the regressors for each input are othogonalised.  This ensures that
% components due to the canonical hrf are not explained away by other basis
% functions or parametric modulators.
%
% Sess(s).ons(u) contains onset times in seconds or scans relative to the
% timing of the first scan
%
% Notes on spm_get_ons, spm_get_bf and spm_Volterra are included below
% for convenience.
%
%                           ----------------
%
% spm_get_ons contructs a struct array containing sparse input
% functions U(i).u specifying occurrence events or epochs (or both).
% These are convolved with a basis set at a later stage to give
% regressors that enter into the design matrix. Interactions of evoked
% responses with some parameter (time or a specified  variate P) enter at
% this stage as additional columns in U(u).u with each trial multiplied
% by the [expansion of the] trial-specific parameter. If parametric
% modulation is modeled, P(p).P contains the original variate and
% P(p).name is its name. The 0th order expansion of this is simply the main
% effect in the first column of U(u).u
%
%                           ----------------
%
% spm_get_bf prompts for basis functions to model hemodynamic
% responses.  The basis functions returned are orthogonalized
% and defined as a function of peri-stimulus time in time-bins.
%
%                           ----------------
%
% For first order expansions spm_Volterra simply convolves the causes
% (e.g. stick functions) in U(u).u by the basis functions in Sess(s).bf
% to create design matrix X.  For second order expansions new entries appear
% in the design matrix that correspond to the hemodynamic interaction among the
% orginal causes (if the events are sufficiently close in time).
% The basis functions for these are two dimensional and are used to
% assemble the second order kernel in spm_graph.m.  Second order effects
% are computed for only the first column of U(u).u.
%
%__________________________________________________________________________
% Copyright (C) 1999-2012 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fMRI_design.m 7739 2019-12-02 14:00:18Z guillaume $

SVNid = '$Rev: 7739 $';

%-Say Hello
%--------------------------------------------------------------------------
spm('FnBanner',mfilename,SVNid);

%-Construct Design matrix {X}
%==========================================================================

%-Microtime onset and microtime resolution
%--------------------------------------------------------------------------
try
    fMRI_T     = SPM.xBF.T;
    fMRI_T0    = SPM.xBF.T0;
catch
    fMRI_T     = spm_get_defaults('stats.fmri.t');
    fMRI_T0    = spm_get_defaults('stats.fmri.t0');
    SPM.xBF.T  = fMRI_T;
    SPM.xBF.T0 = fMRI_T0;
end

%-Time units, dt = time bin {secs}
%--------------------------------------------------------------------------
SPM.xBF.dt     = SPM.xY.RT/SPM.xBF.T;

%-Get basis functions
%--------------------------------------------------------------------------
SPM.xBF        = spm_get_bf(SPM.xBF);

%-Get session specific design parameters
%==========================================================================
Xx    = [];
Xb    = [];
Xname = {};
Bname = {};
for s = 1:length(SPM.nscan)

    %-Number of scans for this session
    %----------------------------------------------------------------------
    k = SPM.nscan(s);

    %-Create convolved stimulus functions or inputs
    %======================================================================

    %-Get inputs, neuronal causes or stimulus functions U
    %----------------------------------------------------------------------
    U = spm_get_ons(SPM,s);

    %-Convolve stimulus functions with basis functions
    %----------------------------------------------------------------------
    [X,Xn,Fc] = spm_Volterra(U, SPM.xBF.bf, SPM.xBF.Volterra);

    %-Resample regressors at acquisition times (32 bin offset)
    %----------------------------------------------------------------------
    if ~isempty(X)
        X = X((0:(k - 1))*fMRI_T + fMRI_T0 + 32,:);
    end

    %-Orthogonalise (within trial type)
    %----------------------------------------------------------------------
    for i = 1:length(Fc)
        if i<= numel(U) && ... % for Volterra kernels
                (~isfield(U(i),'orth') || U(i).orth)
            p = ones(size(Fc(i).i));
        else
            p = Fc(i).p;
        end
        for j = 1:max(p)
            X(:,Fc(i).i(p==j)) = spm_orth(X(:,Fc(i).i(p==j)));
        end
    end

    %-Get user specified regressors
    %======================================================================
    C     = SPM.Sess(s).C.C;
    Cname = SPM.Sess(s).C.name;

    %-Append mean-corrected regressors and names
    %----------------------------------------------------------------------
    X     = [X spm_detrend(C)];
    Xn    = {Xn{:} Cname{:}};

    %-Confounds: Session effects
    %======================================================================
    B     = ones(k,1);
    Bn    = {'constant'};

    %-Session structure array
    %======================================================================
    SPM.Sess(s).U      = U;
    SPM.Sess(s).C.C    = C;
    SPM.Sess(s).C.name = Cname;
    SPM.Sess(s).row    = size(Xb,1) + (1:k);
    SPM.Sess(s).col    = size(Xx,2) + (1:size(X,2));
    SPM.Sess(s).Fc     = Fc;

    %-Append into Xx and Xb
    %======================================================================
    Xx      = blkdiag(Xx,X);
    Xb      = blkdiag(Xb,B);

    %-Append names
    %----------------------------------------------------------------------
    for i = 1:length(Xn)
        Xname{end + 1} = [sprintf('Sn(%i) ',s) Xn{i}];
    end
    for i = 1:length(Bn)
        Bname{end + 1} = [sprintf('Sn(%i) ',s) Bn{i}];
    end

end

%-Place design matrix structure in xX
%==========================================================================
SPM.xX.X    = [Xx Xb];
SPM.xX.iH   = [];
SPM.xX.iC   = 1:size(Xx,2);
SPM.xX.iB   = (1:size(Xb,2)) + size(Xx,2);
SPM.xX.iG   = [];
SPM.xX.name = {Xname{:} Bname{:}};

if nargin < 2 || save_SPM
    %-Save SPM.mat
    %----------------------------------------------------------------------
    fprintf('%-40s: ','Saving fMRI design')                             %-#
    fmt = spm_get_defaults('mat.format');
    s = whos('SPM');
    if s.bytes > 2147483647, fmt = '-v7.3'; end
    save('SPM.mat', 'SPM', fmt);
    fprintf('%30s\n','...SPM.mat saved');
end


text:
function otr = transform(itr,tmx)

%
% FUNCTION
% otr = transform(itr,tmx)
%
% DESCRIPTION
% transform 3D-tensor (Euclidean or Cartesion tensor) of any order (>0) to another coordinate system
%
% PARAMETERS
% otr = output tensor, after transformation; has the same dimensions as the input tensor
% itr = input tensor, before transformation; should be a 3-element vector, a 3x3 matrix, or a 3x3x3x... multidimensional array, each dimension containing 3 elements
% tmx = transformation matrix, 3x3 matrix that contains the direction cosines between the old and the new coordinate system
%

ne = numel(itr);                % number of tensor elements
nd = ndims(itr);                % number of tensor dimensions, i.e. order of tensor
if (ne==3), nd = 1; end         % order of tensor is 1 in case of a 3x1 or 1x3 vector

otr = itr;                      % create output tensor
otr(:) = 0;                     % fill output tensor with zeros; this way a symbolic tensor remains symbolic

iie = zeros(nd,1);              % initialise vector with indices of input tensor element
ioe = zeros(nd,1);              % initialise vector with indices of output tensor element
cne = cumprod(3*ones(nd,1))/3;  % vector with cumulative number of elements for each dimension (divided by three)

for oe = 1:ne,                                 % loop over all output elements
   ioe = mod(floor((oe-1)./cne),3)+1;          % calculate indices of current output tensor element
   for ie = 1:ne,                              % loop over all input elements
      pmx = 1;                                 % initialise product of transformation matrices
      iie = mod(floor((ie-1)./cne),3)+1;       % calculate indices of current input tensor element
      for id = 1:nd,                           % loop over all dimensions
         pmx = pmx * tmx( ioe(id), iie(id) );  % create product of transformation matrices
      end
      otr(oe) = otr(oe) + pmx * itr(ie);       % add product of transformation matrices and input tensor element to output tensor element
   end
end


text:
%% Example: Random Category Registration problem and its semidefinite relaxations
%% Heng Yang, July 06, 2021

clc; clear; close all; restoredefaultpath

%% paths to dependencies
spotpath    = '../spotless';
mosekpath   = '../../mosek';
stridepath  = '../STRIDE';
manoptpath  = '../manopt';
sdpnalpath  = '../../SDPNAL+v1.0';
path.stridepath = stridepath;
path.mosekpath  = mosekpath;
path.manoptpath = manoptpath;
addpath('../utils')
addpath('./solvers')

%% choose if run GNC for STRIDE
rungnc      = true;

%% generate random point cloud registration problem
datatype                 = 'car';
problem.outlierRatio     = 0.1;
problem.noiseSigma       = 0.01;
problem.translationBound = 10.0;
switch datatype
    case 'random'
        %%%%%%%%%%%%%% RANDOM data %%%%%%%%%%%%%%
        problem.N                = 10;
        problem.K                = 3;
        problem                  = gen_category_registration(problem);
    case 'car'
        %%%%%%%%%%%%%% PASCAL Car %%%%%%%%%%%%%%
        problem.path             = './data/car.mat';
        problem                  = gen_catreg_pascal_car(problem);
end
lambda                   = 0.5;

%% generate SDP relaxations
addpath(genpath(spotpath))
SDP     = relax_category_registration_v2(problem,'lambda',lambda,'checkMonomials',false);
rmpath(genpath(spotpath))

%% Solve using STRIDE
% Primal initialization
if rungnc
    solution        = gnc_category_registration(problem,SDP,path,'lambda',lambda);
    [R_est,t_est]   = invert_transformation(solution.R_est,solution.t_est);

    v0       = lift_catreg_v2(R_est(:),...
                              t_est(:),...
                              solution.c_est(:),...
                              solution.theta_est(:),...
                              problem.cBound,...
                              problem.translationBound);
    X0          = rank_one_lift(v0);
    gnc.R_err   = getAngularError(problem.R_gt,solution.R_est);
    gnc.t_err   = getTranslationError(problem.t_gt,solution.t_est);
    gnc.c_err   = getTranslationError(problem.c_gt,solution.c_est);
    gnc.time    = solution.time;
    gnc.info    = solution;

else
    X0       = [];
end

% dual initialization
addpath(genpath(spotpath))
chordalSDP     = chordal_relax_category_registration_v2(problem,'lambda',lambda);
rmpath(genpath(spotpath))
prob = convert_sedumi2mosek(chordalSDP.sedumi.At,...
                            chordalSDP.sedumi.b,...
                            chordalSDP.sedumi.c,...
                            chordalSDP.sedumi.K);
addpath(genpath(mosekpath))
time0   = tic;
param.MSK_IPAR_INTPNT_MAX_ITERATIONS = 20; % set maximum iterations 20
[~,res] = mosekopt('minimize info',prob,param);
time_dualInit = toc(time0);
[Xchordal,ychordal,Schordal,~] = recover_mosek_sol_blk(res,chordalSDP.blk);
S_assm             = catreg_dual_from_chordal_dual_v2(Schordal,problem.N,problem.K);

% STRIDE main algorithm
addpath(genpath(stridepath))
addpath(genpath(manoptpath))

pgdopts.pgdStepSize     = 10;
pgdopts.SDPNALpath      = sdpnalpath;
pgdopts.maxiterPGD      = 5;
% ADMM parameters
pgdopts.tolADMM         = 1e-10;
pgdopts.maxiterADMM     = 1e4;
pgdopts.stopoptionADMM  = 0;

pgdopts.rrOpt           = 1:3;
pgdopts.rrFunName       = 'local_search_catreg_v2';
rrPar.blk = SDP.blk; rrPar.translationBound = problem.translationBound;
rrPar.N = problem.N; rrPar.K = problem.K; rrPar.cBound = problem.cBound;
pgdopts.rrPar           = rrPar;
pgdopts.maxiterLBFGS    = 1000;
pgdopts.maxiterSGS      = 1000;
pgdopts.lbfgseps        = false;
pgdopts.S0              = S_assm;

[outPGD,Xopt,yopt,Sopt]     = PGDSDP(SDP.blk, SDP.At, SDP.b, SDP.C, X0, pgdopts);
rmpath(genpath(manoptpath))

infostride                  = get_performance_catreg_v2(Xopt,yopt,Sopt,SDP,problem,stridepath);
infostride.totaltime        = outPGD.totaltime + time_dualInit;
infostride.time             = [outPGD.totaltime,time_dualInit];
if rungnc
    infostride.gnc  = gnc;
    infostride.totaltime = infostride.totaltime + gnc.time;
    infostride.time      = [gnc.time,infostride.time];
end

fprintf('\n\n\n\n\n')

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

% Signal space
N = 1000;
% Number of measurements
M = 200;
% Sparsity levels
Ks = 4:120;

% Number of dictionaries to be created
num_dict_trials = 100;
% Number of signals to be created for each dictionary
num_signal_trials = 20;

% Number of trials for each K
num_trials = num_dict_trials * num_signal_trials;

omp_success_rates_with_k = zeros(numel(Ks), 1);
omp_average_iterations_with_k = zeros(numel(Ks), 1);
omp_maximum_iterations_with_k = zeros(numel(Ks), 1);

Ls = [2, 4, 6, 8]
num_ls = numel(Ls)
gomp_success_rates_with_k = zeros(numel(Ks), num_ls);
gomp_average_iterations_with_k = zeros(numel(Ks), num_ls);
gomp_maximum_iterations_with_k = zeros(numel(Ks), num_ls);

for K=Ks
    % Trial number
    nt = 0;
    omp_num_successes = 0;
    omp_num_iterations = 0;
    omp_max_iterations = 0;

    gomp_num_successes = zeros(num_ls, 1);
    gomp_num_iterations = zeros(num_ls, 1);
    gomp_max_iterations = zeros(num_ls, 1);

    for ndt=1:num_dict_trials
        % Sensing matrix
        Phi = spx.dict.simple.gaussian_dict(M, N);
        for nst=1:num_signal_trials
            nt = nt + 1;
            % Construct the signal generator.
            gen  = spx.data.synthetic.SparseSignalGenerator(N, K);
            % Generate bi-uniform signals
            x = gen.gaussian;
            % Measurement vectors
            y = Phi.apply(x);

            % OMP solver instance
            solver = spx.pursuit.single.OrthogonalMatchingPursuit(Phi, K);
            % Solve the sparse recovery problem
            omp_result = solver.solve(y);
            % Solution vector
            z = omp_result.z;
            omp_stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
            omp_num_iterations = omp_num_iterations + omp_result.iterations;
            if omp_max_iterations < omp_result.iterations
                omp_max_iterations = omp_result.iterations;
            end
            omp_num_successes = omp_num_successes + omp_stats.success;
            fprintf('K=%d, Trial: %d, OMP: %s, ', ...
                K, nt, spx.io.true_false_short(omp_stats.success));
            for nl=1:num_ls
                L = Ls(nl);
                % GOMP solver instance
                solver = spx.pursuit.single.GOMP(Phi, K);
                % Set the number of atoms to be selected in each iteration
                solver.L = L;
                % Solve the sparse recovery problem
                gomp_result = solver.solve(y);
                % Solution vector
                z = gomp_result.z;
                gomp_stats = spx.commons.sparse.recovery_performance(Phi, K, y, x, z);
                gomp_num_iterations(nl) = gomp_num_iterations(nl) + gomp_result.iterations;
                if gomp_max_iterations(nl) < gomp_result.iterations
                    gomp_max_iterations(nl) = gomp_result.iterations;
                end
                gomp_num_successes(nl) = gomp_num_successes(nl) + gomp_stats.success;
                fprintf(' GOMP-%d:%s', ...
                    L, spx.io.true_false_short(gomp_stats.success));
            end
            fprintf('\n')
        end
    end
    omp_success_rate = omp_num_successes / num_trials;
    omp_average_iterations = omp_num_iterations / num_trials;
    omp_success_rates_with_k(K) = omp_success_rate;
    omp_average_iterations_with_k (K) = omp_average_iterations;
    omp_maximum_iterations_with_k(K) = omp_max_iterations;

    for nl=1:num_ls
        gomp_success_rate = gomp_num_successes(nl) / num_trials;
        gomp_average_iterations = gomp_num_iterations(nl) / num_trials;
        gomp_success_rates_with_k(K, nl) = gomp_success_rate;
        gomp_average_iterations_with_k (K, nl) = gomp_average_iterations;
        gomp_maximum_iterations_with_k(K, nl) = gomp_max_iterations(nl);
    end
end

save('bin/omp_vs_gomp_comparison.mat');

author: indigits

repo: sparse-plex


text:
%% abaqusStruct2inp
% Below is a demonstration of the features of the |abaqusStruct2inp| function

%%
clear; close all; clc;

%% Syntax
% |[T]=abaqusStruct2inp(abaqus_spec,fileName,optionStruct);|

%% Description
% This function provides the basis for coding Abaqus .inp input files

%%
% See also: |abaqusStructTemplate|

%% About Abaqus INP input files
% Abaqus input files (.inp) are plain text files.

%% Coding INP files in MATLAB

%%
% Specify file heading
% Users may provide a heading which can be used to give a description of
% the input file.

% *Heading
% ** Job name: ABAQUS inp file creation demo
% ** Generated by: GIBBON

abaqus_spec.Heading.COMMENT{1}='Job name: ABAQUS inp file creation demo';
abaqus_spec.Heading.COMMENT{2}='Generated by: GIBBON';

%%
% The preprint section
%

% *Preprint, echo=NO, model=NO, history=NO, contact=NO
abaqus_spec.Preprint.ATTR.echo='NO';
abaqus_spec.Preprint.ATTR.model='NO';
abaqus_spec.Preprint.ATTR.history='NO';
abaqus_spec.Preprint.ATTR.contact='NO';

%%
% The part section

V=[ 1 1 1; 1 0 1;  1 1 0; 1 0 0; 0 1 1; 0 0 1; 0 1 0; 0 0 0];

nodeIds=(1:1:size(V,1));
E=[5 6 8 7 1 2 4 3];
elementIds=(1:1:size(E,1));

% Node
abaqus_spec.Part{1}.COMMENT='This section defines the part geometry in terms of nodes and elements';
abaqus_spec.Part{1}.ATTR.name='Cube';
abaqus_spec.Part{1}.Node={nodeIds(:),V};

% Element
abaqus_spec.Part{1}.Element{1}.ATTR.type='C3D8R';
abaqus_spec.Part{1}.Element{1}.VAL={elementIds(:),E};

%Element sets
abaqus_spec.Part{1}.Elset{1}.ATTR.elset='Set-1';
abaqus_spec.Part{1}.Elset{1}.VAL=1;

%Node sets
abaqus_spec.Part{1}.Nset{2}.ATTR.nset='Set-2';
abaqus_spec.Part{1}.Nset{2}.VAL=1:1:size(V,1);

%Node sets using "generate"
abaqus_spec.Part{1}.Nset{1}.ATTR.nset='Set-1';
abaqus_spec.Part{1}.Nset{1}.ATTR.generate='';
abaqus_spec.Part{1}.Nset{1}.VAL=[1 size(V,1) 1];

%Sections
abaqus_spec.Part{1}.Solid_section.ATTR.elset='Set-1';
abaqus_spec.Part{1}.Solid_section.ATTR.material='Elastic';

% abaqus_spec.Part{1}.Shell_section.ATTR.elset='Set-1';
% abaqus_spec.Part{1}.Shell_section.ATTR.material='Material-1';
% % abaqus_spec.Part{1}.Shell_section.VAL=[10 5]; %Use numerical data
% abaqus_spec.Part{1}.Shell_section.VAL='10., 5'; %Use custom text data

%%
% The assembly section

abaqus_spec.Assembly.ATTR.name='Assembly-1';

abaqus_spec.Assembly.Instance{1}.ATTR.name='Cube-1';
abaqus_spec.Assembly.Instance{1}.ATTR.part='Cube';

% abaqus_spec.Assembly.Instance{2}.ATTR.name='Clot-1';
% abaqus_spec.Assembly.Instance{2}.ATTR.part='Clot';

abaqus_spec.Assembly.Nset{1}.ATTR.nset='Set-1';
abaqus_spec.Assembly.Nset{1}.ATTR.instance='Cube-1';
abaqus_spec.Assembly.Nset{1}.ATTR.generate='';
abaqus_spec.Assembly.Nset{1}.VAL=[1 4 1];

abaqus_spec.Assembly.Nset{2}.ATTR.nset='Set-3';
abaqus_spec.Assembly.Nset{2}.ATTR.instance='Cube-1';
abaqus_spec.Assembly.Nset{2}.ATTR.generate='';
abaqus_spec.Assembly.Nset{2}.VAL=[5 7 1];

abaqus_spec.Assembly.Nset{3}.ATTR.nset='Set-4';
abaqus_spec.Assembly.Nset{3}.ATTR.instance='Cube-1';
abaqus_spec.Assembly.Nset{3}.ATTR.generate='';
abaqus_spec.Assembly.Nset{3}.VAL=[2 6 2];

abaqus_spec.Assembly.Nset{4}.ATTR.nset='Set-5';
abaqus_spec.Assembly.Nset{4}.ATTR.instance='Cube-1';
abaqus_spec.Assembly.Nset{4}.VAL=[3 4 7];

abaqus_spec.Assembly.Nset{5}.ATTR.nset='origin';
abaqus_spec.Assembly.Nset{5}.ATTR.instance='Cube-1';
abaqus_spec.Assembly.Nset{5}.VAL=8;

abaqus_spec.Assembly.Nset{6}.ATTR.nset='all';
abaqus_spec.Assembly.Nset{6}.ATTR.instance='Cube-1';
% abaqus_spec.Assembly.Nset{6}.ATTR.generate='';
abaqus_spec.Assembly.Nset{6}.VAL=1:size(V,1);

abaqus_spec.Assembly.Surface{1}.ATTR.type='ELEMENT';
abaqus_spec.Assembly.Surface{1}.ATTR.name='innerSurface';
elemSetsInner={[1 2 3 4],[4 5 6 7]};
elemSetSideInner=[1 2];
valCell=cell(numel(elemSetsInner),1);
for qs=1:1:numel(elemSetsInner)
    surfaceElementSetName=['elementSetInnerSurface',num2str(qs)];
    abaqus_spec.Assembly.Elset{qs}.ATTR.elset=surfaceElementSetName;
    abaqus_spec.Assembly.Elset{qs}.ATTR.internal=''; %Remains hidden uppon import
    abaqus_spec.Assembly.Elset{qs}.ATTR.instance='Cube-1';
    abaqus_spec.Assembly.Elset{qs}.VAL=elemSetsInner{qs};

    sidePick=elemSetSideInner(qs);
    valCell{qs}={[surfaceElementSetName,', S',num2str(sidePick)]};
end
abaqus_spec.Assembly.Surface{1}.VAL=valCell;

%%
% Equation constraint

% *Equation
% 2
% X(n-1), 1, 1.
% XControlN, 1, -1.

%%--> Constraints
abaqus_spec.Assembly.COMMENT='Constraint 1';
abaqus_spec.Assembly.Equation.VAL=['2' newline 'X(n-1), ',  '1, 1.' newline 'XContN, ' '1, -1.'];

%%
% The material section

% *Material, name=Material-1
% *Hyperelastic, neo hooke
%  0.03,1.

abaqus_spec.Material.ATTR.name='Neo-Hooke';
abaqus_spec.Material.Hyperelastic.VAL=[0.03 1];
abaqus_spec.Material.Hyperelastic.ATTR.neo_hooke='';
abaqus_spec.Material.ATTR.name='Elastic';
abaqus_spec.Material.Elastic=[0.5 0.49];

%%
% Orientation

% *Orientation, name=Ori-1, system=RECTANGULAR
% Ori-1-DiscOrient
% 2, 90.
abaqus_spec.Orientation.ATTR.name='Ori-1';
abaqus_spec.Orientation.ATTR.system='RECTANGULAR';
abaqus_spec.Orientation.VAL=['Ori-1-DiscOrient' newline vec2strIntDouble([2 90],'%i')];

%%
% Contact pair

%%--> Contact pair
abaqus_spec.contact_pair{1}.ATTR.interaction='user-1';
abaqus_spec.contact_pair{1}.ATTR.type='node to surface';
abaqus_spec.contact_pair{1}.ATTR.adjust='0.0';
abaqus_spec.contact_pair{1}.VAL='AORTA-ASSEMBLY.OUTER, AORTA1-ASSEMBLY.OUTERX';

% surface interaction
abaqus_spec.surface_interaction{1}.ATTR.name='user-1';
abaqus_spec.surface_interaction{1}.ATTR.user='';
abaqus_spec.surface_interaction{1}.ATTR.properties='7';
abaqus_spec.surface_interaction{1}.ATTR.unsymm='';
abaqus_spec.surface_interaction{1}.ATTR.DEPVAR='13';
abaqus_spec.surface_interaction{1}.CSTM=', '; %Line with just comma

%Wrap to max width of 8 entries
t=vec2strIntDouble([0.2, 0.2, 0.05, 0.05,150,1000,0.0],'%6.7e');
t=strwrap(t,8,', ');
abaqus_spec.surface_interaction{1}.VAL=t;

%%
% The step section
%

abaqus_spec.Step.ATTR.name='Step-1';
abaqus_spec.Step.ATTR.nlgeom='YES';
abaqus_spec.Step.Static=[0.1 1 1e-5 0.1];

% Boundary
abaqus_spec.Step.Boundary{1}.VAL={'Set-1',[1 1],0.1};
abaqus_spec.Step.Boundary{2}.VAL={repmat({'origin'},3,1) [1 1; 2 2; 3 3]};
abaqus_spec.Step.Boundary{3}.VAL={'Set-3',[1,1]};
abaqus_spec.Step.Boundary{4}.VAL={'Set-4',[2,2]};
abaqus_spec.Step.Boundary{5}.VAL={'Set-5',[3,3]};

% Loads
% ** LOADS
% **
% ** Name: Load-1   Type: Pressure
% *Dsload
% Surf-1, P, 0.00533
% **
abaqus_spec.Step.Dsload{1}.VAL={'Surf-1','P',0.00533};

%Output
abaqus_spec.Step.Restart.ATTR.write='';
abaqus_spec.Step.Restart.ATTR.frequency=0;

abaqus_spec.Step.Output{1}.ATTR.field='';
abaqus_spec.Step.Output{1}.ATTR.variable='PRESELECT';
abaqus_spec.Step.Output{2}.ATTR.history='';
abaqus_spec.Step.Output{2}.ATTR.variable='PRESELECT';
abaqus_spec.Step.Node_print.ATTR.nset='all';
abaqus_spec.Step.Node_print.ATTR.frequency = 1;
abaqus_spec.Step.Node_print.VAL='COORD';
abaqus_spec.Step.El_print.VAL='S';

% * NODE SET, set=yoursetname, frequency = 1
% U

%%

%% Creating the INP file
% You can use |abaqusStruct2inp| to write the structure data to a file.

%Create file name for INP file
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');
fileName=fullfile(savePath,'tempModel.inp');
[~,fileNamePart,~]=fileparts(fileName);

%%

[T]=abaqusStruct2inp(abaqus_spec,fileName);

%% View the inp file
%
textView(fileName);

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON