text:
function [U,t,t1] = rest_on_ground(V,F,varargin)
  % [U,t] = rest_on_ground(V,F)
  %
  % Given a solid mesh (with a meaningful centroid) rest the object on the
  % ground (z=0). This effectively simulates dropping the object on the ground
  % in a universe with no momentum, or imagine the object coming to rest on the
  % floor of a viscous ocean.
  %
  % Inputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of mesh triangle indices into rows of V
  %   Optional:
  %     'Centroid' followed by center of mass
  % Outputs:
  %   U  #U by 3 list of new mesh vertex positions
  %   t  3-long list of contact points
  %

  % default values
  cen = [];
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Centroid'}, ...
    {'cen'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  warning('Should reduce V,F to convex hull...');

  %save('bad-rest.mat','V','F')
  if isempty(cen)
    cen = centroid(V,F);
  end
  U = [V;cen];
  t = [];
  stable = false;
  cross2 = @(a,b,c) ...
    [a(:,2).*b(:,3)-a(:,3).*b(:,2), ...
     a(:,3).*b(:,1)-a(:,1).*b(:,3), ...
     a(:,1).*b(:,2)-a(:,2).*b(:,1)];

  iter = 0;
  max_iters = 3*size(F,1);
  t1 = [];
  while ~stable
    iter = iter+1;
    if iter > max_iters
      warning('Failed to converge in %d iterations',max_iters);
      break;
    end
    switch(numel(t))
    case 0
      [z,t] = min(U(:,3));
      U(:,3) = U(:,3)-z;
      assert(all(isreal(U(:))));
    case 1
      cen = U(end,:);
      arm = normalizerow(cen-U(t,:));
      up = [0 0 1];
      ax = normalizerow(cross(arm,up));
      bi = cross(up,ax);
      W = (U-U(t,:))*[bi;up]';
      A = atan2(W(:,2),W(:,1));
      Acen = A(end);
      if Acen>pi/2
        A = pi-A;
        ax = -ax;
      end
      A([t end]) = inf;
      A(A<0) = A(A<0)+2*pi;
      [z,s] = min(A);
      R2 = axisangle2matrix(ax,z);

    %clf;
    %hold on;
    %tsh = tsurf(F,U,'CData',A,fphong,falpha(0.1,0),fsoft);
    %ssh =     sct(U(end,:),'o','MarkerFaceColor','r');
    %ssh =     sct(U(t,:),'o','MarkerFaceColor','g');
    %ssh =     sct(U(s,:),'o','MarkerFaceColor','k');
    %%ssh = sct(U(g,:),'o','MarkerFaceColor','r');
    %qvr(U(t,:),100*bi,'Color','r','LineWidth',3);
    %qvr(U(t,:),100*up,'Color','g','LineWidth',3);
    %qvr(U(t,:),100*ax,'Color','b','LineWidth',3);
    %qvr(U(t,:),100*arm,'Color','c','LineWidth',3);
    %%hold off;
    %%l = light('Position',[0 1e-10 1],'Style','infinite');
    %%add_shadow(tsh,l,'Fade','none');
    %%sh = add_shadow([],l);
    %%sh{1}.MarkerFaceColor = ssh.MarkerFaceColor;
    %axis equal;
    %view(0,0);
    %hold off;

      U = (U-U(t,:))*R2+U(t,:);
      t = [t s];

      assert(all(isreal(U(:))));

    case 2
      cen = U(end,:);
      s = t(2);
      ax = normalizerow(U(s,:)-U(t(1),:));
      W = (U-U(s,:))*[normalizerow(cross(up,ax));up]';
      A = atan2(W(:,2),W(:,1));
      Acen = A(end);
      if Acen>pi/2
        A = pi-A;
        ax = -ax;
      end
      A([t end]) = inf;
      % ignore colinear vertices
      D = normrow(cross2(ax,U-U(t(1),:),2));

      A(D<1e-10) = inf;
      [z,g] = min(A);
      R3 = axisangle2matrix(ax,z);
      if all(t == [1779 8746])
    %clf;
    %hold on;
    %tsh = tsurf(F,U,'FaceColor',blue,falpha(0.1,0),fsoft);
    %ssh =     sct(U(end,:),'o','MarkerFaceColor','r');
    %ssh =     sct(U(t,:),'o','MarkerFaceColor','b');
    %ssh =     sct(U(g,:),'o','MarkerFaceColor','g');
    %%ssh = sct(U(g,:),'o','MarkerFaceColor','r');
    %%qvr(U(t,:),100*arm,'Color','r','LineWidth',3);
    %%qvr(U(t,:),100*up,'Color','g','LineWidth',3);
    %%qvr(U(t,:),100*ax,'Color','b','LineWidth',3);
    %%hold off;
    %%l = light('Position',[0 1e-10 1],'Style','infinite');
    %%add_shadow(tsh,l,'Fade','none');
    %%sh = add_shadow([],l);
    %%sh{1}.MarkerFaceColor = ssh.MarkerFaceColor;
    %axis equal;

      end

      U = (U-U(t(1),:))*R3+U(t(1),:);
      t = [t g];
      assert(all(isreal(U(:))));
    case 3
      if isempty(t1)
        t1 = sort(t);
        if (cross(U(t1(2),:)-U(t1(1),:),U(t1(3),:)-U(t1(1),:))*[0;0;1])>0
          t1 = t1([1 3 2]);
        end
      end
      cen = U(end,:);
      % check if stable
      H = convhull(U(t,1:2));
      stable = inpolygon(cen(1),cen(2),U(t(H),1),U(t(H),2));
      if ~stable
        % which points to remove
        E = t([H(1:end-1) H(2:end)]);
        debug = false;

        if debug
          clf;
          hold on;
          tsh = tsurf(F,U,'FaceColor',blue,falpha(0.1,0),fsoft);
          ssh = sct(cen,'o','MarkerFaceColor','r');
          ssh = sct(U(t,:),'o','MarkerFaceColor','b');
          hold off;
          view(2);
        end

        [sqrD,I,C] = point_mesh_squared_distance(cen(1:2),U(:,1:2),E);
        B = normrow(C-U(E(I,2),1:2))/normrow(U(E(I,1),1:2)-U(E(I,2),1:2));
        B = [B 1-B];
        tI = E(I,:);

        if debug
          hold on;
          plot_edges(U,tI,'LineWidth',3);
          txt(U(t,:),num2str(t'))
          hold off;
        end

        if B(1)<1e-8
          t = tI(2);
        elseif B(2)<1e-8
          t = tI(1);
        else
          t = tI;
        end
        if debug
          hold on;
          ssh = sct(U(t,:),'o','MarkerFaceColor','g');
          hold off;
          pause
        end
      end
      assert(all(isreal(U(:))));
    end

    tsh = [];
    %clf;
    %hold on;
    %tsh = tsurf(F,U,'FaceColor',blue,falpha(1.0,0),fsoft);
    %ssh = sct(U(end,:),'o','MarkerFaceColor','r');
    %ssh = sct(U(t,:),'o','MarkerFaceColor','g');
    %%ssh = sct(U(g,:),'o','MarkerFaceColor','r');
    %%qvr(U(t,:),100*arm,'Color','r','LineWidth',3);
    %%qvr(U(t,:),100*up,'Color','g','LineWidth',3);
    %%qvr(U(t,:),100*ax,'Color','b','LineWidth',3);
    %%hold off;
    %%view(2);
    %view(3);
    %l = light('Position',[0 1e-10 1],'Style','infinite');
    %add_shadow(tsh,l,'Fade','none');
    %%sh = add_shadow([],l);
    %%sh{1}.MarkerFaceColor = ssh.MarkerFaceColor;
    %axis equal;
    %camproj('persp');
    %drawnow;
    if stable
      t = sort(t);
      if (cross(U(t(2),:)-U(t(1),:),U(t(3),:)-U(t(1),:))*[0;0;1])>0
        t = t([1 3 2]);
      end
    end
  end
  % Strip off centroid
  U = U(1:end-1,:);

end

author: alecjacobson

repo: gptoolbox


text:
% util_checkCloseMatrixElements Check if some rows of the matrix A are too
% "similar". Similar means their Euclidean distance being smaller than
% tolerance, which if not given is set automatically heuristically
function [res, distances, tolerance] = util_checkCloseMatrixElements(A, tolerance)

C = nan(size(A,1));
for i = 1:size(A,1)
    for j = i+1:size(A,1)
        C(i,j) = dist2(A(i,:), A(j,:));%/norm(abs(A(i,:))+abs(A(j,:)));
    end
end
nn = C(~isnan(C));
distances = C;
C = C./mean(nn);
if nargin < 2 || isempty(tolerance)
    tolerance = 0.0001;
end

%maxDist = dist2(max(A), min(A));
%C = dist2(A, A);
%C = C + eye(size(C))*maxDist*2; % Don't care about diagonal (same point)
%if nargin < 2 || isempty(tolerance)
%    tolerance = (maxDist/mean(sum(A,2).^2))*0.0001;
%    %[~,j]=sort(sum(A,2).^2);
%    %A = A(j,:);
%end
tmp = find(C <= tolerance);
res = [mod(tmp, size(C,1)) ceil(tmp/size(C,1))];
res(res==0) = size(C,1);
res = unique(sort(res,2), 'rows');

author: SheffieldML

repo: vargplvm


text:
function [cnmfAnalysisOutput] = computeCnmfSignalExtractionOriginal(inputMovie,numExpectedComponents,varargin)
	% Wrapper for CNMF, for use with https://github.com/epnev/ca_source_extraction/commit/8799b13df2b09f30e27fc852e4f5f39ae6f44405
	% Building off of demo_script.m in CNMF github repo
	% Biafra Ahanonu
	% started: 2016.01.20
	% inputs
		% inputMovie - a string or a cell array of strings pointing to the movies to be analyzed (recommended). Else, [x y t] matrix where t = frames.
		% numExpectedComponents - number of expected components
	% outputs
		% cnmfAnalysisOutput - structure containing extractedImages and extractedSignals along with input parameters to the algorithm
	% READ BEFORE RUNNING
		% Get CVX from http://cvxr.com/cvx/doc/install.html
		% Run the below commands in Matlab after unzipping
		% cvx_setup
		% cvx_save_prefs (permanently stores settings)

	[cnmfAnalysisOutput] = ciapkg.signal_extraction.computeCnmfSignalExtractionOriginal(inputMovie,numExpectedComponents,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [newIndex] = cycle(index, cycleLength)
%function [newIndex] = cycle(index, cycleLength)
%function to facilitate indexing on cyclical structures.
%The new index if the index cycled around on a structure of length cycleLength.
%e.g.: cycle(11, 10) = 1, cycle(12, 10) = 2, cycle(9, 10) = 9, etc.
%'index' can be a vector.
%JOD, Feb 02.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

newIndex = zeros(size(index));

for i = 1 : length(index)
  if index(i) > cycleLength
    newIndex(i) = index(i) - cycleLength;
  elseif index(i) < 1
    newIndex(i) = cycleLength + index(i);
  else
    newIndex(i) = index(i);
  end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function  [energyVal] =  energyVADVal(frame_x, Fs)
energyVal = norm(frame_x);

author: singaxiong

repo: SignalGraph


text:
function props = propnames(s)

if isstruct(s)
    props = fieldnames(s);
else
    props = properties(s);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Colormaps that are pretty awesome: DEMO1
% this is an exampe of 4 colormaps that are being considered as the default
% colormap in python's matplotlib lybrary.
%
% All of them look quite good and they don't have any official name, so at
% the moment they are A,B,C,D.
%
% colormaps from https://github.com/bids/colormap
%
% Ander Biguri
%% Clear workspace and get screen data
clear;
clc
close all;

screen=get(0,'ScreenSize') ;
w0=screen(1);
h0=screen(2);
w =screen(3);
h =screen(4);
%% Generate sample data
load flujet
X=X.';

%% Plot original with jet and parula
% Parula
h2=figure('name','Sample data with "parula" colormap');
set(h2,'position',[w0,h0,w/3,h/2])
imagesc(X)

if verLessThan('matlab', '8.4')
   % if parula is the "future"
   colormap(parula());
else
   % if parula is already in Matlab
   colormap('parula');
end
axis image
xlabel('Parula Colormap')
set(gca,'xtick',[],'ytick',[]) % This is axis off without offing the labels

% Jet
h1=figure('name','Sample data with "jet" colormap');
set(h1,'position',[w0,h/2,w/3,h/2])
imagesc(X)
colormap('jet')
xlabel('jet Colormap')
set(gca,'xtick',[],'ytick',[]) % This is axis off without offing the labels
axis image

%% Load new colormaps

m=100;
cm_magma=magma(m);
cm_inferno=inferno(m);
cm_plasma=plasma(m);
cm_viridis=viridis(m);

%% Plot new colormaps
h3=figure('name','Super-cool new colormaps that you can easily use');
set(h3,'position',[w/3,h0,2*w/3,h])
if verLessThan('matlab', '8.4')
% If you are using old Matlab figure engine do  it this way
% (some very lousy colormap problems before)
    subplot(2,2,1,'Position',[0.05 0.55 0.4 0.4])
    subimage(uint8(X/max(X(:))*255),cm_magma)
    xlabel('MAGMA')
    set(gca,'xtick',[],'ytick',[])

    subplot(2,2,2,'Position',[0.55 0.55 0.4 0.4])
    subimage(uint8(X/max(X(:))*255),cm_inferno)
    xlabel('INFERNO')
    set(gca,'xtick',[],'ytick',[])

    subplot(2,2,3,'Position',[0.05 0.05 0.4 0.4])
    subimage(uint8(X/max(X(:))*255),cm_plasma)
    xlabel('PLASMA')
    set(gca,'xtick',[],'ytick',[])

    subplot(2,2,4,'Position',[0.55 0.05 0.4 0.4])
    subimage(uint8(X/max(X(:))*255),cm_viridis)
    xlabel('VIRIDIS')
    set(gca,'xtick',[],'ytick',[])
else
    sp1=subplot(2,2,1,'Position',[0.05 0.55 0.4 0.4]);
    imagesc(X)
    colormap(sp1,cm_magma)
    xlabel('MAGMA')
    set(gca,'xtick',[],'ytick',[])

    sp2=subplot(2,2,2,'Position',[0.55 0.55 0.4 0.4]);
    imagesc(X)
    colormap(sp2,cm_inferno)
    xlabel('INFERNO')
    set(gca,'xtick',[],'ytick',[])

    sp3=subplot(2,2,3,'Position',[0.05 0.05 0.4 0.4]);
    imagesc(X)
    colormap(sp3,cm_plasma)
    xlabel('PLASMA')
    set(gca,'xtick',[],'ytick',[])

    sp4=subplot(2,2,4,'Position',[0.55 0.05 0.4 0.4]);
    imagesc(X)
    colormap(sp4,cm_viridis)
    xlabel('VIRIDIS')
    set(gca,'xtick',[],'ytick',[])
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
% plot_parameters: script to set globally plot parameters
%
% INPUT
%   None.
%
% OUTPUT
%   None.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Version: 001
% Date:    2021/12/14
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

%% COLORS

blue    = [0, 0.4470, 0.7410];
orange  = [0.8500, 0.3250, 0.0980];
green   = [0.4660, 0.6740, 0.1880];
yellow  = [0.9290, 0.6940, 0.1250];
light_blue = [0.3010, 0.7450, 0.9330];

%% TEXT SIZES

font_title  = 25;
font_label  = 15;
font_tick   = 10;
font_legend = 10;

%% LINE WIDTH

lw = 1.75;

%% MARKER SIZE

ms = 10;

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
%
% Initialization for AMICO
%
% NB: DO NOT MODIFY THIS FILE!
%     Make a copy of it, adapt to your paths and rename it to "AMICO_Setup.m"
%

global AMICO_code_path AMICO_data_path

% Path definition: adapt these to your needs
% ==========================================
AMICO_code_path = fileparts(mfilename('fullpath'));
AMICO_data_path = fullfile(AMICO_code_path,'exports');
if(~exist(AMICO_data_path,'dir'))
    mkdir(AMICO_data_path);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function img = plimage2double( plimg )
%PLIMAGE2DOUBLE Summary of this function goes here
%   Detailed explanation goes here
img = reshape(plimg.toDouble(), [plimg.mRow, plimg.mCol, plimg.mLen]);

end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function save_nii_hdr(hdr, fid)

   if ~isequal(hdr.hk.sizeof_hdr,348),
      error('hdr.hk.sizeof_hdr must be 348.');
   end

   write_header(hdr, fid);

   return;					% save_nii_hdr

%---------------------------------------------------------------------
function write_header(hdr, fid)

        %  Original header structures
	%  struct dsr				/* dsr = hdr */
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

   header_key(fid, hdr.hk);
   image_dimension(fid, hdr.dime);
   data_history(fid, hdr.hist);

   %  check the file size is 348 bytes
   %
   fbytes = ftell(fid);

   if ~isequal(fbytes,348),
      msg = sprintf('Header size is not 348 bytes.');
      warning(msg);
   end

   return;					% write_header

%---------------------------------------------------------------------
function header_key(fid, hk)

   fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                      /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */

   fwrite(fid, hk.sizeof_hdr(1),    'int32');	% must be 348.

   % data_type = sprintf('%-10s',hk.data_type);	% ensure it is 10 chars from left
   % fwrite(fid, data_type(1:10), 'uchar');
   pad = zeros(1, 10-length(hk.data_type));
   hk.data_type = [hk.data_type  char(pad)];
   fwrite(fid, hk.data_type(1:10), 'uchar');

   % db_name   = sprintf('%-18s', hk.db_name);	% ensure it is 18 chars from left
   % fwrite(fid, db_name(1:18), 'uchar');
   pad = zeros(1, 18-length(hk.db_name));
   hk.db_name = [hk.db_name  char(pad)];
   fwrite(fid, hk.db_name(1:18), 'uchar');

   fwrite(fid, hk.extents(1),       'int32');
   fwrite(fid, hk.session_error(1), 'int16');
   fwrite(fid, hk.regular(1),       'uchar');	% might be uint8

   % fwrite(fid, hk.hkey_un0(1),    'uchar');
   % fwrite(fid, hk.hkey_un0(1),    'uint8');
   fwrite(fid, hk.dim_info(1),      'uchar');

   return;					% header_key

%---------------------------------------------------------------------
function image_dimension(fid, dime)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%			/*
	%				pixdim[] specifies the voxel dimensions:
	%				pixdim[1] - voxel width
	%				pixdim[2] - voxel height
	%				pixdim[3] - interslice distance
	%				pixdim[4] - volume timing, in msec
	%					..etc
	%			*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

   fwrite(fid, dime.dim(1:8),        'int16');
   fwrite(fid, dime.intent_p1(1),  'float32');
   fwrite(fid, dime.intent_p2(1),  'float32');
   fwrite(fid, dime.intent_p3(1),  'float32');
   fwrite(fid, dime.intent_code(1),  'int16');
   fwrite(fid, dime.datatype(1),     'int16');
   fwrite(fid, dime.bitpix(1),       'int16');
   fwrite(fid, dime.slice_start(1),  'int16');
   fwrite(fid, dime.pixdim(1:8),   'float32');
   fwrite(fid, dime.vox_offset(1), 'float32');
   fwrite(fid, dime.scl_slope(1),  'float32');
   fwrite(fid, dime.scl_inter(1),  'float32');
   fwrite(fid, dime.slice_end(1),    'int16');
   fwrite(fid, dime.slice_code(1),   'uchar');
   fwrite(fid, dime.xyzt_units(1),   'uchar');
   fwrite(fid, dime.cal_max(1),    'float32');
   fwrite(fid, dime.cal_min(1),    'float32');
   fwrite(fid, dime.slice_duration(1), 'float32');
   fwrite(fid, dime.toffset(1),    'float32');
   fwrite(fid, dime.glmax(1),        'int32');
   fwrite(fid, dime.glmin(1),        'int32');

   return;					% image_dimension

%---------------------------------------------------------------------
function data_history(fid, hist)

	% Original header structures
	%struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

   % descrip     = sprintf('%-80s', hist.descrip);     % 80 chars from left
   % fwrite(fid, descrip(1:80),    'uchar');
   pad = zeros(1, 80-length(hist.descrip));
   hist.descrip = [hist.descrip  char(pad)];
   fwrite(fid, hist.descrip(1:80), 'uchar');

   % aux_file    = sprintf('%-24s', hist.aux_file);    % 24 chars from left
   % fwrite(fid, aux_file(1:24),   'uchar');
   pad = zeros(1, 24-length(hist.aux_file));
   hist.aux_file = [hist.aux_file  char(pad)];
   fwrite(fid, hist.aux_file(1:24), 'uchar');

   fwrite(fid, hist.qform_code,    'int16');
   fwrite(fid, hist.sform_code,    'int16');
   fwrite(fid, hist.quatern_b,   'float32');
   fwrite(fid, hist.quatern_c,   'float32');
   fwrite(fid, hist.quatern_d,   'float32');
   fwrite(fid, hist.qoffset_x,   'float32');
   fwrite(fid, hist.qoffset_y,   'float32');
   fwrite(fid, hist.qoffset_z,   'float32');
   fwrite(fid, hist.srow_x(1:4), 'float32');
   fwrite(fid, hist.srow_y(1:4), 'float32');
   fwrite(fid, hist.srow_z(1:4), 'float32');

   % intent_name = sprintf('%-16s', hist.intent_name);	% 16 chars from left
   % fwrite(fid, intent_name(1:16),    'uchar');
   pad = zeros(1, 16-length(hist.intent_name));
   hist.intent_name = [hist.intent_name  char(pad)];
   fwrite(fid, hist.intent_name(1:16), 'uchar');

   % magic	= sprintf('%-4s', hist.magic);		% 4 chars from left
   % fwrite(fid, magic(1:4),           'uchar');
   pad = zeros(1, 4-length(hist.magic));
   hist.magic = [hist.magic  char(pad)];
   fwrite(fid, hist.magic(1:4),        'uchar');

   return;					% data_history

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function data=readbin(fname,nrow,type)
%Open the file
fid=fopen(fname,'rb');
if (fid==-1)
    disp('file not found');
    data=[];
    return;
end

%Determine the type
if (nargin==3)
    if (isempty(type))
        type='double';
    end
else
    type='double';
end

%Read the Data
data=fread(fid,[nrow inf],type);

%Check the last data to see if the nrow is true
sr_a=zeros(1,type);
sr_b=whos('sr_a');
nbyte=sr_b.bytes;

fseek(fid,-nbyte, 'cof');
sr_a=fread(fid,1,type);
if (sr_a~=data(end))
    disp('Achtung: (readbin) "nrow" is probably wrong!');
end

fclose(fid);

repo: ekfmonoslam


text:
function dataType = tapas_uniqc_get_data_type_from_n_voxels(nVoxels)
% Specifies memory-efficient data type for saving from number of voxels
% TODO: incorporate dynamic range of data to be saved as well! add option
% for user to specify data type and bit depth?
%
% dataType = tapas_uniqc_get_data_type_from_n_voxels(nVoxels)
%
% IN
%   nVoxels    [1,n] voxels per dimension to be saved
%               see also MrImageGeometry
% OUT
%
% EXAMPLE
%   tapas_uniqc_get_data_type_from_n_voxels
%
%   See also

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2015-12-06
% Copyright (C) 2015 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% - use double precion (64 bit) for structural images (3D and more than
%   220x220x120 voxel)
% - use int16 for very large data sets (when float32 would exceed 2GB)
% - use single precision (32 bit) for everything in-between

is3D = numel(nVoxels) <= 3 || nVoxels(4) == 1;
isStructural = prod(nVoxels(1:3)) >= 220*220*120;
floatExceeds2GB = prod(nVoxels) >= 2*1024*1024*1024*8/32;

if is3D && isStructural % highest bit resolution for structural images
    dataType   = 'float64';
elseif floatExceeds2GB % int16 for large raw data
    dataType = 'int16';
    warning(['Due to the large number of samples, ', ...
        'data will be converted and saved as 16-bit integer.']);
else
    dataType   = 'float32'; %float32 for everything in between
end

author: translationalneuromodeling

repo: tapas


text:
function S = ft_omri_info_from_header(hdr)

% function S = ft_omri_info_from_header(hdr)
%
% Convenience function to retrieve most important MR information
% from a given header (H) as retrieved from a FieldTrip buffer.
% Will look at both NIFTI-1 and SiemensAP fields, if present, and
% give preference to SiemensAP info.
%
% Returns empty array if no information could be found.

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

SNif = [];
SSap = [];
if isfield(hdr,'nifti_1')
	try
		SNif = mri_info_from_nifti(hdr.nifti_1);
	catch
		warning('Errors occured while inspecting NIFTI-1 header.');
	end
end
if isfield(hdr,'siemensap')
	try
		SSap = mri_info_from_sap(hdr.siemensap);
	catch
		warning('Errors occured while inspecting SiemensAP header.');
	end
end

if ~isempty(SSap)
	S = SSap;
	if ~isempty(SNif)
		if ~isequal(SNif.voxels,SSap.voxels)
			warning('Conflicting information in NIFTI and SiemensAP - trusting SiemensAP...');
    end
		S.mat0 = SNif.mat0;
	end
else
	if ~isempty(SNif)
		S = SNif;
	else
		S = [];
	end
end

function S = mri_info_from_nifti(NH)

S.vx = double(NH.dim(1));
S.vy = double(NH.dim(2));
S.vz = double(NH.dim(3));
S.voxels = [S.vx S.vy S.vz];
S.voxdim = double(NH.pixdim(1:3));
S.size = S.voxels .* S.voxdim;
VoxToWorld = double([NH.srow_x; NH.srow_y; NH.srow_z]);
M = VoxToWorld(1:3,1:3);
P = VoxToWorld(1:3,4);
% correct Mat0 in the same way SPM does (voxel index starts at 1)
S.mat0 = [M (P-M*[1;1;1]); 0 0 0 1];
S.numEchos = 1;	% can't detect this from NIFTI :-(

switch NH.slice_code
	% Long-term TODO: look at slice_start and slice_end for padded slices
	case 1	% NIFTI_SLICE_SEQ_INC
		inds = 1:S.vz;
	case 2  % NIFTI_SLICE_SEQ_DEC
		inds = S.vz:-1:1;
	case 3  % NIFTI_SLICE_ALT_INC
		inds = [(1:2:S.vz) (2:2:S.vz)];
	case 4  % NIFTI_SLICE_ALT_DEC
		inds = [(S.vz:-2:1) ((S.vz-1):-2:1)];
	case 5  % NIFTI_SLICE_ALT_INC2
		inds = [(2:2:S.vz) (1:2:S.vz)];
	case 6  % NIFTI_SLICE_ALT_DEC2
		inds = [((S.vz-1):-2:1) (S.vz:-2:1)];
	otherwise
		warning('Unrecognized slice order - using default');
		inds = 1:S.vz;
end

if NH.slice_duration > 0
	S.TR = double(NH.slice_duration * S.vz);
	% first set up linear
	S.deltaT = (0:(S.vz-1))*double(NH.slice_duration)
	% then re-shuffle
	S.deltaT(inds) = S.deltaT;
else
	% what can we do here?
	S.TR = 2;
	S.deltaT = (0:(S.vz-1))*S.TR/S.vz;
	S.deltaT(inds) = S.deltaT;
end

function S = mri_info_from_sap(SP)

phaseFOV   = SP.sSliceArray.asSlice{1}.dPhaseFOV;
readoutFOV = SP.sSliceArray.asSlice{1}.dReadoutFOV;
sliceThick = SP.sSliceArray.asSlice{1}.dThickness;
distFactor = SP.sGroupArray.asGroup{1}.dDistFact;

S.vx = double(SP.sKSpace.lBaseResolution);
S.vy = S.vx * phaseFOV / readoutFOV;
S.vz = double(SP.sSliceArray.lSize);
S.voxels = [S.vx S.vy S.vz];

% this only takes care of the scaling, not the proper orientation
sx = readoutFOV/S.vx;
sy = phaseFOV/S.vy;	% should always be == sx
sz = sliceThick * (1.0 + distFactor);

S.size = [readoutFOV  phaseFOV  S.vz*sz];
S.mat0 = [sx 0 0 0;  0 sy 0 0; 0 0 sz 0; 0 0 0 1];
S.voxdim = [sx sy sz];

S.numEchos = double(SP.lContrasts);
S.TR = double(SP.alTR) * 1e-6;  % originally in microseconds

switch SP.sSliceArray.ucMode
	case 1	% == NIFTI_SLICE_SEQ_INC
		inds = 1:S.vz;
	case 2  % == NIFTI_SLICE_SEQ_DEC
		inds = S.vz:-1:1;
	case 4  % odd:ALT_INC or even:ALT_INC2
		if mod(S.vz,2) == 1
			inds = [(1:2:S.vz) (2:2:S.vz)];
		else
			inds = [(2:2:S.vz) (1:2:S.vz)];
		end
	otherwise
		warning('Unrecognized slice order - using default');
		inds = 1:S.vz;
end
% first set up linear
S.deltaT = (0:(S.vz-1))*S.TR/S.vz
% then re-shuffle
S.deltaT(inds) = S.deltaT;


text:
%% PTplotSpec - script that computes and plots spectrograms

% ----------------------------------------------------------------------------------
% "THE BEER-WARE LICENSE" (Revision 42):
% <brian.white@queensu.ca> wrote this file. As long as you retain this notice you
% can do whatever you want with this stuff. If we meet some day, and you think
% this stuff is worth it, you can buy me a beer in return. -Brian White
% ----------------------------------------------------------------------------------

if ~isempty(fnameMaster)

prop_max_screen=(max([PTspecfig.Position(3) PTspecfig.Position(4)]));
fontsz=(screensz_multiplier*prop_max_screen);
%% update fonts

f = fields(guiHandlesSpec);
for i = 1 : size(f,1)
    try
        eval(['guiHandlesSpec.' f{i} '.FontSize=fontsz;']);
    catch
        guiHandlesSpec.SpecSelect{1}.FontSize=fontsz;
        guiHandlesSpec.SpecSelect{2}.FontSize=fontsz;
        guiHandlesSpec.SpecSelect{3}.FontSize=fontsz;
        guiHandlesSpec.SpecSelect{4}.FontSize=fontsz;
        guiHandlesSpec.FileSelect{1}.FontSize=fontsz;
        guiHandlesSpec.FileSelect{2}.FontSize=fontsz;
        guiHandlesSpec.FileSelect{3}.FontSize=fontsz;
        guiHandlesSpec.FileSelect{4}.FontSize=fontsz;
        guiHandlesSpec.Sub100HzCheck{1}.FontSize=fontsz;
        guiHandlesSpec.Sub100HzCheck{2}.FontSize=fontsz;
        guiHandlesSpec.Sub100HzCheck{3}.FontSize=fontsz;
        guiHandlesSpec.Sub100HzCheck{4}.FontSize=fontsz;
    end
end

guiHandlesSpec.climMax_input = uicontrol(PTspecfig,'style','edit','string',[num2str(climScale(guiHandlesSpec.checkboxPSD.Value+1, 1))],'fontsize',fontsz,'TooltipString',[TooltipString_scale],'units','normalized','outerposition',[posInfo.climMax_input],...
     'callback','@textinput_call2; climScale(guiHandlesSpec.checkboxPSD.Value+1, 1)=str2num(guiHandlesSpec.climMax_input.String);updateSpec=1;PTplotSpec;');

guiHandlesSpec.climMax_input2 = uicontrol(PTspecfig,'style','edit','string',[num2str(climScale(guiHandlesSpec.checkboxPSD.Value+1, 2))],'fontsize',fontsz,'TooltipString',[TooltipString_scale],'units','normalized','outerposition',[posInfo.climMax_input2],...
     'callback','@textinput_call2; climScale(guiHandlesSpec.checkboxPSD.Value+1, 2)=str2num(guiHandlesSpec.climMax_input2.String);updateSpec=1;PTplotSpec;');

guiHandlesSpec.climMax_input3 = uicontrol(PTspecfig,'style','edit','string',[num2str(climScale(guiHandlesSpec.checkboxPSD.Value+1, 3))],'fontsize',fontsz,'TooltipString',[TooltipString_scale],'units','normalized','outerposition',[posInfo.climMax_input3],...
     'callback','@textinput_call2; climScale(guiHandlesSpec.checkboxPSD.Value+1, 3)=str2num(guiHandlesSpec.climMax_input3.String);updateSpec=1;PTplotSpec;');

guiHandlesSpec.climMax_input4 = uicontrol(PTspecfig,'style','edit','string',[num2str(climScale(guiHandlesSpec.checkboxPSD.Value+1, 4))],'fontsize',fontsz,'TooltipString',[TooltipString_scale],'units','normalized','outerposition',[posInfo.climMax_input4],...
     'callback','@textinput_call2; climScale(guiHandlesSpec.checkboxPSD.Value+1, 4)=str2num(guiHandlesSpec.climMax_input4.String);updateSpec=1;PTplotSpec;');

%%

s1={'';'gyroADC';'debug';'piderr';'setpoint';'axisP';'axisD';'axisDpf';'pidsum'};

datSelectionString=[s1];

clear vars
for i=1:4
    vars(i)=guiHandlesSpec.SpecSelect{i}.Value;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% compute fft %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if guiHandlesSpec.SpecSelect{1}.Value>1 | guiHandlesSpec.SpecSelect{2}.Value>1 | guiHandlesSpec.SpecSelect{3}.Value>1 | guiHandlesSpec.SpecSelect{4}.Value>1
    set(PTspecfig, 'pointer', 'watch')
    if updateSpec==0
        clear s dat ampmat amp2d freq a RC smat amp2d freq2d Throt
        p=0;
         hw = waitbar(0,['please wait... ' ]);

        for k=1:length(vars)
            s=char(datSelectionString(vars(k)));
            for a=1:3,
                if  ( ( ~isempty(strfind(s,'axisD'))) & a==3) | isempty(s)
                    p=p+1;
                    smat{p}=[];%string
                    ampmat{p}=[];%spec matrix
                    freq{p}=[];% freq matrix
                    amp2d{p}=[];%spec 2d
                    freq2d{p}=[];% freq2d
                else
                    eval(['dat{k}(a,:) = T{guiHandlesSpec.FileSelect{k}.Value}.' char(datSelectionString(vars(k))) '_' int2str(a-1) '_(tIND{guiHandlesSpec.FileSelect{k}.Value});';])
                    Throt=T{guiHandlesSpec.FileSelect{k}.Value}.setpoint_3_(tIND{guiHandlesSpec.FileSelect{k}.Value}) / 10;% throttle
                    lograte = A_lograte(guiHandlesSpec.FileSelect{k}.Value);%in kHz
                    p=p+1;
                    waitbar(p/12, hw, ['processing spectrogram... '  int2str(p) ]);
                    smat{p}=s;
                    [freq{p} ampmat{p}]=PTthrSpec(Throt, dat{k}(a,:), lograte, guiHandlesSpec.checkboxPSD.Value); % compute matrices
                    [freq2d{p} amp2d{p}]=PTSpec2d(dat{k}(a,:),lograte, guiHandlesSpec.checkboxPSD.Value); %compute 2d amp spec at same time
               end
            end
        end
        close(hw)
    end
else
    hwarn=warndlg({'Dropdowns set to ''NONE''.'; 'Please select a preset or specific variables to analyze.'});
    pause(3);
    try
        close(hwarn);
    catch
    end
end

if guiHandlesSpec.checkbox2d.Value==0 && ~isempty(ampmat)
    figure(PTspecfig);
    %%%%% plot spec mattrices
    c1=[1 1 1 2 2 2 3 3 3 4 4 4];
    c2=[1 2 3 1 2 3 1 2 3 1 2 3];
    baselineY = [0 -40];
    ftr = fspecial('gaussian',[guiHandlesSpec.smoothFactor_select.Value*5 guiHandlesSpec.smoothFactor_select.Value],4);
    for p=1:size(ampmat,2)
        delete(subplot('position',posInfo.SpecPos(p,:)));
        if ~isempty(ampmat{p})
            delete(subplot('position',posInfo.SpecPos(p,:)));
            h1=subplot('position',posInfo.SpecPos(p,:)); cla
            img = flipud((filter2(ftr, ampmat{p} ))') + baselineY(guiHandlesSpec.checkboxPSD.Value+1);
            imagesc(img);

            lograte=A_lograte(guiHandlesSpec.FileSelect{c1(p)}.Value);

             axLabel={'roll';'pitch';'yaw'};

            if guiHandlesSpec.Sub100HzCheck{c1(p)}.Value==1
                hold on;h=plot([0 100],[size(ampmat{p},2)-round(Flim1/3.33) size(ampmat{p},2)-round(Flim1/3.33)],'y--');set(h,'linewidth',2)
                hold on;h=plot([0 100],[size(ampmat{p},2)-round(Flim2/3.33) size(ampmat{p},2)-round(Flim2/3.33)],'y--');set(h,'linewidth',2)
                % sub100Hz scaling
                xticks=[1 size(ampmat{p},1)/5:size(ampmat{p},1)/5:size(ampmat{p},1)];
                yticks=[(size(ampmat{p},2)-30):6:size(ampmat{p},2)];
                set(h1,'PlotBoxAspectRatioMode','auto','ylim',[size(ampmat{p},2)-30 size(ampmat{p},2)])
                set(h1,'fontsize',fontsz,'CLim',[baselineY(guiHandlesSpec.checkboxPSD.Value+1) climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))],'YTick',[yticks],'yticklabels',[{100} {80} {60} {40} {20} {0}],'XTick',[xticks],'xticklabels',{'0';'20';'40';'60';'80';'100'},'tickdir','out','xminortick','on','yminortick','on');
                a=[];a2=[];a=filter2(ftr, ampmat{p}) + baselineY(guiHandlesSpec.checkboxPSD.Value+1);
                a2 = a(:,(round(Flim1/3.33))+1:(round(Flim2/3.33)));
                meanspec=nanmean(a2(:));
                peakspec=max(max(a(:,(round(Flim1/3.33))+1:(round(Flim2/3.33)))));
                if guiHandlesSpec.ColormapSelect.Value==9 | guiHandlesSpec.ColormapSelect.Value==10
                    h=text(64,(size(ampmat{p},2)-30)+3,['mean=' num2str(meanspec,3)]);
                    set(h,'Color','k','fontsize',fontsz,'fontweight','bold');
                    h=text(64,(size(ampmat{p},2)-30)+1,['peak=' num2str(peakspec,3)]);
                    set(h,'Color','k','fontsize',fontsz,'fontweight','bold');
                else
                    h=text(64,(size(ampmat{p},2)-30)+3,['mean=' num2str(meanspec,3)]);
                    set(h,'Color','w','fontsize',fontsz,'fontweight','bold');
                    h=text(64,(size(ampmat{p},2)-30)+1,['peak=' num2str(peakspec,3)]);
                    set(h,'Color','w','fontsize',fontsz,'fontweight','bold');
                end
                h=text(xticks(1)+1,(size(ampmat{p},2)-30)+1,axLabel{c2(p)});
                set(h,'Color',[1 1 1],'fontsize',fontsz,'fontweight','bold')
            else % full scaling
                xticks=[1 size(ampmat{p},1)/5:size(ampmat{p},1)/5:size(ampmat{p},1)];
                yticks=[1:(size(ampmat{p},2))/10:size(ampmat{p},2) size(ampmat{p},2)];
                set(h1,'fontsize',fontsz,'CLim',[baselineY(guiHandlesSpec.checkboxPSD.Value+1) climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))],'YTick',[yticks],'yticklabels',[{(max(round(yticks * 3.333))/1)-(max(round(yticks * 3.333))/5)*0} {''} {(max(round(yticks * 3.333))/1)-(max(round(yticks * 3.333))/5)*1} {''} {(max(round(yticks * 3.333))/1)-(max(round(yticks * 3.333))/5)*2} {''} {(max(round(yticks * 3.333))/1)-(max(round(yticks * 3.333))/5)*3} {''} {(max(round(yticks * 3.333))/1)-(max(round(yticks * 3.333))/5)*4} {''} {0}],'XTick',[xticks],'xticklabels',{'0';'20';'40';'60';'80';'100'},'tickdir','out','xminortick','on','yminortick','on');
                set(h1,'PlotBoxAspectRatioMode','auto','ylim',[1 size(ampmat{p},2)])
                a=[];a2=[];a=filter2(ftr, ampmat{p}) + baselineY(guiHandlesSpec.checkboxPSD.Value+1);
                a2 = a(:,(size(ampmat{p},2)/10):size(ampmat{p},2));
                meanspec=nanmean(a2(:));
                peakspec=max(max(a(:,(size(ampmat{p},2)/10):size(ampmat{p},2))));
                if guiHandlesSpec.ColormapSelect.Value==9 | guiHandlesSpec.ColormapSelect.Value==10
                    h=text(64,size(ampmat{p},2)*.04,['mean=' num2str(meanspec,3)]);
                    set(h,'Color','k','fontsize',fontsz,'fontweight','bold');
                    h=text(64,size(ampmat{p},2)*.13,['peak=' num2str(peakspec,3)]);
                    set(h,'Color','k','fontsize',fontsz,'fontweight','bold');
                else
                    h=text(64,size(ampmat{p},2)*.04,['mean=' num2str(meanspec,3)]);
                    set(h,'Color','w','fontsize',fontsz,'fontweight','bold');
                    h=text(64,size(ampmat{p},2)*.13,['peak=' num2str(peakspec,3)]);
                    set(h,'Color','w','fontsize',fontsz,'fontweight','bold');
                end
                h=text(xticks(1)+1,size(ampmat{p},2)*.04,axLabel{c2(p)});
                set(h,'Color',[1 1 1],'fontsize',fontsz,'fontweight','bold')
            end

            grid on
            ax = gca;
            ax.GridColor = [1 1 1];
            if guiHandlesSpec.ColormapSelect.Value==9 | guiHandlesSpec.ColormapSelect.Value==10
                ax.GridColor = [0 0 0]; % black on white background
                set(h,'Color',[0 0 0],'fontsize',fontsz,'fontweight','bold')
            end
             ylabel('Frequency (Hz)','fontweight','bold')
             xlabel('% Throttle','fontweight','bold')
        end
    end

    % color bar2 at the top
    try
    delete(hCbar1);delete(hCbar2);delete(hCbar3);delete(hCbar4)
    catch
    end
    if vars(1)>1 % 1=none
        subplot('position',posInfo.SpecPos(1,:));
        hCbar1= colorbar('NorthOutside');
        set(hCbar1,'Position', [posInfo.hCbar1pos]);
    end
    if vars(2)>1 % 1=none
        subplot('position',posInfo.SpecPos(4,:));
        hCbar2= colorbar('NorthOutside');
        set(hCbar2,'Position', [posInfo.hCbar2pos])
    end
    if vars(3)>1 % 1=none
        subplot('position',posInfo.SpecPos(7,:));
        hCbar3= colorbar('NorthOutside');
        set(hCbar3,'Position', [posInfo.hCbar3pos])
    end
    if vars(4)>1 % 1=none
        subplot('position',posInfo.SpecPos(10,:));
        hCbar4= colorbar('NorthOutside');
        set(hCbar4,'Position', [posInfo.hCbar4pos])
    end

    % color maps
    % standard set
    if guiHandlesSpec.ColormapSelect.Value<8,
        colormap(char(guiHandlesSpec.ColormapSelect.String(guiHandlesSpec.ColormapSelect.Value)));
    end
    % new
    if guiHandlesSpec.ColormapSelect.Value==8, colormap(viridis); end
    if guiHandlesSpec.ColormapSelect.Value==9, colormap(linearREDcmap); end
    if guiHandlesSpec.ColormapSelect.Value==10, colormap(linearGREYcmap); end

end

if guiHandlesSpec.checkbox2d.Value==1 && ~isempty(amp2d)
    figure(PTspecfig);
    try
    delete(hCbar1);delete(hCbar2);delete(hCbar3);delete(hCbar4)
    catch
    end
    baselineYlines = [0 -50];
    c1=[1 1 1 2 2 2 3 3 3 4 4 4];
    c2=[1 2 3 1 2 3 1 2 3 1 2 3];
    %%%%% plot 2d amp spec
    for p=1:size(amp2d,2)
         axLabel={'roll';'pitch';'yaw'};

        delete(subplot('position',posInfo.SpecPos(p,:)));
        if ~isempty(amp2d{p})
            h2=subplot('position',posInfo.SpecPos(p,:)); cla
            h=plot(freq2d{p}, smooth(amp2d{p}, log10(size(amp2d{p},1)) * (guiHandlesSpec.smoothFactor_select.Value^2), 'lowess'));hold on
            set(h, 'linewidth', guiHandles.linewidth.Value/2)
            set(h2,'fontsize',fontsz,'fontweight','bold')
            if guiHandlesSpec.specPresets.Value <= 3
                set(h,'Color',[SpecLineCols(c1(p),:,1)])
            end
            if guiHandlesSpec.specPresets.Value > 4 && guiHandlesSpec.specPresets.Value <= 6
                set(h,'Color',[SpecLineCols(c1(p),:,2)])
            end
            if guiHandlesSpec.specPresets.Value > 6
                set(h,'Color',[SpecLineCols(c1(p),:,3)])
            end

            if guiHandlesSpec.Sub100HzCheck{c1(p)}.Value==1
                set(h2,'xtick',[0 20 40 60 80 100], 'yminortick','on')
                axis([0 100 baselineYlines(guiHandlesSpec.checkboxPSD.Value+1) climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))])
                h=plot([round(Flim1) round(Flim1)],[baselineYlines(guiHandlesSpec.checkboxPSD.Value+1) climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))],'k--');
                set(h,'linewidth',1)
                h=plot([round(Flim2) round(Flim2)],[baselineYlines(guiHandlesSpec.checkboxPSD.Value+1) climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))],'k--');
                set(h,'linewidth',1)
            else
                set(h2,'xtick',[0 : ((A_lograte(guiHandlesSpec.FileSelect{k}.Value) / 2) * 1000 / 5) : (A_lograte(guiHandlesSpec.FileSelect{k}.Value) / 2) * 1000],'yminortick','on')
                axis([0 (A_lograte(guiHandlesSpec.FileSelect{k}.Value) / 2) * 1000 baselineYlines(guiHandlesSpec.checkboxPSD.Value+1) climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))])
            end

            xlabel('Frequency (Hz)')
            if guiHandlesSpec.checkboxPSD.Value
                ylabel(['PSD (dB)'])
            else
                ylabel(['Amplitude'])
            end

            h=text(2,climScale(guiHandlesSpec.checkboxPSD.Value+1, c1(p))*.95,axLabel{c2(p)});
            set(h,'Color',[.2 .2 .2],'fontsize',fontsz,'fontweight','bold')

            grid on
        end
    end
end
set(PTspecfig, 'pointer', 'arrow')
updateSpec=0;

end

repo: PIDtoolbox


text:
classdef Layer < handle
  %LAYER Base class for a network layer in a DagNN

  properties (Access = public, Hidden, Transient)
    net
    layerIndex
  end

  methods
    function outputs = forward(obj, inputs, params)
    %FORWARD Forward step
    %  OUTPUTS = FORWARD(OBJ, INPUTS, PARAMS) takes the layer object OBJ
    %  and cell arrays of inputs and parameters and produces a cell
    %  array of outputs evaluating the layer forward.
      outputs = {} ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutpus)
    %BACKWARD  Bacwkard step
    %  [DERINPUTS, DERPARAMS] = BACKWARD(OBJ, INPUTS, INPUTS, PARAMS,
    %  DEROUTPUTS) takes the layer object OBJ and cell arrays of
    %  inputs, parameters, and output derivatives and produces cell
    %  arrays of input and parameter derivatives evaluating the layer
    %  backward.
      derInputs = {} ;
      derOutputs = {} ;
    end

    function reset(obj)
    %RESET Restore internal state
    %  RESET(OBJ)resets the layer OBJ, clearing any internal state.
    end

    function params = initParams(obj)
    %INIT Initialize layer parameters
    %  PARAMS = INIT(OBJ) takes the layer OBJ and returns a cell
    %  array of layer parameters PARAMS with some initial
    %  (e.g. random) values.
      params = {} ;
    end

    function move(obj, device)
    %MOVE Move data to CPU or GPU
    %  MOVE(DESTINATION) moves the data associated to the layer object OBJ
    %  to either the 'gpu' or the 'cpu'. Note that variables and
    %  parameters are moved automatically by the DagNN object
    %  containing the layer, so this operation affects only data
    %  internal to the layer (e.g. the mask in dropout).
    end

    function forwardAdvanced(obj, layer)
    %FORWARDADVANCED  Advanced driver for forward computation
    %  FORWARDADVANCED(OBJ, LAYER) is the advanced interface to compute
    %  the forward step of the layer.
    %
    %  The advanced interface can be changed in order to extend DagNN
    %  non-trivially, or to optimise certain blocks.

      in = layer.inputIndexes ;
      out = layer.outputIndexes ;
      par = layer.paramIndexes ;
      net = obj.net ;

      inputs = {net.vars(in).value} ;

      % give up if any of the inputs is empty (this allows to run
      % subnetworks by specifying only some of the variables as input --
      % however it is somewhat dangerous as inputs could be legitimaly
      % empty)
      if any(cellfun(@isempty, inputs)),         return ; end

      % clear inputs if not needed anymore
      for v = in
        net.numPendingVarRefs(v) = net.numPendingVarRefs(v) - 1 ;
        if net.numPendingVarRefs(v) == 0
          if ~net.vars(v).precious & ~net.computingDerivative & net.conserveMemory
            net.vars(v).value = [] ;
          end
        end
      end

      % call the simplified interface
      outputs = obj.forward(inputs, {net.params(par).value}) ;
      [net.vars(out).value] = deal(outputs{:}) ;

    end

    function backwardAdvanced(obj, layer)
    %BACKWARDADVANCED Advanced driver for backward computation
    %  BACKWARDADVANCED(OBJ, LAYER) is the advanced interface to compute
    %  the backward step of the layer.
    %
    %  The advanced interface can be changed in order to extend DagNN
    %  non-trivially, or to optimise certain blocks.
      in = layer.inputIndexes ;
      out = layer.outputIndexes ;
      par = layer.paramIndexes ;
      net = obj.net ;

      inputs = {net.vars(in).value} ;
      derOutputs = {net.vars(out).der} ;
%       for i = 1:numel(derOutputs)
%         if isempty(derOutputs{i}), return ; end
%       end
        if all(cellfun(@isempty,derOutputs)), return ; end
        if any(cellfun(@isempty, inputs)),
            return ;
        end

      if net.conserveMemory
        % clear output variables (value and derivative)
        % unless precious
        for i = out
          if net.vars(i).precious, continue ; end
          net.vars(i).der = [] ;
          net.vars(i).value = [] ;
        end
      end

      % compute derivatives of inputs and paramerters
      [derInputs, derParams] = obj.backward ...
        (inputs, {net.params(par).value}, derOutputs) ;

      % accumuate derivatives
      for i = 1:numel(in)
        v = in(i) ;
        if net.numPendingVarRefs(v) == 0 || isempty(net.vars(v).der)
          net.vars(v).der = derInputs{i} ;
        elseif ~isempty(derInputs{i})
          net.vars(v).der = net.vars(v).der + derInputs{i} ;
        end
        net.numPendingVarRefs(v) = net.numPendingVarRefs(v) + 1 ;
      end
%       if isempty(derParams),
%           return;
%       end;
      for i = 1:numel(par)
        p = par(i) ;
        if (net.numPendingParamRefs(p) == 0 && ~net.accumulateParamDers) ...
              || isempty(net.params(p).der)
          net.params(p).der = derParams{i} ;
        else
          net.params(p).der = net.params(p).der + derParams{i} ;
        end
        net.numPendingParamRefs(p) = net.numPendingParamRefs(p) + 1 ;
      end

    end

    function rfs = getReceptiveFields(obj)
    %GETRECEPTIVEFIELDS  Get receptive fields.
    %   RFS = GETRECEPTIVEFIELDS(OBJ) gets the receptive fields
    %   of each output varaibles in each input variable.
    %
    %   A *receptive field* is a structure with fields
    %
    %   - size: size of the receptive field
    %   - stride: stride of the receptive field
    %   - offset: offset of the receptive field
    %
    %   It should be interpreted as follows. Given a pixel of
    %   vertical coordinate u in an output variable OUT(y,...) , the first and last
    %   pixels affecting that pixel in an input variable IN(v,...) are:
    %
    %        v_first = stride(1) * (y - 1) + offset(1) - size(1)/2 + 1
    %        v_last  = stride(1) * (y - 1) + offset(1) + size(1)/2 + 1
    %
    %   RFS is a struct array of such structure, with one row for each
    %   input variable and one column for each output variable, expressing
    %   all possible combinations of inputs and outputs.
      rfs = [] ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes = {} ;
    end

    function load(obj, varargin)
    %LOAD Initialize the layer from a paramter structure
    %  LOAD(OBJ, S) initializes the layer object OBJ from the parameter
    %  structure S.  It is the opposite of S = SAVE(OBJ).
    %
    %  LOAD(OBJ, OPT1, VAL1, OPT2, VAL2, ...) uses instead the
    %  option-value pairs to initialize the object properties.
    %
    %  LOAD(OBJ, {OPT1, VAL1, OPT2, VAL2, ...}) is an equivalent form
    %  to the previous call.
      s = dagnn.Layer.argsToStruct(varargin{:}) ;
      for f = fieldnames(s)'
        f = char(f) ;
        obj.(f) = s.(f) ;
      end
    end

    function s = save(obj)
    %SAVE Save the layer configuration to a parameter structure
    %  S = SAVE(OBJ) extracts all the properties of the layer object OBJ
    %  as a structure S. It is the oppostie of LOAD(OBJ, S).
    %
    %  By default, properties that are marked as transient,
    %  dependent, abstract, or private in the layer object are not
    %  saved.
      s = struct ;
      m = metaclass(obj) ;
      for p = m.PropertyList'
        if p.Transient || p.Dependent || p.Abstract, continue ; end
        s.(p.Name) = obj.(p.Name) ;
      end
    end

    function attach(obj, net, index)
    %ATTACH  Attach the layer to a DAG.
    %   ATTACH(OBJ, NET, INDEX). Override this function to
    %   configure parameters or variables.
      obj.net = net ;
      obj.layerIndex = index ;
      for input = net.layers(index).inputs
        net.addVar(char(input)) ;
      end
      for output = net.layers(index).outputs
        net.addVar(char(output)) ;
      end
      for param = net.layers(index).params
        net.addParam(char(param)) ;
        p = net.getParamIndex(char(param)) ;
      end
    end

  end % methods

  methods (Static)
    function s = argsToStruct(varargin)
    %ARGSTOSTRUCT  Convert varadic arguments to structure
    %  S = ARGSTOSTRCUT('opt1', val1, ....) converts the list of options
    %  into a structure.
    %
    %  S = ARGSTOSTRUCT({'opt1', val1, ...}) is equivalent to the
    %  previous call.
    %
    %  S = ARGSTOSTRUCT(S) where S is a structure returns S as is.

      if numel(varargin) == 1 && isstruct(varargin{1})
        s = varargin{1} ;
      else
        if numel(varargin) == 1 && iscell(varargin{1})
          args = varargin{1} ;
        else
          args = varargin ;
        end
        s = cell2struct(args(2:2:end),args(1:2:end),2) ;
      end
    end
  end

end % classdef

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%%***********************************************************
%% lmiexamp1: generate SDP data for the following LMI problem
%%
%%  max  -eta
%%  s.t. B*P + P*B'  <= 0
%%        -P         <= -I
%%         P - eta*I <= 0
%%         P(1,1)     = 1
%%***********************************************************
%% Here is an example on how to use this function to
%% find an optimal P.
%%
%% B = [-1  0  0; 5 -2  0; 1  1 -1];
%% [blk,At,C,b] = lmiexamp1(B);
%% [obj,X,y,Z] = sqlp(blk,At,C,b);
%% P = smat(blk(1,:),y);
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [blk,At,C,b] = lmiexamp1(B);

   n = length(B); n2 = n*(n+1)/2;
   I  = speye(n);
   z0 = sparse(n2,1);
   blktmp{1,1} = 's'; blktmp{1,2} = n;
%%
   blk{1,1} = 's'; blk{1,2} = n;
   blk{2,1} = 's'; blk{2,2} = n;
   blk{3,1} = 's'; blk{3,2} = n;
   blk{4,1} = 'u'; blk{4,2} = 1;
%%
   At{1,1} = [lmifun(B,I),     z0];
   At{2,1} = [lmifun(-I/2,I),  z0];
   At{3,1} = [lmifun(I/2,I),   svec(blktmp,-I,1)];
   At{4,1} = sparse([1, zeros(1,n2)]);
%%
   C{1,1} = sparse(n,n);
   C{2,1} = -speye(n);
   C{3,1} = sparse(n,n);
   C{4,1} = 1;
%%
   b = [zeros(n2,1); -1];
%%**********************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%
% Calculate u-v for general vectors u and v in a product space
%
% Syntax:  >> w = VectorMinus(u,v)
%
%  Input: u, v --- (matrix/string/cell) vectors in a product space given as
%                  (i) mxn matrices of the sime dimension
%                  (ii) polynomials in character strings
%                  (iii) cell array of (i) or (ii)
%
% Output: w --- (matrix/string/cell) the result of u-v stored in the same
%                 way as u and v
%
% Example:
% >> u = {[1 2 3; 4 5 6], '7+8*x*y+9*z^3'};
% >> v = {[1 1 1; 1 1 1],'1+x+x*y+z^3'};
% >> w  = VectorMinus(u,v);
% >> w{:}
% ans =
%
%     0     1     2
%     3     4     5
% ans =
% 6 - x + 7*x*y + 8*z^3
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function sd = loadstoragedata(sd_table, mpc)
%LOADSTORAGEDATA   Load StorageData table into StorageData struct.
%
%   SD = LOADSTORAGEDATA(SD_TABLE)
%   SD = LOADSTORAGEDATA(SD_TABLE, GEN)
%   SD = LOADSTORAGEDATA(SD_TABLE, MPC)
%
%   Loads data from a StorageData table struct, or from a M-file or MAT-file
%   that returns such a struct, and converts it to a StorageData struct.
%   StorageData contains all parameters required by MOST for storage
%   units that is not included in MPC or xGenData.
%
%   If the first argument is the name of an M-file the optional second
%   argument will be passed to it. This can be useful for cases where
%   the user wishes to use data from GEN or GENCOST to set values in the
%   StorageData.
%
%   Inputs:
%       SD_TABLE : a StorageData table struct or the name of an M-file
%                   or MAT-file that returns one, with the following fields
%           .colnames : N dimensional cell array of names corresponding
%               to the columns of the 'data' field. Valid column names
%               are the same as the output field names. All columns are
%               optional and the correpsonding defaults are listed
%               next the output field names below.
%           .data : (NS x N) matrix of data
%           .MinStorageLevel                 (optional)     (scalar)
%           .MaxStorageLevel                 (optional)     (scalar)
%           .OutEff                          (optional)     (scalar)
%           .InEff                           (optional)     (scalar)
%           .LossFactor                      (optional)     (scalar)
%           .rho                             (optional)     (scalar)
%           (values in any of the scalar fields listed above are
%            overridden by any corresponding values in the 'data' table)
%       GEN : (optional) standard GEN matrix for generators corresponding
%             to storage units in SD_TABLE
%       MPC : (optional) MATPOWER case struct containing GEN and GENCOST
%             matrices for generators corresponding to storage units in
%             SD_TABLE
%
%   Output:
%       SD : a StorageData struct, contains the following fields,
%               all of which are (NS x 1) vectors unless otherwise
%               indicated. Defaults are as provided by MD_INIT.
%           .UnitIdx                            (ns x 1)
%           .ExpectedTerminalStorageAim         (ns x 1)
%           .ExpectedTerminalStorageMin         (ns x 1)
%           .ExpectedTerminalStorageMax         (ns x 1)
%           .InitialStorage                     (ns x 1)
%           .InitialStorageLowerBound           (ns x 1)
%           .InitialStorageUpperBound           (ns x 1)
%           .InitialStorageCost                 (ns x 1)
%           .TerminalStoragePrice               (ns x 1)
%           .TerminalChargingPrice0             (ns x 1)
%           .TerminalDischargingPrice0          (ns x 1)
%           .TerminalChargingPriceK             (ns x 1)
%           .TerminalDischargingPriceK          (ns x 1)
%           .MinStorageLevel                    (ns x 1) or (1 x 1)
%           .MaxStorageLevel                    (ns x 1) or (1 x 1)
%           .OutEff                             (ns x 1) or (1 x 1)
%           .InEff                              (ns x 1) or (1 x 1)
%           .LossFactor                         (ns x 1) or (1 x 1)
%           .rho                                (ns x 1) or (1 x 1)

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

%% define named indices into data matrices
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

%% input arg handling
if ischar(sd_table)
    infile = sprintf(' in file: ''%s''', sd_table);
else
    infile = '';
end
if nargin < 2
    args = [];
    gen = [];
else
    args = mpc;
    if isstruct(mpc)
        gen = mpc.gen;
    else
        gen = mpc;
    end
end
ns0 = size(gen, 1);      %% number of storage units in gen matrix

%% load the table
fields = {'colnames', 'data'};
sdt = loadgenericdata(sd_table, 'struct', fields, 'sd_table', args);
if isempty(sdt)
    sdt.colnames = {};
    sdt.data = [];
end
[ns, nc] = size(sdt.data);

%% check consistency of dimensions of data table, colnames and gen
if ns0 > 0 && ns > 0 && ns ~= ns0
    error('loadstoragedata: # of rows in ''data'' table (%d) do not equal rows in GEN (%d)%s', ...
        ns, ns0, infile);
end
if nc ~= length(sdt.colnames)
    error('loadstoragedata: # of columns in ''data'' table (%d) do match entries in ''colnames'' (%d)%s', ...
        nc, length(sdt.colnames), infile);
end

%% initialize storage data with default values from md_init's Storage field
sd = getfield(md_init(), 'Storage');

%% set individual fields
f = {   'MinStorageLevel', ...
        'MaxStorageLevel', ...
        'OutEff', ...
        'InEff', ...
        'LossFactor', ...
        'rho'
    };
for k = 1:length(f)
    if isfield(sd_table, f{k});
        sd.(f{k}) = sd_table.(f{k});
    end
end

%% create cell array with valid fields for checking
valid_fields = {
        'UnitIdx', ...
        'ExpectedTerminalStorageAim', ...
        'ExpectedTerminalStorageMin', ...
        'ExpectedTerminalStorageMax', ...
        'InitialStorage', ...
        'InitialStorageLowerBound', ...
        'InitialStorageUpperBound', ...
        'InitialStorageCost', ...
        'TerminalStoragePrice', ...
        'TerminalChargingPrice0', ...
        'TerminalDischargingPrice0', ...
        'TerminalChargingPriceK', ...
        'TerminalDischargingPriceK', ...
        'MinStorageLevel', ...
        'MaxStorageLevel', ...
        'OutEff', ...
        'InEff', ...
        'LossFactor', ...
        'rho'
    };

%% copy data from table
for j = 1:nc
    %% check for valid colname
    if ~any(strcmp(sdt.colnames{j}, valid_fields))
        error('loadstoragedata: ''%s'' is not a valid StorageData field name%s', sdt.colnames{j}, infile);
    end
    %% copy data
    sd.(sdt.colnames{j}) = sdt.data(:, j);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [handle] = ml_plot_eigenvalues( L )
%PLOT_EIGENVALUES Simple plotting function to visualize eigenvalues
%   The student should convert the Eigenvalue matrix to a vector and
%   visualize the values as a 2D plot.
%   input -----------------------------------------------------------------
%
%       o L      : (N x N), Diagonal Matrix composed of lambda_i
%

lambda = diag(L)';

handle = figure;
plot(lambda, '--b', 'LineWidth', 2)
set(gca,'XTick',[1:1:length(lambda)])
title('Eigenvalues')
ylabel('Eigenvalues')
xlabel('Eigenvector index')
grid on

end

author: epfl-lasa


text:
function SMTout=time_delay_embedding(SMT,CNT,tau,time_interval)
% TDE: time-delay embedding

SMT_in=SMT;
CNT_in=CNT;
interval=time_interval;

SMT_in.x=permute(SMT_in.x,[1,3,2]);
SMT_in_d=prep_segmentation(CNT_in, {'interval', interval-tau});
SMT_in_d.x=permute(SMT_in_d.x,[1,3,2]);
SMT_in.x=[SMT_in.x,SMT_in_d.x];
SMT_in.x=permute(SMT_in.x,[1,3,2]);

SMTout=SMT_in;
end

author: PatternRecognition

repo: OpenBMI


text:
function test_suite = test_distancePointPolyline
%TESTDISTANCEPOINTEDGE  Test suite for function distancePointPolyline
%   output = test_distancePointPolyline(input)
%
%   Example
%   testDistancePointPolyline
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testPointOnPolyline(testCase) %#ok<*DEFNU>
% test if the distance betwwen polyline and a point on polyline is 0

% create polyline
poly = [0 0;10 0;10 10;0 10];

p1 = [0 0];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(0, d, 'AbsTol', .01);

p1 = [10 0];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(0, d, 'AbsTol', .01);

p1 = [5 10];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(0, d, 'AbsTol', .01);

p1 = [10 5];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(0, d, 'AbsTol', .01);

function testPointNotOnPolyline(testCase)

% test various points

% create polyline
poly = [0 0;10 0;10 10;0 10];

p1 = [0 -10];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(10, d, 'AbsTol', .01);

p1 = [20 0];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(10, d, 'AbsTol', .01);

p1 = [0 5];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(5, d, 'AbsTol', .01);

p1 = [5 5];
d = distancePointPolyline(p1, poly);
testCase.assertEqual(5, d, 'AbsTol', .01);

function testManyPoints(testCase)

% test several points in one call...

% create polyline
poly = [0 0;10 0;10 10;0 10];

pts = [0 -10;20 0;0 5; 5 5];
dist = distancePointPolyline(pts, poly);
exp = [10 10 5 5]';

testCase.assertEqual(exp, dist, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [yData,signalData,signalIdx,signalAmps] = ...
                file_embeddingSubSampling(projectionFile,parameters)
%file_embeddingSubSampling finds the potential training set contributions
%from a single file (called by runEmbeddingSubSampling.m)
%
%   Input variables:
%
%       projectionFile -> file to be analyzed (should contain a variable
%                           called 'projections')
%       parameters -> struct containing non-default choices for parameters
%
%
%   Output variables:
%
%       yData -> Nx2 array containing t-SNE embedding
%       signalData -> wavelet data corresponding to yData
%       signalIdx -> idx used
%       signalAmps -> wavelet amplitudes of signalData
%
% (C) Gordon J. Berman, 2014
%     Princeton University

    rtol = parameters.training_relTol;
    perplexity = parameters.training_perplexity;
    numPoints = parameters.training_numPoints;

    fprintf(1,'\t Loading Projections\n');
    load(projectionFile,'projections');

    N = length(projections(:,1));
    numModes = parameters.pcaModes;
    skipLength = floor(N / numPoints);
    if skipLength == 0
        skipLength = 1;
        numPoints = N;
    end
    firstFrame = mod(N,numPoints) + 1;
    signalIdx = firstFrame:skipLength:(firstFrame + (numPoints-1)*skipLength);

    fprintf(1,'\t Calculating Wavelets\n');
    [data,~] = findWavelets(projections,numModes,parameters);
    amps = sum(data,2);

    signalData = bsxfun(@rdivide,data(signalIdx,:),amps(signalIdx));
    signalAmps = amps(signalIdx);

    clear data amps;

    fprintf(1,'\t Calculating Distances\n');
    [D,~] = findKLDivergences(signalData);

    fprintf(1,'\t Running t-SNE\n');
    parameters.relTol = rtol;
    parameters.perplexity = perplexity;
    [yData,~,~,~] = tsne_d(D,parameters);

author: gordonberman

repo: MotionMapper


text:
% Mesh spectrum

function [L,H,d] = mesh_spectrum(S,n,mode)

%[L,H,d] = ct_mesh_spectrum(S,n,mode)
% Compute the mesh laplace matrix and its spectrum
% input,
% S: mesh file, it has to have a pnt and a tri field
% n: number of mesh harmonic functions
% mode: 'full' for the full graph, 'half' if you want to do the first and
% the second half independently (this is useful if your graph is composed
% by two connected components) -> JM note: I don't understand this
% output,
% L: mesh laplacian matrix
% H: matrix containing a mesh harmonic functions per column
% d: spectrum of the negative Laplacian matrix, its units are 1/space^2
% (spatial frequencies are obtained as sqrt(d))

if isempty(mode)
  mode=1;
end

if mode==1
    pnt{1} = S.pos;
    tri{1} = S.tri;
elseif mode==2
    pnt{1} = S.pos(1:end/2,:);
    tri{1} = S.tri(1:end/2,:);
    pnt{2} = S.pos(end/2+1:end,:);
    tri{2} = S.tri(end/2+1:end,:) - size(pnt{1},1);
end

for j = 1:length(pnt)
    if length(pnt)==2&&j == 1
        disp('Computing the spectrum of the the first hemisphere')
    elseif length(pnt)==2&&j == 2
        disp('Computing the spectrum of the the second hemisphere')
    end
    L{j} = mesh_laplacian(pnt{j},tri{j});
    L{j} = (L{j} + L{j}')/2;
    disp('Computing the spectrum of the negative Laplacian matrix')
    [H{j},D] = eigs(L{j},n,'sm');
    d{j} = diag(D);
    disp('Diagonalization completed')
end


text:
function v = setLayerOutputs(obj, layer, outputs)
%SETLAYEROUTPUTS  Set or change the outputs of a layer
%   Example: NET.SETLAYEROUTPUTS('layerName', {'output1', 'output2', ...})

v = [] ;
l = obj.getLayerIndex(layer) ;
for output = outputs
  v(end+1) = obj.addVar(char(output)) ;
end
obj.layers(l).outputs = outputs ;
obj.rebuild() ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
% Copyright (C) 2013 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Quan Wang and Kim L. Boyer.
% Feature Learning by Multidimensional Scaling and its Applications in Object Recognition.
% 2013 26th SIBGRAPI Conference on Graphics, Patterns and Images (Sibgrapi). IEEE, 2013.
%
% For commercial use, please contact the authors.

% The target function to be minimized at each step.

function f=MDS_cost_vector(x,V_dist,X)

X=bsxfun(@minus,X,x);
V=sum(X.*X,2);
f=sqrt(V)-V_dist;

author: Sable

repo: mcbench-benchmarks


text:
function sample = SampleGaussian(mu, sigma)

% sample from the Gaussian distribution specifed by mean value mu and standard deviation sigma
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

 sample = mu + sigma*randn(1,1);

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function vr = rotateVector(v, angle)
%ROTATEVECTOR Rotate a vector by a given angle.
%
%   VR = rotateVector(V, THETA)
%   Rotate the vector V by an angle THETA, given in radians.
%
%   Example
%   rotateVector([1 0], pi/2)
%   ans =
%       0   1
%
%   See also
%   vectors2d, transformVector, createRotation

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

% precomputes angles
cot = cos(angle);
sit = sin(angle);

% compute rotated coordinates
vr = [cot * v(:,1) - sit * v(:,2) , sit * v(:,1) + cot * v(:,2)];


text:
function test_old_ft_multiplotER

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

% this script tests the functionality of ft_multiplotER with respect to the
% different input datatypes. no other functionality is tested.
% the script has been written in order to test a clean up of the code

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/eeg/preproc_neuroscan16');
load(filename);

%there's an unresolved issue with duplicate labels 'FREE'
%FIXME
data.label{1} = 'FREE1';
data.label{2} = 'FREE2';
data.label{3} = 'FREE3';
data.label{4} = 'FREE4';

cfg = [];
cfg.channel = data.label(5:end);
cfg.preproc.demean = 'yes';
cfg.trials = 1:5;
tlck1 = ft_timelockanalysis(cfg, data);
cfg.trials = 6:10;
tlck2 = ft_timelockanalysis(cfg, data);

cfg = [];
cfg.layout = 'biosemi64.lay';
lay = ft_prepare_layout(cfg);

cfg        = [];
cfg.layout = lay;

%plot timelocked-data
ft_multiplotER(cfg, tlck1, tlck2);

%create frequency-data
cfgf = [];
cfgf.method = 'mtmfft';
cfgf.foilim = [0 100];
cfgf.taper  = 'hanning';
cfgf.output = 'pow';
cfgf.trials = 1:5;
cfgf.channel = data.label(5:end); % the first 4 are empty
freq1 = ft_freqanalysis(cfgf, data);
cfgf.trials = 6:10;
freq2 = ft_freqanalysis(cfgf, data);

%plot frequency-data
ft_multiplotER(cfg, freq1, freq2);

%create connectivity-data
cfgf.output = 'fourier';
cfgf.trials = 1:10;
freqx = ft_freqanalysis(cfgf, data);

cfgc2 = [];
cfgc2.method = 'coh';
coh   = ft_connectivityanalysis(cfgc2, freqx);

%plot connectivity-data
cfg.refchannel = 'gui';
cfg.parameter = 'cohspctrm';
ft_multiplotER(cfg, coh);

% FIXME this causes a crash when a new reference is selected and the old one is not unselected
% FIXME it also crashes when more than one cohref is selected

%create connectivity-data with sparse linear indexing
cfgc2.channelcmb = [repmat(freqx.label(5),[numel(freqx.label)-1 1]) freqx.label([1:4 6:end])];
coh2  = ft_connectivityanalysis(cfgc2, freqx);

%plot
cfg.refchannel = 'Cz';
coh2 = ft_checkdata(coh2, 'cmbstyle', 'full');
ft_multiplotER(cfg, coh2);
coh2 = ft_checkdata(coh2, 'cmbstyle', 'sparse');
ft_multiplotER(cfg, coh2);

%create connectivity-data with very sparse linear indexing
cfgc2.channelcmb = cfgc2.channelcmb(1:25,:);
coh3   = ft_connectivityanalysis(cfgc2, freqx);

%plot
ft_multiplotER(cfg, coh3);

%create connectivity-data with even sparser linear indexing
% cfgc2.channelcmb = [repmat(freq2.label(5),[10 1]) freq2.label(21:30)';repmat(freq2.label(10),[10 1]) freq2.label(21:30)'];
% coh4 = ft_connectivityanalysis(cfgc2, freq2);
%
% %plot: this breaks
% cfg.cohrefchannel = 'gui';
% ft_topoplotER(cfg, coh4);
%
% %plot: this works
% cfg.cohrefchannel = coh4.labelcmb(1,1);
% ft_topoplotER(cfg, coh4);
%
% %create connectivity-data with asymmetry
% %the data are probably not full-rank creating a problem for the sf
% %freq3 = ft_selectdata(freq2, 'channel', freq2.label(1:30));
% %%subselection of channels does not help
% freq4 = freq2transfer([], freq2);
%
% cfgc2 = [];
% cfgc2.method = 'granger';
% granger = ft_connectivityanalysis(cfgc2, freq4);
%
% cfg.cohrefchannel = 'gui';
% cfg.zparam = 'grangerspctrm';
% ft_topoplotER(cfg, granger);

author: fieldtrip

repo: fieldtrip


text:
function error_msg = mdp_check_square_stochastic( Z )

% mdp_check_square_stochastic     Check if Z is a square stochastic matrix
% Arguments ---------------------------------------------------------------
%   Z(SxS) = a matrix
% Evaluation --------------------------------------------------------------
%   error_msg = error message or '' if correct

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

error_msg = '';

[s1 s2] = size(Z);
if s1 ~=s2
   error_msg = 'MDP Toolbox ERROR: Matrix must be square';
elseif max(abs(sum(Z,2)-ones(s2,1))) > 10^(-12)
   error_msg = 'MDP Toolbox ERROR: Row sums of the matrix must be 1';
elseif ~isempty(find (Z < 0))
   error_msg = 'MDP Toolbox ERROR: Probabilities must be non-negative';
end;


text:
% compute matrix multiply X * Y
%
function [grad] = B_matrix_multiply(input_layers, curr_layer, future_layers)
X = input_layers{1}.a;
Y = input_layers{2}.a;

% input1 and input2

[Dx,Tx,Nx] = size(X);
[Dy,Ty,Ny] = size(Y);
if Tx ~=Dy
    fprintf('Error: matrix size not match each other\n');
end

future_grad = GetFutureGrad(future_layers, curr_layer);

if Nx==1 && Ny==1
    gradX = conj(future_grad * Y');
    gradY = X' * future_grad;
else
    [validFrameMask, variableLength] = getValidFrameMask(curr_layer);
    if variableLength
        future_grad = PadShortTrajectory(future_grad, validFrameMask, 0);
    end

    if IsInGPU(X(1))
        gradX = gpuArray.zeros(size(X));
        gradY = gpuArray.zeros(size(Y));
    else
        gradX = zeros(size(X));
        gradY = zeros(size(Y));
    end

    if Nx==1 && Ny>1    %use the same X matrix for all Y matrices
        Y = PadShortTrajectory(Y, validFrameMask, 0);
        for i = 1:Ny
            gradX = gradX + conj(future_grad(:,:,i) * Y(:,:,i)');
            gradY(:,:,i) = X' * future_grad(:,:,i);
        end
    elseif Nx>1 && Ny==1    %use the same Y matrix for all X matrices
        X = PadShortTrajectory(X, validFrameMask, 0);
        for i = 1:Nx
            gradX(:,:,i) = conj(future_grad(:,:,i) * Y');
            gradY = gradY + X(:,:,i)' * future_grad(:,:,i);
        end
    elseif Nx>1 && Ny>1
        Y = PadShortTrajectory(Y, validFrameMask, 0);
        X = PadShortTrajectory(X, validFrameMask, 0);

        for i = 1:N
            gradX(:,:,i) = conj(future_grad(:,:,i) * Y(:,:,i)');
            gradY(:,:,i) = X(:,:,i)' * future_grad(:,:,i);
        end
    end
end
grad{1} = gradX;
grad{2} = gradY;
end

author: singaxiong

repo: SignalGraph


text:
classdef GridFlex < uix.GridFlex
    %uiextras.GridFlex  Container with contents arranged in a resizable grid
    %
    %   obj = uiextras.GridFlex() creates a new new grid layout with
    %   draggable dividers between elements. The number of rows and columns
    %   to use is determined from the number of elements in the RowSizes
    %   and ColumnSizes properties respectively. Child elements are
    %   arranged down column one first, then column two etc. If there are
    %   insufficient columns then a new one is added. The output is a new
    %   layout object that can be used as the parent for other
    %   user-interface components. The output is a new layout object that
    %   can be used as the parent for other user-interface components.
    %
    %   obj = uiextras.GridFlex(param,value,...) also sets one or more
    %   parameter values.
    %
    %   See the documentation for more detail and the list of properties.
    %
    %   Examples:
    %   >> f = figure();
    %   >> g = uiextras.GridFlex( 'Parent', f, 'Spacing', 5 );
    %   >> uicontrol( 'Parent', g, 'Background', 'r' )
    %   >> uicontrol( 'Parent', g, 'Background', 'b' )
    %   >> uicontrol( 'Parent', g, 'Background', 'g' )
    %   >> uiextras.Empty( 'Parent', g )
    %   >> uicontrol( 'Parent', g, 'Background', 'c' )
    %   >> uicontrol( 'Parent', g, 'Background', 'y' )
    %   >> set( g, 'ColumnSizes', [-1 100 -2], 'RowSizes', [-1 -2] );
    %
    %   See also: uiextras.Grid
    %             uiextras.HBoxFlex
    %             uiextras.VBoxFlex
    %             uiextras.Empty

    %  Copyright 2009-2014 The MathWorks, Inc.
    %  $Revision: 1062 $ $Date: 2014-10-30 13:30:17 +0000 (Thu, 30 Oct 2014) $

    properties( Hidden, Access = public, Dependent )
        Enable % deprecated
        RowSizes % heights of contents, in pixels and/or weights
        MinimumRowSizes % minimum heights of contents, in pixels
        ColumnSizes % widths of contents, in pixels and/or weights
        MinimumColumnSizes % minimum widths of contents, in pixels
        ShowMarkings
    end

    methods

        function obj = GridFlex( varargin )

            % Call uix constructor
            obj@uix.GridFlex( varargin{:} )

            % Auto-parent
            if ~ismember( 'Parent', varargin(1:2:end) )
                obj.Parent = gcf();
            end

        end % constructor

    end % structors

    methods

        function value = get.Enable( ~ )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

            % Return
            value = 'on';

        end % get.Enable

        function set.Enable( ~, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'on','off'} ) ), ...
                'uiextras:InvalidPropertyValue', ...
                'Property ''Enable'' must be ''on'' or ''off''.' )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

        end % set.Enable

        function value = get.RowSizes( obj )

            % Get
            value = obj.Heights;

        end % get.RowSizes

        function set.RowSizes( obj, value )

            % Set
            obj.Heights = value;

        end % set.RowSizes

        function value = get.MinimumRowSizes( obj )

            % Get
            value = obj.MinimumHeights;

        end % get.MinimumRowSizes

        function set.MinimumRowSizes( obj, value )

            % Set
            obj.MinimumHeights = value;

        end % set.MinimumRowSizes

        function value = get.ColumnSizes( obj )

            % Get
            value = obj.Widths;

        end % get.ColumnSizes

        function set.ColumnSizes( obj, value )

            % Get
            obj.Widths = value;

        end % set.ColumnSizes

        function value = get.MinimumColumnSizes( obj )

            % Get
            value = obj.MinimumWidths;

        end % get.MinimumColumnSizes

        function set.MinimumColumnSizes( obj, value )

            % Get
            obj.MinimumWidths = value;

        end % set.MinimumColumnSizes

        function value = get.ShowMarkings( obj )

            % Get
            value = obj.DividerMarkings;

        end % get.ShowMarkings

        function set.ShowMarkings( obj, value )

            % Set
            obj.DividerMarkings = value;

        end % set.ShowMarkings

    end % accessors

end % classdef

author: goGPS-Project


text:
function a = erfc(a)
%ERFC         Gradient complementary error function erfc(a)
%

% written  05/31/13     S.M. Rump
%

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

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');
  % factorLB <= 2/sqrt(pi) <= factorUB
  INTLAB_STDFCTS_ERF = getappdata(0,'INTLAB_STDFCTS_ERF');
  factorLB = INTLAB_STDFCTS_ERF.TWO_SQRTPIINF;  % round to nearest
  factorUB = INTLAB_STDFCTS_ERF.TWO_SQRTPISUP;  % ~ 1.12

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  a.x = full(a.x);
  ax = exp(-a.x(:).^2);
  a.x = erfc(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      a.dx = intval(-factorUB,-factorLB,'infsup') * ...
          intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
    else
      a.dx = (-factorLB) * sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    if isa(a.x,'intval')
      ax = intval(-factorUB,-factorLB,'infsup') * ax;
      a.dx = a.dx .* ax(:,ones(1,N));
    else
      ax = (-factorLB) * ax;
      a.dx = a.dx .* ax(:,ones(1,N));
    end
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
% 3D molecular dynamics (GPU vectorized)
% by Brian Haidet for AlphaPhoenix
% published 12/15/2018
% CC non-commercial, attribution

% Don't copy this straight into a homework or something rediculous - I'm
% posting this so people can hopefully learn from it! It's not the cleanest
% code, but I added a few more comments in the areas that need
% explaining. and deleted all the extra bits of commented-out code that had
% built up. Feel free to ask more questions on youtube here!
% https://youtu.be/6DlRsPo-dxY

function [] = plotter(iteration)
    load(['c:\MatlabOutput\CsClSimFrames\data\data' num2str(iteration,'%05d') '.mat'])
    %load(['H:\MatlabOutput\CsClSimFrames\data32krun\' num2str(iteration,'%05d') '.mat'])

%PLOTTER Summary of this function goes here
%   Detailed explanation goes here
spheresize=20;
[X,Y,Z]=sphere(spheresize-1);
cdb=zeros(spheresize,spheresize,3);
cdb(:,:,3)=ones(spheresize,spheresize)+((Z-1)/2);
cdr=zeros(spheresize,spheresize,3);
cdr(:,:,1)=ones(spheresize,spheresize)+((Z-1)/2);
X=.7*X;
Y=.7*Y;
Z=.7*Z;

gasFraction=.05;
gasFraction2=.0025;
coldThreshold=.06;
annealNum=1024; %above 2048, annealing instead of instachilling
annealMedianSpeed=.022; %cool until this is the median speed of the simulation
speedLimit=.4;

set(gcf,'color','w');
subplot(3,3,[2,3,5,6,8,9]);
hold on;
xlim([-1-1,l+1+1]);
ylim([-1-1,l+1+1]);
zlim([-1-1,l+1+1]);
axis vis3d;
axis off;
zoom;
setAxes3DPanAndZoomStyle(zoom(gcf),gca,'camera');
view(45+iteration*.5,15);
numParticles=size(r,2);

% ---------------------------- calculate which particles to draw -------------------
% thresholdDist=1.25; %between 2nd and 3rd nearest neighbors
% thresholdCount=14; %number of nearest and 2nd nearest neighbors
% rG=gpuArray(r);
% totForces=gpuArray(zeros(1,numParticles));
% chunkSize=ceil(8192^2/numParticles);
% for i=1:ceil(numParticles/chunkSize) %for each chunk
%     cInds=(i-1)*chunkSize+1 : min((i)*chunkSize,numParticles);
%
%     numPcurr=length(cInds);
%     IP=rG(:,cInds)'*rG(:,:);
%     forces=gpuArray([]);
%     distances=(sqrt(abs(bsxfun(@plus,sum(rG(:,cInds)'.^2,2),sum(rG(:,:).^2,1))-2*IP))<thresholdDist);
%     totForces(cInds)=sum(distances,2);
% end
% drawParticles=totForces<thresholdCount;
% drawPartsIndex=find(drawParticles);
% ----------------------------------------------------------------------------------

% for p=drawPartsIndex
for p=1:numParticles
        delta=mean(r(:,1),2)-[l/2;l/2;l/2];
        %saveData(:,:,f) = saveData(:,:,f) -
        if charge(p)==1
            h=surface(X+r(1,p)-delta(1),Y+r(2,p)-delta(2),Z+r(3,p)-delta(3),cdr);
        else
            h=surface(X+r(1,p)-delta(1),Y+r(2,p)-delta(2),Z+r(3,p)-delta(3),cdb);
        end
        set(h,'edgecolor','none');
end
hold off;
% set(gcf, 'Position', [0 0 1920, 1080]); %<- Set size
% %set(gcf, 'renderer', 'opengl');

subplot(3,3,[4])
[counts,centers] = hist(speedsS*10*30,200);
%semilogy(centers, counts, 'ro', 'MarkerFace', 'r');
semilogy(centers, smooth(counts), 'b-');
xlabel('Distribution of Particle Speeds (units/second)');
xlim([0,max(.15*10*30,speedsS(end)*10*30)]);

title(['CsCl Crystal Structure (equal sized ions)' char(10) ...
    char(10)...
    'Computation Time: ' datestr(compTime/86400, 'DD:HH:MM:SS.FFF') char(10) ...%'GPUTemp: ' num2str(GpuT) ' Celcius' char(10) ...
    char(10) ...
    'Simulation Time: ' datestr(iteration/30/86400, 'HH:MM:SS.FFF') char(10) ...
    'Current Step: ' num2str(maxNum) char(10) ...
    'Number of Particles: ' num2str(numParticles) char(10) ...
    char(10) ...
    'Nucleating until 64 particles' char(10) ...
    'Annealed growth past ' num2str(annealNum) ' particles' char(10) ...
    'Currently cooling crystal with ' num2str(ceil(maxNum/4)) ' sink particles.' char(10) ...
    'Heated walls (2% increase in speed per bounce).' char(10) ...
    'Global speed limit of ' num2str(speedLimit*10*30) ' units/second'])

subplot(3,3,7)
plot(speedsS*10*30, 'b.');%arrayfun(@(it) norm(r(:,it)-rl(:,it)),1:numParticles));
hold on
plot([0,numParticles],[coldThreshold*10*30,coldThreshold*10*30],'r-')
plot([numParticles*(1-gasFraction)],[coldThreshold*10*30],'rx')
plot([numParticles*(1-gasFraction2)],[coldThreshold*10*30],'rx')
plot([0,numParticles],[annealMedianSpeed*10*30,annealMedianSpeed*10*30],'c-')
plot([numParticles*(.5)],[annealMedianSpeed*10*30],'cx')
hold off
ylim([0,max(.15*10*30,speedsS(end)*10*30)]);
xlim([0,numParticles]);
xlabel('All particles (sorted by speed)')
ylabel('Particle velocities')
%drawnow();
% saveas(gcf,['H:\MatlabOutput\CsClSimFrames\frames\Frame' num2str(iteration,'%05i') '.png']);
%fighand = gcf;
end

author: BrianHaidet

repo: AlphaPhoenix


text:
function vTrunc= util_trunc(v, digits, policy)
%vTrunc= util_trunc(v, <digits=4, policy>)
%
% Truncates the number v, digits gives the amount of digits after the
% comma.
%
% policy is 'floor', 'ceil', or 'round'

if ~exist('digits', 'var'), digits=4; end
if ~exist('policy','var'), policy='round'; end

a= 10^digits;
vTrunc= feval(policy, a*v)/a;


text:
function  [alpha, beta, Tau1]    =  Cal_Parameters( im, par, Dict, blk_arr, wei_arr )
[h w ch]   =   size(im);
A          =   Dict.D0;
PCA_idx    =   Dict.cls_idx;
s_idx      =   Dict.s_idx;
seg        =   Dict.seg;

b          =   par.win;
b2         =   b*b*ch;
k          =   0;
s          =   par.step;

N     =  h-b+1;
M     =  w-b+1;
L     =  N*M;
r     =  [1:s:N];
r     =  [r r(end)+1:N];
c     =  [1:s:M];
c     =  [c c(end)+1:M];
X     =  zeros(b*b,L,'single');
for i  = 1:b
    for j  = 1:b
        k    =  k+1;
        blk  =  im(i:h-b+i,j:w-b+j);
        blk  =  blk(:);
        X(k,:) =  blk';
    end
end

m_X       =   zeros(length(r)*length(c),b*b,'single');
X1        =   X';

for i = 1:par.nblk
   v            =   wei_arr(:,i);
   m_X(:,:)     =   m_X(:,:) + X1(blk_arr(:,i),:) .*v(:, ones(1,b2));
end
m_X          =   m_X';

N            =   length(r);
M            =   length(c);
L            =   N*M;
ind          =   zeros(N,M);
ind(r,c)     =   1;

X1           =   X(:, ind~=0);

alpha        =   zeros(b2, L, 'single' );
beta         =   zeros(b2, L, 'single' );
s0           =   zeros(b2, L, 'single' );

idx            =   s_idx(seg(1)+1:seg(2));
L0             =   length(idx);
alpha(:,idx)   =   A*X1(:,idx);
beta(:,idx)    =   A*m_X(:,idx);
for  k  =  1 : L0
    i           =   idx(k);
    a           =   A*( X(:, blk_arr(i, 1:par.nblk)) - repmat( m_X(:, i), 1, par.nblk ));
    s0(:,i)     =   mean(a.^2, 2);
end

for   i  = 2:length(seg)-1
    idx            =   s_idx(seg(i)+1:seg(i+1));
    cls            =   PCA_idx(idx(1));
    P              =   reshape(Dict.PCA_D(:, cls), b2, b2);
    alpha(:,idx)   =   P*X1(:,idx);
    beta(:,idx)    =   P*m_X(:,idx);
    for  j  =  1 : length(idx)
        k           =   idx(j);
        a           =   P*( X(:,blk_arr(k, 1:par.nblk)) - repmat( m_X(:, k), 1, par.nblk ));
        s0(:,k)     =   mean(a.^2, 2);
    end
end
s0       =   max(0, s0-par.nSig^2);
Tau1     =   (par.c1*sqrt(2)*par.nSig^2)./(sqrt(s0) + eps);
return;

author: lbasek

repo: image-denoising-benchmark


text:
function h = beagle(varargin)
% BEAGLE Establish connection with BeagleBoard hardware.
%
% C = beagle creates a connection object, C, for communicating with
% BeagleBoard hardware.
%
% C = beagle(hostname, username, password, builddir) allows specification
% of non-default values for the IP address, user name, password, and the
% build directory on the BeagleBoard hardware.
%
% Methods:
%
% [status, msg] = C.connect tests SSH connection to the BeagleBoard
% hardware. If connection is successful, the return value for status is
% zero. Otherwise an error has occurred during connection. The output
% argument msg provides diagnostic messages if the connection fails.
%
% [status, msg] = C.execute('command') executes the Linux command on the
% BeagleBoard hardware and returns the resulting status and standard output.
%
% C.execute('command', true) executes the Linux command on the BeagleBoard
% hardware and forces the output to the MATLAB command window.
%
% C.run(modelName) runs the previously compiled Simulink model on the
% BeagleBoard hardware.
%
% C.stop(modelName) stops the execution of the Simulink model on the
% BeagleBoard hardware.
%
% C.openShell('ssh') launches a SSH terminal session. Once the terminal
% session is started, you can execute commands on the BeagleBoard hardware
% interactively.
%
% C.openShell('serial') launches a serial terminal session. You must have a
% serial connection between your host computer and the BeagleBoard hardware.
% The parameters of the serial connection must be set to the following:
%
% Baud Rate   : 115200
% Data bits   : 8
% Stop bits   : 1
% Parity      : None
% Flow control: None
%
% Examples:
%
%    C = beagle;
%    C.execute('ls -al ~', true);
%
%  lists the contents of the home directory for the current user on the
%  BeagleBoard hardware.
%
%    [st, msg] = C.connect
%
%  tests the SSH connection to the BeagleBoard hardware.
%
%    C.run('beagleboard_gettingstarted')
%
%  runs the model 'beagleboard_gettingstarted' on the BeagleBoard hardware.
%  The model must have run previously on the BeagleBoard hardware for this
%  method to work properly.
%
%    C.openShell('ssh')
%
%  launches a SSH terminal session. You must enter the user name and the
%  password to login to the Linux shell running over SSH. Once logged in,
%  you can execute interactive shell commands.

%   Copyright 2012 The MathWorks, Inc.

narginchk(0, 4);
if nargin > 0
    h = realtime.internal.getLinuxServicesHandle('BeagleBoard', varargin{:});
else
    h = realtime.internal.getLinuxServicesHandle('BeagleBoard');
end

end


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% problem 11 - Computes the linear convolution of two sequences

% b)
x1=1:4;
x2=7:-1:1;
y=linconv(x1,x2)

% c)
y=conv(x1,x2)


text:
filename='Cantilever_triangle_fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = 1;
constraint = {'volume'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';

nsteps = 15;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
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;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function fz = trifunc (z)

% triaxial objective function

% input

%  z = function argument

% output

%  fz = function value at z

% Orbital Mechanics with Matlab

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

global cx cy cz vcon ucon caxis

den1 = (cz + vcon * z) ^ 2;

den2 = (cz + ucon * z) ^ 2;

fz = (1.0 + (cy / den1) + (cx / den2)) * z ^ 2 - caxis ^ 2;


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: various distances and Multi-Level Parametric Image Registration
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             {'SSD','NCC','MI','NGF'}
%   - transformation       affine2D
% see also E7_PETCT_MLPIR_ext
%==============================================================================

clear, close all, help(mfilename);
setup2DPETCTData;                                         % load data

% a list of distance measures to be used
DM = {'SSD','NCC','MI','NGF'};

OPTpara = FAIRcell2struct(optPara('PIR-GN'));

for dm = 1:length(DM), % run over all distance measures

  % initialize interpolation, using a smooth representation (theta=1e0)
  imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e0);

  % initialize transformation, create initial guess and reference for stopping
  trafo('reset','trafo','affine2D'); wStop = trafo('w0'); w0 = wStop;

  % initialize distance and display options
  distance('reset','distance',DM{dm}); distance('disp')

  % run MLPIR using sufficient amount of details (level=5)
  wSmooth =  MLPIR(ML,'minLevel',5,'plotIter',0,'plotMLiter',0);

  % refine interpolation (theta=1e-3)
  imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e-3);
  level = length(ML); omega = ML{level}.omega; m = ML{level}.m;
  [T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);

  % start PIR, using the result from the smooth problem as starting guess

  % initialize plots
  FAIRplots('set','mode','PIR');
  FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

  % optimize
  xc = getCellCenteredGrid(omega,m);
  Rc = imgModel(R,omega,xc);
  fctn = @(wc) PIRobjFctn(T,Rc,omega,m,0,[],[],xc,wc); fctn([]);
  [wc,his] = GaussNewton(fctn,w0,OPTpara{:});

  % visualize results
  yc = trafo(wc,xc);
  R0 = imgModel(R,omega,xc);
  T0 = imgModel(T,omega,xc);
  Tc = imgModel(T,omega,yc);

  figure(11); clf;
  viewImage(T0,omega,m,'axis','off'); hold on;
  plotGrid(yc,omega,m,'spacing',ceil(m/32),'linewidth',2,'color','w');

  figure(12); clf;
  overlayImage2D(Tc,R0,omega,m); axis off;
end;

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