text:
function sr_parameterstudyNvsTw

disp('~/zmap/src/thomas/seismicrates/sr_parameterstudyNvsTw.m')

% reset random number generator
rand('state',sum(100*clock));
% set default sample size

fYrs=[10]';
fYr1=1980;
fYr2=1980+fYrs;

% create range for N
vN=[25 50 100 150 200]';
% vN=[50 75 100]';

% create range for Tw
vTw=[0.1:0.1:0.9]'.*fYrs;

% fix Tbin by 1/100 of whole period
nTbin1=1/100;
% round time bins
vTbin=ceil((fYr2-fYr1)*nTbin1*365)';

% set default catalog size
nCatSize=5000;

% prepare input parameter matrix  [ vN   vTw   vTbin ]
mVar(:,1)=repmat(vN,size(vTw,1)*size(vTbin,1),1);
mVar(:,2)=reshape(repmat(vTw',size(vN,1)*size(vTbin,1),1),...
    size(vTw,1)*size(vN,1)*size(vTbin,1),1);
mVar(:,3)=repmat(reshape(repmat(vTbin',size(vN,1),1),...
    size(vTbin,1)*size(vN,1),1),size(vTw,1),1);

mCat20=[];mCat00=[];

% bRchange=logical(0)
% rate change?
fRate2=100  ; %  percent of original rate
% nInit=logical(1);
for iCat=1:500
    iCat
    for j=1:size(mVar)
        while isempty(mCat20) | isempty(mCat00)
            vRand=rand(mVar(j,1),1);

            if (fRate2~=100)
                % Fraction of period with 100% rate and fRchange % rate
                fRYr1=1-mVar(j,2)/(fYr2-fYr1);
                fRYr2=mVar(j,2)/(fYr2-fYr1);
                vRand=rand(mVar(j,1),1);
                fRYr1+(fRYr2*fRate2/100);
                vRand1=vRand*(fRYr1+(fRYr2*fRate2/100));
                vSel=(vRand1>fRYr1);
                vRand=vRand1;
                vRand(vSel)=(vRand1(vSel)-fRYr1)*100/fRate2+fRYr1;
                clear vRand1
            else
                vRand=rand(mVar(j,1),1);
            end

            mCat00=vRand*(fYr2-fYr1)+fYr1;
            mCat=rand(nCatSize,1)*(fYr2-fYr1)+fYr1;
            mCat20=mCat00(mCat00>(fYr2-mVar(j,2)));
            if isempty(mCat20) disp('looooop'); end
        end
%         [mLTA(j,iCat), mLTAprob(j,iCat)] =calc_zlta(mCat,mCat00,mCat20,fYr1, fYr2,mVar(j,2),mVar(j,3),mVar(j,1));
%         [mBeta(j,iCat), mBetaprob(j,iCat)] =calc_beta(mCat,mCat00,mCat20,fYr1, fYr2,mVar(j,2),mVar(j,3),mVar(j,1));
        [mResult_(iCat,1,j), mResult_(iCat,2,j)] =calc_zlta(mCat,mCat00,mCat20,fYr1, fYr2,mVar(j,2),mVar(j,3),mVar(j,1));
        [mResult_(iCat,3,j), mResult_(iCat,4,j)] =calc_beta(mCat,mCat00,mCat20,fYr1, fYr2,mVar(j,2),mVar(j,3),mVar(j,1));
        mCat20=[];
        clear mCat00 mCat;
    end
end
%         disp(n);
%     subplot(2,2,1)
%     hold on;plot(1./nTbin1,mLTA);
%     subplot(2,2,2)
%     hold on;plot(1./nTbin1,mLTAprob);
%     subplot(2,2,3)
%     hold on;plot(1./nTbin1,mBeta);
%     subplot(2,2,4)
%     hold on;plot(1./nTbin1,mBetaprob);

% for i=1:size(mLTA,1)
%     mLTA1(:,i)=reshape(mLTA(i,:,:),size(mLTA,2)*size(mLTA,3),1);
%     mLTAprob1(:,i)=reshape(mLTAprob(i,:,:),size(mLTAprob,2)*size(mLTAprob,3),1);
%     mBeta1(:,i)=reshape(mBeta(i,:,:),size(mBeta,2)*size(mBeta,3),1);
%     mBetaprob1(:,i)=reshape(mBetaprob(i,:,:),size(mBetaprob,2)*size(mBetaprob,3),1);
% end

params.mResult_=mResult_;
params.mVar=mVar;
params.vN=vN;
params.vTw=vTw;
params.vTbin=vTbin;

sPrint=sprintf('save NvsTw-Sim%04.0f-R%03.0f.mat params -mat',iCat, fRate2);
eval(sPrint);
disp(sPrint);
% figure;
% subplot(2,1,1)
% errorbar(1./nTbin1,mean(mLTA1),std(mLTA1),'b--x','LineWidth',3);
% hold on;
% errorbar(1./nTbin1,mean(mBeta1),std(mBeta1),'r:o','LineWidth',2);
% xlabel('Bin Size [1/X]');
% ylabel('Z and \beta [ ]');
% legend('z','\beta');
% subplot(2,1,2);
% errorbar(1./nTbin1,mean(mLTAprob1),std(mLTAprob1),'b--x','LineWidth',3);
% YLim([0,1]);
% hold on;
% errorbar(1./nTbin1,mean(mBetaprob1),std(mBetaprob1),'r:o','LineWidth',2);
% xlabel('Bin Size [1/X]');
% ylabel('Probability of Z and \beta [ ]');
% legend('p(z)','p(\beta)');
%

author: CelsoReyes


text:
function [theta,alpha,e]=tls(Z)
% tls
% Computes tls estimate of hyperplane parameters
% [theta, alpha] = tls ( Z)
% Z is of form (p1'; p2'; ..) where pi is p-by-1 (dimension p)
% Computed hyperplane is of form theta' * p = alpha with alpha>0

%* Author: Ranjith Unnikrishnan                                          *
%* Carnegie Mellon University, Vision and Mobile Robotics Laboratory       *
%* THE MATERIAL EMBODIED IN THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"     *
%* AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,      *
%* INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR      *
%* FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL CARNEGIE MELLON  *
%* UNIVERSITY BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,            *
%* SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY         *
%* KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,        *
%* LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF     *
%* THIRD PARTIES, WHETHER OR NOT CARNEGIE MELLON UNIVERSITY HAS BEEN     *
%* ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON        *
%* ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE     *
%* POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.                      *
%

p=size(Z,2);
n=size(Z,1);
centroid=mean(Z,1);
Z=Z-repmat(centroid,n,1);
[U,S,V]=svd(Z'*Z);
theta=U(:,p);
alpha=theta'*centroid';

if(alpha<0)
    theta=-theta;
    alpha=-alpha;
end

e=median( abs(Z*theta - alpha +centroid*theta));
return;

author: zhixy

repo: Laser-Camera-Calibration-Toolbox


text:
function [roiVertInds, roiMapMode] = findROIVertices(roi, nodeInds, v2g)
% Returns a list of indices to each vertex belonging to the current ROI.
% [roiVertInds, roiMapMode] = findROIVertices(~, nodeInds, v2g)

% get ROI map mode
prefs = mrmPreferences;
if isequal(prefs.layerMapMode, 'layer1')
    roiMapMode = 'layer1';
else
    roiMapMode = 'any';
end

switch lower(roiMapMode)
	case 'layer1'
		% The following will give us *a* vertex index for each gray node. However,
		% the same gray node may map to multiple vertices, and we want them
		% all. (Otherwise, the color overlay will have holes.) So, we loop
		% until we've got them all.
		[junk, roiVertInds] = ismember(nodeInds, v2g(1,:));
		roiVertInds = roiVertInds(roiVertInds>0);
		while(any(junk))
			v2g(1,roiVertInds) = 0;
			[junk,tmp] = ismember(nodeInds, v2g(1,:));
			roiVertInds = [roiVertInds; tmp(tmp>0)];
		end

	case 'any'
		% if any of the nodes mapping to a given vertex are in the
		% ROI, include that vertex for drawing the ROI
		I = ismember(v2g, nodeInds);
		roiVertInds = find(sum(I)>0); % find columns w/ at least 1 member

	case 'data'
		% take ROI value from the same nodes as the data mapping,
		% rounding up (e.g., for 'mean' data mapping, will behave
		% like 'any'

	otherwise
		error('Invalid ROI Draw Mode preference.')
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function D = DLSI_updateD(D, E, F, A, lambda, opts)
% function D = DLSI_updateD(D, E, F, A, lambda, opts)
% problem: `D = argmin_D -2trace(ED') + trace(FD'*D) + lambda *||A*D||F^2,`
% subject to: `||d_i||_2^2 <= 1`
% where F is a positive semidefinite matrix
% ========= aproach: ADMM ==============================
% rewrite: `[D, Z] = argmin -2trace(ED') + trace(FD'*D) + lambda ||A*Z||_F^2,`
%     subject to `D = Z; ||d_i||_2^2 <= 1`
% aproach 1: ADMM.
% 1. D = -2trace(ED') + trace(FD'*D) + rho/2 ||D - Z + U||_F^2,
%     s.t. ||d_i||_2^2 <= 1
% 2. Z = argmin lambda*||A*Z|| + rho/2||D - Z + U||_F^2
% 3. U = U + D - Z
% solve 1: D = argmin -2trace(ED') + trace(FD'*D) + rho/2 ||D - W||_F^2
%                       with W = Z - U;
%            = argmin -2trace((E - rho/2*W)*D') +
%               trace((F + rho/2 * eye())*D'D)
% solve 2: derivetaive: 0 = 2A'AZ + rho (Z - V) with V = D + U
% `Z = B*rhoV` with `B = (2*lambda*A'*A + rho I)^{-1}`
% `U = U + D - Z`
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------

    if nargin == 0
        clc;
        d = 300;
        N = 10;
        k = 5;
        k2 = 495;
        load('tmp.mat');
        lambda = 0.01;

        opts.show = 0;
        opts.max_iter = 300;
        opts.verbose = 1;

    end
    if nargin == 6
        opts.lambda = lambda;
    elseif nargin == 5
        opts = lambda;
    end
    %%
    function cost = calcost(D)
        cost =  -2*trace(E*D') + trace(F*D'*D) + lambda*normF2(A*D);
    end
    %%
    iter = 0;
    rho = 1.0;
    Z_old = D;
    U = zeros(size(D));
    I_k = eye(size(D,2));
    % B = inv(2*lambda*A'*A + rho*I_k2); However, this might be very expensive if size(A, 2) is big, which is common
    % Instead, we can use the Sherman–Morrison formula at
    % https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization_(Woodbury_Matrix_Identity)
    X = 2*lambda/rho*A';
    Y = A;
    B1 = X*inv(eye(size(Y, 1)) + Y*X);
    tol = 1e-5;
    optsD.max_iter = 100;
    optsD.tol = 1e-8;
    while iter < opts.max_iter
        iter = iter + 1;
        %% ========= update D ==============================
        W = Z_old - U;
        E2 = E + rho/2 * W;
        F2 = F + rho/2*I_k;
        D = ODL_updateD(D, E2, F2, optsD);
        %% ========= update Z ==============================
        V = D + U;
        % Z_new = rho*B*V; slow
        Z_new = rho*(V - B1*(Y*V)); % fast
        e1 = normF2(D - Z_new);
        e2 = rho*normF2(Z_new - Z_old);
        if (e1 < optsD.tol && e2 < optsD.tol)
            break;
        end
        if opts.verbose
            cost = calcost(D);
            fprintf('iter = %3d | costD = %5.4f | normF2(D - Z) = %5.4f | rho(Z_new - Z_old) = %5.4f\n', iter, cost, e1, e2);
        end
        %% ========= update U ==============================
        U = U + D - Z_new;
        Z_old = Z_new;
    end
%     disp(t1)
%     disp(t2)
    if nargin == 0
        D = [];
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function M=link_communities(W,type_clustering)
%LINK_COMMUNITIES     Optimal overlapping community structure
%
%   M = link_communities(W)
%   M = link_communities(W,'complete');
%
%   The optimal community structure is a subdivision of the network into
%   groups of nodes which have a high number of within-group connections
%   and a low number of between group connections.
%
%   This algorithm uncovers overlapping community structure via
%   hierarchical clustering of network links. This algorith is generalized
%   for weighted/directed/fully-connected networks.
%
%   Input:      W,                  directed (weighted or binary) connection matrix.
%               type_clustering,    type of hierarchical clustering (optional)
%                                       'single'        single-linkage (default)
%                                       'complete'      complete-linkage
%
%   Output:     M,                  nodal community-affiliation matrix
%                                   binary matrix of size CxN [communities x nodes]
%
%   NB: The algorithm can be slow and memory intensive.
%
%   Reference: Ahn, Bagrow and Lehmann (2010) Nature 466, 761–764.
%
%   Mika Rubinov, U Cambridge, 2014-2015

%% initialize

n=size(W,1);                                                        % number of nodes
W(1:n+1:end)=0;
W=W./max(W(:));                                                     % normalize weights

if ~exist('type_clustering','var')
    type_clustering='single';
end

%% get node similarity

W(1:n+1:end) = ( sum(W)/sum(W~=0) + sum(W.')/sum(W.'~=0) )/2;       % mean weight on diagonal
No=sum(W.^2,2);                                                     % out-norm squared
Ni=sum(W.^2,1);                                                     % in-norm squared

Jo=zeros(n);                                                        % weighted in-Jaccard
Ji=zeros(n);                                                        % weighted ou-Jaccard
for b=1:n
    for c=1:n
        Do=W(b,:)*W(c,:).';
        Jo(b,c)=Do./(No(b)+No(c)-Do);

        Di=W(:,b).'*W(:,c);
        Ji(b,c)=Di./(Ni(b)+Ni(c)-Di);
    end
end

%% get link similarity

[A,B]=find( (W|W.') & triu(ones(n),1));
m=length(A);
Ln=zeros(m,2);                                                      % link nodes
Lw=zeros(m,1);                                                      % link weights
for i=1:m
    Ln(i,:) = [A(i) B(i)];                                          % link nodes
    Lw(i) = (W(A(i),B(i))+W(B(i),A(i)))/2;                          % link weight
end

ES=zeros(m,m,'single');                                             % link similarity
for i=1:m
    for j=1:m
        if      Ln(i,1)==Ln(j,1); a=Ln(i,1); b=Ln(i,2); c=Ln(j,2);
        elseif  Ln(i,1)==Ln(j,2); a=Ln(i,1); b=Ln(i,2); c=Ln(j,1);
        elseif  Ln(i,2)==Ln(j,1); a=Ln(i,2); b=Ln(i,1); c=Ln(j,2);
        elseif  Ln(i,2)==Ln(j,2); a=Ln(i,2); b=Ln(i,1); c=Ln(j,1);
        else    continue
        end

        ES(i,j) = (W(a,b)*W(a,c)*Ji(b,c) + W(b,a)*W(c,a)*Jo(b,c))/2;
    end
end
ES(1:m+1:end)=0;

%% perform hierarchical clustering

C=zeros(m,m,'single');                                              % community affiliation matrix
Nc=C; Mc=C; Dc=C;                                                   % communities nodes, links and density
U=1:m;                                                              % initial community assignments
C(1,:)=U;                                                           % as above, in the matrix

for i=1:m-1; fprintf('hierarchy%8d\n',i)                            % hierarchy level

    % compute densities
    for j=1:length(U)                                               % loop over communities
        idx = C(i,:)==U(j);                                         % get link indices
        links = sort(Lw(idx));                                      % sort link weights
        nodes = sort(reshape(Ln(idx,:),2*nnz(idx),1));
        nodes = nodes([true;nodes(2:end)~=nodes(1:end-1)]);         % get unique nodes

        nc = numel(nodes);                                          % community nodes
        mc = sum(links);                                            % community weights
        min_mc = sum(links(1:nc-1));                                % minimal weight
        dc = (mc - min_mc) / (nc.*(nc-1)/2 - min_mc);               % community density

        Nc(i,j)=nc;
        Mc(i,j)=mc;
        Dc(i,j)=dc;
    end

    % cluster
    C(i+1,:)=C(i,:);                                                % copy current partition
    [u1,u2]=find(ES(U,U)==max(max(ES(U,U))));                       % on this line MAXs MUST BE MAXs

    V=U(unique(sortrows(sort([u1 u2],2)),'rows'));                  % get unique links
    for j=1:size(V,1)
        switch type_clustering
            case 'single';      x = max(ES(V(j,:),:),[],1);         % max -> single linkage
            case 'complete';    x = min(ES(V(j,:),:),[],1);         % min -> complete linkage
            otherwise; error('Unknown clustering type.');
        end
        ES(V(j,:),:) = [x;x];                                       % assign distances to whole clusters
        ES(:,V(j,:)) = [x;x].';
        ES(V(j,1),V(j,1)) = 0;                                      % clear diagonal
        ES(V(j,2),V(j,2)) = 0;                                      % clear diagonal

        C(i+1,C(i+1,:)==V(j,2)) = V(j,1);                           % merge communities
        V(V==V(j,2)) = V(j,1);                                      % merge indices
    end

    U=unique(C(i+1,:));                                             % get unique communities
    if numel(U)==1
        break;
    end
end

%%

Dc(isnan(Dc))=0;
[~,i]=max(sum(Dc.*Mc,2));                                           % get maximal density

U=unique(C(i,:));                                                   % unique communities
M=zeros(1,n);                                                       % nodal affiliations
for j=1:length(U)
    M(j,unique( Ln(C(i,:)==U(j),:)) )=1;
end
M=M(sum(M,2)>2,:);

% M2=zeros(n);                                                      % two dimensional nodal affiliation
% for i=1:size(M,1);
%     M2=M2+(M(i,:).'*ones(1,n) & ones(n,1)*M(i,:));
% end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [imsource, D, C, F] = dfviewexcludepreview(outlier, width, height, dsname)
% For use by DFITTOOL

%   $Revision: 1.1.6.5 $
%   Copyright 2003-2004 The MathWorks, Inc.

tempfigure=figure('units','pixels','position',[0 0 width height], ...
      'handlevisibility','off', ...
      'integerhandle','off', ...
      'visible','off', ...
      'paperpositionmode', 'auto', ...
      'color','w');

NONE='(none)';

% We're excluding based on data in one dataset
ds = find(getdsdb,'name',dsname);
if isempty(ds)
   return
end

% Get data w/o NaNs but with no exclusion rule applied
[ydata, cens, freq] = getincludeddata(ds,[]);
ydata = real(ydata);
if isempty(cens)
   cens = zeros(size(ydata));
end
if isempty(freq)
   freq = ones(size(ydata));
end

% Sort y and carry along the rest
[ydata,i] = sort(ydata);
cens = cens(i);
freq = freq(i);

% Create x and y vectors to plot
n = sum(freq);
x = zeros(n,1);
y = zeros(n,1);
g = zeros(n,1);
j = 1;
x(1:freq(1)) = ydata(1);
y(1:freq(1)) = (1:freq(1))';
g(1:freq(1)) = cens(1);
i = freq(1)+1;
for k=2:length(ydata)
   for j=1:freq(k)
      x(i) = ydata(k);
      g(i) = cens(k);
      if (i>1) && (x(i)==x(i-1))
         y(i) = y(i-1) + 1;
      else
         y(i) = 1;
      end
      i = i+1;
   end
end

o = handle(outlier);
ylo = o.YLow;
if isempty(ylo)
	ylo = -Inf;
else
	ylo = str2double(ylo);
end

yhi = o.YHigh;
if isempty(yhi)
	yhi = Inf;
else
	yhi = str2double(yhi);
end

ylotest = o.YLowLessEqual;
yhitest = o.YLowLessEqual;

xlim = [min(x) max(x)];
xlim = xlim + .05 * [-1 1] * diff(xlim);
ylim = [min(y) max(y)];
ylim = ylim + .05 * [-1 1] * diff(ylim);
if ylim(1) == ylim(2)
    ylim = [0 2];
end
ax=axes('position',[.05 .05 .9 .9], ...
      'parent',tempfigure, ...
      'xtick',[],'ytick',[], ...
      'box','on', ...
      'visible','off', 'XLim',xlim,'YLim',ylim);

if ylotest==0
   inbounds = x>=ylo;
else
   inbounds = x>ylo;
end
if yhitest==0
   inbounds = inbounds & x<=yhi;
else
   inbounds = inbounds & x<yhi;
end

figcolor = get(0,'defaultuicontrolbackgroundcolor');

t = inbounds;
l1 = line('XData',x(t),'YData',y(t),...
          'Color','b','Marker','.','LineStyle','none',...
          'Parent',ax);
t = ~inbounds;
l2 = line('XData',x(t),'YData',y(t),...
          'Color',figcolor/2,'Marker','.','LineStyle','none',...
          'Parent',ax);
alllines = [l1 l2];

gr = [.9 .9 .9];

% Create patches to show the area outside the domain and range

allpatches = [];

if (ylo ~= -Inf)
    xlo = max(ylo,xlim(1));
    xp = [xlim(1) xlo xlo xlim(1)];
    yp = [ylim(1) ylim(1) ylim(2) ylim(2)];
    p1=patch(xp,yp,gr,'LineStyle','none','Parent',ax);
    allpatches(end+1)= p1;
end

if (yhi ~= Inf)
    xhi = min(yhi,xlim(2));
    xp = [xlim(2) xhi xhi xlim(2)];
    yp = [ylim(1) ylim(1) ylim(2) ylim(2)];
    p2=patch(xp,yp,gr,'LineStyle','none','Parent',ax);
    allpatches(end+1)= p2;
end

set(ax,'Child',[alllines allpatches]);

x=hardcopy(tempfigure,'-dzbuffer','-r0');
% give the image a black edge
x(1,:,:)=0; x(end,:,:)=0; x(:,1,:)=0; x(:,end,:)=0;
imsource=im2mis(x);
delete(tempfigure);
[D, C, F] = dfviewdata(ds);

author: zouchuhang

repo: LayoutNet


text:
function ng = tetrahedron_grid_count ( n )

%*****************************************************************************80
%
%% TETRAHEDRON_GRID_COUNT counts the grid points inside a tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Output, integer NG, the number of grid points.
%
  ng = ( ( n + 1 ) * ( n + 2 ) * ( n + 3 ) ) / 6;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,VAR]=least_square(v,H,P)

%normal vector
N = (H'*P*H);

%cumpute unknown parameter
x = (N^-1)*H'*P*v;

%convariance of parameter
VAR = N^-1;

return

repo: GINav


text:
function [ a, b ] = p00_lim ( problem )

%*****************************************************************************80
%
%% P00_LIM returns the integration limits for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 September 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Output, real A(2), B(2), the lower and upper limits of integration.
%
  if ( problem == 1 )
    [ a, b ] = p01_lim ( );
  elseif ( problem == 2 )
    [ a, b ] = p02_lim ( );
  elseif ( problem == 3 )
    [ a, b ] = p03_lim ( );
  elseif ( problem == 4 )
    [ a, b ] = p04_lim ( );
  elseif ( problem == 5 )
    [ a, b ] = p05_lim ( );
  elseif ( problem == 6 )
    [ a, b ] = p06_lim ( );
  elseif ( problem == 7 )
    [ a, b ] = p07_lim ( );
  elseif ( problem == 8 )
    [ a, b ] = p08_lim ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_LIM - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem index = %d\n', problem );
    error ( 'P00_LIM - Fatal error!\n' );
  end

  return
end


text:
function [ hyps ] = generateHypsB( lines, vp, tol, omap, gc )
%GENERATEHYPS generate hypothesis, screen out inconsistent ones according
%to OM, GC, and merged map
%   lines: line segment
%   vp: vanishing point, in z,y,x order
%   tol: tolerance when assigning lines to vp
%   omap & gc: evidence
%   hyps: all hyps with >0.5 consistency with OM will be returned;

%% >> ASSIGN LINE SEGMENTS WITH TYPE

%% get line segments align with vanishing point
lines1 = getVanishingLine(lines, vp(1,:), tol);
lines2 = getVanishingLine(lines, vp(2,:), tol);
lines3 = getVanishingLine(lines, vp(3,:), tol);

%% extend line representation with line type
% line: [1x3 n] [planeID] [1x2 uv1] [u/d f/b r/l] (1/-1, 0 not sure, 2 not available)
lines1_attr = zeros(size(lines1,1),9);
lines2_attr = zeros(size(lines2,1),9);
lines3_attr = zeros(size(lines3,1),9);
lines1_attr(:,1:6) = lines1(:,1:6);
lines2_attr(:,1:6) = lines2(:,1:6);
lines3_attr(:,1:6) = lines3(:,1:6);

%% horizontal line, line connecting vp2 and vp3, u=1, d=-1
horizontal = cross(vp(2,:),vp(3,:));
horizontal = horizontal./norm(horizontal);
if horizontal(:,3)<0
    horizontal = -horizontal;
end
% line belong to vp2 and vp3 should have this attribute
% lines1 not available
lines1_attr(:,7) = 2;
% lines2
for i = 1:size(lines2_attr,1)
    l = lines2(i,:);
    u = l(5:6)*2*pi-pi;
    v = computeUVN( l(1:3), u, l(4));
    xyz = uv2xyzN([u' v'], l(4));
    b1 = dot(xyz(1,:), horizontal)>0;
    b2 = dot(xyz(2,:), horizontal)>0;
    if b1&&b2
        lines2_attr(i,7) = 1;
    elseif ~b1&&~b2
        lines2_attr(i,7) = -1;
    else
        lines2_attr(i,7) = 0;
    end
end
% lines3
for i = 1:size(lines3_attr,1)
    l = lines3(i,:);
    u = l(5:6)*2*pi-pi;
    v = computeUVN( l(1:3), u, l(4));
    xyz = uv2xyzN([u' v'], l(4));
    b1 = dot(xyz(1,:), horizontal)>0;
    b2 = dot(xyz(2,:), horizontal)>0;
    if b1&&b2
        lines3_attr(i,7) = 1;
    elseif ~b1&&~b2
        lines3_attr(i,7) = -1;
    else
        lines3_attr(i,7) = 0;
    end
end

%% front and back, line connecting vp1 and vp3, f=1, b=-1
frontNback = cross(vp(1,:),vp(3,:));
frontNback = frontNback./norm(frontNback);
if frontNback(:,2)<0
    frontNback = -frontNback;
end
% line belong to vp1 and vp3 should have this attribute
% lines2 not available
lines2_attr(:,8) = 2;
% lines1
for i = 1:size(lines1_attr,1)
    l = lines1(i,:);
    u = l(5:6)*2*pi-pi;
    v = computeUVN( l(1:3), u, l(4));
    xyz = uv2xyzN( [u' v'], l(4));
    b1 = dot(xyz(1,:), frontNback)>0;
    b2 = dot(xyz(2,:), frontNback)>0;
    if b1&&b2
        lines1_attr(i,8) = 1;
    elseif ~b1&&~b2
        lines1_attr(i,8) = -1;
    else
        lines1_attr(i,8) = 0;
    end
end
% lines3
for i = 1:size(lines3_attr,1)
    l = lines3(i,:);
    u = l(5:6)*2*pi-pi;
    v = computeUVN( l(1:3), u, l(4));
    xyz = uv2xyzN( [u' v'], l(4));
    b1 = dot(xyz(1,:), frontNback)>0;
    b2 = dot(xyz(2,:), frontNback)>0;
    if b1&&b2
        lines3_attr(i,8) = 1;
    elseif ~b1&&~b2
        lines3_attr(i,8) = -1;
    else
        lines3_attr(i,8) = 0;
    end
end

%% left and right, line connecting vp1 and vp2, r=1, l=-1
leftNright = cross(vp(1,:),vp(2,:));
leftNright = leftNright./norm(leftNright);
if leftNright(:,1)<0
    leftNright = -leftNright;
end
% line belong to vp1 and vp2 should have this attribute
% lines3 not available
lines3_attr(:,9) = 2;
% lines1
for i = 1:size(lines1_attr,1)
    l = lines1(i,:);
    u = l(5:6)*2*pi-pi;
    v = computeUVN( l(1:3), u, l(4));
    xyz = uv2xyzN( [u' v'], l(4));
    b1 = dot(xyz(1,:), leftNright)>0;
    b2 = dot(xyz(2,:), leftNright)>0;
    if b1&&b2
        lines1_attr(i,9) = 1;
    elseif ~b1&&~b2
        lines1_attr(i,9) = -1;
    else
        lines1_attr(i,9) = 0;
    end
end
% lines2
for i = 1:size(lines2_attr,1)
    l = lines2(i,:);
    u = l(5:6)*2*pi-pi;
    v = computeUVN( l(1:3), u, l(4));
    xyz = uv2xyzN( [u' v'], l(4));
    b1 = dot(xyz(1,:), leftNright)>0;
    b2 = dot(xyz(2,:), leftNright)>0;
    if b1&&b2
        lines2_attr(i,9) = 1;
    elseif ~b1&&~b2
        lines2_attr(i,9) = -1;
    else
        lines2_attr(i,9) = 0;
    end
end

%% save line segments according to type
typ = [1 1 2; 1 2 -1; 1 -1 2; 1 2 1; ...
       -1 1 2; -1 2 -1; -1 -1 2; -1 2 1; ...
       2 1 -1; 2 1 1; 2 -1 1; 2 -1 -1];
allLines = [lines1_attr; lines2_attr; lines3_attr];
typLine = repmat(struct('lines', [], 'lsLength', [], 'lspool', [], ...
                        'poolLength', []), [size(typ,1) 1]);
for i = 1:size(typ,1)
    valid = (allLines(:,7)==typ(i,1)) & (allLines(:,8)==typ(i,2)) & (allLines(:,9)==typ(i,3));
    typLine(i).lines = allLines(valid,:);
    lsLength = typLine(i).lines(:,6)-typLine(i).lines(:,5);
    lsLength(lsLength<0) = lsLength(lsLength<0) + 1;
    sumLength = sum(lsLength);
    lspool = zeros( ceil(sumLength*100), 1);
    accumuLen = 0;
    poolLength = length(lspool);
    for j = 1:length(lsLength)
        startid = round(accumuLen/sumLength*poolLength) + 1;
        accumuLen = accumuLen + lsLength(j);
        endid = round(accumuLen/sumLength*poolLength);
        lspool(startid:endid) = j;
    end

    typLine(i).lsLength = lsLength;
    typLine(i).lspool = lspool;
    typLine(i).poolLength = poolLength;
end

%% >> BUILD UP RECIPE FROM LINE TYPE TO BOX
% +1: max coords; -1: min coords;
% uf [ 0 +1 +1]  df [ 0 +1 -1]  lf [-1 +1  0];
% ul [-1  0 +1]  dl [-1  0 -1]  rf [+1 +1  0];
% ub [ 0 -1 +1]  db [ 0 -1 -1]  rb [+1 -1  0];
% ur [+1  0 +1]  dr [+1  0 -1]  lb [-1 -1  0];

%% generate recipe: from line type to box by check rank of matrix
allFunc = [ 0 0 0 1 0 -1; ...
            1 0 0 0 0 -1; ...
            0 0 1 0 0 -1; ...
            0 1 0 0 0 -1; ...
            0 0 0 1 -1 0; ...
            1 0 0 0 -1 0; ...
            0 0 1 0 -1 0; ...
            0 1 0 0 -1 0; ...
            1 0 0 -1 0 0; ...
            0 1 0 -1 0 0; ...
            0 1 -1 0 0 0; ...
            1 0 -1 0 0 0];
% combos = combntns(1:size(allFunc,1), 5);
combos = nchoosek(1:size(allFunc,1), 5);
degen = false(size(combos,1),1);
for i = 1:size(combos,1)
    equset = allFunc(combos(i,:), :);
    if rank(equset)<5
        degen(i) = true;
    end
end
cube_recipe = combos(~degen,:);

%% define geometry of cuboid: correlation among corner, edge, and face
% 8 corners, each is the intersection of 3 edges from 3 vps respectively
edgeToPoint = [ 9  2  1; ...
               12  2  3; ...
               11  4  3; ...
               10  4  1; ...
                9  6  5; ...
               12  6  7; ...
               11  8  7; ...
               10  8  5];
% 8 corners, each is assigned with type [u(+1)/d(-1) f(+1)/b(-1) r(+1)/l(-1)]
pointDirect = [ 1  1 -1; ...
                1 -1 -1; ...
                1 -1  1; ...
                1  1  1; ...
               -1  1 -1; ...
               -1 -1 -1; ...
               -1 -1  1; ...
               -1  1  1];
% 12 edges, each connects two corners
PointToEdge = [ 1  4; ...
                2  1; ...
                2  3; ...
                3  4; ...
                5  8; ...
                6  5; ...
                6  7; ...
                7  8; ...
                5  1; ...
                8  4; ...
                7  3; ...
                6  2];
% 6 faces, each face has 4 corners
% ranked as u d f b r l; start from left bottom, clockwise
PointToSurface = [1 2 3 4; ...
                  6 5 8 7; ...
                  5 1 4 8; ...
                  7 3 2 6; ...
                  8 4 3 7; ...
                  6 2 1 5];
SurfaceNm = [1;1;2;2;3;3];

%% >> START TO GENERATE CUBOID HYPOTHESIS, VALIDATE WITH OMAP
maxSample = 200000;
validHyps = repmat(struct('extLine', [], 'hCorner', [], 'srcLine', [], ...
                          'omapScr', 0, 'vpLength', [], 'recipe', [], 'gcScr', -1, 'mgScr', -1, 'opScr', -1 ), maxSample, 1);

%% sample line segments, long line segments are more likely to be selected.
rcpNum = size(cube_recipe,1);
for samid = 1:maxSample
    if rem(samid,1000)==0
%         fprintf('%d/%d hypothesis validated, score: %f/%f\n', samid, maxSample, maxScore, numPoint);
        fprintf('%d/%d hypothesis sampled\r', samid, maxSample);
    end
    rcpid = cube_recipe(randsample(rcpNum, 1),:);

    if any([typLine(rcpid).poolLength]==0)
        validHyps(samid).srcLine = [];
        validHyps(samid).recipe = rcpid;
        continue;
    end

    lineID1 = typLine(rcpid(1)).lspool(randsample(typLine(rcpid(1)).poolLength, 1));
    lineID2 = typLine(rcpid(2)).lspool(randsample(typLine(rcpid(2)).poolLength, 1));
    lineID3 = typLine(rcpid(3)).lspool(randsample(typLine(rcpid(3)).poolLength, 1));
    lineID4 = typLine(rcpid(4)).lspool(randsample(typLine(rcpid(4)).poolLength, 1));
    lineID5 = typLine(rcpid(5)).lspool(randsample(typLine(rcpid(5)).poolLength, 1));

    hlines = zeros(12,9);
    hlines(rcpid(1),:) = typLine(rcpid(1)).lines(lineID1,:);
    hlines(rcpid(2),:) = typLine(rcpid(2)).lines(lineID2,:);
    hlines(rcpid(3),:) = typLine(rcpid(3)).lines(lineID3,:);
    hlines(rcpid(4),:) = typLine(rcpid(4)).lines(lineID4,:);
    hlines(rcpid(5),:) = typLine(rcpid(5)).lines(lineID5,:);

    validHyps(samid).srcLine = hlines;
    validHyps(samid).recipe = rcpid;
end
fprintf('\n');

%% finish up unknown edges, compute corner, check if ls locates in line
for samid = 1:maxSample
    if rem(samid,1000)==0
%         fprintf('%d/%d hypothesis validated, score: %f/%f\n', samid, maxSample, maxScore, numPoint);
        fprintf('%d/%d hypothesis validated\r', samid, maxSample);
    end
    if isempty(validHyps(samid).srcLine)
        validHyps(samid).omapScr = -1;
        continue;
    end
    hlines = validHyps(samid).srcLine;
    hCorner = zeros(8,3);
    lineExist = false(12,1);
    lineExist(validHyps(samid).recipe) = true;
    rcpid = validHyps(samid).recipe;
    % >> compute all the other edges, corner
    bDegenerate = false;
    while any(~lineExist)
        hits = lineExist(edgeToPoint);
        actLine = sum(hits,2)==2;
        i = find(actLine,1)';
        srcLineID = edgeToPoint(i,hits(i,:));
        dstLineID = edgeToPoint(i,~hits(i,:));
        dstVpID = find(~hits(i,:));
        dstCor = cross( hlines(srcLineID(1),1:3), hlines(srcLineID(2),1:3), 2);
        if norm(dstCor)<0.1
            bDegenerate = true;
            break;
        end

        hCorner(i,:) = dstCor./norm(dstCor);
        % check if intersection locates at right position
        v = dot([dstCor;dstCor;dstCor],[horizontal;frontNback;leftNright],2);
        if sign(v(1))~=pointDirect(i,1) ...
         ||sign(v(2))~=pointDirect(i,2) ...
         ||sign(v(3))~=pointDirect(i,3)
            hCorner(i,:) = -hCorner(i,:);
        end

        dstNM = cross(hCorner(i,:), vp(dstVpID,:), 2);
        hlines(dstLineID,1:3) = dstNM./norm(dstNM);
        lineExist(dstLineID) = true;
    end

    if bDegenerate
        validHyps(samid).omapScr = -1;
        continue;
    end

    I = find(actLine); i = I(end);
    srcLineID = edgeToPoint(i,hits(i,:));
    dstCor = cross( hlines(srcLineID(1),1:3), hlines(srcLineID(2),1:3), 2);
    hCorner(i,:) = dstCor./norm(dstCor);
    % check if intersection locates at right position
    v = dot([dstCor;dstCor;dstCor],[horizontal;frontNback;leftNright],2);
    if sign(v(1))~=pointDirect(i,1) ...
     ||sign(v(2))~=pointDirect(i,2) ...
     ||sign(v(3))~=pointDirect(i,3)
        hCorner(i,:) = -hCorner(i,:);
    end

    % >> compute extended line representation
    rcpidR = setdiff(1:12, rcpid);
    areaXY = abs(sum(hlines(rcpidR,1:3).*repmat([0 0 1], [length(rcpidR) 1]),2));
    areaYZ = abs(sum(hlines(rcpidR,1:3).*repmat([1 0 0], [length(rcpidR) 1]),2));
    areaZX = abs(sum(hlines(rcpidR,1:3).*repmat([0 1 0], [length(rcpidR) 1]),2));
    [~, planeIDs] = max([areaXY areaYZ areaZX], [], 2); % 1:XY 2:YZ 3:ZX
    hlines(rcpidR,4) = planeIDs;

    extLine = hlines(:,1:6);
    for i = 1:12
        endPoint = hCorner(PointToEdge(i,:),:);
        uv = xyz2uvN(endPoint,extLine(i,4));
        umin = (min(uv(:,1))+pi)/2/pi;    umax = (max(uv(:,1))+pi)/2/pi;
        if umax-umin>0.5
            extLine(i,5) = umax;
            extLine(i,6) = umin;
        else
            extLine(i,5) = umin;
            extLine(i,6) = umax;
        end
    end

    % check if line segments locates in line
    validated = true;
    for i = rcpid
        lineRange = extLine(i,5:6);
        lsegRange = hlines(i,5:6);
        b = insideRange(lsegRange, lineRange);
        if ~b(1) || ~b(2)
            validated = false;
            break;
        end
    end
    if ~validated
        validHyps(samid).omapScr = -1;
        continue;
    end
    validHyps(samid).srcLine = hlines(rcpid,1:6);
    validHyps(samid).extLine = extLine;
    validHyps(samid).hCorner = hCorner;
end
fprintf('\n');

%% compute the checking direction on omap
[candiSetXYZ, ~] = icosahedron2sphere(6);
[ohei, owid, ~] = size(omap);
candiSetUV = uv2coords(xyz2uvN(candiSetXYZ), owid, ohei);
candiInd = sub2ind([ohei owid], candiSetUV(:,2), candiSetUV(:,1));
% convert gc to omap
gndGC = zeros(ohei, owid, 3);
gndGC(:,:,1) = max(gc(:,:,5),gc(:,:,6));
gndGC(:,:,2) = max(gc(:,:,1),gc(:,:,3));
gndGC(:,:,3) = max(gc(:,:,2),gc(:,:,4));
normGndGC = sum(gndGC(candiInd))+sum(gndGC(candiInd+1*ohei*owid))+sum(gndGC(candiInd+2*ohei*owid));
% omap
gndOmap = omap;
gndOmap(gndOmap>0) = 1;
normGndOmap = sum(gndOmap, 3);
gndOmap = gndOmap./(repmat(normGndOmap+0.0001, [1 1 3]));
normGndOmap = sum(gndOmap(candiInd))+sum(gndOmap(candiInd+1*ohei*owid))+sum(gndOmap(candiInd+2*ohei*owid));
numPoint = size(candiSetXYZ,1);

% merge
gndMerge = zeros(ohei, owid, 3);
gndMerge(1:ohei/2,:,:) = gndOmap(1:ohei/2,:,:);
gndMerge(ohei/2+1:ohei, :, :) = gndGC(ohei/2+1:ohei, :, :);
% normMegMap = sum(gndOmap, 3);
% gndOmap = gndOmap./(repmat(normMegMap+0.0001, [1 1 3]));
normMegMap = sum(gndMerge(candiInd))+sum(gndMerge(candiInd+1*ohei*owid))+sum(gndMerge(candiInd+2*ohei*owid));

% optimal merge
gndOptim = zeros(ohei, owid, 3);
gndOptim(1:685,:,:) = gndOmap(1:685,:,:);
gndOptim(686:ohei, :, :) = gndGC(686:ohei, :, :);
% normMegMap = sum(gndOmap, 3);
% gndOmap = gndOmap./(repmat(normMegMap+0.0001, [1 1 3]));
normOptMap = sum(gndOptim(candiInd))+sum(gndOptim(candiInd+1*ohei*owid))+sum(gndOptim(candiInd+2*ohei*owid));

%% check consistency of orientation on uniformly sampled point on sphere
for samid = 1:maxSample
    hypVp = false(numPoint,3);
    if validHyps(samid).omapScr>-0.5
        hCorner = validHyps(samid).hCorner;
        for i = 1:6
%             suf1 = cross(hCorner(PointToSurface(i,2),:), hCorner(PointToSurface(i,1),:), 2);
%             suf2 = cross(hCorner(PointToSurface(i,3),:), hCorner(PointToSurface(i,2),:), 2);
%             suf3 = cross(hCorner(PointToSurface(i,4),:), hCorner(PointToSurface(i,3),:), 2);
%             suf4 = cross(hCorner(PointToSurface(i,1),:), hCorner(PointToSurface(i,4),:), 2);
%             vad1 = dot(candiSetXYZ, repmat(suf1, [numPoint 1]), 2);
%             vad2 = dot(candiSetXYZ, repmat(suf2, [numPoint 1]), 2);
%             vad3 = dot(candiSetXYZ, repmat(suf3, [numPoint 1]), 2);
%             vad4 = dot(candiSetXYZ, repmat(suf4, [numPoint 1]), 2);
%             vadIN = vad1<0 & vad2<0 & vad3<0 & vad4<0;
            vadIN = insideCone( hCorner(PointToSurface(i,[4 3 2 1]),:), candiSetXYZ, 0);
            hypVp( vadIN, SurfaceNm(i)) = true;
        end
%         response = selVp & hypVp;
%         validHyps(samid).omapScr = sum( sum(response,2)>0 )/size(response,1);
        response = sum(gndOmap(candiInd(hypVp(:,1))+0*ohei*owid)) ...
                 + sum(gndOmap(candiInd(hypVp(:,2))+1*ohei*owid)) ...
                 + sum(gndOmap(candiInd(hypVp(:,3))+2*ohei*owid));
        validHyps(samid).omapScr = response/normGndOmap;
        response = sum(gndGC(candiInd(hypVp(:,1))+0*ohei*owid)) ...
                 + sum(gndGC(candiInd(hypVp(:,2))+1*ohei*owid)) ...
                 + sum(gndGC(candiInd(hypVp(:,3))+2*ohei*owid));
        validHyps(samid).gcScr = response/normGndGC;
        response = sum(gndMerge(candiInd(hypVp(:,1))+0*ohei*owid)) ...
                 + sum(gndMerge(candiInd(hypVp(:,2))+1*ohei*owid)) ...
                 + sum(gndMerge(candiInd(hypVp(:,3))+2*ohei*owid));
        validHyps(samid).mgScr = response/normMegMap;
        response = sum(gndOptim(candiInd(hypVp(:,1))+0*ohei*owid)) ...
                 + sum(gndOptim(candiInd(hypVp(:,2))+1*ohei*owid)) ...
                 + sum(gndOptim(candiInd(hypVp(:,3))+2*ohei*owid));
        validHyps(samid).opScr = response/normOptMap;
    end
end

%% save only validate hypothesis, with high score
omapScr = [validHyps.omapScr];
% [B, IX] = sort(omapScr, 'descend');
% selIX = IX(B>-0.5);
% hyps = validHyps(selIX);

hyps = validHyps(omapScr>-0.5);
end

function b = insideRange(pt, range)
range = range + [-0.02 +0.02];
if range(2)>range(1)
    b = pt>=range(1) & pt<=range(2);
else
    b1 = pt>=range(1) & pt<=1;
    b2 = pt>=0 & pt<=range(2);
    b = b1 | b2;
end

end


text:
function value = p32_i4 ( action, name, value )

%*****************************************************************************80
%
%% P32_I4 sets or gets I4 parameters for problem 32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P32_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P32_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P32_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P32_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P32_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P32_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P32_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P32_I4 - Fatal error!' );

  end

  return
end


text:
function [ Y, gR, gPhi, gTheta ] = AFMT3( X, sigma, mAFMT, nAFMT, pAFMT, extrapVal, interpMethod, algorithm )
%Computes discrete approximation to analytic fourier-mellin transformation.
%Parameters:
%           x                   ...matrix with input data
%           sigma               ...strictly positive number, ensuring
%                                  existence of transformation
%           mAFMT               ...size of mellin-transformation (check)
%
%           nAFMT               ...size of fourier-transformation (check)

%algorithm selects: F-AFMT      ...fast AFMT
%                   D-AFMT      ...discrete AFMT
%
% (c) Lukas Mauch, Thomas Kuestner
% ---------------------------------------------------------------------

switch algorithm
    %% forward-transformation (3d):
    case 'F-AFMT'
        [yLogPol, gR, gPhi, gTheta] = convToLogPol(X, mAFMT, nAFMT, pAFMT, interpMethod);
        scaleFactor = exp(2*pi*(gR./max(max(max(gR)))) * sigma); %maximal radius normalized to 2*pi
        Y = fftnshift(scaleFactor .* yLogPol, 1:3);

    case 'D-AFMT'

    %% backward-transformation:
    case 'iF-AFMT'
        [sizeY, sizeX, sizeZ]  = size(X);

        %Compute center positions oM and oN:
        oY = fix(mAFMT/2)+1;
        oX = fix(nAFMT/2)+1;
        oZ = fix(pAFMT/2)+1;

        %Cartesian sample points:
        [gridCartX, gridCartY, gridCartZ] = meshgrid( -oX+1:nAFMT-oX, -oY+1:mAFMT-oY, -oZ+1:pAFMT-oZ );

        %Log-polar plane sample points:
        R  = sqrt(max(abs(oY-1), abs(mAFMT-oY))^2 + max(abs(oX-1), abs(nAFMT-oX))^2 + ...
             max(abs(oZ-1), abs(pAFMT-oZ))^2);
        [gridR, gridPhi, gridTheta]   = meshgrid( log(R)/sizeX:log(R)/sizeX:log(R), ...
                                       linspace(0,2*pi, sizeY), ...
                                       linspace(-pi/2, pi/2, sizeZ));

        scaleFactor = exp(-2*pi*sigma*(gridR./max(max(max(gridR)))));
        yLogPol = scaleFactor .* ifftnshift(X, 1:3);

        Y = convToCart(yLogPol, gridR, gridPhi, gridTheta, gridCartX, gridCartY, gridCartZ, interpMethod);

    case 'iD-AFMT'

    otherwise
    error('Unknown type of algorithm!')

end

    %% helper functions:
    function [xLogPol, gridR, gridPhi, gridTheta] = convToLogPol( xCart, sizeR, sizePhi, sizeTheta, interpMethod )
        [m, n, p]  = size(xCart);

        %Compute center positions oM and oN:
        oM = fix(m/2)+1;
        oN = fix(n/2)+1;
        oP = fix(p/2)+1;

        %Cartesian sample points:
        [gridX, gridY, gridZ] = meshgrid( -oN+1:n-oN, fliplr(-oM+1:m-oM), -oP+1:p-oP );

        %Log-polar plane sample points:
        maxR = sqrt(max(abs(oN-1), abs(n-oN))^2 + max(abs(oM-1), abs(m-oM))^2 + max(abs(oP-1), abs(p-oP))^2);
        [gridR, gridPhi, gridTheta]   = meshgrid( log(maxR)/sizeR:log(maxR)/sizeR:log(maxR), ...
                                      linspace(0,2*pi,sizePhi), ...
                                      linspace(-pi/2, pi/2, sizeTheta));

        %Perform conversion to log-polar coordinates, using linear interpolation:
        convGridX   = exp(gridR) .* cos(gridPhi) .* cos(gridTheta);
        convGridY   = exp(gridR) .* sin(gridPhi) .* cos(gridTheta);
        convGridZ   = exp(gridR) .* sin(gridTheta);
        xLogPol     = interp3(gridX, gridY, gridZ, xCart, ...
                        convGridX, convGridY, convGridZ, interpMethod, extrapVal);
    end

    function xCart = convToCart( xLogPol, gridR, gridPhi, gridTheta, gridX, gridY, gridZ, interpMethod )
        %Perform conversion to cartesian coordinates:
%         convGridX   = reshape(exp(gridR) .* cos(gridPhi) .* cos(gridTheta),[], 1);
%         convGridY   = reshape(exp(gridR) .* sin(gridPhi) .* cos(gridTheta),[], 1);
%         convGridZ   = reshape(exp(gridR) .* sin(gridTheta),[], 1);
%         F = TriScatteredInterp(convGridX, convGridY, convGridZ, reshape(xLogPol, [], 1), interpMethod);
%         xCart       = F(gridX, gridY, gridZ);

        %interpolation in r-phi plane (worse but fast):
        convGridR       = sqrt(abs(gridX).^2 + abs(gridY).^2 + abs(gridZ).^2);
        convGridPhi     = atan(gridY./gridX);
        convGridTheta   = asin(gridZ./convGridR);
        %unwrap:
        qA = logical(logical((gridX >= 0)) .* logical((gridY<=0)));
        qB = logical(logical((gridX < 0)) .* logical((gridY<=0)));
        qC = logical(logical((gridX < 0)) .* logical((gridY>0)));
        qD = logical(logical((gridX >= 0)) .* logical((gridY>0)));
        convGridPhi(qA) = abs(convGridPhi(qA));
        convGridPhi(qB) = pi-convGridPhi(qB);
        convGridPhi(qC) = pi+abs(convGridPhi(qC));
        convGridPhi(qD) = 2*pi-convGridPhi(qD);
        convGridPhi(isnan(convGridPhi)) = 0;
        convGridTheta(isnan(convGridTheta)) = 0;

        xCart       = interp3(exp(gridR), gridPhi, gridTheta, xLogPol, ...
                        convGridR, convGridPhi, convGridTheta, interpMethod, extrapVal);

    end

end

author: thomaskuestner


text:
function varargout = importERPSS_GUI(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @importERPSS_GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @importERPSS_GUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

%---------------------------------------------------------------------------
function importERPSS_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = [];

%
% Name & version
%
version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   IMPORT ERP from ERPSS (ASCII) GUI'])

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

% help
% helpbutton

set(handles.radiobutton_explicit, 'Value', 1);
set(handles.checkbox_transpose, 'Enable', 'off');

% UIWAIT makes importERPSS_GUI wait for user response (see UIRESUME)
uiwait(handles.gui_chassis);

%---------------------------------------------------------------------------
function varargout = importERPSS_GUI_OutputFcn(hObject, eventdata, handles)
% Get default command line output from handles structure
varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.1)

%---------------------------------------------------------------------------
function radiobutton_explicit_Callback(hObject, eventdata, handles)
if get(hObject, 'Value')
        set(handles.radiobutton_implicit, 'Value', 0);
        set(handles.checkbox_transpose, 'Value', 0);
        set(handles.checkbox_transpose, 'Enable', 'off');
else
        set(handles.radiobutton_explicit, 'Value', 1);
end

%---------------------------------------------------------------------------
function radiobutton_implicit_Callback(hObject, eventdata, handles)
if get(hObject, 'Value')
        set(handles.radiobutton_explicit, 'Value', 0);
        set(handles.checkbox_transpose, 'Enable', 'on');
else
        set(handles.radiobutton_implicit, 'Value', 1);
end
%---------------------------------------------------------------------------
function checkbox_transpose_Callback(hObject, eventdata, handles)

%---------------------------------------------------------------------------
function pushbutton_CANCEL_Callback(hObject, eventdata, handles)
handles.output = [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%---------------------------------------------------------------------------
function pushbutton_OK_Callback(hObject, eventdata, handles)
fname = get(handles.edit_file, 'String');

if isempty(fname)
      msgboxText = 'You have to pick a file first.';
      title = 'ERPLAB: import ERPSS Error';
      errorfound(msgboxText, title);
      return
end
if get(handles.radiobutton_explicit, 'Value')==1
        erpssformat = 0;
else
        erpssformat = 1;
end
if get(handles.checkbox_transpose, 'Value')==1
        dtranspose = 1; % transpose
else
        dtranspose = 0; % do not transpose
end
handles.output = {fname, erpssformat, dtranspose};

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function edit_file_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_file_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function pushbutton_browse_Callback(hObject, eventdata, handles)
[filename, filepath, findex] = uigetfile({'*.txt','text file(*.txt)'; ...
      '*.*',  'All Files (*.*)'}, ...
      'Select an ascii file', ...
      'MultiSelect', 'on');

if ~iscell(filename)
        if isempty(filename)
                disp('User selected cancel.')
                return
        end
        if filename==0
                disp('User selected cancel.')
                return
        end
        filename = cellstr(filename);
end

for a=1:length(filename)
fname{a} = fullfile(filepath, filename{a});
end

lista = get(handles.edit_file, 'String');
lista = [lista ;fname'];
lista = lista(~cellfun(@isempty, lista));
set(handles.edit_file, 'String', lista);

%--------------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)
if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')
      %The GUI is still in UIWAIT, us UIRESUME
      handles.output = '';
      %Update handles structure
      guidata(hObject, handles);
      uiresume(handles.gui_chassis);
else
      % The GUI is no longer waiting, just close it
      delete(handles.gui_chassis);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function spm_clf(F)
% Clear specified figure of objects with visible handles
% FORMAT spm_clf(F)
% F - Figure number, or 'Tag' string of figure(s) to clear
%__________________________________________________________________________
%
% Clears the specified figure, deleting all objects with visible
% handles ('HandleVisibility'=='on').
%
% If the current window is 'Tag'ged interactive, then the figures name
% is cleared and the pointer reset.
%
% F Defaults to the current figure, if there is one.
%
% This is just a gateway to spm_figure('Clear',F).
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_clf.m 3899 2010-05-25 15:36:40Z guillaume $

%-Call spm_figure
%--------------------------------------------------------------------------
if nargin==0
    spm_figure('Clear')
else
    spm_figure('Clear',F)
end

author: spm


text:
function value = p24_exact ( dim_num )

%*****************************************************************************80
%
%% P24_EXACT returns the exact integral for problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  value = 1.0;

  return
end


text:
function varargout = findcoord(DEM)

%FINDCOORD Find coordinates of nonzero elements in GRIDobj
%
% Syntax
%
%     [x,y] = findcoord(DEM)
%     [x,y,val] = findcoord(DEM)
%
% Description
%
%     This function returns the coordinates of nonzero values in DEM.
%
% Input arguments
%
%     DEM     GRIDobj
%
% Output arguments
%
%     x,y     coordinate pair
%     val     value of nonzero element in DEM.Z
%
% See also: find, GRIDobj/coord2ind, GRIDobj/ind2coord, GRIDobj/find
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 10. December, 2014

ix = find(DEM);
[x,y] = ind2coord(DEM,ix);

varargout{1} = x;
varargout{2} = y;
if nargout == 3
    varargout{3} = DEM.Z(ix);
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
%% Set Manopt options (if desired)
Manopt_opts.tolgradnorm = 1e-2;  % Stopping tolerance for norm of Riemannian gradient
Manopt_opts.rel_func_tol = 1e-5;  % Additional stopping criterion for Manopt: stop if the relative function decrease between two successive accepted iterates is less than this value
Manopt_opts.miniter = 1;  % Minimum number of outer iterations (i.e. accepted update steps) to perform
Manopt_opts.maxiter = 300;  % Maximum number of outer iterations (i.e. accepted update steps) to perform
Manopt_opts.maxinner = 500;  % Maximum number of iterations for the conjugate-gradient method used to compute approximate Newton steps
%manopt_options.maxtime = 60*60;  % Maximum computation time to allow, in seconds
%manopt_options.solver = @steepestdescent;  % Select Manopt solver to use: {trustregions (default), conjugategradient, steepestdescent}

%% Set SE-Sync options (if desired)
SE_Sync_opts.r0 = 5;  % Initial maximum-rank parameter at which to start the Riemannian Staircase
SE_Sync_opts.rmax = 10;  % Maximum maximum-rank parameter at which to terminate the Riemannian Staircase
SE_Sync_opts.eig_comp_rel_tol = 1e-4;  % Relative tolerance for the minimum-eigenvalue computation used to test for second-order optimality with MATLAB's eigs() function
SE_Sync_opts.min_eig_lower_bound = -1e-3;  % Minimum eigenvalue threshold for accepting a maxtrix as numerically positive-semidefinite
SE_Sync_opts.Cholesky = false;  % Select whether to use Cholesky or QR decomposition to compute orthogonal projections

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function fh = AWB_fh(A,B,w)
% This is almost an MV2DF, but it does not return derivatives on numeric
% input, w.
%
% Algorithm: Y = A*reshape(w,..)*B

if nargin==0
    test_this();
    return;
end

[m,n] = size(A);
[r,s] =  size(B);

    function y = map_this(w)
        w = reshape(w,n,r);
        y = A*w*B;
    end

    function w = transmap_this(y)
        y = reshape(y,m,s);
        w = A.'*y*B.';
    end

map = @(y) map_this(y);
transmap = @(y) transmap_this(y);

fh = linTrans([],map,transmap);

if exist('w','var') && ~isempty(w)
    fh = fh(w);
end

end

function test_this()
A = randn(2,3);
B = randn(4,5);
f = AWB_fh(A,B);
test_MV2DF(f,randn(3*4,1));
end

author: nesl


text:
function out = isempty( f )
%ISEMPTY   True for empty SEPARABLEAPPROX.
%   ISEMPTY(F) returns 1 if F is an empty SEPARABLEAPPROX object and 0 otherwise.

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

out = isempty( f.cols ) && isempty( f.rows );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% P2.7
clear;  clc;
n1 = 0:20;
x1 = 0.9.^n1;
subplot(3,2,1); stem(n1,x1,'ko');  title('First given sequence x1(n)');
axis([-20 20 0 1]);
[x1f,n1f] = sigfold(x1,n1);                 % obtain x1(-n)
n2 = -20:0;
x2 = 0.8.^(-n2);
subplot(3,2,2); stem(n2,x2,'ko');  title('Second given sequence x2(n)');
axis([-20 20 0 1]);
[x2f,n2f] = sigfold(x2,n2);                 % obtain x2(-n)
[r,nr] = conv_m(x1,n1,x1f,n1f);             % auto-correlation
subplot(3,2,3); stem(nr,r,'ko');  title('Auto-Correlation x1(n)*x1(-n)');
[r,nr] = conv_m(x2,n2,x2f,n2f);             % auto-correlation
subplot(3,2,4); stem(nr,r,'ko');  title('Auto-Correlation x2(n)*x2(-n)');
[r,nr] = conv_m(x1f,n1f,x2,n2);             % cross-correlation
subplot(3,1,3); stem(nr,r,'ko');  title('Cross-Correlation');

author: Sable

repo: mcbench-benchmarks


text:
close all; clear all; clc;
rng('default');
N = 64;
D = 121;
dict = spx.dict.simple.gaussian_dict(N, D);
props = spx.dict.Properties(dict);
G = props.gram_matrix();
F = props.frame_operator();
S = props.singular_values();
lambda = props.gram_eigen_values();
fprintf('Lower frame bound: %f\n', props.lower_frame_bound());
fprintf('Upper frame bound: %f\n', props.upper_frame_bound());
mu = props.coherence();
fprintf('Coherence: %f\n', mu);
% fprintf('Eigen values of Gram matrix:\n');
% fprintf('%f ', lambda);
% fprintf('\n');
mf = spx.graphics.Figures();
mf.new_figure('Eigen values of Gram matrix');
plot(lambda, '.');
title('Eigen values of Gram matrix');
xlabel('Number');
ylabel('Eigen value');

mf.new_figure('Singular values of dictionary');
plot(props.singular_values(), '.');
title('Singular values of dictionary');
xlabel('Number');
ylabel('Singular value');

author: indigits

repo: sparse-plex


text:
filename='Circle_Quadrilateral_Linear_Structured_32';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'rectangle';
widthV = 3;
widthH = 0.5;
cost = {'compliance'};
weights = [1];
constraint = {'volumeConstraint'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 1;
Perimeter_target=3.5;
optimality_final =1e-5;
constr_final =1e-5;

BCscale_factor = 0.3;
HJiter0 = 1;
e2 = 5;

Vfrac_initial = 1;
optimality_initial = 5e-2;
constr_initial = 5e-2;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;

% maxiter = 10;
% maxiter = 1;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

plotting = 0;
printing = 0;
monitoring = 0;
monitoring_interval = 0;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Entropy] = Energy_Entropy_Block(f,winLength,winStep,numOfShortBlocks)

f = f / max(abs(f));
Eol = sum(f.^2);
L = length(f);

if (winLength==0)
    winLength = floor(L);
    winStep = floor(L);
end

numOfBlocks = (L-winLength)/winStep + 1;
curPos = 1;
for (i=1:numOfBlocks)
    curBlock = f(curPos:curPos+winLength-1);
    for (j=1:numOfShortBlocks)
        s(j) = sum(curBlock((j-1)*(winLength/numOfShortBlocks)+1:j*(winLength/numOfShortBlocks)).^2)/Eol;
    end

    Entropy(i) = -sum(s.*log2(s));
    curPos = curPos + winStep;
end


text:
learnLoc = which('startLearning.m');
parentDir = fileparts(learnLoc);
% ADD PATHS
addpath(genpath(parentDir));

% SAVE CONSTANTS (SEE medalConstants.m)
root = parentDir;
data = fullfile(parentDir,'data');
modules = fullfile(parentDir,'modules');
save(fullfile(parentDir,'.config.mat'));

author: dustinstansbury

repo: medal


text:
classdef ADRV9009Tests < HardwareTests

    properties
        uri = 'ip:analog';
        SamplingRateRX = 245.76e6;
        author = 'ADI';
    end

    methods(TestClassSetup)
        % Check hardware connected
        function CheckForHardware(testCase)
            Device = @()adi.ADRV9009.Rx;
            testCase.CheckDevice('ip',Device,testCase.uri(4:end),false);
        end
    end

    methods (Static)
        function estFrequency(data,fs)
            nSamp = length(data);
            FFTRxData  = fftshift(10*log10(abs(fft(data))));
            df = fs/nSamp;  freqRangeRx = (-fs/2:df:fs/2-df).'/1000;
            plot(freqRangeRx, FFTRxData);
        end
    end

    methods (Test)

        function testADRV9009Rx(testCase)
            % Test Rx DMA data output
            rx = adi.ADRV9009.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            [out, valid] = rx();
            rx.release();
            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
        end

        function testADRV9009RxCustomProfile1(testCase)
            % Test Rx custom profiles
            rx = adi.ADRV9009.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            rx.EnableCustomProfile = true;
            rx.CustomProfileFileName = ...
                'Tx_BW200_IR245p76_Rx_BW200_OR245p76_ORx_BW200_OR245p76_DC245p76.txt';
            [out, valid] = rx();
            rxSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',false);
            txSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',true);
            rx.release();
            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(rxSampleRate,int64(245.76e6),...
                'Invalid sample rate after profile write');
            testCase.verifyEqual(txSampleRate,int64(245.76e6),...
                'Invalid sample rate after profile write');
        end

        function testADRV9009RxCustomProfile2(testCase)
            % Test Rx custom profiles
            rx = adi.ADRV9009.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            rx.EnableCustomProfile = true;
            rx.CustomProfileFileName = ...
                'Tx_BW400_IR491p52_Rx_BW200_OR245p76_ORx_BW400_OR491p52_DC245p76.txt';
            [out, valid] = rx();
            rxSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',false);
            txSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',true);
            rx.release();
            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(rxSampleRate,int64(245.76e6),...
                'Invalid sample rate after profile write');
            testCase.verifyEqual(txSampleRate,int64(245.76e6*2),...
                'Invalid sample rate after profile write');
        end
%
%         function testADRV9009RxCustomProfileZC706(testCase)
%             % Test Rx custom profiles
%             rx = adi.ADRV9009.Rx('uri',testCase.uri);
%             rx.EnabledChannels = 1;
%             rx.EnableCustomProfile = true;
%             rx.CustomProfileFileName = ...
%                 'Tx_BW100_IR122p88_Rx_BW100_OR122p88_ORx_BW100_OR122p88_DC122p88.txt';
%             [out, valid] = rx();
%             rxSampleRate = rx.getAttributeLongLong('voltage0',...
%                 'sampling_frequency',false);
%             rx.release();
%             testCase.verifyTrue(valid);
%             testCase.verifyGreaterThan(sum(abs(double(out))),0);
%             testCase.verifyEqual(rxSampleRate,int64(122880000),...
%                 'Invalid sample rate after profile write');
%         end

        function testADRV9009RxWithTxDDS(testCase)
            % Test DDS output
            tx = adi.ADRV9009.Tx('uri',testCase.uri);
            tx.DataSource = 'DDS';
            toneFreq = 30e6;
            tx.DDSFrequencies = repmat(toneFreq,2,4);
            tx.AttenuationChannel0 = -10;
            tx.EnableCustomProfile = true;
            tx.CustomProfileFileName = ...
                'Tx_BW200_IR245p76_Rx_BW200_OR245p76_ORx_BW200_OR245p76_DC245p76.txt';
            tx();
            pause(1);
            rx = adi.ADRV9009.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            rx.kernelBuffersCount = 1;
            for k=1:10
                valid = false;
                while ~valid
                    [out, valid] = rx();
                end
            end

            rxSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',false);

            rx.release();

%             plot(real(out));
%             testCase.estFrequency(out,testCase.SamplingRateRX);
            freqEst = meanfreq(double(real(out)),rxSampleRate);

            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(freqEst,toneFreq,'RelTol',0.01,...
                'Frequency of DDS tone unexpected')

        end

        function testADRV9009RxWithTxDDSTwoChan(testCase)
            % Test DDS output
            tx = adi.ADRV9009.Tx('uri',testCase.uri);
            tx.DataSource = 'DDS';
            toneFreq1 = 30e6;
            toneFreq2 = 85e6;

            tx.DDSFrequencies = [toneFreq1,toneFreq1,toneFreq2,toneFreq2;...
                0,0,0,0];
            tx.DDSScales = [1,1,1,1;...
                0,0,0,0];
            tx.DDSPhases = [90000,0,90000,0;...
                0,0,0,0];

            tx.AttenuationChannel0 = -10;
            tx.EnableCustomProfile = true;
            tx.CustomProfileFileName = ...
                'Tx_BW200_IR245p76_Rx_BW200_OR245p76_ORx_BW200_OR245p76_DC245p76.txt';
            tx();
            pause(1);
            rx = adi.ADRV9009.Rx('uri',testCase.uri);
            rx.EnabledChannels = [1,2];
            rx.kernelBuffersCount = 1;
            for k=1:10
                valid = false;
                while ~valid
                    [out, valid] = rx();
                end
            end

            rxSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',false);

            rx.release();

%             plot(real(out));
%             testCase.estFrequency(out,rxSampleRate);
            freqEst = meanfreq(double(real(out)),rxSampleRate);

            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(freqEst(1),toneFreq1,'RelTol',0.01,...
                'Frequency of DDS tone unexpected')
            testCase.verifyEqual(freqEst(2),toneFreq2,'RelTol',0.01,...
                'Frequency of DDS tone unexpected')
        end

        function testADRV9009RxWithTxData(testCase)
            % Test Tx DMA data output
            amplitude = 2^15; frequency = 20e6;
            swv1 = dsp.SineWave(amplitude, frequency);
            swv1.ComplexOutput = true;
            swv1.SamplesPerFrame = 2^20;
            swv1.SampleRate = 245.76e6;
            y = swv1();

            tx = adi.ADRV9009.Tx('uri',testCase.uri);
            tx.DataSource = 'DMA';
            tx.EnableCyclicBuffers = true;
            tx.AttenuationChannel0 = -10;
            tx.EnableCustomProfile = true;
            tx.CustomProfileFileName = ...
                'Tx_BW200_IR245p76_Rx_BW200_OR245p76_ORx_BW200_OR245p76_DC245p76.txt';
            tx(y);
            rx = adi.ADRV9009.Rx('uri',testCase.uri);
            rx.EnabledChannels = 1;
            rx.kernelBuffersCount = 1;
            for k=1:20
                valid = false;
                while ~valid
                    [out, valid] = rx();
                end
            end

            rxSampleRate = rx.getAttributeLongLong('voltage0',...
                'sampling_frequency',false);

            rx.release();

%             plot(real(out));
%             testCase.estFrequency(out,testCase.SamplingRateRX);
            freqEst = meanfreq(double(real(out)),rxSampleRate);

            testCase.verifyTrue(valid);
            testCase.verifyGreaterThan(sum(abs(double(out))),0);
            testCase.verifyEqual(freqEst,frequency,'RelTol',0.01,...
                'Frequency of ML tone unexpected')
        end
    end

end

author: analogdevicesinc


text:
filename='Cantileverbeam_Tetrahedra_Linear_Structured';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

HJiter0 = 1;
e2 = 100;
N_holes = [12 5 5];
R_holes = 0.2;
phase_holes = [0 0 0];
nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 0.7;
optimality_initial = 5e-2;
constr_initial = 5e-2;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function [mag, freq, h, h2, peak_freq] = fft_calc(dat, TR, varargin)
% Simple function to calculate the FFT power of a
% data vector (dat) as a function of frequency,
% given a sample-to-sample repetition time (TR)
%
% :Usage:
% ::
%
%     [mag, freq, line_handle, nyquist_line_handle, peak_freq] = fft_calc(dat, TR)
%
% Example:
% Create and plot a sin wave:
% ------------------------------------------------------------------
% Fs = 1000;            % Sampling frequency, e.g., Hz (samples/sec)
% T = 1/Fs;             % Sampling period (time between samples, e.g., sec)
% L = 1000;             % Length of signal (in samples)
% t = (0:L-1)*T;        % Time vector
%
% % Sine wave parameters
% theta = 0;            % Phase, in radians
% F = 10;               % Frequency of sin wave (cycles/sec)
%
% y = sin(2 * pi * F * t + theta);
%
% create_figure('sin wave');
% plot(t, y)
% xlabel('Time (sec)')
%
% Downsample the sin wave at new frequency Fs and plot signal and FFT:
% ------------------------------------------------------------------
% Fs = 20;  % New sampling frequency
% downsample_by = round(1000/Fs);
%
% figure;
% plot(t, y, 'k-'); hold on
% set(gcf, 'Position', [0 644        1358         231])
%
% plot(t(1:downsample_by:end), y(1:downsample_by:end), '.-', 'Color', [1 .3 .6], 'LineWidth', 2);
% Nyquist = Fs/2;
% title(sprintf('Sampling rate = %3.1f, Nyquist limit = %3.2f', Fs, Nyquist))
%
% figure; hold on;
% [mag, freq, line_han] = fft_calc(y, 1/1000);
% set(line_han, 'Color', 'k', 'LineWidth', 2, 'LineStyle', '-');
%
% [mag, freq, line_han] = fft_calc(y(1:downsample_by:end), 1/Fs);
% set(line_han, 'Color', [1 .3 .6], 'LineWidth', 2, 'LineStyle', '-');
%
% set(gca, 'XLim', [0 (max(Fs, 1.2 * F))]);

n = length(dat);

if ~isempty(varargin)
    if ischar(varargin{1}) && strcmp(varargin{1}, 'plot')
    else
        doplot = varargin{1};
    end
end

% from matlab example
% power = abs(Y(1:floor(n/2))).^2;
% nyquist = 1/2;
% freq = (1:n/2)/(n/2)*nyquist

nyq = 1 ./ (2 * TR);  % Sampling rate frequency / 2. TR = sample-to-sample time, 1/Fs

timepts = floor(n ./ 2);

freq = (0:timepts-1)/timepts * nyq;
%freq = linspace(0, nyq, timepts)';

mag = fft(dat); %real(abs(fft(dat)));
mag = abs(mag(1:timepts)) .^ 2;  % power

mag = mag ./ sum(mag);

h = plot(freq, mag, 'LineWidth', 2);
xlabel('Frequency (Hz)')
title('Frequency domain')

% nyquist = Fs ./ 2;  % Any signal faster than Nyquist freq will be aliased
hold on
h2 = plot_vertical_line(nyq);
set(h2, 'LineStyle', ':')

peak_freq = freq(mag == max(mag));

end


text:
function test_matlab_ci

% WALLTIME 00:05:00
% MEM 1gb
% DEPENDENCY

disp('--------------------------------------------------------------------');
ver;
disp('--------------------------------------------------------------------');
disp(pwd);
disp('--------------------------------------------------------------------');
disp(ft_version);
disp('--------------------------------------------------------------------');
system('uname -a');
disp('--------------------------------------------------------------------');
system('uptime');
disp('--------------------------------------------------------------------');
system('ls $HOME');
disp('--------------------------------------------------------------------');
system('env');
disp('--------------------------------------------------------------------');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
addpath(genpath('../mutils/My/'));
addpath(genpath('../ptv'));
%%
dat = load('data/dti.mat');
imgs = double(reshape(dat.imgs, [size(dat.imgs, 1), size(dat.imgs, 2), 1, 1, size(dat.imgs,3)]));

opts = [];
opts.check_gradients = 155*0;
opts.pix_resolution = [1,1];
opts.metric = 'nuclear';

opts.grid_spacing = [10, 10];
opts.isoTV = 5e-2;
opts.mean_penalty = 1e-3;

opts.spline_order = 1;
opts.max_iters = 100;
opts.display = 'off';

opts.border_mask = 6;
opts.k_down = 0.7;

tic
[voldef_pl, Tmin_pl,  Kmin_pl] = ptv_register(imgs, [], opts);
toc

%%
savegif('results/dti_register.gif', squeeze([imgs, voldef_pl])*1, 1/20);
%%
dat = load('data/dti2.mat');
imgs = double(abs( reshape(dat.data_to_save(:,:, 2, :, :), 55,116, [])));
imgs = double(reshape(imgs, [size(imgs, 1), size(imgs, 2), 1, 1, size(imgs,3)]));

imgs = imgs ./ max(max(imgs, [], 1), [], 2);

opts = [];
opts.check_gradients = 155*0;
opts.pix_resolution = [1,1];
opts.metric = 'nuclear';

opts.grid_spacing = [6, 6];
opts.isoTV = 5e-3;
opts.mean_penalty = 1e-5*0;

opts.spline_order = 1;
opts.max_iters = 100;
opts.display = 'off';

opts.border_mask = 6;
opts.k_down = 0.7;

tic
[voldef_pl, Tmin_pl,  Kmin_pl] = ptv_register(imgs, [], opts);
toc

%%
savegif('results/dti_register2.gif', squeeze([imgs, voldef_pl])*2, 1/20);

repo: pTVreg


text:
function run_test_single_model
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% Given a model, evaluate the performance.
    baseDir = '../../../';
    addpath([baseDir, filesep, 'codes']);
    addpath([baseDir, filesep, 'tools', filesep, 'bss_eval']);
    addpath([baseDir, filesep, 'tools', filesep, 'bss_eval_3']);
    addpath([baseDir, filesep, 'tools', filesep, 'labrosa']);
    addpath([baseDir, filesep, 'codes', filesep, 'TSP']);

    ModelPath=[baseDir, filesep, 'codes',filesep, 'TSP', filesep, 'model_demo'];

    global SDR;
    SDR.deviter=0;   SDR.devmax=0;   SDR.testmax=0;
    SDR.devsar=0; SDR.devsir=0; SDR.testsar=0; SDR.testsir=0;

    j=7650;

    % Load model
    load model_RNN1_win1_h300_l2_r0_64ms_1000000_softabs_linearout_RELU_logmel_trn0_c1e-10_c0.001_bsz100000_miter10_bf50_c0_d0_7650.mat
    eI.writewav=1;
    eI.bss3=1;
    eI.DataPath=[baseDir, filesep, 'codes', filesep, 'TSP', ...
        filesep, 'Data', filesep];
    eI.saveDir = [baseDir, filesep, 'codes', filesep, 'TSP', ...
        filesep, 'demo', filesep, 'results', filesep];
    eI.CFGPath = [baseDir, filesep, 'tools', filesep, 'htk_features', filesep];
    test_TSP_general_kl_recurrent(eI.modelname, theta, eI, 'done', j);
end

author: posenhuang

repo: deeplearningsourceseparation


text:
% USERDATA  User data for SLAMTB.
%   Edit this script to enter the information you need for SLAM. Variable
%   names and comments should make this file easy to understand. Follow
%   these guidelines:
%
%   * Specify site and estimation details for the current run in ExpOpt.
%   * Specify sampling time and start and end frames in Tim.
%   * Use as many robots and sensors as you wish with Robot{} and Sensor{}.
%   * Assign sensors to robots via Sensor{i}.robot.
%   * Use field Sensor{i}.distortion for radial distortion parameters if
%   desired.
%   * Use the field Opt.map.numLmk and .lmkSize to specify the maximum
%   number of landmarks that the SLAM map must support.
%   * Use Opt.init.initType to select the type of landmarks to use. Try
%   with one in this list:
%       'idpPnt', 'hmgPnt', 'ahmPnt', 'plkLin', 'ahmLin'.
%   * Use World.points and World.segments to create artificial worlds of
%   points or segments. Check functions THICKCLOISTER and HOUSE.
%
%   See further comments within the file for more detailed information.
%
%   NOTE: You can have multiple copies of this file with different names to
%   store different simulation conditions. Just modify the call in SLAMTB
%   to point to the particular 'USERDATA' file you want.
%
%   See also SLAMTB, EULERANGLES, THICKCLOISTER, HOUSE, USERDATAPNT,
%   USERDATALIN.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% Time variables
%   - sampling time, first and last frames
Time = struct(...
  'dt',                   .1,...          % sampling time, seconds
  'firstFrame',           1,...           % first frame #
  'lastFrame',            800);           % last frame #

% Simulated world
%   - Simulation landmark sets, playground dimensions
World = struct(...
  'points',           thickCloister(-6,6,-6,6,1,7),... % 3d point landmarks - see THICKCLOISTER.
  'segments',         []);  % 3D segments - see HOUSE.

% Robot things with their controls
%   - each robot's type and initial configuration, and controls.
%   - motion models (add new model strings if you need more):
%       'constVel'    6D Constant velocity model
%       'odometry'    6D Odometry model
%       'inertial'    6D IMU-based model
%   - See EULERANGLES for orientations specifications.
Robot{1} = struct(...                  % ODOMETRY EXAMPLE
  'id',                 1,...           % robot identifier
  'name',               'Dala',...      % robot name
  'type',               'atrv',...      % type of robot
  'motion',             'odometry',...  % motion model
  'position',           [0;-5;0],...     % robot position in map
  'orientationDegrees', [0;0;0],...     % orientation, in degrees, [roll; pitch; yaw].
  'positionStd',        [0;0;0],...     % position error, std
  'orientationStd',     [0;0;0],...     % orient. error, std, in degrees
  'dx',                 [.08;0;0],...     % position increment
  'daDegrees',          [0;0;.9],...     % angle increment, degrees
  'dxStd',              0.005*[1;1;1],...  % odo linear error std
  'daStd',              0.05*[1;1;1]);      % odo ang error std, degrees

% Robot{2} = struct(...                  % ODOMETRY EXAMPLE
%   'id',                 2,...           % robot identifier
%   'name',               'Dala',...      % robot name
%   'type',               'atrv',...      % type of robot
%   'motion',             'odometry',...  % motion model
%   'position',           [0;5;0],...     % robot position in map
%   'orientationDegrees', [0;0;180],...     % orientation, in degrees, [roll; pitch; yaw].
%   'positionStd',        [0;0;0],...     % position error, std
%   'orientationStd',     [0;0;0],...     % orient. error, std, in degrees
%   'dx',                 [.08;0;0],...     % position increment
%   'daDegrees',          [0;0;.9],...     % angle increment, degrees
%   'dxStd',              0.005*[1;1;1],...  % odo linear error std
%   'daStd',              0.05*[1;1;1]);      % odo ang error std, degrees

% Robot{3} = struct(...                  % CONSTANT VELOCITY EXAMPLE
%   'id',                 3,...           % robot identifier
%   'name',               'Dala',...      % robot name
%   'type',               'atrv',...      % type of robot
%   'motion',             'constVel',...  % motion model
%   'position',           [1;0;0],...     % robot position in map
%   'orientationDegrees', [0;0;45],...    % orientation, in degrees, [roll; pitch; yaw].
%   'positionStd',        [0;0;0],...     % position error, std
%   'orientationStd',     [0;0;0],...     % orient. error, std, degrees
%   'velocity',           [1;0;0],...     % lin. velocity
%   'angularVelDegrees',  [0;0;10],...    % ang. velocity, in degrees
%   'velStd',             [0;0;0],...     % lin. vel. error, std
%   'angVelStd',          [0;0;0],...     % ang. vel. error, std, degrees
%   'dv',                 [0;0;0],...     % veolcity increment
%   'dwDegrees',          [0;0;0],...     % ang. vel. increment, degrees
%   'dvStd',              [0;0;0],...     % vel perturbation std
%   'dwStd',              [0;0;1]);       % ang vel pert. std, degrees

% Sensor things
%   - each sensor's type and parameters, noise, non-measurable prior.
%   - Sensor types (add new type strings if you need more):
%       'pinHole'   Pin-hole camera
%   - See EULERANGLES for orientations specifications.
Sensor{1} = struct(...
  'id',                 1,...           % sensor identifier
  'name',               'Micropix',...  % sensor name
  'type',               'pinHole',...   % type of sensor
  'robot',              1,...           % robot where it is mounted
  'position',           [0;0;.6],...    % position in robot
  'orientationDegrees', [-90;0;-90],... % orientation in robot, [roll; pitch; yaw]
  'positionStd',        [0;0;0],...     % position error std
  'orientationStd',     [0;0;0],...     % orient. error std
  'imageSize',          [640;480],...   % image size
  'pixErrorStd',        1.0,...         % pixel error std
  'intrinsic',          [320;240;320;320],... % intrinsic params [u0 v0 au av]
  'distortion',         [-0.3;0.1],...          % distortion params
  'frameInMap',         false,...       % add sensor frame in slam map?
  'imGrid',               struct(...      % grid for Active Search
    'numCells',         [8;6],...         % number of H and V grid cells
    'skipOuter',        true));           % skip outer cells for initialization?

% Sensor{2} = struct(...
%   'id',                 2,...           % sensor identifier
%   'name',               'Micropix',...      % sensor name
%   'type',               'pinHole',...   % type of sensor
%   'robot',              2,...           % robot where it is mounted
%   'position',           [0;-0.25;.6],...     % position in robot
%   'orientationDegrees', [-90;0;-90],...      % orientation in robot, [roll; pitch; yaw]
%   'positionStd',        [0;0;0],...     % position error std
%   'orientationStd',     [1.5;1.5;1.5],...     % orient. error std
%   'imageSize',          [640;480],...   % image size
%   'pixErrorStd',        1.0,...         % pixel error std
%   'intrinsic',          [320;240;320;320],... % intrinsic params
%   'distortion',         [],...          % distortion params
%   'frameInMap',         true,...         % add sensor frame in slam map?
%   'imGrid',               struct(...      % grid for Active Search
%     'numCells',         [8;6],...         % number of H and V grid cells
%     'skipOuter',        true));           % skip outer cells for initialization?

% Sensor{3} = struct(...
%   'id',                 1,...           % sensor identifier
%   'name',               'Micropix',...  % sensor name
%   'type',               'pinHole',...   % type of sensor
%   'robot',              2,...           % robot where it is mounted
%   'position',           [0;.25;.6],...    % position in robot
%   'orientationDegrees', [-90;0;-90],... % orientation in robot, [roll; pitch; yaw]
%   'positionStd',        [0;0;0],...     % position error std
%   'orientationStd',     [0;0;0],...     % orient. error std
%   'imageSize',          [640;480],...   % image size
%   'pixErrorStd',        1.0,...         % pixel error std
%   'intrinsic',          [320;240;320;320],... % intrinsic params [u0 v0 au av]
%   'distortion',         [-0.3;0.1],...          % distortion params
%   'frameInMap',         false,...       % add sensor frame in slam map?
%   'imGrid',               struct(...      % grid for Active Search
%     'numCells',         [8;6],...         % number of H and V grid cells
%     'skipOuter',        true));           % skip outer cells for initialization?

% Omnidirectional camera model -> MegaPixel Fish Eye lens as example (FoV is ~190 deg )
% k           = [668 438 1 0 0]';
% invProjDist = [4.473e+2 -0.000e+0 -1.068e-3 1.184e-6 -1.856e-9]';
% projDist    = [6.447e+2 -3.410e+2 -2.901e+1 -5.770e+1 1.849e+1 5.415e+0 5.065e+1 -5.614e+1 1.591e+1 0 0]';
% Sensor{1} = struct(...
%   'id',                 1,...           % sensor identifier
%   'name',               'FrontCam',...  % sensor name
%   'type',               'omniCam',...   % type of sensor
%   'robot',              1,...           % robot where it is mounted
%   'position',           [0.2;0;1.2],... % position in robot
%   'orientationDegrees', [-120;0;-90],...% orientation in robot, [roll; pitch; yaw].
%   'positionStd',        [0;0;0],...     % position error std
%   'orientationStd',     [0;0;0],...     % orient. error std
%   'imageSize',          [1280;800],...  % image size
%   'pixErrorStd',        1.0,...         % pixel error std
%   'intrinsic',          k,...           % intrinsic params: [xc yc c d e]';
%   'distortion',         projDist,...    % distortion params -> polynom for projection to cam sensor
%   'invDistortion',      invProjDist,... % distortion params -> polynom for inv proj from cam sensor
%   'frameInMap',         false,...       % add sensor frame in slam map?
%   'imGrid',             struct(...      % grid for Active Search
%     'numCells',         [8;6],...       % number of H and V grid cells
%     'skipOuter',        true));         % skip outer cells for initialization?

% Estimation options
Opt = struct(...
  'map',              struct(...    % options for the map
    'type',           'ekf',...      % type of map {'ekf','graph'}
    'numLmks',        73,...         % number of 3d landmarks
    'lmkSize',        6),...          % Size of landmark state
  'correct',          struct(...    % options for lmk correction
    'reprojectLmks',  true,...       % reproject lmks after active search?
    'reparametrize',  true,...       % reparametrize lmk?
    'nUpdates',       10,...         % max simultaneus updates
    'MD2th',          9,...          % Threshold on Mahalanobis distance squared
    'linTestIdp',     0.1,...        % threshold on IDP linearity test
    'lines',          struct(...     % options for line corrections
      'innType',      'ortDst',...    % innovation type for lines
      'extPolicy',    false,...       % line extending policy ?
      'extSwitch',    10)),...        % extension policy switch point in pixels
  'init',             struct(...    % Options for initialization
    'nbrInits',       [5 1],...      % number of inits [firstFrame, otherFrames]
    'initType',       'idpPnt',...   % Type of lmk to use for init
    'idpPnt',         struct(...     % options for lmk initialization
      'nonObsMean',   .01,...         % mean of non obs
      'nonObsStd',    .5),...         % std of non obs
    'plkLin',         struct(...     % opt. for Plucker and anchored Plucker lines init
      'nonObsMean',   [.1;0],...      % mean of non obs
      'nonObsStd',    [.25;1])),...   % std of non obs
  'obs',              struct(...    % Observation options
    'lines',          struct(...     % lines options
      'minLength',    20)));          % minimum segment length

% Simulation options
%   - random
SimOpt = struct(...
  'random',           struct(...      % random generator options
    'newSeed',        true,...         % select new random seed?
    'fixedSeed',      208948,...            % random seed for non-random runs
    'seed',           []),...          % current seed
  'obs',              Opt.obs);       % Observation options

% Figure options
%   - view, projection, video, ellipses.
%   - figure projections - mapProj:
%       'persp'     Perspective
%       'ortho'     Orthographic
%   - 3D figure views - mapView - see MAPOBSERVER.
%       [a,e,f]     Custom azimuth/elevation/FOV vector. Distance automatic
%       [a,e,f,d]   custom az/el/fov/distance vector.
%   - 3D figure predefined views (edit mapObserver.m to create/edit views):
%       'top'       Top view
%       'side'      Side view
%       'view'      Generic view
%       'normal'    Normal view
%   - objects colors - two options for color specification:
%       'rgbcmykw'  1-char predifined Matlab colors
%       [r g b]     RGB color vector. [0 0 0] is black, [1 1 1] is white.
FigOpt = struct(...
  'renderer',       'zbuffer',...    % renderer
  'rendPeriod',     1,...           % frames to skip for faster processing
  'createVideo',    false,...       % create video sequences?
  'map',            struct(...      % map figure options
    'size',         [320 240],...   % map figure size
    'lims',         struct(...      % playground limits
      'xMin',            -10,...
      'xMax',             10,...
      'yMin',            -10,...
      'yMax',             10,...
      'zMin',            -10,...
      'zMax',             10),...
    'proj',         'persp',...     % projection of the 3d figure
    'view',         'view',...      % viewpoint of the 3d figure [30 45 40 20]
    'orbit',        [0 0],...       % Azimuth and Elevation orbit angle increments - use to animate figure
    'showSimLmk',   true,...        % show simulated landmarks?
    'showEllip',    false,...        % show ellipsoids?
    'colors',       struct(...      % map figure colors
      'border',     [1 1 1],...      %   [r g b]
      'axes',       [0 0 0],...      % with:
      'bckgnd',     [1 1 1],...      %   [0 0 0] black
      'simLmk',     .3*[1 1 1],...   %   [1 1 1] white
      'defPnt',     struct(...       % euclidean point colors
        'mean',     'b',...           % mean dot
        'ellip',    [.7 .7 1]),...    % ellipsoid
      'othPnt',     struct(...       % other point colors
        'mean',     'r',...           % mean dot
        'ellip',    [1 .7 .7]),...    % ellipsoid
      'defLin',     struct(...       % Plucker line colors
        'mean',     [0 .8 0],...      % mean line
        'ellip',    [.6 1 .6]),...    % ellipsoid
      'othLin',     struct(...       % Plucker line colors
        'mean',     [.8 0 0],...      % mean line
        'ellip',    [1 .6 .6]),...    % ellipsoid
      'simu',       'b',...          %   or 'r', 'b', etc.
      'est',        'g',...          % estimated robots and sensors
      'ground',     [.8 .8 .8],...   % simulated robots and sensors
      'label',      [.0 .5 0])),...  % landmark ID labels
  'sensor',         struct(...      % sensor figures options
    'size',         [320 240],...    % sensor figure size
    'showEllip',    false,...        % show ellipses?
    'colors',       struct(...       % Sensor figure colors:
      'border',     .8*[1 1 1],...    %
      'axes',       [0 0 0],...       %
      'bckgnd',     [1 1 1],...       %
      'raw',        .3*[1 1 1],...    %
      'defPnt',     struct(...       % Default point colors
        'updated',  'c',...           % updated
        'predicted','b'),...          % predicted
      'othPnt',     struct(...       % other point colors
        'updated',  'r',...           % updated
        'predicted','m'),...          % predicted
      'defLin',     struct(...       % Default line colors
        'meas',     'b',...           % measurement
        'mean',     'g',...           % mean line
        'ellip',    'y'),...          % ellipsoid
      'othLin',     struct(...       % other line colors
        'meas',     'b',...           % measurement
        'mean',     'm',...           % mean line
        'ellip',    'r'),...          % ellipsoid
      'label',      [.5 .5 .5])));    %

% Experiment options
%   - site name, series gathered, estimation run number
ExpOpt = struct(...
  'root',         '~/SLAM/',...         % root directory
  'site',         'simu',...            % Name of the site
  'dataRun',      1,...                 % Run # on this site
  'estimateRun',  1,...                 % slam run for data and site
  'lmkTypes',     Opt.init.initType,... % types of landmarks used
  'sensingType',  'mono',...            % sensing mode
  'mappingType',  'single');            % mapping mode

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [varargout] = runopf_w_res(varargin)
%RUNOPF_W_RES  Runs an optimal power flow with fixed zonal reserves.
%   RESULTS = RUNOPF_W_RES(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
%   [RESULTS, SUCCESS] = RUNOPF_W_RES(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
%
%   Runs an optimal power flow with the addition of reserve requirements
%   specified as a set of fixed zonal reserves. See RUNOPF for a
%   description of the input and output arguments, which are the same,
%   with the exception that the case file or struct CASEDATA must define
%   a 'reserves' field, which is a struct with the following fields:
%       zones   nrz x ng, zone(i, j) = 1, if gen j belongs to zone i
%                                      0, otherwise
%       req     nrz x 1, zonal reserve requirement in MW
%       cost    (ng or ngr) x 1, cost of reserves in $/MW
%       qty     (ng or ngr) x 1, max quantity of reserves in MW (optional)
%   where nrz is the number of reserve zones and ngr is the number of
%   generators belonging to at least one reserve zone and ng is the total
%   number of generators.
%
%   In addition to the normal OPF output, the RESULTS struct contains a
%   new 'reserves' field with the following fields, in addition to those
%   provided in the input:
%       R       - ng x 1, reserves provided by each gen in MW
%       Rmin    - ng x 1, lower limit on reserves provided by each gen, (MW)
%       Rmax    - ng x 1, upper limit on reserves provided by each gen, (MW)
%       mu.l    - ng x 1, shadow price on reserve lower limit, ($/MW)
%       mu.u    - ng x 1, shadow price on reserve upper limit, ($/MW)
%       mu.Pmax - ng x 1, shadow price on Pg + R <= Pmax constraint, ($/MW)
%       prc     - ng x 1, reserve price for each gen equal to maximum of the
%                         shadow prices on the zonal requirement constraint
%                         for each zone the generator belongs to
%
%   See T_CASE30_USERFCNS for an example case file with fixed reserves,
%   and TOGGLE_RESERVES for the implementation.
%
%   Calling syntax options:
%       results = runopf_w_res(casedata);
%       results = runopf_w_res(casedata, mpopt);
%       results = runopf_w_res(casedata, mpopt, fname);
%       results = runopf_w_res(casedata, mpopt, fname, solvedcase);
%       [results, success] = runopf_w_res(...);
%
%   Example:
%       results = runopf_w_res('t_case30_userfcns');
%
%   See also RUNOPF, TOGGLE_RESERVES, T_CASE30_USERFCNS.

%   MATPOWER
%   Copyright (c) 2008-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

mpc = loadcase(varargin{1});
mpc = toggle_reserves(mpc, 'on');
[varargout{1:nargout}] = runopf(mpc, varargin{2:nargin});

if nargout > 0 && isstruct(varargout{1})
    varargout{1} = toggle_reserves(varargout{1}, 'off');
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
echo on
%--------------------------------------------------
% ToUTM
%   Example conversion from Geodetic to UTM
%   coordinates.  Creates UTM coordinate listing
%   in GeoLab 2.0 NE format.
% 05 Feb 2010
%
% M-files:  dms2rad, refell, ell2utm
%--------------------------------------------------
%clear all

%---------- Define names and geodetic lat,lon of
%---------- points to convert

% PLO      07KC005ECC   N 58 54  12.11752 W111 35  21.20267      104.084
% PLO      07NA001ECC   N 58 59  47.66810 W111 23  32.15643      106.236
% PLO      100A         N 58 50  16.59244 W111 41  23.36311      106.060
% PLO      102A         N 58 52  23.70694 W111 46   9.27700      107.385
% NE       07KC005ECC       6529232.4191      466051.9575            UTM 12
% NE       07NA001ECC       6539528.0046      477460.4120            UTM 12
% NE       100A             6522003.1191      460180.9293            UTM 12
% NE       102A             6525984.6308      455641.7913            UTM 12

name=[
'07KC005ECC   '
'07NA001ECC   '
'100A         '
'102A         '
];
plo=[
58 54  12.11752  111 35  21.20267      104.084
58 59  47.66810  111 23  32.15643      106.236
58 50  16.59244  111 41  23.36311      106.060
58 52  23.70694  111 46   9.27700      107.385
];

%---------- Convert to UTM

latdms=plo(:,1:3);
londms=plo(:,4:6);
lat=dms2rad(latdms);
lon=-dms2rad(londms);
n=length(lat);
[a,b,e2,finv]=refell('NAD27');
[N,E,Zone]=ell2utm(lat,lon,a,e2);

%---------- List results

%[N E Zone]
echo off

disp(' ');
s=['          STATION NAME         N (m)',...
   '            E (m)                 Zone'];
disp(s);
s=['------------------------------------',...
   '---------------------------------------'];
disp(s);
for i=1:n
  s=[' NE       ',name(i,:),'%16.4f %16.4f',...
     '            UTM %3.0f\n'];
  fprintf(s,N(i),E(i),Zone(i));
end
fprintf('\n');


text:
function cl = cluster2subclusters(cl_in,class)
% Take a single cluster cl_in and separate into subclusters based on
% vector of integers class
%
% Class must code unique subclusters subcluster order is only preserved
% if class contains all integers from 1 to nclasses:
%
% i.e., class 3 will only be in subcluster 3 if there are no missing class
% numbers in class
%
% ..
%    tor wager, july 06
% ..

classes = unique(class);    % vector of class numbers

nclust = length(classes);

% Fill in fields
for i = 1:nclust

    % voxels in this cluster
    wh_incluster = find(class == classes(i));
    nvox = length(wh_incluster);

    cl(i).title = sprintf('Subcluster of %3.0f vox from %s',nvox, cl_in.title);

    % copy these fields
    for f = {'threshold' 'M' 'dim' 'voxSize'}
        if isfield(cl_in,f{1}), cl(i).(f{1}) = cl_in.(f{1}); end
    end

    % add new values for these fields
    cl(i).name = '';
    cl(i).numVox = nvox;

    % select in_subcl voxels for these fields
    for f = {'XYZ' 'XYZmm' 'Z'}
        cl(i).(f{1}) = cl_in.(f{1})(:,wh_incluster);
    end

    % center of mass
    cl(i).mm_center = center_of_mass(cl(i).XYZmm,cl(i).Z);

end

author: canlab

repo: CanlabCore


text:
% VOCLABELCOLORMAP Creates a label color map such that adjacent indices have different
% colors.  Useful for reading and writing index images which contain large indices,
% by encoding them as RGB images.
%
% CMAP = VOCLABELCOLORMAP(N) creates a label color map with N entries.
function cmap = VOClabelcolormap(N)

if nargin==0
    N=256;
end
cmap = zeros(N,3);
for i=1:N
    id = i-1; r=0;g=0;b=0;
    for j=0:7
        r = bitor(r, bitshift(bitget(id,1),7 - j));
        g = bitor(g, bitshift(bitget(id,2),7 - j));
        b = bitor(b, bitshift(bitget(id,3),7 - j));
        id = bitshift(id,-3);
    end
    cmap(i,1)=r; cmap(i,2)=g; cmap(i,3)=b;
end
cmap = cmap / 255;

author: quantombone

repo: exemplarsvm