text:
normX = AllFeature';
normX = bsxfun(@rdivide, normX, sqrt(sum(normX.^2,2)));
% [normX, PCAmap] = compute_mapping(normX, 'PCA', 319);
% fid = fopen('G:\face-CASIA-WebFace\sim_label_test.txt');
% C = textscan(fid,'%d');
% fclose(fid);
% sim = C{1};
% fid = fopen('G:\face-CASIA-WebFace\test_data1.txt');
% C = textscan(fid, '%s %d');
% fclose(fid);
% [mappedX, mapping] = JointBayesian(normX, C{2});
num = size(normX,1) / 2;

thresh = zeros(num,1);

% normX = AllFeature';
% normX = mappedX;
for i = 1:num
%     thresh(i) = normX(i,:) * mapping.A * normX(i,:)' + normX(num+i,:) * mapping.A * normX(num+i,:)' - 2 * normX(i,:) * mapping.G * normX(num+i,:)';
    thresh(i) = pdist2(normX(i,:), normX(num+i,:));
end;

num = size(AllFeature1,2);
F1 = AllFeature1';
F1 = bsxfun(@rdivide, F1, sqrt(sum(F1.^2,2)));
% F1 = bsxfun(@minus,F1,PCAmap.mean);
% F1 = F1 * PCAmap.M;
F2 = AllFeature2';
F2 = bsxfun(@rdivide, F2, sqrt(sum(F2.^2,2)));
% F2 = bsxfun(@minus,F2,PCAmap.mean);
% F2 = F2 * PCAmap.M;
% F1 = AllFeature1';
% F2 = AllFeature2';
thresh2 = zeros(num,1);
for i = 1:num
%     thresh2(i) = F1(i,:) * mapping.A * F1(i,:)' + F2(i,:) * mapping.A * F2(i,:)' - 2 * F1(i,:) * mapping.G * F2(i,:)';
    thresh2(i) = pdist2(F1(i,:),F2(i,:));
end;

% model = train(double(sim),sparse(thresh),'-s 2');
% predicted_label = predict(double(sim),sparse(thresh),model);
bestc=256;
lfw_label = ones(6000,1);
lfw_label(3001:6000) = 0;
% predicted_label = predict(double(lfw_label),sparse(thresh2),model);
cmd = [' -t 0 -h 0'];
model = svmtrain(double(sim),thresh,cmd);
[class] = svmpredict(lfw_label,thresh2,model);

author: happynear

repo: FaceVerification


text:
function [Y,W,A] = FDICA(X,W,A,F_i)

%First step: Initialization
% step_size=1e-5;                   %the step distance
step_size = 0.1;
% C_r=0.0002;                             % correlation
max_iteration=1000;                   % the max number of iteration

% if(F_i==33)
%     step_size = step_size/10;
% end
%%%%%%%%%%
frame_N = size(X,2);
% if(F_i==16)
%     a =1;
% end
sign_size = 0;
y_f = W*X;
norm = max(abs(y_f),[],2);
% if(norm>10)
%     norm = repmat(norm,1,size(W,2));
%     W = W./norm;
%     y_f = W*X;
% end

