text:
function [cor_ini, score_all] = sample_opt_pano(cor_id, im_w, im_h, corn, edg, edg2, options)

	score_all = 0;
        [wall_d, x, f] = pano_line_solver(cor_id, im_w, options); % bos-shape case
        cor_id_t = cor_id;
        d_m = 0;
        % initialization
        for j = 1:2:size(cor_id_t,1)
            theta_y = pi*cor_id_t(j+1,2)/im_h-pi/2;
            theta_y_ = pi/2 - pi*cor_id_t(j,2)/im_h;
            d = abs(cot(theta_y));
            d_m = d_m + d/wall_d((j+1)/2);
        end
        d_m = d_m/4;
        cor_ini = [];
        for j = 1:2:size(cor_id_t,1)
            line_d = d_m*wall_d((j+1)/2);
            line_theta_y = acot(line_d);
            line_y = (line_theta_y + pi/2)*im_h/pi;
            cor_ini = [cor_ini;cor_id_t(j,1) line_y];
            cor_ini = [cor_ini;cor_id_t(j,1) line_y];
        end

        score_btn = interp2(corn,cor_ini(2:2:end,1),cor_ini(2:2:end,2));
        score_btn = sum(log(score_btn));

        % add floor score
        cor_all = [cor_ini(2,:);cor_ini(4,:);cor_ini(4,:);cor_ini(6,:);
               cor_ini(6,:);cor_ini(8,:);cor_ini(8,:);cor_ini(2,:)];
    	[ uv ] = coords2uv( cor_all, im_w, im_h );
    	[ xyz ] = uv2xyzN( uv);
    	[ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
    	im = zeros(im_h, im_w, 3);
    	[ panoEdgeC ] = paintParameterLine_my(lines, im_w, im_h, im);
    	score_fl = zeros(1,4);
    	score_fl(1) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 1)));
    	score_fl(2) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 2)));
    	score_fl(3) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 3)));
    	score_fl(4) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 4)));
    	score_fl = sum(log(score_fl));
    	score_btn = score_btn + score_fl;

        % sampling
        % horizontal, bottoms
        d_m_max = d_m + d_m*0.1;
        d_m_min = d_m - d_m*0.1;
        score_best = -Inf;

        for j = d_m_min:d_m*0.02:d_m_max
            cor_opt = [];
            for k = 1:2:size(cor_id_t,1)
                line_d = j*wall_d((k+1)/2);
                line_theta_y = acot(line_d);
                line_y = (line_theta_y + pi/2)*im_h/pi;
                cor_opt = [cor_opt;cor_id_t(k,1) line_y];
            end
            % compute score
            score = interp2(corn,cor_opt(:,1),cor_opt(:,2));
            score = sum(log(score));

            cor_all = [cor_opt(1,:);cor_opt(2,:);cor_opt(2,:);cor_opt(3,:);
               cor_opt(3,:);cor_opt(4,:);cor_opt(4,:);cor_opt(1,:)];
    		[ uv ] = coords2uv( cor_all, im_w, im_h );
    		[ xyz ] = uv2xyzN( uv);
    		[ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
    		[ panoEdgeC ] = paintParameterLine_my(lines, im_w, im_h, im);
    		score_fl = zeros(1,4);
    		score_fl(1) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 1)));
    		score_fl(2) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 2)));
    		score_fl(3) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 3)));
    		score_fl(4) = max(max(edg2(:,:,3).*(panoEdgeC(:,:,1) == 4)));
    		score_fl = sum(log(score_fl));
    		score = score + score_fl;

            if score > score_best
                score_best = score;
                cor_best = cor_opt;
                d_m_best = j;
            end
        end
        if score_best > score_btn
            cor_ini(2:2:end,:) = cor_best;
        else
           score_best = score_btn;
           d_m_best = d_m;
        end
        score_all = score_all + score_best;
        % horizontal, top
        h_m = 0;
        for j = 1:2:size(cor_id_t,1)
            theta_y_ = pi/2 - pi*cor_id_t(j,2)/im_h;
            d = d_m_best*wall_d((j+1)/2);
            h = d*tan(theta_y_);
            h_m = h_m + h;
        end
        h_m = h_m/4;
        for j = 1:2:size(cor_id_t,1)
            line_d = d_m_best*wall_d((j+1)/2);
            line_y_ = atan(h_m/line_d);
            line_y_ = (pi/2-line_y_)*im_h/pi;
            cor_ini(j,2) = line_y_;
        end
        score_tp = interp2(corn,cor_ini(1:2:end,1),cor_ini(1:2:end,2));
        score_tp = sum(log(score_tp));

        % add ceiling score
        cor_all = [cor_ini(1,:);cor_ini(3,:);cor_ini(3,:);cor_ini(5,:);
               cor_ini(5,:);cor_ini(7,:);cor_ini(7,:);cor_ini(1,:)];
    	[ uv ] = coords2uv( cor_all, im_w, im_h );
    	[ xyz ] = uv2xyzN( uv);
    	[ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
    	im = zeros(im_h, im_w, 3);
    	[ panoEdgeC ] = paintParameterLine_my(lines, im_w, im_h, im);
    	score_cl = zeros(1,4);

    	score_cl(1) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 1)));
    	score_cl(2) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 2)));
    	score_cl(3) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 3)));
    	score_cl(4) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 4)));
    	score_cl = sum(log(score_cl));
    	score_tp = score_tp + 0.5*score_cl;

        h_m_max = h_m + h_m*0.1;
        h_m_min = h_m - h_m*0.1;
        score_best = -Inf;
        for j = h_m_min:h_m*0.02:h_m_max
            cor_opt = [];
            for k = 1:2:size(cor_ini,1)
                line_d = d_m_best*wall_d((k+1)/2);
                line_y_ = atan(j/line_d);
                line_y_ = (pi/2-line_y_)*im_h/pi;
                cor_opt = [cor_opt;cor_ini(k,1) line_y_];
            end
            % compute score
            score = interp2(corn,cor_opt(:,1),cor_opt(:,2));
            score = sum(log(score));
            if 1
            % add ceiling score
        	cor_all = [cor_opt(1,:);cor_opt(2,:);cor_opt(2,:);cor_opt(3,:);
               	cor_opt(3,:);cor_opt(4,:);cor_opt(4,:);cor_opt(1,:)];
    		[ uv ] = coords2uv( cor_all, im_w, im_h );
    		[ xyz ] = uv2xyzN( uv);
    		[ lines ] = lineFromTwoPoint( xyz(1:2:end,:), xyz(2:2:end,:) );
    		im = zeros(im_h, im_w, 3);
    		[ panoEdgeC ] = paintParameterLine_my(lines, im_w, im_h, im);
    		score_cl = zeros(1,4);

    		score_cl(1) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 1)));
    		score_cl(2) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 2)));
    		score_cl(3) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 3)));
    		score_cl(4) = max(max(edg2(:,:,2).*(panoEdgeC(:,:,1) == 4)));
    		score_cl = sum(log(score_cl));
    		score = score + 0.5*score_cl;
    	end

            if score > score_best
                score_best = score;
                cor_best = cor_opt;
                h_m_best = j;
            end
        end
        if score_best > score_tp
            cor_ini(1:2:end,:) = cor_best;
        else
            score_best = score_tp;
            h_m_best = h_m;
        end
        score_all = score_all + score_best;
end

author: zouchuhang

repo: LayoutNet


text:
function [FT,VT,CT]=dualLattice(varargin)

% function [FT,VT,CT]=dualLattice(E,V,shrinkFactor,cladOpt)
% ------------------------------------------------------------------------
%
%
%
% 2021/08/16: KMM Fixed error in cladding offset direction. Changed to "out
% of surface offset"
%
% ------------------------------------------------------------------------

%%

switch nargin
    case 2
        E=varargin{1};
        V=varargin{2};
        shrinkFactor=0.2;
        cladOpt=1;
    case 3
        E=varargin{1};
        V=varargin{2};
        shrinkFactor=varargin{3};
        cladOpt=1;
    case 4
        E=varargin{1};
        V=varargin{2};
        shrinkFactor=varargin{3};
        cladOpt=varargin{4};
end

%The default control parameters
cParDefault.shrinkFactor=0.25;
cParDefault.latticeSide=1;
cParDefault.numDigitKeep=5;
cParDefault.meshType='tri';
cParDefault.indBoundary=[];
cParDefault.hexSplit=0;
cParDefault.hexMethod=2;
cParDefault.elementType=[];

%Complement input structure with default
% [cPar]=structComplete(cPar,cParDefault,0);

%%

[Ec,Vc]=patchDetach(E,V,shrinkFactor);

[F]=element2patch(E);
[F_Ec]=element2patch(Ec);

[Fs,indSort]=sort(F,2);
[~,indUni1,indUni2,uniCount,IND_MAP]=unique_map(Fs,'rows');
IND_MAP=sort(IND_MAP,1,'descend');
IND_MAP=full(IND_MAP(1:2,:))';

I=(1:1:size(F,1))';
INDSort=sub2ind(size(F),I(:,ones(size(indSort,2),1)),indSort);

indSort2=indSort;
indSort2=indSort2(indUni1,:); %Unique
indSort2=indSort2(indUni2,:); %Expand

[~,indUnsort]=sort(indSort2,2);
I=(1:1:size(F,1))';
INDUnsort=sub2ind(size(F),I(:,ones(size(indUnsort,2),1)),indUnsort);

F_Ec=F_Ec(INDSort);
F_Ec=F_Ec(INDUnsort);

F=F(indUni1,:);
if numel(shrinkFactor)>1
    if numel(shrinkFactor)==size(E,1) %Shrink factor specified on faces
        [shrinkFactor_V]=faceToVertexMeasure(E,V,shrinkFactor); %Convert to nodal metric
    else
        shrinkFactor_V=shrinkFactor;
    end
    shrinkFactor_E=mean(shrinkFactor_V(F),2); %Convert to edge metric
else
    shrinkFactor_E=shrinkFactor;
end

Fq=NaN(size(IND_MAP,1),6);
Fq(uniCount==2,:)=[F_Ec(IND_MAP(uniCount==2,1),:) F_Ec(IND_MAP(uniCount==2,2),:)];

if any(uniCount==1)
    [Fcc,Vcc]=patchDetach(F,V,shrinkFactor_E);
    Fcc=Fcc(indUni2,:); %Expand
    Vq=[Vc;Vcc];
    F1=fliplr(F_Ec(IND_MAP(uniCount==1,1),:));
    F2=fliplr(Fcc(IND_MAP(uniCount==1,1),:))+size(Vc,1);
    Fq(uniCount==1,:)=[F1 F2];
else
    Vq=Vc;
end

Fq=[Fq(:,[1 2]) Fq(:,[5 4]);...
    Fq(:,[2 3]) Fq(:,[6 5]);...
    Fq(:,[3 1]) Fq(:,[4 6])];

Fq=fliplr(Fq);

if cladOpt==1
    F=F(indUni2,:); %Expand
    [Fq2,Vq2,Fc,~]=dualClad(F(IND_MAP(uniCount==1,1),:),V,shrinkFactor,1);
    [~,~,Nq2]=patchNormal(Fc,Vq2); %Vertex normals

    faceEdgeLengthSet=mean([ sqrt(sum((Vq2(Fc(:,1),:)-Vq2(Fc(:,2),:)).^2,2)) ...
                             sqrt(sum((Vq2(Fc(:,2),:)-Vq2(Fc(:,3),:)).^2,2)) ...
                             sqrt(sum((Vq2(Fc(:,3),:)-Vq2(Fc(:,1),:)).^2,2))],2);

    vertexEdgeLengthSet=faceToVertexMeasure(Fc,Vq2,faceEdgeLengthSet);

    [~,indUni_F2,~]=unique(F2(:)); %Indices of unique nodes

    %Creating offset nodes
    Vq2i=Vq2;
    Vq2i(Fc(indUni_F2),:)=Vq2i(Fc(indUni_F2),:)+Nq2(Fc(indUni_F2),:).*mean(vertexEdgeLengthSet(indUni_F2,:),2);%mean(edgeLengthSet(indPush,:),2); %Vq(F2(indUni_F2),:);

    Vq3=[Vq2;Vq2i]; %Joining node sets

    Fq2_p=fliplr(Fq2)+size(Vq2,1);
    Fq2_2=[Fq2(:,[3 2]) Fq2(:,[2 3])+size(Vq2,1)];%
    Fq2_3=[Fq2(:,[1 4]) Fq2(:,[4 1])+size(Vq2,1)];%
    Fq2_sides=[Fq2_p; Fq2_2; Fq2_3];
    Fq3=[Fq2_sides; Fq2;];
    Cq3=[2*ones(size(Fq2_sides,1),1); 3*ones(size(Fq2,1),1);];

    Cc=4*ones(size(Fc,1),1);

    FT=[Fq(:,[1 2 3]); Fq(:,[3 4 1]); Fq3(:,[1 2 3])+size(Vq,1); Fq3(:,[3 4 1])+size(Vq,1); Fc+size(Vq,1)+size(Vq2,1);]; %Faces
    VT=[Vq; Vq3]; %Vertices
    CT=[ones(size(Fq,1)*2,1); repmat(Cq3,2,1); Cc; ]; %Color data

    [FT,VT]=mergeVertices(FT,VT); %Merge vertices

else
   FT=[Fq(:,[1 2 3]); Fq(:,[3 4 1]); fliplr(F2)]; %Faces
   VT=Vq; %Vertices
   CT=[ones(size(Fq,1)*2,1);2*ones(size(F2,1),1)]; %Color data
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function RPV01 = rpv01(Term,Settle,RPV_Dates,Basis,LIBOR,probfun,b)
%RPV01 Computes the "Risky PV01", ie, the expected present value of the
% premium leg at 1 bp

RPV_Time = yearfrac(Settle,RPV_Dates,Basis);
PaymentTimes = [RPV_Time(1);diff(RPV_Time)];

RPVQ = probfun(RPV_Time,b);
RPVDF = LIBOR.getDiscountFactors(RPV_Dates);

RPV01 = zeros(size(Term));
for spreadidx=1:length(Term)

    RPVtmpidx = RPV_Time < Term(spreadidx);

    PaymentTimestmp = PaymentTimes(RPVtmpidx);
    RPVDFtmp = RPVDF(RPVtmpidx);
    RPVQtmp = RPVQ(RPVtmpidx);

    RPV01(spreadidx) = sum(PaymentTimestmp.*RPVDFtmp.*(RPVQtmp + ...
        1/2*(([1;RPVQtmp(1:end-1)]) - RPVQtmp)));
end

author: Sable

repo: mcbench-benchmarks


text:
function [hmm,fe] = dropstate(hmm,k,X,T)

K = length(hmm.state);
% if isfield(hmm.train,'grouping')
%     Q = length(unique(hmm.train.grouping));
% else
%     Q = 1;
% end
Q = 1;

no_k = setdiff(1:K,k);
K = K - 1;
hmm.state(k) = [];
hmm.K = K;
hmm.prior.Dir2d_alpha = hmm.prior.Dir2d_alpha(no_k,no_k);
hmm.prior.Dir_alpha = hmm.prior.Dir_alpha(no_k);
if Q==1
    hmm.Dir2d_alpha = hmm.Dir2d_alpha(no_k,no_k);
    hmm.Dir_alpha = hmm.Dir_alpha(no_k);
    hmm.P = hmm.P(no_k,no_k);
    hmm.Pi = hmm.Pi(no_k);
else
    hmm.Dir2d_alpha = hmm.Dir2d_alpha(no_k,no_k,:);
    hmm.Dir_alpha = hmm.Dir_alpha(no_k,:);
    hmm.P = hmm.P(no_k,no_k,:);
    hmm.Pi = hmm.Pi(no_k,:);
end

if nargout>1
   if iscell(T)
       fe = 0;
       for i=1:length(T)
           XX_i = cell(1);
           [X_i,XX_i{1},Y_i]  = loadfile(X{i},T{i},hmm.train);
           data = struct('X',X_i,'C',NaN(sum(T{i})-length(T{i})*hmm.train.order,K));
           [Gamma,~,Xi] = hsinference(data,T{i},hmm,Y_i,[],XX_i);
           fe = fe + sum(evalfreeenergy(X_i,T{i},Gamma,Xi,hmm,Y_i,XX_i,[1 1 1 0 0]));
       end
       fe = fe + sum(evalfreeenergy([],[],[],[],hmm,[],[],[0 0 0 1 1]));
   else
       [Gamma,~,Xi] = hsinference(X,T,hmm);
       fe = sum(evalfreeenergy(X,T,Gamma,Xi,hmm));
   end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [ a, ipvt, info ] = zgefa ( a, lda, n )

%*****************************************************************************80
%
%% ZGEFA factors a complex matrix by Gaussian elimination.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%
%  Parameters:
%
%    Input, complex A(LDA,N), the matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); an upper triangular matrix and the multipliers
%    which were used to obtain it.  The factorization can be written A = L*U
%    where L is a product of permutation and unit lower triangular matrices
%    and U is upper triangular.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO,
%    0, normal value.
%    K, if U(K,K) == 0.0.  This is not an error condition for this
%    subroutine, but it does indicate that ZGESL or ZGEDI will divide by zero
%    if called.  Use RCOND in ZGECO for a reliable indication of singularity.
%

%
%  Gaussian elimination with partial pivoting.
%
  info = 0;

  for k = 1 : n - 1
%
%  Find L = pivot index.
%
    l = izamax ( n-k+1, a(k:n,k), 1 ) + k - 1;
    ipvt(k) = l;
%
%  Zero pivot implies this column already triangularized.
%
    if ( zabs1 ( a(l,k) ) == 0.0 )
      info = k;
      continue;
    end
%
%  Interchange if necessary.
%
    if ( l ~= k )
      temp = a(l,k);
      a(l,k) = a(k,k);
      a(k,k) = temp;
    end
%
%  Compute multipliers
%
    a(k+1:n,k) = - a(k+1:n,k) / a(k,k);
%
%  Row elimination with column indexing
%
    for j = k+1 : n
      t = a(l,j);
      if ( l ~= k )
        a(l,j) = a(k,j);
        a(k,j) = t;
      end
      a(k+1:n,j) = a(k+1:n,j) + t * a(k+1:n,k);
    end

  end

  ipvt(n) = n;

  if ( zabs1 ( a(n,n) ) == 0.0 )
    info = n;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function layer = clean_network_layer(layer)
fields = {'a', 'grad', 'grad_W', 'grad_b', 'acc', 'weights', 'grad_W_raw', 'grad2', 'idx', 'X2', 'ft', 'it', 'ot', 'Ct_raw', 'Ct', 'Ct0', 'ht0', 'post', 'validFrameMask'};
for i=1:length(layer)
    for j=1:length(fields)
        if isfield(layer{i}, fields{j})
            layer{i} = rmfield(layer{i}, fields{j});
        end
    end
end
end

author: singaxiong

repo: SignalGraph


text:
kNum = length(prunedStates);

if kNum > 0
    p_I_G_est = zeros(3,kNum);
    p_I_G_imu = zeros(3, kNum);
    kPlot = zeros(1,kNum);

    for k=1:kNum
        state_k = prunedStates{k}.state_k;

        C_CG = quatToRotMat(prunedStates{k}.q_CG);
        C_CI = quatToRotMat(camera.q_CI);
        p_I_G_est(:,k) = prunedStates{k}.p_C_G - C_CG' * C_CI * camera.p_C_I;
        p_I_G_imu(:,k) = msckfState_imuOnly{state_k}.imuState.p_I_G;

        kPlot(k) = state_k;
    end

    figure(1); clf; hold on;
    plot3(p_I_G_est(1,:),p_I_G_est(2,:),p_I_G_est(3,:),'-b');
    plot3(p_I_G_imu(1,:),p_I_G_imu(2,:),p_I_G_imu(3,:),'-r');
    plot3(r_i_vk_i(1,kPlot),r_i_vk_i(2,kPlot),r_i_vk_i(3,kPlot),'-g');
%     if ~isempty(map)
%         scatter3(map(1,:),map(2,:),map(3,:),'or');
%     end
    xlabel('x');ylabel('y');zlabel('z');
    legend('MSCKF','IMU integration','Ground Truth');
%     legend('MSCKF','IMU integration');
    grid on;

    drawnow;
end

repo: MSCKF


text:
function r = coth(a)
%COTH         Taylor hyperbolic cotangent  coth(a)
%
%Thanks to George Corliss for providing the Taylor expansion
%

% written  06/03/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  K = getappdata(0,'INTLAB_TAYLOR_ORDER');

  r = a;
  ct = a.t;
  N = size(a.t,2);
  r.t(1,:) = coth(a.t(1,:));
  ct(1,:) = 1-r.t(1,:).^2;     % 1+a^2
  r.t(2,:) = ct(1,:) .* a.t(2,:);
  for i=2:K
    ct(i,:) = - sum( r.t(1:i,:).*r.t(i:-1:1,:) , 1 );
    r.t(i+1,:) = sum( ct(1:i,:).*a.t(i+1:-1:2,:).*repmat((i:-1:1)',1,N) , 1 )/i;
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [obX,obY,obXM,obYM,lp] = mrTransInplanes(numofanats,obXM,obYM,lp,curInplane)
%NAME:    [obX,obY,obXM,obYM,lp] = mrTransInplanes(numofanats,obXM,obYM,lp,curInplane)
%AUTHOR:  Poirson
%DATE:	  08.04.96
%PURPOSE: One of a set of routines that allows user
%         to set and select a set of oblique planes in saggital slice.
%	  The routines are mrTransInplanes.m, mrRotInplanes.m,
%	  mrClipInplanes.m, mrSelInplane.m, mrSetupInplanes.m
%         mrSpreadInplanes.m
%HISTORY: Started with mrGetOblPlane from G. Boynton 4/6/96
%NOTES:

global sagwin
figure(sagwin)

% One more pair of points for the perpendicular line
nPtPairs = numofanats + 1;

% See if the user has already been working on some inPlanes
if size(obXM,1) == 0
	disp('You must first create a candidate set of Inplanes');
	disp('Choose Set_Up_Inplanes first');
	return
end

% Center point of the inplanes is the middle of perpendicular line
centerX = mean(obXM(nPtPairs,:));
centerY = mean(obYM(nPtPairs,:));

xlim=get(gca,'XLim');
ylim=get(gca,'YLim');
xt = diff(xlim)*(-1.0) * 0.8;
yt=ylim(1)+diff(xlim)*[0.05,0.125,0.20];
msg(1) = text(xt,yt(1),'Translate Buttons:');
msg(2) = text(xt,yt(2),'Left = Center point');
msg(3) = text(xt,yt(3),'Right= Quit');

button = 0;
while(button~=3)
 [tempx,tempy,button]=mrGinput(1,'cross');
 if (button~=3)
  if (button == 1) | (button == 2)
   obXM = obXM + (tempx-centerX);
   obYM = obYM + (tempy-centerY);
   % Recalculate the center
   centerX = tempx;
   centerY = tempy;
  end

  for i=1:length(lp)
    delete(lp(i));
  end

  for i=1:nPtPairs
   if i == curInplane
        lp(i)=line(obXM(i,:),obYM(i,:),'Color','r');
	% Set obX and obY to new values.
   else
        lp(i)=line(obXM(i,:),obYM(i,:),'Color','b');
   end
  end

 end
end

if (curInplane ~= 0)
	obX = obXM(curInplane,:);
	obY = obYM(curInplane,:);
else
	obX = [0,0];
	obY = [0,0];
end

delete(msg(1));delete(msg(2));delete(msg(3));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = (1/(2*m))*sum(power((X*theta - y),2));

% =========================================================================

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 2 Assignments/Linear Regression with Multiple Variables/mlclass-ex1/computeCostMulti.m


text:
function varargout = invpos(varargin)
% INVPOS  Returns model of 1./x

switch class(varargin{1})

    case 'double' % What is the numerical value of this argument (needed for displays etc)
        varargout{1} = 1./varargin{1};

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            t = varargin{2}; % Second arg is the extended operator variable
            X = varargin{3}; % Third arg and above are the args user used when defining t.

            varargout{1} = (cone([2;X-t],X+t));
            varargout{2} = struct('convexity','convex','monotonicity','decreasing','definiteness','positive','model','graph');
            varargout{3} = X;
        else
        end
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Nv, VX, K, EToV] = MeshReader1DGambit(FileName)

% function [Nv, VX, K, EToV] = MeshReader1DGambit(FileName)
% Purpose  : Read in basic grid information to build grid
%
% NOTE     : gambit(Fluent, Inc) *.neu format is assumed

% Open file and read intro
Fid = fopen(FileName, 'rt');
for i=1:6; line = fgetl(Fid); end

% Find number of nodes and number of elements
dims = fscanf(Fid, '%d');
Nv = dims(1); K = dims(2);
for i=1:2; line = fgetl(Fid); end

% read node coordinates
VX = (1:Nv);
for i = 1:Nv
  line = fgetl(Fid);
  tmpx = sscanf(line, '%lf');
  VX(i) = tmpx(2);
end
for i=1:2; line = fgetl(Fid); end

% read element to node connectivity
EToV = zeros(K, 2);
for k = 1:K
  line   = fgetl(Fid);
  tmpcon = sscanf(line, '%lf');
  EToV(k,1) = tmpcon(4); EToV(k,2) = tmpcon(5);
end

% Close file
st = fclose(Fid);
return

repo: Aero-matlab


text:
function vw = zoom3view(vw, resetFlag);
% Zooms the 3D view
%
% view = zoom3view(view, <resetFlag=0>);
%
% Allows the subject to click twice (corners of a rect) on one of
% the views of the brain, gets the new extent of view,  and sets the other
% axes to be consistent with the new zoom extent.
%
% Since you can only select a rectangular zoom range in one orientation
% (axial, coronal, sagittal), leaving the third dimension unconstrained,
% there are two options for how to update the zoom range. If you use the
% left mouse button when selecting the zoom range, the zoom will affect the
% third dimension, causing an isometric zoom range (equal size in each of
% three dimensions). E.g., if you zoom in on an axial slice, the
% superior <-> inferior range will be affected in addition to a <-> p and
% l <-> r directions. If you right or middle click when selecting the range,
% the zoom range will not be isometric (in the example above, the s <-> i
% range will be unchanged; also, a very narrow rectangle will make the
% other directions uneven, while left-clicking will force each direction to
% match the wider edge of the rectangle). It is only the second mouse click
% that determines whether the zoom is isometric or not.
%
% resetFlag: if 1, will reset the zoom to be the whole view size and
% update the axis bounds of the image. If 0 <default>, will get the zoom
% from the user input.
%
%
% 01/04 ras.
% 08/04 ras. Now makes the zoom range cubic by default (equal FOV
% in all directions) if the left mouse button is used, but uses
% an asymmetric zoom (each sub-view may be rectangular) if another
% button is used.
% 06/02 ras: now just updates zoom w/ AXIS command, rather than
% going through a full refresh. Also added resetFlag, so both
% zoom buttons have callbacks to this function (for Volume/Gray views).

dims = viewGet(vw,'Size');

if exist('resetFlag', 'var') & resetFlag==1
    %%%%%reset the zoom
    zoom = [1 dims(1); 1 dims(2); 1 dims(3)];
else
    %%%%%Get from UI
    [X, Y, button] = ginput(2);
    AX = [X(1) X(2) Y(1) Y(2)];

    % use the identity of the last button press
    % to determine if the axes should be isometric
    if button(end)==1
        isoAxes = 1;
    else
        isoAxes = 0;
    end

    % could use getCurSliceOri to get orientation,
    % but this would mean we'd have to allow zooming
    % only on the selected view axes (which is fairly
    % unimportant on a 3-view) -- so figure it out manually
    if isequal(vw.refreshFn, 'volume3View')
        if vw.ui.axiAxesHandle==gca,       ori = 1;
        elseif vw.ui.corAxesHandle==gca,   ori = 2;
        else                                 ori = 3;
        end
    else
        ori = getCurSliceOri(vw);
    end

    % update the zoom for those dimensions the user restricted
    zoom = vw.ui.zoom;
    switch ori
        case 1, % axial view was zoomed in on
            zoom(2,:) = [AX(3) AX(4)];
            zoom(3,:) = [AX(1) AX(2)];
        case 2, % coronal view was zoomed in on
            zoom(1,:) = [AX(3) AX(4)];
            zoom(3,:) = [AX(1) AX(2)];
        case 3, % sagittal view was zoomed in on
            zoom(1,:) = [AX(3) AX(4)];
            zoom(2,:) = [AX(1) AX(2)];
    end

    % make axes isometric, if selected
    if isoAxes==1
        % this is a bit tricky since we want the widest
        % field-of-view from the two selected directions,
        % but ignore the fov-size for the dimension
        % orthogonal to the orientation (e.g., if selected on
        % an axi slice, ignore the axial zoom range):
        fov = diff([X'; Y'],1,2);
        switch ori
            case 1, fov = [0; fov];
            case 2, fov = [fov(1); 0; fov(2)];
            case 3, fov = [fov; 0];
        end

        [maxFov, whichAxis] = max(fov);

        for ori = 1:3
            centerPt = mean(zoom(ori,:));
            zoom(ori,1) = centerPt - maxFov/2;
            zoom(ori,2) = centerPt + maxFov/2;
        end
    end
end

% bounds check
zoom(zoom<1) = 1;
upperBound = min(zoom(:,2), dims(:));
zoom(:,2) = upperBound;

% set the zoom
vw.ui.zoom = zoom;

% update the axes
if isequal(vw.refreshFn, 'volume3View')
    % 3-view
    axis(vw.ui.axiAxesHandle, [zoom(3,:) zoom(2,:)]);
    axis(vw.ui.corAxesHandle, [zoom(3,:) zoom(1,:)]);
    axis(vw.ui.sagAxesHandle, [zoom(2,:) zoom(1,:)]);
else
    % regular gray/volume view
    switch getCurSliceOri(vw)
        case 1, axis(vw.ui.mainAxisHandle, [zoom(3,:) zoom(2,:)]);
        case 2, axis(vw.ui.mainAxisHandle, [zoom(3,:) zoom(1,:)]);
        case 3, axis(vw.ui.mainAxisHandle, [zoom(2,:) zoom(1,:)]);
    end
end

return


text:
% predicting numerical values using Linear Regression
clear all;
format long
disp('===== Linear Regression ====');
disp('Reading featur vector');

figure;
possiblefeaturizations =  {'all','logmultinomial', 'logmultinomial2', 'logmultinomial3','bernouli', 'tfidf','multinomial'};
possiblefeaturizations =  {'bernouli', 'tfidf','multinomial'};

subplot(6,2,1);

for feat = 1:3
      MSEarray =[];
    elapsedarray =[];
    for crossvalidateIter = 1:10
       (fprintf('%d',crossvalidateIter));

        %disp('Splitting up data into training/test sets');
        [num,txt,raw] = xlsread('data\final106.xls');

        % reading the description of each shoe
        descriptions = raw(2:size(raw,1),2);
        style_ratings = num(1:size(num,1),1);
        comfort_ratings = num(1:size(num,1),4);
        overal_ratings = num(1:size(num,1),5);
        shoe_width = num(1:size(num,1),6);
        shoe_size_rating = num(1:size(num,1),9);
        shoe_arch_rating = num(1:size(num,1),10);

%         % only take m data points
%         m=num_data;
%         descriptions = descriptions(1:m);
%         style_ratings = style_ratings(1:m);
%         comfort_ratings = comfort_ratings(1:m);
%         overal_ratings = overal_ratings(1:m);
%         shoe_width = shoe_width(1:m);
%         shoe_size_rating = shoe_size_rating(1:m);
%         shoe_arch_rating =  shoe_arch_rating(1:m);
        valididx = (~isnan(shoe_width) & ~isnan( shoe_arch_rating)) & ~isnan(shoe_size_rating);

        descriptions = descriptions(valididx);
        style_ratings = style_ratings(valididx);
        comfort_ratings = comfort_ratings(valididx);
        overal_ratings = overal_ratings(valididx);
        shoe_width = shoe_width(valididx);
        shoe_size_rating = shoe_size_rating(valididx);
        shoe_arch_rating =  shoe_arch_rating(valididx);
        Indices = crossvalind('Kfold', sum(valididx), 10);

                %featurization = 'bernouli'%'tfidf'%'tfidf'%'multinomial'%'tfidf' %'multinomial'; % 'bernouli', 'tfidf'
        featurization  = possiblefeaturizations{feat};
        featurs = csvread('data\forWeka_featuresonly.csv');
        featurs = featurs(:,2:size(featurs,2));

        featurs = featurs(valididx,:);
        num_data = size(featurs,1); %5000;

        %disp(sprintf('Number of datapoints %d',num_data))
        if strcmp(featurization,'multinomial')
            %just pass
        elseif strcmp(featurization,'bernouli')
            featurs = bernoulli(featurs);
        elseif strcmp(featurization,'tfidf')
            featurs = tfidf(featurs);
        elseif strcmp(featurization,'logmultinomial')
            featurs = log(featurs+1)./log(2);
        elseif strcmp(featurization,'logmultinomial2')
            featurs = round(log(featurs+1)./log(2));
        elseif strcmp(featurization,'logmultinomial3')
            featurs = (log(featurs+1));
        elseif strcmp(featurization,'all')
            featurs = [bernoulli(featurs), tfidf(featurs), (log(featurs+1))] ;

        end

        size_training = floor(.9*num_data);

        trainingset = featurs(Indices~=crossvalidateIter,:);
        testset = featurs(Indices==crossvalidateIter,:);

        responsevals = [style_ratings, comfort_ratings, overal_ratings,shoe_width,shoe_size_rating,shoe_arch_rating];

        responsevals_training = responsevals(Indices~=crossvalidateIter,:);
        responsevals_test = responsevals(Indices==crossvalidateIter,:);

    %disp('Linear Regression');
    %

    tic;

    predictions = [];
    actual = [];
    for i =1:3
        a = responsevals_training(:,i);
        b = responsevals_test(:,i)';
        regresscoeff = regress(a, trainingset);
        C2 = (regresscoeff'*(testset'));
        predictions = [predictions, C2'];
        predictions( predictions>5)=5;
        predictions( predictions<1)=1;
        actual = [actual, b'];
    end

    MSE = mean(sum(((predictions-actual).^2)'))
     elapsed = toc;
        MSEarray = [MSEarray MSE];
        elapsedarray = [elapsedarray elapsed];

    end
    featurization  = possiblefeaturizations{feat}
    fprintf('Avergae MSE for Linear Regression is %0.10f\n', mean(MSEarray));
    subplot(3,2,feat*2-1)
    plot(MSEarray,'-o');
    title(strcat('MSE for Linear Regression (',featurization,') ', sprintf(' avergae MSE = %0.10f\n', mean(MSEarray))));
    xlabel('10 fold cross-validation (iteration no)')
    ylabel('MSE')
    subplot(3,2,feat*2)
    plot(elapsedarray,'-o');
    title(strcat('Elapsed time for Linear Regression (',featurization, ') ' , sprintf(' avergae elapsed time = %0.10f\n', mean(elapsedarray))));
    xlabel('10 fold cross-validation (iteration no)')
    ylabel('Elapsed Time')
    drawnow;
end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [posterior] = tapas_sem_hier_estimate(data, model, inference, pars)
%% Estimates the hgf using mcmc.
%
% Input
%   data        -- Data of the model. It should be a structure array with
%                  fields y, and u of dimensions n times 1, when n is the
%                  number of subjects.
%   model       -- Standar hgf structure object.
%   inference   -- Inference object. (Optional)
%   pars        -- Parameter structure. (Optional)
%
% Output
%  posterior    -- Structure containing the posterior.
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

n = 2;

n = n + 1;
if nargin < n
    inference = struct();
end

n = n + 1;
if nargin < n
    pars = struct();
end

[pars] = tapas_sem_multiv_pars(data, model, pars);
[data] = tapas_sem_multiv_data(data, model, pars);
[model] = tapas_sem_hier_model(data, model, pars);
[inference] = tapas_sem_hier_inference(inference, pars);

[posterior] = tapas_sem_hier_estimate_interface(data, model, inference);

end

author: translationalneuromodeling

repo: tapas


text:
function [ r,flag ] = reward_grid_world_goal(s,x,x_goal,x_bad,r_goal,r_bad,r_step)
%reward
%
%   input ---------------------------------------------------------
%
%       o x      : (1 x 1), current grid index of the agent.
%
%       o x_goal : (1 x 1), goal grid index.
%
%   output ---------------------------------------------------------
%
%       o r      : (1 x 1), reward value
%

if s == x_goal

    r    = r_goal;
    flag = true;

elseif ismember(s,x_bad)

    r    = r_bad;
    flag = true;

else

    r    = r_step;
    flag = false;

end

end

author: epfl-lasa


text:
function femIF = genP1IFEM3DFace(mesh,fem,femI)
%% Usage: Generate Quadrature Information on Interface Faces Used in PPIFEM
%         Each interface face (triangle) is cut into 3 small triangles
%            A1                A1
%           D  E      or      E  D         equivalent
%          A2  A3            A3   A2
%
%          p = [A1;D;E;A2;A3];
%          t = [1 2 3; 2 3 4; 3 4 5]
%
% femIF.tL --- global index of left element
% femIF.tR --- global index of right element
% femIF.basL --- 3*nfi-4-4 matrix basis functions of left elem on small tri
% femIF.basR --- 3*nfi-4-4 matrix basis functions of right elem on small tri
% femIF.gx --- Gaussian x nodes on small tri
% femIF.gy --- Gaussian y nodes on small tri
% femIF.gz --- Gaussian z nodes on small tri
% femIF.gw --- Gaussian weights on small tri (3 Gaussian pt for linear IFEM)
% femIF.area --- Areas of all small tri
% femIF.normal --- unit normal vector of each small tri

% Last Modified by Xu Zhang 08/07/2020
%%
nfI = -min(mesh.fLoc); % number of interface faces
intFID = find(mesh.fLoc<0);
nfIB = size(find(mesh.f_t(intFID,2)==0),1); % # of interface faces on Boundary
nfII = nfI - nfIB; % # of internal interface faces: two tetra share it L&R
tL = zeros(3*nfII,4); tR = zeros(3*nfII,4);
basL = zeros(3*nfII,4,4); basR = zeros(3*nfII,4,4);
gx = zeros(3*nfII,3); gy = zeros(3*nfII,3); gz = zeros(3*nfII,3);
A = zeros(3*nfII,1); normal = zeros(3*nfII,3);

tB = zeros(3*nfIB,4); basB = zeros(3*nfIB,4,4);
gxB = zeros(3*nfIB,3); gyB = zeros(3*nfIB,3); gzB = zeros(3*nfIB,3);
AB = zeros(3*nfIB,1); normalB = zeros(3*nfIB,3);

id = 0; idB = 0;
for i = 1:nfI
    %% Form small triangular partition on the interface face
    fID = intFID(i); % face index
    f_e = mesh.f_e(fID,:); % three surrounding edge index
    idE = find(mesh.eLoc(f_e)<0); % find index of two interface edges
    tmp = [mesh.e(f_e(idE(1)),:), mesh.e(f_e(idE(2)),:)];
    nd1 = sum(tmp) - sum(unique(tmp)); % the node of two interface edges
    nd2 = sum(mesh.e(f_e(idE(1)),:)) - nd1;
    nd3 = sum(mesh.e(f_e(idE(2)),:)) - nd1;
    p = [mesh.p(nd1,:); mesh.eIntP(-mesh.eLoc(f_e(idE(1))),:); ...
        mesh.eIntP(-mesh.eLoc(f_e(idE(2))),:); mesh.p(nd2,:); mesh.p(nd3,:)];
    t = [1 2 3; 2 3 4; 3 4 5];

    %% gx gy gz on a triangle with 3 internal point, accurate upto pd = 2
    X1 = p(t(:,1),:);  X2 = p(t(:,2),:);   X3 = p(t(:,3),:);
    G = zeros(3,9);  w1 = 2/3;  w2 = 1/6; % see gaussPtri.m
    G(:,[1,4,7]) = w1*X1 + w2*(X2+X3);
    G(:,[2,5,8]) = w1*X2 + w2*(X1+X3);
    G(:,[3,6,9]) = w1*X3 + w2*(X1+X2);

    %% triangle area on three-dimension.
%     x1 = X1(:,1);  y1 = X1(:,2);   z1 = X1(:,3);
%     x2 = X2(:,1);  y2 = X2(:,2);   z2 = X2(:,3);
%     x3 = X3(:,1);  y3 = X3(:,2);   z3 = X3(:,3);
%     AT = 1/2*(((x1-x3).*(y2-y1) - (x1-x2).*(y3-y1)).^2 + ...
%         ((y1-y3).*(z2-z1) - (y1-y2).*(z3-z1)).^2 + ...
%         ((z1-z3).*(x2-x1) - (z1-z2).*(x3-x1)).^2).^(1/2);
    AT = TriArea3D(X1,X2,X3);

    %% Left and Right Element
    tIDL = mesh.f_t(fID,1); % element index of left element
    tIDR = mesh.f_t(fID,2); % element index of right element
    if tIDR > 0 % Internal Face
        tIDLi = -mesh.tLoc(tIDL); % intf elem index of left element
        tIDRi = -mesh.tLoc(tIDR); % intf elem index of right element

        %% Determine piece
        nd1ID = mesh.pLoc(nd1); tLp = femI.plusPC(tIDLi); tRp = femI.plusPC(tIDRi);
        if (nd1ID < 0 && tLp == 1) || (nd1ID > 0 && tLp == 2)
            basL(id+1,:,:) = femI.bas2(tIDLi,:,:);
            basL(id+2,:,:) = femI.bas1(tIDLi,:,:);
            basL(id+3,:,:) = femI.bas1(tIDLi,:,:);
        elseif (nd1ID < 0 && tLp == 2) || (nd1ID > 0 && tLp == 1)
            basL(id+1,:,:) = femI.bas1(tIDLi,:,:);
            basL(id+2,:,:) = femI.bas2(tIDLi,:,:);
            basL(id+3,:,:) = femI.bas2(tIDLi,:,:);
        end
        if (nd1ID < 0 && tRp == 1) || (nd1ID > 0 && tRp == 2)
            basR(id+1,:,:) = femI.bas2(tIDRi,:,:);
            basR(id+2,:,:) = femI.bas1(tIDRi,:,:);
            basR(id+3,:,:) = femI.bas1(tIDRi,:,:);
        elseif (nd1ID < 0 && tRp == 2) || (nd1ID > 0 && tRp == 1)
            basR(id+1,:,:) = femI.bas1(tIDRi,:,:);
            basR(id+2,:,:) = femI.bas2(tIDRi,:,:);
            basR(id+3,:,:) = femI.bas2(tIDRi,:,:);
        end

        gx(id+1:id+3,:) = G(:,1:3);
        gy(id+1:id+3,:) = G(:,4:6);
        gz(id+1:id+3,:) = G(:,7:9);

        A(id+1:id+3,:) = AT;
        normal(id+1:id+3,:) = repmat(mesh.f_norm(fID,:),3,1);

        %% tL and tR, use locID, b/c index on interface cell is different
        temp = fem.t(tIDL,:);
        temp1 = temp(femI.locID(tIDLi,:));
        tL(id+1:id+3,:) = repmat(temp1,3,1);

        temp = fem.t(tIDR,:);
        temp2 = temp(femI.locID(tIDRi,:));
        tR(id+1:id+3,:) = repmat(temp2,3,1);
        id = id+3;

    elseif tIDR == 0 % Boundary Face
        tIDLi = -mesh.tLoc(tIDL); % intf elem index of left element

        %% Determine piece: only one element.
        nd1ID = mesh.pLoc(nd1); tLp = femI.plusPC(tIDLi);
        if (nd1ID < 0 && tLp == 1) || (nd1ID > 0 && tLp == 2)
            basB(idB+1,:,:) = femI.bas2(tIDLi,:,:);
            basB(idB+2,:,:) = femI.bas1(tIDLi,:,:);
            basB(idB+3,:,:) = femI.bas1(tIDLi,:,:);
        elseif (nd1ID < 0 && tLp == 2) || (nd1ID > 0 && tLp == 1)
            basB(idB+1,:,:) = femI.bas1(tIDLi,:,:);
            basB(idB+2,:,:) = femI.bas2(tIDLi,:,:);
            basB(idB+3,:,:) = femI.bas2(tIDLi,:,:);
        end

        gxB(idB+1:idB+3,:) = G(:,1:3);
        gyB(idB+1:idB+3,:) = G(:,4:6);
        gzB(idB+1:idB+3,:) = G(:,7:9);

        AB(idB+1:idB+3,:) = AT;
        normalB(idB+1:idB+3,:) = repmat(mesh.f_norm(fID,:),3,1);

        %% tB use locID, b/c index on interface cell is different
        temp = fem.t(tIDL,:);
        temp1 = temp(femI.locID(tIDLi,:));
        tB(idB+1:idB+3,:) = repmat(temp1,3,1);
        idB = idB+3;
    end
end

femIF = struct('tL',tL,'tR',tR,'tB',tB,'basL',basL,'basR',basR,'basB',basB, ...
    'gx',gx,'gy',gy,'gz',gz,'gxB',gxB,'gyB',gyB,'gzB',gzB,'area',A,'areaB',AB,...
    'gw',[1/3;1/3;1/3],'normal',normal,'normalB',normalB);

repo: ifem


text:
y       = load('data/gas.dat')';
time    = 1960:1/4:1986+3/4;
fprintf(1, '\n');

%% Analysis based on Gaussian model %%
stsm        = ssm_stsm('trend', 'dummy', 4);
[stsm logL] = estimate(y, stsm, 0.01, [], 'fmin', 'bfgs', 'disp', 'off');
fprintf(1, '[Gaussian STSM]\n');
fprintf(1, 'loglikelihood:      %g\n', logL);
fprintf(1, 'Irregular variance: %g\n', stsm.param(1));
fprintf(1, 'Trend variance:     %g\n', stsm.param(2));
fprintf(1, 'Level variance:     %g\n', stsm.param(3));
fprintf(1, 'Seasonal variance:  %g\n\n', stsm.param(4));

[alpha irr] = fastsmo(y, stsm);
ycom        = signal(alpha, stsm);
seas        = ycom(2, :);
figure('Name', 'Gaussian model analysis');
subplot(2, 1, 1), plot(time, seas), title('Seasonal'), xlim([1959 1988]), ylim([-0.9 0.8]);
subplot(2, 1, 2), plot(time, irr), title('Gaussian Irregular'), xlim([1959 1988]), ylim([-0.4 0.4]);
if ispc, set(gcf, 'WindowStyle', 'docked'); end
drawnow;

%% Analysis based on t-model %%
stsmt           = [ssm_t ssm_llt ssm_seasonal('dummy', 4)];
randn('state', [3765023265; 2369472656]);
[stsmt logLt]   = estimate(y, stsmt, [stsm.param(1) 4 stsm.param(2:end)], alpha, 'fmin', 'bfgs', 'disp', 'off');
[alphat irrt]   = fastsmo(y, stsmt);
ycomt           = signal(alphat, stsmt);
lvlt            = ycomt(1, :);
seast           = ycomt(2, :);

fprintf(1, '[t-dist error STSM]\n');
fprintf(1, 'loglikelihood:      %g\n', logLt);
fprintf(1, 't-dist variance:    %g\n', stsmt.param(1));
fprintf(1, 't-dist df:          %g\n', stsmt.param(2));
fprintf(1, 'Trend variance:     %g\n', stsmt.param(3));
fprintf(1, 'Level variance:     %g\n', stsmt.param(4));
fprintf(1, 'Seasonal variance:  %g\n', stsmt.param(5));

figure('Name', 'Component comparison between Gaussian and t-dist. models');
subplot(2, 1, 1), plot(time, seas, 'r', 'DisplayName', 'Gaussian'), hold all, plot(time, seast, 'b', 'DisplayName', 't-distribution'), hold off, title('Seasonal component'), xlim([1959 1988]), ylim([-0.9 0.8]);
subplot(2, 1, 2), plot(time, irr, 'r', 'DisplayName', 'Gaussian'), hold all, plot(time, irrt, 'b', 'DisplayName', 't-distribution'), hold off, title('Irregular component'), xlim([1959 1988]), ylim([-0.4 0.4]);
if ispc, set(gcf, 'WindowStyle', 'docked'); end

fprintf(1, '\n');

repo: REST


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function sp = sp_vector_3d_piola_transform (sp1, sp2, sp3, msh)

warning ('geopdes:obsolete', 'The class SP_VECTOR_3D_PIOLA_TRANSFORM is obsolete. Using SP_VECTOR with the transformation as the last argument. Read the help for the usage')

sp = sp_vector ({sp1, sp2, sp3}, msh, 'div-preserving');

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function[h]=lansey(N)
%LANSEY  The Lansey modification of Cynthia Brewer's "Spectral" colormap.
%   __________________________________________________________________
%
%   *|* lansey.png --- The lansey colormap versus parula and jet.
%   Type 'jhelp lansey' to view this image. *|*
%   __________________________________________________________________
%
%   LANSEY(M) returns an M-by-3 matrix containing the Lansey colormap.
%
%   LANSEY with no arguments returns a colormap having the same number of
%   colors as the colormap of the current figure.
%
%   This colormap is Jonathan Lansey's modification of the 11-division
%   version of Cynthia Brewer's "Spectral" colormap.
%
%   LINSPECER itself is available from
%
%       http://www.mathworks.com/matlabcentral/fileexchange/42673.
%
%   This product includes color specifications and designs developed by
%   Cynthia Brewer (http://colorbrewer.org/).
%
%   LANSEY is a simplified version of LINSPECER by Jonathan Lansey,
%   modified and redistributed in accordance with the copyright policies
%   of LINSPECER and ColorBrewer.org, see LANSEY_COPYRIGHT for details.
%
%   To make LANSEY your default colormap, add to your startup.m file the
%   line "set(0,'DefaultFigureColormap',lansey)".
%
%   'lansey --f' generates the figure shown above.
%
%   Usage: h=lansey(M);
%          colormap lansey
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015 J.M. Lilly --- type 'help jlab_license' for details

if nargin>0
    if strcmpi(N, '--f')
       type makefigs_lansey
       makefigs_lansey;
       return
    end
end

if nargin < 1, N = size(get(gcf,'colormap'),1); end

h = cmap2linspecer(colorm(N));
h = cell2mat(h);

function[vOut] = cmap2linspecer(vIn) % changes the format from a double array to a cell array with the right format
vOut = cell(size(vIn,1),1);
for ii=1:size(vIn,1)
    vOut{ii} = vIn(ii,:);
end

% Jonathan Lansay's notes:
% colorm returns a colormap which is really good for creating informative
% heatmap style figures.
% No particular color stands out and it doesn't do too badly for colorblind people either.
% It works by interpolating the data from the
% 'spectral' setting on http://colorbrewer2.org/ set to 11 colors
% It is modified a little to make the brightest yellow a little less bright.
function[cmap] = colorm(varargin)
n = varargin{1};
%frac=1;
frac=.95; % Slight modification from colorbrewer here to make the yellows in the center just a bit darker
cmapp = [158, 1, 66; 213, 62, 79; 244, 109, 67; 253, 174, 97; 254, 224, 139; 255*frac, 255*frac, 191*frac; 230, 245, 152; 171, 221, 164; 102, 194, 165; 50, 136, 189; 94, 79, 162];
x = linspace(1,n,size(cmapp,1));
xi = 1:n;
cmap = zeros(n,3);
for ii=1:3
    cmap(:,ii) = pchip(x,cmapp(:,ii),xi);
end
cmap = flipud(cmap/255);

author: CHLNDDEV


text:
function output = callmpt(interfacedata)

% This file is kept for MPT2 compatability

% Speeds up solving LPs in mpmilp
global mptOptions
if ~isstruct(mptOptions)
    mpt_error
end

% Convert
Matrices = yalmip2mpt(interfacedata);

% Get some MPT options
options = interfacedata.options;
options.mpt.lpsolver = mptOptions.lpsolver;
options.mpt.milpsolver = mptOptions.milpsolver;
options.mpt.verbose = options.verbose;

if options.savedebug
    save mptdebug Matrices
end

if options.mp.unbounded
    Matrices = removeExplorationConstraints(Matrices);
end

[dummy,un] = unique([Matrices.G Matrices.E Matrices.W],'rows');
Matrices.G = Matrices.G(un,:);
Matrices.E = Matrices.E(un,:);
Matrices.W = Matrices.W(un,:);

if isempty(Matrices.binary_var_index)

    showprogress('Calling MPT',options.showprogress);
    solvertime = clock;
    if options.mp.presolve
        [Matrices.lb,Matrices.ub] = mpt_detect_and_improve_bounds(Matrices,Matrices.lb,Matrices.ub,Matrices.binary_var_index,options);
    end

    if any(Matrices.lb(end-Matrices.nx+1:end) == Matrices.ub(end-Matrices.nx+1:end))
        model = [];
    else
        model = mpt_solvenode(Matrices,Matrices.lb,Matrices.ub,Matrices,[],options);
    end
    solvertime = etime(clock,solvertime);

else
    % Pre-solve required on binary problems
    options.mp.presolve = 1;

    solvertime = tic;

    switch options.mp.algorithm
        case 1
            showprogress('Calling MPT via enumeration',options.showprogress);
            model = mpt_enumeration_mpmilp(Matrices,options);
        case 2
            % Still experimental and just for fun. Not working!
            showprogress('Calling MPT via parametric B&B',options.showprogress);
            model = mpt_parbb(Matrices,options);

       case 3
            showprogress('Calling MPT via delayed enumeration',options.showprogress);
            %Matrices = initialize_binary_equalities(Matrices)
            [Matrices.SOS,Matrices.SOSVariables] =  mpt_detect_sos(Matrices);
            [Matrices.lb,Matrices.ub] = mpt_detect_and_improve_bounds(Matrices,Matrices.lb,Matrices.ub,Matrices.binary_var_index,options);
            model = mpt_de_mpmilp(Matrices,options,[]);

        otherwise
    end
    solvertime = toc(solvertime);
end

if isempty(model)
    model = {model};
end

if options.verbose
    if ~isempty(model{1})
        if length(model) == 1
            disp(['-> Generated 1 partition.'])
        else
            disp(['-> Generated ' num2str(length(model)) ' partitions.'])
        end
    end
end

problem = 0;

% Save all data sent to solver?
if options.savesolverinput
    solverinput.Matrices = Matrices;
    solverinput.options  = [];
else
    solverinput = [];
end

% Save all data from the solver?
% This always done
if options.savesolveroutput
    solveroutput.model = model;
    solveroutput.U = interfacedata.used_variables(Matrices.free_var);%(Matrices.free_var <= length( interfacedata.used_variables)));
    solveroutput.x = interfacedata.used_variables(Matrices.param_var);
else
    solveroutput = [];
end

% Standard interface
Primal      = nan*ones(length(interfacedata.c),1);
Dual        = [];
output = createOutputStructure(Primal,Dual,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

function Matrices = initialize_binary_equalities(Matrices)
binary_var_index = Matrices.binary_var_index;
notbinary_var_index = setdiff(1:Matrices.nu,binary_var_index);
% Detect and extract pure binary equalities. Used for simple pruning
nbin = length(binary_var_index);
only_binary = ~any(Matrices.Aeq(:,notbinary_var_index),2);
Matrices.Aeq_bin = Matrices.Aeq(find(only_binary),binary_var_index);
Matrices.beq_bin = Matrices.beq(find(only_binary),:);

function Matrices = removeExplorationConstraints(Matrices);
candidates = find((~any(Matrices.G,2)) & (sum(Matrices.E | Matrices.E,2) == 1));
if ~isempty(candidates)
    Matrices.bndA = -Matrices.E(candidates,:);
    Matrices.bndb = Matrices.W(candidates,:);
    Matrices.G(candidates,:) = [];
    Matrices.E(candidates,:) = [];
    Matrices.W(candidates,:) = [];
end


text:
function F = triangles_from_edges(E)
  % TRIANGLES_FROM_EDGES Given a graph with undirected edges E, find all
  % 3-cliques (triangle)
  %
  % F = triangles_from_edges(E)
  %
  % Inputs:
  %   E  #E by 2 list of undirected edges
  % Outputs:
  %   F  #F by 3 list of unoriented triangles
  %
  % See also: bfs_orient
  %
  n = max(E(:));
  E2V = sparse(repmat(1:size(E,1),2,1)',E,1,size(E,1),n);
  V2V = sparse([E(:,1) E(:,2)],[E(:,2) E(:,1)],1,n,n);
  % If a clique exists with this edge and some vertex then there will be exactly
  % two paths of length one from this edge to this vertex.
  %
  % If there exists exactly two unique paths from an edge to a vertex then both
  % endpoints must be connected to the edge (no other way to get two paths).
  %
  % 3-clique iff #paths from edge to vertex == 2
  [I,J] = find(E2V*V2V==2);
  F = unique(sort([E(I,:) J],2),'rows');
end

author: alecjacobson

repo: gptoolbox


text:
function mObject3 = waveformMorphing(mObject1,mObject2,mRate);
%   Morphing with minimum information
%   (Actually this is not real morphing.
%   It is simply blending two waveform.)
%   mObject3 = waveformMorphing(mObject1,mObject2,mRate);

%   Designed and coded by Hideki Kawahara
%   27/Feb./2005
%   Copyright(c) 2005, Hideki Kawahara

nLength = max(length(mObject1.waveform),length(mObject2.waveform));
if mObject1.samplingFrequency ~= mObject2.samplingFrequency
    mObject3 = [];
    return
end;
x = zeros(nLength,1);
x(1:length(mObject1.waveform)) = (1-mRate)*mObject1.waveform;
x(1:length(mObject2.waveform)) = mRate*mObject2.waveform + x(1:length(mObject2.waveform));

mObject3=createMobject;
mObject3.waveform = x;
mObject3.samplingFrequency = mObject1.samplingFrequency;

author: HidekiKawahara


text:
function expr = str2mathstr(s, varargin)
    % Converts a Matlab string to a Mathematica string
    %
    % @author ayonga @date 2016-11-23
    %
    % Parameters:
    % s: the character array @type char
    % varargin: conversion options. @type varragin
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    p = inputParser;
    p.addParameter('ConvertString',true,@islogical);
    parse(p, varargin{:});

    if p.Results.ConvertString
        expr = ['"', s, '"'];
    else
        expr = s;
    end

end

author: ayonga

repo: frost-dev

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


text:
function [nx] = tapas_ti_gen_state(y, x, u, theta, ptheta)
%% General method to generate state space.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

gen_state = ptheta.method_state;

ns = size(theta, 1);
nc = size(theta, 2);

nx = cell(ns, nc);

for i = 1:ns
    for j = 1:nc
        nx{i, j} = gen_state(y{i, j}, [], u{i, j}, theta{i, j}, ptheta);
    end
end

end


text:
classdef MOCMA < ALGORITHM
% <multi> <real/integer>
% Multi-objective covariance matrix adaptation evolution strategy

%------------------------------- Reference --------------------------------
% C. Igel, N. Hansen, and S. Roth, Covariance matrix adaptation for multi-
% objective optimization, Evolutionary computation, 2007, 15(1): 1-28.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate the initial individuals in CMA-ES
            Population = Problem.Initialization();
            ptarget    = 1/5.5;
            a          = struct('x',num2cell(Population.decs,2)','psucc',ptarget,'sigma',0.5,'pc',0,'C',eye(Problem.D),'Individual',num2cell(Population));

            %% Optimization
            while Algorithm.NotTerminated([a.Individual])
                % Generate new individuals
                for k = 1 : Problem.N
                    a1(k)            = a(k);
                    a1(k).x          = mvnrnd(a(k).x,a(k).sigma^2*a(k).C,1);
                    a1(k).Individual = Problem.Evaluation(a1(k).x);
                end

                % Update the fitness of each individual
                Q           = [a,a1];
                Population  = [Q.Individual];
                % Penalized fitness for handling box constraints
                PopObj      = Population.objs + repmat(1e-6*sum((cat(1,Q.x)-Population.decs).^2,2),1,Problem.M);
                % Calculate the fitness of each individual
                FrontNo     = NDSort(PopObj,inf);
                CrowdDis    = CrowdingDistance(PopObj,FrontNo);
                [~,rank]    = sortrows([FrontNo;-CrowdDis]');
                [~,fitness] = sort(rank);

                % Update the CMA models
                for k = 1 : Problem.N
                    a(k)  = updateStepSize(a(k),fitness(Problem.N+k)<fitness(k),ptarget);
                    a1(k) = updateStepSize(a1(k),fitness(Problem.N+k)<fitness(k),ptarget);
                    a1(k) = updateCovariance(a1(k),(a1(k).x-a(k).x)/a(k).sigma);
                end

                % Individuals for next generation
                Q = [a,a1];
                a = Q(rank(1:Problem.N));
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function test_ft_plot_ortho

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_ortho

dat = randn(64, 64, 64);

figure
ft_plot_ortho(dat);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

for i = 1:1000
   img = getsnapshot(vid);
   imshow(img,[]);

   imwrite(img,[n2s(i,5) '.jpg'],'jpeg');
      drawnow;
end

repo: OpenTLD


text:
function x = jpeg2k2im(y)
%JPEG2K2IM Decodes an IM2JPEG2K compressed image.
%   X = JPEG2K2IM(Y) decodes compressed image Y, reconstructing an
%   approximation of the original image X.  Y is an encoding structure
%   returned by IM2JPEG2K.
%
%   See also IM2JPEG2K.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

narginchk(1,1);       % Check input arguments

% Get decoding parameters: scale, quantization vector, run-length
% table size, zero run code, end-of-data code, wavelet bookkeeping
% array, and run-length table.
n = double(y.n);
q = double(y.q) / 100;
runs = double(y.runs);
zrc = -double(y.zrc);
eoc = zrc - 1;
s = double(y.s);
s = reshape(s,n + 2,2);

% Compute the size of the wavelet transform.
cl = prod(s(1,:));
for i = 2:n + 1
   cl = cl + 3 * prod(s(i,:));
end

% Perform Huffman decoding followed by zero run decoding.
r = huff2mat(y.huffman);

c = [];   zi = find(r == zrc);    i = 1;
for j = 1:length(zi)
   c = [c r(i:zi(j) - 1) zeros(1,runs(r(zi(j) + 1)))];
   i = zi(j) + 2;
end

zi = find(r == eoc);                % Undo terminating zero run
if length(zi) == 1                  % or last non-zero run.
   c = [c r(i:zi - 1)];
   c = [c zeros(1,cl - length(c))];
else
   c = [c r(i:end)];
end

% Denormalize the coefficients.
c = c + (c > 0) - (c < 0);
for k = 1:n
   qi = 3 * k - 2;
   c = wavepaste('h',c,s,k,wavecopy('h',c,s,k) * q(qi));
   c = wavepaste('v',c,s,k,wavecopy('v',c,s,k) * q(qi + 1));
   c = wavepaste('d',c,s,k,wavecopy('d',c,s,k) * q(qi + 2));
end
c = wavepaste('a',c,s,k,wavecopy('a',c,s,k) * q(qi + 3));

% Compute the inverse wavelet transform and level shift.
x = waveback(c,s,'jpeg9.7',n);
x = uint8(x + 128);

author: dipum

repo: dipum-toolbox

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


text:
function [c, r, t, cpulse, acq_codes] = tapas_physio_read_physlogfiles_philips(log_files, cardiac_modality)
% reads out physiological time series and timing vector depending on the
% MR scanner vendor and the modality of peripheral cardiac monitoring (ECG
% or pulse oximetry)
%
%   [c, r, t, cpulse] = tapas_physio_read_physlogfiles_philips(logfile, vendor, cardiac_modality)
%
% IN
%   log_files                   tapas.log_files; see also tapas_physio_new
%           .respiratory
%           .cardiac
%           .sampling_interval
%           .relative_start_acquisition
%   cardiac_modality
%                       'ecg1_filtered'     filtered 1st ECG channel signal
%                                           (Default)
%                       'ecg2_filtered'     filteered 2nd ECG channel
%                                           (sometimes less gradient artifacts)
%                       'ecg1_raw'          raw 1st ECG channel
%                       'ecg2_raw'          raw 2nd ECG channel
%                       'OXY'/'PPU'         pulse plethysmographic unit
%                                           (PPU) signal
%
%                       Note: for ECG, append '_wifi'
%                       for adjusted sampling rate in
%                       wireless Ingenia Scanners
%
%
% OUT
%   c                   cardiac time series (ECG or pulse oximetry)
%   r                   respiratory time series
%   t                   vector of time points (in seconds)
%   cpulse              time events of R-wave peak in cardiac time series (seconds)
%   acq_codes           slice/volume start events marked by number <> 0
%                       for time points in t
%                       10/20 = scan start/end;
%                       1 = ECG pulse; 2 = OXY max; 4 = Resp trigger;
%                       8 = scan volume trigger
%
% EXAMPLE
%   [ons_secs.cpulse, ons_secs.rpulse, ons_secs.t, ons_secs.c] =
%   tapas_physio_read_physlogfiles(logfile, vendor, cardiac_modality);
%
%   See also tapas_physio_main_create_regressors

% Author: Lars Kasper
% Created: 2013-02-16
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

%% read out values
hasCardiac  = ~isempty(log_files.cardiac);
hasResp     = ~isempty(log_files.respiration);

if hasCardiac
    logfile = log_files.cardiac;
else
    logfile = log_files.respiration;
end

if hasCardiac || hasResp
    y = tapas_physio_read_physlogfiles_philips_matrix(logfile);
    acq_codes   = y(:,10);
else
    y = [];
    acq_codes = [];
end

Nsamples    = size(y,1);

dt          = log_files.sampling_interval;

%default: 500 Hz sampling frequency
isWifi      = regexpi(cardiac_modality, '_wifi');

if isWifi
    cardiac_modality = regexprep(cardiac_modality, '_wifi', '', 'ignorecase');
end

if isempty(dt)
    if isWifi
        dt = 1/496;
    else
        dt = 2e-3;
    end
end

t= -log_files.relative_start_acquisition + ((0:(Nsamples-1))*dt)';

% column 3 = ECG, 5 = PPU, 6 = resp,
% 10 = scanner signal: 10/20 = scan start/end; 1 = ECG pulse; 2 = OXY max; 8 = scan event TODO: what is 3 and 9???
% columns 7,8,9: Grad-strengh x,y,z

cpulse = find(acq_codes==1);
if ~isempty(cpulse)
    cpulse = t(cpulse);
end;

if hasResp
    r = y(:,6);
else
    r = [];
end

if hasCardiac

    useDefaultEcgModality = strcmpi(cardiac_modality, 'ecg');

    if useDefaultEcgModality
        cardiacModalityArray = ...
            {'ecg1_filtered', 'ecg2_filtered', 'ecg1_raw', 'ecg2_raw', 'ppu'};
    else
        cardiacModalityArray = {cardiac_modality};
    end

    hasValidCardiacReadout = false;
    iModality = 0;

    % cycle through modalities until non-zero cardiac column found
    while ~hasValidCardiacReadout
        iModality = iModality + 1;
        cardiac_modality = cardiacModalityArray{iModality};
        switch lower(cardiac_modality)
            case {'ecg_raw', 'ecg1_raw'}
                c = y(:,1);
            case {'ecg2_raw'}
                c = y(:,2);
            case {'ecg1', 'ecg_filtered', 'ecg1_filtered'}
                c = y(:,3);
            case { 'ecg2', 'ecg2_filtered'}
                c = y(:,4);
            case {'oxy','oxyge', 'ppu'}
                c = y(:,5);
        end

        hasValidCardiacReadout = any(c);

    end
else
    c = [];
    cpulse = [];
end

end


text:
function extractSubsetFeaturesfromFile(fileName, featureList)
% This function takes a tab delimited file of features (colums) and samples
% (rows) as input and extracts a desired subset of features and their
% corresponding value for each sample.
%
% INPUT
% fileName      file name
% featureList   list of features (must be overlapping with the feature
%               names given in the first row of the file)
%
% OUTPUT
% an output file will be written containing only the subset - the data will
% be comma separated
%
% Ines Thiele, July '22

clear keep colW tline rows fid

fid=fopen(fileName);
tline = fgetl(fid);

cnt = 1;
while ischar(tline)
    %disp(tline);
    % split line
    if cnt == 1 % header
        col = split(tline,'	');
        % get the wanted colum numbers
        colW = find(ismember(col,featureList));
        colW = [1;colW]; % add the first col with data point name
        keep{cnt,1} = col{colW(1)};
        s = ',';
        for i = 2 : length(colW)
            keep{cnt} = strcat(keep{cnt},s,col{colW(i)});
        end
        cnt = cnt + 1;
    else
        %now get only the wanted cols
        rows = split(tline,'	');
        keep{cnt,1} = rows{colW(1)};
        s = ',';
        for i = 2 : length(colW)
            keep{cnt} = strcat(keep{cnt},s,rows{colW(i)});
        end
        cnt = cnt + 1;
    end
    tline = fgetl(fid);
    if mod(cnt,1000)==1
        cnt
    end
end
fclose(fid);

fid =fopen([fileName '_subset' '.txt'], 'w');
for j = 1 : length(keep)-1
    fprintf(fid,strcat(keep{j},'\n'));
end

author: opencobra

repo: cobratoolbox


text:
function r4_normal_01_test ( )

%*****************************************************************************80
%
%% R4_NORMAL_01_TEST tests R4_NORMAL_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4_NORMAL_01_TEST\n' );
  fprintf ( 1, '  R4_NORMAL_01 computes a sequence of \n' );
  fprintf ( 1, '  normally distributed pseudorandom numbers.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ r, seed ] = r4_uniform_01 ( seed );
    fprintf ( 1, '  %6d  %14f\n', i, r );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bnt_data, old_values] = transform_data_into_bnt_format(data,cnodes)
% TRANSFORM_DATA_TO_BNT_FORMAT Ensures discrete variables have values 1,2,..,k
% e.g., if the values of a discrete are [0 1 6], they must be mapped to [1 2 3]
%
% data(i,j) is the value for i-th node in j-th case.
% bnt_data(i,j) is the new value.
% old_values{i} are the original values for node i.
% cnodes is the list of all continous nodes, e.g. [3 5] means the 3rd and 5th node is continuous
%
% Author: yimin.zhang@intel.com
% Last updated: Jan. 22, 2002 by Kevin Murphy.

num_nodes=size(data,1);
num_cases=size(data,2);
old_values=cell(1,num_nodes);

for i=1:num_nodes
  if (myismember(i,cnodes)==1)  %cts nodes no need to be transformed
    %just copy the data
    bnt_data(i,:)=data(i,:);
    continue;
  end
  values = data(i,:);
  sort_v = sort(values);
  %remove the duplicate values in sort_v
  v_set = unique(sort_v);

  %transform the values
  for j=1:size(values,2)
    index = binary_search(v_set,values(j));
    if (index==-1)
      fprintf('value not found in tranforming data to bnt format.\n');
      return;
    end
    bnt_data(i,j)=index;
  end
  old_values{i}=v_set;
end

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

function index=binary_search(vector, value)
% BI_SEARCH do binary search for value in the vector
% Author: yimin.zhang@intel.com
% Last updated: Jan. 19, 2002

begin_index=1;
end_index=size(vector,2);
index=-1;
while (begin_index<=end_index)
  mid=floor((begin_index+end_index)/2);
  if (isstr(vector(mid)))
    % need to write a strcmp to return three result (< = >)
  else
    if (value==vector(mid))
      index=mid;
      return;
    elseif (value>vector(mid))
      begin_index=mid+1;
    else
      end_index=mid-1;
    end
  end
end
return;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym atanh (@var{x})
%% Symbolic atanh function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = atanh (x)
%%   @result{} y = (sym) atanh(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = atanh(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('atanh', x);
end

%!error atanh (sym(1), 2)
%!assert (isequaln (atanh (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1/2;
%! x = sym('1/2');

%!test
%! f1 = atanh(x);
%! f2 = atanh(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = atanh(A);
%! f2 = atanh(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = atanh (d);
%! f = atanh (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function [rootexpdirs,jabfiles,T0,T1] = ethogram_plot(rootexpdirs,jabfiles,nframesplot,varargin)
% ethogram_plot(rootexpdirs,jabfiles,nframesplot,varargin)
% Plots experiment ethograms.
%
% rootexpdirs: cellstr, root directories for experiments to include in
% plot. Eg:
% rootexpdirs = {
%   '/groups/branson/home/kabram/matlab/adamTracking/videos/20130628/Laser On'
%   '/groups/branson/home/kabram/matlab/adamTracking/videos/20130628/Laser Off'
%   };
% jabfiles: cellstr, jab files for classifiers of interest. Eg:
% jabfiles = {
%   '/groups/branson/home/kabram/matlab/adamTracking/data/reach.jab'
%   '/groups/branson/home/kabram/matlab/adamTracking/data/grab.jab'
%   };
% nframesplot: Show this many frames in ethogram.
% varargin: optional P-V pairs:
%  - 'xlsexport', scalar logical.
%  - 'xlsfile', string specifying exported filename for xlsexport.
%  - 'automarks', scalar logical. If true, include automatic predictions
%     for "first lift", etc. In this case, the specified jabfiles will need
%     to include Lift/Handopen/Grab/Sup/Atmouth/Chew.
%  - 'exptags', Either 'readfromjabs' OR (for explicit specification) cell
%     of row cellstrs. If 'readfromjabs', ethogram_plot will look in the
%     jabfiles to find experiment tags. In the case of explicit
%     specification, exptags must have same size as rootexpdirs; exptags{i}
%     is a cellstr specifying tags for all experiments under rootexpdirs{i}.
%
% Example: ethogram_plot(expdirs,jabfiles,nframesplot,'xlsexport',true)
%
% ethogram_plot will run classifyMovie on experiments that do not have the
% appropriate scorefiles, or whose scorefiles are out of date with respect
% to their classifiers.
%
% Double clicking: If you double click on any patch, the corresponding
% experiment and the behavior detector will be opened up in JAABA.
% Shift clicking: Shift-clicking zooms in on an experiment.

[xlsexport,xlsfile,doforce,automarks,smartorder,scoresnotjabs,dontreclassify,exptags,ethcoreargexport] =...
  myparse(varargin,...
  'xlsexport',false,...
  'xlsfile',[],...
  'doforce',false,...
  'automarks',false,...
  'smartorder',true,...
  'scoresnotjabs',false,...
  'dontreclassify',false,...
  'exptags',[],...
  'ethcoreargexport',false);

if nargin < 1,
  rootexpdirs = {};
end
if nargin < 2,
  jabfiles = {};
end
if nargin < 3
  nframesplot = [];
elseif ~isempty(nframesplot)
  assert(isnumeric(nframesplot));
  assert(isscalar(nframesplot) || numel(nframesplot)==2,'Invalid nframesplot');
end

if isempty(rootexpdirs),
  rootexpdirs = uipickfiles('Prompt','Select root directories containing experiments','DirsOnly',true);
  if ~iscell(rootexpdirs) || isempty(rootexpdirs),
    return;
  end
end

% ensure jabfiles
if isempty(jabfiles)
  [tfsuccess,jabfiles] = ExpPP.uiGetJabFiles;
  if ~tfsuccess
    return;
  end
end

% determine exptag status
if ischar(exptags) && strcmp(exptags,'readfromjabs')
  tfExpTagsJab = true;
elseif ~isempty(exptags)
  % It would be odd for a user to initially specified exptags, but not
  % initially specify rootexpdirs. Nonetheless, if this occurs and the user
  % select rootexpdirs (when prompted) that correspond to their exptags,
  % then everything will work and more power to them.
  ExperimentTags.verifyTags(exptags,numel(rootexpdirs));
  tfExpTagsJab = false;
else
  exptags = ExperimentTags.expTags(rootexpdirs);
  tfExpTagsJab = false;
end

% compile expdirs, expdirtags
expdirs = {};
expdirtags = {};
for ndx = 1:numel(rootexpdirs)
  dd = findAllSubDirs(rootexpdirs{ndx});
  nn = numel(dd);
  expdirs(end+1:end+nn) = dd;
  if ~tfExpTagsJab
    expdirtags(end+1:end+nn) = exptags(ndx);
  end
end
nexps = numel(expdirs);
if nexps == 0,
  fprintf('No experiments found.\n');
  T0 = nan;
  T1 = nan;
  return;
end

% compile expdirtags, read-from-jabfile case
if tfExpTagsJab
  % compile exp->tag map
  exp2tagmap = containers.Map();
  for iJab = 1:numel(jabfiles)
    Q = loadAnonymous(jabfiles{iJab});
    jabexpdirs = Q.expDirNames;
    jabexptags = Q.expDirTags;
    if isempty(jabexptags)
      % legacy Jab, no expDirTags
      jabexptags = ExperimentTags.expTags(jabexpdirs);
    else
      jabexptags = ExperimentTags.cleanLegacyNotTags(jabexptags);
    end
    for iExp = 1:numel(jabexpdirs)
      lowerexpname = lower(jabexpdirs{iExp});
      if exp2tagmap.isKey(lowerexpname)
        warning('ethogram_plot:expTags',...
          'Experiment ''%s'' is present in more than one Jabfile. Using tags from one Jab only.',...
          lowerexpname);
      end
      exp2tagmap(lowerexpname) = jabexptags{iExp};
    end
  end

  % look up tags for each element of expdir
  expdirtags = ExperimentTags.expTags(expdirs);
  for iExp = 1:nexps
    lowerexpname = lower(expdirs{iExp});
    if exp2tagmap.isKey(lowerexpname);
      expdirtags{iExp} = exp2tagmap(lowerexpname);
    else
      warning('ethogram_plot:expTags',...
          'Experiment ''%s'' not present in Jabfile(s). No tags will be associated.',...
          lowerexpname);
    end
  end
end

assert(numel(expdirtags)==nexps);

% Reclassify expdirs as appropriate
if ~dontreclassify && ~scoresnotjabs
  for iExp = 1:nexps
    for iJab = 1:numel(jabfiles)
      classifyMovie(expdirs{iExp},jabfiles{iJab},'doforce',doforce,'verbose',1);
    end
  end
end

[scorefns,behaviornames,jabsforscorefiles,trxfilename] = compileJabInfo(jabfiles,scoresnotjabs);
assert(isequal(numel(scorefns),numel(behaviornames),numel(jabsforscorefiles)));

if smartorder
  % order the scorefiles/behaviors/jabnames
  [behaviornames,idx] = orderbypats(behaviornames,ExpPP.BASICBEHAVIORS);
  scorefns = scorefns(idx);
  jabsforscorefiles = jabsforscorefiles(idx);
end

% some experiments may be missing scorefiles at this point
tfGoodExp = cellfun(@(x)ExpPP.isExpDir(x,scorefns),expdirs);
if ~all(tfGoodExp)
  tmpexp = expdirs(~tfGoodExp);
  tmpstr = sprintf('%s\n',tmpexp{:});
  warningNoTrace('ethogram_plot:expsDiscarded','Excluding %d experiments with missing scorefiles: %s',...
    nnz(~tfGoodExp),tmpstr);
end
expdirs = expdirs(tfGoodExp);
expdirtags = expdirtags(tfGoodExp);

[labels,T0,T1] = ...
  formLabelMatrix(expdirs,trxfilename,scorefns,behaviornames,nframesplot);

hfig = figure;

[tfbehspresent,missingbehs] = ExpPP.doNamesSpanAllBasicBehaviors(behaviornames);
reallydoautomarks = automarks && tfbehspresent;
if automarks && ~tfbehspresent
  % user said automarks, but basic behaviors missing
  warning('ethogram_plot:missingBehaviors',...
      'The following behaviors are missing: %s.',...
      civilizedStringFromCellArrayOfStrings(missingbehs));
end
if reallydoautomarks
  dauto = ExpPP.loadexps(expdirs,jabfiles,'includelabels',true);
  assert(numel(dauto)==numel(expdirs)); % ExpPP.loadexps does not guarantee this, but here all expdirs should be valid
  dauto = ExpPP.ensureCustomGroupInit(dauto);
  allUniqueTags = ExperimentTags.allUniqueTags(expdirtags);
  allUniqueTagsValid = matlab.lang.makeValidName(allUniqueTags);
  if numel(unique(allUniqueTagsValid))~=numel(allUniqueTagsValid)
    error('ethogram_plot:tag','Repeated tags.'); % very rare
  end
  Ntag = numel(allUniqueTags);
  for iTag = 1:Ntag
    tagRaw = allUniqueTags{iTag};
    cgrpname = allUniqueTagsValid{iTag};
    tfTag = ExperimentTags.findTag(expdirtags,tagRaw);
    dauto = ExpPP.createCustomGroup(dauto,tfTag,cgrpname);
  end
  ethplotcoreargs = {hfig,labels,T0,T1,expdirs,behaviornames,jabsforscorefiles,'doautomarks',true,'automarkdata',dauto,'scoresnotjabs',scoresnotjabs};
else
  ethplotcoreargs = {hfig,labels,T0,T1,expdirs,behaviornames,jabsforscorefiles,'scoresnotjabs',scoresnotjabs};
end
if ethcoreargexport
  fname = sprintf('ethcoreargs.%s.mat',datestr(now,'yyyymmddTHHMMSS'));
  save(fname,'ethplotcoreargs');
end
[boutmat,line_names] = ethogram_plot_core(ethplotcoreargs{:});

if xlsexport
  if isempty(xlsfile)
    [xlsfile,xlspath] = uiputfile('bouts.txt','Export to file');
    if isequal(xlsfile,0)
      xlsfile = [];
    else
      xlsfile = fullfile(xlspath,xlsfile);
    end
  elseif exist(xlsfile,'file')
    warning('ethogram_plot:fileExists','Overwriting file ''%s''.',xlsfile);
  end

  if ~isempty(xlsfile)
    Ethogram.boutExport(xlsfile,boutmat,line_names,behaviornames);
    fprintf(1,'Exported to tab-delimited file ''%s''.\n',xlsfile);
  end
end

function [scorefns,behaviornames,jabs,trxfilename] = compileJabInfo(jabfiles,scoresnotjabs)
% Compile scorefiles/behaviors from jabfiles.
%
% jabfiles: cellstr of jabfilenames, or cellstr of scorefilenames (archaic)
% if scoresnotjabs is true
%
% scorefns: cellstr of (short) score filenames
% behaviornames: cellstr of behaviornames, one per scorefn
% jabs: cellstr of jabfiles, one per scorefn. eg, jabs{i} is the jabfile
% containing scorefns{i} and behaviornames{i}. jabs may contain repeated
% entries when multi-classifier jabs are present.
% trxfilename: char, unique trxfilename in jabs (all jabs are expected to
% be consistent)

if scoresnotjabs
  scorefns = jabfiles;
  behaviornames = cell(size(scorefns));
  for i = 1:numel(scorefns)
    tok = regexp(scorefns{i},'scores_(.+).mat','tokens');
    assert(~isempty(tok),'Unable to parse score filename.');
    behaviornames{i} = tok{1}{1};
  end
  jabs = scorefns; % just use scorefilenames here; no jabs avail
  trxfilename = 'trx.mat';
else
  % jabfiles specified.
  Njab = numel(jabfiles);
  scorefns = cell(0,1);
  behaviornames = cell(0,1);
  jabs = cell(0,1);
  trxfilename = cell(Njab,1);
  for iJab = 1:Njab
    Q = loadAnonymous(jabfiles{iJab});

    sfns = Q.file.scorefilename;
    if ischar(sfns)
      sfns = {sfns};
    end
    behs = Labels.verifyBehaviorNames(Q.behaviors.names);
    assert(isequal(numel(sfns),numel(behs)),...
      'Inconsistent scores/behaviors in JAB file: %s',jabfiles{iJab});
    Nbeh = numel(behs);
    scorefns = cat(1,scorefns,sfns(:));
    behaviornames = cat(1,behaviornames,behs(:));
    jabs = cat(1,jabs,repmat(jabfiles(iJab),Nbeh,1));

    trxfilename{iJab} = Q.file.trxfilename;
  end

  trxfilename = unique(trxfilename);
  assert(isscalar(trxfilename),'Jabs contain inconsistent trxfilenames.');
  trxfilename = trxfilename{1};
end

function [labels,T0,T1] = formLabelMatrix(expdirs,trxfilename,scorefns,behaviornames,nframesplot)

% behaviornames just used to check against behaviornames in scorefiles

nexps = numel(expdirs);
nbehaviors = numel(behaviornames);
assert(numel(scorefns)==nbehaviors);
labels = false(nexps,0,nbehaviors);
tStarts = nan(1,nexps);
tEnds = nan(1,nexps);

for i = 1:nexps
  trxfile = fullfile(expdirs{i},trxfilename);
  load(trxfile,'trx');
  tStarts(i) = min([trx.firstframe]);
  tEnds(i) = max([trx.endframe]);

  % load score data
  allScores = cell(nbehaviors,1);
  for j = 1:nbehaviors,
    scorefile = fullfile(expdirs{i},scorefns{j});
    allScores{j} = load(scorefile);
    if isfield(allScores{j},'behaviorName') && ...
       ~strcmp(allScores{j}.behaviorName,behaviornames{j})
        warning('ethogram_plot:behaviorMismatch',...
          'Behavior label mismatch in experiment %s. Scorefile: %s. Expected: %s.',...
          expdirs{i},allScores{j}.behaviorName,behaviornames{j});
    end
    allScores{j} = allScores{j}.allScores;
    if iscell(allScores{j})
      assert(isscalar(allScores{j}));
      allScores{j} = allScores{j}{1};
    end
  end

  for j = 1:nbehaviors
    t0 = tStarts(i);
    t1 = min(tEnds(i),numel(allScores{j}.postprocessed{1}));
    if t1 < tEnds(i),
      fprintf('Number of frames for behavior idx%d = %d < number of frames in trajectory file %s = %d\n',...
        j,t1,trxfile,tEnds(i));
      tEnds(i) = t1;
    end

    if t1 > size(labels,2),
      labels = cat(2,labels,nan(nexps,t1-size(labels,2),nbehaviors));
    end

    labels(i,t0:t1,j) = allScores{j}.postprocessed{1}(t0:t1);
  end
end

T0 = min(tStarts);
T1 = max(tEnds)-1;

if ~isempty(nframesplot)
  if numel(nframesplot) == 1,
    T1 = min(T1,T0+nframesplot-1);
  else
    assert(numel(nframesplot)==2);
    T0 = nframesplot(1);
    T1 = min(T1,nframesplot(2));
  end
end
%nframesplot = T1-T0+1;

labels = labels(:,T0:T1,:);

function [cellstr,idx] = orderbypats(cellstr,pats)
% sort cellstr by pats (patterns) as best as possible

Npat = numel(pats);
loc = nan(size(cellstr));
for i = 1:Npat
  tf = regexpmatch(cellstr,pats{i},'caseinsens',true);
  loc(tf) = i;
end
loc(isnan(loc)) = inf; % elements of cellstr that didn't match go at the end

% loc(i) gives the desired relative position of cellstr{i}

[~,idx] = sort(loc);
cellstr = cellstr(idx);

author: kristinbranson

repo: JAABA


text:
%Run only this main program. No need to run the other functions separately
clear; clc;
DEPSK_mrc(1);  %performance of DEPSK by 1 Rx antenna
hold on;
clear; clc;
DEPSK_mrc(2);  %performance of DEPSK by 2 Rx antenna
hold on;
clear; clc;
BPSK_mrc(1);  %performance of BPSK by 1 Rx antenna
hold on;
clear; clc;
BPSK_mrc(2);  %performance of BPSK by 2 Rx antenna

author: Sable

repo: mcbench-benchmarks


text:
function [ypred, ll, mse] = eval_AR_perf(coef, C, y, model)
% Evaluate the performance of an AR model.
%
% Inputs
% coef(:,:,k,m) - coef. matrix to use for k steps back, model m
% C(:,:,m)      - cov. matrix for model m
% y(:,t)        - observation at time t
% model(t)      - which model to use at time t (defaults to 1 if not specified)
%
% Outputs
% ypred(:,t)    - the predicted value of y at t based on the evidence thru t-1.
% ll            - log likelihood
% mse           - mean squared error = sum_t d_t . d_t, where d_t = pred(y_t) - y(t)

[s T] = size(y);
k = size(coef, 3);
M = size(coef, 4);

if nargin<4, model = ones(1, T); end

ypred = zeros(s, T);
ypred(:, 1:k) = y(:, 1:k);
mse = 0;
ll = 0;
for j=1:M
  c(j) = log(normal_coef(C(:,:,j)));
  invC(:,:,j) = inv(C(:,:,j));
end
coef = reshape(coef, [s s*k M]);

for t=k+1:T
  m = model(t-k);
  past = y(:,t-1:-1:t-k);
  ypred(:,t) = coef(:, :, m) * past(:);
  d = ypred(:,t) - y(:,t);
  mse = mse + d' * d;
  ll = ll + c(m) - 0.5*(d' * invC(:,:,m) * d);
end
mse = mse / (T-k+1);


text:
function img = imDrawLine(img, pos1, pos2, varargin)
%IMDRAWLINE Draw a line between two points in the image
%
%   IMG2 = imDrawLine(IMG, P1, P2);
%   IMG2 = imDrawLine(IMG, X1, Y1, X2, Y2);
%
%   [X, Y] = imDrawLine(X1, X2, Y1, Y2) computes an approximation to the
%   line segment joining (X1, Y1) and (X2, Y2) with integer coordinates.
%   X1, X2, Y1, and Y2 should be integers. imDrawLine is reversible; that
%   is, imDrawLine(X1, X2, Y1, Y2) produces the same results as
%   FLIPUD(imDrawLine(X2, X1, Y2, Y1)).
%
%   Example
%   % Draw some lines to make a diamond
%     img = imCreate([16 16], 'uint8');
%     img = imDrawLine(img, [8 1], [16 8]);
%     img = imDrawLine(img, [16 8], [8 16]);
%     img = imDrawLine(img, [8 16], [1 8]);
%     img = imDrawLine(img, [1 8], [8 1]);
%     image(img);
%
%   % Overlay some lines on a grayscale image
%     % read image
%     img = imread('cameraman.tif');
%     % draw white line
%     img = imDrawLine(img, [10 30], [150 110], 255);
%     % draw also a black line
%     img = imDrawLine(img, [10 30], [100 210], 0);
%     % display result
%     imshow(img);
%
%   See also
%     bresenhamLine

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-11-25,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Input argument extraction

% process case of input given as imDrawLine(IMG, X1, Y1, X2, Y2)
if isscalar(pos1)
    if nargin < 5
        error('Please give coordinates either as 2 points, or as 4 coords');
    end
    pos1 = [pos1 pos2];
    pos2 = [varargin{1} varargin{2}];
    varargin(1:2) = [];
end

% extract color
color = 255;
if ~isempty(varargin)
    color = varargin{1};
end

%% Coordinates computation

% coordinates of pixels
[x, y] = bresenhamLine(pos1, pos2, varargin{:});

% remove line pixels outside image
xOut = x < 1 | x > size(img, 2);
yOut = y < 1 | y > size(img, 1);
x(xOut | yOut) = [];
y(xOut | yOut) = [];

%% Write the line into image

if size(img, 3) == 1
    % grayscale image
    for i = 1:length(x)
        img(y(i), x(i)) = color;
    end

else
    % color image
    for i = 1:length(x)
        img(y(i), x(i), :) = color(:);
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage