text:
function [str] = listBiGGModels()
% Reads from url http://bigg.ucsd.edu/api/v2/models
%
% USAGE:
%
%    [str] = listBiGGModels()
%
% OUTPUT:
%    str:    result of url read
%
% .. http://bigg.ucsd.edu/web_api
% .. Get a list of models
str = urlread('http://bigg.ucsd.edu/api/v2/models');

%{
%str = "results_count": 78, "results": [{"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 72, "bigg_id": "e_coli_core", "gene_count": 137, "reaction_count": 95}, {"organism": "Homo sapiens", "metabolite_count": 342, "bigg_id": "iAB_RBC_283", "gene_count": 346, "reaction_count": 469}, {"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 1668, "bigg_id": "iAF1260", "gene_count": 1261, "reaction_count": 2382}, {"organism": "Methanosarcina barkeri str. Fusaro", "metabolite_count": 628, "bigg_id": "iAF692", "gene_count": 692, "reaction_count": 690}, {"organism": "Geobacter metallireducens GS-15", "metabolite_count": 1109, "bigg_id": "iAF987", "gene_count": 987, "reaction_count": 1285}, {"organism": "Escherichia coli APEC O1", "metabolite_count": 1944, "bigg_id": "iAPECO1_1312", "gene_count": 1313, "reaction_count": 2736}, {"organism": "Homo sapiens", "metabolite_count": 738, "bigg_id": "iAT_PLT_636", "gene_count": 636, "reaction_count": 1008}, {"organism": "Escherichia coli BL21(DE3)", "metabolite_count": 1945, "bigg_id": "iB21_1397", "gene_count": 1337, "reaction_count": 2742}, {"organism": "Escherichia coli BW2952", "metabolite_count": 1949, "bigg_id": "iBWG_1329", "gene_count": 1328, "reaction_count": 2742}, {"organism": "Escherichia coli CFT073", "metabolite_count": 1938, "bigg_id": "ic_1306", "gene_count": 1307, "reaction_count": 2727}, {"organism": "Escherichia coli O127:H6 str. E2348/69", "metabolite_count": 1919, "bigg_id": "iE2348C_1286", "gene_count": 1284, "reaction_count": 2704}, {"organism": "Escherichia coli 042", "metabolite_count": 1926, "bigg_id": "iEC042_1314", "gene_count": 1314, "reaction_count": 2715}, {"organism": "Escherichia coli 55989", "metabolite_count": 1953, "bigg_id": "iEC55989_1330", "gene_count": 1330, "reaction_count": 2757}, {"organism": "Escherichia coli ABU 83972", "metabolite_count": 1944, "bigg_id": "iECABU_c1320", "gene_count": 1320, "reaction_count": 2732}, {"organism": "Escherichia coli B str. REL606", "metabolite_count": 1953, "bigg_id": "iECB_1328", "gene_count": 1329, "reaction_count": 2749}, {"organism": "Escherichia coli 'BL21-Gold(DE3)pLysS AG'", "metabolite_count": 1954, "bigg_id": "iECBD_1354", "gene_count": 1354, "reaction_count": 2749}, {"organism": "Escherichia coli BL21(DE3)", "metabolite_count": 1945, "bigg_id": "iECD_1391", "gene_count": 1333, "reaction_count": 2742}, {"organism": "Escherichia coli str. K-12 substr. DH10B", "metabolite_count": 1947, "bigg_id": "iECDH10B_1368", "gene_count": 1327, "reaction_count": 2743}, {"organism": "Escherichia coli DH1", "metabolite_count": 1949, "bigg_id": "iEcDH1_1363", "gene_count": 1363, "reaction_count": 2751}, {"organism": "Escherichia coli DH1", "metabolite_count": 1950, "bigg_id": "iECDH1ME8569_1439", "gene_count": 1439, "reaction_count": 2756}, {"organism": "Escherichia coli E24377A", "metabolite_count": 1974, "bigg_id": "iEcE24377_1341", "gene_count": 1341, "reaction_count": 2764}, {"organism": "Escherichia coli ED1a", "metabolite_count": 1929, "bigg_id": "iECED1_1282", "gene_count": 1279, "reaction_count": 2707}, {"organism": "Escherichia coli O157:H7 str. EC4115", "metabolite_count": 1918, "bigg_id": "iECH74115_1262", "gene_count": 1262, "reaction_count": 2695}, {"organism": "Escherichia coli HS", "metabolite_count": 1965, "bigg_id": "iEcHS_1320", "gene_count": 1321, "reaction_count": 2754}, {"organism": "Escherichia coli IAI1", "metabolite_count": 1970, "bigg_id": "iECIAI1_1343", "gene_count": 1343, "reaction_count": 2766}, {"organism": "Escherichia coli IAI39", "metabolite_count": 1957, "bigg_id": "iECIAI39_1322", "gene_count": 1321, "reaction_count": 2722}, {"organism": "Escherichia coli NA114", "metabolite_count": 1927, "bigg_id": "iECNA114_1301", "gene_count": 1301, "reaction_count": 2719}, {"organism": "Escherichia coli O103:H2 str. 12009", "metabolite_count": 1958, "bigg_id": "iECO103_1326", "gene_count": 1327, "reaction_count": 2759}, {"organism": "Escherichia coli O111:H- str. 11128", "metabolite_count": 1959, "bigg_id": "iECO111_1330", "gene_count": 1328, "reaction_count": 2761}, {"organism": "Escherichia coli O26:H11 str. 11368", "metabolite_count": 1965, "bigg_id": "iECO26_1355", "gene_count": 1355, "reaction_count": 2781}, {"organism": "Escherichia coli IHE3034", "metabolite_count": 1943, "bigg_id": "iECOK1_1307", "gene_count": 1304, "reaction_count": 2730}, {"organism": "Escherichia coli ATCC 8739", "metabolite_count": 1971, "bigg_id": "iEcolC_1368", "gene_count": 1368, "reaction_count": 2769}, {"organism": "Escherichia coli 536", "metabolite_count": 1943, "bigg_id": "iECP_1309", "gene_count": 1309, "reaction_count": 2740}, {"organism": "Escherichia coli O157:H7 str. Sakai", "metabolite_count": 1923, "bigg_id": "iECs_1301", "gene_count": 1301, "reaction_count": 2721}, {"organism": "Escherichia coli S88", "metabolite_count": 1944, "bigg_id": "iECS88_1305", "gene_count": 1305, "reaction_count": 2730}, {"organism": "Escherichia coli SE11", "metabolite_count": 1957, "bigg_id": "iECSE_1348", "gene_count": 1348, "reaction_count": 2769}, {"organism": "Escherichia coli SE15", "metabolite_count": 1951, "bigg_id": "iECSF_1327", "gene_count": 1327, "reaction_count": 2743}, {"organism": "Escherichia coli SMS-3-5", "metabolite_count": 1949, "bigg_id": "iEcSMS35_1347", "gene_count": 1347, "reaction_count": 2747}, {"organism": "Escherichia coli O157:H7 str. TW14359", "metabolite_count": 1920, "bigg_id": "iECSP_1301", "gene_count": 1299, "reaction_count": 2713}, {"organism": "Escherichia coli UMN026", "metabolite_count": 1935, "bigg_id": "iECUMN_1333", "gene_count": 1332, "reaction_count": 2741}, {"organism": "Escherichia coli W", "metabolite_count": 1975, "bigg_id": "iECW_1372", "gene_count": 1372, "reaction_count": 2783}, {"organism": "Escherichia coli KO11FL", "metabolite_count": 1974, "bigg_id": "iEKO11_1354", "gene_count": 1354, "reaction_count": 2779}, {"organism": "Escherichia coli ETEC H10407", "metabolite_count": 1964, "bigg_id": "iETEC_1333", "gene_count": 1333, "reaction_count": 2757}, {"organism": "Escherichia coli O55:H7 str. CB9615", "metabolite_count": 1919, "bigg_id": "iG2583_1286", "gene_count": 1283, "reaction_count": 2705}, {"organism": "Clostridium ljungdahlii DSM 13528", "metabolite_count": 698, "bigg_id": "iHN637", "gene_count": 637, "reaction_count": 785}, {"organism": "Helicobacter pylori 26695", "metabolite_count": 485, "bigg_id": "iIT341", "gene_count": 339, "reaction_count": 554}, {"organism": "Synechocystis sp. PCC 6803", "metabolite_count": 795, "bigg_id": "iJN678", "gene_count": 622, "reaction_count": 863}, {"organism": "Pseudomonas putida KT2440", "metabolite_count": 909, "bigg_id": "iJN746", "gene_count": 746, "reaction_count": 1056}, {"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 1805, "bigg_id": "iJO1366", "gene_count": 1367, "reaction_count": 2583}, {"organism": "Escherichia coli str. K-12 substr. MG1655", "metabolite_count": 761, "bigg_id": "iJR904", "gene_count": 904, "reaction_count": 1075}, {"organism": "Escherichia coli LF82", "metabolite_count": 1940, "bigg_id": "iLF82_1304", "gene_count": 1302, "reaction_count": 2727}, {"organism": "Thermotoga maritima MSB8", "metabolite_count": 570, "bigg_id": "iLJ478", "gene_count": 482, "reaction_count": 652}, {"organism": "Mus musculus", "metabolite_count": 2775, "bigg_id": "iMM1415", "gene_count": 1375, "reaction_count": 3726}, {"organism": "Saccharomyces cerevisiae S288c", "metabolite_count": 1226, "bigg_id": "iMM904", "gene_count": 905, "reaction_count": 1577}, {"organism": "Saccharomyces cerevisiae S288c", "metabolite_count": 1059, "bigg_id": "iND750", "gene_count": 750, "reaction_count": 1266}, {"organism": "Mycobacterium tuberculosis H37Rv", "metabolite_count": 826, "bigg_id": "iNJ661", "gene_count": 661, "reaction_count": 1025}, {"organism": "Escherichia coli O83:H1 str. NRG 857C", "metabolite_count": 1945, "bigg_id": "iNRG857_1313", "gene_count": 1311, "reaction_count": 2736}, {"organism": "Yersinia pestis CO92", "metabolite_count": 1552, "bigg_id": "iPC815", "gene_count": 815, "reaction_count": 1961}, {"organism": null, "metabolite_count": 1706, "bigg_id": "iRC1080", "gene_count": 1086, "reaction_count": 2191}, {"organism": "Shigella flexneri 2a str. 2457T", "metabolite_count": 1914, "bigg_id": "iS_1188", "gene_count": 1188, "reaction_count": 2620}, {"organism": "Staphylococcus aureus subsp. aureus N315", "metabolite_count": 655, "bigg_id": "iSB619", "gene_count": 619, "reaction_count": 743}, {"organism": "Shigella boydii CDC 3083-94", "metabolite_count": 1912, "bigg_id": "iSbBS512_1146", "gene_count": 1147, "reaction_count": 2592}, {"organism": "Shigella boydii Sb227", "metabolite_count": 1910, "bigg_id": "iSBO_1134", "gene_count": 1134, "reaction_count": 2592}, {"organism": "Shigella dysenteriae Sd197", "metabolite_count": 1890, "bigg_id": "iSDY_1059", "gene_count": 1059, "reaction_count": 2540}, {"organism": "Shigella flexneri 2a str. 301", "metabolite_count": 1917, "bigg_id": "iSF_1195", "gene_count": 1195, "reaction_count": 2631}, {"organism": "Shigella flexneri 5 str. 8401", "metabolite_count": 1917, "bigg_id": "iSFV_1184", "gene_count": 1184, "reaction_count": 2622}, {"organism": "Shigella flexneri 2002017", "metabolite_count": 1918, "bigg_id": "iSFxv_1172", "gene_count": 1169, "reaction_count": 2639}, {"organism": "Shigella sonnei Ss046", "metabolite_count": 1938, "bigg_id": "iSSON_1240", "gene_count": 1240, "reaction_count": 2694}, {"organism": "Escherichia coli UM146", "metabolite_count": 1944, "bigg_id": "iUMN146_1321", "gene_count": 1319, "reaction_count": 2736}, {"organism": "Escherichia coli UMNK88", "metabolite_count": 1971, "bigg_id": "iUMNK88_1353", "gene_count": 1353, "reaction_count": 2778}, {"organism": "Escherichia coli UTI89", "metabolite_count": 1942, "bigg_id": "iUTI89_1310", "gene_count": 1310, "reaction_count": 2726}, {"organism": "Escherichia coli W", "metabolite_count": 1975, "bigg_id": "iWFL_1372", "gene_count": 1372, "reaction_count": 2783}, {"organism": "Escherichia coli str. K-12 substr. W3110", "metabolite_count": 1953, "bigg_id": "iY75_1357", "gene_count": 1358, "reaction_count": 2760}, {"organism": "Klebsiella pneumoniae subsp. pneumoniae MGH 78578", "metabolite_count": 1658, "bigg_id": "iYL1228", "gene_count": 1229, "reaction_count": 2262}, {"organism": "Bacillus subtilis subsp. subtilis str. 168", "metabolite_count": 991, "bigg_id": "iYO844", "gene_count": 844, "reaction_count": 1250}, {"organism": "Escherichia coli O157:H7 str. EDL933", "metabolite_count": 1923, "bigg_id": "iZ_1308", "gene_count": 1308, "reaction_count": 2722}, {"organism": "Homo sapiens", "metabolite_count": 2766, "bigg_id": "RECON1", "gene_count": 1905, "reaction_count": 3742}, {"organism": "Salmonella enterica subsp. enterica serovar Typhimurium str. LT2", "metabolite_count": 1802, "bigg_id": "STM_v1_0", "gene_count": 1271, "reaction_count": 2545}]}

%curl 'http://bigg.ucsd.edu/api/v2/models'
%}

author: opencobra

repo: cobratoolbox


text:
function [ joints, s_mat, w_mat ] = InverseKinematics( pose, nsparam, rconf, robot_type )
% INVERSEKINEMATICS
%
% Calculates the inverse kinematics for a KUKA LBR iiwa manipulator.
%
% Input:  pose    - Homogeneous Matrix size(4,4)
%         nsparam - Arm Angle
%         rconf   - Robot Configuration 8-bit number
% Output: joints  - Joint values size(1,7).
%         s_mat   - Shoulder joint matrices As, Bs and Cs
%         w_mat   - Wrist joint matrices Aw, Bw and Cw
%         robot_type - iiwa7 (0) or iiwa14 (~0)

[ arm, elbow, wrist ] = Configuration(rconf);

%Tolerance
tol = 1e-8;

%Robot parameters
%Link length
l = zeros(1,6);
if robot_type == 0
    l = [0.34 0.4 0.4 0.126]; % iiwa7
else
    l = [0.36 0.42 0.4 0.126]; % iiwa14
end

%Denavit-Hartenberg parameters 7 DoF
%DH: [a, alpha,    d, theta]
dh = [0	 -pi/2	 l(1)    0;
      0	  pi/2      0    0;
      0	  pi/2   l(2)	 0;
      0  -pi/2      0    0;
      0  -pi/2   l(3)    0;
      0   pi/2      0    0;
      0      0   l(4)    0];

%Number of joints
nj = size(dh,1);

% Joint values of virtual manipulator
joints = zeros(1,nj);

% Shoulder rotation matrices
s_mat = zeros(3,3,3);

% Wrist rotation matrices
w_mat = zeros(3,3,3);

xend = pose(1:3,4); % end-effector position from base
xs = [0 0 dh(1,3)]'; % shoulder position from base
xwt = [0 0 dh(end,3)]'; % end-effector position from wrist
xw = xend - pose(1:3,1:3)*xwt; % wrist position from base
xsw = xw - xs; % shoulder to wrist vector
usw = unit(xsw);

lbs = l(1);
lse = l(2); % upper arm length (shoulder to elbow)
lew = l(3); % lower arm length (elbow to wrist)

%Check if pose is within arm+forearm reach
assert(norm(xsw) < lse + lew && norm(xsw) > lse - lew, 'Specified pose outside reacheable workspace');

% -- Joint 4 --
% Elbow joint can be directly calculated since it does only depend on the
% robot configuration and the xsw vector
assert(abs((norm(xsw)^2 - lse^2 - lew^2)-(2*lse*lew)) > tol, 'Elbow singularity. Tip at reach limit.');
% Cosine law - According to our robot, joint 4 rotates backwards
joints(4) = elbow * acos((norm(xsw)^2 - lse^2 - lew^2)/(2*lse*lew));

%Added

T34 = dh_calc(dh(4,1),dh(4,2),dh(4,3),joints(4));
R34 = T34(1:3,1:3);

% Shoulder Joints
% First compute the reference joint angles when the arm angle is zero.
[~, R03_o, ~] = ReferencePlane(pose, elbow, robot_type);

skew_usw = skew(usw);
% Following eq. (15), the auxiliary matrixes As Bs and Cs can be calculated
% by substituting eq. (6) into (9).
% R0psi = I3 + sin(psi)*skew_usw + (1-cos(psi))*skew_usw²    (6)
% R03 = R0psi * R03_o                                         (9)
% Substituting (distributive prop.) we get:
% R03 = R03_o*skew_usw*sin(psi) + R03_o*(-skew_usw²)*cos(psi) + R03_o(I3 + skew_usw²)
% R03 =      As       *sin(psi) +        Bs         *cos(psi) +          Cs
As = skew_usw * R03_o;
Bs = -skew_usw^2 * R03_o;
Cs = (usw*usw') * R03_o;

psi = nsparam;
R03 = As*sin(psi) + Bs*cos(psi) + Cs;

% T03 transformation matrix (DH parameters)
%[ cos(j1)*cos(j2)*cos(j3) - sin(j1)*sin(j3), cos(j1)*sin(j2), cos(j3)*sin(j1) + cos(j1)*cos(j2)*sin(j3), 0.4*cos(j1)*sin(j2)]
%[ cos(j1)*sin(j3) + cos(j2)*cos(j3)*sin(j1), sin(j1)*sin(j2), cos(j2)*sin(j1)*sin(j3) - cos(j1)*cos(j3), 0.4*sin(j1)*sin(j2)]
%[                          -cos(j3)*sin(j2),         cos(j2),                          -sin(j2)*sin(j3),  0.4*cos(j2) + 0.34]
%[                                         0,               0,                                         0,                   1]
joints(1) = atan2(arm * R03(2,2), arm * R03(1,2));
joints(2) = arm * acos(R03(3,2));

joints(3) = atan2(arm * -R03(3,3), arm * -R03(3,1));

Aw = R34' * As' * pose(1:3,1:3);
Bw = R34' * Bs' * pose(1:3,1:3);
Cw = R34' * Cs' * pose(1:3,1:3);

R47 = Aw*sin(psi) + Bw*cos(psi) + Cw;

% T47 transformation matrix (DH parameters)
%[ cos(j5)*cos(j6)*cos(j7) - sin(j5)*sin(j7), - cos(j7)*sin(j5) - cos(j5)*cos(j6)*sin(j7), cos(j5)*sin(j6), (63*cos(j5)*sin(j6))/500]
%[ cos(j5)*sin(j7) + cos(j6)*cos(j7)*sin(j5),   cos(j5)*cos(j7) - cos(j6)*sin(j5)*sin(j7), sin(j5)*sin(j6), (63*sin(j5)*sin(j6))/500]
%[                          -cos(j7)*sin(j6),                             sin(j6)*sin(j7),         cos(j6),   (63*cos(j6))/500 + 2/5]
%[                                         0,                                           0,               0,                        1]
joints(5) = atan2(wrist * R47(2,3), wrist * R47(1,3));
joints(6) = wrist * acos(R47(3,3));
joints(7) = atan2(wrist * R47(3,2), wrist * -R47(3,1));

% Grouping Shoulder and Wrist matrices that will be used by the joint
% limit algorithms
s_mat(:,:,1) = As;
s_mat(:,:,2) = Bs;
s_mat(:,:,3) = Cs;
w_mat(:,:,1) = Aw;
w_mat(:,:,2) = Bw;
w_mat(:,:,3) = Cw;

author: xuhuairuogu

repo: V-REP-Simulation-Projects

path: github-repos/MATLAB/xuhuairuogu-V-REP-Simulation-Projects/V-REP-Simulation-Projects-841b944af4ea3a8fb250578d36434515f577f411/admittance control(adapted from an example in the book--A Systematic Approach to Learning Robot Programming with ROS)/iiwa14_kinematics/InverseKinematics.m


text:
function resultingMatrix = blockDiagonalMatrix(param,varargin)
% blockDiagonalMatrix provides an interface for creating block diagonal matrices

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.
  prev_class  = class(varargin{1});

  % Check if the matrices provided for building the block diagonal are all variables
  % of the same class
  for i=1:nargin-1
    if(isa(varargin{i},prev_class))
      prev_class=class(varargin{i});
    else
      error('All the inputs arguments must be of the same class')
    end
  end

  % When running the model using symbolic variables, the creation of the block diagonal has to be carried out
  % differently if Octave is used
  if(isa(varargin{1},'casadi.SX') || isa(varargin{1},'casadi.MX'))
    if(param.isMatlab)
      resultingMatrix = blkdiag(varargin{:});
    else
      size_tot                                                          = [0,0];
      for i=1:nargin-1
        size_tot = size_tot + size(varargin{i});
      end
      resultingMatrix                                                             = SX.zeros(size_tot);
      start_position = [1,1];
      % Create the block diagonal
      for i=1:nargin-1
        resultingMatrix(start_position(1):start_position(1)+size(varargin{i},1)-1,start_position(2):start_position(2)+size(varargin{i},2)-1) = varargin{i};
        start_position(1) = start_position(1) + size(varargin{i},1);
        start_position(2) = start_position(2) + size(varargin{i},2);
      end
    end
  else
    resultingMatrix = blkdiag(varargin{:});
  end

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
i=iroi(im);

% attempt to find the markers using simple vision methods
i2=imorph(i,ones(3,3),'max');	% eliminate the markers
iff=flatfield(i2);		% and compute the flat field
idisp(iff-i)
pause
[k,s]=kmeans(iff-i,3);		% now k-means segmentation
idisp(s)
pause
[l,nblobs]=ilabel(s==3);	% and assume brightest points are markers
idisp(l)
colormap(hot(nblobs));
pause
f=iblobs(l)			% now extract binary image features
hold on
for i=1:12
	if f(i).area < 100,
		plot(f(i).xc,f(i).yc,'*w')
	end
end

% now enter the refinement phase

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)

%function [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)
%
%% READ_RICOH_DATA reads continuous or averaged MEG data
%% generated by the RICOH MEG system and software,
%% and allows the data to be used in FieldTrip.
%%
%% Use as
%%   [dat] = read_ricoh_data(filename, hdr, begsample, endsample, chanindx)
%%
%% This is a wrapper function around the function getRData
%%
%% See also READ_RICOH_HEADER, READ_RICOH_EVENT

% ensure that the required toolbox is on the path
ft_hastoolbox('ricoh_meg_reader', 1);

hdr = hdr.orig; % use the original Ricoh header, not the FieldTrip header

% default is to select all channels
if nargin<5
  chanindx = 1:hdr.channel_count;
end

handles = definehandles;

switch hdr.acq_type
  case handles.AcqTypeEvokedAve
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getRData(filename, start_sample, sample_length);

  case handles.AcqTypeContinuousRaw
    % dat is returned as double
    start_sample  = begsample - 1; % samples start at 0
    sample_length = endsample - begsample + 1;
    epoch_count   = 1;
    start_epoch   = 0;
    dat = getRData(filename, start_sample, sample_length);

    % Unlike Yokogawa system, "AcqTypeEvokedRaw" is not supported for Ricoh system.

  otherwise
    ft_error('unknown data type');
end

if size(dat,1)~=hdr.channel_count
  ft_error('could not read all channels');
elseif size(dat,2)~=(endsample-begsample+1)
  ft_error('could not read all samples');
end

% select only the desired channels
dat = dat(chanindx,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function CPT2 = mk_named_CPT(family_names, names, dag, CPT1)
% MK_NAMED_CPT Permute the dimensions of a CPT so they agree with the internal numbering convention
% CPT2 = mk_named_CPT(family_names, names, dag, CPT1)
%
% This is best explained by example.
% Consider the following directed acyclic graph
%
%     C
%   /   \
%  R     S
%   \   /
%     W
%
% where all arcs point down.
% When we create the CPT for node W, we consider S as its first parent, and R as its
% second, and hence write
%
%      S R W
% CPT1(1,1,:) = [1.0 0.0];
% CPT1(2,1,:) = [0.2 0.8];  % P(W=1 | R=1, S=2) = 0.2
% CPT1(1,2,:) = [0.1 0.9];
% CPT1(2,2,:) = [0.01 0.99];
%
% However, when we create the dag using mk_adj_mat, the nodes get topologically sorted,
% and by chance, node R preceeds node S in this ordering.
% Hence we should have written
%
%      R S W
% CPT2(1,1,:) = [1.0 0.0];
% CPT2(2,1,:) = [0.1 0.9];
% CPT2(1,2,:) = [0.2 0.8]; % P(W=1 | R=1, S=2) = 0.2
% CPT2(2,2,:) = [0.01 0.99];
%
% Since we do not know the order of the nodes in advance, we can write
%   CPT2 = mk_named_CPT({'S', 'R', 'W'}, names, dag, CPT1)
% where 'S', 'R', 'W' are the order of the dimensions we assumed (the child node must be last in this list),
% and names{i} is the name of the i'th node.

n = length(family_names);
family_nums = zeros(1,n);
for i=1:n
  family_nums(i) = stringmatch(family_names{i}, names); % was strmatch
end

fam = family(dag, family_nums(end));
perm = zeros(1,n);
for i=1:n
  %  perm(i) = find(family_nums(i) == fam);
  perm(i) = find(fam(i) == family_nums);
end

CPT2 = permute(CPT1, perm);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [level,bw] = thresh_tool(im,cmap,defaultLevel) %mainfunction
%THRESH_TOOL  Interactively select intensity level for image thresholding.
%   THRESH_TOOL launches a GUI (graphical user interface) for thresholding
%   an intensity input image, IM. IM is displayed in the top left corner. A
%   colorbar and IM's histogram are displayed on the bottom. A line on the
%   histogram indicates the current threshold level. A binary image is
%   displayed in the top right based on the selected level.  To change the
%   level, click and drag the line. The output image updates automatically.
%
%   There are two ways to use this tool.
%
%Mode 1 - nonblocking behavior:
%   THRESH_TOOL(IM) launches GUI tool.  You can continue using the MATLAB
%   Desktop.  Since no results are needed, the function does not block
%   execution of other commands.
%
%   THRESH_TOOL(IM,CMAP) allows the user to specify the colormap, CMAP.  If
%   not specified, the default colormap is used.
%
%   THRESH_TOOL(IM,CMAP,DEFAULTLEVEL) allows the user to specify the
%   default threshold level. If not specified, DEFAULTLEVEL is determined
%   by GRAYTHRESH. Valid values for DEFAULTLEVEL must be consistent with
%   the data type of IM for integer intensity images: uint8 [0,255], uint16
%   [0,65535], int16 [-32768,32767].
%
%   Example
%       x = imread('rice.png');
%       thresh_tool(x)          %no return value, so MATLAB keeps running
%
%Mode 2 - blocking behavior:
%   LEVEL = THRESH_TOOL(...) returns the user selected level, LEVEL, and
%   MATLAB waits for the result before proceeding.  This blocking behavior
%   mode allows the tool to be inserted into an image processing algorithm
%   to support an automated workflow.
%
%   [LEVEL,BW] = THRESH_TOOL(...) also returns the thresholded binary
%   output image, BW.
%
%   Example
%       x = imread('rice.png');
%       lev = thresh_tool(x)    %MATLAB waits for GUI tool to finish
%
%See also COLORMAP, GRAYTHRESH, IM2BW.

%defensive programming
error(nargchk(1,3,nargin))
error(nargoutchk(0,2,nargout))

%validate defaultLevel within range
if nargin>2 %uer specified DEFAULTLEVEL
  dataType = class(im);
  switch dataType
    case 'uint8','uint16','int16'
      if defaultLevel<intmin(dataType) | defaultLevel>intmax(dataType)
        error(['Specified DEFAULTLEVEL outside class range for ' dataType])
      elseif defaultLevel<min(im(:)) | defaultLevel>max(im(:))
        error('Specified DEFAULTLEVEL outside data range for IM')
      end
    case 'double','single'
      %okay, do nothing
    otherwise
      error(['Unsupport image type ' dataType])
  end %switch
end

max_colors=1000;    %practical limit

%calculate bins centers
color_range = double(limits(im));
if isa(im,'uint8') %special case [0 255]
  color_range = [0 255];
  num_colors = 256;
  di = 1;
elseif isinteger(im)
  %try direct indices first
  num_colors = diff(color_range)+1;
  if num_colors<max_colors %okay
    di = 1;                                 %inherent bins
  else %too many levels
    num_colors = max_colors;                %practical limit
    di = diff(color_range)/(num_colors-1);
  end
else %noninteger
  %try infering discrete resolution first (intensities often quantized)
  di = min(diff(sort(unique(im(:)))));
  num_colors = round(diff(color_range)/di)+1;
  if num_colors>max_colors %too many levels
    num_colors = max_colors;                %practical limit
    di = diff(color_range)/(num_colors-1);
  end
end
bin_ctrs = [color_range(1):di:color_range(2)];
FmtSpec = ['%.' num2str(ceil(-log10(di))) 'f'];

%new figure - interactive GUI tool for level segmenting
h_fig = figure;
set(h_fig,'ToolBar','Figure')
if nargin>1 && isstr(cmap) && strmatch(lower(cmap),'gray')
  full_map = gray(num_colors);
elseif nargin>1 && isnumeric(cmap) && length(size(cmap))==2 && size(cmap,2)==3
  full_map = cmap;
else
  full_map = jet(num_colors);
end
setappdata(h_fig,'im',im)
setappdata(h_fig,'FmtSpec',FmtSpec)

%top left - input image
h_ax1 = axes('unit','norm','pos',[0.05 0.35 0.4 0.60]);
rgb = im2rgb(im,full_map);
  function rgb = im2rgb(im,full_map); %nested
    %coerce intensities into gray range [0,1]
    gray = imadjust(im,[],[0 1]);
    %generate indexed image
    num_colors = size(full_map,1);
    ind = gray2ind(gray,num_colors);
    %convert indexed image to RGB
    rgb = ind2rgb(ind,full_map);
  end %im2rgb

image(rgb), axis image
%subimage(im,full_map)
axis off, title('Input Image')

%top right - segmented (eventually)
h_ax2 = axes('unit','norm','pos',[0.55 0.35 0.4 0.60]);
axis off
setappdata(h_fig,'h_ax2',h_ax2)

%next to bottom - intensity distribution
h_hist = axes('unit','norm','pos',[0.05 0.1 0.9 0.2]);
n = hist(double(im(:)),bin_ctrs);
bar(bin_ctrs,n)
axis([color_range limits(n(2:end-1))]) %ignore saturated end scaling
set(h_hist,'xtick',[],'ytick',[])
title('Intensity Distribution')

%very bottom - colorbar
h_cbar = axes('unit','norm','pos',[0.05 0.05 0.9 0.05],'tag','thresh_tool_cbar');
subimage(color_range,[0.5 1.5],1:num_colors,full_map)
set(h_cbar,'ytick',[],'xlim',color_range)
axis normal

v=version;
if str2num(v(1:3))>=7
  %link top axes (pan & zoom)
  linkaxes([h_ax1 h_ax2])
  %link bottom axes (X only - pan & zoom)
  linkaxes([h_hist h_cbar],'x')
end

%colorbar tick locations
set(h_cbar,'xtick',color_range)

%threshold level - initial guess (graythresh)
if nargin>2 %user specified default level
  my_level = defaultLevel;
else %graythresh default
  lo = double(color_range(1));
  hi = double(color_range(2));
  norm_im = (double(im)-lo)/(hi-lo);
  norm_level = graythresh(norm_im); %GRAYTHRESH assumes DOUBLE range [0,1]
  my_level = norm_level*(hi-lo)+lo;
end

%display level as vertical line
axes(h_hist)
h_lev = vline(my_level,'-');
set(h_lev,'LineWidth',2,'color',0.5*[1 1 1],'UserData',my_level)
setappdata(h_fig,'h_lev',h_lev)

%attach draggable behavior for user to change level
move_vline(h_lev,@update_plot);

axes(h_cbar)
y_lim = get(h_cbar,'ylim');

% PLACE TEXT LOCATION ON COLORBAR (Laurens)
%h_text = text(my_level,mean(y_lim),num2str(round(my_level)));
h_text = text(my_level,mean(y_lim),'dummy','HorizontalAlignment','Center');
if nargin<2
  text_color = 0.5*[1 1 1];
else
  text_color = 'm';
end
set(h_text,'FontWeight','Bold','color',text_color,'Tag','cbar_text')
movex_text(h_text,my_level)
%%%%%%%%%%%%%%%%%%%%%%%%

%segmented image
bw = im>my_level;
axes(h_ax2)
hold on
subimage(bw), axis off, axis ij
hold off
title('Segmented')

update_plot

%add reset button (resort to initial guess)
h_reset = uicontrol('unit','norm','pos',[0.0 0.95 .1 .05]);
set(h_reset,'string','Reset','callback',@ResetOriginalLevel)

if nargout>0 %return result(s)
  h_done = uicontrol('unit','norm','pos',[0.9 0.95 0.1 0.05]);
  set(h_done,'string','Done','callback','delete(gcbo)') %better
  %inspect(h_fig)
  set(h_fig,'WindowStyle','modal')
  waitfor(h_done)
  if ishandle(h_fig)
    h_lev = getappdata(gcf,'h_lev');
    level = mean(get(h_lev,'xdata'));
    if nargout>1
      h_im2 = findobj(h_ax2,'type','image');
      bw = logical(rgb2gray(get(h_im2,'cdata')));
    end
    delete(h_fig)
  else
    warning('THRESHTOOL:UserAborted','User Aborted - no return value')
    level = [];
  end
end

end %thresh_tool (mainfunction)

function ResetOriginalLevel(hObject,varargin) %subfunction
h_lev = getappdata(gcf,'h_lev');
init_level = get(h_lev,'UserData');
set(h_lev,'XData',init_level*[1 1])
text_obj = findobj('Type','Text','Tag','cbar_text');
movex_text(text_obj,init_level)
update_plot
end %ResetOriginalLevel (subfunction)

function update_plot %subfunction
im = getappdata(gcf,'im');
h_lev = getappdata(gcf,'h_lev');
my_level = mean(get(h_lev,'xdata'));
h_ax2 = getappdata(gcf,'h_ax2');
h_im2 = findobj(h_ax2,'type','image');
%segmented image
bw = (im>my_level);
rgb_version = repmat(double(bw),[1 1 3]);
set(h_im2,'cdata',rgb_version)
end %update_plot (subfunction)

%function rgbsubimage(im,map), error('DISABLED')

%----------------------------------------------------------------------
function move_vline(handle,DoneFcn) %subfunction
%MOVE_VLINE implements horizontal movement of line.
%
%  Example:
%    plot(sin(0:0.1:pi))
%    h=vline(1);
%    move_vline(h)
%
%Note: This tools strictly requires MOVEX_TEXT, and isn't much good
%      without VLINE by Brandon Kuczenski, available at MATLAB Central.
%<http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=1039&objectType=file>

% This seems to lock the axes position
set(gcf,'Nextplot','Replace')
set(gcf,'DoubleBuffer','on')

h_ax=get(handle,'parent');
h_fig=get(h_ax,'parent');
setappdata(h_fig,'h_vline',handle)
if nargin<2, DoneFcn=[]; end
setappdata(h_fig,'DoneFcn',DoneFcn)
set(handle,'ButtonDownFcn',@DownFcn)

  function DownFcn(hObject,eventdata,varargin) %Nested--%
    set(gcf,'WindowButtonMotionFcn',@MoveFcn)           %
    set(gcf,'WindowButtonUpFcn',@UpFcn)                 %
  end %DownFcn------------------------------------------%

  function UpFcn(hObject,eventdata,varargin) %Nested----%
    set(gcf,'WindowButtonMotionFcn',[])                 %
    DoneFcn=getappdata(hObject,'DoneFcn');              %
    if isstr(DoneFcn)                                   %
      eval(DoneFcn)                                     %
    elseif isa(DoneFcn,'function_handle')               %
      feval(DoneFcn)                                    %
    end                                                 %
  end %UpFcn--------------------------------------------%

  function MoveFcn(hObject,eventdata,varargin) %Nested------%
    h_vline=getappdata(hObject,'h_vline');                  %
    h_ax=get(h_vline,'parent');                             %
    cp = get(h_ax,'CurrentPoint');                          %
    xpos = cp(1);                                           %
    x_range=get(h_ax,'xlim');                               %
    if xpos<x_range(1), xpos=x_range(1); end                %
    if xpos>x_range(2), xpos=x_range(2); end                %
    XData = get(h_vline,'XData');                           %
    XData(:)=xpos;                                          %
    set(h_vline,'xdata',XData)                              %
    %update text                                            %
    text_obj = findobj('Type','Text','Tag','cbar_text');    %
    movex_text(text_obj,xpos)                               %
  end %MoveFcn----------------------------------------------%

end %move_vline(subfunction)

%----------------------------------------------------------------------
function [x,y] = limits(a) %subfunction
% LIMITS returns min & max values of matrix; else scalar value.
%
%   [lo,hi]=LIMITS(a) returns LOw and HIgh values respectively.
%
%   lim=LIMITS(a) returns 1x2 result, where lim = [lo hi] values

if nargin~=1 | nargout>2 %bogus syntax
  error('usage: [lo,hi]=limits(a)')
end

siz=size(a);

if prod(siz)==1 %scalar
  result=a;                         % value
else %matrix
  result=[min(a(:)) max(a(:))];     % limits
end

if nargout==1 %composite result
  x=result;                         % 1x2 vector
elseif nargout==2 %separate results
  x=result(1);                      % two scalars
  y=result(2);
else %no result
  ans=result                        % display answer
end

end %limits (subfunction)

%----------------------------------------------------------------------
function movex_text(h_txt,x_pos) %subfunction
FmtSpec=getappdata(get(get(h_txt,'parent'),'parent'),'FmtSpec');
msg=sprintf(FmtSpec,x_pos);
pos=get(h_txt,'position');
pos(1)=x_pos;
set(h_txt,'Position',pos,'String',msg)
end %movex_text

%--------------------------------------------------------------------------------------------------------------
function hhh=vline(x,in1,in2) %subfunction
% function h=vline(x, linetype, label)
%
% Draws a vertical line on the current axes at the location specified by 'x'.  Optional arguments are
% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line.  The
% label appears in the same color as the line.
%
% The line is held on the current axes, and after plotting the line, the function returns the axes to
% its prior hold state.
%
% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
% legends, but it is not findable by using findobj.  Specifying an output argument causes the function to
% return a handle to the line, so it can be manipulated or deleted.  Also, the HandleVisibility can be
% overridden by setting the root's ShowHiddenHandles property to on.
%
% h = vline(42,'g','The Answer')
%
% returns a handle to a green vertical line on the current axes at x=42, and creates a text object on
% the current axes, close to the line, which reads "The Answer".
%
% vline also supports vector inputs to draw multiple lines at once.  For example,
%
% vline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
%
% draws three lines with the appropriate labels and colors.
%
% By Brandon Kuczenski for Kensington Labs.
% brandon_kuczenski@kensingtonlabs.com
% 8 November 2001

% Downloaded 8/7/03 from MATLAB Central
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=1039&objectType=file

if length(x)>1  % vector input
    for I=1:length(x)
        switch nargin
        case 1
            linetype='r:';
            label='';
        case 2
            if ~iscell(in1)
                in1={in1};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            label='';
        case 3
            if ~iscell(in1)
                in1={in1};
            end
            if ~iscell(in2)
                in2={in2};
            end
            if I>length(in1)
                linetype=in1{end};
            else
                linetype=in1{I};
            end
            if I>length(in2)
                label=in2{end};
            else
                label=in2{I};
            end
        end
        h(I)=vline(x(I),linetype,label);
    end
else
    switch nargin
    case 1
        linetype='r:';
        label='';
    case 2
        linetype=in1;
        label='';
    case 3
        linetype=in1;
        label=in2;
    end

    g=ishold(gca);
    hold on

    y=get(gca,'ylim');
    h=plot([x x],y,linetype);
    if length(label)
        xx=get(gca,'xlim');
        xrange=xx(2)-xx(1);
        xunit=(x-xx(1))/xrange;
        if xunit<0.8
            text(x+0.01*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
        else
            text(x-.05*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
        end
    end

    if g==0
    hold off
    end
    set(h,'tag','vline','handlevisibility','off')
end % else

if nargout
    hhh=h;
end

end %vline (subfunction)

author: Sable

repo: mcbench-benchmarks


text:
function [k, sk] = whitefixedKernCompute(kern, x, x2)

% WHITEFIXEDKERNCOMPUTE Compute the WHITEFIXED kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the fixed parameter white noise
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the inpute matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the fixed parameter white noise
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : whitefixedKernParamInit, kernCompute, kernCreate, whitefixedKernDiagCompute
%
% COPYRIGHT : Nathaniel J. King, 2006
%
% MODIFICATIONS : Neil D. Lawrence, 2009

% KERN

if nargin < 3
  k = whiteKernCompute(kern, x);
else
  k = whiteKernCompute(kern, x, x2);
end
if nargout > 1
  sk = k;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function dat = proc_BeerLambert(dat, varargin)
% PROC_BEERLAMBERT - apply the BeerLambert function, which is useful to analyze
% NIRS data:  relative concentration is calculated as a function of total photon path length
%
% Synopsis:
%   DAT = nirs_LB(DAT, 'Property1',Value1, ...)
%
% Arguments:
%   DAT: data struct with NIRS cnt data split in the x field
%
%   OPT - struct or property/value list of optional properties:
%   'Citation'   - if set, the epsilon (extinction values) is taken from
%                  the specified citation number (see GetExtinctions.m for
%                  details). (default 1 = Gratzer et al)
%   'Epsilon'    - sets extinction coefficients manually (wl1 wl2 vs deoxy
%                  and oxy). In this case, citation is ignored. State in
%                  millimol/liter(?).
%   'Opdist'     - optode (source-detector) distance in cm (default 2.5)
%   'Ival'       - either 'all' (default) or a vector [start end] in samples
%                  specifying the baseline for the LB transformation
%   'DPF'        - differential pathlength factor: probabilistic average
%                  distance travelled by photons, default [5.98 7.15]
%
% Returns:
%   DAT - updated data struct with oxy/deoxy fields specifying absorption
%         values in mmol/l.
%
%
% See also: nirs_* nirsfile_* GetExtinctions
%
% Note: Based on the nirX Nilab toolbox functions u_LBG and u_popLBG.

% matthias.treder@tu-berlin.de 2011, mail@aewald.net 2013
% AE: changed channels label from wl1&wl2 to oxa and deoxy.
% Markus Wenzel 2013 (adapted it to the new toolbox)
% Jan Mehnert February 2014 (ready for public BBCI toolbox) (jan@mehnert.org)
% ToDo: possibility to define individual baseline; online functionability

props={'Citation'   1             'INT'
        'Opdist'    2.5           'DOUBLE'
        'Ival'      'all'         'CHAR|DOUBLE'
        'DPF'       [5.98 7.15]   'DOUBLE'
        'Epsilon'   []            'DOUBLE'
        'Verbose'   0             'BOOL'};

if nargin==0,
    dat= props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);

opt_checkProplist(opt, props);
misc_checkType(dat, 'STRUCT');

if strcmp(opt.Ival,'all')
  opt.Ival = [1 size(dat.x,1)];
end

s1=size(dat.x,1);
s2=size(dat.x,2)/2;

wl1 = dat.x(:,1:end/2) + eps;
wl2 = dat.x(:,end/2+1:end) + eps;

%% Get epsilon
if isempty(opt.Epsilon)
  if ~isfield(dat,'wavelengths')
    error('Wavelengths should be given in the .wavelengths field.')
  end
  [ext,nfo] = procutil_getExtinctions(dat.wavelengths,opt.Citation);
  if opt.Verbose, fprintf('Citation: %s\n',nfo), end
  epsilon = ext(:,1:2);

  % Divide by 1000 to obtain the required unit
  epsilon = epsilon/1000;

else
  epsilon = opt.Epsilon;
end

%% Arrange epsilon so that higher wavelength is on top
[mw,idx] = max(dat.wavelengths);
if idx==2 % higher wavelength is on bottom
  epsilon = flipud(epsilon);
  if opt.Verbose, fprintf('Epsilon matrix was rearranged so that higher WL is on top\n'), end
end

%% Apply Lambert-Beer law
Att_highWL= real(-log10( wl2 ./ ...
    ( repmat(mean(wl2(opt.Ival(1):opt.Ival(2),:),1), [s1,1]))   ));

Att_lowWL= real(-log10( wl1./ ...
    ( repmat(mean(wl1(opt.Ival(1):opt.Ival(2),:),1), [s1,1]))   ));

A=[];
A(:,1)=reshape(Att_highWL,s1*s2,1);
A(:,2)=reshape(Att_lowWL,s1*s2,1);

%----------------------------------
%       3.cc
%----------------------------------
% e=...looks like this
%               oxy-Hb         deoxy-Hb
% higherWL: 830 | e: 0.974       0.693
% lowerWL : 690 | e: 0.35         2.1

e= epsilon/10;

e2=   e.* [opt.DPF' opt.DPF']  .*  opt.Opdist;
c= ( inv(e2)*A'  )';

dat.x = reshape(c(:,1),s1,s2); %in mmol/l
dat.x = [dat.x reshape(c(:,2),s1,s2)]; %in mmol/l

%% Change Channel labels wl1 & wl2 to 'oxy' and 'deoxy'

dat.clab = strrep(dat.clab, 'highWL', 'oxy');
dat.clab = strrep(dat.clab, 'lowWL', 'deoxy');

dat.signal = 'NIRS (oxy, deoxy)';

dat.yUnit = 'mmol/l';

author: bbci


text:
function [data] = ft_determine_coordsys(data, varargin)

% FT_DETERMINE_COORDSYS plots a geometrical object, allowing you to perform
% a visual check on the coordinatesystem, the units and on the anatomical
% labels for the coordinate system axes.
%
% Use as
%   [dataout] = ft_determine_coordsys(datain, ...)
% where the input data structure can be either
%  - an anatomical MRI
%  - a cortical or head surface mesh
%  - an electrode, gradiometer or optode definition
%  - a volume conduction model of the head
% or most other FieldTrip structures that represent geometrical information.
%
% Additional optional input arguments should be specified as key-value pairs
% and can include
%   interactive  = string, 'yes' or 'no' (default = 'yes')
%   axisscale    = scaling factor for the reference axes and sphere (default = 1)
%   clim         = lower and upper anatomical MRI limits (default = [0 1])
%
% This function wil pop up a figure that allows you to check whether the
% alignment of the object relative to the coordinate system axes is correct
% and what the anatomical labels of the coordinate system axes are. You
% should switch on the 3D rotation option in the figure panel to rotate and
% see the figure from all angles. To change the anatomical labels of the
% coordinate system, you should press the corresponding keyboard button.
%
% Recognized and supported coordinate systems are 'ctf', 'bti', '4d', 'yokogawa',
% 'eeglab', 'neuromag', 'itab', 'acpc', 'spm', 'mni', 'fsaverage', 'tal', 'scanras',
% 'scanlps', 'dicom'.
%
% Furthermore, supported coordinate systems that do not specify the origin are 'ras',
% 'als', 'lps', etc. See https://www.fieldtriptoolbox.org/faq/coordsys for more
% details.
%
% See also FT_CONVERT_COORDSYS, FT_DETERMINE_UNITS, FT_CONVERT_UNITS, FT_PLOT_AXES, FT_PLOT_XXX

% Copyright (C) 2015-2021, 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$

dointeractive = ft_getopt(varargin, 'interactive', 'yes');
axisscale     = ft_getopt(varargin, 'axisscale', 1);  % this is used to scale the axmax and rbol
clim          = ft_getopt(varargin, 'clim', [0 1]);   % this is used to scale the orthoplot
fontsize      = ft_getopt(varargin, 'fontsize');      % this is passed to ft_plot_axes

data  = ft_checkdata(data, 'hasunit', 'yes');
dtype = ft_datatype(data);

% the high-level data structures are detected with ft_datatype, but there are
% also some low-level data structures that need to be supproted here
if strcmp(dtype, 'unknown')
  if isfield(data, 'fid') || (isfield(data, 'tri') && isfield(data, 'pos'))
    dtype = 'headshape';
  elseif isfield(data, 'hex') && isfield(data, 'pos')
    dtype = 'mesh';
  elseif isfield(data, 'tet') && isfield(data, 'pos')
    dtype = 'mesh';
  elseif ~strcmp(ft_headmodeltype(data), 'unknown')
    dtype = 'headmodel';
  elseif ~strcmp(ft_senstype(data), 'unknown')
    dtype = 'sens';
  end
elseif strcmp(dtype, 'mesh+label')
  % we don't care about the labels here
  dtype = 'mesh';
end

if isfield(data, 'coordsys') && ~isempty(data.coordsys)
  % ensure that it is in lower case
  data.coordsys = lower(data.coordsys);
  % print the interpretation of the coordinate system
  [labelx, labely, labelz] = coordsys2label(data.coordsys, 2, 0);
  fprintf('The positive x-axis is pointing towards %s\n', labelx);
  fprintf('The positive y-axis is pointing towards %s\n', labely);
  fprintf('The positive z-axis is pointing towards %s\n', labelz);
end

% plot the geometrical object
% the plotting style depends on the data content
figure;
switch dtype
  case 'volume'
    funparam = [];
    if isfield(data, 'anatomy')
      funparam = data.anatomy;
    elseif isfield(data, 'gray')
      funparam = data.gray;
    elseif isfield(data, 'white')
      funparam = data.white;
    elseif isfield(data, 'brick0')
      funparam = data.brick0; % used for an atlas
    elseif isfield(data, 'brick1')
      funparam = data.brick1; % used for an atlas
    else
      % try to determine it automatically
      fn = fieldnames(data);
      for i=1:length(fn)
        if isequal(size(data.(fn{i})), data.dim)
          funparam = data.(fn{i});
          break;
        end
      end
    end

    if isempty(funparam)
      ft_error('don''t know which volumetric parameter to plot');
    end

    % the volumetric data needs to be interpolated onto three orthogonal planes
    % determine a resolution that is close to, or identical to the original resolution
    [corner_vox, corner_head] = cornerpoints(data.dim, data.transform);
    diagonal_head = norm(range(corner_head));
    diagonal_vox  = norm(range(corner_vox));
    resolution    = (diagonal_head+eps)/(diagonal_vox+eps); % this is in units of "data.unit"

    % scale funparam between 0 and 1
    if ~isa(funparam, 'double') % avoid integer datatypes to allow for scaling
      funparam = double(funparam);
    end
    dmin = min(funparam(:));
    dmax = max(funparam(:));
    funparam  = (funparam-dmin)./(dmax-dmin);

    clear ft_plot_slice
    ft_plot_ortho(funparam, 'transform', data.transform, 'unit', data.unit, 'resolution', resolution, 'style', 'intersect', 'clim', clim);
    axis vis3d
    view([110 36]);

  case 'source'
    if isfield(data, 'inside') && ~isfield(data, 'tri')
      % only plot the source locations that are inside the volume conduction model
      ft_plot_mesh(data.pos(data.inside, :));
    else
      ft_plot_mesh(data, 'edgecolor','none', 'facecolor', [0.6 0.8 0.6], 'facealpha', 0.6);
    end
    camlight;

  case 'dip'
    ft_plot_mesh(data, 'edgecolor','none', 'facecolor', 'none');
    camlight;

  case 'headshape'
    ft_plot_headshape(data);
    camlight;

  case {'mesh', 'source+mesh'}
    ft_plot_mesh(data);
    camlight;

  case 'headmodel'
    ft_plot_headmodel(data);
    camlight;

  case {'grad' 'elec' 'opto' 'sens'}
    ft_plot_sens(data, 'label', 'label');
    camlight;

  case {'raw', 'timelock', 'freq', 'mvar', 'freqmvar', 'comp'}
    % the data may contain a gradiometer or electrode definition
    if isfield(data, 'grad')
      ft_plot_sens(data.grad);
    elseif isfield(data, 'elec')
      ft_plot_sens(data.elec, 'label', 'label');
    elseif isfield(data, 'opto')
      ft_plot_sens(data.opto, 'label', 'label');
    end

  case 'unknown'
end % switch dtype{k}

if isfield(data, 'tri')
  % this makes the 3-D object easier to understand
  camlight
  lighting gouraud
end

% plot the 3-D axes, labels, and sphere at the origin
ft_plot_axes(data, 'axisscale', axisscale, 'fontsize', fontsize);

if istrue(dointeractive)
  % ensure the figure is updated prior to asking the question
  % this was needed for FT_ELECTRODEPLACEMENT in combination with MATLAB 2022a
  drawnow

  if ~isfield(data, 'coordsys') || isempty(data.coordsys)
    % default is yes
    fprintf('The coordinate system is not specified.\n')
    value = smartinput('Do you want to change the anatomical labels for the axes [Y, n]? ', 'y');
  else
    % default is no
    fprintf('The coordinate system is specified as "%s".\n', data.coordsys)
    value = smartinput('Do you want to change the anatomical labels for the axes [y, N]? ', 'n');
  end

  if strcmpi(value, 'n')
    return
  end

  % interactively determine orientation
  orientation = '   ';
  while ~any(strcmp(orientation(1), {'r', 'l', 'a', 'p', 's', 'i'}))
    orientation(1) = smartinput('What is the anatomical label for the positive X-axis [r, l, a, p, s, i]? ', '');
  end
  while ~any(strcmp(orientation(2), {'r', 'l', 'a', 'p', 's', 'i'}))
    orientation(2) = smartinput('What is the anatomical label for the positive Y-axis [r, l, a, p, s, i]? ', '');
  end
  while ~any(strcmp(orientation(3), {'r', 'l', 'a', 'p', 's', 'i'}))
    orientation(3) = smartinput('What is the anatomical label for the positive Z-axis [r, l, a, p, s, i]? ', '');
  end

  % interactively determine origin
  origin = ' ';
  while ~any(strcmp(origin, {'a', 'i', 's', 'n'}))
    origin = input('Is the origin of the coordinate system at the a(nterior commissure), i(nterauricular), s(scanner origin), n(ot a landmark)? ', 's');
  end

  % some coordinate systems are identical or very similar, see https://www.fieldtriptoolbox.org/faq/coordsys
  if origin=='a' && strcmp(orientation, 'ras')
    coordsys = 'acpc'; % also used for spm, mni, tal
  elseif origin=='i' && strcmp(orientation, 'als')
    coordsys = 'ctf'; % also used for 4d, bti, eeglab
  elseif origin=='i' && strcmp(orientation, 'ras')
    coordsys = 'neuromag'; % also used for itab
  elseif origin=='s' && strcmp(orientation, 'ras')
    coordsys = 'scanras'; % also used for nifti
  elseif origin=='s' && strcmp(orientation, 'lps')
    coordsys = 'scanlps'; % also used for dicom
  else
    % only use the orientation, not the origin
    coordsys = orientation;
  end

  data.coordsys = coordsys;
end % if interactive


text:
function [mg,vg,pg,pv]=gaussmixg(m,v,w,n)
%GAUSSMIXG global mean, variance and mode of a GMM
%
% Usage: (1) gaussmixg(m,v,w)               % plot the mean and mode positions of a GMM
%        (2) [mg,vg]=gaussmixg(m,v,w)       % find global mean and covariance of a GMM
%        (3) [mg,vg,pg]=gaussmixg(m,v,w)    % find global mean,covariance and mode of a GMM
%        (4) [mg,vg,pg,pv]=gaussmixg(m,v,w) % ... also find log probability of the peak
%
%  Inputs:  M(k,p) = mixture means for pg(p)
%           V(k,p) or V(p,p,k) variances (diagonal or full)
%           W(k,1) = mixture weights
%           N      = maximum number of modes to find [default 1]
%
% Outputs: MG(1,p) = global mean
%          VG(p,p) = global covariance
%          PG(N,p) = sorted list of N modes
%          PV(N,1) = log pdf at the modes PG(N,p) (in decreasing order)
%
%  This routine finds the global mean and covariance matrix of a Gaussian Mixture (GMM). It also
%  attempts to find up to N local maxima using a combination of the fixed point and quadratic
%  Newton-Raphson algorithms from [1]. Currently, N must be less than or equal to the number of
%  mixtures K. In general the PDF surface of a GMM can be very complicated with many local maxima [2]
%  and, as discussed in [1,2], this algorithm is not guaranteed to find the N highest. In [2], it is
%  conjectured that the number of local maxima is <=K for the following cases (a) P=1, (b) all
%  mixture covariance matrices are equal and (c) all mixture covariance matrices are multiples of
%  the identity.
%
% Refs:
%   [1]	M. Á. Carreira-Perpiñán. Mode-finding for mixtures of gaussian distributions.
%       IEEE Trans. Pattern Anal and Machine Intell, 22 (11): 13181323, 2000. doi: 10.1109/34.888716.
%   [2] M. Á. Carreira-Perpiñán and C. K. I. Williams. On the number of modes of a gaussian mixture.
%       In Proc Intl Conf on Scale Space Theories in Computer Vision, volume LNCS 2695, pages 625640,
%       Isle of Skye, June 2003. doi: 10.1007/3-540-44935-3_44.

% Bugs/Suggestions:
% (1) Sometimes the mode is not found, e.g. m=[0 1; 1 0];v=[.01 10; 10 .01];
%     has a true mode near (0,0). Could add to the list of mode candidates
%     all the pairwise intersections of the mixtures.
%     Another is: m=[0 0; 10 0.3]; v=[1 1; 1000 .001];
% (2) When merging candidates, we should keep the one with the highest probability
% (3) could preserve the fixed arrays between calls if p and/or k are unchanged
%
% See also: gaussmix, gaussmixd, gaussmixp, randvec

%      Copyright (C) Mike Brookes 2000-2012
%      Version: $Id: gaussmixg.m 3227 2013-07-04 15:42:04Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Algorithm parameters
nfp=2;          % number of fixed point iterations to do at start
maxloop=60;     % maximum number of iterations
ssf=0.1;        % factor when calculating minimum mode separation
%Sort out arguments
[k,p]=size(m);
if nargin<4
    n=1;
    if nargin<3
        w=ones(k,1);
        if nargin<2
            v=ones(k,p);
        end
    end
end
if ~nargout
    if nargin<4
        n=k;
    end
    nao=4;
else
    nao=nargout;
end

full=numel(size(v))>2 || k==1 && numel(v)>p; % test for full covariance matrices
if full && p==1   % if p=1 then we force diagonal covariance matrices
    v=reshape(v,1,k)';
    full=false;
end
w=w/sum(w);  % force w to sum to unity
% calculate the global mean and covariance
mg=w.'*m;
mz=m-mg(ones(k,1),:);   % means relative to global mean
if nao>2                % need to calculate the modes
    nx=k;               % number of pg values initially
    kk=reshape(repmat(1:nx,k,1),k*nx,1); % [1 1 1 2 2 2 ... nx nx nx]'
    km=repmat(1:k,1,nx)'; % [1:k 1:k ... 1:k]'
    % sort out indexing for all data value pairs; needed to eliminate duplicates
    nxp=nx*(nx-1)/2;
    ja=(1:nxp)';
    jb=floor((3+sqrt(8*ja-3))/2);       % [2 3 3 4 4 4 ... nx]'
    ja=ja-(jb.^2-3*jb+2)/2;             % [1 1:2 1:3 ... 1:nx-1]'
    jc=ones(nxp,1);
    % sort out indexing for vectorized upper triagular matrix
    npu=p*(p+1)/2;                      % number of distinct elements in a symmetrial (p,p) matrix
    kw=1:npu;
    ku=floor((1+sqrt(8*kw-3))/2);       % [1 2 2 3 3 3 ... p]
    kv=kw-(ku.^2-ku)/2;                 % [1 1:2 1:3 ... 1:p]
    zpp=zeros(p,p);
    zpp(kv+p*(ku-1))=kw;
    zpp(ku+p*(kv-1))=kw;
    kw=reshape(zpp,1,[]);               % maps vectorized upper triangular to vectorized full matrix
    kp=repmat(1:p,1,p);                 % row indices for a (p,p) matrix as a row vector
    kq=reshape(repmat(1:p,p,1),1,p^2);  % col indices for a (p,p) matrix as a row vector
    kr=p*kp-p+kq;                       % transpose indexing for a vectorized (p,p) matrix
    kd=1:p+1:p^2;                       % diagonal indices of a (p,p) matrix
    % unity vectors to make efficient replication
    wk=ones(k,1);
    wnx=ones(nx,1);
    if full
        vg=mz.'*(mz.*w(:,ones(1,p)))+reshape(reshape(v,p^2,k)*w,p,p);
        % now determine the mode
        vi=zeros(p*k,p);                    % stack of k inverse cov matrices each size p*p times -0.5
        vim=zeros(p*k,1);                   % stack of k vectors of the form -0.5*inv(vt)*m
        mtk=vim;                             % stack of k vectors of the form m
        lvm=zeros(k,1);
        wpk=repmat((1:p)',k,1);
        for i=1:k    % loop for each mixture
            [uvk,dvk]=eig(v(:,:,i));      % find eigenvalues
            dvk=diag(dvk);
            if(any(dvk<=0))
                error('Covariance matrix for mixture %d is not positive definite',i);
            end
            vik=-0.5*uvk*diag(dvk.^(-1))*uvk';   % calculate inverse including -0.5 factor
            vi((i-1)*p+(1:p),:)=vik;           % vi contains all mixture inverses stacked on top of each other
            vim((i-1)*p+(1:p))=vik*m(i,:)';   % vim contains vi*m for all mixtures stacked on top of each other
            mtk((i-1)*p+(1:p))=m(i,:)';       % mtk contains all mixture means stacked on top of each other
            lvm(i)=log(w(i))-0.5*sum(log(dvk));       % vm contains the weighted sqrt of det(vi) for each mixture
        end
        vif=reshape(permute(reshape(vi,p,k,p),[2 1 3]),k,p^2); % each covariance matrix as a vectorized row
        vimf=reshape(vim,p,k)'; % vi*m as a row for each mixture
        ss=sqrt(min(v(repmat(kd,k,1)+repmat(p^2*(0:k-1)',1,p)),[],1))*ssf/sqrt(p);   % minimum separation of modes [this is a conservative guess]
    else
        vg=mz.'*(mz.*w(:,ones(1,p)))+diag(w.'*v);
        % now determine the mode
        vi=-0.5*v.^(-1);                % vi(k,p) = data-independent scale factor in exponent
        vi2=vi(:,ku).*vi(:,kv);         % vi2(k,npu) = upper triangular Hessian data dependent term
        lvm=log(w)-0.5*sum(log(v),2);   % log of external scale factor (excluding -0.5*p*log(2pi) term)
        vim=vi.*m;
        vim2=vim(:,ku).*vim(:,kv);      % vim2(k,npu) = upper triangular Hessian data independent term
        vimvi=vim(:,kp).*vi(:,kq);      % vimvi(k,p^2) = vectorized Hessian term
        ss=sqrt(min(v,[],1))*ssf/sqrt(p);   % minimum separation of modes [this is a conservative guess]
    end
    pgf=zeros(nx,p);                    % space for fixed point update
    sv=0.01*ss;                         % convergence threshold
    pg=m;  % initialize mode candidates to mixture means
    i=1;   % loop counter
    %     gx=zeros(nx,p);   %%%%%%%%%%%% temp
    while i<=maxloop
        %         pg00=pg0;   %%%%%%%%%%%% temp
        pg0=pg;                         % save previous mode candidates pg(nx,p)
        if full
            py=reshape(sum(reshape((vi*pg'-vim(:,wnx)).*(pg(:,wpk)'-mtk(:,wnx)),p,nx*k),1),k,nx)+lvm(:,wnx);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=(px./ps(wk,:))';             % px(nx,k) = relative mixture probabilities for each data point (rows sum to 1)
            % calculate the fixed point update
            pxvif=px*vif;    % pxvif(nx,p^2)
            pxvimf=px*vimf;  % pxvimf(nx,p)
            for j=1:nx
                pgf(j,:)=pxvimf(j,:)/reshape(pxvif(j,:),p,p);
            end
        else
            py=reshape(sum((pg(kk,:)-m(km,:)).^2.*vi(km,:),2),k,nx)+lvm(:,wnx);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=(px./ps(wk,:))';             % px(nx,k) = relative mixture probabilities for each data point (rows sum to 1)
            % calculate the fixed point update
            pxvim=px*vim;
            pxvi=px*vi;
            pgf=pxvim./pxvi;       % fixed point update for all points
        end
        if i>nfp
            % calculate gradient and Hessian; see [1] equations (4) and (5)
            lp=log(ps)+mx;              % log prob of each data point
            if full
                %                 gx0=gx;   %%%%%%%%%%%% temp
                gx=pxvimf-reshape(sum(repmat(pg,p,1).*reshape(pxvif,[],p),2),nx,p);
                vimpg=repmat(vimf,nx,1)-reshape(permute(reshape(pg*vi',nx,p,k),[3 1 2]),[],p); % vimpg(k*nx,p)
                hx1=2*reshape(sum(reshape(repmat(reshape(px',[],1),1,npu).*vimpg(:,ku).*vimpg(:,kv),k,[]),1),nx,[]);
                hx=pxvif+hx1(:,kw);
            else
                gx=pxvim-pxvi.*pg;               % gradient for each data point (one row per point)
                hx1=px*vim2+(px*vi2).*pg(:,ku).*pg(:,kv);
                hx2=(px*(vimvi)).*pg(:,kq);
                hx=2*(hx1(:,kw)-hx2-hx2(:,kr));
                hx(:,kd)=hx(:,kd)+pxvi;
            end
            hx=reshape(hx',p,p,nx);
            for j=1:nx
                if all(eig(hx(:,:,j))<0)  % if positive definite
                    pg(j,:)=pg(j,:)+gx(j,:)/hx(:,:,j); % do a Newton-Raphson update
                    if full
                        pyj=sum(reshape((vi*pg(j,:)'-vim).*(pg(j,wpk)'-mtk),p,k),1)'+lvm;
                    else
                        pyj=sum((repmat(pg(j,:),k,1)-m).^2.*vi,2)+lvm;
                    end
                    mxj=max(pyj);                % find normalizing factor for each data point to prevent underflow when using exp()
                    pxj=exp(pyj-mxj);            % find normalized probability of each mixture for each datapoint
                    psj=sum(pxj,1);                   % total normalized likelihood of each data point
                    lpj=log(psj)+mxj;              % log prob of updated data point
                    if lpj<lp(j)       % check if the probability has decreased
                        pg(j,:)=pgf(j,:);   % if so, do fixed point update
                    end
                else
                    pg(j,:)=pgf(j,:);   % else do fixed point update
                end
            end
        else
            pg=pgf;       % fixed point update for all points
        end
        if all(all(abs(pg-pg0)<sv(wnx,:))) && i+2<maxloop
            maxloop=min(maxloop,i+2);        % two more loops if converged if converged
        end
        %         [all(pg==pgf,2) pg [i; repmat(NaN,nx-1,1)] pg-pg0]  %debug: [fixed-point x-mode iteration delta-x]
        jd=all(abs(pg(jb,:)-pg(ja,:))<ss(jc,:),2);   % find duplicate modes
        if any(jd)
            jx=sparse([(1:nx)';ja;jb],[(1:nx)';jb;ja],[wnx;jd;jd]);   % neighbour matrix
            kx=any((jx*jx)>0 & ~jx,2);  % find chains that  are not fully connected
            while any(kx)
                kx=any(jx(:,kx),2);     % destroy all links connected to these chains
                jx(kx,:)=0;
                jx(:,kx)=0;
                % jx(kx,kx)=1;
                kx=any((jx*jx)>0 & ~jx,2);
            end
            jx([1:nx+1:nx*nx (ja+(jb-1)*nx)'])=0; % reset the upper triangle + diagonal
            pg(any(jx,2),:)=[];   % delete the duplicates
            % update nx and anything that depends on it
            nx=size(pg,1);
            pgf=zeros(nx,p);                    % space for fixed point update
            wnx=ones(nx,1);
            nxp=nx*(nx-1)/2;
            ja=ja(1:nxp);
            jb=jb(1:nxp);
            kk=reshape(repmat(1:nx,k,1),k*nx,1); % [1 1 1 2 2 2 ... nx nx nx]'
            km=reshape(repmat(1:k,1,nx),k*nx,1); % [1:k 1:k ... 1:k]'
            jc=ones(nxp,1);
        end
        i=i+1;
    end
    %     calculate the log pdf at each mode
    if full
        py=reshape(sum(reshape((vi*pg'-vim(:,wnx)).*(pg(:,wpk)'-mtk(:,wnx)),p,nx*k),1),k,nx)+lvm(:,wnx);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
    else
        py=reshape(sum((pg(kk,:)-m(km,:)).^2.*vi(km,:),2),k,nx)+lvm(:,wnx);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
    end
    [pv,ix]=sort((log(ps)+mx)'-0.5*p*log(2*pi),'descend');
    pg=pg(ix,:);
    if n<numel(pv) % only keep the first n modes
        pg=pg(1:n,:);
        pv=pv(1:n);
    end
elseif nao>1
    if full
        vg=mz.'*(mz.*w(:,ones(1,p)))+reshape(reshape(v,p^2,k)*w,p,p);
    else
        vg=mz.'*(mz.*w(:,ones(1,p)))+diag(w.'*v);
    end
end

if ~nargout
    % now plot the result
    clf;
    pg1=pg(1,:);
    lpm=gaussmixp(mg,m,v,w);
    lpp=pv(1);
    switch p
        case 1
            gaussmixp([],m,v,w);
            hold on;
            ylim=get(gca,'ylim')';
            plot([mg mg]',ylim,'-k',[mg mg; mg mg]+[-1 1; -1 1]*sqrt(vg),ylim,':k');
            plot(pg(1),pv(1),'^k');
            if numel(pg)>1
                plot(pg(2:end),pv(2:end),'xk');
            end
            hold off;
            title(sprintf('Mean+-sd = %.3g+-%.3g LogP = %.3g, Mode\\Delta = %.3g LogP = %.3g',mg,sqrt(vg),lpm,pg1,lpp));
            xlabel('x');
        case 2
            gaussmixp([],m,v,w);
            hold on;
            t=linspace(0,2*pi,100);
            xysd=chol(vg)'*[cos(t); sin(t)]+repmat(mg',1,length(t));
            plot(xysd(1,:),xysd(2,:),':k',mg(1),mg(2),'ok');
            plot(pg(1,1),pg(1,2),'^k');
            if numel(pv)>1
                plot(pg(2:end,1),pg(2:end,2),'xk');
            end
            hold off;
            title(sprintf('Mean = (%.3g,%.3g) LogP = %.3g, Mode:\\Delta = (%.3g,%.3g) LogP = %.3g',mg,lpm,pg1,pv(1)));
            xlabel('x');
            ylabel('y');
        otherwise
            nx=200;
            nc=ceil(sqrt(p/2));
            nr=ceil(p/nc);
            sdx=sqrt(diag(vg))';  % std deviation
            minx=min([mg; pg],[],1)-1.5*sdx;
            maxx=max([mg; pg],[],1)+1.5*sdx;
            ix=2:p; % selected indices
            for i=1:p
                xi=linspace(minx(i),maxx(i),nx)';
                [mm,vm,wm]=gaussmixd(mg(ix),m,v,w,[],ix);
                ym=gaussmixp(xi,mm,vm,wm')+lpm-gaussmixp(mg(i),mm,vm,wm');
                [mp,vp,wp]=gaussmixd(pg(1,ix),m,v,w,[],ix);
                yp=gaussmixp(xi,mp,vp,wp')+lpp-gaussmixp(pg1(i),mp,vp,wp');
                subplot(nr,nc,i);
                plot(xi,ym,'-k',mg(i),lpm,'ok',xi,yp,':b',pg1(i),lpp,'^b');
                axisenlarge([-1 -1 -1 -1.05]);
                hold on
                plot([mg(i) mg(i)],get(gca,'ylim'),'-k');
                hold off
                xlabel(sprintf('x[%d], Mean = %.3g, Mode\\Delta = %.3g',i,mg(i),pg1(i)));
                if i==1
                    title(sprintf('Log Prob: Mean = %.3g, Mode\\Delta = %.3g',lpm,lpp));
                end
                ix(i)=i;
            end
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function theRank=rankTComposition(p,firstElMostSig)
%%RANKTCOMPOSITION Obtain the lexicographic order (the rank) of a particular
%                 composition starting at 0. A composition of n into t
%                 is a way of putting n unlabeled balls into t labeled
%                 slots. Compositions are tuples of t integers >0 that
%                 sum to n. Compositions are given in lexicographic order
%                 where the first element is the least significant unless
%                 otherwise specified.
%
%INPUTS: p An tX1 vector holding the current composition, whose elements
%          sum to n>=1.
% firstElMostSig An optional parameter specifying whether the first element
%          is the most or least significant. The default if omitted or an
%          empty matrix is passed is false (the first element is the least
%          significant).
%
%OUTPUTS: theRank The rank of the composition in a lexicographic ordering
%                 of compositions, counting from zero. There is a total of
%                 binomial(n-1,t-1) compositions. Lexicographic ordering is
%                 defined using p(1) as the least significant element.
%
%This is an implementation of the relation described in Chapter 7.2.1.3 of
%[1] for mapping combinations to compositions. However, a special case is
%if the input is one-dimensional. In such an instance, the rank is always
%0.
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4, Fascicle 3:
%    Generating all Combinations and Partitions, Upper Saddle River, NJ:
%    Addison-Wesley, 2009.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Special case, empty matrix.
if(isempty(p))
   theRank=[];
   return;
end

%Special case, given a scalar.
if(length(p)==1)
    theRank=0;
   return;
end

if(nargin<2||isempty(firstElMostSig))
    firstElMostSig=false;
end

if(firstElMostSig)
   p=flipud(p(:));
end

m=length(p);

%Turn the composition into the equivalent combination.
mCombo=m-1;

c=zeros(mCombo,1);
c(1)=p(1)-1;

for k=2:mCombo
    c(k)=p(k)+c(k-1);
end

%get the rank of the combination.
theRank=rankCombination(c);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y = LPFilter(x,fc)
%
% y = LPFilter(x,fc),
% Overview: Second order zero-phase Lowpass filter. The filter is used to
% estimate the baseline wander in the input ECG signal.
%
% inputs:
% x: vector or matrix of input ECG data (channels x samples)
% fc: -3dB cut-off frequency normalized by the sampling frequency
%
% output:
% y: vector or matrix of filtered ECG data (channels x samples)
%
%
% Open Source ECG Toolbox, version 1.0, November 2006
% Released under the GNU General Public License
% Copyright (C) 2006  Reza Sameni
% Sharif University of Technology, Tehran, Iran -- GIPSA-Lab, INPG, Grenoble, France
% reza.sameni@gmail.com

% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details. You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

k = .7;     % cut-off value
alpha = (1-k*cos(2*pi*fc)-sqrt(2*k*(1-cos(2*pi*fc))-k^2*sin(2*pi*fc)^2))/(1-k);
y = zeros(size(x));
for i = 1:size(x,1)
    y(i,:) = filtfilt(1-alpha,[1 -alpha],x(i,:));
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [Offspring] = DEgenerator2(Problem,Population)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    cv = sum(max(0,Population.cons),2);

    FrontNo = NDSort(Population.objs,Population.cons,1);
    index1  = find(FrontNo==1);
    r       = floor(rand*length(index1))+1;
    best    = index1(r);

    [N,D] = size(Population(1).decs);
    trial = zeros(1*Problem.N,D);

    for i = 1 : Problem.N
        % DE/rand-to-bestf/1/bin--Convergence
        if rand > 0.5
            l = rand;
            if l <= 1/3
            	F = .6;
            elseif l <= 2/3
                F = 0.8;
            else
                F = 1.0;
            end
            l = rand;
            if l <= 1/3
                CR = .1;
            elseif l <= 2/3
                CR = 0.2;
            else
                CR = 1.0;
            end
            indexset    = 1 : Problem.N;
            indexset(i) = [];
            r1  = floor(rand*(Problem.N-1))+1;
            xr1 = indexset(r1);
            indexset(r1) = [];
            r2  = floor(rand*(Problem.N-2))+1;
            xr2 = indexset(r2)  ;
            r3  = floor(rand*(Problem.N-3))+1;
            xr3 = indexset(r3);
            Best_index = Population(best).decs;
            v      = Population(xr1).decs+rand*(Best_index-Population(xr1).decs)+F*(Population(xr2).decs-Population(xr3).decs);
            Lower  = repmat(Problem.lower,N,1);
            Upper  = repmat(Problem.upper,N,1);
            v      = min(max(v,Lower),Upper);
            Site   = rand(N,D) < CR;
            j_rand = floor(rand * D) + 1;
            Site(1, j_rand) = 1;
            Site_  = 1-Site;
            trial(i, :) = Site.*v+Site_.*Population(i).decs;
        % DE/current-to-rand/1--Diversity
        else
            l = rand;
            if l <= 1/3
                F = .6;
            elseif l <= 2/3
                F = 0.8;
            else
                F = 1.0;
            end
            indexset    = 1:Problem.N;
            indexset(i) = [];
            r1  = floor(rand*(Problem.N-1))+1;
            xr1 = indexset(r1);
            indexset(r1) = [];
            r2  = floor(rand*(Problem.N-2))+1;
            xr2 = indexset(r2);
            indexset(r2) = [];
            r3    = floor(rand*(Problem.N-3))+1;
            xr3   = indexset(r3);
            v     = Population(i).decs+rand*(Population(xr1).decs-Population(i).decs)+F*(Population(xr2).decs-Population(xr3).decs);
            Lower = repmat(Problem.lower,N,1);
            Upper = repmat(Problem.upper,N,1);
            trial(i, :) = min(max(v,Lower),Upper);
        end
    end
    Offspring = trial;
    Offspring = Problem.Evaluation(Offspring);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function y=v_teager(x,d,m)
%V_TEAGER calculate v_teager energy waveform Y=(X,D,M)
%
%  Inputs:  x         speech signal
%           d         dimension to apply filter along [default 1st non-singleton]
%           m         Normally Y has the same length as X and the first
%                     and last output samples are extrapolated. Setting m='x'
%                     supresses this extrapolation and Y will be two
%                     samples shorter than X
%
% Outputs:  Y         output signal: y(n)=abs(x(n))^2 - x(n+1)*conj(x(n-1))
%
% Calculates the Teager energy waveform [1]. The following waveforms give
% a constant output (independent of n) where A, B, C are real constants:
%  (a) x(n) = A*sin(B*n+C)    -->   y(n) = (A*sin(B))^2
%  (b) x(n) = A*n + B         -->   y(n) = A^2
%  (c) x(n) = A*exp(j(B*n+C)) -->   y(n) = A^2*(1-exp(2jB))
%  (d) x(n) = A*exp(B*n+C)    -->   y(n) = 0
%
% Reference:
%  [1]	J. Kaiser. On a simple algorithm to calculate the "energy" of a signal.
%       In Proc IEEE Intl Conf Acoustics, Speech and Signal Processing,
%       pages 381-384, vol.1, Apr. 1990. doi: 10.1109/ICASSP.1990.115702.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_teager.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e=size(x);
p=prod(e);
if nargin<2             % if no dimension given, find the first non-singleton
    d=find(e>1,1);
    if ~numel(d)
        d=1;
    end
end
k=e(d);                 % size of active dimension
q=p/k;                  % size of remainder
if d==1
    z=reshape(x,k,q);
else
    z=shiftdim(x,d-1);
    r=size(z);
    z=reshape(z,k,q);
end
if nargin>2 && any(m=='x')
    y=z(2:k-1,:).*conj(z(2:k-1,:))-z(3:k,:).*conj(z(1:k-2,:));
    k=k-2;              % we have lost two elements
elseif k>=4
    y=zeros(k,q);
    y(2:k-1,:)=z(2:k-1,:).*conj(z(2:k-1,:))-z(3:k,:).*conj(z(1:k-2,:));
    y(1,:)=2*y(2,:)-y(3,:);             % linearly interpolate the end points
    y(k,:)=2*y(k-1,:)-y(k-2,:);
elseif k==3
    y=repmat(x(2,:).*conj(x(2,:))-x(3,:).*conj(x(1,:)),3,1);
else
    y=zeros(k,q);
end
if d==1
    e(d)=k;
    y=reshape(y,e);
else
    r(1)=k;
    y=shiftdim(reshape(y,r),length(e)+1-d);
end

author: ImperialCollegeLondon

repo: sap-voicebox