text:
function [y, Fout]=SincResample(xin, npoints, Fin, shift, mode)
% SincResample resamples signals to arbitrary lengths/frequencies and can
% timeshift them.
%
% SincResample can be faster than INTERP or RESAMPLE when X contains many
% relatively short data columns.
%
%
% Example:
% [Y, FOUT]=SincResample(X, NPOINTS, FIN, SHIFT, MODE)
%
% X is a vector or M x N  matrix where each column represents a signal.
% NPOINTS is the number of points required in each column after
% resampling.
% FIN (if supplied) is the sampling rate in X (default 1 samples/sec).
% SHIFT (if supplied) is the time (in seconds) to shift the points at which
% interpolation will be done (default 0). SHIFT should be between ą(1/Fin).
% SHIFT can be used to synchronize signals separated by up to ą(1/Fin)
% in time.
% MODE (if supplied) is a string:
%           'sinc' to use a hamming-windowed sinc function (default)
%           'lanczos' for a 3-lobed Lanczos function
%
% Y is the resulting M x NPOINTS matrix formed by resampling. The time
% of the first sample in each column of Y is shifted by SHIFT seconds
% relative to those in X.
% FOUT is the effective sample rate for output Y (normalised to FIN=1 if
% FIN is not supplied on input).
%
% If X is a bandlimited (filtered) signal containing components at
% DC - 0.5Fs where Fs is the sample rate, SincResample will return
% the signal that would have been seen with a higher sampling rate
% (and the same filter settings). The algorithm depends on the fact that
% the bandlimited signal is completely defined when sampling at Fs.
%
% SincResample returns the data convolved with a set of  time-shifted
% windowed sinc functions, one for each of the samples [1..size(x,1)]
% in the input signal. End-effects are reduced by adding a reflected
% and mirrored copy of the data at the beginning and end of each column
% before resampling.
%
% Note that points at the end of Y (and/or beginning if SHIFT is negative)
% will be extrapolated beyond the boundaries of X.
%
% SincResample is a generalization of the method used in Blanche &
% Swindale (2006).
%
% See also e.g. interp & resample (in the Signal Processing Toolbox) and
% the sinc online documentation.
%
% References:
% T.Blanche & N.V.Swindale (2006)
% Nyquist interpolation improves neuron yield in multiunit recordings
% Journal of Neuroscience Methods 155, 81-91 LinkOut
%
% also
%
% K. Turkowski, Filters for common resampling tasks LinkOut
%
% Toolboxes required: None
%
% Author: Malcolm Lidierth 07/06
% Copyright Š Kings College London 2006-8
%
% Acknowledgements: Steve Alty for comments
%
% Revisions: 10.07  Argument checks tidied (ML)
%            02.08  End-effect minimization now relects and mirrors data
%                   instead of simply reflecting it. This improves the
%                   slopes at the end and reduces the behaviour described at
%                   http://www.dsprelated.com/showmessage/79676/2.php
%                   (ML)
%            02.08  Mode option added (ML)

% Argument checks
if nargin<2
    y=[];
    Fout=0;
    return
end

if nargin<=2
    % Default
    Fin=1;
end

if nargin<=3
    % Zero shift default
    shift=0;
end

if nargin<=4
    % Set sinc as default
    mode='sinc';
end

if Fin==0
    error('SincResample: Fin specified as zero');
end

if abs(shift)>=1/Fin
    error('SincResample: shift must be less than 1/Fin');
end

if rem(npoints,1)~=0
    error('SincResample: npoints must be a whole number');
end

% Main routine
% Prepend and append a reflected and mirrored copy of the data
% to minimize end transients:
% Change 22.02.08: note change from reflected only

if isvector(xin)
    % Force column vector if row input
    xin=xin(:);
end

x=zeros(size(xin,1)*3, size(xin,2));
for k=1:size(xin,2)
    x(2:end-1,k)=[2*xin(1,k)-xin(end:-1:2,k);...
        xin(:,k);...
        2*xin(end,k)-xin(end-1:-1:1,k)];
    x(1,k)=2*x(1,k)-x(2,k);
    x(end,k)=2*x(end,k)-x(end-1,k);
end
clear('xin');
np=npoints*3;

% A column vector of sample times for x  - assume unit spacing and zero
% base
len=size(x,1);
t=(0:len-1)';

% times at which to interpolate data - based on spacing in t.
% ts will be a [npoints by size(x,1)] matrix
ts=(0:(max(t)+1)/np:(np-1)/(np/len))';
ts=(ts(:,ones(size(t)))-t(:,ones(size(ts)))');

% Hamming window
N=size(ts,2);
if N==1 || strcmp(mode, 'lanczos');
    w=1;
else
    th=ts+N-1;
    w=0.54-0.46*cos((2*pi*th/max(th(:)))) ;
end

% Add shift in multiples of a sampling interval i.e. shift/(1/Fin)
ts=ts+(shift*Fin);

switch mode
    case 'lanczos'
        % Generate Lanczos functions and calculate output
        y=lanczos(ts)*x;
    case 'sinc'
        % Generate sinc functions, apply Hamming window if w~=1 and
        % calculate output.
        y=sinc(ts).*w*x;
    otherwise
        error('mode should be ''sinc'' or ''lanczos''');
end

% Throw away what is not needed
y=y(size(y,1)/3+1:(size(y,1)/3)+(size(y,1)/3),:);
Fout=Fin*npoints/(len/3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEBUG Display the data - UNCOMMENT text below to display results
% x=x(size(x,1)/3+1:(size(x,1)/3)+(size(x,1)/3),:);
% debugdisplay(x, y, Fin, Fout, shift)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=sinc(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=find(x==0);
x(i)=1;
y=sin(pi*x)./(pi*x);
y(i)=1;
return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=lanczos(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=3;
y=sinc(x).*sinc(x/n);
y(abs(x)>=n)=0;
return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function debugdisplay(x, y, Fin, Fout, shift) %#ok<DEFNU>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interp requires the Signal Processing Toolbox. The output of interp will
% be shown if Fout is an integer multiple of Fin.
% Note that shift is not applied to interp
% Sample times in seconds
t1= 0 : 1/Fin : (size(x,1)-1)/Fin;
t2= 0 : 1/Fout : (size(y,1)-1)/Fout;
if rem(Fout,Fin)==0
    y3=interp(x(:,1),Fout/Fin);
    t3=linspace(0,max(t2),size(y3,1));
end
t2=t2+shift;
handle=gca(figure);
xlabel(handle,'Time (seconds)');
for i=1:size(x,2)
    if ~ishandle(handle)
        break;
    end
    if rem(Fout,Fin)==0
        y3=interp(x(:,i),Fout/Fin);
        plot(handle,t1,x(:,i),'-bo',t2,y(:,i),'r-x',t3,y3,'g:s');
        xlim([min([0;min(t2)]) max([max(t3);max(t2);max(t1)])]);
    else
        plot(handle,t1,x(:,i),'-bo',t2,y(:,i),'r-x');
        xlim([min([0;min(t2)]) max([max(t2);max(t1)])]);
    end
    xlabel(handle,'Time (seconds)');
    title('Close the figure window to exit this loop');
    if rem(Fout,Fin)==0
        legend('Original','SincResample','MATLAB Interp (no shift)');
    else
        legend('Original','SincResample');
    end
    drawnow();
    pause(.1);
end
return
end


text:
 function cells = blockify_data(ob, data, varargin)
%function cells = blockify_data(ob, data, varargin)

if isempty(ob.handle_blockify_data)
	error(['The Fatrix of type ' ob.caller ' has no handle_blockify_data'])
end
cells = ob.handle_blockify_data(ob, data, varargin{:});

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function hRD = plugin_rd_rx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9371x.common.plugin_rd('ZCU102', 'Rx');
AnalogDevices.adrv9371x.zcu102.add_io(hRD, 'Rx');

author: analogdevicesinc


text:
% ----------------------------------------------------------------------------
% function hfssExportToFile(fid, ReportName, FileName, Type)
%
% Description :
% -------------
% Exports the report data to file, so it can be postprocesed.
%
% Parameters :
% ------------
% fid        - file identifier of the HFSS script file.
% ReportName - name of the report to export.
% FileName   - name of the file.
% Type       - There are four possible types:
%               * txt: Post processor format file
%               * csv: Comma-delimited data file
%               * tab: Tab-separated file
%               * dat: Ansoft plot data file
%
% Note :
% ------
% In HFSS's scripting help (scripting.pdf), the example given is:
%   oDesign.ExportToFile "Plot1", "c:\report1.dat"
% However, it is wrong, since it should use oModule, instead of
% oDesign.
%
% Example :
% ---------
% hfssExportToFile(fid, 'Plot1', 'MyData', 'csv'); % Saves in the same dir.
% hfssExportToFile(fid, 'Plot2', 'C:\MyData2', 'tab');
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% 29-Sept-2012: *Initial release.
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Written by Daniel R. Prado
% danysan@gmail.com / drprado@tsc.uniovi.es
% 29 September 2012
% ----------------------------------------------------------------------------
function hfssExportToFile(fid, ReportName, FileName, Type)

% Arguments processor.
if (nargin < 4)
	error('Insufficient # of arguments !');
end

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'Set oModule = oDesign.GetModule("ReportSetup")\n');

% Parameters.
fprintf(fid, 'oModule.ExportToFile "%s", "%s.%s"\n', ReportName, FileName, Type);

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function x = typeIV(alpha)
%4 times of Bit-wise MAP
alpha_1 = alpha(1 : end/4);
% llr_1 = 2*atanh(prod(tanh(alpha_1/2)));
llr_1 = prod(sign(alpha_1))*min(abs(alpha_1));

alpha_2 = alpha(end/4 + 1 : end/2);
% llr_2 = 2*atanh(prod(tanh(alpha_2/2)));
llr_2 = prod(sign(alpha_2))*min(abs(alpha_2));

alpha_3 = alpha(end/2 + 1 : end*3/4);
% llr_3 = 2*atanh(prod(tanh(alpha_3/2)));
llr_3 = prod(sign(alpha_3))*min(abs(alpha_3));

alpha_4 = alpha(end*3/4 + 1 : end);
% llr_4 = 2*atanh(prod(tanh(alpha_4/2)));
llr_4 = prod(sign(alpha_4))*min(abs(alpha_4));

%Even Parity check
check_bit = (llr_1 + llr_2 + llr_3 + llr_4) < 0;
%Wagner Decoder
x1 = alpha_1 < 0;
if mod(sum(x1), 2) ~= check_bit
    x1(abs(alpha_1) == min(abs(alpha_1))) = mod(x1(abs(alpha_1) == min(abs(alpha_1))) + 1, 2);
end

x2 = alpha_2 < 0;
if mod(sum(x2), 2) ~= check_bit
    x2(abs(alpha_2) == min(abs(alpha_2))) = mod(x2(abs(alpha_2) == min(abs(alpha_2))) + 1, 2);
end

x3 = alpha_3 < 0;
if mod(sum(x3), 2) ~= check_bit
    x3(abs(alpha_3) == min(abs(alpha_3))) = mod(x3(abs(alpha_3) == min(abs(alpha_3))) + 1, 2);
end

x4 = alpha_4 < 0;
if mod(sum(x4), 2) ~= check_bit
    x4(abs(alpha_4) == min(abs(alpha_4))) = mod(x4(abs(alpha_4) == min(abs(alpha_4))) + 1, 2);
end

x_tmp = [x1; x2; x3; x4];

x = x_tmp(:);

end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/mmv_phase_transition_snr_40db_s_32.mat';
N = 256;
S = 32;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 40;
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function params = SetDefaultParams(...
  params,default_windows,default_trans_types,relativeParams)

fns = params(1:2:end-1);
if ismember('windows',fns) || ...
    (ismember('window_offsets',fns) && ...
     (ismember('window_radii',fns) || ...
      (ismember('min_window_radius',fns) && ...
       ismember('max_window_radius',fns) && ...
       ismember('nwindow_radii',fns)))),
else
  params(end+1:end+2) = {'windows',default_windows};
end
if ~ismember('trans_types',fns),
  params(end+1:end+2) = {'trans_types',default_trans_types};
end
params(end+1:end+2) = {'relativeParams',relativeParams};

author: kristinbranson

repo: JAABA


author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_ft_headmovement

% MEM 3gb
% WALLTIME 00:20:00
% DEPENDENCY ft_headmovement

dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectRest.ds');

cfg = [];
cfg.dataset = dataset;
cfg.method  = 'updatesens';
data1 = ft_headmovement(cfg);

cfg = [];
cfg.dataset = dataset;
cfg.method = 'cluster';
cfg.numclusters = 10;
data2 = cell(1,cfg.numclusters);
[data2{:}] = ft_headmovement(cfg);

cfg          = [];
cfg.dataset  = dataset;
cfg.trl      = [(1:1000:100000)' (1000:1000:100000)'];
cfg.trl(:,3) = 0;
cfg.method = 'pertrial_cluster';
cfg.numclusters = 10;
data3 = cell(1,cfg.numclusters);
[data3{:}] = ft_headmovement(cfg);

cfg         = [];
cfg.dataset = dataset;
cfg.trl      = [(1:1000:100000)' (1000:1000:100000)'];
cfg.trl(:,3) = 0;
cfg.method  = 'avgoverrpt';
data4 = ft_headmovement(cfg);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Saving data into ASPAR type 5 format

report_this_filefun(mfilename('fullpath'));
ZG = ZmapGlobal.Data;
if ~exist('tlen', 'var'); tlen = 30 ; end
str = [];
newmatfile = 't1.sum';
newpath = [hodi '/aspar/'];

do = [ ' cd ' newpath ];
eval(do)

% lets addev the mainshock as the fisrt and largest event...

l = ZG.newt2.Date > ZG.maepi.Date(1) & ZG.newt2.Date < mati + days(tlen);
newt3 =  ZG.newt2(l,1:9);
newt3 = [ZG.maepi(1:1:9) ; newt3 ];

lam = (newt3(:,2)-floor(newt3(:,2)))*100*6/10;
lom = (newt3(:,1)-floor(newt3(:,1)))*100*6/10;

s = [floor(newt3(:,3:5))  newt3(:,8:9) floor(newt3(:,2)) lam  floor(abs((newt3(:,1))))  lom  newt3(:,7) newt3(:,6)];
fid = fopen([newpath newmatfile],'w') ;

fprintf(fid,'%7.3f  %7.3f %7f3\n',[min(newt3(:,6)) tmin1 tlen]);

fprintf(fid,'%5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %7.3f %5.0f %7.3f %7.3f %7.3f\n',s');
fclose(fid);
clear s
return

author: CelsoReyes


text:
function lyrics(varargin)

	query = reshape(transpose([strvcat(varargin{:}),[repmat('+',size(strvcat(varargin{:}),1)-1,1);' ']]),1,[]);

	url = ['''http://search.letssingit.com/cgi-exe/am.cgi?a=search&l=archive&s=',query,''''];

	url(url==' ') = [];

	eval(['web ', url ,' -browser'])

end

author: Sable

repo: mcbench-benchmarks


close all; clear all;
run('../vlfeat-0.9.20/toolbox/vl_setup');

addpath('transforms');
addpath('examples');

REW_worktable;
REW_racetracks;
REW_racetracks;
REW_boardingbridge;
REW_theater;
REW_riverbank;
REW_city;
REW_eyot;
REW_lawn;
REW_parterre;
REW_doorwayreliefs;
REW_building;
REW_campussquare;
REW_gym;
REW_tower;
REW_songyalake_square2;
REW_songyalake_square3;
REW_songyalake_square;
REW_palace;
REW_palace2;
REW_palace3;

APAP_railtracks;
APAP_garden;
APAP_conssite;
APAP_train;

ANAP_building;
ANAP_fense;
ANAP_intersection;
ANAP_roundabout;
ANAP_skyline;

DHW_temple;
DHW_carpark;
DHW_forest;

SPHP_bridge;
SPHP_building;
SPHP_campus;
SPHP_garden;
SPHP_park;
SPHP_street;

NISwGSP_01_PisaCathedralInterior;
NISwGSP_02_Uffizi;
NISwGSP_03_SantaMaria;
NISwGSP_03_SantaMaria_all;
NISwGSP_04_SantaMaria2;
NISwGSP_05_SienaCathedral;
NISwGSP_06_SienaCathedralLibrary;
NISwGSP_07_SienaCathedralLibrary2;
NISwGSP_08_SienaCathedralLibrary3;
NISwGSP_09_SienaCathedralLibrary4;
NISwGSP_10_SienaCathedralInterior;
NISwGSP_11_PiazzaCampo;
NISwGSP_12_PalazzoPubblico;
NISwGSP_13_PalazzoPubblico2;
NISwGSP_cise;
NISwGSP_denny;
NISwGSP_grail;
NISwGSP_Memorial_Hall;
NISwGSP_parrington;

CAVE_atrium;
CAVE_eitan_office;
CAVE_times_square;


text:
function [ h, hmin, hmax ] = p13_stepsize ( option )

%*****************************************************************************80
%
%% P13_STEPSIZE returns step sizes for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =     2.000;
  hmin =  0.001;
  hmax = 10.000;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cvx_clearspath

%CVX_CLEARSPATH   Clears the cvx solver path.
%   CVX_CLEARSPATH removes the internal solver directories to the Matlab
%   path. CVX automatically clears the solver path after completion of a
%   model computation, so calling this function should not be necessary.
%   Nevertheless, we provide it for completeness and debugging.

global cvx___
osolv = cvx___.solvers.active;
if osolv,
    tstr = cvx___.solvers.list(osolv).path;
    if ~isempty( tstr ),
        path(strrep(path,tstr,''));
    end
    cvx___.solvers.active = 0;
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function c = plus(a,b)
% PLUS implements a+b, where either a or b is an adiff object.

switch [class(a),class(b)]

case 'adiffdouble'
   c = adiff(a.x+b, repmat(a.dx,[length(b),1]), a.root);

case 'doubleadiff'
   c = adiff(b.x+a, repmat(b.dx,[length(a),1]), b.root);

case 'adiffadiff'
   checkroot(a,b);
   if size(a.dx,1)~=size(b.dx,1)
      if size(a.dx,1)==1
         a.dx = repmat(a.dx,size(b.dx,1),1);
      elseif size(b.dx,1)==1
         b.dx = repmat(b.dx,size(a.dx,1),1);
      end
   end
   c = adiff( a.x+b.x,a.dx+b.dx, a.root);

otherwise
   error(['Can''t add ',class(a),' and ',class(b)]);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout=yticksi(varargin)
%
% For calling details please see v_yticksi.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_yticksi.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: yticksi.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_yticksi(varargin{:});
else
    v_yticksi(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function Ad = adendotd(dense, d, sparAd, Ablk, blkstart) %#ok
% Ad = Adendotd(dense, d, sparAd, Ablk, blkstart)
%
% ADENDOTD  Computes d[k]'*Aj[k] for Lorentz blocks that are to be factored
%  by dpr1fact.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% 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
%

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function prob_test099 ( )

%*****************************************************************************80
%
%% TEST099 tests LOG_NORMAL_CDF, LOG_NORMAL_CDF_INV, LOG_NORMAL_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST099\n' );
  fprintf ( 1, '  For the Lognormal PDF:\n' );
  fprintf ( 1, '  LOG_NORMAL_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  LOG_NORMAL_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  LOG_NORMAL_PDF evaluates the PDF;\n' );

  a = 10.0;
  b = 2.25;

  check = log_normal_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST099 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = log_normal_sample ( a, b, seed );

    pdf = log_normal_pdf ( x, a, b );

    cdf = log_normal_cdf ( x, a, b );

    x2 = log_normal_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function jed = epoch_to_jed_roman ( )

%*****************************************************************************80
%
%% EPOCH_TO_JED_ROMAN returns the epoch of the Roman calendar as a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%     15 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the epoch.
%
  jed = 1446389.5;

  return
end


text:
function value = torus_square_area_3d ( r1, r2 )

%*****************************************************************************80
%
%% TORUS_SQUARE_AREA_3D returns the area of a square torus in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      R1 - R2 <= SQRT ( X**2 + Y**2 ) <= R1 + R2,
%      -R2 <= Z <= R2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, the two radii that define the torus.
%
%    Output, real TORUS_SQUARE_AREA_3D, the area of the torus.
%
  value = 16.0E+00 * pi * r1 * r2;

  return
end


text:
%% Optimizing the operation of a Hydroelectric Dam
%
% This demo maximizes the revenue generated from a hydroelectric dam using
% Quadratic Programming.  We compute the optimal flow rate to the turbine
% as well as the optimal spill flow rate.
%
% This problem is constrained such that:
%
%  1)  0 < Turbine Flow < 25,000 CFS
%  2)  0 < Spill Flow
%  3)  Max change in Turbine and Spill Flow < 500 CFS
%  4)  Combined Turbine and Spill Flow > 1000 CFS
%  5)  50000 < Reservoir Storage < 100000 Acre-Feet
%  6)  The final reservoir level must equal the initial level (90000 AF)
%
% Two solvers are used to solve this problem: FMINCON and QUADPROG
% FMINCON is used first as it takes less time to set-up the problem for
% ths solver.  Supplying the Gradient and Hessian to the FMINCON algorithm
% significantly speeds up the solution of the problem.
% QUADPROG ionterior-point-convex is used third to provide a faster
% algorithm for solving this type of problem.  This solver is also better
% suited for handling a larger version of the same problem.
%
% Copyright (c) 2012, MathWorks, Inc.

%% 1 - Load data
% Two variables are loaded from data set:
%   inFlow - flow rate into the reservoir measured hourly (CFS)
%   price - price of electricity measured hourly ($/MWh)
load('FlowAndPriceData.mat');

%% 2 - Define Starting Point and Constants
% We wish to find the vector x of optimal decision variables.
% For this problem:
% x = [turbineFlow(1) ... turbineFlow(N) spillFlow(1) ... spillFlow(N)]

% Our initial starting point is turbineFlow = inFlow and spillFlow = 0
x0 = [inFlow; zeros(size(inFlow))];

stor0 = 90000; % initial vol. of water stored in the reservoir (Acre-Feet)

k1 = 0.00003; % K-factor coefficient
k2 = 9; % K-factor offset

MW2kW = 1000; % MW to kW
C2A = 1.98347/24; % Convert from CFS to AF/HR

C = [C2A k1 k2 MW2kW]; % Vector of constants

N = length(inFlow);

%% 3 - Define Linear Inequality Constraints and Bounds
DefineConstraints;

%% 4 - Objective function
% Use symbolic math for creating the objective function and then convert
% into a MATLAB function that can be evaluated.
%% 4.1 - Derive Objective Function using Symbolic Math
% First define the decision variables as X, which is a vector consisting
% of [TurbineFlow, SpillFlow] and all the parameters as symbolic objects
X = sym('x',[2*N,1]);   % turbine flow and spill flow
F = sym('F',[N,1]);     % flow into the reservoir
P = sym('P',[N,1]);     % price
s0 = sym('s0');         % initial storage
c = sym({'c1','c2','c3','c4'}.','r');   % constants

%%
% Total out flow equation
TotFlow = X(1:N)+X(N+1:end);

%%
% Storage Equations
S = cell(N,1);
S{1} = s0 + c(1)*(F(1)-TotFlow(1));

for ii = 2:N
    S{ii} = S{ii-1} + c(1)*(F(ii)-TotFlow(ii));
end

%%
% K factor equation
k = c(2)*([s0; S(1:end-1)]+ S)/2+c(3);

% MWh equation
MWh = k.*X(1:N)/c(4);

%%
% Revenue equation
R = P.*MWh;

%%
% Total Revenue (Minus sign changes from maximization to minimization)
totR = -sum(R);

%% 4.2 - Create Function File that can be Evaluated in MATLAB
% Create a function that when evaluated returns the value at a point

% Perform substitutions so that X is the only remaining variable
totR = subs(totR,[c;s0;P;F],[C';stor0;price;inFlow]);

% Generate a MATLAB file from the equation
matlabFunction(totR,'vars',{X},'file','objFcn');

%% 5 - Optimize with FMINCON
options = optimset('MaxFunEvals',Inf,'MaxIter',5000,...
    'Algorithm','interior-point','Display','iter');

tic
[x1, fval1] = fmincon(@objFcn,x0,A,b,Aeq,beq,LB,UB,[],options);
toc

% Visualize Results
plotResults(price, x1, N);

%% 6 - Optimize with FMINCON with supplied Gradient and Hessian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Let's see if we can speed up the optimization by supplying the gradient
% and Hessian.  Passing these to the algorithm allows it to make more
% accurate decisions because it no longer has to approximate the gradient
% and hessian (approximation also takes longer as it requires many
% objective function evaluations).

%% 6.1 - Gradient
% Use the GRADIENT function from Symbolic Math to calculate the gradient
gradObj = gradient(totR,X);

% Create a function that when evaluated returns the gradient at a point
matlabFunction(gradObj,'vars',{X},'file','grad');

% The optimization routine expects two return arguments for a user-supplied
% gradient.  The first argument is the value of the objective function and
% the second argument is the gradient.  The DEAL function is used to return
% multiple output arguments through a function handle.
ValAndGrad = @(x) deal(objFcn(x),grad(x));

%% 6.2 - Hessian
% Use the HESSIAN function from Symbolic Math to calculate the Hessian
hessObj = hessian(totR,X);

% Create a function that when evaluated returns the Hessian at a point
matlabFunction(hessObj,'vars',{X},'file','Hess');

% The Hessian function has two input arguments (the current point and the
% Lagrange Multipliers) and returns the value of the Hessian at the current
% point.
hfcn = @(x,lambda) Hess(x);

%% 6.3 - Optimize
options = optimset('Disp','iter','Algorithm','interior-point',...
                   'MaxFunEvals',500,'MaxIter',5000,...
                   'GradObj','on','Hessian','user-supplied',...
                   'HessFcn',hfcn);

tic
[x2,fval2] = fmincon(ValAndGrad,x0,A,b,Aeq,beq,LB,UB,[],options);
toc

% Visualize Results
plotResults(price, x2, N);

%% 7 - Optimize using the QUADPROG Solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QUADPROG minimizes functions of the type (1/2) * x' * H * x  +  f' * x
% First perform some substitutions into the revenue equation
Rsub = subs(R,[c;s0;P;F],[C';stor0;price;inFlow]);

% Total revenue (minus sign to flip maximization to minimization)
Rtot = -sum(Rsub);

%% 7.2 - Calculate and evaluate the linear component
fsym = gradient(Rtot,X);
f = double(subs(fsym,X,zeros(size(X))));

%% 7.1 - Calculate and evaluate the quadratic matrix
H = 0.5.*double(hessian(Rtot,X));

%% 7.3 - Optimize
qpoptions = optimset('Algorithm','interior-point-convex','Disp','iter');

tic
[x3,fval3] = quadprog(H,f,A,b,Aeq,beq,LB,UB,[],qpoptions);
toc

% Visualize Results
plotResults(price, x3, N);


text:
function angdataw=angtimewarp(evLatency, newLatency, angdata)
% ANGTIMEWARP - Given two event marker vectors, computes a
%                 warping of the input angular time series so that its
%                 evlatencies match newlatencies. Values of the warped
%                 timeserie that falls between two frames in the original
%                 timeserie will be linearly interpolated under the
%                 assumption that phase change is minimal between two
%                 successive time points.
% Usage:
%   >> warpAngs = angtimewarp(evlatency, newlatency, angData)
%
% Necessary inputs:
%   evlatency  - [vector] time markers on the original time-series, in
%                frames. Markers must be ordered by increasing
%                latency. If you want to warp the entire time series,
%                make sure frame 1 and the last frame are in the vector.
%   newlatency - [vector] desired time marker latencies. The original
%                time series will be warped so that its time markers (see
%                evlatency) match the ones in newlatency. newlatency
%                frames must be sorted by ascending latencies in frames.
%                Both vectors have to be the same length.
%   angData    - [vector] original angular time series (in radians).
%                Angles should be between -pi and pi.
%
% Optional outputs:
%   warpAngs   - [vector] warped angular time-course, with values between
%                -pi and pi
%
% Example:
%   >> angs = 2*pi*rand(1,10)-pi;
%   >> warpangs = angtimewarp([1 5 10], [1 6 10], angs)
%
% Authors: Jean Hausser, SCCN/INC/UCSD, 2006
%
% See also: TIMEWARP, PHASECOHER, ERPIMAGE, NEWTIMEF

  if min(sort(evLatency) == evLatency) == 0
    error('evlatency should be sorted');
    return;
  end
  if min(sort(newLatency) == newLatency) == 0
    error('newlatency should be sorted');
    return;
  end
  if length(evLatency) ~= length(newLatency)
    error('evlatency and newlatency must have the same length.');
    return;
  end
  if length(evLatency) < 2 || length(newLatency) < 2
    error(['There should be at least two events in evlatency and ' ...
          'newlatency, that is "begin" and "end"' ]);
    return;
  end
  if evLatency(1) ~= 1
    disp(['Assuming old and new time series beginnings are ' ...
          'synchronized']);
    disp(['Make sure you defined an end event for both old and new time ' ...
          'series !']);
    evLatency(end+1)=1;
    newLatency(end+1)=1;
    evLatency = sort(evLatency);
    newLatency = sort(newLatency);
  end

  t = 1:max(evLatency);

  for k=1:length(evLatency)-1
    for i=evLatency(k):evLatency(k+1)-1
      tp(i) = (t(i)-evLatency(k)) * ...
              (newLatency(k+1) - newLatency(k))/...
              (evLatency(k+1) - evLatency(k)) + ...
              newLatency(k);
    end
  end

  %Check what's going on at tp(max(newLatency)), should equal t(max(evLatency))
  tp(max(evLatency)) = max(newLatency);
  ts = tp-min(newLatency)+1;

  angdataw = zeros(1, max(newLatency)-min(newLatency)+1);

  k = 0;
  for i=1:length(angdataw)
    while i > ts(k+1)
      k = k+1;
    end

    if k == 0

      angdataw(1) = angdata(1);
    else

      %Linear interp
      angdataw(i) = angdata(k)*(1 - (i-ts(k))/(ts(k+1)-ts(k))) + ...
                    angdata(k+1)*(1 - (ts(k+1)-i)/(ts(k+1)-ts(k)));

%       %Correction because angles have a ring structure
%       theta1 = [angdata(k) angdata(k+1) angdataw(i)];
%       theta2 = theta1 - min(angdata(k), angdata(k+1));
%       theta2max = max(theta2(1), theta2(2));
%       if ~ ( (theta2max <= pi & theta2(3) <= theta2max) | ...
%              (theta2max >= pi & theta2(3) >= theta2max) | ...
%              theta2(3) == theta2(1) | theta2(3) == theta2(2) )
%         angdataw(i) = angdataw(i) + pi;
%       end
%       if angdataw(i) > pi %Make sure we're still on [-pi, pi]
%         angdataw(i) = angdataw(i) - 2*pi;
%       end
    end
  end
  angdataw = wrap2pi(angdataw);

  function a = wrap2pi(a, a_center )
% function a = wrap(a,a_center)
%
% Wraps angles to a range of 2*pi.
% Inverse of Matlab's "unwrap", and better than wrapToPi ( which has
% redundant [-pi,pi])
% Optional input "a_center" defines the center angle.  Default is 0, giving
% angles from (-pi,pi], chosen to match angle(complex(-1,0)).  Maximum
% possible value is pi.

% T.Hilmer, UH
% 2010.10.18 version 2
%   removed code from version 1. Have not bug-checked second input
%   "a_center"

if nargin < 2, a_center = 0; end

% new way
a = mod(a,2*pi); % [0 2pi)

% shift
j = a > pi - a_center;
a(j) = a(j) - 2*pi;
j = a < a_center - pi;
a(j) = a(j) + 2*pi;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function X = stack2vecs(S)
%STACK2VECS Convert stack to observation-by-features matrix.
% Usage:
%   X = stack2vecs(S)
%
% Args:
%   S: 4-d array of size [h, w, c, n]
%
% Returns:
%   X: 2-d array of size [n, (h * w * c)]
%
% See also: vecs2stack

if ndims(S) == 3; S = permute(S,[1 2 4 3]); end

X = reshape(S,[],size(S,4))';

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function [isSameModel,differences] = compareModels(modelA,modelB,printLevel)
%compares modelA with modelB, looking for differences between the
%structures
%
% INPUT
% modelA:       structure
% modelB:       structure
% printLevel:
%
% OUTPUT
% isSameModel:   true if identical models, false otherwise
% differences:   structure listing differences between models
%               *.reason: gives a text stack of why the difference occurred
%                         as well as a field
%               *.where: contains the indices and subfields of the structure
%                        where the comparison failed.
%
%
% Note: This function depends on structeq.m and celleq.m

if~exist('printLevel','var')
    printLevel=1;
end

differences=[];
result = 0;
i=0;
while result==0
    [result, why] = structeq(modelA,modelB);
    if result==0
        i=i+1;
        if printLevel>0
            fprintf('%s\n',why.Reason)
            fprintf('%s\n',why.Where)
        end
        differences(i).reason = why.Reason;
        differences(i).where = why.Where;
        %fieldName = strrep(why.Where,'(1)','');
        fieldName = strtok(why.Where,'(1)');
        eval(['modelB' fieldName ' = modelA' fieldName ';']);
    end
end
isSameModel = i==0;

end

author: opencobra

repo: cobratoolbox


text:
classdef IN1KMOP6 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'transformer';
            config.args.objs = 'err&params&flops';
            config.args.normalized_objectives = false;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [W, info] = MTFLCd_ADMM_WSolver...
    (Xdiag, yvect, Thvect, Zvect, rho, lambda1, lambda2, opts)
%
% Multi-Task Feature Learning with Calibration - ADMM
% Subproblem: W
% diagnoalized version.
%
% Objective
%   min_W {   rho/2 sum_i^m ||theta_i/rho - z_i + y+i - X_iwi||
%           + lambda2/2 ||W||_F^2 + lambda1 ||W||_{1,2}
%          }
%
% INPUT
%  X - cell array of {n_i by d matrices} by m
%  y - cell array of {n_i by 1 vectors}  by m
%  Th - cell array of {n_i by 1 vectors} by m
%  Z  - cell array of {n_i by 1 vectors} by m
%  rho - parameter of the augmented Lagrange
%  lambda1 - regularization parameter of the l2,1 norm penalty
%  lambda2 - regularization parameter of the Fro norm penalty
%  opts - MANDATORY optimization options.
%
% OUTPUT
%  W - task weight d by m.
%
% Author: Jiayu

%% Initialization

W0 = opts.init; % the outer loop must pass in the last solution

funcVal = zeros(opts.maxIter,1);

TrZy = Thvect/rho - Zvect + yvect;

%% Computation

bFlag = 0;

W     = W0;
W_old = W0;

gamma = 1; gamma_inc = 2;

t =1; t_old = 1;
for iter = 1: opts.maxIter
    alpha = (t_old  -1 )/t;
    V = W + alpha * (W - W_old);

    [fV, gV] = smoothObj(V);

    for lsIter = 1:100
        W = proj(V - gV/gamma, lambda1/gamma);
        f = smoothObj(W);

        delta_W = W - V;
        r_sum = sum(sum(delta_W.^2));

        if(r_sum <= 1e-20), bFlag = 1; break; end

        if f<= fV + sum(sum(delta_W .* gV)) + gamma/2 * r_sum
            break;
        end

        gamma = gamma * gamma_inc;
    end

    W_old = W;

    funcVal(iter) = f + lambda1 * sum(sqrt(sum(W.^2, 2)));

    if(bFlag), break; end

    % convergence
    if(iter>1)
        if (abs( funcVal(iter) - funcVal(iter-1) ) ...
                <= opts.tol* abs(funcVal(iter-1)))
            break;
        end
    end

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

info.funcVal = funcVal(1:iter);

%% Nested functions
    function [X] = proj(D, t) % l2.1 norm projection.
        X = repmat(max(0, 1 - t./sqrt(sum(D.^2,2))),1,size(D,2)).*D;
    end

    function [f, g] = smoothObj(W)
        f = 0.5 * lambda2 *sum(sum(W.^2));
        tpVect = TrZy - Xdiag * W(:);
        f = f +  0.5 * rho * sum(tpVect.^2);
        if nargout >= 2
            g = reshape(lambda2 * W(:) - rho * (Xdiag' * tpVect), size(W));
        end
    end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function Next = SurrogateAssistedSelection(Problem,net,p0,p1,Ref,Input,wmax,tr)
% Surrogate-assisted selection for selecting promising solutions

%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    Next  = OperatorGA(Problem,[Input;Ref.decs],{1,15,1,5});
    Label = net.predict(Next);
    a     = tr;
    b     = 1 - tr;
    i     = 0;
    if p0<0.4 || (p1<a&&p0<b)
        while i < wmax
            [~,index] = sort(Label,'descend');
            Input     = Next(index(1:length(Ref)),:);
            Next      = OperatorGA(Problem,[Input;Ref.decs],{1,15,1,5});
            Label     = net.predict(Next);
            i = i+size(Next,1);
        end
        Next = Next(Label>0.9,:);
    elseif p0>b && p1<a
        Next = [];
    elseif p1 > b
        while i<wmax
            [~,index] = sort(Label);
            Input     = Next(index(1:length(Ref)),:);
            Next      = OperatorGA(Problem,[Input;Ref.decs],{1,15,1,5});
            Label     = net.predict(Next);
            i = i+size(Next,1);
        end
        Next = Next(Label<0.1,:);
    else
        Next = Next(randi(end),:);
    end
end


text:
function [dims, isndarray, maxlevel, count] = nestbracket2dim(str,brackets,testndarray)
%
% [dims, isndarray, maxlevel, count] = nestbracket2dim(str,brackets)
%
% Extracting the dimension vector of a JSON string formatted array
% by analyzing the pairs of opening/closing bracket tokenss; this function
% only returns valid dimension information when the array is an N-D array
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      str: a string-formatted JSON array using square-brackets for enclosing
%           elements and comma as separators between elements
%      brackets: (optional), a string of length 2, with the first character
%               being the opening token and the 2nd being the closing token.
%               if not given, brackets is set to '[]' to find matching square-brackets;
%               for example, '{}' looks for a matching closing curly-bracket in
%               the string key(pos(startpos,:end))
%      testndarray: (optional), 1 to test if the input string contains an
%               ND array, i.e. with uniform element lengths (recursively)
%
% output:
%      dims: the speculated dimension vector with the length matching the maximum
%            depth of the embedded bracket pairs. When the input string encodes an
%            N-D array, the dims vector contains all integers; however, returning
%            an all-integer dims vector does not mean the array is
%            rectangular. if testndarray is set to 1, dims returns isndarray
%      isndarray: 1 to indicate the input string contains an ND array,
%            otherwise, 0
%      maxlevel: return the depth of the enclosed brackets in the string, i.e. the
%            length of the dims vector.
%      count: the relative depth from the level 0 - scanning from the left
%            to right of the string, an opening token increases the level by 1
%            and a closing token decreases the level by 1; a zero indicates
%            the positions of a matching bracket of the same level.
%
% example:
%      str='[[ [1,2,3], [4,2,1]], [ [10,1,0], [2,5,10]] ]'; % an N-D array
%      [dim,dep]=nestbracket2dim(str)
%      str='[[ [1,2,3], [4,2,1]], [ [10,1,0], [2,5]] ]'; % an invalid N-D array
%      [dim,dep]=nestbracket2dim(str)

% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin<2)
    brackets='[]';
end
str=str(str==brackets(1) | str==brackets(2) | str==',');
count=cumsum(str==brackets(1)) - cumsum(str==brackets(2));
isndarray=testuniform(count, max(count));
if(nargin>2 && testndarray)
    dims=isndarray;
    return;
end
maxlevel=max(count);
dims=histc(count,1:maxlevel);
dims(1:end-1)=dims(1:end-1)*0.5;
dims(2:end)=dims(2:end)./dims(1:end-1);
dims=fliplr(dims);

function isndarray=testuniform(count, maxval)
isndarray=false;
if(length(count)>2)
    idx=find(count(2:end)==count(1),1);
    if(idx==length(count)-2 && count(1)<maxval)
        isndarray=testuniform(count(2:end-1), maxval);
        return;
    end
    if(~isempty(idx) && mod(length(count)-1,idx+1)==0)
        count2d=reshape(count(1:end-1),idx+1,[]);
        if(all(diff(count2d')==0))
            isndarray=true;
        end
    end
end


text:
function value = i4_min ( a, b )

%*****************************************************************************80
%
%% I4_MIN returns the minimum of two I4's.
%
%  Discussion:
%
%    An I4 is an integer value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, B, values to compare.
%
%    Output, integer VALUE, the minimum of A and B.
%
  a = floor ( a );
  b = floor ( b );

  if ( a < b )
    value = a;
  else
    value = b;
  end

  return
end


text:
function [Isub,domSub,femSub] = femSubdivide(dom,fem,Nsub,fig)
%+========================================================================+
%|                                                                        |
%|              OPENFEM - LIBRARY FOR FINITE ELEMENT METHOD               |
%|           openFem is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018.          |
%| 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                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : femSubdivide.m                                |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & François Alouges            |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Subdivide finite element space for domain     |
%|  `---'  |                decomposition method                          |
%|         |                ! ONLY WORKS FOR NON DIRICHLET FEM !          |
%+========================================================================+

% Number of subdomains is power of 2
if (Nsub ~= 2^floor(log2(Nsub)))
    error('femSubdivide.m : number of subdomain is not a power of 2.')
end

% Mesh and dof associated to finite element space
mesh          = fem.msh;
[dof,elt2dof] = fem.dof;                                   % to be improved

% Dof subdivision with binary tree
Nleaf  = ceil(length(fem)/Nsub);
bitree = tree(msh(dof),'binary',Nleaf);
Isub   = bitree{end}.ind;

% Security
if (length(Isub) ~= Nsub)
    error('femSubdivide.m : unavailable case')
end
if (norm(sort(cell2mat(Isub)) - (1:length(fem))','inf') > 1e-12)
    error('femSubdivide.m : unavailable case')
end

% Output initialization
domSub = cell(Nsub,1);
femSub = cell(Nsub,1);

% For each subindices
for i = 1:Nsub
    % Indices of valid dof and extended mesh
    I       = ismember(elt2dof,Isub{i});
    meshSub = mesh.sub(sum(I,2)>0);

    % Quadrature
    domSub{i}     = dom;
    domSub{i}.msh = meshSub;

    % Finite element
    femSub{i}     = fem;
    femSub{i}.msh = meshSub;

    % Dirichlet condition for boundary dof
    dir = setdiff( msh(femSub{i}.dof) , msh(dof(Isub{i},:)) );
    if (size(dir,1)>0)
        femSub{i} = dirichlet(femSub{i},dir);
    end

    % Security
    if (length(femSub{i}) ~= length(Isub{i}))
        error('femSubdivide.m : unavailable case');
    end

    % Graphical representation
    if fig
        tmp = mesh.sub(sum(I,2)==size(elt2dof,2));
        figure(fig)
        hold on
        plot(tmp,'b')
        plot(tmp.bnd,'r')
        plot(setdiff(meshSub,tmp),'w')
        axis equal
    end
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function gj_2d_ani(before,after,x,xy_ans,varargin)

params = inputParser;
params.CaseSensitive = false;
params.addParameter('record', false);
params.addParameter('record_filename','record');
params.parse(varargin{:});

%Extract values from the inputParser
h_record = params.Results.record;
record_filename = params.Results.record_filename;

if h_record
    v = VideoWriter([record_filename,'.mp4'],'MPEG-4');
    v.FrameRate = 60;
    v.Quality = 100;
    open(v);
end

n_steps = 200;

diff = after-before;
% figure;
for i_step = 1:n_steps
    tempA = before+diff/n_steps*i_step;

    y1 = (-tempA(1,1)*x+tempA(1,3))/(tempA(1,2)+eps);
    y2 = (-tempA(2,1)*x+tempA(2,3))/(tempA(2,2)+eps);

    plot(x,y1);
    hold on;
    plot(x,y2);
    plot(xy_ans(1,1), xy_ans(2,1),'r.','markersize',20)

    mArrow2(xy_ans(1,1), xy_ans(2,1), xy_ans(1,1)+tempA(1,1),xy_ans(2,1)+tempA(1,2),{'color',[0 0.4470 0.7410]});
    mArrow2(xy_ans(1,1), xy_ans(2,1), xy_ans(1,1)+tempA(2,1),xy_ans(2,1)+tempA(2,2),{'color',[0.85 0.325 0.098]});
    grid on;
    hold off;
    xlim([-5 5])
    ylim([-5 5])

    set(gcf,'color','w')

    if h_record
        F=getframe(gcf);
        writeVideo(v,F);
    end

    pause(0.01);

end

end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function [in3] = l2in3(l)
% Convert volume from liters to cubic inches.
% Chad Greene 2012
in3 = l*61.023744095;


text:
function [weight, exp_w] = anticor_kernel(data, W, exp_ret, exp_w)
% This program output the final portfolio the BAH(Anticor) strategy
% BAH(Anticor) has one folds of experts
%
% function [weight, exp_w] = anticor_kernel(data, W, exp_ret, exp_w)
%
% weight: final portfolio, used for next rebalance
% exp_w: experts weights in the first fold
%
% data: market sequence vectors
% W: maximum window size, the number of experts (W-1)
% exp_ret: experts' return in the first fold
% exp_w: experts' weights in the first fold
%
% Example: [weight, exp_w] = anticor_kernel(data, W, exp_ret, exp_w)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=2:W,
    exp_w(k-1, :) = anticor_expert(data, exp_w(k-1, :)', k);
end

% Combine portfolios according to q(k, l) and previous expert return
numerator = 0;
denominator = 0;
for k=2:W,
    numerator = numerator + exp_ret(k-1, 1)*exp_w(k-1, :);
    denominator = denominator + exp_ret(k-1, 1);
end

weight = numerator'/denominator;

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [ wgft_coefficients ] = WGFT_wgft(g,V,signal,varargin)

% WGFT coefficients are output in an NxN matrix. The ith column corresponds
% to translation to vertex i. The kth row corresponds to frequency
% \lambda_{k-1}

N=size(V,1);
wgft_coefficients=zeros(N,N);

if ~isempty(varargin)
    param=varargin{1};
    for i=1:N
        for k=1:N
            g_ik=WGFT_atom(g,V,i,k,param);
            wgft_coefficients(k,i)=g_ik'*signal;
        end
    end
else
    for i=1:N
        for k=1:N
            g_ik=WGFT_atom(g,V,i,k);
            wgft_coefficients(k,i)=g_ik'*signal;
        end
    end
end

end

repo: gspbox


text:
% The box [x1 y1 x2 y2] here is different for box in matlab [x y w h]
function o = box_overlap(a, b)

% Compute the symmetric intersection over union overlap between a set of
% bounding boxes in a and a single bounding box in b.
%
% a  a matrix where each row specifies a bounding box
% b  a single bounding box

x1 = max(a(:,1), b(1));
y1 = max(a(:,2), b(2));
x2 = min(a(:,3), b(3));
y2 = min(a(:,4), b(4));

w = x2-x1+1;
h = y2-y1+1;
inter = w.*h;
aarea = (a(:,3)-a(:,1)+1) .* (a(:,4)-a(:,2)+1);
barea = (b(3)-b(1)+1) * (b(4)-b(2)+1);
% intersection over union overlap
o = inter ./ (aarea+barea-inter);
% set invalid entries to 0 overlap
o(w <= 0) = 0;
o(h <= 0) = 0;

author: ShapeNet

repo: RenderForCNN


text:
% mean ROI tcs across all sessions
%
% kgs 1/2007
clear all
% for PC:
expDir = ('y:\RAID\projects\Kids\fmri\adolescents');
% for linux:
%expDir='/biac1/kgs/projects/Kids/fmri/adolescents';

sessions = { 'ar_12yo_011908' 'is_16yo_121607' 'kwl_14yo_011208' 'dw_14yo_102707'};
roi=      {'lPPA_IOvsOJ_p3' 'lPPA_IOvsOJ_p3' 'lPPA_IOvsOJ_p3' 'lPPA_IOvsOJ_p3'};

%% data type:
dt =    [4 3 4 3]; % motion corrected both within and between
%  scan number
scan = [7 8 8 8]; %
nsbj=length(sessions);

% get some basic params from first session e.g. eventParams
cd( fullfile(expDir, sessions{1}) );
fprintf('*** Session 1 *** \n', sessions{1});
% initialize hidden inplane
hI = initHiddenInplane(dt(1), scan(1));

% params is a struct in which i set the following fields the following fields:
params.dataType=dt;
params.scan=scan;
params.studyDir=expDir;
params.viewType='INPLANE';

% enforce all analysis params to be the same
params.eventParams  = er_getParams(hI);

params.eventParams.ampType='difference'; % params.eventParams.ampType can be one of the following:
                                        % 'difference' or 'betas' or 'relamps' or 'deconvolved';
                                        % this field specifies which amptype will be plotted in the meanamps plot
                                        % and which amps will be calculated in the tc.allMeanAmps and tc.meanAmps
params.eventParams.peakPeriod=[4:12];

params.methodFlag=4; % method to combine Tcs across subjects; Across-sessions analysis ("random effects")
                     % Note that the function tc_acrossSessions
                     % first generates the subjectTcs array which is a
                     % structure array of individual subjects tc
                     % and then calls tc_combineTcs to calculate the across
                     % subject average; Check tc_combineTcs to learn about
                     % methods to combine data across subjects

params.openUI=1; % 1 to open TC UI with average time course across subjects
                 % will open a left panel with mean amps according to params.eventParams.ampType field
                 % and a right panel with the event trigged average time
                 % course (except if
                 % params.eventParams.ampType='deconvolved' in that case
                 % will show the mean deconvolved time course

% set glm and analysis params to be the same across sessions
params.eventParams.detrend = 1; % high-pass filter (low-pass is blurring, below)
params.eventParams.inhomoCorrection = 1; % divide by mean
params.eventParams.temporalNormalization = 0; % no matching 1st temporal frame

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set HRF function to be the GLM SPM HRF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
params.eventParams.glmHRF = 3; %
% number of events per block
params.eventParams.eventsPerBlock = 6; % see above

params.sessionPlots={'betas' 'meanAmps' 'meantcs'}; % single subject plots - see all options in tc_sessionPlots
params.savePlotsPath='Images';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop across sessions
% returns mean tc across subject in the tc variable
% and structure array containing single subject tcs in subjectTcs
% calls tc_combineTcs and tc_sessionPlots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[tc subjectTcs] = tc_acrossSessions(sessions, roi, params);

author: vistalab

repo: vistasoft


text:
function mrUtilInvertAndSaveFig(fn,trans)
%
% mrUtilInvertAndSaveFig(fileName,[trans=true])
%
%
%

if(~exist('trans','var')|isempty(trans))
    trans = true;
end

tcol = [1 2 3];
fc = [1 1 1];
if(trans)
    set(gcf,'InvertHardCopy','off','color',tcol/255);
else
    set(gcf,'InvertHardCopy','off','color',[0 0 0]);
end
set(gca,'ycolor',fc,'xcolor',fc,'color', fc);
set(get(gca,'XLabel'),'color',fc);
mrUtilPrintFigure([fn '.eps']);
%print2im([fn '.png'], gcf, true, 120);
unix(['pstoimg -antialias -aaliastext -density 300 -type png -crop a -out ' fn '.png ' fn '.eps']);
if(trans)
    [im,cm] = imread([fn '.png']);
    transInd = cm(:,1)==tcol(1)/255 & cm(:,2)==tcol(2)/255 & cm(:,3)==tcol(3)/255;
    imwrite(im,cm,[fn '.png'],'Transparency',double(~transInd));
end

return


text:
function C=cofactorMatrix(A)
%%COFACTORMATRIX Compute the cofactor matrix of a square matrix A. The cofactor
%     matrix of a matrix A is one such that C(i,j)=(-1)^(i+j)det(A_{i,j}), where
%     A_{i,j} is the matrix A with row i and column j removed.
%
%INPUTS: A An nXn square matrix. The elements can be real or complex. The
%          matrix can be singular.
%
%OUTPUTS: C The nXn cofactor matrix of the matrix A.
%
%The adjugate of a matrix is defined in Ch. 4 of [1] as the transpose of the
%cofactor matrix C. This function just calls adj(A) to get the adjugate and then
%takes the transpose of the result.
%
%REFERENCES:
%[1] K. M. Abadir and J. R. Magnus, "Matrix algebra," New York, 2005.
%
%August 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

C=adj(A)';

end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [header_size, version, f_height, f_width, bytes_per_chunk, ...
          max_n_frames, data_format] = fmf_read_header( filename )
% [header_size,version,nr,nc,bytes_per_chunk,nframes,data_format] = fmf_read_header(filename)
%
% reads FlyMovieFormat header data from FILENAME
%
% HEADER_SIZE is in bytes for the appropriate file format (fmf VERSION)
% F_HEIGHT and F_WIDTH are the number of pixels in a frame
% BYTES_PER_CHUNK is the number of bytes per frame in the file
% MAX_N_FRAMES is the number of frames in the file
%
% JAB 7/1/04
%
% Fixed bug in reading number of frames.
%
% KMB 05/22/07

header_size = 28;

fp = fopen( filename, 'r' );

% read header
version = double( fread( fp, 1, 'uint32' ) );
if (version ~= 1) && (version ~= 3),
  error( 'version not supported -- FMF versions 1 and 3 only' );
end

if version == 1,
  data_format = 'MONO8';
  bits_per_pixel = 8;
end

if version == 3,
  format_len = double( fread( fp, 1, 'uint32' ) );
  data_format = strcat( char( fread( fp, format_len))');
  bits_per_pixel = double( fread( fp, 1, 'uint32' ) );
  header_size = header_size + format_len + 8;
end

f_height = double( fread( fp, 1, 'uint32' ) );
f_width = double( fread( fp, 1, 'uint32' ) );
bytes_per_chunk = double( fread( fp, 1, 'uint64' ) );
frame_count_location = ftell(fp);
max_n_frames = double( fread( fp, 1, 'uint64' ) );

if max_n_frames <= 0,
  fprintf('Resorting to computing number of frames from size of file.\n');
  pcurr = ftell(fp);
  fseek(fp,0,'eof');
  nbytestot = ftell(fp);
  fseek(fp,pcurr,'bof');
  max_n_frames = floor((nbytestot - header_size) / bytes_per_chunk);
end

fclose( fp );

author: kristinbranson

repo: JAABA


text:
clear all; close all; clc;
se1=strel('square', 3)
se2=strel('line', 10, 45)

author: HuangCongQing


text:
function crackFVMP2
%% CRACK Problem
%
% crack solves the Poisson equation $-\Delta u =f$ in $\Omega$ and $u =
% g_D$ on $\partial \Omega$ in a crack domain
% $\Omega=\{|x|+|y|<1\}\backslash \{0<= x <=1, y=0\}$
%  using adaptive finite element method (AFEM). We choose f=1 and g_D such
%  that the exact solution is $u = r^{\beta}\sin(\beta\theta)-0.25r^2,
%  \beta = 1/2$ in the polar coordinate.
%
% EXAMPLE
%
%    crack
%
% See also  crack_performance, Lshape
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all;clc;clear;
%% Parameters
maxN = 2e3;     theta = 0.5;    maxIt = 50;
N = zeros(maxIt,1); err = zeros(maxIt,3);
%%  Generate an initial mesh
node = [1,0; 0,1; -1,0; 0,-1; 0,0; 1,0];        % nodes
elem = [5,1,2; 5,2,3; 5,3,4; 5,4,6];            % elements
elem = label(node,elem);                        % label the mesh
bdEdge = setboundary(node,elem,'Dirichlet');    % Dirichlet boundary condition
showmesh(node,elem);                            % plot mesh
findelem(node,elem);                            % plot element indices
findnode(node,1:5);                             % plot node indices
findedge(node,[1 5; 5 6],1);                    % plot the crack edge
text(node(6,1)+0.05,node(6,2)+0.075,int2str(6), ...
     'FontSize',14,'FontWeight','bold');
%%
% node 1 and node 6 are the same point (1,0)

%%  Get a fine mesh by uniform bisection
for k = 1:2
    [node,elem,bdEdge] = uniformbisect(node,elem,bdEdge);
end
clf; showmesh(node,elem);
%% Set up PDE data
pde.f = @f;
pde.exactu = @exactu;
pde.g_D = @exactu;
pde.Du=@Du;
pde.d=[];
%%  Adaptive Finite Element Method
% *SOLVE* -> *ESTIMATE* -> *MARK* -> *REFINE*
for k = 1:maxIt
    % Step 1: SOLVE
   [u,A,~,~,edge] = PoissonP2FVM(node,elem,[],@pde.f,@pde.g_D);
    % Plot mesh and solution
    figure(1);  showresult(node,elem,u(1:size(node,1)),[-7,12]);
    % Step 2: ESTIMATE
   % eta = estimaterecovery(node,elem,u); % recovery type % need eta^2 for mark
    eta = recoverFlux2(node,elem,pde,bdEdge,u);
    % Record error and number of vertices
   uI = pde.exactu([node; (node(edge(:,1),:)+node(edge(:,2),:))/2]);
   err(k,1) = sqrt((u-uI)'*A*(u-uI));
   err(k,2) = getL2error(node,elem,@exactu,u);
   err(k,3) = getH1error(node,elem,@pde.Du,u);
    N(k) = size(node,1);
    if (N(k)>maxN), break; end
    % Step 3: MARK
    markedElem = mark(elem,eta,theta);
    % Step 4: REFINE
    [node,elem,bdEdge] = bisect(node,elem,markedElem,bdEdge);
end

%% Plot convergence rates
N= N(1:k);
erruIuh= err(1:k,1); errL2 = err(1:k,2); errH1 = err(1:k,3);
figure(2);
r1 = showrate(N,erruIuh,5,'-*');
hold on;
r2 = showrate(N,errL2,5,'k-+');
r3 = showrate(N,errH1,10,'r-+');
legend('||Du_I-Du_h||',['cN^{' num2str(r1) '}'], ...
       '||u-u_h||',['cN^{' num2str(r2) '}'], ...
       '||Du-Du_h||',['cN^{' num2str(r3) '}'], ...
       'LOCATION','Best');
%%
% Using AFEM, we obtain optimal rate of convergence the error in the energy
% norm ($N^{-0.5}$) and in the $L^2$ norm ($N^{-1}$).
end % End of function CRACK

%%  Data of CRACK
function z = f(p)   % load data (right hand side function)
z = ones(size(p,1),1);
end

function z = exactu(p) % exact solution
r = sqrt(sum(p.^2,2));
z = sqrt(0.5*(r-p(:,1)))-0.25*r.^2;
end

function z = Du(p) % derivative of the exact solution
r = sqrt(sum(p.^2,2));
z(:,1) = (p(:,1)./r-1)./sqrt(8*(r-p(:,1)))-0.5*p(:,1); % u_x
z(:,2) = p(:,2)./r./sqrt(8*(r-p(:,1)))-0.5*p(:,2);     % u_y
end

repo: ifem


text:
function [ face_pointer, face_data ] = xyzf_data_read ( input_filename, ...
  face_num, face_data_num )

%*****************************************************************************80
%
%% XYZF_DATA_READ reads the data in an XYZF file.
%
%  Discussion:
%
%    This routine assumes that the file contains exactly three kinds of
%    records:
%
%    COMMENTS which begin with a '#' character in column 1;
%    BLANKS which contain nothing but 'whitespace';
%    FACE ITEMS, which are indices of points on a face.
%
%    The routine ignores comments and blanks and returns
%    the number of face items.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer FACE_DATA_NUM, the number of face items.
%
%    Output, integer FACE_POINTER(FACE_NUM+1), pointers to the
%    first face item for each face.
%
%    Output, integer FACE_DATA(FACE_DATA_NUM), the face items.
%
  input_unit = fopen ( input_filename, 'r' );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XYZF_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file "%s".\n', input_filename );
    error ( 'XYZF_DATA_READ - Error!' );
    return;
  end

  face = 0;
  face_pointer(1) = 1;
%
%  Use FGETL, not FGETS, so we don't pick up the NEWLINE character.
%
  while ( face < face_num )

    text = fgetl ( input_unit );

    if ( text(1) == '#' || s_len_trim ( text ) == 0 )
      continue
    end

    face = face + 1;

    n = s_word_count ( text );
    face_pointer(face+1) = face_pointer(face) + n;

    ilo = face_pointer(face);
    ihi = face_pointer(face+1) - 1;

    [ temp, count ] = sscanf ( text, '%d' );

    face_data(ilo:ihi) = temp(1:count);

  end
%
%  Close the file.
%
  fclose ( input_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function i4_ceiling_test ( )

%*****************************************************************************80
%
%% I4_CEILING_TEST tests I4_CEILING.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2014
%
%  Author:
%
%    John Burkardt
%
  r8_lo = -100.0;
  r8_hi =  100.0;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_CEILING_TEST\n' );
  fprintf ( 1, '  I4_CEILING evaluates the "ceiling" of an R8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      R8    I4_CEILING(R8)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ r8, seed ] = r8_uniform_ab ( r8_lo, r8_hi, seed );
    i4 = i4_ceiling ( r8 );
    fprintf ( 1, '  %8.4f            %4d\n', r8, i4 );
  end

  return
end


text:
function ds_info=my_gen_ds_info_gtea(ds_config)

ds_dir=fullfile('../datasets', 'gtea');

train_idx_file=fullfile(ds_dir, 'ImageSets/Segmentation/train.txt');
fid=fopen(train_idx_file);
train_file_names=textscan(fid, '%s');
train_file_names=train_file_names{1};
fclose(fid);

val_idx_file=fullfile(ds_dir, 'ImageSets/Segmentation/val.txt');
fid=fopen(val_idx_file);
val_file_names=textscan(fid, '%s');
val_file_names=val_file_names{1};
fclose(fid);

train_num=length(train_file_names);
img_names=cat(1, train_file_names, val_file_names);
img_num=length(img_names);

img_files=cell(img_num, 1);
mask_files=cell(img_num, 1);

for t_idx=1:img_num
    file_name=img_names{t_idx};
    img_files{t_idx}=[file_name(1:end-4) '.jpg'];
    mask_file_name = strrep(file_name,'JPEGImages','SegmentationClass');
    mask_files{t_idx}=[mask_file_name(1:end-4) '.png'];

end

train_idxes=1:train_num;
val_idxes=train_num+1:img_num;

ds_info=[];

ds_info.img_names=img_names;
ds_info.img_files=img_files;
ds_info.mask_files=mask_files;

ds_info.train_idxes=uint32(train_idxes');
ds_info.test_idxes=uint32(val_idxes');

disp('here::')
img_dir=fullfile(ds_dir, 'JPEGImages')
mask_dir=fullfile(ds_dir, 'SegmentationClass')

data_dirs=[];
data_dirs{1}=img_dir;
data_dirs{2}=mask_dir;

data_dir_idxes_img=zeros([img_num 1], 'uint8')+1;
data_dir_idxes_mask=zeros([img_num 1], 'uint8')+2;

ds_info.data_dir_idxes_img=data_dir_idxes_img;
ds_info.data_dir_idxes_mask=data_dir_idxes_mask;
ds_info.data_dirs=data_dirs;
ds_info.ds_dir=ds_dir;

ds_info.class_info=gen_class_info_ego();

ds_info.ds_name='gtea';

ds_info=process_ds_info_classification(ds_info, ds_config);

end

author: aurooj

repo: Hand-Segmentation-in-the-Wild


text:
function c=melcepst(s,fs,w,nc,p,n,inc,fl,fh)
%MELCEPST Calculate the mel cepstrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH)
%
%
% Simple use: (1) c=melcepst(s,fs)          % calculate mel cepstrum with 12 coefs, 256 sample frames
%			  (2) c=melcepst(s,fs,'e0dD')   % include log energy, 0th cepstral coef, delta and delta-delta coefs
%
% Inputs:
%     s	  speech signal
%     fs  sample rate in Hz (default 11025)
%     w   mode string (see below)
%     nc  number of cepstral coefficients excluding 0'th coefficient [default 12]
%     p   number of filters in filterbank [default: floor(3*log(fs)) =  approx 2.1 per ocatave]
%     n   length of frame in samples [default power of 2 < (0.03*fs)]
%     inc frame increment [default n/2]
%     fl  low end of the lowest filter as a fraction of fs [default = 0]
%     fh  high end of highest filter as a fraction of fs [default = 0.5]
%
%		w   any sensible combination of the following:
%
%               'R'  rectangular window in time domain
%				'N'	 Hanning window in time domain
%				'M'	 Hamming window in time domain (default)
%
%               't'  triangular shaped filters in mel domain (default)
%               'n'  hanning shaped filters in mel domain
%               'm'  hamming shaped filters in mel domain
%
%				'p'	 filters act in the power domain
%				'a'	 filters act in the absolute magnitude domain (default)
%
%               '0'  include 0'th order cepstral coefficient
%				'E'  include log energy
%				'd'	 include delta coefficients (dc/dt)
%				'D'	 include delta-delta coefficients (d^2c/dt^2)
%
%               'z'  highest and lowest filters taper down to zero (default)
%               'y'  lowest filter remains at 1 down to 0 frequency and
%			   	     highest filter remains at 1 up to nyquist freqency
%
%		       If 'ty' or 'ny' is specified, the total power in the fft is preserved.
%
% Outputs:	c     mel cepstrum output: one frame per row. Log energy, if requested, is the
%                 first element of each row followed by the delta and then the delta-delta
%                 coefficients.
%

% BUGS: (1) should have power limit as 1e-16 rather than 1e-6 (or possibly a better way of choosing this)
%           and put into VOICEBOX
%       (2) get rdct to change the data length (properly) instead of doing it explicitly (wrongly)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: melcepst.m 3497 2013-09-26 16:10:51Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 fs=11025; end
if nargin<3 w='M'; end
if nargin<4 nc=12; end
if nargin<5 p=floor(3*log(fs)); end
if nargin<6 n=pow2(floor(log2(0.03*fs))); end
if nargin<9
   fh=0.5;
   if nargin<8
     fl=0;
     if nargin<7
        inc=floor(n/2);
     end
  end
end

if isempty(w)
   w='M';
end
if any(w=='R')
   z=enframe(s,n,inc);
elseif any (w=='N')
   z=enframe(s,hanning(n),inc);
else
   z=enframe(s,hamming(n),inc);
end
f=rfft(z.');
[m,a,b]=melbankm(p,n,fs,fl,fh,w);
pw=f(a:b,:).*conj(f(a:b,:));
pth=max(pw(:))*1E-20;
if any(w=='p')
   y=log(max(m*pw,pth));
else
   ath=sqrt(pth);
   y=log(max(m*abs(f(a:b,:)),ath));
end
c=rdct(y).';
nf=size(c,1);
nc=nc+1;
if p>nc
   c(:,nc+1:end)=[];
elseif p<nc
   c=[c zeros(nf,nc-p)];
end
if ~any(w=='0')
   c(:,1)=[];
   nc=nc-1;
end
if any(w=='E')
   c=[log(max(sum(pw),pth)).' c];
   nc=nc+1;
end

% calculate derivative

if any(w=='D')
  vf=(4:-1:-4)/60;
  af=(1:-1:-1)/2;
  ww=ones(5,1);
  cx=[c(ww,:); c; c(nf*ww,:)];
  vx=reshape(filter(vf,1,cx(:)),nf+10,nc);
  vx(1:8,:)=[];
  ax=reshape(filter(af,1,vx(:)),nf+2,nc);
  ax(1:2,:)=[];
  vx([1 nf+2],:)=[];
  if any(w=='d')
     c=[c vx ax];
  else
     c=[c ax];
  end
elseif any(w=='d')
  vf=(4:-1:-4)/60;
  ww=ones(4,1);
  cx=[c(ww,:); c; c(nf*ww,:)];
  vx=reshape(filter(vf,1,cx(:)),nf+8,nc);
  vx(1:8,:)=[];
  c=[c vx];
end

if nargout<1
   [nf,nc]=size(c);
   t=((0:nf-1)*inc+(n-1)/2)/fs;
   ci=(1:nc)-any(w=='0')-any(w=='E');
   imh = imagesc(t,ci,c.');
   axis('xy');
   xlabel('Time (s)');
   ylabel('Mel-cepstrum coefficient');
	map = (0:63)'/63;
	colormap([map map map]);
	colorbar;
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%%*******************************************************************
%% HSDNTdirfun: compute (dX,dZ), given dy, for the NT direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

    function [par,dX,dy,dZ] = HSDNTdirfun(blk,At,par,Rd,EinvRc,xx);

    global solve_ok

    dX = cell(size(blk,1),1); dZ = cell(size(blk,1),1); dy = [];
    if (any(isnan(xx)) | any(isinf(xx)))
       solve_ok = 0;
       fprintf('\n  HSDNTdirfun: solution contains NaN or inf.');
       return;
    end
%%
    m = par.m;
    dy2 = xx(1:m+2);
%%
    for p=1:size(blk,1)
       pblk = blk(p,:);
       if strcmp(pblk{1},'l')
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),[],[],dy2));
          tmp   = par.dd{p}.*dZ{p};
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'q')
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),[],[],dy2));
          tmp = par.dd{p}.*dZ{p} + qops(pblk,qops(pblk,dZ{p},par.ee{p},1),par.ee{p},3);
          dX{p} = EinvRc{p} - tmp;
       elseif strcmp(pblk{1},'s')
          dZ(p) = ops(Rd(p),'-',Atyfun(pblk,At(p,:),par.permA(p,:),par.isspAy(p),dy2));
          tmp   = Prod3(pblk,par.W{p},dZ{p},par.W{p},1);
          dX{p} = EinvRc{p}-tmp;
       end
    end
    dy  = dy2(1:m);
    par.dtau = dy2(m+1);
    par.dtheta = dy2(m+2);
    par.dkap = (par.mu./par.tau - par.kap) - par.kap*(par.dtau/par.tau);
%%*******************************************************************

repo: REST


text:
function x_yuv = rgb2yuv(x_rgb)
r = x_rgb(:,:,1);
g = x_rgb(:,:,2);
b = x_rgb(:,:,3);
wr = 0.2126;
wb = 0.0722;
wg = 1-wr-wb;
umax = 0.436;
vmax = 0.615;
Y = wr*r +wg*g+wb*b;
x_yuv(:,:,1) = Y;
x_yuv(:,:,2) = umax*((b-Y)/(1-wb));
x_yuv(:,:,3) = vmax* ((r-Y)/(1-wr));
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function y=chk_param(val,name,varargin)

% y=chk_param(val,name,varargin)
%
% Check that 'val' satisfies all the requirements given in 'varargin'.
% If any of the requirements are not met, the parameter's 'name' is printed
% along with information on which requirement failed.
%
% Example usage:
%   error(chk_param(speed,'speed','>',0,'<=',100));
% The above will throw an error if speed is <=0 or >100. Otherwise,
% nothing will be printed.
%
% Supported parameter checks are:
%   'string'
%   TBD

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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/>.

global CHK_PARAM_BYPASS
if (CHK_PARAM_BYPASS==1)
  y=[];
  return
end

failed=0;
t=1;
while (t<=length(varargin))
  switch lower(varargin{t})
    case 'scalar'
      if (~isscalar(val)) failed=1; break; end
    case 'vector'
      if (~isvector(val)) failed=1; break; end
    case 'horizontal'
      if (size(val,1)~=1) failed=1; break; end
    case 'vertical'
      if (size(val,2)~=1) failed=1; break; end
    case 'string'
      if (~ischar(val)) failed=1; break; end
    case 'numeric'
      if (~isnumeric(val)) failed=1; break; end
    case 'struct'
      if (~isstruct(val)) failed=1; break; end
    case 'real'
      if (~isreal(val)) failed=1; break; end
    case '>'
      t=t+1;
      if (t>length(varargin)) error('> requires an argument'); end
      if (~all(val(:)>varargin{t})) failed=1; break; end
    case '<'
      t=t+1;
      if (t>length(varargin)) error('< requires an argument'); end
      if (~all(val(:)<varargin{t})) failed=1; break; end
    case '<='
      t=t+1;
      if (t>length(varargin)) error('<= requires an argument'); end
      if (~all(val(:)<=varargin{t})) failed=1; break; end
    case '>='
      t=t+1;
      if (t>length(varargin)) error('>= requires an argument'); end
      if (~all(val(:)>=varargin{t})) failed=1; break; end
    case 'rows>'
      t=t+1;
      if (t>length(varargin)) error('rows> requires an argument'); end
      if (size(val,1)<=varargin{t}) failed=1; break; end
    case 'rows<'
      t=t+1;
      if (t>length(varargin)) error('rows< requires an argument'); end
      if (size(val,1)>=varargin{t}) failed=1; break; end
    case 'rows<='
      t=t+1;
      if (t>length(varargin)) error('rows<= requires an argument'); end
      if (size(val,1)>varargin{t}) failed=1; break; end
    case 'rows>='
      t=t+1;
      if (t>length(varargin)) error('rows>= requires an argument'); end
      if (size(val,1)<varargin{t}) failed=1; break; end
    case 'rows=='
      t=t+1;
      if (t>length(varargin)) error('rows>= requires an argument'); end
      if (size(val,1)~=varargin{t}) failed=1; break; end
    case 'cols>'
      t=t+1;
      if (t>length(varargin)) error('cols> requires an argument'); end
      if (size(val,2)<=varargin{t}) failed=1; break; end
    case 'cols<'
      t=t+1;
      if (t>length(varargin)) error('cols< requires an argument'); end
      if (size(val,2)>=varargin{t}) failed=1; break; end
    case 'cols<='
      t=t+1;
      if (t>length(varargin)) error('cols<= requires an argument'); end
      if (size(val,2)>varargin{t}) failed=1; break; end
    case 'cols>='
      t=t+1;
      if (t>length(varargin)) error('cols>= requires an argument'); end
      if (size(val,2)<varargin{t}) failed=1; break; end
    case 'cols=='
      t=t+1;
      if (t>length(varargin)) error('cols== requires an argument'); end
      if (size(val,2)~=varargin{t}) failed=1; break; end
    case 'length>'
      t=t+1;
      if (t>length(varargin)) error('length> requires an argument'); end
      if (length(val)<=varargin{t}) failed=1; break; end
    case 'length<'
      t=t+1;
      if (t>length(varargin)) error('length< requires an argument'); end
      if (length(val)>=varargin{t}) failed=1; break; end
    case 'length<='
      t=t+1;
      if (t>length(varargin)) error('length<= requires an argument'); end
      if (length(val)>varargin{t}) failed=1; break; end
    case 'length>='
      t=t+1;
      if (t>length(varargin)) error('length>= requires an argument'); end
      if (length(val)<varargin{t}) failed=1; break; end
    case 'length=='
      t=t+1;
      if (t>length(varargin)) error('length== requires an argument'); end
      if (length(val)~=varargin{t}) failed=1; break; end
    case 'integer'
      if (any(val(:)~=floor(val(:)))) failed=1; break; end
    otherwise
      if (ischar(varargin{t}))
        error(sprintf('Unrecognized requirement passed to chk_param: %s',varargin{t}));
      else
        error(sprintf('Unrecognized requirement passed to chk_param index %i',t));
      end
  end
  t=t+1;
end

if (failed)
  y=['parameter ''' name ''' must satisfy the following requirements:'];
  for t=1:length(varargin)
    if (t~=1)
      y=[y ', '];
    else
      y=[y ' '];
    end
    if (ischar(varargin{t}))
      y=[y varargin{t}];
    else
      y=[y num2str(varargin{t})];
    end
  end
else
  y=[];
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function [inpts,volpts]=mrClearPoints(inpts,volpts)
%function [inpts,volpts]=mrClearPoints(inpts,volpts)
%
%Clear all alignment points

%6/16/96	gmb	wrote it.

ynstr=input('Clear all points (y/n)? ','s');

if ynstr(1)=='y'
	inpts=[];
	volpts=[];
end

npoints = size(volpts,1);
disp([num2str(npoints),' pairs of points selected.']);


text:
function x=isleap(Year)
%ISLEAP True for leap year.
%     ISLEAP(Year) returns 1 if Year is a leap year and 0 otherwise.
%     ISLEAP is only set for gregorian calendar, so Year >= 1583
%
% Syntax: 	ISLEAP(YEAR)
%
%     Inputs:
%           YEAR - Year of interest (default = current year).
%           You can input a vector of years.
%     Outputs:
%           Logical vector.
%
%      Example:
%
%           Calling on Matlab the function: isleap
%
%           Answer is: 0
%
%
%           Calling on Matlab the function: x=isleap([2007 2008])
%
%           Answer is:
%           x = 0 1
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%           Modified after Simon Jan suggestions
% To cite this file, this would be an appropriate format:
% Cardillo G. (2007) Isleap: a simple routine to test if a year is a leap
% year.
% http://www.mathworks.com/matlabcentral/fileexchange/14172

%Input Error handling
switch nargin
    case 0
         c=clock; Year=c(1); clear c
    case 1
        if ~isvector(Year) || ~all(isnumeric(Year)) || ~all(isfinite(Year)) || isempty(Year)
            error('Warning: Year values must be numeric and finite')
        end
        if ~isequal(Year,round(Year))
            error('Warning: Year values must be integer')
        end
        L=Year-1583;
        if L(L<0)
            error('Warning: Every value of Year must be >1582')
        end
    otherwise
        error('stats:Isleap:TooMuchInputs','Year must be a scalar or a vector.');
end

% The Gregorian calendar has 97 leap years every 400 years:
% Every year divisible by 4 is a leap year.
% However, every year divisible by 100 is not a leap year.
% However, every year divisible by 400 is a leap year after all.
% So, 1700, 1800, 1900, 2100, and 2200 are not leap years,
% but 1600, 2000, and 2400 are leap years.
x = ~mod(Year, 4) & (mod(Year, 100) | ~mod(Year, 400));
return

author: Sable

repo: mcbench-benchmarks


text:
function colored_noise_test01 ( n, q_d, alpha, seed_init )

%*****************************************************************************80
%
%% COLORED_NOISE_TEST01 calls F_ALPHA with particular parameters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of the sequence
%    to generate.
%
%    Input, real Q_D, the variance of the sequence.
%
%    Input, real ALPHA, the exponent of the power law.
%
%    Input, integer SEED_INIT, the initial seed for the
%    random number generator.
%
  output_filename = sprintf ( 'alpha_%4.2f.txt', alpha );
%
%  Report parameters.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COLORED_NOISE_TEST01:\n' );
  fprintf ( 1, '  Generating %d sample points.\n', n );
  fprintf ( 1, '  1/F^ALPHA noise has ALPHA = %f\n', alpha );
  fprintf ( 1, '  Variance is %f\n', q_d );
  fprintf ( 1, '  Initial random number seed = %d\n', seed_init );

  seed = seed_init;

  [ x, seed ] = f_alpha ( n, q_d, alpha, seed );
%
%  Print no more than 10 entries of the data.
%
  r8vec_print_part ( n, x, 10, '  Noise sample:' );
%
%  Write the data to a file.
%
  output_unit = fopen ( output_filename, 'wt' );

  for i = 1 : n
    fprintf ( output_unit, '%f\n', x(i) );
  end

  fclose ( output_unit );

  fprintf ( 1, '  Data written to file "%s".\n', output_filename );

  return
end


text:
function r = isintval(a)
%ISINTVAL     Returns 1 if  a  is intval
%
%   r = isintval(a)
%

% written  09/02/00     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  r = logical(1);

repo: OpenMAS


text:
function varargout = gui_receivedataset(varargin)
% GUI_RECEIVEDATASET MATLAB code for gui_receivedataset.fig
%      GUI_RECEIVEDATASET, by itself, creates a new GUI_RECEIVEDATASET or raises the existing
%      singleton*.
%
%      H = GUI_RECEIVEDATASET returns the handle to a new GUI_RECEIVEDATASET or the handle to
%      the existing singleton*.
%
%      GUI_RECEIVEDATASET('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUI_RECEIVEDATASET.M with the given input arguments.
%
%      GUI_RECEIVEDATASET('Property','Value',...) creates a new GUI_RECEIVEDATASET or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before gui_receivedataset_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to gui_receivedataset_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help gui_receivedataset

% Last Modified by GUIDE v2.5 19-Nov-2010 21:07:06

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @gui_receivedataset_OpeningFcn, ...
                   'gui_OutputFcn',  @gui_receivedataset_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before gui_receivedataset is made visible.
function gui_receivedataset_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to gui_receivedataset (see VARARGIN)

% update the popup contents (and datasets variables)
update_datasets(hObject);

% select dataset
handles = guidata(hObject);
if isempty([handles.datasets{2:2:end}])
    % no dataset specified: create a new one
    dataset = gui_loadset;
    % cancelled: cancel this dialog, too
    if isempty(dataset)
        return; end
    % (re-)create the dataset list
    update_datasets(hObject);
end

% preferably select the lastdata in the pulldown menu
handles = guidata(hObject);
dataids = strtok(cellstr(get(handles.popupmenu1,'String')));
selected_id = find(strcmp(dataids,'lastdata'),1);
if isempty(selected_id)
    selected_id = 1; end
set(handles.popupmenu1,'Value',selected_id);
guidata(hObject, handles);
% invoke the popup menu to propagate selection
popupmenu1_Callback(handles.popupmenu1,{},guidata(hObject));

% Choose default command line output for gui_receivedataset
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes gui_receivedataset wait for user response (see UIRESUME)
uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = gui_receivedataset_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
try
    varargout{1} = handles.output;
    close(handles.figure1);
catch
    varargout = {};
end

function update_datasets(hObject)
handles = guidata(hObject);
% retrieve and attach the known datasets...
[handles.datasets,printed,handles.indexable_datasets] = gui_listdatasets;
guidata(hObject,handles);
set(handles.popupmenu1,'String',printed);

function edit1_Callback(hObject, eventdata, handles)
gui_check_scalar(hObject,'25');

function edit2_Callback(hObject, eventdata, handles)
gui_check_varname(hObject,'laststream');

function pushbutton1_Callback(hObject, eventdata, handles)
% get dataset name
choices = get(handles.popupmenu1,'String');
selection = get(handles.popupmenu1,'Value');
dataset = strtok(choices{selection});
% start process
run_receivedataset(get(handles.edit2,'String'),dataset,str2num(get(handles.edit1,'String')));
% resume GUI
uiresume(handles.figure1);

function pushbutton2_Callback(hObject, eventdata, handles)
% cancel pressed
uiresume(handles.figure1);

function pushbutton3_Callback(hObject, eventdata, handles)
% help pressed
disp('coming soon...');

% --- auto-generated junk ---

function popupmenu1_Callback(hObject, eventdata, handles)
function popupmenu1_CreateFcn(hObject, eventdata, handles)
function edit1_CreateFcn(hObject, eventdata, handles)
function edit2_CreateFcn(hObject, eventdata, handles)


text:
function varargout = createRhombododecahedron()
%CREATERHOMBODODECAHEDRON Create a 3D mesh representing a rhombododecahedron.
%
%   [V, E, F] = createRhombododecahedron
%   V is a 14-by-3 array with vertex coordinate,
%   E is a 12-by-2 array containing indices of neighbour vertices,
%   F is a 8-by-3 array containing vertices array of each face.
%
%   [V, F] = createRhombododecahedron;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createRhombododecahedron;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%   [v, e, f] = createRhombododecahedron;
%   drawMesh(v, f);
%
%
%   See also
%   meshes3d, drawMesh

%   ---------
%   author : David Legland
%   e-mail: david.legland@inra.fr
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/02/2005.
%

%   HISTORY
%   04/01/2007: remove unused variables

nodes = [0 0 2;...
    1 -1 1;1 1 1;-1 1 1;-1 -1 1;...
    2 0 0;0 2 0;-2 0 0;0 -2 0;...
    1 -1 -1;1 1 -1;-1 1 -1;-1 -1 -1;...
    0 0 -2];

edges = [...
    1 2;1 3;1 4;1 5;...
    2 6;2 9;3 6;3 7;4 7;4 8;5 8;5 9;...
    6 10;6 11;7 11;7 12;8 12;8 13;9 10;9 13; ...
    10 14;11 14;12 14;13 14];

faces = [...
    1 2 6 3;...
    1 3 7 4;...
    1 4 8 5;...
    1 5 9 2;...
    2 9 10 6;...
    3 6 11 7;...
    4 7 12 8;...
    5 8 13 9;...
    6 10 14 11;...
    7 11 14 12;...
    8 12 14 13;...
    9 13 14 10];

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef ParadigmTFR < ParadigmDataflowSimplified
    % Time/Frequency Regression. This is a new approach that is not yet published.
    %
    % This method learns a linear classifier (or regressor) of second-order dynamics in the EEG over
    % a particular set of time and frequency bins. The key trick lies in a set of extra assumptions
    % that make the (very high-dimensional) solution statistically tractable. The algorithm amounts
    % to a single convex optimization problem dependent on a small number of hyper-parameters that can
    % be optimized using grid search. For a reasonably exhaustive grid the solution is the globally
    % optimal second-order dynamics solution.
    %
    % Notes:
    %  This implementation is not yet optimally tuned -- the spectral estimation should be replaced by a
    %  similar approach to multi-taper CSP.
    %
    % Examples:
    %   % learn a spectral classifier for features within -2 to 2 seconds relative to some marker in some standard EEG frequency bands
    %   % here using area under curve (AUC) to optimize the regularization parameter (assuming that the number of exemplars per class is imbalanced)
    %   % and setting the set of relative regularization term weights that should be searched over to {[1 1 1 1]}, i.e., not searching over those, for speed.
    %   myapproach = {'TFR' 'SignalProcessing',{'EpochExtraction',[-2 2]},'Prediction',{'FeatureExtraction',{'Times',[-1.5,-1,-0.5,0,0.5,1,1.5],'Frequencies',[4 7 11 14 20 30]}, ...
    %       'MachineLearning',{'Learner',{'proximal','LambdaSearch',{'ParameterMetric','auc'},'TermWeights',{[1 1 1 1]}}}}};
    %
    % Name:
    %   Time/Frequency Regression
    %
    %                            Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                            2013-03-08

    methods

        function defaults = preprocessing_defaults(self)
            % define the default pre-processing parameters of this paradigm
            defaults = {'FIRFilter',{[0.5 1],'highpass'}, 'EpochExtraction',[-3 3], 'Resampling',100};
        end

        function defaults = machine_learning_defaults(self)
            % global sharing approach
            defaults = {'proximal', 'lambdaSearch',{'lambdas',2.^(5:-0.1:-4),'foldmargin',0},'regularizers',{ ...
                    'term1',{'trace'}, ...                                                                      % the weights are a small sum of second-order spatial filters per T/F resel
                    'term2',{'trace','LinearOperator','@(x)reshape(x,a*b,c*d)'}, ...                            % across all of time/frequency we learn combinations of few latent spatial filters
                    'term3',{'l2', 'LinearOperator','@(x)vec(diff(x,[],4))','NonorthogonalTransform',true}, ... % temporal smoothness
                    'term4',{'l2', 'LinearOperator','@(x)vec(diff(x,[],3))','NonorthogonalTransform',true}, ... % spectral smoothness
                }, 'regweights',{[1 1 1 1],[1 1 0.5 0.5],[1 1 2 2],[1 1 0.25 0.25],[1 1 4 4],[1 2 1 2],[2 1 1 2],[1 2 2 1],[2 1 2 1]}};
        end

        function model = feature_adapt(self,varargin)
            % adapt a feature representation using the CSP algorithm
            args = arg_define(varargin, ...
                arg_norep('signal'), ...
                arg({'freqpoints','Frequencies'},[4 7 11 14 20 30],[0 0.2 200 1000],'Frequency points to consider.'),...
                arg({'timepoints','Times'},-2.5:0.5:2.5,[],'Time points to consider.'), ...
                arg({'from_edges','FromEdges'},true,[],'Run up to epoch edges. This is both in time and frequency.'), ...
                arg({'vectorize_features','VectorizeFeatures'},false,[],'Vectorize the features. For compatibility with basic classifiers.'));
            model.args = rmfield(args,'signal');
            model.chanlocs = args.signal.chanlocs;
        end

        function [features,shape] = feature_extract(self,signal,featuremodel)
            % first perform the time/freq decomposition
            args = featuremodel.args;
            if args.from_edges
                args.timepoints = [signal.xmin args.timepoints signal.xmax];
                args.freqpoints = [0 args.freqpoints signal.srate/2];
            end
            time2idx = @(t) min(signal.pnts,max(1,1+round((t-signal.xmin)*signal.srate)));
            half_hann = @(a,b) 0.5*(1-cos(pi*((((a+1):b) - a) / (b - a))));
            shape = [size(signal.data,1),size(signal.data,1),length(args.freqpoints)-2,length(args.timepoints)-2];
            features = zeros([shape signal.trials]);
            for t=2:length(args.timepoints)-1
                % calculate time window function
                last_tp = time2idx(args.timepoints(t-1));
                cur_tp = time2idx(args.timepoints(t));
                next_tp = time2idx(args.timepoints(t+1));
                wndrange = last_tp+1:next_tp;
                wndfunc = [half_hann(last_tp,cur_tp) 1-half_hann(cur_tp,next_tp)];
                % window and fourier-transform the data
                X = fft(bsxfun(@times,wndfunc,signal.data(:,wndrange,:)),[],2);
                % multiply out the cross-spectral covariance
                tmp = zeros(size(X,1),size(X,1),ceil(size(X,2)/2),size(X,3));
                freq2idx = @(f) min(size(tmp,3),max(1,1+round(f*length(wndrange)/signal.srate)));
                for f=1:size(tmp,3)
                    for n=1:size(tmp,4)
                        tmp(:,:,f,n) = 2*real(X(:,f,n)*X(:,f,n)');
                    end
                end
                % average around frequency centers
                for f=2:length(args.freqpoints)-1
                    last_fp = freq2idx(args.freqpoints(f-1));
                    cur_fp = freq2idx(args.freqpoints(f));
                    next_fp = freq2idx(args.freqpoints(f+1));
                    wndrange = last_fp+1:next_fp;
                    wndfunc = [half_hann(last_fp,cur_fp) 1-half_hann(cur_fp,next_fp)];
                    features(:,:,f-1,t-1,:) = mean(bsxfun(@times,tmp(:,:,wndrange,:),reshape(wndfunc,1,1,[])),3);
                end
            end
            % do final vectorization if desired
            if featuremodel.args.vectorize_features
                features = reshape(features,[],signal.trials)'; end
        end

        function layout = dialog_layout_defaults(self)
            layout = {'SignalProcessing.Resampling.SamplingRate','SignalProcessing.FIRFilter.Frequencies', ...
                '', 'SignalProcessing.EpochExtraction','', ...
                'Prediction.FeatureExtraction.Frequencies', 'Prediction.FeatureExtraction.Times', ...
                'Prediction.FeatureExtraction.FromEdges','', ...
                'Prediction.MachineLearning.Learner.LossType', '', ...
                'Prediction.MachineLearning.Learner.LambdaSearch.Lambdas', ...
                'Prediction.MachineLearning.Learner.TermWeights', ...
                'Prediction.MachineLearning.Learner.LambdaSearch.ParameterMetric'};
        end

    end
end


text:
function f = p42_f ( n, x )

%*****************************************************************************80
%
%% P42_F evaluates the objective function for problem 42.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2002
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    MJD Powell,
%    An Efficient Method for Finding the Minimum of a Function of
%    Several Variables Without Calculating Derivatives,
%    Computer Journal,
%    Volume 7, Number 2, pages 155-162, 1964.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  if ( x(2) == 0.0 )
    term = 0.0;
  else
    arg = ( x(1) + 2.0 * x(2) + x(3) ) / x(2);
    term = exp ( - arg^2 );
  end

  f = 3.0 ...
    - 1.0 / ( 1.0 + ( x(1) - x(2) )^2 ) ...
    - sin ( 0.5 * pi * x(2) * x(3) ) ...
    - term;

  return
end