text:
function rgb = hsi2rgb(hsi)
%HSI2RGB Converts an HSI image to RGB.
%   RGB = HSI2RGB(HSI) converts an HSI image to RGB, where HSI is
%   assumed to be of class double with:
%     HSI(:,:,1) = hue image, assumed to be in the range [0,1] by
%                    having been divided by 2*pi.
%     HSI(:,:,2) = saturation image, in the range [0,1].
%     HSI(:,:,3) = intensity image, in the range [0,1].
%
%   The components of the output image are:
%     RGB(:,:,1) = red.
%     RGB(:,:,2) = green.
%     RGB(:,:,3) = blue.
%
%   All the conversion equations are explained in Chapter 7 of DIPUM3E.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Extract the individual HSI component images.
H = hsi(:,:,1) * 2 * pi;
S = hsi(:,:,2);
I = hsi(:,:,3);

% Implement the conversion equations.
R = zeros(size(hsi,1),size(hsi,2));
G = zeros(size(hsi,1),size(hsi,2));
B = zeros(size(hsi,1),size(hsi,2));

% RG sector (0 <= H < 2*pi/3).
idx = find( (0 <= H) & (H < 2*pi/3));
B(idx) = I(idx).*(1 - S(idx));
R(idx) = I(idx).*(1 + S(idx).*cos(H(idx))./...
                                          cos(pi/3 - H(idx)));
G(idx) = 3*I(idx) - (R(idx) + B(idx));

% BG sector (2*pi/3 <= H < 4*pi/3).
idx = find( (2*pi/3 <= H) & (H < 4*pi/3) );
R(idx) = I(idx).*(1 - S(idx));
G(idx) = I(idx).*(1 + S(idx).*cos(H(idx) - 2*pi/3)./...
                    cos(pi - H(idx)));
B(idx) = 3*I(idx) - (R(idx) + G(idx));

% BR sector.
idx = find( (4*pi/3 <= H) & (H <= 2*pi));
G(idx) = I(idx).*(1 - S(idx));
B(idx) = I(idx).*(1 + S(idx).*cos(H(idx) - 4*pi/3)./...
                                           cos(5*pi/3 - H(idx)));
R(idx) = 3*I(idx) - (G(idx) + B(idx));

% Combine all three results into an RGB image. Clip to [0,1] to
% compensate for floating-point arithmetic rounding effects.
rgb = cat(3,R,G,B);
rgb = max(min(rgb,1),0);

author: dipum

repo: dipum-toolbox

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


text:
function [prob,sol,fmin] = qp_prob(varargin)
%QP_PROB  Return an OPTI QP
%
%   prob = qp_prob(no) return a pre-built optiprob of a saved QP.
%
%   [prob,sol,fmin] = qp_prob(no) returns the optimum solution and
%   function eval at the optimum
%
%   no = qp_prob() returns the number of problems available for testing.

%   (C) 2011 Jonathan Currie (I2C2)

%Check if just returning no problems
if(nargin < 1)
    prob = 3; sol = []; fmin = [];
    return;
else
    no = varargin{1};
end

%Big switch yard
switch(no)
    case 1
        H = eye(3);
        f = -[2 3 1]';
        A = [1 1 1;3 -2 -3; 1 -3 2];
        b = [1;1;1];
        prob = optiprob('H',H,'f',f,'ineq',A,b);
        sol = [1/3;4/3;-2/3];
        fmin = -2.83333333301227;

    case 2
        H = [1 -1; -1 2];
        f = -[2 6]';
        A = [1 1; -1 2; 2 1];
        b = [2; 2; 3];
        lb = [0;0];
        prob = optiprob('H',H,'f',f,'ineq',A,b,'lb',lb);
        sol = [1/3;4/3];
        fmin = -8.22222220552525;

    case 3
        H = [1 -1; -1 2];
        f = -[2 6]';
        A = [1 1; -1 2; 2 1];
        b = [2; 2; 3];
        Aeq = [1 1.5];
        beq = 2;
        lb = [0;0];
        ub = [10;10];
        prob = optiprob('H',H,'f',f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub);
        sol = [0.34482763667623;1.10344824221585];
        fmin = -6.41379310344827;

    otherwise
        error('Problem not available or not implemented yet');
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function appendDataToHdf5(filename, datasetName, inputData, varargin)
	% Appends 3D data to existing dataset in hdf5 file, assumes x,y are the same, only appending z.
	% Biafra Ahanonu
	% started: 2014.01.07
	% inputs
		%
	% outputs
		%

	ciapkg.hdf5.appendDataToHdf5(filename, datasetName, inputData,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% FEATURE_WRAPPER Wrapper for feature functions
%
% Usage
%    feature = FEATURE_WRAPPER(x, objects, feature_fun, options)
%
% Input
%    x (numeric): The file data.
%    object (struct): The objects contained in the data.
%    feature_fun (function handle): The real feature function handle, takes as
%       input one signal (or multiple, arranged as columns of a matrix) and
%       outputs the corresponding feature vectors. These are arranged with
%       feature dimension along the first axis, time/space along second and
%       third axes and signal index along the fourth axis (if more than one
%       signal are input).
%    options (struct): Various parameters for the wrapper, such as:
%       options.input_sz (int): The size of the input vectors to be given to
%          feature_fun. If empty, takes the rectangle specified by the objects
%          structure, otherwise takes the rectangle of size input_sz centered
%          on the object (default empty).
%       options.output_sz (int): The desired size of the data covered by the
%          feature vector. If empty, keeps the output from feature_fun, other-
%          wise extracts the time/space rectangle of size output_sz centered
%          on the original data, taking into account any subsampling by
%          feature_fun (default empty).
%       options.obj_normalize (int): The normalization of the input vectors
%          before being given to feature_fun. Can be empty, 1, 2, or Inf
%          (default []).
%       options.collapse (boolean): If true, collapses the time/space dimen-
%          sion into one vector, otherwise leaves this dimension intact
%          (default false).
%
% Output
%    feature (numeric): An PxNxK array where P is the feature dimension, N is
%       the space/time dimension and K is the signal index, if multiple
%       objects are given as input.
%
% See also
%    PREPARE_DATABASE

function t = feature_wrapper(x,objects,fun,options)
	if nargin < 4
		options = struct();
	end

	options = fill_struct(options, 'input_sz', []);
	options = fill_struct(options, 'output_sz', []);
	options = fill_struct(options, 'obj_normalize', []);
	options = fill_struct(options, 'collapse', 0);

	if isempty(options.input_sz)
		sz=objects(1).u2-objects(1).u1+ones(size(objects(1).u1));

        if length(sz)==1
            sz=[sz 1];
        end
        buf = zeros([sz, ...
			length(objects)]);

		u1 = [objects.u1];
		u2 = [objects.u2];
	else
		if length(options.input_sz) == 1
			options.input_sz = [options.input_sz 1];
		end

		buf = zeros([options.input_sz,length(objects)]);

		if options.input_sz(2)==1
			u1 = round(([objects.u1]+[objects.u2]+1)/2-options.input_sz(1)/2);
			u2 = u1+options.input_sz(1)-1;
		else
			u1 = round(([objects.u1]+[objects.u2]+1)/2-options.input_sz/2);
			u2 = u1+options.input_sz-1;
		end
	end

	% extract objects with bounding box
	for l = 1:length(objects)
		n_dim = sum(size(x)>1);
		switch n_dim
			case 1
				ind = max(u1(l),1):min(u2(l),length(x));
				buf(:,1,l) = [zeros(max(0,1-u1(l)),1); x(ind); zeros(max(0,u2(l)-length(x)),1)];
			case 2
				buf = x;
				% TODO : extract bounding box
		end
	end

	% normalize
	if ~isempty(options.obj_normalize)
		if options.obj_normalize == 1
			n = sum(abs(buf),1);
		elseif options.obj_normalize == 2
			n = sqrt(sum(abs(buf).^2,1));
		elseif options.obj_normalize == Inf
			n = max(abs(buf),[],1);
		end

		buf = bsxfun(@times,buf,1./n);
	end

	t = fun(buf);

	if ~isempty(options.output_sz)
		if length(options.output_sz) == 1
			options.output_sz = [options.output_sz 1];
		end

		N = [size(t,2) size(t,3)];
		extent = floor(options.output_sz./(2*options.input_sz).*N);

		if options.input_sz(2) > 1
			t = t(:,N(1)/2+1-extent(1):N(1)/2+1+extent(1), ...
				N(2)/2+1-extent(2):N(2)/2+1+extent(2),:);
		else
			t = t(:,N(1)/2+1-extent(1):N(1)/2+1+extent(1),1,:);
		end
	end

	t = reshape(t,[size(t,1) size(t,2)*size(t,3) size(t,4)]);

	if options.collapse
		t = reshape(t,[size(t,1)*size(t,2) 1 size(t,3)]);
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function fix_lines(fname, fname2)
%FIX_LINES  Improves the line style of eps files generated by print
%
% Examples:
%   fix_lines fname
%   fix_lines fname fname2
%
% This function improves the style of lines in eps files generated by
% MATLAB's print function, making them more similar to those seen on
% screen. Grid lines are also changed from a dashed style to a dotted
% style, for greater differentiation from dashed lines.
%
% The function also places embedded fonts after the postscript header, in
% versions of MATLAB which place the fonts first (R2006b and earlier), in
% order to allow programs such as Ghostscript to find the bounding box
% information.
%
% IN:
%   fname - Name or path of source eps file.
%   fname2 - Name or path of destination eps file. Default: same as fname.

% Copyright: (C) Oliver Woodford, 2008-2010

% The idea of editing the EPS file to change line styles comes from Jiro
% Doke's FIXPSLINESTYLE (fex id: 17928)
% The idea of changing dash length with line width came from comments on
% fex id: 5743, but the implementation is mine :)

% Thank you to Sylvain Favrot for bringing the embedded font/bounding box
% interaction in older versions of MATLAB to my attention.
% Thank you to D Ko for bringing an error with eps files with tiff previews
% to my attention.
% Thank you to Laurence K for suggesting the check to see if the file was
% opened.

% Read in the file
fh = fopen(fname, 'r');
if fh == -1
    error('File %s not found.', fname);
end
try
    fstrm = fread(fh, '*char')';
catch ex
    fclose(fh);
    rethrow(ex);
end
fclose(fh);

% Move any embedded fonts after the postscript header
if strcmp(fstrm(1:15), '%!PS-AdobeFont-')
    % Find the start and end of the header
    ind = regexp(fstrm, '[\n\r]%!PS-Adobe-');
    [ind2 ind2] = regexp(fstrm, '[\n\r]%%EndComments[\n\r]+');
    % Put the header first
    if ~isempty(ind) && ~isempty(ind2) && ind(1) < ind2(1)
        fstrm = fstrm([ind(1)+1:ind2(1) 1:ind(1) ind2(1)+1:end]);
    end
end

% Make sure all line width commands come before the line style definitions,
% so that dash lengths can be based on the correct widths
% Find all line style sections
ind = [regexp(fstrm, '[\n\r]SO[\n\r]'),... % This needs to be here even though it doesn't have dots/dashes!
       regexp(fstrm, '[\n\r]DO[\n\r]'),...
       regexp(fstrm, '[\n\r]DA[\n\r]'),...
       regexp(fstrm, '[\n\r]DD[\n\r]')];
ind = sort(ind);
% Find line width commands
[ind2 ind3] = regexp(fstrm, '[\n\r]\d* w[\n\r]');
% Go through each line style section and swap with any line width commands
% near by
b = 1;
m = numel(ind);
n = numel(ind2);
for a = 1:m
    % Go forwards width commands until we pass the current line style
    while b <= n && ind2(b) < ind(a)
        b = b + 1;
    end
    if b > n
        % No more width commands
        break;
    end
    % Check we haven't gone past another line style (including SO!)
    if a < m && ind2(b) > ind(a+1)
        continue;
    end
    % Are the commands close enough to be confident we can swap them?
    if (ind2(b) - ind(a)) > 8
        continue;
    end
    % Move the line style command below the line width command
    fstrm(ind(a)+1:ind3(b)) = [fstrm(ind(a)+4:ind3(b)) fstrm(ind(a)+1:ind(a)+3)];
    b = b + 1;
end

% Find any grid line definitions and change to GR format
% Find the DO sections again as they may have moved
ind = int32(regexp(fstrm, '[\n\r]DO[\n\r]'));
if ~isempty(ind)
    % Find all occurrences of what are believed to be axes and grid lines
    ind2 = int32(regexp(fstrm, '[\n\r] *\d* *\d* *mt *\d* *\d* *L[\n\r]'));
    if ~isempty(ind2)
        % Now see which DO sections come just before axes and grid lines
        ind2 = repmat(ind2', [1 numel(ind)]) - repmat(ind, [numel(ind2) 1]);
        ind2 = any(ind2 > 0 & ind2 < 12); % 12 chars seems about right
        ind = ind(ind2);
        % Change any regions we believe to be grid lines to GR
        fstrm(ind+1) = 'G';
        fstrm(ind+2) = 'R';
    end
end

% Isolate line style definition section
first_sec = strfind(fstrm, '% line types:');
[second_sec remaining] = strtok(fstrm(first_sec+1:end), '/');
[remaining remaining] = strtok(remaining, '%');

% Define the new styles, including the new GR format
% Dot and dash lengths have two parts: a constant amount plus a line width
% variable amount. The constant amount comes after dpi2point, and the
% variable amount comes after currentlinewidth. If you want to change
% dot/dash lengths for a one particular line style only, edit the numbers
% in the /DO (dotted lines), /DA (dashed lines), /DD (dot dash lines) and
% /GR (grid lines) lines for the style you want to change.
new_style = {'/dom { dpi2point 1 currentlinewidth 0.08 mul add mul mul } bdef',... % Dot length macro based on line width
             '/dam { dpi2point 2 currentlinewidth 0.04 mul add mul mul } bdef',... % Dash length macro based on line width
             '/SO { [] 0 setdash 0 setlinecap } bdef',... % Solid lines
             '/DO { [1 dom 1.2 dom] 0 setdash 0 setlinecap } bdef',... % Dotted lines
             '/DA { [4 dam 1.5 dam] 0 setdash 0 setlinecap } bdef',... % Dashed lines
             '/DD { [1 dom 1.2 dom 4 dam 1.2 dom] 0 setdash 0 setlinecap } bdef',... % Dot dash lines
             '/GR { [0 dpi2point mul 4 dpi2point mul] 0 setdash 1 setlinecap } bdef'}; % Grid lines - dot spacing remains constant

if nargin < 2
    % Overwrite the input file
    fname2 = fname;
end

% Save the file with the section replaced
fh = fopen(fname2, 'w');
if fh == -1
    error('Unable to open %s for writing.', fname2);
end
try
    fwrite(fh, fstrm(1:first_sec), 'char*1');
    fwrite(fh, second_sec, 'char*1');
    fprintf(fh, '%s\r', new_style{:});
    fwrite(fh, remaining, 'char*1');
catch ex
    fclose(fh);
    rethrow(ex);
end
fclose(fh);
return

author: Cloud-CV

repo: object-proposals


text:
function x = inv_digamma(y,niter)
% INV_DIGAMMA    Inverse of the digamma function.
%
% inv_digamma(y) returns x such that digamma(x) = y.

% a different algorithm is provided by Paul Fackler:
% http://www.american.edu/academic.depts/cas/econ/gaussres/pdf/loggamma.src

% Newton iteration to solve digamma(x)-y = 0
x = exp(y)+1/2;
i = find(y <= -2.22);
x(i) = -1./(y(i) - digamma(1));

% never need more than 5 iterations
if nargin < 2
  niter = 5;
end
for iter = 1:niter
  x = x - (digamma(x)-y)./trigamma(x);
end
return

% test
y = -3:0.01:0.1;
x = digamma(inv_digamma(y));
max(abs(x-y))
max(abs(x-y)./inv_digamma(y))

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function cMap=wcbp(varargin)

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

cPink=[220 000 163];
cBlue=[0 0 204];
cCyan=[0 160 160];

cMap=[255*ones(1,3);...
    120*ones(1,3);...
    cCyan;...
    cBlue;...
    cPink/1.5;...
    cPink;]./255;

[cMap]=resampleColormap(cMap,n);

%%
% _*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:
%%	Author:  Tianli Yu	https://netfiles.uiuc.edu/tyu/www/tianli/
%%	Date:	Monday Aug. 01 2005  11:58 pm
%%
%%  SampleTest: This is a sample script to test the TriangulateSpherePoints.m
%%

	VertexM = rand(800, 3) - 0.5;
	VNormV = (sum(VertexM .^ 2, 2)) .^ .5;
	VertexM = VertexM ./ repmat(VNormV, 1, 3);
	FacetM = TriangulateSpherePoints(VertexM);
	figure;
	patch('vertices', VertexM, 'faces', FacetM, 'FaceColor', [1, 1, 0], 'EdgeColor', [1 .2 .2]);
	axis equal;

author: Sable

repo: mcbench-benchmarks


text:
function [tdata, validIdx] = translate_movie(data, dv, ops)

%% Parameters
[ly, lx, nFrames] = size(data);
if nargin < 3
  ops = [];
end
subpixel = getOr(ops, {'subPixel' 'SubPixel'}, 1);
if isfinite(subpixel)
  dv = round(subpixel*dv)./subpixel;
end
useGPU = getOr(ops, 'useGPU', false);

dv = permute(dv, [2 3 1]);
fy = ifftshift((-fix(ly/2):ceil(ly/2) - 1)/ly)';% freq along first dimension
fx = ifftshift((-fix(lx/2):ceil(lx/2) - 1)/lx); % freq along second dimension

if useGPU
  batchSize = 32;
  fx = gpuArray(fx);
  fy = gpuArray(fy);
  dv = gpuArray(dv);
else
  batchSize = 3;
end
tdata = zeros([ly, lx, nFrames], 'like', data);
%% Work through data in batches
nBatches = ceil(nFrames/batchSize);
for bi = 1:nBatches
  fi = (bi - 1)*batchSize + 1:min(bi*batchSize, nFrames);
  if useGPU
    batchData = gpuArray(single(data(:,:,fi)));
  else
    batchData = single(data(:,:,fi));
  end
  phaseShift = bsxfun(@times,...
    exp(-1j*2*pi*bsxfun(@times, fy, dv(1,:,fi))),... y rotation
    exp(-1j*2*pi*bsxfun(@times, fx, dv(2,:,fi)))); % x rotation
  tdata(:,:,fi) = gather(real(ifft2(fft2(batchData).*phaseShift)));
end
if nargout > 1
  dvMax = max(0, ceil(max(gather(dv), [], 3)));
  dvMin = min(0, floor(min(gather(dv), [], 3)));
  validX = (1 + dvMax(2)):(lx + dvMin(2));
  validY = (1 + dvMax(1)):(ly + dvMin(1));
  validIdx = {validY validX};
end

end

author: cortex-lab


text:
function [ obj ] = remove( obj, idx )
% Remove data (fMRI time series, confounds, DCM network structure, ... )
% and estimation results from HUGE object.
%
% INPUTS:
%   obj - A tapas_Huge object.
%
% OPTIONAL INPUTS:
%   idx - Only remove data of the subjects indicated in 'idx'. 'idx' must
%         be a vector containing numeric or logical array indices.
%
% OUTPUTS:
%   obj - A tapas_Huge object with data and results removed.
%
% EXAMPLES:
%   [obj] = REMOVE(obj)    Remove results and data of all subjects.
%
%   [obj] = REMOVE(obj,1:5)    Remove results and data for first 5
%                              subjects.
%
%   [obj] = REMOVE(obj,'all') is the same as [obj] = REMOVE(obj)
%
%   [obj] = REMOVE(obj,0) is the same as [obj] = REMOVE(obj)
%
% See also tapas_Huge.IMPORT
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

%% check input
if nargin < 2 || (ischar(idx) && strcmpi(idx, 'all'))
    idx = 0;
end

if obj.N == 0 || isempty(idx)
    return
end

if isscalar(idx) && idx == 0
    idx = 1:obj.N;
end

assert((all(idx >= 1) && all(rem(idx,1) == 0)) || islogical(idx),...
    'TAPAS:HUGE:invalidIndex',...
    'idx must be real positive integers or logicals');

assert(max(idx) <= obj.N, 'TAPAS:HUGE:invalidIndex',...
    'idx exceeds number of subjects');

if islogical(idx)
    assert(numel(idx) == obj.N, 'TAPAS:HUGE:invalidIndex',...
        'idx must have length %u', obj.N);
end

%% remove data
obj.data(idx)   = [];
obj.inputs(idx) = [];

if isnumeric(obj.options.nvp.startingvaluedcm) && ...
        ~isempty(obj.options.nvp.startingvaluedcm)
    obj.options.nvp.startingvaluedcm(idx,:) = [];
end

if ~isempty(obj.model)
    obj.model.d(idx, :) = [];
    obj.model.theta_c(idx, :) = [];
    obj.model.theta_h(idx, :) = [];
    obj.model.lambda(idx, :) = [];
    obj.model.pi = sum(obj.model.d)/sum(obj.model.d(:));
end

obj.N = numel(obj.data);
if ~obj.N
    obj = reset_properties(obj);
end

% reset prior, posterior and trace
obj.prior = [];
obj.posterior = [];
obj.trace = [];

end

function [ obj ] = reset_properties( obj )

tmp = obj.default_options();
obj.options.nvp.priorclustermean        = tmp.nvp.priorclustermean;
obj.options.nvp.priorclustervariance    = tmp.nvp.priorclustervariance;
obj.options.nvp.priordegree             = tmp.nvp.priordegree;
obj.options.nvp.priorvarianceratio      = tmp.nvp.priorvarianceratio;
obj.options.nvp.startingvaluedcm        = tmp.nvp.startingvaluedcm;
obj.options.nvp.startingvaluegmm        = tmp.nvp.startingvaluegmm;
obj.options.nvp.confoundsvariant        = tmp.nvp.confoundsvariant;

obj.R = 0;
obj.L = 0;
obj.M = 0;
obj.dcm = [];
obj.labels = struct();
obj.model = [];
obj.idx = struct();

end

author: translationalneuromodeling

repo: tapas


text:
function pdf = beta_pdf ( x, a, b )

%*****************************************************************************80
%
%% BETA_PDF evaluates the Beta PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = X**(A-1) * (1-X)**(B-1) / BETA(A,B).
%
%    A = B = 1 yields the Uniform distribution on [0,1].
%    A = B = 1/2 yields the Arcsin distribution.
%        B = 1 yields the power function distribution.
%    A = B -> Infinity tends to the Normal distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    0.0 <= X <= 1.0.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A,
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < 0.0 | 1.0 < x )
    pdf = 0.0;
  else
    pdf = x^( a - 1.0 ) * ( 1.0 - x )^( b - 1.0 ) / beta ( a, b );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%--------------------------------------------------------------------------
% HDL Workflow Script
% Generated with MATLAB 9.3 (R2017b) at 14:05:38 on 07/05/2018
% This script was generated using the following parameter values:
%     Filename  : '/tmp/MathWorks_tools/targeting_models/modem-qpsk/FixedPoint/demos/ADI_DMA_TT/hdlworkflow.m'
%     Overwrite : true
%     Comments  : true
%     Headers   : true
%     DUT       : 'combinedTxRx_ADIDMA/Combined TX and RX'
% To view changes after modifying the workflow, run the following command:
% >> hWC.export('DUT','combinedTxRx_ADIDMA/Combined TX and RX');
%--------------------------------------------------------------------------

%% Load the Model
load_system('combinedTxRx_ADIDMA');

%% Restore the Model to default HDL parameters
%hdlrestoreparams('combinedTxRx_ADIDMA/Combined TX and RX');

%% Model HDL Parameters
%% Set Model 'combinedTxRx_ADIDMA' HDL parameters
hdlset_param('combinedTxRx_ADIDMA', 'HDLSubsystem', 'combinedTxRx_ADIDMA/Combined TX and RX');
hdlset_param('combinedTxRx_ADIDMA', 'ReferenceDesign', 'adrv9361z7035 box lvds Base System (Vivado 2017.4)');
hdlset_param('combinedTxRx_ADIDMA', 'SynthesisTool', 'Xilinx Vivado');
hdlset_param('combinedTxRx_ADIDMA', 'SynthesisToolChipFamily', 'Zynq');
hdlset_param('combinedTxRx_ADIDMA', 'SynthesisToolDeviceName', 'xc7z035i');
hdlset_param('combinedTxRx_ADIDMA', 'SynthesisToolPackageName', 'fbg676');
hdlset_param('combinedTxRx_ADIDMA', 'SynthesisToolSpeedValue', '-2L');
hdlset_param('combinedTxRx_ADIDMA', 'TargetDirectory', 'hdl_prj/hdlsrc');
hdlset_param('combinedTxRx_ADIDMA', 'TargetLanguage', 'Verilog');
hdlset_param('combinedTxRx_ADIDMA', 'TargetPlatform', 'AnalogDevices adrv9361z7035 box lvds (modem)');
hdlset_param('combinedTxRx_ADIDMA', 'Workflow', 'IP Core Generation');

% Set SubSystem HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX', 'ProcessorFPGASynchronization', 'Free running');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/bytesIn', 'IOInterface', 'IP Data IN [0:63]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/bytesIn', 'IOInterfaceMapping', '[0:63]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/validIn', 'IOInterface', 'IP Valid Tx Data IN');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/validIn', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/FromRadioR', 'IOInterface', 'AD9361 ADC Data I0 [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/FromRadioR', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/FromRadioI', 'IOInterface', 'AD9361 ADC Data Q0 [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/FromRadioI', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/FRLoopBw', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/FRLoopBw', 'IOInterfaceMapping', 'x"100"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/EQmu', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/EQmu', 'IOInterfaceMapping', 'x"104"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Scope Select', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Scope Select', 'IOInterfaceMapping', 'x"108"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/DebugSelector', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/DebugSelector', 'IOInterfaceMapping', 'x"10C"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TxDMASelectIn', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TxDMASelectIn', 'IOInterfaceMapping', 'x"110"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/BypassEQ', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/BypassEQ', 'IOInterfaceMapping', 'x"114"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/EnableDecode', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/EnableDecode', 'IOInterfaceMapping', 'x"118"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/PDThreshold', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/PDThreshold', 'IOInterfaceMapping', 'x"11C"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TransmitToggle', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TransmitToggle', 'IOInterfaceMapping', 'x"120"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TransmitAlways', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TransmitAlways', 'IOInterfaceMapping', 'x"124"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/PacketSource', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/PacketSource', 'IOInterfaceMapping', 'x"128"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Loopback', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Loopback', 'IOInterfaceMapping', 'x"12C"');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/BypassCoding', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/BypassCoding', 'IOInterfaceMapping', 'x"130"');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/bytesOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/bytesOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/validBytesOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/validBytesOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/startPacketOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/startPacketOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/real', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/real', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/imag', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/IQ Interface Mapper/imag', 'IOInterfaceMapping', '');

% Set SubSystem HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL', 'ProcessorFPGASynchronization', 'Free running');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/re', 'IOInterface', 'Rx data I1 In [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/re', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/Im', 'IOInterface', 'Rx data Q1 In [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/Im', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/Enable', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/Enable', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/FRLoopBW', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/FRLoopBW', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/EQmu', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/EQmu', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/Scope Select', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/Scope Select', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/DebugSelector', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/DebugSelector', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/EQBypass', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/EQBypass', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/EnableDecode', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/EnableDecode', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/PDThreshold', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/PDThreshold', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/BypassCoding', 'IOInterface', 'Rx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/BypassCoding', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/bytesOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/bytesOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/validOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/validOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/sync', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/sync', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/payloadLenOut', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/payloadLenOut', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/dataRe', 'IOInterface', 'Rx data I1 Out [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/dataRe', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/dataIm', 'IOInterface', 'Rx data Q1 Out [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/dataIm', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/validIQ', 'IOInterface', 'Rx data Valid Out');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/validIQ', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/DebugSelection', 'IOInterface', 'Rx data Valid Out');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/DebugSelection', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/RxTransferComplete', 'IOInterface', 'Rx data Valid Out');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Receiver HDL/RxTransferComplete', 'IOInterfaceMapping', '[0]');

% Set SubSystem HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL', 'ProcessorFPGASynchronization', 'Free running');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/bytesIn', 'IOInterface', 'Tx data I1 In [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/bytesIn', 'IOInterfaceMapping', '[0:15]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/validIn', 'IOInterface', 'Tx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/validIn', 'IOInterfaceMapping', '[0]');

% Set Inport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/BypassEncode', 'IOInterface', 'Tx data Valid In');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/BypassEncode', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/Packet Generation/trueData', 'IOInterface', 'Tx data Valid Out');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/Packet Generation/trueData', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/real', 'IOInterface', 'Tx data I1 Out [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/real', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/imag', 'IOInterface', 'Tx data Q1 Out [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/imag', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/validOut', 'IOInterface', 'Tx data Valid Out');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/validOut', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/Need Data', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/Need Data', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/TxTransferComplete', 'IOInterface', 'Tx data Valid Out');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/Transmitter HDL/TxTransferComplete', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ToRadioR', 'IOInterface', 'AD9361 DAC Data I0 [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ToRadioR', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ToRadioI', 'IOInterface', 'AD9361 DAC Data Q0 [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ToRadioI', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ToRadioValid', 'IOInterface', 'IP Valid Tx Data OUT');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ToRadioValid', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/NeedData', 'IOInterface', 'IP Load Tx Data OUT');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/NeedData', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/BytesOutRx', 'IOInterface', 'IP Data OUT [0:63]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/BytesOutRx', 'IOInterfaceMapping', '[0:63]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/payloadLenOutRx', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/payloadLenOutRx', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ValidOutRx', 'IOInterface', 'IP Data Valid OUT');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/ValidOutRx', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/syncRx', 'IOInterface', 'ADC DMA sync');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/syncRx', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/reRx', 'IOInterface', 'IP Debug 1 OUT [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/reRx', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/imRx', 'IOInterface', 'IP Debug 2 OUT [0:15]');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/imRx', 'IOInterfaceMapping', '[0:15]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/validRx', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/validRx', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TxDMASelectOut', 'IOInterface', 'Tx Mux Sel');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TxDMASelectOut', 'IOInterfaceMapping', '[0]');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/debugSelectionAXI', 'IOInterface', 'AXI4-Lite');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/debugSelectionAXI', 'IOInterfaceMapping', 'x"140"');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TxTransferComplete', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/TxTransferComplete', 'IOInterfaceMapping', '');

% Set Outport HDL parameters
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/RxTransferComplete', 'IOInterface', 'No Interface Specified');
hdlset_param('combinedTxRx_ADIDMA/Combined TX and RX/RxTransferComplete', 'IOInterfaceMapping', '');

%% Workflow Configuration Settings
% Construct the Workflow Configuration Object with default settings
hWC = hdlcoder.WorkflowConfig('SynthesisTool','Xilinx Vivado','TargetWorkflow','IP Core Generation');

% Specify the top level project directory
hWC.ProjectFolder = 'hdl_prj';
hWC.ReferenceDesignToolVersion = '2016.4';
hWC.IgnoreToolVersionMismatch = false;

% Set Workflow tasks to run
hWC.RunTaskGenerateRTLCodeAndIPCore = true;
hWC.RunTaskCreateProject = true;
hWC.RunTaskGenerateSoftwareInterfaceModel = false;
hWC.RunTaskBuildFPGABitstream = true;
hWC.RunTaskProgramTargetDevice = false;

% Set properties related to 'RunTaskGenerateRTLCodeAndIPCore' Task
hWC.IPCoreRepository = '';
hWC.GenerateIPCoreReport = true;

% Set properties related to 'RunTaskCreateProject' Task
hWC.Objective = hdlcoder.Objective.None;
hWC.AdditionalProjectCreationTclFiles = '';
hWC.EnableIPCaching = false;

% Set properties related to 'RunTaskGenerateSoftwareInterfaceModel' Task
hWC.OperatingSystem = '';

% Set properties related to 'RunTaskBuildFPGABitstream' Task
hWC.RunExternalBuild = true;
hWC.TclFileForSynthesisBuild = hdlcoder.BuildOption.Custom;
hWC.CustomBuildTclFile = 'adi_build.tcl';

% Set properties related to 'RunTaskProgramTargetDevice' Task
hWC.ProgrammingMethod = hdlcoder.ProgrammingMethod.Download;

% Validate the Workflow Configuration Object
hWC.validate;

%% Run the workflow
hdlcoder.runWorkflow('combinedTxRx_ADIDMA/Combined TX and RX', hWC);

author: analogdevicesinc


text:
function commit = git_hash()
	[temp,commit] = system('git show --format=%H');
end


text:
function eclassifier = mcmcTrainEdgeClassifier(efeatures, adjlist, imsegs, labels)

ntrees = 20;
nnodes = 8;
ndata = 50000;

% train {ground, vertical, sky} classifier
[edata, elab] = formatData(efeatures, adjlist, labels, ndata);
mean(elab==1)
eclassifier = train_boosted_dt_2c(edata, [], elab, ntrees, nnodes, 0);

function [data, lab] = formatData(features, adjlist, labels, ndata)
% concatenate data and select ndata random datapoints

nimages = numel(features);

[tmp, nvars] = size(features{1});

% count edges
ne = 0;
for f = 1:nimages
    ne = ne + size(features{f}, 1);
end

data = zeros(ne, nvars);
lab = zeros(ne, 1);

% concatenate data
c = 0;
for f = 1:nimages
    cf = size(features{f}, 1);
    data(c+1:c+cf, :) = features{f};
    s1 = adjlist{f}(:, 1);
    s2 = adjlist{f}(:, 2);
    lab(c+1:c+cf) = (labels{f}(s1)==labels{f}(s2))*2-1;
    c = c + cf;
end

% select random ndata points
if ne > ndata
    rind = randperm(ne);
    rind = rind(1:ndata);
    data = data(rind, :);
    lab = lab(rind);
end

author: zouchuhang

repo: LayoutNet


text:
function a = gk323 ( m, n )

%*****************************************************************************80
%
%% GK323 returns the GK323 matrix.
%
%  Discussion:
%
%    This matrix is occasionally known as the "Todd" matrix.
%
%  Formula:
%
%    A(I,J) = abs ( I - J )
%
%  Example:
%
%    N = 5
%
%     0  1  2  3  4
%     1  0  1  2  3
%     2  1  0  1  2
%     3  2  1  0  1
%     4  3  2  1  0
%
%  Rectangular Properties:
%
%    A is integral: int ( A ) = A.
%
%    A is a special case of the Fiedler matrix.
%
%  Square Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    det ( A ) = (-1)^(N-1) * 2^(N-2) * ( N - 1 ).
%
%    A has a dominant positive eigenvalue, and N-1 real negative eigenvalues.
%
%    If N = 2 mod 4, then -1 is an eigenvalue, with an eigenvector
%    of the form ( 1, -1, -1, 1, 1, -1, -1, 1, ... ).
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   09 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.23,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 51,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns
%    of the matrix.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n
      a(i,j) = abs ( i - j );
    end
  end

  return
end


text:
% Function to derive the Linear Prediction residual signal
%
% Octave compatible
%
% Description
%  Function to derive the Linear Prediction residual signal
%
% Inputs
%  x               : [samples] [Nx1] Input signal
%  L               : [samples] [1x1] window length (e.g., 25ms =>  25/1000*fs)
%  shift           : [samples] [1x1] window shift (e.g., 5ms => 5/1000*fs)
%  order           : [samples] [1x1] Order of Linear Prediction
%
% Outputs
%  res             : [samples] [Mx1] Linear Prediction residual
%  LPCcoeff        : [samples] [order+1xM] Linear Prediction coefficients
%
% Example
%  [res,LPCcoeff] = lpcresidual(x,L,shift,order);
%
% Copyright (c) 2011 University of Mons, FNRS
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Thomas Drugman, TCTS Lab.
%
% $Id <info set by the versioning system> $

function [res,LPCcoeff] = lpcresidual(x,L,shift,order)

%% Initial settings
x=x(:);
shift = round(shift);
order = round(order);

start=1;
L = round(L);
stop=start+L;

% Allocate space
res=zeros(1,length(x));
LPCcoeff=zeros(order+1,round(length(x)/shift));

%% Do processing
n=1;
win = hanning(stop-start+1);
while stop<length(x)
     segment=x(start:stop);
     segment=segment.*win;

     A=lpc(segment,order);
     LPCcoeff(:,n)=A(:);

     inv=filter(A,1,segment);

     inv=inv*sqrt(sum(segment.^2)/sum(inv.^2));

     res(start:stop)=res(start:stop)+inv'; % Overlap and add

     % Increment
     start=start+shift;
     stop=stop+shift;
     n=n+1;
end

res=res/max(abs(res)); % Normalise amplitude

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [ imgOut ] = hyperOrthorectify( imgIn, altitude, hpbw )
%HYPERORTHORECTIFY Orthorectifies areal observed data.
%   Orthorectifies areal observed data using nearest neighbor interpolation.
%
% Inputs
%   imgIn       Input image (m x n) or (m x n x p)
%   altitude    Sensor altitude (meters)
%   hpbw        Half power beam width (radians).
% Outputs
%   imgOut      Orthorectified image.

% Input parameters
if (ndims(imgIn) == 2)
    [h, w] = size(imgIn);
    p = 1;
elseif (ndims(imgIn) == 3)
    [h, w, p] = size(imgIn);
end

radPerPix = hpbw/w;
x = tan(hpbw/2)*altitude;  % m
gsd = altitude*radPerPix;  % m
n = x/gsd;

imgOut = zeros(h, floor(n)*2, p);
for k=1:p
    for j=1:h
        for i=-floor(n):1:floor(n)-1
            boresiteDistance = gsd*i;
            theta = atan(boresiteDistance/altitude);
            imagePix = round(theta / radPerPix);
            imgOut(j, floor(n)+i+1, k) = imgIn(j, (w/2)+imagePix+1, k);
        end
    end
end

author: davidkun

repo: HyperSpectralToolbox


text:
function varargout = pivotplot(varargin)
%PIVOTPLOT   Semilogy plot of pivot values.
%   PIVOTPLOT( F ) semilogy plot related to the pivots values taken during
%   the construction of the CHEBFUN2 F.
%
%   H = PIVOTPLOT( F ) returns a handle H to the figure.
%
%   PIVOTPLOT( F, S ) allows further plotting options, such as linestyle,
%   linecolor, etc. If S contains a string 'LOGLOG', the pseudo sig will be
%   displayed on a log-log scale.

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

[varargout{1:nargout}] = pivotplot@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function hyperDemo_detectors
% HYPERDEMO_DETECTORS Demonstrates target detector algorithms
clear; clc; dbstop if error; close all;
%--------------------------------------------------------------------------
% Parameters
resultsDir = 'results\\';
dataDir = 'data\\AVIRIS\\';
%--------------------------------------------------------------------------

mkdir(resultsDir);

% Read part of AVIRIS data file that we will further process
M = hyperReadAvirisRfl(sprintf('%s\\f970620t01p02_r03_sc02.a.rfl', dataDir), [1 100], [1 614], [1 224]);
M = hyperNormalize(M);

% Read AVIRIS .spc file
lambdasNm = hyperReadAvirisSpc(sprintf('%s\\f970620t01p02_r03.a.spc', dataDir));

% Isomorph
[h, w, p] = size(M);
M = hyperConvert2d(M);

% Resample AVIRIS image.
desiredLambdasNm = 400:(2400-400)/(224-1):2400;
M = hyperResample(M, lambdasNm, desiredLambdasNm);

% Remove low SNR bands.
goodBands = [10:100 116:150 180:216];
M = M(goodBands, :);
p = length(goodBands);

% Demonstrate difference spectral similarity measurements
M = hyperConvert3d(M, h, w, p);
target = squeeze(M(11, 77, :));
figure; plot(desiredLambdasNm(goodBands), target); grid on;
    title('Target Signature; Pixel (32, 257)');

M = hyperConvert2d(M);

% RX Anomly Detector
r = hyperRxDetector(M);
r = hyperConvert3d(r.', h, w, 1);
figure; imagesc(r); title('RX Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\rx detector.png', resultsDir));

% Constrained Energy Minimization (CEM)
r = hyperCem(M, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(abs(r)); title('CEM Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\cem detector.png', resultsDir));

% Adaptive Cosine Estimator (ACE)
r = hyperAce(M, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(r); title('ACE Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\ace detector.png', resultsDir));

% Signed Adaptive Cosine Estimator (S-ACE)
r = hyperSignedAce(M, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(r); title('Signed ACE Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\signed ace detector.png', resultsDir));

% Matched Filter
r = hyperMatchedFilter(M, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(r); title('MF Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\mf detector.png', resultsDir));

% Generalized Likehood Ratio Test (GLRT) detector
r = hyperGlrt(M, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(r); title('GLRT Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\cem detector.png', resultsDir));

% Estimate background endmembers
U = hyperAtgp(M, 5);

% Hybrid Unstructured Detector (HUD)
r = hyperHud(M, U, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(abs(r)); title('HUD Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\hud detector.png', resultsDir));

% Adaptive Matched Subspace Detector (AMSD)
r = hyperAmsd(M, U, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(abs(r)); title('AMSD Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\amsd detector.png', resultsDir));
figure; mesh(r); title('AMSD Detector Results');

% Orthogonal Subspace Projection (OSP)
r = hyperOsp(M, U, target);
r = hyperConvert3d(r, h, w, 1);
figure; imagesc(abs(r)); title('OSP Detector Results'); axis image;
    colorbar;
hyperSaveFigure(gcf, sprintf('%s\\osp detector.png', resultsDir));


text:
function [gains]=lookup_gain_in_table(G,a_post,a_priori,a_post_range,a_priori_range,step);
% function [gains]=lookup_gain_in_table(G,a_post,a_priori,a_post_range,a_priori_range,step);
% This function selects the right gain value from the table G, given
% vectors with a priori and a posteriori SNRs
%
% INPUT variables:
% G: Matrix with gain values for speech DFT or magnitude estimation,
% evaluated at all combinations of a priori and a posteriori SNR in the
% input variables Rksi and Rgam.
%
%
% a_priori: Array of "a priori" SNR (SNRprior) values for which values
% have to be selected from the gain table   NOTE: The values must be in dBs.
% a_post: Array of "a posteriori" SNR (SNRpost) values for which values
% have to be selected from the gain table  NOTE: The values must be in dBs.
%
% a_post_range: The range of "a posteriori" SNR values
%
% a_priori_range: The range of "a priori" SNR values
%
% step: step is the stepsize in db's in the table
%
% OUTPUT variables:
% gains: Matrix with gain values that are selected from the gain table G
%
%
% Copyright 2007: Delft University of Technology, Information and
% Communication Theory Group. The software is free for non-commercial use.
% This program comes WITHOUT ANY WARRANTY.
%
% Last modified: 22-11-2007.

a_prioridb=round(10*log10(a_priori)/step)*step;
a_postdb=round(10*log10(a_post)/step)*step;
[Ia_post]=min(max(min(a_post_range),a_postdb), max(a_post_range));
Ia_post=Ia_post-min(a_post_range)+1;
Ia_post=Ia_post/step;
[Ia_priori]=min(max(min(a_priori_range),a_prioridb), max(a_priori_range));
Ia_priori=Ia_priori-min(a_priori_range)+1;
Ia_priori=Ia_priori/step;

gains=G(Ia_priori+(Ia_post-1)*length(G(:,1)));


text:
function [X,Y,beta_c,noise] = GenDataLinear(m,beta,v)

X         = randn(m, length(beta));
X         = X./repmat(sqrt(sum(X.^2,1)),size(X,1),1);
noise     = randn(m,1)/sqrt(v);
beta_c    = beta(:);
Y         = X*beta_c+noise;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
% FILE NAME:    Recon3Rxns2HMR.m
%
% PURPOSE: This script attempts to assign proper HMR id(s) to each
%          Recon3D rxn, and output the association into a defined
%          array structure.
%

% Load rxn association from HMR2 and Recon3D to MNX
load('Recon3Rxns2MNX.mat');    %Recon3D rxn association to MNX
load('ihumanRxns2MNX.mat');    %HMR rxn association to MNX

% Add assocation based on BiGG database
Recon3D.BiGG2HMR=cell(numel(Recon3D.rxns),1);
Recon3D.BiGG2HMR(:)={''};
index=find(~cellfun(@isempty,Recon3D.rxnBiGGID));
[a, b]=ismember(Recon3D.rxnBiGGID(index),ihuman.HMR2BiGG);
I=find(a);
Recon3D.BiGG2HMR(index(I))=ihuman.rxns(b(I));
numel(find(~cellfun(@isempty,Recon3D.BiGG2HMR)))  % ans = 4522
% Check consistency
indexBiGG=find(~cellfun(@isempty,Recon3D.BiGG2HMR));
indexHMR=find(~cellfun(@isempty,Recon3D.rxnHMRID));
overlap=intersect(indexBiGG,indexHMR);
ind=find(~cellfun(@isequal,Recon3D.BiGG2HMR(overlap),Recon3D.rxnHMRID(overlap)));
fprintf('%s: %s-%s\n',num2str(overlap(ind)),Recon3D.BiGG2HMR{overlap(ind)},Recon3D.rxnHMRID{overlap(ind)});
% 11741: HMR_8671-HMR_6533
% Keep the direct association based on Recon3 rxn id and remove the other one
Recon3D.BiGG2HMR{11741}='';
% These BiGG association seem to be problematic

% Add assocation based on MNX database

% Prepare the rxnCompIdx field to represent the rxn compartment info,
% because the rxnComps field cannot accurately descirbe cases, such as
% the reactions involve multiple compartments (e.g. exchagne/transport)

load('Recon3DRaven.mat');    %Load Recon3D in RAVEN format
Recon3D.rxnCompIdx=addRxnCompIdx(Recon3DRaven);
Recon3D.comps=Recon3DRaven.comps;
Recon3D.compNames=Recon3DRaven.compNames;
Recon3D.mets=Recon3DRaven.mets;
ihuman.rxnCompIdx=addRxnCompIdx(ihuman);

% Get additional Recon3D-HMR2 association using MNX association and
% compartment info, and then assign mapping results to rxnHMRID field
load('mergedModel.mat');
for i=1:numel(Recon3D.rxns)
		% Consider the rxns with MNX association and without HMR association
		if isempty(Recon3D.rxnHMRID{i}) && ~isempty(Recon3D.rxnMNXID{i})
				%attempt to get the HMR assoc
				for j=1:numel(mergedModel.rxns)
						if ismember(Recon3D.rxnMNXID{i},mergedModel.confirmedMNXID{j})
								rxns=[{};mergedModel.rxns{j}];
								if ~isempty(mergedModel.duplicateRxns{j}) && ....
								~ismember(mergedModel.rxns{j},{'HMR_8771','HMR_1592'})
										rxns=[rxns;transpose(strsplit(mergedModel.duplicateRxns{j},';'))];
								end

								% Specify rxn assoc using compartment index
								for k=1:numel(rxns)
										index=find(strcmp(rxns{k},ihuman.rxns));
										if isequal(ihuman.rxnCompIdx{index},Recon3D.rxnCompIdx{i})

												if isempty(Recon3D.rxnHMRID{i})
														Recon3D.rxnHMRID{i}=ihuman.rxns{index};
												else
														%Allow multiple HMR associations to one Recon3D rxn
														Recon3D.rxnHMRID{i}=strcat(Recon3D.rxnHMRID{i},';',ihuman.rxns{index});
												end
										end
								end

						end
				end
		end
end
numel(find(~cellfun(@isempty,Recon3D.rxnHMRID)))  % ans = 4998

save('Recon3Rxns2HMR.mat','Recon3D');    %Save the information for merging
save('ihumanRxns2MNX.mat','ihuman');     %2018-05-28

% flag the reactions with MNX association but without HMR association 2018-05-30
Recon3D.withMNXnoHMR=cell(numel(Recon3D.rxns),1);
Recon3D.withMNXnoHMR(:)={''};
for i=1:numel(Recon3D.rxns)
		if isempty(Recon3D.rxnHMRID{i}) && ~isempty(Recon3D.rxnMNXID{i})
				for j=1:numel(mergedModel.rxns)
						% locate these reactions extended to other compartments
						if ismember(Recon3D.rxnMNXID{i},mergedModel.confirmedMNXID{j}) & ~ismember(mergedModel.rxns{j},{'HMR_8771','HMR_1592'})
								Recon3D.withMNXnoHMR{i}=Recon3D.rxns{i};
								Recon3D.BiGG2HMR{i}='';   % clean off these HMR assoc
						end
				end
		end
end
numel(find(~cellfun(@isempty,Recon3D.withMNXnoHMR)))  % ans = 259
save('Recon3Rxns2HMR.mat','Recon3D');   % 2018-05-30

%===============
%% sub functions
% add compIdx field to accurately represent the rxn compartment info
function compIdx = addRxnCompIdx(model)
		compIdx=cell(numel(model.rxns),1);
		compIdx(:)={''};
		for i=1:numel(model.rxns)
				compIdx{i}=unique(model.metComps(find(model.S(:,i))));
		end
end

author: SysBioChalmers

repo: Human-GEM


text:
function check = logistic_check ( a, b )

%*****************************************************************************80
%
%% LOGISTIC_CHECK checks the parameters of the Logistic CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LOGISTIC_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
function [headmodel, cfg] = ft_prepare_singleshell(cfg, mri)

% FT_PREPARE_SINGLESHELL is deprecated, please use FT_PREPARE_HEADMODEL and
% FT_PREPARE_MESH
%
% See also FT_PREPARE_HEADMODEL

% TODO the spheremesh option should be renamed consistently with other mesh generation cfgs
% TODO shape should contain pnt as subfield and not be equal to pnt (for consistency with other use of shape)

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

warning('FT_PREPARE_SINGLESHELL is deprecated, please use FT_PREPARE_HEADMODEL with cfg.method = ''singleshell'' instead.')

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar mri
ft_preamble provenance mri

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'spheremesh', 'numvertices'});
cfg = ft_checkconfig(cfg, 'deprecated', 'mriunits');

% set the defaults
if ~isfield(cfg, 'threshold'),     cfg.threshold = 0.5;     end % relative
if ~isfield(cfg, 'numvertices'),   cfg.numvertices = [];    end % approximate number of vertices in sphere

% the data is specified as input variable or input file
hasmri = exist('mri', 'var');

if hasmri
  headmodel.bnd = ft_prepare_mesh(cfg, mri);
else
  headmodel.bnd = ft_prepare_mesh(cfg);
end

headmodel.type = 'singleshell';

% ensure that the geometrical units are specified
headmodel = ft_determine_units(headmodel);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug

ft_postamble previous   mri
ft_postamble provenance headmodel
ft_postamble history    headmodel


text:
function show_particles(X, Y_k)

figure(1)
image(Y_k)
title('+++ Showing Particles +++')

hold on
plot(X(2,:), X(1,:), '.')
hold off

drawnow


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

%% -*- texinfo -*-
%% @defun  assert_have_python_and_sympy (pyexec)
%% @defunx assert_have_python_and_sympy (pyexec, verbose)
%% Check that Python and SymPy are installed and working.
%%
%% @seealso{sympref}
%% @end defun

function assert_have_python_and_sympy (pyexec, verbose)

  minsympyver = '1.4';

  if (nargin < 2)
    verbose = false;
  end

  if (verbose)
    disp ('')
    disp ('Symbolic package diagnostics')
    disp ('============================')
    disp ('')
    disp ('Python and SymPy are needed for most features of the Symbolic package.')
    disp ('')
    fprintf ('The Python interpreter is currently: "%s".\n', pyexec)
    disp ('')
    disp ('Computers may have more than one Python interpreter installed.  If you')
    disp ('need to, you can select a different one using the PYTHON environment')
    disp ('variable (see "help sympref").  For example, to use Python 3, try')
    disp ('    setenv PYTHON python3')
    disp ('    sympref reset')
    disp ('')
    fprintf ('Attempting to run %s -c "print(\\"Python says hello\\")"\n\n', pyexec);
  end

  [status, output] = system ([pyexec ' -c "print(\"Python says hello\")"']);
  if (verbose)
    status
    output
  end

  if ( (status ~= 0) || (~ strcmp(strtrim(output), 'Python says hello')) )
    if (~ verbose)
      error ('OctSymPy:nopython', ...
             ['Cannot run the Python executable "%s"\n' ...
              '    Try "sympref diagnose" for more information.'], ...
             pyexec)
    end
    disp ('')
    disp ('Unfortunately, that command failed!')
    disp ('We expected to see "status = 0" and "output = Python says hello".')
    disp ('')
    disp ('  * Is there an error message above?')
    disp ('')
    disp ('  * Do you have Python installed?')
    disp ('')
    disp ('  * Please try using "setenv" as described above.')
    disp ('')
    disp ('  * Most systems search the PATH environment when looking for commands.')
    disp ('    Your path seems to be:')
    disp ('')
    disp (getenv ('PATH'))

    # if (ispc () && (~isunix ()))
    #   disp ('')
    #   disp ('  * Did you install the Symbolic bundle for Windows that includes Python?')
    #   disp ('    Please refer to the installation notes for Windows users at')
    #   disp ('    <https://github.com/cbm755/octsympy/wiki/Notes-on-Windows-installation>')
    # end
    return
  end

  if (verbose)
    disp ('Good, Python ran correctly.')
    disp ('')
    disp ('')
    disp ('Python version')
    disp ('--------------')
    disp ('')
    disp ('Let''s check what version of Python we are calling...')
    disp ('')
    fprintf ('Attempting to run %s -c "import sys; print(sys.version)"\n', pyexec);
    [status, output] = system ([pyexec ' -c "import sys; print(sys.version)"']);
    status
    output
    fprintf ('Running %s -c "import sys; print(sys.executable)"\n', pyexec);
    [status, output] = system ([pyexec ' -c "import sys; print(sys.executable)"']);
    output
    fprintf ('Checking if %s is running in a Cygwin-like POSIX environment...', pyexec);
    if python_env_is_cygwin_like (pyexec)
      fprintf (' yes\n\n');
    else
      fprintf (' no\n\n');
    end
    disp ('');
    show_system_info ();
  end

  if (verbose)
    disp ('')
    disp ('SymPy Python Library')
    disp ('--------------------')
    disp ('')
    disp ('SymPy is a Python library used by Symbolic for almost all features.')
    disp ('')
    fprintf ('Attempting to run %s -c "import sympy; print(sympy.__version__)"\n\n', pyexec);
  end

  [status, output] = system([pyexec ' -c "import sympy; print(sympy.__version__)"']);
  if (verbose)
    status
    output
  end

  if (status ~= 0)
    if (~ verbose)
      error ('OctSymPy:nosympy', ...
            ['Python cannot import SymPy: have you installed SymPy?\n' ...
             '    Try "sympref diagnose" for more information.'])
    end
    disp ('')
    disp ('Unfortunately status was non-zero: probably Python cannot import sympy.')
    disp ('')
    disp ('  * Is there an error message above?')
    disp ('')
    disp ('  * Do you have SymPy installed?  If not, please try to install it and')
    disp ('    try again.')
    disp ('')
    disp ('  * If you do have SymPy installed, maybe it''s installed for a different')
    disp ('    Python interpreter than the one we found?  Please try "setenv" as')
    disp ('    described above to change your python interpreter.')
    return
  end

  spver = strtrim(output);

  if (verbose)
    fprintf ('SymPy must be at least version %s; you have version %s.\n', ...
             minsympyver, spver);
  end

  if (~ exist ('OCTAVE_VERSION', 'builtin'))
    % no compare_versions on matlab, just assume its ok (!)
    if (verbose)
      disp ('We cannot easily compare versions on your system: please verify the above.')
    end
  else
    if (compare_versions (spver, minsympyver, '<'))
      if (~ verbose)
        error('OctSymPy:oldsympy', ...
              ['SymPy version %s found but is too old (%s required)\n' ...
               '    Try "sympref diagnose" for more information.'], ...
              spver, minsympyver)
      end
      disp ('**** Your SymPy is too old! ****')
      disp ('Installed newer version already?  Perhaps it was for a different Python?')
      disp ('Try "setenv" as described above to change your python interpreter.')
      return
    end
  end

  if (verbose)
    disp ('')
    disp ('Good, a working version of SymPy is installed.')
    disp ('')
    disp ('')
    disp ('Python XML Parsing and DOM Support')
    disp ('----------------------------------')
    disp ('')
    disp ('The XML DOM library is used by Symbolic for passing values to and from Python.')
    disp ('Some older versions of Python formatted XML output differently.  As long as you')
    disp ('have any reasonably recent version of Python, this should pass.')
    disp ('')
  end

  script = ['import xml.dom.minidom as minidom; ' ...
            'doc = minidom.parseString(\"<item>value</item>\"); ' ...
            'print(doc.toprettyxml(indent=\"\"))'];

  if (verbose)
    fprintf ('Attempting to run %s -c "%s"\n\n', pyexec, script);
  end

  if (ispc () && (~isunix ()))
    script = strrep (script, '<', '^<');
    script = strrep (script, '>', '^>');
  end

  [status, output] = system([pyexec ' -c "' script '"']);
  if (verbose)
    status
    output
  end

  if (status ~= 0)
    if (~ verbose)
      error ('OctSymPy:noxmldom', ...
             ['Something wrong, perhaps with XML DOM: is this an unusual Python setup?\n' ...
              '    Try "sympref diagnose" for more information.'])
    end
    disp ('')
    disp ('Unfortunately status was non-zero: perhaps Python cannot import xml.dom?')
    disp ('')
    disp ('  * Is there an error message above?')
    disp ('')
    disp ('  * Are you using an unusual Python installation?  If not, please try to')
    disp ('    reinstall it and try again.')
    return
  end

  xmlout = strtrim (output);

  if (isempty (strfind (xmlout, '<item>value</item>')))
    if (~ verbose)
      error ('OctSymPy:oldxmldom', ...
             ['Python XML output is not compatible, probably an old version of Python?\n' ...
              '    Try "sympref diagnose" for more information.'])
    end
    disp ('**** Your Python is too old! ****')
    disp ('')
    disp ('The XML output shown above must have "<item>value</item>" all on one line.')
    disp ('If your Python interpreter is an older release, please try upgrading to a more')
    disp ('recent version.  Symbolic should work with Python 2 (version 2.7.3 or newer)')
    disp ('and Python 3 (version 3.2.3 or newer).')
    return
  end

  if (verbose)
    fprintf ('\nYour kit looks good for running the Symbolic package.  Happy hacking!\n\n')
  end
end

repo: octsympy


text:
function w = get_w(sv)

% w = get_w()
% returns value of w from svm
% AE: i have changed this procedure to include non linear feature
% selection with RFE. When non-linear, it gives  a score vector for each
% feature (used in RFE) which is not the margin but is sorted as the margin
if strcmp(sv.ker{1},'linear'),

  w=sv.alpha'*get_x(sv.Xsv);

else
    Xtmp = get_x(sv.Xsv);
    alphatmp = sv.alpha;

    if strcmp(sv.ker{1},'rbf'),

      % compute the kernel matrix for all components
      K = Xtmp*Xtmp';
      Kdn = sum(Xtmp.^2,2);
      Kn = sum(Xtmp.^2,2);
      K = ones(size(Xtmp,1),1)*Kdn' + Kn*ones(1,size(Xtmp,1)) - 2*K;
      K = exp(-K/(2*sv.ker{2}));
      % compute the margin when one component is removed
       for i = 1:size(Xtmp,2),
           Ki = Xtmp(:,i)*ones(1,size(Xtmp,1)) - ones(size(Xtmp,1),1)*Xtmp(:,i)';
           Ki = Ki.^2;
           Ki = Ki/(2*sv.ker{2});
           Ki = exp(Ki);
           w(i) = alphatmp'*(K.*Ki)*alphatmp;
       end;

    elseif strcmp(sv.ker{1},'poly'),

        % compute the margin when one component is removed
        Ktmp = Xtmp*Xtmp';
        for i = 1:size(Xtmp,2),
           Ki = Xtmp(:,i)*Xtmp(:,i)';
           K_i = (Ktmp - Ki+1).^(sv.ker{2});
           w(i) = alphatmp'*K_i*alphatmp;
	end;
    end;% if strcmp(...,'rbf')
    w=max(w)-w;  %% make largest the smallest --- wrong way round!

end;% if strcmp(...,'linear')

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Yosra Boukari & Houssem Haddar (c) 2017-2018.|
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtIpbHelmholtz0.m                            |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & Yosra Boukari               |
%|  ( # )  |                & Houssem Haddar                              |
%|  / 0 \  |   CREATION   : 14.03.2017                                    |
%| ( === ) |   LAST MODIF : 14.03.2018                                    |
%|  `---'  |   SYNOPSIS   : Completion data for spherical helmholtz       |
%|         |                scatering                                     |
%|         |   ref1 : A Convergent Data Completion Algorithm Using Surface|
%|         |   Integral Equation, Inverse Problems, IOP Publishing ...    |
%|         |   ref2: Poly cours Terasse p. 194, Calderon operators        |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Wave number and frequency
k = 1
f = (k*340)/(2*pi)

% Data noise
noise = 0%1e-2

% Incident plane wave with normlized direction
X0         = [0 0 -1];
X0         = X0./norm(X0);
PW         = @(X) exp(1i*k*X*X0');
gradxPW{1} = @(X) 1i*k*X0(1) .* PW(X);
gradxPW{2} = @(X) 1i*k*X0(2) .* PW(X);
gradxPW{3} = @(X) 1i*k*X0(3) .* PW(X);

% Green kernel function --> G(x,y) = exp(ik|x-y|)/|x-y|
Gxy = @(X,Y) femGreenKernel(X,Y,'[exp(ikr)/r]',k);

% grady Green kernel function --> G(x,y) = grady[exp(ik|x-y|)/|x-y|]
dyGxy{1} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]1',k);
dyGxy{2} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]2',k);
dyGxy{3} = @(X,Y) femGreenKernel(X,Y,'grady[exp(ikr)/r]3',k);

% Green kernel function --> G(x,y) = gradx[exp(ik|x-y|)/|x-y|]
dxGxy{1} = @(X,Y) femGreenKernel(X,Y,'gradx[exp(ikr)/r]1',k);
dxGxy{2} = @(X,Y) femGreenKernel(X,Y,'gradx[exp(ikr)/r]2',k);
dxGxy{3} = @(X,Y) femGreenKernel(X,Y,'gradx[exp(ikr)/r]3',k);

% Finite element
gss = 3;
typ = 'P1';
tol = 1e-3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SCATERING PROBLEM %%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fictive mesh for diffraction problem (sphere)
mesh = mshSphere(4e2,1);

% % Fictive mesh for diffraction problem (cube)
% mesh = mshCube(4e2,[1 1 1]);
% mesh = mesh.bnd;

% Verify wave number
stp = mesh.stp;
if k > 1/stp(2)
    error('completionHelmholtz.m : wave number k is too high for mesh')
end

% Quadrature
sigma = dom(mesh,gss);

% Finite element
u = fem(mesh,typ);

% Boundary element operator on fictive domain (Hypersingular)
H = 1/(4*pi) .* (k^2 * integral(sigma,sigma,ntimes(u),Gxy,ntimes(u)) ...
    - integral(sigma,sigma,nxgrad(u),Gxy,nxgrad(u)));
H = H + 1/(4*pi) .* (k^2 * regularize(sigma,sigma,ntimes(u),'[1/r]',ntimes(u)) ...
    - regularize(sigma,sigma,nxgrad(u),'[1/r]',nxgrad(u)));

% Solve neumann problem on fictive boundary : - [H] mu = - dnP0
RHS = - integral(sigma,ntimes(u),gradxPW);
mu  = (-H) \ RHS;

% Radiation on boundary :  p(x) =  - [D + I/2] mu
Id  = integral(sigma,u,u);
D   = 1/(4*pi) .* integral(sigma,sigma,u,dyGxy,ntimes(u));
D   = D + 1/(4*pi) .* regularize(sigma,sigma,u,'grady[1/r]',ntimes(u));
sol = - Id \ ( (D + 0.5*Id) * mu);

% Comparare to analytic solution
ref = sphereHelmholtz('dom','neu',1,k,1.0001*mesh.vtx);
errScatMesh = norm(ref - sol)/norm(ref)

% Graphical representation
figure
plot(mesh)
hold on
plot(mesh,abs(sol+PW(u.dof)))
xlabel('X');   ylabel('Y');   zlabel('Z');
title('Total field')
axis equal
colorbar

%%%%%%%%%%%%%%%%%%%%%%%%%%%% SCATERING RADIATION %%%%%%%%%%%%%%%%%%%%%%%%%%
% Interior boundary
meshInt = mshSphere(5e2,1.25);
meshInt = swap(meshInt);

% Exterior boundary
meshExt = mshSphere(1e3,1.5);

% Fusion des maillages
meshTot = union(meshExt,meshInt);

% Verify wave number
stp = mesh.stp;
if k > 1/stp(2)
    error('completionHelmholtz.m : wave number k is too high for mesh')
end

% Quadrature
sigmaTot = dom(meshTot,gss);

% Finite element
uTot = fem(meshTot,typ);

% Mass matrix
IdTot = integral(sigmaTot,uTot,uTot);

% Solution of the scatering problem : p(x) = - [D] mu
p = - 1/(4*pi) .* integral(sigmaTot,sigma,uTot,dyGxy,ntimes(u)) * mu;

% Extract Galerkin
p = IdTot \ p;

% Compare to analytical solution
ref = sphereHelmholtz('dom','neu',1,k,uTot.dof);
errScat = norm(ref-p)/norm(ref)

% Scatered speed : dnp(x) = - [H] mu
dnp = - 1/(4*pi) .* (k^2 * integral(sigmaTot,sigma,ntimes(uTot),Gxy,ntimes(u)) ...
    - integral(sigmaTot,sigma,nxgrad(uTot),Gxy,nxgrad(u))) * mu;

% Extract Galerkin
dnp = IdTot \ dnp;

% % Graphical representation
% figure
% plot(mesh)
% hold on
% plot(mesh,abs(p+PW(uTot.dof)))
% plotNrm(mesh,'r')
% xlabel('X');   ylabel('Y');   zlabel('Z');
% title('Interior radiation')
% axis equal
% alpha(0.5)
% colorbar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OPERATOR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Single layer --> \int_Sx \int_Sy psi(x)' G(x,y) psi(y) dx dy
S = 1/(4*pi) .* integral(sigmaTot,sigmaTot,uTot,Gxy,uTot);
S = S + 1/(4*pi) .* regularize(sigmaTot,sigmaTot,uTot,'[1/r]',uTot);

% Double layer --> \int_Sx \int_Sy psi(x)' dny G(x,y) psi(y) dx dy
D = 1/(4*pi) .* integral(sigmaTot,sigmaTot,uTot,dyGxy,ntimes(uTot));
D = D + 1/(4*pi) .* regularize(sigmaTot,sigmaTot,uTot,'grady[1/r]',ntimes(uTot));

% Double layer --> \int_Sx \int_Sy psi(x)' dnx G(x,y) psi(y) dx dy
Dt = D.';

% Hypersingular --> k^2 * \int_Sx \int_Sy n.psi(x) G(x,y) n.psi(y) dx dy
%                   - \int_Sx \int_Sy nxgrad(psi(x)) G(x,y) nxgrad(psi(y)) dx dy
N = 1/(4*pi) .* (k^2 * integral(sigmaTot,sigmaTot,ntimes(uTot),Gxy,ntimes(uTot)) ...
    - integral(sigmaTot,sigmaTot,nxgrad(uTot),Gxy,nxgrad(uTot)));
N = N + 1/(4*pi) .* (k^2 * regularize(sigmaTot,sigmaTot,ntimes(uTot),'[1/r]',ntimes(uTot)) ...
    - regularize(sigmaTot,sigmaTot,nxgrad(uTot),'[1/r]',nxgrad(uTot)));

% Calderon operators
Z = sparse(length(uTot),length(uTot));
I = [IdTot Z ; Z IdTot];
H = [-D S ; -N Dt];

% Calderon operator for interior problem (projector)
A = H + 0.5*I;
errCalderon = norm(A*(I\A)-A,'fro')/norm(A,'fro')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INVERSE PROBLEM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Indices des dof
N    = size(meshTot.vtx,1);
Iext = 1:size(meshExt.vtx,1);
Iint = Iext(end)+1:N;

% Excitation
f = p(Iext);
g = dnp(Iext);

% Solutions
pSol   = p(Iint);
dnpSol = dnp(Iint);

% Indice operateur
Iext = [Iext , N + Iext];
Iint = [Iint , N + Iint];

% Extraction des sous matrices
A11 = A(Iext,Iext);
A12 = A(Iext,Iint);
A21 = A(Iint,Iext);
A22 = A(Iint,Iint);

% Matrices identite
I11 = I(Iext,Iext);
I22 = I(Iint,Iint);

% Linear system
LHS = [A12 ; A22 - I22];

% Right hand side
RHS = [ I11 - A11 ; -A21] * [f;g];

% Solve linear system with gaussian noise
X = LHS \ (RHS .* (1 + noise*(1 + randn(size(RHS,1),1))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TICHONOV %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% data  = F1;
% datan = B * F1n;
% [U,S,V]=svd(A);
% global s2 delta rhs2
% a1 = length(fe{2});
% rhs=(U')*datan;
% rhsr =rhs(1:a1+a1);
% s2 = diag(S).^2;
% rhs2 = abs(rhsr).^2;
% delta=norm(datan-B*(data))/norm(datan);
% alpha0=delta*min(s2)/(1-delta);
% alpha1=delta*max(s2)/(1-delta);
%
% fmor = @(x)sum(((x ./(x+s2)).^2 - delta^2).*rhs2);
% alphaM = fzero(fmor,[alpha0 alpha1]);
%
% X2 = V*((diag(S)./(alphaM+s2)).*rhsr);
% errTichonov = norm(X-X2)./norm(X2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ERROR ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare external pression
sol = X(1:end/2);
ref = pSol;
errSolution = norm(ref-sol)/norm(ref)

% Graphical representation
figure(13)
plot(meshInt,abs(sol+PW(meshInt.vtx)))
xlabel('X');   ylabel('Y');   zlabel('Z');
title('Reconstructed total field')
axis equal
colorbar

% Graphical representation
figure(14)
plot(meshInt,abs(sol-ref)./abs(ref))
xlabel('X');   ylabel('Y');   zlabel('Z');
title('Relative error on reconstructed field')
axis equal
colorbar

% Compare normal derivative
sol = X(end/2+1:end);
ref = dnpSol;
errDerivative = norm(ref-sol)/norm(ref)

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [xPred, PPred, yPred, S, Pxy] = KalmanFilterX_Predict(x,P,F,Q,H,R,u,B,O)
% KALMANFILTERX_PREDICT Perform the discrete-time KF state and measurement
% prediction steps, under the assumption of additive process noise.
%
% Parameters
% ----------
% x: column vector
%   The (xDim x 1) state estimate at the previous time-step.
% P: matrix
%   The (xDim x xDim) state covariance matrix at the previous
%   time-step.
% F: matrix
%   An (xDim x xDim) state transition matrix.
% Q: matrix
%   The (xDim x xDim) process noise covariance matrix.
% H: matrix
%   A (xDim x yDim) measurement matrix.
% R: matrix
%   The (yDim x yDim) measurement noise covariance matrix.
% u: column vector, optional
%   A optional (xDim x 1) control input.
%   If omitted, no control input is used.
% B: matrix, optional
%   An optional (xDim x xDim) control gain matrix.
%   If omitted, B is assumed to be 1.
% O: matrix, optional
%   An optional (xDim x xDim) control noise covariance
%   matrix. If omitted, Q is assumed to be 0.
%
% Returns
% -------
% xPred: column vector
%   The (xDim x 1) predicted state estimate.
% PPred: matrix
%   The (xDim x xDim) predicted state covariance matrix.
% yPred: column vector
%   The (yDim x 1) predicted measurement estimate.
% Pxy: matrix
%   The (xDim x yDim) cross-covariance matrix.
% S: matrix
%   The (yDim x yDim) innovation covariance matrix.
%
%October 2017 Lyudmil Vladimirov, University of Liverpool.

    switch(nargin)
        case(6)
            u = 0;
            B = 0;
            O = 0;
        case(7)
            B = 1;
            O = 0;
        case(8)
            O = 0;
    end

   [xPred, PPred] = KalmanFilterX_PredictState(x,P,F,Q,u,B,O);
   [yPred, S, Pxy] = KalmanFilterX_PredictObs(xPred,PPred,H,R);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
classdef dme_xfmr3p < mp.dm_element
%MP.DME_XFMR3P  MATPOWER data model class for 3-phase transformer data

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

    properties
        fbus    %% bus index vector for "from" port (port 1) (all branches)
        tbus    %% bus index vector for "to" port (port 2) (all branches)
        r       %% series resistance (p.u.) for branches that are on
        x       %% series reactance (p.u.) for branches that are on
        base_kva
        base_kv
%         g_to    %% shunt conductance (p.u.) at "to" end for branches that are on
%         b_fr    %% shunt susceptance (p.u.) at "from" end for branches that are on
%         b_to    %% shunt susceptance (p.u.) at "to" end for branches that are on
%         tm      %% transformer off-nominal turns ratio for branches that are on
%         ta      %% xformer phase-shift angle (radians) for branches that are on
%         rate_a  %% long term flow limit (p.u.) for branches that are on
    end     %% properties

    methods
        function name = name(obj)
            name = 'xfmr3p';
        end

        function label = label(obj)
            label = '3-ph Transformer';
        end

        function label = labels(obj)
            label = '3-ph Transformers';
        end

        function name = cxn_type(obj)
            name = 'bus3p';
        end

        function name = cxn_idx_prop(obj)
            name = {'fbus', 'tbus'};
        end

        function names = main_table_var_names(obj)
            names = horzcat( main_table_var_names@mp.dm_element(obj), ...
                {'bus_fr', 'bus_to', 'r', 'x', 'base_kva', 'base_kv', ...
                 'pl1_fr', 'ql1_fr', 'pl2_fr', 'ql2_fr', 'pl3_fr', 'ql3_fr', ...
                 'pl1_to', 'ql1_to', 'pl2_to', 'ql2_to', 'pl3_to', 'ql3_to' ...
                 });
        end

%         function vars = export_vars(obj)
%             vars = {'pl1_fr', 'ql1_fr', ...
%                     'pl2_fr', 'ql2_fr', ...
%                     'pl3_fr', 'ql3_fr', ...
%                     'pl1_to', 'ql1_to', ...
%                     'pl2_to', 'ql2_to', ...
%                     'pl3_to', 'ql3_to' };
%         end

        function obj = initialize(obj, dm)
            initialize@mp.dm_element(obj, dm);  %% call parent

            %% get bus mapping info
            b2i = dm.elements.bus3p.ID2i;   %% bus num to idx mapping

            %% set bus index vectors for port connectivity
            obj.fbus = b2i(obj.tab.bus_fr);
            obj.tbus = b2i(obj.tab.bus_to);
        end

        function obj = update_status(obj, dm)
            %% get bus status info
            bs = dm.elements.bus3p.tab.status;  %% bus status

            %% update status of branches connected to isolated/offline buses
            obj.tab.status = obj.tab.status & bs(obj.fbus) & ...
                                              bs(obj.tbus);

            %% call parent to fill in on/off
            update_status@mp.dm_element(obj, dm);
        end

        function obj = build_params(obj, dm)
            obj.r = obj.tab.r(obj.on);
            obj.x = obj.tab.x(obj.on);
            obj.base_kva = obj.tab.base_kva(obj.on);
            obj.base_kv  = obj.tab.base_kv(obj.on);
        end

        function obj = pretty_print(obj, dm, section, out_e, mpopt, fd, pp_args)
            switch section
                case 'det'
                    %% compute currents/powers for pp_args
                    s_fr = [  obj.tab.pl1_fr + 1j * obj.tab.ql1_fr ...
                            obj.tab.pl2_fr + 1j * obj.tab.ql2_fr ...
                            obj.tab.pl3_fr + 1j * obj.tab.ql3_fr    ];
                    s_to = [  obj.tab.pl1_to + 1j * obj.tab.ql1_to ...
                            obj.tab.pl2_to + 1j * obj.tab.ql2_to ...
                            obj.tab.pl3_to + 1j * obj.tab.ql3_to    ];
                    t = dm.elements.bus3p.tab;    %% bus3p table
                    vm = [ t.vm1 t.vm2 t.vm3 ] .* (t.base_kv/sqrt(3) * [1 1 1]);
                    va = [ t.va1 t.va2 t.va3 ];
                    v_ = vm .* exp(1j * va * pi/180);
                    i_fr = conj( s_fr ./ v_(obj.fbus, :));
                    i_to = conj( s_to ./ v_(obj.tbus, :));
                    c_fr = struct('cm', abs(i_fr), 'ca', angle(i_fr) * 180/pi);
                    c_to = struct('cm', abs(i_to), 'ca', angle(i_to) * 180/pi);

                    pp_args.xfmr3p = {'c', 'f', abs(i_fr), angle(i_fr) * 180/pi};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);

                    pp_args.xfmr3p = {'c', 't', abs(i_to), angle(i_to) * 180/pi};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);

                    pp_args.xfmr3p = {'s', 'f', real(s_fr), imag(s_fr)};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);

                    pp_args.xfmr3p = {'s', 't', real(s_to), imag(s_to)};
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);
                otherwise
                    pretty_print@mp.dm_element(obj, dm, section, out_e, mpopt, fd, pp_args);
            end
        end

        function TorF = pp_have_section_sum(obj, mpopt, pp_args)
            TorF = true;
        end

        function obj = pp_data_sum(obj, dm, rows, out_e, mpopt, fd, pp_args)
            %% call parent
            pp_data_sum@mp.dm_element(obj, dm, rows, out_e, mpopt, fd, pp_args);

            %% print generation summary
            t = obj.tab;
            ploss = [ t.pl1_fr t.pl2_fr t.pl3_fr ] + ...
                    [ t.pl1_to t.pl2_to t.pl3_to ];
            qloss = [ t.ql1_fr t.ql2_fr t.ql3_fr ] + ...
                    [ t.ql1_to t.ql2_to t.ql3_to ];
            fprintf(fd, '  %-29s %12.1f kW %12.1f kVAr\n', 'Total 3-ph transformer loss', ...
                sum(sum(ploss(obj.on, :))), sum(sum(qloss(obj.on, :))) );
        end

        function TorF = pp_have_section_det(obj, mpopt, pp_args)
            TorF = true;
        end

        function h = pp_get_headers_det(obj, dm, out_e, mpopt, pp_args)
            cs = pp_args.xfmr3p{1};
            ft = pp_args.xfmr3p{2};
            if cs == 'c' && ft == 'f'
                h1 = pp_get_headers_det@mp.dm_element(obj, dm, out_e, mpopt, pp_args);
            else
                h1 = {};
            end
            if cs == 'c'
                if ft == 'f'
                    h2 = {'-->  Current Injections at "From" Bus' };
                else    %% ft == 't'
                    h2 = {'', ...
                        '<--  Current Injections at "To" Bus' };
                end
            else    %% cs == 's'
                if ft == 'f'
                    h2 = {'', ...
                        '-->  Power Injections at "From" Bus' };
                else    %% ft == 't'
                    h2 = {'', ...
                        '<--  Power Injections at "To" Bus' };
                end
            end
            switch cs
                case 'c'
                    h = [ h1 h2 ...
                        {   '  3-ph    3-ph Bus  3-ph Bus          Phase A Current  Phase B Current  Phase C Current', ...
                            'Xfrm ID    From ID   To ID    Status   (A)    (deg)     (A)    (deg)     (A)    (deg)', ...
                            '--------  --------  --------  ------  ------  ------   ------  ------   ------  ------' } ];
                    %%       1234567 123456789 123456789 -----1 123456.89 12345.7 12345.78 12345.7 12345.78 12345.7
                case 's'
                    h = [ h1 h2 ...
                        {   '  3-ph    3-ph Bus  3-ph Bus          Phase A Power    Phase B Power    Phase C Power', ...
                            'Xfmr ID    From ID   To ID    Status   (kW)   (kVAr)    (kW)   (kVAr)    (kW)   (kVAr)', ...
                            '--------  --------  --------  ------  ------  ------   ------  ------   ------  ------' } ];
                    %%       1234567 123456789 123456789 -----1 1234567.9 12345.7 123456.8 12345.7 123456.8 12345.7
            end
        end

        function str = pp_data_row_det(obj, dm, k, out_e, mpopt, fd, pp_args)
            switch pp_args.xfmr3p{1}    %% cs
                case 'c'
                    cm = pp_args.xfmr3p{3}(k, :);
                    ca = pp_args.xfmr3p{4}(k, :);
                    str = sprintf('%7d %9d %9d %6d %9.2f %7.1f %8.2f %7.1f %8.2f %7.1f', ...
                        obj.tab.uid(k), obj.tab.bus_fr(k), obj.tab.bus_to(k), ...
                        obj.tab.status(k), ...
                        cm(:, 1), ca(:, 1), ...
                        cm(:, 2), ca(:, 2), ...
                        cm(:, 3), ca(:, 3) );
                case 's'
                    p = pp_args.xfmr3p{3}(k, :);
                    q = pp_args.xfmr3p{4}(k, :);
                    str = sprintf('%7d %9d %9d %6d %9.1f %7.1f %8.1f %7.1f %8.1f %7.1f', ...
                        obj.tab.uid(k), obj.tab.bus_fr(k), obj.tab.bus_to(k), ...
                        obj.tab.status(k), ...
                        p(:, 1), q(:, 1), ...
                        p(:, 2), q(:, 2), ...
                        p(:, 3), q(:, 3) );
            end
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function precisions = precision_plot(positions, ground_truth, title, show)
%PRECISION_PLOT
%   Calculates precision for a series of distance thresholds (percentage of
%   frames where the distance to the ground truth is within the threshold).
%   The results are shown in a new figure if SHOW is true.
%
%   Accepts positions and ground truth as Nx2 matrices (for N frames), and
%   a title string.
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/

	max_threshold = 50;  %used for graphs in the paper

	precisions = zeros(max_threshold, 1);

	if size(positions,1) ~= size(ground_truth,1),
% 		fprintf('%12s - Number of ground truth frames does not match number of tracked frames.\n', title)

		%just ignore any extra frames, in either results or ground truth
		n = min(size(positions,1), size(ground_truth,1));
		positions(n+1:end,:) = [];
		ground_truth(n+1:end,:) = [];
	end

	%calculate distances to ground truth over all frames
	distances = sqrt((positions(:,1) - ground_truth(:,1)).^2 + ...
				 	 (positions(:,2) - ground_truth(:,2)).^2);
	distances(isnan(distances)) = [];

	%compute precisions
	for p = 1:max_threshold,
		precisions(p) = nnz(distances <= p) / numel(distances);
	end

	%plot the precisions
	if show == 1,
		figure('NumberTitle','off', 'Name',['Precisions - ' title])
		plot(precisions, 'k-', 'LineWidth',2)
		xlabel('Threshold'), ylabel('Precision')
	end

end

repo: Context-Aware-CF-Tracking