text:
function [synCat] = syn_catalog(nEvents, fBValue, fMc, fInc, fMinLat, fMaxLat, fMinLon, fMaxLon, fMinDepth, fMaxDepth, fMinTime, fMaxTime)
    % SYN_CATALOG Creates a synthetic catalog.
    %
    % [mNewCatalog] = SYNCATALOG(nNumberEvents, fBValue, fMc, fInc, fMinLat, fMaxLat,
    %                                      fMinLon, fMaxLon, fMinDepth, fMaxDepth, fMinTime, fMaxTime)
    % ------------------------------------------------------------------------------------------------
    %
    % Input parameters:
    %   nEvents     Number of events in the catalog
    %   fBValue           b-value of the catalog
    %   fMc               Magnitude of completeness of the catalog (= minimum magnitude)
    %   fInc              Magnitude increment (usually 0.1)
    %   fMinLat           Minimum latitude value of events in the catalog
    %   fMaxLat           Maximun latitude value of events in the catalog
    %   fMinLon           Minimum longitude value of events in the catalog
    %   fMaxLon           Maximun longitude value of events in the catalog
    %   fMinDepth         Minimum depth of events in the catalog (positive value)
    %   fMaxDepth         Maximun depth of events in the catalog (positive value)
    %   fMinTime          Minimum date/time of events in the catalog (decimal year: e.g. 1983.5)
    %   fMaxTime          Maximun date/time of events in the catalog (decimal year: e.g. 1987.9)
    %
    % Output parameters:
    %   synCat       Synthetic catalog
    %
    % Danijel Schorlemmer, Thomas van Stiphout, Celso Reyes
    % 2002, 2007, 2019

    % If other options for catalog creation are requested, then create_synthetic_catalog subfunction
    % (below) can be activated and wired up.  Default behavior is "homogenous"

    report_this_filefun();

    tb = table;

    rng('shuffle');
    getRandomVals = @(minv,maxv) rand(nEvents,1) * (maxv - minv) + minv;

    tb.Date = dateshift(getRandomVals(fMinTime, fMaxTime),'start','minute');
    tb.Magnitude = syn_create_magnitudes(nEvents, fBValue, fMc, fInc);
    tb.Longitude = getRandomVals(fMinLon, fMaxLon);
    tb.Latitude  = getRandomVals(fMinLat, fMaxLat);
    tb.Depth     = getRandomVals(fMinDepth, fMaxDepth);

    synCat = ZmapCatalog.fromTable(tb);
    synCat.Name = 'synthetic';
    synCat.sort('Date','ascend');
    return

    %% code salvaged from thomas/synthetic/syn_catalog

    function synCat = create_synthetic_catalog(synMode)
        % allowing some choices in the catalog creation method
        % nSynMode          Type of synthetic catalog 0:homogeneous distr; 1:based
        % on real catalog 2:hypocenter based on real catalog, magnitude and focal
        % time is randomly computed
        % mCatalog          declustered catalog needed only for nSynMode=1
        %
        % Output parameters:
        %   mNewCatalog       Synthetic catalog
        %
        % Danijel Schorlemmer
        % April 17, 2002
        % Updates
        % Mai  9, 2007  van Stiphout, Thomas     replaced datevec.m by decyear2mat.m

        if isnumeric(synMode)
            switch synMode
                case 0
                    synMode = "homogeneous"; % homogenous distribution
                case 1
                    synMode = "catalog-based"; % based on real catalog"
                case 2
                    synMode = "hypocenter-based"; % magnitude and focal time is randomly computed.
            end
        end

        rng('shuffle');
        getRandomVals = @(minv,maxv) rand(nEvents,1) * (maxv - minv) + minv;

        tb = table;

        switch synMode
            case "homogeneous"
                tb.Date = dateshift(getRandomVals(fMinTime, fMaxTime),'start','minute');
                tb.Magnitude = syn_create_magnitudes(nEvents, fBValue, fMc, fInc);
                tb.Longitude = getRandomVals(fMinLon,fMaxLon);
                tb.Latitude  = getRandomVals(fMinLat,fMaxLat);
                tb.Depth     = getRandomVals(fMinDepth, fMaxDepth);

            case"catalog-based"
                % take from real catalog hypocenters and magnitudes  ....
                dates = dateshift(getRandomVals(fMinTime, fMaxTime),'start','minute');
                [tb.Date, idx] = sort(dates);
                tb.Magnitude = mCatalog.Magnitude(idx);
                tb.Longitude = mCatalog.Longitude(idx);
                tb.Latitude = mCatalog.Latitude(idx);
                tb.Depth = mCatalog.Depth(idx);

            case "hypocenter-based"

                hypoShift = 5; % initial shift of hypocenter from reference eevent.
                dates = dateshift(getRandomVals(fMinTime, fMaxTime),'start','minute');
                [tb.Date, idx] = sort(dates);
                tb.Magnitude = syn_create_magnitudes(nEvents, fBValue, fMc, fInc);

                vTmp = km2deg(rand(nNumberEvents,2).* hypoShift - hypoShift/2);

                tb.Longitude = mCatalog.Longitude(idx)+vTmp;
                tb.Latitude = mCatalog.Latitude(idx)+vTmp;
                tb.Depth = mCatalog.Depth(idx);

            otherwise
                error('unknown synthetic catalog creation mode, should be one of "homogenous", "catalog-based", "hypocenter-based"')
        end

        synCat = ZmapCatalog.fromTable(tb);
        synCat.Name = 'synthetic';
    end
end

author: CelsoReyes


text:
function [a, b, c, d] = corners(M)
    a = M(1, 1);
    b = M(1, end);
    c = M(end, 1);
    d = M(end, end)
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function K = validateScsCone(K)
if (~isfield(K, 'f'))
    K.f = 0;
end
if (~isfield(K, 'l'))
    K.l = 0;
end
if (~isfield(K, 'ep'))
    K.ep = 0;
end
if (~isfield(K, 'ed'))
    K.ed = 0;
end
if (~isfield(K, 'q'))
    K.q = [];
end
if (~isfield(K, 's'))
    K.s = [];
end
if (~isfield(K, 'p'))
    K.p = [];
end
end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [bit_errs_ss,i_out,q_out]=despread(tiq,sumiq,cs_t,sn_t,seq1,seq2,data,N,fcarr)
%Bob Gess, June 2008, for EE473
%Portions copied and modified from QPSK_TX_IQ_RX by JC and posted on
%Mathworks Dec 2005

%This module multplies the carrier frequency by the prn sequence and then
%demodulates the incoming signal by the composite signal (prn and carrier)
%The data is also filtered to produce the I-channel, Q-channel.

%The program aso displays the number of bit errors and the bit error rate
%in the command window at the end of the simulation.

%fcarr=20e3;         % Carrier frequency(Hz)
%N=1e3;		        % Number of data bits(bit rate)
fs=8*10e3;		    % Sampling frequency
Fn=fs/2;            % Nyquist frequency
Ts=1/fs;	        % Sampling time = 1/fs
T=1/N;		        % Bit time
randn('state',0);   % Keeps PRBS from changing on reruns
td=[0:Ts:(N*T)-Ts]';% Time vector(data)(transpose)

%prn and carrier
si=cs_t.*seq1;
sq=sn_t.*seq2;

%demodulate and despread
ibit=sumiq.*si;

qbit=sumiq.*sq;

combined=ibit+qbit;

%simple low pass filter

rc1=1/N;                    %time constant
ht1=(1/rc1).*exp(-tiq/rc1);%impulse response
flt_ibit=filter(ibit,1,ht1)/fs;
f5=length(flt_ibit)     %Verifies progress of program

%filter for q channel
rc=1/N;                 %time constant
ht=(1/rc).*exp(-tiq/rc);%impulse response
flt_qbit=filter(qbit,1,ht)/fs;
f6=length(flt_qbit)     %Verifies progress of program

bit1=sign(flt_ibit);%+/-1
bit2=sign(flt_qbit);%+/-1

%Second filtering stage (optional, but does give slight improvement)
bit3=filter(bit1,1,ht1)/fs;
f7=length(bit3)         %Verifies progress of program
bit4=filter(bit2,1,ht)/fs;
f8=length(bit4)         %Verifies progress of program

bit5=sign(bit3);%+/-1
bit6=sign(bit4);%+/-1

shaped_ibit=bit5 >0;%0 and 1
shaped_qbit=bit6 >0;%0 and 1

i_out=[shaped_ibit]';%transpose
q_out=[shaped_qbit]';%transpose

%compare received data to transmitted data by sampling at end of time
%cycle for data bit.  This compensates for delay created by filtering (plus
%more accurately simulates how processor would capture the data).
x=1;
for i=(fs/N-2):(fs/N):39998

    serial(x)=(i_out(i)+i_out(i+1)+i_out(i+2))/3;
    x=x+1;
    serial(x)=(q_out(i)+q_out(i+1)+q_out(i+2))/3;
    x=x+1;
end
data_1=data >0;
serial_1=serial >.5;
bit_errs_ss=xor(data_1,serial_1);
ss_bit_errors=sum(bit_errs_ss)
ss_ber=ss_bit_errors/N

%========================================================================
%Take FFT of modulated spread carrier
%========================================================================
y=sumiq;
NFFY=2.^(ceil(log(length(y))/log(2)));
FFTY=fft(y,NFFY);%pad with zeros
NumUniquePts=ceil((NFFY+1)/2);
FFTY=FFTY(1:NumUniquePts);
MY=abs(FFTY);
MY=MY*2;
MY(1)=MY(1)/2;
MY(length(MY))=MY(length(MY))/2;
MY=MY/length(y);
f1=(0:NumUniquePts-1)*2*Fn/NFFY;

%========================================================================
%Take FFT of despread combined signal (for demo purposes)
%========================================================================
y=combined;
NFFY=2.^(ceil(log(length(y))/log(2)));
FFTY=fft(y,NFFY);%pad with zeros
NumUniquePts=ceil((NFFY+1)/2);
FFTY=FFTY(1:NumUniquePts);
DY=abs(FFTY);
DY=DY*2;
DY(1)=DY(1)/2;
DY(length(DY))=DY(length(DY))/2;
DY=DY/length(y);
f1=(0:NumUniquePts-1)*2*Fn/NFFY;

%=====================================================================
%Plots
%======================================================================

figure(8)
subplot(4,1,1)
plot(tiq,ibit)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('I channel Output Waveform')

subplot(4,1,2)
plot(tiq,qbit)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Q channel Output Waveform')

subplot(4,1,3)
plot(tiq,si)
axis([0 5/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('i demod sig Waveform')

subplot(4,1,4)
plot(tiq,sq)
axis([0 5/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('q demod sig Waveform')

figure(9)

subplot(2,1,1)
plot(tiq,i_out)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('SS I Channel Data')

subplot(2,1,2)
plot(tiq,q_out)
axis([0 20/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('SS Q Channel Data')

figure(12)
subplot(2,1,1); plot(f1,MY);xlabel('');ylabel('AMPLITUDE');
axis([0 40000 0 1]);%zoom in/out
title('Frequency domain received');
grid on

subplot(2,1,2); plot(f1,20*log10(abs(MY).^2));xlabel('FREQUENCY(Hz)');ylabel('DB');
axis([0 40000 -80 0]);%zoom in/out
grid on
title('Received QPSK spread signal')

figure(12)
subplot(2,1,1); plot(f1,MY);xlabel('');ylabel('AMPLITUDE');
axis([0 40000 0 1]);%zoom in/out
title('Frequency domain received');
grid on

subplot(2,1,2); plot(f1,20*log10(abs(MY).^2));xlabel('FREQUENCY(Hz)');ylabel('DB');
axis([0 40000 -80 0]);%zoom in/out
grid on
title('Received QPSK spread signal')

figure(14)
subplot(2,1,1); plot(f1,DY);xlabel('');ylabel('AMPLITUDE');
axis([0 4*N 0 1]);%zoom in/out
title('Frequency domain despread combined signal');
grid on

subplot(2,1,2); plot(f1,20*log10(abs(DY).^2));xlabel('FREQUENCY(Hz)');ylabel('DB');
axis([0 4*N -80 0]);%zoom in/out
grid on
title('Received despread combined signal')
end

author: Sable

repo: mcbench-benchmarks


text:
% Generate exe file for Simu_Matlab

data       = coder.load('GEN_initData.mat');
N          = data.N;
xDim       = data.dim.x;
uDim       = data.dim.u;

% global variable
global ParNMPCGlobalVariable
globalVariable = {'ParNMPCGlobalVariable',coder.Constant(ParNMPCGlobalVariable)};
% config
cfg = coder.config('lib');
cfg.FilePartitionMethod = 'SingleFile';
cfg.TargetLang = 'C';
stackUsageMax = (xDim+uDim)*N/360*200000;
cfg.StackUsageMax = stackUsageMax;
cfg.BuildConfiguration = 'Faster Runs';
cfg.SupportNonFinite = false;
cfg.GenerateExampleMain = 'GenerateCodeAndCompile';
cfg.GenCodeOnly = false; % true to generate code only
myCCompiler = mex.getCompilerConfigurations(cfg.TargetLang,'Selected');
clear NMPC_Solve % must be cleared before code generation
if ~strcmp(myCCompiler.Manufacturer,'Microsoft')
    cfg.PostCodeGenCommand = 'buildInfo.addLinkFlags(''-fopenmp'')';
end
% generate exe
codegen -config cfg Simu_Matlab -globals globalVariable
%% Run
!Simu_Matlab

author: deng-haoyang

repo: ParNMPC


text:
function [model,rxnNames] = addDemandReaction(model,metaboliteNameList, printLevel)
% Adds demand reactions for a set of metabolites
% The reaction names for the demand reactions will be `DM_[metaboliteName]``
%
% USAGE:
%
%    model = addDemandReaction(model, metaboliteNameList)
%
% INPUTS:
%    model:                 COBRA model structure
%    metaboliteNameList:    List of metabolite names (cell array)
%    addMissingMetabolites  adds missing metabolites to model (default:
%                           true)
%
% OPTIONAL INPUT:
%    printLevel:            If > 0 will print out the reaction formulas (Default: 1).
%
% OUTPUTS:
%    model:                 COBRA model structure with added demand reactions
%    rxnNames:              List of added reactions
%
% .. Authors:
%       - Markus Herrgard 5/8/07
%       - Ines Thiele 03/09 - Corrected reaction coefficient for demand reaction
%       - Thomas Pfau, June 2018 - Change to use addMultipleReactions and adding printLevel

if (~iscell(metaboliteNameList))
    metaboliteNameList = {metaboliteNameList};
end

if nargin < 3 %No PrintLevel
    printLevel = 0;
end

missingMets = setdiff(metaboliteNameList,model.mets);
if ~isempty(missingMets)
    warning('The following Metabolites have been added to the model, as they were not in the model before:\n%s', strjoin(missingMets,'\n'));
    model = addMultipleMetabolites(model,missingMets);
end
rxnNames = rowvector(strcat('DM_',metaboliteNameList));
nMets = length(metaboliteNameList);
stoich = -1 * speye(nMets);
subSystems = repmat({'Demand'},nMets,1);
ubs = repmat(1000,nMets,1);
lbs = zeros(nMets,1);
model = addMultipleReactions(model,rxnNames,metaboliteNameList,stoich,...
    'lb',lbs,'ub',ubs,'subSystems',subSystems,...
    'printLevel',printLevel);
end

author: opencobra

repo: cobratoolbox


text:
function [correlations, rms_errors, patch_experts, visi_index, centres, normalisation_options] ...
    = TrainPatchExperts(trainingLoc, view, scale, sigma, patch_type, ratio_neg, num_samples, varargin)

    % Training patch experts for all of the landmarks for a particular view
    % and scale

    % Get the training parameters
    normalisation_options = Parse_settings( sigma, patch_type, ratio_neg, num_samples, varargin{:});

    % Load the training data
    if(sum(strcmp(varargin,'data_loc')))
        ind = find(strcmp(varargin,'data_loc')) + 1;
        data_loc = varargin{ind};
        data_loc = sprintf(['%s/' data_loc '%s_%s.mat'], trainingLoc, num2str(scale), num2str(view));
    else
        data_loc = sprintf('%s/wild_%s_%s.mat', trainingLoc, num2str(scale), num2str(view));
    end

    load(data_loc);
    examples = all_images;
    visi_index = visiIndex;

    clear 'all_images'

    numPoints = size(landmark_locations,2);

    correlations = zeros(1, numPoints);
    rms_errors = zeros(1, numPoints);

    patch_experts = zeros(1, numPoints, normalisation_options.patchSize(1)*normalisation_options.patchSize(2)+2);

    for j=1:numPoints

        % can only do mirroring if there is no yaw
        if((numPoints == 68 || numPoints == 29 )&& centres(2) == 0)
            % Do not redo a mirror feature (just flip them)
            if(numPoints == 68)
                mirrorInds = [1,17;2,16;3,15;4,14;5,13;6,12;7,11;8,10;18,27;19,26;20,25;21,24;22,23;...
                      32,36;33,35;37,46;38,45;39,44;40,43;41,48;42,47;49,55;50,54;51,53;60,56;59,57;...
                      61,65;62,64;68,66];
            else
                mirrorInds = [1,2; 3,4; 5,7; 6,8; 9,10; 11,12; 13,15; 14,16; 17,18; 19,20; 23,24];
            end
            mirror_idx = j;
            if(any(mirrorInds(:,1)==j))
                mirror_idx = mirrorInds(mirrorInds(:,1)==j,2);
            elseif(any(mirrorInds(:,2)==j))
                mirror_idx = mirrorInds(mirrorInds(:,2)==j,1);
            end

            if(mirror_idx~=j & correlations(1,mirror_idx) ~= 0)

                correlations(1,j) = correlations(1,mirror_idx);
                rms_errors(1, j) = rms_errors(1,mirror_idx);

                % Bias and scaling are the same
                patch_experts(1, j, 1:2) = patch_experts(1, mirror_idx, 1:2);

                % flip the support weights
                w = patch_experts(1, mirror_idx, 3:end);

                w = reshape(w, normalisation_options.patchSize(1), normalisation_options.patchSize(2));
                w = fliplr(w);
                w = reshape(w, normalisation_options.patchSize(1)*normalisation_options.patchSize(2),1);

                patch_experts(1, j, 3:end) = w;

                fprintf('Feature %d done\n', j);

                continue;

            end
        end

        if(visi_index(j))
            % instead of loading the patches compute them here:
            numSamples = normalisation_options.numSamples;

            [curr_examples, currentLabels, unnormed_samples] = ExtractTrainingSamples(examples, landmark_locations, sigma, normalisation_options.normalisationRegion, normalisation_options.patchSize, numSamples, j, normalisation_options);

            [ patchExpert, correlation, rmsError ] = CreatePatchExpert( curr_examples, currentLabels, unnormed_samples, normalisation_options);

            correlations(1,j) = correlation;
            rms_errors(1, j) = rmsError;
            patch_experts(1, j, :) = patchExpert(:);

            fprintf('Feature %d done\n', j);
        end
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function u = tanh(a)
%TANH         Slope hyperbolic tangent tanh(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  u = a;

  u.r = tanh(a.r);
  indexc = 1:INTLAB_SLOPE.NUMVAR;
  indexr = 2:INTLAB_SLOPE.NUMVAR+1;
  Xxs = hull(a.r(:,indexc),a.r(:,indexr));
  Index = 1:size(a.r.inf,1);

  index = all( a.r.sup<=0 , 2);
  if any(index)
    aindex.r = a.r(index,:);
    aindex.s = a.s(index,:);
    u.s(index,:) = slopeconvexconcave('tanh','1./sqr(cosh(%))',aindex,1);
    Index(index) = 0;
  end

  index = all( a.r.inf>=0 , 2);
  if any(index)
    aindex.r = a.r(index,:);
    aindex.s = a.s(index,:);
    u.s(index,:) = slopeconvexconcave('tanh','1./sqr(cosh(%))',aindex,0);
    Index(index) = 0;
  end

  if any(Index)
    Index( Index==0 ) = [];
    u.s(Index,:) = a.s(Index,:) ./ sqr(cosh(Xxs(Index)));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%% Nonconvex QP
clc
%Objective & Constraints
H = sparse([0 -1; -1 0]);
f = [0;0];
%Constraints
lb = [-0.5;-0.5];
ub = [1;1];

%Build & Solve
Opt = opti('qp',H,f,'bounds',lb,ub,'options',optiset('solver','ooqp','display','final','warnings','critical'))
[x,fval,exitflag,info] = multisolve(Opt)
%Plot
figure(1)
plot(Opt)
figure(2)
multiplot(Opt)

%% Nonconvex QP w Inequalities and Maximization
clc
%Objective & Constraints
H = sparse([0 -1; -1 0]);
f = [0;0];
%Constraints
A = [1 1; 1 -1]; b = [1.5;1.2];
lb = [-0.5;-0.5];
ub = [1;1];

%Build & Solve
Opt = opti('qp',H,f,'ineq',A,b,'bounds',lb,ub,'sense',-1,'options',optiset('solver','ooqp','display','iter','warnings','all'))
[x,fval,exitflag,info] = multisolve(Opt)
%Plot
figure(1)
plot(Opt)
figure(2)
multiplot(Opt)

%% Nonconvex QP2
clc
%Objective & Constraints
H = sparse([0 0 -1;0 -1 0;-1 0 0]);
f = [0;0;0];
%Constraints
lb = [0;0;0]; ub = [4;4;4];

%Build & Solve
Opt = opti('qp',H,f,'bounds',lb,ub,'options',optiset('solver','auto','display','final'))
[x,fval,exitflag,info] = multisolve(Opt)
multiplot(Opt)

%% NLS1
clc
%Function
fun = @(x,xdata) x(1)*exp(x(2)*xdata);
grad = @(x,xdata) [ exp(x(2).*xdata), x(1).*xdata.*exp(x(2).*xdata)];
%Fitting Data
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
%Setup Options
opts = optiset('solver','mkltrnls','display','final');
%Build & Solve
Opt = opti('fun',fun,'grad',grad,'data',xdata,ydata,'bounds',[400;-1],[600;1],'options',opts)
[x,fval,exitflag,info] = multisolve(Opt)
%Plot
figure(1)
plot(Opt)
figure(2)
multiplot(Opt,1)

%% NLS2
clc
%Function
i = (1:40)';
fun = @(x) x(1)*exp(-x(2)*i) + x(3);
%Fitting Data
ydata=[5.8728, 5.4948, 5.0081, 4.5929, 4.3574, 4.1198, 3.6843, 3.3642, 2.9742, 3.0237, 2.7002, 2.8781,...
       2.5144, 2.4432, 2.2894, 2.0938, 1.9265, 2.1271, 1.8387, 1.7791, 1.6686, 1.6232, 1.571, 1.6057,...
       1.3825, 1.5087, 1.3624, 1.4206, 1.2097, 1.3129, 1.131, 1.306, 1.2008, 1.3469, 1.1837, 1.2102,...
       0.96518, 1.2129, 1.2003, 1.0743];
%Setup Options
opts = optiset('solver','mkltrnls','display','final');
%Build & Solve
x0=[1.0; 0.0; 0.0];
Opt = opti('fun',fun,'ydata',ydata,'bounds',[0;0;0],[10;5;5],'options',opts,'x0',x0)
[xs,fvals] = multisolve(Opt,x0)

plot(Opt)

%% NLP1
clc
% Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
% Linear Constraints
A = [-1 1];  b = -1;
Aeq = [1.1 1]; beq = 5;
lb = [0;0]; ub = [4;4];
% Initial Guess
x0 = [2;2];

% Build OPTI Object
Opt = opti('obj',obj,'ndec',2,'bounds',lb,ub,'ineq',A,b,'eq',Aeq,beq,'options',optiset('display','final'))
% Solve with Violation Penalty of 1e2
[x,fval,ef,info] = multisolve(Opt,x0)
%Plot
multiplot(Opt,1)

%% NLP2
clc
% Objective
 fun = @(x) norm([x(1) - sin(2*x(1) + 3*x(2)) - cos(3*x(1) - 5*x(2));
                  x(2) - sin(x(1) - 2*x(2)) + cos(x(1) + 3*x(2))]);

% Bounds
 lb = [-4;-4];
 ub = [4;4];

% Initial Solution
 x0 = [-4;-4];

% Build optiprob structure (intermediate structure)
 Opt = opti('fun',fun,'bounds',lb,ub,'x0',x0,'options',optiset('solver','filterSD','display','iter'));

 [x,fval,ef,info] = multisolve(Opt,x0)
%Plot
figure(1)
plot(Opt,[],1)
figure(2)
multiplot(Opt,1)

%% NLP3
 % Objective
 fun = @(x) 20 + x(1)^2 + x(2)^2 - 10*(cos(2*pi*x(1)) + cos(2*pi*x(2)));

 A = [1 1]; b = -1.1;

% Bounds
 lb = [-5;-5];
 ub = [5;5];

%Initial Guess
 x0 = [0;0];

 % Build optiprob structure (intermediate structure)
 Opt = opti('fun',fun,'ineq',A,b,'bounds',lb,ub,'x0',x0,'options',optiset('solver','filterSD','display','iter'));

 [x,fval,ef,info] = multisolve(Opt,x0)
%Plot
figure(1)
plot(Opt,[-2 0 -2 0])
figure(2)
multiplot(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [U,output] = cpd3_sd(T,U0,options)
%CPD3_SD CPD by simultaneous diagonalization.
%   [U,output] = cpd3_sd(T,U0) computes the factor matrices U{1}, U{2} and
%   U{3} belonging to a canonical polyadic decomposition of the third-order
%   tensor T. The algorithm is initialized with the factor matrices U0{n}.
%   The number of rank-one terms R should satisfy 1 <= R <= max(size(T)).
%   The structure output returns additional information:
%
%      output.SDName   - The name of the selected simultaneous
%                        diagonalization algorithm.
%      output.SDOutput - Information returned by the selected simultaneous
%                        diagonalization algorithm.
%
%   cpd3_sd(T,U0,options) may be used to set the following options:
%
%      options.SD =            - A function handle to the desired
%      [@cpd_als|@cpd_minf|...   simultaneous diagonalization algorithm.
%       @cpd_nls|{@cpd3_sgsd}]
%      options.SDOptions       - An options structure passed to the
%                                selected SD algorithm. See also
%                                help [options.SD].

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. De Lathauwer, "A Link between the Canonical Decomposition in
%       Multilinear Algebra and Simultaneous Matrix Diagonalization,"
%       SIAM J. Matrix Anal. Appl., Vol. 28, 2006, pp. 642-666.

% Incomplete/sparse tensors not supported yet.
T = fmt(T,true);
if isstruct(T)
    error('This method does not yet support incomplete/sparse tensors.');
end

% Check the tensor T.
N = ndims(T);
if N ~= 3, error('cpd3_sd:T','ndims(T) should be 3.'); end

% Check the initial factor matrices U0.
U0 = U0(:).';
R = size(U0{1},2);
size_tens = size(T);
if any(cellfun('size',U0,2) ~= R)
    error('cpd3_sd:U0','size(U0{n},2) should be the same for all n.');
end
if R > max(size_tens)
    error('cpd3_sd:U0','size(U0{n},2) should be <= max(size(T)).');
end
if any(cellfun('size',U0,1) ~= size_tens)
    error('cpd3_sd:U0','size(T,n) should equal size(U0{n},1).');
end

% Check the options structure.
xsfunc = @(f)isa(f,'function_handle')&&exist(func2str(f),'file');
if nargin < 3, options = struct; end
if ~isfield(options,'SD')
    options.SD = @cpd3_sgsd;
elseif ~xsfunc(options.SD)
    error('cpd3_sd:SD','Algorithm does not exist.');
end

% Prepermute the tensor so the longest mode is the first one, and for
% efficiency, the last one is the shortest one.
[size_tens,idx] = sort(size_tens,'descend');
if any(idx ~= 1:3)
    T = permute(T,idx);
    U0 = U0(idx);
    iperm(idx) = 1:N;
end

% Best rank-R approximation of the mode-1 unfolding of T.
% T(1) = U{1}*kr(U([3 2]))^T = F*W^-H * W^H*E^H
% U{1} = F*W^-H; kr(U([3 2])) = conj(E*W)
T1 = reshape(T,size_tens(1),[]);
[~,~,E] = svd(T1,'econ');
E = E(:,1:R);

% Compute the Gramian P'*P for the rank-one detecting device P. The Gramian
% is computed as a fourth-order tensor in which element with indices
% (r,s,t,u) is equal to <Er,Et>*<Es,Eu> + <Er,Eu>*<Es,Et> - <Eu'*Es,Er'*Et>
% - <Et'*Es,Er'*Eu>, where En denotes the matricized n-th column of E.
% First compute all products P34(r,s,t,u) = conj(<Er'*Es,Et'*Eu>).
EE = reshape(E,size_tens(2),[]);
P34 = zeros(size_tens(3)^2,R^2);
for r = 1:R
    X = EE(:,(r-1)*size_tens(3)+(1:size_tens(3)))'*EE;
    P34(:,(r-1)*R+(1:R)) = reshape(X,size_tens(3)^2,R);
end
P34 = reshape(P34'*P34,[R R R R]);
% Compute all four terms P1 to P4.
rr = false(1,R*(R+1)/2);
rr([1 1+cumsum(R:-1:2)]) = true;
P12 = ones(R*(R+1)/2,1);
P12(rr) = 2;
P12 = diag(P12);
P34 = -reshape(permute(P34,[1 4 2 3]),[R^2 R^2]) ...
      -reshape(permute(P34,[1 4 3 2]),[R^2 R^2]);
% Compose the matrix P'*P, without redundant columns/rows.
lt = tril(true(R,R));
PP = P12+P34(lt(:),lt(:));

% Compute the kernel of the adjusted Gramian and build the tensor M,
% consisting of symmetric matrices Mk = W Lk W^T.
[~,~,K] = svd(PP,'econ');
K = K(:,end:-1:end-R+1);
K(~rr,:) = 0.5*K(~rr,:);
M = zeros(R,R,R);
Mlt = zeros(R,R);
for r = 1:R
    Mlt(lt) = K(:,r);
    M(:,:,r) = Mlt+tril(Mlt,-1).';
end

% Simultaneously diagonalize the frontal slices of M.
W0 = E\conj(kr(U0(end:-1:2)));
C0 = reshape(M,[R*R R]).'*conj(kr(W0,W0))/conj((W0'*W0).^2);
if ~isfield(options,'SDOptions'), options.SDOptions = struct; end
[W,output.SDOutput] = options.SD(M,{W0,W0,C0},options.SDOptions);
output.SDName = func2str(options.SD);

% Reconstruct the factor matrices.
U = cell(1,N);
krCB = conj(E*W{1});
for r = 1:R
    [u,s,v] = svd(reshape(krCB(:,r),size_tens([2 3])),'econ');
    U{2} = [U{2} s(1)*u(:,1)];
    U{3} = [U{3} conj(v(:,1))];
end
U{1} = T1*conj(kr(U([3 2])))/conj((U{3}'*U{3}).*(U{2}'*U{2}));

% Inverse permute the factor matrices.
if exist('iperm','var')
    U = U(iperm);
end

author: andrewssobral

repo: lrslibrary


text:
function[varargout]=trajunwrap(varargin)
%TRAJUNWRAP  Unwraps Lagrangian trajectories from a periodic domain.
%
%   CX=TRAJUNWRAP(CX,L) unwraps the complex-valued Lagrangian trajectory CX
%   from within a doubly periodic domain having edges at +/- L.
%
%   CX specifies the partile position in complex-valued form CX=X+iY.
%
%   CX may be of any dimension, but time is assumed to increase along rows.
%   The output array will be the same size as the input array.
%
%   Set L=[LX LY] to specify different periodicities in the X and Y
%   directions. Use L=[LX NaN] for only zonal periodicity, or L=[NaN Ly]
%   for only meridional periodicity.
%
%   [CX1,CX2,...,CXN]=TRAJUNWRAP(CX1,CX2,...CXN,L) also works for multiple
%   arrays of trajectories with the same periodicity.
%
%   TRAJUNWRAP(CX1,CX2,...CXN,L); with no output arguments overwrites the
%   original named input variables.
%
%   See also TRAJWRAP.
%
%   Usage: cx=trajunwrap(cx,L);
%          cx=trajunwrap(cx,[Lx Ly],'nan');
%          cx=trajunwrap(cx,[Lx Ly]);
%          [cx1,cx2]=trajunwrap(cx1,cx2,[Lx Ly]);
%          trajunwrap(cx1,cx2,[Lx Ly]);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2014--2015 J.M. Lilly --- type 'help jlab_license' for details

na=nargin;
L=varargin{end};
if length(L)==1
    Lx=L;
    Ly=L;
else
   Lx=L(1);
   Ly=L(2);
end

for i=1:na-1
    z=varargin{i};
    if iscell(z)
        for j=1:length(z)
            varargout{i}{j}=trajunwrap_one(Lx,Ly,z{j});
        end
    else
       varargout{i}=trajunwrap_one(Lx,Ly,z);
    end
end

eval(to_overwrite(na-1));

function[zw]=trajunwrap_one(Lx,Ly,z)
x=real(z);
y=imag(z);
if ~isnan(Lx)
    x=frac(Lx,pi)*unwrap(frac(pi,Lx)*real(z),[],1);
end
if ~isnan(Ly)
    y=frac(Ly,pi)*unwrap(frac(pi,Lx)*imag(z),[],1);
end

zw=x+sqrt(-1)*y;

% if ~isempty(strfind(str,'nan'))
%     zwshift=vshift(zw,1,1);
%     bool=(abs(real(zw-zwshift))>Lx/2)|(abs(imag(zw-zwshift))>Ly/2);
%     bool(1,:,:,:,:,:,:)=true;
%     bool(end,:,:,:,:,:,:)=true;
%     zw(bool)=nan;
% end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function H = freqwin(name,L,bw,varargin)
%FREQWIN Frequency response window
%   Usage: H = freqwin(name,L,bw);
%
%   `freqwin(name,L,bw)` returns a frequency window *name* of length *L*
%   with the mainlobe -6dB (half height) bandwidth *bw*. It is intended to
%   represent frequency response of a band-pass filter/window with
%   bandwidth *bw*. The bandwidth is given in frequencies normalized
%   to the Nyquist rate.
%
%   The function is not periodically wrapped should it be nonzero outside
%   of the *L* samples (as opposed to e.g. |pgauss|).
%
%   The following windows can be used in place of *name*:
%
%     'gauss'        Gaussian window
%
%     'gammatone'    Gammatone window. The default order is 4. Custom order
%                    can be set by `{'gammatone',order}`.
%
%     'butterworth'  Butterworth window. The default order is 4. Custom
%                    order can be set by `{'butterworth',order}`.
%
%   `freqwin(name,L,bw,fs)` does the same as above except *bw* is expected
%   to be in Hz given sampling frequency *fs*.
%
%   `freqwin` understands the following key-value pairs and flags at the end of
%   the list of input parameters:
%
%     'fs',fs      If the sampling frequency *fs* is specified then the *bw*
%                  is expected to be in Hz.
%
%     'shift',s    Shift the window by $s$ samples. The value can be a
%                  fractional number.
%
%     'wp'         Output is whole point even. This is the default. It
%                  corresponds to a shift of $s=0$.
%
%     'hp'         Output is half point even, as most Matlab filter
%                  routines. This corresponds to a shift of $s=-.5$
%
%   Additionally, the function accepts flags to normalize the output. Please see
%   the help of |setnorm|. Default is to use `'peak'` normalization.
%
%
%   See also: firwin, setnorm, plotfft

% AUTHORS: Nicki Holighaus

complainif_notenoughargs(nargin,3,upper(mfilename));

if ~isscalar(L)
    error('%s: L must be a scalar',upper(mfilename));
end

if ~isscalar(bw)
    error('%s: bw must be a scalar',upper(mfilename));
end

freqwintypes = arg_freqwin(struct);
freqwintypes = freqwintypes.flags.wintype;

if ~iscell(name), name = {name}; end

if ~ischar(name{1}) || ~any(strcmpi(name{1},freqwintypes))
  error('%s: First input argument must the name of a supported window.',...
        upper(mfilename));
end;

winArgs = name(2:end);
winName = lower(name{1});

definput.import={'setnorm'};
definput.importdefaults={'null'};
definput.flags.centering={'wp','hp','shift'};
definput.keyvals.shift = 0;
definput.keyvals.fs = 2;
definput.keyvals.atheight = 10^(-3/10);
[flags,kv,fs]=ltfatarghelper({'fs'},definput,varargin,'freqwin');

if flags.do_wp, kv.shift=0; end
if flags.do_hp, kv.shift=0.5; end

if ( kv.shift >= .5 || kv.shift < -.5 )
    error('%s: Parameter shift must be in ]-.5,.5].',upper(mfilename));
end

if ( bw > fs || bw < eps )
     error('%s: Parameter bw must be in ]0,fs].',upper(mfilename));
end

step = fs/L;
bwrelheight = kv.atheight;

H = (-kv.shift+[0:1:ceil(L/2)-1,-floor(L/2):-1]');

switch winName
    case 'gauss'
        H = exp(4*H.^2*log(bwrelheight)/(bw/step)^2);

    case 'butterworth'
        definputbutter.keyvals.order=4;
        [~,~,order]=ltfatarghelper({'order'},definputbutter,winArgs);
        H = 1./(sqrt(1 + (H/(bw/step/2)).^(2*order)));

    case 'gammatone'
        definputgamma.keyvals.order=4;
        [~,~,order]=ltfatarghelper({'order'},definputgamma,winArgs);
        if order <= 1
            error('%s: order must be bigger than 1. (passed %.2f)',...
                  upper(mfilename),order);
        end

        gtInverse = @(yn) sqrt(yn^(-2/order)-1);
        dilation = bw/2/gtInverse(bwrelheight)/step;
        peakpos = (order-1)/(2*pi*dilation);
        peakmod = exp(2*pi*1i*H*peakpos);
        H = (1+1i*H/dilation).^(-order).*peakmod;
    otherwise
        error('%s: SENTINEL. Unknown window.',upper(mfilename));
end

H=setnorm(H,flags.norm);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function momentspace = tsss_config_momentspace
% configuration file for cropping
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: tsss_config_momentspace.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
% D
%--------------------------------------------------------------------------
D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

%--------------------------------------------------------------------------
% condthresh
%--------------------------------------------------------------------------
condthresh         = cfg_entry;
condthresh.tag     = 'condthresh';
condthresh.name    = 'Condition number threshold';
condthresh.help    = {'Threshold on condition number applied for basis regularisation'};
condthresh.strtype = 'r';
condthresh.num     = [1 1];
condthresh.val     = {80};

%--------------------------------------------------------------------------
% none
%--------------------------------------------------------------------------

none = cfg_const;
none.tag = 'none';
none.name = 'None';
none.val  = {0};

%--------------------------------------------------------------------------
% addchannels
%--------------------------------------------------------------------------

addchannels      = cfg_choice;
addchannels.tag  = 'addchannels';
addchannels.name = 'Extra channels to add';
addchannels.values  = {none, spm_cfg_eeg_channel_selector};
addchannels.val  = {none};

%--------------------------------------------------------------------------
% momentspace
%--------------------------------------------------------------------------
momentspace          = cfg_exbranch;
momentspace.tag      = 'momentspace';
momentspace.name     = 'TSSS space conversion';
momentspace.val      = {D, condthresh, addchannels};
momentspace.help     = {'Engage the TSSS space virtual montage'}';
momentspace.prog     = @eeg_momentspace;
momentspace.vout     = @vout_eeg_momentspace;
momentspace.modality = {'EEG'};

%==========================================================================
function out = eeg_momentspace(job)
% construct the S struct
S   = [];
S.D = char(job.D);
S.condthresh = job.condthresh;

D = spm_eeg_load(S.D);

if isfield(job.addchannels, 'channels')
    S.addchannels = D.chanlabels(D.selectchannels(spm_cfg_eeg_channel_selector(job.addchannels.channels)));
else
    S.addchannels = {};
end

out.D        = tsss_spm_momentspace(S);
out.Dfname   = {fullfile(out.D)};

%==========================================================================
function dep = vout_eeg_momentspace(job)
% return dependencies
dep(1)            = cfg_dep;
dep(1).sname      = 'Moment space MEG data';
dep(1).src_output = substruct('.','D');
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'Moment space MEG datafile';
dep(2).src_output = substruct('.','Dfname');
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function pass = test_isPeriodicTech(pref)
% Test isPeriodicTech().

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

f = chebfun3(@(x,y,z) x);
pass(1) = ~isPeriodicTech(f);

f = chebfun3(@(x,y,z) cos(pi*x), 'trig');
pass(2) = isPeriodicTech(f);

% Test for complex function values:
f = chebfun3(@(x,y,z) exp(1i*pi*x), 'trig');
pass(3) = isPeriodicTech(f);

% On different domain:
f = chebfun3(@(x,y,z) cos(x).*sin(z), [ -pi, pi, -pi, pi, -pi, pi ], 'trig');
pass(4) = isPeriodicTech(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function C = train_ridgeRegressionShrink(xTr, yTr, varargin)
% TRAIN_RIDGEREGRESSIONSHRINK - Regularized Regression with automatic shrinkage selection
%
%Synopsis:
%   C = train_ridgeRegressionShrink(XTR, YTR)
%   C = train_ridgeRegressionShrink(XTR, YTR, OPTS)
%
%Arguments:
%   XTR: DOUBLE [NxM] - Data matrix, with N feature dimensions, and M training points/examples.
%   YTR: INT [CxM] - Regression targets. C by M matrix of regression
%                     labels, with C representing the number of target signals and M the number of training examples/points.
%   OPT: PROPLIST - Structure or property/value list of optional
%                   properties. Options are also passed to clsutil_shrinkage.
%     'ExcludeInfs' - BOOL (default 0): If true, training data points with value 'inf' are excluded from XTR
%
%Returns:
%   C: STRUCT - Trained regressor structure, with the hyperplane given by
%               fields C.w and C.b.  C includes the fields:
%    'w' : weight matrix
%    'b' : FLOAT bias
%
%Description:
%   TRAIN_RLDA trains a regularized linear regressionon data X with
%   regression targets given in YTR. The shrinkage parameter is selected by the
%   function clsutil_shrinkage.
%
%
%Examples:
%   train_RLDA(X, regression_target)
%   train_RLDA(X, regression_target, 'Target', 'D')
%
%See also:
%   APPLY_SEPARATINGHYPERPLANE, CLSUTIL_SHRINKAGE,
%   TRAIN_LDA, TRAIN_RDAREJECT

% Sven Daehne

props= {'ExcludeInfs'      0                             'BOOL'
       };

% get props list of the subfunction
props_shrinkage= clsutil_shrinkage;

if nargin==0,
  C= opt_catProps(props, props_shrinkage);
  return
end

opt = opt_proplistToStruct(varargin{:});
opt = opt_setDefaults(opt, props);
opt_checkProplist(opt, props, props_shrinkage);

if opt.ExcludeInfs,
  ind = find(sum(abs(xTr),1)==inf);
  xTr(:,ind) = [];
  yTr(:,ind) = [];
end

% compute the regularized covariance matrix
opt_shrinkage = opt_substruct(opt, props_shrinkage(:,1));
[C_cov, C.gamma] = clsutil_shrinkage(xTr, opt_shrinkage);

% compute the weight vector(s)
C_cov = C_cov * size(xTr,2);
Xn = xTr - repmat(mean(xTr,2), [1 size(xTr,2)]);
C.w = C_cov \ Xn * yTr';

% compute the bias term(s)
mu_x = mean(xTr,2);
mu_y = mean(yTr,2);
C.b = mu_y - C.w'*mu_x;

author: bbci


text:
%BDRNEU Compute Neumann (flux) boundary condition vector.
%
%   [ F, T ] = BDRNEU( PROB, I_CUB ) Computes a load vector F with
%   Neumann/flux boundary condition contributions for the boundary
%   expressions defined in PROB.BDR.N.
%
%   Boundary coefficients are specified as (1 x n_dvar) cell arrays in
%   the BDR.N field where the entry for each boundary is a (n_bc_groups
%   x n_bdr) nested cell array containing the coefficients
%   (n_bc_groups is 1, except for special element types such as
%   Hermite basis functions). The coefficient entries can be specified
%   either as constant numeric values, or string expressions which
%   will be evaluated during the simulation.
%
%       Input       Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       prob        struct                 Problem data struct
%       i_cub       scalar {2}             Numerical integration rule
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       f           (neq,1)                Neumann boundary condition load vector
%       t           scalar                 Time spent in function

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function hrOut = getLVD_HR(hText)
%getMA_HR Summary of this function goes here
%   Detailed explanation goes here
    try
        hr = {repmat('# ',1,250)};
        hrOut = textwrap(hText,hr);
        hrOut = strtrim(hrOut{1});
        hrOut = strrep(hrOut,' ','#');
        hrOut = hrOut(1:end-4);
    catch
        hrOut = '##################################################################################';
    end
end


text:
function linplus_test06 ( )

%*****************************************************************************80
%
%% TEST06 tests R83_NP_FA, R83_NP_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  For a tridiagonal matrix that can be factored\n' );
  fprintf ( 1, '    with no pivoting,\n' );
  fprintf ( 1, '  R83_NP_FA factors;\n' );
  fprintf ( 1, '  R83_NP_SL solves a factored system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  a = r83_random ( n, seed );

  r83_print ( n, a, '  The tridiagonal matrix:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r83_mxv ( n, a, x );
  x(1:n) = 0.0E+00;
%
%  Factor the matrix.
%
  [ a_lu, info ] = r83_np_fa ( n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST06 - Fatal error!\n' );
    fprintf ( 1, '  The test matrix is singular.\n' );
    return;
  end
%
%  Solve the linear system.
%
  job = 0;
  x = r83_np_sl ( n, a_lu, b, job );

  r8vec_print ( n, x, '  Solution:' );
%
%  Set the desired solution
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side, using the factored matrix.
%
  job = 1;
  b = r83_np_ml ( n, a_lu, x, job );
%
%  Solve the linear system.
%
  job = 1;
  x = r83_np_sl ( n, a_lu, b, job );

  r8vec_print ( n, x, '  Solution to tranposed system:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [mri, hdr] = read_ctf_mri(filename);

% READ_CTF_MRI reads header and image data from CTF format MRI file
%
% [mri, hdr] = read_ctf_mri(filename)
%
% See also READ_CTF_MEG4, READ_CTF_RES4

% Copyright (C) 2003 Robert Oostenveld
%
fid = fopen(filename,'rb', 'ieee-be');

if fid<=0
  error(sprintf('could not open MRI file: %s\n', filename));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ THE IMAGE HEADER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

warning off
% general header information
hdr.identifierString = setstr(fread(fid,32,'char'))'; % CTF_MRI_FORMAT VER 2.2
hdr.imageSize = fread(fid,1,'int16'); % always = 256
hdr.dataSize = fread(fid,1,'int16'); % 1 or 2(bytes)
hdr.clippingRange = fread(fid,1,'int16'); % max.integer value of data
hdr.imageOrientation = fread(fid,1,'int16'); % eg., 0 = left on left, 1 = left on right
hdr.mmPerPixel_sagittal = fread(fid,1,'float'); % voxel dimensions in mm
hdr.mmPerPixel_coronal = fread(fid,1,'float'); % voxel dimensions in mm
hdr.mmPerPixel_axial = fread(fid,1,'float'); % voxel dimensions in mm

% HeadModel_Info specific header items
hdr.HeadModel.Nasion_Sag = fread(fid,1,'int16'); % fid.point coordinate(in voxels) for nasion - sagittal
hdr.HeadModel.Nasion_Cor = fread(fid,1,'int16'); % nasion - coronal
hdr.HeadModel.Nasion_Axi = fread(fid,1,'int16'); % nasion - axial
hdr.HeadModel.LeftEar_Sag = fread(fid,1,'int16'); % left ear - sagittal
hdr.HeadModel.LeftEar_Cor = fread(fid,1,'int16'); % left ear - coronal
hdr.HeadModel.LeftEar_Axi = fread(fid,1,'int16'); % left ear - axial
hdr.HeadModel.RightEar_Sag = fread(fid,1,'int16'); % right ear - sagittal
hdr.HeadModel.RightEar_Cor = fread(fid,1,'int16'); % right ear - coronal
hdr.HeadModel.RightEar_Axi = fread(fid,1,'int16'); % right ear - axial
fread(fid,2,'char'); % padding to 4 byte boundary
hdr.HeadModel.defaultSphereX = fread(fid,1,'float'); % sphere origin x coordinate(in mm)
hdr.HeadModel.defaultSphereY = fread(fid,1,'float'); % sphere origin y coordinate(in mm)
hdr.HeadModel.defaultSphereZ = fread(fid,1,'float'); % sphere origin z coordinate(in mm)
hdr.HeadModel.defaultSphereRadius = fread(fid,1,'float'); % default sphere radius(in mm)

% Image_Info specific header items
hdr.Image.modality = fread(fid,1,'int16'); % 0 = MRI, 1 = CT, 2 = PET, 3 = SPECT, 4 = OTHER
hdr.Image.manufacturerName = setstr(fread(fid,64,'char'))';
hdr.Image.instituteName = setstr(fread(fid,64,'char'))';
hdr.Image.patientID = setstr(fread(fid,32,'char'))';
hdr.Image.dateAndTime = setstr(fread(fid,32,'char'))';
hdr.Image.scanType = setstr(fread(fid,32,'char'))';
hdr.Image.contrastAgent = setstr(fread(fid,32,'char'))';
hdr.Image.imagedNucleus = setstr(fread(fid,32,'char'))';

fread(fid,2,'char'); % padding to 4 byte boundary

hdr.Image.Frequency = fread(fid,1,'float');
hdr.Image.FieldStrength = fread(fid,1,'float');
hdr.Image.EchoTime = fread(fid,1,'float');
hdr.Image.RepetitionTime = fread(fid,1,'float');
hdr.Image.InversionTime = fread(fid,1,'float');
hdr.Image.FlipAngle = fread(fid,1,'float');
hdr.Image.NoExcitations = fread(fid,1,'int16');
hdr.Image.NoAcquisitions = fread(fid,1,'int16');
hdr.Image.commentString = setstr(fread(fid,256,'char'))';
hdr.Image.forFutureUse = setstr(fread(fid,64,'char'))';

% continuation general header
hdr.headOrigin_sagittal = fread(fid,1,'float'); % voxel location of head origin
hdr.headOrigin_coronal = fread(fid,1,'float'); % voxel location of head origin
hdr.headOrigin_axial = fread(fid,1,'float'); % voxel location of head origin
% euler angles to align MR to head coordinate system(angles in degrees !)
hdr.rotate_coronal = fread(fid,1,'float'); % 1. rotate in coronal plane by this angle
hdr.rotate_sagittal = fread(fid,1,'float'); % 2. rotate in sagittal plane by this angle
hdr.rotate_axial = fread(fid,1,'float'); % 3. rotate in axial plane by this angle
hdr.orthogonalFlag = fread(fid,1,'int16'); % if set then image is orthogonal
hdr.interpolatedFlag = fread(fid,1,'int16'); % if set than image was interpolated
hdr.originalSliceThickness = fread(fid,1,'float'); % original spacing between slices before interpolation
transformMatrix = fread(fid,[4 4],'float')', % transformation matrix head->MRI[column][row]
fread(fid,204,'char'); % unused, padding to 1028 bytes
warning on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ THE IMAGE DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if hdr.dataSize==1
  mri = uint8(fread(fid, 256*256*256, 'uint8'));
elseif hdr.dataSize==2
  mri = uint16(fread(fid, 256*256*256, 'uint16'));
else
  error('unknown datasize in CTF mri file');
end
mri = reshape(mri, [256 256 256]);
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO POST-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% reorient the image data to obtain corresponding image data and transformation matrix
mri = permute(mri, [3 1 2]);	% this was determined by trial and error

% reorient the image data and the transformation matrix along the left-right direction
% remember that the fiducials in voxel coordinates also have to be flipped (see down)
mri = flipdim(mri, 1);
flip = [-1 0 0 256
         0 1 0 0
         0 0 1 0
         0 0 0 1    ];
transformMatrix = flip*transformMatrix;

% re-compute the homogeneous transformation matrices (apply voxel scaling)
scale = eye(4);
scale(1,1) = hdr.mmPerPixel_sagittal;
scale(2,2) = hdr.mmPerPixel_coronal;
scale(3,3) = hdr.mmPerPixel_axial;
hdr.transformHead2MRI = transformMatrix*inv(scale);
hdr.transformMRI2Head = scale*inv(transformMatrix);

% determine location of fiducials in MRI voxel coordinates
% flip the fiducials in voxel coordinates to correspond to the previous flip along left-right
hdr.fiducial.mri.nas = [256 - hdr.HeadModel.Nasion_Sag hdr.HeadModel.Nasion_Cor hdr.HeadModel.Nasion_Axi];
hdr.fiducial.mri.lpa = [256 - hdr.HeadModel.LeftEar_Sag hdr.HeadModel.LeftEar_Cor hdr.HeadModel.LeftEar_Axi];
hdr.fiducial.mri.rpa = [256 - hdr.HeadModel.RightEar_Sag hdr.HeadModel.RightEar_Cor hdr.HeadModel.RightEar_Axi];

% compute location of fiducials in MRI and HEAD coordinates

% hdr.fiducial.head.nas = warp3d(hdr.transformMRI2Head, hdr.fiducial.mri.nas, 'homogenous');
% hdr.fiducial.head.lpa = warp3d(hdr.transformMRI2Head, hdr.fiducial.mri.lpa, 'homogenous');
% hdr.fiducial.head.rpa = warp3d(hdr.transformMRI2Head, hdr.fiducial.mri.rpa, 'homogenous');

author: PatternRecognition

repo: OpenBMI


text:
% COT    Cotangent of argument in radians.
%    COT(X) is the cotangent of the elements of X.
%
%    Class support for input X:
%       float: double, single
%
%    See also ACOT, COTD.
%
%    Reference page in Doc Center
%       doc cot
%
%    Other functions named cot
%
%       codistributed/cot    gpuArray/cot    sym/cot    ts/cot
%

author: jmaih


text:
function mom = cpot_to_mpot(can)
% CPOT_TO_MPOT Convert a canonical potential to moment form.
% mom = cpot_to_mpot(can)

[logp, mu, Sigma] = canonical_to_moment(can.g, can.h, can.K);
mom = mpot(can.domain, can.sizes, logp, mu, Sigma);

%%%%%%%

function [logp, mu, Sigma] = canonical_to_moment(g, h, K)
% CANONICAL_TO_MOMENT Convert canonical characteristics to moment form.
% [logp, mu, Sigma] = canonical_to_moment(g, h, K)

if det(K)==0
  Sigma = inf*size(K);
else
  Sigma = inv(K);
end
mu = Sigma*h;
n = length(mu);
if isempty(mu)
  logp = g - 0.5*(log(det(K)) - n*log(2*pi));
else
  logp = g - 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [err, out_fn] = Norm (InBrik, N_runs, run, run_norm, run_mask, foutid_log)
%
%   [err, norm_fn] = Norm(InBrik, N_runs, run)
%
%Purpose:
%
%   Normalize the original dataset so that the beta values in 3dDeconvolve would be automatically percent signal change.
%
%Input Parameters:
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Gang Chen
%     Date : Mon Oct 20 11:26:38 EDT 2003
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'Norm';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

[InBrikPrefix, InBrikView, InBrikExt] =  AfniPrefix (InBrik);
rm_cmmd = 'rm -f ';

%Originally 3dClipLevel was used to generate a cliplevel for each run,
%then a step function chops off those voxel outside of the brain.
%But due to the concern of losing the sensitive amygdala and orbitofrontal cortex,
%it might be better to run 3dAutomask, which is roughly equivalent to
%3dClipLevel but would protect amygdala and orbitofrontal cortex regions
%from being masked.

%create command for executing 3dAutomask
if (run_mask.do == 1),
   automask_fn = sprintf('%s_msk', InBrikPrefix);
   automask = sprintf('3dAutomask -dilate %i -prefix %s %s%s', run_mask.N_dilate, automask_fn, InBrikPrefix, InBrikView);
   fprintf(2, '\nRunning: %s\n', automask);
   fprintf(foutid_log, '%s\n', automask);
   [s,w] = system(automask);
   if (s),
      fprintf(2,'Error during 3dAutomask: \n\n%s\n\n',w);
      while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
   end
	rm_cmmd = sprintf('%s %s%s%s', rm_cmmd, automask_fn, InBrikView, InBrikExt);

	%If no normalization, do masking
   if (run_norm == 0),
	   out_fn = sprintf('%s_mask', InBrikPrefix);
      mask_signal = sprintf('3dcalc -fscale -a  %s%s  \\\n', InBrikPrefix, InBrikView);
 	   mask_signal = sprintf('%s    -b %s%s \\\n', mask_signal, automask_fn, InBrikView);
	   mask_signal = sprintf('%s    -expr  "a*b" \\\n', mask_signal);
      mask_signal = sprintf('%s    -prefix %s', mask_signal, out_fn);
	   fprintf(2, '\nRunning: %s\n', mask_signal);
	   fprintf(foutid_log, '%s\n', mask_signal);
	   [s,w] = system(mask_signal);
      if (s),
         fprintf(2,'Error during 3dcalc: \n\n%s\n\n',w);
         while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
      end
   end
	rm_cmmd = sprintf('%s %s%s%s', rm_cmmd, automask_fn, InBrikView, InBrikExt);
end

if (run_norm == 1),   %for normalization
   out_fn = sprintf('%s_norm', InBrikPrefix);
   concat = sprintf('3dTcat -prefix %s ', out_fn);

   for (i=1:1:N_runs),

   %store each mean dataset in a temporary file, and delete it before leaving this function.
      mean_fn = sprintf('%s_mean%i', InBrikPrefix, i);
	   get_mean = sprintf('3dTstat -mean -prefix %s %s%s''[%i..%i]''', mean_fn, InBrikPrefix, InBrikView, run(i).first, run(i).last);
      fprintf(2, '\nRunning: %s\n', get_mean);
	   fprintf(foutid_log, '%s\n', get_mean);
	   [s, w] = system(get_mean);
	   if (s),
         fprintf(2,'Error at running 3dTstat: \n\n\n%s\n', w);
         while (1); fprintf(2,'Halted: Ctrl+c to exit');pause; end
      end

	%normalizing
%	norm_signal = sprintf('3dcalc -datum float -a  %s%s''[%i..%i]''  \\\n', InBrikPrefix, InBrikView, run(i).first, run(i).last);
	   norm_signal = sprintf('3dcalc -fscale -a  %s%s''[%i..%i]''  \\\n', InBrikPrefix, InBrikView, run(i).first, run(i).last);
   	norm_signal = sprintf('%s    -b %s%s \\\n', norm_signal, mean_fn, InBrikView);
		if (run_mask.do == 1),
		   norm_signal = sprintf('%s    -c %s%s \\\n', norm_signal, automask_fn,InBrikView);
		%	norm_signal = sprintf('%s    -expr  "(a/b*100)*step(b-clip_level)" \\\n', norm_signal);
   	   norm_signal = sprintf('%s    -expr  "(a/b*100)*c" \\\n', norm_signal);
		else norm_signal = sprintf('%s    -expr  "a/b*100" \\\n', norm_signal);
		end
	   norm_out_fn = sprintf('%s_norm%i', InBrikPrefix, i);
   	norm_signal = sprintf('%s    -prefix %s', norm_signal, norm_out_fn);
		fprintf(2, '\nRunning: %s\n', norm_signal);
	   fprintf(foutid_log, '%s\n', norm_signal);
	   [s,w] = system(norm_signal);
      if (s),
         fprintf(2,'Error during 3dcalc: \n\n%s\n\n',w);
         while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
      end

	%Concatenate the above files together
	   concat = sprintf('%s %s%s%s', concat, norm_out_fn, InBrikView, InBrikExt);

	%create remove command for deleting those individual percent signal change and mean file for each run.
	   rm_cmmd = sprintf('%s %s%s%s %s%s%s', rm_cmmd, norm_out_fn, InBrikView, InBrikExt, mean_fn, InBrikView, InBrikExt);

   end

%Concatenate all runs back together
   fprintf(2, '\nRunning: %s\n', concat);
   fprintf(foutid_log, '%s\n', concat);
   [s,w] = system(concat);
   if (s),
      fprintf(2,'Error during running 3dTcat: \n\n%s\n\n',w);
      while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
   end
end

%Don't forget to remove those temporary files above.
fprintf(2, '\nRemoving: %s\n', rm_cmmd);
fprintf(foutid_log, '%s\n', rm_cmmd);
[s,w] = system(rm_cmmd);
if (s),
    fprintf(2,'Error during rm: \n\n%s\n\n',w);
    while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dens = dmvt_ls (y, mu, Sigma, lambda, nu)
    % This function computes Density/CDF of Skew-T with scale location.
    % y must be a matrix where each row has
    %    a multivariate vector of dimension
    %    ncol(y) = p , nrow(y) = sample size.
    % mu, lambda: must be of the vector type of
    %             the same dimension equal to ncol(y) = p
    % lambda: 1 x p
    % Sigma: Matrix p x p
    nrow = @(x) size(x,1);
    ncol = @(x) size(x,2);
    n = nrow(y);
    p = ncol(y);
    mahalanobis_d = (pdist2(y,mu,'mahalanobis',Sigma)).^2;
    denst = (gamma((p + nu)./2)./(gamma(nu./2) .* pi.^(p./2))) .* ...
            nu.^(-p./2) .* det(Sigma).^(-1/2) .* ...
            (1 + mahalanobis_d./nu).^(-(p + nu)./2);
    dens = 2 .* (denst) .* tcdf(sqrt((p + nu)./(mahalanobis_d + nu)) .* sum(repmat(lambda / matrix_sqrt(Sigma), n, 1) .* (y - mu), 2), nu + p);
end

author: ramintoosi

repo: ROSS

save_path: github-repos/MATLAB/ramintoosi-ROSS


text:
function p = runICPPairList(p, g, VH)

msg('S', {g.procICP{:} 'FIND OVERLAPPING POINT CLOUDS'}, 'LogLevel', 'basic');

% List of all possible point cloud pairList
p.PairList = nchoosek(1:g.nPC, 2);

idx2del = [];

% Check intersection of voxel hulls for each pair
for i = 1:size(p.PairList,1)

    % Indices of point clouds of actual pair
    idxPC1 = p.PairList(i,1);
    idxPC2 = p.PairList(i,2);

    % Intersection of voxel hulls
    C = intersect(VH{idxPC1}(:, 1:3), VH{idxPC2}(:, 1:3), 'rows');

    % Number of intersecting voxels
    noIntersectingVoxel(i) = size(C,1);

    % Delete pair?
    if (noIntersectingVoxel(i) < p.MinNoIntersectingVoxel) || ... % to less intersecting voxels?
       all(ismember([idxPC1 idxPC2], p.IdxFixedPointClouds))   % both point clouds are fixed
        idx2del = [idx2del; i];
    end

end

% Delete pairs
p.PairList(idx2del, :) = [];
noIntersectingVoxel(idx2del) = []; % for report

% Switch columns so that fixed point clouds are always in first column (i.e. query points are selected in fixed point clouds)
for i = 1:numel(p.IdxFixedPointClouds)
    idx2switch = p.PairList(:,2) == p.IdxFixedPointClouds(i);
    p.PairList(idx2switch,:) = [p.PairList(idx2switch,2) p.PairList(idx2switch,1)]; % switch columns
end

% Report
if isempty(p.PairList)
    error('Input point clouds are not overlapping!');
else
    msg('T', '-------------------------------------------------------------------------------', 'LogLevel', 'basic');
    msg('T', 'LIST OF POINT CLOUD PAIRS:', 'LogLevel', 'basic');
    msg('T', sprintf('%12s %14s %14s %14s', 'pair', 'point_cloud', 'point_cloud', 'overlapping'), 'LogLevel', 'basic');
    msg('T', sprintf('%12s %14s %14s %14s', ''    , ''           , ''           , 'voxel'      ), 'LogLevel', 'basic');
    for i = 1:size(p.PairList,1)
        msg('T', sprintf('%12s %14s %14s %14d', ['[' num2str(i) ']'], ['[' num2str(p.PairList(i,1)) ']'], ['[' num2str(p.PairList(i,2)) ']'], noIntersectingVoxel(i)), 'LogLevel', 'basic');
    end
    msg('T', '-------------------------------------------------------------------------------', 'LogLevel', 'basic');
end

% Warning if point cloud is not overlapping with other point clouds
for i = 1:g.nPC
    if sum(p.PairList(:) == i) == 0
        msg('X', {g.procICP{:} 'FIND OVERLAPPING POINT CLOUDS'}, sprintf('attention, point cloud [%d] is not overlapping with other point clouds!', i), 'LogLevel', 'basic');
    end
end
msg('E', {g.procICP{:} 'FIND OVERLAPPING POINT CLOUDS'}, 'LogLevel', 'basic');

end

author: pglira


text:
function [ i, j ] = ij_next_gt ( i, j, n )

%*****************************************************************************80
%
%% IJ_NEXT_GT returns the next matrix index, with the constraint that I < J.
%
%  Discussion:
%
%    For N = 5, the sequence of indices returned is:
%
%      (1,2), (1,3), (1,4), (1,5), (2,3), (2,4), (2,5), (3,4), (3,5), (4,5).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the current pair of indices.
%
%    Input, integer N, the maximum value for I and J.
%    A value of N less than 2 is nonsense.
%
%    Output, integer I, J, the next pair of indices.  If either index
%    is illegal on input, the output value of (I,J) will be (1,2).
%
  if ( n < 2 )
    i = 0;
    j = 0;
    return
  end

  if ( i < 1 | n < i | j < 1 | n < j | j <= i )
    i = 1;
    j = 2;
    return
  end

  if ( j < n )
    j = j + 1;
  elseif ( i < n - 1 )
    i = i + 1;
    j = i + 1;
  else
    i = 0;
    j = 0;
  end

  return
end


text:
function [dat,marker, hdr, HDR]=Load_Emotiv(filename)

[Temp_hdr Temp_data]=edfread(filename);

[hdr, dat ,marker]=Convert_Emotiv(Temp_hdr,Temp_data);

% hdr.Label=Temp_hdr.label';
%
% hdr.NS=length(hdr.Label);
%
% [HDR] = leadidcodexyz(hdr)


text:
function boundary = p12_boundary_nearest ( m, n, point )

%*****************************************************************************80
%
%% P12_BOUNDARY_NEAREST returns a nearest boundary point in problem 12.
%
%  Discussion:
%
%    Actually, we only need "reasonably accurate" values when the point
%    to be tested is outside the region, and not to far from it.  Since
%    the politics of interior propinquity are surprisingly complicated for
%    this region, we will plump for a simple, approximate scheme.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the coordinates of the points.
%
%    Output, real BOUNDARY(M,N), points on the boundary
%    that are nearest to each point.
%
  x1 = 0.000;
  x2 = 0.500;
  x2p5 = 0.5625;
  x3 = 0.625;
  x3p5 = 0.8125;
  x4 = 1.000;
  y1 = 0.000;
  y2 = 0.250;
  y2p5 = 0.31250;
  y3 = 0.375;
  y4 = 1.000;

  t1(1:2,1:3) = [ ...
    0.0000, 0.0000; ...
    0.5000, 0.0000; ...
    0.2500, 0.2500 ]';

  q1(1:2,1:4) = [ ...
    0.0000, 0.0000; ...
    0.2500, 0.2500; ...
    0.2500, 0.7500; ...
    0.0000, 1.0000 ]';

  t2(1:2,1:3) = [ ...
    0.5000, 0.0000; ...
    0.5000, 0.2500; ...
    0.2500, 0.2500 ]';

  q2(1:2,1:4) = [ ...
    0.2500, 0.2500; ...
    0.5000, 0.2500; ...
    0.5000, 0.3125; ...
    0.2500, 0.3125 ]';

  q3(1:2,1:4) = [ ...
    0.2500, 0.3125; ...
    0.5000, 0.2125; ...
    0.5000, 0.3750; ...
    0.2500, 0.375 ]';

  q4(1:2,1:4) = [ ...
    0.5000, 0.3750; ...
    0.5000, 1.0000; ...
    0.2500, 0.7500; ...
    0.2500, 0.375 ]';

  t3(1:2,1:3) = [ ...
    0.5000, 1.0000; ...
    0.0000, 1.0000; ...
    0.2500, 0.7500 ]';

  q5(1:2,1:4) = [ ...
    0.5000, 0.0000; ...
    0.5625, 0.0000; ...
    0.5625, 0.0625; ...
    0.5000, 0.1250 ]';

  q6(1:2,1:4) = [ ...
    0.6250, 0.0000; ...
    0.6250, 0.1250; ...
    0.6625, 0.0625; ...
    0.6625, 0.0000 ]';

  t4(1:2,1:3) = [ ...
    0.5000, 0.1250; ...
    0.5625, 0.0625; ...
    0.6250, 0.1250 ]';

  t5(1:2,1:3) = [ ...
    0.5000, 0.3750; ...
    0.6250, 0.3750; ...
    0.5625, 0.4375 ]';

  q7(1:2,1:4) = [ ...
    0.5000, 0.3750; ...
    0.5625, 0.4375; ...
    0.5625, 1.0000; ...
    0.5000, 1.0000 ]';

  q8(1:2,1:4) = [ ...
    0.6250, 0.3750; ...
    0.6250, 1.0000; ...
    0.5626, 1.0000; ...
    0.5625, 0.4375 ]';

  t6(1:2,1:3) = [ ...
    0.6250, 0.0000; ...
    1.0000, 0.0000; ...
    0.8125, 0.1875 ]';

  q9(1:2,1:4) = [ ...
    0.6250, 0.0000; ...
    0.8125, 0.1875; ...
    0.8125, 0.2500; ...
    0.6250, 0.2500 ]';

  q10(1:2,1:4) = [ ...
    0.6250, 0.2500; ...
    0.8125, 0.2500; ...
    0.8125, 0.3125; ...
    0.6250, 0.3125 ]';

  q11(1:2,1:4) = [ ...
    0.6250, 0.3125; ...
    0.8125, 0.3125; ...
    0.8125, 0.3750; ...
    0.6250, 0.3750 ]';

  q12(1:2,1:4) = [ ...
    0.6250, 0.3750; ...
    0.8125, 0.3750; ...
    0.8125, 0.8125; ...
    0.6250, 1.0000 ]';

  q13(1:2,1:4) = [ ...
    1.0000, 0.0000; ...
    1.0000, 1.0000; ...
    0.8125, 0.8125; ...
    0.8125, 0.1875 ]';

  t7(1:2,1:3) = [ ...
    0.6250, 1.0000; ...
    0.8125, 0.8125; ...
    1.0000, 1.0000 ]';

  for j = 1 : n
%
%  Column 1
%
    if ( point(1,j) <= x1 & point(2,j) <= y1 )

      boundary(1,j) = x1;
      boundary(2,j) = y1;

    elseif ( point(1,j) <= x1 & point(2,j) <= y4 )

      boundary(1,j) = x1;
      boundary(2,j) = point(2,j);

    elseif ( point(1,j) <= x1 & y4 <= point(2,j) )

      boundary(1,j) = x1;
      boundary(2,j) = y4;
%
%  Column 2
%
    elseif ( point(1,j) <= x2 & point(2,j) <= y1 )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y1;

    elseif ( triangle_contains_point_2d ( t1, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y1;

    elseif ( quad_contains_point_2d ( q1, point(1:2,j) ) )

      boundary(1,j) = x1;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( t2, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( q2, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = y2;

    elseif ( quad_contains_point_2d ( q3, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = y3;

    elseif ( quad_contains_point_2d ( q4, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = point(2,j);

    elseif ( triangle_contains_point_2d ( t3, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y4;

    elseif ( point(1,j) <= x2 & y4 <= point(2,j) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y4;
%
%  Under the H crossbar.
%
    elseif ( point(1,j) <= x2p5 & point(2,j) <= y1 )

      boundary(1,j) = x2;
      boundary(2,j) = y1;

    elseif ( point(1,j) <= x3 & point(2,j) <= y1 )

      boundary(1,j) = x3;
      boundary(2,j) = y1;

    elseif ( quad_contains_point_2d ( q5, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( q6, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = point(2,j);

    elseif ( triangle_contains_point_2d ( t4, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y2;
%
%  In the H crossbar.
%
    elseif ( point(1,j) <= x3 & y2 <= point(2,j) & point(2,j) <= y2p5 )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y2;

    elseif ( point(1,j) <= x3 & y2p5 <= point(2,j) & point(2,j) <= y3 )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y3;
%
%  Above the H crossbar.
%
    elseif ( triangle_contains_point_2d ( t5, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y3;

    elseif ( quad_contains_point_2d ( q7, point(1:2,j) ) )

      boundary(1,j) = x2;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( q8, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = point(2,j);

    elseif ( point(1,j) <= x2p5 & y4 <= point(2,j) )

      boundary(1,j) = x2;
      boundary(2,j) = y4;

    elseif ( point(1,j) <= x3 & y4 <= point(2,j) )

      boundary(1,j) = x3;
      boundary(2,j) = y4;
%
%  Column 4
%
    elseif ( point(1,j) <= x4 & point(2,j) <= y1 )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y1;

    elseif ( triangle_contains_point_2d ( t6, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y1;

    elseif ( quad_contains_point_2d ( q9, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( q10, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = y2;

    elseif ( quad_contains_point_2d ( q11, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = y3;

    elseif ( quad_contains_point_2d ( q12, point(1:2,j) ) )

      boundary(1,j) = x3;
      boundary(2,j) = point(2,j);

    elseif ( quad_contains_point_2d ( q13, point(1:2,j) ) )

      boundary(1,j) = x4;
      boundary(2,j) = point(2,j);

    elseif ( triangle_contains_point_2d ( t7, point(1:2,j) ) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y4;

    elseif ( point(1,j) <= x4 & y4 <= point(2,j) )

      boundary(1,j) = point(1,j);
      boundary(2,j) = y4;
%
%  Column 5
%
    elseif ( x4 <= point(1,j) & point(2,j) <= y1 )

      boundary(1,j) = x4;
      boundary(2,j) = y1;

    elseif ( x4 <= point(1,j) & point(2,j) <= y4 )

      boundary(1,j) = x4;
      boundary(2,j) = point(2,j);

    elseif ( x4 <= point(1,j) & y4 <= point(2,j) )

      boundary(1,j) = x4;
      boundary(2,j) = y4;
    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'P12_BOUNDARY_NEAREST - Fatal error!\n' );
      fprintf ( 1, '  Orphan point = ( %f, %f )\n', point(1:2,j) );
      error ( 'P12_BOUNDARY_NEAREST - Fatal error!' );

    end

  end

  return
end


text:
function area = polygonArea3d(poly, varargin)
%POLYGONAREA3D Area of a 3D polygon.
%
%   AREA = polygonArea3d(POLY)
%   POLY is given as a N-by-3 array of vertex coordinates. The resulting
%   area is positive.
%   Works also for polygons given as a cell array of polygons.
%
%   Example
%     % area of a simple 3D square
%     poly = [10 30 20;20 30 20;20 40 20;10 40 20];
%     polygonArea3d(poly)
%     ans =
%        100
%
%     % Area of a 3D mesh
%     [v f] = createCubeOctahedron;
%     polygons = meshFacePolygons(v, f);
%     areas = polygonArea3d(polygons);
%     sum(areas)
%     ans =
%         18.9282
%
%   See also
%     polygons3d, triangleArea3d, polygonArea, polygonCentroid3d

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2012-02-24, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

% Check multiple polygons
if iscell(poly) || sum(sum(isnan(poly))) > 0
    % split the polygons into a cell array
    polygons = splitPolygons3d(poly);
    nPolys = length(polygons);

    % compute area of each polygon
    area = zeros(nPolys, 1);
    for  i = 1:nPolys
        area(i) = polygonArea3d(polygons{i});
    end

    return;
end

% put the first vertex at origin (reducing computation errors for polygons
% far from origin)
v0 = poly(1, :);
poly = bsxfun(@minus, poly, v0);

% indices of next vertices
N = size(poly, 1);
iNext = [2:N 1];

% compute cross-product of each elementary triangle formed by origin and
% two consecutive vertices
cp = cross(poly, poly(iNext,:), 2);

% choose one of the triangles as reference for the normal direction
vn = vectorNorm3d(cp);
[tmp, ind] = max(vn); %#ok<ASGLU>
cpRef = cp(ind,:);

% compute the sign of the area of each triangle
% (need to compute the sign explicitely, as the norm of the cross product
% does not keep orientation within supporting plane)
sign_i = sign(dot(cp, repmat(cpRef, N, 1), 2));

% compute area of each triangle, using sign correction
area_i = vectorNorm3d(cp) .* sign_i;

% sum up individual triangles area
area = sum(area_i) / 2;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ orig ] = plx_orig_header( fname )

% PLX_ORIG_HEADER Extracts the header informations of plx files using the
% Plexon Offline SDK, which is available from
% http://www.plexon.com/assets/downloads/sdk/ReadingPLXandDDTfilesinMatlab-mexw.zip
%
% Use as
%   [orig] = plx_orig_header(filename)
%
% Copyright (C) 2012 by Thomas Hartmann
%
% This code can be redistributed under the terms of the GPL version 3 or
% newer.

% get counts...
[tscounts, wfcounts, evcounts, contcounts] = plx_info(fname, 0);

orig.TSCounts = tscounts;
orig.WFCounts = wfcounts;

% get event channels...
[dummy, evchans] = plx_event_chanmap(fname);

orig.EVCounts = zeros(512, 1);

for i=1:length(evchans)
  try
    orig.EVCounts(evchans(i)+1) = evcounts(i);
  catch
  end
end %for

% get more infos...
[OpenedFileName, Version, Freq, Comment, Trodalness, NPW, PreTresh, SpikePeakV, SpikeADResBits, SlowPeakV, SlowADResBits, Duration, DateTime] = plx_information(fname);
orig.MagicNumber = 1480936528;
orig.Version = Version;
orig.ADFrequency = Freq;
orig.Comment = Comment;
orig.NumEventChannels = length(evchans);
orig.NumSlowChannels = length(contcounts);
[orig.NumDSPChannels, dummy] = plx_chanmap(fname);
orig.NumPointsWave = NPW;
orig.NumPointsPreThr = PreTresh;
[orig.Year, orig.Month, orig.Day, orig.Hour, orig.Minute, orig.Second] = datevec(DateTime);
orig.LastTimestamp = Duration * Freq;
orig.Trodalness = Trodalness;
orig.DataTrodalness = Trodalness;
orig.BitsPerSpikeSample = SpikeADResBits;
orig.BitsPerSlowSample = SlowADResBits;
orig.SpikeMaxMagnitudeMV = SpikePeakV;
orig.SlowMaxMagnitudeMV = SlowPeakV;

% gather further info for additional headers...
[dummy, chan_names] = plx_chan_names(fname);
[dummy, dspchans] = plx_chanmap(fname);
[dummy, gains] = plx_chan_gains(fname);
[dummy, filters] = plx_chan_filters(fname);
[dummy, thresholds] = plx_chan_thresholds(fname);
[dummy, evnames] = plx_event_names(fname);
[dummy, ad_names] = plx_adchan_names(fname);
[dummy, adchans] = plx_ad_chanmap(fname);
[dummy, ad_freqs] = plx_adchan_freqs(fname);
[dummy, ad_gains] = plx_adchan_gains(fname);

% do channelheaders...
for i=1:orig.NumDSPChannels
  head = [];
  head.Name = ['mua_' chan_names(i, :)];
  head.SIGName = ['mua_' chan_names(i, :)];
  head.Channel = dspchans(i);
  head.SIG = i;
  head.Gain = gains(i);
  head.Filter = filters(i);
  head.Threshold = thresholds(i);

  orig.ChannelHeader(i) = head;
end %for

% do eventheaders...
for i=1:orig.NumEventChannels
  head = [];
  head.Name = evnames(i, :);
  head.Channel = evchans(i);

  orig.EventHeader(i) = head;
end %for

% do slowchannelheaders...
for i=1:orig.NumSlowChannels
  head = [];
  head.Name = ['lfp_' ad_names(i, :)];
  head.Channel = adchans(i);
  head.ADFreq = ad_freqs(i);
  head.Gain = ad_gains(i);

  orig.SlowChannelHeader(i) = head;
end %for

end


text:
%% Fit GLG and GFM models to image
% The GLG model and the GFM model are both fitted to the wavelet
% coefficients of the 'lena' image.
% The fitted densities are compared to the observed distribution of
% coefficients

%% Load the image and compute the wavelet transform

X = imread('lena.png');
I = im2double( X );

dwtmode('per', 'nodisp');
[C, S] = wavedec2(I, 3, 'db4');
tree = dwt2_to_cell( C, S );

%% Fit GLG model

M = compute_moments2D( tree );
theta_glg = moment_estimation( M );
FitInfo_glg = GLG_EM_wrapper( theta_glg, tree );

%% Fit GFM model

theta_gfm = GFM_init( tree, 2, 0.01, 0.001 );
FitInfo_gfm = GFM_EM_wrapper( theta_gfm, tree, 'max_itr', 50 );

%% Plot distributions

% Level and direction of transform
l = 2;
d = 2;

y = tree{l+1}{d}(:);
x = min(y):5e-3:max(y);
n = length(y);
[heights, locations] = hist(y, x);

% Normalized histogram
width = locations(2) - locations(1);
heights = heights / (n*width);
bar( locations, heights, 'hist' )

% Density from GLG model
w = linspace(min(y), max(y), 200);
glg = GLG_marginal_density( w, FitInfo_glg.model(l,1,d), FitInfo_glg.model(l,2,d) );
glg_h = line(w, glg, 'color', 'r');

% Density from GFM model
gfm = GFM_marginal_density( w, FitInfo_gfm.model{l,1,d}, FitInfo_gfm.model{l,3,d} );
gfm_h = line(w, gfm, 'color', 'g');

% Format figure
xlim( 0.5*[-1 1] );
legend( [glg_h, gfm_h], 'GLG density', 'GFM density' );

author: Sable

repo: mcbench-benchmarks


text:
function [relative_yaw_angle_rad, bankingangle_rad, index, info] = TrackPosition(ref_banking, States)
%__________________________________________________________________________
%% Documentation
%
% Author:       Dominik St�rk (dominik.staerk@tum.de)
%
% Start Date:   17.09.2020
%
% Description:
% In reference to the current track position and yaw angle this function
% interpolates the coresponding banking angle and heading offset. The
% data is provided by the 'ref_banking' entry in the track�s data
% dictionary.
%
% Inputs:
%   ref_banking
%   States
%
% Outputs:
%   rel_yaw_agngle  - yaw angle offset to raceline yaw angle dataset
%   bankingangle_rad     - banking angle
%   index           - current index in raceline dataset
%   info            - more info of the current track position

%% init
x_m             = States(1);
y_m             = States(2);
psi_rad         = States(3);

x_ref_m         = ref_banking.x_ref_m;
y_ref_m         = ref_banking.y_ref_m;

%% get banking angle and relative yaw angle
% get the index of the current track position t2 out of ref_banking in
% track.sldd
[vdp, xy_loc_m, index] = getIndex(x_m, y_m, x_ref_m, y_ref_m);

% get raceline vectors to interpolate
[Scenario32, Scenario21, length32_m, length21_m, progress_m, offset_m, sgn] = ...
                getRacelineVectors(xy_loc_m, vdp, index, x_ref_m, y_ref_m);

% interpolate the actual banking angle, relative yaw angle offset and track
% width to the left and right corresponding to the active Scenario
if Scenario32 == true
   DX = [0, length32_m];
   if index == vdp
   DY_eps = [ref_banking.bankingangle_rad(vdp),...
            ref_banking.bankingangle_rad(1)];
   DY_psi = [ref_banking.psi_racetraj_rad(vdp),...
            ref_banking.psi_racetraj_rad(1)];
   DY_w_l = [ref_banking.width_left_m(vdp),...
            ref_banking.width_left_m(1)];
   DY_w_r = [ref_banking.width_right_m(vdp),...
            ref_banking.width_right_m(1)];
   else
   DY_eps = [ref_banking.bankingangle_rad(index),...
            ref_banking.bankingangle_rad(index+1)];
   DY_psi = [ref_banking.psi_racetraj_rad(index),...
            ref_banking.psi_racetraj_rad(index+1)];
   DY_w_l = [ref_banking.width_left_m(index),...
            ref_banking.width_left_m(index+1)];
   DY_w_r = [ref_banking.width_right_m(index),...
            ref_banking.width_right_m(index+1)];
   end
   bankingangle_rad      = InterpLin(DX,DY_eps,progress_m);
   psi_racetraj_rad = InterpLin(DX,DY_psi,progress_m);
   track_w_l_m      = InterpLin(DX,DY_w_l,progress_m);
   track_w_r_m      = InterpLin(DX,DY_w_r,progress_m);

elseif Scenario21 == true
   DX = [0, length21_m];
   if index == 1
   DY_eps = [ref_banking.bankingangle_rad(vdp),...
            ref_banking.bankingangle_rad(1)];
   DY_psi = [ref_banking.psi_racetraj_rad(vdp),...
            ref_banking.psi_racetraj_rad(1)];
   DY_w_l = [ref_banking.width_left_m(vdp),...
            ref_banking.width_left_m(1)];
   DY_w_r = [ref_banking.width_right_m(vdp),...
            ref_banking.width_right_m(1)];
   else
   DY_eps = [ref_banking.bankingangle_rad(index-1),...
            ref_banking.bankingangle_rad(index)];
   DY_psi = [ref_banking.psi_racetraj_rad(index-1),...
            ref_banking.psi_racetraj_rad(index)];
   DY_w_l = [ref_banking.width_left_m(index-1),...
            ref_banking.width_left_m(index)];
   DY_w_r = [ref_banking.width_right_m(index-1),...
            ref_banking.width_right_m(index)];
   end
   bankingangle_rad      = InterpLin(DX,DY_eps,progress_m);
   psi_racetraj_rad = InterpLin(DX,DY_psi,progress_m);
   track_w_l_m      = InterpLin(DX,DY_w_l,progress_m);
   track_w_r_m      = InterpLin(DX,DY_w_r,progress_m);
else
   bankingangle_rad      = ref_banking.bankingangle_rad(index);
   psi_racetraj_rad = ref_banking.psi_racetraj_rad(index);
   track_w_l_m      = ref_banking.width_left_m(index);
   track_w_r_m      = ref_banking.width_right_m(index);
end

% return the relative yaw angle offset
rot_psi_rad  = [cos(psi_rad),-sin(psi_rad);sin(psi_rad),...
                  cos(psi_rad)] * [0;1];
ref             = [0,-1;1,0]*rot_psi_rad;
rot             = [cos(psi_racetraj_rad),-sin(psi_racetraj_rad);...
                  sin(psi_racetraj_rad),cos(psi_racetraj_rad)];

% fix to avoid error due to numerical inaccuracy
beta = acos(calc_acosInput(rot*[0;1],ref));

relative_yaw_angle_rad = beta - pi/2;

% return data to calculate the track profile, if ndtm is selected
info = [sgn, offset_m, track_w_l_m, track_w_r_m];

author: TUMFTM


text:
%IMREADMULTI  Loads a multi-page image from a file
%
%     imgs = cv.imreadmulti(filename)
%     imgs = cv.imreadmulti(filename, 'OptionName',optionValue, ...)
%
% ## Input
% * __filename__ Name of file to be loaded.
%
% ## Output
% * __imgs__ A cell-array of images holding each page, if more than one. If
%   If the image cannot be read (because of IO errors, improper permissions,
%   unsupported or invalid format), the function throws an error.
%
% ## Options
% * __Unchanged__ If set, return the loaded image as is (with alpha channel,
%   otherwise it gets cropped). Both the depth and number of channels are
%   unchanged as determined by the decoder. default false
% * __AnyDepth__ If set, return 16-bit/32-bit image when the input has the
%   corresponding depth, otherwise convert it to 8-bit. default false
% * __AnyColor__ If set, the image is read in any possible color format.
%   default true
% * __Color__ If set, always convert image to the 3 channel BGR color image.
%   default false
% * __Grayscale__ If set, always convert image to the single channel grayscale
%   image. default false
% * __GDAL__ If set, use the gdal driver for loading the image. default false
% * __IgnoreOrientation__ If set, do not rotate the image according to EXIF's
%   orientation flag. default false
% * __Flags__ Advanced option to directly set the flag specifying the depth
%   and color type of a loaded image. Note that setting this integer flag
%   overrides all the other flag options. Not set by default:
%   * `>0`: Return a 3-channel color image. Note that in the current
%           implementation the alpha channel, if any, is stripped from the
%           output image. For example, a 4-channel RGBA image is loaded as
%           RGB if `Flags >= 0`.
%   * `=0`: Return a grayscale image
%   * `<0`: Return the loaded image as is (with alpha channel if present).
% * __FlipChannels__ in case the output is color image, flips the color order
%   from OpenCV's BGR/BGRA to MATLAB's RGB/RGBA order. default true
%
% The function cv.imreadmulti loads a multi-page image from the specified file
% into a cell-array of matrices.
%
% See cv.imread for details.
%
% See also: cv.imread, imread, multibandread
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv