text:
function [ir,x0,delay] = ir_nfchoa(X,head_orientation,xs,src,sofa,conf)
%IR_NFCHOA binaural simulation of NFCHOA
%
%   Usage: [ir,x0,delay] = ir_nfchoa(X,head_orientation,xs,src,sofa,conf)
%
%   Input parameters:
%       X                - listener position / m
%       head_orientation - orientation of the listener with [phi theta] /
%                          (rad, rad)
%       xs               - virtual source position / m
%       src              - source type: 'pw' -plane wave
%                                       'ps' - point source
%       sofa             - impulse response data set for the secondary sources
%       conf             - configuration struct (see SFS_config)
%
%   Output parameters:
%       ir               - impulse response for the desired HOA synthesis
%                          (nx2 matrix)
%       x0               - secondary sources
%       delay            - delay added by driving function / s
%
%   IR_NFCHOA(X,head_orientation,xs,src,L,sofa,conf) calculates a binaural room
%   impulse response for a virtual source at xs for a virtual NFCHOA array and
%   a listener located at X.
%
%   See also: ssr_brs_nfchoa, ir_nfchoa, ir_point_source, auralize_ir

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 6;
nargmax = 6;
narginchk(nargmin,nargmax);
if conf.debug
    isargposition(X);
    isargxs(xs);
    isargvector(head_orientation);
    isargpositivescalar(L);
    isargchar(src);
    isargstruct(config);
end

%% ===== Variables ======================================================
% Loudspeaker positions
x0 = secondary_source_positions(conf);

