text:
function T = slgbfe(X, G, Gc, dy, fm, varargin)
%SLGBFE Performs Graph-based Feature Extraction Learning
%
% $ Syntax $
%   - T = slgbfe(X, G, Gc, dy, fm, ...)
%
% $ Arguments $
%   - X:        The sample matrix
%   - G:        The graph to be optimized
%   - Gc:       The constraint graph
%   - dy:       The dimension of feature space
%   - fm:       The formulation type
%   - T:        The learned transform matrix (dx x dy)
%               the transform is done by y = T' * x
%
% $ Description $
%   - T = slgbfe(X, G, Gc, dy, fm, ...) performs graph-based feature
%     extraction learning. It is to solve the following optimization.
%
%       min/max  T'X M(G) X'T,    s.t. T'X M(Gc) X'T = I
%
%     The concrete formulation depends on the formulation type given in
%     fm = {fg, fc}. For fg, it has the following three types:
%       - 'minW':   do minimization with M(G) = W
%       - 'maxW':   do maximization with M(G) = W
%       - 'minL':   do minimization with M(G) = L = D - W
%       - 'maxL':   do maximization with M(G) = L = D - W
%     For fc, it has the following three types:
%       - 'O':      constraint T be orthogonal: T'*T = I (ignore Gc)
%       - 'I':      constraint T'* X * X' * T = I (ignore Gc)
%       - 'WC':     constraint with M(Gc) = W of Gc
%       - 'LC':     constraint with M(Gc) = L of Gc: D - W
%     In the aforementioned formulation, W is the adjacency matrix, while
%     L is the Laplacian matrix. When Gc is ignored (as in 'O' and 'I'),
%     you can just input Gc as [].
%
%     You can further specify the following properties to control the
%     learning process:
%       - 'whparams':  The parameters for doing whitening of M(Gc), please
%                      refer to the function slwhiten_from_cov. The params
%                      are given in a cell array as {method, ...}.
%                      default = {}
%       - 'skip':      The number of components to be skipped. default = 0
%
% $ Remarks $
%   - The implementation is based on slgembed.
%
%   - The function will not centralize the samples, if it is needed please
%     centralize them before invoking.
%
% $ History $
%   - Created by Dahua Lin, on Sep 17, 2006
%

%% parse and verify input arguments

if nargin < 5
    raise_lackinput('slgbfe', 5);
end

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

if ~iscell(fm) || length(fm) ~= 2
    error('sltoolbox:invalidarg', 'fm should be a length-2 cell array');
end
fg = fm{1};
fc = fm{2};

gi = slgraphinfo(G, {[n, n]});
W = sladjmat(G, ...
    'valtype', 'numeric', ...
    'sparse', strcmp(gi.form, 'adjmat') && issparse(G));

if strcmp(fc, 'WC') || strcmp(fc, 'LC')
    if isempty(Gc)
        error('sltoolbox:invalidarg', ...
            'When fc is WC or LC, Gc should not be empty');
    end
    slgraphinfo(Gc, {'adjmat', [n, n]});
    if isnumeric(Gc)
        Wc = Gc;
    else
        Wc = double(Gc);
    end
else
    Wc = [];
end

opts.whparams = {};
opts.skip = 0;
opts = slparseprops(opts, varargin{:});

%% Construct problem

% enforce symmetry
W = (W + W') * (1/2);
if ~isempty(Wc)
    Wc = (Wc + Wc') * (1/2);
end

% construct re-formulated G: R
switch fg
    case 'maxW'
        R = X * W * X';
        rfg = 'maxW';
    case 'minW'
        R = X * W * X';
        rfg = 'minW';
    case 'maxL'
        R = X * make_Lmat(W) * X';
        rfg = 'maxW';
    case 'minL'
        R = X * make_Lmat(W) * X';
        rfg = 'minW';
    otherwise
        error('sltoolbox:invalidarg', 'Invalid fg name: %s', fg);
end

% construct re-formulated Gc: Rc
switch fc
    case 'O'
        Rc = [];
        rfc = 'I';
    case 'I'
        Rc = X * X';
        rfc = 'WC';
    case 'WC'
        Rc = X * Wc * X';
        rfc = 'WC';
    case 'LC'
        Rc = X * make_Lmat(Wc) * X';
        rfc = 'WC';
    otherwise
        error('sltoolbox:invalidarg', 'Invalid fc name: %s', fc);
end

%% solve problem

Y = slgembed(R, Rc, dy, {rfg, rfc}, ...
    'inv', opts.whparams, ...
    'skip', opts.skip);
T = Y';

%% Computational routines

function L = make_Lmat(W)

vD = sum(W, 1)';
if issparse(vD)
    vD = full(vD);
end

n = size(W, 1);
if issparse(W)
    D = sparse((1:n)', (1:n)', vD, n, n, n);
    L = D - W;
else
    L = -W;
    dinds = (1:n)'*(n+1)-n;
    L(dinds) = L(dinds) + vD;
end

author: lmthang


text:
clc;
close all;
clearvars;
rng('default');

n = 8;
A1 = gallery('binomial',n);
A2 = A1;

r = 4;
Q_factors = zeros(n, r);
for i=1:r
    [v, beta] = spx.la.house.gen(A2(i:n, i));
    Q_factors(i+1:n, i) = v(2:n-i+1);
    A2(i:n, i:n) = spx.la.house.premul(A2(i:n, i:n), v, beta)
    v_full = [zeros(i-1, 1); v];
end
% We compute the factors W and Y for Q = Q_1 * Q_2 * .. * Q_r
% such that Q = I - WY'
[W, Y] = spx.la.house.wy(Q_factors);
% We pre-multiply with I - WY'
A3 = spx.la.house.wy_premul(A2, W, Y);
A1
A3
max(max(abs(A1 - A3)))

author: indigits

repo: sparse-plex


text:
% Solution_Curve_ODE.
% Copyright Springer 2013. Stephen Lynch.
% See Figure 1.2.
deqn=@(t,c) .00713*(4-c)^2*(1-c/2);
[t,ca]=ode45(deqn,[0 400],0);
plot(t,ca(:,1));
axis([0 400 0 3]);
fsize=15;
set(gca,'xtick',0:100:400,'FontSize',fsize);
set(gca,'ytick',0:1:3,'FontSize',fsize);
xlabel('t','FontSize',fsize);
ylabel('c(t)','FontSize',fsize);

author: Sable

repo: mcbench-benchmarks


text:
function jv = toJava(v,collection)
%toJava       Convert to Java representation.
% toJava(v) returns a Java object representing the matlab value v.
% toJava(v,1) requires the result to be a collection, not a simple object
% such as java.lang.Double.
%
% Conversions:
% scalar           -> java.lang.Double
% numeric array    -> javaArray of Double
% character array  -> java.lang.String or array of same
% empty array      -> java.util.Vector
% cell array       -> javaArray of Objects
% empty cell array -> java.util.BitSet
% structure        -> java.util.HashTable
% class            -> struct with '_class' property
%
% See also fromJava.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

if nargin < 2
  collection = 0;
end
% in case v already is a Java object
jv = v;
if isnumeric(v)
  if isempty(v)
    jv = java.util.Vector;
  elseif length(v) == 1
    if collection
      jv = javaArray('java.lang.Double',1);
      jv(1) = java.lang.Double(v);
    else
      jv = java.lang.Double(v);
    end
  else
    sz = size(v);
    jv = javaArray('java.lang.Double',sz);
    v = v(:);
    for i = 1:prod(sz)
      index = substruct('()',num2cell(ind2subv(sz,i)));
      subsasgn(jv,index,toJava(v(i)));
    end
  end
  return
end
if ischar(v)
  n = rows(v);
  if n == 1 & ~collection
    jv = java.lang.String(v);
  else
    jv = javaArray('java.lang.String',n);
    for i = 1:n
      jv(i) = java.lang.String(deblank(v(i,:)));
    end
  end
end
if iscell(v)
  if isempty(v)
    jv = java.util.BitSet;
  else
    sz = size(v);
    jv = javaArray('java.lang.Object',sz);
    v = v(:);
    for i = 1:prod(sz)
      index = substruct('()',num2cell(ind2subv(sz,i)));
      subsasgn(jv,index,toJava(v{i}));
    end
  end
  return
end
if isobject(v) | isstruct(v)
  jv = java.util.Hashtable;
  if isobject(v)
    jv.put('_class',class(v));
    v = struct(v);
  end
  % record the order of the fieldnames
  fields = java.util.Vector;
  for f = fieldnames(v)'
    fields.addElement(char(f));
    jv.put(char(f),toJava(getfield(v,char(f))));
  end
  jv.put('_fields',fields);
  return
end

author: andrewssobral

repo: lrslibrary


text:
%% Dynamic Matrix Control Tutorial
% Dynamic Matrix Control (DMC) was the first Model Predictive Control (MPC)
% algorithm introduced in early 1980s. Nowadays, DMC is available in almost
% all commercial industrial distributed control systems and process
% simulation software packages. This tutorial intends to explain the
% features of DMC using the dmc function developed by the author.

%% Example: A Water Heater
% Consider a water heater as shown in the following figure, where the cold
% water is heated by means of a gas burner. The aim of DMC is by
% manipulating valve to control the gas flow so that the outlet temperature
% is at desired level.
%%
%
% <<WaterHeater.PNG>>
%

%% The Step Response Model
% The DMC algorithm works with a step response model. The first step to
% design a DMC controller is to perform a step test on the plant to
% generate a step response model. The step response of the water heater is
% obtained through such a test and given as follows:
p.sr = [0;0;0.271;0.498;0.687;0.845;0.977;1.087;1.179;1.256;...
        1.320;1.374;1.419;1.456;1.487;1.513;1.535;1.553;1.565;1.581;...
        1.592;1.600;1.608;1.614;1.619;1.632;1.627;1.630;1.633;1.635];

% The response is converged after 30 steps.

%% DMC without Setpoint Prediction
% Setup the DMC
N=120;  % Total simulation length (samples)
p.p=10; % Prediction horizon
p.m=5;  % Moving horizon
p.la=1; % Control weight
% Reference (setpoint)
R=[ones(30,1);zeros(30,1);ones(30,1);zeros(30,1)];
p.y=0;  % Initial output
p.v=[]; % empty past input to indicate initialization
% buffer of input to cope with time delay
u=zeros(3,1);
% Initialization of variables for results
Y=zeros(N,1);
U=zeros(N,1);
% DMC Simulation
for k=1:120
    p.a=0;
    p.r=R(k);   % DMC only knows current setpoint
    if k>60     % change smoothing factor for second half simulation
        p.a=0.7;
    end
    p=dmc(p);
    Y(k)=p.y;
    U(k)=p.u;
    u=[u(2:3);p.u];
    p.y=0.8351*p.y+0.2713*u(1); % actual plant output
end
% DMC results
subplot(211)
plot(1:N,Y,'b-',1:N,R,'r--',[60 60],[-0.5 1.5],':','linewidth',2)
title('solid: output, dashed: reference')
text(35,1,'\alpha=0')
text(95,1,'\alpha=0.7')
axis([0 120 -0.5 1.5])
subplot(212)
[xx,yy]=stairs(1:N,U);
plot(xx,yy,'-',[60 60],[-0.5 1.5],':','linewidth',2)
axis([0 120 -0.5 1.5])
title('input, \lambda=1')
xlabel('time, min')

%% DMC with Setpoint Prediction
% Setup the DMC
N=120;  % Total simulation length (samples)
p.p=10; % Prediction horizon
p.m=5;  % Moving horizon
p.la=1; % Control weight
% Reference (setpoint)
R=[ones(30,1);zeros(30,1);ones(30,1);zeros(30,1)];
p.y=0;  % Initial output
p.v=[]; % empty past input to indicate initialization
% buffer of input to cope with time delay
u=zeros(3,1);
% Initialization of variables for results
Y=zeros(N,1);
U=zeros(N,1);
% DMC Simulation
for k=1:120
    p.a=0;
    p.r=R(k:min(N,k+p.p)); % DMC knows future setpoint
    if k>60     % change smoothing factor for second half simulation
        p.a=0.7;
    end
    p=dmc(p);
    Y(k)=p.y;
    U(k)=p.u;
    u=[u(2:3);p.u];
    p.y=0.8351*p.y+0.2713*u(1); % actual plant output
end
% DMC results
subplot(211)
plot(1:N,Y,'b-',1:N,R,'r--',[60 60],[-0.5 1.5],':','linewidth',2)
title('solid: output, dashed: reference')
text(35,1,'\alpha=0')
text(95,1,'\alpha=0.7')
axis([0 120 -0.5 1.5])
subplot(212)
[xx,yy]=stairs(1:N,U);
plot(xx,yy,'-',[60 60],[-0.5 1.5],':','linewidth',2)
axis([0 120 -0.5 1.5])
title('input, \lambda=1')
xlabel('time, min')

%% DMC with Different Control Weight
% Setup the DMC
N=120;  % Total simulation length (samples)
p.p=10; % Prediction horizon
p.m=5;  % Moving horizon
p.la=0.1; % Control weight
% Reference (setpoint)
R=[ones(30,1);zeros(30,1);ones(30,1);zeros(30,1)];
p.y=0;  % Initial output
p.v=[]; % empty past input to indicate initialization
% buffer of input to cope with time delay
u=zeros(3,1);
% Initialization of variables for results
Y=zeros(N,1);
U=zeros(N,1);
% DMC Simulation
for k=1:120
    p.a=0;
    p.r=R(k:min(N,k+p.p)); % DMC knows future setpoint
    if k>60     % change smoothing factor for second half simulation
        p.a=0.7;
    end
    p=dmc(p);
    Y(k)=p.y;
    U(k)=p.u;
    u=[u(2:3);p.u];
    p.y=0.8351*p.y+0.2713*u(1); % actual plant output
end
% DMC results
subplot(211)
plot(1:N,Y,'b-',1:N,R,'r--',[60 60],[-0.5 1.5],':','linewidth',2)
title('solid: output, dashed: reference')
text(35,1,'\alpha=0')
text(95,1,'\alpha=0.7')
axis([0 120 -0.5 1.5])
subplot(212)
[xx,yy]=stairs(1:N,U);
plot(xx,yy,'-',[60 60],[-0.5 1.5],':','linewidth',2)
axis([0 120 -1 2])
title('input, \lambda=0.1')
xlabel('time, min')

%% Reference
% Camacho, E.F. and Bordons, C., Model Predictive Control, Springer-Verlag,
% 1999.


text:
function result = nii_read_volume(info)
% function for reading volume of NifTi ( .nii ) volume file
% nii_read_header(file-info)
%
% volume = nii_read_volume(file-header)
%
% examples:
% 1: info = nii_read_header()
%    V = nii_read_volume(info);
%    imshow(squeeze(V(:,:,round(end/2))),[]);
%
% 2: V = nii_read_volume('test.nii');

if(~isstruct(info)) info=nii_read_header(info); end

% Open v3d file
fid=fopen(info.Filename,'rb');

  % Seek volume data start
  datasize=prod(info.Dimensions)*(info.BitVoxel/8);
  fsize=info.Filesize;
  fseek(fid,fsize-datasize,'bof');

  % Read Volume data
  switch(info.DataTypeStr)
      case 'INT8'
        V = int8(fread(fid,datasize,'int8'));
      case 'UINT8'
        V = uint8(fread(fid,datasize,'uint8'));
      case 'INT16'
        V = int16(fread(fid,datasize,'int16'));
      case 'UINT16'
        V = uint16(fread(fid,datasize,'uint16'));
      case 'INT32'
        V = int32(fread(fid,datasize,'int32'));
      case 'UINT32'
        V = uint32(fread(fid,datasize,'uint32'));
      case 'INT64'
        V = int64(fread(fid,datasize,'int64'));
      case 'UINT64'
        V = uint64(fread(fid,datasize,'uint64'));
      case 'FLOAT'
        V = single(fread(fid,datasize,'single'));
      otherwise
        V = uint8(fread(fid,datasize,'uint8'));
  end
fclose(fid);

% Reshape the volume data to the right dimensions
result = reshape(V,info.Dimensions(1:3));

author: tomdoel

repo: pulmonarytoolkit


text:
function npartitions = partn_enum ( n, nmax )

%*****************************************************************************80
%
%% PARTN_ENUM enumerates the partitions of N with maximum element NMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    Normally N must be positive, but for this routine any
%    N is allowed.
%
%    Input, integer NMAX, the maximum element in the partition.
%    Normally, 1 <= NMAX <= N is required,
%    but for this routine any value of NMAX is allowed.
%
%    Output, integer NPARTITIONS is the number of partitions of N
%    with maximum element NMAX.
%
  if ( n <= 0 )

    npartitions = 0;

  elseif ( nmax <= 0 || n < nmax )

    npartitions = 0;

  else

    offset = 1;

    p = npart_table ( n, nmax );

    npartitions = p(n+offset,nmax+offset);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function alpha_pq = moment_normalized ( n, x, y, p, q )

%*****************************************************************************80
%
%% MOMENT_NORMALIZED computes a normalized moment of a polygon.
%
%  Discussion:
%
%    Alpha(P,Q) = Integral ( x, y in polygon ) x^p y^q dx dy / Area ( polygon )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carsten Steger,
%    On the calculation of arbitrary moments of polygons,
%    Technical Report FGBV-96-05,
%    Forschungsgruppe Bildverstehen, Informatik IX,
%    Technische Universitaet Muenchen, October 1996.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real X(N), Y(N), the vertex coordinates.
%
%    Input, integer P, Q, the indices of the moment.
%
%    Output, real ALPHA_PQ, the normalized moment Alpha(P,Q).
%
  nu_pq = moment ( n, x, y, p, q );
  nu_00 = moment ( n, x, y, 0, 0 );

  alpha_pq = nu_pq / nu_00;

  return
end


text:
function [PowLoss] = PowerLoss(modres, modfit, moddf, tc, TR, Run, alpha)
% function [PowLoss] = PowerLoss(modres, modfit, moddf, tc, TR, Run, alpha)
%
% Estimates Power-loss due to mis-modeling.
%
% INPUT:
%
% modres - residuals
% modfit - model fit
% moddf  - model degrees of freedom
% tc     - time course
% TR     - time resolution
% Runs   - expermental design
% alpha  - alpha value
%
% OUTPUT:
%
% PowLoss - Estimated power loss
%
%

len = length(tc);               % length of time course
%T = round(30./TR);              % length of estimated HRF
T = 30;
tstar = tinv(1-alpha,moddf);    % t-threshold

% Fit FIR model to find 'baseline' power.
[h, fit, e] = Fit_sFIR(tc,TR,Run,T,1);
s = (1/(len-T))*e'*e;
t = 1/sqrt(s*inv(fit'*fit));
basePow = 1- nctcdf(tstar,(len-T),t);

% Compute model power.
sig = (1/moddf)*modres'*modres;
ts = 1/sqrt(sig*inv(modfit'*modfit));
modPow = 1- nctcdf(tstar,moddf,ts);

% Compute 'power loss'
PowLoss = basePow - modPow;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = sqdread(filename,varargin)
% [DATA,INFO] = SQDREAD(FILENAME,'CHANNELS',[C1 C2 C3,...],'SAMPLES',[S1 SN],'FORMAT','double');    % READ SQD-FILE
%
% SQDREAD Get Channel data and info saved in sqd file 'filename' using
% specifications from other options like 'Info','Samples', 'Channels',
% 'Format'. When the format is 'double', the MEG data has units of picotesla (pT).
% Usage:
% info = sqdread(filename,'Info');
%   where 'info' belongs to class sqdhandle. Use GET and SET methods to review or
%   change object properties. It contains experiment parameter values and
%   data acquisition parameters.
%
% [data,info] = sqdread(filename);
%   Get all the raw data saved in 'filename'. Data is arranged as following:
%	|Chan(0)Sample(1)   Chan(1)Sample(1)	... Chan(end-1)Sample(1)  |
%	|Chan(0)Sample(2)   Chan(1)Sample(2)	...	        .	        |
%	|   .			        .			                .	        |
%	|   .			        .			                .           |
%	|Chan(0)Sample(end) Chan(1)Sample(end)	... Chan(end-1)Sample(end)|
%
% [data,info] = sqdread(filename,'Channels',[0 1 2 3 4],'Samples',[2001 3000],...
%	'Format','int16');
%   Get data specified by Channels,Samples, and Format fields of inputs, ie
%   in the above example, data is returned for Channel numbers [0,1,2,3,4] and Samples
%   from 2001 to 3000. The data is returned as 'INT16'. Data will be arranged
%   as shown above in a increasing order. You may choose above options in any
%   order. If option is not specified, default values are assumed.
%   Default values for :
%       'Channels' - all available channels
%       'Samples'  - all available samples
%       'Format'   - 'double'
%
% Examples:
% If the path for the sqd-toolbox has been set properly and all of the files present
% (see README), there should be a sample raw sqd-file "testsqd.sqd" for testing. The
% sqd-file "testsqd.sqd" has data for 192 channels, 10 seconds of data and sampled
% at 1000Hz. The data is sinusoids of varying frequencies for the first 157 channels
% and random triggers for the rest. For the examples below, "testsqd.sqd" is used as
% the test sqd-file.
%
% 1. To read only the information regarding the data acquisition parameters:
%   a. info = sqdread('testsqd.sqd','info');
%       This will return you "info" which contains all the acquisition information.
%       Please make sure that when you specify the filename, it either includes the
%       filepath or you are in the directory where the file exists.
%   b. get(info);
%       This will tell you the parameters included in "info"
%   c. To get value of a parameter "fieldname" of "info".Choose fieldname
%      from parameters returned from get(info).
%      Ex.
%       info.PatientInfo        % this returns the PatientInfo = Matt Cheely,...
%       info.InputGain          % this returns the Amplifier InputGain = 2
%       info.SamplesAvailable   % this returns number of samples available per channel = 10000
%      In general, to get value for parameter "fieldname"
%       info.fieldname
%       get(info,'fieldname')
%
% 2. To read all the data:
%   data = sqdread('testsqd.sqd');
%       This will return you a 2-Dimensional matrix array "data" whose size will
%   be 10000 x 192 ie (Number of Samples x Number of Channels)
%
% 3. To read data and info together:
%   [data,info] = sqdread('testsqd.sqd');
%
% 4. To read in data for only specific channels and/or samples
%   a. To read in all the data for only channel N
%      N = 0;                                       % First channel
%      data = sqdread('testsqd.sqd','Channels',N);  % Read in the data
%      plot(data);                                  % To look at the signal
%      size(data)                                   % should be 10000 x 1
%                                                   % ie  (NumSamples x NumChannels)
%   b. To read in the data for all channels but only from Sample1 to SampleN
%      Sample1 = 1;                                 % First Sample
%      SampleN = 1000;                              % 1000th sample = 1second of data
%      data = sqdread('testsqd.sqd','Samples',[Sample1 SampleN]); % Read in the data
%      size(data)                                   % should be 1000 x 192
%                                                   % ie (NumSamples x NumChannels)
%   c. To read in data for channels N1,N2 and samples from Sample1 to SampleN
%      N1 = 0;                                      % First Channel
%      N2 = 99;                                     % 100th channel
%      Sample1 = 1;                                 % First Sample
%      SampleN = 1000;                              % 1000th sample = 1second of data
%      % Read in the data
%      data = sqdread('testsqd.sqd','Channels',[N1,N2],'Samples',[Sample1 SampleN]);
%      size(data)                                   % should be 1000 x 2
%                                                   % ie (NumSamples x NumChannels)
%
%
%   SEE ALSO,
%   SQDHANDLE,GETDATA,GET,SET

% This function is just a wrapper around the methods of CLASS SQDHANDLE
% It uses the filename to make a sqdhandle_obj and then calls its methods
% to return the information asked for. Very little error correction is done
% here. It is assumed that the methods of SQDHANDLE will take care of it.

% Get filename
% If cannot find file, ask user input
if nargin<1,
    filenameqn = input(...
        'Do you wish to choose an sqd-file to open? [Return,Yes|No]',...
        's');
    if isempty(filenameqn)|strncmpi(filenameqn,'y',1),
        [f,p] = uigetfile({'*.sqd;*.SQD'},'Get SQD-FILE');
        if ((~isstr(f)) | ~min(size(f))), return;end;  % If Load cancelled
        filename = [p f];
    else
        error('Valid sqd-filename required');
    end;
end;

if ischar(filename)
    info = sqdhandle(filename);                     % Get sqdhandle_obj from file
    switch nargin-1
        case 0                                      % Ex. data = sqdread(file.sqd);
            varargout{1} = getdata(info);           % Get data from sqdhandle_obj
            if nargout>=2                           % Ex. [data,info] = sqdread(file.sqd);
                varargout{2} = info;                % Return sqdhandle_obj as info
            end;
        case 1
            if strcmpi(varargin{1},'info')          % Ex. info = sqdread(file.sqd,'info');
                varargout{1} = info;                % Return sqdhandle_obj as info
            else
                error('Incorrect inputs')
            end;
        otherwise
            % Ex. data = sqdread(file.sqd,'Channels',[1 2],'Samples',[1 1000]);
            varargout{1} = getdata(info,varargin{:});
            if nargout>=2
                varargout{2} = info;
            end;
    end;
else
   error('First argument must be a valid filename');
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function niter = loopy_converged(engine)
% LOOPY_CONVERGED Did loopy belief propagation converge? 0 means no, eles we return the num. iterations.
% function niter = loopy_converged(engine)
%
% We use a simple heuristic: we say convergence occurred if the number of iterations
% used was less than the maximum allowed.

if engine.niter == engine.max_iter
  niter = 0;
else
  niter = engine.niter;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function addnoise_asl(cleanfile, noisefile, outfile, snr)
% ----------------------------------------------------------------------
%   This function adds noise to a file at a specified SNR level. It uses
%   the active speech level to compute the speech energy. The
%   active speech level is computed as per ITU-T P.56 standard [1].
%
%   Usage:  addnoise_asl(cleanFile.wav, noiseFile.wav, noisyFile.wav, SNR)
%
%         cleanFile.wav  - clean input file in .wav format
%         noiseFile.wav  - file containing the noise signal in .wav format
%         noisyFile.wav  - resulting noisy file
%         SNR            - desired SNR in dB
%
%   Note that if the variable IRS below (line 38) is set to 1, then it applies the IRS
%   filter to bandlimit the signal to 300 Hz - 3.2 kHz. The default IRS
%   value is 0, ie, no IRS filtering is applied.
%
%  Example call:
%       addnoise_asl('sp04.wav','white_noise.wav','sp04_white_5db.wav',5);
%
%
%  References:
%   [1] ITU-T (1993). Objective measurement of active speech level. ITU-T
%       Recommendation P. 56
%
%   Author: Yi Hu and Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
% ----------------------------------------------------------------------

if nargin ~=4
    fprintf('USAGE: addnoise_asl(cleanFile.wav, noiseFile.wav, noisyFile.wav, SNR) \n');
    fprintf('For more help, type: help addnoise_asl\n\n');
    return;
end

IRS=0;  % if 1 apply IRS filter simulating telephone handset bandwidth (300 Hz -3.2 kHz)

% wavread gives floating point column data
[clean, srate, nbits]= wavread(cleanfile);
% filter clean speech with irs filter
if IRS==1, clean= apply_IRS( clean, srate, nbits); end;

[Px, asl, c0]= asl_P56 ( clean, srate, nbits);
% Px is the active speech level ms energy, asl is the active factor, and c0
% is the active speech level threshold.

x=clean;
x_len= length( x); % length of speech signal

[noise, srate1, nbits1]= wavread( noisefile);
if (srate1~= srate)| (nbits1~= nbits)
    error( 'the formats of the two files dont match!');
end
noise_len= length( noise);
if (noise_len<= x_len)
    error( 'the noise length has to be greater than speech length!');
end

rand_start_limit= noise_len- x_len+ 1;
% the start of the noise segment can vary between [1 rand_start_limit]
rand_start= round( (rand_start_limit- 1)* rand( 1)+ 1);
% random start of the noise segment
noise_segment= noise( rand_start: rand_start+ x_len- 1);

if IRS==1, noise_segment= apply_IRS( noise_segment, srate, nbits); end;

% this is the randomly selected noise segment that will be added to the
% clean speech x
Pn= noise_segment'* noise_segment/ x_len;
% we need to scale the noise segment samples to obtain the desired snr= 10*
% log10( Px/ (sf^2 * Pn))
sf= sqrt( Px/Pn/ (10^ (snr/ 10))); % scale factor for noise segment data
noise_segment= noise_segment * sf;

noisy = x+ noise_segment;

if ( (max( noisy)>= 1) | (min( noisy)< -1))
    error( 'Overflow occurred!\n');
end;

wavwrite( noisy, srate, nbits, outfile);

fprintf( 1, '\n NOTE: For comparison, the SNR based on long-term RMS level is %4.2f dB.\n\n', 10*log10((x'*x)/ ...
     (noise_segment'*noise_segment)));

%------------------------------------------------------------------------
function data_filtered= apply_IRS( data, Fs, nbits);

n= length( data);

% now find the next power of 2 which is greater or equal to n
pow_of_2= 2^ (ceil( log2( n)));

align_filter_dB= [0, -200; 50, -40; 100, -20; 125, -12; 160, -6; 200, 0;...
    250, 4; 300, 6; 350, 8; 400, 10; 500, 11; 600, 12; 700, 12; 800, 12;...
    1000, 12; 1300, 12; 1600, 12; 2000, 12; 2500, 12; 3000, 12; 3250, 12;...
    3500, 4; 4000, -200; 5000, -200; 6300, -200; 8000, -200];

[number_of_points, trivial]= size( align_filter_dB);
overallGainFilter= interp1( align_filter_dB( :, 1), align_filter_dB( :, 2), ...
    1000);

x= zeros( 1, pow_of_2);
x( 1: n)= data;

x_fft= fft( x, pow_of_2);

freq_resolution= Fs/ pow_of_2;

factorDb( 1: pow_of_2/2+ 1)= interp1( align_filter_dB( :, 1), ...
    align_filter_dB( :, 2), (0: pow_of_2/2)* freq_resolution)- ...
    overallGainFilter;
factor= 10.^ (factorDb/ 20);

factor= [factor, fliplr( factor( 2: pow_of_2/2))];
x_fft= x_fft.* factor;

y= ifft( x_fft, pow_of_2);

data_filtered= y( 1: n)';

function [asl_ms, asl, c0]= asl_P56 ( x, fs, nbits)
% this implements ITU P.56 method B.
% 'speechfile' is the speech file to calculate active speech level for,
% 'asl' is the active speech level (between 0 and 1),
% 'asl_rms' is the active speech level mean square energy.

% x is the column vector of floating point speech data

x= x(:); % make sure x is column vector
T= 0.03; % time constant of smoothing, in seconds
H= 0.2; % hangover time in seconds
M= 15.9;
% margin in dB of the difference between threshold and active speech level
thres_no= nbits- 1; % number of thresholds, for 16 bit, it's 15

I= ceil( fs* H); % hangover in samples
g= exp( -1/( fs* T)); % smoothing factor in envelop detection
c( 1: thres_no)= 2.^ (-15: thres_no- 16);
% vector with thresholds from one quantizing level up to half the maximum
% code, at a step of 2, in the case of 16bit samples, from 2^-15 to 0.5;
a( 1: thres_no)= 0; % activity counter for each level threshold
hang( 1: thres_no)= I; % hangover counter for each level threshold

sq= x'* x; % long-term level square energy of x
x_len= length( x); % length of x

% use a 2nd order IIR filter to detect the envelope q
x_abs= abs( x);
p= filter( 1-g, [1 -g], x_abs);
q= filter( 1-g, [1 -g], p);

for k= 1: x_len
    for j= 1: thres_no
        if (q(k)>= c(j))
            a(j)= a(j)+ 1;
            hang(j)= 0;
        elseif (hang(j)< I)
            a(j)= a(j)+ 1;
            hang(j)= hang(j)+ 1;
        else
            break;
        end
    end
end

asl= 0;
asl_rms= 0;
if (a(1)== 0)
    return;
else
    AdB1= 10* log10( sq/ a(1)+ eps);
end

CdB1= 20* log10( c(1)+ eps);
if (AdB1- CdB1< M)
    return;
end

AdB(1)= AdB1;
CdB(1)= CdB1;
Delta(1)= AdB1- CdB1;

for j= 2: thres_no
    AdB(j)= 10* log10( sq/ (a(j)+ eps)+ eps);
    CdB(j)= 20* log10( c(j)+ eps);
end

for j= 2: thres_no
    if (a(j) ~= 0)
        Delta(j)= AdB(j)- CdB(j);
        if (Delta(j)<= M)
            % interpolate to find the asl
            [asl_ms_log, cl0]= bin_interp( AdB(j), ...
                AdB(j-1), CdB(j), CdB(j-1), M, 0.5);
            asl_ms= 10^ (asl_ms_log/ 10);
            asl= (sq/ x_len)/ asl_ms;
            c0= 10^( cl0/ 20);
            break;
        end
    end
end

function [asl_ms_log, cc]= bin_interp(upcount, lwcount, ...
    upthr, lwthr, Margin, tol)

if (tol < 0)
    tol = -tol;
end

% Check if extreme counts are not already the true active value
iterno = 1;
if (abs(upcount - upthr - Margin) < tol)
    asl_ms_log= upcount;
    cc= upthr;
    return;
end
if (abs(lwcount - lwthr - Margin) < tol)
    asl_ms_log= lwcount;
    cc= lwthr;
    return;
end

% Initialize first middle for given (initial) bounds
midcount = (upcount + lwcount) / 2.0;
midthr = (upthr + lwthr) / 2.0;

% Repeats loop until `diff' falls inside the tolerance (-tol<=diff<=tol)
while ( 1)

    diff= midcount- midthr- Margin;
    if (abs(diff)<= tol)
        break;
    end

    % if tolerance is not met up to 20 iteractions, then relax the
    % tolerance by 10%

    iterno= iterno+ 1;

    if (iterno>20)
      tol = tol* 1.1;
    end

    if (diff> tol)   % then new bounds are ...
        midcount = (upcount + midcount) / 2.0;
        % upper and middle activities
        midthr = (upthr + midthr) / 2.0;
        % ... and thresholds
    elseif (diff< -tol)	% then new bounds are ...
        midcount = (midcount + lwcount) / 2.0;
        % middle and lower activities
        midthr = (midthr + lwthr) / 2.0;
        % ... and thresholds
    end

end
%   Since the tolerance has been satisfied, midcount is selected
%   as the interpolated value with a tol [dB] tolerance.

asl_ms_log= midcount;
cc= midthr;

author: jtkim-kaist

repo: Speech-enhancement


text:
% run_exiftool_demo
% Demonstration of functions
%    getexif
%    putexif
% Peter Burns, 28 May 2013
%              22 July fixed minor mistake calling putexif

[f,p] = uigetfile({'*.jpg';'*.tif'}, ...
'Select image file');

% Read Exif data
exifdata = getexif([p,f]);
disp(exifdata)

dat = imread([p,f]);
imagesc(dat), axis image

%
% Do processing image data here ...
% (For this demo. we save the original data)
%

[f1,p1] = uiputfile({'*.jpg';'*.tif'},'Save as');

% Save image with original Exif data
status = putexif(dat,[p1,f1], [p,f]);

% Check
[exifdata1, nf1] = getexif([p,f]);
disp(exifdata1)
[exifdata2, nf2] = getexif([p1,f1]);
disp(exifdata2)


text:
function semantic_segmentation_rgb(fileName, sensorType, varName)

    %% Check if input is valid
    if ~isstring(fileName) && ~ischar(fileName)
        error("File name must a string or char array\n")
    end

    if ~strcmp(sensorType, "semantic_segmentation_rgb")
       error("Wrong sensor selected\n")
    end

    if ~isvarname(varName)
        error("Invalid variable name\n")
    end

    %% Create a python file
    file = fopen(strcat(fileName, '.py'), 'w');

    % Automatically genrates a python file containing the sensor call back
    % bindings
    fprintf(file, 'import numpy as np\n');
    fprintf(file, 'from collections import defaultdict\n');
    fprintf(file, '\n');
    fprintf(file, '# Converts the tags to human observable colors\n');
    fprintf(file, '_mapDict = defaultdict(lambda: (255, 255, 255), {\n');
    fprintf(file, '    0: (0, 0, 0),\n');
    fprintf(file, '    0: (0, 0, 0),\n');
    fprintf(file, '    1: (70, 70, 70),\n');
    fprintf(file, '    2: (100, 40, 40),\n');
    fprintf(file, '    3: (55, 90, 80),\n');
    fprintf(file, '    4: (220, 20, 60),\n');
    fprintf(file, '    5: (153, 153, 153),\n');
    fprintf(file, '    6: (157, 234, 50),\n');
    fprintf(file, '    7: (128, 64, 128),\n');
    fprintf(file, '    8: (244, 35, 232),\n');
    fprintf(file, '    9: (107, 142, 35),\n');
    fprintf(file, '    10: (0, 0, 142),\n');
    fprintf(file, '    11: (102, 102, 156),\n');
    fprintf(file, '    12: (220, 220, 0),\n');
    fprintf(file, '    13: (70, 130, 180),\n');
    fprintf(file, '    14: (81, 0, 81),\n');
    fprintf(file, '    15: (150, 100, 100),\n');
    fprintf(file, '    16: (230, 150, 140),\n');
    fprintf(file, '    17: (180, 165, 180),\n');
    fprintf(file, '    18: (250, 170, 30),\n');
    fprintf(file, '    19: (110, 190, 160),\n');
    fprintf(file, '    20: (170, 120, 50),\n');
    fprintf(file, '    21: (45, 60, 150),\n');
    fprintf(file, '    22: (145, 170, 100),\n');
    fprintf(file, '})\n');
    fprintf(file, '\n');
    fprintf(file, 'def _setup_mapping():\n');
    fprintf(file, '    global _mapDict\n');
    fprintf(file, '    global _mapping_array\n');
    fprintf(file, '\n');
    fprintf(file, '    tags = np.array(list(_mapDict.keys()))\n');
    fprintf(file, '    colors = np.array(list(_mapDict.values()))\n');
    fprintf(file, '\n');
    fprintf(file, '    _mapping_array = np.zeros((tags.max()+1, 3), dtype="uint8")\n');
    fprintf(file, '    _mapping_array[tags] = colors\n');
    fprintf(file, '\n');

    fprintf(file, 'def bindSensor(sensor):\n');
    fprintf(file, '    _setup_mapping()\n');
    fprintf(file, '    sensor.listen(lambda _image: do_something(_image))\n');
    fprintf(file, '\n');
    fprintf(file, 'def do_something(_image):\n');
    fprintf(file, '    global %s\n', varName);
    fprintf(file, '    global _mapping_array\n');
    fprintf(file, '\n');
    fprintf(file, '    data = np.frombuffer(_image.raw_data, dtype=np.dtype("uint8"))\n');
    fprintf(file, '\n');
    fprintf(file, '    # Get the red channel which has the tags\n');
    fprintf(file, '    data = np.reshape(data, (_image.height, _image.width, 4))\n');
    fprintf(file, '    data = data[:, :, 2]\n');
    fprintf(file, '\n');
    fprintf(file, '    # Map the tags to their respective RGB colors\n');
    fprintf(file, '    data = _mapping_array[data]\n');
    fprintf(file, '\n');
    fprintf(file, '    # Convert the data into MATLAB cast compatible type\n');
    fprintf(file, '    %s = np.ascontiguousarray(data)\n', varName);

    fclose(file);
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function value = h_hofstadter ( n )

%*****************************************************************************80
%
%% H_HOFSTADTER computes the Hofstadter H sequence.
%
%  Discussion:
%
%    H(N) = 0                          if N = 0
%         = N - H ( H ( H ( N - 1 ) ), otherwise.
%
%    H(N) is defined for all nonnegative integers.
%
%  Table:
%
%     N  H(N)
%    --  ----
%
%     0     0
%     1     1
%     2     1
%     3     2
%     4     3
%     5     4
%     6     4
%     7     5
%     8     5
%     9     6
%    10     7
%    11     7
%    12     8
%    13     9
%    14    10
%    15    10
%    16    11
%    17    12
%    18    13
%    19    13
%    20    14
%    21    14
%    22    15
%    23    16
%    24    17
%    25    17
%    26    18
%    27    18
%    28    19
%    29    20
%    30    20
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Douglas Hofstadter,
%    Goedel, Escher, Bach,
%    Basic Books, 1979.
%
%  Parameters:
%
%    Input, integer N, the argument of the function.
%
%    Output, integer VALUE, the value of the function.
%
  if ( n <= 0 )
    value = 0;
  else
    value = n - h_hofstadter ( h_hofstadter ( h_hofstadter ( n-1 ) ) );
  end

  return
end


text:
function []=panel6irfdisp(N,n,Units,endo,irf_estimates,strshocks_estimates,IRFperiods,IRFt,stringdates1,T,decimaldates1,pref)

% IRFs
if pref.plot
% plot the figure
irf=figure('Tag','BEARresults');
set(irf,'Color',[0.9 0.9 0.9]);
   if IRFt==1
   set(irf,'name',['impulse response functions (no structural identifcation)']);
   elseif IRFt==2
   set(irf,'name',['impulse response functions (structural identification by Choleski ordering)']);
   elseif IRFt==3
   set(irf,'name',['impulse response functions (structural identification by triangular factorisation)']);
   end
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
      % loop over units (for shocks)
      for kk=1:N
         % loop over shocks
         for ll=1:n
         % increment count
         count=count+1;
         % then plot
         subplot(N*n,N*n,count)
         temp=irf_estimates{(ii-1)*n+jj,(kk-1)*n+ll};
         hold on
         Xpatch=[(1:IRFperiods) (IRFperiods:-1:1)];
         Ypatch=[temp(1,:) fliplr(temp(3,:))];
         IRFpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
         set(IRFpatch,'facealpha',0.5);
         set(IRFpatch,'edgecolor','none');
         plot(temp(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
         plot([1,IRFperiods],[0 0],'k--');
         hold off
         minband=min(temp(1,:));
         maxband=max(temp(3,:));
         space=maxband-minband;
         Ymin=minband-0.2*space;
         Ymax=maxband+0.2*space;
         set(gca,'XLim',[1 IRFperiods],'YLim',[Ymin Ymax],'FontName','Times New Roman');
            % top labels
            if jj==1
            title([Units{kk,1} '\_' endo{ll,1}],'FontWeight','normal');
            end
            % side labels
            if ll==1
            ylabel([Units{ii,1} '\_' endo{jj,1}],'FontWeight','normal');
            end
         end
      end
   end
end
% top supertitle
ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
set(get(ax,'Title'),'Visible','on')
title('Shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
% side supertitle
ylabel('Response of:','FontSize',12,'FontName','Times New Roman','FontWeight','normal');
set(get(ax,'Ylabel'),'Visible','on');
end
% save on Excel
% create the cell that will be saved on excel
irfcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},IRFperiods+3,1);
horzspace=repmat({''},2,5*N*n+N);
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % initiate the cell of results
   rowcell={};
      % loop over units (for shocks)
      for kk=1:N
      endocell={};
         % loop over shocks
         for ll=1:n
         % create a header
         header=[{['response of ' Units{ii,1} '_' endo{jj,1} ' to ' Units{kk,1} '_' endo{ll,1} ' shocks']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
         % complete the cell
         tempcell=[[header;num2cell((1:IRFperiods)') num2cell((irf_estimates{(ii-1)*n+jj,(kk-1)*n+ll})')] vertspace];
         % concatenate to the previous parts of unitcell
         endocell=[endocell tempcell];
         end
      rowcell=[rowcell vertspace endocell];
      end
   % concatenate to the previous parts of irfcell
   irfcell=[irfcell;horzspace;rowcell];
   end
irfcell=[irfcell;horzspace];
end
% trim
irfcell=irfcell(3:end-2,2:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),irfcell,'IRF','B2');
end

% structural shocks

% generate plot, if there is any structural decomposition
if IRFt~=1
    if pref.plot
% plot
strshocks=figure('Tag','BEARresults');
set(strshocks,'Color',[0.9 0.9 0.9]);
set(strshocks,'name','structural shocks')
% initiate the count
count=0;
   % loop over units
   for ii=1:N
      % loop over endogenous variables
      for jj=1:n
      % increment count
      count=count+1;
      % then plot
      subplot(N,n,count)
      hold on
      Xpatch=[decimaldates1' fliplr(decimaldates1')];
      Ypatch=[strshocks_estimates{jj,1,ii}(1,:) fliplr(strshocks_estimates{jj,1,ii}(3,:))];
      Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
      set(Fpatch,'facealpha',0.6);
      set(Fpatch,'edgecolor','none');
      strs=plot(decimaldates1,strshocks_estimates{jj,1,ii}(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
      plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
      hold off
      minband=min(strshocks_estimates{jj,1,ii}(1,:));
      maxband=max(strshocks_estimates{jj,1,ii}(3,:));
      space=maxband-minband;
      Ymin=minband-0.2*space;
      Ymax=maxband+0.2*space;
      set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'YLim',[Ymin,Ymax],'FontName','Times New Roman');
         % top labels
         if count<=n
         title(endo{count,1},'FontWeight','normal');
         end
         % side labels
         if jj==1
         ylabel(Units{ii,1},'FontWeight','normal');
         end
      end
   end
    end
% save on Excel
% create the cell that will be saved on excel
strshockcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},3,5*n);
   % loop over units
   for ii=1:N
   % initiate the cell of results
   unitcell={};
      % loop over endogenous variables (horizontal dimension)
      for jj=1:n
      % create a header
      header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
      % complete the cell
      endocell=[[header;stringdates1 num2cell(strshocks_estimates{jj,:,ii}')] vertspace];
      % concatenate to the previous parts of unitcell
      unitcell=[unitcell endocell];
      end
   % concatenate to the previous parts of afcell
   strshockcell=[strshockcell;horzspace;unitcell];
   end
% trim
strshockcell=strshockcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strshockcell,'shocks','B2');
end
end

author: european-central-bank

repo: BEAR-toolbox


text:
classdef DMSingleImageMetaInfo
    % DMSingleImageMetaInfo. A structure for holding information related to a
    %     single image file
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    properties
        ImagePath
        ImageFilename
        ImageFileFormat
        Modality
        Date
        Time

        PatientId
        StudyUid
        SeriesUid
        ImageUid

        PatientName
        StudyDescription
        SeriesDescription

        Version
    end

    methods
        function obj = DMSingleImageMetaInfo(path, filename, image_type, modality, date, time, ...
                patient_id, study_uid, series_uid, image_uid, ...
                patient_name, study_description, series_description)

            if nargin > 0
                obj.ImagePath = path;
                obj.ImageFilename = filename;
                obj.ImageFileFormat = image_type;
                obj.Modality = modality;
                obj.Date = date;
                obj.Time = time;

                obj.PatientId = patient_id;
                obj.StudyUid = study_uid;
                obj.SeriesUid = series_uid;
                obj.ImageUid = image_uid;

                obj.PatientName = patient_name;
                obj.StudyDescription = study_description;
                obj.SeriesDescription = series_description;
                obj.Version = 1;
            end
        end
    end
end


text:
function drawBox3D(h,occlusion,corners,face_idx,orientation)

% set styles for occlusion and truncation
occ_col    = {'w','g','y','r'};

% draw projected 3D bounding boxes
if ~isempty(corners)
  for f=1:4
    line([corners(1,face_idx(f,:)),corners(1,face_idx(f,1))]+1,...
         [corners(2,face_idx(f,:)),corners(2,face_idx(f,1))]+1,...
         'parent',h(2).axes, 'color',occ_col{occlusion+2},...
         'LineWidth',6);
    line([corners(1,face_idx(f,:)),corners(1,face_idx(f,1))]+1,...
         [corners(2,face_idx(f,:)),corners(2,face_idx(f,1))]+1,...
         'parent',h(2).axes,'color','k','LineWidth',2);
  end
end

% draw orientation vector
if ~isempty(orientation)
  line([orientation(1,:),orientation(1,:)]+1,...
       [orientation(2,:),orientation(2,:)]+1,...
       'parent',h(2).axes,'color','w','LineWidth',6);
  line([orientation(1,:),orientation(1,:)]+1,...
       [orientation(2,:),orientation(2,:)]+1,...
       'parent',h(2).axes,'color','k','LineWidth',2);
end

repo: MSCKF


text:
function evaluate_clustering(name, embedding_dimension, distance)

if nargin < 1
    name = 'liftedstructsim_softmax_pair_m128_multilabel';
end

if nargin < 2
    embedding_dimension = 64;
end

if nargin < 3
    distance = 'euclidean';
end

K = 11316;
filename = sprintf('idx_kmeans_googlenet_%s_embed%d_%s.mat', name, embedding_dimension, distance);
if exist(filename, 'file')
    object = load(filename);
    idx = object.idx;
else
    % compute similarity
    load(sprintf('validation_googlenet_feat_matrix_%s_embed%d_baselr_1E4_gaussian2k.mat', ...
                                        name, embedding_dimension), 'fc_embedding');
    X = double(fc_embedding');

    % kmeans clustering
    fprintf('2d kmeans %d\n', K);
    opts = struct('maxiters', 1000, 'mindelta', eps, 'verbose', 1);
    [center, sse] = vgg_kmeans(X, K, opts);
    [idx_kmeans, d] = vgg_nearest_neighbour(X, center);

    % construct idx
    num = size(X, 2);
    idx = zeros(num, 1);
    for i = 1:K
        index = find(idx_kmeans == i);
        [~, ind] = min(d(index));
        cid = index(ind);
        idx(index) = cid;
    end

    fprintf('Number of clusters: %d\n', length(unique(idx)));
    save(sprintf('idx_kmeans_googlenet_%s_embed%d_%s.mat', name, embedding_dimension, distance), 'idx');
end

% evaluation
num_validation_classes = K;

% load ground truth from filenames
[image_ids, class_ids, superclass_ids, path_list] = ...
    textread('/cvgl/group/Ebay_Dataset/Ebay_test.txt', '%d %d %d %s',...
    'headerlines', 1);

num = numel(class_ids);
item_ids = cell(num,1);
for i = 1:num
    class_id = class_ids(i);
    item_ids{i} = num2str(class_id);
end
assert(length(unique(item_ids)) == num_validation_classes);

% Given cluster assignment and the class names
%   Compute the three clustering metrics.
[NMI, RI, F1] = compute_clutering_metric(idx, item_ids);

fprintf('[method: %s, distance: %s] NMI: %.3f, RI: %.3f, F1: %.3f\n\n', ...
    name, distance, NMI, RI, F1);

author: rksltnl


text:
function M = rmPlotGUI_update(M)
% Refresh the Retinotopy Model Plotting GUI.
%
% M = rmPlotGUI_update(M);
%
% Update function for rmPlotGUI. Visualizes the receptive field and time
% series (data, prediction, and residual) for a given voxel for a
% retinotopy model fit.
%
% All the information about the model should be in the struct M, created
% with rmPlotGUI_getModel.
%
% ras, 09/06.
% ras, 11/08: broekn off into separate subfunction.
if ~exist('M','var') || isempty(M), M = get(gcf, 'UserData'); end

% get current voxel from the GUI
voxel = get(M.ui.voxel.sliderHandle, 'Value');

% get needed params
coords = M.coords(:,voxel);
if isequal(M.roi.viewType, 'Gray')  % convert coords into an index
    coords = M.coords(voxel);
end

M.modelNum = get(M.ui.model, 'Value');

%% compute RF, get tSeries for this voxel
[pred, RF, rfParams, variance_explained, blanks] = rmPlotGUI_makePrediction(M, coords);

% check if the prediction for this voxel is empty. If so, don't plot it:
if all( isnan(pred) | isinf(pred) )
	a = axes(M.ui.tsAxes);  cla(a);
	AX = axis;
	text(a, AX(2) + .5*(AX(2)-AX(1)), AX(3) + .5*(AX(4)-AX(3)), ...
		 '(No data available for this voxel)', ...
		 'FontSize', 14, ...
		 'HorizontalAlignment', 'center');
	return
end

% store the current prediction
M.prediction = pred;
M.RF = RF;

%%%%% VISUALIZE
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (1) plot time series
%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes(M.ui.tsAxes); cla; hold on;

hTSeries = plot(M.x, M.tSeries(:,voxel), 'k--', 'LineWidth', 1.5);
hFit = plot(M.x, pred(:, 1), 'b', 'LineWidth', 1);
hResidual = plot(M.x, M.tSeries(:,voxel)-pred(:,1), 'r:', 'LineWidth', 1);

xlim([min(M.x) max(M.x)]);

% Determine appropriate limits for y-axis. To do so, we concatenate the
% visible series. These include some or all of raw, predicted, and residual
% tseries. The user selects which are visible by checkbox in the GUI.
allPlotted = [];
if get(M.ui.tsCheck,   'Value'), allPlotted = cat(1, allPlotted, M.tSeries(:,voxel)); end
if get(M.ui.predCheck, 'Value'), allPlotted = cat(1, allPlotted, pred(:,1)); end
if get(M.ui.resCheck,  'Value'), allPlotted = cat(1, allPlotted, M.tSeries(:,voxel)-pred(:,1)); end
if max(allPlotted) > min(allPlotted), ylim([min(allPlotted) max(allPlotted)]); end

h = axis;
for n=1:numel(M.sepx),
    plot([1 1].*M.sepx(n), [h(3) h(4)], 'k:', 'LineWidth', 2);
end;
xlabel('Time (sec)');
ylabel('BOLD signal change (%)');

% set a button-down function which will allow the user to double-click on
% any time point and view the stimuls aperture at that time point.
bdf = ['M = get(gcf, ''UserData''); ' ...
	   'pt = get(gca, ''CurrentPoint''); ' ...
	   'pt = round(pt(1) / M.params.stim(1).framePeriod); ' ...
	   'st = get(gcf, ''SelectionType''); ' ...
	   'if isequal(st, ''open''), ' ...
	   '  set(M.ui.overlayStimCheck, ''Value'', 1); ' ...
	   '  mrvSliderSet(M.ui.time, ''Value'', pt); ' ...
	   '  rmPlotGUI(''update''); ' ...
	   'end; ' ...
	   'clear M st pt'];
set(M.ui.tsAxes, 'ButtonDownFcn', bdf);

% Indicate the time points with no stimulus
plot(M.x(blanks), zeros(sum(blanks),1), 'rx')

% set the user data of each checkbox to point to the PLOT curves
% we just created. The callbacks for these checkboxes will then toggle
% the visibility of each curve.
set(M.ui.tsCheck,   'UserData', hTSeries);
set(M.ui.predCheck, 'UserData', hFit);
set(M.ui.resCheck,  'UserData', hResidual);
if get(M.ui.tsCheck,   'Value')==0, set(hTSeries,  'Visible', 'off'); end
if get(M.ui.predCheck, 'Value')==0, set(hFit,      'Visible', 'off'); end
if get(M.ui.resCheck,  'Value')==0, set(hResidual, 'Visible', 'off'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (2) set text fields
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% on Windows, we can save space expressing the "degrees" units string
if ispc, degStr = char(176);
else     degStr = 'deg'; end

% print out percent variance (R^2) explained: This calculation is incorrect
% because the predictions contain the trends.
% (ras) so how do we make it correct?
txt = sprintf('Variance explained: %.2f%%', variance_explained*100);
set(M.ui.r2Text, 'String', txt);

% print out coords of voxel
txt = ['Coords: [' num2str(M.params.roi.coords(:,voxel)') ']'];
set(M.ui.coordsText, 'String', txt);

txt = sprintf('x=%.1f, y=%.1f (%s)',rfParams(1,1), rfParams(1,2), degStr);
set(M.ui.xyText, 'String', txt);

[th r] = cart2pol(rfParams(1,1), rfParams(1,2));
th = round( (pi/2 - th) * (180/pi) ); % rad -> deg CW from up
txt = sprintf('r=%.1f %s, theta=%i %s', r, degStr, th, degStr);
set(M.ui.rthText, 'String', txt);

if size(rfParams,1)<5 && rfParams(3)==rfParams(5)
	% sigma major and minor are the same: circular pRF
	% (on the rare instance where an oval Gaussian was fit with exactly the
	% same major and minor axes, I think this should still be sufficiently
	% clear)
	txt = sprintf('sigma=%.1f%s', rfParams(3), degStr);
elseif size(rfParams, 1)==2
    txt = sprintf('sigma_1=%.2f, sigma_2=%.2f %s', rfParams(1,3), rfParams(2,3), degStr);
else
	% oval Gaussian: report both major and minor sigma
	txt = sprintf('sigma_1=%.1f, sigma_2=%.1f %s', rfParams(3), rfParams(5), degStr);
end
set(M.ui.sigmaText, 'String', txt);

txt = sprintf('beta=%.1f %% %s/sec', rfParams(4), degStr);
set(M.ui.betaText, 'String', txt);

% from rmPlot: cache the data for the current time series plot as well
M.currTsData.x    = M.x;
M.currTsData.pred = pred;
M.currTsData.raw  = M.tSeries(:,voxel);
M.currTsData.sepx = M.sepx;
M.currTsData.pRF  = RF;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (3) show receptive field
%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(M.ui.peakCheck, 'UserData', hResidual);
rescaleToPeak = isequal( get(M.ui.peakCheck, 'Checked'), 'on' );
if get(M.ui.overlayStimCheck, 'Value')==0
    if min(M.params.analysis.Y) == max(M.params.analysis.Y)
        subplot(M.ui.rfAxes);
        plot(M.params.analysis.X,RF,'bo-','MarkerFaceColor','b');
        axis square

    else % default

        % show the RF only, don't overlay the stimulus image
        if rescaleToPeak==0,
            peak = [];
        else
            TR = M.params.stim(1).framePeriod;  % TODO: deal w/ diff't TRs in diff't scans
            peak = rfParams(4) .* M.params.analysis.sampleRate^2 .*  TR;
        end;

        rfPlot(M.params, RF, M.ui.rfAxes,peak);

        % jw add for exponentiated sigma
        modelName = rmGet(M.model{1}, 'desc');
        if strcmpi(modelName, 'fitprf')
            rfParamsExp    = rfParams;
            rfParamsExp(3) = rfParams(5); % put the exponentiated sigma into the sigma slot
            RFsEXP = rmPlotGUI_makeRFs(modelName, rfParamsExp, M.params.analysis.X, M.params.analysis.Y);
            rfPlot(M.params, RFsEXP, M.ui.rfAxes,peak);
            M.currTsData.pRF = RFsEXP; %store the exponentiated pRF in the User Data
        end
    end

    % make sure the time slider is not visible
    mrvSliderSet(M.ui.time, 'Visible', 'off');

else
	%% overlay the stimulus image on the RF
	TR = M.params.stim(1).framePeriod;  % TODO: deal w/ diff't TRs in diff't scans

    % Our time series starts at t = 0, not t = TR, hence we need to
    % subtract 1 from the x-axis index before multiplying by the TR to get
    % the correct time
    %t = get(M.ui.time.sliderHandle, 'Value') * TR;
    t = (get(M.ui.time.sliderHandle, 'Value') -1 )* TR;

	% show the time point on the time series axes
	axes(M.ui.tsAxes);
	delete( findobj('Parent', gca, 'Tag', 'CurStimPoint') );
	AX = axis;
	hLine = line([t t], AX(3:4), 'LineWidth', 2, 'LineStyle', '--', 'Color', 'm');
	set(hLine, 'Tag', 'CurStimPoint');

	% get the stimulus image for this frame
    %   becausae the t-series starts at t = 0, and the frame number starts at
    %   f = 1, the forumala to derive the frame from the time is
    %       f = (t+tr)/tr;
    %   and not
    %       f = t/tr;
    tr = M.params.stim(1).framePeriod; % MR frame
	f = round((t+tr)/ tr);
	if f < 1 || f > size(M.params.analysis.allstimimages, 1)
		fprintf('Can''t visualize stimulus: selected time point outside scan range.\n');
		return
	end
	[stimImage RF] = getCurStimImage(M, f, RF);
    if max(stimImage(:)) > 1 && max(stimImage(:)) < 256
        stimImage = stimImage/255;
    end
    if min(stimImage(:))<0
        stimImage = max(stimImage,0); % clip
    end

	% overlay and display
	axes(M.ui.rfAxes); cla
	RF_img(:,:,1) = stimImage;
	RF_img(:,:,2) = RF;
	RF_img(:,:,3) = RF;

    [x,y] = prfSamplingGrid(M.params);
    x = x(:); y = y(:);
    imagesc(x, -y, RF_img); hold on;
    plot([min(x) max(x)], [0 0], 'k-');
    plot( [0 0], [min(y) max(y)], 'k-');

    axis image  xy;
    ylabel('y (deg)');
    xlabel('x (deg)');

    % make sure the time slider is visible
	mrvSliderSet(M.ui.time, 'Visible', 'on');
end

% turn on grid if the menu option is selected
if isequal( get(M.ui.grid, 'Checked'), 'on' ),	grid on
else                                            grid off; end

%% store the updated M struct in the figure's user data
M.prevVoxel = voxel;  % also update the most recent voxel
set(M.fig, 'UserData', M);

return;
%--------------------------------------

%--------------------------------------
function [stimImage RF] = getCurStimImage(M, f, RFvals)
% Get a stimulus image matching the sampling positions as the RF.
% Also returns the RF resampled into a square grid.
x = prfSamplingGrid(M.params);

% account for the different stimuli that are shown next to each other
% f originally refers to the frame in the combined time series across scans:
% we want to break this down into scan n, frame f within that scan.
n = 1;
nStimScans = numel(M.params.stim);
while n <= nStimScans,
    tmp = f + M.params.stim(n).prescanDuration;
    if tmp > size(M.params.stim(n).images_org,2),
        f = tmp - size(M.params.stim(n).images_org,2);
        n = n + 1;
    else
        f = tmp;
        break;
    end
end

% stim image
stimImage     = NaN(size(x));
stimImage(M.params.stim(1).instimwindow) = M.params.stim(n).images_org(:,f);
stimImage     = reshape(stimImage, size(x));

% RF
RF     = NaN(size(x));
RF(M.params.stim(1).instimwindow) = normalize(RFvals, 0, 1);
RF     = reshape(RF, size(x));

return
% -------------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [header_mat,data_mat] = mhdrload(file)
%MHDRLOAD Load data from an ASCII file containing multiple text
% headers throughout the file.
% [header, data] = MHDRLOAD('filename.ext') reads a data file
% called 'filename.ext', which contains a text header.  There
% is no default extension; any extensions must be explicitly
% supplied.
%
% The first output, HEADER, is the header information, returned
% as a text array.
% The second output, DATA, is the data matrix.  This data matrix
% has the same dimensions as the data in the file, one row per
% line of ASCII data in the file.  If the data is not regularly
% spaced (i.e., each line of ASCII data does not contain the
% same number of points), the data is returned as a column
% vector.
%
% Limitations:  No lines of the text header can begin with
% a number.  The header must come before the data.
%
% MODIFIED from hdrload.m: Dec 20, 2002  Jeff Daniels, NSWCCD - ARD
% UPDATED September 20, 2006, J. Daniels
%
% See also LOAD, SAVE, SPCONVERT, FSCANF, FPRINTF, STR2MAT, HDRLOAD.
% See also the IOFUN directory.
%
% EXAMPLE:
%   If example_data.txt is:
%   Recorded Data: 12/15/2001
%   header 1
%   rows = 2 cols = 2
%   12 23
%   34 21
%   header 2
%   rows = 3 cols = 3
%   19 73 13
%   33 32 47
%   34 12 68
%
%   MHDRLOAD returns:
% 	header(:,:,1) =
%
% 	Recorded Data: 12/15/2001
% 	header 1
% 	rows = 2 cols = 2
%
% 	header(:,:,2) =
%
% 	header 2
% 	rows = 3 cols = 3
%
% 	data(:,:,1) =
%
%        12    23     0
%        34    21     0
%         0     0     0
%
% 	data(:,:,2) =
%
%        19    73    13
%        33    32    47
%        34    12    68

% check number and type of arguments
if nargin < 1
    error('Function requires one input argument');
elseif ~ischar(file)
    error('Input argument must be a string representing a filename');
end

% Open the file.  If this returns a -1, we did not open the file
% successfully.
fid = fopen(file);
if fid==-1
    error('File not found or permission denied.');
end

% Initialize loop variables
% We store the number of lines in the header, and the maximum length
% of any one line in the header.  These are used later in assigning
% the 'header' output variable.
no_lines = 0;
max_line = 0;

% We also store the number of columns in the data we read.  This way
% we can compute the size of the output based on the number of
% columns and the total number of data points.
ncols = 0;

% Finally, we initialize the data to [].
data = [];
j=1;

% Start processing.
line = fgetl(fid);
if ~ischar(line)
    disp('Warning: file contains no header and no data')
end;

while line~=-1
    [data, ncols, errmsg, nxtindex] = sscanf(line, '%f');
    % One slight problem, pointed out by Peter VanderWal: If the first
    % character of the line is 'e', then this will scan as 0.00e+00.
    % We can trap this case specifically by using the 'next index'
    % output: in the case of a stripped 'e' the next index is one,
    % indicating zero characters read.  See the help entry for 'sscanf'
    % for more information on this output parameter.
    % We loop through the file one line at a time until we find some
    % data.  After that point we stop checking for header information.
    % This part of the program takes most of the processing time, because
    % fgetl is relatively slow (compared to fscanf, which we will use
    % later).
    while isempty(data) | nxtindex==1
		no_lines = no_lines+1;
		max_line = max([max_line, length(line)]);
		% Create unique variable to hold this line of text information.
		% Store the last-read line in this variable.
		eval(['line', num2str(no_lines), '=line;'])
		line = fgetl(fid);
		if ~ischar(line)
			disp('Warning: file contains no data')
		break
		end;
		[data, ncols, errmsg, nxtindex] = sscanf(line, '%f');
    end % while

    % Create header output from line information. The number of lines and
	% the maximum line length are stored explicitly, and each line is
	% stored in a unique variable using the 'eval' statement within the
	% loop. Note that, if we knew a priori that the headers were 10 lines
	% or less, we could use the STR2MAT function and save some work.
	% First, initialize the header to an array of spaces.
	header_mat(1:no_lines,1:max_line,j) = char(' '*ones(no_lines, max_line));
	for i = 1:no_lines
		if length(eval(['line' num2str(i)]))
          		varname = eval(['line' num2str(i)]);
        	else
            		varname = ['line',num2str(i)];
        	end
		% Note that we only assign this line variable to a subset of this
		% row of the header array.  We thus ensure that the matrix sizes in
		% the assignment are equal.
		%eval(['header(i, 1:length(' varname ')) = ' varname ';']);
		header_mat(i,1:length(varname),j)=varname;
	end

    % Now that we have read in the first line of data, we can skip the
    % processing that stores header information, and just read in the
    % rest of the data.
    data = [data; fscanf(fid, '%f')];
    % Resize output data, based on the number of columns (as returned
    % from the sscanf of the first line of data) and the total number of
    % data elements. Since the data was read in row-wise, and MATLAB
    % stores data in columnwise format, we have to reverse the size
    % arguments and then transpose the data.  If we read in irregularly
    % spaced data, then the division we are about to do will not work.
    % Therefore, we will trap the error with an EVAL call; if the reshape
    % fails, we will just return the data as is.
    eval('data_mat(1:length(data)/ncols,1:ncols,j) = reshape(data, ncols, length(data)/ncols)'';', '')
    line = fgetl(fid);
    j=j+1;
    no_lines = 0;
end

fclose(fid);
% And we're done!


text:
function out= cell_flaten(C, dim)
%CELL_FLATEN - Flaten all layers of a cell array
%
%Synopsis:
% CFLAT= cell_flaten(C)
% CFLAT= cell_flaten(C, DIM)
%
%Arguements:
%  C: Arbitrary cell array
%  DIM: Dimension in which CFLAT should live, default 2.
%
%Returns:
% CFLAT: Cell array, where each entry in CFLAT is a non-cell object.

% blanker@cs.tu-berlin.de

if nargin<2,
  dim= 2;
end

if ~iscell(C),
  out= {C};
  return;
end

out= {};
for ci= 1:length(C),
  out= cat(dim, out, cell_flaten(C{ci}));
end

author: bbci


text:
function [X, meta] = xifread(filename, debug)

%xifread version 0.3 - Reads HDILab/ATL xif ultrasound data files.
%Written by Bo Lind - bolind@ikol.dk
%
%Usage: [X, meta] = xifread('myfile.xif');
%or:    [X, meta] = xifread('myfile.xif', 1); for verbose mode
%
%Where meta is an (optional) struct with the following fields:
%
%scanlines   Number of scanlines
%samples     Number of samples per scanline
%frames      Number of frames
%framerate   Frame rate
%width       Width in mm
%height      Height in mm
%
%Note, X is a three-dimensional matrix of type uint8, to save memory.
%Convert to double if needed.

%Check number of arguments. One is required, the filename, if number two
%evaluates to true, verbose/debug mode info will be printed.
if (nargin < 1)
    error('xifread: Bad number of arguments.');
elseif (nargin == 1)
    debug = 0;
elseif (nargin > 1)
    debug = 1;
end

%Open file.
FID = fopen(filename, 'r');

%Abort if error opening file (does not exist, etc.)
if (FID == -1)
    error('Error opening file: %s, aborting...', filename);
end

%Run through the XML header, reading in scanlines, samples etc.
s = [];
while (size(strfind(s, '#SOH')) == 0)
    s = fgets(FID);
        if (strfind(s, 'echoNumLines'))
        s = fgets(FID);
        s = fgets(FID);
        meta.scanlines = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'echoNumDisplaySamples'))
        s = fgets(FID);
        s = fgets(FID);
        meta.samples = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'numFrames'))
        s = fgets(FID);
        s = fgets(FID);
        meta.frames = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'frameRate'))
        s = fgets(FID);
        s = fgets(FID);
        meta.framerate = str2num(cell2mat(regexp(s, '(\d+)', 'match')));
    end
    if (strfind(s, 'echoScan_linearWidth'))
        s = fgets(FID);
        s = fgets(FID);
        %Different regexp, as we need to match a decimal number.
        meta.width = str2num(cell2mat(regexp(s, '(\d+\.\d+)', 'match')));
    end
    if (strfind(s, 'echoDepth_twodDepth'))
        s = fgets(FID);
        s = fgets(FID);
        %Same as above.
        meta.height = str2num(cell2mat(regexp(s, '(\d+\.\d+)', 'match')));
    end
end

if (debug)
    fprintf('numLines: %d\n', meta.scanlines);
    fprintf('numDisplaySamples: %d\n', meta.samples);
    fprintf('numFrames: %d\n', meta.frames);
    fprintf('frameRate: %d\n', meta.framerate);
    fprintf('echoScan_linearWidth: %d\n', meta.width);
    fprintf('echoDepth_twodDepth: %d\n', meta.height);
end

%Read throuh C-style header.
s = [];
while (size(strfind(s, '#EOH')) == 0)
    s = fgets(FID);
end

%Now filepointer is at start of data. We simply read them in.
%Initializing X first, to allocate memory in a sane fashion.
X = uint8(zeros(uint16([meta.samples, meta.scanlines, meta.frames])));
for i = 1:meta.frames
    X(:,:,i) = fread(FID, [meta.samples, meta.scanlines], 'uchar');
end

fclose(FID);


text:
function x = prox_sum_square(v, lambda)
% PROX_SUM_SQUARE    Proximal operator of sum-of-squares.
%
%   prox_sum_square(v,lambda) is the proximal operator of
%   (1/2)||.||_2^2 with parameter lambda.

    x = (1/(1 + lambda))*v;
end

author: cvxgrp

repo: proximal

save_path: github-repos/MATLAB/cvxgrp-proximal


text:
function init_axes2d(axes_handle, grid2d, withinterp, withpml, isswapped)
chkarg(~isempty(axes_handle) && ishandle(axes_handle), '"axes_handle" should be handle.');
chkarg(istypesizeof(grid2d, 'Grid2d'), '"grid2d" should be instance of Grid2d.');
chkarg(istypesizeof(withinterp, 'logical'), '"withinterp" should be logical.');
chkarg(istypesizeof(withpml, 'logical'), '"withpml" should be logical.');
chkarg(istypesizeof(isswapped, 'logical'), '"isswapped" should be logical.');

hold(axes_handle, 'on');

if withinterp
	lplot = grid2d.lplot(GT.prim, withinterp, withpml);
else
	lplot = grid2d.lvoxelbound(GT.prim, withpml);
end

if ~isswapped
	h = Dir.h; v = Dir.v;
else
	h = Dir.v; v = Dir.h;
end

axis(axes_handle, [lplot{h}([1 end]), lplot{v}([1 end])]);
daspect(axes_handle, [1 1 1]);  % axis(axes_handle, 'image');
box(axes_handle, 'on');

hr = rotate3d;
setAllowAxesRotate(hr, axes_handle, false);

% set(axes_handle, 'FontSize', 18, 'FontWeight', 'Bold');
str = [char(grid2d.axis(h)), '  (', char(hex2dec('00D7'))];
if grid2d.unitvalue < 1e5 && grid2d.unitvalue > 1e-3
	str = [str, num2str(grid2d.unitvalue), ')'];
else
	str = [str, num2str(grid2d.unitvalue, '%.2e'), ')'];
end

xlabel(axes_handle, str, 'Rotation', 0, 'FontSize', 15);
% s = [char(grid2d.axis(v)), '  (', char(hex2dec('00D7')), num2str(grid2d.unitvalue, '%e'), ')'];
str = char(grid2d.axis(v));
ylabel(axes_handle, str, 'Rotation', 0, 'FontSize', 15);

set(axes_handle, 'TickDir', 'out');

% L0 = scalar3d.gi.L0;
% dlu = scalar3d.gi.display_length_unit;
% dlu_scale = scalar3d.gi.dlu_scale;
% xlabel(strcat(AxisName(Xx), ' (', char(hex2dec('00D7')), ...  % 00D7 is a unicode character for the multiplication sign.
% 	num2str(L0*dlu_scale), dlu, ')'));
% ylabel(strcat(AxisName(Yy), ' (', char(hex2dec('00D7')), ...  % 00D7 is a unicode character for the multiplication sign.
% 	num2str(L0*dlu_scale), dlu, ')'));
% zlabel(strcat(AxisName(Zz), ' (', char(hex2dec('00D7')), ...  % 00D7 is a unicode character for the multiplication sign.
% 	num2str(L0*dlu_scale), dlu, ')'));
% title(plot_title);

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [U] = spm_voice_dct(N,K,n,s)
% Logarithmically sampled discrete cosine transform matrix
% FORMAT [U] = spm_voice_dct(N,K,n,[s])
%
%  N         - dimension
%  K         - order
%  n         - log scaling parameter (typically 4)
%  s         - optional linear scaling [default 1]
%
% This routine returns a discrete cosine transform matrix sampled
% logarithmically according to a scaling parameter.
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_dct.m 7750 2019-12-05 17:54:29Z spm $

% defaults
%--------------------------------------------------------------------------
if nargin < 4, s = 1; end

% discrete cosine transform matrix
%--------------------------------------------------------------------------
nu = exp(-n*(0:(N - 1))/N);                  % log spacing
nu = 1 - nu;                                 % ensure range lies between
nu = s*N*nu/nu(end);                         % 0 and s*N
U  = spm_dctmtx(N,K,nu);                     % DCT matrix

author: spm


text:
function pass = test_isPeriodicTech(pref)
% Test isPeriodicTech().

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1000*pref.cheb2Prefs.chebfun2eps;

f = chebfun2v(@(x,y) x, @(x,y) y);
pass(1) = ~isPeriodicTech(f);

f1 = chebfun2(@(x,y) cos(pi*x), [ -1, 1, -pi, pi ], 'trig');
f2 = chebfun2(@(x,y) sin(y), [ -1, 1, -pi, pi ], 'trig');
f = [f1; f2];
pass(2) = isPeriodicTech(f);
f = [f1; f1; f2];
pass(3) = isPeriodicTech(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef CommunityDetection < PROBLEM
% <single> <label> <large/none> <expensive/none>
% The community detection problem with label based encoding
% dataNo --- 1 --- Number of dataset

%------------------------------- Reference --------------------------------
% Y. Tian, S. Yang, and X. Zhang, An evolutionary multiobjective
% optimization based fuzzy method for overlapping community detection, IEEE
% Transactions on Fuzzy Systems, 2020, 28(11): 2841-2855.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% The datasets are taken from
% http://www-personal.umich.edu/~mejn/netdata/
% No.   Name        Nodes   Edges
% 1     Karate      34      78
% 2     Dolphin     62      159
% 3     Polbook     105     441
% 4     Football    115     613

    properties(Access = private)
        Adj;    % Adjacency matrix of the network
        G;      % The graph object
    end
    methods
    	%% Default settings of the problem
        function Setting(obj)
            % Load data
            dataNo    = obj.ParameterSet(1);
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'Dataset_CD.mat'),'Dataset');
            str     = {'Karate','Dolphin','Polbook','Football'};
            obj.Adj = Dataset.(str{dataNo});
            obj.G   = graph(obj.Adj);
            % Parameter setting
            obj.M = 1;
            obj.D = size(obj.Adj,2);
            obj.lower    = 1     + zeros(1,obj.D);
            obj.upper    = obj.D + zeros(1,obj.D);
            obj.encoding = 3     + zeros(1,obj.D);
        end
        %% Repair invalid solutions
        function PopDec = CalDec(obj,PopDec)
            for i = 1 : size(PopDec,1)
                P = zeros(1,obj.D);
                while ~all(P)
                    x = find(~P,1);
                    P(PopDec(i,:)==PopDec(i,x)) = max(P) + 1;
                end
                PopDec(i,:) = P;
            end
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = zeros(size(PopDec,1),1);
            for i = 1 : size(PopObj,1)
                PopObj(i) = 1 - Modularity(obj.Adj,Decoding(PopDec(i,:)));
            end
        end
        %% Display a population in the decision space
        function DrawDec(obj,Population)
            [~,best] = min(Population.objs);
            C = Decoding(Population(best).dec);
            h = plot(Draw([]),obj.G,'MarkerSize',6,'EdgeColor',[.3 .3 .3]);
            tempStream = RandStream('mlfg6331_64','Seed',2);
            for i = 1 : length(C)
                highlight(h,C{i},'NodeColor',rand(tempStream,1,3));
            end
        end
    end
end

function Community = Decoding(Dec)
    Community = {};
    while any(Dec)
        current      = find(Dec==Dec(find(Dec,1)));
        Community    = [Community,current];
        Dec(current) = 0;
    end
end

function Q = Modularity(Adj,C)
    Q = 0;
    M = sum(sum(Adj))/2;
    for i = 1 : length(C)
        Q = Q + sum(sum(Adj(C{i},C{i})))/2/M - (sum(sum(Adj(C{i},:)))/2/M)^2;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [correct_class_decision_values normalized_rank_results rank_confusion_matrix] = get_rank_and_decision_value_results(cv, YTe, classifier_labels, decision_values, create_rank_confusion_matrix)
% This helper method calculates the normalized rank results, the decision values for the correct class, and the rank confusion matrix information.
%  This method should be run every time the classifier is tested.
%
%  Input parameters:
%    YTe:  the labels for the current test points
%    decision_values:  The decision values for the current test point
%    create_rank_confusion_matrix:  whether a confusion matrix should be created
%
%  The returned values are in the structure NORMALIZED_RANK_AND_DECISION_VALUE_RESULTS and have the fields
%
%   correct_class_decision_values: a [num_test_points x 1] vector containing the decision value for the ith test point (for the correct class of the ith point)
%
%   normalized_rank_results:  a [num_test_points x 1] vector containing the normalized rank value for the ith test point
%
%   rank_confusion_matrix: [num_predicted_classes x num_actual classes] current rank confusion matrix.
%      the i, j entry of this matrix tells how high up on the rank of predictions for the jth class
%      was entry i.
%

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.edu)
%
%     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/>.

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

rank_confusion_matrix = [];

YTe_unique_values = unique(YTe);

% % %[vals sorted_inds] = sort(decision_values', 'descend');   % NEED TO BE CAREFUL THERE IS NOT TIE IN THE MAX DECISION VALUE (could create a bias in results)
[vals sorted_inds] = sort((decision_values + eps .* rand(size(decision_values)))', 'descend');   % adding a small amount of noise so that there is not a tie in the decision values
the_ranks_all_test_points = sorted_inds';   % each row has all the ranked order results for a given test point (i.e., all ranks for YTe(i))

% go through all test points and find the ranking and the decision value of the real label
for iTestPoint = 1:length(YTe)
    curr_rank_results(iTestPoint) = find(YTe_unique_values(the_ranks_all_test_points(iTestPoint, :)) == YTe(iTestPoint));  % find rank of the real label YTe
    correct_class_decision_values(iTestPoint) = decision_values(iTestPoint, (classifier_labels == YTe(iTestPoint)));   % find the decision value for real label
end

normalized_rank_results = 1 - (((curr_rank_results) - 1)./(length(YTe_unique_values) - 1));

% get information to create the rank confusion matrix
if create_rank_confusion_matrix == 1
    for iUniqueYTe = 1:length(YTe_unique_values)   % assuming that each time unique(YTe) has all possible label values (if it doesn't then this code will not work)

            curr_rank_inds = find(YTe == YTe_unique_values(iUniqueYTe));
            [rank_cm_vals rank_cm_inds] = sort(the_ranks_all_test_points(curr_rank_inds, :)');   % rank_cm_inds contains predicted classes from the most likely (first index value), to least likely (last index value)

            rank_confusion_matrix(:, iUniqueYTe) = sum(1 -((rank_cm_inds -1)./(length(YTe_unique_values) - 1)), 2);  % add all test points together for a total rank confusion matrix

    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ out ] = proj_Euclidean_ball( x,c,r )
%PROJ_EUCLIDEAN_BALL computes the orthogonal projection onto the Euclidean ball {x:||x-c||<=r}
%
%  Usage:
%  out = PROJ_EUCLIDEAN_BALL(x,[c],[r])
%  ===========================================
%  Input:
%  x - point to be projected (vector/matrix)
%  c - center of the ball (vector/matrix) [default: c=0]
%  r - positive radius (a positive scalar) [default: r=1]
%  ===========================================
%  Assumptions:
%  l2 or Frobenius norm
%  ===========================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

%reading the user x and setting defalut values when required.
if (nargin < 1)
     error ('usage: proj_Euclidean_ball(input,[c],[r])') ;
end

if (nargin < 3)
    %setting default value of r to 1
    r = 1;
end

if ((nargin < 2) || (isempty (c)))
    %setting default value of c to zeros
    c  = zeros(size(x)) ;
end

out=  c+r *(x-c) / max(norm(x-c,'fro'),r) ;
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [w, momentum] = sgd(w, momentum, grad, opts, lr)
%SGD
%   Example SGD solver, with momentum, for use with CNN_TRAIN and
%   CNN_TRAIN_DAG.
%
%   The convergence of SGD depends heavily on the learning rate (set in the
%   options for CNN_TRAIN and CNN_TRAIN_DAG).
%
%   If called without any input argument, returns the default options
%   structure.
%
%   Solver options: (opts.train.solverOpts)
%
%   `momentum`:: 0.9
%      Parameter for Momentum SGD; set to 0 for standard SGD.
%
%   Note: for backwards compatibility, the parameter can also be set in
%   opts.train.momentum.

% Copyright (C) 2016 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin == 0 % Return the default solver options
  w = struct('momentum', 0.9);
  return;
end
if isempty(momentum)
  momentum = 0 ;
end

momentum = opts.momentum * momentum - grad ;
w = w + lr * momentum ;

repo: CREST-Release


text:
% This script reads word2vec bin file

addpath('./hglmm_fv_v1.6/utilities');
addpath('./hglmm_fv_v1.6/fv');
f_name = 'GoogleNews-vectors-negative300.bin';

% some settings
sqrt_flag = false;
norm_flag = false;

[w_names, w_features] = read_word2vec_binfile(f_name, sqrt_flag, norm_flag);

save('GoogleNews_words.mat','w_names');
save('GoogleNews_vectors.mat','w_features','-v7.3');

author: layumi

repo: Image-Text-Embedding


text:
function st = makecall(args, origin, structargs)

    % check if signed in and grab username, key, domain
    [un, key, domain] = signin;
    if isempty(un) || isempty(key)
        error('Plotly:CredentialsNotFound',...
             ['It looks like you haven''t set up your plotly '...
              'account credentials yet.\nTo get started, save your '...
              'plotly username and API key by calling:\n'...
              '>>> saveplotlycredentials(username, api_key)\n\n'...
              'For more help, see https://plot.ly/MATLAB or contact '...
              'chris@plot.ly.']);
    end

    platform = 'MATLAB';

    args = m2json(args);
    kwargs = m2json(structargs);
    url = [domain '/clientresp'];
    payload = {'platform', platform, 'version', plotly_version, 'args', args, 'un', un, 'key', key, 'origin', origin, 'kwargs', kwargs};

    if (is_octave)
        % use octave super_powers
        resp = urlread(url, 'post', payload);
    else
        % do it matlab way
        resp = urlread(url, 'Post', payload);
    end

    st = loadjson(resp);

    response_handler(resp);

end

author: plotly


text:
function varargout = drawGraphEdges(varargin)
%DRAWGRAPHEDGES Draw edges of a graph.
%
%   drawGraphEdges(NODES, EDGES)
%   Draws a graph specified by a set of nodes (array N-by-2 or N-by-3,
%   corresponding to coordinate of each node), and a set of edges (an array
%   Ne-by-2, containing to the first and the second node of each edge).
%
%   drawGraphEdges(..., SEDGES)
%   Specifies the draw mode for each element, as in the classical 'plot'
%   function.
%   Default drawing is a blue line for edges.
%
%
%   H = drawGraphEdges(...)
%   Returns handle to the set of edges.
%
%   See also
%     graphs, drawGraph, fillGraphFaces
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-11-24
% Copyright 2005-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

%% Input argument processing

% initialisations
e = [];

% check input arguments number
if nargin == 0
    help drawGraphEdges;
    return;
end

% extract handle of axis to draw on
if isAxisHandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% First extract the graph structure
var = varargin{1};
if iscell(var)
    % TODO: should consider array of graph structures.
    % graph is stored as a cell array: first cell is nodes, second one is
    % edges, and third one is faces
    n = var{1};
    if length(var) > 1
        e = var{2};
    end
    varargin(1) = [];

elseif isstruct(var)
    % graph is stored as a structure, with fields 'nodes', 'edges'
    n = var.nodes;
    e = var.edges;
    varargin(1) = [];

else
    % graph is stored as set of variables: nodes + edges
    n = varargin{1};
    e = varargin{2};
    varargin(1:2) = [];
end

% check if there are edges to draw
if size(e, 1) == 0
    return;
end

% setup default drawing style if not specified
if isempty(varargin)
    varargin = {'-b'};
end

%% main drawing processing

if size(n, 2) == 2
    % Draw 2D edges
    x = [n(e(:,1), 1) n(e(:,2), 1)]';
    y = [n(e(:,1), 2) n(e(:,2), 2)]';
    he = plot(ax, x, y, varargin{:});

elseif size(n, 2) == 3
    % Draw 3D edges
    x = [n(e(:,1), 1) n(e(:,2), 1)]';
    y = [n(e(:,1), 2) n(e(:,2), 2)]';
    z = [n(e(:,1), 3) n(e(:,2), 3)]';
    he = plot3(ax, x, y, z, varargin{:});

end

%% format output arguments

if nargout == 1
    varargout = {he};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% DEMCMU35VARGPLVMANIMATE Load results for the dyn. GPLVM on CMU35 data and
% produce animation and plots. The plots (and errors) are in the data
% space. For plots / errors  in the scaled space (for comparison with FGPLVM)
% see demCmu35VargplvmPlotsScaled.m
% DESC Load results for the dyn. GPLVM on CMU35 data and produce animation
% and plots.
% COPYRIGHT :  Andreas C. Damianou, Michalis K. Titsias, 2011
%
% SEEALSO : demCmu35VargplvmLoadChannels, demCmu35gplvmVargplvm3.m,
% demCmu35VargplvmPlotsScaled.m
% VARGPLVM

randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'cmu35gplvm';
if ~exist('experimentNo') experimentNo = 33; end
% Options for the following: 'Legs', 'Body'
if ~exist('predictPart') predictPart = 'Legs'; end
% Sampling:
if ~exist('doSampling') doSampling = 0; end
if ~exist('showSkel') showSkel = 1; end
if ~exist('displayPlots')  displayPlots = 1; end
% -1 is a flag meaning that plotRange == missingInd
if ~exist('plotRange') plotRange = -1; end

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);
dataSetName = 'cmu35gplvm';

capName = dataSetName;
capName(1) = upper(capName(1));
fileName=['dem' capName 'Vargplvm' num2str(experimentNo)];

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);
origBias = mean(Y);
origScale = 1./sqrt(var(Y));

legInd = [8:14];
bodyInd = [21:50];

startInd = 63;
dt=0.05;
timeStampsTest = ([0:size(Ytest,1)-1].*dt)';

load(fileName)
model.dynamics.t_star = timeStampsTest;

% Choose one of the following lines as appropriate to obtain the results for the leg or body reconstruction
% respectively.
if strcmp(predictPart,'Legs')
    load(['demCmu35Vargplvm' num2str(experimentNo) 'PredLegs.mat']); missingInd = legInd;
elseif strcmp(predictPart, 'Body')
    load(['demCmu35Vargplvm' num2str(experimentNo) 'PredBody.mat']) ;missingInd = bodyInd;
end

YtestOrig = Ytest;
YtestGplvm = Ytest;
YtestGplvm(startInd:end, missingInd) = NaN;
indexMissingData = startInd:size(YtestGplvm);
% In case only barmu and lambda are loaded from the results then the next step is necessary to find
% the predictive dencity Ypred.
[x, varx, modelUpdated] = vargplvmDynamicsUpdateModelTestVar(model, barmu, lambda, YtestGplvm);
Testmeans = x(indexMissingData, :);
Testcovars = varx(indexMissingData, :);
[mu, sigma] = vargplvmPosteriorMeanVar(modelUpdated, Testmeans, Testcovars);
Ypred = mu;

%%
origYtest = Ytest;
err = origYtest(startInd:end, missingInd) - Ypred(:, missingInd);
err = err*180/pi;
angleErrorGplvm = sqrt(mean(mean((err.*err))))
load cmu35TaylorScaleBias
err = err.*repmat(scale(1, missingInd+9), size(origYtest, 1)-startInd+1, 1);
taylorErrorGplvm = sum(sum(err.*err))/length(missingInd)

%%
% NN
temp = ones(1, size(Ytest, 2));
temp(missingInd) = 0;
presentInd = find(temp);
YtrainNn = Y(:, presentInd);
YtestNn = origYtest(startInd:end, presentInd);
dists = dist2(YtrainNn, YtestNn);
[void, bestIndNn] = min(dists);

%%

if displayPlots
    colordef white
    if plotRange == -1 % flag
        plotRange = missingInd;
    end
    for plotNo = plotRange
        figNo = plotNo - min(plotRange) + 1;
        figure(figNo)
        clf
        lin = plot(1:size(origYtest, 1), origYtest(:, plotNo), '-')
        hold on
        lin = [lin; plot(1:size(origYtest, 1), [origYtest(1:startInd-1, plotNo); Y(bestIndNn, plotNo)], ':')];
        lin = [lin; plot(1:size(origYtest, 1), [origYtest(1:startInd-1, plotNo); Ypred(:, plotNo)], '--')];
        ax = gca;
        xlabel('Frame')
        ylabel('Normalized joint angle')
        set(ax, 'fontname', 'arial');
        set(ax, 'fontsize', 20);
        set(lin, 'lineWidth', 2);
        pos = get(gcf, 'paperposition')
        origpos = pos;
        pos(3) = pos(3)/2;
        pos(4) = pos(4)/2;
        set(gcf, 'paperposition', pos);
        fontsize = get(gca, 'fontsize');
        set(gca, 'fontsize', fontsize/2);
        lineWidth = get(gca, 'lineWidth');
        set(gca, 'lineWidth', lineWidth*2);
        set(gcf, 'paperposition', origpos);
    end
end

%%

if showSkel
    skel = acclaimReadSkel('35.asf');
    [tmpchan, skel] = acclaimLoadChannels('35_01.amc', skel);
    channels1 = demCmu35VargplvmLoadChannels(Ytest,skel);
    channels2 = demCmu35VargplvmLoadChannels(Ypred,skel);
    channels3 = demCmu35VargplvmLoadChannels(Y(bestIndNn, :),skel);
    %skelPlayData(skel, channels, 1/25);

    startInd = 63;
    skelPlayData2(skel, 1/15, channels1(startInd:end,:), channels2, channels3,{'Ytest','YpredGPLVM','YpredNN'});
end

%%

if doSampling
 %   timeStampsTestOrig = timeStampsTest;
 %   timeStampsTest = (timeStampsTest(end):dt:(timeStampsTest(end)+100*dt))';
    model.dynamics.t_star = timeStampsTest;
    seqTest = 1; % WALK
    %seqTest = 17; % RUN

    if seqTest ~= 1
        datFrom = model.dynamics.seq(seqTest-1)+1;
    else
        datFrom = 1;
    end
    datEnd=model.dynamics.seq(seqTest);

    model.dynamics.t_star = model.dynamics.t_star(1:(datEnd+1-datFrom));%%%
    %     clear model
    %     load demCmu35VargplvmOne
    %     %load demCmu35VargplvmOneRun
    %     ySamp, xSamp] = vargpTimeDynamicsSample(model,1);
    %     channels4 = demCmu35VargplvmLoadChannels(ySamp,skel);
    %     skelPlayData(skel, channels4, 1/25);

    [ySamp, xSamp] = vargpTimeDynamicsSampleSeq(model,1,seqTest);
    channels4 = demCmu35VargplvmLoadChannels(ySamp,skel);
    skelPlayData(skel, channels4, 1/10);

    %     [Testmeans2 Testcovars2] = vargplvmPredictSeq(model.dynamics, timeStampsTest,1);
    %     Varmu2 = vargplvmPosteriorMeanVar(model, Testmeans2, Testcovars2);
    %     channels5 = demCmu35VargplvmLoadChannels(Varmu2,skel);
    %     skelPlayData(skel, channels5, 1/25);

end
%%
%lvmVisualise(model, [], 'skelVisualise', 'skelModify', channels, skel);
%lvmVisualise(model, [], 'TEMPskelVisualise', 'skelModify', channels1,
%skel);

author: SheffieldML

repo: vargplvm


text:
function [handle,hs] = ml_plot_value_func(X,f,dims,options)
%ML_PLOT_VALUE_FUNC Plots a value function associated with a
% generative density model, a regressor which provides a likelihood or any
% other function.
%
%
%   input -----------------------------------------------------------------
%
%       o X    : (N x D), original dataset on which f was learned. It is
%                         needed to be able to determine the plotting boundaries.
%
%       o dims : (1 x 2), input dimensions of data X to be used with f.
%
%
%       o f  : function handle, z = f(X)
%
%

%% Process input and set default options

num_samples = 50;
title_name  = 'your title';
handle      = [];
surf_type   = 'surf';
regr_type   = 'LR';
color       = 'k';
points_size = 15;
bFigure     = true;
bColorbar   = false;
Cmap        = 'hot';

if isfield(options,'Cmap'),         Cmap         = options.Cmap;        end
if isfield(options,'color'),        color        = options.color;       end
if isfield(options,'bFigure'),      bFigure      = options.bFigure;     end
if isfield(options,'bColorbar'),    bColorbar    = options.bColorbar;   end
if isfield(options,'points_size'),  points_size  = options.points_size; end
if isfield(options,'title'),        title_name   = options.title;       end
if isfield(options,'surf_type'),    surf_type    = options.surf_type;   end
if isfield(options,'regr_type'),    regr_type    = options.regr_type;   end

%% Get boundary of original data.
X = X(:,dims);
min_x = min(X(:,1));
max_x = max(X(:,1));
min_y = min(X(:,2));
max_y = max(X(:,2));
[X,Y]  = meshgrid(linspace(min_x,max_x,num_samples),linspace(min_y,max_y,num_samples));

%% Evaluate f

switch regr_type
    case 'LR'
        Data = [X(:),Y(:)];
    case 'GMR'
        Data = [X(:),Y(:)]';
end
z      = f(Data);

%% Plot

if bFigure
    handle = figure;
    set(gcf,'color','w');

end

if strcmp(surf_type,'surf')

    hs = surf(X,Y,reshape(z,size(X)));

elseif strcmp(surf_type,'scatter')

    hs = scatter3(X(:),Y(:),z(:),points_size, color,'filled');

elseif strcmp(surf_type,'pcolor')

    hs = pcolor(X,Y,reshape(z,size(X))); shading interp;

end

%% Plot Colorbar
if bColorbar
    colorbar;
end
%% Plot attributes, name, lables, scaling, etc...

if bFigure
    title(title_name, 'Interpreter','Tex','FontName','Times', 'FontWeight','Light','FontSize',15);
    xlabel('$\xi_x$','Interpreter','LaTex','FontName','Times', 'FontWeight','Light','FontSize',15);
    ylabel('$\xi_y$','Interpreter','LaTex','FontName','Times', 'FontWeight','Light','FontSize',15);
    zlabel('$\kappa$','Interpreter','LaTex','FontName','Times', 'FontWeight','Light', 'FontSize',15);
    axis tight
end

colormap (Cmap)

end

author: epfl-lasa


text:
classdef(Abstract) AbstractGradientCalculationMethod < matlab.mixin.SetGet
    %AbstractGradientCalculationMethod Summary of this class goes here
    %   Detailed explanation goes here

    methods
        tf = useBuiltInMethod(obj);

        tf = shouldComputeSparsity(obj)

        computeGradientSparsity(obj, fun, x0, fAtX0, useParallel)

        g = computeGrad(obj, fun, x0, fAtX0, useParallel);

        J = computeJacobian(obj, cFun, x0, cAtX0, useParallel);

        openOptionsDialog(obj)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%
% BEZIER
%
% ger en bezierkurza mellan p1 och p2, med styrpunkter b och c. n stycken
% punkter, inkluderar p1, men inte p2.

function r=bezier(p1,b,c,p2,n)

t=(0:1/n:1-(1/n)/2)';
r=(1-t).^3*p1+3*t.*(1-t).^2*b+3*t.^2.*(1-t)*c+t.^3*p2;

author: Sable

repo: mcbench-benchmarks