for n_i=1:max_iteration
%         y_fun = power(1+exp(-y_r),-1)+j*power(1+exp(-y_i),-1);
%         W_o1 = step_size*(diag(diag(y_fun*y_f'/frame_N))-(y_fun*y_f'/frame_N))*W_o{K_i}+W_o{K_i};
    y_f = W*X;
    y_r = real(y_f);
    y_i = imag(y_f);
    y_fun = tanh(y_r)+1i*tanh(y_i);
    sign = abs(det(eye(size(X,1))-(y_fun*y_f'/frame_N)));
    A(n_i,F_i-1) = sign;
    %%%%%This part trys to adjust the step_size which is unnecessary
    if(sign>100 && sign_size ==0)
        step_size = step_size/10;
        sign_size =1;
    end
    if(sign>1000 && sign_size ==1)
        step_size = step_size/10;
        sign_size =2;
    end
    if(sign>10000 && sign_size ==2)
        step_size = step_size/10;
        sign_size =3;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    W1 = step_size*(eye(size(X,1))-(y_fun*y_f'/frame_N))*W+W;
    W = W1;

    if(sign<1e-4)
        break;
    end
%     A(n_i,F_i-1) = abs(det(eye(size(X,1))-(y_fun*y_f'/frame_N)));
%     A(n_i,K_i) = abs(det(eye(size(s,2))-(y_fun*y_f'/frame_N)));
%     if(abs(A(n_i,K_i))<0.001) break; end;
end
A(end,F_i-1) = n_i;
Y = y_f;
% Y = Project_back(W,y_f);
return;

repo: Voice-Separation-and-Enhancement


text:
function f = tone2freq(T)
% MUSIC.TONE2FREQ converts a musical semitone to a frequency.
%    F = MUSIC.TONE2FREQ(T) converts the musical semitones in T to frequencies.
%
%    Example
%       f = music.tone2freq(0:2);  % returns [261.63  277.19  293.67]
%
%    See also music.tone2interval, music.tone2note, music.freq2tone.

%    Author: E. Johnson
%    Copyright 2010 The MathWorks, Inc.

fC4 = 261.625565300599;  % Middle C (C4) is 261.63 Hz

f = fC4 .* 2 .^ (T / 12);

author: Sable

repo: mcbench-benchmarks


text:
function value = r8vec_length ( dim_num, x )

%*****************************************************************************80
%
%% R8VEC_LENGTH returns the Euclidean length of an R8VEC
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real X(DIM_NUM), the vector.
%
%    Output, real VALUE, the Euclidean length of the vector.
%
  value = sqrt ( sum ( ( x(1:dim_num) ).^2 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y0 = invlink(y,fali)
% inverse link function for locfit.
% y is a vector of raw fitted values.
% fali is the integer [family link] vector from locfit.
% output is the inv. link.

link = fali(2);

switch(link)
    case 3   % identity
        y0 = y;
    case 4   % log
        y0 = exp(y);
    case 5   % logit - should invert carefully!
        y0 = 1 - 1./(1+exp(y));
    case 6   % inverse
        y0 = 1/y;
    case 7   % sqrt
        y0 = y*abs(y);
    case 8   % arcsin
        y0 = sin(y)*sin(y);
    otherwise
        disp('invlink: Unknown link function');
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef MimSaveImage < MimGuiPlugin
    % MimSaveImage. Gui Plugin for exporting the image currently in the
    % visualisation window to a file
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimSaveImage is a Gui Plugin for the MIM Toolkit.
    %     The gui will create a button for the user to run this plugin.
    %     Running this plugin will raise a Save dialog allowing the user to
    %     choose a filename and format, and then save the image currently in the
    %     visualisation panel.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Export Image'
        SelectedText = 'Export Image'
        ToolTip = 'Save the current image to files'
        Category = 'Export'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1

        Icon = 'export_image.png'
        Location = 20
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.SaveBackgroundImage();
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded();
        end

        function is_selected = IsSelected(gui_app)
            is_selected = false;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function  [X,Y,vals,labI]=mp_utm(optn,varargin)
% MP_UTM   Universal Transverse Mercator projection
%           This function should not be used directly; instead it is
%           is accessed by various high-level functions named M_*.

% mp_utm.m, Peter Lemmond (peter@whoi.edu)

% created mp_utm.m 13Aug98 from mp_tmerc.m, v1.2d distribution, by:
%
% Rich Pawlowicz (rich@ocgy.ubc.ca) 2/Apr/1997
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
% Mathematical formulas for the projections and their inverses are taken from
%
%      Snyder, John P., Map Projections used by the US Geological Survey,
%      Geol. Surv. Bull. 1532, 2nd Edition, USGPO, Washington D.C., 1983.
%

% 10/Dec/98 - PL added various ellipsoids.
% 17/May/12 - clarified hemisphere setting

global MAP_PROJECTION MAP_VAR_LIST

% define a structure of various ellipsoids. each has a name, and
% a vector consisting of equatorial radius and flattening. the first
% two are somewhat special cases.

MAP_ELLIP = struct ( 'normal', [1.0, 0], ...
    'sphere', [6370997.0, 0], ...
    'grs80' , [6378137.0, 1/298.257], ...
    'grs67' , [6378160.0, 1/247.247], ...
    'wgs84' , [6378137.0, 1/298.257], ...
    'wgs72' , [6378135.0, 1/298.260], ...
    'wgs66' , [6378145.0, 1/298.250], ...
    'wgs60' , [6378165.0, 1/298.300], ...
    'clrk66', [6378206.4, 1/294.980], ...
    'clrk80', [6378249.1, 1/293.466], ...
    'intl24', [6378388.0, 1/297.000], ...
    'intl67', [6378157.5, 1/298.250]);

name={'UTM'};

switch optn

  case 'name'

    X=name;

  case {'usage','set'}

    m_names=fieldnames(MAP_ELLIP);

    X=char({['     ''' varargin{1} ''''],...
	'     <,''lon<gitude>'',[min max]>',...
	'     <,''lat<itude>'',[min max]>',...
	'     <,''zon<e>'',value>',...
	'     <,''hem<isphere>'',[1|0] (0 for N)>',...
	'     <,''ell<ipsoid>'', one of',...
    reshape(sprintf('         %6s',m_names{:}),15,length(m_names))',...
        '                               >',...
	'     <,''rec<tbox>'', ( ''on'' | ''off'' )>'});

  case 'get'

    X=char([' Projection: ' MAP_PROJECTION.name '  (function: ' ...
	  MAP_PROJECTION.routine ')'],...
	[' longitudes: ' num2str(MAP_VAR_LIST.ulongs) ],...
	[' latitudes: ' num2str(MAP_VAR_LIST.ulats) ],...
	[' zone: ' num2str(MAP_VAR_LIST.zone) ],...
	[' hemisphere: ' num2str(MAP_VAR_LIST.hemisphere) ],...
	[' ellipsoid: ' MAP_VAR_LIST.ellipsoid ], ...
	[' Rectangular border: ' MAP_VAR_LIST.rectbox ]);

  case 'initialize'

    MAP_VAR_LIST=[];
    MAP_PROJECTION.name=varargin{1};
    MAP_VAR_LIST.ulongs = [-72 -68];
    MAP_VAR_LIST.ulats = [40 44];
    MAP_VAR_LIST.zone = 0;		% will be computed if not there
    MAP_VAR_LIST.hemisphere = -1;
    MAP_VAR_LIST.ellipsoid = 'wgs84';
    MAP_VAR_LIST.rectbox='off';
    k=2;

    while k<length(varargin)
      switch varargin{k}(1:3)
	case 'lon'
	  MAP_VAR_LIST.ulongs=varargin{k+1}(:)';
	case 'lat'
	  MAP_VAR_LIST.ulats=varargin{k+1}(:)';
	case 'zon'
	  MAP_VAR_LIST.zone=varargin{k+1};
	case 'hem'
	  MAP_VAR_LIST.hemisphere=varargin{k+1};
        case 'ell'
	  MAP_VAR_LIST.ellipsoid=varargin{k+1};
        case 'rec'
          MAP_VAR_LIST.rectbox=varargin{k+1};
	otherwise
	  disp(['Unknown option: ' varargin{k}]);
      end
      k=k+2;
    end

    % set the zone and hemisphere if not specified

    if (~MAP_VAR_LIST.zone)
      MAP_VAR_LIST.zone = 1 + fix((mod(mean(MAP_VAR_LIST.ulongs)+180,360))/6);
    end

    if MAP_VAR_LIST.hemisphere==-1
      MAP_VAR_LIST.hemisphere=mean(MAP_VAR_LIST.ulats)<0;
    end

    % check for a valid ellipsoid. if not, use the normalized sphere

    if ~isfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid)
       MAP_VAR_LIST.ellipsoid = 'normal';
    end

    % Get X/Y and (if rectboxs are desired) recompute lat/long limits.

    mu_util('xylimits');
    if strcmp(MAP_VAR_LIST.rectbox,'on'), mu_util('lllimits'); end

  case 'll2xy'

    long=varargin{1};
    lat=varargin{2};
    vals=zeros(size(long));

    % Clip out-of-range values (lat/long boxes)

    if  strcmp(MAP_VAR_LIST.rectbox,'off') && ~strcmp(varargin{4},'off')
        vals=vals | long<=MAP_VAR_LIST.longs(1)+eps*10 | long>=MAP_VAR_LIST.longs(2)-eps*10 | ...
	              lat<=MAP_VAR_LIST.lats(1)+eps*10 |   lat>=MAP_VAR_LIST.lats(2)-eps*10;
        [long,lat]=mu_util('clip',varargin{4},long,MAP_VAR_LIST.longs(1),long<MAP_VAR_LIST.longs(1),lat);
        [long,lat]=mu_util('clip',varargin{4},long,MAP_VAR_LIST.longs(2),long>MAP_VAR_LIST.longs(2),lat);
        [lat,long]=mu_util('clip',varargin{4},lat,MAP_VAR_LIST.lats(1),lat<MAP_VAR_LIST.lats(1),long);
        [lat,long]=mu_util('clip',varargin{4},lat,MAP_VAR_LIST.lats(2),lat>MAP_VAR_LIST.lats(2),long);
    end

    % do the forward transformation

    [X,Y] = mu_ll2utm(lat,long,MAP_VAR_LIST.zone,MAP_VAR_LIST.hemisphere, ...
	getfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid));

    % Clip out-of-range values (rectboxes)

    if strcmp(MAP_VAR_LIST.rectbox,'on')  && ~strcmp(varargin{4},'off')
        vals= vals | X<=MAP_VAR_LIST.xlims(1)+eps*10 | X>=MAP_VAR_LIST.xlims(2)-eps*10 | ...
                     Y<=MAP_VAR_LIST.ylims(1)+eps*10 | Y>=MAP_VAR_LIST.ylims(2)-eps*10;
        [X,Y]=mu_util('clip',varargin{4},X,MAP_VAR_LIST.xlims(1),X<MAP_VAR_LIST.xlims(1),Y);
        [X,Y]=mu_util('clip',varargin{4},X,MAP_VAR_LIST.xlims(2),X>MAP_VAR_LIST.xlims(2),Y);
        [Y,X]=mu_util('clip',varargin{4},Y,MAP_VAR_LIST.ylims(1),Y<MAP_VAR_LIST.ylims(1),X);
        [Y,X]=mu_util('clip',varargin{4},Y,MAP_VAR_LIST.ylims(2),Y>MAP_VAR_LIST.ylims(2),X);
    end

  case 'xy2ll'

    [Y,X] = mu_utm2ll(varargin{1}, varargin{2}, MAP_VAR_LIST.zone, ...
	MAP_VAR_LIST.hemisphere, getfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid));

  case 'xgrid'

    [X,Y,vals,labI]=mu_util('xgrid',MAP_VAR_LIST.longs,MAP_VAR_LIST.lats,varargin{1},31,varargin{2:3});

  case 'ygrid'

    [X,Y,vals,labI]=mu_util('ygrid',MAP_VAR_LIST.lats,MAP_VAR_LIST.longs,varargin{1},31,varargin{2:3});

  case 'box'

    [X,Y]=mu_util('box',31);

end

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

function [x,y] = mu_ll2utm (lat,lon, zone, hemisphere,ellipsoid)
%mu_ll2utm		Convert geodetic lat,lon to X/Y UTM coordinates
%
%	[x,y] = mu_ll2utm (lat, lon, zone, hemisphere,ellipsoid)
%
%	input is latitude and longitude vectors, zone number,
%		hemisphere(N=0,S=1), ellipsoid info [eq-rad, flat]
%	output is X/Y vectors
%
%	see also	mu_utm2ll, utmzone

% some general constants

DEG2RADS    = 0.01745329252;
RADIUS      = ellipsoid(1);
FLAT        = ellipsoid(2);
K_NOT       = 0.9996;
FALSE_EAST  = 500000;
FALSE_NORTH = 10000000;

% check for valid numbers

if (max(abs(lat)) > 90)
  error('latitude values exceed 89 degree');
  return;
end

if ((zone < 1) || (zone > 60))
  error ('utm zones only valid from 1 to 60');
  return;
end

% compute some geodetic parameters

lambda_not  = ((-180 + zone*6) - 3) * DEG2RADS;

e2  = 2*FLAT - FLAT*FLAT;
e4  = e2 * e2;
e6  = e4 * e2;
ep2 = e2/(1-e2);

% some other constants, vectors

lat = lat * DEG2RADS;
lon = lon * DEG2RADS;

sinL = sin(lat);
tanL = tan(lat);
cosL = cos(lat);

T = tanL.*tanL;
C = ep2 * (cosL.*cosL);
A = (lon - lambda_not).*cosL;
A2 = A.*A;
A4 = A2.*A2;
S = sinL.*sinL;

% solve for N

N = RADIUS ./ (sqrt (1-e2*S));

% solve for M

M0 = 1 - e2*0.25 - e4*0.046875 - e6*0.01953125;
M1 = e2*0.375 + e4*0.09375 + e6*0.043945313;
M2 = e4*0.05859375 + e6*0.043945313;
M3 = e6*0.011393229;
M = RADIUS.*(M0.*lat - M1.*sin(2*lat) + M2.*sin(4*lat) - M3.*sin(6*lat));

% solve for x

X0 = A4.*A/120;
X1 = 5 - 18*T + T.*T + 72*C - 58*ep2;
X2 = A2.*A/6;
X3 = 1 - T + C;
x = N.*(A + X3.*X2 + X1.* X0);

% solve for y

Y0 = 61 - 58*T + T.*T + 600*C - 330*ep2;
Y1 = 5 - T + 9*C + 4*C.*C;

y = M + N.*tanL.*(A2/2 + Y1.*A4/24 + Y0.*A4.*A2/720);

% finally, do the scaling and false thing. if using a unit-normal radius,
% we don't bother.

x = x*K_NOT + (RADIUS>1) * FALSE_EAST;

y = y*K_NOT;
if (hemisphere)
  y = y + (RADIUS>1) * FALSE_NORTH;
end

return

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

function [lat,lon] = mu_utm2ll (x,y, zone, hemisphere,ellipsoid)
%mu_utm2ll		Convert X/Y UTM coordinates to geodetic lat,lon
%
%	[lat,lon] = mu_utm2ll (x,y, zone, hemisphere,ellipsoid)
%
%	input is X/Y vectors, zone number, hemisphere(N=0,S=1),
%		ellipsoid info [eq-rad, flat]
%	output is lat/lon vectors
%
%	see also	mu_ll2utm, utmzone

% some general constants

DEG2RADS    = 0.01745329252;
RADIUS      = ellipsoid(1);
FLAT        = ellipsoid(2);
K_NOT       = 0.9996;
FALSE_EAST  = 500000;
FALSE_NORTH = 10000000;

if ((zone < 1) || (zone > 60))
  error ('utm zones only valid from 1 to 60');
  return;
end

% compute some geodetic parameters

e2  = 2*FLAT - FLAT*FLAT;
e4  = e2 * e2;
e6  = e4 * e2;
eps = e2 / (1-e2);
em1 = sqrt(1-e2);
e1  = (1-em1)/(1+em1);
e12 = e1*e1;

lambda_not  = ((-180 + zone*6) - 3) * DEG2RADS;

% remove the false things

x = x - (RADIUS>1)*FALSE_EAST;
if (hemisphere)
  y = y - (RADIUS>1)*FALSE_NORTH;
end

% compute the footpoint latitude

M = y/K_NOT;
mu = M/(RADIUS * (1 - 0.25*e2 - 0.046875*e4 - 0.01953125*e6));
foot = mu + (1.5*e1 - 0.84375*e12*e1)*sin(2*mu) ...
    + (1.3125*e12 - 1.71875*e12*e12)*sin(4*mu) ...
    + (1.57291666667*e12*e1)*sin(6*mu) ...
    + (2.142578125*e12*e12)*sin(8*mu);

% some other terms

sinF = sin(foot);
cosF = cos(foot);
tanF = tan(foot);

N = RADIUS ./ sqrt(1-e2*(sinF.*sinF));
T = tanF.*tanF;
T2 = T.*T;
C = eps * cosF.*cosF;
C2 = C.*C;
denom = sqrt(1-e2*(sinF.*sinF));
R = RADIUS * em1*em1 ./ (denom.*denom.*denom);
D = x./(N*K_NOT);
D2 = D.*D;
D4 = D2.*D2;

% can now compute the lat and lon

lat = foot - (N.*tanF./R) .* (0.5*D2 - (5 + 3*T + 10*C - 4*C2 - 9*eps).*D4/24 ...
    + (61 + 90*T + 298*C + 45*T2 - 252*eps - 3*C2) .* D4 .* D2/720);

lon = lambda_not + (D - (1 + 2*T +C).*D2.*D/6 + ...
    (5 - 2*C + 28*T - 3*C2 + 8*eps + 24*T2).*D4.*D./120)./cosF;

% convert back to degrees;

lat=lat/DEG2RADS;
lon=lon/DEG2RADS;

return

author: goGPS-Project


text:
function respStruct = mrMeshTcpParseResponseData(respByteStream)
%
% respStruct = mrMeshTcpParseResponseData(respByteStream)
%
% Parses the mrMesh-encoded byteStream into a matlab structure.
%
% HISTORY
% 2007.04.12 RFD wrote it.

respStruct = [];
if(isempty(respByteStream))
    return;
end

done = false;
while(~done)
    % max param name length is 200
    numBytesLeft = length(respByteStream);
    eq = strfind(respByteStream(1:min(204,numBytesLeft)),' = ');
    if(isempty(eq))
        done = true;
    else
        fn = char(respByteStream(1:eq(1)-1));
        if(fn(end)==']')
            % then it's an array
            sep = strfind(fn,'[');
            dims = str2num(fn(sep(1):end));
            fn = fn(1:sep(1)-1);
            % read all the returned bytes as doubles (8 bytes per element)
            dataEnd = eq+3+prod(dims)*8+1;
            data = respByteStream(eq+3:dataEnd);
            % the data block should begin and end with a ' (ascii=39)
            if(data(1)~=39||data(end)~=39)
                done = true;
                warning('response parsing teminated abnormally.');
            else
                data = typecast(data(2:end-1),'double');
                if(length(dims)>1)
                    data = reshape(data,dims);
                end
                respStruct = setfield(respStruct, fn, data);
                respByteStream = respByteStream(dataEnd(1)+2:end);
            end
        else
            % scalar
            % MaxDoubleAsString = 40
            dataEnd = strfind(respByteStream(1:min(245,numBytesLeft)),sprintf('\n'));
            data = char(respByteStream(eq+3:dataEnd(1)));
            respStruct = setfield(respStruct, fn, str2double(data));
            respByteStream = respByteStream(dataEnd(1)+1:end);
        end
    end
end
if isempty(respStruct)
    respStruct=respByteStream;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% POP_IMPORTEGIMAT - import EGI Matlab segmented file
%
% Usage:
%   >> EEG = pop_importegimat(filename, srate, latpoint0);
%
% Inputs:
%   filename    - Matlab file name
%   srate       - sampling rate
%   latpoint0   - latency in sample ms of stimulus presentation.
%                 When data files are exported using Netstation, the user specify
%                 a time range (-100 ms to 500 ms for instance). In this
%                 case, the latency of the stimulus is 100 (ms). Default is 0 (ms)
%
% Output:
%   EEG        - EEGLAB dataset structure
%
% Authors: Arnaud Delorme, CERCO/UCSD, Jan 2010

% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [EEG com] = pop_importegimat(filename, srate, latpoint0, dataField);

    EEG = [];
    com = '';
    if nargin < 3, latpoint0 = 0; end
    if nargin < 4, dataField = 'Session'; end
    if nargin < 1
        % ask user
        [filename, filepath] = uigetfile('*.mat', 'Choose a Matlab file from Netstation -- pop_importegimat()');
        if filename == 0 return; end
        filename = fullfile(filepath, filename);
        tmpdata = load('-mat', filename);
        fieldValues = fieldnames(tmpdata);
        sessionPos = strmatch('Session', fieldValues);
        posFieldData = 1;
        if ~isempty(sessionPos), posFieldData = sessionPos; end

        if ~isfield(tmpdata, 'samplingRate'), srate = 250; else srate = tmpdata.samplingRate; end
        % epoch data files only
        promptstr    = { { 'style' 'text'       'string' 'Sampling rate (Hz)' } ...
                         { 'style' 'edit'       'string'  int2str(srate) } ...
                         { 'style' 'text'       'string' 'Sample latency for stimulus (ms)' } ...
                         { 'style' 'edit'       'string' '0' } ...
                         { 'style' 'text'       'string' 'Field containing data' } ...
                         { 'style' 'popupmenu'  'string' fieldValues 'value' posFieldData } ...
                         };
        geometry = { [1 1] [1 1] [1 1] };
        result       = inputgui( 'geometry', geometry, 'uilist', promptstr, ...
                                 'helpcom', 'pophelp(''pop_importegimat'')', ...
                                 'title', 'Import a Matlab file from Netstation -- pop_importegimat()');

        if length(result) == 0 return; end;
        srate = str2num(result{1});
        latpoint0 = str2num(result{2});
        dataField = fieldValues{result{3}};
        if isempty(latpoint0), latpoint0 = 0; end
    end

    EEG = eeg_emptyset;
    fprintf('Reading EGI Matlab file %s\n', filename);
    tmpdata = load('-mat', filename);
    if isfield(tmpdata, 'samplingRate') % continuous file
        srate = tmpdata.samplingRate;
    end

    fieldValues = fieldnames(tmpdata);
    if all(cellfun(@(x)isempty(findstr(x, 'Segment')), fieldValues))
        EEG.srate = srate;
        indData = strmatch(dataField, fieldValues);
        EEG.data = tmpdata.(fieldValues{indData(1)});
        EEG = eeg_checkset(EEG);
        EEG = readegilocs(EEG);

        com = sprintf('EEG = pop_importegimat(''%s'');', filename);
    else
        % get data types
        % --------------
        allfields = fieldnames(tmpdata);
        for index = 1:length(allfields)
            allfields{index} = allfields{index}(1:findstr(allfields{index}, 'Segment')-2);
        end
        datatypes = unique_bc(allfields);
        datatypes(cellfun(@isempty, datatypes)) = [];

        % read all data
        % -------------
        counttrial = 1;
        EEG.srate  = srate;
        latency = (latpoint0/1000)*EEG.srate+1;
        for index = 1:length(datatypes)
            tindex = 1;
            for tindex = 1:length(allfields)
                if isfield(tmpdata, sprintf('%s_Segment%d', datatypes{index}, tindex))
                    datatrial = getfield(tmpdata, sprintf('%s_Segment%d', datatypes{index}, tindex));
                    if counttrial == 1
                        EEG.pnts = size(datatrial,2);
                        EEG.data = repmat(single(0), [size(datatrial,1), size(datatrial,2), 1000]);
                    end
                    EEG.data(:,:,counttrial) = datatrial;

                    EEG.event(counttrial).type    = datatypes{index};
                    EEG.event(counttrial).latency = latency;
                    EEG.event(counttrial).epoch   = counttrial;

                    counttrial = counttrial+1;
                    latency = latency + EEG.pnts;
                end
            end
        end
        fprintf('%d trials read\n', counttrial-1);
        EEG.data(:,:,counttrial:end) = [];

        EEG.setname  = filename(1:end-4);
        EEG.nbchan   = size(EEG.data,1);
        EEG.trials   = counttrial-1;
        if latpoint0 ~= 1
            EEG.xmin     = -latpoint0/1000;
        end
        EEG = eeg_checkset(EEG);

        % channel location
        % ----------------
        if all(EEG.data(end,1:10) == 0)
            disp('Deleting empty data reference channel (reference channel location is retained)');
            EEG.data(end,:)   = [];
            EEG.nbchan        = size(EEG.data,1);
            EEG = eeg_checkset(EEG);
        end
        EEG = readegilocs(EEG);

        com = sprintf('EEG = pop_importegimat(''%s'', %3.2f, %3.2f, %d);', filename, srate, latpoint0, dataField);
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function index = subv2ind(siz,sub)
%SUBV2IND   Linear index from subscript vector.
% SUBV2IND(SIZ,SUB) returns an equivalent single index corresponding to a
% subscript vector for an array of size SIZ.
% If SUB is a matrix, with subscript vectors as rows, then the result is a
% column vector.
%
% This is the opposite of IND2SUBV, so that
%   SUBV2IND(SIZ,IND2SUBV(SIZ,IND)) == IND.
%
% See also IND2SUBV, SUB2IND.

%index = subv2indTest(siz,sub);
prev_cum_size = [1 cumprod(siz(1:end-1))];
%index = (sub-1)*prev_cum_size' + 1;
index = sub*prev_cum_size' - sum(prev_cum_size) + 1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% testWriteSegy : Script to test WriteSegy and WriteSegyStructure
%
n=200;
seisdata=peaks(n);
imagesc(seisdata);colorbar;
xlabel('Traces');
ylabel('Time')
title('The data that should be written to disk')

% use WriteSegy to write an IEEE Revision 1 formatted file
%
WriteSegy('test.segy',seisdata);

% Read the file we just created with all of its header values
%
[Data,SegyTraceHeaders,SegyHeader,HeaderInfo]=ReadSegy('test.segy');

% NOW WE CAN WRITE ALL THE FORMATS SUPPORTED BY SegyMAT :

% WRITE IBM FLOATING POINT REV 0
% THIS IS THE ONLY IMPLEMENTED FORMAT FOR REV 0
SegyHeader.SegyFormatRevisionNumber=0;
SegyHeader.DataSampleFormat=1;
WriteSegyStructure('data_IBM_REV0.segy',SegyHeader,SegyTraceHeaders,Data);

% WRITE IBM FLOATING POINT REV 1
SegyHeader.DataSampleFormat=1;
WriteSegyStructure('data_IBM_REV1.segy',SegyHeader,SegyTraceHeaders,Data);

% WRITE 4BYTE INT POINT REV 1
SegyHeader.DataSampleFormat=2;
WriteSegyStructure('data_4byteINT.segy',SegyHeader,SegyTraceHeaders,Data);

% WRITE 2BYTE INT POINT REV 1
SegyHeader.DataSampleFormat=3;
WriteSegyStructure('data_2byteINT.segy',SegyHeader,SegyTraceHeaders,Data);

% WRITE 1BYTE INT POINT REV 1
SegyHeader.DataSampleFormat=8;
WriteSegyStructure('data_1byteINT.segy',SegyHeader,SegyTraceHeaders,Data);

% WRITE IEEE POINT REV 1
SegyHeader.DataSampleFormat=5;
WriteSegyStructure('data_IEEE.segy',SegyHeader,SegyTraceHeaders,Data);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function planC = importDICOM(source,destn,recursiveFlag)
% function importDICOM(source,destn)
%
% Function to import DICOM from source and write CERR file to destn.
%
% APA, 8/14/2018

if ~exist('recursiveFlag','var')
    recursiveFlag = false;
end

zipFlag = 'No';

% % Define file and directory filter patterns
% fileFilterC = {'scout','corTmp','sag','report','_nac','mip',...
%     'coronal','sagittal','screen_save'};
% dirFilterC = {'ot'};

% Get all directories and files
tic,
[filesInCurDir,dirsInCurDir] = rdir(source);
toc

% Convert to lower case
dirsToImportC = lower({dirsInCurDir.fullpath});
if isempty(dirsToImportC)
    dirsToImportC = {source};
end

% filter directories containing no files
indV1 = false(1,length(dirsToImportC));
tic,
for dirNum = 1:length(dirsToImportC)
    dirS = dir(dirsToImportC{dirNum});
    if ~all([dirS.isdir])
        indV1(dirNum) = 1;
    end
end
toc
dirsToImportC = dirsToImportC(indV1);

%% Import DICOM to CERR
% Read options file
pathStr = getCERRPath;
optName = [pathStr,'CERROptions.json'];
optS = opts4Exe(optName);

% Import all the dirs
for dirNum = 1:length(dirsToImportC)
    try
        init_ML_DICOM
        %hWaitbar = waitbar(0,'Scanning Directory Please wait...');
        sourceDir = dirsToImportC{dirNum};
        %patient = scandir_mldcm_babs(sourceDir, hWaitbar, 1);
        excludePixelDataFlag = true;
        patient = scandir_mldcm(sourceDir, excludePixelDataFlag, recursiveFlag);
        %close(hWaitbar);
        dcmdirS = struct(['patient_' num2str(1)],patient.PATIENT(1));
        for j = 2:length(patient.PATIENT)
            dcmdirS.(['patient_' num2str(j)]) = patient.PATIENT(j);
        end
        patNameC = fieldnames(dcmdirS);
        mergeScansFlag = 'No';
        combinedDcmdirS = struct('STUDY',dcmdirS.(patNameC{1}).STUDY);
        combinedDcmdirS = struct('STUDY',dcmdirS.(patNameC{1}).STUDY,'info',dcmdirS.(patNameC{1}).info);
        for i = 2:length(patNameC)
            for j = 1:length(dcmdirS.(patNameC{i}).STUDY.SERIES)
                combinedDcmdirS.STUDY.SERIES(end+1) = ...
                    dcmdirS.(patNameC{i}).STUDY.SERIES(j);
            end
        end
        % Pass the java dicom structures to function to create CERR plan
        try
            planC = dcmdir2planC(combinedDcmdirS,mergeScansFlag,optS);
        end

        indexS = planC{end};

        % build the filename for storing planC
        if sourceDir(end) == filesep
            [~,folderNam] = fileparts(sourceDir(1:end-1));
        else
            [~,folderNam] = fileparts(sourceDir);
        end

        %mrn = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.PatientID;
        %studyDscr = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.StudyDescription;
        %seriesDscr = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.SeriesDescription;
        %modality = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.Modality;
        % outFileName = [mrn,'~',studyDscr,'~',seriesDscr,'~',modality];

        outFileName = folderNam;

        %outFileName = mrn;   % store file names as MRNs

        %[~,outFileName] = fileparts(sourceDir);  % store file names as DICOM directory names
        %fullOutFileName = fullfile(destn,fileName);

        %Check for duplicate name of fullOutFileName
        if exist('destn','var') && ~isempty(destn)
            dirOut = dir(destn);
            allOutNames = {dirOut.name};
            if any(strcmpi([outFileName,'.mat'],allOutNames))
                fullOutFileName = [outFileName,'_duplicate_',num2str(rand(1))];
            end
            if strcmpi(zipFlag,'Yes')
                saved_fullFileName = fullfile(destn,[outFileName,'.mat.bz2']);
            else
                saved_fullFileName = fullfile(destn,[outFileName,'.mat']);
            end

            if ~exist(fileparts(saved_fullFileName),'dir')
                mkdir(fileparts(saved_fullFileName))
            end
            save_planC(planC,[], 'passed', saved_fullFileName);
        end

    catch

        disp(['Cannot convert ',source])

    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function plotyy(obj1, obj2, varargin)
   p = inputParser;
   p.addParameter('snum',max([obj1.dnum(1) obj2.dnum(1)]));
   p.addParameter('enum', min([obj1.dnum(end) obj2.dnum(end)]));
   p.addParameter('fun1','plot');
   p.addParameter('fun2','plot');
   p.parse(varargin{:});
   Args = p.Results;

    [ax, ~, ~] = plotyy(obj1.dnum, obj1.data, obj2.dnum, obj2.data, Args.fun1, Args.fun2);
    datetick('x');
    set(ax(2), 'XTick', [], 'XTickLabel', {});
    set(ax(1), 'XLim', [Args.snum Args.enum]);
end

author: geoscience-community-codes

repo: GISMO


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

% Cleaning
clear all
close all
clc

% Library path
run('../../addpathGypsilab')

% Dimensions
N = 200;

% Accuracy
tol = 1e-3;

% Particles receptors X
mesh  = mshSphere(N,1);
omega = dom(mesh,3);
phi0  = fem(mesh,'P0');
phi1  = fem(mesh,'P1');

% Wave number or frequency (Hz)
k = 5;

% Green kernel -> exp(1i*k*r)/r
green = '[exp(ikr)/r]';
Gxy   = @(X,Y) femGreenKernel(X,Y,green,k) + ...
    sqrt(N).*(X(:,1)==Y(:,1)).*(X(:,2)==Y(:,2)).*(X(:,3)==Y(:,3))  ;

% Particles charges (multiples)
V0 = (-1+2*rand(length(phi0),2)) + (-1+2i*rand(length(phi0),2));
V1 = (-1+2*rand(length(phi1),2)) + (-1+2i*rand(length(phi1),2));

% Spatial representation of particles
figure
plot(mesh)
axis equal

% All forms
tic
Ah = integral(omega,omega,phi0,Gxy,phi0,tol);
Bv = integral(omega,omega,phi0,green,k,phi1,tol);
Cs = integral(omega,phi1,phi0);
Df = integral(omega,omega,phi1,Gxy,phi1);
toc

%%% Single Builder
disp('~~~~~~~~~~~~~ SINGLE BUILDERS ~~~~~~~~~~~~~')
Mb  = bmm(Ah);
sol = Mb * V0;
ref = Ah * V0;
norm(ref-sol,'inf')/norm(ref,'inf')

Mb  = bmm(Bv);
sol = Mb * V1;
ref = Bv * V1;
norm(ref-sol,'inf')/norm(ref,'inf')

Mb  = bmm(Cs);
sol = Mb * V0;
ref = Cs * V0;
norm(ref-sol,'inf')/norm(ref,'inf')

Mb  = bmm(Df);
sol = Mb * V1;
ref = Df * V1;
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Multi Builder
disp('~~~~~~~~~~~~~ MULTI BUILDERS ~~~~~~~~~~~~~')
Mb  = bmm({Ah,Cs';Cs,Df});
Mr  = [full(Ah) Cs' ; Cs Df];
Vb  = [V0;V1];
sol = Mb * Vb;
ref = Mr * Vb;
norm(ref-sol,'inf')/norm(ref,'inf')

figure
spy(Mb)

disp(' ')

%%% Full conversion
disp('~~~~~~~~~~~~~ FULL CONVERSION ~~~~~~~~~~~~~')
sol = full(Mb) * Vb;
ref = Mr * Vb;
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Sparse conversion
disp('~~~~~~~~~~~~~ SPARSE CONVERSION ~~~~~~~~~~~~~')
sol = sparse(Mb) * Vb;
ref = Mr * Vb;
norm(ref-sol,'inf')/norm(ref,'inf')

figure
spy(sparse(Mb))

disp(' ')

%%% Transposition
disp('~~~~~~~~~~~~~ TRANSPOSITION ~~~~~~~~~~~~~')
sol = Vb.' * Mb.';
ref = (Mr * Vb).';
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Transposition
disp('~~~~~~~~~~~~~ CONJUGATE TRANSPOSITION ~~~~~~~~~~~~~')
sol = Vb' * Mb';
ref = (Mr * Vb)';
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Concatenation
disp('~~~~~~~~~~~~~ CONCATENATION ~~~~~~~~~~~~~')
Ab  = Mb;
Bb  = bmm({Bv;Df});
Cb  = bmm({Ah,Bv});
Db  = bmm(Cs');
Mb2 = [Ab , Bb ; Cb , Db];
sol = Mb2 * [Vb;V1];
ref = [Ab*Vb + Bb*V1;Cb*Vb+Db*V1];
norm(ref-sol,'inf')/norm(ref,'inf')

spy(Mb2)

disp(' ')

%%% Scalar product
disp('~~~~~~~~~~~~~ SCALAR PRODUCT ~~~~~~~~~~~~~')
sol = (-(3.*Mb.*2i)) * Vb;
ref = (-3*2i) * (Mr*Vb);
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Addition
disp('~~~~~~~~~~~~~ ADDITION ~~~~~~~~~~~~~')
sol = (Mb+Mb) * Vb;
ref = 2*(Mr*Vb);
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Multiplication
disp('~~~~~~~~~~~~~ MULTIPLICATION ~~~~~~~~~~~~~')
sol = (Mb.'*Mb) * Vb;
ref = Mr.'*(Mr*Vb);
norm(ref-sol,'inf')/norm(ref,'inf')

figure
spy(Mb.'*Mb)

disp(' ')

%%% LU factorization
disp('~~~~~~~~~~~~~ LU FACTORISATION ~~~~~~~~~~~~~')
[Lb,Ub] = lu(Mb);
sol = Lb * (Ub * Vb);
ref = Mb * Vb;
norm(ref-sol,'inf')/norm(ref,'inf')

figure
subplot(1,2,1); spy(Lb)
subplot(1,2,2); spy(Ub)

disp(' ')

%%% Solve LU
disp('~~~~~~~~~~~~~ SOLVE LU SYSTEM ~~~~~~~~~~~~~')
sol = Mb \ Vb;
ref = Mr \ Vb;
norm(ref-sol,'inf')/norm(ref,'inf')

disp(' ')

%%% Inversion
disp('~~~~~~~~~~~~~ INVERSION ~~~~~~~~~~~~~')
Mbm1 = inv(Mb);
sol  = Mbm1 * Vb;
ref  = Mr \ Vb;
norm(ref-sol,'inf')/norm(ref,'inf')

figure
spy(Mbm1)

disp(' ')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function value = imDynamicRange(img, bRobust, type)
%
%
%        value = imDynamicRange(img, bRobust, type)
%
%
%        Input:
%           -img: the input image
%           -bRobust: if bRobust > 0 --> robust statistics for min and max
%                     luminance values. bRboust becomes the percentile!
%           -type: 'Classic', 'Michelson', and 'Weber'
%
%        Output:
%           -value(1): dynamic range of img
%           -value(2): dynamic range of img in f-stops (if 'Classic')
%           -value(3): dynamic range of img in log10 space space (if
%           'Classic')
%
%     Copyright (C) 2011-20  Francesco Banterle
%
%     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/>.
%

if(~exist('bRobust', 'var'))
    bRobust = 0;
end

if(~exist('type', 'var'))
    type = 'Classic';
end

if(bRobust >= 0.5)
   bRobust = 0.01;
end

L = lum(img);

if(bRobust > 0.0)
    minL = MaxQuart(L, bRobust);
    maxL = MaxQuart(L, 1 - bRobust);
else
    minL = min(L(:));
    maxL = max(L(:));
end

if(minL < 1e-6)
    warning('minL is less than 1e-6 cd/m^2');
end

if(minL <= 0.0)
    warning('minL is 0.0 is set to the first value greater than zero.');

    minL = min(min(L(L > 0)));
end

switch type

    case 'Classic'
        value(1) = maxL / minL;
        value(2) = log2(value(1));
        value(3) = log10(value(1));

    case 'Michelson'
        value = (maxL - minL) / (maxL + minL);

    case 'Weber'
        value = (maxL - minL) / minL;
end

end

author: banterle


text:
bg = 0.7;
im = bg*ones(50,50);

im = ipaste(im, ones(20,20), [5,5]);
im = ipaste(im, ones(10,15), [30,30]);
c = kcircle(10);
c(c==0) = bg;
im = ipaste(im, c, [5,27]);

author: petercorke

repo: machinevision-toolbox-matlab


text:
function res = bf_output_image_mv(BF, S)
% Computes multivariate test on a number of frequency bands
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes, modified from Vladimir Litvak's example code
% $Id: bf_output_image_mv.m 7703 2019-11-22 12:06:29Z guillaume $

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

%% no covariance matrix creation -so need to check bands are within this window

if nargin == 0
    all = cfg_const;
    all.tag = 'all';
    all.name = 'All';
    all.val  = {1};

    condlabel = cfg_entry;
    condlabel.tag = 'condlabel';
    condlabel.name = 'Condition label';
    condlabel.strtype = 's';
    condlabel.val = {''};

    conditions = cfg_repeat;
    conditions.tag = 'conditions';
    conditions.name = 'Conditions';
    conditions.help = {'Specify the labels of the conditions to be included in the inversion'};
    conditions.num  = [1 Inf];
    conditions.values  = {condlabel};
    conditions.val = {condlabel};

    whatconditions = cfg_choice;
    whatconditions.tag = 'whatconditions';
    whatconditions.name = 'What conditions to include?';
    whatconditions.values = {all, conditions};
    whatconditions.val = {all};

    %     design = cfg_const;
    %     design.tag = 'design';
    %     design.name = 'design';
    %     design.help = {'Use default settings for the inversion'};
    %     design.val  = {1};

    design = cfg_files;
    design.tag = 'design';
    design.name = 'design matrix';
    design.filter = 'mat';
    design.num = [1 Inf];
    design.help = {'Select the design matrix'};

    woi = cfg_entry;
    woi.tag = 'woi';
    woi.name = 'Time windows of interest';
    woi.strtype = 'r';
    woi.num = [Inf 2];
    woi.help = {'Time windows (in ms)'};

    foi = cfg_entry;
    foi.tag = 'foi';
    foi.name = 'Frequency bands of interest';
    foi.strtype = 'r';
    foi.num = [Inf 2];
    foi.help = {'Freq bands (in Hz)'};

    datafeatures = cfg_menu;
    datafeatures.tag = 'datafeatures';
    datafeatures.name = 'Features';
    datafeatures.labels =get_data_features;
    datafeatures.values =get_data_features;

    datafeatures.help = {'Data features of interest'};
    datafeatures.val =datafeatures.values(end);

    contrast = cfg_entry;
    contrast.tag = 'contrast';
    contrast.name = 'contrast';
    contrast.strtype = 'i';
    contrast.num = [1 Inf];
    contrast.val = {1};

    result         = cfg_menu;
    result.tag     = 'result';
    result.name    = 'What to output';
    result.help    = {'Specify output type'};
    result.labels  = {
        'chi square'
        'BIC'
        'r square'
        }';
    result.values  = result.labels;
    result.val = {'chi square'};

    modality         = cfg_menu;
    modality.tag     = 'modality';
    modality.name    = 'Modality';
    modality.help    = {'Specify modality'};
    modality.labels  = {
        'MEG'
        'MEGPLANAR'
        'EEG'
        }';
    modality.values  = {
        'MEG'
        'MEGPLANAR'
        'EEG'
        }';
    modality.val = {'MEG'};

    custom = cfg_branch;
    custom.tag = 'custom';
    custom.name = 'Custom';
    custom.help = {'Define custom settings for the inversion'};
    custom.val  = {whatconditions, contrast, woi};

    isdesign = cfg_choice;
    isdesign.tag = 'isdesign';
    isdesign.name = 'Design matrix parameters';
    isdesign.help = {'Choose whether to load custom design'};
    isdesign.values = {design, custom};
    isdesign.val = {design};

    sametrials = cfg_menu;
    sametrials.tag = 'sametrials';
    sametrials.name = 'Trials same as for filters';
    sametrials.labels = {'yes', 'no'};
    sametrials.values = {true, false};
    sametrials.val = {false};
    sametrials.help = {'Take the same trials as used for filter computation',...
        'This is useful for bootstrap.'};

    image_mv      = cfg_branch;
    image_mv.tag  = 'image_mv';
    image_mv.name = 'Mv image';
    image_mv.val  = { isdesign, datafeatures, foi,  result, sametrials, modality};

    res = image_mv;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

D = BF.data.D;

if isfield(S.isdesign,'custom'),
    %% gui specified conditions and contrast

    woitmp = S.isdesign.custom.woi;

    % DP - The indsample function does not work for matrices, so I have
    % looped through. Otherwise, we are left with only one windows of
    % interest.
    for wi = 1:size(woitmp,1)
        woiind(wi,:)=D.indsample(woitmp(wi,:)/1000);
        woi(wi,:)=D.time(woiind(wi,:)); %% in seconds
    end

    duration=unique(woiind(:,2)-woiind(:,1))./D.fsample; %% in sec
    if numel(duration)>1,
        error('both windows need to be the same length');
    end;

    % DP - this can result in a tiny bit of residual, presumably due to
    % numerical imprecision. Not sure. Doesn't seem necessary anyway since
    % the number of samples is equal.
%     duration=unique(woi(:,2)-woi(:,1));
%     if numel(duration)>1,
%         error('both windows need to be the same length');
%     end;
%     woi = S.isdesign.custom.woi;
%     woiind=D.indsample(woi/1000);
%     woi=D.time(woiind); %% in seconds

%     duration=unique(woiind(:,2)-woiind(:,1))./D.fsample; %% in sec
%     if numel(duration)>1,
%         error('both windows need to be the same length');
%     end;

    duration=unique(woi(:,2)-woi(:,1));
    if numel(duration)>1,
        error('both windows need to be the same length');
    end;

    %
    whatconditions=S.isdesign.custom.whatconditions;

    if isfield(whatconditions, 'all')
        if S.sametrials
            trials{1} = BF.features.trials;
        else
            trials{1} = 1:D.ntrials;
        end
        clabel{1}='all';
    else
        for i = 1:numel(whatconditions.condlabel)
            if isempty(D.indtrial(whatconditions.condlabel{i}, 'GOOD'))
                error('No trials matched the selection.');
            end

            clabel{i}=whatconditions.condlabel{i};

            if S.sametrials
                trials{i} = BF.features.trials(strmatch(clabel{i}, D.conditions(BF.features.trials)));
            else
                trials{i} = D.indtrial(whatconditions.condlabel{i}, 'GOOD');
            end

        end
        if isempty(trials)
            error('No trials matched the selection, check the specified condition labels');
        end
    end

    %%check for number of trials //// added by ANNA 30/04/2013
    for i=1:numel(trials)
        num_trials(i) = length(trials{i});
    end

    if min(num_trials)~= max(num_trials)
        warning ('Number of trials are not the same accross conditions- throwing away');
        num_trials
    end
    nt = min(num_trials); % ANNA -throw away the extra trials

    col=0;
    X=[];
    Xtrials=[];
    Xstartlatencies=[];
    xlabel=[];
    for j=1:size(woi, 1),
        for i=1:numel(trials)
            col=col+1;
            %nt=numel(trials{i}); % ANNA look above for the warning
            Xtmp=[zeros(size(X,1),1); ones(nt,1)];
            if col>1,
                X=[X;zeros(nt,size(X,2))]; %1)]; ANNA
            end;
            X=[X Xtmp];
            tlist = spm_vec(trials{i}); % ANNA list of trials
            Xtrials=[Xtrials ; tlist(1:nt)]; % ANNA list of trials - this must be a vector
            Xstartlatencies=[Xstartlatencies; ones(nt,1).*woi(j,1)];%% again this has to be a vector
            xlabel=strvcat(xlabel,[clabel{i} ',' num2str(woi(j,1))]);
        end;
    end;
%     % ORIG
%     col=0;
%     X=[];
%     Xtrials=[];
%     Xstartlatencies=[];
%     xlabel=[];
%     for j=1:size(woi, 1),
%         for i=1:numel(trials)
%             col=col+1;
%             nt=numel(trials{i});
%             Xtmp=[zeros(size(X,1),1); ones(nt,1)];
%             if col>1,
%                 X=[X;zeros(nt,1)];
%             end;
%             X=[X Xtmp];
%             Xtrials=[Xtrials ;[trials{i}]];
%             Xstartlatencies=[Xstartlatencies; ones(nt,1).*woi(j,1)];
%             xlabel=strvcat(xlabel,[clabel{i} ',' num2str(woi(j,1))]);
%         end;
%     end;

    allsamples=[D.indsample(Xstartlatencies); D.indsample(Xstartlatencies+ones(size(Xstartlatencies)).*duration)]';
    contrast=S.isdesign.custom.contrast';

else %%  conditions and contrast  specified in a file

    if ~exist(cell2mat(S.isdesign.design)),
        error('Cannot load design matrix');
    end;
    disp('loading design matrix');
    a=load(cell2mat(S.isdesign.design));
    X=a.design.X; %% design matrix
    contrast=a.design.contrast;
    ntrials=size(X,1);
    if (size(a.design.Xstartlatencies,1)~=ntrials)||(size(a.design.Xtrials,1)~=ntrials)
        error('start latencies and Xtrials and X should have a value per row of the design');
    end;

    Xtrials=a.design.Xtrials; %% indices of trials to use
    for j=1:ntrials,
        allsamples(j,1)=D.indsample(a.design.Xstartlatencies(j));
        allsamples(j,2)=D.indsample(a.design.Xstartlatencies(j)+a.design.Xwindowduration);
    end;

end;

if size(Xtrials)~=size(Xstartlatencies)
    error('Xtrials and start latencies must be the same length');
end;
if size(Xtrials,1)~=size(X,1)
    error('X must have same number of rows as trials and start latencies');
end;
if size(contrast,1)~=size(X,2),
    error('contrast needs to match number of columns in design');
end;

Fgraph  = spm_figure('GetWin','Graphics'); spm_figure('Clear',Fgraph);

subplot(4,1,1);
imagesc(X);
title('X');
set(gca,'Xtick',1:col);
set(gca,'Xticklabel',xlabel);
subplot(4,1,2);

imagesc(contrast);
title('c');
subplot(4,1,3);

imagesc(Xtrials);
title('trials');
subplot(4,1,4);

imagesc(Xstartlatencies);
title('start latency');

chanind = BF.features.(S.modality).chanind;
U       = BF.features.(S.modality).U;

Nchans=size(U,2); %% effective number of channels

nsamples = unique(allsamples(:,2)-allsamples(:,1));

if length(nsamples) > 1
    error('All time windows should be equal lentgh')
end

alltrials = spm_vec(Xtrials);
ntrials   = length(alltrials);

%% now identify frequency bands of interest

Nbands=size(S.foi,1);
windowduration=(nsamples/D.fsample);
dctfreq = (0:nsamples-1)/2/windowduration;           % DCT frequencies (Hz)
dctT      = spm_dctmtx(nsamples,nsamples);

freqstr=[];
allfreqind=[];
for fband=1:Nbands, %% allows one to break up spectrum and ignore some frequencies
    freqrange=S.foi(fband,:);
    j      = find( (dctfreq >= freqrange(1)) & (dctfreq<=freqrange(2)) );
    featureind{fband}=j;
    allfreqind=sort(unique([allfreqind j]));
    freqstr=[freqstr sprintf('%3.1f-%3.1f,',dctfreq(min(j)),dctfreq(max(j)))];
end; % for fband=1:Nbands

if isempty(allfreqind),
    error('No valid frequency range found');
end;
% Tfull      = dctT(:,allfreqind); %% A filter for all bands (not necessarily continuous)
%% end of freq band section

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', ntrials , 'Computing covariance'); drawnow;
if ntrials  > 100, Ibar = floor(linspace(1, ntrials ,100));
else Ibar = 1:ntrials; end

%% load in data and make up simple design matrix
%ncond=numel(samples); %% number of conditions= columns in design matrix
%Nt=ntrials*ncond; %% total number of time windows under consideration
%X=zeros(Nt,ncond);

flatdata=zeros(ntrials*nsamples,Nchans);
%% want flatdata in form Nchans,Nt*Nsamples

%count=0;
for i = 1:ntrials
    %for j = 1:numel(samples)
    %   count=count+1;
    %   X(count,j)=1;
    Y  = U'*squeeze(D(chanind, allsamples(i,1):allsamples(i,2)-1, alltrials(i)));
    Y  = detrend(Y'); %% detrend and throw away low freq drift
    flatdata((i-1)*nsamples+1:i*nsamples,:) =Y;

    %end
    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

W = BF.inverse.(S.modality).W;
nvert = numel(W);
S.regressout=[]; %% turn off for now
regressout=S.regressout;

%% set up the data features
weights=-1; %% set up flag
Yfull=get_data_features(flatdata,nsamples,ntrials,weights,dctT,S.datafeatures,featureind,regressout); %% set up data structures

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nvert, 'Scanning grid points'); drawnow;
if nvert > 100, Ibar = floor(linspace(1, nvert,100));
else Ibar = 1:nvert; end

outval = nan(1, nvert);

for i = 1:nvert

    if ~isnan(W{i})

        w    = W{i};

        %% returns columns of a matrix with rows as different observations
        [Yfull,vedata]=get_data_features(flatdata,nsamples,ntrials,w,dctT,S.datafeatures,featureind,regressout); %% extract the data features

        Yfull=Yfull-repmat(mean(Yfull),size(Yfull,1),1); %% remove dc level from each column/feature
        Yfull=Yfull./repmat(std(Yfull),size(Yfull,1),1); %% normalize features to have unit variance by default
        [chival,BIC,cva] = output_image_mv_cva(X,Yfull,contrast); %% run the multivariate test

        switch S.result
            case 'chi square'
                resultstr='chisq';
                outval(i) = chival(1);
            case 'r square'
                outval(i) = cva.ccorr.^2;
                resultstr='rsq';
            case 'BIC'
                outval(i)=BIC(1);
                resultstr='BIC';
        end;

    end

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

image(1).val   = outval;

image(1).label = ['mv' resultstr S.datafeatures freqstr  spm_file(D.fname, 'basename')];

res = image;

author: spm


text:
classdef TestingOrthogonalCorrectors < handle

    properties (Access = public)

    end

    properties (Access = private)

    end

    properties (Access = private)
       mesh
       orientation
       interpolator
       orthogonalCorrector
    end

    methods (Access = public)

        function obj = TestingOrthogonalCorrectors()
          %  obj.createToyMesh();
          %  obj.createToyOrientation();
            obj.createBenchmarkMesh();
            obj.createBenchmarkOrientation();
            obj.plotOrientationVector();
            obj.createInterpolator();
            obj.createOrthogonalCorrector();
        end

    end

    methods (Access = private)

        function createToyMesh(obj)
            x = [66.90 128.89 115.25 76.73 26.84 157.24 ...
                168.58 141.74 98.65 45.74 3.78 2.65 31.37 ...
                83.53 137.21 174.63 147.03 99.03 38.55];
            y = [89.20 89.58 120.20 130.40 119.44 113.39 ...
                153.08 154.59 160.26 164.80 152.32 220.74 ...
                207.13 202.977 190.88 192.77 239.64 241.90 242.28];
            s.coord(:,1) = x;
            s.coord(:,2) = y;
            s.connec = delaunay(s.coord);
            m = Mesh(s);
            m.plot();
            obj.mesh = m;
        end

        function createToyOrientation(obj)
            alpha = pi/180*[190 140 300 150 185 ...
                            120 75  45  150 130 ...
                            160+180 0   190  5  -5   ...
                            170 0   190 190];
            a(:,1) = cos(alpha);
            a(:,2) = sin(alpha);
            obj.orientation = a;
        end

       function createBenchmarkMesh(obj)
           h = 0.03;
           xmin = 0.50;
           xmax = 2.0;
           ymin = 0.25;
           ymax = 1.75;
           xv = xmin:h:xmax;
           yv = ymin:h:ymax;
           [X,Y] = meshgrid(xv,yv);
           s.coord(:,1) = X(:);
           s.coord(:,2) = Y(:);
           s.connec = delaunay(s.coord);
            m = Mesh(s);
            obj.mesh = m;
        end

        function createBenchmarkOrientation(obj)
            s1 = 0.32;
            s2 = -0.8;
            x1 = obj.mesh.coord(:,1);
            x2 = obj.mesh.coord(:,2);
            v(:,1) = cos(pi*(x1 + s1*x2));
            v(:,2) = cos(pi*(x2 + s2*x1));
            beta = atan2(v(:,2),v(:,1));
            alpha = beta/2;
            obj.orientation(:,1) = cos(alpha);
            obj.orientation(:,2) = sin(alpha);
        end

        function plotOrientationVector(obj)
            %figure()
            a = obj.orientation;
            x = obj.mesh.coord(:,1);
            y = obj.mesh.coord(:,2);
            ax = a(:,1);
            ay = a(:,2);
            quiver(x,y,ax,ay);
        end

        function createInterpolator(obj)
            s.meshCont    = obj.mesh;
            s.meshDisc    = obj.mesh.createDiscontinousMesh;
            s.orientation = obj.orientation;
            s = SymmetricContMapCondition(s);
            sC = s.computeCondition();
            obj.interpolator = sC;
        end

        function createOrthogonalCorrector(obj)
            s.mesh            = obj.mesh;
            s.interpolator    = obj.interpolator;
            s.correctorValue  = obj.computeCorrector();
            o = OrthogonalCorrectorComputer(s);
            oC = o.compute();
            o.plot();
            obj.orthogonalCorrector = oC;
        end

        function cV = computeCorrector(obj)
            s.mesh               = obj.mesh;
            s.orientation        = obj.orientation;
            s.singularityCoord   = obj.computeSingularities();
            c = CorrectorComputer(s);
            cV = c.compute();
            c.plot()
        end

        function sCoord = computeSingularities(obj)
            s.mesh        = obj.mesh;
            s.orientation = obj.orientation;
            sF = SingularitiesFinder(s);
            isS = sF.computeSingularElements();
            sF.plot();
            coordB = obj.mesh.computeBaricenter();
            coordB = transpose(coordB);
            sCoord =  coordB(isS,:);
            sCoord = sCoord(1,:);
        end

    end

end


text: %PLUS. Mapping overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function c = ttv(a,v,dims)
%TTV Tensor times vector for ktensor.
%
%   Y = TTV(X,A,N) computes the product of Kruskal tensor X with a
%   (column) vector A.  The integer N specifies the dimension in X
%   along which A is multiplied.  If size(A) = [I,1], then X must have
%   size(X,N) = I.  Note that ndims(Y) = ndims(X) - 1 because the N-th
%   dimension is removed.
%
%   Y = TTV(X,{A1,A2,...}) computes the product of tensor X with a
%   sequence of vectors in the cell array.  The products are computed
%   sequentially along all dimensions (or modes) of X. The cell array
%   contains ndims(X) vectors.
%
%   Y = TTV(X,{A1,A2,...},DIMS) computes the sequence tensor-vector
%   products along the dimensions specified by DIMS.
%
%   See also TENSOR/TTV, KTENSOR, KTENSOR/TTM.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%%%%%%%%%%%%%%%%%%%%%%
%%% ERROR CHECKING %%%
%%%%%%%%%%%%%%%%%%%%%%

% Check the number of arguments
if (nargin < 2)
    error('TTV requires at least two arguments.');
end

% Check for 3rd argument
if ~exist('dims','var')
    dims = [];
end

% Check that 2nd argument is cell array. If not, recall with v as a
% cell array with one element.
if ~iscell(v)
    c = ttv(a,{v},dims);
    return;
end

% Get sorted dims and index for multiplicands
[dims,vidx] = tt_dimscheck(dims,ndims(a),numel(v));

% Check that each multiplicand is the right size.
for i = 1:numel(dims)
    if ~isequal(size(v{vidx(i)}),[size(a,dims(i)) 1])
        error('Multiplicand is wrong size');
    end
end

% Figure out which dimensions will be left when we're done
remdims = setdiff(1:ndims(a),dims);

% Collapse dimensions that are being multiplied out
newlambda = a.lambda;
for i = 1:numel(dims)
    newlambda = newlambda .* ( a.u{dims(i)}' * v{vidx(i)} );
end

% Create final result
if isempty(remdims)
    c = sum(newlambda);
else
    c = ktensor(newlambda,a.u{remdims});
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function kmlStruct = kml2struct(kmlFile)
% kmlStruct = kml2struct(kmlFile)
%
% Import a .kml file as a vector array of shapefile structs, with Geometry, Name,
% Description, Lon, Lat, and BoundaryBox fields.  Structs may contain a mix
% of points, lines, and polygons.
%
% .kml files with folder structure will not be presented as such, but will
% appear as a single vector array of structs.
%
%

[FID msg] = fopen(kmlFile,'rt');

if FID<0
    error(msg)
end

txt = fread(FID,'uint8=>char')';
fclose(FID);

expr = '<Placemark.+?>.+?</Placemark>';

objectStrings = regexp(txt,expr,'match');

Nos = length(objectStrings);

for ii = 1:Nos
    % Find Object Name Field
    bucket = regexp(objectStrings{ii},'<name.*?>.+?</name>','match');
    if isempty(bucket)
        name = 'undefined';
    else
        % Clip off flags
        name = regexprep(bucket{1},'<name.*?>\s*','');
        name = regexprep(name,'\s*</name>','');
    end

    % Find Object Description Field
    bucket = regexp(objectStrings{ii},'<description.*?>.+?</description>','match');
    if isempty(bucket)
        desc = '';
    else
        % Clip off flags
        desc = regexprep(bucket{1},'<description.*?>\s*','');
        desc = regexprep(desc,'\s*</description>','');
    end

    geom = 0;
    % Identify Object Type
    if ~isempty(regexp(objectStrings{ii},'<Point', 'once'))
        geom = 1;
    elseif ~isempty(regexp(objectStrings{ii},'<LineString', 'once'))
        geom = 2;
    elseif ~isempty(regexp(objectStrings{ii},'<Polygon', 'once'))
        geom = 3;
    end

    switch geom
        case 1
            geometry = 'Point';
        case 2
            geometry = 'Line';
        case 3
            geometry = 'Polygon';
        otherwise
            geometry = '';
    end

    % Find Coordinate Field
    bucket = regexp(objectStrings{ii},'<coordinates.*?>.+?</coordinates>','match');
    % Clip off flags
    coordStr = regexprep(bucket{1},'<coordinates.*?>(\s+)*','');
    coordStr = regexprep(coordStr,'(\s+)*</coordinates>','');
    % Split coordinate string by commas or white spaces, and convert string
    % to doubles
    coordMat = str2double(regexp(coordStr,'[,\s]+','split'));
    % Rearrange coordinates to form an x-by-3 matrix
    [m,n] = size(coordMat);
    coordMat = reshape(coordMat,3,m*n/3)';

    % define polygon in clockwise direction, and terminate
    [Lat, Lon] = poly2ccw(coordMat(:,2),coordMat(:,1));
    if geom==3
        Lon = [Lon;NaN];
        Lat = [Lat;NaN];
    end

    % Create structure
    kmlStruct(ii).Geometry = geometry;
    kmlStruct(ii).Name = name;
    kmlStruct(ii).Description = desc;
    kmlStruct(ii).Lon = Lon;
    kmlStruct(ii).Lat = Lat;
    kmlStruct(ii).BoundingBox = [[min(Lon) min(Lat);max(Lon) max(Lat)]];
end

author: Sable


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

function [plug] = getPlugNr(appliance, house, dataset)

    % returns the plug nr of a specified (appliance, house, dataset)-triple

    appliance_house_matrix = getApplianceHouseMatrix(dataset);
    plug = num2str(appliance_house_matrix(appliance, house), '%02d');

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [c]=findcntr(H)
    %
    %    [Center]=findcntr(Handle)
    %    Returns the geometric center of Handle
    %
    %    This can be used as the ORIGIN input in ROTATE to rotate
    %    an object about its geometric center.
    %
    %    See also ROTATE, TRNSLATE, SCALHAND
    %
    %    Richard G. Cobb    3/96
    %
    c=[0 0 0];

    tmp=(max(get(H,'Xdata'))+min(get(H,'Xdata')))/2;
    if tmp ~= []
        c(1)=tmp;
    else
        set(H,'Xdata',zeros(size(get(H,'Ydata'))))
    end
    tmp=(max(get(H,'Ydata'))+min(get(H,'Ydata')))/2;
    if tmp ~= []
        c(2)=tmp;
    else
        set(H,'Ydata',zeros(size(get(H,'Xdata'))))
    end

    tmp=(max(get(H,'Zdata'))+min(get(H,'Zdata')))/2;
    if tmp ~= []
        c(3)=tmp;
    else
        set(H,'Zdata',zeros(size(get(H,'Ydata'))))
    end

    %eof

author: CelsoReyes


text:
function moverotateobj(action)
%moverotateobj(action)
%
%Called when an object is clicked on
%whos ButtonDownFcn = moverotateobj('BtnDown')
%
%This is can be done by using
%   set(h,'ButtonDownFcn','moverotateobj(''BtnDown'')')
%where h is handle vector or typing
%moverotateobj('On') does this automatically
%moverotateobj('Off') takes it off
%
%moverotateobj by its self turns the property 'On'
%if it is 'Off' or turns it 'Off' if it is 'On'
%
%Copy-Left, Alejandro Sanchez-Barba, 2005

if nargin==0
    hFig = gcf;
    ud = get(hFig,'UserData');
    if isfield(ud,'manual')
        if strcmp(ud.manual,'On')
            action = 'Off';
        else %off
            action = 'On';
        end
    else
        action = 'On';
    end
end
eval(action)
return

%-----------------------------------------------------------
function BtnDown
hFig = gcf;
Button = get(hFig,'SelectionType');
ud = get(hFig,'UserData');
ud.Button = Button;
ud.hAx = get(hFig,'CurrentAxes');
ud.CurrentPoint = get(ud.hAx,'CurrentPoint');
ud.hObj = get(hFig,'CurrentObject');
set(ud.hObj,'Selected','On')
ud.ObjType = get(ud.hObj,'Type');
switch ud.Button
    case 'normal' %Move
        setptr(hFig,'closedhand');
    case 'alt' %Rotate
        setptr(hFig,'rotate');
        ud.hLine = line('XData',[],'YData',[], ...
            'Parent',ud.hAx,'Color','k');
        ud.Orig = ud.CurrentPoint;
    otherwise
        return
end %switch
ud.xmode = get(ud.hAx,'XLimMode');
ud.ymode = get(ud.hAx,'YLimMode');
set(ud.hAx,'XLimMode','manual')
set(ud.hAx,'YLimMode','manual')
set(hFig,'UserData',ud)
set(hFig,'WindowButtonMotionFcn','moverotateobj(''BtnMotion'')',...
   'WindowButtonUpFcn', 'moverotateobj(''BtnUp'')');
return

%-----------------------------------------------------------
function BtnMotion
hFig = gcf;
ud = get(hFig,'UserData');
hObj = ud.hObj;
ObjType = ud.ObjType;
switch ObjType
    case 'line'
        x = get(hObj,'XData');
        y = get(hObj,'YData');
    case 'text'
        pos = get(hObj,'Position');
        x = pos(1);
        y = pos(2);
    otherwise
        error('Unknown Object');
end %switch
hAx = ud.hAx;
LastPoint = ud.CurrentPoint;
CurrentPoint = get(hAx,'CurrentPoint');
Button = ud.Button;
switch Button
    case 'normal' %Move
        %Calculate Difference
        dx = CurrentPoint(1,1) - LastPoint(1,1);
        dy = CurrentPoint(1,2) - LastPoint(1,2);
        x = x + dx;
        y = y + dy;
        %Update position on Plot
        switch ObjType
            case 'line'
                set(hObj,'XData',x,'YData',y)
            case 'text'
                set(hObj,'Position',[x,y])
            otherwise
                error('Unknown Object')
        end %switch
    case 'alt' %Rotate
        xline = [ud.Orig(1,1),CurrentPoint(1,1)];
        yline = [ud.Orig(1,2),CurrentPoint(1,2)];
        set(ud.hLine,'XData',xline,'YData',yline)
        u1 = CurrentPoint(1,1) - ud.Orig(1,1);
        v1 = CurrentPoint(1,2) - ud.Orig(1,2);
        th1 = atan2(v1,u1);
        u2 = LastPoint(1,1) - ud.Orig(1,1);
        v2 = LastPoint(1,2) - ud.Orig(1,2);
        th2 = atan2(v2,u2);
        theta = th1 - th2;
%         disp([th1,th2,theta]*180/pi)
        switch ObjType
            case 'line'
                rot = [cos(theta) -sin(theta); ...
                    sin(theta) cos(theta)];
                vec = [x(:)-ud.Orig(1,1),y(:)-ud.Orig(1,2)];
                newvec = (rot*vec')';
                x = newvec(:,1) + ud.Orig(1,1);
                y = newvec(:,2) + ud.Orig(1,2);
                set(hObj,'XData',x,'YData',y)
            case 'text'
                thetaold = get(hObj,'Rotation');
                thetanew = thetaold + theta*180/pi;
                set(hObj,'Rotation',thetanew)
            otherwise
                error('Unknown Object')
        end %switch
end %switch
ud.CurrentPoint = CurrentPoint;
set(hFig,'UserData',ud)
return

%-----------------------------------------------------------
function BtnUp
hFig = gcf;
setptr(hFig,'arrow');
set(hFig,'WindowButtonMotionFcn','');
ud = get(hFig,'UserData');
ud.xmode = get(ud.hAx,'XLimMode');
ud.ymode = get(ud.hAx,'YLimMode');
set(ud.hObj,'Selected','Off')
if strcmp(ud.xmode,'auto')
    set(ud.hAx,'XLimMode','auto')
end
if strcmp(ud.ymode,'auto')
    set(ud.hAx,'YLimMode','auto')
end
if strcmp(ud.Button,'alt')
    set(ud.hLine,'XData',[],'YData',[])
end
return

%-----------------------------------------------------------
function On
hFig = gcf;
ud = get(hFig,'UserData');
ud.hAx = get(hFig,'CurrentAxes');
ud.manual = 'On';
set(hFig,'UserData',ud)
h = get(ud.hAx,'Children');
set(h,'ButtonDownFcn','moverotateobj(''BtnDown'')')
return

%-----------------------------------------------------------
function Off
hFig = gcf;
ud = get(hFig,'UserData');
ud.hAx = get(hFig,'CurrentAxes');
ud.manual = 'Off';
set(hFig,'UserData',ud)
h = get(ud.hAx,'Children');
set(h,'ButtonDownFcn','')
set(hFig,'WindowButtonMotionFcn','',...
            'WindowButtonUpFcn','');
return


text:
function varargout = optimN(varargin)
% Full multigrid matrix solver stuff (circulant boundaries)
%_______________________________________________________________________
%
% FORMAT v = optimN('fmg',A, b, param)
% v     - the solution n1*n2*n3*n4
% A     - parameterisation of 2nd derivatives
%         n1*n2*n3*(n4*(n4+1)/2)
%         The first n4 volumes are the diagonal elements, which are
%         followed by the off-diagonals (note that 2nd derivs are%
%         symmetric).  e.g. if n4=3, then the ordering would be
%         (1,1),(2,2),(3,3),(1,2),(1,3),(2,3)
% b     - parameterisation of first derivatives n1*n2*n3*n4
% param - 6 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Voxel sizes
%         - [5][6][7] Regularisation parameters
%           - For membrane and bending energy, the parameters
%             are lambda, unused and id.
%         - [8] Number of Full Multigrid cycles
%         - [9] Number of relaxation iterations per cycle
%
%           Note that more cycles and iterations may be needed
%           for bending energy than for membrane energy.
%
% Solve equations using a Full Multigrid method.  See Press et al
% for more information.
% v = inv(A+H)*b
% A, b and v are all single precision floating point.
% H is a large sparse matrix encoded by param(1:7).
% The tensor field encoded by A MUST be positive-definite. If it is not,
% then anything could happen (see references about Fisher scoring for
% help on ensuring that second derivatives are positive definite).
%
%_______________________________________________________________________
%
% FORMAT m = optimN('vel2mom', v, param)
% v     - velocity (flow) field n1*n2*n3*n4.
% param - 4 parameters (settings)
%         - [1] Regularisation type, can take values of
%           - 1 Membrane energy
%           - 2 Bending energy
%         - [2][3][4] Voxel sizes
%         - [5][6][7] Regularisation parameters
%           - For membrane and bending energy, the parameters
%             are lambda, unusaed and id.
% m       - `momentum' field n1*n2*n3*n4.
%
% Convert a flow field to a momentum field by m = H*v, where
% H is the large sparse matrix encoding some form of regularisation.
% v and m are single precision floating point. This function has uses
% beyond only image registration.
%
%_______________________________________________________________________
%
% Note that the boundary conditions are circulant throughout.
% For Neumann boundary conditions (zero gradients at the boundaries)
% use optimNn.
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: optimN.m 4758 2012-05-29 15:34:08Z john $

[varargout{1:nargout}] = optim_compat(0,varargin{:});

%error('Not compiled for %s in MATLAB %s  (see make.m)\n', computer, version);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function x = cs_utsolve (U,b)                                               %#ok
%CS_UTSOLVE solve a sparse lower triangular system U'*x=b.
%   x = cs_utsolve(U,b) computes x = U'\b, U must be upper triangular with a
%   zero-free diagonal.  b must be a full vector.
%
%   Example:
%       Prob = UFget ('HB/arc130') ; A = Prob.A ; n = size (A,1) ;
%       b = rand (n,1);
%       [L U p q] = cs_lu (A) ;
%       x = cs_ltsolve (L, cs_utsolve (U, b(q))) ;   % x = L' \ (U' \ b(q)) ;
%       x (p) = x ;
%       norm (A'*x-b)
%
%   See also CS_LSOLVE, CS_LTSOLVE, CS_USOLVE, MLDIVIDE.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_utsolve mexFunction not found') ;


text:
function [SE3] = SE3MatrixFromComponents( x, y, z, r, p, yaw )

% SE3MatrixFromComponents - build a 4x4 matrix representing an SE(3) transform
%
% [SE3] = SE3MatrixFromComponents( x, y, z, r, p, yaw )
%
% INPUTS:
%   x, y, z: translation
%   r, p, yaw: rotation in Euler angle representation
%
% OUTPUTS:
%   SE3: 4x4 matrix representing the SE(3) transform

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2016 University of Oxford
% Authors:
%  Geoff Pascoe (gmp@robots.ox.ac.uk)
%  Will Maddern (wm@robots.ox.ac.uk)
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % Allow passing of a single 6-element vector
  if nargin == 1
    y = x(2);
    z = x(3);
    r = x(4);
    p = x(5);
    yaw = x(6);
    x = x(1);
  end

  % Convert euler angles to rotation matrices
  R_x = [
    1, 0, 0;
    0, cos(r), -sin(r);
    0, sin(r), cos(r) ];

  R_y = [
    cos(p), 0, sin(p);
    0, 1, 0;
    -sin(p), 0, cos(p) ];

  R_z = [
    cos(yaw), -sin(yaw), 0;
    sin(yaw), cos(yaw), 0;
    0, 0, 1];

  R = R_z * R_y * R_x;

  SE3 = [R [x; y; z]; zeros(1,3), 1];

end

author: ori-mrg

repo: robotcar-dataset-sdk


text:
function out = chebfunroot()
%CHEBFUNROOT   Root directory of Chebfun installation.
%   S = CHEBFUNROOT() returns a string that is the name of the directory where
%   Chebfun is installed.

out = fileparts(which('chebfunroot'));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [p,t] = circlemesh(x,y,r,h)
%% CIRCLEMESH generates a mesh for a circle
%
% [node,elem] = circlemesh(x,y,r,h) generates a quasi-uniform mesh with
% size h of the circle centered at (x,y) with readius r.
%
% Example
%
%   [node,elem] = circlemesh(0,0,1,0.1);
%   showmesh(node,elem);
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

[p,t] = odtmesh2d(@fd,@huniform,h,[-1,-1;1,1],[],0,x,y,r);

    function s = fd(p,x,y,r)
    s = sqrt(sum((p(:,1)-x).^2+(p(:,2)-y).^2,2))-r;
    end

    function h = huniform(p,varargin)
    h = ones(size(p,1),1);
    end
end

repo: ifem


text:
function varargout = mat2xml(matfile, xmlfile)
%MAT2XML Convert a MAT-file into an XML file
%  MAT2XML(MATFILE, XMLFILE) loads MAT-file MATFILE and converts it
%  into an XML file saved in XMLFILE.
%  T = MAT2XML(MATFILE) also converts the MAT-file MATFILE to XML
%  but returns it instead as an XMLTree object T.
%
%  See also XML2MAT, XMLTREE, SAVE.

%  Copyright 2003 Guillaume Flandin.
%  Revision: 1.0 $  $Date: 2003/07/01 21:43 $

error(nargchk(1,2,nargin));

s = load(matfile);
names = fieldnames(s);

t = xmltree(['<matfile name="' matfile '"/>']);

for i=1:length(names)
	[t, uid] = add(t,root(t),'element',names{i});
	t = attributes(t,'add',uid,'type',class(getfield(s,names{i})));
	t = attributes(t,'add',uid,'size',sub_num2str(size(getfield(s,names{i}))));
	t = sub_var2xml(t,getfield(s,names{i}),uid);
end

if nargin == 2
	save(t,xmlfile);
else
	varargout{1} = t;
end

%=======================================================================
function t = sub_var2xml(t,v,uid)

	switch class(v)
		case 'double'
			t = add(t,uid,'chardata',sub_num2str(v));
		case 'sparse'
			% TODO % better names for sparse elements
			[i,j,s] = find(v);
			[t, uid2] = add(t,uid,'element','row');
			t = attributes(t,'add',uid2,'size',sub_num2str(size(i)));
			t = add(t,uid2,'chardata',sub_num2str(i));
			[t, uid2] = add(t,uid,'element','col');
			t = attributes(t,'add',uid2,'size',sub_num2str(size(j)));
			t = add(t,uid2,'chardata',sub_num2str(j));
			[t, uid2] = add(t,uid,'element','val');
			t = attributes(t,'add',uid2,'size',sub_num2str(size(s)));
			t = add(t,uid2,'chardata',sub_num2str(s));
		case 'struct'
			names = fieldnames(v);
			for j=1:prod(size(v))
				for i=1:length(names)
					[t, uid2] = add(t,uid,'element',names{i});
					t = attributes(t,'add',uid2,'index',num2str(j));
					t = attributes(t,'add',uid2,'type',...
						class(getfield(v(j),names{i})));
					t = attributes(t,'add',uid2,'size', ...
						sub_num2str(size(getfield(v(j),names{i}))));
					t = sub_var2xml(t,getfield(v(j),names{i}),uid2);
				end
			end
		case 'cell'
			for i=1:prod(size(v))
				[t, uid2] = add(t,uid,'element','cell');
				% TODO % special handling of cellstr ?
				t = attributes(t,'add',uid2,'index',num2str(i));
				t = attributes(t,'add',uid2,'type',class(v{i}));
				t = attributes(t,'add',uid2,'size',sub_num2str(size(v{i})));
				t = sub_var2xml(t,v{i},uid2);
			end
		case 'char'
			% TODO % char values should be in CData
			t = add(t,uid,'chardata',v);
		case {'int8','uint8','int16','uint16','int32','uint32'}
			[t, uid] = add(t,uid,'element',class(v));
			% TODO % Handle integer formats (cannot use sprintf or num2str)
		otherwise
			if ismember('serialize',methods(class(v)))
				% TODO % is CData necessary for class output ?
				t = add(t,uid,'cdata',serialize(v));
			else
				warning(sprintf(...
				'[MAT2XML] Cannot convert from %s to XML.',class(v)));
			end
	end

%=======================================================================
function s = sub_num2str(n)
	% TODO % use format ?
	if isempty(n)
		s = '[]';
	else
		s = ['[' sprintf('%g ',n(1:end-1))];
		s = [s num2str(n(end)) ']'];
	end

author: buzsakilab

repo: buzcode


text:
function z = classifierHead(x,poolerWeights,classifierWeights)
% classifierHead   The standard classification head for a BERT model.
%
%   Z = classifierHead(X,poolerWeights,classifierWeights) applies
%   bert.layer.pooler and bert.layer.classifier to X with poolerWeights and
%   classifierWeights respectively. Both poolerWeights and
%   classifierWeights must be structs with fields 'kernel' and 'bias'.

% Copyright 2021 The MathWorks, Inc.
z = bert.layer.pooler(x,poolerWeights);
z = bert.layer.classifier(z,classifierWeights);
end

author: matlab-deep-learning

repo: transformer-models


text:
%%% [video_aux] = crop_4dvideo(video,r,s)
%
function [video_aux] = crop_4dvideo(video,r,s)
  if(s == -1)
    s = size(video,4);
  end

  k = 1;
  for i = r : s
    video_aux(:,:,:,k) = video(:,:,:,i);
    k = k + 1;
  end
end

author: andrewssobral

repo: lrslibrary


text:
function [ output_args ] = interfaceError( fname,fid )
%

if ~exist(fname,'file')
  e = MException('mtex:fileNotFound',...
    ['File ''' fname ''' not found']);
  throwAsCaller(e);
end

[st] = dbstack(1);
i = find(strncmp('load',{st.name},4),1,'first');
name = st(i).name;

type = regexp(name,'(?<=load)(\w*)(?=_)','match');
interface = regexp(lower(name),'(?<=[a-z]*_)\w*','match');

if nargin > 1
  try
    fclose(fid);
  catch
  end
end

fname = regexprep(fname,'\','/');

e = MException('mtex:wrongInterface',...
  [type{1} ' format ''' upper(interface{1}) ''' does not match the data\n file: ''' fname '''']);
throwAsCaller(e);

author: mtex-toolbox

repo: mtex

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


text:
classdef Two_Arch2 < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Two-archive algorithm 2
% CAsize --- --- Convergence archive size
% p      --- --- The parameter of fractional distance

%------------------------------- Reference --------------------------------
% H. Wang, L. Jiao, and X. Yao, Two_Arch2: An improved two-archive
% algorithm for many-objective optimization, IEEE Transactions on
% Evolutionary Computation, 2015, 19(4): 524-541.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [CAsize,p] = Algorithm.ParameterSet(Problem.N,1/Problem.M);

            %% Generate random population
            Population = Problem.Initialization();
            CA = UpdateCA([],Population,CAsize);
            DA = UpdateDA([],Population,Problem.N,p);

            %% Optimization
            while Algorithm.NotTerminated(DA)
                [ParentC,ParentM] = MatingSelection(CA,DA,Problem.N);
                Offspring         = [OperatorGA(Problem,ParentC,{1,20,0,0}),OperatorGA(Problem,ParentM,{0,0,1,20})];
                CA = UpdateCA(CA,Offspring,CAsize);
                DA = UpdateDA(DA,Offspring,Problem.N,p);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% SPARFWSLV Solves block sparse upper-triangular system.
%    y = sparfwslv(L,b) yields the same result as
%              y = L.L\b(L.perm,:)
%    However, SPARFWSLV is faster than the built-in operator "\",
%    because it uses dense linear algebra and loop-unrolling on
%    supernodes.
%
%    For sparse b, one should use
%    y = sparfwslv(L,b,symbfwblk(L.L,L.xsuper, b));
%
%    Typical use, with X sparse m x m positive definite and b is m x n:
%            L = sparchol(symbchol(X),X);
%            L.d(L.dep) = inf;
%            y = sparbwslv(L,sparfwslv(L,b) ./ L.d);
%    Then y solves X*y=b.
%
% See also symbchol, sparchol, sparbwslv, mrdivide, mldivide.

function y = sparfwslv(L,b, ysymb)
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

% ----------------------------------------
% Solve L.L * y = b
% ----------------------------------------
if nargin > 2
    y = fwblkslv(L,b,ysymb);
else
    y = fwblkslv(L,b);
end

repo: REST


text:
function writeCbModel(model,format,fileName,compSymbolList,compNameList,sbmlLevel,sbmlVersion)
%writeCbModel Write out COBRA models in various formats
%
% writeCbModel(model,format,fileName,compSymbolList,compNameList,sbmlLevel,sbmlVersion)
%
%INPUTS
% model             Standard COBRA model structure
% format            File format to be used ('text','xls' or 'sbml')
%
%OPTIONAL INPUTS
% fileName          File name for output file (optional, default opens
%                   dialog box)
% compSymbolList    List of compartment symbols
% compNameList      List of compartment names corresponding to compSymbolList
% sbmlLevel         SBML Level (default = 2)
% sbmlVersion       SBML Version (default = 1)
%
% Markus Herrgard 2/5/07
% Ines Thiele 01/10 - Added more options for field to write in xls format
% Richard Que (3/17/10) Added ability to specify compartment names and
%                       symbols

if ~exist('compSymbolList','var') || isempty(compSymbolList)
    compSymbolList = {'c','m','v','x','e','t','g','r','n','p','l','y'};
    compNameList = {'Cytoplasm','Mitochondrion','Vacuole','Peroxisome','Extracellular','Pool','Golgi','Endoplasmic_reticulum','Nucleus','Periplasm','Lysosome','Glycosome'};
end

if nargin < 6
    sbmlLevel = 2;
    sbmlVersion = 1;
end

[nMets,nRxns] = size(model.S);

formulas = printRxnFormula(model,model.rxns,false,false,false,1,false);

%% Open a dialog to select file name
if (nargin < 3 & ~strcmp(format,'sbml'))
    switch format
        case 'xls'
            [fileNameFull,filePath] = uiputfile({'*.xls'});
        case {'text','txt'}
            [fileNameFull,filePath] = uiputfile({'*.txt'});
        case 'xml'
            [fileNameFull,filePath] = uiputfile({'*.xml'});
        otherwise
            [fileNameFull,filePath] = uiputfile({'*'});
    end
    if (fileNameFull)
        [t1,t2,t3,t4,tokens] = regexp(fileNameFull,'(\w*)\.(\w*)');
        fileName = [filePath tokens{1}{1}];
        switch tokens{1}{2}
            case 'xls'
                format = 'xls';
            case 'txt'
                format = 'text';
                fileName = [fileName '.txt'];
            case 'xml'
                format = 'sbml';
%                 fprintf('Note that you will be asked to supply the file name again (this is a feature, not a bug)');
            otherwise
                format = 'unknown';
        end
    else
        return;
    end
end
switch format
    %% Text file
    case {'text','txt'}
        fid = fopen(fileName,'w');
        fprintf(fid,'Rxn name\t');
        if (isfield(model,'rxnNames'))
            fprintf(fid,'Rxn description\t');
        end
        fprintf(fid,'Formula\t');
        if (isfield(model,'grRules'))
            fprintf(fid,'Gene-reaction association\t');
        end
        fprintf(fid,'Reversible\tLB\tUB\tObjective\n');
        for i = 1:nRxns
            fprintf(fid,'%s\t',model.rxns{i});
            if (isfield(model,'rxnNames'))
                fprintf(fid,'%s\t',model.rxnNames{i});
            end
            fprintf(fid,'%s\t',formulas{i});
            if (isfield(model,'grRules'))
                fprintf(fid,'%s\t',model.grRules{i});
            end
            fprintf(fid,'%d\t%6.2f\t%6.2f\t%6.2f\n',model.rev(i),model.lb(i),model.ub(i),model.c(i));
        end
        fprintf(fid,'Metabolite name\tMetabolite description\tMetabolite formula\n');
        for i = 1:nMets
            fprintf(fid,'%s',model.mets{i});
            if isfield(model,'metNames')
                fprintf(fid,'\t%s',model.metNames{i});
            end
            if isfield(model,'metFormulas')
                fprintf(fid,'\t%s',model.metFormulas{i});
            end
            fprintf(fid,'\n');
        end
        fclose(fid);
        %% Excel file
    case 'xls'
        tmpData{1,1} = 'Abbreviation';
        tmpData{1,2} = 'Description';
        baseInd = 3;
        tmpData{1,baseInd} = 'Reaction';
        tmpData{1,baseInd+1} = 'GPR';
        tmpData{1,baseInd+2} = 'Genes';
        tmpData{1,baseInd+3} = 'Proteins';
        tmpData{1,baseInd+4} = 'Subsystem';
        tmpData{1,baseInd+5} = 'Reversible';
        tmpData{1,baseInd+6} = 'Lower bound';
        tmpData{1,baseInd+7} = 'Upper bound';
        tmpData{1,baseInd+8} = 'Objective';
        tmpData{1,baseInd+9} = 'Confidence Score';
        tmpData{1,baseInd+10} = 'EC Number';
        tmpData{1,baseInd+11} = 'Notes';
        tmpData{1,baseInd+12} = 'References';
        for i = 1:nRxns
            tmpData{i+1,1} = chopForExcel(model.rxns{i});
            if (isfield(model,'rxnNames'))
                tmpData{i+1,2} = chopForExcel(model.rxnNames{i});
            else
                tmpData{i+1,2} =  '';
            end

            tmpData{i+1,baseInd} = chopForExcel(formulas{i});
            if (isfield(model,'geneNameRules'))
                tmpData{i+1,baseInd+1} = chopForExcel(model.geneNameRules{i});
            elseif (isfield(model,'grRules'))
                tmpData{i+1,baseInd+1} = chopForExcel(model.grRules{i});
            else
                tmpData{i+1,baseInd+1} = '';
            end
            if (isfield(model,'geneNames'))
                geneNames = model.geneNames(model.rxnGeneMat(i,:) == 1);
                tmpData{i+1,baseInd+2} = constructGeneStr(geneNames);
            elseif (isfield(model,'genes'))
                geneNames = model.genes(model.rxnGeneMat(i,:) == 1);
                tmpData{i+1,baseInd+2} = constructGeneStr(geneNames);
            else
                tmpData{i+1,baseInd+2} = '';
            end
            if (isfield(model,'proteins'))
                tmpData{i+1,baseInd+3} = chopForExcel(model.proteins{i});
            else
                tmpData{i+1,baseInd+3} = '';
            end
            if (isfield(model,'subSystems'))
                tmpData{i+1,baseInd+4} = chopForExcel(char(model.subSystems{i}));
            else
                tmpData{i+1,baseInd+4} = '';
            end
            tmpData{i+1,baseInd+5} = model.rev(i)*1.0;
            tmpData{i+1,baseInd+6} = model.lb(i);
            tmpData{i+1,baseInd+7} = model.ub(i);
            tmpData{i+1,baseInd+8} = model.c(i);
            if (isfield(model,'confidenceScores'))
                tmpData{i+1,baseInd+9} =  chopForExcel(num2str(model.confidenceScores{i}));
            else
                tmpData{i+1,baseInd+9} = '';
            end
            if (isfield(model,'rxnECNumbers'))
                tmpData{i+1,baseInd+10} = chopForExcel(model.rxnECNumbers{i});
            else
                tmpData{i+1,baseInd+10} = '';
            end
            if (isfield(model,'rxnNotes'))
                tmpData{i+1,baseInd+11} = chopForExcel(char(model.rxnNotes{i}));
            else
                tmpData{i+1,baseInd+11} = '';
            end
            if (isfield(model,'rxnReferences'))
                tmpData{i+1,baseInd+12} = chopForExcel(char(model.rxnReferences{i}));
            else
                tmpData{i+1,baseInd+12} = '';
            end
        end
        %keyboard
        xlswrite(fileName,tmpData,'Reaction List');
        if isfield(model,'metNames')
            tmpMetData{1,1} = 'Abbreviation';
            tmpMetData{1,2} = 'Description';
            tmpMetData{1,3} = 'Neutral formula';
            tmpMetData{1,4} = 'Charged formula';
            tmpMetData{1,5} = 'Charge';
            tmpMetData{1,6} = 'Compartment';
            tmpMetData{1,7} = 'KEGG ID';
            tmpMetData{1,8} = 'PubChem ID';
            tmpMetData{1,9} = 'ChEBI ID';
            tmpMetData{1,10} = 'InChI String';
            tmpMetData{1,11} = 'SMILES';
            for i = 1:nMets
                tmpMetData{i+1,1} = chopForExcel(model.mets{i});
                tmpMetData{i+1,2} = chopForExcel(model.metNames{i});
                if isfield(model,'metFormulasNeutral')
                    tmpMetData{i+1,3} = chopForExcel(model.metFormulasNeutral{i});
                else
                    tmpMetData{i+1,3} = '';
                end
                if isfield(model,'metFormulas')
                    tmpMetData{i+1,4} = chopForExcel(model.metFormulas{i});
                else
                    tmpMetData{i+1,4} = '';
                end
                if isfield(model,'metCharge')
                    tmpMetData{i+1,5} = chopForExcel(model.metCharge(i));
                else
                    tmpMetData{i+1,5} = '';
                end
                if isfield(model,'metCompartment')
                    tmpMetData{i+1,6} = chopForExcel(model.metCompartment{i});
                else
                    tmpMetData{i+1,6} = '';
                end
                if isfield(model,'metKEGGID')
                    tmpMetData{i+1,7} = chopForExcel(model.metKEGGID{i});
                else
                    tmpMetData{i+1,7} = '';
                end
                if isfield(model,'metPubChemID')
                    if iscell(model.metPubChemID(i))
                    tmpMetData{i+1,8} = chopForExcel(model.metPubChemID{i});
                    else
                    tmpMetData{i+1,8} = chopForExcel(model.metPubChemID(i));
                    end
                else
                    tmpMetData{i+1,8} = '';
                end
                if isfield(model,'metChEBIID')

                    tmpMetData{i+1,9} = chopForExcel(model.metChEBIID(i));
                else
                    tmpMetData{i+1,9} = '';
                end
                if isfield(model,'metInchiString')
                    tmpMetData{i+1,10} = chopForExcel(model.metInchiString{i});
                else
                    tmpMetData{i+1,10} = '';
                end
                if isfield(model,'metSmiles')
                    tmpMetData{i+1,11} = chopForExcel(model.metSmiles{i});
                else
                    tmpMetData{i+1,11} = '';
                end
            end
            xlswrite(fileName,tmpMetData,'Metabolite List');
        else
            xlswrite(fileName,model.mets,'Metabolite List');
        end
        %% SBML
    case 'sbml'
        sbmlModel = convertCobraToSBML(model,sbmlLevel,sbmlVersion,compSymbolList,compNameList);
        if exist('fileName','var')&&~isempty(fileName)
            OutputSBML(sbmlModel,fileName);
        else
            OutputSBML(sbmlModel);
        end
        %% Unknown
    otherwise
        error('Unknown file format');
end

%% Chop strings for excel output
function strOut = chopForExcel(str)

if (length(str) > 5000)
    strOut = str(1:5000);
    fprintf('String longer than 5000 characters - truncated for Excel output\n%s\n',str);
else
    strOut = str;
end

%% Construct gene name string
function geneStr = constructGeneStr(geneNames)

geneStr = '';
for i = 1:length(geneNames)
    geneStr = [geneStr ' ' geneNames{i}];
end
geneStr = strtrim(geneStr);

author: opencobra

repo: cobratoolbox