text:
function eddy = mod_new_eddy(STATS, amplitude, lat, lon, thresh, sa, cyc, geospeed, extrema_num, detect)
%NEW_EDDY Initializes new eddy objects, run with no arguments to create an
%empty matrix
    if nargin
        eddy = struct('Stats', STATS, ...
            'Lat', lat, ...
            'Lon', lon, ...
            'Amplitude', amplitude, ...
            'ThreshFound', thresh, ...
            'SurfaceArea', sa, ...
            'Date', NaN, ...
            'Cyc', cyc, ...
            'MeanGeoSpeed', geospeed, ...
            'NumberOfExtrema', extrema_num, ...
            'DetectedBy', detect);
    else
        eddy = struct('Stats', {}, ...
            'Lat', {}, ...
            'Lon', {}, ...
            'Amplitude', {}, ...
            'ThreshFound', {}, ...
            'SurfaceArea', {}, ...
            'Date', {}, ...
            'Cyc', {}, ...
            'MeanGeoSpeed', {}, ...
            'NumberOfExtrema', {}, ...
            'DetectedBy', {});
    end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
classdef MatNatSession < MatNatBase
    % MatNatSession A Matlab class representing an XNAT session
    %
    %     Licence
    %     -------
    %     Part of MatNat. https://github.com/tomdoel/matnat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (SetAccess = protected)
        Label
        Id
        ProjectId
        SubjectLabel
    end

    properties (Access = private)
        RestClient
        ScanMap
    end

    methods
        function obj = MatNatSession(restClient)
            obj.RestClient = restClient;
        end

        function scanMap = getScanMap(obj)
            % Returns a map of scan IDs to MatNatScan objects

            obj.populateScanMapIfNecessary;
            scanMap = obj.ScanMap;
        end

        function number_of_scans = CountScans(obj)
            % Returns the number of scans in this session

            obj.populateScanMapIfNecessary;
            number_of_scans = obj.ScanMap.Count;
        end
    end

    methods (Access = private)
        function populateScanMapIfNecessary(obj)
            if isempty(obj.ScanMap)
                obj.ScanMap = obj.RestClient.getScanMap(obj.ProjectId, obj.SubjectLabel, obj.Label);
            end
        end
    end

    methods (Static)
        function obj = createFromServerObject(restClient, serverObject, projectId, subjectLabel)
            % Creates a MatNatSession based on the information
            % structure returned from the XNAT server

            obj = MatNatSession(restClient);
            obj.ProjectId = projectId;
            obj.SubjectLabel = subjectLabel;
            obj.Label = MatNatBase.getOptionalProperty(serverObject, 'label');
            obj.Id = MatNatBase.getOptionalProperty(serverObject, 'ID');
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function test_pull2099

% WALLTIME 00:20:00
% MEM 3gb
% DEPENDENCY ft_multiplotER ft_singleplotER ft_topoplotER lineattributes_common ft_databrowser

%%
load(dccnpath('/home/common/matlab/fieldtrip/data/test/pull2099.mat'));

% try a bunch of visualizations
cfg = [];
cfg.layout = 'CTF275_helmet.mat';
ft_multiplotER(cfg, t);

cfg = [];
cfg.linecolor = 'spatial';
cfg.layout = 'butterfly';
ft_multiplotER(cfg, t); % ->monochrome because spatial information of the sensors is not known (grad missing)

cfg = [];
cfg.linecolor = 'spatial';
cfg.layout = 'CTF275_helmet.mat';
cfg.linewidth = 2;
ft_multiplotER(cfg, t);

cfg = [];
cfg.linecolor = 'spatial';
cfg.viewmode = 'butterfly';
cfg.layout = 'CTF275_helmet.mat';
ft_multiplotER(cfg, t);

cfg = [];
cfg.viewmode = 'butterfly';
cfg.layout = 'CTF275_helmet.mat';
cfg.colorgroups = 'labelchar3';
ft_multiplotER(cfg, t);

cfg = [];
cfg.channel = 'MRO';
cfg.showlocations = 'yes';
cfg.linecolor = [0.8 0.4 0.3];
cfg.layout = 'CTF275_helmet.mat';
ft_singleplotER(cfg, t);

cfg = [];
cfg.channel = 'MRO';
cfg.showlocations = 'yes';
cfg.layout = 'CTF275_helmet.mat';
ft_singleplotER(cfg, t1, t2, t3, t4);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Ri = get_Ri(r_i, theta_i, sigma_r, sigma_theta)
% computes the covariance of the point p_i = r_i*[cos(theta_i); sin(theta_i)
% in cartesian coordinates

% first we compute the matrix square root of R, which consists of the
% Jacobian * sqrt(R_in_polar_coords)
Ri = [cos(theta_i) -r_i*sin(theta_i)
      sin(theta_i)  r_i*cos(theta_i) ] * diag([sigma_r sigma_theta]);

% and then we get the entire Ri
Ri = Ri * Ri';

author: rpng

repo: ocekf-slam

save_path: github-repos/MATLAB/rpng-ocekf-slam


text:
function [A,C] = fd_bilinear_coefficients(mn,mx,side,V)
  % FD_BILINEAR_COEFFICIENTS  Given a grid from minimum corner mn to maximum
  % corner mx, with `side` nodes on each side, compute a matrix `A` so that
  % `V=A*C` if `C` are the node centers.
  %
  % [A,C] = fd_bilinear_coefficients(mn,mx,side,V)
  %
  % Inputs:
  %   mn  3d position of minimum corner
  %   mx  3d position of maximum corner
  %   side  number of nodes in each dimension [x y]
  %   V  #V by 2 list of query points
  % Outputs:
  %   A  x*y by #V matrix
  %   C  x*y by 3 list of node center positions
  %

  [X,Y] = meshgrid( ...
    linspace(mn(1),mx(1),side(1)), ...
    linspace(mn(2),mx(2),side(2)));
  C = [X(:) Y(:)];
  r = (mx-mn)./(side-1);

  I = floor(bsxfun(@times,bsxfun(@rdivide,bsxfun(@minus,V,mn),mx-mn),side-1));
  A = sparse(size(V,1),prod(side));
  for x = 1:2
    for y = 1:2
      Ixy = sub2ind(side([2 1]),I(:,2)+x,I(:,1)+y);
      Sxy = C(Ixy,:);
      Axy = prod(bsxfun(@minus,r,abs(V-Sxy)),2)/prod(r);
      A = A + sparse((1:size(V,1))',Ixy,Axy,size(V,1),size(C,1));
    end
  end

end

author: alecjacobson

repo: gptoolbox


text:
% Copyright 2017 Lime Microsystems Ltd.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%    http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
% Page indication channel
% requires STTD?
% 288 bits/frame in 15 slots, last 12 bits not used
% PI relates to system frame number SFN
% table 4.3 and section 4.4.5 of Richardson W-CDMA book (ISBN 978-0-521-18782-4)
function y=WCDMADLpich(sfn,dly) % delay is 30x256 chips earlier than sccpch
	y=[((1-2*randint(1,144))+i*(1-2*randint(1,144))),zeros(1,6)]/sqrt(2); % last 6x2 bits unused
	y=shift( y, dly );
end

author: myriadrf


text:
function writeBB(filename,D,type)
  % WRITEBB  write medit data file
  %
  % writeBB(filename,D,type)
  %
  % Inputs:
  %   filename  path to *.bb file
  %   D  #D by dim  data matrix
  %   type  type identifier as per
  %     http://www.ann.jussieu.fr/frey/publications/RT-0253.pdf page 34
  %     2 means defined at vertices
  %     1 means defined at elements

  f = fopen(filename,'w');
  % print header
  % only vector data is supported
  fprintf(f,'%d %d %d %d\n',2,1,size(D,1),type);
  fprintf(f,'%g\n',D');
  fclose(f);

end


text:
function [M_n,Threshold_graph,H_est_time,RX_Payload_1_no_Equalizer,RX_Payload_2_no_Equalizer,RX_Payload_1_no_pilot,RX_Payload_2_no_pilot,BER] = OFDM_RX(RX,Parameters_struct)
%% Debug mode
Debug_mode = 'off';
if strcmp(Debug_mode,'on')
   clearvars -except Debug_mode;close all;clc;
   Global_Parameters;
   load('RX');
end
%% j Parameter
j = 1i;
%% Root Raised Cosine filter
rolloff = 0.5;
L_RRC = 6;
OVR = 2;
RRC = rcosdesign(rolloff,L_RRC,OVR,'sqrt'); % [1x13]
RX_signal = conv(RX,RRC); % [1x3012]
%% Packet Detection
D = 16;
L = 32;
C_n = zeros(1,length(RX)-D+1-L);
P_n = zeros(1,length(RX)-D+1-L);
C_k = zeros(1,L);
P_k = zeros(1,L);

for n=1:length(RX)-D+1-L
    for k=1:L
        C_k(k) = RX(n+k-1)*complex(RX(n+k-1+D));
        P_k(k) = abs(RX(n+k-1+D))^2;
    end
    C_n(n) = sum(C_k);
    P_n(n) = sum(P_k);
end
M_n = (abs(C_n).^2)./(P_n.^2);
%% Packet_select
Threshold = 0.75;
loc = find(M_n>Threshold);
temp_1 = [loc,0];
temp_2 = [0,loc];
temp_3 = temp_1-temp_2;
Packet_Front = find(temp_3>300);
Packet_Front_idx = loc(Packet_Front);
Length_over_Threshold = 230;

for x=1:length(Packet_Front_idx)-1
    if M_n(Packet_Front_idx(x)+Length_over_Threshold)>Threshold;
        idx = Packet_Front_idx(x)+L_RRC+1;
    end % if Loop
end % for Loop
Threshold_graph = Threshold*ones(1,length(M_n));
Threshold_graph(idx-L_RRC-1) = 1.15;
%% Downsampling
Frame_DWN_sampling = RX_signal(idx:OVR:OVR*480+idx-1); % [1x480] Frame length
%% Coarse CFO Estimation
Short_preamble_slot_length = 16;
z = Frame_DWN_sampling(Short_preamble_slot_length*5+1:Short_preamble_slot_length*6)*Frame_DWN_sampling(Short_preamble_slot_length*6+1:Short_preamble_slot_length*7)'; % [1x16]*[16x1]
f_Coarse_est = (-1/(2*pi*Short_preamble_slot_length*Parameters_struct.Ts))*angle(z);
Frame_After_Coarse = Frame_DWN_sampling.*exp(-j*2*pi*f_Coarse_est*Parameters_struct.Ts*(0:480-1)); % [1x480]
%% Fine CFO Estimation
z = Frame_After_Coarse(Short_preamble_slot_length*12+1:Short_preamble_slot_length*16)*Frame_After_Coarse(Short_preamble_slot_length*16+1:Short_preamble_slot_length*20)'; % [1x64]*[64x1]=[1x1]
f_Fine_est = (-1/(2*pi*64*Parameters_struct.Ts))*angle(z);
Frame_After_Fine = Frame_After_Coarse.*exp(-j*2*pi*f_Fine_est*Parameters_struct.Ts*(0:480-1)); % [1x160]
%% Symbol Timing Estimation
%% Channel Estimation
Long_preamble_1 = Frame_After_Fine(Short_preamble_slot_length*12+1:Short_preamble_slot_length*16); % [1x64]
Long_preamble_2 = Frame_After_Fine(Short_preamble_slot_length*16+1:Short_preamble_slot_length*20); % [1x64]
Long_preamble_1_After_FFT = fftshift(fft(Long_preamble_1)); % [1x64]
Long_preamble_2_After_FFT = fftshift(fft(Long_preamble_2)); % [1x64]
H_est = 0.5*(Long_preamble_1_After_FFT+Long_preamble_2_After_FFT).*conj(Parameters_struct.Long_preamble_slot_Frequency); % [1x64]
H_est_time = ifft(ifftshift(H_est)); % [1x64]
%% One tap Equalizer
RX_Payload_1_time = Frame_After_Fine(320+1:400); % [1x80]
RX_Payload_1_no_CP = RX_Payload_1_time(17:end); % [1x64]
RX_Payload_1_Frequency = fftshift(fft(RX_Payload_1_no_CP)); % [1x64]
RX_Payload_1_Frequency_Equalizer = RX_Payload_1_Frequency./H_est; % [1x64]

RX_Payload_2_time = Frame_After_Fine(400+1:480); % [1x80]
RX_Payload_2_no_CP = RX_Payload_2_time(17:end); % [1x64]
RX_Payload_2_Frequency = fftshift(fft(RX_Payload_2_no_CP)); % [1x64]
RX_Payload_2_Frequency_Equalizer = RX_Payload_2_Frequency./H_est; % [1x64]
%% De-Mapping
RX_Payload_1_no_Equalizer = [RX_Payload_1_Frequency(7:11),RX_Payload_1_Frequency(13:25),RX_Payload_1_Frequency(27:32),RX_Payload_1_Frequency(34:39),RX_Payload_1_Frequency(41:53),RX_Payload_1_Frequency(55:59)]; % [1x48]
RX_Payload_1_no_pilot = [RX_Payload_1_Frequency_Equalizer(7:11),RX_Payload_1_Frequency_Equalizer(13:25),RX_Payload_1_Frequency_Equalizer(27:32),RX_Payload_1_Frequency_Equalizer(34:39),RX_Payload_1_Frequency_Equalizer(41:53),RX_Payload_1_Frequency_Equalizer(55:59)]; % [1x48]
RX_Payload_1_Final = pskdemod(RX_Payload_1_no_pilot,4,pi/4); % [1x48]

RX_Payload_2_no_Equalizer = [RX_Payload_2_Frequency(7:11),RX_Payload_2_Frequency(13:25),RX_Payload_2_Frequency(27:32),RX_Payload_2_Frequency(34:39),RX_Payload_2_Frequency(41:53),RX_Payload_2_Frequency(55:59)]; % [1x48]
RX_Payload_2_no_pilot = [RX_Payload_2_Frequency_Equalizer(7:11),RX_Payload_2_Frequency_Equalizer(13:25),RX_Payload_2_Frequency_Equalizer(27:32),RX_Payload_2_Frequency_Equalizer(34:39),RX_Payload_2_Frequency_Equalizer(41:53),RX_Payload_2_Frequency_Equalizer(55:59)]; % [1x48]
RX_Payload_2_Final = pskdemod(RX_Payload_2_no_pilot,4,pi/4); % [1x48]
%% BER calculation
Error_bits = sum([abs(sign(Parameters_struct.data_Payload_1-RX_Payload_1_Final)),abs(sign(Parameters_struct.data_Payload_2-RX_Payload_2_Final))]);
BER = Error_bits/(length(Parameters_struct.data_Payload_1)+length(Parameters_struct.data_Payload_2));
%% Plot
if strcmp(Debug_mode,'on')
    subplot(2,4,1),plot(RX,'.');title('RX-Raw');axis([-1.5 1.5 -1.5 1.5]);axis square;
    %--------------------------------------------------------------------------------%
    subplot(2,4,2),plot(real(RX));title('I');axis([1 3000 -1.5 1.5]);axis square;
    subplot(2,4,3),plot(imag(RX));title('Q');axis([1 3000 -1.5 1.5]);axis square;
    %--------------------------------------------------------------------------------%
    [Spectrum_waveform,Welch_Spectrum_frequency] = pwelch(RX,[],[],[],1/Parameters_struct.Ts,'centered','power');
    subplot(2,4,4),plot(Welch_Spectrum_frequency,pow2db(Spectrum_waveform));
    title('Welch Power Spectral Density');axis square;
    %--------------------------------------------------------------------------------%
    subplot(2,4,5),plot(1:length(M_n),M_n,1:length(M_n),Threshold_graph);title('Packet Detection');axis([1,length(M_n),0,1.2]);axis square;
    subplot(2,4,6),plot(abs(H_est_time));title('Channel Estimation');axis([1 64 0 7]);axis square;xlabel('Time');
    %--------------------------------------------------------------------------------%
    subplot(2,4,7),plot(RX_Payload_1_no_Equalizer,'*');
    hold on
    subplot(2,4,7),plot(RX_Payload_2_no_Equalizer,'*');
    title('Before Equalizer');axis([-8 8 -8 8]);axis square;
    hold off
    %--------------------------------------------------------------------------------%
    subplot(2,4,8),plot(RX_Payload_1_no_pilot,'*');
    hold on
    subplot(2,4,8),plot(RX_Payload_2_no_pilot,'*');
    title({'Demodulation';['BER = ',num2str(BER)]});axis([-1.5 1.5 -1.5 1.5]);axis square;
    hold off
    set(gcf,'Units','centimeters','position',[1 2 49 24]);
end % Plot end
%% End function
end

author: MeowLucian


text:
function [equality center] = doquality(Xcx,X)
% :Usage:
% ::
%
%     [equality] = doquality(Xcx,X)
%
% :Inputs:
%
%   **Xcx:**
%        binary indicator matrix of cluster assignments,
%
%   **X:**
%        stimulus coordinates in group space
%
%        also: takes group spaces with zeros;

if size(Xcx>1);
    for i = 1:size(Xcx,2)    % for each class
        tmp = mean(X(find(Xcx(:,i)),:),1);     % get center of this class
        center(i,:) = tmp;
        % dist of all points in X from each center
        % number of cols of X is the number of dimensions; sums squared
        % vals across dims, takes sqrt to get Euclidean distance in N-d
        % space
        d(:,i) = sum((repmat(tmp,size(X,1),1) - X).^2,2).^0.5;
    end

    % rows of d are objects, columns are classes, values are dist from
    % class center
    for i=1:size(d,1)                   % i is the object (point)

        myclass = find(Xcx(i,:)==1);    % index of which class it is

        edist(i) = d(i,myclass);        %distance to center of own class

        otherclass = find(Xcx(i,:)==0);     % indices of columns for other classes

        otherdist(i) = min(d(i,otherclass));  %distance to center of cluster

    end

    % for each point, quality is distance to nearest neighbor - dist to own
    % cl / mak of those two
    equality = (otherdist - edist) ./ max([edist;otherdist]);

else
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function putlabel(x,y,names,z,znames)

%PUTLABEL plots user-specified labels to the observations in a two- or
%  a three-dimensional figure.
%  If no labels are given, the indices are plotted.
%
% Required input arguments:
%    x : x-coordinates of the data
%    y : y-coordinates of the data
%
% Optional input arguments:
%   names : labels to be added on the plot. They must be listed in a
%           column vector.
%       z : z-coordinates of the data
%  znames : labels to be added on the 3D-plot. They must be listed in a
%           columnvector.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% I/O: putlabel(x,y,names,z,znames)
%
% Written by S. Verboven on 01/10/2002
% Last update on 18/02/2004

xrange=get(gca,'Xlim');
range=xrange(2)-xrange(1);
if nargin<3
   for i=1:length(x)
      text(x(i)+range/50,y(i),num2str(i));
   end
else
   if nargin<4
      for i=1:length(x)
         text(x(i)+range/50,y(i),names(i,:));
      end
   else
      if nargin<5
         for i=1:length(x)
            text(x(i)+range/50,y(i),z(i),num2str(i));
         end
      else
         for i=1:length(x)
            text(x(i)+range/50,y(i),z(i),znames(i,:));
         end
      end
   end
end

author: marianux

repo: ecg-kit

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


text:
% spdiag - sparse diagonal matrix
%
% Copyright(c) 2009 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt

function D = spdiag(d)

if isempty(d)
  D = [];
  return;
end

if size(d,1)<size(d,2)
  d = d';
end

D = spdiags(d,0,size(d,1),size(d,1));

repo: REST


text:
% IndexToAssignment Convert index to variable assignment.
%
%   A = IndexToAssignment(I, D) converts an index, I, into the .val vector
%   into an assignment over variables with cardinality D. If I is a vector,
%   then the function produces a matrix of assignments, one assignment
%   per row.
%
%   See also AssignmentToIndex.m and FactorTutorial.m

function A = IndexToAssignment(I, D)

D = D(:)'; % ensure that D is a row vector
A = mod(floor(repmat(I(:) - 1, 1, length(D)) ./ repmat(cumprod([1, D(1:end - 1)]), length(I), 1)), ...
        repmat(D, length(I), 1)) + 1;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function rts_retrieveDataAndStore(tcpipClient, numData, dataType)
%rts_retrieveDataAndStore Summary of this function goes here
%   Detailed explanation goes here
    readasync(tcpipClient);

    bytesDbl = 8;
    dataPts = [];
    frameName = [];
    frameNum = [];
    transmitDateTime = [];

    try
        if(tcpipClient.BytesAvailable >= 88) %88 bytes is the size of the frame header
            while(tcpipClient.BytesAvailable > 88)
                [frameName, frameNum, transmitDateTime] = rts_readMinorFrameHeader(tcpipClient);
                if(strcmpi(dataType,'d'))
                    if(numData == -1)
                        dataPts = fread(tcpipClient,floor(tcpipClient.BytesAvailable/bytesDbl),'double');
                    else
                        numDataToRet = min(numData, floor(tcpipClient.BytesAvailable/bytesDbl));
                        if(numDataToRet > 0)
                            dataPts = fread(tcpipClient,numDataToRet,'double');
                        else
                            dataPts = [];
                        end
                    end
                elseif(strcmpi(dataType,'s'))
                    dataPts = deblank(fscanf(tcpipClient,'%c',numData));
                end
            end

            retDateTime = datevec(now);
            data = {dataPts, frameName, frameNum, transmitDateTime', retDateTime};
            set(tcpipClient,'UserData',data);
        end
    catch ME
        disp(ME.message);
    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% I'm sure this one is not mine; I think it was Fran who wrote it.
% Last checked: June 2021

clc
% define parameters
c0 = 1;

% define a 1D mesh
L = 1;
Nr = 1000;
mesh1 = createMeshSpherical1D(Nr, L);
x = mesh1.cellcenters.x;
xf = mesh1.facecenters.x;

% define the boundaries
BC = createBC(mesh1); % all Neumann
BC.right.a=0;
BC.right.b=1;
BC.right.c=0;
[Mbc, RHSbc] = boundaryCondition(BC);

% define diffusion term
D = createCellVariable(mesh1,1);
Dave = harmonicMean(D);
Mdiff = diffusionTermSpherical1D(Dave);

% initial values
c_old = createCellVariable(mesh1, c0);
c = c_old;

% solver
dt = 1E-3;
t = dt;
N = 30;

err = 0;
tvec = dt;

it = 1;

while t < 1
    % define the transient term
    [Mt, RHSt] = transientTerm(c_old, dt);

    % assemble the matrix and vector
    M = Mt-Mdiff+Mbc;
    RHS = RHSt+RHSbc;

    % solve PDE
    c = solvePDE(mesh1,M,RHS);
    c_old = c;

    % calculate analytic solution
    f = zeros(Nr,1);
    for n=1:N
        f = f + 2*(-1)^(n+1)/(n*pi)*exp(-(n*pi)^2*t)*sin(n*pi*x)./x;
    end

    tvec = [tvec; t];
    err = [err; trapz(x,(f-c.value(2:Nr+1)).^2)];

    if mod(it,10) == 0
    figure(1);
    plot(x, f, 'b', [0; x; 1], c.value,'ro');
    title(sprintf('t=%.2e',t));
    xlabel('r')
    ylabel('c')
    legend('analytical','numerical')
    ylim([0,1]);
    drawnow;
    end
    t = t+dt;
    it = it + 1;
end

figure(2)
plot(tvec,sqrt(err))
xlabel('t')
ylabel('absolute error')

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function B = ndim_expand(A, v)
%NDIM_EXPAND expand an array in a new dimension by multiplying it with a vector
%	B = NDIM_EXPAND(A, v)
%
%	A - multidimensional array (can be vector or matrix)
%	v - vector
%	B - add one new dimension to A as: [A.*v(1) | A.*v(2) | ... | A.*v(n)]
%
%	eg. ndim_expand(ones(2,3), [2 3 4])

s = size(A);
d = length(s);
if d==2 && s(2)==1
	v = v(:);
	if s(1)==1
		B = A*v;
	else
		B = A*v';
	end
else
	d = d + 1;
	n = length(v);
	B = v(1)*A;
	for i = 2:n
		B = cat(d, B, v(i)*A);
	end
end

author: Sable

repo: mcbench-benchmarks


text:
% DEMCLASSIFICATIONSVARGPLVM
% COPYRIGHT Andreas C. Damianou, 2012
% SEEALSO: demClassificationSvargplvm2, demClassificationSvargplvm3,  demClassificationSvargplvm4
% VARGPLVM

% clear;scale2var1 = 1;dataPerClass = 10;initial_X = 'together';itNo =
% [];indPoints = 30; initVardistIters = 300;latentDim = 6;demClassificationSvargplvm

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);
if ~exist('onlyTest'), onlyTest = false; end
if exist('diaryFile'),    diary(diaryFile), end
if ~exist('doPredictions'),    doPredictions = false; end
if ~exist('testDataPerClass'), testDataPerClass = 100; end
if ~exist('printPlots'), printPlots = true; end

addpath(genpath('../../bc-vargplvm/matlab'));
addpath(genpath('../../vargplvm/matlab'));

% This script initialises the options structure 'globalOpt'.
svargplvm_init;

if onlyTest
    load(['demOilSvargplvm' num2str(experimentNo)]);
    globalOpt = model.globalOpt;
end

globalOpt.dataSetNames = {'oilData', 'oilLabels'};
globalOpt.dataSetName = 'oil';

if exist('timeStamps')
    [globalOpt, YtrAll] = bc_loadData(globalOpt, timeStamps);
else
    [globalOpt, YtrAll] = bc_loadData(globalOpt);
end
Yall{1} = YtrAll.Y;
lblsTrain = YtrAll.lbls;
labelsTrain = transformLabels(lblsTrain);
% ---- We want labels with -1 and 1, not 0 and 1 (anti-correlation)
YtrAll.lbls(YtrAll.lbls == 0) = -1;
Yall{2} = YtrAll.lbls;

globalOptTest = globalOpt;
globalOptTest.dataPerClass = testDataPerClass;
globalOptTest.dataSetName = 'oilTest';
if exist('timeStamps')
    [globalOptTest, YtsAll] = bc_loadData(globalOptTest, timeStamps);
else
    [globalOptTest, YtsAll] = bc_loadData(globalOptTest);
end
lblsTest = YtsAll.lbls;
labelsTest = transformLabels(lblsTest);
Yts{1} = YtsAll.Y;
YtsAll.lbls(YtsAll.lbls == 0) = -1;
Yts{2} = YtsAll.lbls;

if isfield(globalOpt, 'normaliseData') && globalOpt.normaliseData
    Yall{1} = utils_normaliseData(Yall{1});
    Yts{1} = utils_normaliseData(Yts{1});
end

if onlyTest
    model = svargplvmRestorePrunedModel(model, Yall);
end

numberOfDatasets = length(Yall);

% globalOpt.baseKern = {'rbfardjit', 'linard2'};

%globalOpt.baseKern = {'rbfardjit', 'rbfardjit'};

globalOpt.indPoints = min(globalOpt.indPoints, size(Yall{1},1));

%-- Load datasets
for i=1:numberOfDatasets
    Y = Yall{i};
    dims{i} = size(Y,2);
    N{i} = size(Y,1);
    indTr = globalOpt.indTr;
    if indTr == -1
        indTr = 1:N{i};
    end
    if ~exist('Yts')
        indTs = setdiff(1:size(Y,1), indTr);
        Yts{i} = Y(indTs,:);
    end
    Ytr{i} = Y(indTr,:);

    t{i} = linspace(0, 2*pi, size(Y, 1)+1)'; t{i} = t{i}(1:end-1, 1);
    timeStampsTraining{i} = t{i}(indTr,1); %timeStampsTest = t(indTs,1);
end

for i=2:numberOfDatasets
    if N{i} ~= N{i-1}
        error('The number of observations in each dataset must be the same!');
    end
end

%%
if ~onlyTest
    % Free up some memory
    clear('Y')

    options = svargplvmOptions(Ytr, globalOpt, labelsTrain);

    if ~isempty(globalOpt.dynamicsConstrainType)
        for i=1:numberOfDatasets
            % Set up dynamcis (i.e. back-constraints) model
            optionsDyn{i}.type = 'vargpTime';
            optionsDyn{i}.inverseWidth=30;
            %   optionsDyn.vardistCovars = vardistCovarsMult;
            optionsDyn{i}.initX = globalOpt.initX;
            optionsDyn{i}.constrainType = globalOpt.dynamicsConstrainType;

            if exist('timeStampsTraining')
                optionsDyn{i}.t = timeStampsTraining;
            end
            if exist('labelsTrain') && ~isempty(labelsTrain)
                optionsDyn{i}.labels = labelsTrain;
            end
        end
    else
        optionsDyn= [];
    end

    model = svargplvmModelCreate(Ytr, globalOpt, options, optionsDyn);
    if exist('diaryFile')
        model.diaryFile = diaryFile;
    end

    model.globalOpt = globalOpt;
    model.options = options;

    %-- Define what level of parallelism to use (w.r.t submodels or/and w.r.t
    % datapoints).
    %{
fprintf('# Parallel computations w.r.t the submodels!\n');
model.parallel = 1;
model = svargplvmPropagateField(model,'parallel', 1);
%
fprintf('# Parallel computations w.r.t the datapoints!\n');
model.vardist.parallel = 1;
for i=1:model.numModels
    model.comp{i}.vardist.parallel = 1;
end
    %}

            %%%% TEMP
    if exist('whiteVar')
        fprintf('aaa\n\n\n')
        model.dynamics.kern.comp{2}.variance = whiteVar;
    end
    %%%%

    % Force kernel computations
    params = svargplvmExtractParam(model);
    model = svargplvmExpandParam(model, params);

    %%
    fprintf('# Median of vardist. covars: %d \n',median(median(model.vardist.covars)));
    fprintf('# Min of vardist. covars: %d \n',min(min(model.vardist.covars)));
    fprintf('# Max of vardist. covars: %d \n',max(max(model.vardist.covars)));

    model = svargplvmOptimiseModel(model);

    svargplvmShowScales(model)
    figure
end

%%
    capName = model.globalOpt.dataSetName;
    capName(1) = upper(capName(1));
    errors = fgplvmNearestNeighbour(model.comp{1}, lblsTrain);
    model2 = vargplvmReduceModel(model.comp{1}, 2);
    errors2 = fgplvmNearestNeighbour(model2, lblsTrain);
    fprintf('# Visualisation errors in all dims/in 2D:  %d / %d\n', errors, errors2)
    if printPlots
        vargplvmPrintPlot(model2, lblsTrain, [capName 'Vargplvm'], model.globalOpt.experimentNo);

        %%
        % For 3D plots:
         labels = labelsTrain; dims = [1 2 3];
         plot3k({model.X(:,dims(1)) model.X(:,dims(2)) model.X(:,dims(3))}, 'ColorData', labels, 'Marker', {'x',6});
    end

%%
obsMod = 1; % one of the involved sub-models (the one for which we have the data)
infMod = setdiff(1:2, obsMod);

sharedDims = svargplvmFindSharedDims(model);

%---------------------------- PREDICTIONS ---------------
if ~doPredictions
    return
end

%%

if ~exist('testOnTraining')
    testOnTraining=0;
end

%--------------------%
svargplvmPredictions %---- Script returning: ZpredMuAll and mini(the indices for NN)
%--------------------%

%--- For classification we are interested in labels
ZpredMuAll(ZpredMuAll > 0) = 1;
ZpredMuAll(ZpredMuAll < 0) = -1;
NNpred = Ytr{infMod}(mini,:);
realLabels = lblsTest(testInd,:);

% Find the error in the labels ------------
[labelErrors, labelSharingErrors, wrongIndices] = findLabelErrors(realLabels, ZpredMuAll);
[labelErrorsNN, void, wrongIndicesNN] = findLabelErrors(realLabels,NNpred);

fprintf('# Gplvm label errors: %d\n', labelErrors);
fprintf('# Gplvm label sharing errors: %d\n', labelSharingErrors);
fprintf('# NN label errors:    %d\n', labelErrorsNN);

% Confusion matrix
cmat_a = confMatrix( XpredAll, model.X, labelsTrain, labelsTest, 1);
cmat_rel = cmat_a ./ repmat(sum(cmat_a, 2), 1, length(unique(labelsTrain)));
cmat_rel

author: SheffieldML

repo: vargplvm


text:
function [S,Y,YS,lbfgs_start,lbfgs_end,Hdiag,skipped] = lbfgsAdd(y,s,S,Y,YS,lbfgs_start,lbfgs_end,Hdiag,useMex)
ys = y'*s;
skipped = 0;
corrections = size(S,2);
if ys > 1e-10
	if lbfgs_end < corrections
		lbfgs_end = lbfgs_end+1;
		if lbfgs_start ~= 1
			if lbfgs_start == corrections
				lbfgs_start = 1;
			else
				lbfgs_start = lbfgs_start+1;
			end
		end
	else
		lbfgs_start = min(2,corrections);
		lbfgs_end = 1;
	end

	if useMex
		lbfgsAddC(y,s,Y,S,ys,int64(lbfgs_end));
	else
		S(:,lbfgs_end) = s;
		Y(:,lbfgs_end) = y;
	end
	YS(lbfgs_end) = ys;

	% Update scale of initial Hessian approximation
	Hdiag = ys/(y'*y);
else
	skipped = 1;
end

repo: pTVreg


text:
% Quadratic equation function
%  determines coefficients of ax^2 + bx + c
%  input x,y arrays
%  by: Dr. Sherif Omran
%
%
function [a,b,c]=Quadratic(x,y)

p1=x(2)-x(3);
p2=x(3)-x(1);
p3=x(1)-x(2);
p4=x(3)^2-x(2)^2;
p5=x(1)^2-x(3)^2;
p6=x(2)^2-x(1)^2;
p7=x(2)^2*x(3)-x(2)*x(3)^2;
p8=x(1)*x(3)^2-x(1)^2*x(3);
p9=x(1)^2*x(2)-x(1)*x(2)^2;

delta=x(1)^2*(x(2)-x(3))-x(1)*(x(2)^2-x(3)^2)+1*(x(2)^2*x(3)-x(2)*x(3)^2);

a=(1/delta)*((x(2)-x(3))*y(1)+(x(3)-x(1))*y(2)+(x(1)-x(2))*y(3));
b=(1/delta)*(p4*y(1)+p5*y(2)+p6*y(3));
c=(1/delta)*(p7*y(1)+p8*y(2)+p9*y(3));

return;


text:
% plot the filtered estimates for navigation states, imu errors and imu
% error covariance stds.

kf = readdata(filresfile, 1+18);
if (~isOutNED)
%    save_google_kml(kf(:,2:4), kmlfilename);
end
if(useGPS)
inillh_ant2=options.inillh_ant;

posdata=loadAllGPSData(gpsfile, [kf(1,1), kf(end,1)]);
maxl=size(posdata,1);
if(isOutNED)
    for i=1:maxl
        rovXYZ = posdata(i,2:4);
        posdata(i,2:4)=posdiff_v001(rovXYZ',inillh_ant2);
    end
end
end
nextFig=1;
f(nextFig) = figure;
if (isOutNED)
    %         plot(kf(:,3),kf(:,2),'g.')
    %         hold on
    %         plot(posdata(:,3),posdata(:,2),'r+')
    %         grid
    %         xlabel('East [m]')
    %         ylabel('North[m]')
    plot3(kf(:,1)-kf(1,1),kf(:,3),kf(:,2),'g.')
    hold on
    if(useGPS)
    plot3(posdata(:,1)-kf(1,1),posdata(:,3),posdata(:,2),'r+')
    end
    grid
    axis equal
    xlabel('Time [s]')
    ylabel('East [m]')
    zlabel('North[m]')
else
    plot(kf(:,3),kf(:,2),'g.')
    hold on
    if(useGPS)
    plot(posdata(:,3),posdata(:,2),'r+')
    end
    grid
    axis equal
    xlabel('Longitude [radian]')
    ylabel('Lattitude [radian]')
end
title('green KF trajectory and the red GPS reference for GPS antenna');
saveas(f(nextFig),[resdir,'red truth and track'],'fig');

nextFig=nextFig+1;
f(nextFig) = figure;
plot(kf(:,1)-kf(1,1),kf(:,4),'g.')
hold on
if(useGPS)
plot(posdata(:,1)-kf(1,1),posdata(:,4),'r+');
end
grid
xlabel('Time [s]')
ylabel('Height/m')
title('Height of antenna by KF(green) and reference(red)');
saveas(f(nextFig),[resdir 'red truth and height'],'fig');

isPrintToFile=0;
if(options.useCam)
    campose=load(options.camPoseFile);
    if(~isempty(campose))
        nextFig=nextFig+1;
        f(nextFig) = figure;
        plot(campose(:,2)-campose(1,2),campose(:,6:8),'marker','.')
        grid
        xlabel('Time [s]')
        ylabel('m')
        legend('x','y','z')
        title('Ts2c');
        if isPrintToFile
            print(f(nextFig), '-dtiff', [plotDir 'Ts2c']);
        end

        nextFig=nextFig+1;
        f(nextFig) = figure;
        plot(campose(:,2)-campose(1,2),campose(:, 3:5),'marker','+');
        grid
        xlabel('Time [s]')
        ylabel('degree')
        legend('Roll','Pitch','Yaw/Heading')
        title('Cs2c');
        if isPrintToFile
            print(f(nextFig), '-dtiff', [plotDir 'qs2c']);
        end
        nextFig=nextFig+1;
        f(nextFig) = figure;
        plot(campose(:,2)-campose(1,2),campose(:,12:14),'marker','.')
        grid
        xlabel('Time [s]')
        ylabel('m')
        legend('x','y','z')
        title('Xs02e');
        if isPrintToFile
            print(f(nextFig), '-dtiff', [plotDir 'Xs02e']);
        end

        nextFig=nextFig+1;
        f(nextFig) = figure;
        plot(campose(:,2)-campose(1,2),campose(:, 9:11),'marker','+');
        grid
        xlabel('Time [s]')
        ylabel('degree')
        legend('Roll','Pitch','Yaw/Heading')
        title('qs02e');
        if isPrintToFile
            print(f(nextFig), '-dtiff', [plotDir 'qs02e']);
        end
    end
end
%close all
filresfile=[resdir, 'filresult.bin'];
kf = readdata(filresfile, 1+18);
plotNEDKalmanFilterResult(kf, resdir);
imuresfile=[resdir, 'imuresult.bin'];
err = readdata(imuresfile, 1+24);
ploterr_v001(err,resdir);

repo: ekfmonoslam


text:
% Compute the log intensity for the inverse link function g(f) = exp(-exp(-f)).
%
% The function is used in GLM likelihoods such as likPoisson, likGamma, likBeta
% and likInvGauss.
%
% Copyright (c) by Hannes Nickisch, 2013-10-16.

function [lg,dlg,d2lg,d3lg] = glm_invlink_expexp(f)
  lg = -exp(-f);
  if nargout>1
    dlg = -lg;
    if nargout>2
      d2lg = lg;
      if nargout>2
        d3lg = -lg;
      end
    end
  end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.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.

function CstTDMMeshSettings(mws, CellsPerWavelenth, MinimumCell)

%'@ set mesh properties (Hexahedral)
Mesh = invoke(mws,'Mesh');
invoke(Mesh,'MeshType','PBA');
invoke(Mesh,'SetCreator','High Frequency');

MeshSettings = invoke(mws,'MeshSettings');
invoke(MeshSettings,'SetMeshType','Hex');
invoke(MeshSettings,'Set','Version','1%');

%'MAX CELL - WAVELENGTH REFINEMENT
invoke(MeshSettings,'Set','StepsPerWaveNear',num2str(CellsPerWavelenth));
invoke(MeshSettings,'Set','StepsPerWaveFar',num2str(CellsPerWavelenth));
invoke(MeshSettings,'Set','WavelengthRefinementSameAsNear','1');
%'MAX CELL - GEOMETRY REFINEMENT
invoke(MeshSettings,'Set','StepsPerBoxNear',num2str(CellsPerWavelenth));
invoke(MeshSettings,'Set','StepsPerBoxFar',num2str(CellsPerWavelenth));
invoke(MeshSettings,'Set','MaxStepNear',num2str(CellsPerWavelenth));
invoke(MeshSettings,'Set','MaxStepFar',num2str(CellsPerWavelenth));
invoke(MeshSettings,'Set','ModelBoxDescrNear','maxedge');
invoke(MeshSettings,'Set','ModelBoxDescrFar','maxedge');
invoke(MeshSettings,'Set','UseMaxStepAbsolute','0');
invoke(MeshSettings,'Set','GeometryRefinementSameAsNear','1');
%'MIN CELL
invoke(MeshSettings,'Set','UseRatioLimitGeometry','1');
invoke(MeshSettings,'Set','RatioLimitGeometry',num2str(MinimumCell));
invoke(MeshSettings,'Set','MinStepGeometryX','0');
invoke(MeshSettings,'Set','MinStepGeometryY','0');
invoke(MeshSettings,'Set','MinStepGeometryZ','0');
invoke(MeshSettings,'Set','UseSameMinStepGeometryXYZ','1');

MeshSettings = invoke(mws,'MeshSettings');
invoke(MeshSettings,'SetMeshType','Hex');
invoke(MeshSettings,'Set','FaceRefinementOn','0');
invoke(MeshSettings,'Set','FaceRefinementPolicy','2');
invoke(MeshSettings,'Set','FaceRefinementRatio','2');
invoke(MeshSettings,'Set','FaceRefinementStep','0');
invoke(MeshSettings,'Set','FaceRefinementNSteps','2');
invoke(MeshSettings,'Set','EllipseRefinementOn','0');
invoke(MeshSettings,'Set','EllipseRefinementPolicy','2');
invoke(MeshSettings,'Set','EllipseRefinementRatio','2');
invoke(MeshSettings,'Set','EllipseRefinementStep','0');
invoke(MeshSettings,'Set','EllipseRefinementNSteps','2');
invoke(MeshSettings,'Set','FaceRefinementBufferLines','3');
invoke(MeshSettings,'Set','EdgeRefinementOn','1');
invoke(MeshSettings,'Set','EdgeRefinementPolicy','1');
invoke(MeshSettings,'Set','EdgeRefinementRatio','2');
invoke(MeshSettings,'Set','EdgeRefinementStep','0');
invoke(MeshSettings,'Set','EdgeRefinementBufferLines','3');
invoke(MeshSettings,'Set','RefineEdgeMaterialGlobal','0');
invoke(MeshSettings,'Set','RefineAxialEdgeGlobal','0');
invoke(MeshSettings,'Set','BufferLinesNear','3');
invoke(MeshSettings,'Set','UseDielectrics','1');
invoke(MeshSettings,'Set','EquilibrateOn','0');
invoke(MeshSettings,'Set','Equilibrate','1.5');
invoke(MeshSettings,'Set','IgnoreThinPanelMaterial','0');

MeshSettings = invoke(mws,'MeshSettings');
invoke(MeshSettings,'SetMeshType','Hex');
invoke(MeshSettings,'Set','SnapToAxialEdges','1');
invoke(MeshSettings,'Set','SnapToPlanes','1');
invoke(MeshSettings,'Set','SnapToSpheres','1');
invoke(MeshSettings,'Set','SnapToEllipses','1');
invoke(MeshSettings,'Set','SnapToCylinders','1');
invoke(MeshSettings,'Set','SnapToCylinderCenters','1');
invoke(MeshSettings,'Set','SnapToEllipseCenters','1');

Discretizer = invoke(mws,'Discretizer');
invoke(Discretizer,'MeshType','PBA');
invoke(Discretizer,'PBAType','Fast PBA');
invoke(Discretizer,'AutomaticPBAType','True');
invoke(Discretizer,'FPBAAccuracyEnhancement','enable');
invoke(Discretizer,'ConnectivityCheck','False');
invoke(Discretizer,'ConvertGeometryDataAfterMeshing','True');
invoke(Discretizer,'UsePecEdgeModel','True');
invoke(Discretizer,'GapDetection','False');
invoke(Discretizer,'FPBAGapTolerance','1e-3');
invoke(Discretizer,'SetMaxParallelMesherThreads','Hex','12');
invoke(Discretizer,'SetParallelMesherMode','Hex','Maximum');
invoke(Discretizer,'PointAccEnhancement','0');
invoke(Discretizer,'UseSplitComponents','True');
invoke(Discretizer,'EnableSubgridding','False');
invoke(Discretizer,'PBAFillLimit','99');
invoke(Discretizer,'AlwaysExcludePec','False');
end

repo: CST-MATLAB-API


text:
function[]=flipmap
%FLIPMAP Flips the current colormap upside-down.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2002, 2004 J.M. Lilly --- type 'help jlab_license' for details
ZZmap=colormap;
colormap(flipud(ZZmap));
clear ZZmap

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function RF = rfGaussian2d(X,Y,sigmaMajor,sigmaMinor,theta, x0,y0)
% rfGaussian2d - Create a two dimensional Gaussian receptive field
%
%  RF = rfGaussian2d(X,Y,sigmaMajor,sigmaMinor,theta,x0,y0);
%
%    X,Y        : Sample positions in deg
%    sigmaMajor : standard deviation longest direction
%    sigmaMinor : standard deviation shortest direction
%                 [default: sigmaMinor = sigmaMajor]
%    theta      : angle of sigmaMajor (radians, 0=vertical)
%                 [default = 0];
%    x0         : x-coordinate of center of RF [default = 0];
%    y0         : y-coordinate of center of RF [default = 0];
%
%
% Example:
% To make one rf:
%    fieldRange = 20;  % Deg
%    sampleRate = 0.2; % Deg
%    x = [-fieldRange:sampleRate:fieldRange];
%    y = x;
%    [X,Y] = meshgrid(x,y);
%    sigma = 5;  % Deg
%    rf = rfGaussian2d(X,Y,sigma);
%

% Programming notes:
%  Maybe we should always make sure we use an odd number of samples and
%  have a sample at the origin

% According to profile half of the executing time for the function
% is spend at notDefined, and this functions is called lots (it
% is the slowest part in rmMain together with rfMakePrediction),   so
% if it seems that all arguments are given just skip it.
% ras 08/06: well, we can do what Bob does: not call it at all anyway.
if nargin ~= 7,
    if ~exist('X', 'var') || isempty(X),
        error('Must define X grid');
    end;

    if ~exist('Y', 'var') || isempty(Y),
        error('Must define Y grid');
    end;

    if ~exist ('sigmaMajor', 'var') || isempty(sigmaMajor),
        error('Must scale on major axis');
    end;

    if ~exist ('sigmaMinor', 'var') || isempty(sigmaMinor),
        sigmaMinor = sigmaMajor;
    end;

    if ~exist ('theta', 'var') || isempty(theta), theta = false; end;
    if ~exist ('x0', 'var') || isempty(x0),       x0 = 0;    end;
    if ~exist ('y0', 'var') || isempty(y0),       y0 = 0;    end;
end;

% Allow sigma, x,y to be a matrix so that the final output will be
% size(X,1) by size(x0,2). This way we can make many RFs at the same time.
% Here I assume that all parameters are given.
if numel(sigmaMajor)~=1,
    sz1 = numel(X);
    sz2 = numel(sigmaMajor);

    X   = repmat(X(:),1,sz2);
    Y   = repmat(Y(:),1,sz2);

    sigmaMajor = repmat(sigmaMajor(:)',sz1,1);
    sigmaMinor = repmat(sigmaMinor(:)',sz1,1);

    if any(theta(:)),
        theta = repmat(theta(:)',sz1,1);
    end;

    x0 = repmat(x0(:)',sz1,1);
    y0 = repmat(y0(:)',sz1,1);
end;

% Translate grid so that center is at RF center
X = X - x0;   % positive x0 moves center right
Y = Y - y0;   % positive y0 moves center up

% Rotate grid around the RF center, positive theta rotates the
% grid to the right. No need for this if theta is 0.
if any(theta(:)),
    Xold = X;
    Yold = Y;
    X = Xold .* cos(theta) - Yold .* sin(theta);
    Y = Xold .* sin(theta) + Yold .* cos(theta);
end;

% make gaussian on current grid
RF = exp( -.5 * ((Y ./ sigmaMajor).^2 + (X ./ sigmaMinor).^2));

% Normalize the Gaussian.
% The idea is that if you stimulate the entire RF you will
% always get the same activation, independent of the RF parameters.
% RF = RF ./ (sigmaMajor.*2.*pi.*sigmaMinor);
% Decided not to do this. The fitting will deal with this.
% Thus amplitude will always be the same.

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ h ] = plot_q_mountain_car(Q)
%PLOT_Q_MOUNTAIN_CAR Summary of this function goes here
%   Detailed explanation goes here

% (N x 2)
s     = BuildStateList();
[v,I] = max(Q,[],2);

vs = -v;

xs      = linspace(-1.5,0.6,100);
ys      = linspace(-0.07,0.07,100);
[xi,yi] = meshgrid(xs,ys);

% Q (N x 2)

h = figure; hold on;
subplot(2,2,1);
zi = griddata(s(:,1),s(:,2),vs,xi,yi);
surf(xi,yi,zi); shading interp;
view(2);
title('Mountain car -max_u Q(x,u)');
xlabel('Position');
ylabel('Velocity');
colorbar;
ca = caxis;
xlim([-1.5,0.6]);
ylim([-0.07,0.07]);
axis square;
box on;

subplot(2,2,2);
plot_policy_moutain_car(Q);
xlim([-1.5,0.6]);
ylim([-0.07,0.07]);

subplot(2,2,3);
zi = griddata(s(:,1),s(:,2),-Q(:,1),xi,yi);
surf(xi,yi,zi); shading interp;
view(2);
title('-Q(x,-1)');
xlabel('Position');
ylabel('Velocity');
colorbar;
caxis(ca);
xlim([-1.5,0.6]);
ylim([-0.07,0.07]);
axis square;
box on;

subplot(2,2,4);
zi = griddata(s(:,1),s(:,2),-Q(:,2),xi,yi);
surf(xi,yi,zi); shading interp;
view(2);
title('-Q(x,1)');
xlabel('Position');
ylabel('Velocity');
colorbar;
caxis(ca);
xlim([-1.5,0.6]);
ylim([-0.07,0.07]);
axis square;
box on;

end

author: epfl-lasa


text:
function rel_els = find_rel_subjs(up, group_type)
%FIND_REL_SUBJS identifies the subjects for each group (YHVs, EHVs and
%OHVs)
%
%	Inputs:
%		db_data         demographic data extracted from the Vortal database
%       up              universal parameters structure
%
%	Outputs:
%       rel_els         a structure of subject els for each group
%

if strcmp(group_type, 'group')
    groups = up.paramSet.groups;
    group_names = unique(groups);
elseif strcmp(group_type, 'sub_group')
    groups = up.paramSet.sub_groups;
    group_names = unique(groups);
end

for group_no = 1 : length(group_names)
    eval(['rel_els.' group_names{group_no} ' = find(strcmp(groups, group_names{group_no} ));']);
end

%% Old Code
%
% loaded_info = load(up.paths.db_data);
%
% temp.group = loaded_info.db_data.group(loaded_info.rel_els);
%
% if length(unique(temp.group)) == 2
%     rel_els.yhv = find(temp.group == 1);
%     rel_els.ehv = find(temp.group == 2);
%     rel_els.global = 1:length(temp.group);
% elseif length(unique(temp.group)) == 3
%     rel_els.bw = find(temp.group == 1);
%     rel_els.am = find(temp.group == 2);
%     rel_els.fm = find(temp.group == 3);
%     rel_els.global = 1:length(temp.group);
% else
%     fprint('Unexpected Number of Groups');
% end

end

author: peterhcharlton

repo: RRest


text:
function imdb = setup_imdb_0131(datasetDir, varargin)

suffix = '.jpg';
allframes_subdir = 'allFrames';

imdb.imageDir = datasetDir;
imdb.meta.sets = {'train', 'val', 'test'};

seqs =dir(datasetDir);
seqs = {seqs([seqs.isdir]).name};
seqs = setdiff(seqs, {'.', '..'});

classes = dir(fullfile(datasetDir, seqs{1}));
classes = {classes([classes.isdir]).name};
imdb.meta.classes = setdiff(classes, {'.', '..', allframes_subdir});

imdb.images.name = {};
imdb.images.time = [];
imdb.images.class = [];
for i = 1:numel(seqs),
  seqName = seqs{i};
  imdb.images.time = [imdb.images.time ; parse_seq_name(seqName,30/(24*3600),0,-1)]; % TODO fps wrong at 1st frame
  files = dir(fullfile(datasetDir, seqName, allframes_subdir, ['*' suffix]));
  fileNames = sort({files.name});
	fileNames = fileNames(2:end); % TODO fps wrong at 1st frame
  class_tbl = -1*ones(1,numel(fileNames));
  for c = 1:numel(imdb.meta.classes),
    className = imdb.meta.classes{c};
    files = dir(fullfile(datasetDir, seqName, className, ['*' suffix]));
    class_tbl(ismember(fileNames,{files.name})) = c;
  end
  assert(~any(class_tbl==-1));
  imdb.images.class = [imdb.images.class class_tbl];
  imdb.images.name = [imdb.images.name ...
    cellfun(@(s) fullfile(seqName,allframes_subdir,s), fileNames,'UniformOutput',false)];
end

imdb.images.id = 1:numel(imdb.images.name);

end

function tstr = parse_seq_name(str, off_1, off_2, off_3)

if ~exist('off_1','var') || isempty(off_1), off_1 = 0; end;
if ~exist('off_2','var') || isempty(off_2), off_2 = 0; end;
if ~exist('off_3','var') || isempty(off_3), off_3 = 0; end;

fmt_in = 'yymmddHHMMSS';
fmt_out = 'yyyy-mm-dd HH:MM:SS.FFF';

sep_locs = strfind(str,'_');
assert(numel(sep_locs)==2, 'Wrong format: %s', str);

tnum_1 = datenum(str(1:sep_locs(1)-1), fmt_in) + off_1;
tnum_2 = datenum(str(sep_locs(1)+1:sep_locs(2)-1), fmt_in) + off_2;
N = str2num(str(sep_locs(2)+1:end)) + off_3;
tnums = linspace(tnum_1, tnum_2, N);

tstr = datestr(tnums, fmt_out);

end

repo: ResNet-Matconvnet


text:
function [fp, fm] = partition(f)
% PARTITION     Partition a spherefun into its even/periodic
%   odd/anti-periodic parts.
%
%   [fp, fm] = partition(f) partitions f into two spherefuns fp & fm with
%   the following properties:
%
%   fp has a CDR decomposition such that C is even and R is pi periodic
%   fm has a CDR decomposition such that C is odd and R is pi anti-periodic
%
% See also COMBINE

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

if ( ~isa(f, 'spherefun') )
    error('SPHEREFUN:partition:unknown', ['Undefined function '...
        'partition'' for input argument of type %s.'], class(f));
end

if ( isempty(f) )
    fp = spherefun();
    fm = spherefun();
    return
end

% Do the even-pi-periodic case first
id = f.idxPlus;
if ( isempty(id) )
    fp = spherefun();
else
    fp = f;
    fp.cols = fp.cols(:, id);
    fp.rows = fp.rows(:, id);
    fp.pivotValues = fp.pivotValues(id);
    fp.pivotLocations = fp.pivotLocations(id, :);
    fp.idxPlus = 1:length(id);
    fp.idxMinus = [];
end

% Now do the odd case
id = f.idxMinus;
if ( isempty(id) )
    fm = spherefun();
else
    fm = f;
    fm.cols = fm.cols(:, id);
    fm.rows = fm.rows(:, id);
    fm.pivotValues = fm.pivotValues(id);
    fm.pivotLocations = fm.pivotLocations(id, :);
    fm.idxMinus = 1:length(id);
    fm.idxPlus = [];
    fm.nonZeroPoles = 0;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function str = showdoc(item, indent)

% function str = showdoc(item, indent)
% Display help text for a cfg_menu and all of its options.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showdoc.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

str = showdoc(item.cfg_item, indent);
str{end+1} = sprintf('One of the following options must be selected:');
% Display cfg_menu labels
for k = 1:numel(item.labels)
    str{end+1} = sprintf('* %s', item.labels{k});
end;

author: spm


text:
function []=olsfeval(data_endo_c,stringdates3,Fstartdate,Fcenddate,Fcperiods,Fcomp,n,forecast_estimates,names,endo,pref)

% function []=olsfeval(data_endo_c,stringdates3,Fstartdate,Fcenddate,Fcperiods,Fcomp,n,forecast_estimates,names,endo,datapath)
% calculates, display and saves forecast evaluation results for a OLS VAR model
% inputs:  - matrix 'data_endo_c': matrix of endogenous data for the forecast evaluation period (i.e. period for which forecast is estimated and actual data exists)
%          - cell 'stringdates3': date strings for the forecast evaluation period (i.e. period for which forecast is estimated and actual data exists)
%          - string 'Fstartdate': start date of the forecasts
%          - string 'Fcenddate': end date of the forecat evaluation (i.e. period for which forecast is estimated and actual data exists)
%          - integer 'Fcperiods': number of periods for which forecast evaluation can be conducted (i.e.for which forecast is estimated and actual data exists)
%          - integer 'Fcperiods': number of periods for which forecast evaluation can be conducted (i.e.for which forecast is estimated and actual data exists)
%          - integer 'Fcomp': 0-1 value indicating if forecast evaluation is possible
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - cell 'forecast_estimates': lower bound, point estimates, and upper bound for the unconditional forecasts
%          - cell 'names': cell containing the excel spreadsheet labels (names and dates)
%          - cell 'endo': list of endogenous variables of the model
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% first, note that forecast evaluation can only be conducted if there is some observable data after the beginning of the forecast
if Fcomp==1

   % preliminary task: obtain a matrix of forecasts over the common periods
   for ii=1:n
   forecast_c(:,ii)=forecast_estimates{ii,1}(2,1:Fcperiods)';
   end
   % then compute the matrix of forecast errors
   ferrors=data_endo_c-forecast_c;

   % compute first the sequential RMSE, defined in (a.8.11)

   % square the forecast error matrix entrywise
   sferrors=ferrors.^2;
   % sum entries sequentially
   sumsferrors=sferrors(1,:);
   for ii=2:Fcperiods
   sumsferrors(ii,:)=sumsferrors(ii-1,:)+sferrors(ii,:);
   end
   % divide by the number of forecast periods and take square roots to obtain RMSE
   for ii=1:Fcperiods
   RMSE(ii,:)=((1/ii)*sumsferrors(ii,:)).^0.5;
   end

   % compute then the sequential MAE, defined in (a.8.12)

   % take the absolute value of the forecast error matrix
   absferrors=abs(ferrors);
   % sum entries sequentially
   sumabsferrors=absferrors(1,:);
   for ii=2:Fcperiods
   sumabsferrors(ii,:)=sumabsferrors(ii-1,:)+absferrors(ii,:);
   end
   % divide by the number of forecast periods to obtain MAE
   for ii=1:Fcperiods
   MAE(ii,:)=(1/ii)*sumabsferrors(ii,:);
   end

   % compute the sequential MAPE, defined in (a.8.13)

   % divide entrywise by actual values and take absolute values
   absratioferrors=abs(ferrors./data_endo_c);
   % sum entries sequentially
   sumabsratioferrors=absratioferrors(1,:);
   for ii=2:Fcperiods
   sumabsratioferrors(ii,:)=sumabsratioferrors(ii-1,:)+absratioferrors(ii,:);
   end
   % divide by 100*(number of forecast periods) to obtain MAPE
   for ii=1:Fcperiods
   MAPE(ii,:)=(100/ii)*sumabsratioferrors(ii,:);
   end

   % compute the Theil's inequality coefficient, defined in (a.8.14)

   % first compute the left term of the denominator
   % square entrywise the matrix of actual data
   sendo=data_endo_c.^2;
   % sum entries sequentially
   sumsendo=sendo(1,:);
   for ii=2:Fcperiods
   sumsendo(ii,:)=sumsendo(ii-1,:)+sendo(ii,:);
   end
   % divide by the number of forecast periods and take square roots
   for ii=1:Fcperiods
   leftterm(ii,:)=((1/ii)*sumsendo(ii,:)).^0.5;
   end
   % then compute the right term of the denominator
   % square entrywise the matrix of forecast values
   sforecasts=forecast_c.^2;
   % sum entries sequentially
   sumsforecasts=sforecasts(1,:);
   for ii=2:Fcperiods
   sumsforecasts(ii,:)=sumsforecasts(ii-1,:)+sforecasts(ii,:);
   end
   % divide by the number of forecast periods and take square roots
   for ii=1:Fcperiods
   rightterm(ii,:)=((1/ii)*sumsforecasts(ii,:)).^0.5;
   end
   % finally, compute the U stats
   Ustat=RMSE./(leftterm+rightterm);

% if forecast evaluation is not possible, do not do anything
elseif Fcomp==0
end

% now, print the results and display them

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'at');

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

Fevalinfo='Forecast evaluation:';
fprintf('%s\n',Fevalinfo);
fprintf(fid,'%s\n',Fevalinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% if forecast evaluation is not possible, return a message to signal it

if Fcomp==0

finfo1=['Forecast evaluation cannot be conducted.'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);
finfo2=['Forecasts start in ' Fstartdate ', while observable data is available only until ' names{end,1} '.'];
fprintf('%s\n',finfo2);
fprintf(fid,'%s\n',finfo2);
finfo3=['To obtain forecast evaluation, the forecast start date must be anterior to the end of the data set.'];
fprintf('%s\n',finfo3);
fprintf(fid,'%s\n',finfo3);

% if forecast evaluation is possible, display the results
elseif Fcomp==1

finfo1=['Evaluation conducted over ' num2str(Fcperiods) ' periods (from ' Fstartdate ' to ' Fcenddate ').'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);

   % loop over endogenous variables
   for ii=1:n

   fprintf('%s\n','');
   fprintf(fid,'%s\n','');

   endoinfo=['Endogenous: ' endo{ii,1}];
   fprintf('%s\n',endoinfo);
   fprintf(fid,'%s\n',endoinfo);

   temp='fprintf(''%12s';
      for jj=1:Fcperiods-1
      temp=[temp ' %10s'];
      end
   temp=[temp ' %10s\n'','''''];
      for jj=1:Fcperiods
      temp=[temp ',''' stringdates3{jj,1} ''''];
      end
   temp=[temp ');'];
   eval(temp);
   temp='fprintf(fid,''%12s';
      for jj=1:Fcperiods-1
      temp=[temp ' %10s'];
      end
   temp=[temp ' %10s\n'','''''];
      for jj=1:Fcperiods
      temp=[temp ',''' stringdates3{jj,1} ''''];
      end
   temp=[temp ');'];
   eval(temp);

   label='RMSE:       ';
   values=RMSE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='MAE:        ';
   values=MAE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='MAPE:       ';
   values=MAPE(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

   label='Theil''s U:  ';
   values=Ustat(1:Fcperiods,ii)';
   temp='fprintf(''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);
   temp='fprintf(fid,''%12s';
   for jj=1:Fcperiods-1
   temp=[temp ' %10.3f'];
   end
   temp=[temp ' %10.3f\n'''];
   temp=[temp ',label,values);'];
   eval(temp);

    end

end

fclose(fid);

author: european-central-bank

repo: BEAR-toolbox


text:
function f=loglikGaPExpWfixPriorH(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

[Hkt_r, cx, cy, Wxk] = loglikGaPreadparamWfix(x,varargin);
Hkt=exp(Hkt_r); %nonnegativity constrains

[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
P=Wxkbg*Hktbg; %current approximation

%Poisson contribution
t1=Vxt.*log(P) - P;
%Gamma contribution
t2=(alpha-1)*log(Hkt)-1/beta*Hkt-alpha*log(beta)-log(gamma(alpha));

f=sum(t1(:))+sum(t2(:));
f=-f; %conjugate gradient is mimimizing!
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function test_bug3124

% WALLTIME 00:30:00
% MEM 3gb
% DEPENDENCY ft_sourceanalysis

%%

load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/headmodel.mat'))
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer/data_all.mat'))

%%

cfg = [];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.output = 'fourier';
freq = ft_freqanalysis(cfg, data_all);

%%
% with 1 cm grid -> 3000 grid points, 1484 inside
% with 2 cm grid -> 186 grid points

cfg = [];
cfg.headmodel = ft_convert_units(headmodel, 'cm');
cfg.resolution = 2;
cfg.channel = 'MEG';

sourcemodel = ft_prepare_leadfield(cfg, freq);

%%
% make a manual selection of a single frequency

fcfg = [];
fcfg.frequency = 10;
freq1 = ft_selectdata(fcfg, freq);

fcfg = [];
fcfg.foilim = 10;
freq2 = ft_selectdata(fcfg, freq);

%%
% the result should be at 10 Hz in all cases

cfg = [];
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
cfg.frequency = 10;
source0 = ft_sourceanalysis(cfg, freq);

cfg = [];
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
source1 = ft_sourceanalysis(cfg, freq1);

cfg = [];
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
source2 = ft_sourceanalysis(cfg, freq2);

assert(isequal(source0.freq, [10]));
assert(isequal(source1.freq, [10]));
assert(isequal(source2.freq, [10]));

assert(isequaln(source0.avg.pow, source1.avg.pow));
assert(isequaln(source0.avg.pow, source2.avg.pow));

%%
% make a manual selection of a range

fcfg = [];
fcfg.frequency = [9 11];
freq1 = ft_selectdata(fcfg, freq);

fcfg = [];
fcfg.foilim = [9 11];
freq2 = ft_selectdata(fcfg, freq);

%%
% the result should be at 10 Hz in all cases, i.e. averaged from 9-11

cfg = [];
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
cfg.frequency = [9 11];
source0 = ft_sourceanalysis(cfg, freq);

cfg = [];
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
source1 = ft_sourceanalysis(cfg, freq1);

cfg = [];
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
source2 = ft_sourceanalysis(cfg, freq2);

assert(isequal(source0.freq, [10]));
assert(isequal(source1.freq, [10]));
assert(isequal(source2.freq, [10]));

assert(isequaln(source0.avg.pow, source1.avg.pow));
assert(isequaln(source0.avg.pow, source2.avg.pow));

author: fieldtrip

repo: fieldtrip


text:
function fx = edcheb ( x, coef, npl )

%*****************************************************************************80
%
%% EDCHEB evaluates the derivative of a Chebyshev series at a point.
%
%  Discussion:
%
%    This routine evaluates the derivative of a Chebyshev series
%    at a point in [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Roger Broucke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Broucke,
%    Algorithm 446:
%    Ten Subroutines for the Manipulation of Chebyshev Series,
%    Communications of the ACM,
%    October 1973, Volume 16, Number 4, pages 254-256.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) X, the evaluation point.
%    -1 <= X <= +1.
%
%    Input, real ( kind = 8 ) COEF(NPL), the Chebyshev series.
%
%    Input, integer ( kind = 4 ) NPL, the number of terms in the
%    Chebyshev series.
%
  xjp2 = 0.0;
  xjpl = 0.0;
  bjp2 = 0.0;
  bjpl = 0.0;
  n = npl - 1;

  for k = 1 : n
    j = npl - k;
    dj = j;
    xj = 2.0 * coef(j+1) * dj + xjp2;
    bj = 2.0 * x * bjpl - bjp2 + xj;
    bf = bjp2;
    bjp2 = bjpl;
    bjpl = bj;
    xjp2 = xjpl;
    xjpl = xj;
  end

  fx = 0.5 * ( bj - bf );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function P = compute_image_pyramid2(img, sig, nL, ratio)
%%  COMPUTE_IMAGE_PYRAMID computes nL level image pyramid of the input image IMG using filter F
% downsample each level directly from the original images
%
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-10 $
%   $Revision $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

P   = cell(nL,1);
tmp = img;
P{1}= tmp;

% Get version information (from http://www.mathworks.com/matlabcentral/fileexchange/17285)
v = sscanf (version, '%d.%d.%d') ; v = 10.^(0:-1:-(length(v)-1)) * v ;

for m = 2:nL

    % Gaussian filtering
    s   = sig / ratio^(m-2);
    f   = fspecial('gaussian', 2*round(1.5*s) +1, s);
    tmp = imfilter(img, f, 'corr', 'symmetric', 'same');

    sz  = round([size(img,1) size(img,2)]*ratio^(m-1));

    %tmp = imresize(img, sz, 'bilinear', 'Antialiasing', true);   % better than 'nearest'
    % IMRESIZE changes default algorithm since version 7.4 (R2007a)
    if v > 7.3
        tmp = imresize(tmp, sz, 'bilinear', 'Antialiasing', false);   % better than 'nearest'
    else
        tmp = imresize(tmp, sz, 'bilinear', 0); % Disable antialiasing, old version for cluster
    end;

    P{m} = tmp;
end;

author: kristinbranson

repo: JAABA


text:
function f = f10_f0 ( n, x, y )

%*****************************************************************************80
%
%% F10_F0 returns the value of function 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real F(N,1), the function values.
%
  t1(1:n,1) = sqrt ( ( 80.0 * x(1:n,1) - 40.0 ).^2 + ( 90.0 * y(1:n,1) - 45.0 ).^2 );
  t2(1:n,1) = exp ( - 0.04 * t1(1:n,1) );
  t3(1:n,1) = cos ( 0.15 * t1(1:n,1) );

  f(1:n,1) = t2(1:n,1) .* t3(1:n,1);

  return
end


text:
function [rgb] = bg_rgba2rgb(bg, rgba, varargin)

% BG_RGBA2RGB overlays a transparency masked colored image on a colored background,
% and represents the result as an RGB matrix.
%
% Use as:
%   rgb = bg_rgba2rgb(bg, rgba)
%
% or
%   rgb = bg_rgba2rgb(bg, rgba, cmap, clim, alpha, amap, alim);
%
% When 2 input arguments are supplied:
%   bg   = Nx3 matrix of background rgb-coded color-values, or MxNx3
%   rgba = Nx4 matrix of rgb + alpha values, or MxNx4
%
% When 7 input arguments are supplied:
%   bg   = Nx3 matrix, Nx1 vector, 1x3 vector, MxN, or MxNx3.
%   rgba = Nx1 vector with 'functional values', or MxN.
%   cmap = Mx3 colormap, or MATLAB-supported name of colormap
%   clim = 1x2 vector denoting the color limits
%   alpha = Nx1 vector with 'alpha values', or MxN
%   amap = Mx1 alphamap, or MATLAB -supported name of alphamap ('rampup/down', 'vup/down')
%   alim = 1x2 vector denoting the opacity limits

if numel(varargin)==0
  siz1 = size(bg);
  siz2 = size(rgba);
  assert(isequal(siz1(1:end-1),siz2(1:end-1)), 'number of data points should be the same');
  assert(siz1(end)==3 && siz2(end)==4,         'inconsistent input size');

  bg   = reshape(bg,   [], 3);
  rgba = reshape(rgba, [], 4);
  rgb  = do_conversion(bg, rgba);
  rgb  = reshape(rgb, siz1);
else
  % this requires the data to be converted into rgb values irst, and
  % needs a cmap + clim, and an alpha, alphamap and alim
  if numel(varargin)~=5
    error('if a vector of color data is supplied, more input arguments are required');
  end

  % colormap handling
  cmap = varargin{1};
  if ischar(cmap)
    cmap = ft_colormap(cmap);
  elseif size(cmap,2)~=3
    error('invalid specification of colormap, should be nx3');
  end
  clim  = varargin{2}; if isempty(clim), clim(1) = min(rgba(:)); clim(2) = max(rgba(:)); end
  alpha = varargin{3};
  amap  = varargin{4};
  alim  = varargin{5}; if isempty(alim), alim(1) = min(alpha(:)); alim(2) = max(alpha(:)); end

  % deal with the color data
  dat     = rgba;
  siz     = size(dat);
  dat     = reshape(dat, [], 1);
  alpha   = reshape(alpha, [], 1); % assume to be same siza as input data!
  finvals = isfinite(dat);

  rgba              = zeros(size(dat,1),4);
  rgba(finvals,1:3) = dat2rgb(dat(finvals), cmap, clim);

  finvals          = isfinite(alpha);
  rgba(finvals,4)  = alpha2a(alpha(finvals), amap, alim);
  rgba(~finvals,4) = 0;

  % deal with the background: allow it to be 1x3, i.e. a single color per
  % pixel.
  if isequal(size(bg),[1 3])
    bg = permute(repmat(bg(:), [1 siz]), [1+(1:numel(siz)) 1]);
  else
    siz_bg = size(bg);
    if isequal(siz_bg, siz)
      % make bg rgb
      if all(bg(:)<=1) && all(bg(:)>=0)
        % don't scale
      else
        bg_min = min(bg(:));
        bg_max = max(bg(:));
        bg     = (bg-bg_min)./(bg_max-bg_min);
      end
      bg = repmat(bg, [ones(1,ndims(bg)) 3]);
    else
      %FIXME
    end
  end
  bg  = reshape(bg, [], 3);

  if numel(siz)==2 && siz(2)==1, siz = siz(1); end
  rgb = do_conversion(bg, rgba);
  rgb = reshape(rgb, [siz 3]);
end

function rgb = do_conversion(bg, rgba)

rgb = zeros(size(rgba,1),3);
a_  = 1-rgba(:,4);
a   = rgba(:,4);

for k = 1:3
    rgb(:,k) = bg(:,k).*a_  + rgba(:,k).*a;
end

function rgb = dat2rgb(dat, cmap, clim)

dat(end+1) = clim(1);
dat(end+1) = clim(2); % add the extremes to be sure that they are included

dat(dat<clim(1)) = clim(1);
dat(dat>clim(2)) = clim(2);

% scale between 0 and 1
dat = dat-min(dat);
dat = dat/max(dat);

ind = round(dat.*(size(cmap,1)-1))+1;
rgb = cmap(ind(1:end-2),:);

function a = alpha2a(alpha, amap, alim)

alpha(end+1) = alim(1);
alpha(end+1) = alim(2);

alpha(alpha<alim(1)) = alim(1);
alpha(alpha>alim(2)) = alim(2);

if ischar(amap)
    switch amap
        case 'rampup'
            a = alpha - min(alpha);
            a = a./max(a);
        case 'rampdown'
            a = alpha - min(alpha);
            a = a./max(a);
            a = 1-a;
        case 'vup'
            a = alpha - min(alpha);
            a = a./max(a);
            a = 1 - 2.*abs(a - 0.5);
        case 'vdown'
            a = alpha - min(alpha);
            a = a./max(a);
            a = 2.*abs(a - 0.5);
        otherwise
            error('unknown alphamap specified');
    end
else
    amap = amap(:);

    % assume it to be a vector that has M elements, scaled between 0 and 1
    a = alpha-min(alpha);
    a = a/max(a);
    ind = round(a.*(size(amap,1)-1))+1;
    a   = amap(ind(1:end),:);
end
a = a(1:end-2);


function dataout=modeDRAEC_BSS(nummics, numrefs, datain)
%
% Perform dr and aec together, then bss.
% nummics:              no. of mic channels
% numrefs:              no. of reference channels
% datain:               input data
% dataout:              output data
%

%% perform stft
addpath('stft2');
% fft size
fftsize=512;
stftshift=fftsize/2;

M=nummics;
R=numrefs;
N=M;

Xtf=cell(M+R, 1);
for m=1:M+R
    Xtf{m}=stft(datain(:, m), stftshift, fftsize, false);
end
[K, T]=size(Xtf{1});

Ytf=cell(M, 1);
for m=1:M
    Ytf{m}=zeros(K, T);
end

%% params
% dr filter length
DR_FLEN=5;
% forgetting factor for dr and aec
DRAEC_FORGET=0.999;
% forgetting factor for bss
BF_FORGET=0.999;
% the shape parameter of the source prior
GAMMA=0.2;
%
% used to keep stable
%
VAR_BIAS=0.01;
STABLE_EPS=1e-3;
DRAEC_DIAGLOAD=1e-6;
BF_DIAGLOAD=1e-6;

%% space for dr and aec
% current mic data backup
Miccurrent=zeros(K, M);
% ref and delayed mic data
draecfsize=R+M*DR_FLEN;
Refmicdelay=zeros(K, draecfsize);

% mic-ref correlation
Cxr=cell(K, 1);
for k=1:K
    Cxr{k}=zeros(M, draecfsize);
end

% reference auto correlation
Crr=cell(K, 1);
for k=1:K
    Crr{k}=zeros(draecfsize, draecfsize);
end

% reverb and echo path
REPath=cell(K, 1);
for k=1:K
    REPath{k}=zeros(M, draecfsize);
end

%% space for bss
% the weighted correlation matrices
C1=cell(K, 1);
C2=cell(K, 1);
for k=1:K
    C1{k}=STABLE_EPS*eye(M, M);
    C2{k}=STABLE_EPS*eye(M, M);
end

% demixing matrices
Demix=cell(K, 1);
for k=1:K
    Demix{k}=eye(N, M);
end

%% perform iteration
for tau=1:T
    %% perform dr and aec together
    % direct nearend and early reverberation
    Early=zeros(K, M);

    %
    % shift in new data
    %
    % shift old mic data back
    Refmicdelay=circshift(Refmicdelay, M, 2);
    % shift in delayed mic data
    Refmicdelay(:, R+1:R+M)=Miccurrent;

    % mic data backup
    for m=1:M
        Miccurrent(:, m)=Xtf{m}(:, tau);
    end

    % shift in reference data
    for r=1:R
        Refmicdelay(:, r)=Xtf{M+r}(:, tau);
    end

    for k=1:K
        % calculate late reverberation and echo
        ref=Refmicdelay(k, :).';
        late=REPath{k}*ref;

        % direct nearend and early reverberation
        mic=Miccurrent(k, :).';
        early=mic-late;
        % output data
        Early(k, :)=early.';

        %
        % calculate nonlinearity
        %
        xsq=abs(mic).^2;
        ysq=abs(early).^2;

        phi=0;
        for m=1:M
            if ysq(m)<=xsq(m)
                phi=phi+ysq(m);
            else
                phi=phi+xsq(m);
            end
        end

        phi=(1-DRAEC_FORGET)*(phi+VAR_BIAS)^((GAMMA-2)/2);

        % update mic ref correlation
        Cxr{k}=DRAEC_FORGET*Cxr{k}+phi*(mic*ref');

        % update ref auto-correlation
        Crr{k}=DRAEC_FORGET*Crr{k}+phi*(ref*ref');

        % update echo and reverb path
        REPath{k}=Cxr{k}/(Crr{k}+DRAEC_DIAGLOAD*eye(draecfsize, draecfsize));
    end

    %% perform bss
    Bssout=zeros(K, M);

    %
    % calculate nonlinearity
    %
    phi1=0;
    phi2=0;

    for k=1:K
        x=Early(k, :).';
        y=Demix{k}*x;
        % output data
        Bssout(k, :)=y.';

        phi1=phi1+abs(y(1))^2;
        phi2=phi2+abs(y(2))^2;
    end

    phi1=(1-BF_FORGET)*(phi1+VAR_BIAS)^((GAMMA-2)/2);
    phi2=(1-BF_FORGET)*(phi2+VAR_BIAS)^((GAMMA-2)/2);

    % update the demixing matrices
    for k=1:K
        %
        % accumulate the weighted correlation
        %
        x=Early(k, :).';
        C1{k}=BF_FORGET*C1{k}+phi1*(x*x');
        C2{k}=BF_FORGET*C2{k}+phi2*(x*x');

        % solve gev problem
        D=heig2(BF_DIAGLOAD, C2{k}, C1{k});
        Demix{k}=D;
    end

    for m=1:M
        Ytf{m}(:, tau)=Bssout(:, m);
    end
end

%% perform istft and output signal
dataout=zeros(dataLength(T, stftshift, fftsize ), N);
for n=1:N
    dataout(:, n)=istft(Ytf{n}, stftshift, false);
end

end


text:
function zz=lpcfq2zz(f,q)
%LPCFQ2ZZ Convert frequencies and q factors to z-plane poles ZZ=(F,Q)
%all input values are in normalized Hz
% roots are at exp(2*pi*f*(-1/(2q) +- j)
% if f has more columns than q, remaining columns are real roots at -f

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcfq2zz.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,pf]=size(f);
if nargin < 2
   pq=0;
else
   pq=size(q,2);
end;
zz=zeros(nf,pf+pq);
if pq
   ii=1:pq;
   zz(:,2*ii-1)=exp(pi*f(:,ii).*(2i-q.^(-1)));
   zz(:,2*ii)=conj(zz(:,2*ii-1));
end
if pf>pq
   ii=1+pq:pf;
   zz(:,ii+pq)= exp(-2*pi*f(:,ii));
end

author: decouples


text:
function [ fea, out ] = ex_planestress2( varargin )
%EX_PLANESTRESS2 NAFEMS benchmark challenge 1 plane stress example.
%
%   [ FEA, OUT ] = EX_PLANESTRESS2( VARARGIN ) NAFEMS benchmark example to calculate von Mieses stress
%   thin plate under plane stress assumption with loads all around.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       E           scalar {210e9}         Modulus of elasticity
%       nu          scalar {0.3}           Poissons ratio
%       igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
%       hmax        scalar {1/20}          Max grid cell size
%       sfun        string {sflag1}        Shape function for displacements
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'E',        210e9; ...
  'nu',       0.3; ...
  'igrid',    0; ...
  'hmax',     1/20; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
gobj1 = gobj_rectangle( 0, 1, 0, 1, 'R1' );
fea.geom.objects = { gobj1 };
fea.sdim = { 'x' 'y' };

% Grid generation.
if( opt.igrid==1 )
  fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
else
  nx = 1/opt.hmax;
  fea.grid = rectgrid( nx, nx );
end

% Boundary conditions.
dtol = 0.1;
if( opt.igrid==1 )
  lbdr = findbdr( fea, ['x<',num2str(dtol)] );     % Left boundary number.
  rbdr = findbdr( fea, ['x>',num2str(1-dtol)] );   % Right boundary number.
  tbdr = findbdr( fea, ['y>',num2str(1-dtol)] );   % Top boundary number.
  bbdr = findbdr( fea, ['y<',num2str(dtol)] );     % Bottom boundary number.
else
  lbdr = 4;
  rbdr = 2;
  tbdr = 3;
  bbdr = 1;
end

% Add plane stress physics mode.
fea = addphys(fea,@planestress);
fea.phys.pss.eqn.coef{1,end} = { opt.nu };
fea.phys.pss.eqn.coef{2,end} = { opt.E  };
fea.phys.pss.sfun            = { opt.sfun opt.sfun };

% Set boundary condition types.
bctype = mat2cell( zeros(2,4), [1 1], [1 1 1 1] );
fea.phys.pss.bdr.coef{1,5}   = bctype;

% Set loads on boundary.
bccoef = mat2cell( zeros(2,4), [1 1], [1 1 1 1] );
bccoef{1,tbdr} = '-x';
bccoef{2,tbdr} = 'x';
bccoef{1,bbdr} = '-(1-x)';
bccoef{2,bbdr} = '1-x';
bccoef{1,lbdr} = '1-y';
bccoef{2,lbdr} = '-(1-y)';
bccoef{1,rbdr} = 'y';
bccoef{2,rbdr} = '-y';
fea.phys.pss.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea       = parsephys(fea);             % Check and parse physics modes.
fea       = parseprob(fea);             % Check and parse problem struct.
fea.sol.u = solvestat(fea,'fid',fid);   % Call to stationary solver.

% Postprocessing.
s_vm = fea.phys.pss.eqn.vars{1,2};
if ( opt.iplot>0 )
  figure
  postplot( fea, 'surfexpr', s_vm, 'isoexpr', s_vm )
  title('von Mieses stress')
end

% Error checking.
dtol = sqrt(eps);
sm  = evalexpr( s_vm, [0.5;0.5], fea );
s01 = evalexpr( s_vm, [dtol;1],    fea );
s10 = evalexpr( s_vm, [1;dtol],    fea );
s00 = evalexpr( s_vm, [dtol;dtol], fea );
s11 = evalexpr( s_vm, [1;1],       fea );

sol = [ sm s01 s10 s00 s11 ];
ref = [  0   0   0   2   2 ];
out.err    = norm(sol-ref)/norm(ref);
out.pass   = out.err < 0.1;

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function triangle_dunavant_rule_test02 ( )

%*****************************************************************************80
%
%% TRIANGLE_DUNAVANT_RULE_TEST02 tests DUNAVANT_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_DUNAVANT_RULE_TEST02\n' );
  fprintf ( 1, '  DUNAVANT_RULE returns the points and weights\n' );
  fprintf ( 1, '  of a Dunavant rule for the triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we simply check that the weights\n' );
  fprintf ( 1, '  sum to 1.\n' );

  rule_num = dunavant_rule_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of available rules = %d\n', rule_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '      Rule    Order    Sum of weights\n' );
  fprintf ( 1, '\n' );

  for rule = 1 : rule_num

    order_num = dunavant_order_num ( rule );

    [ xy, w ] = dunavant_rule ( rule, order_num );

    w_sum = sum ( w(1:order_num) );

    fprintf ( 1, '  %8d  %8d  %14f\n', rule, order_num, w_sum );

  end

  return
end


text:
function varargout = pm_get_defaults(defstr, varargin)

% FORMAT defval = pm_get_defaults(defstr)
% Return the defaults value associated with identifier "defstr".
% Currently, this is a '.' subscript reference into the global
% "defaults" variable defined in spm_defaults.m.
%
% FORMAT spm_get_defaults(defstr, defval)
% Sets the defaults value associated with identifier "defstr". The new
% defaults value applies immediately to:
% * new modules in batch jobs
% * modules in batch jobs that have not been saved yet
% This value will not be saved for future sessions of SPM. To make
% persistent changes, edit spm_defaults.m.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton
% $Id: pm_get_defaults.m 4572 2011-11-23 17:35:10Z chloe $
% ---------------------------------------------------------------------
global pm_def
if isempty(pm_def)
    pm_defaults;
    pm_def.sessname='session';
    pm_def.pedir=2;
elseif ~isfield(pm_def,'pedir');
    pm_def.pedir=2;
elseif ~isfield(pm_def,'sessname');
    pm_def.sessname='session';
end

[default_file_path, tmpname] = fileparts(mfilename('fullpath'));
pm_def.defaultsfilename{1} = sprintf('%s%s%s',default_file_path,filesep,'pm_defaults.m');
if ischar(pm_def.MFLAGS.TEMPLATE)
    pm_def.MFLAGS.TEMPLATE = {pm_def.MFLAGS.TEMPLATE};
end
% tmpTEMPLATE = pm_def.MFLAGS.TEMPLATE;
% pm_def.MFLAGS.TEMPLATE=[];
% pm_def.MFLAGS.TEMPLATE{1}=tmpTEMPLATE;

% construct subscript reference struct from dot delimited tag string
tags = textscan(defstr,'%s', 'delimiter','.');
subs = struct('type','.','subs',tags{1}');

if nargin == 1
    varargout{1} = subsref(pm_def, subs);
else
    pm_def = subsasgn(pm_def, subs, varargin{1});
end


text:
%REFINEDETECTEDMARKERS  Refind not detected markers based on the already detected and the board layout
%
%     [detectedCorners, detectedIds, rejectedCorners] = cv.refineDetectedMarkers(img, board, detectedCorners, detectedIds, rejectedCorners)
%     [detectedCorners, detectedIds, rejectedCorners, recoveredIdxs] = cv.refineDetectedMarkers(...)
%     [...] = cv.refineDetectedMarkers(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __img__ input image (8-bit grayscale or color).
% * __board__ layout of markers in the board.
% * __detectedCorners__ cell array of already detected marker corners
%   `{{[x,y],..}, ..}`.
% * __detectedIds__ vector of already detected marker identifiers (0-based).
% * __rejectedCorners__ cell array of rejected candidates during the marker
%   detection process `{{[x,y],..}, ..}`.
%
% ## Output
% * __detectedCorners__ output refined marker corners.
% * __detectedIds__ output refined marker identifiers
% * __rejectedCorners__ output refined rejected corners.
% * __recoveredIdxs__ Optional array that returns the indexes of the recovered
%   candidates in the original `rejectedCorners` array.
%
% ## Options
% * __CameraMatrix__ optional input 3x3 floating-point camera matrix
%   `A = [fx 0 cx; 0 fy cy; 0 0 1]`.
% * __DistCoeffs__ optional vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4]` of 4, 5, 8 or 12 elements.
% * __MinRepDistance__ minimum distance between the corners of the rejected
%   candidate and the reprojected marker in order to consider it as a
%   correspondence. default 10.0
% * __ErrorCorrectionRate__ rate of allowed erroneous bits respect to the
%   error correction capability of the used dictionary. -1 ignores the error
%   correction step. default 3.0
% * __CheckAllOrders__ Consider the four posible corner orders in the
%   `rejectedCorners` array. If it set to false, only the provided corner
%   order is considered (default true).
% * __DetectorParameters__ marker detection parameters.
%
% This function tries to find markers that were not detected in the basic
% cv.detectMarkers function. First, based on the current detected marker and
% the board layout, the function interpolates the position of the missing
% markers. Then it tries to find correspondence between the reprojected
% markers and the rejected candidates based on the `minRepDistance` and
% `errorCorrectionRate` parameters. If camera parameters and distortion
% coefficients are provided, missing markers are reprojected using
% cv.projectPoints function. If not, missing marker projections are
% interpolated using global homography, and all the marker corners in the
% board must have the same Z coordinate.
%
% See also: cv.detectMarkers, cv.estimatePoseBoard, cv.findHomography,
%  cv.perspectiveTransform
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%FEATSELI Trainable mapping for individual feature selection
%
%   [W,R] = FEATSELI(A,CRIT,K,T)
%   [W,R] = A*FEATSELI([],CRIT,K,T)
%   [W,R] = A*FEATSELI(CRIT,K,T)
%   [W,R] = FEATSELI(A,CRIT,K,N)
%   [W,R] = A*FEATSELI([],CRIT,K,N)
%   [W,R] = A*FEATSELI(CRIT,K,N)
%
% INPUT
%   A    Training dataset
%   CRIT Name of the criterion or untrained mapping
%        (default: 'NN', i.e. the LOO 1-Nearest Neighbor error)
%   K    Number of features to select (default: sort all features)
%   T    Tuning dataset (optional)
%   N    Number of cross-validations (optional)
%
% OUTPUT
%   W    Feature selection mapping
%   R    Matrix with criterion values
%
% DESCRIPTION
% Individual selection of K features using the dataset A. CRIT sets the
% criterion used by the feature evaluation routine FEATEVAL. If the dataset
% T is given, it is used as test set for FEATEVAL. For K = 0 all features are
% selected, but reordered according to the criterion. The result W can be
% used for selecting features using B*W.
% The selected features are stored in W.DATA and can be found by +W.
% In R, the search is reported step by step as:
%
% 	R(:,1) : number of features
% 	R(:,2) : criterion value
% 	R(:,3) : added / deleted feature
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, FEATEVAL, FEATSELO, FEATSELB, FEATSELF,
% FEATSEL, FEATSELP,FEATSELM

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: featseli.m,v 1.5 2009/07/01 07:48:05 davidt Exp $

function [w,r] = featseli(varargin)

  varargin = shiftargin(varargin,{'char','prmapping'});
  argin = setdefaults(varargin,[],'NN',0,[]);
  if mapping_task(argin,'definition')
    w = define_mapping(argin,'untrained','Individual FeatSel');
    return
  end

  [a,crit,ksel,t] = deal(argin{:});

	[m,k,c] = getsize(a); featlist = getfeatlab(a);

	% If KSEL is not given, return all features.

	if (ksel == 0), ksel = k; end

	isvaldfile(a,1,2); % at least 1 object per class, 2 classes
	a = testdatasize(a);
	if isdataset(t), iscomdset(a,t); end

	critval = zeros(k,1);

	% Evaluate each feature in turn.

	s = sprintf('Evaluation of %i features: ',k);
	prwaitbar(k,s);
  a = setprior(a,getprior(a));
	if (isempty(t))
		for j = 1:k
			prwaitbar(k,j,[s int2str(j)]);
			critval(j) = feateval(a(:,j),crit);
		end
  elseif isdataset(t)
		for j = 1:k
			prwaitbar(k,j,[s int2str(j)]);
			critval(j) = feateval(a(:,j),crit,t(:,j));
    end
  else
		for j = 1:k
			prwaitbar(k,j,[s int2str(j)]);
			critval(j) = feateval(a(:,j),crit,t);
    end
	end
	prwaitbar(0);

	% Sort the features by criterion value (maximum first).
	[critval_sorted,J] = sort(-critval);
	r = [[1:k]', -critval_sorted, J];
	J = J(1:ksel)';

	% Return the mapping found.
	w = featsel(k,J);
  w = setmapping_type(w,'trained');
  w = setsize(w,[k length(J)]);
	if ~isempty(featlist)
		w = setlabels(w,featlist(J,:));
	end
	w = setname(w,'Individual FeatSel');

return

author: marianux

repo: ecg-kit


text:
function m = shiftleft(m,r)
%SHIFTLEFT    Shift array m left by r bits
%
%   m = shiftleft(m,r)
%
%first r bits of m(:,1) are zero, otherwise lost
%

% written  12/30/98     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_LONG_LOGBETA = getappdata(0,'INTLAB_LONG_LOGBETA');

  Ones = ones(1,size(m,2));

  % first r bits of m
  factor = 2.^(r-INTLAB_LONG_LOGBETA)*Ones ;
  ms = floor( m .* factor );

  % array m shifted left by r bits, last r bits zero
  m =  ( m - ms ./ factor ) .* ( 2.^r*Ones );

  % add last r bits
  m(:,1:end-1) = m(:,1:end-1) + ms(:,2:end);

repo: OpenMAS


text:
function [ Accum_Grad ] = average_gradients_in_frames( Frames )
%AVERAGE_GRADIENTS_IN_FRAMES Summary of this function goes here
%   Detailed explanation goes here
n_frames=length(Frames);
Accum_Grad=[];
for layer=1:length(Frames{1})
    if isfield(Frames{1}(layer),'dzdw')&& ~isempty(Frames{1}(layer).dzdw)
        Accum_Grad(layer).dzdw=0;
        for f=1:n_frames
            Accum_Grad(layer).dzdw=Accum_Grad(layer).dzdw+Frames{f}(layer).dzdw;
        end
        Accum_Grad(layer).dzdw=Accum_Grad(layer).dzdw./n_frames;
    end

    if isfield(Frames{1}(layer),'dzdb')
        if ~isempty(Frames{1}(layer).dzdb)
            Accum_Grad(layer).dzdb=0;
            for f=1:n_frames
                Accum_Grad(layer).dzdb=Accum_Grad(layer).dzdb+Frames{f}(layer).dzdb;
            end
            Accum_Grad(layer).dzdb=Accum_Grad(layer).dzdb./n_frames;
        else
            Accum_Grad(layer).dzdb=[];
        end
    end
end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function mObject3 = timeAlignedDirectSTRAIGHTmorphing(mObject1,mObject2,mRate,mixMethod);
%   Morphing based on time-aligned mixing of STRAIGHT parameters
%   mObject3 = timeAlignedDirectSTRAIGHTmorphing(mObject1,mObject2,mRate,mixMethod);

%   Designed and coded by Hideki Kawahara
%   28/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara

mObject3 = checkForSimilarity(mObject1,mObject2);
if length(mObject3) ==0;return;end;
dtFrame = mObject1.frameUpdateInterval;
endLocation1 = (length(mObject1.F0)-1)*dtFrame; % in ms
endLocation2 = (length(mObject2.F0)-1)*dtFrame; % in ms
timeAnchor1 = [0;mObject1.anchorTimeLocation;endLocation1];
timeAnchor2 = [0;mObject2.anchorTimeLocation;endLocation2];
locationOn1 = (0:length(mObject1.F0)-1)*dtFrame;
locationOn2 = (0:length(mObject2.F0)-1)*dtFrame;
mapFrom1to2 = interp1(timeAnchor1,timeAnchor2,locationOn1);
[nr1,nc1] = size(mObject1.spectrogram);
[nr2,nc2] = size(mObject2.spectrogram);

%---- mixing on mObject1's time axis
nAxis1 = length(locationOn1);
nAxis2 = length(locationOn2);
morphedF0 = zeros(nAxis1,1);
morphedAp = zeros(nr1,nAxis1);
morphedSgram = zeros(nr1,nAxis1);
weightSumF0 = zeros(nAxis1,1);
for ii=1:nAxis1
    mappedIndexOn2 = mapFrom1to2(ii)/dtFrame+1;
    iFloor = floor(mappedIndexOn2);
    iFraction = mappedIndexOn2-iFloor;
    dAp = iFraction*(mObject2.aperiodicityIndex(:,min(iFloor+1,nAxis2))-mObject2.aperiodicityIndex(:,iFloor));
    morphedAp(:,ii) = (1-mRate)*mObject1.aperiodicityIndex(:,ii)+mRate*(mObject2.aperiodicityIndex(:,iFloor)+dAp);
    switch mixMethod
        case 'linear'
            dSgram = iFraction*(mObject2.spectrogram(:,min(iFloor+1,nAxis2))-mObject2.spectrogram(:,iFloor));
            morphedSgram(:,ii) = (1-mRate)*mObject1.spectrogram(:,ii)+mRate*(mObject2.spectrogram(:,iFloor)+dSgram);
        case 'log'
            dSgram = iFraction*(log(mObject2.spectrogram(:,min(iFloor+1,nAxis2)))-log(mObject2.spectrogram(:,iFloor)));
            tmp = (1-mRate)*log(mObject1.spectrogram(:,ii))+mRate*(log(mObject2.spectrogram(:,iFloor))+dSgram);
            morphedSgram(:,ii) = exp(tmp);
    end;
    if mObject1.F0(ii)>0
        morphedF0(ii) = (1-mRate)*log(mObject1.F0(ii));
        weightSumF0(ii) = (1-mRate);
    end;
    if (mObject2.F0(iFloor)>0) & (mObject2.F0(min(iFloor+1,nAxis2))>0)
        dF0 = iFraction*(log(mObject2.F0(min(iFloor+1,nAxis2)))-log(mObject2.F0(iFloor)));
        morphedF0(ii) = mRate*(log(mObject2.F0(iFloor))+dF0)+morphedF0(ii);
        weightSumF0(ii) = weightSumF0(ii)+mRate;
    end;
end;
morphedF0(weightSumF0>0) = exp(morphedF0(weightSumF0>0)./weightSumF0(weightSumF0>0));

%----- mapping back onto morphed time axis
timeAnchorMorph = (1-mRate)*timeAnchor1 + mRate*timeAnchor2;
locationOnMorph = (0:(timeAnchorMorph(end)/dtFrame))*dtFrame;
mapFormMorphTo1 = interp1(timeAnchorMorph,timeAnchor1,locationOnMorph);
nAxisMorph = length(locationOnMorph);
morphedApOnMorph = zeros(nr1,nAxisMorph);
morphedSgramOnMorph = zeros(nr1,nAxisMorph);
morphedF0onMorph = zeros(nAxisMorph,1);
for ii=1:nAxisMorph
    mappedIndexOn1 = mapFormMorphTo1(ii)/dtFrame+1;
    iFloor = floor(mappedIndexOn1);
    iFraction = mappedIndexOn1-iFloor;
    morphedApOnMorph(:,ii) = morphedAp(:,iFloor) ...
        +iFraction*(morphedAp(:,min(iFloor+1,nAxis1))-morphedAp(:,iFloor));
    morphedSgramOnMorph(:,ii) = morphedSgram(:,iFloor) ...
        +iFraction*(morphedSgram(:,min(iFloor+1,nAxis1))-morphedSgram(:,iFloor));
    if (morphedF0(iFloor)>0) & (morphedF0(min(iFloor+1,nAxis1))>0)
        dF0 = iFraction*(morphedF0(min(iFloor+1,nAxis1))-morphedF0(iFloor));
        morphedF0onMorph(ii) = morphedF0(iFloor)+dF0;
    end;
end;
mObject3.F0 = morphedF0onMorph;
mObject3.aperiodicityIndex = morphedApOnMorph;
mObject3.spectrogram = morphedSgramOnMorph;
mObject3.anchorTimeLocation = timeAnchorMorph(2:end-1);
mObject3.anchorFrequency = (1-mRate)*mObject1.anchorFrequency+mRate*mObject2.anchorFrequency;
%mObject3 = morphedAp; % This line is a dummy.

%%% ------ Internal function to check for object's similarity
function mObject3 = checkForSimilarity(mObject1,mObject2)
mObject3 = [];
if mObject1.samplingFrequency ~= mObject2.samplingFrequency;mObject3 = [];return;end;
if mObject1.frameUpdateInterval ~= mObject2.frameUpdateInterval;mObject3 = [];return;end;
if length(mObject1.anchorTimeLocation) ~= length(mObject2.anchorTimeLocation);mObject3 = [];return;end;
nAnchor = length(mObject1.anchorTimeLocation);
for ii=1:nAnchor % check for similarity of anchor structure
    frequencyAnchor1 = mObject1.anchorFrequency;
    frequencyAnchor2 = mObject2.anchorFrequency;
    if (sum(frequencyAnchor1>0) ~= sum(frequencyAnchor2>0)) | ...
            (sum(frequencyAnchor1<0) ~= sum(frequencyAnchor2<0))
        return;
    end;
end;
mObject3 = createMobject;

author: HidekiKawahara


text:
clear all;
clc

% File folder and name of int16 input signal
int16FileFolder = '/Users/';
int16FileName = 'inputsignal';

fileName = [int16FileFolder int16FileName '.int16'];

nMics = 128; % # mics in array
fs = 44.1e3; %sampling rate

data = readInt16(nMics, fs, 0, fileName);

% Get microphone positions
%Here the microphone positions needs to be loaded. That is mic #1
%corresponds to the first row in the timeSignal matrix, mic #2 to the
%second row and so on
Xm = zeros(1, nMics);
Ym = zeros(1, nMics);
Zm = zeros(1, nMics);
Wm = ones(1, numel(Xm))/numel(Xm);

% Scanning range (multiple frequencies, broadband)
fmin = 3e3;
fmax = 5e3;
c = 340; %speed of sound

t_start = 0.0;
T = 0.035; %time weighting
block = T*fs;

% Scanning parameters
distance = 2.0; %distance to source
maxX = 1; %max scanning extent x
maxY = 1; %max scanning extent y
deltaX = 0; %offset scanning grid
deltaY = 0; %offset scanning grid
anglRes = 1; %angle resolution
filterFrequencies = [fmin, fmax];

tic
[anglex,angley, pow] = sweepPow2(Xm, Ym, Zm, Wm, data(:, (t_start*fs)+1:end), fs, filterFrequencies, distance, maxX, maxY, anglRes, block, deltaX, deltaY);
toc

figure;
imagesc(angley, anglex, db(pow'))

author: jorgengrythe

repo: beamforming