%% ===== BRIR ===========================================================
% Calculate driving function
[d,~,delay] = driving_function_imp_nfchoa(x0,xs,src,conf);
% Generate the impulse response for NFCHOA
ir = ir_generic(X,head_orientation,x0,d,sofa,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function [L, stored_at] = quadtree_laplacian(C,W,CH,D,A)
% QUADTREE_LAPLACIAN
% Builds a finite difference Laplacian on a quadtree following the scheme
% suggested by Bickel et al. "Adaptative Simulation of Electrical
% Discharges". This code is *purposefully* not optimized beyond
% asymptotics for simplicity in understanding its functionality and
% translating it to other programming languages beyond prototyping.
%
% L = octree_laplacian(C,W,CH,D,A)
%
% Inputs:
%   C #nodes by 3 matrix of cell centers
%   W #nodes vector of cell widths (**not** half widths)
%   CH #nodes by 4 matrix of child indeces (-1 if leaf node)
%   D #nodes vector of tree depths
%   A #nodes by #nodes sparse adjacency matrix, where a value of a in the
%       (i,j) entry means that node j is to the a-th direction of i
%       (a=1: left;  a=2: right;  a=3: bottom;  a=4: top).
%
% Outputs:
%   L #num_children by #num_children sparse Laplacian matrix
%   stored_at #num_children by 3 matrix of child cell centers, where the
%       values of L are stored
%
%
% Example:
%
% % Build an octree
% P = 0.5*[cos(th),sin(th)];
% P = [P;[-1,-1];[1,1]];
% [C,W,CH,PAR,D,A] = initialize_quadtree(P,'MaxDepth',8,'Graded',true);
% % This is for plotting
% [V,Q] = bad_quad_mesh_from_quadtree(C,W,CH);
% % Call function to construct Laplacian
% [L, stored_at] = quadtree_laplacian(C,W,CH,D,A);
% % Dummy Laplacian function
% gt_fun = stored_at(:,1).^2.0;
% laplacian_fun = 2 + 0.*stored_at(:,1);
% % Find boundary
% Vmin = min(stored_at,[],1);
% Vmax = max(stored_at,[],1);
% is_boundary = (stored_at(:,1)<=(Vmin(1)+0.2)) + (stored_at(:,2)<=(Vmin(2)+0.2)) + ...
%     (stored_at(:,1)>=(Vmax(1)-0.2)) + (stored_at(:,2)>=(Vmax(2)-0.2));
% bb = find(is_boundary);
% bc = gt_fun(bb);
% % Solve as energy minimization
% u = min_quad_with_fixed(0.5*L,-laplacian_fun,bb,bc);
% % Plot solution
% tsurf(Q,V,falpha(1,1),'FaceVertexCData',u)
% hold on
% sct(stored_at(bb,:)) % Visualize boundary conditions
% set(gcf,'Color','w'); grid off; axis equal; colorbar; caxis([0 1]);
%
%
% See also: initialize_quadtree.m

% We will store Laplacian values at
% child cell indeces
children = find(CH(:,1)==-1);
% map from all cells to children
cell_to_children = -ones(size(W,1),1);
cell_to_children(children) = 1:length(children);

% Vectors for constructing the Laplacian
I = [];
J = [];
vals = [];

for i=1:size(children,1)
    new_I = [];
    new_J = [];
    new_vals = [];
    l = [1,1,1,1];
    new_dirs = [];
    child = children(i);
    d = D(child);
    num_dirs = 0;
    % Let's build d u(child)/dx^2 ~ u(child+W(child)*[1,0])/hr(hl+hr) -
    % 2u(child)/hlhr + u(child-W(child)*[1,0])/hr(hl+hr)
    % So, let's look for the value to the j direction. To do this, we seek the
    % lowest-depth neighbor to the j direction. As a reminder the octree
    % adjacency convention is i->j (1:left-2:right-3:bottom-4:top)
    for j=1:4
        j_neighbors = find(A(child,:)==j);
        if ~isempty(j_neighbors)
            depths_j_neighbors = D(j_neighbors);
            [max_depth_j, max_depth_j_neighbor] = max(depths_j_neighbors);
            max_depth_j_neighbor = j_neighbors(max_depth_j_neighbor);
            % There are two options:
            % One: the leaf node to our j direction has lower or equal depth to
            % us
            if max_depth_j<=d
                l(j) = (W(child) + W(max_depth_j_neighbor))/2;
                % then it's easy, just add this node
                new_I = [new_I;i];
                % THIS HAS TO BE A CHILD !
                assert(cell_to_children(max_depth_j_neighbor)>0);
                new_J = [new_J;cell_to_children(max_depth_j_neighbor)];
                new_vals = [new_vals;-1]; % Todo fix this
                new_dirs = [new_dirs;j];
            else
                % In this case, assuming the grid is graded, there should
                % be two j-neighbors at depth d+1
                nn = j_neighbors(D(j_neighbors)==(d+1));
                assert(length(nn)==2,"Are you sure you are inputting a graded quadtree?")
                assert(all(CH(nn,1)==[-1;-1]))
                % Then we simply average both
                l(j) = (W(child) + W(nn(1)))/2;
                new_I = [new_I;i;i];
                new_J = [new_J;cell_to_children(nn(1));cell_to_children(nn(2))];
                new_vals = [new_vals;-0.5;-0.5];
                new_dirs = [new_dirs;j;j];
            end
            num_dirs = num_dirs + 1;
        end
    end
    new_I = [new_I;i];
    new_J = [new_J;i];
    new_vals = [new_vals; 1.0];
    new_dirs = [new_dirs;5]; % just a hack
    % At this point, we have to divide by the edge-lengths
    new_vals(new_dirs==1) = new_vals(new_dirs==1)/(l(1)*(l(1)+l(2)));
    new_vals(new_dirs==2) = new_vals(new_dirs==2)/(l(2)*(l(1)+l(2)));
    new_vals(new_dirs==3) = new_vals(new_dirs==3)/(l(3)*(l(3)+l(4)));
    new_vals(new_dirs==4) = new_vals(new_dirs==4)/(l(4)*(l(3)+l(4)));
    new_vals(new_dirs==5) = 1/(l(1)*l(2)) + 1/(l(3)*l(4));

    % And add them to the big sparse Laplacian construction vectors
    I = [I;new_I];
    J = [J;new_J];
    vals = [vals;new_vals];
end

% THE LAPLACIAN IS NEGATIVE SEMI DEFINITE!
L = -2*sparse(I,J,vals,length(children),length(children));
stored_at = C(children,:);

end

author: alecjacobson

repo: gptoolbox


text:
function [x]=gdist(a,x)
%[Y] = GDIST(A, X) Guitar Distortion
%
%   GDIST creates a distortion effect like that of
%   an overdriven guitar amplifier. This is a Matlab
%   implementation of an algorithm that was found on
%   www.musicdsp.org.
%
%   A = The amount of distortion.  A
%       should be chosen so that -1<A<1.
%   X = Input.  Should be a column vector
%       between -1 and 1.
%
%coded by: Steve McGovern, date: 09.29.04
%URL: http://www.steve-m.us

k = 2*a/(1-a);
x = (1+k)*(x)./(1+k*abs(x));

author: Sable

repo: mcbench-benchmarks


text:
function [HDR]=openxml(arg1,CHAN,arg4,arg5,arg6)
% OPENXML reads XML files and tries to extract biosignal data
%
% This is an auxilary function to SOPEN.
% Use SOPEN instead of OPENXML.
%

%
% HDR = openxml(HDR);
%
% HDR contains the Headerinformation and internal data
%
% see also: SOPEN, SREAD, SSEEK, STELL, SCLOSE, SWRITE, SEOF

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

%	$Id$
%	Copyright 2006,2007,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

if ischar(arg1);
        HDR.FileName = arg1;
        HDR.FILE.PERMISSION = 'r';
else
        HDR = arg1;
end;

%if strncmp(HDR.TYPE,'XML',3),
%        if any(HDR.FILE.PERMISSION=='r'),
                fid = fopen(HDR.FileName,HDR.FILE.PERMISSION,'ieee-le');
                s = char(fread(fid,[1,1024],'char'));
                if all(s(1:2)==[255,254]) & all(s(4:2:end)==0)
                        HDR.TYPE='XML-UTF16';
                elseif ~isempty(findstr(char(s),'?xml version'))
                        HDR.TYPE='XML-UTF8';
		end;
		fseek(fid,0,'bof');
                if strcmp(HDR.TYPE,'XML-UTF16'),
                        magic = char(fread(fid,1,'uint16'));
                        HDR.XML = char(fread(fid,[1,inf],'uint16'));
                elseif strcmp(HDR.TYPE,'XML-UTF8'),
                        HDR.XML = char(fread(fid,[1,inf],'uint8'));
                end;
                fclose(fid);
                HDR.FILE.FID = fid;

                if 0, ~exist('xmlstruct')
                        warning('XML toolbox missing')
                end;
                if 1,
                        HDR.XMLstruct = xmlstruct(HDR.XML,'sub');
                        HDR.XMLlist   = xmlstruct(HDR.XML);
                end;

                try
                        XML = xmltree(HDR.XML);
                        XML = convert(XML);
                        HDR.XML  =  XML;
			HDR.TYPE = 'XML';
                catch
                        fprintf(HDR.FILE.stderr,'ERROR SOPEN (XML): XML-toolbox missing or invalid XML file.\n');
                        return;
                end;

                tmp = fieldnames(HDR.XML);
                if any(strmatch('PatientDemographics',tmp)) & any(strmatch('TestDemographics',tmp)) & any(strmatch('RestingECGMeasurements',tmp)) & any(strmatch('Diagnosis',tmp)) & any(strmatch('Waveform',tmp)),
                        % GE-Marquette FDA-XML   MAC5000

                        HDR.Patient.ID = HDR.XML.PatientDemographics.PatientID;
                        tmp = HDR.XML.PatientDemographics.Gender;
                        HDR.Patient.Sex = 2*(upper(tmp(1))=='F') + (upper(tmp(1))=='M');
                        HDR.Patient.Name = [HDR.XML.PatientDemographics.PatientLastName,', ',HDR.XML.PatientDemographics.PatientFirstName];

                        tmp = HDR.XML.TestDemographics.AcquisitionDate;
                        tmp(tmp=='-') = ' ';
                        HDR.T0([3,2,1])=str2double(tmp);
                        tmp = HDR.XML.TestDemographics.AcquisitionTime;
                        tmp(tmp==':') = ' ';
                        HDR.T0(4:6) = str2double(tmp);

                        HDR.NS = str2double(HDR.XML.Waveform.NumberofLeads);
                        HDR.SampleRate = str2double(HDR.XML.Waveform.SampleBase);
                        HDR.Filter.LowPass  = str2double(HDR.XML.Waveform.LowPassFilter)*ones(1,HDR.NS);
                        HDR.Filter.HighPass = str2double(HDR.XML.Waveform.HighPassFilter)*ones(1,HDR.NS);
                        HDR.Filter.Notch    = str2double(HDR.XML.Waveform.ACFilter)*ones(1,HDR.NS);

                        HDR.NRec = 1;
                        HDR.SPR = 1;
                        for k = 1:HDR.NS,
                                CH = HDR.XML.Waveform.LeadData{k};
                                HDR.AS.SPR(k)  = str2double(CH.LeadSampleCountTotal);
                                HDR.SPR        = lcm(HDR.SPR,HDR.AS.SPR(k));
                                HDR.Cal(k)     = str2double(CH.LeadAmplitudeUnitsPerBit);
                                HDR.PhysDim{k} = CH.LeadAmplitudeUnits;
                                HDR.Label{k}   = CH.LeadID;

                                t = radix64d(CH.WaveFormData);
                                t = 256*t(2:2:end) + t(1:2:end);
                                t = t - (t>=(2^15))*(2^16);
                                HDR.data(:,k) = t(:);
                        end;
                        HDR.TYPE = 'native';

                elseif any(strmatch('component',tmp))
                          % FDA-XML Format
                        tmp   = HDR.XML.component.series.derivation;
                        if isfield(tmp,'Series');
                                tmp = tmp.Series.component.sequenceSet.component;
                        else    % Dovermed.CO.IL version of format
                                tmp = tmp.derivedSeries.component.sequenceSet.component;
                        end;
                        HDR.NS = length(tmp)-1;
                        HDR.NRec = 1;
                        HDR.Cal = 1;
                        HDR.PhysDim = {' '};
                        HDR.SampleRate = 1;
                        HDR.TYPE = 'XML-FDA';     % that's an FDA XML file

                elseif any(strmatch('dataacquisition',tmp)) & any(strmatch('reportinfo',tmp)) & any(strmatch('patient',tmp)) & any(strmatch('documentinfo',tmp)),
                        % SierraECG  1.03  *.open.xml from PHILIPS
                        HDR.SampleRate = str2double(HDR.XML.dataacquisition.signalcharacteristics.samplingrate);
                        HDR.NS  = str2double(HDR.XML.dataacquisition.signalcharacteristics.numberchannelsvalid);
                        HDR.Cal = str2double(HDR.XML.reportinfo.reportgain.amplitudegain.overallgain);
                        HDR.PhysDim = {'uV'};
                        HDR.Filter.HighPass = str2double(HDR.XML.reportinfo.reportbandwidth.highpassfiltersetting);
                        HDR.Filter.LowPass  = str2double(HDR.XML.reportinfo.reportbandwidth.lowpassfiltersetting);
                        HDR.Filter.Notch    = str2double(HDR.XML.reportinfo.reportbandwidth.notchfiltersetting);

                        t = HDR.XML.reportinfo.reportformat.waveformformat.mainwaveformformat;
                        k = 0;
                        HDR.Label={};
                        while ~isempty(t),
                                [s,t] = strtok(t,' ');
                                k = k+1;
                                HDR.Label{k, 1} = [s,' '];
                        end;
                        HDR.Patient.Id     = str2double(HDR.XML.patient.generalpatientdata.patientid);
                        tmp    = HDR.XML.patient.generalpatientdata.age;
                        if isfield(tmp,'years'),
                                HDR.Patient.Age    = str2double(tmp.years);
                        end
                        if isfield(tmp,'dateofbirth')
                                tmp = tmp.dateofbirth;
                                tmp(tmp=='-')=' ';
                                HDR.Patient.Birthday([6,5,4]) = str2double(tmp);
                        end;

                        tmp    = HDR.XML.patient.generalpatientdata.sex;
                        HDR.Patient.Sex    = strncmpi(tmp,'Male',1) + strncmpi(tmp,'Female',1)*2;
                        HDR.Patient.Weight = str2double(HDR.XML.patient.generalpatientdata.weight.kg);
                        HDR.Patient.Height = str2double(HDR.XML.patient.generalpatientdata.height.cm);

                        HDR.VERSION = HDR.XML.documentinfo.documentversion;
                        HDR.TYPE = HDR.XML.documentinfo.documenttype;

                elseif any(strmatch('component',tmp)) & any(strmatch('reportinfo',tmp)) & any(strmatch('patient',tmp)) & any(strmatch('documentinfo',tmp)),
                           % FDA-XML Format
                        tmp   = HDR.XML.component.series.derivation;
                        if isfield(tmp,'Series');
                                tmp = tmp.Series.component.sequenceSet.component;
                        else    % Dovermed.CO.IL version of format
                                tmp = tmp.derivedSeries.component.sequenceSet.component;
                        end;
                        HDR.NS = length(tmp)-1;
                        HDR.NRec = 1;
                        HDR.Cal = 1;
                        HDR.PhysDim = {' '};
                        HDR.SampleRate = 1;
                        HDR.TYPE = 'XML-FDA';     % that's an FDA XML file

                elseif any(strmatch('StripData',tmp)) & any(strmatch('ClinicalInfo',tmp)) & any(strmatch('PatientInfo',tmp)) & any(strmatch('ArrhythmiaData',tmp)),
                        % GE Case8000 stress ECG

                	HDR.SampleRate = str2double(HDR.XML.StripData.SampleRate);
                	tmp = HDR.XML.ClinicalInfo.ObservationDateTime;
                	HDR.T0 = [str2double(tmp.Year), str2double(tmp.Month), str2double(tmp.Day), str2double(tmp.Hour), str2double(tmp.Minute), str2double(tmp.Second)];

                	HDR.Patient.Id = HDR.XML.PatientInfo.PID;
                	HDR.Patient.Name = 'X'; % [HDR.XML.PatientInfo.Name, ', ', HDR.XML.PatientInfo.GivenName];
                	HDR.Patient.Age = str2double(HDR.XML.PatientInfo.Age);
                	tmp = HDR.XML.PatientInfo.Gender;
                	HDR.Patient.Sex = any(tmp(1)=='Mm') + any(tmp(1)=='Ff')*2;
                	HDR.Patient.Height = str2double(HDR.XML.PatientInfo.Height);
                	HDR.Patient.Weight = str2double(HDR.XML.PatientInfo.Weight);
                	tmp = HDR.XML.PatientInfo.BirthDateTime;
                	HDR.Patient.Birthday = [str2double(tmp.Year), str2double(tmp.Month), str2double(tmp.Day),0,0,0];

                	tmp = HDR.XML.StripData.Strip;
                	HDR.NS = length(tmp{1}.WaveformData);
               		tmax   = str2double(tmp{end}.Time.Minute)*60 + str2double(tmp{end}.Time.Second)+10;
                	data   = repmat(NaN,tmax*HDR.SampleRate,HDR.NS);
                	for k  = 1:length(tmp);
	                	t = HDR.SampleRate*(str2double(tmp{k}.Time.Minute)*60 + str2double(tmp{k}.Time.Second));
	                	for k2 = 1:HDR.NS,
		                	x = str2double(tmp{k}.WaveformData{k2});
		                	data(t+1:t+length(x),k2)=x(:);
	                	end;
                	end;
                	tmp = HDR.XML.ArrhythmiaData.Strip;
                	for k  = 1:length(tmp);
	                	t = HDR.SampleRate*(str2double(tmp{k}.Time.Minute)*60 + str2double(tmp{k}.Time.Second));
	                	for k2 = 1:HDR.NS,
		                	x = str2double(tmp{k}.WaveformData{k2});
		                	data(t+1:t+length(x),k2)=x(:);
	                	end;
                	end;
                	HDR.data = data - 2^12*(data>2^11);
                	HDR.TYPE = 'native';
                	HDR.NRec = 1;
                	HDR.SPR  = size(HDR.data,1);
                	HDR.Calib = sparse(2:HDR.NS,1:HDR.NS,1);
                	HDR.FLAG.UCAL = 1;

                else
                        fprintf(HDR.FILE.stderr,'Warning SOPEN (XML): File %s is not supported.\n',HDR.FileName);
                        return;
                end

                try
                	tmp=HDR.XML.componentOf.timepointEvent.componentOf.subjectAssignment.subject.trialSubject.subjectDemographicPerson.name;
                	HDR.Patient.Name = sprintf('%s, %s',tmp.family, tmp.given);
		catch
                end;

                HDR.Calib = sparse(2:HDR.NS+1,1:HDR.NS,HDR.Cal);
                HDR.FILE.OPEN = 1;
                HDR.FILE.POS  = 0;
%        end;
%end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Auxillary functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = radix64d(x);
% RADIX64D - decoding of radix64 encoded sequence

% 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

%	$Id$
%	(C) 2006 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

global BIOSIG_GLOBAL

if ~isfield(BIOSIG_GLOBAL,'R64E');
	BIOSIG_GLOBAL.R64E = ['A':'Z','a':'z','0':'9','+','/'];
	BIOSIG_GLOBAL.R64D = zeros(256,1)-1;
	for k = 1:length(BIOSIG_GLOBAL.R64E),
		BIOSIG_GLOBAL.R64D(BIOSIG_GLOBAL.R64E(k)) = k-1;
	end;
end;

% http://www.faqs.org/rfcs/rfc2440.html

t = BIOSIG_GLOBAL.R64D(x);
t = [t(t>=0); zeros(3,1)];
N = floor(length(t)/4);
t = reshape(bitand(t(1:N*4),2^6-1),4,N);

y(1,:) = bitshift(t(1,:),2)         + bitshift(t(2,:),-4);
y(2,:) = bitshift(mod(t(2,:),16),4) + bitshift(t(3,:),-2);
y(3,:) = bitshift(mod(t(3,:),4),6)  + t(4,:);

y = y(:)';
y = y(1:end-sum(x=='=')); % remove possible pad characters

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef FDA3 < PROBLEM
% <multi> <real> <large/none> <dynamic>
% Benchmark dynamic MOP proposed by Farina, Deb, and Amato
% taut --- 10 --- Number of generations for static optimization
% nt   --- 10 --- Number of distinct steps

%------------------------------- Reference --------------------------------
% M. Farina, K. Deb, and P. Amato, Dynamic multiobjective optimization
% problems: Test cases, approximations, and applications, IEEE Transactions
% on Evolutionary Computation, 2004, 8(5): 425-442.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        taut;       % Number of generations for static optimization
        nt;         % Number of distinct steps
        Optimums;   % Point sets on all Pareto fronts
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.taut,obj.nt] = obj.ParameterSet(10,10);
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = [0,-ones(1,obj.D-1)];
            obj.upper    = [1, ones(1,obj.D-1)];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate solutions
        function Population = Evaluation(obj,varargin)
            PopDec     = obj.CalDec(varargin{1});
            PopObj     = obj.CalObj(PopDec);
            PopCon     = obj.CalCon(PopDec);
            % Attach the current number of function evaluations to solutions
            Population = SOLUTION(PopDec,PopObj,PopCon,zeros(size(PopDec,1),1)+obj.FE);
            obj.FE     = obj.FE + length(Population);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            t = floor(obj.FE/obj.N/obj.taut)/obj.nt;
            PopObj(:,1) = PopDec(:,1).^(10.^(2*sin(0.5*pi*t)));
            G = abs(sin(0.5*pi*t));
            g = 1 + G + sum((PopDec(:,2:end)-G).^2,2);
            h = 1 - (PopObj(:,1)./g).^0.5;
            PopObj(:,2) = g.*h;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            % Generate point sets on all Pareto fronts
            t = floor(0:obj.maxFE/obj.N/obj.taut)/obj.nt;
            G = abs(sin(0.5.*pi.*t));
            G = unique(round(G*1e6)/1e6);
            x = linspace(0,1,N)';
            obj.Optimums = {};
            for i = 1 : length(G)
                g = 1 + G(i);
                obj.Optimums(i,:) = {G(i),[x,g*(1-sqrt(x/g))]};
            end
            % Combine all point sets
            R = cat(1,obj.Optimums{:,2});
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            t      = floor(Population.adds/obj.N/obj.taut)/obj.nt;
            G      = abs(sin(0.5.*pi.*t));
            G      = round(G*1e6)/1e6;
            change = [0;find(G(1:end-1)~=G(2:end));length(G)];
            Scores = zeros(1,length(change)-1);
            allG   = cell2mat(obj.Optimums(:,1));
            for i = 1 : length(change)-1
                subPop    = Population(change(i)+1:change(i+1));
                Scores(i) = feval(metName,subPop,obj.Optimums{find(G(change(i)+1)==allG,1),2});
            end
            score = mean(Scores);
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            t      = floor(Population.adds/obj.N/obj.taut)/obj.nt;
            G      = abs(sin(0.5.*pi.*t));
            G      = round(G*1e6)/1e6;
            change = [0;find(G(1:end-1)~=G(2:end));length(G)];
            allG   = cell2mat(obj.Optimums(:,1));
            tempStream = RandStream('mlfg6331_64','Seed',2);
            for i = 1 : length(change)-1
                color = rand(tempStream,1,3);
                Draw(Population(change(i)+1:change(i+1)).objs,'o','MarkerSize',5,'Marker','o','Markerfacecolor',sqrt(color),'Markeredgecolor',color,{'\it f\rm_1','\it f\rm_2',[]});
                Draw(obj.Optimums{find(G(change(i)+1)==allG,1),2},'-','LineWidth',1,'Color',color);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [model, L] = ppcaVb(X, q, prior)
% Perform variatioanl Bayeisan inference for probabilistic PCA model.
% Input:
%   X: d x n data matrix
%   q: dimension of target space
% Output:
%   model: trained model structure
%   L: variantional lower bound
% Reference:
%   Pattern Recognition and Machine Learning by Christopher M. Bishop
% Written by Mo Chen (sth4nth@gmail.com).
[m,n] = size(X);
if nargin < 3
    a0 = 1e-4;
    b0 = 1e-4;
    c0 = 1e-4;
    d0 = 1e-4;
else
    a0 = prior.a;
    b0 = prior.b;
    c0 = prior.c;
    d0 = prior.d;
end

if nargin < 2
    q = m-1;
end
tol = 1e-6;
maxIter = 500;
L = -inf(1,maxIter);

mu = mean(X,2);
Xo = bsxfun(@minus, X, mu);
s = dot(Xo(:),Xo(:));
I = eye(q);
% init parameters
a = a0+m/2;
c = c0+m*n/2;
Ealpha = 1e-4;
Ebeta = 1e-4;
EW = rand(q,m);
EWo = bsxfun(@minus,EW,mean(EW,2));
EWW = EWo*EWo'/m+EW*EW';
for iter = 2:maxIter
%     q(z)
    LZ = I+Ebeta*EWW;
    V = inv(chol(LZ));                   % inv(LZ) = V*V';
    EZ = LZ\EW*Xo*Ebeta;
    EZZ = n*(V*V')+EZ*EZ';
    KLZ = n*sum(log(diag(V)));           % KLZ = 0.5*n*log(det(inv(LZ)));
%     q(w)
    LW = diag(Ealpha)+Ebeta*EZZ;
    V = inv(chol(LW));                   % inv(LW) = V*V';
    EW = LW\EZ*Xo'*Ebeta;
    EWW = m*(V*V')+EW*EW';
    KLW = m*sum(log(diag(V)));           % KLW = 0.5*n*log(det(inv(LW)));
%     q(alpha)
    b = b0+diag(EWW)/2;
    Ealpha = a./b;
    KLalpha = -sum(a*log(b));
%     q(beta)
    WZ = EW'*EZ;
    d = d0+(s-2*dot(Xo(:),WZ(:))+dot(EWW(:),EZZ(:)))/2;
    Ebeta = c/d;
    KLbeta = -c*log(d);
%     q(mu)
%     Emu = Ebeta/(lambda+n*Ebeta)*sum(X-WZ,2);

%     lower bound
    L(iter) = KLalpha+KLbeta+KLW+KLZ;
    if L(iter)-L(iter-1) < tol*abs(L(iter-1)); break; end
end
L = L(2:iter);

model.Z = EZ;
model.W = EW;
model.apha = Ealpha;
model.beta = Ebeta;
model.a = a;
model.b = b;
model.c = c;
model.d = d;
model.mu = mu;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%   CLASS GUI_Remote_Resources
% =========================================================================
%
% DESCRIPTION
%   class to manage the about window of goGPSz
%
% EXAMPLE
%   ui = GUI_Remote_Resources.getInstance();
%
% FOR A LIST OF CONSTANTs and METHODS use doc Core_UI

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef GUI_Remote_Resources < GUI_Unique_Win

    properties (Constant)
        WIN_NAME = 'goGPS_Inspect_RR';
    end

    properties (Constant, Access = 'protected')
        BG_COLOR = Core_UI.DARK_GREY_BG;
    end

    %% PROPERTIES SINGLETON POINTERS
    % ==================================================================================================================================================
    properties % Utility Pointers to Singletons
        log
        state
    end

    %% PROPERTIES GUI
    % ==================================================================================================================================================
    properties
        w_edt       % Handle of the main window
        win         % Handle to this window
        j_rrini     % ini resources file java pointer
    end

    %% PROPERTIES STATUS
    % ==================================================================================================================================================
    properties (GetAccess = private, SetAccess = private)
    end

    %% METHOD CREATOR
    % ==================================================================================================================================================

    methods (Static, Access = private)
        function this = GUI_Remote_Resources(w_edt)
            % GUI_MAIN object creator
            this.init();
            this.openGUI();
            if nargin == 1
                this.w_edt = w_edt;
            end
        end
    end

    methods (Static, Access = public)
        function this = getInstance(w_edt)
            % Get the persistent instance of the class
            persistent unique_instance_gui_rr__

            if nargin == 0
                w_edt = [];
            end

            if isempty(unique_instance_gui_rr__)
                this = GUI_Remote_Resources(w_edt);
                unique_instance_gui_rr__ = this;
            else
                this = unique_instance_gui_rr__;
                this.init();
                this.openGUI();
            end
        end

        function closeGUI()
            fh_list = get(groot, 'Children');
            fig_handle = [];

            % bad code writing style but fast
            for f = 1 : numel(fh_list)
                try
                    if isfield(fh_list(f).UserData, 'name') && strcmp(fh_list(f).UserData.name, GUI_Remote_Resources.WIN_NAME)
                        fig_handle = fh_list(f);
                        break
                    end
                catch
                end
            end
            if ~isempty(fig_handle)
                delete(fig_handle);
            end
        end

    end

    %% METHODS INIT
    % ==================================================================================================================================================
    methods
        function init(this)
            this.log = Core.getLogger();
            this.state = Core.getState();
        end

        function openGUI(this)
            % Main Window ----------------------------------------------------------------------------------------------

            old_win = this.getUniqueWinHandle();
            if ~isempty(old_win)
                log = Core.getLogger();
                log.addMarkedMessage('Resetting the old Inspector window');
                win = old_win;

                if strcmp(this.win.Visible, 'off')
                    if isunix && not(ismac())
                        win.Position(1) = round((win.Parent.ScreenSize(3) - win.Position(3)) / 2);
                        win.Position(2) = round((win.Parent.ScreenSize(4) - win.Position(4)) / 2);
                    else
                        win.OuterPosition(1) = round((win.Parent.ScreenSize(3) - win.OuterPosition(3)) / 2);
                        win.OuterPosition(2) = round((win.Parent.ScreenSize(4) - win.OuterPosition(4)) / 2);
                    end
                end
                this.win = win;
            else
                win = figure( 'Name', 'Remote Resources Inspector', ...
                    'Visible', 'on', ...
                    'DockControls', 'off', ...
                    'MenuBar', 'none', ...
                    'ToolBar', 'none', ...
                    'NumberTitle', 'off', ...
                    'Position', [0 0 1040, 640], ...
                    'Resize', 'on');
                win.UserData.name = this.WIN_NAME;

                this.win = win;

                if isunix && not(ismac())
                    win.Position(1) = round((win.Parent.ScreenSize(3) - win.Position(3)) / 2);
                    win.Position(2) = round((win.Parent.ScreenSize(4) - win.Position(4)) / 2);
                else
                    win.OuterPosition(1) = round((win.Parent.ScreenSize(3) - win.OuterPosition(3)) / 2);
                    win.OuterPosition(2) = round((win.Parent.ScreenSize(4) - win.OuterPosition(4)) / 2);
                end

                try
                    main_vb = uix.VBox('Parent', win, ...
                        'Padding', 5, ...
                        'BackgroundColor', Core_UI.DARKER_GREY_BG);
                catch
                    log = Core.getLogger;
                    log.setOutMode(1,[],0); % to plot a Warning I need to disable GUI and enable
                    log.addError('Please install GUI Layout Toolbox (https://it.mathworks.com/matlabcentral/fileexchange/47982-gui-layout-toolbox)');
                    open('GUI Layout Toolbox 2.3.4.mltbx');
                    log.newLine();
                    log.addWarning('After installation re-run goGPS');
                    close(win);
                    return;
                end
                top_bh = uix.HBox('Parent', main_vb);

                logo_GUI.BG_COLOR = Core_UI.DARK_GREY_BG;
                left_tbv = uix.VBox('Parent', top_bh, ...
                    'BackgroundColor', logo_GUI.BG_COLOR, ...
                    'Padding', 5);

                % Logo/title box -------------------------------------------------------------------------------------------

                logo_g = uix.Grid('Parent', left_tbv, ...
                    'Padding', 5, ...
                    'BackgroundColor', logo_GUI.BG_COLOR);

                logo_ax = axes( 'Parent', logo_g);
                logo_g.Widths = 64;
                logo_g.Heights = 64;
                [logo, transparency] = Core_UI.getLogo();
                logo(repmat(sum(logo,3) == 0,1,1,3)) = 0;
                logo = logo - 20;
                image(logo_ax, ones(size(logo)), 'AlphaData', transparency);
                logo_ax.XTickLabel = [];
                logo_ax.YTickLabel = [];
                axis off;

                Core_UI.insertEmpty(left_tbv, logo_GUI.BG_COLOR);
                left_tbv.Heights = [82 -1];

                % Title Panel -----------------------------------------------------------------------------------------------
                right_tvb = uix.VBox('Parent', top_bh, ...
                    'Padding', 5, ...
                    'BackgroundColor', logo_GUI.BG_COLOR);

                top_bh.Widths = [106 -1];

                title = uix.HBox('Parent', right_tvb, ...
                    'BackgroundColor', logo_GUI.BG_COLOR);
                txt = this.insertBoldText(title, 'Remote Resources in use', 10, Core_UI.LBLUE, 'left');
                txt.BackgroundColor = logo_GUI.BG_COLOR;

                Core_UI.insertEmpty(right_tvb, logo_GUI.BG_COLOR);

                % Disclaimer Panel -----------------------------------------------------------------------------------------------
                txt = this.insertText(right_tvb, {'Inspect the remoute resources locations', '(read-only)'},  8, [], 'left');

                right_tvb.Heights = [20 5 -1];

                string_bh = uix.VBox('Parent', main_vb, ...
                    'Padding', 10, ...
                    'BackgroundColor', GUI_Remote_Resources.BG_COLOR);

                %Core_UI.insertEmpty(string_bh);

                this.j_rrini = this.insertRRBox(string_bh);
                this.j_rrini.setText('RR');

                % Manage dimension -------------------------------------------------------------------------------------------
                main_vb.Heights = [80 -1];
            end
            this.updateResourcePopUpsState();
            this.win.Visible = 'on';
            this.bringOnTop();
        end

        function updateResourcePopUpsState(this)
            % Getting current remote resource manager
            r_man = Remote_Resource_Manager.getInstance();

            state = Core.getCurrentSettings;
            % read current center
            [center_list, center_ss] = r_man.getCenterList();
            cur_center = Core.getState.getCurCenter;
            if isempty(cur_center)
                cur_center = {'default'};
            end
            value = 1;
            while (value < numel(center_list)) && ~strcmp(center_list{value}, cur_center)
                value = value + 1;
            end

            % display resources tree of the current center
            if ~isempty(value)
                try
                    str = r_man.centerToString(state.getRemoteCenter());
                    str = strrep(['% ' str], char(10), [char(10) ' ']);
                catch
                    str = sprintf('[!!] Resource file missing:\n"%s"\nnot found\n\ngoGPS may not work properly', state.getRemoteSourceFile);
                end
                this.j_rrini.setText(str);
            end
        end

    end

    %% METHODS INSERT
    % ==================================================================================================================================================
    methods (Static)
        function txt = insertBoldText(parent, title, font_size, color, alignment)
            if nargin < 4 || isempty(color)
                color = Core_UI.WHITE;
            end
            if nargin < 5 || isempty(alignment)
                alignment = 'center';
            end
            txt = uicontrol('Parent', parent, ...
                'Style', 'Text', ...
                'String', title, ...
                'ForegroundColor', color, ...
                'HorizontalAlignment', alignment, ...
                'FontSize', Core_UI.getFontSize(font_size), ...
                'FontWeight', 'bold', ...
                'BackgroundColor', GUI_Remote_Resources.BG_COLOR);
        end

        function txt = insertText(parent, title, font_size, color, alignment)
            if nargin < 4 || isempty(color)
                color = Core_UI.WHITE;
            end
            if nargin < 5 || isempty(alignment)
                alignment = 'center';
            end
            txt = uicontrol('Parent', parent, ...
                'Style', 'Text', ...
                'String', title, ...
                'ForegroundColor', color, ...
                'HorizontalAlignment', alignment, ...
                'FontSize', Core_UI.getFontSize(font_size), ...
                'BackgroundColor', GUI_Remote_Resources.BG_COLOR);
        end

        function insertLink(parent, url, prefix, prefix_color, link_label, font_size, color, alignment)
            if nargin < 4 || isempty(color)
                color = Core_UI.WHITE;
            end

            % Create and display the text label
            label_str = sprintf('%s %s', ...
                round(prefix_color(1) * 255), round(prefix_color(2) * 255), round(prefix_color(3) * 255), prefix, ...
                round(color(1) * 255), round(color(2) * 255), round(color(3) * 255), link_label);
            j_label = javaObjectEDT('javax.swing.JLabel', label_str);
            bg_color = num2cell(GUI_Remote_Resources.BG_COLOR);
            j_label.setBackground(java.awt.Color(bg_color{:}));
            %% DEPRECATE!!!
            warning on
            [hj_label, h_container] = javacomponent(j_label, [10,10,250,20], parent);
            warning off

            % Modify the mouse cursor when hovering on the label
            hj_label.setCursor(java.awt.Cursor.getPredefinedCursor(java.awt.Cursor.HAND_CURSOR));

            % Set the label's tooltip
            hj_label.setToolTipText(['Visit the ' url ' website']);

            % Set the mouse-click callback
            set(hj_label, 'MouseClickedCallback', @(h,e)web(['http://' url], '-browser'))
        end

        function j_rrini = insertRRBox(container)
            j_rrini = com.mathworks.widgets.SyntaxTextPane;
            %codeType = j_rrini.M_MIME_TYPE;  % j_rrini.contentType='text/m-MATLAB'
            %j_rrini.setContentType(codeType);
            str = ' Remote Resources in use: cd ';
            j_rrini.setText(str);
            % Create the ScrollPanel containing the widget
            j_scroll_settings = com.mathworks.mwswing.MJScrollPane(j_rrini);
            % Inject edit box with the Java Scroll Pane into the main_window
            %% DEPRECATE!!!
            warning off
            [panel_j, panel_h] = javacomponent(j_scroll_settings, [1 1 1 1], container);
            warning on
            j_rrini.setEditable(0);
        end

    end

    %% METHODS EVENTS
    % ==================================================================================================================================================
    methods (Access = public)
    end
end

author: goGPS-Project


text:
%% numOrder
% Below is a demonstration of the features of the |numOrder| function

%%
clear; close all; clc;

%% Syntax
% |[m]=numOrder(A);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function outpoints = tal2mni(inpoints)
% Converts coordinates to MNI brain best guess
% from Talairach coordinates
% FORMAT outpoints = tal2mni(inpoints)
% Where inpoints is N by 3 or 3 by N matrix of coordinates
%  (N being the number of points)
% outpoints is the coordinate matrix with MNI points
% Matthew Brett 2/2/01

% ensure that SPM is available, needed for spm_matrix
hasspm = ft_hastoolbox('spm8up', 3) || ft_hastoolbox('spm2', 1);

dimdim = find(size(inpoints) == 3);
if isempty(dimdim)
  ft_error('input must be a N by 3 or 3 by N matrix')
end
if dimdim == 2
  inpoints = inpoints';
end

% Transformation matrices, different zooms above/below AC
rotn  = spm_matrix([0 0 0 0.05]);
upz   = spm_matrix([0 0 0 0 0 0 0.99 0.97 0.92]);
downz = spm_matrix([0 0 0 0 0 0 0.99 0.97 0.84]);

inpoints = [inpoints; ones(1, size(inpoints, 2))];
% Apply inverse translation
inpoints = inv(rotn)*inpoints;

tmp = inpoints(3,:)<0;  % 1 if below AC
inpoints(:, tmp) = inv(downz) * inpoints(:, tmp);
inpoints(:, ~tmp) = inv(upz) * inpoints(:, ~tmp);
outpoints = inpoints(1:3, :);
if dimdim == 2
  outpoints = outpoints';
end

author: spm


text:
function [model, hasEffect, constrRxnNames, deletedGenes] = deleteModelGenes(model, geneList, downRegFraction)
% Deletes one or more genes and constrain the reactions
% affected to zero and appends '_deleted' to the gene(s)
%
% USAGE:
%
%    [model, hasEffect, constrRxnNames, deletedGenes] = deleteModelGenes(model, geneList, downRegFraction)
%
% INPUT:
%    model:             COBRA model with the appropriate constrains for a
%                       particular condition
%
% OPTIONAL INPUTS:
%    geneList:          List of genes to be deleted (Default =  all genes in model)
%    downRegFraction:   Fraction of the original bounds that the reactions
%                       corresponding to downregulated genes will be assigned
%                       (Default = 0 corresponding to a full deletion)
%
% OUTPUTS:
%    model:             COBRA model with the selected genes deleted
%    hasEffect:         True if the gene deletion has an effect on the model
%    constrRxnNames:    Reactions that are associated to the genes in `deletedGenes`
%    deletedGenes:      The list of genes removed from the model.
%
% .. Authors:
%       - Markus Herrgard 8/28/06
%       - Josh Lerman and Richard Que 04/21/10 - Added an error if non-existent gene.
%       - Richard Que (04/22/2010) - '_deleted' is appended to deleted gene names

if (nargin < 2)
    geneList = model.genes;
end

if (nargin < 3)
    downRegFraction = 0;
end

if (~iscell(geneList))
    geneName = geneList;
    clear geneList;
    geneList{1} = geneName;
end

if (~isfield(model,'genes'))
    error('Gene-reaction associations not included with the model');
end

%RxnGeneMat is required for this function, so we will have to build it if
%it does not exist
if ~isfield(model,'rxnGeneMat')
    model = buildRxnGeneMat(model);
end

hasEffect = false;
constrRxnNames = {};
%deletedGenes is a cell array for returning the genes that are
%eliminated from the model.
deletedGenes = {};

% Find gene indices in model
[isInModel,geneInd] = ismember(geneList,regexprep(model.genes,'_deleted',''));

if (all(isInModel))

  %If there are any zero elements in geneInd remove them from the
  %geneList and the geneInd because they correspond to genes that
  %are not in the model.
  deletedGenes = geneList( find( geneInd ) );
  geneInd = geneInd( find( geneInd ) );

  %mark genes for deletion
  model.genes(geneInd) = strcat(model.genes(geneInd),'_deleted');

    % Find rxns associated with this gene
    rxnInd = find(any(model.rxnGeneMat(:,geneInd),2));
    if (~isempty(rxnInd))
        x = true(size(model.genes));
        % set genes marked "_deleted" to false
        x(~cellfun('isempty',(regexp(model.genes,'_deleted')))) = false;
        constrainRxn = false(length(rxnInd),1);
        % Figure out if any of the reaction states is changed
        for j = 1:length(rxnInd)
            if (isfield(model, 'rules') && ~isempty(model.rules{rxnInd(j)})) %To avoid errors if the rule is empty
                if (~eval(model.rules{rxnInd(j)}))
                    constrainRxn(j) = true;
                end
            end
        end
        % Constrain flux through the reactions associated with these genes
        if (any(constrainRxn))
            constrRxnNames = model.rxns(rxnInd(constrainRxn));
            if (nargin > 2)
                model = changeRxnBounds(model,constrRxnNames,downRegFraction*model.lb(findRxnIDs(model,constrRxnNames)),'l');
                model = changeRxnBounds(model,constrRxnNames,downRegFraction*model.ub(findRxnIDs(model,constrRxnNames)),'u');
            else
                % Full deletion
                model = changeRxnBounds(model,constrRxnNames,0,'b');
            end
            hasEffect = true;
        end
    end
else
    error(['Gene',' ',geneList{~isInModel}, ' not in model!']);
end

author: opencobra

repo: cobratoolbox


text:
function s = prtUtilSimpleStruct(varargin)
%s = prtUtilSimpleStruct(varargin)
%s = prtUtilSimpleStruct('asdf',randn(100,1),'fff',randn(100,1))
%s = prtUtilSimpleStruct('asdf',randn(100,1),'fff',randn(100,1),'asfsadf',prtUtilCellPrintf('%s',num2cell(1:100)'))
%s = prtUtilSimpleStruct(origStruct,...)

% fields = varargin(1:2:end);
% vals = varargin(2:2:end);

if mod(nargin,2) %odd # inputs
    origStruct = varargin{1};

    if isempty(origStruct)
        origStruct = struct;
    end

    if ~isa(origStruct,'struct')
        error('prtUtilSimpleStruct:firstInputStruct','prtUtilSimpleStruct with odd # of inputs expects first input to be a structure');
    end
    varargin = varargin(2:end);
else
    origStruct = struct;
end
structInput = cell(1,length(varargin));

for i = 1:2:length(varargin)
    cField = varargin{i};
    cVal = varargin{i+1};
    assert(isvarname(cField),'can not create struct with a field named %s',cField);

    if iscellstr(cVal)
        %ok
    else
        cVal = mat2cell(cVal,ones(size(cVal,1),1),size(cVal,2));
    end

    structInput(i:i+1) = {cField,cVal};
end
newStruct = struct(structInput{:});

f1 = fieldnames(origStruct);
f2 = fieldnames(newStruct);

%If the new structure data specifies a pre-existing field in origStruct,
%remove the field from origStruct before trying to merge the structures
uNames = unique(cat(1,f1(:),f2(:)));
for i = 1:length(uNames)
    if any(strcmpi(f2,uNames{i})) && any(strcmpi(f1,uNames{i}))
        origStruct = rmfield(origStruct,uNames{i});
    end
end

s = mergestruct(origStruct,newStruct);

function sout = mergestruct(varargin)
%MERGESTRUCT Merge structures with unique fields.
% This is from http://blogs.mathworks.com/loren/2009/10/15/concatenating-structs/
%   Copyright 2009 The MathWorks, Inc.

% Start with collecting fieldnames, checking implicitly
% that inputs are structures.
fn = [];
for k = 1:nargin
    try
        fn = [fn ; fieldnames(varargin{k})];
    catch MEstruct
        throw(MEstruct)
    end
end

% Make sure the field names are unique.
if length(fn) ~= length(unique(fn))
    error('mergestruct:FieldsNotUnique',...
        'Field names must be unique');
end

% Now concatenate the data from each struct.  Can't use
% structfun since input structs may not be scalar.
c = [];
for k = 1:nargin
    try
        c = [c ; struct2cell(varargin{k})];
    catch MEdata
        throw(MEdata);
    end
end

% Construct the output.
sout = cell2struct(c, fn, 1);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%

function setAlwaysOnTopTimer(jd, time, focus, ignoreIfPc)
if nargin<4
    ignoreIfPc=true;
    if nargin<3
        focus=[];
        if nargin<2
            time=1;
        end
    end
end
if isempty(jd)
    return;
end
jd.toFront;
if ispc && ignoreIfPc
    return;
end
%disp(['Setting window on top for ' num2str(time) ' seconds']);
javaMethodEDT( 'setAlwaysOnTop', jd, true);
tmr=timer;
tmr.StartDelay=time;
tmr.TimerFcn=@(h,e)dismiss;
start(tmr);

    function dismiss
        javaMethodEDT('setAlwaysOnTop', jd, false);
        if ~isempty(focus)
            javaMethodEDT('requestFocus',focus);
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function X = VBA_spm_inv(A,TOL)
% inverse for ill-conditioned matrices
% FORMAT X = spm_inv(A,TOL)
%
% A   - matrix
% X   - inverse
%
% TOL - tolerance: default = max(eps(norm(A,'inf'))*max(m,n),exp(-32))
%
% This routine simply adds a small diagonal matrix to A and calls inv.m
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_inv.m 4360 2011-06-14 16:46:37Z ged $

% check A
%--------------------------------------------------------------------------
[m,n] = size(A);
if isempty(A), X = sparse(n,m); return, end

% tolerance
%--------------------------------------------------------------------------
if nargin == 1
    TOL  = max(eps(norm(A,'inf'))*max(m,n),exp(-32));
end

% inverse
%--------------------------------------------------------------------------
X     = inv(A + speye(m,n)*TOL);

author: MBB-team

repo: VBA-toolbox


text:
classdef PTKContext
    % PTKContext. An enumeration used to specify a region of interest
    %
    %     Contexts are used in relation to PTKImageTemplate to request a template
    %     image for a particular region of interest.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    enumeration
        OriginalImage,  % The full image
        LungROI,        % A region containing the lung and airways
        LeftLungROI,    % A region containing the left lung
        RightLungROI,   % A region containing the right lung

        Lungs,          % Left and right lung masks after removal of airways

        LeftLung,       % Left lung mask
        RightLung,      % Right lung mask

        RightUpperLobe,
        RightMiddleLobe,
        RightLowerLobe,
        LeftUpperLobe,
        LeftLowerLobe,

        % Right lung segments
        R_AP,
        R_P,
        R_AN,
        R_L,
        R_M,
        R_S,
        R_MB,
        R_AB,
        R_LB,
        R_PB,

        % Left lung segments
        L_APP,
        L_APP2,
        L_AN,
        L_SL,
        L_IL,
        L_S,
        L_AMB,
        L_LB,
        L_PB

    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function data = simulate_hinge(setup)

%This function simulates the hinge phase of the dynamics

%Set up the integration algorithm
Tspan = setup.Tspan;
IC = [setup.IC.th;
    setup.IC.dth];
eventFunc = @(t,z)events_hinge(t,z,setup.P);
options = odeset(...
    'RelTol',setup.tol,...
    'AbsTol',setup.tol,...
    'Vectorized','on',...
    'MaxStep',setup.odeMaxStep,...
    'Events',eventFunc);
userfun = @(t,z)dynamics_hinge(t,z,setup.P);

%Check that the event function is initially satisfied:
eventVal = feval(eventFunc,Tspan(1),IC);
if (sum(eventVal<0)>0) %IC violates constraints
    tspan = [Tspan(1), Tspan(1) + eps];
    nTime = 2;
    t = tspan;
    Z = IC*ones(1,2);
    [~, C, E] = dynamics_hinge(t,Z,setup.P);

    triggered = eventVal<0;
    index = 1:length(triggered);
    code = index(triggered);
    %Check that the number of possible transitions is one
    if length(code) >1
        disp('WARNING - multiple possible transitions!');
    end

else  %IC satisfied constraints

    %Run simulation
    sol = feval(setup.solver,userfun,Tspan,IC,options);

    %Format for post processing
    tspan = [sol.x(1),sol.x(end)];
    nTime = ceil(setup.dataFreq*diff(tspan));
    t = linspace(tspan(1),tspan(2),nTime);
    Z = deval(sol,t);
    [~, C, E] = dynamics_hinge(t,Z,setup.P);

    code = sol.ie;

end

%Store in a nice format for plotting
data.time = t;
data.state.th = Z(1,:);
data.state.dth = Z(2,:);
data.state.x = zeros(1,nTime) + setup.IC.x;
data.state.dx = zeros(1,nTime);
data.state.y = zeros(1,nTime) + setup.IC.y;
data.state.dy = zeros(1,nTime);
data.contact.h = C(1,:);
data.contact.v = C(2,:);
data.energy.potential = E(1,:);
data.energy.kinetic = E(2,:);
data.P = setup.P;

%Get transitions for finite state machine
data.phase = 'HINGE';
if isempty(code)
    data.exit = 'TIMEOUT';
else
    switch code(end)
        case 1
            data.exit = 'FALL_NEG';
        case 2
            data.exit = 'FALL_POS';
        case 3
            data.exit = 'LIFT';
        case 4
            data.exit = 'SLIP_NEG';
        case 5
            data.exit = 'SLIP_POS';
        otherwise
            error('Invalid exit condition in simulate hinge')
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [geom_data, valid_ids] = Read_geom_files_dynamic(users, params_data_dir)

    geom_data = [];
    valid_ids = [];

    load('../../pca_generation/pdm_68_aligned_wild.mat');

    user1 = {'train_001', 'train_002', 'train_003', 'train_004', 'train_005',...
            'train_006', 'train_007', 'train_008', 'train_009', 'train_010',...
            'train_011', 'train_012', 'train_013', 'train_014', 'train_015',...
            'train_016', 'train_017', 'train_018'};

    user2 = {'train_019', 'train_020', 'train_021', 'train_022', 'train_023',...
             'train_024', 'train_025', 'train_026', 'train_027', 'train_028',...
             'train_029', 'train_030', 'train_031', 'train_032'};

    user3 = {'train_033', 'train_034', 'train_035', 'train_036', 'train_037',...
             'train_038', 'train_039', 'train_040', 'train_041'};

    user4 = {'train_042', 'train_043', 'train_044', 'train_045', 'train_046',...
             'train_047', 'train_048', 'train_049', 'train_050', 'train_051',...
             'train_052', 'train_053', 'train_054', 'train_055', 'train_056'};

    user5 = {'train_057', 'train_058', 'train_059', 'train_060', 'train_061',...
             'train_062', 'train_063'};

    user6 = {'train_064', 'train_065', 'train_066', 'train_067', 'train_068',...
             'train_069', 'train_070', 'train_071', 'train_072', 'train_073',...
             'train_074', 'train_075', 'train_076', 'train_077', 'train_078', 'train_079'};

    user7 = {'train_080', 'train_081', 'train_082', 'train_083', 'train_084',...
             'train_085', 'train_086', 'train_087'};

    users_group = cat(1, {user1}, {user2}, {user3}, {user4}, {user5}, {user6}, {user7});
    user_inds = [];

    for i=1:numel(users)

        geom_file = [params_data_dir, '/au_training_' users{i} '.csv'];

        if(i == 1)
            tab = readtable(geom_file);
            column_names = tab.Properties.VariableNames;
            valid_ind = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'success'));
            shape_inds = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'p_'));
        end

        res = dlmread(geom_file, ',', 1, 0);

        % Check the confidence of detection
        valid = logical(res(:, valid_ind));
        res = res(:, shape_inds);

        % Do not consider global parameters
        res = res(:, 7:end);

        for k=1:numel(users_group)
            if(~isempty(strmatch(users{i}, users_group{k})))
                user_inds = cat(1, user_inds, k * ones(size(valid,1),1));
            end
        end

        actual_locs = res * V';
        res = cat(2, actual_locs, res);

        valid_ids = cat(1, valid_ids, valid);

        geom_data = cat(1, geom_data, res);

    end

    if(numel(users) > 0)

        % Perform normalization here
        u_id = unique(user_inds)';
        for u=u_id
            geom_data(user_inds==u, :) = bsxfun(@plus, geom_data(user_inds==u, :), -median(geom_data(user_inds==u & valid_ids, :)));
        end
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [mu,h] = apgarchcore(parameters, data, ar, ma, x, p, q, y, m, z, v, T)
%{
-----------------------------------------------------------------------
 PURPOSE:
 Estimation of conditional mean and variance of the
 APGARCH: Asymmetric Power GARCH: Ding, Granger and Engle (1993)
-----------------------------------------------------------------------
 USAGE:
 [mu,h] = apgarchcore(parameters, data, p, q, m, T)

 INPUTS:
 parameters:	vector of parameters
 data:         (T x 1) vector of data
 ar:        positive scalar integer representing the order of AR
 am:        positive scalar integer representing the order of MA
 x:         (T x N) vector of factors for the mean process
 p:         positive scalar integer representing the order of ARCH
 q:         positive scalar integer representing the order of GARCH
 y:         (T x N) vector of factors for the volatility process, must be positive!

 OUTPUTS:
 mu:           conditional mean
 h:            conditional variance
-----------------------------------------------------------------------
 Author:
 Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
 Date:     10/2010
 Update 1: 08/2011: ARMA-X Support
-----------------------------------------------------------------------
%}
% Verifying that the vector of parameters is a column vector
[r,c] = size(parameters);
if c>r
    parameters = parameters';
    [r,c] = size(parameters);
end

% Asymmetric  and Nonlinear coefficients
delta = parameters(3+z+p+q+v);
asym=parameters(4+z+p+q:3+z+2*p+q+v);

% Initial parameters
mu = [];
h = [];
mu(1:m,1) = parameters(1);
h(1:m,1) = var(data);

% Dimension of factors
if isscalar(x)
    xy=zeros(size(data,1),1);
else
    xy=x;
end
if isscalar(y)
    yy=zeros(size(data,1),1);
else
    yy=y;
end
% Estimation of Conditional Mean and Variance
for t = (m+1):T;
   mu(t,1) = parameters(1:1+z)'*[1; data(t-(1:ar)); data(t-(1:ma))-mu(t-(1:ma),1); xy(t,:)'*ones((isscalar(x) < 1))];
   h(t,1) = (parameters(2+z:2+z+p+q+v)'*[1; (abs(data(t-(1:p)) - mu(t-(1:p))) - asym(1:p)*(data(t-(1:p)) - mu(t-1))).^delta; (h(t-(1:q))).^delta; yy(t,:)'*ones((isscalar(y) < 1))]).^(1/delta);
end

end


text:
function LeftStanceConstraints(nlp, bounds, varargin)
    domain = nlp.Plant;

    ip = inputParser;
    ip.addParameter('LoadPath',[],@ischar);
    ip.parse(varargin{:});

    %% virtual constraints
    opt.constraint.virtual_constraints(nlp, bounds, ip.Results.LoadPath);
    %% foot clearance
    [right_foot_frame] = sys.frames.RightFoot(domain);
    opt.constraint.foot_clearance(nlp, bounds, right_foot_frame);

    %% swing toe position
    opt.constraint.step_distance(nlp, bounds);

    %% swing foot orientation
    %     opt.constraint.foot_orientation(nlp, bounds, 'Right');
    %% swing foot velocity
    opt.constraint.impact_velocity(nlp, bounds, right_foot_frame);

    %% feet distance
    %     opt.constraint.feet_distance(nlp, bounds);

    opt.constraint.yaw_start(nlp, bounds);

    opt.constraint.knee_angle(nlp, bounds);

    opt.constraint.average_velocity(nlp, bounds);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function f = sqrt(f, pref)
%SQRT   Square root of a CHEBFUN.
%   SQRT(F) returns the square root of a CHEBFUN F.
%
% See also POWER.

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

% Trivial case: (f is empty)
if ( isempty(f) )
    return
end

if ( nargin < 2 )
    pref = chebfunpref();
end

% Simply call POWER()
f = power(f, 0.5, pref);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [Coor, EffIdx, ind] = CvtPtsToVoxelFun( data, minVal, Size, gridRes )
Coor = [];
EffIdx = [];
ind = [];
Dim = size(data, 1);
for i = 1 : 1 : Dim
    Coor(i, :) = ceil(( data(i, :) - minVal(i) )/gridRes);
end
if Dim == 3
    EffIdx = find( Coor(1, :) >= 1 & Coor(1, :) <= Size(1) & ...
        Coor(2, :) >= 1 & Coor(2, :) <= Size(2) & ...
        Coor(3, :) >= 1 & Coor(3, :) <= Size(3) );
    ind = sub2ind(Size, Coor(1, EffIdx), Coor(2, EffIdx), Coor(3, EffIdx));
end
if Dim == 2
    EffIdx = find( Coor(1, :) >= 1 & Coor(1, :) <= Size(1) & ...
        Coor(2, :) >= 1 & Coor(2, :) <= Size(2) );
    ind = sub2ind(Size, Coor(1, EffIdx), Coor(2, EffIdx));
end
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function [ratio,initdl] = fitint(ratio);
%fitint   computes contraction/expansion ratio
%   [ratio,initdl] = fitint(ratio);
%   input
%          ratio      initial guess
%   output
%          ratio      converged value
%          initdl     associated first increment
%
%   called by subint
%   IFISS function: DJS; 1 April 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      global global_N global_INTL global_LASTDL
      if(global_N==1)
      initdl=global_INTL;
      ratio=global_LASTDL/initdl;
      return
      else
      x=ratio;
% call nonlinear equation solver
%     fprintf('\n      x            f   \n')
      ratio=fzero('fint',x,optimset('Display','off'));
      initdl=global_LASTDL/(ratio^global_N);
      end
% check solution validity
      if ratio <= 1,
         fprintf('\n\n\nInfeasible ratio computed for stretch grid.\n');
         fprintf('Try changing the initial value for "ratio" in function "subint.m"\n');
         error(' ');
      elseif ratio >= 2
         fprintf('\n\n*** warning \n computed ratio is large. \n')
      end
      return

author: johannesgerer

repo: jburkardt-m


text:
%kkmsquish 'Compress Map to One Column by Means of Average, RMS, or MAX'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kmsquish.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% MultiChoice: t 'Squishing scheme', default: 1: 'squishing scheme (1)average, (2)RMS, (3)Maximum'
%    Choices are:
%   1: 'Average of map columns'
%   2: 'RMS value of map columns'
%   3: 'Maximum value of mape columns'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkmsquish(i, {'i','';'t',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kmsquish - Compress Map to One Column by Means of Average, RMS, or MAX
%
%  DESCRIPTION
% kmsquish compresses a map down to one column, effectively
% converting the map to a "greyscale" map. An example of
% this operation is encountered when trying to print a
% pseudocolor RGB image on a B/W laser printer: the RGB
% color information must be somehow converted to greyscale
% in order for the printer to know how to dither the image
% during printing.
%
% kmsquish is very happy to operate on simple 3-column
% maps, such as would be encountered when processing a
% color image in any of the standard coordinate systems
% (ntscRGB, etc). However, kmsquish is not limited to 3 column
% maps - it will operate on maps with any number of cloumns.
%
% The conversion from multiple columns to a single column
% is done by computing a grey level for each row of the map,
% one at a time.  For t=1, the grey level for a row is the
% average of the columns in the row. For t=2, the grey
% level is the RMS value of the columns in the row. For t=3,
% the grey level is the maximum value found in the colums of
% the row. Clearly, these schemes are not optimal in any sense,
% but one of them will usually give an acceptable result.
%
% Ordinarily (map=0), the map is the only part of the input
% object that is modified. The value data is not touched; it
% still points to various rows in the map just as before the
% compression. If map=1, the value data is pulled thru the
% map and the map segment itself is removed. If mask data
% is present, then this operation proceeds as described in
% the man page for kmapdata(1).
%
%
%
%  EXAMPLES
%
% Squishing and mapping an image that is to be printed.
% kmsquish -i ball.xv.color -o ball.ready.to.print
% here, -m pand -t do not need to be specified since the defaults are
% set correctly.
%
%  "SEE ALSO"
% igamut(1), kmapdata(1)
%
%  RESTRICTIONS
% Complex data is not supported.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkmsquish(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkmsquish(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'t', 1;'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','MultiChoice','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kmsquish"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [sens] = ft_read_sens(filename, varargin)

% FT_READ_SENS read sensor positions from various manufacturer specific files. See
% further down for the list of file types that are supported.
%
% Use as
%   grad = ft_read_sens(filename, ...)  % for gradiometers
%   elec = ft_read_sens(filename, ...)  % for electrodes
%
% Additional options should be specified in key-value pairs and can be
%   'fileformat'     = string, see the list of supported file formats (the default is determined automatically)
%   'senstype'       = string, can be 'eeg', 'meg' or 'nirs', specifies which type of sensors to read from the file (default = 'eeg')
%   'coordsys'       = string, 'head' or 'dewar' (default = 'head')
%   'coilaccuracy'   = scalar, can be empty or a number (0, 1 or 2) to specify the accuracy (default = [])
%
% The electrode, gradiometer and optode structures are defined in more detail
% in FT_DATATYPE_SENS.
%
% Files from the following acquisition systems and analysis platforms file formats
% are supported.
%   asa_elc besa_elp besa_pos besa_sfp yokogawa_ave yokogawa_con yokogawa_raw 4d
%   4d_pdf 4d_m4d 4d_xyz ctf_ds ctf_res4 itab_raw itab_mhd netmeg neuromag_fif
%   neuromag_mne neuromag_mne_elec neuromag_mne_grad polhemus_fil polhemus_pos
%   zebris_sfp spmeeg_mat eeglab_set localite_pos artinis_oxy3 artinis_oxyproj matlab
%
% See also FT_READ_HEADER, FT_DATATYPE_SENS, FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD,

% Copyright (C) 2005-2018 Robert Oostenveld
%
% 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$

% optionally get the data from the URL and make a temporary local copy
filename = fetch_url(filename);

% get the options
fileformat     = ft_getopt(varargin, 'fileformat', ft_filetype(filename));
senstype       = ft_getopt(varargin, 'senstype');         % can be eeg or meg, default is automatic when []
coordsys       = ft_getopt(varargin, 'coordsys', 'head'); % this is used for ctf and neuromag_mne, it can be head or dewar
coilaccuracy   = ft_getopt(varargin, 'coilaccuracy');     % empty, or a number between 0 to 2

realtime = any(strcmp(fileformat, {'fcdc_buffer', 'ctf_shm', 'fcdc_mysql'}));

if realtime
  % skip the rest of the initial checks to increase the speed for realtime operation
else
  % test whether the file exists
  if ~exist(filename, 'file')
    ft_error('file ''%s'' does not exist', filename);
  end
end

switch fileformat
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % read the content from various files that contain EEG electrode positions
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case 'asa_elc'
    sens = read_asa_elc(filename);

  case 'artinis_oxy3'
    ft_hastoolbox('artinis', 1);
    hdr = read_artinis_oxy3(filename, false);
    sens = hdr.opto;

  case 'artinis_oxyproj'
    ft_hastoolbox('artinis', 1);
    hdr = read_artinis_oxyproj(filename);
    sens = hdr.opto;

  case 'polhemus_pos'
    sens = read_polhemus_pos(filename);

  case 'besa_elp'
    ft_error('unknown fileformat for electrodes or gradiometers');
    % the code below does not yet work
    fid = fopen_or_error(filename);
    % the ascii file contains: type, label, angle, angle
    tmp = textscan(fid, '%s%s%f%f');
    fclose(fid);
    sel = strcmpi(tmp{1}, 'EEG');  % type can be EEG or POS
    sens.label = tmp{2}(sel);
    az = tmp{3}(sel) * pi/180;
    el = tmp{4}(sel) * pi/180;
    r  = ones(size(el));
    [x, y, z] = sph2cart(az, el, r);
    sens.chanpos = [x y z];

  case 'besa_pos'
    tmp = importdata(filename);
    if ~isnumeric(tmp)
      ft_error('unexpected file format for fileformat=besa_pos')
    end
    [nchan,nrow] = size(tmp);
    if nrow==3
      sens.pnt = tmp;
    elseif nrow==9
      pnt1 = tmp(:,1:3);  % bottom coil
      pnt2 = tmp(:,4:6);  % top coil
      ori  = tmp(:,7:9);  % orientation of bottom coil
      sens.pnt = [pnt1; pnt2];
      sens.ori = [ori; ori];
      sens.tra = [eye(nchan) -eye(nchan)];
    else
      ft_error('unexpected file format for fileformat=besa_pos')
    end
    [p, f, x] = fileparts(filename);
    elpfile = fullfile(p, [f '.elp']);
    elafile = fullfile(p, [f '.ela']);
    if exist(elpfile, 'file')
      ft_warning('reading channel labels from %s', elpfile);
      % read the channel names from the accompanying ELP file
      lbl = importdata(elpfile);
      sens.label = strrep(lbl.textdata(:,2) ,'''', '');
    elseif exist(elafile, 'file')
      ft_warning('reading channel labels from %s', elafile);
      % read the channel names from the accompanying ELA file
      lbl = importdata(elafile);
      lbl = strrep(lbl, 'MEG ', ''); % remove the channel type
      lbl = strrep(lbl, 'EEG ', ''); % remove the channel type
      sens.label = lbl;
    else
      % the file does not have channel labels in it
      ft_warning('creating fake channel names for besa_pos');
      for i=1:nchan
        sens.label{i} = sprintf('%03d', i);
      end
    end

  case 'besa_sfh'
    sfh = readBESAsfh(filename);
    sens.label   = sfh.SurfacePointsLabels(:);
    sens.elecpos = sfh.SurfacePointsCoordinates(:,1:3);
    sel = true(sfh.NrOfPoints, 1);
    for i=1:sfh.NrOfPoints
      tok = tokenize(sens.label{i}, '_');
      sens.label{i} = tok{2};
      sel(i) = ~strcmp(tok{1}, 'Fid');
    end
    sens.label   = sens.label(sel);
    sens.elecpos = sens.elecpos(sel,:);

  case 'besa_sfp'
    [lab, pos] = read_besa_sfp(filename);
    sens.label   = lab;
    sens.elecpos = pos;

  case 'bioimage_mgrid'
    sens = read_bioimage_mgrid(filename);

  case {'ctf_ds', 'ctf_res4', 'ctf_old', 'neuromag_fif', 'neuromag_mne', '4d', '4d_pdf', '4d_m4d', '4d_xyz', 'yokogawa_ave', 'yokogawa_con', 'yokogawa_raw', 'ricoh_ave', 'ricoh_con', 'itab_raw' 'itab_mhd', 'netmeg'}
    % gradiometer information is always stored in the header of the MEG dataset, hence uses the standard fieldtrip/fileio ft_read_header function
    hdr = ft_read_header(filename, 'headerformat', fileformat, 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    % sometimes there can also be electrode position information in the header
    if isfield(hdr, 'elec') && isfield(hdr, 'grad')
      if isempty(senstype)
        % set the default
        ft_warning('both electrode and gradiometer information is present, returning the electrode information by default');
        senstype = 'eeg';
      end
      switch lower(senstype)
        case 'eeg'
          sens = hdr.elec;
        case 'meg'
          sens = hdr.grad;
        otherwise
          ft_error('incorrect specification of senstype');
      end
    elseif isfield(hdr, 'grad')
      sens = hdr.grad;
    elseif isfield(hdr, 'elec')
      sens = hdr.elec;
    else
      ft_error('neither electrode nor gradiometer information is present');
    end

  case {'curry_dat', 'curry_cdt'}

    hdr = ft_read_header(filename);

    if ~isempty(hdr.orig.sensorpos)
      sens.elecpos = hdr.orig.sensorpos';
      sens.label   = hdr.label(1:size(sens.elecpos, 1));
    end

  case 'fcdc_buffer'
    % the online header should have a binary blob with the sensor information
    hdr = ft_read_header(filename, 'headerformat', fileformat);
    if isfield(hdr, 'orig') && isfield(hdr.orig, 'neuromag_header')
      hdr = decode_fif(hdr.orig.neuromag_header);
    end
    if isfield(hdr, 'orig') && isfield(hdr.orig, 'ctf_res4')
      hdr = decode_res4(hdr.orig.ctf_res4);
    end
    if isempty(senstype)
      % set the default
      ft_warning('both electrode and gradiometer information is present, returning the electrode information by default');
      senstype = 'eeg';
    end
    switch lower(senstype)
      case 'eeg'
        sens = hdr.elec;
      case 'meg'
        sens = hdr.grad;
      otherwise
        ft_error('incorrect specification of senstype');
    end

  case 'neuromag_mne_grad'
    % the file can contain both, force reading the gradiometer info
    % note that this functionality overlaps with senstype=eeg/meg
    hdr = ft_read_header(filename, 'headerformat', 'neuromag_mne', 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    sens = hdr.grad;

  case 'neuromag_mne_elec'
    % the file can contain both, force reading the electrode info
    % note that this functionality overlaps with senstype=eeg/meg
    hdr = ft_read_header(filename, 'headerformat', 'neuromag_mne', 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    sens = hdr.elec;

  case {'spmeeg_mat', 'eeglab_set'}
    % this is for EEG formats where electrode positions can be stored with the data
    hdr = ft_read_header(filename, 'coordsys', coordsys, 'coilaccuracy', coilaccuracy);
    if isfield(hdr, 'grad')
      sens = hdr.grad;
    elseif isfield(hdr, 'elec')
      sens = hdr.elec;
    else
      ft_error('no electrodes or gradiometers found in the file')
    end

  case 'polhemus_fil'
    % these are created at the FIL in London with a polhemus tracker
    [sens.fid.pnt, sens.pnt, sens.fid.label] = read_polhemus_fil(filename, 0);
    % the file does not have channel labels in it
    ft_warning('no channel names in polhemus file, using numbers instead');
    for i=1:size(sens.pnt, 1)
      sens.label{i} = sprintf('%03d', i);
    end

  case 'matlab'
    % MATLAB files can contain all sensor arrays
    matfile = filename;   % this solves a problem with the MATLAB compiler v3
    var = whos('-file', filename);
    sel = intersect({'elec', 'grad', 'opto', 'sens', 'elc'}, {var(:).name});
    if numel(sel)==1
      % read the specific variable
      sens = loadvar(matfile, sel{1});
    else
      % read whatever variable is in the file, this will error if the file contains multiple variables
      sens = loadvar(matfile);
    end

  case 'zebris_sfp'
    % these are created by a Zebris tracker, at CRC in Liege at least.
    [sens.fid.pnt, sens.chanpos, sens.fid.label, sens.label] = read_zebris(filename, 0);
    % convert to columns
    sens.label = sens.label(:);
    sens.fid.label = sens.fid.label(:);

  case '4d_el_ascii'
    fid = fopen_or_error(filename, 'rt');
    c = textscan(fid, '%s%s%f%f%f');
    l = c{:,1}; % label
    s = c{:,2}; % status, it can be 'Collected' or empty
    x = c{:,3};
    y = c{:,4};
    z = c{:,5};
    % shift the columns with one where the status is not specified
    sel = isnan(z);
    z(sel) = y(sel);
    y(sel) = x(sel);
    x(sel) = str2double(s(sel));
    s(sel) = {''};
    fclose(fid);
    if false
      % return all positions, including the ones that do not correspond to
      % electrodes per see, such as the fiducials and localizer coils
      sens          = [];
      sens.label    = l;
      sens.elecpos  = [x y z];
    else
      % split the electrodes and fiducials
      % this is consistent with zebris_sfp and with the output of ft_read_headshape
      sens            = [];
      sens.label      = l(~sel);
      sens.elecpos    = [x(~sel) y(~sel) z(~sel)];
      sens.fid.label  = l(sel);
      sens.fid.pnt    = [x(sel) y(sel) z(sel)];
    end

  case {'localite_pos','localite_ins'}
    if ~usejava('jvm') % Using xml2struct requires java
      fid = fopen_or_error(filename);

      % Read marker-file and store contents in cells of strings
      tmp = textscan(fid,'%s');

      fclose(fid);

      % Search for cells that contain coordinates
      selx = strncmp('data0',tmp{1},5);
      sely = strncmp('data1',tmp{1},5);
      selz = strncmp('data2',tmp{1},5);
      sellab = strncmp('description',tmp{1},5);

      % Extract cells that contain coordinates
      xtemp  = tmp{1}(selx);
      ytemp  = tmp{1}(sely);
      ztemp  = tmp{1}(selz);
      labtemp = tmp{1}(sellab);

      % Determine which channels are set. In localite channels that are not set
      % automatically receive coordinates [0, 0, 0] and should therefore
      % be discarded.
      settemp = tmp{1}(strncmp('set',tmp{1},3));
      selset = strncmp('set="f',settemp,6);

      % Remove channels that are not set
      xtemp(selset) = [];
      ytemp(selset) = [];
      ztemp(selset) = [];
      labtemp(selset) = [];

      % Convert cells that contain coordinates from string to double
      x = [];
      y = [];
      z = [];
      lbl = [];

      for i=1:numel(xtemp)
        x(i,1) = str2double(xtemp{i}(8:end-1));
        y(i,1) = str2double(ytemp{i}(8:end-1));
        z(i,1) = str2double(ztemp{i}(8:end-3));
        lbl{i,1} = labtemp{i}(14:end-1);
      end

      % Create and fill sens structure
      sens = [];
      sens.elecpos = [x y z];
      sens.chanpos = sens.elecpos;
      sens.label = lbl;
    else
      tmp = xml2struct(filename);

      sens = [];

      % Loop through structure obtained from xml-file and store
      % coordinate information into sens structure of channels that are
      % set.
      for i=1:numel(tmp)
        if strcmp(tmp(i).Marker.set, 'true')
          sens.elecpos(i,1) = str2double(tmp(i).Marker.ColVec3D.data0);
          sens.elecpos(i,2) = str2double(tmp(i).Marker.ColVec3D.data1);
          sens.elecpos(i,3) = str2double(tmp(i).Marker.ColVec3D.data2);
          sens.label{i} = tmp(i).Marker.description;
        end
      end

      sens.chanpos = sens.elecpos;
    end

  case 'easycap_txt'
    % Read the file and store all contents in cells of strings
    fid = fopen_or_error(filename);
    tmp = textscan(fid,'%s%s%s%s');
    fclose(fid);

    sens = [];
    if all(cellfun(@isempty, tmp{4}))
      % it contains theta and phi
      sens.label = tmp{1}(2:end);
      theta = cellfun(@str2double, tmp{2}(2:end));
      phi   = cellfun(@str2double, tmp{3}(2:end));
      radians = @(x) pi*x/180;
      ft_warning('assuming a head radius of 85 mm');
      x = 85*cos(radians(phi)).*sin(radians(theta));
      y = 85*sin(radians(theta)).*sin(radians(phi));
      z = 85*cos(radians(theta));
      sens.unit = 'mm';
      sens.elecpos = [x y z];
      sens.chanpos = [x y z];
    else
      % it contains X, Y, Z
      sens.label = tmp{1}(2:end);
      x = cellfun(@str2double, tmp{2}(2:end));
      y = cellfun(@str2double, tmp{3}(2:end));
      z = cellfun(@str2double, tmp{4}(2:end));
      sens.elecpos = [x y z];
      sens.chanpos = [x y z];
    end

  case 'neuromag_iso'
    ft_hastoolbox('mne', 1);
    FIFF = fiff_define_constants();
    [fid, tree, dir] = fiff_open(filename);
    isotrak = fiff_dir_tree_find(tree, FIFF.FIFFB_ISOTRAK);
    sel = find([isotrak.dir.kind]==FIFF.FIFF_DIG_POINT);
    sens = [];
    sens.elecpos = nan(numel(sel),3);
    sens.chanpos = nan(numel(sel),3);
    coordsys     = nan(numel(sel),1);
    for i=sel
      tag = fiff_read_tag(fid,isotrak.dir(i).pos);
      sens.elecpos(i,:) = tag.data.r;
      sens.chanpos(i,:) = tag.data.r;
      coordsys(i)       = tag.data.coord_frame;
    end
    fclose(fid);

    if all(coordsys==FIFF.FIFFV_COORD_DEVICE)
      sens.coordsys = 'device';
    elseif all(coordsys==FIFF.FIFFV_COORD_ISOTRAK)
      sens.coordsys = 'isotrak';
    elseif all(coordsys==FIFF.FIFFV_COORD_HPI)
      sens.coordsys = 'hpi';
    elseif all(coordsys==FIFF.FIFFV_COORD_HEAD)
      sens.coordsys = 'head';
    else
      sens.coordsys = 'unknown';
    end

    ft_warning('creating fake channel names for neuromag_iso');
    for i=1:size(sens.chanpos,1)
      sens.label{i} = sprintf('%d', i);
    end

  case 'neuromag_cal'
    dat = cell(1,14);
    [dat{:}] = textread(filename, '%s%f%f%f%f%f%f%f%f%f%f%f%f%f');
    label = dat{1};
    x = dat{2};
    y = dat{3};
    z = dat{4};
    rot = cat(2, dat{5:13});
    cal = dat{14}; % ??
    % construct the sensor structucture
    % it seems that the channel positions are expressed in dewar cordinates
    % it would be possible to use coil_def.dat to construct the coil positions
    sens.label = label;
    sens.chanpos = [x y z];

  otherwise
    ft_error('unknown fileformat for electrodes or gradiometers');
end % switch fileformat

% ensure that the sensor description is up-to-date
% this will also add chantype and units to the sensor array if missing
sens = ft_datatype_sens(sens);

% ensure that the output is consistent with the type requested by the user
if strcmpi(senstype, 'meg')
  assert(isfield(sens,'coilpos'), 'cannot read gradiometer information from %s', filename);
elseif strcmpi(senstype, 'eeg')
  assert(isfield(sens,'elecpos'), 'cannot read electrode information from %s', filename);
else
  % it is empty if not specified by the user, in that case either one is fine
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cam.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function f = cam(x,y)
% Four-hump camel
function f = cam(x,y)
if nargin == 1
  x1 = x(1);
  x2 = x(2);
else
  x1 = x;
  x2 = y;
end
f=(4-2.1.*x1.^2+x1.^4./3).*x1.^2+x1.*x2+(-4+4.*x2.^2).*x2.^2;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function lambdamax = tvdiplmax(y)
% Calculate the value of lambda so that if lambda >= lambdamax, the TVD
% functional solved by TVDIP is minimized by the trivial constant
% solution x = mean(y). This can then be used to determine a useful range
% of values of lambda, for example.
%
% Usage:
% lambdamax = tvdiplmax(y)
%
% Input arguments:
% - y          Original signal to denoise, size N x 1.
%
% Output arguments:
% - lambdamax  Value of at which x = mean(y) is the output of the TVDIP
%              function.
%
% (c) Max Little, 2010. Based around code originally written by
% S.J. Kim, K. Koh, S. Boyd and D. Gorinevsky. If you use this code for
% your research, please cite:
% M.A. Little, Nick S. Jones (2010)
% "Sparse Bayesian Step-Filtering for High-Throughput Analysis of Molecular
% Machine Dynamics", in 2010 IEEE International Conference on Acoustics,
% Speech and Signal Processing, 2010, ICASSP 2010 Proceedings.
%
% This code is released under the terms of GNU General Public License as
% published by the Free Software Foundation; version 2 or later.

error(nargchk(1,1,nargin));
y = y(:);
N = length(y);
M = N - 1;

% Construct sparse operator matrices
I1 = speye(M,M);
O1 = spalloc(M,1,M);
D = [I1 O1]-[O1 I1];

DDT = D*D';
Dy  = D*y;

lambdamax = max(abs(DDT\Dy));

author: beckel

repo: nilm-eval

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


text:
function [A,varargout] = selectMaxbyColumn(A,varargin)
% find maximum in each row

% find maximum values
ind = A == repmat(max(A,[],1),size(A,1),1);

% select only the first maximum
ind = ind & ind == cumsum(ind,1);

% return results
A = A(ind);

for i = 1:nargout-1
  varargout{i} = varargin{i}(ind);
end

author: mtex-toolbox

repo: mtex

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


text:
function [a, s]=getarrowslice
% Return the Arrow and Slice based on the GCO

  if isempty(getappdata(gco,'controlarrow')) & ...
        isempty(getappdata(gco,'isosurface'))
    a = gco;
    s = getappdata(a,'arrowslice');
    if isempty(s)
      s=getappdata(a,'arrowiso');
    end
  else
    s = gco;
    if ~isempty(getappdata(s,'isosurface'))
      s=getappdata(s,'isosurface');
    end
    a = getappdata(s,'controlarrow');
  end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% Multi-axle steering simple
% This script presents open-loop multi-axle steering of an simple vehicle.
%
% <<MultiaxleSteeringSimple.gif>>
%
%% Simulation models and parameters
% First, all classes of the package are imported with

clear ; close all ; clc

import VehicleDynamicsLateral.*

%%
% Choosing tire and vehicle model. In this case, the parameters are defined
% by the user.

% Choosing tire
TireModel           = TirePacejka();

% Choosing vehicle
VehicleModel        = VehicleSimpleNonlinear();
VehicleModel.tire   = TireModel;

%% Simulation parameters
% Choosing simulation time span
%

T       = 4;                            % Total simulation time     [s]
resol   = 80;                           % Resolution
TSPAN   = 0:T/resol:T;                  % Time span                 [s]

%% Open-loop steering input
% Single period sine wave with:
%

T_period            = 2;                % Steering single period    [s]
freq                = 1/T_period;       % Steering frequency        [Hz]
delta_freq          = freq*2*pi;        % Steering frequency        [rad/s]

% End input index
% sine_index          = find(TSPAN > T_period);

% Defining steering input
VehicleModel.deltaf = 15*pi/180*sin(delta_freq*TSPAN);
VehicleModel.deltar = 15*pi/180*sin(delta_freq*TSPAN);

figure
subplot(2,1,1)
    hold on ; grid on ; box on
    plot(TSPAN,VehicleModel.deltaf*180/pi,'r','linewidth',2)
    ylabel('Delta F [deg]')
subplot(2,1,2)
    hold on ; grid on ; box on
    plot(TSPAN,VehicleModel.deltar*180/pi,'r','linewidth',2)
    ylabel('Delta R [deg]')
    xlabel('Time [s]')

%%
% To define a simulation object (simulator) the arguments must be the
% vehicle object and the time span.

simulator = Simulator(VehicleModel, TSPAN);

%%
% Initial conditions
% Changing initial conditions of the simulation object

simulator.V0 = 60/3.6;              % Initial velocity              [m/s]

%% Run simulation
% To simulate the system we run the Simulate method of the simulation
% object.

simulator.Simulate();

%% Results
%

g = Graphics(simulator);
g.Frame();
g.Animation();
% g.Animation('html/MultiaxleSteeringSimple');       % Uncomment to save animation gif

%% See Also
%
% <../../../index.html Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral