text:
function t = uminus(t)
%UMINUS Unary minus (-) for sptenmat.
%
%   See also SPTENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

t.vals = -t.vals;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%demo_SOFAstrings - Script for testing the string array feature of SOFA.

% #Author: Piotr Majdak
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
% #Author: Michael Mihocic: updated with variable ReceiverDescriptions instead of Ears (03.08.2022)
%
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or – as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% 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.

%% Test Strings as application-specific variable
% Load some arbritrary HRTFs
hrtf = SOFAload(fullfile(SOFAdbPath, 'database','ari','dtf_nh2.sofa'));
% Add a string array
str={};
for ii=1:hrtf.API.M
  str{ii,1}=['String' num2str(round(rand(1,1)*10000))];
end
% SOFAaddVariable(Obj,Name,Dim,Value)
hrtf2 = SOFAaddVariable(hrtf,'Test','MS',str);

% Add a new string with dimensions [RS]
strn={'left ear'; 'right ear'};
hrtf2 = SOFAaddVariable(hrtf2, 'ReceiverDescriptions', 'RS', strn);

% Save as SOFA
SOFAsave('stringtest_applicationvar.sofa',hrtf2);
% Reload the file
hrtf = SOFAload('stringtest_applicationvar.sofa');
% compare the strings
if prod(strcmp(hrtf.Test,hrtf2.Test))
    disp('SimpleFreeFieldHRIR: String Load-Reload: OK');
    delete('stringtest_applicationvar.sofa');
else
    error('String comparison showed differences');
end
clear

%% Test with conventions GeneralString (non-standardized convention, just for testing)
% Create an empty object
Obj = SOFAgetConventions('GeneralString');
% Create numeric data with M=15, R=2, N=10
Obj.Data.Double=rand(15,2,10);
% Create string arrays
str2={}; str={};
for ii=1:15
  id = num2str(round(rand(1,1)*1000000));
  str{ii,1}=['X' id];
  str2{ii,1}=['Left' id];
  str2{ii,2}=['Right' id];
end
Obj.String2 = str2;      % String1=[MRS]
Obj.Data.String1 = str;  % Data.String1=[MS]
Obj.Data.String2 = str2; % Data.String2=[MRS]

% Add a new string with dimensions [RS]
strn={'left ear'; 'right ear'};
Obj = SOFAaddVariable(Obj, 'ReceiverDescriptions', 'RS', strn);

% Update dimensions
Obj = SOFAupdateDimensions(Obj);
% Save as SOFA
SOFAsave('stringtest_generalstring.sofa',Obj);
% Reload the file
Obj2 = SOFAload('stringtest_generalstring.sofa');
% Compare the strings
if ~prod(strcmp(Obj2.Data.String2,Obj.Data.String2))
    error('Data.String2: Comparison showed differences');
end
if ~prod(strcmp(Obj2.String2,Obj.String2))
    error('String2: Comparison showed differences');
end
if ~prod(strcmp(Obj2.Data.String1,Obj.Data.String1))
    error('Data.String1: Comparison showed differences');
end
if ~prod(strcmp(Obj2.ReceiverDescriptions,Obj.ReceiverDescriptions))
    error('ReceiverDescriptions: Comparison showed differences');
end
disp('GeneralString: String1, String2, Data, ReceiverDescriptions: Load-Reload: OK');
clear
delete('stringtest_generalstring.sofa');

author: sofacoustics

repo: SOFAtoolbox


text:
function p1 = gibbs(p1, x)
%GIBBS      Gibbs sampling
%
%           Description
%           P1 = GIBBS(P1, X) Samples one sample for parameters
%           contained in structure P1 at level n and and gets a
%           parameter structure X from level n-1 in hierarchy.
%           Returns a structure P1 containing new samples for
%           parameters in structure P1.
%
%           If P1 contains higher level parameters, e.g P1 =
%           net.p.w{1} and though there is a field P1.p, the
%           function calls itself to sample for those higher level
%           parameters also. So function takes care of the whole
%           hierarchical structure of ARD if called for the lowest
%           level parameters.
%
% Copyright (c) 1999-2003 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

%Matlab 6.5
%if any(strcmp(fieldnames(p1),'p')) & ~isempty(p1.p) %isfield(p1,'p')
%  for a = fieldnames(p1.a)';
%    a=a{:};
%    if any(strcmp(fieldnames(p1.p),a)) %isfield(p1.p,a);
%      p2=p1.p.(a);
%      % call the spesified distribution sampler
%      p1.a.(a)=feval(['cond_' p1.f '_' p2.f],a,p1.a,p2.a,x);
%      % sample the upper level parameters
%      p1.p.(a)=gibbs(p2, p1.a.(a));
%    end
%  end
%end

% Pre matlab 6.5
S.type='.';
if any(strcmp(fieldnames(p1),'p')) & ~isempty(p1.p) %isfield(p1,'p')
  for a = fieldnames(p1.a)';
    a=a{:};
    S.subs=a;
    if any(strcmp(fieldnames(p1.p),a)) %isfield(p1.p,a);
      p2=subsref(p1.p,S);
      p1.a=subsasgn(p1.a, S, feval(['cond_' p1.f '_' p2.f],a,p1.a,p2.a,x));
      p1.p=subsasgn(p1.p, S, gibbs(p2, subsref(p1.a,S)));
    end
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function LLF = cauchy_log(parameters,data)
%---------------------------------------------------
% PURPOSE:
%  Log Likelihood function of the Cauchy-Lorentz Distribution,
%  to estimate the location and scale parameters of a data set
%---------------------------------------------------
% USAGE:
% LLF = cauchy_log(x, data)
%
% INPUTS:
% parameters: a vector of parameters of the form [m; s]
% data:       the data set
%
% OUTPUTS:
% LLF:      The log-likelihood function
%---------------------------------------------------
% Author:
% Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
% Date: 06/2010
%---------------------------------------------------
[a,b]=size(parameters);
if b>a
   parameters=parameters';
end
n=length(data);
LLF  = -n*log(pi) +n*log(parameters(2)) - sum(log(parameters(2)^2 + (data-parameters(1)).^2));
end

author: Sable

repo: mcbench-benchmarks


text:
function varargout = vscale(varargin)
%VSCALE   Vertical scale of a SPHEREFUN.
%
% VSCL = VSCALE(F) returns the vertial scale of a SPHEREFUN as determined
% by evaluating on a coarse tensor-product grid.

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

[varargout{1:nargout}] = vscale@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [rInt,vInt] = Hill2ECI_Vectorized(rTgt,vTgt,rHill,vHill)
%% Purpose:
% Convert those position (rHill) and velocity (vHill) values back into an
% ECI coordinate frame of reference using the reference satellite
% (rTgt,vTgt) position and velocity data.
%
%% Inputs:
%rTgt                       [3 x N]                 ECI Position vector of
%                                                   reference frame (km)
%
%vTgt                       [3 x N]                 ECI Velocity vector of
%                                                   reference frame (km/s)
%
%rHill                      [3 x N]                 Hill's relative
%                                                   position vector (km)
%
%vHill                      [3 x N]                 Hill's relative
%                                                   velocity vector (km/s)
%
%
%
%% Outputs:
%rInt                       [3 x N]
%
%vInt                       [3 x N]
%
%
% References:
% Vallado 2007.
% Programed by Darin C Koblick 11/30/2012
%% Begin Code Sequence
%Declare Local Functions
rTgtMag = sqrt(sum(rTgt.^2,1));
vTgtMag = sqrt(sum(vTgt.^2,1));
matrixMultiply = @(x,y)permute(cat(2,sum(permute(x(1,:,:),[3 2 1]).*permute(y,[2 1]),2), ...
                                     sum(permute(x(2,:,:),[3 2 1]).*permute(y,[2 1]),2), ...
                                     sum(permute(x(3,:,:),[3 2 1]).*permute(y,[2 1]),2)),[2 1]);
%Find the RSW matrix from the target ECI positions
RSW = ECI2RSW(rTgt,vTgt); rIntMag = rTgtMag + rHill(1,:);
%Compute rotation angles to go from tgt to interceptor
lambda_int = rHill(2,:)./rTgtMag; phi_int = sin(rHill(3,:)./rTgtMag);
CLI = cos(lambda_int); SLI = sin(lambda_int); CPI = cos(phi_int); SPI = sin(phi_int);
%find rotation matrix to go from rsw to SEZ of inerceptor
RSW_SEZ = zeros(3,3,size(rTgt,2));
RSW_SEZ(1,1,:) = SPI.*CLI;  RSW_SEZ(1,2,:) = SPI.*SLI; RSW_SEZ(1,3,:) = -CPI;
RSW_SEZ(2,1,:) = -SLI;      RSW_SEZ(2,2,:) = CLI;      RSW_SEZ(3,1,:) = CPI.*CLI;
                            RSW_SEZ(3,2,:) = CPI.*SLI; RSW_SEZ(3,3,:) = SPI;
%Find velocity component positions by using angular rates in SEZ frame
vIntSEZ = cat(1,-rIntMag.*vHill(3,:)./rTgtMag, ...
                 rIntMag.*(vHill(2,:)./rTgtMag + vTgtMag./rTgtMag).*CPI, ...
                 vHill(1,:));
vInt = matrixMultiply(permute(RSW,[2 1 3]), ...
       matrixMultiply(permute(RSW_SEZ,[2 1 3]), ...
       vIntSEZ));
%Find the position components
rIntRSW = bsxfun(@times,rIntMag,cat(1,CPI.*CLI, ...
                                      CPI.*SLI, ...
                                      SPI));
rInt = matrixMultiply(permute(RSW,[2 1 3]),rIntRSW);
end


text:
classdef tuneAD9361AGC < AD9361TRx & gen80211aTestWaveform & demod80211aTestWaveform
    properties (Access = public)
        %%
        % Toggle SIM_STUDY to control whether a simulation study or
        % a physical experiments-based study using a radio is conducted
        % SIM_STUDY - true/false
        SIM_STUDY = true
        % SIM_MODE controls the simulation mode as follows:
        %   value  |     channel     |   radio
        %---------------------------------------
        %     0    |   AWGN          |    No
        %     1    |   AWGN          |    Yes
        %     2    |   AWGN+Fading   |    No
        %     3    |   AWGN+Fading   |    Yes
        SIM_MODE = 3
        % GAIN_MODE controls whether the simulated gain applied to each
        % packet is constant or varying (applicable when SIM_STUDY is
        % set to true and SIM_MODE is set to 1 or 3)
        % 0 - uniform
        % 1 - non-uniform
        % If GAIN_MODE is set to 1, but SIM_MODE is set to 0 or 2, no AGC
        % is applied to the received signal and therefore, more bad packets
        % might be received than what is expected.
        GAIN_MODE = 1
        % SNR (dB) of the received test waveform (applicable when SIM_STUDY
        % is set to true)
        snr
    end

    methods
        %%
        % sim_settings needs to have the following fields:
        % SIM_STUDY - see above
        % SIM_MODE - see above
        % GAIN_MODE - see above
        % snr - SNR (dB)
        %
        % wlan_settings needs to have the following fields:
        % fc - carrier frequency
        % numPackets - number of WLAN frames
        % mcs - WLAN's MCS Index
        % seed - 0/xyz (0 - random seed; xyz - seed value)
        % See gen80211aTestWaveform.m for more details.
        %
        % ad9361_settings needs to have the following fields:
        % AGC_MODE - AGC mode (manual/slow/fast/receiver control)
        % LOG_ADC_OUTPUT - true/false (Toggle to log ADC block - It is
        % suggested that this quantity is set to false unless the number of
        % simulated packets is small. ADC runs at a higher rate than the
        % sample clock. So, the size of the log file will be significantly
        % larger if set to true.)
        % SAVE_LOG_DATA - true/false (Toggle to save log data to file)
        %
        % See AD9361TRx.m for more details.
        function obj = tuneAD9361AGC(sim_settings, wlan_settings, varargin)
            % generate IEEE 802.11a compliant test waveform
            obj = obj@gen80211aTestWaveform(sim_settings, wlan_settings);

            % apply AGC settings to the model and run simulations
            obj = obj@AD9361TRx(sim_settings, varargin);

            % demodulate IEEE 802.11a compliant test waveform
            obj = obj@demod80211aTestWaveform(sim_settings);

            obj.SIM_STUDY = sim_settings.SIM_STUDY;
            obj.SIM_MODE = sim_settings.SIM_MODE;
            obj.GAIN_MODE = sim_settings.GAIN_MODE;
            obj.snr = sim_settings.snr;
        end
    end
