text:
% fdr() - compute false detection rate mask
%
% Usage:
%   >> [p_fdr, p_masked] = fdr( pvals, alpha);
%
% Inputs:
%   pvals   - vector or array of p-values
%   alpha   - threshold value (non-corrected). If no alpha is given
%             each p-value is used as its own alpha and FDR corrected
%             array is returned.
%   fdrtype - ['parametric'|'nonParametric'] FDR type. Default is
%             'parametric'.
%
% Outputs:
%   p_fdr    - pvalue used for threshold (based on independence
%              or positive dependence of measurements)
%   p_masked - p-value thresholded. Same size as pvals.
%
% Author: Arnaud Delorme, SCCN, 2008-
%         Based on a function by Tom Nichols
%
% Reference: Bejamini & Yekutieli (2001) The Annals of Statistics

% Copyright (C) 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [pID, p_masked] = fdr(pvals, q, fdrType);

if nargin < 3, fdrType = 'parametric'; end;
if isempty(pvals), pID = []; return; end;
p = sort(pvals(:));
V = length(p);
I = (1:V)';

cVID = 1;
cVN = sum(1./(1:V));

if nargin < 2
    pID = ones(size(pvals));
    thresholds = exp(linspace(log(0.1),log(0.000001), 100));
    for index = 1:length(thresholds)
        [tmp p_masked] = fdr(pvals, thresholds(index));
        pID(p_masked) = thresholds(index);
    end;
else
    if strcmpi(fdrType, 'parametric')
        pID = p(max(find(p<=I/V*q/cVID))); % standard FDR
    else
        pID = p(max(find(p<=I/V*q/cVN)));  % non-parametric FDR
    end;
end;
if isempty(pID), pID = 0; end;

if nargout > 1
    p_masked = pvals<=pID;
end;

author: PatternRecognition

repo: OpenBMI


text:
classdef MimRemapBlue < MimGuiPlugin
    % MimRemapBlue. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimRemapBlue is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Right upper'
        SelectedText = 'Right upper'
        ToolTip = 'Marks airway as right upper'
        Category = 'Airway label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{1}
        Location = 31
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 1;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Map');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 1;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%This Matlab script can be used to reproduce Figure 6.7 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Number of BSs
L = 16;

%Number of UEs per BS
K = 10;

%Number of BS antennas
M = 100;

%Define the range of pilot reuse factors
fRange = [1 2];

%Select the number of setups with random UE locations
nbrOfSetups = 20;

%Select the number of channel realizations per setup
nbrOfRealizations = 200;

%Select the range of hardware qualities. The values at the same position in
%the different vectors are considered simultaneously.
kappatUE = [0.99 0.99 0.99 0.99 0.99 0.95 0.97 0.98 0.99 0.995 1];
kapparUE = [0.95 0.97 0.98 0.995 1 0.99 0.99 0.99 0.99 0.99 0.99];
kappatBS = kapparUE;
kapparBS = kappatUE;

%% Propagation parameters

%Communication bandwidth
B = 20e6;

%Total uplink transmit power per UE (mW)
p = 100;

%Total downlink transmit power per UE (mW)
rho = 100;

%Noise figure at the BS (in dB)
noiseFigure = 7;

%Compute noise power
noiseVariancedBm = -174 + 10*log10(B) + noiseFigure;

%Select length of coherence block
tau_c = 200;

%Use the approximation of the Gaussian local scattering model
accuracy = 2;

%Angular standard deviation in the local scattering model (in degrees)
ASDdeg = 10;

%Prepare to save simulation results
sumSE_MR = zeros(length(kappatUE),length(fRange),nbrOfSetups);
sumSE_RZF = zeros(length(kappatUE),length(fRange),nbrOfSetups);
sumSE_MMMSE = zeros(length(kappatUE),length(fRange),nbrOfSetups);

%% Go through all setups
for n = 1:nbrOfSetups

    %Output simulation progress
    disp([num2str(n) ' setups out of ' num2str(nbrOfSetups)]);

    %Compute channel statistics for one setup
    [R,channelGaindB] = functionExampleSetup(L,K,M,accuracy,ASDdeg);

    %Compute the normalized average channel gain, where the normalization
    %is based on the noise power
    channelGainOverNoise = channelGaindB - noiseVariancedBm;

    %Go through all pilot reuse factors
    for s = 1:length(fRange)

        %Extract pilot reuse factor
        f = fRange(s);

        %Output simulation progress
        disp([num2str(f) ' reuse factor out of ' num2str(length(fRange))]);

        %Go through all hardware qualities
        for r = 1:length(kappatUE)

            %Generate channel realizations with estimates and estimation
            %error correlation matrices, using LMMSE estimation
            [Hhat,C,tau_p,~,H] = functionChannelEstimates_impairments(R,channelGainOverNoise,nbrOfRealizations,M,K,L,p,f,kappatUE(r),kapparBS(r));

            %Compute SE using Monte-Carlo realizations
            [SE_MR,SE_RZF,SE_MMMSE] = functionComputeSE_DL_impairments(H,Hhat,C,tau_c,tau_p,nbrOfRealizations,M,K,L,p,rho,kappatBS(r),kapparUE(r));

            %Save results
            sumSE_MR(r,s,n) = mean(sum(SE_MR,1));
            sumSE_RZF(r,s,n) = mean(sum(SE_RZF,1));
            sumSE_MMMSE(r,s,n) = mean(sum(SE_MMMSE,1));

            %Delete large matrices
            clear Hhat C H;

        end

    end

    %Delete large matrices
    clear R;

end

%Plot the simulation results
figure;
hold on; box on;

plot(kapparUE([1 2 3 9 4 5]),max(mean(sumSE_MMMSE([1 2 3 9 4 5],:,:),3),[],2),'rd-','LineWidth',1);
plot(kapparUE([1 2 3 9 4 5]),max(mean(sumSE_RZF([1 2 3 9 4 5],:,:),3),[],2),'r-.','LineWidth',1);
plot(kapparUE([1 2 3 9 4 5]),max(mean(sumSE_MR([1 2 3 9 4 5],:,:),3),[],2),'bs-','LineWidth',1);

xlabel('DL hardware quality');
ylabel('Average sum SE [bit/s/Hz/cell]');
legend('M-MMSE','RZF','MR','Location','SouthEast');
xlim([0.95 1]);
ylim([0 40]);

%Plot the simulation results
figure;
hold on; box on;

plot(kappatUE(6:11),max(mean(sumSE_MMMSE(6:11,:,:),3),[],2),'rd-','LineWidth',1);
plot(kappatUE(6:11),max(mean(sumSE_RZF(6:11,:,:),3),[],2),'r-.','LineWidth',1);
plot(kappatUE(6:11),max(mean(sumSE_MR(6:11,:,:),3),[],2),'bs-','LineWidth',1);

xlabel('UL hardware quality');
ylabel('Average sum SE [bit/s/Hz/cell]');
legend('M-MMSE','RZF','MR','Location','SouthEast');
xlim([0.95 1]);
ylim([0 40]);

author: emilbjornson

repo: massivemimobook


text:
function A = HALS_spatial_thresh(Y, A, C, active_pixel, maxIter, sn)
%% run HALS by fixating all spatial components
% input:
%   Y:  d*T, fluorescence data
%   A:  d*K, spatial components
%   C:  K*T, temporal components
%   active_pixel, mask for pixels to be updated

% output:
%   A: d*K, updated spatial components

% Author: Pengcheng Zhou, Carnegie Mellon University, adapted from Johannes
% Friedrich's NIPS paper "Fast Constrained Non-negative Matrix
% Factorization for Whole-Brain Calcium Imaging Data

%% options for HALS
if nargin<5;    maxIter = 1;    end   %maximum iteration number
if nargin<4;    active_pixel=true(size(A));
elseif isempty(active_pixel)
    active_pixel = true(size(A));
else
    active_pixel = logical(active_pixel);
end     %determine nonzero pixels
if ~exist('sn', 'var')||isempty(sn)
    sn = GetSn(Y);
end
sn = reshape(sn, [], 1);
%% initialization
A(~active_pixel) = 0;
K = size(A, 2);     % number of components
Cmean = mean(C,2);
Ymean = mean(Y,2);
T = size(C,2);
U = double(Y*C'-T*(Ymean*Cmean'));
V = double(C*C'-T*(Cmean*Cmean'));
cc = diag(V);   % squares of l2 norm all all components
cc_thr = 3.0./sqrt(cc);

%% updating
for miter=1:maxIter
    for k=1:K
        if cc(k)==0
            continue;
        end
        tmp_ind = active_pixel(:, k);
        if sum(tmp_ind)==0
            A(:, k) = 0;
            continue;
        end
        ak = A(tmp_ind, k)+(U(tmp_ind, k)-A(tmp_ind,:)*V(:, k))/cc(k);
        ak(ak<sn(tmp_ind)*cc_thr(k)) = 0;
        A(tmp_ind, k) = ak;
    end
end

author: zhoupc


text:
function linpack_z_test37 ( )

%*****************************************************************************80
%
%% TEST37 tests ZTRSL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST37\n' );
  fprintf ( 1, '  For a double precision complex (C)\n' );
  fprintf ( 1, '  triangular matrix (TR),\n' );
  fprintf ( 1, '  ZTRSL solves a linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the matrix.
%
  seed = 123456789;

  for i = 1 : n
    for j = 1 : i
      [ a(i,j), seed ] = c8_uniform_01 ( seed );
    end
    a(i,i+1:n) = 0.0;
  end
%
%  Set the desired solution
%
  for i = 1 : n
    x(i) = complex ( i, 10 * i );
  end
%
%  Compute the corresponding right hand side.
%
  b(1:n) = a(1:n,1:n) * transpose ( x(1:n) );
%
%  Solve the lower triangular system.
%
  job = 0;
  [ b, info ] = ztrsl ( a, lda, n, b, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed                     Exact\n' );
  fprintf ( 1, '  Solution                     Solution\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  (%8f  %8f)  (%8f  %8f)\n', ...
      real ( b(i) ), imag ( b(i) ), real ( x(i) ), imag ( x(i) ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% distance to wall
function [data,units] = compute_dist2wall_rect(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  x = trx(fly).x_mm;
  y = trx(fly).y_mm;

  % The ordering of the arguments is to account for weird y directions of images.
  dtop = getDist(trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly), trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly),...
    x,y);
  dleft = getDist(trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly), trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly),...
    x,y);
  dright = getDist(trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly), trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly),...
    x,y);
  dbottom = getDist(trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly), trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly),...
    x,y);

  data{i} = min([dtop;dleft; dright; dbottom],[],1);

end
units = parseunits('mm');

function d = getDist(p1_x,p1_y,p2_x,p2_y,p_x,p_y)

dp1p2 = (p1_x-p2_x).^2 + (p1_y-p2_y).^2;

dotpr_x = (p1_x - p_x).*(p1_x-p2_x);
dotpr_y = (p1_y - p_y).*(p1_y-p2_y);

t = (dotpr_x+dotpr_y)/dp1p2;
proj_x = p1_x + t.*(p2_x-p1_x);
proj_y = p1_y + t.*(p2_y-p1_y);
d = sqrt(  (p_x-proj_x).^2 + (p_y-proj_y).^2);

% Find whether the mice is in the interior of the square.
side = sign((p2_x-p1_x)*(p_y-p1_y)-(p2_y-p1_y)*(p_x-p1_x));
d = side.*d;

author: kristinbranson

repo: JAABA


text:
function [pnt] = elec1020_fraction(cnt1, cnt2, fraction)

% ELEC1020_FRACTION

% Copyright (C) 2003, 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/>.

ncnt = size(cnt1,1);

% determine the total length of the contour
tot_l = 0;
for i=1:ncnt
  tot_l = tot_l + pntdist(cnt1(i,:), cnt2(i,:));
end

frac_l = fraction * tot_l;

% propagate along the contour untill we get at the desired fraction
sum_l = 0;
for i=1:ncnt
  seg_l = pntdist(cnt1(i,:), cnt2(i,:));
  if (sum_l+seg_l)>=frac_l
    % the desired point lies on this segment
    la = frac_l - sum_l;
    vec = cnt2(i,:)-cnt1(i,:);
    pnt = cnt1(i,:) + la * vec/norm(vec);
    return
  else
    sum_l = sum_l + seg_l;
    sum_f = sum_l/tot_l;
  end
end

pnt = [nan nan nan];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function [fnormhat,t]=instfreq(x,t,L,trace);
%INSTFREQ Instantaneous frequency estimation.
%	[FNORMHAT,T]=INSTFREQ(X,T,L,TRACE) computes the instantaneous
%	frequency of the analytic signal X at time instant(s) T, using the
%	trapezoidal integration rule.
%	The result FNORMHAT lies between 0.0 and 0.5.
%
%	X : Analytic signal to be analyzed.
%	T : Time instants	        (default : 2:length(X)-1).
%	L : If L=1, computes the (normalized) instantaneous frequency
%	    of the signal X defined as angle(X(T+1)*conj(X(T-1)) ;
%	    if L>1, computes a Maximum Likelihood estimation of the
%	    instantaneous frequency of the deterministic part of the signal
%	    blurried in a white gaussian noise.
%	    L must be an integer       	(default : 1).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                (default : 0).
%	FNORMHAT : Output (normalized) instantaneous frequency.
%	T : Time instants.
%
%	Examples :
%	 x=fmlin(70,0.05,0.35,25); [instf,t]=instfreq(x); plot(t,instf)
%	 N=64; SNR=10.0; L=4; t=L+1:N-L; x=fmsin(N,0.05,0.35,40);
%	 sig=sigmerge(x,hilbert(randn(N,1)),SNR);
%	 plotifl(t,[instfreq(sig,t,L),instfreq(x,t)]); grid;
%	 title ('theoretical and estimated instantaneous frequencies');
%
%	See also  KAYTTH, SGRPDLAY.

%	F. Auger, March 1994, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%	------------------- CONFIDENTIAL PROGRAM --------------------
%	This program can not be used without the authorization of its
%	author(s). For any comment or bug report, please send e-mail to
%	f.auger@ieee.org

if (nargin == 0),
 error('At least one parameter required');
end;
[xrow,xcol] = size(x);
if (xcol~=1),
 error('X must have only one column');
end

if (nargin == 1),
 t=2:xrow-1; L=1; trace=0.0;
elseif (nargin == 2),
 L = 1; trace=0.0;
elseif (nargin == 3),
 trace=0.0;
end;

if L<1,
 error('L must be >=1');
end
[trow,tcol] = size(t);
if (trow~=1),
 error('T must have only one row');
end;

if (L==1),
 if any(t==1)|any(t==xrow),
  error('T can not be equal to 1 neither to the last element of X');
 else
  fnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi);
 end;
else
 H=kaytth(L);
 if any(t<=L)|any(t+L>xrow),
  error('The relation L<T<=length(X)-L must be satisfied');
 else
  for icol=1:tcol,
   if trace, disprog(icol,tcol,10); end;
   ti = t(icol); tau = 0:L;
   R = x(ti+tau).*conj(x(ti-tau));
   M4 = R(2:L+1).*conj(R(1:L));

   diff=2e-6;
   tetapred = H * (unwrap(angle(-M4))+pi);
   while tetapred<0.0 , tetapred=tetapred+(2*pi); end;
   while tetapred>2*pi, tetapred=tetapred-(2*pi); end;
   iter = 1;
   while (diff > 1e-6)&(iter<50),
    M4bis=M4 .* exp(-j*2.0*tetapred);
    teta = H * (unwrap(angle(M4bis))+2.0*tetapred);
    while teta<0.0 , teta=(2*pi)+teta; end;
    while teta>2*pi, teta=teta-(2*pi); end;
    diff=abs(teta-tetapred);
    tetapred=teta; iter=iter+1;
   end;
   fnormhat(icol,1)=teta/(2*pi);
  end;
 end;
end;


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_SCARA(robot, T)
%   Solves the inverse kinematic problem for the SCARA example robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_SCARA returns 2 possible solutions, thus,
%   Q is a 4x4 matrix where each column stores 4 feasible joint values.
%
%
%   Example code:
%
%   robot=load_robot('example', 'scara');
%   q = [0 0 0 0];
%   T = directkinematic(robot, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(robot, T);
%   %check that all of them are feasible solutions!
%   %and every Ti equals T
%   for i=1:2,
%        Ti = directkinematic(robot, qinv(:,i))
%   end
%
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

fprintf('\nComputing inverse kinematics for the %s robot', robot.name);

%initialize q
q=zeros(4,2);

%Evaluate the DH table to obtain geometric parameters
d = eval(robot.DH.d);
a = eval(robot.DH.a);

%Store geometric parameters
L1=abs(d(1));
L2=abs(a(1));
L3=abs(a(2));

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Distance of the point to the origin of S0
R = sqrt(Px^2+Py^2);

%Compute angles
gamma = real(acos((L2^2+R^2-L3^2)/(2*R*L2)));
beta = atan2(Py,Px);
delta = real(acos((L2^2+L3^2-R^2)/(2*L2*L3)));

%find the last rotation for the two possible configurations
q4_1= find_last_rotation(robot,[beta+gamma delta-pi L1-Pz 0], T);
q4_2= find_last_rotation(robot,[beta-gamma pi-delta L1-Pz 0], T);

%Arrange all possible solutions
q=[beta+gamma beta-gamma;
    delta-pi pi-delta;
    L1-Pz    L1-Pz;
    q4_1 q4_2];

% Compute the last rotation
function q4 = find_last_rotation(robot, q, T)

U = T(1:3,1);

%Recompute the DH table according to q1, q2 and q3
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%now compute the position/orientation of the system S3
H=eye(4);
for i=1:3,
    H=H*dh(theta(i), d(i), a(i), alpha(i));
end

X3=H(1:3,1);
Y3=H(1:3,2);

coseno=X3'*U;
seno=U'*Y3;
%compute the last rotation
q4=atan2(seno,coseno);

repo: ARTE


text:
report_this_filefun(mfilename('fullpath'));

ti = [];
ret = [];
figure_w_normalized_uicontrolunits(bmapc)
[px,py]  = ginput(1)
hold on
plot(px,py,'w+','era','normal')
drawnow
for ra = 1:0.25:15
    ra
    retpoint
    l = isnan(ret);
    ret(l) = [];

    [i,j] = find(ret == min(min(ret)));
    if isempty(ret) == 1; ret = nan;end

    tru = (teb-t0b)/(10^(bvg(i,8)-(bvg(i,1)-1.3*bvg(i,9))*6));
    trl = (teb-t0b)/(10^(bvg(i,8)-(bvg(i,1)+1.3*bvg(i,9))*6));

    ti = [ti ; min(ret) ra bvg(i,1)  bvg(i,9)  bvg(i,6) bvg(i,8)  -tru+min(ret) trl-min(ret)];
end

figure
rect = [0.20, 0.70, 0.6, 0.25];
axes('position',rect)
%pl = plot(ti(:,2),ti(:,3))
errorbar(ti(:,2),ti(:,3),ti(:,4),ti(:,4))
%set(pl,'LineWidth',2)
hold on
pl = plot(ti(:,2),ti(:,3),'o')
set(pl,'LineWidth',2)
ylabel('b-value')
xl = get(gca,'Xlim');
xl = [0 16];
set(gca,'Xlim',xl);
matdraw

rect = [0.20, 0.40, 0.6, 0.25];
axes('position',rect)
errorbar(ti(:,2),ti(:,6),2*ti(:,4),2*ti(:,4))
%pl = plot(ti(:,2),ti(:,6))
%set(pl,'LineWidth',2)
hold on
pl = plot(ti(:,2),ti(:,6),'o')
set(pl,'LineWidth',2)
ylabel('a-value')
set(gca,'Xlim',xl);

rect = [0.20, 0.10, 0.6, 0.25];
axes('position',rect)
errorbar(ti(:,2),ti(:,1),ti(:,7),ti(:,8))
%pl = plot(ti(:,2),ti(:,1))
set(pl,'LineWidth',2)
set(gca,'Xlim',xl,'Ylim',[0 150]);
hold on
pl = plot(ti(:,2),ti(:,1),'o')
set(pl,'LineWidth',2)
ylabel('Tr [yrs.]')

xlabel('Radius in [km]')

author: CelsoReyes


text:
function result = ball_f1_nd ( func, n, xc, r )

%*****************************************************************************80
%
%% BALL_F1_ND approximates an integral inside a ball in ND.
%
%  Integration region:
%
%    Points X(1:N) such that:
%
%      Sum ( X(1:N) - XC(1:N) )**2 <= R**2.
%
%  Discussion:
%
%    An (N+1)*2**N point 5-th degree formula is used, Stroud number SN:5-6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F at the N-vector X, of the form
%      function value = func ( n, x )
%
%    Input, integer N, the dimension of the space.
%
%    Input, real XC(N), the center of the ball.
%
%    Input, real R, the radius of the ball.
%
%    Output, real RESULT, the approximate integral of the function.
%
  if ( r == 0.0E+00 )
    result = 0.0E+00;
    return
  end

  u2 = ( 1.0E+00 - 2.0E+00 * sqrt ( 1.0E+00 / ( n + 4 ) ) ) / ( n + 2 );
  u = sqrt ( u2 );
  x(1:n) = xc(1:n) - r * u;

  w = 1.0E+00 / ( ( n + 1 ) * 2^n );

  quad = 0.0E+00;
  ihi = 2^n;

  for i = 1 : ihi

    itemp = i - 1;

    for j = 1 : n

      u = ( xc(j) - x(j) ) / r;

      if ( mod ( itemp, 2 ) == 1 )
        x(j) = xc(j) - abs ( x(j) - xc(j) );
      else
        x(j) = xc(j) + abs ( x(j) - xc(j) );
      end

      itemp = floor ( itemp / 2 );

    end

    quad = quad + w * feval ( func, n, x );

  end

  temp = sqrt ( n + 4 );

  t = sqrt ( 2.0E+00 * ( n + 1 ) / ( n + 2 ) ) / ( n * temp );

  y = ( 1.0E+00 + 2.0E+00 / ( n * temp ) ) / ( n + 2 );
  v = sqrt ( y - t );
  u = sqrt ( y + ( n - 1 ) * t );

  khi = 2^n;

  for i = 1 : n

    x(1:n) = xc(1:n) - r * v;

    x(i) = xc(i) - r * u;

    for k = 1 : khi

      ktemp = k - 1;

      for j = 1 : n

        if ( mod ( ktemp, 2 ) == 1 )
          x(j) = xc(j) - abs ( x(j) - xc(j) );
        else
          x(j) = xc(j) + abs ( x(j) - xc(j) );
        end

        ktemp = floor ( ktemp / 2 );

      end

      quad = quad + w * feval ( func, n, x );

    end

    x(i) = xc(i) - r * v;

  end

  volume = ball_volume_nd ( n, r );
  result = quad * volume;

  return
end


text:
function hf_out = lhs_operation(hf, samplesf, reg_filter, sample_weights)

% This is the left-hand-side operation in Conjugate Gradient

% Get sizes
num_features = length(hf);
filter_sz = zeros(num_features,2);
for k = 1:num_features
    filter_sz(k,:) = [size(hf{k},1), size(hf{k},2)];
end
[~, k1] = max(filter_sz(:,1));  % Index for the feature block with the largest spatial size
block_inds = 1:num_features;
block_inds(k1) = [];
output_sz = [size(hf{k1},1), 2*size(hf{k1},2)-1];

% Compute the operation corresponding to the data term in the optimization
% (blockwise matrix multiplications)
%implements: A' diag(sample_weights) A f

% sum over all features and feature blocks
sh = mtimesx(samplesf{k1}, permute(hf{k1}, [3 4 1 2]), 'speed');    % assumes the feature with the highest resolution is first
pad_sz = cell(1,1,num_features);
for k = block_inds
    pad_sz{k} = (output_sz - [size(hf{k},1), 2*size(hf{k},2)-1]) / 2;

    sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) = ...
        sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) + mtimesx(samplesf{k}, permute(hf{k}, [3 4 1 2]), 'speed');
end

% weight all the samples
sh = bsxfun(@times,sample_weights,sh);

% multiply with the transpose
hf_out = cell(1,1,num_features);
hf_out{k1} = permute(conj(mtimesx(sh, 'C', samplesf{k1}, 'speed')), [3 4 2 1]);
for k = block_inds
    hf_out{k} = permute(conj(mtimesx(sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end), 'C', samplesf{k}, 'speed')), [3 4 2 1]);
end

% compute the operation corresponding to the regularization term (convolve
% each feature dimension with the DFT of w, and the tramsposed operation)
% add the regularization part
% hf_conv = cell(1,1,num_features);
for k = 1:num_features
    reg_pad = min(size(reg_filter{k},2)-1, size(hf{k},2)-1);

    % add part needed for convolution
    hf_conv = cat(2, hf{k}, conj(rot90(hf{k}(:, end-reg_pad:end-1, :), 2)));

    % do first convolution
    hf_conv = convn(hf_conv, reg_filter{k});

    % do final convolution and put toghether result
    hf_out{k} = hf_out{k} + convn(hf_conv(:,1:end-reg_pad,:), reg_filter{k}, 'valid');
end

end

author: martin-danelljan

repo: ECO


text:
%VERTCAT Mapping overload
%
%  Vertical concatenation of mappings is performed by PARALLEL

author: marianux

repo: ecg-kit

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


text:
% Nonnegativity Constrained Least Squares with Multiple Righthand Sides
%      using Block Principal Pivoting method
%
% This software solves the following problem: given A and B, find X such that
%              minimize || AX-B ||_F^2 where X>=0 elementwise.
%
% Reference:
%      Jingu Kim and Haesun Park. Fast Nonnegative Matrix Factorization: An Activeset-like Method and Comparisons,
%      SIAM Journal on Scientific Computing, 33(6), pp. 3261-3281, 2011.
%
% Written by Jingu Kim (jingu.kim@gmail.com)
%            School of Computational Science and Engineering,
%            Georgia Institute of Technology
%
% Note that this algorithm assumes that the input matrix A has full column rank.
% This code comes with no guarantee or warranty of any kind.
% Please send bug reports, comments, or questions to Jingu Kim.
%
% Modified Feb-20-2009
% Modified Mar-13-2011: numChol and numEq
%
% <Inputs>
%        A : input matrix (m x n) (by default), or A'*A (n x n) if isInputProd==1
%        B : input matrix (m x k) (by default), or A'*B (n x k) if isInputProd==1
%        isInputProd : (optional, default:0) if turned on, use (A'*A,A'*B) as input instead of (A,B)
%        init : (optional) initial value for X
% <Outputs>
%        X : the solution (n x k)
%        Y : A'*A*X - A'*B where X is the solution (n x k)
%        success : 0 for success, 1 for failure.
%                  Failure could only happen on a numericall very ill-conditioned problem.
%        numChol : number of unique cholesky decompositions done
%        numEqs : number of systems of linear equations solved

function [ X,Y,success,numChol,numEq ] = nnlsm_blockpivot( A, B, isInputProd, init )
    if nargin<3, isInputProd=0;, end
    if isInputProd
        AtA = A;, AtB = B;
    else
        AtA = A'*A;, AtB = A'*B;
    end

    if size(AtA,1)==1
        X = AtB/AtA; X(X<0) = 0;
        Y = AtA*X - AtB;
        numChol = 1; numEq = size(AtB,2); success = 1;
        return
    end

    [n,k]=size(AtB);
    MAX_BIG_ITER = n*5;
    % set initial feasible solution
    X = zeros(n,k);
    if nargin<4
        Y = - AtB;
        PassiveSet = false(n,k);
        numChol = 0;
        numEq = 0;
    else
        PassiveSet = (init > 0);
        [ X,numChol,numEq] = normalEqComb(AtA,AtB,PassiveSet);
        Y = AtA * X - AtB;
    end
    % parameters
    pbar = 3;
    P = zeros(1,k);, P(:) = pbar;
    Ninf = zeros(1,k);, Ninf(:) = n+1;

    NonOptSet = (Y < 0) & ~PassiveSet;
    InfeaSet = (X < 0) & PassiveSet;
    NotGood = sum(NonOptSet)+sum(InfeaSet);
    NotOptCols = NotGood > 0;

    bigIter = 0;, success=0;
    while(~isempty(find(NotOptCols)))
        bigIter = bigIter+1;
        if ((MAX_BIG_ITER >0) && (bigIter > MAX_BIG_ITER))   % set max_iter for ill-conditioned (numerically unstable) case
            success = 1;, break
        end

        Cols1 = NotOptCols & (NotGood < Ninf);
        Cols2 = NotOptCols & (NotGood >= Ninf) & (P >= 1);
        Cols3Ix = find(NotOptCols & ~Cols1 & ~Cols2);
        if ~isempty(find(Cols1))
            P(Cols1) = pbar;,Ninf(Cols1) = NotGood(Cols1);
            PassiveSet(NonOptSet & repmat(Cols1,n,1)) = true;
            PassiveSet(InfeaSet & repmat(Cols1,n,1)) = false;
        end
        if ~isempty(find(Cols2))
            P(Cols2) = P(Cols2)-1;
            PassiveSet(NonOptSet & repmat(Cols2,n,1)) = true;
            PassiveSet(InfeaSet & repmat(Cols2,n,1)) = false;
        end
        if ~isempty(Cols3Ix)
            for i=1:length(Cols3Ix)
                Ix = Cols3Ix(i);
                toChange = max(find( NonOptSet(:,Ix)|InfeaSet(:,Ix) ));
                if PassiveSet(toChange,Ix)
                    PassiveSet(toChange,Ix)=false;
                else
                    PassiveSet(toChange,Ix)=true;
                end
            end
        end
        [ X(:,NotOptCols),tempChol,tempEq ] = normalEqComb(AtA,AtB(:,NotOptCols),PassiveSet(:,NotOptCols));
        numChol = numChol + tempChol;
        numEq = numEq + tempEq;
        X(abs(X)<1e-12) = 0;            % One can uncomment this line for numerical stability.
        Y(:,NotOptCols) = AtA * X(:,NotOptCols) - AtB(:,NotOptCols);
        Y(abs(Y)<1e-12) = 0;            % One can uncomment this line for numerical stability.

        % check optimality
        NotOptMask = repmat(NotOptCols,n,1);
        NonOptSet = NotOptMask & (Y < 0) & ~PassiveSet;
        InfeaSet = NotOptMask & (X < 0) & PassiveSet;
        NotGood = sum(NonOptSet)+sum(InfeaSet);
        NotOptCols = NotGood > 0;
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function varargout = qr(f, ignored)
%QR Orthogonal-triangular decomposition of a SEPARABLEAPPROX.
%
% [Q, R] = QR( F ), where F is a separableApprox, produces an unitary column
% quasimatrix Q and a upper-triangular row quasimatrix R so that F = Q * R. This
% is computed by a continuous analogue of QR.
%
% [Q, R] = QR( F, 0 ) is the same as QR( F ).
%
% [Q, R, E] = QR( F ) and [Q, R, E] = QR( F, 'vector') produces a vector E that
% stores the pivoting locations.
%
% For more information about this decomposition:
% A. Townsend and L. N. Trefethen, Continuous analogues of matrix
% factorizations, Proc. Royal Soc. A., 2015.
%
% See also LU, and CHOL.

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

% Empty check:
if ( isempty( f ) )
    varargout = cell(1, nargout);
    return
end

% As always start with the CDR decomposition:
[C, D, R] = cdr( f );

% Balance out the scaling, becareful about signs:
sgns = sign( diag( D ) );
C = C * diag(sgns) * sqrt( abs(D) );   % Put the signs into Q
R = R * sqrt( abs(D) );

% QR of the column
[Q, RC] = qr( C );

% Form R:
R = RC * R.';   % still an upper-triangular quasimatrix! :)

% Output to user:
if ( nargout <= 1 )
    varargout = { Q };
elseif ( nargout == 2 )
    varargout = { Q, R };
elseif ( nargout == 3 )
    E = f.pivotLocations;   % complete pivot Locations from GE.
    varargout = { Q, R, E( :, 1 ) };  % Only pass the ones in x.
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fx = p07_fun ( n, x )

%*****************************************************************************80
%
%% P07_FUN evaluates the integrand for problem 7.
%
%  Discussion:
%
%    The integrand is singular at x = 0.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    1 / sqrt ( x )
%
%  Antiderivative:
%
%    2 * sqrt ( x )
%
%  Exact Integral:
%
%    2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner,
%    Comparison of Numerical Quadrature Formulas,
%    in Mathematical Software, edited by John R Rice,
%    Academic Press, 1971
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = 1.0 ./ sqrt ( x );
%
%  Replace "Inf" value at 0 by 0.
%
  i = find ( x == 0.0 );
  fx(i) = 0.0;

  return
end


text:
function[varargout]=mat2col(varargin)
%MAT2COL  Compress NAN-padded matrix data into long columns.
%
%   C=MAT2COL(M), where the columns of matrix M are data vectors of
%   nonuniform length with NANs filling out the gaps in each column,
%   appends together all the columns into a vector C in which isolated
%   NANs mark the transitions between columns.
%
%   [C1,C2,...,CN]=MAT2COL(M1,M2,...,MN) also works for multiple input
%   matrices of the same size.  In this case the locations of NANs in
%   M1 are used as the key for appending all the MN into columns.
%
%   MAT2COL, COL2MAT, and COLBREAKS together form a system for moving
%   data with segments of nonuniform length rapidly back and forth
%   between a column format and a padded-matrix format.
%
%   MAT2COL(M1,M2,...,MN); with no output arguments overwrites the
%   original input variables.
%
%   See also COL2MAT, COLBREAKS, COL2CELL, CELL2COL.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2016 J.M. Lilly --- type 'help jlab_license' for details

%no loops!

if strcmpi(varargin{1},'--t')
    return
end

if isempty(varargin{1})
   for i=1:nargout
       varargout{i}=[];
   end
   return
end

i=0;
while i<nargin
	i=i+1;
	mat=varargin{i};

	if i==1
		if ~all(isnan(mat(end,:)))
			mat(size(mat,1)+1,:)=nan*ones(size(mat(1,:)));
		end
		mat=mat(:);

		%new matrix should include all data points
		bool=false(size(mat));
		index=find(~isnan(mat));
		bool(index)=ones(size(index));
		index2=find(diff([bool;0])==-1)+1;
		if ~isempty(index2)
			bool(index2)=ones(size(index2));
		end
		index=find(bool);
	end
	mat=mat(index);
	col=mat;
	varargout{i}=col;
end

if nargout>nargin
        varargout{end+1}=index;
	%eval(['c' int2str(nargout) '=index;']);
end

if nargout==0
   %for i=1:nargin
   %   eval(['varargout{',int2str(i),'}=c',int2str(i),';'])
   %end
  eval(to_overwrite(nargin));
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [A]=comp2real(V)
%these function give vector of 'V' that have complex or non_complex number
%and separate real & imaginary part and attach those in one Matrix.
%this function put together real & imaginary part of any number in any row
m=length(V);
q=[];
for i=1:m
    q(i,1)=real(V(i));
    q(i,2)=imag(V(i));
end
A=q;

author: Sable

repo: mcbench-benchmarks


text:
function bw = Bw_Img(I)

if ndims(I) == 3
    I = rgb2gray(I);
end
bw = im2bw(I, graythresh(I));
bw = ~bw;

author: decouples


text:
function g = computeForwardKinematics(obj)
    % computes the forward kinematics transformation matrix of the
    % coordinate frame
    %
    % Return values:
    % g: the forward transformation matrix under the coordinate
    % system @type SymExpression

    frame = obj;
    while ~isprop(frame, 'TwistPairs') %isempty(frame.TwistPairs)
        frame = frame.Reference;
        if isempty(frame)
            error('The coordinate system is not fully defined.');
        end
    end

    twists = frame.TwistPairs;

    g = eval_math_fun('ForwardKinematics',[twists,{obj.gst0}]);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function y = fftn(varargin)
% For two inputs, if second input has a zero it gives an empty result with the
% size specified by the second input, like Matlab does, instead of giving
% an error
if nargin==2 && any(varargin{2}==0)
    y = ones(varargin{2});
else
    y = builtin('fftn', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function  [F_struc,K,binary_variables] = expandsemivar(F_struc,K,semicont_variables)

model.F_struc = F_struc;
model.K = K;
model.lb = -inf(size(model.F_struc,2)-1,1);
model.ub = -model.lb
model = presolve_bounds_from_modelbounds(model,1);

if any(isinf(model.lb(semicont_variables))) ||  any(isinf(model.ub(semicont_variables)))
    error('There are unbounded semi-continuous variables.');
end

% m new binaries required.
m = length(semicont_variables);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function obj = changeCaxis(obj)
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)

	userInput = inputdlg('CAXIS min max');str2num(userInput{1});
	S = findobj(gcf,'Type','Axes');
	% C = cell2mat(get(S,'Clim'));
	C = str2num(userInput{1});
	% C = [-1 7];
	set(S,'CLim',C);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function d = dot_outer(g1,g2,varargin)
% outer inner product between two quaternions
%
% Input
%  g1, g2 - @quaternion
%
% Output
%  d - double
%
% Description
% cos angle(g1,g2)/2 = dot(g1,g2)

if ~isempty(g1) && ~isempty(g2)

  q1 = [g1.a(:) g1.b(:) g1.c(:) g1.d(:)];
  q2 = [g2.a(:) g2.b(:) g2.c(:) g2.d(:)];

  d = q1 * q2.';

else
    d = [];
end

author: mtex-toolbox

repo: mtex

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


text:
clear all; close all;
N=0:4
X=fftshift(N)
Y=fftshift(fftshift(N))
Z=ifftshift(fftshift(N))

author: HuangCongQing


text:
function [data] = ft_megplanar(cfg, data)

% FT_MEGPLANAR computes planar MEG gradients gradients for raw data or average
% event-related field data. It can also convert frequency-domain data that was computed
% using FT_FREQANALYSIS, as long as it contains the complex-valued fourierspcrm and not
% only the powspctrm.
%
% Use as
%    [interp] = ft_megplanar(cfg, data)
% where the input data corresponds to the output from FT_PREPROCESSING,
% FT_TIMELOCKANALYSIS or FT_FREQANALYSIS (with output='fourier').
%
% The configuration should contain
%   cfg.planarmethod   = string, can be 'sincos', 'orig', 'fitplane', 'sourceproject' (default = 'sincos')
%   cfg.channel        =  Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.trials         = 'all' or a selection given as a 1xN vector (default = 'all')
%
% The methods orig, sincos and fitplane are all based on a neighbourhood interpolation.
% For these methods you need to specify
%   cfg.neighbours     = neighbourhood structure, see FT_PREPARE_NEIGHBOURS
%
% In the 'sourceproject' method a minumum current estimate is done using a large number
% of dipoles that are placed in the upper layer of the brain surface, followed by a
% forward computation towards a planar gradiometer array. This requires the
% specification of a volume conduction model of the head and of a source model. The
% 'sourceproject' method is not supported for frequency domain data.
%
% A dipole layer representing the brain surface must be specified with
%   cfg.inwardshift = depth of the source layer relative to the head model surface ,
%                     (default = 2.5 cm, which is appropriate for a skin-based head model)
%   cfg.spheremesh  = number of dipoles in the source layer (default = 642)
%   cfg.tolerance   = tolerance ratio for leadfield matrix inverse based on a truncated svd,
%                     reflects the relative magnitude of the largest singular value
%                     to retain (default = 1e-3)
%   cfg.headshape   = a filename containing headshape, a structure containing a
%                     single triangulated boundary, or a Nx3 matrix with surface
%                     points
% If no headshape is specified, the dipole layer will be based on the inner compartment
% of the volume conduction model.
%
% Optionally, you can modify the leadfields by reducing the rank, i.e. remove the weakest orientation
%   cfg.reducerank    = 'no', or number (default = 3 for EEG, 2 for MEG)
%   cfg.backproject   = 'yes' or 'no',  determines when reducerank is applied whether the
%                       lower rank leadfield is projected back onto the original linear
%                       subspace, or not (default = 'yes')
%
% The volume conduction model of the head should be specified as
%   cfg.headmodel     = structure with volume conduction model, see FT_PREPARE_HEADMODEL
%
% The following cfg fields are optional:
%   cfg.feedback
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_COMBINEPLANAR, FT_PREPARE_NEIGHBOURS

% Copyright (C) 2004-2019, Robert Oostenveld
% Copyright (C) 2020-      Robert Oostenveld and Jan-Mathijs Schoffelen
%
% 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$

% 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 data
ft_preamble provenance data

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

% store the original input representation of the data, this is used later on to convert it back
isfreq = ft_datatype(data, 'freq');
istlck = ft_datatype(data, 'timelock');  % this will be temporary converted into raw

% check if the input data is valid for this function, this converts the data if needed
data = ft_checkdata(data, 'datatype', {'raw' 'freq'}, 'feedback', 'yes', 'hassampleinfo', 'yes', 'ismeg', 'yes', 'senstype', {'ctf151', 'ctf275', 'bti148', 'bti248', 'itab153', 'yokogawa160', 'yokogawa64'});

if isfreq
  if ~isfield(data, 'fourierspctrm'), ft_error('freq data should contain Fourier spectra'); end
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels', 'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed',    {'hdmfile', 'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'vol',     'headmodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'grid',    'sourcemodel'});
cfg = ft_checkconfig(cfg, 'renamed',    {'pruneratio', 'tolerance'});

% set the default configuration
cfg.channel      = ft_getopt(cfg, 'channel',      'all');
cfg.tolerance    = ft_getopt(cfg, 'tolerance',    1e-3);
cfg.trials       = ft_getopt(cfg, 'trials',       'all', 1);
cfg.planarmethod = ft_getopt(cfg, 'planarmethod', 'sincos');
cfg.feedback     = ft_getopt(cfg, 'feedback',     'text');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamedval',  {'headshape', 'headmodel', []});

if ~strcmp(cfg.planarmethod, 'sourceproject')
  % this is limited to reading neighbours from disk and/or selecting channels
  % the user should call FT_PREPARE_NEIGHBOURS directly for the actual construction
  tmpcfg = keepfields(cfg, {'neighbours', 'channel', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  cfg.neighbours = ft_prepare_neighbours(tmpcfg);
end

% put the low-level options pertaining to the dipole grid in their own field
cfg = ft_checkconfig(cfg, 'renamed', {'tightgrid', 'tight'}); % this is moved to cfg.sourcemodel.tight by the subsequent createsubcfg
cfg = ft_checkconfig(cfg, 'renamed', {'sourceunits', 'unit'}); % this is moved to cfg.sourcemodel.unit by the subsequent createsubcfg

% put the low-level options pertaining to the sourcemodel in their own field
cfg = ft_checkconfig(cfg, 'createsubcfg', {'sourcemodel'});
% move some fields from cfg.sourcemodel back to the top-level configuration
cfg = ft_checkconfig(cfg, 'createtopcfg', {'sourcemodel'});

% select trials of interest
tmpcfg = keepfields(cfg, {'trials', 'channel', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'}); % don't keep tolerance, it is used differently here
data = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

if strcmp(cfg.planarmethod, 'sourceproject')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Do an inverse computation with a simplified distributed source model
  % and compute forward again with the axial gradiometer array replaced by
  % a planar one.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % method specific configuration options
  cfg.headshape   = ft_getopt(cfg, 'headshape',   []);
  cfg.inwardshift = ft_getopt(cfg, 'inwardshift', 2.5); % this number assumes that all other inputs are in cm
  cfg.pruneratio  = ft_getopt(cfg, 'pruneratio',  1e-3);
  cfg.spheremesh  = ft_getopt(cfg, 'spheremesh',  642);

  if isfreq
    ft_error('the method ''sourceproject'' is not supported for frequency data as input');
  end

  % PREPARE_HEADMODEL will match the data labels, the gradiometer labels and the
  % volume model labels (in case of a localspheres model) and result in a gradiometer
  % definition that only contains the gradiometers that are present in the
  % data. This should exclude the non-MEG channels, so the user-defined
  % cfg.channel should be overruled
  tmpcfg = cfg;
  tmpcfg.channel = ft_channelselection('MEG', cfg.channel);
  [headmodel, axial.grad, tmpcfg] = prepare_headmodel(tmpcfg, data);

  % construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
  leadfieldopt = {};
  leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(cfg, 'reducerank'));
  leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(cfg, 'backproject'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(cfg, 'normalize'));
  leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(cfg, 'normalizeparam'));
  leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(cfg, 'weight'));

  % copy all options that are potentially used in FT_PREPARE_SOURCEMODEL
  tmpcfg           = keepfields(cfg, {'sourcemodel', 'mri', 'headshape', 'symmetry', 'smooth', 'threshold', 'spheremesh', 'inwardshift', 'xgrid' 'ygrid', 'zgrid', 'resolution', 'tight', 'warpmni', 'template', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
  tmpcfg.headmodel = headmodel;
  tmpcfg.grad      = axial.grad;
  % determine the dipole layer that represents the surface of the brain
  sourcemodel = ft_prepare_sourcemodel(tmpcfg);

  % compute the forward model for the axial gradiometers
  ft_info('computing forward model for %d dipoles\n', size(sourcemodel.pos,1));
  lfold = ft_compute_leadfield(sourcemodel.pos, axial.grad, headmodel, leadfieldopt{:});

  % construct the planar gradient definition and compute its forward model
  % this will not work for a localspheres model, compute_leadfield will catch
  % the error
  planar.grad = constructplanargrad([], axial.grad);
  lfnew = ft_compute_leadfield(sourcemodel.pos, planar.grad, headmodel, leadfieldopt{:});

  % compute the interpolation matrix
  transform = lfnew * ft_inv(lfold, 'method', 'tsvd', 'tolerance', cfg.tolerance);

  planarmontage = [];
  planarmontage.tra      = transform;
  planarmontage.labelold = axial.grad.label;
  planarmontage.labelnew = planar.grad.label;

  % apply the linear transformation to the data
  interp  = ft_apply_montage(data, planarmontage, 'keepunused', 'yes');

  % also apply the linear transformation to the gradiometer definition
  interp.grad = ft_apply_montage(data.grad, planarmontage, 'balancename', 'planar', 'keepunused', 'yes');

  % ensure there is a type string describing the gradiometer definition
  if ~isfield(interp.grad, 'type')
    interp.grad.type = [ft_senstype(data.grad) '_planar'];
  else
    interp.grad.type = [interp.grad.type '_planar'];
  end

  %   % interpolate the data towards the planar gradiometers
  %   for i=1:Ntrials
  %     ft_info('interpolating trial %d to planar gradiometer\n', i);
  %     interp.trial{i} = transform * data.trial{i}(dataindx,:);
  %   end % for Ntrials
  %
  %   % all planar gradiometer channels are included in the output
  %   interp.grad  = planar.grad;
  %   interp.label = planar.grad.label;
  %
  %   % copy the non-gradiometer channels back into the output data
  %   other = setdiff(1:Nchan, dataindx);
  %   for i=other
  %     interp.label{end+1} = data.label{i};
  %     for j=1:Ntrials
  %       interp.trial{j}(end+1,:) = data.trial{j}(i,:);
  %     end
  %   end
  %
else

  sens = ft_determine_units(data.grad);
  chanposnans = any(isnan(sens.chanpos(:))) || any(isnan(sens.chanori(:)));
  if chanposnans
    if isfield(sens, 'chanposold')
      % temporarily replace chanpos and chanorig with the original values
      sens.chanpos = sens.chanposold;
      sens.chanori = sens.chanoriold;
      sens.label   = sens.labelold;
      sens = rmfield(sens, {'chanposold', 'chanoriold', 'labelold'});
    else
      ft_error('The channel positions (and/or orientations) contain NaNs; this prohibits correct behavior of the function. Please replace the input channel definition with one that contains valid channel positions');
    end
  end
  cfg.channel = ft_channelselection(cfg.channel, sens.label);
  cfg.channel = ft_channelselection(cfg.channel, data.label);

  % ensure channel order according to cfg.channel (there might be one check
  % too much in here somewhere or in the subfunctions, but I don't care.
  % Better one too much than one too little - JMH @ 09/19/12
  cfg = struct(cfg);
  [neighbsel] = match_str({cfg.neighbours.label}, cfg.channel);
  cfg.neighbours = cfg.neighbours(neighbsel);
  cfg.neighbsel = channelconnectivity(cfg);

  assert(any(cfg.neighbsel(:)), 'no neighbours found')
  ft_info('average number of neighbours is %.2f\n', mean(sum(cfg.neighbsel)));

  Ngrad = length(sens.label);
  distance = zeros(Ngrad,Ngrad);

  for i=1:size(cfg.neighbsel,1)
    j=find(cfg.neighbsel(i, :));
    d = sqrt(sum((sens.chanpos(j,:) - repmat(sens.chanpos(i, :), numel(j), 1)).^2, 2));
    distance(i,j) = d;
    distance(j,i) = d;
  end

  ft_info('minimum distance between neighbours is %6.2f %s\n', min(distance(distance~=0)), sens.unit);
  ft_info('maximum distance between gradiometers is %6.2f %s\n', max(distance(distance~=0)), sens.unit);

  % The following does not work when running in deployed mode because the
  % private functions that compute the planar montage are not recognized as
  % such and won't be compiled, unless explicitly specified.

  % % generically call megplanar_orig megplanar_sincos or megplanar_fitplane
  %fun = ['megplanar_'  cfg.planarmethod];
  %if ~exist(fun, 'file')
  %  ft_error('unknown method for computation of planar gradient');
  %end
  %planarmontage = eval([fun '(cfg, data.grad)']);

  switch cfg.planarmethod
    case 'sincos'
      planarmontage = megplanar_sincos(cfg, sens);
    case 'orig'
      % method specific info that is needed
      cfg.distance  = distance;

      planarmontage = megplanar_orig(cfg, sens);
    case 'fitplane'
      planarmontage = megplanar_fitplane(cfg, sens);
    otherwise
      fun = ['megplanar_' cfg.planarmethod];
      if ~exist(fun, 'file')
        ft_error('unknown method for computation of planar gradient');
      end
      planarmontage = eval([fun '(cfg, data.grad)']);
  end

  % apply the linear transformation to the data
  interp = ft_apply_montage(data, planarmontage, 'keepunused', 'yes', 'feedback', cfg.feedback);

  % also apply the linear transformation to the gradiometer definition
  interp.grad = ft_apply_montage(sens, planarmontage, 'balancename', 'planar', 'keepunused', 'yes');

  % ensure there is a type string describing the gradiometer definition
  if ~isfield(interp.grad, 'type')
    % put the original gradiometer type in (will get _planar appended)
    interp.grad.type = ft_senstype(sens);
  end
  interp.grad.type = [interp.grad.type '_planar'];

  % add the chanpos info back into the gradiometer description
  tmplabel = interp.grad.label;
  for k = 1:numel(tmplabel)
    if ~isempty(strfind(tmplabel{k}, '_dV')) || ~isempty(strfind(tmplabel{k}, '_dH'))
      tmplabel{k} = tmplabel{k}(1:end-3);
    end
  end
  [ix,iy] = match_str(tmplabel, sens.label);
  interp.grad.chanpos(ix,:) = sens.chanpos(iy,:);

  % if the original chanpos contained nans, make sure to put nans in the
  % updated one as well, and move the updated chanpos values to chanposold
  if chanposnans
    interp.grad.chanposold = sens.chanpos;
    interp.grad.chanoriold = sens.chanori;
    interp.grad.labelold   = sens.label;
    interp.grad.chanpos    = nan(size(interp.grad.chanpos));
    interp.grad.chanori    = nan(size(interp.grad.chanori));
  end
end

if istlck
  % convert the raw structure back into a timelock structure
  interp = ft_checkdata(interp, 'datatype', 'timelock');
end

% copy the trial specific information into the output
if isfield(data, 'trialinfo')
  interp.trialinfo = data.trialinfo;
end

% copy the sampleinfo field as well
if isfield(data, 'sampleinfo')
  interp.sampleinfo = data.sampleinfo;
end

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

% rename the output variable to accomodate the savevar postamble
data = interp;

ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data


text:
%  The soacm code is modified on the basis of an implementation of C-V model, which is:
%  created on 04/26/2004
%  author: Chunming Li
%  email: li_chunming@hotmail.com
%  Copyright (c) 2004-2006 by Chunming Li
function phi = evol_soacm(ort,img,phi_init,mu,nu,lambda_1,lambda_2,delta_t,epsilon)

    img = mirror_expand(img);

    phi = mirror_expand(phi_init);
    phi = mirror_ensure(phi);
    delta_h = Delta(phi,epsilon);
    cur = curvature(phi); % compute curvature of phi function
    [gray_c1,gray_c2] = binaryfit(phi,img,epsilon);

    if ort == -1
        t1 = lambda_1 * (img-gray_c1).^2;
        t2 = lambda_2 * (img-gray_c2).^2;
        phi = phi + delta_t * delta_h .* ( mu*cur - nu - t1 + t2);
    else
        t3 = (img-ort).^2;
        phi = phi + delta_t * delta_h .* ( mu*cur - nu - t3);
    end

    phi = mirror_shrink(phi);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function enforceFields(structure, field_names)
%ENFORCEFIELDS   Checks structure field names for existance.
%
% DESCRIPTION:
%       enforceFields checks a MATLAB structure for the existance of a set
%       of required field name defined by a cell array. If a field name
%       within the cell array is not found within the structure, an error
%       is thrown with the name of the missing field given in the error
%       message. The function is useful for enforcing a set of required
%       fields within a structure.
%
% USAGE:
%       enforceFields(structure, field_names)
%
% INPUTS:
%       structure   - MATLAB structure
%       field_names - cell array of allowable field names
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 13th October 2009
%       last update - 13th October 2009
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also fieldnames

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

names = fieldnames(structure);
for field_index = 1:length(field_names)
    field_found = false;
    for names_index = 1:length(names)
        if strcmp(field_names{field_index}, names{names_index})
            field_found = true;
        end
    end
    if ~field_found
    	error(['The field ' inputname(1) '.' field_names{field_index} ' must be defined']);
    end
end

repo: Aero-matlab


text:
function a = r8vec_indicator ( n )

%*****************************************************************************80
%
%% R8VEC_INDICATOR sets an R8VEC to the indicator vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Output, real A(N), the vector with entries (1, 2, ..., N ).
%
  a = ( ( 1 : n ) )';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [B0,BV] = spm_MDP_DP(MDP)
% dynamic programming using active inference
% FORMAT [B0,BV] = spm_MDP_DP(MDP)
%
% MDP.A(O,N)      - Likelihood of O outcomes given N hidden states
% MDP.B{M}(N,N)   - transition probabilities among hidden states (priors)
% MDP.C(N,1)      - prior preferences (prior over future states)
%
% MDP.V(T - 1,P)  - P allowable policies (control sequences)
%
% B0      - optimal state action policy or transition matrix
% BV      - corresponding policy using value iteration
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_DP.m 6598 2015-11-11 19:48:30Z karl $

% set up and preliminaries
%==========================================================================

% generative model and initial states
%--------------------------------------------------------------------------
T     = size(MDP.V,1) + 1;        % number of outcomes
Ns    = size(MDP.B{1},1);         % number of hidden states
Nu    = size(MDP.B,2);            % number of hidden controls
p0    = exp(-8);                  % smallest probability

% likelihood model (for a partially observed MDP implicit in G)
%--------------------------------------------------------------------------
try
    A = MDP.A + p0;
catch
    A = speye(Ns,Ns) + p0;
end
A     = A*diag(1./sum(A));        % normalise
lnA   = log(A);                   % log probabilities
H     = sum(A.*lnA)';             % negentropy of observations

% transition probabilities (priors)
%--------------------------------------------------------------------------
for j = 1:Nu
    B{j}   = MDP.B{j} + p0;
    B{j}   = B{j}*diag(1./sum(B{j}));
    sB{j}  = B{j};
    rB{j}  = spm_softmax(log(B{j})');
    lnB{j} = log(B{j});
end

% terminal probabilities over outcomes (priors)
%--------------------------------------------------------------------------
try
    C = MDP.C;
catch
    C = zeros(No,1);
end

% asume constant preferences over states
%--------------------------------------------------------------------------
if size(C,2) ~= T
    C = C(:,end)*ones(1,T);
end
C = A'*diag(1./sum(A,2))*spm_softmax(C);
C = log(C);

% policies, states and their expectations
%--------------------------------------------------------------------------
V     = MDP.V;
Np    = size(V,2);                % number of allowable policies

% policy iteration
%==========================================================================
for s = 1:Ns

    % Variational iterations (hidden states)
    %======================================================================
    x     = zeros(Ns,T,Np) + 1/Ns;
    for k = 1:Np

        % gradient descent on free energy
        %------------------------------------------------------------------
        for i = 1:16

            % hiddens states (x)
            %--------------------------------------------------------------
            x(:,1,k) = 0;
            x(s,1,k) = 1;

            for j = 2:T

                % current state
                %----------------------------------------------------------
                xj   = x(:,j,k);
                qx   = log(xj);
                v    = 0;

                % evaluate free energy and gradients (v = dFdx)
                %----------------------------------------------------------
                if j > 1, v = v + qx - log(sB{V(j - 1,k)}*x(:,j - 1,k)); end
                if j < T, v = v      - log(rB{V(j    ,k)}*x(:,j + 1,k)); end

                % update
                %----------------------------------------------------------
                x(:,j,k) = spm_softmax(qx - v/4);
                F(j,k)   = xj'*v;

            end

            % convergence
            %--------------------------------------------------------------
            if i > 1
                dF = F0 - sum(F(:,k));
                if dF > 1/128, F0 = F0 - dF; else, break, end
            else
                F0 = sum(F(:,k));
            end

        end
    end

    % value of policies (Q)
    %======================================================================
    Q     = zeros(Np,1);
    for k = 1:Np

        % path integral of expected free energy
        %------------------------------------------------------------------
        for j = 2:T
            v    = C(:,j) - log(x(:,j,k)) + H;
            Q(k) = Q(k) + v'*x(:,j,k);

        end
    end

    % optimal transition from this state
    %======================================================================
    [u,k]   = max(Q);
    B0(:,s) = lnB{V(1,k)}(:,s);

end

if nargout < 2, return, end

% value iteration
%==========================================================================
V     = zeros(Ns,1);
C     = C(:,end);
g     = 1 - 1/T;
for i = 1:32
    for s  = 1:Ns

        % value of actions (Q)
        %------------------------------------------------------------------
        Q     = zeros(Nu,1);
        for k = 1:Nu
            Q(k)   = B{1,k}(:,s)'*(C + g*V);
        end

        % optimal transition from this state
        %------------------------------------------------------------------
        [u,k]   = max(Q);
        BV(:,s) = B{1,k}(:,s);

    end

    % optimal transition from this state
    %----------------------------------------------------------------------
    dV   = BV'*(C + g*V) - V;
    V    = V + dV;

    % convergence
    %----------------------------------------------------------------------
    if norm(dV) < 1e-2, break, end

end

author: spm


text:
%COLORNAME Map between color names and RGB values
%
% RGB = COLORNAME(NAME) is the RGB-tristimulus value (1x3) corresponding to
% the color specified by the string NAME.  If RGB is a cell-array (1xN) of
% names then RGB is a matrix (Nx3) with each row being the corresponding
% tristimulus.
%
% XYZ = COLORNAME(NAME, 'xyz') as above but the XYZ-tristimulus value
% corresponding to the color specified by the string NAME.
%
% XY = COLORNAME(NAME, 'xy') as above but the xy-chromaticity coordinates
% corresponding to the color specified by the string NAME.
%
% NAME = COLORNAME(RGB) is a string giving the name of the color that is
% closest (Euclidean) to the given RGB-tristimulus value (1x3).  If RGB is
% a matrix (Nx3) then return a cell-array (1xN) of color names.
%
% NAME = COLORNAME(XYZ, 'xyz') as above but the color is the closest (Euclidean)
% to the given XYZ-tristimulus value.
%
% NAME = COLORNAME(XYZ, 'xy') as above but the color is the closest (Euclidean)
% to the given xy-chromaticity value with assumed Y=1.
%
% Notes::
% - Color name may contain a wildcard, eg. "?burnt"
% - Based on the standard X11 color database rgb.txt.
% - Tristimulus values are in the range 0 to 1

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function r = colorname(a, varargin)

    opt.color = {'rgb', 'xyz', 'xy'};
    opt = tb_optparse(opt, varargin);

    persistent  rgbtable;

    mvtb_present = exist('tristim2cc');

    % ensure that the database is loaded
    if isempty(rgbtable)
        % load mapping table from file
        fprintf('loading rgb.txt\n');
        f = fopen('private/rgb.txt', 'r');
        k = 0;
        rgb = [];
        names = {};
        xy = [];

        while ~feof(f),
            line = fgets(f);
            if line(1) == '#',
                continue;
            end

            [A,count,errm,next] = sscanf(line, '%d %d %d');
            if count == 3
                k = k + 1;
                rgb(k,:) = A' / 255.0;
                names{k} = lower( strtrim(line(next:end)) );
                if mvtb_present
                    xy = tristim2cc( colorspace('RGB->XYZ', rgb) );
                end
            end
        end
        s.rgb = rgb;
        s.names = names;
        if mvtb_present
            s.xy = xy;
        end
        rgbtable = s;
    end

    if isstr(a)
        % map name to rgb/xy
        if a(1)  == '?'
            % just do a wildcard lookup
            r = namelookup(rgbtable, a(2:end), opt);
        else
            r = name2rgb(rgbtable, a, opt);
        end
    elseif iscell(a)
        % map multiple names to rgb
        r = [];
        for name=a,
            rgb = name2rgb(rgbtable, name{1}, opt.xy);
            if isempty(rgb)
                warning('Color %s not found', name{1});
            end
            r = [r; rgb];
        end
    else
        % map values to strings
        switch opt.color
            case 'rgb'
                if numel(a) == 3
                    r = rgb2name(rgbtable, a(:)');
                elseif numcols(a) ~= 3
                    error('RGB data must have 3 columns');
                else
                    r = {};
                    for i=1:numrows(a)
                        r{i} = rgb2name(rgbtable, a(i,:));
                    end
                end

            case 'xyz'
                if numel(a) == 3
                    rgb = colorspace('XYZ->RGB', a(:)');
                    r = rgb2name(rgbtable, rgb);
                elseif numcols(a) ~= 3
                    error('XYZ data must have 3 columns');
                else
                    rgb = colorspace('XYZ->RGB', a);

                    r = {};
                    for i=1:numrows(a)
                        r{i} = rgb2name(rgbtable, rgb(i,:));
                    end
                end

            case 'xy'
                if numel(a) == 2
                    Y = 1;  XYZ = 1/a(2);
                    X = a(1) * XYZ;
                    Z = (1-a(1)-a(2)) * XYZ;
                    rgb = colorspace('XYZ->RGB', [X Y Z]);
                    r = rgb2name(rgbtable, rgb);
                elseif numcols(a) ~= 2
                    error('xy data must have 2 columns');
                else
                    Y = ones(numrows(a),1);  XYZ = 1./a(:,2);
                    X = a(:,1) .* XYZ;
                    Z = (1-a(:,1)-a(:,2)) .* XYZ;
                    rgb = colorspace('XYZ->RGB', [X Y Z]);

                    r = {};
                    for i=1:numrows(a)
                        r{i} = rgb2name(rgbtable, rgb(i,:));
                    end
                end
        end

    end
end

function r = namelookup(table, s)
    s = lower(s);   % all matching done in lower case

    r = {};
    count = 1;
    for k=1:length(table.names),
        if ~isempty( findstr(table.names{k}, s) )
            r{count} = table.names{k};
            count = count + 1;
        end
    end
end

function out = name2rgb(table, s, opt)

    s = lower(s);   % all matching done in lower case

    for k=1:length(table.names),
        if strcmp(s, table.names(k)),
            rgb = table.rgb(k,:);
            switch opt.color
                case 'rgb'
                    out = rgb;
                case 'xy'
                    XYZ = colorspace('RGB->XYZ', r);
                    out = tristim2cc(XYZ);
                case 'xyz'
                    out = colorspace('RGB->XYZ', rgb);
            end
            return;
        end
    end
    out = [];
end

function r = rgb2name(table, v)
    d = table.rgb - ones(numrows(table.rgb),1) * v;
    n = colnorm(d');
    [z,k] = min(n);
    r = table.names{k};
end

function r = xy2name(table, v)
    d = table.xy - ones(numrows(table.xy),1) * v;
    n = colnorm(d');
    [z,k] = min(n);
    r = table.names{k};
end

author: Allopart

repo: rbpf-gmapping


text:
function drawCartPoleTraj(t,p1,p2,nFrame)
% drawCartPoleTraj(t,p1,p2,nFrame)
%
% INPUTS:
%   t =  [1,n] = time stamp for the data in p1 and p2
%   p1 = [2,n] = [x;y] = position of center of the cart
%   p2 = [2,n] = [x;y] = position of tip of the pendulum
%   nFrame = scalar integer = number of "freeze" frames to display
%

clf; hold on;

Cart_Width = 0.15;
Cart_Height = 0.05;

Pole_Width = 4;  %pixels

%%%% Figure out the window size:

[xLow, xUpp, yLow, yUpp] = getBounds(p1,p2);

xLow = xLow - 0.7*Cart_Width;
xUpp = xUpp + 0.7*Cart_Width;

yLow = yLow - 0.7*Cart_Height;
yUpp = yUpp + 0.7*Cart_Height;

Limits = [xLow,xUpp,yLow,yUpp];

%%%% Get color map for the figure
map = colormap;
tMap = linspace(t(1),t(end),size(map,1))';

%%%% Plot Rails
plot([Limits(1) Limits(2)],-0.5*Cart_Height*[1,1],'k-','LineWidth',2)

%%%% Draw the trace of the pendulum tip  (continuously vary color)

nTime = length(t);
for i=1:(nTime-1)
    idx = i:(i+1);
    x = p2(1,idx);
    y = p2(2,idx);
    c = interp1(tMap,map,mean(t(idx)));
    plot(x,y,'Color',c);
end

%%%% Compute the frames for plotting:
tFrame = linspace(t(1), t(end), nFrame);
cart = interp1(t',p1',tFrame')';
pole = interp1(t',p2',tFrame')';

for i = 1:nFrame

    % Compute color:
    color = interp1(tMap,map,tFrame(i));

    %Plot Cart
    x = cart(1,i) - 0.5*Cart_Width;
    y = -0.5*Cart_Height;
    w = Cart_Width;
    h = Cart_Height;
    hCart = rectangle('Position',[x,y,w,h],'LineWidth',2);
    set(hCart,'FaceColor',color);
    set(hCart,'EdgeColor',0.8*color);

    %Plot Pendulum
    Rod_X = [cart(1,i), pole(1,i)];
    Rod_Y = [cart(2,i), pole(2,i)];
    plot(Rod_X,Rod_Y,'k-','LineWidth',Pole_Width,'Color',color)

    %Plot Bob and hinge
    plot(pole(1,i),pole(2,i),'k.','MarkerSize',40,'Color',color)
    plot(cart(1,i),cart(2,i),'k.','MarkerSize',60,'Color',color)

end

%These commands keep the window from automatically rescaling in funny ways.
axis(Limits);
axis('equal');
axis manual;
axis off;

end

function [xLow, xUpp, yLow, yUpp] = getBounds(p1,p2)
%
% Returns the upper and lower bound on the data in val
%

val = [p1,p2];
xLow = min(val(1,:));
xUpp = max(val(1,:));
yLow = min(val(2,:));
yUpp = max(val(2,:));

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [hd_estimates] = hd_new_for_signres(const,exo,beta,k,n,p,D,m,T,X,Y,data_exo,IRFt,signreslabels_shocks)
%computes the historical decomposition for the series
%ouput hd_estimates = cell array where columns capture variables, and rows
%the contributions of shocks, exogenous, constant and initial conditions to
%these variables. The second to last rows capture the unexplained part for models
%that are not fully identified, while the last row captures the part of
%the fluctuation that should be explained by the VAR, after accounting for
%exogenous/deterministic components.

%row 1 to n = contribution of shock x the movement in variable y hd_estimates(x,y)
%row n+1 = contribution of the constant
%row n+2 = contribution of initial conditions (past shocks)
%row n+3 = unexplained part (for partially identified model
%row n+4 = part that was left to explain by the structural shocks after
%accounting for exogenous, constant and initial conditions
%% Preliminaries

% preliminaries for historical decomposition
%1. Determine how many contributions we are going to calculate
contributors = n + 1 + 1 + length(exo); %variables + constant + initial conditions + exogenous
hd_estimates2=cell(contributors+2,n); %shocks+constant+initial values+exogenous+unexplained+to be explained by shocks only
% number of identified shocks
if IRFt==2 || IRFt==3
    identified=n; % fully identified
elseif IRFt==4 || IRFt==6 %if the model is identified by sign restrictions or sign restrictions (+ IV)
    identified=size(signreslabels_shocks,1); % count the labels provided in the sign res sheet (+ IV)
elseif IRFt==5
    identified=1; % one IV shock
elseif IRFt==1
    identified=0;
end

%===============================================
Bfull=reshape(beta,k,n);                    %get the Bfull matrix
B=Bfull(1:n*p,:);                           %drop the coefficients for all exogenous variables from the matrix
Bcomp = [B'; eye(n*(p-1)) zeros(n*(p-1),n)];%put into companion form
EPS=Y-X*Bfull;                               %get reduced form residuals
ETA=(D\EPS');                        %get structural shocks

%% Compute historical decompositions
%===============================================
% Contribution of each shock
    aux_D = zeros(n*p,n); %auxilary D matrix, that is consistent with companion matrix
    aux_D(1:n,:) = D; %set the first entrys equal to Dinv
    Selec = [eye(n) zeros(n,(p-1)*n)]; %selection matrix
    HDestimates_store = zeros(p*n,T+1,n); %cell aray to store results
    HDestimates = zeros(n,T+1,n);
    for j=1:n % for each variable
        ETAcomp = zeros(n,T+1); %structural shock matrix that is consistent with companion form
        if j <= identified %if j is an identified shock
        ETAcomp(j,2:end) = ETA(j,:); %fill in the entry for shock n, leave the first entry blank
        end
        for i = 2:T+1
            HDestimates_store(:,i,j) = aux_D*ETAcomp(:,i) + Bcomp*HDestimates_store(:,i-1,j); %recursively sum over shock impulse at period i on variable j and the previous period
            HDestimates(:,i,j) =  Selec*HDestimates_store(:,i,j); %select the entry corresponding to the current period
        end
    end

% contribution of the initial values
    HDinitial_storage   = zeros(p*n,T+1);
    HDinitial_estimates = zeros(n, T+1);
    Xnoexo = X(:,1:n*p);
    HDinitial_storage(:,1) = Xnoexo(1,:)'; %set the initial values to the first row of X (n*P)+exo
    HDinitial_estimates(:,1) = Selec*HDinitial_storage(:,1); %select the initial values for the first n variables (i.e. the values at Y_{t-1}
    for i = 2:T+1 %loop over periods and compute the impact of the initial conditions recursively
        HDinitial_storage(:,i) = Bcomp*HDinitial_storage(:,i-1); %compute the impact of those values (which in principle consist of past shocks)
        HDinitial_estimates(:,i) = Selec*HDinitial_storage(:,i);
    end

%  Contribution of the Constant
    HDconstant_storage = zeros(p*n,T+1);
    HDconstant_estimates = zeros(n, T+1);
    Coefficients = zeros(p*n,1);
    if const==1
        Coefficients(1:n,:) = Bfull(n*p+1,:);
        for i = 2:T+1 %loop over periods
            HDconstant_storage(:,i) = Coefficients + Bcomp*HDconstant_storage(:,i-1);
            HDconstant_estimates(:,i) = Selec * HDconstant_storage(:,i);
        end
    end

 % Contribution of exogenous variables, no lags of exogenous variables just
 % the contemporaneous values
 if m > 1
    HDexo_storage = zeros(p*n,T+1);
    HDexo_estimates = zeros(n,T+1);
%     Coefficients_exo = zeros(p*n,(m-const)*(1+p));
    Coefficients_exo = zeros(p*n,m-const);
    data_exocut=X(:,n*p+const+1:end); % take exogenous data from X
%     Coefficients_exo(1:n,:) = Bfull(n*p+const+1:end,:)'; %get the corresponding coefficients
    Coefficients_exo(1:n,:) = Bfull(n*p+const+1:end,:)'; %get the corresponding coefficients
        for i = 2:T+1
            HDexo_storage(:,i) = Coefficients_exo*data_exocut(i-1,:)' + Bcomp*HDexo_storage(:,i-1);
            HDexo_estimates(:,i) = Selec * HDexo_storage(:,i);
        end
  end

 %% put these values into the corresponding cell for hd_estimates such that
 % for variable x (hd_estimates(x,n+1)) = HDinitial_estimates(x,:)
 % for variable x (hd_estimates(x,n+1)) = HDconstant_estimates(x,:)
%reorganize storage
        for jj=1:n %for variables
            for kk=1:T+1 %for periods
                for ii=1:n %for shock contributions
                    hd_estimates2{ii,jj}(1,kk) = HDestimates(jj,kk,ii);
                end
                hd_estimates2{n+1,jj}(1,kk)=HDinitial_estimates(jj,kk);
                if const==1
                hd_estimates2{n+2,jj}(1,kk)=HDconstant_estimates(jj,kk);
                else
                hd_estimates2{n+2,jj}(1,kk)=0;
                end
                if m>1
                hd_estimates2{n+3,jj}(1,kk)=HDexo_estimates(jj,kk);
                end
            end
        end

 HDsum = zeros(T+1,n); %if we sum over all variables this should give Y
 for jj=1:n %loop over variables (columns)
 sumvariable = zeros(1,T+1);
 for kk=1:T+1 %loop over periods
     sumperiod=0;
  for ii=1:contributors %loop over contributors (rows)
      value = hd_estimates2{ii,jj}(1,kk);
      sumperiod = sumperiod+value;
  end
  sumvariable(1,kk) = sumperiod;
 end
 HDsum(:,jj)=sumvariable(1,:);
 end

 %determine the unexplained part (if model is not fully identified)
 aux = zeros(1,n);
 unexplained = Y-HDsum(2:end,:);
 unexplained = [aux; unexplained];
for jj=1:n
    hd_estimates2{contributors+1,jj}=unexplained(:,jj)';
end

%finally substract the sum of the contribution of the
%exogenous, constant,initial conditions from Y to get the
%part that was left to be explained by the shocks (for plotting reasons)
 Exosum = zeros(T+1,n); %if we sum over all variables this should give Y
 for jj=1:n %loop over variables (columns)
 sumvariable = zeros(1,T+1);
 for kk=1:T+1 %loop over periods
     sumperiod=0;
  for ii=n+1:contributors %loop over contributors (rows)
      value = hd_estimates2{ii,jj}(1,kk);
      sumperiod = sumperiod+value;
  end
  sumvariable(1,kk) = sumperiod;
 end
 Exosum(:,jj)=sumvariable(1,:);
 end

 %determine the part that was left to be explained by the shocks
 aux = zeros(1,n);
 tobeexplained = Y - Exosum(2:end,:);
 tobeexplained = [aux; tobeexplained];
for jj=1:n
    hd_estimates2{contributors+2,jj}=tobeexplained(:,jj)';
end

hd_estimates = cell(contributors+2,n);

%drop the initial entry for each cell in hd_estimates
for jj=1:n
    for ii=1:contributors+2
        hd_estimates{ii,jj}(2,:)=hd_estimates2{ii,jj}(2:end);
    end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [i,j,k] = PickRandomLatticeSite_3D_QPOTTS(x,y,z)
Sx = size(x,1);
i=floor(1+rand*Sx);
j=floor(1+rand*Sx);
k=floor(1+rand*Sx);

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34985-monte-carlo-simulation-of-three-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 3D square-lattice - microstructure/PickRandomLatticeSite_3D_QPOTTS.m


text:
function data_list = read_data_list(base_path, classes, data_size, train_or_test, debug)
% read pre-computed volumetric filenames
% base_path: root data folder
% classes: categories to be fetched
% data_size: the size of the volumetric representation
% train_or_test: select training data or testing data('train', 'test')
% debug: true for debug.

if nargin < 5
    debug = 0;
end

if debug
    maxNum = 20;
elseif strcmp(train_or_test, 'train')
    maxNum = 80*12;
else
    maxNum = 20*12;
end

num_classes = length(classes);
data_list = cell(num_classes,1);
for c = 1 : num_classes
    fprintf('reading the %s category\n', classes{c});
    category_path = [base_path '/' classes{c} '/' num2str(data_size) '/' train_or_test];
    files = dir(category_path);

    cat_ind = 0;
    data_list{c} = [];
    for i = 1 : length(files)
        idx = i;
        if strcmp(files(idx).name, '.') || strcmp(files(idx).name, '..') || files(idx).isdir == 1 || ~strcmp(files(idx).name(end-2:end), 'mat')
            continue;
        end
        filename = [category_path '/' files(idx).name];

        cat_ind = cat_ind + 1;
        data_list{c}(cat_ind,1).filename = filename;
        data_list{c}(cat_ind,1).label = c;
        if cat_ind == maxNum
          break;
        end
    end
    fprintf('fetched %d number of instances\n', cat_ind);
end

author: zhirongw


text:
function p = gaussProb(X, mu, Sigma)
% Multivariate Gaussian distribution, pdf
% X(i,:) is i'th case
% *** In the univariate case, Sigma is the variance, not the standard
% deviation! ***

% This file is from pmtk3.googlecode.com

d = size(Sigma, 2);
X  = reshape(X, [], d);  % make sure X is n-by-d and not d-by-n
X = bsxfun(@minus, X, rowvec(mu));
logp = -0.5*sum((X/(Sigma)).*X, 2);
logZ = (d/2)*log(2*pi) + 0.5*logdet(Sigma);
logp = logp - logZ;
p = exp(logp);

end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function angVal=angleBetweenLineAndPlane(a,n)
%%ANGLEBETWEENLINEANDPLANE Given a line specified in parameteric form as
%       x=x0+a*t, where t is the scalar parametric parameter, and a plane
%       given in parameteric form as x'*n=d, return the (acute) angle
%       between the line and the plane.
%
%INPUTS: a The xDimX1 parameteric slope vector of the line.
%        n An xDimX1 normal to the plane. If this is omitted or an empty
%          matrix is passed, then n is a vector that is all zero except for
%          the last element, which is 1. In 3D that is equivalent to using
%          the x-y plane.
%
%OUTPUTS: angVal The angle between the line and the plane given in radians.
%
%This just uses the identity that
%sin(angVal)=abs(dot(a,n)/(norm(a)*norm(n)))
%
%EXAMPLE:
%An easy test for the function is to use the x-y plane and just have a line
%in the x-z plane defined at a specified angle to the x axis. In that
%instance, the result of this function can be quickly compared ot the
%truth.
% n=[0;0;1];
% theta=pi/3;
% a=[cos(theta);0;sin(theta)];
% AbsErr=abs(angleBetweenLineAndPlane(a,n)-theta)
%In this instance the absolute error is 0.
%
%July 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

angVal=asin(abs(dot(a,n)/(norm(a)*norm(n))));

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [o, p, k] = R2opk(R)

% % phi --------------------------------------------------------------------------
%
% try
%
%     p = asing( R(1,3));
%     p = ang0400(p);
%
% catch
%
%     p = NaN;
%
% end
%
% % omega ------------------------------------------------------------------------
%
% try
%
%     o = asing(-R(2,3) / cosg(p));
%
% catch
%
%     o = NaN;
%
% end
%
% % kappa ------------------------------------------------------------------------
%
% try
%
%     k1_1 = asing(-R(1,2) / cosg(p)); % pos or neg
%
%     if k1_1 >= 0, k1_2 = ang0400(200 - k1_1);      end
%     if k1_1 <  0, k1_2 = ang0400(200 + abs(k1_1)); end
%
%     k1_1 = ang0400(k1_1); % always pos
%
%     k2_1 = acosg(R(1,1) / cosg(p)); % always pos
%     k2_2 = ang0400(-k2_1);
%
%     if abs(k1_1 - k2_1) < 1, k = mean([k1_1 k2_1]); end
%     if abs(k1_1 - k2_2) < 1, k = mean([k1_1 k2_2]); end
%     if abs(k1_2 - k2_1) < 1, k = mean([k1_2 k2_1]); end
%     if abs(k1_2 - k2_2) < 1, k = mean([k1_2 k2_2]); end
%
%     if ~exist('k'), k = NaN; end
%
% catch
%
%     k = NaN;
%
% end

% From omphika.m
o = atan2(-R(2,3), R(3,3)) * 200/pi;
p = asin(R(1,3))           * 200/pi;
k = atan2(-R(1,2), R(1,1)) * 200/pi;

end

author: pglira


text:
function y = atan(x)
%ATAN         Implements  atan(x)  for intervals
%
%   y = atan(x)
%
%interval standard function implementation
%

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  complex allowed, sparse input,
%                                  major revision, improved accuracy
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/06/07     S.M. Rump  approximate std fcts removed
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  if issparse(x)
    [ix,jx,sx] = find(x);
    [m,n] = size(x);
    y = sparse(ix,jx,atan(full(sx)),m,n);
    return
  end

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

  if x.complex
    y = log( 2./(1-j*x) - 1) / (2*j);
    if rndold
      setround(rndold)
    end
    return
  end

  y = x;

  xinf = x.inf(:);
  xsup = x.sup(:);

  IndexInfPos = ( xinf>=0 );
  len1 = sum(IndexInfPos);
  IndexSupNeg = ( xsup<=0 );

  Y = atan_pos( [ xinf(IndexInfPos) ; -xsup(IndexSupNeg) ] , -1 );
  y.inf(IndexInfPos) = Y(1:len1);
  y.sup(IndexSupNeg) = -Y( len1+1 : end );

  IndexInfNeg = ( xinf<0 );
  len1 = sum(IndexInfNeg);
  IndexSupPos = ( xsup>0 );

  Y = atan_pos( [ -xinf(IndexInfNeg) ; xsup(IndexSupPos) ] , 1 );
  y.inf(IndexInfNeg) = -Y(1:len1);
  y.sup(IndexSupPos) = Y( len1+1 : end );

  setround(rndold)

repo: OpenMAS


text:
function loc_3d_torus(NumOfSamples)

% Create Localized 3-d data from uniform distribution.
% Data are bount to the neighborhood near a torus surface.
% Torus surface equation:
%
%         (c - sqrt(x^2+y^2))^2 + z^2 = a^2
%
% This torus is created by revoloution around z axis of a circle centered
% at (c,0) with radious a.

clc;
input_dims = 3;

a = 1;
a1 = .9;
c = 5;
% Initialize Data
Data = [];

while size(Data,2)<NumOfSamples

In1 = 20*(rand(input_dims,1)-.5);

% Apply Mask to input Data
if abs(In1(3,1))<= sqrt(a^2 - (c-sqrt(In1(1,1)^2 + In1(2,1)^2))^2) && ...
                                                                       abs(In1(3,1))>=sqrt(a1^2 - (c-sqrt(In1(1,1)^2 + In1(2,1)^2))^2)
    Data = [Data In1];
end
end

a = randperm(NumOfSamples);
Data = Data(:,a);

save local_torus Data;
plot3(Data(1,:),Data(2,:),Data(3,:),'.','MarkerSize',1);
grid on;
clear all;

author: Sable


text:
function obj = Calc_Mannings_Landcover(obj,data,type,varargin)
% obj = Calc_Mannings_Landcover(obj,data,type,varargin)
% Input a msh class object and interpolates a land-cover database file
% (netcdf format) onto the msh while converting to mannings values through
% a look-up table
%
% type:
%  - 'nlcd':
%    NLCD 2006 Land Cover (2011 Edition) (1.1Gb) from
%    https://www.mrlc.gov/nlcd06_data.php
%  - 'ccap':
%    CCAP NOAA land cover:
%    https://coast.noaa.gov/data/digitalcoast/pdf/ccap-class-scheme-regional.pdf
%
% varargin: Accepts the same options as for msh.interp to control how
%           data is interpolated; see 'help msh.interp'
%
%  Author:            KR & WP July 19, 2018
%  Update for CCAP:   WP May 5, 2021
%  Update for storing incremental interpolation results: KJ & OK July 9, 2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 3 || isempty(type)
   type = 'nlcd';
end

if strcmp(type,'nlcd')
    disp('Info: Using NLCD table')
    % NLCD table
    load nlcd
    nlcd_mannings(end+1)=nlcd_mannings(11); % <--make NaN = default value
    varargin{end+1}='lut';
    varargin{end+1}=nlcd_mannings;
elseif strcmp(type,'ccap')
    disp('Info: Using CCAP table')
    % CCAP table
    load ccap
    ccap_mannings(end+1)=ccap_mannings(21); % <-- make NaN = default value
    varargin{end+1}='lut';
    varargin{end+1}=ccap_mannings;
else
    error('Land-cover database not supported')
end
% The mannings name and default value
attrname = 'mannings_n_at_sea_floor';
default_val = varargin{end}(end); %end of the lut is default
dmy = msh();  dmy.p = obj.p; dmy.t = obj.t;
% Convert to Mannings and interpolate how the user wants
dmy = interp(dmy,data,varargin{:});
Man = dmy.b';
Man( isnan(Man) | Man==0 ) = default_val; %points outside of lut can still be NaN
%% Make into f13 struct
if isempty(obj.f13)
    % Add add mannings as first entry in f13 struct
    obj.f13.AGRID = obj.title;
    obj.f13.NumOfNodes = length(obj.p);
    obj.f13.nAttr = 1;
    NA = 1;
else
    broken = 0;
    for NA = 1:obj.f13.nAttr
        if strcmp(attrname,obj.f13.defval.Atr(NA).AttrName)
            broken = 1;
            % overwrite existing tau0
            break
        end
    end
    if ~broken
        % add mannings to list
        obj.f13.nAttr = obj.f13.nAttr + 1;
        NA = obj.f13.nAttr;
    end
end

% Default Values
obj.f13.defval.Atr(NA).AttrName = attrname;
% We can just put in the options here
obj.f13.defval.Atr(NA).Unit = 'unitless';
valpernode = 1;
obj.f13.defval.Atr(NA).ValuesPerNode = valpernode ;
obj.f13.defval.Atr(NA).Val = default_val ;

% User Values
obj.f13.userval.Atr(NA).AttrName = attrname ;
numnodes = length(find(Man ~= default_val));
obj.f13.userval.Atr(NA).usernumnodes = numnodes ;
% Print out list of nodes for each
K = find(Man ~= default_val);
if isfield(obj.f13.userval.Atr(NA),'Val')
    obj.f13.userval.Atr(NA).Val = [obj.f13.userval.Atr(NA).Val'; K' , Man(K)']';
else
    obj.f13.userval.Atr(NA).Val = [K ; Man(K)];
end
%EOF
end

author: CHLNDDEV


text:
% Equalizer design example
% "Filter design" lecture notes (EE364) by S. Boyd
% (figures are generated)
%
% Designs a frequency-domain and time-domain FIR equalizer for
% a single-input single-output (SISO) channel.
%
% Frequency-domain equalization uses a Chebychev criteria and
% is specified in terms of frequency response functions.
% It is a convex problem (which can be formulated as an SOCP):
%
%   minimize   max |G(w)H(w) - G_des(w)|     for w in [0,pi]
%
% where H is the frequency response function and our variable
% is the filter impulse response h. Function G is the unequalized
% frequency response and G_des is the desired freq response.
%
% Time-domain equalization immediately designs the impulse
% response function by specifying the problem in time (it's an LP):
%
%   minimize   max_{t neq D} |g_tilde(t)|
%       s.t.   g_tilde(D) = 1
%
% where g_tilde is the impulse response of equalized system,
% and D is the delay of the system.
%
% Written for CVX by Almir Mutapcic 02/02/06

%********************************************************************
% problem specs
%********************************************************************
% sample channel with impulse response g
g =.5*[ 0.6526;  0.2157; -0.2639;  1.8024; -0.6430; ...
        0.1096; -0.7190;  0.4206; -0.0193;  0.6603;];

% problem parameters
n  = 30;              % filter order
D  = 10;              % overall delay

%********************************************************************
% frequency domain equalization
%********************************************************************
% number of freq samples (rule-of-thumb)
m  = 15*(length(g) + n);

w = linspace(0,pi,m)';
G = exp( -j*kron(w,[0:length(g)-1]) )*g;
A = exp( -j*kron(w,[0:n-1]) );

% desired frequency response is a pure delay (equalized channel)
Gdes = exp(-j*D*w);

% formulate and solve the Chebyshev design problem
cvx_begin
  variable hf(n,1)
  minimize( max( abs( G.*(A*hf) - Gdes ) ) )
cvx_end

% check if problem was successfully solved
disp(['Frequency equalization problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
end

%********************************************************************
% time-domain equalization
%********************************************************************
% define the convolution matrix
Tconv = toeplitz([g; zeros(n-1,1)],[g(1) zeros(1,n-1)]);

% create array of all times without t=D
times_not_D = [1:D D+2:size(Tconv,1)];

% formulate and solve the time equalization problem
cvx_begin
  variable t
  variable ht(n,1)

  minimize( max( abs( Tconv(times_not_D,:)*ht ) ) )
  subject to
    Tconv(D+1,:)*ht == 1;
cvx_end

% check if problem was successfully solved
if ~strfind(cvx_status,'Solved')
  disp(['Frequency equalization problem is ' cvx_status])
  return
end

%********************************************************************
% equalizer plots
%********************************************************************
% plot g
figure(1)
plot([0:length(g)-1],g,'o',[0:length(g)-1],g,'b:')
xlabel('t')
ylabel('g(t)')

figure(2)
H = exp(-j*kron(w,[0:length(g)-1]))*g;
% magnitude
subplot(2,1,1);
plot(w,20*log10(abs(H)))
axis([0,pi,-20,20])
xlabel('w')
ylabel('mag G(w) in dB')
% phase
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w')
ylabel('phase G(w)')

% freq equalizer
figure(3)
plot([0:n-1],hf,'o',[0:n-1],hf,'b:')
xlabel('t')
ylabel('h(t)')

% plot g_tilde
figure(4)
gt=conv(g,hf);
plot([1:length(gt)]-1,gt,'o',[1:length(gt)]-1,gt,'b:')
xlabel('t')
ylabel('g tilde(t)')
axis([0,length(gt)-1,-.2 1.2])

figure(5)
H = exp(-j*kron(w,[0:length(gt)-1]))*gt;
% amplitude
subplot(2,1,1)
plot(w,20*log10(abs(H)))
axis([0,pi,-20,20])
xlabel('w')
ylabel('mag G tilde(w) in dB')
% phase
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w')
ylabel('phase G tilde(w)')

% time equalizer
figure(6)
plot([0:n-1],ht,'o',[0:n-1],ht,'b:')
xlabel('t')
ylabel('h(t)')

% plot g_tilde
figure(7)
gt=conv(g,ht);
plot([1:length(gt)]-1,gt,'o',[1:length(gt)]-1,gt,'b:')
xlabel('t')
ylabel('g tilde(t)')

figure(8)
H = exp(-j*kron(w,[0:length(gt)-1]))*gt;
% magnitude
subplot(2,1,1)
plot(w,20*log10(abs(H)))
axis([0,pi,-20,20])
xlabel('w')
ylabel('mag G tilde(w) in dB')
% phase
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel('w')
ylabel('phase G tilde(w)')

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function test_bug2909

% WALLTIME 00:10:00
% MEM 3gb
% DEPENDENCY ft_read_cifti

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2909'));

%%
% this is the data mentioned by JM on http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2865
filename1 = 'RSN-networks.32k_fs_LR.dlabel.nii';
filename2 = 'RSN-networks.32k_fs_LR.dlabel.nii';
surfL = 'Q1-Q6_R440.L.midthickness.32k_fs_LR.surf.gii';
surfR = 'Q1-Q6_R440.R.midthickness.32k_fs_LR.surf.gii';

cii1 = ft_read_cifti(filename1, 'cortexleft', surfL, 'cortexright', surfR, 'mapname', 'field');
cii2 = ft_read_cifti(filename2, 'cortexleft', surfL, 'cortexright', surfR, 'mapname', 'array');

% this is the file I got from Mike Harms over email
filename3 = 'rfMRI_REST1_LR_Atlas_stats_ALLCOL_MEAN.dscalar.nii';
filename4 = 'rfMRI_REST1_LR_Atlas_stats_ALLCOL_MEAN.dscalar.nii';

cii3 = ft_read_cifti(filename3, 'debug', true, 'maplabel', 'field');
cii4 = ft_read_cifti(filename4, 'debug', true, 'maplabel', 'array');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% try writing them and subsequently reading them again

filename1b = tempname;
ft_write_cifti(filename1b, cii1, 'parameter', 'rsn_consensus_communities___pcn11__power_neuron11');
d = dir([filename1b '*.nii']);
filename1b = fullfile(fileparts(filename1b), d.name); % determine the name of the cii file
cii1b = ft_read_cifti(filename1b, 'cortexleft', surfL, 'cortexright', surfR);
assert(isequal(cii1.rsn_consensus_communities___pcn11__power_neuron11, cii1b.dscalar));

%% this one does not work yet
% filename2b = tempname;
% ft_write_cifti(filename2b, cii2, 'parameter', 'x32k_fs_lr');
% d = dir([filename2b '*.nii']);
% filename1b = fullfile(fileparts(filename2b), d.name); % determine the name of the cii file
% cii2b = ft_read_cifti(filename1b, 'cortexleft', surfL, 'cortexright', surfR);

author: fieldtrip

repo: fieldtrip