text:
function C= train_QDA(xTr, yTr, varargin)
%TRAIN_QDA - Quadratic Discriminant Analysis
%
%Synopsis:
%   C = train_QDA(XTR, YTR)
%   C = train_QDA(XTR, YTR, PRIOR)
%
%Arguments:
% XTR: DOUBLE [NxM] - Data matrix, with N feature dimensions, and M
%                     training points/examples.
% YTR: INT [CxM] - Class membership labels of points in X_TR. C by M matrix
%                  of training labels, with C representing the number of
%                  classes and M the number of training examples/points.
%                  Y_TR(i,j)==1 if the point j belongs to class i.
% PRIOR: DOUBLE - (default ones(nClasses, 1)/nClasses): Empirical class
%                 priors
%
%Returns:
%   C: STRUCT - Trained classifier structure, with the hyperplane given by
%               fields C.w, C.b and C.sq
%
%See also:
%   APPLY_QDA

if size(yTr,1)==1,
  nClasses= 2;
  clInd{1}= find(yTr==-1);
  clInd{2}= find(yTr==1);
  N= [length(clInd{1}) length(clInd{2})];
else
  nClasses= size(yTr,1);
  clInd= cell(nClasses,1);
  N= zeros(nClasses, 1);
  for ci= 1:nClasses,
    clInd{ci}= find(yTr(ci,:));
    N(ci)= length(clInd{ci});
  end
end

if nargin==2
  priorP= ones(nClasses,1)/nClasses;
else
  priorP= varargin{1};
end
if isequal(priorP, '*')
  priorP = N/sum(N);
end

d= size(xTr,1);
C.w= zeros(d, nClasses);
C.b= zeros(1, nClasses);
C.sq= zeros(d, d, nClasses);
for ci= 1:nClasses,
  cli= clInd{ci};
  C.w(:,ci)= mean(xTr(:,cli), 2);
  yc= xTr(:,cli) - C.w(:,ci)*ones(1,N(ci));
  Sq= yc*yc';
  C.sq(:,:,ci)= Sq / (N(ci)-1);
end
C.b= zeros(1, nClasses);
for ci= 1:nClasses,
  S= C.sq(:,:,ci);
  S= pinv(S);
  C.sq(:,:,ci)= -0.5*S;
  C.b(ci)=  -0.5*C.w(:,ci)' * S*C.w(:,ci) + ...
            0.5*log(max([det(S),realmin])) + log(priorP(ci));
  C.w(:,ci)= S*C.w(:,ci);
end
C.b=C.b';

if nClasses==2,
  sq(:,:)= C.sq(:,:,2) - C.sq(:,:,1);
  C.sq= sq;
  C.w= C.w(:,2)-C.w(:,1);
  C.b= C.b(2) - C.b(1);
end

author: bbci


text:
function [p_, sdp_, c_, sdc_, dk_, sdk_, aa_, bb_]=mypval2(var1, mati)

    % clpvla.m                            A.Allmann
    % function to calculate the parameters of the modified Omori Law
    %
    %
    %
    % this function is a modification of a program by Paul Raesenberg
    % that is based on Programs by Carl Kisslinger and Yoshi Ogata
    %
    % function finds the maximum liklihood estimates of p,c and k, the
    % parameters of the modifies Omori equation
    % it also finds the standard deviations of these parameters
    %
    % Input: Earthquake Catalog of an Cluster Sequence
    %
    % Output: p c k values of the modified Omori Law with respective
    %         standard deviations
    %         A and B values of the Gutenberg Relation based on k
    %
    % Create an input window for magnitude thresholds and
    % plot cumulative number versus time to allow input of start and end
    % time

    global backcat ttcat
    global clu1 pyy tmp1 tmp2 tmp3 tmp4 difp
    global xt cumu cumu2
    global p c dk tt pc loop nn pp nit t err1x err2x ieflag isflag
    global cstep pstep tmpcat ts tend eps1 eps2
    global sdc sdk sdp  aa bb pcheck loopcheck
    global autop tmeqtime tmvar
    %if var1 == 3
    tmvar=[];
    %input of parameters(Magnitude,Time)
    ZG.newt2=ttcat;              %function operates with single cluster
    autop=0;
    if var1==4
        autop=1;
    end
    %calculate start -end time of overall catalog
    t0b = min(ZG.newt2.Date);
    teb = max(ZG.newt2.Date);
    tdiff=days(teb-t0b);       %time difference in days

    par3=tdiff/100;

    par5=par3;
    if par5>.5
        par5=par5/5;
    end

    % calculate cumulative number versus time and bin it
    %
    n = ZG.newt2.Count;
    [cumu, xt] = hist(ZG.newt2.Date,(t0b:days(par3):teb));
    [cumu, xt] = hist(days(ZG.newt2.Date-ZG.newt2.Date(1)),(0:par5:tdiff)); %WHY BOTH?
    difp= [0 diff(cumu)];
    cumu2 = cumsum(cumu);

    % find start time of time series
    %
    nn=find(cumu==max(cumu));
    nnn=nn(1,1)-2;
    tmvar=1;           %temperal variable
    if par3>=1
        tmp3=t0b+nnn*days(par3);
    else
        tmp3=nnn*par5;
    end
    tmp3 = mati;
    tmp2=min(ttcat.Magnitude);
    tmp1=max(ttcat.Magnitude);

    tmp3=max(0,tmp3);

    tmp4=teb;

    %%end

    %cumputation part after parameter input
    %elseif var1==8  | var1==6 | var1==7
    %set the error test values
    eps1=.0005;
    eps2=.0005;

    %set the parameter starting values
    PO=1.1;
    CO=0.1;

    %set the initial step size
    pstep=.05;
    cstep=.05;
    pp=PO;
    pc=CO;
    nit=0;
    ieflag=0;
    isflag=0;
    pcheck=false;
    err1x=0;
    err2x=0;
    ts=0.0000001;
    if autop ~= 1             %input was manual

        %Build timecatalog

        mains=find(ttcat.Magnitude == max(ttcat.Magnitude),1);
        mains=ttcat.subset(mains);         %biggest shock in sequence (first one only!)
        assert(mains.Count == 1);
        if par3<0.001
            daycriteria1= ttcat.Date>=days(tmp3)+ttcat.Date(1);
            daycriteria2= ttcat.Date<=days(tmp4)+ttcat.Date(1);
            tmpcat = ttcat.subset(daycriteria1&daycriteria2);
            tmp6=days(tmp3)+ttcat.Date(1);
        else
            tmpcat=ttcat.subset(tmp3<=ttcat.Date & ttcat.Date<=tmp4);
            tmp6=tmp3;
        end
        tmpcat=tmpcat.subset(tmp2<=tmpcat.Magnitude & tmpcat.Magnitude<=tmp1);
        if var1 ==6 || var1==7
            tmpcat=tmpcat.subset(tmpcat.Date>mains.Date);
            tmpcat=[mains; tmpcat];
            ts=(tmp6-mains.Date); %# days
            ts=max(0.0000001,ts)
        end
        tmeqtime=clustime(tmpcat);
        tmeqtime=tmeqtime-tmeqtime(1);     %time in days relative to first eq
        tmeqtime=tmeqtime(2:length(tmeqtime));

        %automatic estimate works with whole sequence
    else
        tmeqtime=clustime(ttcat);
        tmeqtime=tmeqtime-tmeqtime(1);
        tmeqtime=tmeqtime(2:length(tmeqtime));

    end

    tp1 = input('tp1=   ');
    tp2 = input('tp2=    ');
    ts = tp1;
    l = tmeqtime >= tp1 & tmeqtime <= tp2;
    tmeqtime = tmeqtime(l);

    tend=tmeqtime(length(tmeqtime)); %end time

    %Loop begins here
    nn=length(tmeqtime);
    loop=0;
    loopcheck=0;
    tt=tmeqtime(nn)+1;
    t=tmeqtime;
    if pc < 0 ; pc = 0.0; end
    if pc <= ts; pc = ts + 0.05;end

    MIN_CSTEP = 0.000001;
    MIN_PSTEP = 0.00001;

    ploop_c_and_p_calcs(MIN_CSTEP, MIN_PSTEP, true,'kpc'); %call of function who calculates parameters

    if loopcheck<500
        %round values on two digits
        p=round(p, -2);
        sdp=round(sdp, -2);
        c=round(c, -3);
        sdc=round(sdc, -3);
        dk=round(dk, -2);
        sdk= round(sdk, -2);
        aa=round(aa, -2);
        bb=round(bb, -2);

        disp(['p = ' num2str(p)  ' +/- ' num2str(sdp)]);
        disp(['c = ' num2str(c)  ' +/- ' num2str(sdc)]);
        disp(['k = ' num2str(dk)  ' +/- ' num2str(sdk)]);
        disp(['b = ' num2str(bb)  ' +/- ' num2str(sdp)]);
        disp(['a = ' num2str(aa)  ' +/- ' num2str(sdp)]);
    else
        disp('No result');
        %p = nan;
        %c = nan;
        %k = nan;
        %bb = nan;
        %aa = nan;

        % p = nan; sdp = nan;
    end
    if autop~=1
        if par3>=1
            tdiff = round(tmpcat(length(tmpcat(:,1)),3)-tmpcat(1,3));
        else
            tdiff = (tmpcat(length(tmpcat(:,1)),3)-tmpcat(1,3))*365;
        end
        % set arrays to zero
        %
        if par3>=1
            cumu = 0:1:(tdiff/days(par3))+1;
            cumu2 = 0:1:(tdiff/days(par3))-1;
        else
            par5=par3/5;
            cumu = 0:par5:tdiff+2*par3;
            cumu2 =  0:par5:tdiff-1;
        end
        cumu = zeros(size(cumu));
        cumu2 = zeros(size(cumu2));
        %
    end
    tmvar=[];
    %end;

    p_=p;
    sdp_=sdp;
    c_=c;
    sdc_=sdc;
    dk_=dk;
    sdk_=sdk;
    aa_=aa;
    bb_=bb;
end

author: CelsoReyes


text:
function sliceThicknessV = deduceSliceWidths(planC,scanNum)
%function sliceThicknessV = deduceSliceWidths(planC)
%
%This function deduces slice thicknesses when the treatment planning system
%provides only z Values.
%
%
%Created:  30 Apr 03, JOD.
% copyright (c) 2001-2006, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).
indexS = planC{end};

zValuesV = [planC{indexS.scan}(scanNum).scanInfo(:).zValue];

sliceThicknessV = ones(size(zValuesV)) * nan;

for i = 2 : length(zValuesV) - 1
  nextDelta = abs(zValuesV(i+1) - zValuesV(i));
  sliceThicknessV(i) = nextDelta;
end
if length(zValuesV) > 1
    sliceThicknessV(1) = 2 * (abs(zValuesV(2) - zValuesV(1)) - 0.5 * sliceThicknessV(2));
    sliceThicknessV(end) = 2 * (abs(zValuesV(end) - zValuesV(end - 1)) - 0.5 * sliceThicknessV(end - 1));
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%This Matlab script can be used to reproduce Figure 7.13 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%% Define parameters
lambda = 1; %Wavelength (normalized)
M_H = 8; %Number of antenna per horizontal row
M_V = 4; %Number of rows
d_H = 0.5*lambda; %Horizontal antenna spacing
d_V = 0.5*lambda; %Vertical antenna spacing

%Define the antenna geometry
M = M_H*M_V; %Total number of antennas
U = zeros(3,M); %Matrix containing the position of the antennas

i = @(m) mod(m-1,M_H); %Horizontal index
j = @(m) floor((m-1)/M_H); %Vertical index

for m = 1:M
    U(:,m) = [0; i(m)*d_H; j(m)*d_V]; %Position of the mth element
end

%% Compute the spatial signature for various directions
varphi = linspace(-pi/2,pi/2,1000);
theta= [0,-pi/3];
P = zeros(length(varphi),length(theta));
v = ones(M,1);

for i = 1:length(varphi)

    for j = 1:length(theta)

        P(i,j) = abs(v'*functionSpatialSignature3DLoS(U,varphi(i),theta(j),lambda))/M;

    end

end

P = P/max(max(P));
P = 10*log10(P);

%% Plot the simulation results
figure;
hold on; box on;
plot(varphi/pi,P(:,1)','-k','LineWidth',1)
plot(varphi/pi,P(:,2)','r--','LineWidth',1)

ylim([-30,0])
ax = gca;
set(ax, 'XTick', [-0.5, -0.25, 0, 0.25, 0.5]);
xlabel('Azimuth angle $\varphi$ in multiples of $\pi$','Interpreter','latex');
ylabel('Normalized array response [dB]');
legend({'$\theta = 0$', '$\theta = -\pi/3$   '},'Interpreter','latex');

author: emilbjornson

repo: massivemimobook


text:
function w = ymdf_to_weekday_hebrew ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_WEEKDAY_HEBREW returns the weekday of a Hebrew YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, integer W, is the week day number of the date, with
%    1 for Sunday, through 7 for Saturday.
%
  jed = ymdf_to_jed_hebrew ( y, m, d, f );

  [ w, f2 ] = jed_to_weekday ( jed );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function net = make_net(opts)
% MAKE_NET constructs the network up to the output.
% Loss functions are not added.
%
% Options include:
% Architecture to use for initial branches of network.
% Whether to use scores or bounding box regression.
% Whether to use xcorr or concat layer to join.
% Architecture to put after concat (1x1 convs).

switch opts.branch.type
    case 'alexnet'
        [branch1, branch2] = make_branches_alexnet(opts);
    otherwise
        error('Unknown branch type')
end
[repr_sz, repr_stride] = output_size(branch1, [opts.exemplarSize*[1 1], 3, 1]);
switch opts.join.method
    case 'xcorr'
        [join, final] = make_join_xcorr(opts);
    case 'corrfilt'
        [join, final] = make_join_corr_filt(opts, repr_sz, repr_stride);
    otherwise
        error('Unknown join method')
end

net = make_siamese(branch1, branch2, join, final, ...
                   {'exemplar', 'instance'}, 'score', ...
                   'share_all', opts.share_params);

end

function [branch1, branch2] = make_branches_alexnet(opts)
    branch_opts.last_layer = 'conv5';
    branch_opts.num_out     = [96, 256, 384, 384, 256];
    branch_opts.num_in      = [ 3,  48, 256, 192, 192];
    branch_opts.conv_stride = [ 2,   1,   1,   1,   1];
    branch_opts.pool_stride = [ 2,   2];
    branch_opts.batchNormalization = true;
    branch_opts = vl_argparse(branch_opts, {opts.branch.conf});

    branch_opts.exemplarSize = opts.exemplarSize * [1 1];
    branch_opts.instanceSize = opts.instanceSize * [1 1];
    branch_opts.weightInitMethod = opts.init.weightInitMethod;
    branch_opts.scale            = opts.init.scale;
    branch_opts.initBias         = opts.init.initBias;

    f = @() make_branch_alexnet(branch_opts);
    branch1 = f();
    branch2 = f();
end

function [join, final] = make_join_xcorr(opts)
    join_opts.finalBatchNorm = true;
    join_opts.adjustGainInit = 1;
    join_opts.adjustBiasInit = 0;
    % Learning rates ignored if batch-norm is enabled.
    join_opts.adjustGainLR = 0;
    join_opts.adjustBiasLR = 1;
    join_opts = vl_argparse(join_opts, {opts.join.conf});

    join = dagnn.DagNN();
    join.addLayer('xcorr', XCorr(), {'in1', 'in2'}, {'out'});

    % Create adjust layer.
    final.layers = {};
    convOpts = {'CudnnWorkspaceLimit', 1024*1024*1024};
    if join_opts.finalBatchNorm
        % Batch-norm layer only.
        final.layers{end+1} = struct(...
            'type', 'bnorm', 'name', 'adjust_bn', ...
            'weights', {{single(join_opts.adjustGainInit), ...
                         single(join_opts.adjustBiasInit), ...
                         zeros(1, 2, 'single')}}, ...
            'learningRate', [2 1 0.3], ...
            'weightDecay', [0 0]);
    else
        % Linear layer only.
        final.layers{end+1} = struct(...
            'type', 'conv', 'name', 'adjust', ...
            'weights', {{single(join_opts.adjustGainInit), ...
                         single(join_opts.adjustBiasInit)}}, ...
            'learningRate', [join_opts.adjustGainLR, join_opts.adjustBiasLR], ...
            'weightDecay', [1 0], ...
            'opts', {convOpts});
    end
end

function [join, final] = make_join_corr_filt(opts, in_sz, in_stride)
    join_opts.finalBatchNorm = false;
    join_opts.const_cf = false;
    join_opts.lambda = nan;
    join_opts.window = 'not-set';
    join_opts.window_lr = 0;
    join_opts.bias = false;
    join_opts.adjust = true;
    join_opts.sigma = 0;
    join_opts.target_lr = 0;

    join_opts.adjustGainInit = 1;
    join_opts.adjustBiasInit = 0;
    % Learning rates ignored if batch-norm is enabled.
    join_opts.adjustGainLR = 0;
    join_opts.adjustBiasLR = 1;

    join_opts = vl_argparse(join_opts, {opts.join.conf});
    convOpts = {'CudnnWorkspaceLimit', 1024*1024*1024};

    join = dagnn.DagNN();
    % Apply window before correlation filter.
    join.addLayer('cf_window', MulConst(), ...
                  {'in1'}, {'cf_example'}, {'window'});
    p = join.getParamIndex('window');
    join.params(p).value = single(make_window(in_sz, join_opts.window));
    join.params(p).learningRate = join_opts.window_lr;

    % Establish whether there is a bias parameter to the XCorr.
    cf_outputs = {'tmpl'};
    xcorr_inputs = {'tmpl_cropped', 'in2'};

    % learnt alphas instead of CF for Adaptation Experiment
    if join_opts.const_cf
        join.addLayer('circ', ConvCircScalar(), ...
                      {'cf_example'}, cf_outputs, {'circf'});
        p = join.getParamIndex('circf');
        join.params(p).value = init_weight(opts.init, in_sz(1), in_sz(2), 1, 1, 'single');
    else
    % Add a correlation filter before the XCorr in branch 1.
        if join_opts.bias
            % Connect correlation filter bias to xcorr bias.
            cf_outputs = [cf_outputs, {'bias'}];
            xcorr_inputs = [xcorr_inputs, {'bias'}];
        end
        join.addLayer('cf', ...
                      CorrFilter('lambda', join_opts.lambda, ...
                                 'bias', join_opts.bias), ...
                      {'cf_example'}, cf_outputs, {'cf_target'});
        % Set correlation filter target.
        p = join.getParamIndex('cf_target');
    end

    join.addLayer('crop_z', ...
                    CropMargin('margin', 16), ...
                    cf_outputs, xcorr_inputs{1});

    % Cross-correlate template with features of other image.
    join.addLayer('xcorr', XCorr('bias', join_opts.bias), ...
                  xcorr_inputs, {'out'});

    assert(join_opts.sigma > 0);
    join.params(p).value = single(gaussian_response(in_sz, join_opts.sigma/in_stride));
    join.params(p).learningRate = join_opts.target_lr;

    % Add scalar layer to calibrate corr-filt scores for loss function.
    final.layers = {};
    if join_opts.adjust
        if  join_opts.finalBatchNorm
            % Batch-norm layer only.
            final.layers{end+1} = struct(...
                'type', 'bnorm', 'name', 'adjust_bn', ...
                'weights', {{single(join_opts.adjustGainInit), ...
                             single(join_opts.adjustBiasInit), ...
                             zeros(1, 2, 'single')}}, ...
                'learningRate', [2 1 0.3], ...
                'weightDecay', [0 0]);
        else
            final.layers{end+1} = struct(...
                'type', 'conv', 'name', 'adjust', ...
                'weights', {{single(1), single(-0.5)}}, ...
                'learningRate', [1, 2], ...
                'weightDecay', [0 0], ...
                'opts', {convOpts});
        end
    end
end

function num_out = output_dim(opts)
    % TODO: Restructure options so that this code does not need to know about
    % different types of loss function?
    switch opts.loss.type
        case {'simple', 'structured'}
            num_out = 1;
        case 'regression'
            assert(opts.instanceSize==opts.exemplarSize, 'Exemplar and Instance should have the same size.');
            assert(opts.negatives==0 && opts.hardNegatives==0, 'No negative pairs for the moment.');
            num_out = 4;
        otherwise
            error('unknown loss');
    end
end

function [out_sz, out_stride] = output_size(net, in_sz)
    % Assume that net has 1 input and 1 output.
    if isa(net, 'dagnn.DagNN')
        input = only(net.getInputs());
        output = only(net.getOutputs());
        sizes = net.getVarSizes({input, in_sz});
        out_sz = sizes{net.getVarIndex(output)}(1:3);
        rfs = net.getVarReceptiveFields(input);
        out_stride = rfs(net.getVarIndex(output)).stride;
    else
        info = vl_simplenn_display(net, 'inputSize', in_sz);
        out_sz = info.dataSize(1:3, end);
        out_stride = info.receptiveFieldStride(:, end);
    end
    out_sz = reshape(out_sz, 1, []);
    assert(all(out_stride == out_stride(1)));
    out_stride = out_stride(1);
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function [C, Ceq] = multipleShooting2(z,pack,dyn)
%
% This function enforces the dynamics of the cart-pole system
% using 4th-order runge kutta:
%
%
% INPUTS:
% 	z = [N,1] = vector of decision variables
% 	pack = struct with matrix dimensions
% 	dyn = struct with physical parameters
%
% OUTPUTS:
% 	C = []
% 	Ceq = column vector of defects. Trajectory is feasible iff norm(Ceq)->0
%
% NOTES:
%  	[t,x,u] = unPackDecVar(z,pack);
%   size(t) == [1,1];
% 	size(x) == [nState, nSegment+1]
% 	size(u) == [nControl, nSubStep*nSegment+1]
%
% * The trajectory is represented with piecewise-linear control
% 	and piece-wise quartic state (implicityly).
%

[t,x,u] = unPackDecVar(z,pack);
nx = size(x,2);
nu = size(u,2);

nSubStep = (nu-1)/(nx-1);
dt = t/(nu-1);

xLow = x(:,1:(nx-1));  %Each segment starts here
xUpp = x(:,2:nx);   %And goes here

% Time-step along each trajectory segment in parallel
xStar = xLow;
uIdx = 0:nSubStep:(nu-1);
for i=1:nSubStep
	uStarLow = u(:,i + uIdx);
	uStarUpp = u(:,i+1 + uIdx);
	uStarMid = 0.5*(uStarLow+uStarUpp);
    k1 = cartPoleDynamics(xStar,uStarLow,dyn);
    k2 = cartPoleDynamics(xStar+0.5*k1,uStarMid,dyn);
    k3 = cartPoleDynamics(xStar+0.5*k2,uStarMid,dyn);
    k4 = cartPoleDynamics(xStar+k3,uStarUpp,dyn);
    xStar = xStar + (dt/6)*(k1 + 2*k2 + 2*k3 + k4);
end

% Defect constraint:
defect = xStar-xUpp;

% Pack up and return:
C = [];
Ceq = reshape(defect,numel(defect),1);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% TEST_RING_MIXED_BC_G_NMNN: data function for Neumann boundary condition.

function g = test_ring_mixed_bc_g_nmnn (x, y, ind)

  [theta, r] = cart2pol (x,y);
  switch (ind)
    case 1
      g = -cos (theta) .* exp (x) .* (sin (x.*y) + y.*cos (x.*y)) - ...
          sin (theta) .* exp(x) .* x .* cos (x.*y);
    case 2
      g = cos (theta) .* exp (x) .* (sin (x.*y) + y.*cos (x.*y)) + ...
          sin (theta) .* exp(x) .* x .* cos (x.*y);
    case 3
      g = -x .* exp(x) .* cos (x.*y);
    case 4
      g = -exp(x) .* (sin (x.*y) + y .* cos (x.*y));
    otherwise
      error ('g_nmnn: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function pde = jumpmgdata1
%% Data of JUMPMG1
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'g_D',@g_D,'d',@omega);

    function s = f(p) % load data (right hand side function)
    s = zeros(size(p,1),1);
    end

    function s = g_D(p) % Dirichlet boundary condition
    s = zeros(size(p,1),1);
    x = p(:,1);
    idx = (abs(x-1)<eps);
    s(idx) = 1;
    s(~idx) = 0;
    end

    function c = omega(p) % diffusion constant
    global epsilon
    c = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2); z = p(:,3);
    idx = (x>0) & (x<1) & (y>0) & (y<1) & (z>0) & (z<1);
    c(idx) = 1/epsilon;
    c(~idx) = 1;
    end
end

repo: ifem


text:
function assert(condition,message)

if ~condition
  if nargin < 2
    error('Assertion failed');
  end
  error(message)
end

author: Cloud-CV

repo: object-proposals


text:
function raeq = eqxra (tjd, k)

% this function computes the intermediate right ascension
% of the equinox at julian date tjd, using an analytical expression
% for the accumulated precession in right ascension.  for the
% true equinox the result is the equation of the origins.

% input

%  tjd = tdb julian date

%  k = equinox selection code

%      set k = 0 for mean equinox
%      set k = 1 for true equinox (equation of the origins)

% output

%  raeq = intermediate right ascension of the equinox, in hours (+ or -)

% NOTE: this function computes the accumulated precession in right
%       ascension using the analytic equation in capitaine et al. (2003),
%       astronomy and astrophysics 412, 567-586, eq. (42).

% ported from NOVAS 3.0

%%%%%%%%%%%%%%%%%%%%%%%

% t0 = tdb julian date of epoch j2000.0 (tt)

t0 = 2451545.0d0;

t = (tjd - t0) / 36525.0d0;

% for the true equinox, obtain the equation of the equinoxes in time seconds

if (k == 1)

    [a, a, ee, a, a] = etilt (tjd);

    eqeq = ee;

else

    eqeq = 0.0d0;

end

% precession in ra in arcseconds taken from capitaine et al. (2003),
% astronomy and astrophysics 412, 567-586, eq. (42)

precra = 0.014506d0 + ...
    (((( -0.0000000368d0 * t ...
    - 0.000029956d0) * t ...
    - 0.00000044d0) * t ...
    + 1.3915817d0) * t ...
    + 4612.156534d0) * t;

raeq = -(precra / 15.0d0 + eqeq) / 3600.0d0;

author: Sable

repo: mcbench-benchmarks


text:
function varargout = size(cS,varargin)
% overloads size

varargout{1} = length(cS);

author: mtex-toolbox

repo: mtex

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


text:
function uI = edgeinterpolate2(exactu,node,edge,face,face2edge)
%% EDGEINTERPOLATE2 interpolate to the quadratic (1st type) edge finite element space.
%
% uI = edgeinterpolate(u,node,edge) interpolates a given function
% u into the lowesr order edge finite element spaces. The coefficient
% is given by the line integral int_e u*t ds. Simpson rule is used to
% evaluate this line integral.
%
% The input u could be a funtional handel or an array of length N (linear
% element) or N+NE (quadratic element).
%
% Example
%
%   node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%   elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%   maxIt = 3;
%   HcurlErr = zeros(maxIt,1);
%   L2Err = zeros(maxIt,1);
%   N = zeros(maxIt,1);
%   for k = 1:maxIt
%       [node,elem] = uniformbisect3(node,elem);
%       [elem2dof,T] = dof3NE2(elem);
%       pde = Maxwelldata2;
%       uI = edgeinterpolate2(pde.exactu,node,T.edge,T.face,T.face2edge);
%       HcurlErr(k) = getHcurlerror3NE2(node,elem,pde.curlu,uI);
%       L2Err(k) = getL2error3NE2(node,elem,pde.exactu,uI);
%       N(k) = length(uI);
%   end
%   figure(1)
%   r1 = showrate(N,HcurlErr,1,'r-+');
%   hold on
%   r2 = showrate(N,L2Err,1,'b-+');
%   legend('||u-u_I||_{curl}',['N^{' num2str(r1) '}'],...
%          '||u-u_I||',['N^{' num2str(r2) '}'],'LOCATION','Best');
%
% See also edgeinterpolate, edgeinterpolate1
%
% Maxwell doc Section: Dirichlet
% boundary condition
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Evaluate function at vertices and middle points
if isnumeric(exactu)
    uQ = exactu;
else
    mid = (node(edge(:,1),:) + node(edge(:,2),:))/2;
    uQ = exactu([node; mid]);
end

%% Edge dof coefficients
N = size(node,1);	NE = size(edge,1);    NF = size(face,1);
edgeVec = node(edge(:,2),:)-node(edge(:,1),:);
uI = zeros(2*(NE+NF),1);
uI(1:NE,1) = dot(edgeVec,(uQ(edge(:,1),:)+uQ(edge(:,2),:)+4*uQ(N+1:N+NE,:))/6,2);
uI(NE+1:2*NE,1) = dot(edgeVec,0.5*(uQ(edge(:,1),:)-uQ(edge(:,2),:)),2);

%% Face dof coefficients
eik = node(face(:,3),:)-node(face(:,1),:);
eij = node(face(:,2),:)-node(face(:,1),:);
uquadpts = uQ(N+face2edge(:,1),:)+uQ(N+face2edge(:,2),:)+uQ(N+face2edge(:,3),:);
lf = [4*dot(eik,uquadpts,2) 4*dot(eij,uquadpts,2)];
face2edgeDofValue = [uI(face2edge(:,1:3)) uI(face2edge(:,1:3)+NE)];
localMatrix = [4 8 4 -4  0 4; 8 4 -4 0 -4 4]';
lf = lf - face2edgeDofValue*localMatrix;
uI((2*NE+1):end) = [2*lf(:,1) - lf(:,2); 2*lf(:,2) - lf(:,1)]/3;

repo: Aero-matlab


text:
function [vSel, vM]=misd(mCat)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: vMain=misd(mCat)
% Author: van Stiphout, Thomas
% Email: vanstiphout@sed.ethz.ch
% Created: 14. Feb. 2007
% Changed: -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variables:
% mCat              Catalog to be declustered
% iYr1                Starting year (read from mCat)
% iYr2                Ending year (read from mCat)
% fXmeff            Magnitude cutoff
% fRfact             rfact
% fTau0              Tau0 (is equal fTaumin)
% fTaumin          Taumin
% fTaumax          Taumax
% fP1                  P1
% fXk                  xk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% prepare variables
% fTaumin=fTaumin*24*60;
% fTau0=fTaumin;
% fTaumax=fTaumax*24*60;
% sYr1=num2str(floor(min(mCat(:,3))));
% sYr2=num2str(ceil(max(mCat(:,3))));
% sYr1=sYr1(3:4);
% sYr2=sYr2(3:4);
%
% % does catalog contain 12 columns?
% if (size(mCat,2)<12)
%     mTmp=nan(size(mCat,1),(12-size(mCat,2)));
%     mCat=[mCat mTmp];
% end
!rm ~/zmap/src/thomas/decluster/marsan/*.out
!rm ~/zmap/src/thomas/decluster/marsan/misd_input.dat
!rm ~/zmap/src/thomas/decluster/marsan/CORRECT.DAT
% !rm ~/zmap/src/thomas/decluster/marsan/misdMcorr.dat
% export to readalbe format for tmp
export4misd(mCat,'misd_input.dat');
!mv misd_input.dat  ~/zmap/src/thomas/decluster/marsan/.
calc_misdMcorr(mCat,'CORRECT.DAT');
!mv CORRECT.DAT  ~/zmap/src/thomas/decluster/marsan/.
% calc_misdMcorr(mCat,'misdMcorr.dat');
% !mv misdMcorr.dat  ~/zmap/src/thomas/decluster/marsan/.
% save ~/zmap/src/thomas/decluster/marsan/misd_input.dat mCat -ascii

% write input file (input.cmn)
% infile=fopen('input.cmn','w');
% fprintf(infile,'CA.hypo71\n4\n%2s\n%2s\n%3.1f\n%06.3f\n%08.3f\n%010.2f\n%010.2f\n%05.2f\n%05.2f',...
%     sYr1,sYr2,fXmeff,fRfact,fTau0,fTaumin,fTaumax,fP1,fXk);
% fclose(infile);
% run misd algorithm of Marsan2007
sPath=pwd;
cd ~/zmap/src/thomas/decluster/marsan
tic
unix('./misd inputmisd.in');
toc
cd(sPath);
% extract vector with 1 and 0's of events in declustered catalog
% unix(['awk -f ~/zmap/src/thomas/decluster/reasen/clu2list.awk cluster.ano > v.dat ' ]);
% import vector
vM=load('~/zmap/src/thomas/decluster/marsan/W0.out');
vSel=(vM>rand(size(vM,1),1));
% % clean up
% !rm cluster.*


text:
function F = CalHV(points,bounds,k,nSample)
% Calculate the hypervolume-based fitness value of each solution

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

% This function is modified from the code in
% http://www.tik.ee.ethz.ch/sop/download/supplementary/hype/

    [N,M] = size(points);
    if M > 2
        % Use the estimated method for three or more objectives
        alpha = zeros(1,N);
        for i = 1 : k
            alpha(i) = prod((k-[1:i-1])./(N-[1:i-1]))./i;
        end
        Fmin = min(points,[],1);
        S    = unifrnd(repmat(Fmin,nSample,1),repmat(bounds,nSample,1));
        PdS  = false(N,nSample);
        dS   = zeros(1,nSample);
        for i = 1 : N
            x        = sum(repmat(points(i,:),nSample,1)-S<=0,2) == M;
            PdS(i,x) = true;
            dS(x)    = dS(x) + 1;
        end
        F = zeros(1,N);
        for i = 1 : N
            F(i) = sum(alpha(dS(PdS(i,:))));
        end
        F = F.*prod(bounds-Fmin)/nSample;
    else
        % Use the accurate method for two objectives
        pvec  = 1:size(points,1);
        alpha = zeros(1,k);
        for i = 1 : k
            j = 1:i-1;
            alpha(i) = prod((k-j)./(N-j))./i;
        end
        F = hypesub(N,points,M,bounds,pvec,alpha,k);
    end
end

function h = hypesub(l,A,M,bounds,pvec,alpha,k)
% The recursive function for the accurate method

    h     = zeros(1,l);
    [S,i] = sortrows(A,M);
    pvec  = pvec(i);
    for i = 1 : size(S,1)
        if i < size(S,1)
            extrusion = S(i+1,M) - S(i,M);
        else
            extrusion = bounds(M) - S(i,M);
        end
        if M == 1
            if i > k
                break;
            end
            if alpha >= 0
                h(pvec(1:i)) = h(pvec(1:i)) + extrusion*alpha(i);
            end
        elseif extrusion > 0
            h = h + extrusion*hypesub(l,S(1:i,:),M-1,bounds,pvec(1:i),alpha,k);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function readHTKresult(input, output);

INPUT = fopen(input,'r');
OUTPUT = fopen(output,'w');

% Read in the data for one test
for i=1:70
    fsearch2('Date:',INPUT);
    tmp = textscan(INPUT,'%s%s%d%s%d',1);
    data.date.week(i) = tmp{1};
    data.date.month(i) = tmp{2};
    data.date.day(i) = tmp{3};
    data.date.time(i) = tmp{4};
    data.date.year(i) = tmp{5};

    fsearch2('labels/',INPUT);
    tmp = textscan(INPUT,'%2c',1);
    data.label(i) = tmp;

    fsearch2('Correct=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.sent_corr(i) = tmp{1};
    fsearch2('H=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.sent_H(i) = tmp{1};
    fsearch2('S=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.sent_S(i) = tmp{1};
    fsearch2('N=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.sent_N(i) = tmp{1};

    fsearch2('Corr=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_corr(i) = tmp{1};
    fsearch2('Acc=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_acc(i) = tmp{1};
    fsearch2('H=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_H(i) = tmp{1};
    fsearch2('D=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_D(i) = tmp{1};
    fsearch2('S=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_S(i) = tmp{1};
    fsearch2('I=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_I(i) = tmp{1};
    fsearch2('N=',INPUT);
    tmp = textscan(INPUT,'%f',1);
    data.word_N(i) = tmp{1};
end

% write the data
for i=1:70
    if i==1, fprintf(OUTPUT,'\nTEST SET A\n');
    elseif i==29, fprintf(OUTPUT,'\nTEST SET B\n');
    elseif i==57, fprintf(OUTPUT,'\nTEST SET C\n'); end
    if mod(i,7) == 1, fprintf(OUTPUT,'%s\n',data.label{i}); end
%     fprintf(OUTPUT,'%s %s %d %s %d ',...
%         data.date.week{i},data.date.month{i},data.date.day(i),data.date.time{i},data.date.year(i));
    fprintf(OUTPUT,'%s %d %s - ',data.date.month{i},data.date.day(i),data.date.time{i});
    fprintf(OUTPUT,'%2.2f	%2.2f   	%d	%d   	%d   	%d  	%d  - ',data.word_corr(i),data.word_acc(i),data.word_H(i),...
        data.word_D(i),data.word_S(i),data.word_I(i),data.word_N(i));
    fprintf(OUTPUT,'%2.2f %d %d %d \n',data.sent_corr(i),data.sent_H(i),data.sent_S(i),data.sent_N(i));
end
fprintf(OUTPUT,'\n');

fclose(INPUT);
fclose(OUTPUT);

author: singaxiong

repo: SignalGraph


text:
function jed = moon_phase_to_jed ( n, phase )

%*****************************************************************************80
%
%% MOON_PHASE_TO_JED calculates the JED of a moon phase.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2013
%
%  Reference:
%
%    William Press, Brian Flannery, Saul Teukolsky, William Vetterling,
%    Numerical Recipes: The Art of Scientific Computing,
%    Cambridge University Press.
%
%  Parameters:
%
%    Input, integer N, specifies that the N-th such phase
%    of the moon since January 1900 is to be computed.
%
%    Input, integer PHASE, specifies which phase is to be computed.
%    0=new moon,
%    1=first quarter,
%    2=full,
%    3=last quarter.
%
%    Output, real JED, the Julian Ephemeris Date on which the
%    requested phase occurs.
%
  degrees_to_radians = pi / 180.0;
%
%  First estimate.
%
  j = 2415020 + 28 * n + 7 * phase;
%
%  Compute a correction term.
%
  c = n + phase / 4.0;

  t = c / 1236.85;

  xtra = 0.75933 + 1.53058868 * c + ( 0.0001178 - 0.000000155 * t ) * t * t;

  as = degrees_to_radians * ( 359.2242 + 29.105356 * c );

  am = degrees_to_radians * ( 306.0253 + 385.816918 * c + 0.010730 * t * t );

  if ( phase == 0 || phase == 2 )

    xtra = xtra + ( 0.1734 - 0.000393 * t ) * sin ( as ) - 0.4068 * sin ( am );

  elseif ( phase == 1 || phase == 3 )

    xtra = xtra + ( 0.1721 - 0.0004 * t ) * sin ( as ) - 0.6280 * sin ( am );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'MOON_PHASE_TO_JED - Fatal error!\n' );
    fprintf ( 1, '  Illegal PHASE option = %d\n', phase );
    error ( 'MOON_PHASE_TO_JED - Fatal error!' );

  end

  jed = j + xtra;

  return
end


text:
function linpack_d_test07 ( )

%*****************************************************************************80
%
%% TEST07 tests DGBFA and DGBSL.
%
%  Discussion:
%
%    DGBFA and DGBSL are for general banded matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  ml = 25;
  mu = 25;
  lda = 2*ml+mu+1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST07\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  DGBFA factors the matrix,\n' );
  fprintf ( 1, '  DGBSL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to matrix A and right hand side B.
%
%  We want to try a problem with a significant bandwidth.
%
  m = ml + mu + 1;
  fprintf ( 1, '  The bandwidth of the matrix is %d\n', m );

  for j = 1 : n

    ilo = max ( 1, j - mu );
    ihi = min ( n, j + ml );

    temp = 0.0;
    for i = ilo : ihi
      a(i-j+m,j) = -1.0;
      temp = temp - 1.0;
    end

    temp = temp + 1.0;
    a(m,j) = 4.0 - temp;
    b(j) = 4.0;

  end
%
%  Force B to be a column vector.
%
  b = b';
%
%  Factor the matrix A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, ipivot, info ] = dgbfa ( a, lda, n, ml, mu );

  if ( info ~= 0 )
   fprintf ( 1, '  Error!  DGBFA returns INFO = %d\n', info );
    return
  end
%
%  Call DGBSL to solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  job = 0;
  b = dgbsl ( a, lda, n, ml, mu, ipivot, b, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last 5 entries of the solution:\n' );
  fprintf ( 1, '  (All should be 1):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end


text:
function centroids = computeNewCentroids(X, idx, K)

% computeNewCentroids computes the new centroids of each cluster based on
% the mean value of the all the points belonging to that cluster.

% Initialize variables
[m n] = size(X);
centroids = zeros(K, n);

for i=1:K
    temp = find(idx==i);
    Xtemp = X(temp,:); % Get all points belonging to that cluster
    centroids(i,:) = (sum(Xtemp,1))./length(Xtemp); % Assign new centroid based on mean
end


text:
%% Finding contours in an image
% We learn how to find contours of objects in our image.
%
% In this sample you will learn how to use the following OpenCV functions:
%
% * <matlab:doc('cv.findContours') cv.findContours>
% * <matlab:doc('cv.drawContours') cv.drawContours>
%
% Sources:
%
% * <https://docs.opencv.org/3.1.0/df/d0d/tutorial_find_contours.html>
% * <https://github.com/opencv/opencv/blob/3.1.0/samples/cpp/tutorial_code/ShapeDescriptors/findContours_demo.cpp>
%

function varargout = findContours_demo_gui(im)
    % load source image
    if nargin < 1
        src = imread(fullfile(mexopencv.root(),'test','HappyFish.jpg'));
    elseif ischar(im)
        src = imread(im);
    else
        src = im;
    end

    % Convert image to gray and blur it
    if size(src,3) == 3
        src = cv.cvtColor(src, 'RGB2GRAY');
    end
    src = cv.blur(src, 'KSize',[3 3]);

    % create the UI
    h = buildGUI(src);
    if nargout > 0, varargout{1} = h; end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    thresh = round(get(h.slid, 'Value'));
    set(h.txt, 'String',sprintf('Canny thresh: %3d',thresh));

    % Detect edges using canny
    canny_output = cv.Canny(h.src, [thresh thresh*2], 'ApertureSize',3);

    % Find contours
    [contours, hierarchy] = cv.findContours(canny_output, ...
        'Mode','Tree', 'Method','Simple');

    % Draw contours
    drawing = zeros([size(canny_output) 3], 'uint8');
    for i=1:numel(contours)
        clr = randi([0 255], [1 3], 'uint8');
        drawing = cv.drawContours(drawing, contours, ...
            'Hierarchy',hierarchy, 'ContourIdx',i-1, 'MaxLevel',0, ...
            'Color',clr, 'Thickness',2, 'LineType',8);
    end

    % show result
    set(h.img, 'CData',drawing);
    drawnow;
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    thresh = 100;
    max_thresh = 255;
    sz = size(img);
    sz(2) = max(sz(2), 250);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Contours Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2) sz(1)+29]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 30 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img = imshow(img, 'Parent',h.ax);
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax);
        h.img = imshow(img);
    end
    h.txt = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[5 5 130 20], 'String',sprintf('Canny thresh: %3d',thresh));
    h.slid = uicontrol('Parent',h.fig, 'Style','slider', 'Value',thresh, ...
        'Min',0, 'Max',max_thresh, 'SliderStep',[1 10]./(max_thresh-0), ...
        'Position',[135 5 sz(2)-135-5 20]);

    % hook event handlers, and trigger default start
    set(h.slid, 'Callback',{@onChange,h}, ...
        'Interruptible','off', 'BusyAction','cancel');
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef HSVDirectionKey < directionColorKey
  % converts crystal or specimen directions to rgb values
  %
  % The priciple idea is to take the fundamental sector, apply white to the
  % center and red, blue and green to the vertices. This works well if all
  % the edges of the fundamental sector are reflections, i.e. for for m, mm2,
  % mmm, 3m, 4mm, 4/mmm, 6mm, -62m, 6/mmm, -43m, m-3m.
  % In almost all other cases the fundamental sector can be divided by an
  % additional reflection into two subsectors which are colored one with
  % white and one with black center.
  % There are three cases, -1, -3, -4, where this does not work. Actually one
  % can show that in this cases it is impossible to have a smooth one to one
  % relation between the color space and the fundamental sector.

  properties
    colorStretching = 1;
    whiteCenter = vector3d(1,0,0)
    grayValue = [0.2 0.5]
    grayGradient = 0.5 % 0.25
    maxAngle = inf
  end

  properties %(Access = private)
    refl = [];
    rot = rotation.id;
    alpha = 0;
  end

  methods

    function dM = HSVDirectionKey(varargin)

      dM = dM@directionColorKey(varargin{:});
      if ismember(dM.sym.id,[2,18,26])
        warning(['Not a topological correct colormap! Please use the point group ' char(dM.sym.properGroup)]);
      end

      dM.updatesR;
    end

    function rgb = direction2color(dM,h,varargin)

      h = vector3d(h);
      h.antipodal = false;
      h = h.project2FundamentalRegion(dM.sym);

      wC = vector3d(project2FundamentalRegion(dM.whiteCenter,dM.sym));

      if dM.maxAngle < inf

        if angle(vector3d.Z,wC) < 1*degree
          owC = orth(wC);
        else
          owC = vector3d.Z;
        end

        h = project2FundamentalRegion(h,dM.sym,wC);

        rho = angle(h,owC,wC);
        radius = max(0,1 - angle(h,wC) ./ dM.maxAngle);

      else

        switchWB = false;

        % copy to the reduced sector
        h_sR = h;
        whiteOrBlack = true(size(h));
        for i = 1:length(dM.refl)
          ind = dot(h_sR,dM.refl(i)) < 1e-5;
          h_sR(ind) = reflection(dM.refl(i)) * h_sR(ind);
          whiteOrBlack = whiteOrBlack & ~ind;

          if dot(wC,dM.refl(i))<1e-5
            wC = reflection(dM.refl(i)) * wC;
            switchWB = ~switchWB;
          end
        end

        % which are white
        whiteOrBlack = xor(whiteOrBlack,switchWB);

        % compute angle of the points "sh" relative to the center point "center"
        % this should be between 0 and 1
        if isa(dM.sym,'crystalSymmetry')
          ref = vector3d(dM.sym.aAxisRec);
        else
          ref = xvector;
        end
        [radius,rho] = polarCoordinates(dM.sR,h_sR,wC,ref,'maxAngle',dM.maxAngle);

        % white center
        radius(whiteOrBlack) = 0.5+radius(whiteOrBlack)./2;

        % black center
        radius(~whiteOrBlack) = (1-radius(~whiteOrBlack))./2;

      end

      % stretch colors
      radius = radius*(1+dM.alpha)-dM.alpha;

      % compute the color vector on the sphere
      v = vector3d('rho',rho,'theta',radius.*pi);

      % post processing of the color vector
      % by default we have white at the z, black at the -z, red
      % at x and green and blue at 120 and 240 degree accordingly
      % post rotate the color
      v = dM.colorPostRotation * dM.rot * v;
      [th,rh] = polar(v);

      % stretching of the colors
      %th = (th ./ pi) .* pi;

      % correct white -> color gradient
      ind = th > pi/2;

      % the white region
      th(ind) = ((2 * dM.grayGradient(1) * th(ind) ./pi + ...
        (1 - dM.grayGradient(1)) * (1-cos(th(ind)))) ./ 2) .^ dM.colorStretching;

      % the black region
      th(~ind) = ((2 * dM.grayGradient(end) * th(~ind) ./pi + ...
        (1 - dM.grayGradient(end)) * (1-cos(th(~ind)))) ./ 2) .^ dM.colorStretching;

      gray = th;
      gray(ind) = 1 - 2*dM.grayValue(1)*abs(gray(ind) - 0.5);
      gray(~ind) = 1 - 2*dM.grayValue(end)*abs(gray(~ind) - 0.5);

      gray = get_option(varargin,'grayValue',gray);

      % compute rgb values
      rgb = ar2rgb(mod(rh./ 2 ./ pi,1),th,gray);

      rgb(isnan(h.x),:) = NaN;

      %rgb = rgb2gray(rgb);

    end

  end

  methods (Access=protected)

    function updatesR(oM)
      % spherical region to be colorized

      oM.sR = oM.sym.fundamentalSector;
      r30 = rotation.byAxisAngle(zvector,[30,-30]*degree);

      % symmetry dependent settings
      switch oM.sym.id
        case 0
          sR = oM.sym.Laue.fundamentalSector;
          oM.refl = setdiff(sR.N,oM.sR.N);
        case 1                                                   % 1
          if isa(oM.sym,'crystalSymmetry')
            oM.refl = oM.sym.rot.axis;
          else
            oM.refl = vector3d.Z;
          end
        case {3,9}                                               % 211, 112
          oM.refl = -rotate(oM.sR.N,rotation.byAxisAngle(oM.sym.rot(2).axis,90*degree));
        case 6                                                   % 121
          oM.refl = rotate(oM.sR.N,rotation.byAxisAngle(-oM.sym.rot(2).axis,90*degree));
        case {5}, oM.refl = rotate(oM.sR.N(2),-90*degree);       % 2/m11
        case {8}, oM.refl = rotate(oM.sR.N(2),90*degree); %      % 12/m1
        case {11,12}, oM.refl = rotate(oM.sR.N(2),-90*degree); % 222
        case 17, oM.refl = -rotate(sum(oM.sR.N),90*degree);      % 3
        case 18, oM.refl = -rotate(sum(oM.sR.N(2:3)),90*degree); % -3
        case 19
          oM.refl = r30 .* oM.sR.N(end-1:end);                   % 321
          if angle(oM.refl(1),oM.refl(2)) < 1*degree
            oM.refl = inv(r30) .* oM.sR.N(end-1:end);
          end
        case 21, oM.refl =  rotate(sum(oM.sR.N(2:3)),90*degree);    % -31m, -3m1
        case 22, oM.refl =  -rotate(sum(oM.sR.N(2:3)),90*degree);   % 312
        case 24, oM.refl =  -rotate(sum(oM.sR.N(2:3)),90*degree);   % -31m, -3m1
        case {25,27,28}, oM.refl = rotate(oM.sR.N(end),-45*degree); % 4,4/m,422
        case 26, oM.refl = rotate(oM.sR.N(end),-90*degree);      % -4
        case 30, oM.refl = rotate(oM.sR.N(2),45*degree);            % -42m
        case 31, oM.refl = -rotate(oM.sR.N(2),45*degree);        % -4m2
        case {33,35,36}, oM.refl = rotate(oM.sR.N(end),-30*degree); % 6,6/m, 622,
        case 34, oM.refl = rotate(oM.sR.N(end),-60*degree);            % -6
        case {41}, oM.refl = sum(oM.sR.N(3:4))- sum(oM.sR.N(1:2));  % 23
        case {42,43}, oM.refl = oM.sR.N(end-2) - oM.sR.N(end-1);      % 432, m-3
      end

      % reduce fundamental sector by reflectors for black-white colorcoding
      oM.sR.N = [oM.sR.N(:);oM.refl(:)];
      oM.sR.alpha = [oM.sR.alpha(:);zeros(length(oM.refl),1)];
      oM.whiteCenter = oM.sR.center;

    end
  end

  methods (Static = true)
   function rot = green2white % rotate green 2 white
      rot = rotation.byAxisAngle(xvector,90*degree);
    end

    function rot = blue2green % switch blue and green
      rot = reflection(yvector);
    end

    function rot = black2white % switch white and black
      rot = reflection(zvector);
    end
  end
end


text:
function [y,fs,hd,ffx]=readsfs(ff,ty,sub,mode,nmax,nskip,xpath)
%READSFS  Read a .SFS format sound file [Y,FS,HD,FFX]=(FF,TY,SUB,MODE,NMAX,NSKIP,XPATH)
%
% Usage:  [s,fs]=readsfs(filename,1); % read the last speech item in a file
%
% Input Parameters:
%
%  FF gives the name of the file or alternatively
%                 can be the ffx output from a previous call to readsfs
%  TY gives the type of data item required e.g.:
%              0 Main header, 1 Speech data, 2 Laryngograph, 5 Annotation
%  SUB specifies the instance of type TY in the file: 0 for first (default), -1 for last or else
%      it can specify the start of the processing history field as a string (e.g. 'hqtx')
%  MODE		specifies the following (*=default):
%
%           File I/O: 'f'    Do not close file on exit
%             'd'    Look in data directory: voicebox('dir_data')
%           Int Format: 'i'  Force integer data to be at least 16 bits
%                              (some sfs files have a header error which falsely indicates 8-bit data)
%           Create item: 'c' Create item if necessary
%
%  NMAX     maximum number of samples to read (or -1 for unlimited [default])
%  NSKIP    number of samples or frames to skip from start of file
%               (or -1 to continue from previous read when FFX is given instead of a filename [default])
%  XPATH    (used with 'c' option) gives the full name of the program needed to generate the data or
%           the directory containing it.
%
% Output Parameters:
%
%  Y        data matrix or cell matrix whose format depends on TY:
%        TY=0: empty
%			TY=5: cell array {nf,3} = {position length annotation}
%        TY=1,2: column vector containing data
%        TY=11: data array with one row per frame
%  FS       sample frequency in Hz
%  HD     cell matrix whose format depends on TY:
%        TY=0: cell{14,1}
%              {1} row vector
%                {1}(1) = serial_date (see DATENUM() for format)
%                {1}(2) = file_number
%                {1}(3) = machine_type
%              {2} = File type (= 'UC2')
%              {3} = username of creator
%              {4} = site of creator
%              {5} = source
%              {6} = database
%              {7} = speaker name
%              {8} = session code
%              {9} = session date (as a string)
%             {10} = name of token
%             {11} = token repetition code
%             {12} = recording conditions
%             {13} = archiving details
%             {14} = general comments
%        TY>0: cell{4,1}
%              {1} = (1,14) array:
%                 {1}(1)  = processdate (see DATENUM() for format)
%                 {1}(2)  = datatype: 1=speech, 2=lx, 3=tx cycle lengths, 4=fx freq
%                                     5=annotations, 6=phonetic, 7=synthesiser, 8=words
%                                     9=grey-level, 10=voicing, 11=energy, 12=formants
%                                     13=energy, 14=lpc, 15=markov, 16=acoustic, 17=?,
%                                     18=geometry, 19=aerodynamics, 20=articulatory
%                                     21=source, 22=physiological, 23=rational filter
%                                     24=poles/zeros, 25=glottal flow, 26=excitation model
%                                     27=nose, 28=calibration, 29=area
%                 {1}(3)  = subtype
%                 {1}(4)  = floating: 1=float, 0=int, -1=structure+
%                 {1}(5)  = datasize in bytes
%                 {1}(6)  = framesize in units of datasize
%                 {1}(7)  = numframes
%                 {1}(8)  = length in bytes of data
%                 {1}(9)  = frameduration=1/sample_rate
%                 {1}(10) = datapresent: 0=deleted, 1=present, 2=link
%                 {1}(11) = time offset
%                 {1}(12) = windowsize
%                 {1}(13) = overlap
%                 {1}(14) = lxsync
%              {2} = processing history
%              {3} = parameter field
%              {4} = comment
%
%  FFX     cell array containing:
%              {1} = filename
%              {2} = (1,4) = [fid byte_order item_row values_read]
%              {3} = (nitem,5) = one row per item [type subtype length position byteorder]
%              {4} = {nitem,3} = cell: one row per item {processing parameters comment} text strings
%
%
% The SFS (Speech Filing System) is a package mostly written by Mark Huckvale
% and is available for UNIX and PC systems from
% http://www.phon.ucl.ac.uk/resource/sfs/

% Features yet to be implemented:
%
%   (1) If no output parameters are specified, header information will be printed.
%   (2) following link items
%   (3) MODE options:
%                Scaling: 's'    Auto scale to make data peak = +-1
%                         'r'    Raw unscaled data (integer values)
%                         'q'    Scaled to make 0dBm0 be unity mean square
%                         'p' *	Scaled to make +-1 equal full scale
%                Errors   'r'    Return if file is non-existant

%	   Copyright (C) Mike Brookes 1998
%      Version: $Id: readsfs.m 1876 2012-05-25 15:43:29Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<7
    xpath=voicebox('sfsbin');       % path for sfs programs
end
EXESUF=voicebox('sfssuffix');                                              % suffix for executable O/S dependent
if nargin<4 mode='p';
else mode = [mode(:).' 'p'];
end
if nargout==0
    if nargin<2 | ty<=0
        [yy,fs,hd,ffx]=readsfs(ff,0,0,mode);
        fprintf(1,'File: %s\n',ffx{1});
        fprintf(1,'Database: %s, Speaker: %s, Date: %s, Token: %s\n',hd{6},hd{7},hd{9},hd{10});
        lst=ffx{3};
        for it=2:size(lst,1);
            [yy,fs,hd,ffx]=readsfs(ffx,lst(it,1),lst(it,2),mode,0);
            nf=hd{1}(7);
            fd=hd{1}(9);
            fprintf(1,'%3d.%02d %ss @ %sHz (%d frames): %s\n',lst(it,1),lst(it,2),sprintsi(nf*fd),sprintsi(1/fd),nf,ffx{4}{it,1});
        end
    end
else
    it=[];
    xfid=[];                  % xfid will be non-empty second time around
    while (isempty(it))                 % may go round this loop twice
        if ischar(ff)           % If ff is a string we must read file
            if any(mode=='d')
                ff=fullfile(voicebox('dir_data'),ff);
            end
            fid=fopen(ff,'rb','b');
            if fid == -1 error(sprintf('Can''t open %s for input',ff)); end

            t=fread(fid,512,'uchar').';
            if (t(1:3)~='UC2')
                error(sprintf('%s is not an SFS file type UC2',ff))
            end
            itemlist = [0 1 0 0 t(512)];

            proglist={};
            for i=2:200
                pos = ftell(fid);
                [t,n]=fread(fid,512,'uchar');
                if (n < 512) break; end
                mm=pow2(1,8*([0 1 2 3]+(t(512)==0)*[3 1 -1 -3]));
                itemlist(i,:)=[mm*[t(389:392) t(393:396) t(413:416)] pos t(512)];
                if itemlist(i,1)>29
                    error(sprintf('%d is not a valid SFS item type',itemlist(i,1)))
                end
                proglist{i,1}=char(zerotrim(t(1:256)'));
                proglist{i,2}=char(zerotrim(t(257:384)'));
                proglist{i,3}=char(zerotrim(t(437:456)'));
                fseek(fid,itemlist(i,3),'cof');
            end
            ffx={ff; [fid 0 0 0]; itemlist; proglist};
        else
            ffx=ff;
            ff=ffx{1};
            fid=ffx{2}(1);
            if fid<0
                fid=fopen(ffx{1},'rb',char('b'+(ffx{2}(2)~=0)*('l'-'b')));
            end
        end

        % now try to find the requested item

        list=ffx{3};
        if nargin<2 ty=0; end
        if nargin<3 sub=0; end
        if ty<0 ty=list(1,1); end
        if ischar(sub)
            lsub=length(sub);
            proglist=ffx{4};

            for itt=size(proglist,1):-1:2
                if list(itt,1)==ty & length(proglist{itt,1})>=lsub
                    if strcmpi(sub,proglist{itt,1}(1:lsub))
                        it=itt;
                    end
                end
            end
            if (isempty(it))
                if any(mode=='c') & isempty(xfid)      % try to create item if we haven't tried before
                    xfid=-1;
                    if nargin>=7
                        xname=xpath;
                        xfid=fopen(xname);
                    end
                    if xfid<0
                        if any('/\'==xpath(end))        % would be better to use fullfile()
                            xname=[xpath sub EXESUF];
                        else
                            xname=[xpath '/' sub EXESUF];
                        end
                        xfid=fopen(xname);
                    end
                    if xfid<0
                        error(sprintf('Cannot find executable program %s',sub));
                    else
                        fclose(xfid);
                        fclose(fid); % close this file
                        doscom=['cmd /c "' xname '" ' ffx{1}];
                        %fprintf(1,'Executing: %s\n',doscom);
                        if dos(doscom) % run the program
                            error(sprintf('Error running DOS command: %s',doscom));
                        end
                        ff=ffx{1};          % force reread of header information
                    end
                else
                    error(sprintf('Cannot find item %d.%s in file: %s',ty,sub,ff));
                end
            end
        else % numeric subitem specification
            if sub>0
                it = find(list(:,1)==ty & list(:,2)==sub);
            elseif sub==0
                it = min(find(list(:,1)==ty));
            else
                it = max(find(list(:,1)==ty));
            end
            if (isempty(it))
                error(sprintf('Cannot find item %d.%d in file: %s',ty,sub,ff));
            end
        end
    end % loop up to two times while (isempty(it))
    lit=list(it,:);
    if ffx{2}(3)~=it
        ffx{2}(3)=it;
        ffx{2}(4)=0;
    end

    % read the selected item with the correct byte order

    if lit(5)~=ffx{2}(2)
        fclose(fid);
        fid=fopen(ffx{1},'rb',char('b'+(lit(5)~=0)*('l'-'b')));
        ffx{2}(1:2)=[fid lit(5)];
        if fid == -1 error(sprintf('Can''t open %s for input',ff)); end
    end
    fseek(fid,lit(4),'bof');

    y=[];
    fs=0;
    if ~lit(1)                              % read main header
        mb=fread(fid,512,'uchar').';
        if nargout>2
            mc=[1 4; 9 28; 29 32; 37 56; 57 76; 77 96; 97 116; 117 136; 137 296; 297 304; 305 312; 313 332; 333 412];
            hd=cell(14,1);
            hd{1}=[pow2(1,8*([2 3 0 1]+(mb(512)==0)*[1 -1 1 -1]))*[mb(5:8); mb(33:36)].' mb(512)];
            hd{1}(1)=hd{1}(1)/86400+719529;  % convert date format
            for i=1:13
                hd{i+1}=char(zerotrim(mb(mc(i,1):mc(i,2))));
            end
        end
    else

        % read the item header

        hd=cell(4,1);
        hdr=zeros(1,14);

        [str,n]=fread(fid,256,'uchar');
        if (n<256) error(sprintf('Error reading item %d.%d in file: %s',ty,sub,file)); end
        hd{2}=char(zerotrim(str'));

        [str,n]=fread(fid,128,'uchar');
        if (n<128) error(sprintf('Error reading item %d.%d in file: %s',ty,sub,file)); end
        hd{3}=char(zerotrim(str'));

        hdr(1:8)=fread(fid,8,'long');
        hdr(9)=fread(fid,1,'double');
        if hdr(9) fs=1/hdr(9); end
        hdr(10)=fread(fid,1,'long');
        hdr(11)=fread(fid,1,'double');

        [str,n]=fread(fid,20,'uchar');
        if (n<20) error(sprintf('Error reading item %d.%d in file: %s',ty,sub,file)); end
        hd{4}=char(zerotrim(str'));

        [hdr(12:14),n]=fread(fid,3,'long');
        if (n<3) error(sprintf('Error reading item %d.%d in file: %s',ty,sub,file)); end
        fseek(fid,44,'cof');
        hd{1}=hdr;
        hd{1}(1)=hd{1}(1)/86400+719529;  % convert date format

        % now read the actual data

        if nargin<6 nskip=ffx{2}(4);
        elseif nskip<0 nskip=ffx{2}(4);
        end

        ksamples=hdr(7)-nskip;
        if nargin>4
            if nmax>=0
                ksamples=min(nmax,ksamples);
            end
        end

        if ksamples>0
            ffx{2}(4)=nskip+ksamples;
            fsz=hdr(6);
            if(hdr(10)==1)		% data present
                if(hdr(4)>=0)		% non-structured
                    ds=hdr(5);
                    if(hdr(4)>0)
                        if(ds==4) fmt='float';
                        elseif (ds==8) fmt='double';
                        else error('error in sfs file'); end
                    else
                        if(ds==1 & all(mode~='i')) fmt='uchar';
                        elseif(ds<=2) fmt='short'; fsz=ceil(fsz*ds/2);
                        elseif(ds==4) fmt='long';
                        else error('error in sfs file'); end
                    end
                    fseek(fid,lit(4)+512+nskip*fsz*ds,'bof');
                    nd=fsz*ksamples;
                    [y,n]=fread(fid,nd,fmt);
                    if (n<nd) error(sprintf('Error reading item %d.%d in file: %s',ty,sub,file)); end
                    y = reshape(y,fsz,ksamples)';
                else
                    if (hdr(2)==5)
                        y = cell(ksamples,3);
                        for ifr=1:nskip
                            lf=fread(fid,1,'uchar');
                            fseek(fid,lf,'cof');
                        end

                        for ifr=1:ksamples
                            lf=fread(fid,1,'uchar');
                            tdat=fread(fid,2,'long');
                            y(ifr,:)={tdat(1) tdat(2) char(fread(fid,lf-9,'uchar').')};
                            lf=fread(fid,1,'uchar');
                        end
                    else
                        error(sprintf('Cannot convert item %d.%d in file: %s',ty,sub,file));
                    end
                end
            end
        end
    end
    if all(mode~='f') fclose(fid); ffx{2}(1)=-1; end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function test_qr
%TEST_QR test various QR factorization methods
%
% Example:
%   test_qr
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

index = UFget ;
[ignore f] = sort (max (index.nrows,index.ncols)) ;

% f = 276
% f = 706
f = f (1:100) ;

for i = f

    % Prob = UFget (i,index)
    Prob = UFget (i) ;
    disp (Prob) ;
    A = Prob.A ;
    [m n] = size (A) ;
    if (m < n)
        A = A' ;
    end
    [m n] = size (A) ;
    if (sprank (A) < n | ~isreal (A))                                       %#ok
        continue ;
    end

    [V,beta,p,R1,q] = cs_qr(A) ;
    A = A (p,q) ;
    parent = etree (A, 'col') ;                                             %#ok

    R0 = qr (A) ;
    R2 = qr_givens (full (A)) ;
    R3 = qr_givens_full (full (A)) ;

    subplot (2,2,1) ; cspy (R0) ; title ('matlab') ;
    subplot (2,2,2) ; cspy (R3) ; title ('qr-full') ;
    subplot (2,2,3) ; cspy (R2) ; title ('qr-givens') ;
    subplot (2,2,4) ; cspy (R1) ; title ('cs-qr') ;

    e0 = norm (A'*A-R0'*R0,1) / norm (A,1) ;
    e1 = norm (A'*A-R1'*R1,1) / norm (A,1) ;
    e2 = norm (A'*A-R2'*R2,1) / norm (A,1) ;
    e3 = norm (A'*A-R3'*R3,1) / norm (A,1) ;
    fprintf ('error %6.2e %6.2e %6.2e %6.2e\n', e0, e1, e2, e3) ;
    drawnow
    if (e1 > e0*1e3 | e2 > e0*1e3)                                          %#ok
        pause
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%
% Generate mixtures of Gaussian distributions

rot = @(t)[cos(t) sin(t);-sin(t) cos(t)];
cov = @(t,s)rot(-t)*diag(s)*rot(t);

% m=mean, C=cov
randncov = @(m,C,N)sqrtm(C)*randn(2,N)+repmat(m(:),[1 N]);
gaussian = @(m,C,X)1/sqrt(det(C))*exp( inv(C) );

C = cov(.5,[1 .1]);
m = [.5 1];
N = 1000;
X = randncov(m,C,N);

ms = 15;

clf;
plot(X(1,:), X(2,:), '.', 'MarkerSize', ms, 'color', [1 0 0]);
axis equal;

author: optimaltransport


text:
%This is a more real-life example

load DataTest; %load some simulation data for acoustic pressure coming from a focused device

h1=figure;
subplot(1,2,1);

 [X,Y,Z]=meshgrid(Rx,Ry,Rz);
 pf = abs(uf);

[Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(pf,1,0,10,10); %We find the first 10 minima and maxima

ufa=20*log10(abs(uf)/max(abs(uf(:))));

is=isosurface(X,Y,Z,ufa,-6); %we display the focused pressure at -6dB
patch(is,'facecolor','red','edgecolor','none');
view(24,28);
lighting gouraud;
camlight;
daspect([1 1 1]);
title(sprintf('(%3.1f,%3.1f,%3.1f) mm',CoordinatesToEval(nc,1)*1000, CoordinatesToEval(nc,2)*1000,CoordinatesToEval(nc,3)*1000));
xlabel('x (mm)');
ylabel('y (mm)');
zlabel('z (mm)');
for nm=1:5
   %we display the location of the maxima, note that for my needs I
   %switched X<->Y
   posmax=[Rx(MaxPos(nm,2)),Ry(MaxPos(nm,1)),Rz(MaxPos(nm,3))];
   line([posmax(1)-1 posmax(1)+1],[posmax(2) posmax(2)],[posmax(3) posmax(3)],'linewidth',1,'color','b');
   line([posmax(1) posmax(1)],[posmax(2)-1 posmax(2)+1],[posmax(3) posmax(3)],'linewidth',1,'color','b');
   line([posmax(1) posmax(1)],[posmax(2) posmax(2)],[posmax(3)-1 posmax(3)+1],'linewidth',1,'color','b');
   TT=['\bf \leftarrow ' char(nm+96)];
   Align='Cap';
   FSize=12;
   text(posmax(1),posmax(2),posmax(3),TT,'FontSize',FSize,'color','b','VerticalAlignment',Align,'Tag','IgnoreReFormatting','HorizontalAlignment','center');
end
grid on;
subplot(1,2,2);

xlim([0 10]);
ylim([0 10]);
daspect([1 1 1]);

TT={};
TT{1}='First 5 Maxima';
for nm=1:5
    posmax=[Rx(MaxPos(nm,2)),Ry(MaxPos(nm,1)),Rz(MaxPos(nm,3))];
    TT{nm+1}= [char(nm+96) sprintf(': %4.3f%',Maxima(nm)) ', X=' sprintf('%2.1f, ',posmax(1)) ', Y=' sprintf('%2.1f',posmax(2)) ', Z=' sprintf('%2.1f',posmax(3))];
end
Align='Cap';
FSize=12;

t1=text(-1.5,7,TT,'FontSize',FSize,'color','k','VerticalAlignment',Align,'BackgroundColor','w','Edgecolor','k','Margin',4);
set(gca,'visible','off');
set(t1,'visible','on');

%at end, you should see three iosurfaces, where in each isosurface we can
%see that global maxima is inside the -6dB region (as expected for my
%needs) and we can see that the other 4 maxima are located oustide this
%region (what for me is also expected).


text:
function shepard_sr_x3_configure()
    global config;
    % general configurations for both training and testing
    config.input_size = [48 48];
    config.chs = 1;
    config.forward_pass_scheme = {'conv_v_mask_norm', 'conv_v', 'conv_v', 'conv_f', 'out'};
    config.mask_for_SR = 'true';
    config.misc.mask_type = 9; % 3x super-resolution

    config.nonlinearity = 'relu';   % 'relu', 'tanh', 'sigmoid'
    config.output_activation = 'nil';   % 'softmax', 'inherit', 'nil'
    config.cost_function = 'L2 norm'; % 'cross entropy', 'L2 norm'
    config.kernel_size = [8 8; 9 9; 1 1; 8 8];
    config.conv_hidden_size = [16 512 512];
    config.full_hidden_size = [];
    config.output_size = [40 40 1];
    config.batch_size = 10;
    config.compute_device = 'GPU';

    % the following items are only for training
    config.learning_rate = 0.001;
    config.weight_range = 0.03;
    config.decay = 5e-7 / 10;
    config.normalize_init_weights = 0;
    config.dropout_full_layer = 0;
    config.optimization = 'adagrad';
end

author: jimmy-ren


text:
function r4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R4MAT_PRINT_SOME prints out a portion of an R4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%7d :', i );

      for j = j2lo : j2hi
        fprintf ( 1, '%12g  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function S = slmax(A, d)
%SLMAX Compute the maximum of values in subarrays
%
% $ Syntax $
%   - S = slmax(A)
%   - S = slmax(A, d)
%   - S = slmax(A, [d1 d2 ... dk])
%
% $ Arguments $
%   - A:        the input array
%   - d:        the dimensions along which the maximum is searched
%   - S:        the resultant max matrix
%
% $ Description $
%   - S = slmax(A) finds the maximums along column vectors of A. It is
%     equivalent to S = max(A)
%
%   - S = slmax(A, d) finds the maximums along dimension d. It is
%     equivalent to S = max(A, [], d)
%
%   - S = slmax(A, [d1 d2 ... dk]) finds the maximum values along dimension
%     d1, d2, ... dk.
%
% $ History $
%   - Created by Dahua Lin on Nov 19th, 2005
%

%% parse and verify input arguments
if nargin < 2 || isempty(d)
    d = 1;
end

%% compute
if isscalar(d)
    S = max(A, [], d);
else
    k = length(d);
    S = A;
    for i = 1 : k
        S = max(S, [], d(i));
    end
end

author: lmthang


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  battistn@tcnj.edu
% Date Created: May 27th, 2015
% Institution Created: UNC-CH
% Date Modified: July 13, 2021
% Institution Modified: TCNJ
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (battistn@tcnj.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in (x,y) positions of the immersed boundary from .vtk
%           format
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Fdata = read_Force_Scalar_Data_From_vtk(path,simNums,strChoice)

cd(path);

filename = [strChoice '.' num2str(simNums) '.vtk'];  % desired FORCE-LAG-DATA.xxxx.vtk file

fileID = fopen(filename);
if ( fileID== -1 )
    error('\nCANNOT OPEN THE FILE!');
end

str = fgets(fileID); %-1 if eof
if ~strcmp( str(3:5),'vtk');
    error('\nNot in proper VTK format');
end

% read in the header info %
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);
str = fgets(fileID);

% Check whether VTK file has time info. This is a VTK file with time,
%       need to read the next 3 lines to have read in appropriate
%       number of header lines.
if ~strcmp( str(1), 'P')
	str = fgets(fileID);
    str = fgets(fileID);
    str = fgets(fileID);
end

% stores # of Lagrangian Pts. as stated in .vtk file
numLagPts = sscanf(str,'%*s %f %*s',1);

for i=1:numLagPts+5 % +5 to get to FORCE DATA
    str = fgets(fileID);
end

% read in the vertices %
[mat,count] = fscanf(fileID,'%f',numLagPts);
if count ~= numLagPts
   error('\nProblem reading in Lagrangian Pts. Force Data');
end

mat = reshape(mat, 1, count/1); % Reshape vector -> matrix (every 3 entries in vector make into matrix row)
forces = mat';              % Store vertices in new matrix

fclose(fileID);               % Closes the data file.

Fdata = forces(:,1);         % magnitude of the force

cd ..;                        % Change directory back to ../hier_IB2d_data/ directory

clear mat str filename fileID count;

author: nickabattista


text:
function result = warpAffine3B(in,A,badVal,B,interpMethod)
%
% function result = warpAffine3B(in,A,badVal,B)
%
% in: input volume, 3D array
% A: 3x4 affine transform matrix or a 4x4 matrix with [0 0 0 1]
%    for the last row.
% badVal: if a transformed point is outside of the volume, badVal is used
% B:  number of voxels to put in the border (default =0, no border)
%
% result: output volume, same size as in
%
% 10/99, on - added Border parameter to permit nearest neighbor interpolation at the edges
%

if ( ~exist('badVal') | isempty(badVal) )
  badVal=NaN;
end

if ( ~exist('B') | isempty(B) )
  B = 0;
end

if(ieNotDefined('interpMethod'))
    interpMethod = '*linear';
end

if (size(A,1)>3)
  A=A(1:3,:);
end

% original size
[NyO NxO NzO] = size(in);

% if B~=0, put a border by replicating edge voxels
if B~=0
  % put border at each slice
  for k =1:size(in,3)
    inB(:,:,k) = putborde(in(:,:,k),B,B,3);
 end
 % repeat the first and last slices
 in = cat(3,repmat(inB(:,:,1),[1 1 B]), inB, repmat(inB(:,:,end),[1 1 B]));
end

% coordinates corresponding to the volume with borders
[xB,yB,zB]=meshgrid(1-B:size(in,2)-B,1-B:size(in,1)-B,1-B:size(in,3)-B);

% Compute coordinates corresponding to input volume
% and transformed coordinates for result
[xgrid,ygrid,zgrid]=meshgrid(1:NxO,1:NyO,1:NzO);
coords=[xgrid(:)'; ygrid(:)'; zgrid(:)'];
homogeneousCoords=[coords; ones(1,size(coords,2))];
warpedCoords=A*homogeneousCoords;

% Compute result using interp3
result = interp3(xB,yB,zB,in,warpedCoords(1,:),warpedCoords(2,:),warpedCoords(3,:),interpMethod);
result = reshape(result,[NyO NxO NzO]);

% replace NaNs with badval
if(~isnan(badVal))
  NaNIndices = find(isnan(result));
  result(NaNIndices)=badVal*ones(size(NaNIndices));
end
return;

%%% Debug

slice=[1 2 3; 4 5 6; 7 8 9]';
slice=[1 1 1; 3 3 3; 5 5 5]';
input=ones(3,3,4);
for z=1:4
  input(:,:,z)=slice;
end

A= [1 0 0 .5;
    0 1 0 0;
    0 0 1 0;
    0 0 0 1];

A= [1 0 0 .5;
    0 1 0 .5;
    0 0 1 0];

res=warpAffine3(input,A)
resB=warpAffine3B(input,A,NaN,1)

for z=1:4
  input(:,:,z)=z*ones(3,3);
end
A= [1 0 0 0;
    0 1 0 0;
    0 0 1 .5;
    0 0 0 1];
res=warpAffine3(input,A)
resB=warpAffine3B(input,A,NaN,1)

input=rand(5,5,5);
res=warpAffine3(input,eye(4));
resB=warpAffine3B(input,eye(4),NaN,1);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% This file is the main file of using the SINDy-PI method to
% infer the ODE of the double pendulum.
%
% Date: 2019/05/06
% Coded By: K

%% Close all, clear all, clc
close all;clear all; clc;
[status,msg] = mkdir('Results');
addpath('Functions')
set(0,'defaulttextInterpreter','latex')
%% Simulate the double pendulum and gather the simulation data

%Load the system parameters, those parameters are based on the actual system
load("EstimatedValueDou.mat")

%Pendulum mass
m1=parsEsDou(1);m2=parsEsDou(2);

%Pendulum center of mass
a1=parsEsDou(3);a2=parsEsDou(4);

%Pendulum arm inertial
I1=parsEsDou(5);I2=parsEsDou(6);

%Gravity constant and the first pendulum arm length
g=9.81;L=0.2667;

%Dapming ratio
k1=0;k2=0;

%Define the simulation time length
Tf=10;dt=0.001;tspan=0:dt:Tf;
T_test=3;tspan_test=0:dt:T_test;

%Define the inital state of the pendulum: theta1, theta2, dtheta1, dtheta2
state0=[pi+1.2;pi-0.6;0;0];
state0_test=[pi-1;pi-0.4;0.3;0.4];

% Define noise level and add gaussian noise to the data
noise=0;

% Define whehter you have control, if you have it, please define it
Control=0;u=0;

%Define whether you want to shuffel the final data
Shuffle=0;

% Run the ODE files and gather the simulation data
[dData,Data]=Get_Sim_Data(@(t,y)DouPenODE(t,y,m1,m2,a1,a2,L,I1,I2,k1,k2,g),state0,u,tspan,noise,Control,Shuffle);

[dData_test,Data_test]=Get_Sim_Data(@(t,y)DouPenODE(t,y,m1,m2,a1,a2,L,I1,I2,k1,k2,g),state0_test,u,tspan,noise,Control,Shuffle);

%% Plot the data
figure(1)
plot(tspan,Data(:,3),'linewidth',3,'color','black')
box('off')
axis('on')
set(gca,'FontSize',24)

figure(2)
plot(tspan,Data(:,4),'linewidth',3,'color','black')
box('off')
axis('on')
set(gca,'FontSize',24)

% figure(3)
% plot(Data(:,3),dData(:,3),'linewidth',3,'color','black')
% box('off')
% axis('off')
%
% figure(4)
% plot(Data(:,4),dData(:,4),'linewidth',3,'color','black')
% box('off')
% axis('off')

%% Now perform sparse regression of non-linear dynamics

% Get the number of states we have
[dtat_length,n_state]=size(Data);

% Define the control input(Should be zero in our example)
n_control=0;

% Choose whether you want to display actual ODE or not
disp_actual_ode=1;

% If the ODEs you want to display is the actual underlyting dynamics of the
% system, please set actual as 1
actual=1;

% Print the actual ODE we try to discover
Print_ODEs(@(t,y)DouPenODE(t,y,m1,m2,a1,a2,L,I1,I2,k1,k2,g),n_state,n_control,disp_actual_ode,actual);

% Create symbolic states
dz=sym('dz',[n_state,1]);

% Now we first create the parameters of the function right hand side
Highest_Poly_Order_Guess=1;
Highest_Trig_Order_Guess=1;
Highest_U_Order_Guess=0;

% Then create the right hand side library parameters
Highest_Poly_Order=1;
Highest_Trig_Order=4;
Highest_U_Order=0;
Highest_dPoly_Order=1;

%% Define parameters for the sparese regression
lam=[1e-4;5e-4;1e-3;2e-3;3e-3;4e-3;5e-3;6e-3;7e-3;8e-3;9e-3;1e-2;2e-2;3e-2;4e-2;5e-2;...
    6e-2;7e-2;8e-2;9e-2;1e-1;2e-1;3e-1;4e-1;5e-1;6e-1;7e-1;8e-1;9e-1;1;1.5;2;2.5;3;3.5;4;4.5;5;...
    6;7;8;9;10;20;30;40;50;100;200];

N_iter=20;
disp=0;
NormalizeLib=0;

for iter=1:n_state
    fprintf('\n \n Calculating the %i expression...\n',iter)

    % According to the previous parameter generate the left hand side guess
    [LHS_Data,LHS_Sym]=GuessLib(Data,dData(:,iter),iter,u,Highest_Poly_Order_Guess,Highest_Trig_Order_Guess,Highest_U_Order_Guess);

    %Generate the corresponding data
    [SINDy_Data,SINDy_Struct]=SINDyLib(Data,dData(:,iter),iter,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order);

    % Run the for loop and try all the left hand guess
    for i=1:length(LHS_Sym)
        if iter==1 && i==1
            Xi=cell(n_state,length(LHS_Sym),length(lam));
            ODE=cell(n_state,length(LHS_Sym),length(lam));
            ODEs=cell(n_state,length(LHS_Sym),length(lam));
        end

        % Print the left hand side that we are testing
        fprintf('\t Testing the left hand side as %s:\n',char(LHS_Sym{i}))

        % Exclude the guess from SINDy library
        [RHS_Data,RHS_Struct]=ExcludeGuess(SINDy_Data,SINDy_Struct,LHS_Sym{i});

        parfor j=1:length(lam)
            % Select the sparse threashold
            lambda=lam(j);

            % Perform the sparse regression problem
            [Xi{iter,i,j},ODE{iter,i,j}]=sparsifyDynamics(RHS_Data,LHS_Data(:,i),LHS_Sym{i},lambda,N_iter,RHS_Struct,disp,NormalizeLib);

            % Perform sybolic calculation and solve for dX
            digits(6)
            ODE_Guess=vpa(solve(LHS_Sym{i}==ODE{iter,i,j},dz(iter)));

            % Print the discovered ODE
            fprintf(strcat('\t The corresponding ODE we found is: ',char(dz(iter,1)),'=',char((ODE_Guess)),'\n \n'));

            % Store the result
            ODEs{iter,i,j}=ODE_Guess;
        end
    end
end

%% Now generate the ODE function file and test the accuracy of the
% identified system

fprintf('\v Start calculating the best model that could represent the training data...\n \n')
for iter=1:n_state
    % Print which expression are you working on
    fprintf('\t Calculating the best model for the %d expression...\n',iter)

    for i=1:length(LHS_Sym)
        % Print the process
        fprintf('\t Calculating the score of previously found ODE on the test data, %d %% finished. \n',round((i/length(LHS_Sym))*100))

        for j=1:length(lam)
            % If the previous ODE is 0, set the score as NaN, else calculate
            % it.
            if isempty(ODEs{iter,i,j})
                ODE_Not_Exist=1;
                Score(iter,i,j)=NaN;
            else
                % Generate the ODE file
                Generate_ODE_RHS(ODEs{iter,i,j},n_state,n_control);
                % Calculate the accuracy of the file
                Score(iter,i,j)=Get_Score(dData_test(:,iter),Data_test,u,Control,tspan_test,state0_test,Shuffle);
            end
        end

        % Get the best lambda
        [minVal1(iter,i),minIndex1(iter,i)]=min(Score(iter,i,:));

    end

    % Get the best score and use this ODE file
    [minVal2(iter,1),minIndex2(iter,1)]=min(minVal1(iter,:));

    % Store the best ODE
    ODE_Best(iter,1)=ODEs{iter,minIndex2(iter,1),minIndex1(iter,minIndex2(iter,1))};

    % Print the Result
    fprintf('\n\n\n\t The SINDy-PI discovered Best ODE for the %d expression is:\n',iter)
    fprintf('\t %s = %s \n\n\n',char(dz(iter)),char((ODE_Best(iter,1)))')
end

%% Now generate this best guess ODE and print its result
disp_best=1;
if disp_best==1
    fprintf('\n\n\n\v The SINDy-PI discovered Best ODE for the whole system is:\n')
    digits(4)
    for iter=1:n_state
        fprintf(strcat('\v ******\v\t',char(dz(iter)),'=',char(simplify(ODE_Best(iter,1))),'\n'));
    end
end

% Also Print the actual ODE for comparison
digits(4)
fprintf('\n\n\n')
Print_ODEs(@(t,y)DouPenODE(t,y,m1,m2,a1,a2,L,I1,I2,k1,k2,g),n_state,n_control,disp_actual_ode,actual);

%% Get the simulation result
% Generate the ODE file
fprintf('\n\n\n\v Generating the Best Model for comparision...\n')
Generate_ODE_RHS(ODE_Best(:,1),n_state,n_control);

% Define a new test data for the comparison
Noise_test=0;
state0_test=[pi+0.3;pi-0.5;0;0];
[dData_Es,Data_Es]=Get_Sim_Data(@(t,z)Sindy_ODE_RHS(t,z,u),state0_test,u,tspan_test,Noise_test,Control,Shuffle);
[dData_test,Data_test]=Get_Sim_Data(@(t,y)DouPenODE(t,y,m1,m2,a1,a2,L,I1,I2,k1,k2,g),state0_test,u,tspan_test,Noise_test,Control,Shuffle);

%% Save the result
File_Name=strcat('Results/DoublePendulum_NoiseLevel_',num2str(noise),'.mat');
save(File_Name,'Score','ODEs','ODE_Best','state0_test','Data_Es','dData_Es',...
    'Data_test','dData_test','Xi','tspan_test')

%% Plot the simulation data
% Show process
fprintf('\n Simulation finished, plotting the result...... \n')

close all
% Create the new directory to save the plot
[fld_status, fld_msg, fld_msgID]=mkdir('Figures');

%%
close all
figure(1)
plot(tspan_test,Data_test(:,1),'linewidth',4.5,'Color','black')
hold on
plot(tspan_test,Data_Es(:,1),'linewidth',4.5,'linestyle','--','color','blue')
% title('Validation')
% xlabel('Time $(t)$')
% ylabel('$\theta_1$')
% legend('Actual Dynamics','Approximated Dynamics')
set(gca,'XTickLabel',[])
set(gca,'FontSize',34);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
%print('-depsc2', '-loose', 'Figures/DoublePendulum_Theta1.eps');

%
figure(2)
plot(tspan_test,Data_test(:,2),'linewidth',4.5,'Color','black')
hold on
plot(tspan_test,Data_Es(:,2),'linewidth',4.5,'linestyle','--','color','blue')
% title('Validation')
% xlabel('Time $(t)$')
% ylabel('$\theta_2$')
% legend('Actual Dynamics','Approximated Dynamics')
set(gca,'FontSize',34);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/DoublePendulum_Theta2.eps');

%%
figure(3)
plot(Data_test(:,1),dData_test(:,1),'linewidth',4.5,'Color','green')
hold on
plot(Data_Es(:,1),dData_Es(:,1),'linewidth',4.5,'linestyle','--','color','blue')
% title('Phase Plot: $\theta_1$ vs $\dot{\theta_1}$')
% xlabel('$\theta_1$')
% ylabel('$\dot{\theta_1}$')
% legend('Actual Dynamics','Approximated Dynamics')
set(gca,'FontSize',24);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/DoublePendulum_PhasePlot_Theta1_vs_Theta2.eps');

%
figure(4)
plot(Data_test(:,2),dData_test(:,2),'linewidth',4.5,'Color','green')
hold on
plot(Data_Es(:,2),dData_Es(:,2),'linewidth',4.5,'linestyle','--','color','blue')
% title('Phase Plot: $\theta_2$ vs $\dot{\theta_2}$')
% xlabel('$\theta_2$')
% ylabel('$\dot{\theta_2}$')
% legend('Actual Dynamics','Approximated Dynamics')
set(gca,'FontSize',24);
grid on

set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
print('-depsc2', '-loose', 'Figures/DoublePendulum_PhasePlot_dTheta1_vs_dTheta2.eps');

author: dynamicslab

repo: SINDy-PI


text:
function [facecolor, linestyle, marker] = getFaceColorForSRMethod(srMethodName)

    srMethods = {'nn', 'bicubic', 'ebsr', 'scsr', 'nbsrf', 'aplus', 'srcnn', 'drcn', 'vdsr', 'sesr', 'nuisr', 'wnuisr', 'hysr', 'dbrsr', 'l1btv', 'bepsr', 'irwsr', 'bvsr', 'srb', 'vsrnet'};

    % Get face color for the SR method.
    colorMap = jet(30);
    facecolor = [colorMap(1:10,:); colorMap(21:30,:)];
    facecolor = facecolor( strcmp(srMethodName, srMethods), :);

    if nargout > 1
        % Get line style for the SR method.
        linestyle = cellstr(char('-', ':', '-.', '--', '-', ':', '-.', '--', '-', ':', '-.', '--', '-', ':', '-.', '--', '-', ':', '-.', '--'));
        linestyle = linestyle{ strcmp(srMethodName, srMethods) };
    end

    if nargout > 2
        % Get marker for the SR method.
        marker = ['o', 'x', '+', '*', 's', 'd', 'v', '^', '<', '>', 'o', 'x', '+', '*', 's', 'd', 'v', '^', '<', '>'];
        marker = marker(strcmp(srMethodName, srMethods));
    end

author: thomas-koehler

repo: SupER


text:
function msm_to_mm ( output_filename, a, format, type, symmetry )

%*****************************************************************************80
%
%% MSM_TO_MM writes a MATLAB Sparse Matrix to a Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
%    Input, string FORMAT, is 'COORDINATE' for a sparse matrix, or 'ARRAY'
%    for a dense matrix.
%
%    Input, string TYPE, is 'REAL', 'INTEGER', 'COMPLEX', or 'PATTERN'.
%
%    Input, string SYMMETRY is 'GENERAL', 'SYMMETRIC', 'SKEW-SYMMETRIC'
%    or 'HERMITIAN'.
%
  if ( s_eqi ( format, 'coordinate' ) )
    msm_to_mm_coordinate ( output_filename, a, type, symmetry );
  elseif ( s_eqi ( format, 'array' ) )
    msm_to_mm_array ( output_filename, a, type, symmetry );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized choice for FORMAT = "%s".\n', format );
    error ( 'MSM_TO_MM - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [im,regions,gt] = os_get_gt_regions(imdb, imageId)
% Get the image and GT regions for imageId

% read the image
ii = find(imdb.images.id == imageId) ;
im = imread(fullfile(imdb.imageDir, imdb.images.name{ii})) ;

if size(im,3) == 1, im = repmat(im, [1 1 3]) ; end
height = size(im,1) ;
width = size(im,2) ;

% read all corresponding segments
sel = find(imdb.segments.imageId == imdb.images.id(ii)) ;
area = zeros(1, numel(sel)) ;
support = cell(1, numel(sel)) ;
for j = 1:numel(sel)
  if isfield(imdb.segments, 'mask') && ...
      exist(fullfile(imdb.maskDir, imdb.segments.mask{sel(j)}), 'file')
    mask = logical(imread(fullfile(imdb.maskDir, imdb.segments.mask{sel(j)}))) ;
  else
    mask = true(height, width);
  end

  if size(mask,3) > 1, mask = mask(:,:,1) ; end
  area(j) = sum(mask(:)) ;
  support{j} = mask ;
end

% subtract smaller regions from larger ones, in case they overlap
occupied = false(height, width) ;
[~,perm] = sort(area) ;
for i = 1:numel(perm)
  j = perm(i) ;
  support{j} = support{j} & ~occupied ;
  occupied = occupied | support{j} ;
end

% create region structure
regions.basis = zeros(height, width) ;
regions.labels = cell(1, numel(support)) ;
regions.area = zeros(1, numel(support)) ;
regions.segmentIndex = zeros(1, numel(support)) ;
for j = 1:numel(support)
  regions.basis(support{j}) = j ;
  regions.labels(j) = {j} ;
  regions.area(j) = sum(support{j}(:)) ;
  regions.segmentIndex(j) = sel(j) ;
end

% drop empty regions (very uncommon but crashes the code)
ok = false(1, numel(regions.labels)) ;
for r = 1:numel(regions.labels)
  ok(r) = any(ismember(regions.basis(:), regions.labels{r})) ;
end
if ~all(ok)
  warning('Dropping some empty gt regions') ;
end
regions.labels = regions.labels(ok) ;

% now create a segmentation mask with class labels for the regions
if (size(imdb.segments.label, 1) == 1)
    gt = zeros(height, width) ;
    for r = 1:numel(regions.labels)
      mask = ismember(regions.basis, regions.labels{r}) ;
      gt(mask) = imdb.segments.label(regions.segmentIndex(r)) ;
    end
else
    numLabels = size(imdb.segments.label, 1);

    gt = zeros(numLabels, height, width) ;
    for ll = 1 : numLabels
        tmp = zeros(height, width);
        for r = 1:numel(regions.labels)
            mask = ismember(regions.basis, regions.labels{r}) ;
            tmp(mask) = (imdb.segments.label(ll, regions.segmentIndex(r)) == 1);
            gt(ll, :, :) = tmp;
        end
    end
end

if 0
  figure(3) ; clf ;
  [~,gt_] = ismember(gt, find(imdb.meta.inUse)) ;
  subplot(1,2,1) ; imagesc(im) ; title(imdb.images.name{ii}) ;
  subplot(1,2,2) ; image(gt_ + 1);
  colormap([ [1,1,1]; distinguishable_colors(sum(imdb.meta.inUse))]) ;
end

end

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function trl = sampleinfo2trl(data)

% SAMPLEINFO2TRL constructs the trial definition from the sampleinfo, the time axes
% and optionally from the trialinfo
%
% Use as
%   trl = sampleinfo2trl(data)
%
% See also ARTIFACT2BOOLVEC, ARTIFACT2EVENT, ARTIFACT2TRL, BOOLVEC2ARTIFACT, BOOLVEC2EVENT, BOOLVEC2TRL, EVENT2ARTIFACT, EVENT2BOOLVEC, EVENT2TRL, TRL2ARTIFACT, TRL2BOOLVEC, TRL2EVENT

% get the begin and end sample of each trial
begsample = data.sampleinfo(:,1);
endsample = data.sampleinfo(:,2);

% recreate the offset
offset = zeros(numel(data.trial), 1);
for i=1:numel(data.trial)
  offset(i) = time2offset(data.time{i}, data.fsample);
end

if isfield(data, 'trialinfo') && istable(data.trialinfo)
  trl = table(begsample, endsample, offset);
  trl = horzcat(trl, data.trialinfo);
elseif isfield(data, 'trialinfo') && isnumeric(data.trialinfo)
  trl = [begsample endsample offset data.trialinfo];
else
  trl = [begsample endsample offset];
end