end

author: analogdevicesinc


text:
function changes = get_program_changes(file_name)
% Gets the program changes in a midi file.
%
%   changes = get_program_changes(file_name)
%
% INPUTS:
%   file_name - the name of the midi file
%
% OUTPUTS:
%   changes - an Nx5 array of the patch changes.  Each row
%     is a single patch changes.  The columns are:
%        (1) the patch (instrument) number
%        (2) the time of the command in beats
%        (3) the time of the command in seconds
%        (4) the track number
%        (5) the channel
%
% 2010-05-03 Christine Smit csmit@ee.columbia.edu

import edu.columbia.ee.csmit.MidiKaraoke.read.*;
import java.io.File;
import javax.sound.midi.*;

midiFile = File(file_name);
seq = MidiSystem.getSequence(midiFile);

% get the number of ticks/quarter note, which I assume is the
% 'beat' in the nm
ticksPerQuarterNote = seq.getResolution();

changesInMidi =  ProgramChangeViewParser.parse(seq);

changes = changesInMidi.getProgramChangesDoubles();

% add 1 to the patch numbers
changes(:,1) = changes(:,1) +1;
% convert the ticks to beats...
changes(:,2) = changes(:,2)./ticksPerQuarterNote;
% add 1 to the track number
changes(:,4) = changes(:,4)+1;
% add 1 to the channel numbers
changes(:,5) = changes(:,5)+1;

end


text:
function results = rpca_for_tensor( data, params )
% params.mode = 0: apply RPCA to each frontal slice of the tensor. only
% supported for IsTC = false.
% params.mode = i: apply RPCA to the i-th mode unfolding of the tensor
% Uses IALM by Yi Ma's group
% For partial observations case, RPCA does not support mode 0 yet.

tic;
if ~params.IsTC
switch params.mode
    case 0
        iters = [];
        for k = 1:size(data.T,3)
            T = double( data.T(:,:,k) );
            lambda = 1 / sqrt(max(size(T)));
            [Xk Ek iter] = inexact_alm_rpca(T, lambda*params.rRatio, params.opt_tol, params.max_iter, 1/(params.mu1fac*std(T(:))) );
            results.X(:,:,k) = Xk;
            results.E(:,:,k) = Ek;
            iters(k) = iter;
        end
        results.iter = mean(iters);
    otherwise
%     N = length(size(double(data.T)));
        Tmat = tenmat(data.T,params.mode);
        T = double( Tmat );
        lambda = 1 / sqrt(max(size(T)));
        [X E iter] = inexact_alm_rpca(T, lambda*params.rRatio, params.opt_tol, params.max_iter, 1/params.mu1 );
        results.X = tensor(tenmat( X, Tmat.rdims, Tmat.cdims, Tmat.tsize) );
        results.E = tensor(tenmat( E, Tmat.rdims, Tmat.cdims, Tmat.tsize) );
        results.iter = iter;
end

else
    % RPCA with missing data
    results = rpca_tc( data, params );
end

results.cpu = toc;
results.IsTC = params.IsTC;

end

% ======================================================================= %
function results = rpca_tc( data, params )
% solves
%   min_{X,E} ||X||_* + \lambda_1||E||_1
%   s.t.      Ac(Y+E) = Tc
%             Y = X
%
% data.matInd: indices in Omega w.r.t. i-th mode of T

Tmat = tenmat(data.T,params.mode);
Xmat = double(tenmat(data.X,params.mode));
[ n, m ] = size(double(Tmat));

data = tenInd2matInd_core( data, params.mode );
Omega = data.matInd;
b = double(Tmat(Omega));
X = zeros( n, m );
Y = X;
E = zeros( n, m );

%%%%%%%% Fix lambda value here!! %%%%%%%%%%
% lambda = params.lambda;
lambda = params.rRatio / sqrt(max(m,n));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu = params.mu1;

V = cell(1,2);
V{1} = zeros( size(b) );
V{2} = zeros( size(X) );
bnorm2 = norm(b)^2;

% assign operators
Aprod = @(X_param)A_select( X_param, Omega );
Atprod = @(b_param)At_select( b_param, Omega, size(E) );

for iter = 1:params.max_iter
    % solve X_i's
    P = Y - mu*V{2};
    X = matrix_shrinkage( P, mu );

    % solve E
    p = b + mu*V{1} - Aprod(Y);     P = Atprod(p);
    E(:) = shrinkage_v( P(:), mu*lambda );

    % solve Y
    Yp = Y;
    Y = X + mu*V{2};
    Y(Omega) = (Y(Omega) - E(Omega) + b + mu*V{1}) / 2;

    % compute optimality stats
    tdiff{1} = Aprod( Y + E ) - b;
    tdiff{2} = Y - X;
    Ydiff = Y - Yp;
    YdiffOmega = Ydiff(Omega);
    YpOmega = Yp(Omega);

    pres = norm( [tdiff{1}; tdiff{2}(:)] ) / norm( [b; X(:)] );
    dres = norm( [YdiffOmega(:); Ydiff(:)] ) / norm( [YpOmega(:); Y(:)] );

    rel_err = norm(Y-Xmat, 'fro') / norm(Xmat,'fro');

    % print
    if params.verbose
        fprintf('Iter: %d,   pinf: %3.2e,   dinf: %3.2e,    rel_err: %3.2e\n', iter, pres, dres, rel_err );
    end

    if max(pres, dres) < params.opt_tol
%     if pres < params.opt_tol
        break;
    end

    % update Lagrange multipliers
    V{1} = V{1} - tdiff{1}/mu;
    V{2} = V{2} - tdiff{2}/mu;

%     mu = max(mu / 1.5, params.mu_min);
end

results.X = tensor( tenmat( X, Tmat.rdims, Tmat.cdims, Tmat.tsize ) );
results.E = tensor( tenmat( E, Tmat.rdims, Tmat.cdims, Tmat.tsize ) );
results.V = V;
results.b = b;
results.iter = iter;
results.mu = mu;
results.lambda = lambda;

end

function AX = A_select( X, Omega )
% AX is a vector

AX = X(Omega);

end

function X = At_select( b, Omega, sz )

X = zeros( sz );
X(Omega) = b;
end

repo: lrslibrary


text:
report_this_filefun(mfilename('fullpath'));

% TODO delete this, probably. -CGR
% Now lets plot the color-map of the z-value
%
figure

set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
    'FontWeight','bold','LineWidth',1.5,...
    'Box','on','SortMethod','childorder')

rect = [0.18,  0.10, 0.7, 0.75];
rect1 = rect;

% set values greater tresh = nan
%
re4 = ret;l = r > tresh; re4(l) = zeros(1,length(find(l)))*nan;

% plot image
%
orient portrait
set(gcf,'PaperPosition', [2. 1 7.0 5.0])

l = isnan(re4);

axes('position',rect);hold on

pco1 = pcolor(gx,gy,(re4))
caxis([ 25 150]);
shading interp
axis([ min(gx) max(gx) min(gy) max(gy)]);
axis image;
hold on
h = cool(64);h = [  h]; colormap(h)

set(gca,'Color',[0.9 0.9 0.9])
set(gca,'YTickLabels',[  10 8  6 4  2])

%xlabel('Distance in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
%ylabel('depth in [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

if exist('maex', 'var')
    pl = plot(maex,-maey,'*y');
    set(pl,'MarkerSize',6,'LineWidth',2)
end
overlay

set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
    'FontWeight','normal','LineWidth',1.5,...
    'Box','on','TickDir','out')
h1 = gca;hzma = gca;

% Create a colobar
%
% h5 = colorbar('horiz');
%set(h5,'Pos',[0.25 0.05 0.5 0.05],...
%'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

% Make the figure visible
%
set(gca,'visible','on','FontSize',10,'FontWeight','bold',...
    'FontWeight','normal','LineWidth',1.0,...
    'Box','on','TickDir','out')

axes(h1)

author: CelsoReyes


text:
function [xform, fiterr] = affineSolve(A,B);
% Solve the equation coordsA = xform * coordsB for the xform.
%
% Usage: [xform, fiterr] = affineSolve(coordsA,coordsB);
%
% Problem: You have coordinates of the same points in two
% coordinate spaces, A and B. You want to find a 4x4 affine
% transform matrix that optimally maps between these points.
%
% Sub-problem: Simply dividing coordsA / coordsB returns a
% 3x3 xform, not a 4x4 affine.
%
% Solution: use this. It uses singular value decomposition to
% find the rotation, scales (and skews, if the points are not
% really matching), and subtraction to find the translations,
% and builds a 4x4 affine xform which minimizes the error
% in the above equation.
%
% coordsA and coordsB should be 3xN coordinate matrices with the
% same number of columns. Each column of coordsA should reflect
% the (row, column, slice) of a point in coordinate space A (the
% NEW coordinate space being xformed to), and each column in coordsB
% should be the (row, column, slice) of the SAME POINT in coordinate
% space B (the OLD coordinate space being xformed from).
%
% The optional output argument fiterr is the root mean squared error
% between coordsA and coordsB given the best set of rotations and
% tanslations. It is a measure of the goodness-of-fit of the xform
% (since, if the points don't really correspond, not all equations
% can be perfectly solved).
%
% ras, 10/2005.
if nargin<2, help(mfilename); error('Not enough args.'); end

% size + number checks:
if ~isnumeric(A) | ~isnumeric(B), error('Need numeric args.'); end
if size(A,2)~=size(B,2)
    error('Coords should have same # of columns.');
end
if size(A,1)<3 | size(B,1)<3,
    error('Coords need to be specified in 3 dimensions.')
end

% just curious: does just doing / division work better?
n = size(A, 2);
xform = [A; ones(1, n)] / [B; ones(1, n)];

% % express points relative to center of all points:
% nPoints = size(A,2);
% centeredA = A - repmat(mean(A,2),[1 nPoints]);
% centeredB = B - repmat(mean(B,2),[1 nPoints]);
%
% % solve for rotation + scales (& skews, if points are bad):
% H = zeros(3,3);
% for i = 1:nPoints
%     H = H + (centeredB(:,i) * centeredA(:,i)');
% end
% [U S V] = svd(H);
% rot = V*(U');
%
% % solve for translation
% rotatedB = rot * B;
% trans = mean(A,2) - mean(rotatedB,2);
%
% % build 4x4 affine xform matrix
% xform = [rot trans; 0 0 0 1];
%
% % force skew to equal 0
% [trans rot scale skew] = affineDecompose(xform);
% skew = [0 0 0];
% xform = affineBuild(trans,rot,scale,skew);
%
% % return error of fit if requested
% if nargout > 1
% 	xB = xform * [B; ones(1, size(B, 2))];
% 	xB = xB(1:3,:);
% 	fiterr = mean(abs(A - xB).^2);
% end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function events = sta_lta(wave,varargin)

%STA_LTA: Short-Time-Average/Long-Time-Average event detector.
%
%USAGE: events = sta_lta(wave,prop_name_1,prop_val_1,...)
%
%DIAGRAM:
%                       /\      /\/\        /\
%              /\  /\/\/  \/\  /    \/\    /  \/\
%                \/          \/        \/\/      \
%                                           |-STA-| -->
%         --> |-----------------LTA---------------| -->
%
%INPUTS: wave     - A waveform object containing events (maybe)
%        varargin - User-defined parameter name/value pairs (below)
%
%VALID PROP_NAME:
%
%  'edp'      - Event Detection Parameters (Detailed in 'VALID PROP_VAL')
%  'lta_mode' - LTA window behavior following a 'trigger on' flag
%  'skip'     - Times for STA/LTA to skip over (if known), examples are
%               data gaps, calibration pulses, periods of excessive noise
%  'eot'      - Event output type (Detailed in 'VALID PROP_VAL')
%  'fix'      - Output event length fixed or variable?
%  'pad'      - Pad output events by a fixed amount of time before start
%               time (trigger on) and after stop time (trigger off).
%
%VALID PROP_VAL:
%
%  'edp'(1x6 numeric) ... default: [1 8 2 1.6 0 3]
%     --> [l_sta l_lta th_on th_off min_sep min_dur]
%        --> l_sta    - STA window length (s)
%        --> l_lta    - LTA window length (s)
%        --> th_on    - STA/LTA trigger on threshold
%        --> th_off   - STA/LTA trigger off threshold
%        --> min_sep  - Minimum seperation between trigger off time of an
%                       event and trigger on time of the next event (s)
%        --> min_dur  - Minimum event duration to be recorded(s)
%
%	'lta_mode' (string) ... default: 'continuous'
%     --> 'frozen' - LTA window fixed in place after trigger is turned on
%                    while STA window continues forward.
%     --> 'continuous' - LTA window continues w/ STA window after trigger
%                        is turned on (Same behavior as before trigger)
%     --> 'grow' - LTA left edge fixed, right edges continues w/ STA right
%                  edge after trigger on, thus size of LTA window is
%                  growing until trigger off, at which point LTA window
%                  returns to original size.
%
%	'skip' (nx2 numeric) ... default: [] (don't skip anything)
%     --> skip_val - List of start/stop times to skip
%
%	'eot' (string) ... default: 'sst'
%     --> 'st'  - Return event start times (nx1 numeric)
%                 (a.k.a. trigger on times)
%     --> 'sst' - Return event start/stop times (nx2 numeric)
%                 (a.k.a. trigger on/off times)
%     --> 'ssd' - Return event start/stop datapoints (nx2 integer)
%                 (referenced from first datapoint in 'wave')
%     --> 'wfa' - Return event waveform array (1xn waveform)
%
%   'fix' (1x1 numeric) ... default: 0
%     --> 'fix_val' - If 0, output events will be variable length
%                     (trigger on time to trigger off time)
%                     If a positive numeric value N, output events will be
%                     fixed length (trigger on time to fixed length N)
%                     Units of N is seconds if non-zero
%
%   'pad' (1x2 numeric) ... default: [0 0] (i.e. no padding)
%     --> 'pad_val' - If [0 0] events will have no time padding
%                     [1 2] - pad by 1s before, 2s after event
%                     this will pad events regarless of the value of fix:
%                     i.e. fix = 7, pad = [1 2], returned events are 10s
%
%OUTPUTS: events - events in format specified by 'return'

% Author: Dane Ketner, Alaska Volcano Observatory
% $Date$
% $Revision$

%% Check waveform variable
if isa(wave,'waveform')
   Fs = get(wave,'freq');         % Sampling frequency
   l_v = get(wave,'data_length'); % Length of time series
   tv = get(wave, 'timevector');  % Time vector of waveform
   if isempty(wave)
       disp('Input waveform empty, no events detected')
       events = [];
       return
   end
else
   error('STA_LTA: First Argument must be a waveform object')
end

%% Set all default parameters
l_sta = 1*Fs;     % STA window length
l_lta = 8*Fs;     % LTA window length
th_on = 2;        % Trigger on when sta_to_lta exceeds this theshold
th_off = 1.6;     % Trigger off when sta_to_lta drops below threshold
min_sep = 0*Fs;   % Skip ahead after end of event
min_dur = 3*Fs;   % Any triggers shorter than min_dur are discarded

lta_mode = 'continuous'; % Post trigger-on LTA behavior
skip = [];               % List of start/stop times to skip
eot = 'sst';             % Event output type
fix = 0*Fs;              % Fix event length
pad = [0 0]*Fs;          % Pad event start/stop times

%% Check varargin size
nv = numel(varargin);
if ~rem(nv,2) == 0
   error(['STA_LTA: Arguments after wave must appear in ',...
          'property_name/property_value pairs'])
end

%% User-defined parameters (varargin)
if nv > 0
    for p = 1:2:nv-1
        v1 = varargin{p};
        v2 = varargin{p+1};
        switch lower(v1)
            case 'edp'
                if isnumeric(v2) && numel(v2) == 6
                    l_sta = v2(1)*Fs;    % STA window length
                    l_lta = v2(2)*Fs;    % LTA window length
                    th_on = v2(3);       % Trigger on theshold
                    th_off = v2(4);      % Trigger off threshold
                    min_sep = v2(5)*Fs;  % Skip ahead after end of event
                    min_dur = v2(6)*Fs;  % Minimum event duration
                else
                    error('STA_LTA: Wrong format for input ''edp''')
                end
            case 'lta_mode'
                switch lower(v2)
                    case {'freeze','frozen'}
                        lta_mode = 'frozen';
                    case {'continue','continuous'}
                        lta_mode = 'continuous';
                    case {'grow','growing'}
                        lta_mode = 'grow';
                    otherwise
                      error('STA_LTA: Wrong format for input ''lta_mode''')
                end

            case 'skip'
                if isnumeric(v2) && size(v2,2) == 2
                    skip = v2;
                else
                    error('STA_LTA: Wrong format for input ''skip''')
                end
            case 'eot'
                switch lower(v2)
                    case 'st'
                        eot = v2;
                    case 'sst'
                        eot = v2;
                    case 'ssd'
                        eot = v2;
                    case 'wfa'
                        eot = v2;
                    otherwise
                        error('STA_LTA: Wrong format for input ''eot''')
                end
            case 'fix'
                if isnumeric(v2) && numel(v2) == 1
                    fix = v2*Fs;
                else
                    error('STA_LTA: Wrong format for input ''fix''')
                end
            case 'pad'
                if isnumeric(v2) && size(v2,1) == 1 && size(v2,2) == 2
                    pad = v2*Fs;
                else
                    error('STA_LTA: Wrong format for input ''pad''')
                end
            otherwise
                error(['STA_LTA: ''edp'', ''lta_mode'', ''skip'' ',...
                    '''eot'', ''fix'', and ''pad'' are the only ',...
                    'property names recognized'])
        end
    end
end

%% Initialize waveform data
wave = set2val(wave,skip,NaN);   % NaN skip times
wave = zero2nan(wave,10);        % NaN any data gaps
v = get(wave,'data');            % Waveform data
abs_v = abs(v);                  % Absolute value of time series

%% Initialize flags and other variables
lta_calc_flag = false;       % has the full LTA window been calculated?
ntrig = 0;               % number of triggers
trig_array = zeros(1,2); % array of trigger times: [on,off;on,off;...]

%% Loops over data
% i is the primary reference point (right end of STA/LTA window)
i = l_lta+1;
while i <= l_v % START STA_LTA MAIN LOOP

%% Skip data gaps (NaN values in LTA window)?
   if any(isnan(abs_v(i-l_lta:i)))
      gap = true;
      lta_calc_flag = false; % Force full claculations after gap
      while (gap) && (i < l_v)
         i = i+1;
         if ~any(isnan(abs_v(i-l_lta:i)))
            gap = false;
         end
      end
   end

%% Calculate STA & LTA Sum (Do Full Calculation?)
   if ~lta_calc_flag
      lta_sum = 0;
      sta_sum = 0;
      for j = i-l_lta:i-1              % Loop to compute LTA & STA
         lta_sum = lta_sum + abs_v(j); % Sum LTA window
         if (i - j) <= l_sta           % Sum STA window (right side of LTA)
            sta_sum = sta_sum + abs_v(j);
         end
      end
      lta_calc_flag = true;
   else

%% Calculate STA & LTA Sum (Single new data point if not Full)
      lta_sum = lta_sum - abs_v(i-l_lta-1) + abs_v(i-1);
      sta_sum = sta_sum - abs_v(i-l_sta-1) + abs_v(i-1);
   end

%% Calculate STA & LTA
   lta = lta_sum/l_lta;
   sta = sta_sum/l_sta;

%% Calculate STA/LTA Ratio
   sta_to_lta = sta/lta;

%% Trigger on? (Y/N)
   if (sta_to_lta > th_on)
      j = i;   % Set secondary reference point = primary
      g = 0;   % l_lta growth, only used if LTA growing
      while (sta_to_lta > th_off)
         j = j+1;
         if j < l_v
            sta_sum = sta_sum - abs_v(j-l_sta-1) + abs_v(j-1);
            switch lta_mode
               case 'frozen'
                  % LTA is good just the way it is
               case 'continuous'
                  % Add new data point, remove oldest data point
                  lta_sum = lta_sum - abs_v(j-l_lta-1) + abs_v(j-1);
               case 'grow'
                  % Add new data point, increase
                  lta_sum = lta_sum + abs_v(j-1);
                  l_lta = l_lta + 1;
                  g = g+1;
            end
            sta = sta_sum/l_sta;
            lta = lta_sum/l_lta;
            sta_to_lta = sta/lta;
            if any(isnan(abs_v(j-l_sta:j))) % NaN gaps to skip?
               sta_to_lta = 0; % Force trigger off (data gap in STA window)
            end
         else
            sta_to_lta = 0; % Force trigger off (end of data)
         end
      end
      duration = (j-i); % span from trigger on to trigger off
      l_lta = l_lta-g;

%% Triggered period long enough? (Y/N)
      if duration > min_dur % If duration < min_dur then skip it
         trig_t = i-l_sta;  % Beginning of STA window during trigger on
         end_t  = j;        % End of STA window during trigger off
         ntrig = ntrig + 1; % Event counter
         trig_array(ntrig,:) = [trig_t, end_t];
      end
      i = j + min_sep;   % Skip ahead by minimum event seperation
      lta_calc_flag = false; % Reset LTA calc flag to force new computation
   end
   i = i + 1;
end % END STA_LTA MAIN LOOP

%% Return events
if (trig_array(1,1)==0)&&(trig_array(1,2)==0)
    disp('No events detected')
    events = [];
    return
end
if fix>0
    trig_array(:,2) = trig_array(:,1) + fix;
end
trig_array(:,1) = trig_array(:,1) - pad(1);
trig_array(:,2) = trig_array(:,2) + pad(2);
% If events are padded, make sure the first and last events are not out of
% range (beginning before or ending after the span of input 'wave')
trig_array(trig_array<1) = 1;      % Set to first datapoint
trig_array(trig_array>l_v) = l_v;  % Set to last datapoint
switch lower(eot)
    case {'st'}
        events = tv(trig_array(:,1));
    case {'sst'}
        events = tv(trig_array);
    case {'ssd'}
        events = trig_array;
    case {'wfa'}
        events = [];
        for n = 1:size(trig_array,1)
            events = [events; extract(wave,'INDEX',trig_array(n,1),...
                                                   trig_array(n,2))];
        end
end

author: geoscience-community-codes

repo: GISMO


text:
% c = str2cell(v)
% c = str2cell(v,delim[s])
%
%	creates a cell array <c> from input array <v>
%	seperated by delimiter[s] <del>
%
% v	: a string/vector of a ML supported data class
% delim	: array of delimiter[s]				[def: isspace]
%	  v-class	syntax
%	  ----------------------
%	  char		'xyz'
%	  other		[1 2 pi]
%
% note:
%	an input matrix (NxM) will be turned into a
%	1x(N*M) row vector WITHOUT warning!
%
% examples:
%	s='this is a ;;; test case; +-0 ;;;_;;; or; +;+ 1';
%	c=str2cell(s,' +-;0_')
%		'this'
%		'is'
%		'a'
%		'test'
%		'case'
%		'or'
%		'1'
%
%	d=[pi pi 1 2 pi 4 5 6 pi];
%	c=str2cell(d,[pi 5])
%		[1x2 double]	% contents: 1 2
%		[         4]
%		[         6]

% created:
%	us	13-Nov-2000
% modified:
%	us	20-Dec-2003		/ TMW
%	us	06-Jun-2005 20:05:27

%--------------------------------------------------------------------------------
function	c=str2cell(s,varargin)

	if	nargin < 1
		help str2cell
		return;
	end

% make sure we have a row vec
		s=s(:).';

% check input
	if	nargin < 2
	if	ischar(s)
		ix=~isspace(s);
	else
		ix=ones(size(s));
	end
	else
		ix=~ismember(s,varargin{1});
	if	~isempty(find(isnan(varargin{1})))
		ix=ix&~isnan(s);
	end
	end

% find indices of delim[s]
		dx=diff(ix);
		k=sort([findstr(dx,-1) findstr(dx,1)+1]);
	if	~isempty(k)
	if	ix(1) == 1
		k=[1 k];
	end
	if	ix(end) == 1
		k=[k length(ix)];
	end
	else
		k=1:length(s);
		k=sort([k k]);
	end
		k=reshape(k,2,[]).';

% create cell[s]
		nk=size(k,1);
		c=cell(nk,1);
	for	i=1:nk
		c{i}=s(k(i,1):k(i,2));
	end
		return;


text:
function geometry_test1745 ( )

%*****************************************************************************80
%
%% TEST1745 tests R8MAT_SOLVE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  nrhs = 2;

  a = [ ...
     1.0,  4.0,  7.0; ...
     2.0,  5.0,  8.0; ...
     3.0,  6.0,  0.0; ...
    14.0, 32.0, 23.0; ...
     7.0, 16.0,  7.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1745\n' );
  fprintf ( 1, '  R8MAT_SOLVE solves linear systems.\n' );
  fprintf ( 1, '\n' );
%
%  Print out the matrix to be inverted.
%
  r8mat_print ( n, n+nrhs, a, '  The linear system:' );
%
%  Solve the systems.
%
  [ a, info ] = r8mat_solve ( n, nrhs, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The input matrix was singular.\n' );
    fprintf ( 1, '  The solutions could not be computed.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The computed solutions:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    for j = 1 : nrhs
      fprintf ( 1, '  %12f', a(i,n+j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function model = gpdisimOptimise(model, display, iters);

% GPDISIMOPTIMISE Optimise the GPSIM model.
% FORMAT
% DESC optimises the Gaussian process single input motif model for
% a given number of iterations.
% ARG model : the model to be optimised.
% ARG display : whether or not to display while optimisation
% proceeds, set to 2 for the most verbose and 0 for the least
% verbose.
% ARG iters : number of iterations for the optimisation.
% RETURN model : the optimised model.
%
% SEEALSO : scg, conjgrad, gpdisimCreate, gpdisimGradient, gpdisimObjective
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007

% SHEFFIELDML

if nargin < 3
  iters = 2000;
  if nargin < 2
    display = 1;
  end
end

params = gpdisimExtractParam(model);

options = optOptions;
if display
  options(1) = 1;
  if length(params) <= 100
    options(9) = 1;
  end
end
options(14) = iters;

if isfield(model, 'optimiser')
  optim = str2func(model.optimiser);
else
  optim = str2func('conjgrad');
end

params = optim('gpdisimObjective', params,  options, ...
               'gpdisimGradient', model);

model = gpdisimExpandParam(model, params);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [tg,theta] = tgmo2(tmap,ntex,radius,norient,varargin)
% function [tg,theta] = tgmo2(tmap,ntex,radius,norient,...)
%
% Compute the texture gradient at a single scale and multiple
% orientations.
%
% INPUT
%	tmap		Texton map, values in [1,ntex].
%	ntex		Number of textons.
%	radius		Radius of disc for texture gradient.
%	norient		Number of orientation at which to compute
%			the texture gradient.
%	'smooth'	Smoothing method, one of
%			{'gaussian','savgol','none'}, default 'none'.
%	'sigma'		Sigma for smoothing, default to radius.
%	'tsim'		Texton similarity matrix.  If not
%			provided, then use chi-squared.
%
% OUTPUT
%	tg		Size [h w norient] array of tg images.
%	theta		Vector of disc orientations (which are
%			orthogonal to the texture gradient).
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

% process options
smooth = 'none';
sigma = radius;
usechi2 = true;
for i = 1:2:numel(varargin),
  opt = varargin{i};
  if ~ischar(opt), error('option names not a string'); end
  if i==numel(varargin), error(sprintf('option ''%s'' has no value',opt)); end
  val = varargin{i+1};
  switch opt,
   case 'smooth',
    switch val,
     case {'none','gaussian','savgol'}, smooth=val;
     otherwise, error(sprintf('invalid option smooth=''%s''',val));
    end
   case 'sigma', sigma=val;
   case 'tsim', tsim=val; usechi2=false;
   otherwise, error(sprintf('invalid option ''%s''',opt));
  end
end

radius = max(1,radius);
norient = max(1,norient);
theta = (0:norient-1)/norient*pi;

% check texton labels
if any(tmap~=round(tmap)),
  error('texton labels not integral');
end
if min(tmap(:)) < 1 | max(tmap(:))>ntex,
  error(sprintf('texton labels out of range [1,%d]',ntex));
end

% radius of discrete disc
wr = floor(radius);

% count number of pixels in a disc
[u,v] = meshgrid(-wr:wr,-wr:wr);
gamma = mod(atan2(v,u),2*pi);
mask = (u.^2 + v.^2 <= radius^2);
mask(wr+1,wr+1) = 0; % mask out center pixel to remove bias
count = sum(mask(:));

% determine which pie slice each pixel falls into
% 0=masked [1,2*norient]=slice
slice = 1 + floor(gamma/(pi/norient));
slice = slice .* mask;

[h,w] = size(tmap);
tg = zeros(h,w,norient);
fwrite(2,'[');
for x = 1:w,
  fwrite(2,'.');
  for y = 1:h,
    pie = zeros(ntex,2*norient);
    for u = -wr:wr,
      xi = x + u;
      if xi<1 | xi>w, continue; end
      for v = -wr:wr,
        yi = y + v;
        if yi<1 | yi>h, continue; end
        s = slice(v+wr+1,u+wr+1);
        if s==0, continue; end % masked out
        t = tmap(yi,xi);
        pie(t,s) = pie(t,s) + 1;
      end
    end
    pie = pie .* (2/count); % normalize
    % initialize left/right histograms
    lhist = sum(pie(:,1:norient),2);
    rhist = sum(pie(:,norient+1:end),2);
    % spin the disc to compute tg at each orientation
    for i = 1:norient,
      if usechi2,
        chi = (lhist-rhist).^2 ./ (lhist+rhist+eps);
        tg(y,x,i) = 0.5*sum(chi);
      else
        lrdiff = abs(lhist-rhist);
        tg(y,x,i) = lrdiff' * tsim * lrdiff;
      end
      if i<norient,
        inc = pie(:,norient+i) - pie(:,i);
        lhist = lhist + inc;
        rhist = rhist - inc;
      end
    end
  end
end
fprintf(2,']\n');

for i = 1:norient,
  switch smooth,
   case 'gaussian',
    f = oeFilter([sigma .5],3,theta(i)+pi/2);
    tg(:,:,i) = applyFilter(f,tg(:,:,i));
   case 'savgol',
    a = fitparab(tg(:,:,i),sigma,sigma/4,theta(i));
    tg(:,:,i) = max(0,a);
  end
end

author: Cloud-CV

repo: object-proposals


text:
function run_REDaddon_sourcery(db, ops0)
ops = build_ops3(db, ops0);
% red channel addon to already processed data
mimgG = [];

if sum(ismember(ops.expts, ops.expred))>0
    mimgR = regRedChannelExpts(ops);
else
    if ops.nchannels_red==2
        % return green mean image from short red/green recording
        [mimgR, mimgG] = regRedGreenChannel(ops);
    else
        mimgR = regRedChannelOnly(ops);
    end
end

save(fullfile(ops.ResultsSavePath,'redchannel.mat'), 'mimgR','mimgG');

add_red_channel_sourcery(ops);

author: cortex-lab


text:
function [ fea, out ] = ex_laplace1( varargin )
%EX_LAPLACE1 2D Laplace equation example on a unit square.
%
%   [ FEA, OUT ] = EX_LAPLACE1( VARARGIN ) Laplace equation on a unit square
%   with exact solution 2*y/((1+x)^2+y^2). Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
%       hmax        scalar {1/10}          Max grid cell size
%       refsol      string {2*y/((1+x)^2+y^2)}   Reference solution
%       sfun        string {sflag1}        Shape function
%       iphys       scalar 0/{1}           Use physics mode to define problem    (=1)
%                                          or directly define fea.eqn/bdr fields (=0)
%       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 = { ...
  'igrid',    0; ...
  'hmax',     0.1; ...
  'refsol',   '2*y/((1+x)^2+y^2)'; ...
  'sfun',     'sflag1'; ...
  'iphys',    1; ...
  'icub',     2; ...
  'iplot',    1; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
gobj = gobj_rectangle();
fea.geom.objects = { gobj };

% Grid generation.
if ( opt.igrid==-1 )
  fea.grid = rectgrid(round(1/opt.hmax));
  fea.grid = quad2tri( fea.grid );
elseif ( opt.igrid==0 )
  fea.grid = rectgrid(round(1/opt.hmax));
else
  fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
end
n_bdr = max(fea.grid.b(3,:));           % Number of boundaries.

% Problem definition.
fea.sdim  = { 'x' 'y' };                % Coordinate names.
if ( opt.iphys==1 )

  fea = addphys(fea,@poisson);          % Add Poisson equation physics mode.
  fea.phys.poi.sfun = { opt.sfun };     % Set shape function.
  fea.phys.poi.eqn.coef{3,4} = { 0 };   % Set source term coefficient to zero.
  fea.phys.poi.bdr.coef{1,end} = repmat({opt.refsol},1,n_bdr);   % Set Dirichlet boundary coefficient to reference solution.
  fea = parsephys(fea);                 % Check and parse physics modes.

else

  fea.dvar  = { 'u' };                  % Dependent variable name.
  fea.sfun  = { opt.sfun  };            % Shape function.

  % Define equation system.
  fea.eqn.a.form = { [2 3;2 3] };       % First row indicates test function space   (2=x-derivative + 3=y-derivative),
                                        % second row indicates trial function space (2=x-derivative + 3=y-derivative).
  fea.eqn.a.coef = { 1 };               % Coefficient used in assembling stiffness matrix.

  fea.eqn.f.form = { 1 };               % Test function space to evaluate in right hand side (1=function values).
  fea.eqn.f.coef = { 0 };               % Coefficient used in right hand side.

  % Define boundary conditions.
  fea.bdr.d     = cell(1,n_bdr);
 [fea.bdr.d{:}] = deal(opt.refsol);     % Assign reference solution to all boundaries (Dirichlet).

  fea.bdr.n     = cell(1,n_bdr);        % No Neumann boundaries ('fea.bdr.n' empty).

end

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

% Postprocessing.
s_err = ['abs(',opt.refsol,'-u)'];
if ( opt.iplot>0 )
  figure
  subplot(2,1,1)
  postplot(fea,'surfexpr','u')
  title('Solution u')
  subplot(2,1,2)
  postplot(fea,'surfexpr',s_err)
  title('Error')
end

% Error checking.
if ( size(fea.grid.c,1)==4 )
  xi = [0;0];
else
  xi = [1/3;1/3;1/3];
end
err = evalexpr0(s_err,xi,1,1:size(fea.grid.c,2),[],fea);
ref = evalexpr0('u',xi,1,1:size(fea.grid.c,2),[],fea);
err = sqrt(sum(err.^2)/sum(ref.^2));

if( ~isempty(fid) )
  fprintf(fid,'\nL2 Error: %f\n',err)
  fprintf(fid,'\n\n')
end

out.err  = err;
out.pass = out.err<0.01;
if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function [sFile, ChannelMat] = in_fopen_intan(DataFile)
% IN_FOPEN_INTAN Open Intan RHS/RHD recordings
%
% DESCRIPTION:
%    Intan has 3 different ways of saving the raw file (Indicated by the AcqType variable throughout the code)
%       1. One .rhd/.rhs file for saving everything
%       2. A separate file for each channel (.dat), a .rhd/.rhs header, a separate file for each pin of the parallel port (.dat)
%       3. A separate file for each type of channel (one file for all data-files, one file for aux etc.)
%    This code supports 1 and 2.
%
%    The events are read from a parallel port, and all Digital outputs are saved in seperate files as well.
%    Information for the datatypes can be found at:
%        http://intantech.com/files/Intan_RHD2000_data_file_formats.pdf
%        http://intantech.com/files/Intan_RHS2000_data_file_formats.pdf

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Konstantinos Nasiotis 2018
%          Francois Tadel, 2019

%% ===== GET FILES =====
% Get base dataset folder
[hdr.BaseFolder, isItInfo, hdr.FileExt] = bst_fileparts(DataFile);
[filePath, fileComment] = bst_fileparts(hdr.BaseFolder);

% Check the type of file
if strcmp(isItInfo,'info')
    AcqType = 2; % One separate file per channel
else
    AcqType = 1; % One file to rule them all
end

% Get extra files
if (AcqType == 2)
    % Event files (.dat) % They are collected from a parallel port and saved as board-DIN-*.dat files
    EventFiles = dir(bst_fullfile(hdr.BaseFolder, '*DIN*.dat'));
    if isempty(EventFiles)
        disp(['BST> Warning: board-DIN-*.dat files not found in folder: ' 10 hdr.BaseFolder]);
        EventFiles = [];
    end
    % Recordings (amp*.dat)
    ampFiles = dir(bst_fullfile(hdr.BaseFolder, 'amp*.dat'));
    if isempty(ampFiles)
        error(['Could not find any amp*.dat recordings in folder: ' 10 hdr.BaseFolder]);
    end
end

%% ===== READ DATA HEADERS =====
hdr.chan_headers = {};
hdr.chan_files = {};

% Read the header
switch (hdr.FileExt)
    case '.rhd'
        newHeader = read_Intan_RHD2000_bst_2022(DataFile,1,1,1,100);
    case '.rhs'
        newHeader = read_Intan_RHS2000_bst_2018(DataFile,1,1,1,100);
end
newHeader.AcqType = AcqType; % This will be used later in in_fread_intan.m

% Check for the magic Number
if newHeader.magic_number ~= hex2dec('d69127ac') && newHeader.magic_number ~= hex2dec('c6912702')
    error('Magic Number Incorrect. The Intan header was not loaded properly');
end

% Check for some important fields
if ~isfield(newHeader, 'frequency_parameters') || ~isfield(newHeader, 'amplifier_channels') || ~isfield(newHeader, 'sample_rate')
    error('Missing fields in the file header of file');
end

% Check if all the channels' files are present
if AcqType==2
    if length(ampFiles) ~= newHeader.num_amplifier_channels
        error('Missing channel files. Check if the .dat files from all channels are present');
    end
    % Check if there are missing timestamps in the file
    % Check from the time.dat file how many samples exist, and compare with
    % every amp file that was collected
    fileinfo = dir(bst_fullfile(hdr.BaseFolder, 'time.dat'));
    num_samples_time = fileinfo.bytes/4; % int32 = 4 bytes
    for iChannel = 1:newHeader.num_amplifier_channels
        num_samples_channel = ampFiles(iChannel).bytes/2; % int32 = 4 bytes
        if (num_samples_time ~= num_samples_channel)
            error(['There are some missing blocks of recordings in file: ' ampFiles(iChannel).name]);
        end
    end
    % Read the time.dat file to extract first and last timestamp
    fid = fopen(bst_fullfile(hdr.BaseFolder, 'time.dat'), 'r');
    Time = fread(fid, num_samples_time, 'int32');
    fclose(fid);
    Time = Time / newHeader.frequency_parameters.amplifier_sample_rate; % sample rate from header file

    hdr.chan_files   = ampFiles;  % I ONLY LOAD THE AMP FILES HERE. THE AUXILIARY FILES HAVE DIFFERENT SAMPLING FREQUENCY. FIX THIS ON A LATER VERSION
    hdr.ChannelCount = length(hdr.chan_files);
    hdr.NumSamples   = num_samples_time;

else
    Time             = linspace(0,newHeader.nSamples/newHeader.sample_rate,newHeader.nSamples);
    num_samples_time = newHeader.sample_rate;
    hdr.ChannelCount = newHeader.num_amplifier_channels;
    hdr.NumSamples   = newHeader.nSamples;
end

% Extract information needed for opening the file
hdr.FirstTimeStamp        = Time(1);
hdr.LastTimeStamp         = Time(end);
hdr.SamplingFrequency     = newHeader.sample_rate;
% Save all file names
hdr.chan_headers = newHeader;
hdr.DataFile     = DataFile;

%% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder = 'l';
sFile.filename  = hdr.BaseFolder;
sFile.format    = 'EEG-INTAN';
sFile.device    = 'Intan';
sFile.header    = hdr;
sFile.comment   = fileComment;
sFile.condition = fileComment;
% Consider that the sampling rate of the file is the sampling rate of the first signal
sFile.prop.sfreq = hdr.SamplingFrequency;
sFile.prop.times = [0, hdr.NumSamples - 1] ./ sFile.prop.sfreq;
sFile.prop.nAvg  = 1;
% No info on bad channels
sFile.channelflag = ones(hdr.ChannelCount, 1);

%% ===== CREATE EMPTY CHANNEL FILE =====
ChannelMat = db_template('channelmat');
ChannelMat.Comment = 'Intan channels';
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, hdr.ChannelCount]);
% For each channel
for i = 1:hdr.ChannelCount
    if AcqType==2
        [fPath,fName,fExt] = bst_fileparts(hdr.chan_files(i).name);
        ChannelMat.Channel(i).Name = fName;
    else
        ChannelMat.Channel(i).Name = newHeader.amplifier_channels(i).custom_channel_name;
    end
    ChannelMat.Channel(i).Loc     = [0; 0; 0];
    ChannelMat.Channel(i).Type    = 'EEG';
    ChannelMat.Channel(i).Orient  = [];
    ChannelMat.Channel(i).Weight  = 1;
    ChannelMat.Channel(i).Comment = [];
end

%% ===== READ EVENTS =====
areThereEvents = 0;

if AcqType==2
    if ~isempty(EventFiles)
        % Read blocks of BINARY INPUT channels to create events
        parallel_input = false(num_samples_time, length(EventFiles));

        for iDIN = 1:length(EventFiles)
            fid = fopen(bst_fullfile(hdr.BaseFolder, EventFiles(iDIN).name), 'r');
            temp = fread(fid, num_samples_time, 'uint16');
            parallel_input(:,iDIN) = logical(temp);
            fclose(fid);
        end
        events_vector = bst_bi2de(parallel_input,'right-msb'); % Collapse the parallel ports to a vector that shows the events
        areThereEvents = 1;
    end
else
    if isfield(newHeader,'board_dig_in_data')
        events_vector = bst_bi2de(newHeader.board_dig_in_data,'right-msb'); % Collapse the parallel ports to a vector that shows the events
        areThereEvents = 1;
    end
end

if areThereEvents
    if bst_get('UseSigProcToolbox')
        [event_labels, event_samples] = findpeaks(events_vector);
    else
        [event_samples, event_labels] = peakseek(events_vector, min(events_vector));
        event_labels = event_labels';
        event_samples = event_samples';
    end

    % Create events list
    if ~isempty(event_labels)

        % Get list of events
        event_labels_unique = unique(event_labels); % These are still numbers

        % Initialize list of events
        events = repmat(db_template('event'), 1, length(event_labels_unique));
        % Format list
        for iEvt = 1:length(event_labels_unique)
            % Fill the event fields
            events(iEvt).label      = num2str(event_labels_unique(iEvt));
            events(iEvt).color      = rand(1,3);
            events(iEvt).reactTimes = [];
            events(iEvt).select     = 1;
            events(iEvt).times      = event_samples(event_labels == event_labels_unique(iEvt))' ./ sFile.prop.sfreq;
            events(iEvt).epochs     = ones(1, length(events(iEvt).times));    % Epoch: set as 1 for all the occurrences
            events(iEvt).channels   = [];
            events(iEvt).notes      = [];
        end
        % Import this list
        sFile = import_events(sFile, [], events);
    end
end
end

function de = bst_bi2de(bi, flg)
    if nargin < 2 || isempty(flg)
        flg = 'right-msb';
    end

    % Initialize array of powers of two
    [numElems, numDecimals] = size(bi);
    pows = zeros(numDecimals, 1);
    for iDec = 1:numDecimals
        pows(iDec) = 2 ^ (iDec - 1);
    end

    % Flip if required
    if strcmp(flg, 'left-msb')
        pows = flipud(pows);
    end

    de = bi * pows;
end

author: brainstorm-tools


text:
function pass = test_ivp_chebmatrix_syntax
% Test CHEBMATRIX syntax for IVPs, using full Brusselator IVP:

%% Setup
dom = [0 5];
tol = 1e-14;

%% With CHEBOP, multiple variable syntax for reference
N = chebop(@(t,u,v,w) [
    diff(u)-1-u^2*v + (w+1)*u;
    diff(v) - u*w+u^2*v;
    diff(w) - 1.5 + w*u], dom);
%% Test with left and right BCs:
N.lbc = [1;3;4];
uL = N\0;
N.lbc = [];
pass(1) = norm(feval(uL,dom(1)) - [1;3;4]) < tol;
%%
N.lbc = [];
N.rbc = [.7; 0.9; 1.2];
uR = N\0;
pass(2) = norm(feval(uR,dom(end)) - [.7; 0.9; 1.2]) < tol;
%% CHEBMATRIX syntax

% Define a new CHEBOP
M = chebop(@(t,u) [
    diff(u{1})-1-u{1}^2*u{2} + (u{3}+1)*u{1};
    diff(u{2}) - u{1}*u{3}+u{1}^2*u{2};
    diff(u{3}) - 1.5 + u{3}*u{1}], dom);

%% CHEBMATRIX left boundary conditions
M.lbc = @(u)[u{1}-1; u{2}-3; u{3}-4];
vL1 = M\0;
% Results should be identical to above
pass(3) = norm(vL1 - uL) == 0;

%% Vector LBC
M.lbc = [1; 3; 4];
vL2 = M\0;
pass(4) = norm(vL2 - uL) == 0;

%% CHEBMATRIX right boundary conditions
M.lbc = [];
M.rbc = @(u)[u{1}-.7; u{2}-0.9; u{3}-1.2];
vR1 = M\0;
% Results should be identical to above
pass(5) = norm(vR1 - uR) == 0;

%% Vector RBC
M.rbc = [.7; .9; 1.2];
vR2 = M\0;
pass(6) = norm(vR2 - uR) == 0;

end


text:
function scale = p15_scale ( neqn )

%*****************************************************************************80
%
%% P15_SCALE returns scale factors for problem p15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end


text:
function [indic, xlevels] = condf2indic(X, varargin)
% Create N x k indicator matrix of ones/zeros from N x 1 list of numeric
% values (X)
%
% :Usage:
% ::
%
%     [indic, xlevels] = condf2indic(X, ['integers'])
%
% :Inputs:
%
%   **'integers':**
%       Treat X as integer vector. Check that all are integers, and return
%       one column per integer in order, excluding zeros. Insert NaNs in
%       columns whose corresponding integers are missing.  Returns an n x k
%       matrix indic, where n is size(X, 1) and k is max(X).
%
%       'integers' may be followed by another variable specifying the
%       number of columns to use. this is useful because integers at the
%       end of a list in X may be missing; this will append columns for
%       those.
%
% :Outputs:
%
%   **indic:**
%        is returned as single precision type, because it can then be
%        used as a design matrix in a GLM
%
%   **xlevels:**
%        are the values of X corresponding to columns of indic
%
% :Examples:
% condf = [1 1 1 2 2 2 4 4 4 6 6 1 1 1 8 8]';
% [indic, xlevels] = condf2indic(condf); indic, xlevels
% [indic, xlevels] = condf2indic(condf, 'integers'); indic, xlevels
% [indic, xlevels] = condf2indic(condf, 'integers', 10); indic, xlevels

% ..
%    tor wager, nov 2007. Edited Feb 2018 to expand functionality to
%    'integers' case, for use with atlas object.
%    Aug 2018: fixed some inconsistencies in 'integers' option.
% ..

% PRELIM CALCULATIONS

X = double(X); % for objects/special types

[xlevels, ~, indx] = unique(X);

max_integers = length(xlevels);

% if X contains only zero values, return
if ~any(X)
    indic=[];
    return;
end

% INPUTS

dointegers = false;

if any(strcmp(varargin, 'integers'))

    dointegers = true;

    wh_input = strcmp(varargin, 'integers');

    if length(varargin) > find(wh_input)      % we have entered max integers

        max_integers = varargin{wh_input + 1};

    else

        max_integers = max(xlevels);

    end

end

% INITIALIZE - cols equal to num unique entries, or all integers

indic = single(false(length(indx), max_integers));

if dointegers

    % Check integers
    u = unique(X);
    if ~all(u == round(u))
        warning('Some condition function values are not integers.');
    end

    % Remove zero-valued integer
    wh = find(xlevels == 0);
    if ~isempty(wh)
        xlevels(wh) = [];
        indic(:, wh) = [];
        badIndx = indx(find(indx == wh, 1, 'first'));
        indx(indx == wh) = 0;
        indx(indx > badIndx) = indx(indx > badIndx) - 1;
    end

    for i = 1:length(xlevels)
        indic(indx == i, xlevels(i)) = 1;
    end

    xlevels = (1:max_integers)';
%
%     % if codes are 1...n, xlevels is index values...but if there are
%     % missing integers, need to insert into indic
%
%     is_missing = true(max(X), 1);
%     is_missing(xlevels) = false;
%
%     % Insert NaNs for any missing index values.
%
%     indic = naninsert(is_missing, indic')';

else
    % no integers, do not preserve empty columns

    for i = 1:length(xlevels)
        indic(indx == i, i) = 1;
    end

end % integers

end % main function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% BV_READ  Reads BrainVision data over existing TCP socket.
%     DATA = BV_READ(H)
%     In:
%     h : handle to a BrainVision connection (previously opened via bv_open)
%     Out:
%     data : a block of new data, [#Channels x #Samples]

% Author: Hal Greenwald, The MITRE Corporation, 29-NOV-2011
% This software incorporates portions of the BrainVision Recorder RDA Client.  Used with permission of Brain Products GmbH.

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

function data = bv_read(h)
data = [];
if (~h.initialized)
    return;
end

header_size = 24;
try
    % check for existing data in socket buffer
    tryheader = pnet(h.handle, 'read', header_size, 'byte', 'network', 'view', 'noblock');
    while ~isempty(tryheader)

        % Read header of RDA message
        hdr = ReadHeader(h.handle);

        % Perform some action depending on the type of the data package
        switch hdr.type
            case 1       % Start message
                % Reset block counter to check overflows
                h.lastBlock = -1;
                pnet(h.handle, 'read', hdr.size - header_size);

            case 4       % 32Bit Data block
                % Read data and markers from message
                [datahdr, eegdata, markers] = ReadDataMessage(h.handle, hdr, h.channelCount);
                data = cat(2, data, eegdata);

                % check tcpip buffer overflow
                if h.lastBlock ~= -1 && datahdr.block > h.lastBlock + 1
                    disp(['******* Overflow with ' int2str(datahdr.block - h.lastBlock) ' blocks ******']);
                end
                h.lastBlock = datahdr.block;

                % Edit the following lines to process marker information
%                 for m = 1:datahdr.markerCount
%                     disp(markers(m));
%                 end

            case 3       % Stop message
                disp('Received STOP signal from Recorder');
                pnet(h.handle, 'read', hdr.size - header_size);
                bv_close(h);
                h.initialized = false;
                data = [];

            otherwise    % ignore all unknown types, but read the package from buffer
                pnet(h.handle, 'read', hdr.size - header_size);
        end
        if h.initialized
            tryheader = pnet(h.handle, 'read', header_size, 'byte', 'network', 'view', 'noblock');
        else
            tryheader = [];
        end
    end
catch er
    disp(er.message);
end

if ~isempty(data)
    % scale the data to uV (consistent with the pop_loadbv plugin)
    data = bsxfun(@times,data,h.resolutions');
end

%% ***********************************************************************
% Read the message header
function hdr = ReadHeader(con)
% con    tcpip connection object

% define a struct for the header
hdr = struct('uid',[],'size',[],'type',[]);

% read id, size and type of the message
% swapbytes is important for correct byte order of MATLAB variables
hdr.uid = pnet(con,'read', 16);
hdr.size = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));
hdr.type = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));

%% ***********************************************************************
% Read a data message
function [datahdr, data, markers] = ReadDataMessage(con, hdr, channelCount)
% con       tcpip connection object
% hdr       message header
% datahdr   data header with information on datalength and number of markers
% data      data as one dimensional arry
% markers   markers as array of marker structs

% Define data header struct and read data header
datahdr = struct('block',[],'points',[],'markerCount',[]);

datahdr.block = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));
datahdr.points = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));
datahdr.markerCount = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));

% Read data in float format
data = swapbytes(pnet(con,'read', channelCount * datahdr.points, 'single', 'network'));
data = reshape(data,channelCount,datahdr.points);

% Define markers struct and read markers
markers = struct('size',[],'position',[],'points',[],'channel',[],'type',[],'description',[]);
for m = 1:datahdr.markerCount
    marker = struct('size',[],'position',[],'points',[],'channel',[],'type',[],'description',[]);

    % Read integer information of markers
    marker.size = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));
    marker.position = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));
    marker.points = swapbytes(pnet(con,'read', 1, 'uint32', 'network'));
    marker.channel = swapbytes(pnet(con,'read', 1, 'int32', 'network'));

    % type and description of markers are zero-terminated char arrays
    % of unknown length
    c = pnet(con,'read', 1);
    while c ~= 0
        marker.type = [marker.type c];
        c = pnet(con,'read', 1);
    end

    c = pnet(con,'read', 1);
    while c ~= 0
        marker.description = [marker.description c];
        c = pnet(con,'read', 1);
    end

    % Add marker to array
    markers(m) = marker;
end

author: PatternRecognition

repo: OpenBMI


text:
function [beamGeometryInitS,tags] = dicomrt_d2c_beam(beamGeometryInitS,indexS,study,tags)
% dicomrt_d2c_beam(beamGeometryInitS,indexS,plan,tags)
%
% Convert DICOM beam data in CERR format.
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check input data
[study,type,dummylabel,PatientPosition]=dicomrt_checkinput(study);

if strcmpi('RTPLAN',type)~=1
    error('dicomrt_d2c_scan: input data does not have the right format. Exit now!');
else
    type='BEAM GEOMETRY';
end

% Get DICOM-RT toolbox dataset info
study_pointer=study{1,1};
temp_study_header=study_pointer{1};

% Get parameters
beams=fieldnames(temp_study_header.BeamSequence);
nbeams=size(beams,1);

for i=1:nbeams; % loop through the number of slices
    beamGeometryInitS(i).imageNumber                = tags.nimages +i;
    beamGeometryInitS(i).imageType                  = type;
    beamGeometryInitS(i).caseNumber                 = 1;
    beamGeometryInitS(i).patientName                = [temp_study_header.PatientName.FamilyName ,' ', ...
                temp_study_header.PatientName.GivenName];
    beamGeometryInitS(i).beamNumber                 = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'BeamNumber');
    % RadiationType
    RadiationType=getfield(temp_study_header,'BeamSequence',char(beams(i)),'RadiationType');
    if strcmpi(RadiationType,'PHOTON')==1
        beamGeometryInitS(i).beamModality               = 'X-RAY';
    else % 'OTHER' is not a possibility
        beamGeometryInitS(i).beamModality               = RadiationType;
    end
    beamGeometryInitS(i).beamEnergyMeV              = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'ControlPointSequence','Item_1','NominalBeamEnergy');
    beamGeometryInitS(i).beamDescription            = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'BeamName');
    beamGeometryInitS(i).RxDosePerTxGy              = getfield(temp_study_header,'FractionGroupSequence','Item_1',...
        'ReferencedBeamSequence',char(beams(i)),'BeamDose');
    beamGeometryInitS(i).numberOfTx                 = getfield(temp_study_header,'FractionGroupSequence','Item_1',...
        'NumberOfFractionsPlanned');

    beamGeometryInitS(i).fractionGroupID            = temp_study_header.RTPlanDescription;
    beamGeometryInitS(i).beamType                   = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'BeamType');
    beamGeometryInitS(i).planIDOfOrigin             = temp_study_header.StudyID;
    % collimatorType
    xjaw = getfield(temp_study_header,'BeamSequence',char(beams(i)),'ControlPointSequence','Item_1',...
        'BeamLimitingDevicePositionSequence','Item_1','LeafJawPositions');
    yjaw =  getfield(temp_study_header,'BeamSequence',char(beams(i)),'ControlPointSequence','Item_1',...
        'BeamLimitingDevicePositionSequence','Item_2','LeafJawPositions');
    if -xjaw(1)==xjaw(2) & -yjaw(1)==yjaw(2)
        beamGeometryInitS(i).collimatorType             = 'SYMMETRIC';
    elseif -xjaw(1)==xjaw(2) & -yjaw(1)~=yjaw(2)
        beamGeometryInitS(i).collimatorType             = 'ASYMMETRIC_Y';
    elseif -xjaw(1)==xjaw(2) & -yjaw(1)~=yjaw(2)
        beamGeometryInitS(i).collimatorType             = 'ASYMMETRIC_X';
    else
        beamGeometryInitS(i).collimatorType             = 'ASYMMETRIC';
    end
    % apertureType
    % set default
    nblocks = 0;
    ncomp   = 0;
    mlcmod  = 0;
    nwedges = 0;
    % retrieve info
    nblocks = getfield(temp_study_header,'BeamSequence',char(beams(i)),'NumberOfBlocks');
    ncomp   = getfield(temp_study_header,'BeamSequence',char(beams(i)),'NumberOfCompensators');
    mlcmod  = 0;
    nwedges = getfield(temp_study_header,'BeamSequence',char(beams(i)),'NumberOfWedges');
    % check
    if nblocks ~= 0
        beamGeometryInitS(i).apertureType               = 'BLOCKS';
    elseif ncomp ~= 0
        beamGeometryInitS(i).apertureType               = 'TRANSMISSION_MAP';
    elseif mlcmod ~= 0
        beamGeometryInitS(i).apertureType               = 'MLC_XY';
    else
        beamGeometryInitS(i).apertureType               = 'COLLIMATOR';
    end
    beamGeometryInitS(i).apertureDescription        = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'BeamName');
    beamGeometryInitS(i).collimatorAngle            = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'ControlPointSequence','Item_1','BeamLimitingDeviceAngle');
    beamGeometryInitS(i).gantryAngle                = 360-getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'ControlPointSequence','Item_1','GantryAngle');
    beamGeometryInitS(i).couchAngle                 = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'ControlPointSequence','Item_1','TableTopEccentricAngle');
    beamGeometryInitS(i).headInOut                  = tags.hio;
    beamGeometryInitS(i).nominalIsocenterDistance   = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'SourceAxisDistance').*0.1;
    beamGeometryInitS(i).numberRepresentation       = 'CHARACTER';
    % Optional
    beamGeometryInitS(i).apertureID                 = '';
    if nwedges ~= 0
        beamGeometryInitS(i).wedgeAngle                 = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
            'WedgeSequence','Item_1','WedgeAngle');
        beamGeometryInitS(i).wedgeRotationAngle         = 360-getfield(temp_study_header,'BeamSequence',char(beams(i)),...
            'WedgeSequence','Item_1','WedgeOrientation');
    else
        beamGeometryInitS(i).wedgeAngle                 = '';
        beamGeometryInitS(i).wedgeRotationAngle         = '';
    end
    beamGeometryInitS(i).arcAngle                   = '';
    beamGeometryInitS(i).machineID                  = getfield(temp_study_header,'BeamSequence',char(beams(i)),...
        'TreatmentMachineName');
    beamGeometryInitS(i).beamWeight                 = '';
    beamGeometryInitS(i).weightUnits                = '';
    beamGeometryInitS(i).compensator                = '';
    beamGeometryInitS(i).compensatorFormat          = '';
    beamGeometryInitS(i).file                       = '';
end

% Update nimages
tags.nimages = tags.nimages + nbeams;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% [B3_] = tensor_fibers_tube(A)
%
function [B3_] = tensor_fibers_tube(A)
  [n1, n2, n3] = size(A);
  for i = 1:n1
    for j = 1:n2
      B3_{i,j} = A(i,j,:);
    end
  end
end

author: andrewssobral


text:
function [R,L_lb,L_cpd] = rankest(T,options)
%RANKEST Estimate rank.
%   rankest(T) plots an L-curve of number of rank-one terms in a canonical
%   polyadic decomposition. The x-axis corresponds to the number of
%   rank-one terms, and the y-axis corresponds to the relative error of the
%   CPD in that many rank-one terms. Additionally, the corner R of the
%   resulting L-curve is estimated.
%
%   R = rankest(T) does not plot anything and instead returns the number of
%   rank-one terms R corresponding to the corner of the L-curve.
%
%   [R,L_lb,L_cpd] = rankest(T) also returns the L-curve ranks L_lb(:,1)
%   and L_cpd(:,1) and the corresponding lower bound on the truncation
%   error L_lb(:,2) and relative error of the CPD approximation L_cpd(:,2),
%   respectively.
%
%   rankest(T,options) may be used to set the following options:
%
%      options.MaxR =           - Maximum number of rank-one terms to try.
%      numel(T)/max(size_tens)
%      options.MinR = 1         - Minimum number of rank-one terms to try.
%      options.MinRelErr = 1e-2 - Determines an upper threshold for the
%                                 number of rank-one terms R to try.
%      options.Solver = @cpd    - The solver used to compute the CPD for
%                                 each number of rank-one terms R. Called
%                                 as options.Solver(T,R, ...
%                                 options.SolverOptions), where
%                                 options.SolverOptions is an options
%                                 structure passed to the solver.
%      options.XMultiplier = 1  - The importance of the number of rank-one
%                                 terms in determining the L-curve corner,
%                                 relative to the importance of the
%                                 relative error of the approximation.
%
%   See also mlrankest.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] J.L. Castellanos, S. Gomez, V. Guerra, "The triangle method for
%       finding the corner of the L-curve," Applied Numerical Mathematics,
%       Vol. 43, No. 4, 2002, pp. 359-373.

% Incomplete/sparse tensors not supported yet.
T = fmt(T,true);
if isstruct(T), size_tens = T.size;
else size_tens = size(T); end

% Check the options structure.
if nargin < 2, options = struct; end
if ~isfield(options,'MaxR')
    options.MaxR = prod(size_tens)/max(size_tens);
end
if ~isfield(options,'MinR'), options.MinR = 1; end
if ~isfield(options,'MinRelErr'), options.MinRelErr = 1e-2; end
if ~isfield(options,'Solver'), options.Solver = @cpd; end
if ~isfield(options,'SolverOptions'), options.SolverOptions = struct; end
if ~isfield(options.SolverOptions,'Compression')
    options.SolverOptions.Compression = false;
end
if ~isfield(options.SolverOptions,'Display')
    options.SolverOptions.Display = false;
end
if ~isfield(options,'XMultiplier'), options.XMultiplier = 1; end

% Compute lower bound on truncation error.
frobT = frob(T);
if ~isstruct(T)
    [~,~,sv] = mlsvd(T);
    Rmax = max(cellfun(@length,sv))-1;
    lb = max(cell2mat(cellfun( ...
             @(s)[sqrt(flipud(cumsum(flipud(s(:).^2)))).'/frobT ...
             zeros(1,Rmax+1-length(s))],sv(:),'UniformOutput',false)));
    lb = lb(2:end);
else
    Rmax = 0;
    lb = [];
end

% Determine the range of rank-one terms to test.
if ~isempty(lb)
    R = find(lb(options.MinR:end)<=options.MinRelErr,1,'first')+ ...
        options.MinR-1;
    if isempty(R), R = max(options.MinR,Rmax+1); end
    R = R:options.MaxR;
else
    R = 1:options.MaxR;
end

% Compute the relative error for each number of rank-one terms R(r).
relerr = zeros(1,length(R));
for r = 1:length(R)

    % Compute the CPD in R rank-one terms.
    [U,output] = options.Solver(T,R(r),options.SolverOptions);
    if isfield(output,'Refinement') && ...
       isfield(output.Refinement,'fval')
        relerr(r) = sqrt(2*output.Refinement.fval(end))/frobT;
    elseif isfield(output,'Algorithm') && ...
           isfield(output.Algorithm,'fval')
        relerr(r) = sqrt(2*output.Algorithm.fval(end))/frobT;
    elseif isfield(output,'fval')
        relerr(r) = sqrt(2*output.fval(end))/frobT;
    else
        relerr(r) = frob(cpdres(T,U))/frobT;
    end

    % Stop is options.MinRelErr has been reached.
    if relerr(r) <= options.MinRelErr
        R = R(1:r);
        relerr = relerr(1:r);
        break;
    end

    % Update the plot.
    if r > 1 && r < length(R) && nargout == 0, Lcurve(true); end

end

% Compute L-curve corner using the triangle method [1].
logx = options.XMultiplier*R(:);
logy = log10(relerr(:));
if R(1) ~= 1 && ~isempty(lb)
    logx = [options.XMultiplier;logx];
    logy = [log10(lb(1));logy];
end
ab = cat(3,bsxfun(@minus,logx,logx.'),bsxfun(@minus,logy,logy.'));
ac = cat(3,logx(end)-logx.',logy(end)-logy.');
area = bsxfun(@times,ab(:,:,1),ac(:,:,2)) - ...
       bsxfun(@times,ab(:,:,2),ac(:,:,1));
cosa = bsxfun(@times,ab(:,:,1),ac(:,:,1)) + ...
       bsxfun(@times,ab(:,:,2),ac(:,:,2));
cosa = bsxfun(@rdivide,cosa./sqrt(ab(:,:,1).^2+ab(:,:,2).^2), ...
                       sqrt(ac(:,:,1).^2+ac(:,:,2).^2));
cosa(area >= 0 | tril(true(size(cosa))) | cosa <= cos(7*pi/8)) = -1;
[a,opt] = max(max(cosa(:,2:end)));
if isempty(a), opt = 1; end
if a == -1, opt = size(cosa,2)-1; end

% Display output.
L_lb = [(1:Rmax).' lb(:)];
L_cpd = [R(:) relerr(:)];
if nargout == 0, Lcurve(false); end
R = R(opt);

function Lcurve(update)
    if R(r) <= 1, return; end
    style = {'Marker','+','MarkerSize',2.5};
    if ~update
        semilogy(R(opt),relerr(opt),'rs','LineStyle','none'); hold on;
    end
    semilogy(R(1:length(relerr)),relerr,'r',style{:}); hold on;
    semilogy(1:length(lb),lb,style{:});
    semilogy([1 R(r)],options.MinRelErr*[1 1],'k:');
    text(1,options.MinRelErr,'MinRelErr','VerticalAlignment','Top');
    hold off;
    ylim([min(options.MinRelErr/(10^0.5),min(relerr(1:r))) lb(1)]);
    xlim([1 R(r)]);
    ylabel('frob(cpdres(T,U))/frob(T)'); xlabel('R');
    xt = get(gca,'XTick'); set(gca,'XTick',xt(mod(xt,1) == 0));
    if update
        if ~isempty(lb)
            legend(['CPD error (trying R = ' int2str(R(r+1)) '...)'], ...
                'Lower bound on error','Location','NE');
        else
            legend(['CPD error (trying R = ' int2str(R(r+1)) '...)'], ...
                'Location','NE');
        end
    else
        if ~isempty(lb)
            legend(['L-curve corner at R = ' int2str(R(opt))], ...
                'CPD error','Lower bound on error','Location','NE');
        else
            legend(['L-curve corner at R = ' int2str(R(opt))], ...
                'CPD error','Location','NE');
        end
    end
    set(datacursormode(gcf),'UpdateFcn',@datacursor);
    drawnow;
end

function txt = datacursor(~,event_obj)
    pos = get(event_obj,'Position');
    txt = {['R: ' int2str(pos(1))], ...
           ['relative error: ' num2str(pos(2))]};
end

end


text:
function pass = test_linearize_init_fails(~)
%TEST_LINEARIZE_INIT_FAILS   Test that we get sensible error messages when we
%are dealing with an initial guess that makes the chebop fail to evaluate.

%% First order problem which is OK, as it's solved as an IVP
N = chebop(0, 1);
N.op = @(u) diff(u) - sqrt(u);
N.lbc = 1;
u = N\1;
pass(1) = norm(N(u)-1) < 1e-10;

%% First order problem which is OK, not OK when try to solve as a BVP
N = chebop(0, 1);
N.op = @(u) diff(u) - sqrt(u);
N.bc = @(u) u(0) - 1;
try
    u = N\1;
catch ME
    pass(2) = strcmp(ME.identifier, ...
        'CHEBFUN:CHEBOP:linearize:invalidInitialGuess');
end

%% Nor when solvebvp is called directly
N = chebop(0, 1);
N.op = @(u) diff(u) - sqrt(u);
N.lbc = 1;
try
    u = solvebvp(N, 1);
catch ME
    pass(3) = strcmp(ME.identifier, ...
        'CHEBFUN:CHEBOP:linearize:invalidInitialGuess');
end

%% Second order problem
N = chebop(0, 1);
N.op = @(u) diff(u,2) - sqrt(u);
N.bc = 1;
try
    u = N\1;
catch ME
    pass(4) = strcmp(ME.identifier, ...
        'CHEBFUN:CHEBOP:linearize:invalidInitialGuess');
end

%% Second order problem, different function causing issues
N = chebop(0, 1);
N.op = @(u) diff(u,2) - 1./u;
N.bc = 1;
try
    u = N\1;
catch ME
    pass(5) = strcmp(ME.identifier, ...
        'CHEBFUN:CHEBOP:linearize:invalidInitialGuess');
end

%% Coupled system, OK when solved as an IVP
N = chebop(0, 1);
N.op = @(x, u, v) [diff(u,2) - sqrt(v); diff(v,2) + 1./u];
N.lbc = @(u,v) [u-1; diff(u)-.1; v-2; diff(v) + .2];
uv = N\[1; 2];
pass(6) = norm(N(uv) - [1;2]) < 1e-8;

%% Coupled system, OK when solved as an IVP
N = chebop(0, 1);
N.op = @(x, u, v) [diff(u,2) - sqrt(v); diff(v,2) + 1./u];
N.bc = @(x,u,v) [u(0) - 1; u(1) + 1; v(0) - 2; v(0) - 3];
try
    uv = N\1;
catch ME
    pass(7) = strcmp(ME.identifier, ...
        'CHEBFUN:CHEBOP:linearize:invalidInitialGuess');
end


text:
function [ V ] = GetValues( Q , knn , p )
%GETVALUE Summary of this function goes here
%   Detailed explanation goes here

V = Q(knn,:)' * p ;

end

author: epfl-lasa


text:
function [nodes, elements, nodedata, elementdata] = read_gmsh_binary(filename)

% READ_GMSH_BINARY reads a gmsh .msh binary file. Current support is only
% for version 2. There are some ASCII-readers floating around on the net,
% but they do not seem to work with the primary use case in FieldTrip (and
% the test data that I have available), which is SimNibs generated data.
%
% See also MESH_LOAD_GMSH4

% Copyright (C) 2023, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen(filename, 'rb');
txt = fgetl(fid);

if startsWith(txt, '$MeshFormat')
  % the file starts, as expected with the file header
  txt    = fgetl(fid);
  format = sscanf(txt, '%f %d %d');
  version  = format(1);
  isbinary = format(2)==1;
  dtype    = format(3);

  % read until the end of this chunk
  while 1
    % in a binary file, there's an endianness check byte, which is
    % currently not handled here.
    txt = fgetl(fid);
    if startsWith(txt, '$EndMeshFormat')
      break;
    end
  end
else
  ft_error('File %s cannot be read, header information is missing', filename);
end

if floor(version)~=2 || ~isbinary
  if ~isbinary
    str = 'not';
  else
    str = '';
  end
  ft_error('Currently only binary version 2 files can be read, the version is: %1.1f, and the file is %s binary', version, str);
end
if dtype~=8
  ft_error('Currently only 8-byte precision data is supported');
end

nodes       = [];
elements    = [];
nodedata    = [];
elementdata = [];

% read the contents of the file
while 1
  txt = fgetl(fid);
  if txt==-1
    % end of file
    fclose(fid);
    break;
  end
  % Nodes and Elements are expected to be minimally present, needed to
  % describe the mesh topology
  if startsWith(txt, '$Nodes')
    nodes = getnodes(fid);
  elseif startsWith(txt, '$Elements')
    elements = getelements(fid);
  elseif startsWith(txt, '$PhysicalNames')
    keyboard
  elseif startsWith(txt, '$Entities')
    keyboard
  elseif startsWith(txt, '$PartitionedEntities')
    keyboard
  elseif startsWith(txt, '$Periodic')
    keyboard
  elseif startsWith(txt, '$GhostElements')
    keyboard
  elseif startsWith(txt, '$Parametrizations')
    keyboard
  elseif startsWith(txt, '$NodeData')
    keyboard
  elseif startsWith(txt, '$ElementData')
    elementdata = getelementdata(fid);
  elseif startsWith(txt, '$ElementNodeData')
    keyboard
  elseif startsWith(txt, '$InterpolationScheme')
    keyboard
  else
    keyboard

  end
end

%%%%%%%%%%%%%%
% subfunctions

% $Nodes
function [nodes] = getnodes(fid)

% the chunk starts with an ascii-line with the number of elements
%
% the format of the nodes is 4-byte indx followed by 3x a float for the
% coordinates

txt = fgetl(fid);
N   = sscanf(txt, '%f');
ptr = ftell(fid);

fprintf('Reading position information for %d nodes\n', N);
nodes.indx = fread(fid, N, 'uint32', 24);
fseek(fid, ptr+4, 'bof');
nodes.nodes = fread(fid, 3*N, '3*double', 4);
nodes.nodes = reshape(nodes.nodes, 3, [])';
fseek(fid, ptr+28*N, 'bof');
txt = fgetl(fid);
assert(isequal(txt, '$EndNodes'), 'Reading of the nodes unexpectedly failed');

% $Elements
function [elements] = getelements(fid)

% the chunk starts with an ascii-line mentioning the overall number of elements, then 3 integers
% per sub-chunk of elements, that indicate the type, the number of elements
% in the sub-chunk, and the number of tags

txt = fgetl(fid);
N   = sscanf(txt, '%f');

fprintf('Reading element information for %d elements\n', N);
Ntotal = 0;
while N>Ntotal

  hdr  = fread(fid, 3, 'uint32');
  type = hdr(1);
  Nsub = hdr(2);
  ntag = hdr(3);

  [nnode, type_str] = type2nnode(type);

  fprintf('Reading %d %s with %d tags\n', Nsub, type_str, ntag);
  tmp = fread(fid, Nsub.*(nnode+ntag+1), '*uint32');
  tmp = reshape(tmp, nnode+ntag+1, [])';

  elements.(type_str)                     = tmp(:, (ntag+2):end);
  elements.(sprintf('%s_indx', type_str)) = tmp(:, 1);
  elements.(sprintf('%s_tag',  type_str)) = tmp(:, 1 + (1:ntag));
  Ntotal = Ntotal + Nsub;
end
txt = fgetl(fid);
assert(isequal(txt, '$EndElements'), 'Reading of the elements unexpectedly failed');

% $ElementData
function [elementdata] = getelementdata(fid)

% stringTags
txt = fgetl(fid);
N   = sscanf(txt, '%d');
for k = 1:N
  stringTag{k,1} = fgetl(fid);
end

% realTags
txt = fgetl(fid);
N   = sscanf(txt, '%d');
for k = 1:N
  txt = fgetl(fid);
  realTag(k,1) = sscanf(txt, '%d');
end

% integerTags
txt = fgetl(fid);
N   = sscanf(txt, '%d'); % N is kind of expected to be 4
for k = 1:N
  txt = fgetl(fid);
  integerTag(k,1) = sscanf(txt, '%d');
end
nfc = integerTag(2); % number of field components
N   = integerTag(3); % number of elements

fprintf('Reading elementdata for %d elements\n', N);
ptr = ftell(fid);
elementdata.indx = fread(fid, N, 'uint32', nfc*8);
fseek(fid, ptr+4, 'bof');
elementdata.value = fread(fid, nfc*N, sprintf('%d*double', nfc), 4);
elementdata.value = reshape(elementdata.value, nfc, [])';
fseek(fid, ptr+(4+nfc*8)*N, 'bof');

txt = fgetl(fid);
assert(isequal(txt, '$EndElementData'), 'Reading of the elementdata unexpectedly failed');

% helper function to get the number of nodes per element + a human interpretable name
function [nnode, type_str] = type2nnode(type)

switch type
  case 1
    % 2-node line
    nnode = 2; type_str = 'lines';
  case 2
    % 3-node triangle
    nnode = 3; type_str = 'triangles';
  case 3
    % 4-node quadrangle
    nnode = 4; type_str = 'quadrangles';
  case 4
    % 4-node tetrahedron
    nnode = 4; type_str = 'tetrahedra';
  case 5
    % 8-node hexahedron
    nnode = 8; type_str = 'hexahedra';
  case 6
    % 6-node prism
    nnode = 6; type_str = 'prisms';
  case 7
    % 5-node pyramid
    nnode = 5; type_str = 'pyramids';
  case 8
    % 3-node seconde order line (2 nodes associated with the vertices and 1 with the edge)
    nnode = 3; type_str = 'line';
  case 9
    % 6-node second order triangle (3 nodes associated with the vertices and 3 with the edges)
    nnode = 6; type_str = 'triangles';
  case 10
    % 9-node second order quadrangle (4 nodes associated with the vertices, 4 with the edges and 1 with the face)
    nnode = 9; type_str = 'quadrangles';
  case 11
    % 10-node second order tetrahedron (4 nodes associated with the vertices and 6 with the edges)
    nnode = 10; type_str = 'tetrahedra';
  case 12
    % 27-node second order hexahedron (8 nodes associated with the vertices, 12 with the edges,  6 with the faces and 1 with the volume)
    nnode = 27; type_str = 'hexahedra';
  case 13
    % 18-node second order prism (6 nodes associated with the vertices, 9 with the edges and 3 with the quadrangular faces)
    nnode = 18; type_str = 'prisms';
  case 14
    % 14-node seconde order pyramid (5 nodes associated with the vertices, 8 with the edges and 1 with the quadrangular face)
    nnode = 14; type_str = 'pyramids';
  case 15
    % 1-node point
    nnode = 1; type_str = 'points';
  otherwise
    %% FIXME add types 16-31, and 92 93
    keyboard
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function f_star = processOutputQP_oneFoot(analyticalSolution,primalSolution,QPStatus,feetContactStatus,Config)

    % PROCESSOUTPUTQP_ONEFOOT evaluates the output of the WBToolbox QP block.
    %                         In case the QP block exited with an error, a
    %                        "default", user defined solution to the QP problem
    %                         is provided instead of the one coming from the QP block.
    %
    % FORMAT: f_star = processOutputQP_oneFoot(analyticalSolution,primalSolution,QPStatus,feetContactStatus,Config)
    %
    % INPUT:  - analyticalSolution = the alternative user defined solution to the QP
    %                                problem to be used when the QP block fails;
    %         - primalSolution = the solution to the QP problem provided by the
    %                            WBToolbox QP block;
    %         - QPStatus = the status check returned by the QP block.
    %
    %         - feetContactStatus = [2 x 1] a vector describing the feet contact
    %                               status. Format: [leftFoot; rightFoot];
    %         - Config = a structure containing the robot-related configuration
    %                    parameters;
    %
    % OUTPUT: - f_star = the final solution to the QP problem that will be used in
    %                    the controller.
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

    CONTACT_THRESHOLD   = 0.1;
    QP_STATUS_THRESHOLD = 0.01;

    % if the robot is balancing on one foot, extend the QP solution to
    % have the correct dimension of f_star for matrix multiplication
    if feetContactStatus(1) > (1 - CONTACT_THRESHOLD)

        % left foot balancing
        updated_primalSolution     = [primalSolution; zeros(6,1)];
        updated_analyticalSolution = [analyticalSolution; zeros(6,1)];

    else
        % right foot balancing
        updated_primalSolution     = [zeros(6,1); primalSolution];
        updated_analyticalSolution = [zeros(6,1); analyticalSolution];
    end

    if Config.USE_QP_SOLVER && abs(QPStatus)< QP_STATUS_THRESHOLD

        f_star = updated_primalSolution;
    else
        f_star = updated_analyticalSolution;
    end
end

author: robotology

repo: whole-body-controllers


text:
function bpa = tapas_bayesian_parameter_average(varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function calculates the Bayesian parameter average for the individual estimates handed  to
% it.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% USAGE:
%     est1 = tapas_fitModel(responses1, inputs);
%     est2 = tapas_fitModel(responses2, inputs);
%     ...
%     estn = tapas_fitModel(responsesn, inputs);
%
%     bpa = tapas_bayesian_parameter_average(est1, est2,..., estn);
%
% INPUT ARGUMENTS:
%     varargin           Estimate structures generated by tapas_fitModel(...). Note that all estimates
%                        must have been made under the same priors.
%
% OUTPUT:
%     bpa.u              Input to agent (i.e., the inputs array from the arguments)
%     bpa.c_prc          Configuration settings for your chosen perceptual model
%                        (see the configuration file of that model for details)
%     bpa.c_obs          Configuration settings for your chosen observation model
%                        (see the configuration file of that model for details)
%     bpa.optim          A place for the optimization algorithm to dump infos of interest to it
%     bpa.p_prc          Bayesian average of estimates of perceptual parameters
%                        (see the configuration file of your chosen perceptual model for details)
%     bpa.p_obs          Bayesian average of estimates of observation parameters
%                        (see the configuration file of your chosen observation model for details)
%     bpa.traj:          Trajectories of the environmental states tracked by the perceptual model
%                        (see the configuration file of that model for details)
%
%
% PLOTTING OF RESULTS:
%     To plot the trajectories of the inferred perceptual states (as implied by the averaged
%     parameters), there is a function <modelname>_plotTraj(...) for each perceptual model. This
%     takes the structure returned by bpa(...) as its only argument.
%
%     Additionally, the function tapas_fit_plotCorr(...) plots the posterior correlation of the
%     averaged parameters. It takes the structure returned by bpa(...) as its only
%     argument. Note that this function only works if the optimization algorithm makes the
%     posterior correlation available in est.optim.Corr for all of the estimate structures handed
%     to bpa(...).
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Number of estimates to average
n = size(varargin,2);

% Inputs
u = varargin{1}.u;

% Determine the models involved
prc_model = varargin{1}.c_prc.model;
obs_model = varargin{1}.c_obs.model;

% Get priors
prc_priormus = varargin{1}.c_prc.priormus;
prc_priorsas = varargin{1}.c_prc.priorsas;
obs_priormus = varargin{1}.c_obs.priormus;
obs_priorsas = varargin{1}.c_obs.priorsas;

% Check whether everything matches up
for i = 2:n
    if ~strcmp(prc_model,varargin{i}.c_prc.model)
        error('tapas:hgf:bpa:PrcModNoMatch', 'Perceptual models do not match.');
    end

    if ~strcmp(obs_model,varargin{i}.c_obs.model)
        error('tapas:hgf:bpa:ObsModNoMatch', 'Observation models do not match.');
    end

    if ~isequalwithequalnans(prc_priormus,varargin{i}.c_prc.priormus) || ~isequalwithequalnans(prc_priorsas,varargin{i}.c_prc.priorsas)
        error('tapas:hgf:bpa:PrcPriorsNoMatch', 'Perceptual priors do not match.');
    end

    if ~isequalwithequalnans(obs_priormus,varargin{i}.c_obs.priormus) || ~isequalwithequalnans(obs_priorsas,varargin{i}.c_obs.priorsas)
        error('tapas:hgf:bpa:ObsPriorsNoMatch', 'Observation priors do not match.');
    end

    if ~isequalwithequalnans(u(:),varargin{i}.u(:))
        disp(['Warning: inputs for argument number ' num2str(i) ' do not match those for first argument.']);
    end
end

% Record configuration
bpa       = struct;
bpa.u     = u;
bpa.ign   = [];
bpa.c_prc = varargin{1}.c_prc;
bpa.c_obs = varargin{1}.c_obs;

% Determine indices of parameters that have been optimized (i.e., those that are not fixed or NaN)
opt_idx = [bpa.c_prc.priorsas, bpa.c_obs.priorsas];
opt_idx(isnan(opt_idx)) = 0;
opt_idx = find(opt_idx);

% Prior precision
priorsas = [prc_priorsas, obs_priorsas];
H0 = diag(1./priorsas(opt_idx));

% Posterior precision and covariance
H = (1-n).*H0;

for i=1:n
    H = H + varargin{i}.optim.H;
end

Sigma = inv(H);
Sigma = tapas_nearest_psd(Sigma);
Corr = tapas_Cov2Corr(Sigma);

% Record results
bpa.optim.H     = H;
bpa.optim.Sigma = Sigma;
bpa.optim.Corr  = Corr;

% Prior mean
priormus = [prc_priormus, obs_priormus]';
mu0 = priormus(opt_idx);

% Posterior mean
mu = (1-n).*H0*mu0;

for i=1:n
    mui = [varargin{i}.p_prc.ptrans, varargin{i}.p_obs.ptrans]';
    mui = mui(opt_idx);
    mu = mu + varargin{i}.optim.H*mui;
end

mu = Sigma*mu;

% Replace optimized values in priormus with averaged values
ptrans = priormus';
ptrans(opt_idx) = mu';

% Separate perceptual and observation parameters
n_prcpars = length(bpa.c_prc.priormus);
ptrans_prc = ptrans(1:n_prcpars);
ptrans_obs = ptrans(n_prcpars+1:end);

% Transform MAP parameters back to their native space
[dummy, bpa.p_prc]   = bpa.c_prc.transp_prc_fun(bpa, ptrans_prc);
[dummy, bpa.p_obs]   = bpa.c_obs.transp_obs_fun(bpa, ptrans_obs);
bpa.p_prc.p      = bpa.c_prc.transp_prc_fun(bpa, ptrans_prc);
bpa.p_obs.p      = bpa.c_obs.transp_obs_fun(bpa, ptrans_obs);

% Store transformed MAP parameters
bpa.p_prc.ptrans = ptrans_prc;
bpa.p_obs.ptrans = ptrans_obs;

% Store representations at MAP estimate
bpa.traj = bpa.c_prc.prc_fun(bpa, bpa.p_prc.p);

% Print results
disp(' ')
disp('Results:');
disp(bpa.p_prc)
disp(bpa.p_obs)

end

author: translationalneuromodeling

repo: tapas