text:
function[mat,xmid,ymid]=twodhist(varargin)
%TWODHIST  Two-dimensional histogram.
%   __________________________________________________________________
%
%   *|* twodhist.png --- Histogram of the global drifter dataset.
%   Type 'jhelp twodhist' to view this image. *|*
%   __________________________________________________________________
%
%   MAT=TWODHIST(X,Y,XBIN,YBIN) where X and Y are arrays of the same
%   length, creates a two-dimensional histogram MAT with bin edges
%   specified by XBIN and YBIN.
%
%   If XBIN and YBIN are length N and M, respectively, then MAT is of
%   size M-1 x N-1.  XBIN and YBIN must be monotonically increasing.
%
%   [MAT,XMID,YMID]=TWODHIST(...) optionally returns the midpoints XMID
%   and YMID of the bins.
%
%   TWODHIST, TWODSTATS, and TWODMED are three related functions for
%   computing statistics as a function two variables using very fast
%   algorithms that avoid any loops through efficient use of indexing.
%
%   X and Y can also be cell arrays of numerical arrays, in which case
%   all data values are concatented prior to finding the histogram.
%   __________________________________________________________________
%
%   Automatic bin calculation
%
%   TWODHIST can compute appropriate bins internally.
%
%   [MAT,XMID,YMID]=TWODIST(X,Y,N) uses N bins in the X and Y directions,
%   linearly spaced between the minimum and maximum values, and returns the
%   bin midpoints in XMID and YMID.  MAT is N-1 x N-1.
%
%   [MAT,XMID,YMID]=TWODIST(X,Y,[XMIN XMAX],[YMIN YMAX],N) similarly uses N
%   bins, linearly spaced between the designated X and Y values.
%   __________________________________________________________________
%
%   Algorithms
%
%   By default, TWODHIST now works with an internal call to Matlab's
%   HISTCOUNTS2 function, available as of Matlab 2015b.  This is much
%   faster than the previous algorithm.
%
%   If HISTCOUNTS2 is not available, TWODHIST uses loopless algorithm that
%   is in turn much faster than an explicit loop.  TWODHIST(...,'jlab')
%   uses this algorithm, while TWODHIST(...,'slow') uses the explicit loop.
%   These options are mostly used for testing purposes.
%   __________________________________________________________________
%
%   Parallelization
%
%   TWODHIST(...,'parallel') parallelizes the computation using the fast
%   algorithm together with SPMD.  This requires that Matlab's Parallel
%   Computing Toolbox be installed.  While TWODHIST is already very fast,
%   parallelization may be useful for extremely large datasets.
%   __________________________________________________________________
%
%   See also TWODMED, TWODSTATS.
%
%   'twodhist --f' generates the sample figure shown above.
%   'twodhist --t' runs some tests.
%
%   Usage: [mat,xmid,ymid]=twodhist(x,y,N);
%          [mat,xmid,ymid]=twodhist(x,y,[xmin xmax],[ymin ymax],N);
%          mat=twodhist(x,y,xbin,ybin);
%          [mat,xmid,ymid]=twodhist(x,y,xbin,ybin,'parallel');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 J.M. Lilly --- type 'help jlab_license' for details

%   Additional output
%
%   [MAT,XMID,YMID,INDEX]=TWODHIST(X,Y,...) also returns INDEX, an array of
%   the same size as X and Y giving the index into the matrix MAT
%   corresponding to each (X,Y) data point.
%   __________________________________________________________________

if strcmpi(varargin,'--t')
   twodhist_test;return
elseif strcmpi(varargin,'--f')
   type makefigs_twodhist
   makefigs_twodhist;
   return
end

parstr='serial';
str='histcounts2';

xdata=varargin{1};
ydata=varargin{2};

if iscell(xdata)
    [xdata,ydata]=cell2col(xdata,ydata);
end
if ~isreal(xdata)||~isreal(ydata)
    error('X and Y must be real-valued.');
end
if ~aresame(size(xdata),size(ydata))
     error('X and Y should have the same size.')
end
vcolon(xdata,ydata);

bool=isfinite(xdata)&isfinite(ydata);
xdata=xdata(bool);
ydata=ydata(bool);

for i=1:2
    if ischar(varargin{end})
        if strcmpi(varargin{end}(1:3),'ser')||strcmpi(varargin{end}(1:3),'par')
            parstr=varargin{end};
        else
            str=varargin{end};
        end
        varargin=varargin(1:end-1);
    end
end

if strcmpi(str(1:3),'his')&&(exist('histcounts2')~=2)
    str='jlab';
    disp('HISTCOUNTS2 not found; reverting to former JLAB algorithm.')
end

if strcmpi(parstr(1:3),'par')
    if exist('parpool')~=2
        disp('Sorry, parallel algorithm requires the Parallel Computing Toolbox.')
        disp('Defaulting to the serial algorithm.')
        str='serial';
    end
end

if length(varargin)==4
    xbin=varargin{3};
    ybin=varargin{4};
elseif length(varargin)==3
    N=varargin{3};
    xbin=linspace(minmin(xdata),maxmax(xdata),N);
    ybin=linspace(minmin(ydata),maxmax(ydata),N);
elseif length(varargin)==5
    N=varargin{5};
    xbin=linspace(varargin{3}(1),varargin{3}(2),N);
    ybin=linspace(varargin{4}(1),varargin{4}(2),N);
end

xbin=xbin(:);
ybin=ybin(:);
if any(diff(xbin)<0)
  error('XBIN must be monotonically increasing')
end
if any(diff(ybin)<0)
  error('YBIN must be monotonically increasing')
end

%Exclude points which are obviously outside of the domain
bool=xdata<xbin(end)&xdata>xbin(1)&ydata<ybin(end)&ydata>ybin(1);
xdata=xdata(bool);
ydata=ydata(bool);

if ~isempty(xdata)
    if strcmpi(parstr(1:3),'ser')
        mat=twodhist_one(xdata,ydata,xbin,ybin,str);
    elseif strcmpi(parstr(1:3),'par')
        %Parallel Algorithm
        disp('TWODHIST employing parallel algorithm.')
        pool=gcp;
        Nworkers=pool.NumWorkers;
        N=length(xdata);
        M=floor(N/Nworkers);
        spmd
            %Determine the data to send to each worker
            if labindex<Nworkers
                spmdindex=(labindex-1)*M+1:labindex*M;
            else
                %A few leftover time series go to the last worker
                spmdindex=(labindex-1)*M+1:N;
            end
            mati=twodhist_one(xdata(spmdindex),ydata(spmdindex),xbin,ybin,str);
        end
        mat=zeros(size(mati{1}));
        for i=1:length(mati)
            mat=mat+mati{i};
        end
    end
else
    disp('Warning: No valid data in specified region.')
    mat=0*oprod(ybin(1:end-1),xbin(1:end-1));
end

if nargout>1
  xmid=(xbin+vshift(xbin,1,1))./2;
  xmid=xmid(1:end-1);
end
if nargout>2
  ymid=(ybin+vshift(ybin,1,1))./2;
  ymid=ymid(1:end-1);
end

function[mat]=twodhist_one(xdata,ydata,xbin,ybin,str)

if strcmpi(str(1:3),'his')
    mat=histcounts2(xdata,ydata,xbin,ybin)';
elseif strcmpi(str(1:3),'jla')
    mat=twodhist_jlab(xdata,ydata,xbin,ybin);
elseif strcmpi(str(1:3),'slo')
    mat=twodhist_slow(xdata,ydata,xbin,ybin);
end

function[mat]=twodhist_jlab(xdata,ydata,xbin,ybin)

[xnum,xi,xmid]=bindata(xbin,xdata);
[ynum,yi,ymid]=bindata(ybin,ydata);

mat=zeros([length(ybin)-1,length(xbin)-1]);
index=nan*zeros(size(xdata));

nani=(~isnan(xnum)&~isnan(ynum));

if sum(nani(:))>0
    index(nani)=sub2ind([length(ybin)-1,length(xbin)-1],ynum(nani),xnum(nani));
    [indexsorted,sorter]=sort(index(nani));
    [L,ia]=blocklen(indexsorted);
    mat(indexsorted(ia))=L(ia);
end

function[mat]=twodhist_slow(xdata,ydata,xbin,ybin)
mat=zeros(length(ybin),length(xbin));
[xbinb,ybinb]=vshift(xbin,ybin,1,1);
for i=1:length(xbin)
    for j=1:length(ybin)
        mat(j,i)=length(find(xdata>xbin(i)&xdata<=xbinb(i)&...
            ydata>ybin(j)&ydata<=ybinb(j)));
    end
end
mat=mat(1:end-1,:);
mat=mat(:,1:end-1);

function[]=twodhist_test
L=10000;
xdata=3*abs(rand(L,1));
ydata=3*abs(rand(L,1));
xbin=(0:.1:2);
ybin=(0:.2:2);
tic;
mat1=twodhist(xdata,ydata,xbin,ybin);
dt1=toc;
tic
mat2=twodhist(xdata,ydata,xbin,ybin,'slow');
dt2=toc;
bool=aresame(mat1,mat2,1e-10);
reporttest('TWODHIST fast vs. slow algorithm',bool)
%disp(['TWODHIST fast algorithm was ' num2str(dt2./dt1) ' times faster than direct algorithm.'])

if exist('parpool')==2
    mat3=twodhist(xdata,ydata,xbin,ybin,'parallel');
    bool=aresame(mat1,mat3,1e-10);
    reporttest('TWODHIST fast vs. parallel algorithm',bool)
end

xdata=-3*abs(rand(L,1));
ydata=-3*abs(rand(L,1));
xbin=(-2:.1:0);
ybin=(-2:.2:0);
mat1=twodhist(xdata,ydata,xbin,ybin);
mat2=twodhist(xdata,ydata,xbin,ybin,'slow');
bool=aresame(mat1,mat2,1e-10);
reporttest('TWODHIST fast vs. slow algorithm, negative bins',bool)

xdata=randn(L,1);
ydata=randn(L,1);
xbin=(-2:.1:2);
ybin=(-2:.2:2);
mat1=twodhist(xdata,ydata,xbin,ybin);
mat2=twodhist(xdata,ydata,xbin,ybin,'slow');
bool=aresame(mat1,mat2,1e-10);
reporttest('TWODHIST fast vs. slow algorithm, crossing zero',bool)

L=10000;
xdata=3*abs(rand(L,1));
ydata=3*abs(rand(L,1));
xbin=(0:.1:2);
ybin=(0:.2:2);
tic;
mat1=twodhist(xdata,ydata,xbin,ybin);
dt1=toc;
tic
mat2=twodhist(xdata,ydata,xbin,ybin,'slow');
dt2=toc;
bool=aresame(mat1,mat2,1e-10);
reporttest('TWODHIST fast vs. slow algorithm',bool)

if exist('histcounts2')==2
    xdata=randn(L,1);
    ydata=randn(L,1);
    xbin=(-2:.1:2);
    ybin=(-2:.2:2);
    tic;mat1=twodhist(xdata,ydata,xbin,ybin,'histcounts2');toc;
    tic;mat2=twodhist(xdata,ydata,xbin,ybin,'jlab');toc;
    bool=aresame(mat1,mat2,1e-10);
    reporttest('TWODHIST JLAB vs. HISTCOUNTS algorithm, crossing zero',bool)
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
 function out = subsref(ob, args)
%function out = subsref(ob, args)
%	evaluate call of the form ob.ref or ob(ref,:)
%	fix: needs a lot of work to handle ordered subsets!

%
%	enable handling of multiple subscripts, e.g., ob.field()
%
if length(args) > 1
	arg = args(1);
else
	arg = args;
end

%
%	ob.?
%
if arg.type == '.'
	out = getfield(struct(ob), arg.subs);

