text:
%%*******************************************************************
%% Prod2: compute the block diagonal matrix A*B
%%
%%    C = Prod2(blk,A,B,options);
%%
%% INPUT:  blk = a cell array describing the block structure of A and B
%%         A,B = square matrices or column vectors.
%%
%%   options = 0  if no special structure
%%             1  if C is symmetric
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

   function  C = Prod2(blk,A,B,options);

   global spdensity

   if (nargin == 3); options = 0; end;
   iscellA = iscell(A); iscellB = iscell(B);
%%
   if (~iscellA & ~iscellB)
      if (size(blk,1) > 1);
         error('Prod2: blk and A,B are not compatible');
      end;
      if strcmp(blk{1},'s')
         numblk = length(blk{2});
         isspA = issparse(A); isspB = issparse(B);
         if (numblk > 1)
            if ~isspA; A=sparse(A); isspA=1; end
            if ~isspB; B=sparse(B); isspB=1; end
         end
         %%use_matlab = (options==0 & ~isspA & ~isspB) | (isspA & isspB);
         use_matlab = (~isspA & ~isspB) | (isspA & isspB);
         if (use_matlab)
            C = A*B;
            if (options==1); C = 0.5*(C+C'); end;
         else
            C = mexProd2(blk,A,B,options);
         end
         checksparse = (numblk==1) & (isspA | isspB);
         if (checksparse)
            n2 = sum(blk{2}.*blk{2});
            if (mexnnz(C) <= spdensity*n2);
               if ~issparse(C); C = sparse(C); end;
            else
               if issparse(C); C = full(C); end;
            end
         end
      elseif (strcmp(blk{1},'q') | strcmp(blk{1},'l') | strcmp(blk{1},'u'))
         C = A.*B;
      end
   else
      error(['Prod2: A,B must be matrices']);
   end
%%*******************************************************************

repo: REST


text:
function [Y, s] = standardize(X)
% Unitize the vectors to be unit length
%   By default dim = 1 (columns).
% Written by Mo Chen (sth4nth@gmail.com).
X = bsxfun(@minux,X,mean(X,2));
s = sqrt(mean(sum(X.^2,1)));
Y = X/s;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
 function tdir = test_dir(tdir)
%function tdir = test_dir(tdir)
%
% Return path to test directory.
% Users can/should customize this file.

persistent pdir
if ~isvar('pdir') || isempty(pdir)
	try
		pdir.name = getenv('USER');
	catch
		pdir.name = 'user';
	end
	pdir.name = ['test,'  pdir.name];
	pdir.top = '/tmp/';
	pdir.path = [pdir.top pdir.name '/']; % trailing '/' is essential!
end

if nargin
	pdir.path = tdir;
end

if ~exist(pdir.path, 'dir')
	printf('Directory "%s" not found, so it will be created', pdir.path)
%	printf('run the following command!:')
%	'!mkdir /tmp/test'
	try
		mkdir(pdir.top, pdir.name) % compat with matlab 6.x
	catch
		try
			mkdir(pdir.path) % matlab 7.x accepts full path
		catch
			error 'give up'
		end
	end
end

tdir = pdir.path;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function colorTest()

%   colorTest() creates a test image showing the color encoding scheme

%   According to the c++ source code of Daniel Scharstein
%   Contact: schar@middlebury.edu

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 20:22:10 (Wed, 31 Oct 2006) $

% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%
%                         All Rights Reserved
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

%% test color pattern of Daniel's c++ code

truerange = 1;
height = 151;
width  = 151;
range = truerange * 1.04;

s2 = round(height/2);

[x y] = meshgrid(1:width, 1:height);

u = x*range/s2 - range;
v = y*range/s2 - range;

img = computeColor(u/truerange, v/truerange);

img(s2,:,:) = 0;
img(:,s2,:) = 0;

figure;
imshow(img);
title('test color pattern');
pause; close;

% test read and write flow
F(:,:,1) = u;
F(:,:,2) = v;
writeFlowFile(F, 'colorTest.flo');
F2 = readFlowFile('colorTest.flo');

u2 = F2(:,:,1);
v2 = F2(:,:,2);

img2 = computeColor(u2/truerange, v2/truerange);

img2(s2,:,:) = 0;
img2(:,s2,:) = 0;

figure; imshow(img2);
title('saved and reloaded test color pattern');
pause; close;

% color encoding scheme for optical flow
img = computeColor(u/range/sqrt(2), v/range/sqrt(2));

img(s2,:,:) = 0;
img(:,s2,:) = 0;

figure;
imshow(img);
title('optical flow color encoding scheme');
pause; close;

author: kristinbranson

repo: JAABA


text:
function ft_destroy_buffer

% FT_DESTROY_BUFFER stops the thread with the TCP server attached to
% the local MATLAB instance and removes all data from memory.
%
% Use as
%   ft_destroy_buffer
%
% See also FT_CREATE_BUFFER

% Copyright (C) 2010, 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$

% clearing the mex file from memory will cause the function registered with
% mexAtExit to be executed. This function will then stop the threads and
% release all allocated memory to the operating system
clear buffer

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function R = wishrand(S,nu);
%WISHRND Random matrices from Wishard distribution.
%   R = WISHRAND(S,N) returns a matrix of random numbers chosen
%   from the central Wishard distribution with parameters S and NU.
%
%   S is a symmetric positive definite scale matrix
%   NU is degrees of freedom
%
%   Note: E[R]=S
%
%   References: Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin (2013).
%                 Bayesian Data Analysis, third edition.
%               Gentle (2003), Random Number Generation and Monte Carlo
%                 Methods, 2nd ed, Springer, p. 199, Algorithm 5.8.
%               Smith & Hocking (1972), Algorithm AS 53: Wishard Variate
%                 Generator, Applied Statistics, 21(3), pp. 341-345.
%
%	See also INVWISHRAND
%
% Copyright (c) 1999-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 2
  error('Requires two input arguments.');
end;

[d d2] = size(S);
if d ~= d2
  error('Matrix S must be square');
end

% Algorithm 5.8 step 0.
[T,p]=chol(S);
if p > 0
  error('Matrix S must be positive definite.');
end

if (nu >= d) && (nu == round(nu))
  % distribution is proper and degrees of freedom is integer
  % brute-force may be used, which is surprisingly faster in Matlab
  % at least up to d>10, nu>1000
  % Algorithm described e.g. in (Gelman et al., 2013)
  Y = T'*randn(d,nu);
  R = Y*Y'./nu;
else
  % distribution is not proper or degrees of freedom is not integer
  % (Gentle, 2003) Algorithm 5.8
  % Algorithm 5.8 step 1.
  Z=zeros(d);
  for j=2:d
    Z(1:j,j)=randn(j,1);
  end
  % Algorithm 5.8 step 2. Note that there is error in the book.
  % Book says nu-i, while it should be nu+1-i
  % See errata <http://www.scs.gmu.edu/~jgentle/rngbk/errata.htm>
  % gamrand below is same as chi2rnd(nu-(1:d)+1), but much faster using c-code
  y=gamrand((nu+1-(1:d)),nu+1-(1:d));
  % Algorithm 5.8 step 3.
  Z2=Z.^2;
  sy=sqrt(y);
  B=zeros(d);
  B(1,1)=y(1);
  % In Matlab 6.5 following loop is accelerated and thus quite fast
  for j=2:d
    B(j,j)=y(j)+sum(Z2(1:(j-1),j));
    B(1,j)=Z(1,j).*sy(1);
    B(j,1)=B(1,j);
    for i=2:(j-1)
      B(i,j)=Z(i,j).*sy(i)+sum(Z(1:(i-1),i).*Z(1:(i-1),j));
      B(j,i)=B(i,j);
    end
  end
  % Algorithm 5.8 step 4.
  R=T'*B*T/nu;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
classdef nnhuberloss < nntest
  methods (Test)

    function basic(test)
      x = test.randn([100 1]) ;
      t = test.randn([100 1]) ;
      y = vl_nnhuberloss(x, t) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nnhuberloss(x, t, dzdy) ;
      test.der(@(x) vl_nnhuberloss(x, t), x, dzdy, dzdx, 1e-3*test.range) ;
    end

    function basicWithInstanceWeights(test)
      x = test.randn([100 1]) ;
      t = test.randn([100 1]) ;
      w = test.randn([100 1]) ;
      y = vl_nnhuberloss(x, t, 'instanceWeights', w) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nnhuberloss(x, t, dzdy, 'instanceWeights', w) ;
      test.der(@(x) vl_nnhuberloss(x, t, 'instanceWeights', w), x, ...
                                        dzdy, dzdx, 1e-3*test.range) ;
    end

  end
end

repo: MFT


text:
function [avgArray] = calculateAvgArray(tau, sPeriod, readings)
%This function uses raw readings and returns an array of readings for the
%input Tau argument. readings can either be freq or time

avgCount = tau / sPeriod; %gets count for averaging
loopCount = numel(readings) / avgCount; %get number of loop iterations needed
loopCount = floor(loopCount); %convert loopCount to integer in case it is a non int

avgArray = zeros(1,loopCount); %allocate array
iter = 0;
iter = int32(iter);

%loop to build array
for i = 1:loopCount
    temp = mean(readings((iter+1):(iter+avgCount)));
    avgArray(i) = temp;
    iter = iter + avgCount;
end

author: Sable

repo: mcbench-benchmarks


text:
function B = blockrot90(A, p, q, k)
%BLOCKROT90 Perform 90 degree block rotation on arrays.
%
%   BLOCKROT90(A,P,Q) is the 90 degree anticlockwise rotation of all P-by-Q
%   submatrices of the matrix A.
%
%   BLOCKROT90(A,P,Q,K) is the K*90 degree rotation of all P-by-Q submatrices
%   of A, K = +-1,+-2,...
%
%   BLOCKROT90 is vectorized along higher dimensions, so the same operations is
%   performed on all 2-D slices A(:,:,I,J,K,...).
%
%   See also ROT90.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 13:25:20 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   nargsin = nargin;
   error(nargchk(3, 4, nargsin));

   z = size(A);
   if nargsin < 4
      k = 1;
   else
      if any(size(k) ~= 1)
         error('K must be a scalar.');
      end
      k = k - 4*floor(k/4);     % map k to {0,1,2,3}
   end

   r = z(1)/p;
   s = z(2)/q;
   if (r ~= round(r)) | (s ~= round(s))
      error('A can not be subdivided into P-by-Q submatrices.');
   end

   u = [ p r q prod(z)/(z(1)*q) ];
   v = [ q*r p*s z(3:end) ];

   switch k
      case 1                    % 90 degrees anticlockwise
         B = reshape(A, u);
         B = reshape(permute(B(:,:,q:-1:1,:), [ 3 2 1 4 ]), v);
      case 2                    % 180 degrees
         B = reshape(A, u);
         B = reshape(B(p:-1:1,:,q:-1:1,:), z);
      case 3                    % 90 degrees clockwise
         B = reshape(A, u);
         B = reshape(permute(B(p:-1:1,:,:,:), [ 3 2 1 4 ]), v);
      otherwise                 % 0 degrees
         B = A;
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function w=weight2(n,m) % n-length of weighting vector, m - is alpha
% quadratic linguistic quantifier: Q(r)=(1/(1-alpha*(r)^0.5))
re=[];
for h=1:n
    re=[re ((1-m*(h/n)^0.5)^-1 - (1-m*((h-1)/n)^0.5)^-1)];
end
reT=sum(re);
w=re/reT;
end


text:
function [ state ] = getPin13State( t )
%% This function is used to get the state of the input pin 13
%     on the media flange, of the KUKA iiwa 7 R 800.

%% Syntax:
% [ state ] = getPin13State( t )

%% About:
% This function is used to get the state of the input pin number 13
%  on the media flange, of the KUKA iiwa 7 R 800.

%% Arreguments:
% t: is the TCP/IP connection.
% state:  represents the state of the input-pin. It is an integer equal to
% 1 when the pin state is true, or it is an integer equal to 0 when the
% pin state is false. Otherwise, if the pin is not connected the return
% value is empty.

% Copyright, Mohammad SAFEEA, 9th of May 2017

theCommand='getPin13';
fprintf(t, theCommand);
message=fgets(t);
state=str2num(message);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [seedinfo,outstate] = ncorr_gui_seedpreview(reference,current,roi,num_region,pos_seed,radius,cutoff_diffnorm,cutoff_iteration,enabled_stepanalysis,subsettrunc,num_img,total_imgs,pos_parent)
% This is a GUI for displaying the seeds before accepting them. Its
% possible the seeds are not correct or do not converge, so I included this
% utility as a check.
%
% Inputs -----------------------------------------------------------------%
%   reference - ncorr_class_img; used for displaying the background
%   image and calculations
%   current -  ncorr_class_img(s); used for displaying the background
%   image and calculations
%   roi - ncorr_class_roi; ROI corresponding to the reference image
%   num_region - integer; number corresponding to the region
%   pos_seed - integer array; 1st column is the x-position of the seed; 2nd
%   is the y-position. Each row corresponds to a thread in ascending order
%   (i.e. the first row is the first thread).
%   radius - integer; radius of subset
%   cutoff_diffnorm - double; cutoff of norm of the difference vector
%   cutoff_iteration - integer; cutoff for number of iterations
%   enabled_stepanalysis - logical; if true, then process as many seeds as
%   possible. If false, process all the seeds.
%   subsettrunc - logical; if true, then enable subset truncation
%   num_img - integer; number of reference image being analyzed
%   total_imgs - integer; total number of images being analyzed
%   pos_parent - integer array; this is the position of the parent figure
%   which determines where to position this figure
%
% Outputs ----------------------------------------------------------------%
%   seedinfo - struct; contains struct('paramvector',{},'num_region',{},'num_thread',{},'computepoints',{})
%   which is information on the location of the seed, the deformation
%   parameters for that seed, the region number, and the thread assigned for
%   computation.
%   outstate - integer; returns either out.cancelled, out.failed, or
%   out.success.
%
% Note that if step analysis is enabled, outstate will return success if at
% least one image is seeded. If no seeds are found, it will return failed.
% However, if step analysis is disabled, outstate will return success only
% if all images are seeded correctly, or else it will return failed.

    % Data ---------------------------------------------------------------%
    % Initialize Outputs
    outstate = out.cancelled;
    seedinfo = struct('paramvector',{},'num_region',{},'num_thread',{},'computepoints',{}); % paramvector = [x y u v du/dx du/dy dv/dx dv/dy corrcoef]
    % Get GUI handles
    handles_gui = init_gui();
    % Run c-tor
    feval(ncorr_util_wrapcallbacktrycatch(@constructor,handles_gui.figure));

    % Callbacks ----------------------------------------------------------%
    function constructor()
        % Calculate seeds. ncorr_alg_seedanalysis will return failed if no
        % images are seeded correctly in the case that step analysis is
        % enabled. If step analysis is disabled, it will return failed if
        % any of the seeds are processed incorrectly.
        [seedinfo_prelim,convergence_prelim,outstate_seeds] = ncorr_alg_seedanalysis(reference, ...
                                                                                current, ...
                                                                                roi, ...
                                                                                num_region, ...
                                                                                pos_seed, ...
                                                                                radius, ...
                                                                                cutoff_diffnorm, ...
                                                                                cutoff_iteration, ...
                                                                                enabled_stepanalysis, ...
                                                                                subsettrunc, ...
                                                                                num_img, ...
                                                                                total_imgs);

        if (outstate_seeds == out.success)
            % NOTE: Seeds are ordered in seedinfo_prelim WRT to thread
            % number when processed by ncorr_alg_seedanalysis (i.e.
            % seedinfo_prelim(i) == thread number i), so this isnt
            % explicitly checked when displaying subsets corresponding
            % to different threads.

            % Check to see if any seeds have a relatively high
            % correlation coefficient and warn the user about it.
            paramvectorcollection = vertcat(seedinfo_prelim.paramvector); % The ninth column is the correlation coefficients
            if (any(paramvectorcollection(:,9) > 0.4))
                h_error = errordlg('High correlation coefficient detected. One of the seed points is possibly incorrect. Proceed with caution.','Error','modal');
                uiwait(h_error);
            end

            % Precompute cirrois
            cirrois = struct('mask',{},'region',{},'boundary',{},'x',{},'y',{},'radius',{});
            for i = 0:size(pos_seed,1)-1
                cirrois(i+1) = roi.get_cirroi(pos_seed(i+1,1),pos_seed(i+1,2),radius,subsettrunc);
            end

            % Precompute the subsets, so they do not need to be
            % recalculated when switching between images/threads.
            % Initialize
            preview_subsets_ref = zeros(cirrois(1).radius*2+1,cirrois(1).radius*2+1,length(cirrois));
            preview_subsets_cur = zeros(cirrois(1).radius*2+1,cirrois(1).radius*2+1,length(cirrois),size(seedinfo_prelim,3));

            % Interpolate
            gs_buffer_ref = reference.get_gs();
            for i = 0:size(seedinfo_prelim,3)-1
                % Getting bcoefficients could be time consuming if a
                % lot of high resolution images are being processed.
                bcoef_buffer_cur = current(i+1).get_bcoef();
                for j = 0:size(seedinfo_prelim,1)-1
                    for k = 0:size(cirrois(j+1).region.noderange,1)-1
                        x_ref = k+(cirrois(j+1).x-cirrois(j+1).radius);
                        for l = 0:2:cirrois(j+1).region.noderange(k+1)-1
                            % Transform x_ref and y_ref and then
                            % interpolate
                            vec_y_ref = (cirrois(j+1).region.nodelist(k+1,l+1):cirrois(j+1).region.nodelist(k+1,l+2))';
                            vec_x_ref = repmat(x_ref,length(vec_y_ref),1);

                            % Fill reference subset - only do this once
                            if (i == 0)
                                preview_subsets_ref(vec_y_ref-(cirrois(j+1).y-cirrois(j+1).radius)+1,k+1,j+1) = gs_buffer_ref(vec_y_ref+1,x_ref+1);
                            end

                            % Transform ref coordinates
                            vec_x_cur = vec_x_ref + seedinfo_prelim(j+1,1,i+1).paramvector(3) + seedinfo_prelim(j+1,1,i+1).paramvector(5)*(vec_x_ref-cirrois(j+1).x) + seedinfo_prelim(j+1,1,i+1).paramvector(6)*(vec_y_ref-cirrois(j+1).y);
                            vec_y_cur = vec_y_ref + seedinfo_prelim(j+1,1,i+1).paramvector(4) + seedinfo_prelim(j+1,1,i+1).paramvector(7)*(vec_x_ref-cirrois(j+1).x) + seedinfo_prelim(j+1,1,i+1).paramvector(8)*(vec_y_ref-cirrois(j+1).y);

                            % Interpolate
                            vec_interp_cur = ncorr_alg_interpqbs([vec_x_cur vec_y_cur],bcoef_buffer_cur,0,0,current(i+1).border_bcoef);

                            % Note that any values interpolated out of
                            % the range of the bcoefficient matrix in
                            % ncorr_alg_interpqbs is set to NaN.
                            % Set these NaNs to zero for display
                            % purposes.
                            vec_interp_cur(isnan(vec_interp_cur)) = 0;

                            % Fill current subset
                            preview_subsets_cur(vec_y_ref-cirrois(j+1).y+cirrois(j+1).radius+1,k+1,j+1,i+1) = vec_interp_cur;
                        end
                    end
                end
            end

            % Set Data
            setappdata(handles_gui.figure,'cirrois',cirrois);
            setappdata(handles_gui.figure,'num_cur',size(seedinfo_prelim,3)-1);
            setappdata(handles_gui.figure,'num_thread',size(pos_seed,1)-1);
            setappdata(handles_gui.figure,'preview_subsets_ref',preview_subsets_ref);
            setappdata(handles_gui.figure,'preview_subsets_cur',preview_subsets_cur);
            setappdata(handles_gui.figure,'seedinfo_prelim',seedinfo_prelim);
            setappdata(handles_gui.figure,'convergence_prelim',convergence_prelim);

            % Update
            update_axes('set');

            % Set Visible
            set(handles_gui.figure,'Visible','on');
        elseif (outstate_seeds == out.failed)
            h_error = errordlg('Seed placement failed; please replace seeds.','Error','modal');
            uiwait(h_error);

            outstate = out.failed;

            % Exit
            close(handles_gui.figure);
        else
            % Analysis was cancelled - exit.
            close(handles_gui.figure);
        end
    end

    function callback_edit_threadnum(hObject,eventdata) %#ok<INUSD>
        % Get data
        num_thread = getappdata(handles_gui.figure,'num_thread');

        num_thread_prelim = str2double(get(handles_gui.edit_threadnum,'string'));
        if (ncorr_util_isintbb(num_thread_prelim,1,size(pos_seed,1),'Thread number') == out.success)
            num_thread = num_thread_prelim-1;
        end

        % Set data
        setappdata(handles_gui.figure,'num_thread',num_thread);

        % Update
        update_axes('set');
    end

    function callback_edit_imgnum(hObject,eventdata) %#ok<INUSD>
        % Get data
        num_cur = getappdata(handles_gui.figure,'num_cur');
        seedinfo_prelim = getappdata(handles_gui.figure,'seedinfo_prelim');

        num_cur_prelim = str2double(get(handles_gui.edit_imgnum,'string'));
        if (ncorr_util_isintbb(num_cur_prelim,1,size(seedinfo_prelim,3),'Current image number') == out.success)
            num_cur = num_cur_prelim-1;
        end

        % Set data
        setappdata(handles_gui.figure,'num_cur',num_cur);

        % Update
        update_axes('set');
    end

    function callback_button_finish(hObject,eventdata) %#ok<INUSD>
        % Get data
        seedinfo_prelim = getappdata(handles_gui.figure,'seedinfo_prelim');

        % Set output
        for i = 0:size(seedinfo_prelim,3)-1
            for j = 0:size(seedinfo_prelim,1)-1
                seedinfo(j+1,1,i+1) = seedinfo_prelim(j+1,1,i+1);
            end
        end
        outstate = out.success;

        % Exit
        close(handles_gui.figure);
    end

    function callback_button_cancel(hObject,eventdata) %#ok<INUSD>
        close(handles_gui.figure);
    end

    function callback_button_left_thread(hObject,eventdata) %#ok<INUSD>
        % Get data
        num_thread = getappdata(handles_gui.figure,'num_thread');

        % Check for overshoot
        if (num_thread > 0)
            num_thread = num_thread-1;
        end

        % Set data
        setappdata(handles_gui.figure,'num_thread',num_thread);

        % Update
        update_axes('set');
    end

    function callback_button_right_thread(hObject,eventdata) %#ok<INUSD>
        % Get data
        num_thread = getappdata(handles_gui.figure,'num_thread');

        % Check for overshoot
        if (num_thread < size(pos_seed,1)-1)
            num_thread = num_thread+1;
        end

        % Set data
        setappdata(handles_gui.figure,'num_thread',num_thread);

        % Update
        update_axes('set');
    end

    function callback_button_left_cur(hObject,eventdata) %#ok<INUSD>
        % Get data
        num_cur = getappdata(handles_gui.figure,'num_cur');

        % Check for overshoot
        if (num_cur > 0)
            num_cur = num_cur-1;
        end

        % Set data
        setappdata(handles_gui.figure,'num_cur',num_cur);

        % Update
        update_axes('set');
    end

    function callback_button_right_cur(hObject,eventdata) %#ok<INUSD>
        % Get data
        num_cur = getappdata(handles_gui.figure,'num_cur');
        seedinfo_prelim = getappdata(handles_gui.figure,'seedinfo_prelim');

        % Check for overshoot
        if (num_cur < size(seedinfo_prelim,3)-1)
            num_cur = num_cur+1;
        end

        % Set data
        setappdata(handles_gui.figure,'num_cur',num_cur);

        % Update
        update_axes('set');
    end

    function update_axes(action)
        % Get data
        num_cur = getappdata(handles_gui.figure,'num_cur');
        num_thread = getappdata(handles_gui.figure,'num_thread');
        preview_subsets_ref = getappdata(handles_gui.figure,'preview_subsets_ref');
        preview_subsets_cur = getappdata(handles_gui.figure,'preview_subsets_cur');
        cirrois = getappdata(handles_gui.figure,'cirrois');
        seedinfo_prelim = getappdata(handles_gui.figure,'seedinfo_prelim');
        convergence_prelim = getappdata(handles_gui.figure,'convergence_prelim');

        if (strcmp(action,'set'))
            gs_buffer_ref = reference.get_gs();
            gs_buffer_cur = current(num_cur+1).get_gs();

            % Paint reference subset and subset locations
            preview_subsetloc_ref = gs_buffer_ref;
            preview_subsetloc_cur = gs_buffer_cur;

            % Iterate over num_thread last so that it highlights this
            % subset over the others
            for i = [0:num_thread-1 num_thread+1:size(pos_seed,1)-1 num_thread]
                for j = 0:size(cirrois(i+1).region.noderange,1)-1
                    x_ref = j+(cirrois(i+1).x-cirrois(i+1).radius);
                    for k = 0:2:cirrois(i+1).region.noderange(j+1)-1
                        for l = cirrois(i+1).region.nodelist(j+1,k+1):cirrois(i+1).region.nodelist(j+1,k+2)
                            y_ref = l;

                            % Fill subset locations
                            if (i == num_thread)
                                % User brighter highlight to indicate this is
                                % the seed corresponding to the current
                                % threads
                                preview_subsetloc_ref(y_ref+1,x_ref+1) = gs_buffer_ref(y_ref+1,x_ref+1)+1.75*reference.max_gs;
                            else
                                % Do regular highlight
                                preview_subsetloc_ref(y_ref+1,x_ref+1) = gs_buffer_ref(y_ref+1,x_ref+1)+reference.max_gs;
                            end

                            % defvector_init(1) = u
                            % defvector_init(2) = v
                            % defvector_init(3) = du_dx
                            % defvector_init(4) = du_dy
                            % defvector_init(5) = dv_dx
                            % defvector_init(6) = dv_dy
                            % y_cur = y_ref + v + dv/dx*deltax + dv/dy*deltay
                            % x_cur = x_ref + u + du/dx*deltax + du/dy*deltay

                            x_cur = x_ref + seedinfo_prelim(i+1,1,num_cur+1).paramvector(3) + seedinfo_prelim(i+1,1,num_cur+1).paramvector(5)*(x_ref-cirrois(i+1).x) + seedinfo_prelim(i+1,1,num_cur+1).paramvector(6)*(y_ref-cirrois(i+1).y);
                            y_cur = y_ref + seedinfo_prelim(i+1,1,num_cur+1).paramvector(4) + seedinfo_prelim(i+1,1,num_cur+1).paramvector(7)*(x_ref-cirrois(i+1).x) + seedinfo_prelim(i+1,1,num_cur+1).paramvector(8)*(y_ref-cirrois(i+1).y);

                            % x_cur and y_cur could be out of bounds, check
                            % before using them
                            if (floor(x_cur) >= 0 && floor(y_cur) >= 0 && ceil(x_cur) < size(gs_buffer_cur,2) && ceil(y_cur) < size(gs_buffer_cur,1))
                                % Fill subset locations. This is a crude
                                % approximation of the current subset. Use both
                                % floor and ceil of transformed locations to
                                % highlight the area of the current subset.
                                if (i == num_thread)
                                    % User brighter highlight to indicate this is
                                    % the seed corresponding to the thread of
                                    % interest
                                    preview_subsetloc_cur(floor(y_cur)+1,floor(x_cur)+1) = gs_buffer_cur(floor(y_cur)+1,floor(x_cur)+1)+1.75*current(num_cur+1).max_gs;
                                    preview_subsetloc_cur(ceil(y_cur)+1,ceil(x_cur)+1) = gs_buffer_cur(ceil(y_cur)+1,ceil(x_cur)+1)+1.75*current(num_cur+1).max_gs;
                                else
                                    % Do regular highlight
                                    preview_subsetloc_cur(floor(y_cur)+1,floor(x_cur)+1) = gs_buffer_cur(floor(y_cur)+1,floor(x_cur)+1)+current(num_cur+1).max_gs;
                                    preview_subsetloc_cur(ceil(y_cur)+1,ceil(x_cur)+1) = gs_buffer_cur(ceil(y_cur)+1,ceil(x_cur)+1)+current(num_cur+1).max_gs;
                                end
                            end
                        end
                    end
                end
            end

            % Set Images
            imshow(preview_subsets_ref(:,:,num_thread+1),[reference.min_gs reference.max_gs],'Parent',handles_gui.axes_refsubset);
            set(handles_gui.axes_refsubset,'Visible','off');
            imshow(preview_subsetloc_ref,[reference.min_gs 2*reference.max_gs],'Parent',handles_gui.axes_refsubsetloc);
            set(handles_gui.axes_refsubsetloc,'Visible','off');

            imshow(preview_subsets_cur(:,:,num_thread+1,num_cur+1),[current(num_cur+1).min_gs current(num_cur+1).max_gs],'Parent',handles_gui.axes_cursubset);
            set(handles_gui.axes_cursubset,'Visible','off');
            imshow(preview_subsetloc_cur,[current(num_cur+1).min_gs 2*current(num_cur+1).max_gs],'Parent',handles_gui.axes_cursubsetloc);
            set(handles_gui.axes_cursubsetloc,'Visible','off');

            % Set Texts
            if (convergence_prelim(num_thread+1,1,num_cur+1).num_iterations < cutoff_iteration)
                set(handles_gui.text_num_iterations_num,'String',num2str(convergence_prelim(num_thread+1,1,num_cur+1).num_iterations,'%6.0f'),'ForegroundColor', 'k');
            else
                set(handles_gui.text_num_iterations_num,'String',[num2str(convergence_prelim(num_thread+1,1,num_cur+1).num_iterations,'%6.0f') ' (Max)' ],'ForegroundColor', 'r');
            end
            set(handles_gui.text_gradnorm_num,'String',num2str(convergence_prelim(num_thread+1,1,num_cur+1).diffnorm));
            set(handles_gui.text_corrcoef_num,'String',num2str(seedinfo_prelim(num_thread+1,1,num_cur+1).paramvector(9)));
            set(handles_gui.text_thread,'String', ['Thread: ' num2str(num_thread+1)]);
            set(handles_gui.text_cur,'String', ['Name: ' current(num_cur+1).name(1:min(end,22))]);
            set(handles_gui.text_ref,'String', ['Name: ' reference.name(1:min(end,22))]);

            % Set buttons
            % Set left/right buttons for current image
            set(handles_gui.edit_imgnum,'String',num2str(num_cur+1));
            if (size(seedinfo_prelim,3) == 1)
                set(handles_gui.button_right_cur,'Enable','off');
                set(handles_gui.button_left_cur,'Enable','off');
                set(handles_gui.edit_imgnum,'Enable','off');
            elseif (num_cur == 0)
                set(handles_gui.button_right_cur,'Enable','on');
                set(handles_gui.button_left_cur,'Enable','off');
                set(handles_gui.edit_imgnum,'Enable','on');
            elseif (num_cur == size(seedinfo_prelim,3)-1)
                set(handles_gui.button_right_cur,'Enable','off');
                set(handles_gui.button_left_cur,'Enable','on');
                set(handles_gui.edit_imgnum,'Enable','on');
            else
                set(handles_gui.button_right_cur,'Enable','on');
                set(handles_gui.button_left_cur,'Enable','on');
                set(handles_gui.edit_imgnum,'Enable','on');
            end

            % Set left/right buttons for thread number
            set(handles_gui.edit_threadnum,'String',num2str(num_thread+1));
            if (size(pos_seed,1) == 1)
                set(handles_gui.button_right_thread,'Enable','off');
                set(handles_gui.button_left_thread,'Enable','off');
                set(handles_gui.edit_threadnum,'Enable','off');
            elseif (num_thread == 0)
                set(handles_gui.button_right_thread,'Enable','on');
                set(handles_gui.button_left_thread,'Enable','off');
                set(handles_gui.edit_threadnum,'Enable','on');
            elseif (num_thread == size(pos_seed,1)-1)
                set(handles_gui.button_right_thread,'Enable','off');
                set(handles_gui.button_left_thread,'Enable','on');
                set(handles_gui.edit_threadnum,'Enable','on');
            else
                set(handles_gui.button_right_thread,'Enable','on');
                set(handles_gui.button_left_thread,'Enable','on');
                set(handles_gui.edit_threadnum,'Enable','on');
            end
        end
    end

    function handles_gui = init_gui()
    % GUI controls -------------------------------------------------------%
        % Figure
        handles_gui.figure = figure( ...
            'Tag', 'figure', ...
            'Units', 'characters', ...
            'Position', ncorr_util_figpos(pos_parent,[56.1 167]), ...
            'Name', 'Seed Preview', ...
            'MenuBar', 'none', ...
            'NumberTitle', 'off', ...
            'Color', get(0,'DefaultUicontrolBackgroundColor'), ...
            'handlevisibility','off', ...
            'DockControls','off', ...
            'WindowStyle','modal', ...
            'Resize','off', ...
            'Visible','off', ...
            'IntegerHandle','off', ...
            'Interruptible','off');

        % Panels
        handles_gui.group_menu = uibuttongroup( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_menu', ...
            'Units', 'characters', ...
            'Position', [2 48.7 35 6.7], ...
            'Title', 'Menu', ...
            'Interruptible','off');

        handles_gui.group_ref = uibuttongroup( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_reference', ...
            'Units', 'characters', ...
            'Position', [39.0 0.75 62 54.6], ...
            'Title', 'Reference', ...
            'Interruptible','off');

        handles_gui.group_cur = uibuttongroup( ...
            'Parent', handles_gui.figure, ...
            'Tag', 'group_current', ...
            'Units', 'characters', ...
            'Position', [103 0.75 62 54.6], ...
            'Title', 'Current', ...
            'Interruptible','off');

        % Axes
        handles_gui.axes_refsubsetloc = axes( ...
            'Parent', handles_gui.group_ref, ...
            'Tag', 'axes_refsubsetloc', ...
            'Units', 'characters', ...
            'Visible', 'off', ...
            'Position', [2.1 35.7 56.5 17], ...
            'Interruptible','off');

        handles_gui.axes_cursubsetloc = axes( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'axes_cursubsetloc', ...
            'Units', 'characters', ...
            'Visible', 'off', ...
            'Position', [2.1 35.7 56.5 17], ...
            'Interruptible','off');

        handles_gui.axes_refsubset = axes( ...
            'Parent', handles_gui.group_ref, ...
            'Tag', 'axes_refsubset', ...
            'Units', 'characters', ...
            'Visible', 'off', ...
            'Position', [2.1 13.1 56.5 17], ...
            'Interruptible','off');

        handles_gui.axes_cursubset = axes( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'axes_cursubset', ...
            'Units', 'characters', ...
            'Visible', 'off', ...
            'Position', [2.1 13.1 56.5 17], ...
            'Interruptible','off');

        % Static Texts
        handles_gui.text_refsubloc = uicontrol( ...
            'Parent', handles_gui.group_ref, ...
            'Tag', 'text_cur', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 33.6 56.5 1.2], ...
            'String', 'Location of Reference Subset', ...
            'Interruptible','off');

        handles_gui.text_cursubloc = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_cursubloc', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 33.6 56.5 1.2], ...
            'String', 'Approximate Location of Current Subset', ...
            'Interruptible','off');

        handles_gui.text_refsub = uicontrol( ...
            'Parent', handles_gui.group_ref, ...
            'Tag', 'text_refsub', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 11.0 56.5 1.2], ...
            'String', 'Reference Subset', ...
            'Interruptible','off');

        handles_gui.text_cursub = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_cursub', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 11.0 56.5 1.2], ...
            'String', 'Transformed Current Subset', ...
            'Interruptible','off');

        handles_gui.text_num_iterations = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_num_iterations', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 8.4 38.5 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', '# of Gauss Newton Iterations:', ...
            'Interruptible','off');

        handles_gui.text_gradnorm = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_gradnorm', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 6.6 38.5 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', 'Norm of difference vector:', ...
            'Interruptible','off');

        handles_gui.text_corrcoef = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_corrcoef', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 4.8 38.5 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', 'Correlation Coefficient:', ...
            'Interruptible','off');

        handles_gui.text_num_iterations_num = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_num_iterations_num', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [41.5 8.4 13 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', '', ...
            'Interruptible','off');

        handles_gui.text_gradnorm_num = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_gradnorm_num', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [41.5 6.6 13 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', '', ...
            'Interruptible','off');

        handles_gui.text_corrcoef_num = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_corrcoef_num', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [41.5 4.8 13 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', '', ...
            'Interruptible','off');

        handles_gui.text_thread = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_thread', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 3.0 56.9 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', 'Thread:', ...
            'Interruptible','off');

        handles_gui.text_cur = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'text_cur', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 1.2 56.9 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', 'Name:', ...
            'Interruptible','off');

        handles_gui.text_ref = uicontrol( ...
            'Parent', handles_gui.group_ref, ...
            'Tag', 'text_ref', ...
            'Style', 'text', ...
            'Units', 'characters', ...
            'Position', [2.1 1.2 56.9 1.2], ...
            'HorizontalAlignment', 'Left', ...
            'String', 'Name:', ...
            'Interruptible','off');

        % Edit
        handles_gui.edit_threadnum = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'edit_threadnum', ...
            'Style', 'edit', ...
            'Units', 'characters', ...
            'Position', [44.3 3.1 7 1.2], ...
            'String', '', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_edit_threadnum,handles_gui.figure), ...
            'Enable', 'off', ...
            'Interruptible','off');

        handles_gui.edit_imgnum = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'edit_imgnum', ...
            'Style', 'edit', ...
            'Units', 'characters', ...
            'Position', [44.3 1.2 7 1.2], ...
            'String', '', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_edit_imgnum,handles_gui.figure), ...
            'Enable', 'off', ...
            'Interruptible','off');

        % Pushbuttons
        handles_gui.button_finish = uicontrol( ...
            'Parent', handles_gui.group_menu, ...
            'Tag', 'button_finish', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [2.1 3 29.7 1.7], ...
            'String', 'Finish', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_finish,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_cancel = uicontrol( ...
            'Parent', handles_gui.group_menu, ...
            'Tag', 'button_cancel', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [2.1 1 29.7 1.7], ...
            'String', 'Cancel', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_cancel,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_left_thread = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'button_left_thread', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [37.3 2.9 6 1.5], ...
            'String', '<', ...
            'Enable', 'off', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_left_thread,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_right_thread = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'button_right_thread', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [52.3 2.9 6 1.5], ...
            'String', '>', ...
            'Enable', 'off', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_right_thread,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_left_cur = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'button_left_cur', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [37.3 1.1 6 1.5], ...
            'String', '<', ...
            'Enable', 'off', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_left_cur,handles_gui.figure), ...
            'Interruptible','off');

        handles_gui.button_right_cur = uicontrol( ...
            'Parent', handles_gui.group_cur, ...
            'Tag', 'button_right_cur', ...
            'Style', 'pushbutton', ...
            'Units', 'characters', ...
            'Position', [52.3 1.1 6 1.5], ...
            'String', '>', ...
            'Enable', 'off', ...
            'Callback', ncorr_util_wrapcallbacktrycatch(@callback_button_right_cur,handles_gui.figure), ...
            'Interruptible','off');
    end

    % Pause until figure is closed ---------------------------------------%
    waitfor(handles_gui.figure);
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
%GETSYSTEMTIMEPRECISEASFILETIME Returns system time with high precision.
% Usage:
%   t = GetSystemTimePreciseAsFileTime
%
% Returns:
%   t: seconds from 00:00:00 UTC, 1/1/1601 with 0.1 ns tick precision
%
% Example:
% >> a = GetSystemTimePreciseAsFileTime
% a =
%    1.3113e+10
% >> b = GetSystemTimePreciseAsFileTime
% b =
%    1.3113e+10
% >> b - a
% ans =
%     4.0448
%
% Reference:
%   https://msdn.microsoft.com/en-us/library/windows/desktop/hh706895(v=vs.85).aspx
%   https://msdn.microsoft.com/en-us/library/windows/desktop/dn553408(v=vs.85).aspx

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
% GP_COV_SUM - Covariance function as a sum of two covariance functions.
%
% Usage:
%
%   COVFUNC = GP_COV_SUM(COVFUNC1, COVFUNC2)
%
% The returned covariance function is called as
%
%   K = COVFUNC(THETA)
%   [K, DK] = COVFUNC(THETA)
%
% where THETA is a vector of parameters. The parameters of COVFUNC1 and
% COVFUNC2 are concatenated into a single parameter vector.
%
% If no arguments are given, dimensionality information is returned:
%
%   [N_THETA, N1, N2] = COVFUNC()
%
% N_THETA : Number of non-fixed parameters
% N1      : Number of rows in the covariance matrix
% N2      : Number of columns in the covariance matrix
%
% See also GP_COV, GP_COV_PRODUCT.

% Last modified 2010-01-25
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function func = gp_cov_sum(varargin)

func = @get_covariance;

if nargin < 2
  error('Must give at least two covariance functions')
end

covfuncs = varargin ;

% Number of covariance functions
n_funcs = nargin ;

% Number of parameters for each covariance function
n_theta = zeros(n_funcs,1) ;
% Dimensionalities of each covariance function
M = zeros(n_funcs,1) ;
N = zeros(n_funcs,1) ;
% Extract the values from the covariance functions
for i = 1:n_funcs
  [n_theta(i), M(i), N(i)] = varargin{i}();
end

% Indices of the hyperparameters for each covariance function
ind_theta = cell(n_funcs,1);
for i = 1:n_funcs
  ind_theta{i} = (1+sum(n_theta(1:(i-1)))):(sum(n_theta(1:i))) ;
end

%[n_theta1, M1, N1] = covfunc1();
%[n_theta2, M2, N2] = covfunc2();

if ~all(M == M(1)) || ~all(N == N(1))
  error('Can''t sum covariance matrices with different dimensionalities');
end
M = M(1) ;
N = N(1) ;

  function varargout = get_covariance(theta)

  if nargout == 0
    nout = 1 ;
  else
    nout = nargout ;
  end

  varargout = cell(nout,1);

  out = cell(nout,1);
  varargout = cell(nout,1);

  % Return only dimension information if requested
  if nargin == 0
    if nout >= 1
      varargout{1} = sum(n_theta); % number of parameters
      if nout >= 2
        varargout{2} = M; % dimensionalities
        if nout >= 3
          varargout{3} = N; % dimensionalities
        end
      end
    end
    return
  end

  varargout{1} = 0 ;
  varargout{2} = cell(sum(n_theta),1) ;

  for i = 1:n_funcs
    % Compute the covariance matrix
    [out{:}] = covfuncs{i}(theta(ind_theta{i}));
    varargout{1} = varargout{1} + out{1};

    % Compute the derivative
    if nout >= 2
      varargout{2}(ind_theta{i}) = out{2}(:);
    end
  end

  end

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function ppFlipper(data)
% ppFlipper(data)
%
% write data on the parallelport;
% It assumes the parallelport to be at IO_ADDR.
%

% kraulem 07/07

global IO_ADDR
ppWriteStay(IO_ADDR, data);
return;

author: PatternRecognition

repo: OpenBMI


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model/simulation parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sim=1;
OS_RATE = 8;
SNR = 1;
fc = 10e3/20e6; % sample rate is 20 MHz, top is 10 kHz offset
muFOC = floor(.01*2^12)/2^12;
muTOC = floor(.01*8*2^12)/2^12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize LUTs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
make_srrc_lut;
make_train_lut;
make_trig_lut;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (sim)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Emulate microprocessor packet creation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % data payload creation
    messageASCII = 'Hello World';
    message = double(unicode2native(messageASCII));
    % add on length of message to the front with four bytes
    msgLength = length(message);
    messageWithNumBytes =[ ...
        mod(msgLength,2^8) ...
        mod(floor(msgLength/2^8),2^8) ...
        mod(floor(msgLength/2^16),2^8) ...
        1 ... % message ID
        message];
    % add two bytes at the end, which is a CRC
    messageWithCRC = CreateAppend16BitCRC(messageWithNumBytes);
    ml = length(messageWithCRC);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FPGA radio transmit core
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    data_in = 0;
    empty_in = 1;
    tx_en_in = 0;
    numBytesFromFifo = 0;
    num_samp = ml*8*2*2*3;
    x = zeros(1,num_samp);
    CORE_LATENCY = 4;
    data_buf = zeros(1,CORE_LATENCY);
    empty_buf = ones(1,CORE_LATENCY);
    clear_buf = zeros(1,CORE_LATENCY);
    tx_en_buf = zeros(1,CORE_LATENCY);
    for i1 = 1:num_samp
        % first thing the processor does is clear the internal tx fifo
        if i1 == 1
            clear_fifo_in = 1;
        else
            clear_fifo_in = 0;
        end

        if i1 == 5 % wait a little bit then begin to load the fifo
            empty_in = 0;
            numBytesFromFifo = 0;
        end

        data_buf = [data_buf(2:end) data_in];
        empty_buf = [empty_buf(2:end) empty_in];
        clear_buf = [clear_buf(2:end) clear_fifo_in];
        tx_en_buf = [tx_en_buf(2:end) tx_en_in];
        [i_out, q_out, re_byte_out, tx_done_out, d1, d2, d3] = ...
            qpsk_tx(data_buf(1),empty_buf(1),clear_buf(1),tx_en_buf(1));
        x_out = complex(i_out,q_out)/2^11;
        x(i1) = x_out;

        %%% Emulate read FIFO AXI interface
        if re_byte_out == 1 && numBytesFromFifo < length(messageWithCRC)
            data_in = messageWithCRC(numBytesFromFifo+1);
            numBytesFromFifo = numBytesFromFifo + 1;
        end
        % processor loaded all bytes into FIFO so begin transmitting
        if numBytesFromFifo == length(messageWithCRC)
            empty_in = 1;
            tx_en_in = 1;
        end
    end

    index = find(abs(x) > sum(SRRC))+24*8; % constant is pad bits
    offset = index(1)+6+length(TB_i)*OS_RATE;
    idx = offset:8:(offset+8*ml*4-1);
    y = x(idx); % four symbos per byte of data
    sc = zeros(1,(msgLength+6)*8);
    sc(1:2:end) = real(y);
    sc(2:2:end) = imag(y);
    sh = sign(sc);
    sb = (sh+1)/2;
    d = zeros(1,ml);
    for i1 = 1:ml
        si = sb(1+(i1-1)*8:i1*8);
        d(i1) = bi2de(si);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Emulate channel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % pad on either side with zeros
    p = complex(zeros(1,100),zeros(1,100));
    xp = [p x p]; % pad

    % Apply frequency offset and receive/over-the-air AWGN
    y = xp.*exp(1i*2*pi*fc*(0:length(xp)-1));
    rC = y/max(abs(y))*.1*2^11; % this controls receive gain
    r = awgn(rC,SNR,0,1);
    r1 = rC;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load Chipscope samples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (~sim)
    fid = fopen('rx.prn');
    M = textscan(fid,'%d %d %d %d','Headerlines',1);
    fclose(fid);
    is = double(M{3});
    qs = double(M{4});
    r = complex(is,qs);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main receiver core
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r_out = zeros(1,length(r));
s_f = zeros(1,length(r));
s_t = zeros(1,length(r));
s_c = zeros(1,length(r));
f_est = zeros(1,length(r));
t_est = zeros(1,length(r));
s_p = zeros(1,length(r));
s_o = zeros(1,length(r));
byte_count = 0;
if (sim)
    bytes = zeros(1,ml);
end
for i1 = 1:length(r)+200
    if i1 > length(r)
        r_in = 0;
    else
        r_in = r(i1);
    end
    i_in = round(real(r_in));
    q_in = round(imag(r_in));
    [r_out(i1), s_f(i1), s_c(i1), s_t(i1), t_est(i1), f_est(i1), ...
        byte, en, s_p(i1), s_o(i1)] = ...
        qpsk_rx(i_in, q_in, floor(muFOC*2^12), floor(muTOC*2^12));

    if en == 1
        byte_count = byte_count + 1;
        bytes(byte_count) = byte;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
subplot(2,2,1)
scatter(real(r),imag(r))
title('Pre FOC Signal');
subplot(2,2,3)
plot(real(r_out));
title('Pre FOC Signal (real part)');
subplot(2,2,2)
scatter(real(s_t),imag(s_t))
title('Post TOC Signal');
subplot(2,2,4)
plot(real(s_p));
title('Message bits');

figure(3)
axis([0 length(s_p) -1.25 1.25]);
plot(s_o);
title('Correlation magnitude');

numRecBytes = bytes(1)+bytes(2)+bytes(3);
msgBytes = bytes((1+4):(numRecBytes+4));

if (sim)
    if sum(msgBytes-message) == 0
        disp('Received message correctly');
    else
        disp('Received message incorrectly');
    end
    native2unicode(bytes);
    native2unicode(msgBytes)
end

if (~sim)
    native2unicode(msgBytes)
end


text:
%% ODE_BATCH_LOCAL uses the BATCH command to run the ODE code locally.
%
%  Discussion:
%
%    The ODE code is a function, so first we must write a script
%    called ODE_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2012
%
%  Author:
%
%    John Burkardt
%
  workers = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ODE_BATCH_LOCAL:\n' );
  fprintf ( 1, '  Run ODE_SCRIPT locally.\n' );
%
%  BATCH defines the job and sets it for execution.
%
  my_job = batch ( 'ode_script', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'ode_fun', 'ode_system' }, ...
    'Configuration', 'local', ...
    'matlabpool', workers );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  LOAD makes the script's workspace available.
%
%  bVals = B values
%  kVals = K values
%  peakVals = peak values
%
  load ( my_job );
%
%  Display the results.
%
  ode_display ( bVals, kVals, peakVals );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ODE_BATCH_LOCAL:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun chebyshevT (@var{n}, @var{x})
%% Numerically evaluate Chebyshev polynomials of the first kind.
%%
%% Evaluates the Chebyshev polynomial of the first kind of degree
%% @var{n} at the point @var{x}, in double precision.  Both inputs
%% can be arrays but their sizes must be either the same or scalar.
%%
%% Example:
%% @example
%% @group
%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
%% chebyshevT (18, 0.9)
%%   @result{} ans = -0.2614
%% @end group
%% @end example
%%
%% Using this function may be preferable to evaluating the Chebyshev
%% polynomial in monomial form because the latter can give poor
%% accuracy due to numerical instability.
%% See the example in @pxref{@@double/chebyshevU}.
%%
%% This function may be slow for large numbers of inputs.
%% This is because it is not a native double-precision implementation
%% but rather the numerical evaluation of the Python @code{mpmath} function
%% @code{chebyshevt}.
%%
%% @seealso{@@sym/chebychevT, @@double/chebyshevU}
%% @end defun

function y = chebyshevT (n, x)
  if (nargin ~= 2)
    print_usage ();
  end

  if (isequal (size (n), size (x)) || isscalar(n))
    y = zeros (size (x));
  elseif (isscalar (x))
    y = zeros (size (n));
  else
    error ('chebyshevT: inputs N and X must have compatible sizes')
  end

  cmd = { 'Ln = _ins[0]'
          'Lx = _ins[1]'
          'if len(Ln) == 1 and len(Lx) != 1:'
          '    Ln = Ln*len(Lx)'
          'if len(Ln) != 1 and len(Lx) == 1:'
          '    Lx = Lx*len(Ln)'
          'c = [complex(mpmath.chebyt(n, x)) for n,x in zip(Ln, Lx)]'
          'return c,' };
  c = pycall_sympy__ (cmd, num2cell (n(:)), num2cell (x(:)));
  for i = 1:numel (c)
    y(i) = c{i};
  end
end

%!error chebyshevT (1)
%!error chebyshevT (1, 2, 3)

%!error <sizes> chebyshevT ([1 2], [1 2 3])
%!error <sizes> chebyshevT ([1 2], [1; 2])

%!test
%! y = sym(11)/10;
%! t = sym(2);
%! x = 1.1;
%! s = 2;
%! A = chebyshevT (s, x);
%! B = double (chebyshevT (t, y));
%! assert (A, B, -2*eps);

%!test
%! % maple
%! A = -0.304681164165948269030369;
%! B = chebyshevT (18.1, 0.9);
%! assert (A, B, -10*eps)

%!test
%! % maple, complex inputs
%! % ChebyshevT(12.1+3.1*I, 0.5+0.2*I);
%! A = 0.637229289490379273451 - 0.475324703778957991318*1i;
%! B = chebyshevT (12.1+3.1*i, 0.5+0.2i);
%! assert (A, B, -5*eps);

%!test
%! % maple, matrix inputs
%! A = [0.59523064198266880000  0.57727442996887552000];
%! B = chebyshevT ([16 17], [0.9 0.7]);
%! assert (A, B, -10*eps);

%!test
%! % x matrix, s scalar
%! y = [1 2 sym(pi); exp(sym(1)) 5 6];
%! t = sym(2);
%! x = double (y);
%! s = 2;
%! A = chebyshevT (s, x);
%! B = double (chebyshevT (t, y));
%! assert (A, B, -eps);

%!test
%! % s matrix, x scalar
%! t = [1 2 sym(pi); exp(sym(1)) 5 6];
%! y = sym(2);
%! s = double (t);
%! x = 2;
%! A = chebyshevT (s, x);
%! B = double (chebyshevT (t, y));
%! assert (A, B, -eps);

%!xtest
%! % https://github.com/fredrik-johansson/mpmath/issues/469
%! assert (chebyshevT (4, inf), inf)
%! assert (chebyshevT (4, -inf), inf)
%! assert (chebyshevT (3, inf), inf)
%! assert (chebyshevT (3, -inf), -inf)

repo: octsympy


text:
function binnedmatrix = toy_simulation(Network_opts,Assembly_opts)

% creating matrix
binnedmatrix=poissrnd(Network_opts.meanspikebin,[Network_opts.nneurons Network_opts.nbins]);

allassemblyneurons = unique([Assembly_opts.assembly_neurons{:}]);
nonassemblyneurons = setdiff(1:Network_opts.nneurons,allassemblyneurons);
nassemblies = length(Assembly_opts.assembly_neurons);

% add activations of all assemblies in random bins
for assemblyindex=1:nassemblies

    assemblysize = length(Assembly_opts.assembly_neurons{assemblyindex});

    % drawing activations
    activation_bins = randi([1,Network_opts.nbins],1,Assembly_opts.number_of_activations);

    % introducting activations
    binnedmatrix(Assembly_opts.assembly_neurons{assemblyindex},activation_bins) = ...
        poissrnd(Assembly_opts.meanspikerate_activations,...
        assemblysize,Assembly_opts.number_of_activations);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef S2KernelHandle < S2Kernel
  % defines a kernel function as a function of the rotational angle

  properties
    fun = @(x) 1;
  end

  methods

    function psi = S2KernelHandle(fun,varargin)

      % extract parameter and halfwidth
      if nargin == 0, return;end

      psi.fun = fun;
      psi.A = calcFourier(psi,getMTEXpref('maxS2Bandwidth'),varargin{:});

    end

    function value = eval(psi,t)
      value = psi.fun(t);
    end

  end
end

author: mtex-toolbox

repo: mtex

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


text:
function[Q1,Q2,Q3,Q4,Q]=divgeom(varargin)
%DIVGEOM  Geometric decomposition of eddy vorticity flux divergence.
%
%   [F1,F2,F3,F4]=DIVGEOM(DX,DY,K,L,THETA) returns the geometric decomposition
%   of eddy vorticity flux divergence associated with variance ellipses
%   having kinetic energy K, anisotropy L, and orientation THETA.
%
%   K, L, and THETA are matrices of the same size.  These are defined on an
%   X-Y grid with x oriented in *columns* and y oriented in *rows*.  DX and
%   DY are the sampling intervals in the X and Y directions, respectively.
%
%   Note that K and L are related to ellipse parameters KAPPA and LAMBDA
%   used elsewhere in JLAB by K=KAPPA^2 and L=LAMBDA*KAPPA^2.
%
%   F1, F2, F3, and F4 are four different contributions to the eddy
%   vorticity flux divergence, as follows:
%
%       F1     Quadratic variations in the linear energy L
%       F2     Product of linear variations in THETA and L
%       F3     Quadratic variations in the orientation THETA
%       F4     Product of linear variations in orientation THETA
%
%   For details, see Waterman and Lilly (2015), Geometric decomposition of
%   eddy-mean flow feedbacks in barotropic systems, J. Phys. Oceanogr.
%
%   [F1,F2,F3,F4,F]=DIVGEOM(...) also returns the total eddy flux
%   divergence F, calculated directly, with F1+F2+F3+F4 = F apart from
%   numerical error.
%
%   By default, DIVGEOM calculates derivates with repeated applications of
%   a first central difference.  DIVGEOM(...,'arakawa') alternately uses a
%   modified first central difference appropriate for models that employ an
%   Awakawa advection scheme.
%   __________________________________________________________________
%
%   Usage: [f1,f2,f3,f4]=divgeom(dx,dy,K,L,theta);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--f')
    %   'divgeom --f' generates a sample figure. XX not currently working
    %divgeom_figure,return
end

%Divgeom does run tests, but these are hidden because it involves mat-files
%not distributed as a part of JLAB
%if strcmpi(varargin{1},XXX)
%    divgeom_test,return
%end

%   Do I do this?
%   Filtering
%
%   DIVGEOM can optionally filter the second-order derivative terms, F1,
%   F3, and F to reduce small-scale noise.
%
%   [F1,F2,F3,F4,F]=DIVGEOM(...,N) smooths F1, F3, and F with an N point
%   boxcar filter in both the X and Y directions.  N should be odd.
%   __________________________________________________________________
%

diffstr='cartesian';
str='second';

for i=1:2
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'car')||strcmpi(varargin{end}(1:3),'ara')
            diffstr=varargin{end};
        elseif strcmpi(varargin{end}(1:3),'dir')||strcmpi(varargin{end}(1:3),'sec')
            str=varargin{end};
        end
        varargin=varargin(1:end-1);
    end
end

dx=varargin{1};
dy=varargin{2};
K=varargin{3};
L=varargin{4};
theta=varargin{5};

if length(varargin)==6
    N=varargin{6};
else
    N=0;
end

%[x,y] = meshgrid([1:size(K,2)]*dx,[1:size(K,1)]*dx);

%/************************************************
%Compute gradients
Lx=divgeom_vdiff(dx,L,2,diffstr);
Ly=divgeom_vdiff(dy,L,1,diffstr);

thetax=divgeom_vdiff(dx,frac(1,2)*unwrap(2*theta,[],2),2,diffstr);
thetay=divgeom_vdiff(dy,frac(1,2)*unwrap(2*theta,[],1),1,diffstr);

cos2x=divgeom_vdiff(dx,cos(2*theta),2,diffstr);
cos2y=divgeom_vdiff(dy,cos(2*theta),1,diffstr);

sin2x=divgeom_vdiff(dx,sin(2*theta),2,diffstr);
sin2y=divgeom_vdiff(dy,sin(2*theta),1,diffstr);

M=L.*cos(2*theta);
N=L.*sin(2*theta);

Nx=divgeom_vdiff(dx,N,2,diffstr);
Ny=divgeom_vdiff(dy,N,1,diffstr);

QM=divgeom_mixederiv(dx,dy,divgeom_vdiff(dx,M,2,diffstr),divgeom_vdiff(dy,M,1,diffstr),diffstr);
QN=divgeom_vdiff(dx,Nx,2,diffstr)-divgeom_vdiff(dy,Ny,1,diffstr);

Q=QN-QM;
%\************************************************

%
% cx=L;
% cx=vdiff(dx,cx,2)+sqrt(-1)*vdiff(dx,cx,1);
% cx=vdiff(dx,cx,2)+sqrt(-1)*vdiff(dx,cx,1);
% Q1=-imag(rot(-2*theta).*cx);

Q1a=-cos(2*theta).*divgeom_mixederiv(dx,dy,Lx,Ly,diffstr);
Q1b= sin(2*theta).*(divgeom_vdiff(dx,Lx,2,diffstr)-divgeom_vdiff(dy,Ly,1,diffstr));
Q1=Q1a+Q1b;

%vsize(L,theta,thetax,thetay)
Q3a= 2*L.*cos(2*theta).*(divgeom_vdiff(dx,thetax,2,diffstr)-divgeom_vdiff(dy,thetay,1,diffstr));
Q3b= 2*L.*sin(2*theta).*divgeom_mixederiv(dx,dy,thetax,thetay,diffstr);
Q3=Q3a+Q3b;

if findstr(str,'dir')
    Q2a= 4*cos(2*theta).*(thetax.*Lx-thetay.*Ly);
    Q2b= 4*sin(2*theta).*(thetax.*Ly+thetay.*Lx);
    Q2=Q2a+Q2b;

    Q4a= 4*L.*cos(2*theta).*(2*thetax.*thetay);
    Q4b= -4*L.*sin(2*theta).*(thetax.^2-thetay.^2);
    Q4=Q4a+Q4b;
else
    dkdsin2=L.*(divgeom_vdiff(dx,sin2x,2,diffstr)-divgeom_vdiff(dy,sin2y,1,diffstr));
    dldcos2=L.*divgeom_mixederiv(dx,dy,divgeom_vdiff(dx,cos(2*theta),2,diffstr),...
        divgeom_vdiff(dy,cos(2*theta),1,diffstr),diffstr);

    Q2a=QN-Q1b-dkdsin2;
    Q2b=-(QM+Q1a-dldcos2);
    Q2=Q2a+Q2b;

    Q4a=-dldcos2-Q3b;
    Q4b=dkdsin2-Q3a;
    Q4=Q4a+Q4b;
end

function[df]=divgeom_vdiff(dx,f,dim,schemestr);
%This is basically just to implement what I think is the way the Arakawa
%advection scheme takes derivatives.  Taken from PSI2FIELDS

f0ca=f(:,1,:,:);
f0cb=f(:,end,:,:);
f0ra=f(1,:,:,:);
f0rb=f(end,:,:,:);

if strcmpi(schemestr(1:3),'ara')
    if dim==1
        dim2=2;
    elseif dim==2
        dim2=1;
    end
    f=frac(1,4)*(2*f + vshift(f,1,dim2) + vshift(f,-1,dim2));
    if dim==1
        %size(f),dim,size(frac(1,3)*(2*f + vshift(f,1,dim2)))
        %f(:,1,:,:)=0;
        %f(:,end,:,:)=0;
        f(:,1,:,:)=frac(1,3)*(2*f0ca + vshift(f0ca,1,dim2));
        f(:,end,:,:)=frac(1,3)*(2*f0cb + vshift(f0cb,-1,dim2));
        %f(:,end,:,:)=f0(:,end,:,:);
    elseif dim==2
        %f(1,:,:,:)=0;
        %f(end,:,:,:)=0;
        f(1,:,:,:)=frac(1,3)*(2*f0ra + vshift(f0ra,1,dim2));
        f(end,:,:,:)=frac(1,3)*(2*f0rb + vshift(f0rb,-1,dim2));
        %f(end,:,:,:)=f0(end,:,:,:);
    end
end
% dim
% size(f)
% edgestr
df=vdiff(dx,f,dim);

function[gxy]=divgeom_mixederiv(dx,dy,fx,fy,diffstr)
%Correction for two possible forms of mixed derivative

gxy= divgeom_vdiff(dy,fx,1,diffstr);
gyx= divgeom_vdiff(dx,fy,2,diffstr);
gxy=gxy+gyx;

%bool=abs(gyx)<abs(gxy);
%gxy(bool)=gyx(bool);
%gxy=2*gxy;

function[]=divgeom_test
load highresjet
use highresjet

[q1,q2,q3,q4,q]=divgeom(dx,K,L,theta);
reporttest('DIVGEOM Q1--Q4 add to total Q',aresame(q1+q2+q3+q4,q,1e-16))

[q1d,q2d,q3d,q4d,qd]=divgeom(dx,K,L,theta,'direct');

rat=abs(frac(q2-q2d,q2+q2d));
rat=sort(rat(:));
index=find(log10(rat)>-1,1,'first');

%No more than 15% of the data has the error ratio > 1/10
reporttest('DIVGEOM Q2 direct and implicit forms agree',index/length(rat)>0.85);

rat=abs(frac(q4-q4d,q4+q4d));
rat=sort(rat(:));
index=find(log10(rat)>-1,1,'first');

%No more than 15% of the data has the error ratio > 1/10
reporttest('DIVGEOM Q4 direct and implicit forms agree',index/length(rat)>0.85);

function[]=divgeom_figure
load jetellipses_highres
use jetellipses

[q1,q2,q3,q4,q]=divgeom(x(2)-x(1),x(2)-x(1),kappabar,lambdabar,thetabar,5);

figure
subplot(2,2,1),jpcolor(x,y,q1),axis equal,axis tight
subplot(2,2,2),jpcolor(x,y,q2),axis equal,axis tight
subplot(2,2,3),jpcolor(x,y,q3),axis equal,axis tight
subplot(2,2,4),jpcolor(x,y,q4),axis equal,axis tight
for i=1:4
    subplot(2,2,i),%caxis([-.2 .2]/16)
    hold on,%contour(x,y,lambdabar,[0.1 0.1]);
    %linestyle k
end
packfig(2,2)

%[q1b,q2b,q3b,q4b,qb]=divgeom(x(2)-x(1),kappabar,lambdabar,thetabar,5,'arakawa');

figure
subplot(2,1,1),jpcolor(x,y,q1+q2+q3+q4),axis equal,axis tight
subplot(2,1,2),jpcolor(x,y,q),axis equal,axis tight
for i=1:2
    subplot(2,1,i),caxis([-.2 .2]/16)
    hold on,contour(x,y,lambdabar,[0.1 0.1]);
    linestyle k
end
packfig(2,1)

function[]=divgeom_figure2

filename='/Users/lilly/Data/early/QGBetaPlaneTurbulenceFloats_experiment_04.nc';
ii=425:675;
[u,v]=vzeros(length(ii),1024,1001);
for k=1:1001
    k
    [utemp,vtemp]=FieldsFromTurbulenceFile(filename,k, 'u','v');
    u(:,:,k)=utemp(ii,:);
    v(:,:,k)=vtemp(ii,:);
end

u2=u;
v2=v;

for k=1:251
    k
    u2(k,:,:)=anatrans(u(k,:,:),3,'periodic');
    v2(k,:,:)=anatrans(v(k,:,:),3,'periodic');
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% Function that optimizes Pose given an initial solution of the rotation only
%
% IMPORTANT: THIS FUNCTION ONLY DEALS WITH THE ORTHOGRAPHIC CASE RIGHT NOW
%
% This function optimizes the pose estimation of one camera with respect to
% another given only the 2 corresponding views.  The data is supposed to be
% centered and only the rotation has to be estimated. I it is here
% performed using simple gradient descent.
%
% The following code explains how the code to compute the
% error/gradient/hessian was created using matlab symbolic toolbox.
%
% Vincent's Structure From Motion Toolbox      Version 2.0\n
% Copyright (C) 2009 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

%%%%%  Define the rotation matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms a b c d real

R = [ a^2+b^2-c^2-d^2 2*b*c-2*a*d 2*a*c+2*b*d; ...
 2*a*d+2*b*c a^2-b^2+c^2-d^2 2*c*d-2*a*b; ...
 2*b*d-2*a*c 2*a*b+2*c*d a^2-b^2-c^2+d^2 ]/(a^2+b^2+c^2+d^2);
R = R(1:2,:);

%%%%%  Compute its derivatives with respect to the quaternions %%%%%%%%%%%%%%%%%%%%%%
l = R;
l = [ diff(l,a) diff(l,b) diff(l,c) diff(l,d) ];
l = simple( l );

%%% Now, use the following:
%%% You need an expression of the form : '[ a[1][1] = x, a[1][2] =y^2 ]'
%%% not the overall brackets and the indexing starting from 1 (as it is
%%% maple. The quotes: simply coz it has to be a string
%%% so you create an array of assigments basically
maple restart;
com = 'res := [ ';
k = 1;
for i = 1 : 12
  for j = 1 : 2
    com = [ com 'dR[' num2str(k) ']=' char(l(j,i)) ', '];
	k = k + 1;
  end
end
com = [ com(1:end-2) '];' ];

%%%%%  Below is the code to generate C-code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Mupad try, but it sucks ...
%  reset(symengine);
%  out = evalin(symengine, [ 'opts := generate::optimize([R=expr(' char(Rp') ') ]);' ] );
%
%  out = evalin(symengine, [ 'opts := generate::optimize([R=expr(' char(Rp') '), dR=expr(' char(l(:,3:10)) '), ddR=expr(' char(l(:,11:end)) ')]);' ] );
%  evalin(symengine, [ ':=rhs(opts[-1]);' ] );
%
%
%
%  out = [ char(evalin(symengine, [ 'generate::C(opts[1..-2]);' ] )) ...
%    char( evalin(symengine, [ 'generate::C(rhs(opts[-1]))' ] )) ];
%  out = strrep( out, '"', '' );
%  out = strrep( out, [ var '[0][' ], [ var '[' ] );

%%% now, save your crazy array into the res variable
maple( com );

%  %%% cost of the assigments with no optimization
%  maple('codegen[cost](res)')
%  % maple( 'opt1 := [codegen[optimize](l)]' );
%  %%% cost of the assigment with the normal optimization
%  in = [ 'l:=' char( l ) ]; maple( in ); maple('codegen[cost](codegen[optimize](l))')

%%% code for better optimization and corresponding cost
maple( 'opt2 := [codegen[optimize](res,tryhard)]' );
maple('codegen[cost](opt2)')

%%% Convert to C code
out = maple( 'codegen[C](opt2)' );

%%% clean the C code and convert it to Matlab code (replace [] by ()
out = strrep( out, ';   ', ';\n' );
out = strrep( out, '~', '' );

tout = ''
existing = cell(1,0);
tmp = regexp( out,'(t\d*)', 'tokens' )
for i = 1 : length(tmp)
  ii = tmp(i);
  ii = ii{1}{1};
  doExist = false;
  for j = 1 : size(existing,2)
    if strcmp(existing{j},ii)
	  doExist = true;
	  break;
	end
  end
  if ~doExist
    existing{1,end+1} = ii;
    tout = [ tout, ', ', ii ];
  end
end
out = [ 'double ', tout(2:end), ';\n', out ];
fprintf(out)

author: vrabaud


text:
% test_cgal_check_self_intersect.m

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% nice mesh without self intersections

xyz = [
    -1     1    0
     0     1    0
    -0.5   0.5  0
    -1     0    0
     0     0    0
     1     0    0
     0    -1    0
    ];

x = xyz(:, 1);
y = xyz(:, 2);
z = xyz(:, 3);

tri = [
    1 2 3
    1 4 3
    3 2 5
    3 4 5
    2 5 6
    4 5 7
    5 7 6
    ];

% interior vertices that are nearest neighbors
I = 3;
J = 5;

% plot mesh, with each face in a different colour
hold off
gplot3d(dmatrix_mesh(tri), xyz)
hold on
plot3(x, y, z, 'o')
axis equal xy
view(2)

% check triangle intersections
cross = cgal_check_self_intersect(tri, xyz);
cross'
% expected: all zeros

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% move vertex I down so that triangle 1 overlaps with triangles 4 and 6

y(I) = -1/4.5;
xyz = [x y z];

% plot mesh, with each face in a different colour
hold off
gplot3d(dmatrix_mesh(tri), xyz)
hold on
plot3(x, y, z, 'o')
axis equal xy
view(2)

% check triangle intersections
cross = cgal_check_self_intersect(tri, xyz);
cross'
% expected: 2     2     2     4     0     4     0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% spherical mesh

% uniformly distribute particles across the surface of the unit sphere
% using Anton Semechko's implementation of Reisz s-energy minimisation
% [xyz, tri] = ParticleSampleSphere('N', 20);

% use always the same points for reproducibility
xyz = [
    0.1865    0.1323   -0.9735
    0.8520    0.0405   -0.5219
    0.8562    0.4596    0.2358
    0.5982   -0.1356    0.7897
    0.2115    0.5448    0.8115
    0.4822    0.7204   -0.4986
   -0.9091    0.3984   -0.1217
    0.1588    0.9741    0.1609
   -0.5982    0.1356   -0.7897
   -0.5127    0.6897    0.5113
   -0.3189    0.8227   -0.4706
   -0.8520   -0.0405    0.5219
   -0.1866   -0.1323    0.9735
    0.5127   -0.6897   -0.5113
   -0.1588   -0.9741   -0.1609
   -0.2115   -0.5448   -0.8115
   -0.4822   -0.7204    0.4986
   -0.8562   -0.4596   -0.2358
    0.3189   -0.8227    0.4706
    0.9091   -0.3984    0.1217
];

tri = [
     7    10    11
     8    11    10
     1    11     6
     1     6     2
     8     6    11
    20     2     3
    20     3     4
     2     6     3
     8     3     6
     8    10     5
     8     5     3
     4     3     5
    13     5    10
    13     4     5
    20    14     2
     1     2    14
    20     4    19
    20    19    14
    15    14    19
    13    19     4
     7    11     9
     1     9    11
     7    12    10
    13    10    12
    15    16    14
     1    14    16
     1    16     9
    15    18    16
     9    16    18
     7     9    18
     7    18    12
    15    19    17
    15    17    18
    12    18    17
    13    17    19
    13    12    17
     ];

% check triangle intersections
cross = cgal_check_self_intersect(tri, xyz);
cross'

% plot mesh, colouring the facets according to whether they are overlapping
% or not
subplot(1, 2, 1)
hold off
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3), cross);
axis equal

%% create some triangle overlap

% select one vertex
I = 17;

% select a nearest neighbour
J = 19;

% plot them on the mesh
hold on
plot3(xyz(I, 1), xyz(I, 2), xyz(I, 3), 'ro')
plot3(xyz(J, 1), xyz(J, 2), xyz(J, 3), 'r*')

% move the nearest neighbor in a way such that I is still within the convex
% hull, but the nearest neighbour causes overlap
lat(J) = 13/180*pi;
lon(J) = -120/180*pi;

[xyz(J, 1), xyz(J, 2), xyz(J, 3)] = sph2cart(lon(J), lat(J), 1);

% check triangle intersections
cross = cgal_check_self_intersect(tri, xyz);
cross'

% triangles that contain the vertex we moved
find(sum(tri == J, 2)>0)
find(cross)

% plot mesh, colouring the facets according to whether they are overlapping
% or not
subplot(1, 2, 2)
hold off
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3), double(cross~=0));
axis equal
hold on
plot3(xyz(I, 1), xyz(I, 2), xyz(I, 3), 'ro')
plot3(xyz(J, 1), xyz(J, 2), xyz(J, 3), 'r*')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% two 3D triangles orthogonal and intersecting, intersection=segment

xyz = [
    -1 -1 0   % point for horizontal triangle
    1 -1 0    % point for horizontal triangle
    -1 1 0    % point for horizontal triangle
    -.5 0 1   % point for vertical triangle
    .5 0 1    % point for vertical triangle
    -.5 0 -1  % point for vertical triangle
    ];

tri = [
    1 2 3     % horizontal
    4 5 6     % vertical
    ];

% check triangle intersections
cross = cgal_check_self_intersect(tri, xyz);
cross'

% plot mesh, colouring the facets according to whether they are overlapping
% or not
subplot(1, 1, 1)
hold off
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3), cross);
axis equal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% two 3D triangles orthogonal and intersecting, intersection=point

xyz = [
    -1 -1 0   % point for horizontal triangle
    1 -1 0    % point for horizontal triangle
    -1 1 0    % point for horizontal triangle
    -.5 0 2   % point for vertical triangle
    .5 0 2    % point for vertical triangle
    -.5 0 0  % point for vertical triangle
    ];

tri = [
    1 2 3     % horizontal
    4 5 6     % vertical
    ];

% check triangle intersections
cross = cgal_check_self_intersect(tri, xyz);
cross'

% plot mesh, colouring the facets according to whether they are overlapping
% or not
subplot(1, 1, 1)
hold off
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3), cross);
axis equal

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function dlY = spatialDropout(dlX,dropoutFactor,fmt,doTraining)

if doTraining
    maskSize = size(dlX);
    maskSize(fmt=='S') = 1;

    dropoutScaleFactor = single(1 - dropoutFactor);
    dropoutMask = (rand(maskSize,'like',dlX) > dropoutFactor) / dropoutScaleFactor;

    dlY = dlX .* dropoutMask;
else
    dlY = dlX;
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function result = MaskRecon(Data, v)

for i = 1 : size(Data, 1)
    dis(i) = norm(v-Data(i, :));
end
[mindis, ind] = min(dis);
if ind < 11
    result = ind-1;
else
    result = ind-11;
end

author: decouples


text:
%kktan 'Compute Tangent or Arc Tangent Function'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ktan.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% Toggle: tan 'Tangent', default: 0: 'Calculate tangent of input'
% Toggle: atan 'Arc Tangent', default: 0: 'Calculate arc tangent of input'
%
% Example: o = kktan(i, {'i','';'o','';'tan',0;'atan',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ktan - Compute Tangent or Arc Tangent Function
%
%  DESCRIPTION
% The "Tan & ArcTan" operator returns the tangent, or arc tangent, value
% of each data point in the input data object, \fBInput".  The forward or
% inverse trigonometric operation is specified by selecting the \fBTangent" or
% \fBArc Tangent" option on the user interface.
%
% Executing "Tan & ArcTan" runs the program \fIkarith1\fP
% with the -tan or -atan flag.
%
%  "\fBData Type"" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "\fBMap Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_1input
%
%  "\fBValidity Mask"" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "\fBExplicit Location and Time Data"" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "\fBFailure Modes""
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith1
%
%  RESTRICTIONS
% Operations on complex data are not supported at this time.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kktan(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kktan(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'tan', 0;'atan', 0};
maxval={0,0,0,0};
minval={0,0,0,0};
istoggle=[0,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith1"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% FUNCTION init_opts
% initialization options for multi-task learning library
%
% If one of the ncessary opts are empty then it will be set to default
% values as specified in this file.
%
% Table of Options.  * * indicates default value.
% FIELD                DESCRIPTION
%% Optimization options
%
%  .max_iter               Maximum iteration step number
%                           *1000*
%  .tol                    Tolerance
%                           *10e-3*
%  .tFlag                  Termination condition
%                           0 => change of absolute function value:
%                             abs( funcVal(i)- funcVal(i-1) ) <= .tol
%                         * 1 => change of relative function value:
%                             abs( funcVal(i)- funcVal(i-1) )
%                              <= .tol * funcVal(i-1)
%                           2 => absolute function value:
%                             funcVal(end)<= .tol
%                           3 => Run the code for .maxIter iterations
%% Starting Point
%
% .W0               Starting point of W.
%                   Initialized according to .init.
%
% .C0               Starting point for the intercept C (for Logistic Loss)
%                   Initialized according to .init.
%
% .init             .init specifies how to initialize W, C.
%                         0 => .W0, C0 are set by a guess value infered from data
%                         1 => .W0 and .C0 are defined
%                       * 2 => .W0= zeros(.), .C0=0 *
%
%% Parallel Computing
%
% .pFlag            Enable Map-Reduce (needs Parallel Toolbox support).
%                        *false*
%
% .pSeg_num         set the number of total parallel segmentations. if the
%                   number is non-positive, then current matlab pool number
%                   will be used. [supported in the incoming version]
%                        * pool size *
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 21, 2012.
%

function opts = init_opts (opts)

%% Default values
DEFAULT_MAX_ITERATION = 1000;
DEFAULT_TOLERANCE     = 1e-4;
MINIMUM_TOLERANCE     = eps * 100;
DEFAULT_TERMINATION_COND = 1;
DEFAULT_INIT = 2;
DEFAULT_PARALLEL_SWITCH = false;

%% Starting Point
if isfield(opts,'init')
    if (opts.init~=0) && (opts.init~=1) && (opts.init~=2)
        opts.init=DEFAULT_INIT; % if .init is not 0, 1, or 2, then use the default 0
    end

    if (~isfield(opts,'W0')) && (opts.init==1)
        opts.init=DEFAULT_INIT; % if .W0 is not defined and .init=1, set .init=0
    end
else
    opts.init = DEFAULT_INIT; % if .init is not specified, use "0"
end

%% Tolerance
if isfield(opts, 'tol')
    % detect if the tolerance is smaller than minimum
    % tolerance allowed.
    if (opts.tol <MINIMUM_TOLERANCE)
        opts.tol = MINIMUM_TOLERANCE;
    end
else
    opts.tol = DEFAULT_TOLERANCE;
end

%% Maximum iteration steps
if isfield(opts, 'maxIter')
    if (opts.maxIter<1)
        opts.maxIter = DEFAULT_MAX_ITERATION;
    end
else
    opts.maxIter = DEFAULT_MAX_ITERATION;
end

%% Termination condition
if isfield(opts,'tFlag')
    if opts.tFlag<0
        opts.tFlag=0;
    elseif opts.tFlag>3
        opts.tFlag=3;
    else
        opts.tFlag=floor(opts.tFlag);
    end
else
    opts.tFlag = DEFAULT_TERMINATION_COND;
end

%% Parallel Options
if isfield(opts, 'pFlag')
    if opts.pFlag == true && ~exist('matlabpool', 'file')
        opts.pFlag = false;
        warning('MALSAR:PARALLEL','Parallel Toolbox is not detected, MALSAR is forced to turn off pFlag.');
    elseif opts.pFlag ~= true && opts.pFlag ~= false
        % validate the pFlag.
        opts.pFlag = DEFAULT_PARALLEL_SWITCH;
    end
else
    % if not set the pFlag to default.
    opts.pFlag = DEFAULT_PARALLEL_SWITCH;
end

if opts.pFlag
    % if the pFlag is checked,
    % check segmentation number.
    if isfield(opts, 'pSeg_num')
        if opts.pSeg_num < 0
            opts.pSeg_num = matlabpool('size');
        else
            opts.pSeg_num = ceil(opts.pSeg_num);
        end
    else
        opts.pSeg_num = matlabpool('size');
    end
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
%% Unimodal ODF Shapes
%
%%
% Also take a look at the page <SO3Kernels.html SO3Kernels>.
%
% In order to control the shape of unimodal ODF
% The classes @SO3Kernel and @S2Kernel are needed in MTEX to define the
% specific form of unimodal and fibre symmetric ODFs. It has to be passed
% as an argument when calling the methods <uniformODF.html uniformODF> or
% <fibreODF.html fibreODF>.
%
% A kernel is defined by specifying its name and its free parameter.
% Alternatively one can also specify the halfwidth of the kernel. Below you
% find a list of some important SO3Kernel functions supported by MTEX.

psi{1} = SO3AbelPoissonKernel(0.79);
psi{2} = SO3DeLaValleePoussinKernel(13);
psi{3} = SO3BumpKernel(35*degree);
psi{4} = SO3DirichletKernel(3);
psi{5} = SO3vonMisesFisherKernel(7.5);
psi{6} = SO3GaussWeierstrassKernel(0.07);
psi{7} = fibreVonMisesFisherKernel(7.2);
psi{8} = SO3SquareSingularityKernel(0.72);

%%
% Lets visualize these kernel functions as one dimensional sections through
% the orientation space

% the kernel on SO(3)
close;
figure('position',[100,100,1000,450])
hold all
for i = 1:numel(psi)
  plot(psi{i},'DisplayName',class(psi{i}));
end
hold off
legend(gca,'show','Location','eastoutside')

%%
% one dimensional sections through the corresponding PDF

close; figure('position',[100,100,1000,450])
hold all
for i = 1:numel(psi)
  plot(psi{i}.radon,'symmetric','DisplayName',class(psi{i}),'linewidth',2);
end
hold off

ylim([-5,20])

legend(gca,'show','Location','eastoutside')
%%
% the Fourier coefficients of the kernels

close; figure('position',[100,100,500,450])
hold all
for i = 1:numel(psi)
  plotSpektra(psi{i},'bandwidth',32,'linewidth',2,'DisplayName',class(psi{i}));
end
hold off
legend(gca,'show')


text:
function M = esvm_perform_calibration(grid, val_set, models, params)
% 1. Perform LABOO calibration procedure and 2. Learn a combination
% matrix M which multiplexes the detection results (by compiling
% co-occurrence statistics on true positives)
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

%% Perform calibration
betas = esvm_perform_platt_calibration(grid, val_set, models, ...
                                       params);

%% Estimate the co-occurrence matrix M
if ~(isfield(params,'SKIP_M') && params.SKIP_M==1)
  M = esvm_estimate_M(grid, models, params);
end

%concatenate results
M.betas = betas;

author: quantombone

repo: exemplarsvm


text:
function testGaborFilter
%% by dzhg: zhgdai@126.com
%% This code is to search the direction of the vessel for each point of
%% the Image. I apply Gabor Filter to this work and the corresponding
%% output Gabor filter will be choosed to be the direction of the vessel.
%% if you have any question about this code, please contact with me by
%% zhgdai@126.com without hesitate.
%% This edition is to apply adaptative size of gabor filter to detect the
%% direction of the vessel.

I = imread('00005.bmp');

if isgray(I) ==0
    I = rgb2gray(I);
end

[ImageWidth,ImageHeight] = size(I);
sizeinterval = 10;
Begin_s = 10;
End_s = 40;
angleinterval = 18;

%% papers parameters
s = Begin_s:sizeinterval:End_s;
theta = 0:pi/angleinterval:(angleinterval-1)/angleinterval*pi;

f = 1./s;
Sx = s/pi;
Sy = s/pi;
LengthTheta = length(theta);
LengthS = length(s);

GaborFilter = zeros(ImageWidth,ImageHeight,(LengthTheta-1)*LengthS);
Index = 1;
for j = 1:LengthS
    for i = 1:LengthTheta-1
    [G,gabout] = gabordzhg(I,Sx(j),Sy(j),f(j),theta(i));
    GaborModel{Index} = G;
    GaborFilter(:,:,Index) = gabout;
    Index = Index + 1;
    end
end

imshow(I);
hold on;
[a,b] = ginput;
% a = 275;
% b = 245;

GaborFilterValue = GaborFilter(floor(b),floor(a),1:end);
savevar = [];
for i = 1:LengthS
    tempvar = var(GaborFilterValue((i-1)*(angleinterval-1)+1:i*(angleinterval-1)));
    savevar = [savevar;tempvar];
end
[maxvar,maxindex] = max(savevar);
[minvalue,minindex] = min(GaborFilterValue((maxindex-1)*(angleinterval-1)+1:maxindex*(angleinterval-1)));
imshow(GaborModel{(maxindex-1)*(angleinterval-1)+minindex},[]);

author: Sable

repo: mcbench-benchmarks