%
%	ob(?)
%
elseif arg.type == '()'
	subs = arg.subs;

	if ~isempty(ob.index1) | ~isempty(ob.index2)
		error 'subscripts after subscripts perhaps not done'
	end

	%
	%	G(:,:)
	%
	if length(subs) == 2 & subs{1} == ':' & subs{2} == ':'
		out = ob(:,[1:ob.dims(2)]);

	%
	%	Gt(:,ii) - for ordered subsets
	%
	elseif length(subs) == 2 & subs{1} == ':' ...
		& ndims(subs{2}) == 2 & ob.is_transpose
		nb = ob.nb;
		ii = subs{2};
		ia = 1 + round((ii(1,:) - 1) / nb);
		if any(ii(:) ~= col(outer_sum(1:nb,(ia-1)*nb)))
			error 'bad subref'
		end
		ob.ia_start = ia(1) - 1;
		ob.ia_inc = ia(2) - ia(1);
		if any(diff(ia) ~= ob.ia_inc)
			error 'nonuniform angles'
		end
		ob.is_subset = true;
		out = ob;

	%
	%	G(:,j)
	%
	elseif length(subs) == 2 & subs{1} == ':'
		jj = col(subs{2});
		if islogical(jj)
			if length(jj) ~= ob.dims(2)
				error 'bad column logical length'
			end
		elseif any(jj < 1 | jj > ob.dims(2))
			error 'bad column index'
		end

		%
		%	handle G(:,mask(:)) as a special case
		%
		if length(jj) == ob.nx*ob.ny & all(jj == ob.mask(:)) ...
				& ~ob.is_masked & ~ob.is_transpose
			ob.is_masked = true;
			ob.dims(2) = sum(ob.mask(:));
			out = ob;

		else
			out = zeros(ob.dims(1),length(jj));
			for nn=1:length(jj)
				x = zeros(ob.dims(2),1);
				x(jj(nn)) = 1;
				out(:,nn) = ob * x;
			end
			out = sparse(out);
		end

	%
	%	G(i,:)
	%
	elseif length(subs) == 2 & subs{2} == ':'
		ii = col(subs{1});
		if islogical(ii)
			if length(ii) ~= ob.dims(1)
				error 'bad row logical length'
			end
		elseif any(ii < 1 | ii > ob.dims(1))
			error(sprintf('bad row index %s', int2str(ii)))
		end

		out = zeros(ob.dims(2), length(ii));
		for nn=1:length(ii)
			y = zeros(ob.dims(1), 1);
			y(ii(nn)) = 1;
			out(:,nn) = ob' * y;
		end
		out = sparse(out);

	else
		arg.subs
		error todo
	end

else
	error(sprintf('type %s notdone', arg.type))
end

if length(args) > 1
	out = subsref(out, args(2:end));
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%DIAGNOSTIC TOOLS (in the diag-folder):
%
% Convergence diagnostics
%   PSRF     - Potential Scale Reduction Factor
%   CPSRF    - Cumulative Potential Scale Reduction Factor
%   MPSRF    - Multivariate Potential Scale Reduction Factor
%   CMPSRF   - Cumulative Multivariate Potential Scale Reduction Factor
%   IPSRF    - Interval-based Potential Scale Reduction Factor
%   CIPSRF   - Cumulative Interval-based Potential Scale Reduction Factor
%   KSSTAT   - Kolmogorov-Smirnov goodness-of-fit hypothesis test
%   HAIR     - Brooks' hairiness convergence diagnostic
%   CUSUM    - Yu-Mykland convergence diagnostic for MCMC
%   SCORE    - Calculate score-function convergence diagnostic
%   GBINIT   - Initial iterations for Gibbs iteration diagnostic
%   GBITER   - Estimate number of additional Gibbs iterations
%
% Time series analysis
%   ACORR      - Estimate autocorrelation function of time series using xcorr
%   ACORR2     - Estimate autocorrelation function of time series using fft
%   ACORRTIME  - Estimate autocorrelation evolution of time series (simple)
%   GEYER_ICSE - Compute autocorrelation time tau using Geyer's
%                initial convex sequence estimator
%                (requires Optimization toolbox)
%   GEYER_IMSE - Compute autocorrelation time tau using Geyer's
%                initial monotone sequence estimator
%
% Survival model criteria
%   AUCS       - Compute area under curve for survival model
%   AUCT       - Compute area under curve for survival model at given time
%   EXT_AUC    - Compute Extended AUC proposed by Chambless et al (2011)
%   HCS        - Compute Harrell's C for survival model at given time
%   HCT        - Compute Harrel's C for survival model at several time points
%   IDIS       - Integrated Discrimination Improvement between two models
%   RSQR       - R^2 statistic given probabilities at time point T
%
% Kernel density estimation etc.:
%   KERNEL1  - 1D Kernel density estimation of data
%   KERNELS  - Kernel density estimation of independent components of data
%   KERNELP  - 1D Kernel density estimation, with automatic kernel width
%   NDHIST   - Normalized histogram of N-dimensional data
%   HPDI     - Estimates the Bayesian HPD intervals
%
% Misc:
%   CUSTATS   - Calculate cumulative statistics of data
%   GRADCHEK  - Checks a user-defined gradient function using finite
%               differences.
%   DERIVATIVECHECK - Compare user-supplied derivatives to
%                     finite-differencing derivatives.
%

author: gpstuff-dev

repo: gpstuff

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


text:
function value = p01_autonomous ( )

%*****************************************************************************80
%
%% P01_AUTONOMOUS reports whether problem P01 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Tests the solvers on a simple unconstrained quadratic function

%{
    Solve:

    minimize_x  c'x + x'Dx/2

    as an example of using TFOCS without the "SCD" interface

%}

randn( 'state', sum('quadratic test') );
N       = 100;
c       = randn(N,1);
D       = randn(N,N);
D       = D * D' + .5*eye(N);
x_star  = - D \ c;          % the true answer
x0      = zeros(N,1);

% Here's what you could do...
% f       = @(x) c'*x + x'*D*x/2;
% grad_f  = @(x) c + D*x;
% smoothF = @(x) wrapper_objective( f, grad_f, x );

% Here's a simpler way:
smoothF = smooth_quad(D,c);

x_error = @(x) norm(x-x_star,Inf);

opts = [];
opts.errFcn     = { @(f,x) x_error(x)};
opts.restart = 100;

[ x, out, optsOut ] = tfocs( smoothF, [], [],x0, opts );
% Check that we are within allowable bounds
if out.err(end) < 1e-5
    disp('Everything is working');
else
    error('Failed the test');
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function Tn = tangent_vector(P,Cov)
Tn = logm(P^-0.5*Cov*P^-0.5);

author: alexandrebarachant

repo: covariancetoolbox


text:
function symb_pvec = sdpvar2str(pvec)
%SDPVAR2STR Converts an SDPVAR object to MATLAB string representation
%
% S = SDPVAR2STR(P)
%
% S : String
% P : SDPVAR object

for pi = 1:size(pvec,1)
    for pj = 1:size(pvec,2)
        p = pvec(pi,pj);

        if isa(p,'double')
            symb_p = num2str(p);
        else
            LinearVariables = depends(p);
            x = recover(LinearVariables);
            exponent_p = full(exponents(p,x));
            names = cell(length(LinearVariables),1);
            for i = 1:length(LinearVariables)
                names{i}=['x(' num2str(LinearVariables(i)) ')'];
            end

            symb_p = '';
            if all(exponent_p(1,:)==0)
                symb_p = num2str(getbasematrix(p,0));
                exponent_p = exponent_p(2:end,:);
            end

            for i = 1:size(exponent_p,1)
                coeff = getbasematrixwithoutcheck(p,i);
                switch full(coeff)
                    case 1
                        coeff='+';
                    case -1
                        coeff = '-';
                    otherwise
                        if coeff >0
                            coeff = ['+' num2str2(coeff)];
                        else
                            coeff=[num2str2(coeff)];
                        end
                end
                if strcmp(symb_p,'') & (strcmp(coeff,'+') | strcmp(coeff,'-'))
                        symb_p = [symb_p coeff symbmonom(names,exponent_p(i,:))];
                    else
                symb_p = [symb_p coeff '*' symbmonom(names,exponent_p(i,:))];
            end
            end
            if symb_p(1)=='+'
                symb_p = symb_p(2:end);
            end
        end

        symb_p = strrep(symb_p,'*^0','');
        symb_p = strrep(symb_p,'^0','');
        symb_p = strrep(symb_p,'+*','+');
        symb_p = strrep(symb_p,'-*','-');
        symb_pvec{pi,pj} = symb_p;
    end
end

function s = symbmonom(names,monom)
s = '';
for j = 1:length(monom)
    if monom(j)~=0
        if strcmp(s,'')
            s = [s names{j}];
        else
            s = [s '*' names{j}];
        end
    end
    if monom(j)~=1
        s = [s '^' num2str(monom(j))];
    end
%      if monom(j)>1
%        s = [s '^' num2str(monom(j))];
%    end

end

function s = num2str2(x)
        s = num2str(x);
        if isequal(s,'1')
            s = '';
        end
        if isequal(s,'-1')
            s = '-';
        end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function view_varmap(lab1,valueMap, newgri, sor)
    % view_maxz plots the maxz LTA values calculated
    % with maxzlta.m or other similar values as a color map
    % needs valueMap, gx, gy, stri
    %
    % define size of the plot etc.
    %
    ZG=ZmapGlobal.Data;
    % Prmap = nan(size(valueMap));

    % Find out if figure already exists
    %
    bmap=findobj('Type','Figure','-and','Name','variance-value-map');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(bmap)
        bmap = figure_w_normalized_uicontrolunits( ...
            'Name','variance-value-map',...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

        lab1 = 'b-value:';

        create_my_menu();

        re4 = valueMap;

        colormap(jet)
        ZG.tresh_km = nan;
        minpe = nan;
        Mmin = nan;

    end   % This is the end of the figure setup

    % Now lets plot the color-map of the z-value
    %
    figure(bmap);
    delete(findobj(bmap,'Type','axes'));
    % delete(sizmap);
    ax = reset(gca);
    cla(ax)
    set(ax,'NextPlot','replace')
    watchon;
    set(ax,'visible','off',...
        'FontSize',ZmapGlobal.Data.fontsz.s,...
        'FontWeight','normal',...
        'LineWidth',1,...
        'Box','on','SortMethod','childorder')

    % find max and min of data for automatic scaling
    ZG.maxc = fix(max(valueMap(:))+1);
    ZG.minc = fix(min(valueMap(:)))-1;

    % set values gretaer ZG.tresh_km = nan
    %
    re4 = valueMap;

    % plot image
    %
    orient(ax,'landscape')

    ax.Position = [0.12,  0.10, 0.8, 0.8];
    set(gca,'NextPlot','add')
    pcolor(ax, gx,gy,re4);

    axis(ax, [min(gx) max(gx) min(gy) max(gy)])
    axis(ax,'image')
    set(gca,'NextPlot','add')

    shading(ax, ZG.shading_style);

    % make the scaling for the recurrence time map reasonable

    fix_caxis.ApplyIfFrozen(ax);

    title(ax, sprintf('%s ;  %g to %g',name, t0b, teb),...
        'FontSize',ZmapGlobal.Data.fontsz.s,...
        'Color','k',...
        'FontWeight','normal')

    xlabel(ax, 'Longitude [deg]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel(ax, 'Latitude [deg]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)

    % plot overlay
    %
    set(ax,'NextPlot','add')
    zmap_update_displays();

    set(ax,'NextPlot','add')
    quiver(ax, newgri(:,1),newgri(:,2), -cos(sor(:,SA*2)*pi/180), sin(sor(:,SA*2)*pi/180),0.8,'.', ...
        'LineWidth',1, 'Color','k')

    set(ax,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,...
        'FontWeight','normal',...
        'LineWidth',1,...
        'Box','on',...
        'TickDir','out');

    h1 = ax;
    hzma = ax;

    % Create a colorbar
    %
    h5 = colorbar('horiz');
    set(h5,'Pos',[0.35 0.06 0.4 0.02],...
        'FontWeight','normal',...
        'FontSize',ZmapGlobal.Data.fontsz.s,...
        'TickDir','out')

    axes('position', [0.00,  0.0, 1 1])
    axis('off')
    %  Text Object Creation
    text(...
        'Units','normalized',...
        'Position',[ 0.33 0.06 0 ],...
        'HorizontalAlignment','right',...
        'FontSize',ZmapGlobal.Data.fontsz.s,....
        'FontWeight','normal',...
        'String','Variance');

    % Make the figure visible
    %
    set(gca,'FontSize',ZmapGlobal.Data.fontsz.s,...
        'FontWeight','normal','LineWidth',1,...
        'Box','on','TickDir','out');
    set(gcf,'color','w');
    figure(bmap);
    axes(h1)
    watchoff(bmap)

    %% ui functions
    function create_my_menu()
        add_menu_divider();

        options = uimenu('Label',' Select ');
        uimenu(options,'Label','Refresh ','MenuSelectedFcn',@cb_refresh)
        uimenu(options,'Label','Select EQ in Circle',...
            'MenuSelectedFcn',@select_in_circle)
        uimenu(options,'Label','Select EQ in Circle - Constant R',...
            'MenuSelectedFcn',@cb_select_in_constr)

        uimenu(options,'Label','Select EQ in Polygon -new ',...
            'MenuSelectedFcn',@cb_select_in_polygon_new)

        op1 = uimenu('Label',' Maps ');

        uimenu(op1,'Label','Variance map',...
            'MenuSelectedFcn',@cb_variancemap)
        uimenu(op1,'Label','Resolution map',...
            'MenuSelectedFcn',@cb_resolutionmap)
        uimenu(op1,'Label','Plot map on top of topography ',...
            'MenuSelectedFcn',@cb_plot_on_topography)

        uimenu(op1,'Label','Histogram ','MenuSelectedFcn',@(~,~)zhist())

        add_display_menu(1)
    end

    %% callback functions

    function cb_refresh(~,~)
        view_varmap(lab1, re4);
    end

    function select_in_circle(~,~)
        h1 = gca;
        met = 'ni';
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        circle;
        watchon;
        doinvers_michael;
        watchoff;
    end

    function cb_select_in_constr(~,~)
        h1 = gca;
        met = 'ra';
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        circle;
        watchon;
        doinvers_michael;
        watchoff;
    end

    function cb_select_in_polygon_new(~,~)
        cufi = gcf;
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        selectp;
        watchon;
        doinvers_michael;
        watchoff;
    end

    function cb_variancemap(~,~)
        view_varmap('b-value',r);
    end

    function cb_resolutionmap(~,~)
        view_varmap('Radius', rama);
    end

    function cb_plot_on_topography(~,~)
        dramap_z('stress2','w',valueMap);
    end

end

author: CelsoReyes


text:
function [img, seg, names, counts] = LM2segments(D, imagesize, HOMEIMAGES, HOMELMSEGMENTS)
%
% Transforms all the labelme labels into segmentation masks.
% It takes into account the occlusions, removing labeled pixels that are
% occluded by other objects. Depth ordering is achieved using the function
% LMsortlayers.m
%
% For unlabeled pixels, seg(...) = 0
%
%  [img, seg, names, counts] = LM2segments(D, imagesize, HOMEIMAGES, HOMELMSEGMENTS)
% 'img' and 'seg' are matrices
%
% To precompute the segmentations:
%
%   LM2segments(D(1), [], HOMEIMAGES, HOMELMSEGMENTS) % removing the second
%     argument, makes the segmentation mask to have the same size than the
%     images.
%
% To read the precomputed segmentations:
%   [img, seg, names] = LM2segments(D(1), [], HOMEIMAGES, HOMELMSEGMENTS);
%
%
% HOMEANNOTATIONS = 'http://labelme.csail.mit.edu/Annotations'
% HOMEIMAGES = 'http://labelme.csail.mit.edu/Images'

if nargin == 2 && nargout<3
    % SPECIAL BEHAVIOR
    %
    % seg = LM2segments(annotation, imagesize);
    imgtmp = zeros(imagesize);
    annotationtmp = LMsortlayers(D.annotation, imgtmp);

    % Get segmentation
    [S_instances, classes] = LMobjectmask(annotationtmp, size(imgtmp));
    %classesndx = [annotationtmp.object.namendx];
    %area = squeeze(sum(sum(S_instances,1),2));

    %j = find(area>2); % remove small objects
    %S_instances = S_instances(:,:,j);
    %classesndx = classesndx(j);

    % Assign labels taking into account occlusions!
    Mclasses = zeros([imagesize(1) imagesize(2)]);
    for k = size(S_instances,3):-1:1;
        S_instances(:,:,k) = k*S_instances(:,:,k);
        Mclasses = Mclasses+(Mclasses==0).*S_instances(:,:,k);
    end
    img  = uint16(Mclasses);
    seg = classes;

    return
end

if nargin==4
    precomputed = 1;

    if length(D) == 1
        fileseg = fullfile(HOMELMSEGMENTS, D(1).annotation.folder, [D(1).annotation.filename(1:end-4) '.mat']);
        % read the image and return
        if exist(fileseg, 'file')
            load(fileseg)
            seg = S;
            img = [];
            return
        end
    end
else
    precomputed = 0;
    HOMELMSEGMENTS = '';
end

Nimages = length(D);

% Create list of objects:
n=0;
while 1
    n=n+1;
    if n >length(D)
        [D, names, counts] = LMcreateObjectIndexField(D);
        break
    end

    if isfield(D(n).annotation, 'object')
        if ~isfield(D(n).annotation.object, 'namendx')
            [D, names, counts] = LMcreateObjectIndexField(D);
        else
            [names, counts, imagendx, objectndx, objclass_ndx] = LMobjectnames(D, 'name');
        end
        break
    end
end

% [names, counts, imagendx, objectndx, objclass_ndx] = LMobjectnames(D, 'name');
% for k = 1:length(objclass_ndx)
%     D(imagendx(k)).annotation.object(objectndx(k)).namendx = objclass_ndx(k);
% end
Nobjectclasses = length(names);

if nargout>0
    if Nimages > 1
        % Initalize output variables
        seg = zeros([imagesize(1) imagesize(2) Nimages], 'uint16');
        img = zeros([imagesize(1) imagesize(2) 3 Nimages], 'uint8');
    end
end

if Nimages > 1
    figure
end

for ndx = 1:Nimages
    % Load image
    imgtmp = LMimread(D, ndx, HOMEIMAGES);
    annotation = D(ndx).annotation;
    if size(imgtmp,3)==1; imgtmp = repmat(imgtmp, [1 1 3]); end
    [nrows ncols cc] = size(imgtmp);

    % Scale image so that image box fits tight in image
    if ~isempty(imagesize)
        scaling = max(imagesize(1)/nrows, imagesize(2)/ncols);
        [annotationtmp, imgtmp] = LMimscale(annotation, imgtmp, scaling, 'bicubic');

        % Crop image to final size
        [nr nc cc] = size(imgtmp);
        sr = floor((nr-imagesize(1))/2);
        sc = floor((nc-imagesize(2))/2);
        [annotationtmp, imgtmp] = LMimcrop(annotationtmp, imgtmp, [sc+1 sc+imagesize(2) sr+1 sr+imagesize(1)]);

        Mclasses = zeros([imagesize(1) imagesize(2)]);
    else
        annotationtmp = annotation;

        Mclasses = zeros([size(imgtmp,1) size(imgtmp,2)]);
    end

    if isfield(annotationtmp, 'object')
        % Sort layers
        annotationtmp = LMsortlayers(annotationtmp, imgtmp);

        % Get segmentation
        [S_instances, classes] = LMobjectmask(annotationtmp, size(imgtmp));
        classesndx = [annotationtmp.object.namendx];
        area = squeeze(sum(sum(S_instances,1),2));

        j = find(area>2); % remove small objects
        S_instances = S_instances(:,:,j);
        classesndx = classesndx(j);

        % Assign labels taking into account occlusions!
        for k = size(S_instances,3):-1:1;
            S_instances(:,:,k) = classesndx(k)*S_instances(:,:,k);
            Mclasses = Mclasses+(Mclasses==0).*S_instances(:,:,k);
        end
    end

    if nargout>0
        % Store values
        seg(:,:,ndx)   = uint16(Mclasses);
        img(:,:,:,ndx) = imgtmp;
    end

    % Save gist if a HOMELMSEGMENTS file is provided
    if precomputed
        I = imgtmp;
        S = uint16(Mclasses);
        S_instances = uint16(S_instances);
        mkdir(fullfile(HOMELMSEGMENTS, D(ndx).annotation.folder))
        fileseg = fullfile(HOMELMSEGMENTS, D(ndx).annotation.folder, [D(ndx).annotation.filename(1:end-4) '.mat']);
        if ~isempty(imagesize)
            save (fileseg, 'I', 'S', 'names', 'S_instances')
        else
            save (fileseg, 'S', 'names', 'S_instances')
        end
    end

    % Visualization
    if Nimages > 1
        subplot(121)
        image(imgtmp); axis('equal'); axis('tight');
        title(sprintf('%d (out of %d)', ndx, Nimages))
        subplot(122)
        image(mod(Mclasses+1,256)); axis('equal'); axis('tight');
        colormap([0 0 0; hsv(min(Nobjectclasses+1,256))])
        drawnow
    end
end
%
% if Nimages > 1
%     if nargout > 0
%         % plot stats
%         figure
%         subplot(121)
%         loglog(sort(counts, 'descend'))
%         xlabel('count rank')
%         ylabel('Number of instances')
%         axis('tight')
%         pixelcounts = hist(single(seg(:)), 0:single(max(seg(:))));
%         unlabeled = pixelcounts(1);
%         pixelcounts = pixelcounts(2:end); % remove unlabeled pixels;
%         subplot(122)
%         loglog(sort(pixelcounts, 'descend'))
%         xlabel('area rank')
%         ylabel('Number of pixels')
%         axis('tight')
%         title(sprintf('%d categories', length(names)))
%
%         Ns = min(20, size(seg,3));
%         figure
%         montage(reshape(uint8(mod(seg(:,:,1:Ns),256)), [imagesize(1) imagesize(2) 1 Ns]))
%         colormap(gray(min(Nobjectclasses ,256)))
%     else
%         figure
%         loglog(sort(counts, 'descend'))
%         xlabel('count rank')
%         ylabel('Number of instances')
%         axis('tight')
%     end
% end

author: CSAILVision

repo: LabelMeToolbox


text:
function  kern = sdrbfKernParamInit(kern)

% SDRBFKERNPARAMINIT SDRBF kernel initialization
% FORMAT
% DESC
% Initializes the switching dynamical latent force model structure with
% some initial parameters for the RBF underlying kernels. The initial
% parameters are passed through an option in kern.
% RETURN kern : the kernel structure with the default parameters placed in.
% ARG kern : the kernel structure which requires initialisation.
%
% SEEALSO : kernCreate, kernParamInit, rbfKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

% Create the basic structure based on the lfm kernel.
kern = rbfKernParamInit(kern);

if isfield(kern, 'options') && isfield(kern.options, 'isNormalised') && ...
        kern.options.isNormalised,
    kern.isNormalised = 1;
end

if isfield(kern, 'options') && isfield(kern.options, 'nIntervals')
    kern.nIntervals = kern.options.nIntervals;
    if isfield(kern.options, 'nlfPerInt')
        kern.nlfPerInt = kern.options.nlfPerInt;
    else
        kern.nlfPerInt = 1;
    end
    kern.inverseWidth = ones(kern.nlfPerInt, kern.nIntervals); % Total number of inverse widths.
    % An option for the initialization of the switching times
    if isfield(kern.options, 'switchingTimes')
        if kern.nIntervals == length(kern.options.switchingTimes)
            kern.switchingTimes = kern.options.switchingTimes;
        else
            error('The number of intervals does not match the information of the swicthing time points')
        end
    else
        partition = linspace(-0.1,1, kern.nIntervals + 1);
        kern.switchingTimes = partition(1:end-1);
    end
else
    kern.nIntervals = 1;
    kern.switchingTimes = -0.1;
    warning('LFM:Instead:SDLFM', 'Use the LFM kernel instead.')
end

% Number of parameters computed as: inverseWidth (= nlf*nIntervals) and the
% switching points (= nIntervals). The variance, which is a parameter in
% the original rbf kernel, is not counted as a parameter. It's left here,
% so that it can be used when making sdrbf kernel computations and gradients.

kern.nParams = kern.nlfPerInt*kern.nIntervals + kern.nIntervals;
dimParam = [kern.nlfPerInt*kern.nIntervals kern.nIntervals];
kern.inverseWidthIndx = 1:dimParam(1);
kern.switchingTimesIndx = dimParam(1)+1:sum(dimParam(1:2));
%kern.transforms.index = 1:(kern.nlfPerInt*kern.nIntervals);
kern.transforms.index = [1:(kern.nlfPerInt*kern.nIntervals) (kern.nlfPerInt*kern.nIntervals)+2:kern.nParams];
kern.transforms.type = optimiDefaultConstraint('positive');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function DEM = ADEM_eyeblink(OPTION)
% Simulation of eyeblink conditioning
% FORMAT DEM = ADEM_eyeblink(OPTION)
%
%  OPTION:
%     case{'EYEBLINK'}   : spontaneous eye blinking
%     case{'AIRPUFF'}    : unconditioned eyeblink response to air puff
%     case{'STARTLE'}    : unconditioned  startle response to a sound
%     case{'TRACE'}      : trace conditioning to the sound
%     case{'DELAY'}      : delay conditioning to the sound
%     case{'EXTINCTION'} : extinction of trace conditioning to sound
%
%__________________________________________________________________________
%
% This demonstration routine illustrates Pavlovian learning under active
% inference. It uses the eyeblink conditioning paradigm to model startle
% responses and the subsequent acquisition of an eyeblink - under delay and
% trace learning. The various options above determine the nature of the
% simulation (or edit the OPTION below). The generative model, in this
% example, starts with a heteroclinic cycle with an inset. The cycle per se
% generates autonomous eyeblinks periodically, while the inset is
% activated by a conditioned stimulus (CS). The subsequent unstable fixed
% points play the role of an echo-state and enables the learning or
% association of a high-level hidden cause with subsequent unconditioned
% responses (UR).
%
% In active inference, an unconditioned response corresponds to a prior
% belief that a hidden state will generate action and the unconditioned
% stimulus (US). Pavlovian conditioning is the learning of the Association
% between a conditioned stimulus (CS) and the unconditioned response. The
% dynamics entailed by the heteroclinic cycle enable trace conditioning,
% which may be related to hippocampal function.
%
% In this example, there are two levels with the hidden states at the first
% level modelling beliefs about eyeblinks, unconditioned responses and
% unconditioned stimuli.  Proprioceptive predictions are generated by
% beliefs about ensuing eyeblinks and unconditioned responses (which
% also predict the conditioned stimulus. Hidden states at the second level
% embody a sense of time through Lotka-Volterra dynamics. Successive epochs
% of time are passed to the first level via a softmax transform. Learning
% corresponds to Hebbian plasticity (that minimises free energy) in the
% connections between the state unit encoding expectations about a UR and
% expectations about the CS (for delay conditioning) and heteroclinic
% states (for trace conditioning): see the functions at the end of this
% routine.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_eyeblink.m 7679 2019-10-24 15:54:07Z spm $

% paradigm and stimuli
%==========================================================================

% options:
%--------------------------------------------------------------------------
if ~nargin, OPTION = 'STARTLE'; end

P.x    = 0;                                          % connection strengths
P.v    = [0 0 0];                                    % connection strengths
PC     = 0;                                          % prior covariance
V      = [8 4 0];                                    % sensory precision

switch OPTION

    case{'EYEBLINK'}  % spontaneous eye blinking
        N      = 256;                                % number of time bins
        NT     = 1;                                  % number of trials
        C(:,1) = sparse(1,N);                        % CS: loud sound
        C(:,2) = sparse(1,N);                        % US: air puff

    case{'AIRPUFF'}   % unconditional response to US
        N      = 256;                                % number of time bins
        NT     = 1;                                  % number of trials
        CS     = 128;                                % onset of CS (bins)
        C(:,1) = sparse(1,N);                        % CS: loud sound
        C(:,2) = exp(-((1:N) - CS).^2/(2*4^2));      % US: air puff

    case{'STARTLE'}   % startle response to CS
        N      = 256;                                % number of time bins
        NT     = 1;                                  % number of trials
        CS     = 128;                                % onset of CS (bins)
        C(:,1) = exp(-((1:N) - CS).^2/(2*4^2));      % CS: loud sound
        C(:,2) = sparse(1,N);                        % US: air puff

    case{'DELAY'}     % delay learning of association
        N      = 64;                                 % number of time bins
        NT     = 16;                                 % number of trials
        CS     = 48;                                 % onset of CS (bins)
        C(:,1) = exp(-((1:N) - CS).^2/(2*8^2));      % CS: loud sound
        C(:,2) = exp(-((1:N) - CS - 4).^2/(2*4^2));  % US: air puff
        PC     = 1;                                  % enable learning

    case{'TRACE'}     % trace learning of association (over time)
        N      = 128;                                % number of time bins
        NT     = 32;                                 % number of trials
        CS     = 48;                                 % onset of CS (bins)
        C(:,1) = exp(-((1:N) - CS).^2/(2*4^2));      % CS: loud sound
        C(:,2) = exp(-((1:N) - CS - 32).^2/(2*4^2)); % US: air puff
        PC     = 1;                                  % enable learning

    case{'EXTINCTION'} % normalextinction
        N      = 128;                                % number of time bins
        NT     = 8;                                  % number of trials
        CS     = 48;                                 % onset of CS (bins)
        C(:,1) = exp(-((1:N) - CS).^2/(2*4^2));      % CS: loud sound
        C(:,2) = sparse(1,N);                        % US: air puff
        P.x    = 0;                                  % connection strengths
        P.v    = [0 0 0.6];                          % connection strengths
        PC     = 1;                                  % enable learning
        V      = [8 0 0];                            % pirotoxin lesion
        V      = [8 2 0];                            % sensory precision

    otherwise
end

% generative process
%==========================================================================
M(1).E.n  = 2;
M(1).E.d  = 2;
M(1).E.s  = 1;
M(1).E.nN = 2;
M(1).E.nE = 1;

% level 1
%--------------------------------------------------------------------------
G(1).g  = @(x,v,a,P) G1g(x,v,a,P);
G(1).V  = exp(16);                                % error precision
G(1).U  = exp(4);                                 % motor gain

% level 2
%--------------------------------------------------------------------------
G(2).v  = [0;0];                                  % stimulus (CS and US)
G(2).a  = 0;                                      % action
G(2).V  = exp(16);

% generative model
%==========================================================================

% positions associated with each state (on unit circle)
%--------------------------------------------------------------------------
x.CS    = 0;
x.UR    = 0;
x.EB    = 0;

% level 1
%--------------------------------------------------------------------------
M(1).x  = x;
M(1).f  = @(x,v,P) M1f(x,v,P);
M(1).g  = @(x,v,P) M1g(x,v,P);
M(1).pE = P;
M(1).pC = PC;
M(1).W  = diag(exp([4 0 4]));                     % cerebellar (IPN) lesion
M(1).W  = exp(4);                                 % error precision
M(1).V  = exp(V);                                 % sensory attenuation

% for use with spm_ALAP (for state-dependent sensory attenuation)
%--------------------------------------------------------------------------
M(1).E.method.x = 0;
M(1).ph = @(x,v,h,M) [8 4 4] - 2*(x.UR + x.EB);   % sensory attenuation

% level 2
%--------------------------------------------------------------------------
P = [
    +1   -0.5  0.2  0.2  0.2  0.2;
    +0.5  1   -0.5  0    0    0  ;
    -0.5  0.5  1   -0.5  0    0  ;
    -0.5  0    0.5  1    0    0.5;
    -0.5  0    0    0.5  1   -0.5;
    -0.5  0    0   -0.5  0.5  1  ] - 1;

M(2).x  = [0 0 0 4 0 0]' - 4;
M(2).f  = @(x,v,P) spm_lotka_volterra(x,v,P);
M(2).g  = @(x,v,P) spm_softmax(x);
M(2).pE = P;
M(2).V  = exp(0);                                 % hippocampal lesion
M(2).V  = exp(4);                                 % error precision
M(2).W  = exp(4);                                 % error precision

% level 3
%--------------------------------------------------------------------------
M(3).v  = 0;                                      % inputs (null)
M(3).V  = exp(16);

% ADEM
%==========================================================================
DEM.U  = sparse(N,1);
DEM.C  = C;
DEM.G  = G;
DEM.M  = M;
DEM.db = 0;
for i  = 1:NT

    % integrate active inference scheme
    %----------------------------------------------------------------------
    DEM  = spm_ADEM(DEM);

    if i == 1
        spm_figure('GetWin','Before learning');
    else
        spm_figure('GetWin','After learning');
    end
    spm_DEM_qU(DEM.qU,DEM.pU);

    if NT == 1, break, end

    % latency and vigour of CR
    %----------------------------------------------------------------------
    a       = DEM.qU.a{2};
    a       = abs(a(1:min(end,(CS + 32))));
    if max(a) > 1/4;
        q = sum((1:length(a)).*a/sum(a));
    else
        q = NaN;
    end

    CR(1,i) = q;
    CR(2,i) = max(a);

    % Baysian belief update
    %----------------------------------------------------------------------
    p           = DEM.qP.P{1};
    DEM.M(1).pE = p;
    qP(:,i)     = spm_vec(p);

    % plot
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 1'); clf

    subplot(2,2,1),cla
    plot((1:i),CR(1,:) - CR(1,1)), hold on
    plot((1:i),qP',':b'),          hold on
    plot((1:i),CR(2,:),'--'),      hold off
    set(gca,'XLim',[0 (NT + 1)])
    title('latency  and vigour','FontSize',16);
    xlabel('trials')
    ylabel('response')
    axis square

    subplot(2,2,2),cla
    imagesc(qP)
    title('plasticity','FontSize',16);
    set(gca,'XLim',[0.5 (NT + .5)])
    xlabel('trials')
    ylabel('connection')
    axis square

end

% functions of generative process and model (at the first level)
%==========================================================================
function s = G1g(x,v,a,P)                        % sensory process
s.CS = v(1);                                     % CS: amplitude sound
s.US = v(2);                                     % US: air puff
s.R  = a;                                        % proprioception

function s = M1g(x,v,P)                          % sensory prediction
s.CS = x.CS;                                     % prediction of CS
s.US = x.UR;                                     % prediction of US
s.R  = x.UR + x.EB;                              % prediction of UR

function f = M1f(x,v,P)                          % dynamics of states
f.CS = v(1) - x.CS;                              % dynamics of CS
f.UR = P.v*v(1:3) + P.x*x.CS - x.UR;             % dynamics of UR construct
f.EB = v(6) - x.EB;                              % dynamics of eyeblinks

author: spm


text:
function [qp]=zn2fr(input)
% [qp]=zn2fr(input)
% Ziegler-Nichols PID controller for processes of 2nd order.
% This function computes parameters of the controller (q0, q1, q2, p1, p2).
% Controller is based on forward rectangular method of discretization.
% Transfer function of the controller is as follows:
%
%            q0 + q1*z^-1 + q2*z^-2     q0 + q1*z^-1 + q2*z^-2
% G(z^-1) = ------------------------ = ------------------------
%                  1 - z^-1              1 + p1*z^-1 + p2*z^-2
%
% where p1=-1 and p2=0.
%
% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2
% Gs(z^-1) = -----------------------
%             1 + a1*z^-1 + a2*z^-2
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... b2
%                  input(5) ... sample time T0
% Output: qp ... controller parameters
%                qp(1) ... q0
%                qp(2) ... q1
%                qp(3) ... q2
%                qp(4) ... -1   (p1 of the controller)
%                qp(5) ... 0    (p2 of the controller)

a1 = input(1);
b1 = input(2);
a2 = input(3);
b2 = input(4);
T0 = input(5);

% compute ultimate gain and frequency
[Kpu, Tu] =  ultim([b1 b2],[a1 a2],T0);

Kp = 0.6*Kpu;
Ti = Tu/2;
Td = Tu/8;

q0 = Kp*(1+Td/T0);
q1 = -Kp*(1-T0/Ti+2*Td/T0);
q2 = Kp*(Td/T0);
p1 = -1;
p2 = 0;

qp=[q0; q1; q2; p1; p2];

author: Sable

repo: mcbench-benchmarks


text:
function display(X)
%DISPLAY (overloaded)

switch(X.typeflag)
    case {0,9,40}
        n = X.dim(1);
        m = X.dim(2);
        if (n*m==1)

            vars = depends(X);
            if any(ismember(vars,yalmip('extvariables')))
                classification = 'Nonlinear scalar ';
            else
                linearbilinearquadraticsigmonial = is(X,'LBQS');
                if linearbilinearquadraticsigmonial(1)
                    classification = 'Linear scalar ';
                elseif linearbilinearquadraticsigmonial(4)
                    classification = 'Signomial scalar ';
                elseif linearbilinearquadraticsigmonial(2)
                    classification = 'Bilinear scalar ';
                elseif linearbilinearquadraticsigmonial(3)
                    classification = 'Quadratic scalar ';
                else
                    classification = 'Polynomial scalar ';
                end
            end

            if ~isreal(X.basis)
                classification = [classification '(complex'];
            else
                classification = [classification '(real'];
            end

            if is(X,'compound')
                if ~isequal(X.extra.opname,'')
                    classification = [classification ', models ''' X.extra.opname ''''];
                end
            end

            if ~islinear(X)
                variables = getvariables(X);
                monomtable = yalmip('monomtable');
                s = sum((monomtable(variables,:)),2);
               % s = sum(full(monomtable(variables,:)),2);
                if ((nnz(getbasematrix(X,0))==0) ) && all(s==s(1))
                    classification = [classification ', homogeneous'];
                end
            end
            if any(ismember(depends(X),yalmip('intvariables')))
                classification = [classification ', integer'];
            else
                if any(ismember(depends(X),yalmip('binvariables')))
                    classification = [classification ', binary'];
                else
                    if any(ismember(depends(X),yalmip('uncvariables')))
                        classification = [classification ', uncertain'];
                    else
                        if any(ismember(depends(X),yalmip('semicontvariables')))
                        classification = [classification ', semi-continuous'];
                        end
                    end
                end
            end

            nvars = length(depends(X));
            if nvars == 1
                classification = [classification ', ' num2str(nvars) ' variable)'];
            else
                classification = [classification ', ' num2str(nvars) ' variables)'];
            end

            if any(ismember(depends(X),yalmip('parvariables')))
                classification = [classification ', parametric'];
            end

            try
                v = value(X);
            catch
                v = nan;
            end
             if ~isnan(v)
                classification = [classification '\nCurrent value: ' num2str(v)];
            end

            B = getbase(X);
            [ii,jj,ss1] = find(real(getbase(B)));
            [ii,jj,ss2] = find(imag(getbase(B)));
            ss = [ss1;ss2];
            DynamicalRange = [num2str( min(abs(ss))) ' to ' num2str( max(abs(ss)))];
            classification = [classification '\nCoefficients range: ' DynamicalRange];
            fprintf([classification '\n']);
        else

            vars = depends(X);
            if any(ismember(vars,yalmip('extvariables')))
                classification = 'Nonlinear matrix variable ';
            else
                if islinear(X)
                    classification = 'Linear matrix variable ';
                elseif is(X,'sigmonial')
                    classification = 'Signomial matrix variable ';
                elseif is(X,'bilinear')
                    classification = 'Bilinear matrix variable ';
                elseif is(X,'quadratic')
                    classification = 'Quadratic matrix variable ';
                else
                    classification = 'Polynomial matrix variable ';
                end
            end

            if isreal(X.basis)
                if issymmetric(X)
                    info = ' (symmetric, ';
                else
                    info = ' (full, ';
                end;
                info = [info 'real'];
            else
                if issymmetric(X)
                    info = ' (symmetric, ';
                elseif ishermitian(X)
                    info = ' (hermitian, ';
                else
                    info = ' (full, ';
                end;
                info = [info 'complex'];
            end;

            if X.typeflag==9
                info = [info ', KYP'];
            end

            if X.typeflag==40
                info = [info ', Generalized KYP'];
            end

            if any(ismember(depends(X),yalmip('intvariables')))
                info = [info ', integer'];
            else
                if any(ismember(depends(X),yalmip('binvariables')))
                    info = [info ', binary'];
                else
                    if any(ismember(depends(X),yalmip('uncvariables')))
                        info = [info ', uncertain'];
                    end
                end
            end

            if any(ismember(depends(X),yalmip('parvariables')))
                info = [info ', parametric'];
            end

            xvars = depends(X);
            nvars = length(xvars);
            if nvars == 1
                info = [info ', ' num2str(nvars) ' variable)'];
            else
                info = [info ', ' num2str(nvars) ' variables)'];
            end

            n = X.dim(1);
            m = X.dim(2);

            if (n<100) && (n==m)
                x = recover(xvars);
                if ~any(any(isnan(value(x))))
                    doubleX = double(X);
                    try
                    eigX = eig(doubleX);
                    info = [info '\nEigenvalues between [' num2str(min(eigX)) ',' num2str(max(eigX)) ']'];
                    catch
                    end
                end
            elseif n~=m
                x = recover(xvars);
                if ~any(any(isnan(value(x))))
                    doubleX = value(X);
                    try
                        info = [info '\nValues in range [' num2str(min(min(doubleX))) ',' num2str(max(max(doubleX))) ']'];
                    catch
                    end
                end
            end

            B = getbase(X);
            [ii,jj,ss1] = find(real(getbase(B)));
            [ii,jj,ss2] = find(imag(getbase(B)));
            ss = [ss1;ss2];
            if any(ss)
                DynamicalRange = [num2str( min(abs(ss))) ' to ' num2str( max(abs(ss)))];
            else
                DynamicalRange = ['0 to 0'];
            end

            info = [info '\nCoefficient range: ' DynamicalRange];
            fprintf([classification num2str(n) 'x' num2str(m) info '\n']);
        end;
    case 1
        disp('Relational object');
    case 3
        disp('Relational object');
    case 4
        disp('Relational object');
    case 5
        disp('Cone object');
    case 6
        disp('logdet object');
    case 7
        disp('Integrality constraint');
    case 10
        disp('Eigenvalue object');
    case 11
        disp('SOS object');
    case 20
        disp('Power cone object');

    otherwise
end


text:
function tetrahedron_arbq_rule_test03 ( degree, n, header )

%*****************************************************************************80
%
%% TETRAHEDRON_ARBQ_RULE_TEST03 gets a rule and creates GNUPLOT input files.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    10 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired total polynomial degree exactness
%    of the quadrature rule.  0 <= DEGREE <= 15.
%
%    Input, integer N, the number of nodes to be used by the rule.
%
%    Input, string HEADER, an identifier for the filenames.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_ARBQ_RULE_TEST03\n' );
  fprintf ( 1, '  Get a quadrature rule for the tetrahedron.\n' );
  fprintf ( 1, '  Set up GNUPLOT graphics input.\n' );
  fprintf ( 1, '  Polynomial exactness degree DEGREE = %d\n', degree );
%
%  Retrieve a symmetric quadrature rule.
%
   [ x, w ] = tetrahedron_arbq ( degree, n );
%
%  Create files for input to GNUPLOT.
%
  tetrahedron_arbq_gnuplot ( n, x, header );

  return
end


text:
% Make some synthetic quaternions
syms   q0 q1 q2 q3   p0 p1 p2 p3   real

q = [q0 q1 q2 q3]';
p = [p0 p1 p2 p3]';

vgg_assert_equal('vgg_quat_mul(p,q)', 'vgg_quat_matrix(p) * q', 1e-8,1);

% Generate a random unit quaternion.
% This represents a 3D rotation.
q = randn(4,1);
q = q / norm(q);

% Now convert it to a rotation matrix,
% and convert back to ensure these routines
% are self-inverse
R = vgg_quat_rotation_matrix(q);
qback = vgg_quat_from_rotation_matrix(R);

vgg_assert_equal(q, qback, -1e-8,1);

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function Z = xor(X,Y)
%XOR Logical EXCLUSIVE OR for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Z = tenfun(@xor,X,Y);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Infrared Small-Target Detection Using Multiscale Gray Difference Weighted Image Entropy
clearvars;
close all;
clc;
wflag = 1;
for kk = 1
    fold = '.\data\';% 27 images
    try
        img = imread([fold, num2str(kk), '.jpg']);
    catch
        img = imread([fold, num2str(kk), '.bmp']);
    end
    img = img(:,:,1);
    re = mgdwe(img);
    bw = bwfunc(re);
    if wflag
        refold = '.\';
        if exist(refold,'dir')==0
            mkdir(refold);
        end
        imwrite(uint8(img), [refold, num2str(kk), '1.png']);
        imwrite(uint8( mat2gray(re) * 255), [refold, num2str(kk), '2.png']);
        imwrite(bw, [refold, num2str(kk), '3.png']);
    end
end

author: daxjuanxiong

repo: infrared-small-target-detection


text:
function hand_plot ( )

%*****************************************************************************80
%
%% HAND_PLOT plots the hand data.
%
%  Discussion:
%
%     This program assumes that the file 'HAND_NODES.TXT' is available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625.
%

%
%  Read the data.
%
  xy = load ( 'hand_nodes.txt' );
%
%  Make XY an array of column vectors.
%
  xy = xy';
%
%  Repeat the first column at the end so the polygon closes.
%
  xy = [ xy, [ xy(:,1) ] ];
%
%  Clear the graphics frame.
%
  clf

  plot ( xy(1,:), xy(2,:), 'Color', 'r', 'LineWidth', 2 );
  hold on;
  plot ( xy(1,:), xy(2,:), 'b.', 'MarkerSize', 15 );
  axis equal
  grid on
  title ( 'Hand data and straight line interpolant' )

  hold off

  return
end


text:
function prob_test1306 ( )

%*****************************************************************************80
%
%% PROB_TEST1306 tests QUASIGEOMETRIC_MEAN, *_SAMPLE, *_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2009
%
%  Author:
%
%    John Burkardt
%
  sample_num = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PROB_TEST1306\n' );
  fprintf ( 1, '  For the Quasigeometric PDF:\n' );
  fprintf ( 1, '  QUASIGEOMETRIC_MEAN computes the mean;\n' );
  fprintf ( 1, '  QUASIGEOMETRIC_SAMPLE samples;\n' );
  fprintf ( 1, '  QUASIGEOMETRIC_VARIANCE computes the variance.\n' );

  a = 0.4825;
  b = 0.5893;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );

  check = quasigeometric_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PROB_TEST1306 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = quasigeometric_mean ( a, b );
  variance = quasigeometric_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : sample_num
    [ x(i), seed ] = quasigeometric_sample ( a, b, seed );
  end

  mean = i4vec_mean ( sample_num, x );
  variance = i4vec_variance ( sample_num, x );
  xmax = max ( x(1:sample_num) );
  xmin = min ( x(1:sample_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', sample_num );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %6d\n', xmax );
  fprintf ( 1, '  Sample minimum =  %6d\n', xmin );

  return
end


text:
function [rgb] = bg_rgba2rgb(bg, rgba, varargin)

% BG_RGBA2RGB overlays a transparency masked colored image on a colored background,
% and represents the result as an RGB matrix.
%
% Use as:
%   rgb = bg_rgba2rgb(bg, rgba)
%
% or
%   rgb = bg_rgba2rgb(bg, rgba, cmap, clim, alpha, amap, alim);
%
% When 2 input arguments are supplied:
%   bg   = Nx3 matrix of background rgb-coded color-values, or MxNx3
%   rgba = Nx4 matrix of rgb + alpha values, or MxNx4
%
% When 7 input arguments are supplied:
%   bg   = Nx3 matrix, Nx1 vector, 1x3 vector, MxN, or MxNx3.
%   rgba = Nx1 vector with 'functional values', or MxN.
%   cmap = Mx3 colormap, or MATLAB-supported name of colormap
%   clim = 1x2 vector denoting the color limits
%   alpha = Nx1 vector with 'alpha values', or MxN
%   amap = Mx1 alphamap, or MATLAB -supported name of alphamap ('rampup/down', 'vup/down')
%   alim = 1x2 vector denoting the opacity limits

if numel(varargin)==0
  siz1 = size(bg);
  siz2 = size(rgba);
  assert(isequal(siz1(1:end-1),siz2(1:end-1)), 'number of data points should be the same');
  assert(siz1(end)==3 && siz2(end)==4,         'inconsistent input size');

  bg   = reshape(bg,   [], 3);
  rgba = reshape(rgba, [], 4);
  rgb  = do_conversion(bg, rgba);
  rgb  = reshape(rgb, siz1);
else
  % this requires the data to be converted into rgb values irst, and
  % needs a cmap + clim, and an alpha, alphamap and alim
  if numel(varargin)~=5
    error('if a vector of color data is supplied, more input arguments are required');
  end

  % colormap handling
  cmap = varargin{1};
  if ischar(cmap)
    cmap = colormap(cmap);
  elseif size(cmap,2)~=3
    error('invalid specification of colormap, should be nx3');
  end
  clim  = varargin{2}; if isempty(clim), clim(1) = min(rgba(:)); clim(2) = max(rgba(:)); end
  alpha = varargin{3};
  amap  = varargin{4};
  alim  = varargin{5}; if isempty(alim), alim(1) = min(alpha(:)); alim(2) = max(alpha(:)); end

  % deal with the color data
  dat     = rgba;
  siz     = size(dat);
  dat     = reshape(dat, [], 1);
  alpha   = reshape(alpha, [], 1); % assume to be same siza as input data!
  finvals = isfinite(dat);

  rgba              = zeros(size(dat,1),4);
  rgba(finvals,1:3) = dat2rgb(dat(finvals), cmap, clim);

  finvals          = isfinite(alpha);
  rgba(finvals,4)  = alpha2a(alpha(finvals), amap, alim);
  rgba(~finvals,4) = 0;

  % deal with the background: allow it to be 1x3, i.e. a single color per
  % pixel.
  if isequal(size(bg),[1 3])
    bg = permute(repmat(bg(:), [1 siz]), [1+(1:numel(siz)) 1]);
  else
    siz_bg = size(bg);
    if isequal(siz_bg, siz)
      % make bg rgb
      if all(bg(:)<=1) && all(bg(:)>=0)
        % don't scale
      else
        bg_min = min(bg(:));
        bg_max = max(bg(:));
        bg     = (bg-bg_min)./(bg_max-bg_min);
      end
      bg = repmat(bg, [ones(1,ndims(bg)) 3]);
    else
      %FIXME
    end
  end
  bg  = reshape(bg, [], 3);

  if numel(siz)==2 && siz(2)==1, siz = siz(1); end
  rgb = do_conversion(bg, rgba);
  rgb = reshape(rgb, [siz 3]);
end

function rgb = do_conversion(bg, rgba)

rgb = zeros(size(rgba,1),3);
a_  = 1-rgba(:,4);
a   = rgba(:,4);

for k = 1:3
    rgb(:,k) = bg(:,k).*a_  + rgba(:,k).*a;
end

function rgb = dat2rgb(dat, cmap, clim)

dat(end+1) = clim(1);
dat(end+1) = clim(2); % add the extremes to be sure that they are included

dat(dat<clim(1)) = clim(1);
dat(dat>clim(2)) = clim(2);

% scale between 0 and 1
dat = dat-min(dat);
dat = dat/max(dat);

ind = round(dat.*(size(cmap,1)-1))+1;
rgb = cmap(ind(1:end-2),:);

function a = alpha2a(alpha, amap, alim)

alpha(end+1) = alim(1);
alpha(end+1) = alim(2);

alpha(alpha<alim(1)) = alim(1);
alpha(alpha>alim(2)) = alim(2);

if ischar(amap)
    switch amap
        case 'rampup'
            a = alpha - min(alpha);
            a = a./max(a);
        case 'rampdown'
            a = alpha - min(alpha);
            a = a./max(a);
            a = 1-a;
        case 'vup'
            a = alpha - min(alpha);
            a = a./max(a);
            a = 1 - 2.*abs(a - 0.5);
        case 'vdown'
            a = alpha - min(alpha);
            a = a./max(a);
            a = 2.*abs(a - 0.5);
        otherwise
            error('unknown alphamap specified');
    end
else
    amap = amap(:);

    % assume it to be a vector that has M elements, scaled between 0 and 1
    a = alpha-min(alpha);
    a = a/max(a);
    ind = round(a.*(size(amap,1)-1))+1;
    a   = amap(ind(1:end),:);
end
a = a(1:end-2);


text:
function bnet = learn_params_dbn(bnet, data)
% LEARN_PARAM_DBN Estimate params of a DBN for a fully observed model
% bnet = learn_params_dbn(bnet, data)
%
% data(i,t) is the value of node i in slice t (can be a cell array)
% We currently assume there is a single time series
%
% We set bnet.CPD{i} to its ML/MAP estimate.
%
% Currently we assume each node in the first 2 slices has its own CPD (no param tying);
% all nodes in slices >2 share their params with slice 2 as usual.

[ss T] = size(data);

% slice 1
for j=1:ss
  if adjustable_CPD(bnet.CPD{j})
    fam = family(bnet.dag,j);
    bnet.CPD{j} = learn_params(bnet.CPD{j}, data(fam,1));
  end
end

% slices 2:T
% data2(:,t) contains [data(:,t-1); data(:,t)].
% Then we extract out the rows corresponding to the parents in the current and previous slice.
data2 = [data(:,1:T-1);
	 data(:,2:T)];
for j=1:ss
  j2 = j+ss;
  if adjustable_CPD(bnet.CPD{j2})
    fam = family(bnet.dag,j2);
    bnet.CPD{j2} = learn_params(bnet.CPD{j2}, data2(fam,:));
  end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
%% getVoronoiStats
% Below is a demonstration of the features of the |getVoronoiStats| function

%%
clear; close all; clc;

%% Syntax
% |voronoiStats=getVoronoiStats(Fv_cell,Vv,np,numBins);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function  [cct, ta1, tetael] = srcintf1(la, nn1, ct4, src4, src5, ctc0, ctad, nprimtilt, lmax, K, stfname)
% srcintf1.m integrates the source including polarisation over conformal time using the approach of Seljak :
% integration over the photon past light cone, see M Zaldarriaga et al., ApJ nr 494, 491 (1998) for K > 0
% the function needs an l-range (la), a k-range (nn1) a ctime-range (ct4), the present conformal time (ctc0),
% decoupling time (ctad), the source in conformal time for the temperature and polarisation, resp
% src4(1..lct4, 1..lnn1), src5(1..lct4, 1..lnn1), the index of the primary power spectrum (nprimtilt),
% the maximum la of the anistropy spectrum (lmax) and the curvature energy content from the Friedman equation (K),
% the startfile name of the ultra-spherical function (stfname).
%
% output: ctt is the square of the temperature spectrum (1) resp the E-polarisation spectrum (2)
% resp the cross-correlation of temperature and polarisation (3), tetatl, tetael temperature resp, polarisation
% temperature for all k1 values (for the last la).
%
% D Vangheluwe 13 june 2005
% remark 1a: attention :nn1 is an integer array :wavenumber array k1 is made from it by taking k1 = nn1 * sqrt(K)
% remark 1: we use ultra-spherical bessel values which have two parameters and must be found by integration:
%   for the calculation of the ultra-spherical bessel functions see: cmb/usphint.m
% remark 2: for the ode equation of the u-function, see equation (36) of Zaldarriaga & Seljak.
% remark 3: take attention the value of nn1 should lie within the range of kb1 values : there is no check!!
%   see spline of start  values x10 and pbd0.
% remark 4: how do we get the ultra-spherical bessel function values? see development and test routine usphpar.m
%  and usphpar1.m
% remark 5: this function is a fast version of srcpint.m (the source need not to be in the loop
% if there is enough memory, probably 80Mbytes will be used)
% remark6 on 13 june: this routine has been obtained by modifying srcpintf.m according to usphpar1.m

if (K < 0),
   error('the space curvature constant K should be > 0 for this routine');
   return;
end
kk0 = sign(K);
lla = size(la, 2);
k1 = nn1 * sqrt(K);
lk1 = size(k1, 2);
lct4 = size(ct4, 2);
k1step = k1(2) - k1(1);
ct4step = ct4(2) - ct4(1);
k1max = k1(lk1);
k1min = k1(1);
ct4sym = 0.5 * pi/sqrt(K);

% if necessary adapt the source to the period of the ultra-spherical function for K >0
ct0 = ct4(1);
ict1 = 1 : lct4;
if ct4sym < ct4(end)
   ctc4 = ctc0 - ct4;
   ict1 = find(ctc4 < ct4sym);
   lct1 = size(ict1,2);
%   ct0 = ctc0 - ct4sym;
   ct0 = ct4(ict1(1));
% map the region ctc4 >= ct4sym onto ict1 by mirroring it wrt the point: ctc4 = ct4sym
   ict2 = find(ctc4 < ct4sym  &  ctc4 > 2*ct4sym - ctc0 + ct4step);
% values of ct4 to be used to calculate (by interpolation) the mirrored source values:
   ct44 = 2*(ctc0 - ct4sym) - ct4(ict2);
% define a (small) region (index icts) where the source changes rapidely and should be splined for reconstruction:
% ct4(icts(1..end)) is projected onto ct44s(end..1) :notice that the index is inverted, the phaseshift between
% the two amounts to : 2*rem(ctc0 - ct4sym - ct4(1), ct4step).
% The same phaseshift and index inversion exists also between ct4(ict1(1..end) and ct44(end..1)
% icts and ict3 have the same size and take attention that they can be empty for ctc0 - ctad +100 < ct4sym!!
   ctcsmin = 2*ct4sym - ctc0 + (ctad - 100);
   ctcsmax = 2*ct4sym - ctc0 + (ctad + 100);
   icts = find((ct4 < ctad + 100)  &  (ct4 >= ctad - 100)  &  (ctc0 - ct4 >= ct4sym));
   ict3 = find(ctc4(ict2) < ctcsmax  &  ctc4(ict2) >= ctcsmin);
   if ~isempty(ict3),    ct44s = ct44(ict3); end
end

% load the table and find the start values for the integration of the ode for u-functions
% the startfile can be obtained by running usphst.m
%  stv1 = load('usphst.dat');
stv1 = load(stfname);
lmax0 = stv1.ust.la(end);

nrsteps = 10000;   % default 10000
%lmax = 1500;
xmax = k1(lk1) * ctc0;   % =default 3000
xstep = 2*lmax0/nrsteps;
kx = 1e-12 : xstep : xmax;

% allocate the data in order to speed up the routine
kctc0 = zeros(1, lct4);
kctc = zeros(1, lct4);
jl = zeros(lct4, 1);
htable = zeros(lct4, 1);
cmbtable = zeros(1, lk1);

% calculate the curvature length times wavenumber : a parameter for the ultra_spherical bessel function
kb = sqrt(abs(K)) ./k1;

%#########
%kb(1)
%kb(end)
if all(kb == 0), kbzero = 1; x10 = 1e-12 * ones(1, lk1);
else  % not all kb are zero

   kbzero = 0;
   ikb = find(kb > 10);
   if ~isempty(ikb),
      kb(ikb) = 10 * ones(1, size(ikb, 2));
      message('boundary kb > 10 reached')
   end

% split the structure stv1 from the startfile : x1: start values where pbi_beta = 1e-6,
% pbd : dphi_beta/dx at the start values, all data exact within 1e-6.
% the la-values resp kb-vector should be the same as in the program usphst.m
   la1 = stv1.ust.la;
   lla1 = size(la1, 2);
   kb2 = stv1.ust.kb;
   x1v1 = stv1.ust.x1;
   pbdv1 = stv1.ust.pbd;

end %all(kb == 0)

%############################# start
for il = 1:lla

    l = la(il)
    laa = sqrt(l * (l + 1));

% make a range of wavenumber and kb starting at 1/(l+1)
    clear('k2')
%    ik2 = (l+1) : lk1;
    ik2 = find(nn1 >= l+1);
    k2 = k1(ik2);
    lk2 = size(k2, 2);
    nn2 = nn1(ik2);
    cmbtable = zeros(1,lk2);

% the range of indices in src for k2 will be ik2 : take care, do not neglect!!:
    k2min = k2(1);
    k2max = k2(lk2);
    tetatl = zeros(1, lk2);
    tetael = zeros(1, lk2);
    kb = sqrt(abs(K)) ./k2;
    odd_la = xor(rem(nn2,2), rem(l,2));

% the kb-range is made more progressive towards 1/la
    nrkbsteps = 100;
    kbmax = 1/laa;
    ekbmax = -7;
    ekbmin = log(kbmax - 1e-5)/log(10);
    ekbstep = (ekbmax - ekbmin)/nrkbsteps;
    ekb1 = ekbmin : ekbstep : ekbmax;
    kb1 = kbmax - 10 .^ekb1;
    lkb1 = size(kb1, 2);
    if l == la1(1)
        if any(kb1 ~= kb2)  error('kb vector is unequal to the start vector'); return; end
    end

% find the start values (hs) for the ode integration step as a function of la and kb (is a surface):
% the long sought magic formula : la1 is the la-vector from usphst1.m!!
    hs = 0.5 * x1v1 .* repmat((la1 .^-0.9)', 1, lkb1);
    if ~kbzero
       il1 = find(la1 == l);
       if isempty(il1)
           error('la not in the range of la1')
           break;
       else
           x10 = spline(kb1, x1v1(il1,:), kb);
           pbd0 = spline(kb1, pbdv1(il1,:), kb);
           hstart = spline(kb1, hs(il1,:), kb);
       end

% calculate the start values for the u-function : u0 = [u, du/dx] with u = r(x) * phi_beta(x) and
%  du/dx = phi_beta(x) * dr(x)/dx + r(x) * dphi_beta(x)/dx :
       r10 = sin(kb .* x10) ./kb;
       rd10 = cos(kb .* x10);
       u0 = [r10 *1e-6; r10 .* pbd0 + 1e-6 * rd10];

% define the step and set 'odeint' to a constant number of steps for all cases (kb)
       toi = 1e-4;
       hmax = 0.3;   %we take hmax = 0.3 as the default value
       maxstp = 150;   % default maxstp = 150
       x1e = ones(1, lk2) * xmax;
% solve the u-function instead of the phi_beta function : xs2 does not have a constant step
       [xs2, u2, nok, nbad, nfev] = odeintp(u0, x10, x1e, toi, hstart, 0, hmax, maxstp, @uspheq1, l, kb, kk0);

% prepare a spline of the solution for a constant step in xs (the values are found with ppval1):
       y2der = splinep(xs2, u2);
% calculate the number of constant steps, xstep in the solution xs2: ilast is the last one
%       ilast = ceil((xs2(end,:) - x10)/xstep);
% set the number of overlapping steps for the matching to the asymptotic solution (default= 40):
       noverflow = 85;
       ncsteps = 0;

    end  % kbzero

% calculate also the spherical bessel function as we may need it for kb < 1e-5
    table_bv = sphbes(l, kx);
% set some constants needed for the polarisation formula
    gl = sqrt((l + 2) * (l + 1) * l * (l - 1));

% interpolate and integrate (sum with the Simpson rule) over conformal time :
    for i = 1:lk2
%    for i = 47:47

       clear('htable')
       if kb(i) > 1e-5 * (1500/l)

% define a vector of argument values for the ultra-spherical bessel function (x1): reverse ct4
%          clear('xs', 'x1', 'x2', 'yspl2','iover')
          x1 = x10(i) : xstep : xmax;
% define some parameters for the calculation ; x1=x1sym up to the symmetry point, x1* kb = pi/2 needed for K>0 :
          x1sym = min(0.5*pi/kb(i), xmax);
% in the next 10 lines make a table over x1 (ul) of ultra-spherical bessel values:

          if  xs2(end,i) >=  x1sym

% asymptotic expansion is not needed here as we found already the complete solution: include one point xs2 > x1sym+xstep
% adding xs2step + xstep to x1sym in the find logic statement of xs2 makes sure that xs(end) > x1(ilast) > x1sym
             xs2step = max(diff(xs2(:,i)));
             ix2 = find(xs2(:,i) <= (x1sym + xs2step + xstep));
% include one point xs2 > x1sym in the interpolation of the ode solution:
%             ix2 = [ix2', (ix2(end) + 1)];
             [xs, yspl2] = ppval1(xs2(ix2,i), u2(ix2,i), y2der(ix2,i), xstep);
%             ilast = floor((x1sym - x10(i))/xstep) + 1;
             ilast = floor((x1sym - x10(i))/xstep) + 2;

          else  % all steps in odeint are used (full interpolation) : an asymptotic solution is needed

             [xs, yspl2] = ppval1(xs2(:,i), u2(:,i), y2der(:,i), xstep);
             ilast = floor((xs(end) - x10(i))/xstep) + 1 - ncsteps;
% define the overflow index (iover): default 85 (or 40) steps should at least include one zero and one maximum of u2
             iover = (ilast - noverflow) : ilast;

% find the asymptotic values and the phase correction in the overflow region:
             dphi = phdif(x1(iover), usphas(l, kb(i), x1(iover), 0, 0, kk0), yspl2(iover));

          end

% calculate the ultra-spherical bessel values
          kctc0 = ctc0 * k2(i) * ones(1, lct4);
          kctc = k2(i) * ct4;
          xctc = kctc0 - kctc;
% xctc > x1(1) is necessary as x1(1)=x10 is not the origin of ct, but the point where odeint started
          if   xs2(end,i) >=  x1sym
% debug action 12-6-2005 : interpolation over xs limits our range to xs(end) at most if it happens xs(end) < x1sym
%             ixode = find(xctc > x1(1)  &  xctc <= x1sym);
             ixode = find(xctc > x1(1)  &  xctc <= min(x1sym, xs(end)));
          else
             ixode = find(xctc > x1(1)  &  xctc <= x1(ilast));
          end
          ul = zeros(1, lct4);
          if  ~isempty(ixode)
             r1 = sin(kb(i) * xs)/kb(i);

% interpolate the ultra-spherical bessel function when we are in the k-range of the ode solution
             ul(ixode) = interpl(xs, yspl2 ./ r1, xctc(ixode));
% calculate the necessary values of the asymptotic expansion of the ultra-spherical bessel function
             if  (xs2(end,i) <  x1sym)   % asymptotic expansion needed for K>0
                ixas = find(xctc > x1(ilast)  &  xctc <= x1sym);
                ul(ixas) = usphas1(l, kb(i), xctc(ixas), dphi, 1, kk0);
             end

% if necessary adapt the source to the period of the ultra-spherical function for K >0 by mirroring:
% we should have: find(k2(i)* (ctc0 - ct4(ict1)) < x1sym) == find(xctc < x1sym) == ict1
             src41 = src4(:,ik2(i));
             src51 = src5(:,ik2(i));
             if  ct4sym < ct4(end)
% reconstruct the source by linear interpolation : not so accurate but sufficient where changes are slow
                src42 = interpl(ct4, src4(:,ik2(i))', ct44)';
                src52 = interpl(ct4, src5(:,ik2(i))', ct44)';
% for efficieny reason reconstruct the source by splines only in the region (ict3) where it changes fast
                if ~isempty(ict3)
                  src42(ict3) = spline(ct4(icts), src4(icts,ik2(i)), ct44s)';
                  src52(ict3) = spline(ct4(icts), src5(icts,ik2(i)), ct44s)';
                end
                if odd_la(i) == 1
                   src41(ict2) = src41(ict2) + src42;
                   src51(ict2) = src51(ict2) + src52;
                else
                   src41(ict2) = src41(ict2) - src42;
                   src51(ict2) = src51(ict2) - src52;
                end
             end

debug = 0;
if debug == 1
ct4step
nn2(i)
ctc0-ct4sym
ul(ict1(1))
src41(ict1(1))
ct0
odd_la(i)
xs2(end,i)
ct4(ict1(1))
ct4(end) - ct4sym
ct4(end) - x10(i)/k2(i)

figure(1)
plot(ct4(ict1), src41(ict1), ct4, ul, 'k--', ct4, src4(:,ik2(i)), 'b--', ct4(ict2(ict3)), src42(ict3),'r')
%plot(ct4(ict1), src41(ict1), ct4, ul, 'k--', ct4, src4(:,ik2(i)), 'b--')
%axis([6500, 7100, -0.05, 0.05])
end

% integrate the source function over ct4, following formula (40) of Zaldarriaga et al.(1998)
             htable = src41(ict1) .* ul(ict1)';
             tetatl(i) = simpsint(ct0, ct4(end), htable) + (ct0 - ctc0 + ct4sym) * src41(ict1(1)) * ul(ict1(1));
             htable = src51(ict1) .* ul(ict1)';
             tetael(i) = simpsint(ct0, ct4(end), htable) + (ct0 - ctc0 + ct4sym) * src51(ict1(1)) * ul(ict1(1));
          else
             tetatl(i) = 0;
             tetael(i) = 0;
          end

       else  %if kb(i) <= 1e-5 *(1500/l)

% interpolate the spherical bessel function
          kctc0 = ctc0 * k2(i) * ones(1, lct4);
          kctc = k2(i) * ct4;
          xctc = kctc0 - kctc;
          jl = interpl(kx, table_bv, xctc)';
% integrate the source function following (12) and (13) of Seljak resp (18) of Zaldarriaga and Seljak
% the integration interval is limited to ct4 where the source <> 0
          htable = src4(:,ik2(i)) .* jl;
          tetatl(i) = simpsint(ct4(1), ct4(end), htable);
          htable = src5(:,ik2(i)) .* jl;
          tetael(i) = simpsint(ct4(1), ct4(end), htable);

       end  %if kb(i)

   end  % forloop k2

%ta1 = tetatl;
f2 = 1;
%if (K ~= 0),  f2 = coth(pi*k2/sqrt(abs(K))); end
factor = sqrt((k2 .^2 - 4*K) ./ (k2 .^2 - K));
%factor = 1;
% perform the integration over k2min-k2max, following (9) of Seljak resp (19) of Zaldarriaga and Seljak
   cmbtable = f2 .* factor .* (tetatl .^2) .* k2 ./ (k2 .^2 - K) .^ (1.5 - 0.5*nprimtilt);
ta1 = cmbtable;
%    cct.ctt(il) = l*(l + 1) * sqrt(K) * sum(cmbtable');
    cct.ctt(il) = l*(l + 1) * k1step * sum(cmbtable');
%   cct.ctt(il) = l*(l + 1) * simpsint(k2min, k2max, cmbtable')
   cmbtable = factor .* (tetael .^2) .* k2 ./ (k2 .^2 - K) .^ (1.5 - 0.5*nprimtilt);
   cct.cee(il) = gl^2 * l*(l + 1) * k1step * sum(cmbtable');
   cmbtable = factor .* (tetael .* tetatl) .* k2 ./ (k2 .^2 - K) .^ (1.5 - 0.5*nprimtilt);
   cct.cte(il) = l*(l + 1) * gl * k1step * sum(cmbtable');
end  % forloop la


text:
function [pnt,jac] = nrbdeval(nurbs, dnurbs, tt)
% Evaluation of the derivative NURBS curve or surface.
%
%     [pnt, jac] = nrbdeval(crv, dcrv, tt)
%     [pnt, jac] = nrbdeval(srf, dsrf, {tu tv})
%
% INPUTS:
%
%   crv    - original NURBS curve.
%
%   srf    - original NUBRS surface
%
%   dcrv   - NURBS derivative represention of crv
%
%   dsrf   - NURBS derivative represention of surface
%
%   tt     - parametric evaluation points
%            If the nurbs is a surface then tt is a cell
%            {tu, tv} are the parametric coordinates
%
%   pnt  - evaluated points.
%   jac  - evaluated first derivatives (Jacobian).
%
% Examples:
%
%   // Determine the first derivatives a NURBS curve at 9 points for 0.0 to
%   // 1.0
%   tt = linspace(0.0, 1.0, 9);
%   dcrv = nrbderiv(crv);
%   [pnts,jac] = nrbdeval(crv, dcrv, tt);

%  D.M. Spink
%  Copyright (c) 2000.

if ~isstruct(nurbs)
  error('NURBS representation is not structure!');
end

if ~strcmp(nurbs.form,'B-NURBS')
  error('Not a recognised NURBS representation');
end

[cp,cw] = nrbeval(nurbs, tt);

if iscell(nurbs.knots)

  % NURBS structure represents a surface
  temp = cw(ones(3,1),:,:);
  pnt = cp./temp;

  [cup,cuw] = nrbeval(dnurbs{1}, tt);
  tempu = cuw(ones(3,1),:,:);
  jac{1} = (cup-tempu.*pnt)./temp;

  [cvp,cvw] = nrbeval(dnurbs{2}, tt);
  tempv = cvw(ones(3,1),:,:);
  jac{2} = (cvp-tempv.*pnt)./temp;

else

  % NURBS is a curve
  temp = cw(ones(3,1),:);
  pnt = cp./temp;

  % first derivative
  [cup,cuw] = nrbeval(dnurbs,tt);
  temp1 = cuw(ones(3,1),:);
  jac = (cup-temp1.*pnt)./temp;

end


text:
function pass = test_eulerTricomi( prefs )
% Check that Euler--Tricomi equation is working.
% Alex Townsend, August 2013.

% First, check that all the particular solutions are being
% calculated correctly.

if ( nargin < 1 )
    prefs = chebfunpref();
end
tol = 100*prefs.techPrefs.chebfuneps;

exact = chebfun2(@(x,y) 1+0*x);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(1) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) y);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(2) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) x);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(3) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) x.*y);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(4) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) 3*y.^2+x.^3);
%exact = chebfun2(@(x,y) y.^2-x.^2);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0;
pass(5) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) 3*x.^2+y.^3);
% exact = chebfun2(@(x,y) 3*x.^2+y.^3);
N = chebop2(@(x,y,u) y.*diff(u,2,2) - diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(6) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) 3*x.^2+y.^3);
% exact = chebfun2(@(x,y) 3*x.^2+y.^3);
N = chebop2(@(x,y,u) -y.*diff(u,2,2) + diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(7) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) y.^3+x.^3.*y);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(8) = ( norm( u - exact ) < tol );

exact = chebfun2(@(x,y) 6*x.*y.^2+x.^4);
N = chebop2(@(x,y,u) diff(u,2,2) - x.*diff(u,2,1));
N.lbc = exact(-1,:); N.rbc = exact(1,:);
N.ubc = exact(:,1); N.dbc = exact(:,-1);
u = N \ 0 ;
pass(9) = ( norm( u - exact ) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function s = minus(f,g)
%-   Subtraction of two SINGFUN objects.
%   F - G subtracts G from F, where F and G are SINGFUN objects or scalars.
%
% See also PLUS.

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

% Call PLUS():
s = plus(f, -g);

end


text:
function pf = loadPoleFigure_dubna(fname,varargin)
% load dubna cnv file
%
% Syntax
%   pf = loadPoleFigure_dubna(fname)
%
% Input
%  fname - file name
%
% Options
%  STANDARD_GRID - use reguar 5 degree grid
%  DUBNA_GRID    - use original rotated Dubna grid
%
% Output
%  pf    - @PoleFigure
%
% See also
% loadPoleFigure dubna_demo ImportPoleFigureData

% ensure right extension
[pathstr, name, ext] = fileparts(fname); %#ok<ASGLU>
if ~any(strcmpi(ext,{'.cnv','.cns'}))
  interfaceError(fname);
end

% load data
try
  d = load(fname);
catch
  interfaceError(fname);
end

% ensure sufficients data
if ~isa(d,'double') || mod(numel(d),72)~=0 ||...
    numel(d)<3*72 || numel(d)>40*72
  interfaceError(fname);
end

d = reshape(d.',72,[]);
h = string2Miller(fname);

switch lower(ext)
  case '.cnv'
    r  = DubnaGrid(size(d,2));
  case '.cns'
    r = regularS2Grid('points',size(d),'maxtheta',5*degree*(size(d,2)-1),'antipodal');
end

pf = PoleFigure(h,r,d,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
classdef nnnonorm < nntest

  properties (TestParameter)
    rows = {2 4}
    cols = {2 5}
    numDims = {4 8}
    batchSize = {1 2 3}
  end

  methods (Test)
    function basic(test, rows, cols, numDims, batchSize)
      H = rows ;
      W = cols ;
      C = numDims ;
      bs = batchSize ;
      x = test.randn(H, W, C, bs) ;
      g = test.randn(1, 1, C, 1) / test.range ;
      b = test.randn(1, 1, C, 1) / test.range ;

      y = vl_nnnonorm(x, g, b) ;
      dzdy = test.randn(size(y)) ;
      [dzdx,dzdg,dzdb] = vl_nnnonorm(x, g, b, dzdy) ;
      test.der(@(x) vl_nnnonorm(x, g, b), x, dzdy, dzdx, test.range * 1e-3) ;
      test.der(@(g) vl_nnnonorm(x, g, b), g, dzdy, dzdg, test.range * 1e-3) ;
      test.der(@(b) vl_nnnonorm(x, g, b), b, dzdy, dzdb, test.range * 1e-3) ;
    end
  end
end

repo: MFT


text:
%IM_INVERT Fixed mapping for images inversion
%
% A = IM_INVERT(A)
% A = A*IM_INVERT
%
% DESCRIPTION
% Inverts image A by subtracting it from its maximum. Note that binary
% images can better be inverted by A = ~A. In that case also A = 1-A can
% be done.
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES

% Copyright: D. de Ridder, R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_invert(a)

  if nargin < 1 | isempty(a)
    b = prmapping(mfilename,'fixed');
    b = setname(b,'Image inverse');
  elseif isa(a,'prdataset') % allows datafiles too
    isobjim(a);
    b = filtim(a,mfilename);
  elseif isa(a,'double') | isa(a,'dip_image') % here we have a single image
		if isa(a,'dip_image'), a = double(a); end
		b = max(max(max(a)))-a;
  else
    error('Illegal input')
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% Author: Julio Zaragoza
% The University of Adelaide, Adelaide, South Australia

% Example: using the CACC function for discretizing continuos data.
%
% Note that such data appeared in the paper:
%  'A Discretization Algorithm based on Class-Attribute Contingency Coefficient' (CACC),
%  by Sheng-Jung Tsai, Shien-I Lee and Wei-Pang Yang, Information Sciences, Elsevier, 2008
%  and my results are the same as those reported in the paper.

contdata = [3,1;    % here 1 (in the second column) stands for class 'Care'
            5,1;    % 2 stands for class 'Edu' and 3 stands for class 'Work'.
            6,1;
            15,2;
            17,2;
            21,2;
            35,3;
            45,3;
            46,3;
            51,2;
            56,2;
            57,2;
            66,1;
            70,1;
            71,1];
[ discdata,discvalues,discscheme ] = cacc(contdata)


text:
function [ w, success ] = quad_form( x, Q, v, w )

%QUAD_FORM quadratic form.
%   QUAD_FORM(x,Q) is real(x'*Q*x) = x'*((Q+Q')/2)*x.
%   QUAD_FORM(x,Q,v,w) is real(x'*(Q*x+v)+w).
%
%   x must be a row or column vector, and Q must either be a scalar or
%   a square matrix with the same number of rows as x. If supplied, v must
%   be a scalar or a vector of the same size as x, and w must be a scalar.
%   If x is a row vector, then x (and v) are Hermitian transpoed before
%   evaluation.
%
%   NOTE: The use of QUAD_FORM can often be replaced by a call to NORM. For
%   example, if Q is positive definite, then the constraint
%       quad_form(x,Q) <= 1
%   is equivalent to
%       norm(sqrtm(Q)*x) <= 1
%   Generally speaking, the NORM version will be more reliable and more
%   accurate, so we encourage you to make similar conversions whenever
%   possible. We *strongly* discourage the QP-era practice of converting
%   NORM expressions into quadratic forms.
%
%   Disciplined convex programming information:
%       QUAD_FORM(x,Q,v,w) is neither convex nor concave in x and (Q,v)
%       jointly, so at least one of the two must be constant.
%
%       If (Q,v) is constant, then QUAD_FORM is convex if Q is positive
%       semidefinite, and concave if Q is negative semidefinite. An error
%       is generated if Q is indefinite (unless x is also constant).
%       QUAD_FORM is nonmonotonic in x, so x must be affine.
%
%       If x is constant, then QUAD_FORM is affine in Q, v, and w. The
%       signs of x will govern whether the elements of Q, v, and w may
%       be convex, concave, or affine.

tol = 16 * eps;
tolLDL = 4 * eps;
if nargin < 4,
    w = 0;
    if nargin < 3,
        v = 0;
    end
end

sx = size( x );
if length( sx ) ~= 2 || all( sx ~= 1 ),
    cvx_throw( 'The first argument must be a vector.' );
end
nx = prod( sx );

sQ = size( Q );
if length( sQ ) ~= 2 || sQ( 1 ) ~= sQ( 2 ),
    cvx_throw( 'The second argument must be a scalar or a square matrix.' );
elseif sQ( 1 ) ~= 1 && sQ( 1 ) ~= nx,
    cvx_throw( 'The size of Q is incompatible with the size of x.' );
end

if nargin < 3,
    v = sparse( sx(1), sx(2) );
elseif ~isequal( size( v ), sx ),
    cvx_throw( 'The size of v is incompatible with the size of x.' );
end

if nargin < 4,
    w = 0;
elseif numel(w) ~= 1 || ~isreal(w),
    cvx_throw( 'The fourth argument must be a real scalar.' );
end

if sx(1) ~= nx,
    x = x';
    v = v';
end

success = true;
cvx_optval = [];
if cvx_isconstant( x ),
    x = cvx_constant( x );
    cvx_optval = real( x' * Q * x ) + sum( real( v' * x ) ) + w;
    return
elseif ~cvx_isconstant( Q ) || ~cvx_isconstant( v ),
    cvx_throw( 'Either x or (Q,v) must be constant.' );
elseif ~cvx_isaffine( x ),
    cvx_throw( 'First argument must be affine.' );
end

Q = cvx_constant( Q );
Q = 0.5 * ( Q + Q' );
v = cvx_constant( v );
dQ = diag( Q );

if nnz( Q ) == 0,
    cvx_optval = real( v' * x ) + w;
    return
end

if sQ( 1 ) == 1,
    x = x + 0.5 * ( v / Q );
    w = w - 0.25 * ( v' * v ) / Q;
    if ~isreal( x ), x = abs( x ); end
    w = real( Q ) * sum_square( x ) + w;
    return
end

while true,

    %
    % Remove zero rows and columns from Q. If a diagonal element of Q is
    % zero but there are elements on that row or column that are not,
    % then we know that neither Q nor -Q is PSD.
    %

    dQ = diag( Q );
    trQ = sum( dQ );
    if ~all( dQ ),
        nnzQ = nnz( Q );
        tt = dQ ~= 0;
        Q = Q( tt, tt );
        if nnz( Q ) ~= nnzQ,
            success = false;
            break
        end
        dQ = dQ( tt );
        w = w + real( v( ~tt, : )' * cvx_subsref( x, ~tt, ':' ) );
        v = v( tt, : );
        x = cvx_subsref( x, tt, ':' );
        sx = length( x );
    end

    %
    % Determine the sign of the elements of Q. If they are not all of
    % the same sign, then neither Q nor -Q is PSD. Note that trQ has
    % preserved the sign of our quadratic form, so setting Q=-Q here
    % in the concave case does not cause a problem.
    %

    dQ = dQ > 0;
    if ~all( dQ ),
        if any( dQ ),
            success = false;
        else
            Q = -Q;
        end
    end

    %
    % We've had to modify this portion of the code because MATLAB has
    % removed support for the CHOLINC function.
    %
    % First, try a Cholesky. If it successfully completes its
    % factorization without fail, we accept it without question. If
    % it terminates early, we perform a numerical test to see if the
    % result still approximates the square root to good precision.
    %
    % If Cholesky fails, then we assume the matrix is either rank
    % deficient or indefinite. For sparse matrices, we perform an LDL
    % factorization, and remove the contributions of any 2x2 blocks,
    % negative 1x1 blocks, and near-zero 1x1 blocks on the diagonal.
    % If there are no such blocks, we accept it without question; if
    % so, we perform the same numerical test. If the test fails, we
    % assume, for sparse matrices, at least, that the matrix is
    % indefinite.
    %
    % If the matrix is dense, our final test is an eigenvalue
    % decomposition, the most expensive but the most accurate.
    %

    spQ = nnz(Q) <= 0.1 * nx * nx;
    if spQ,
        Q = sparse( Q );
        [ R, p, prm ] = chol( Q, 'upper', 'vector' );
        if any( diff(prm) ~= 1 ),
            R( :, prm ) = R; %#ok
        end
    else
        Q = full( Q );
        [ R, p ] = chol( Q, 'upper' );
        if p > 1,
            R = [ R , R' \ Q(1:p-1,p:end) ]; %#ok
        end
    end
    valid = p == 0;
    if ~valid,
        tolQ = tol * norm( Q, 'fro' );
        if p > 1,
            valid = norm( Q - R' * R, 'fro' ) < tolQ;
        end
    end
    if ~valid && spQ,
        [ R, DD, prm ] = ldl( sparse( Q ), 'upper', 'vector' );
        if nnz( R ) > max( sx, 0.1 * sx * ( sx + 1 ) / 2 ), spQ = false; end
        tt = diag(DD,1) == 0;
        tt = [ tt ; true ] & [ true ; tt ] & ( diag(DD) > tolLDL * trQ );
        DD = diag(DD);
        R  = bsxfun( @times, sqrt( DD(tt,:) ), R(tt,:) );
        if any( diff(prm) ~= 1 ), R( :, prm ) = R; end
        valid = all( tt ) || norm( Q - R' * R, 'fro' ) < tolQ;
    end
    if ~valid && ~spQ,
        [ V, D ] = eig( full( Q ) );
        if nnz( V ) <= max( length(V), 0.1 * numel(V) ), V = sparse(V); end
        D = diag( D );
        if any( D(2:end) < D(1:end-1) ),
            [D,ndxs] = sort(D);
            V = V(:,ndxs);
        end
        valid = D(1) > -tol * D(end);
        if valid,
            nzero = nnz( cumsum(D) < tol * abs(trQ) );
            V = V(:,nzero+1:end);
            D = sqrt(D(nzero+1:end));
            R = diag(sparse(D)) * V';
        end
    end
    if ~valid,
        success = false;
        break;
    end

    %
    % Scale so that the mean eigenvalue of (1/alpha)*R'*R is one.
    % Hopefully this will minimize scaling issues.
    %

    alpha = trQ / size(R,1);
    w = w + alpha * sum_square_abs( ( R * x ) / sqrt(alpha) ) + real( v' * x );
    break;

end

if ~success && nargout == 1,
    cvx_throw( 'Disciplined convex programming error:\n    The second argument must be positive or negative semidefinite.' );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
% When you have a light and object described by a mesh, you want to know
% the shadow of the object.
% This function calculates a mesh named "shadowvolume" (see wikipedia) from
% a triangulated object mesh and a light.
%
% [SVvertices,SVfaces]=patchshadowvolume(OBJvertices,OBJfaces,L);
%
% Inputs,
%    OBJvertices, OBJfaces : The triangulated patch vertices and faces
%                       of the object causing a shadow
%    L: The light must be a 1x4 array with x,y,z,d,
%		 with d=0 for parallel light, then x,y,z is the light direction
%		 and d=1 for point light, then x,y,z is the light position
%
% Outputs,
%    SVvertices, SVfaces : The triangulated shadow volume
%
% Function is written by D.Kroon University of Twente (March 2010)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft