text:
%==========================================================================
%
%  By Gennady Zilberman, Ben-Gurion University of the Negev
%  genazilberman@yahoo.com
%
%   This software performs Forward Channel simulation of IS-95 standrd
%   system.
%
%   The simulation environment: Static AWGN channel.
%   Link Rate = 9600 KBps
%
%==========================================================================

global Zi Zq Zs show R Gi Gq

clear j;
show = 0;
SD = 0;         % ---- Soft/Hard Decision Receiver

%==========================================================================
% ---------------------------------------- MAIN SIMULATION DEFINITIONS --------------------------------
%==========================================================================
BitRate = 9600;
ChipRate = 1228800;
N = 184;  		         % 9.6 KBps rate -> 184 netto data bits in each 20 msec packet
MFType = 1;		     % Matched filter type - Raised Cosine
R = 5;			            % Analog Signal Simulation rate

% ------------------------ Viterbi Polynom -------------------
G_Vit = [1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1];
K = size(G_Vit, 2); 		% number of cheap per data bit
L = size(G_Vit, 1); 		% number of cheap per data bit

% ------------------------ Walsh Matrix -------------------
WLen = 64;
Walsh = reshape([1;0]*ones(1, WLen/2), WLen , 1);
%Walsh = zeros(WLen ,1);

%------- Spreading PN polynomials ---------
%Gi = [ 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1]';
%Gq = [ 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0 1]';

Gi_ind = [15, 13, 9, 8, 7, 5, 0]';
Gq_ind = [15, 12, 11, 10, 6, 5, 4, 3, 0]';

Gi = zeros(16, 1);
Gi(16-Gi_ind) = ones(size(Gi_ind));
Zi = [zeros(length(Gi)-1, 1); 1];     % Initial State for I-Channel PN generator

Gq = zeros(16, 1);
Gq(16-Gq_ind) = ones(size(Gq_ind));
Zq = [zeros(length(Gq)-1, 1); 1]; 	  % Initial State for Q-Channel PN generator

%------- Scrambler Generation Polynom ---------
Gs_ind = [42, 35, 33, 31, 27, 26, 25, 22, 21, 19, 18, 17, 16, 10, 7, 6, 5, 3, 2, 1, 0]';
Gs = zeros(43, 1);
Gs(43-Gs_ind) = ones(size(Gs_ind));
Zs = [zeros(length(Gs)-1, 1); 1];  		% Initial State for Long Sequence Generator

%------- AWGN definitions -------------
EbEc = 10*log10(ChipRate/BitRate); 		% Bit/Chip rate loss
EbEcVit = 10*log10(L);
EbNo = [-2 : 0.5 : 6.5]; 		% measured EbNo range (dB)

%==========================================================================
% ----------------------------------------------- MAIN SIMULATION LOOP ---------------------------------
%==========================================================================
ErrorsB = []; ErrorsC = []; NN = [];
if (SD == 1)
   fprintf('\n SOFT Decision Viterbi Decoder\n\n');
else
   fprintf('\n HARD Decision Viterbi Decoder\n\n');
end

for i=1:length(EbNo)
   fprintf('\nProcessing %1.1f (dB)', EbNo(i));
   iter = 0;	ErrB = 0; ErrC = 0;
   while (ErrB <300) & (iter <150)
      drawnow;

      %----------------------------------- Data Transmit ---------------------------------
      TxData = (randn(N, 1)>0);          %- Gamble the Tx data

      [TxChips, Scrambler] = PacketBuilder(TxData, G_Vit, Gs);		% 19.2 Kcps rate
      [x PN MF] = Modulator(TxChips, MFType, Walsh);  % 1.2288 Mcps rate

      %-------------------------------- AWGN Channel ------------------------------------
      noise = 1/sqrt(2)*sqrt(R/2)*( randn(size(x)) + j*randn(size(x)))*10^(-(EbNo(i) - EbEc)/20);
      r = x+noise;

      %------------------------------------ Receiver ---------------------------------
      RxSD = Demodulator(r, PN, MF, Walsh); 		% Soft Decision   19.2 Kcps
      RxHD = (RxSD>0); 		% Define Hard Decisions received chips
      if (SD)
        % RxSD = sign(RxSD); 		% SOFT Decision = HARD Decision
        % RxSD = sign(RxSD).*(log2(abs(RxSD)+1)); 		% SOFT decision = 3 bits representation
        [RxData Metric]= ReceiverSD(RxSD, G_Vit, Scrambler); 		% Get Received Data 9.6 KBps (Soft Decision)
      else
        [RxData Metric]= ReceiverHD(RxHD, G_Vit, Scrambler); 		% Get Received Data 9.6 KBps (Hard Decision)
      end

      if(show)
         subplot(311); plot(RxSD, '-o'); title('Soft Decisions');
         subplot(312); plot(xor(TxChips, RxHD), '-o'); title('Chip Errors');
         subplot(313); plot(xor(TxData, RxData), '-o');
         title(['Data Bit Errors. Metric = ', num2str(Metric)]);
         pause;
      end

      if(mod(iter, 50)==0)
         fprintf('.');
         save TempResults ErrB ErrC N iter
      end

      ErrB = ErrB + sum(xor(RxData, TxData));  % Data bits Error
      ErrC = ErrC + sum(xor(RxHD, TxChips));  	  % Chip Error (before Viterbi)
      iter = iter+ 1;
   end
   %---------------------------------------- Save the data of current  iteration ---------------------
   ErrorsB = [ErrorsB; ErrB];
   ErrorsC = [ErrorsC; ErrC];
   NN = [NN; N*iter];
   save SimData *
end

%---------------------------------------- Calculate Error's probability ----------------------------
PerrB = ErrorsB./NN;
%PerrB1 = ErrorsB1./NN1;
PerrC = ErrorsC./NN;
Pbpsk= 1/2*erfc(sqrt(10.^(EbNo/10)));
PcVit= 1/2*erfc(sqrt(10.^((EbNo-EbEcVit)/10)));
Pc =   1/2*erfc(sqrt(10.^((EbNo-EbEc)/10)));

%---------------------------------------- Show simulation Results -------- ---------------------
figure;
semilogy(EbNo(1:length(PerrB)), PerrB, 'r-o'); hold on;
%semilogy(EbNo(1:length(PerrB1)), PerrB1, 'k-o'); hold on;
semilogy(EbNo(1:length(PerrC)), PerrC, 'm-o'); grid on;
semilogy(EbNo, Pbpsk, 'b-.x'); xlabel('EbNo (dB)');
%semilogy(EbNo, PcVit, 'k-.x'); ylabel('BER');
semilogy(EbNo, Pc, 'g-.x');

title('Error probability as function of EbNo');

legend('Pb of System (HD)', 'Pb of System (SD)', 'Pc before Viterbi of System', ...
   'Pb of BPSK with no Viterbi (theory)', 'Pc on Receiver (theory)');

legend('Pb of System', 'Pc before Viterbi of System', ...
   'Pb of BPSK with no Viterbi (theory)', 'Pc before Viterbi (theory)', 'Pc on Receiver (theory)');

repo: mcbench-benchmarks


text:
function gx = g_ip(x,P,u,in)

gx = P(1) + real((u+P(3)).^(1-P(2)));

if isnan(gx)
    gx
end

author: MBB-team

repo: VBA-toolbox


text:
%KWAVEGRID   Class definition for k-Wave grid structure.
%
% DESCRIPTION:
%       See makeGrid for function arguments.
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 22nd July 2010
%       last update - 29th August 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% 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/>.

classdef kWaveGrid

    % define the private properties (these parameters are stored and cannot
    % be altered by the user after the object is created, but they can be
    % accessed)
    properties (GetAccess = public, SetAccess = private)

        % grid dimensions in grid points
        Nx = 0;
        Ny = 0;
        Nz = 0;

        % grid point spacing in metres
      	dx = 0;
        dy = 0;
        dz = 0;

        % wavenumber grids
        kx_vec = 0;
        ky_vec = 0;
        kz_vec = 0;
        k = 0;

        % maximum supported spatial frequency
        kx_max = 0;
        ky_max = 0;
        kz_max = 0;
        k_max = 0;

        % number of dimensions
        dim = 0;

        % non-uniform grid
        nonuniform = false;

    end

    % define the public properties (these parameters are stored and can be
    % both accessed and altered by the user after the object is created)
    properties (GetAccess = public, SetAccess = public)

        % time array
        t_array = 'auto';

    end

    % define the dependent properties (these parameters are not stored but
    % re-computed each time they are needed)
    properties(Dependent = true, GetAccess = public, SetAccess = private)

        % 3D plaid Cartesian spatial grids
        x;
        y;
        z;

        % 3D plaid wavenumber grids
        kx;
        ky;
        kz;

        % 1D spatial grids
        x_vec;
        y_vec;
        z_vec;

        % Cartesian grid size
        x_size;
        z_size;
        y_size;

        % time step and number of time points
        dt;
        Nt;

        % total number of grid points
        total_grid_points;

    end

    % define the staggered grid properties (these can be accessed but not
    % altered by the user after the object is created, the parameters are
    % also hidden from view)
    properties (Hidden = true, GetAccess = public, SetAccess = private)

        % position vectors for the staggered grid points in [0, 1]
        xn_vec = 0;
        yn_vec = 0;
        zn_vec = 0;
        xn_vec_sgx = 0;
        yn_vec_sgy = 0;
        zn_vec_sgz = 0;

        % transformation gradients between uniform and staggered grids
        dxudxn = 0;
        dyudyn = 0;
        dzudzn = 0;
        dxudxn_sgx = 0;
        dyudyn_sgy = 0;
        dzudzn_sgz = 0;

    end

    % define the dependent staggered grid properties (these parameters are
    % not stored but re-computed each time they are needed, the parameters
    % are also hidden from view)
    properties(Dependent = true, Hidden = true, GetAccess = public, SetAccess = private)

        % 3D plaid non-uniform spatial grids
        xn;
        yn;
        zn;

    end

    % constructor function
    methods
        function kgrid = kWaveGrid(varargin)

            % assign the input values to the grid object
            switch nargin

                case 6

                    % 3D uniform grid
                    kgrid.Nx = varargin{1};
                    kgrid.dx = varargin{2};
                    kgrid.Ny = varargin{3};
                    kgrid.dy = varargin{4};
                    kgrid.Nz = varargin{5};
                    kgrid.dz = varargin{6};

                    % set the number of dimensions
                    kgrid.dim = 3;

                case 4

                    % 2D uniform grid
                    kgrid.Nx = varargin{1};
                    kgrid.dx = varargin{2};
                    kgrid.Ny = varargin{3};
                    kgrid.dy = varargin{4};

                    % set the number of dimensions
                    kgrid.dim = 2;

                case 2

                    % 1D uniform grid
                    kgrid.Nx = varargin{1};
                    kgrid.dx = varargin{2};

                    % set the number of dimensions
                    kgrid.dim = 1;

                otherwise
                    error('Incorrect number of input arguments');
            end

            switch kgrid.dim
                case 1

                    % assign the grid parameters for the x spatial direction
                    kgrid.kx_vec = kgrid.makeDim(kgrid.Nx, kgrid.dx);

                    % define the scalar wavenumber based on the wavenumber components
                    kgrid.k = abs(kgrid.kx_vec);

                    % define maximum supported frequency
                    kgrid.kx_max = max(abs(kgrid.kx_vec(:)));
                    kgrid.k_max = kgrid.kx_max;

                case 2

                    % assign the grid parameters for the x and z spatial directions
                    kgrid.kx_vec = kgrid.makeDim(kgrid.Nx, kgrid.dx);
                    kgrid.ky_vec = kgrid.makeDim(kgrid.Ny, kgrid.dy);

                    % define the wavenumber based on the wavenumber components
                    % (using bsxfun saves memory by avoiding ndgrid)
                    kgrid.k = zeros(kgrid.Nx, kgrid.Ny);
                    kgrid.k = bsxfun(@plus, (reshape(kgrid.kx_vec, [], 1, 1).^2), kgrid.k);
                    kgrid.k = bsxfun(@plus, (reshape(kgrid.ky_vec, 1, [], 1).^2), kgrid.k);
                    kgrid.k = sqrt(kgrid.k);

                    % define maximum supported frequency
                    kgrid.kx_max = max(abs(kgrid.kx_vec(:)));
                    kgrid.ky_max = max(abs(kgrid.ky_vec(:)));
                    kgrid.k_max = min([kgrid.kx_max, kgrid.ky_max]);

                case 3

                    % assign the grid parameters for the x ,y and z spatial directions
                    kgrid.kx_vec = kgrid.makeDim(kgrid.Nx, kgrid.dx);
                    kgrid.ky_vec = kgrid.makeDim(kgrid.Ny, kgrid.dy);
                    kgrid.kz_vec = kgrid.makeDim(kgrid.Nz, kgrid.dz);

                    % define the wavenumber based on the wavenumber components
                    % (using bsxfun saves memory by avoiding ndgrid)
                    kgrid.k = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz);
                    kgrid.k = bsxfun(@plus, (reshape(kgrid.kx_vec, [], 1, 1).^2), kgrid.k);
                    kgrid.k = bsxfun(@plus, (reshape(kgrid.ky_vec, 1, [], 1).^2), kgrid.k);
                    kgrid.k = bsxfun(@plus, (reshape(kgrid.kz_vec, 1, 1, []).^2), kgrid.k);
                    kgrid.k = sqrt(kgrid.k);

                    % define maximum supported frequency
                    kgrid.kx_max = max(abs(kgrid.kx_vec(:)));
                    kgrid.ky_max = max(abs(kgrid.ky_vec(:)));
                    kgrid.kz_max = max(abs(kgrid.kz_vec(:)));
                    kgrid.k_max = min([kgrid.kx_max, kgrid.ky_max, kgrid.kz_max]);

            end
        end
    end

    % functions for dependent variables that only run when queried
    methods

        function x = get.x(obj)
            % calculate x based on kx
            x = obj.x_size*obj.kx*obj.dx/(2*pi);
        end

        function y = get.y(obj)
            % calculate y based on ky
            y = obj.y_size*obj.ky*obj.dy/(2*pi);
        end

        function z = get.z(obj)
            % calculate z based on kz
            z = obj.z_size*obj.kz*obj.dz/(2*pi);
        end

        function kx = get.kx(obj)
            % duplicate kx vector
            switch obj.dim
                case 1
                    kx = obj.kx_vec;
                case 2
                    kx = repmat(obj.kx_vec, [1, obj.Ny]);
                case 3
                    kx = repmat(obj.kx_vec, [1, obj.Ny, obj.Nz]);
            end
        end

        function ky = get.ky(obj)
            % rotate ky vector to y direction then duplicate
            switch obj.dim
                case 1
                    ky = 0;
                case 2
                    ky = repmat(obj.ky_vec.', [obj.Nx, 1]);
                case 3
                    ky = repmat(obj.ky_vec.', [obj.Nx, 1, obj.Nz]);
            end
        end

        function kz = get.kz(obj)
            % permute kz vector to z direction then duplicate
            switch obj.dim
                case 1
                    kz = 0;
                case 2
                    kz = 0;
                case 3
                    kz = repmat(permute(obj.kz_vec, [2 3 1]), [obj.Nx, obj.Ny, 1]);
            end
        end

        function x_vec = get.x_vec(obj)
            % calculate x_vec based on kx_vec
            x_vec = obj.x_size*obj.kx_vec*obj.dx/(2*pi);
        end

        function y_vec = get.y_vec(obj)
            % calculate y_vec based on ky_vec
            y_vec = obj.y_size*obj.ky_vec*obj.dy/(2*pi);
        end

        function z_vec = get.z_vec(obj)
            % calculate z_vec based on kz_vec
            z_vec = obj.z_size*obj.kz_vec*obj.dz/(2*pi);
        end

        function x_size = get.x_size(obj)
            % calculate x_size based on Nx and dx
            x_size = obj.Nx*obj.dx;
        end

        function y_size = get.y_size(obj)
            % calculate y_size based on Ny and dy
            y_size = obj.Ny*obj.dy;
        end

        function z_size = get.z_size(obj)
            % calculate z_size based on Nz and dz
            z_size = obj.Nz*obj.dz;
        end

        function dt = get.dt(obj)
            % extract dt from the time array
            if strcmp(obj.t_array, 'auto')
                dt = 'auto';
            else
                dt = obj.t_array(2) - obj.t_array(1);
            end
        end

        function Nt = get.Nt(obj)
            % extract dt from the time array
            if strcmp(obj.t_array, 'auto')
                Nt = 'auto';
            else
                Nt = length(obj.t_array);
            end
        end

        function N = get.total_grid_points(obj)
            % calculate the total number of points in the grid
            switch obj.dim
                case 1
                    N = obj.Nx;
                case 2
                    N = obj.Nx*obj.Ny;
                case 3
                    N = obj.Nx*obj.Ny*obj.Nz;
            end
        end
    end

    % functions that can only be accessed by class members
    methods (Access = 'protected', Static = true)

        % subfunction to create the grid parameters for a single spatial direction
        function kx_vec = makeDim(Nx, dx)

            % define the discretisation of the spatial dimension such that there is
            % always a DC component
            if rem(Nx, 2) == 0
                % grid dimension has an even number of points
                nx = ((-Nx/2:Nx/2-1)/Nx).';
            else
                % grid dimension has an odd number of points
                nx = ((-(Nx-1)/2:(Nx-1)/2)/Nx).';
            end

            % force middle value to be zero in case 1/Nx is a recurring
            % number and the series doesn't give exactly zero
            nx(floor(Nx/2) + 1) = 0;

            % define the wavenumber vector components
            kx_vec = (2*pi/dx).*nx;
        end

    end

    % functions for non-uniform grids
    methods

        % subfunction to return plaid xn matrix
        function xn = get.xn(obj)
            if obj.nonuniform
                % duplicate xn vector
                switch obj.dim
                    case 1
                        xn = obj.xn_vec;
                    case 2
                        xn = repmat(obj.xn_vec, [1, obj.Ny]);
                    case 3
                        xn = repmat(obj.xn_vec, [1, obj.Ny, obj.Nz]);
                end
            else
                xn = 0;
            end
        end

        % subfunction to return plaid yn matrix
        function yn = get.yn(obj)
            if obj.nonuniform
                % rotate yn vector to y direction then duplicate
                switch obj.dim
                    case 1
                        yn = 0;
                    case 2
                        yn = repmat(obj.yn_vec.', [obj.Nx, 1]);
                    case 3
                        yn = repmat(obj.yn_vec.', [obj.Nx, 1, obj.Nz]);
                end
            else
                yn = 0;
            end
        end

        % subfunction to return plaid zn matrix
        function zn = get.zn(obj)
            if obj.nonuniform
                % permute kz vector to z direction then duplicate
                switch obj.dim
                    case 1
                        zn = 0;
                    case 2
                        zn = 0;
                    case 3
                        zn = repmat(permute(obj.zn_vec, [2 3 1]), [obj.Nx, obj.Ny, 1]);
                end
            else
                zn = 0;
            end
        end

        % subfunction to set non-uniform grid parameters in specified
        % dimension
        function obj = setNUGrid(obj, dim, xn_vec, dxudxn, xn_vec_sgx, dxudxn_sgx)

            % check the dimension to set the nonuniform grid is appropriate
            if dim > obj.dim
                error(['Cannot set nonuniform parameters for dimension ' num2str(dim) ' of ' num2str(obj.dim) '-dimensional grid']);
            end

            % force non-uniform grid spacing to be column vectors, and the
            % gradients to be in the correct direction for use with bsxfun
            switch dim
                case 1
                    obj.xn_vec = reshape(xn_vec, [], 1);
                    obj.dxudxn = reshape(dxudxn, [], 1);
                    obj.xn_vec_sgx = reshape(xn_vec_sgx, [], 1);
                    obj.dxudxn_sgx = reshape(dxudxn_sgx, [], 1);
                case 2
                    obj.yn_vec = reshape(xn_vec, [], 1);
                    obj.dyudyn = reshape(dxudxn, 1, []);
                    obj.yn_vec_sgy = reshape(xn_vec_sgx, [], 1);
                    obj.dyudyn_sgy = reshape(dxudxn_sgx, 1, []);
                case 3
                    obj.zn_vec = reshape(xn_vec, [], 1);
                    obj.dzudzn = reshape(dxudxn, 1, 1, []);
                    obj.zn_vec_sgz = reshape(xn_vec_sgx, [], 1);
                    obj.dzdzn_sgz = reshape(dxudxn_sgx, 1, 1, []);
            end

            % set non-uniform flag
            obj.nonuniform = true;

        end
    end
end

repo: Aero-matlab


text:
function [mvy2,mvx2,mvz2]=expand_motion_field(mvy,mvx,mvz,newdim,offsets)

if length(offsets) == 1
	offsets = [0 0 offsets];
end

mvy2=zeros(newdim,class(mvy));
mvx2 = mvy2;
mvz2 = mvy2;

dim = size(mvy);
ub = dim+offsets;

ys = 1:newdim(1); ys2 = ys-offsets(1); ys2(ys2<1) = 1; ys2(ys2>dim(1)) = dim(1);
xs = 1:newdim(2); xs2 = xs-offsets(2); xs2(xs2<1) = 1; xs2(xs2>dim(2)) = dim(2);
zs = 1:newdim(3); zs2 = zs-offsets(3); zs2(zs2<1) = 1; zs2(zs2>dim(3)) = dim(3);

mvy2(ys,xs,zs) = mvy(ys2,xs2,zs2);
mvx2(ys,xs,zs) = mvx(ys2,xs2,zs2);
mvz2(ys,xs,zs) = mvz(ys2,xs2,zs2);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function dcmobj = build_module_template(moduleName)
% function tags = build_module_template(moduleName)
%"build_module_template"
%   Create a Java DICOM object containing fields used to represent the
%   specified module.
%
%   Valid moduleNames are currently 'patient'
%
%JRA 06/06/06
%NAV 07/19/16 updated to dcm4che3
% AI 11/20/2020 Added module subsets used in dcmdir_add to group images.
%
%Usage:
%   dcmobj = general_equipment_module_template;
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Create a dcmobj to hold the structure and data.
%%dcmobj = org.dcm4che2.data.BasicDicomObject;
  %dcmobj = org.dcm4che3.data.Attributes;

%Done = 0;

switch lower(moduleName)
    %Get the top level tags used in this module.
    case 'patient'
        tagS = patient_module_tags;
    case 'patient_subset'
        tagS = patient_module_tags_subset;
    case 'general_study'
        tagS = general_study_module_tags;
    case 'general_study_subset'
        tagS = general_study_module_tags_subset;
    case 'general_series'
        tagS = general_series_module_tags;
    case 'general_series_subset'
        tagS = general_series_module_tags_subset;
    case 'frame_of_reference'
        tagS = frame_of_reference_module_tags;
    case 'general_equipment'
        tagS = general_equipment_module_tags;
    case 'content_identification'
        tagS = content_identification_module_tags;
    case 'general_image'
        tagS = general_image_module_tags;
    case 'image_plane'
        tagS = image_plane_module_tags;
    case 'image_pixel'
        tagS = image_pixel_module_tags;
    case 'ct_image'
        tagS = CT_image_module_tags;
    case 'pt_image'
        tagS = PT_image_module_tags;
    case 'mr_image'
        tagS = MR_image_module_tags;
    case 'mr_image_subset'
        tagS = MR_image_module_tags_subset;
    case 'sop_common'
        tagS = SOP_common_module_tags;
    case 'structure_set'
        tagS = structure_set_module_tags;
    case 'roi_contour'
        tagS = ROI_contour_module_tags;
    case 'rt_series'
        tagS = RT_series_module_tags;
    case 'rt_roi_observations'
        tagS = RT_ROI_observations_module_tags;
    case 'rt_dose'
        tagS = RT_dose_module_tags;
    case 'rt_dvh'
        tagS = RT_DVH_module_tags;
    case 'multi_frame'
        tagS = multi_frame_module_tags;
    case 'gsps'
        tagS = gsps_module_tags;
    otherwise
        error('Unrecognized or unsupported module template requested.')
end

%tags = hex2dec({tagS.tag});

%Create all top level tags in this object.
dcmobj = javaObject('org.dcm4che3.data.Attributes');
dcmobj = createEmptyFields(dcmobj, tagS);


text:
clear all;

[X,Y,Z] = peaks(100);
X = single(X);
Y = single(Y);
Z = single(Z);
Z_index = uint8((Z - min(Z(:))) * (255 / (max(Z(:)) - min(Z(:)))));
Z_color = uint8(ind2rgb(Z_index, hsv(256)*256));
Z_color_multiframe = reshape([Z_color(:)*0.2 Z_color(:)*0.6 Z_color(:)], 100, 100, 3, 3);
Z_color_noisy = uint8(single(Z_color) + rand(100, 100, 3).*50);

% 8-bit, grayscale image
saveastiff(uint8(Z_index), 'Z_uint8.tif');

% Lossless LZW compression
options.comp = 'lzw';
saveastiff(uint8(Z_index), 'Z_uint8_LZW.tif', options);
options.comp = 'no';

% Ask a question if the file is already exist
options.ask = true;
saveastiff(uint8(Z_index), 'Z_uint8_LZW.tif', options);
options.ask = false;

% Allow message printing.
options.message = true;
saveastiff(uint8(Z_index), 'Z_uint8_LZW.tif', options);
options.message = false;

% 16-bit, grayscale image
saveastiff(uint16(Z_index), 'Z_uint16.tif');

% 32-bit single, grayscale image
saveastiff(Z, 'Z_single.tif');

% RGB color image
options.color = true;
saveastiff(Z_color, 'Z_rgb.tif', options);
options.color = false;

% Save each R, G and B chanels of the color image, separately.
saveastiff(Z_color, 'Z_rgb_channel.tif');

% Save the multi-frame RGB color image
options.color = true;
saveastiff(Z_color_multiframe, 'Z_rgb_multiframe.tif', options);
options.color = false;

% Save the noise-added RGB color image
options.color = true;
saveastiff(Z_color_noisy, 'Z_rgb_noisy.tif', options);
options.color = false;

% 32-bit single, 50x50x50 volume data
saveastiff(single(rand(50, 50, 50)), 'volume_50x50x50.tif');

% Append option is ignored if path dose not exist.
options.append = true;
saveastiff(Z_index, 'Z_uint8_append.tif', options);
options.append = false;

% You can append any type of image to an existing tiff file.
options.append = true;
saveastiff(Z, 'complex_frame.tif');
saveastiff(single(rand(10, 10, 3)), 'complex_frame.tif', options);
options.color = true;
saveastiff(Z_color_multiframe, 'complex_frame.tif', options);
options.color = false;
options.append = false;

% Load multiframe tiff
multiframe = loadtiff('volume_50x50x50.tif');
complexframe = loadtiff('complex_frame.tif');

author: Sable


text:
function output = F1measure(P,R)
% Calculates the F1 measure with alpah = .5
%
% © Siamak Faridani, UC Berkeley, 2/2/2012
%
% v 1.0
%
% F = 1/(alpha(1/P)+(1-alpha)(1/R))
% F1 is when alpha is 0.5
%
% (see page 268 of Manning and Schutze)
%
% Inputs
%    P: precision
%    R: Recall
% Outputs
%    F1-measure

output = 2*(P*R)/(P+R);

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [PSCI,xSCI]=covarianceIntersectSCI(CovHyp,covType,xHyp,u,numSamp)
%%COVARIANCEINTERSECTSCI Perform sampling covariance intersection. This is
%                   a method of fusing the first two moments of estimates
%                   when the correlation between the estimates is unknown.
%                   Given the z(inverse) covariance matrices of the
%                   estimates, this function returns a covariance matrix.
%                   Alternatively, if the estimates themselves are given,
%                   this function can also return the merged estimate.
%                   Unlike covarianceIntersect, this function has a random
%                   component and does not specifically optimize any
%                   particular function. Compare this function to
%                   ellipsoidIntersect.
%
%INPUTS: CovHyp The xDimXxDimXN covariance matrices or inverse covariance
%           matrices of the values to be merged.
%   CovType An optional input specifying whether covariance matrices or
%           inverse covariance matrices are in CovHyp. Possible values are:
%           0 (The default if omitted or an empty matrix is passed) CovHyp
%             contains covariance matrices.
%           1 CovHyp contains inverse covariance matrices.
%      xHyp The optional xDimXN set of vectors to merge. These are only
%           needed if xSCI is requested on the output.
%         u A parameter between 0 and 1 that affects the performance of the
%           algorithm. The default if omitted or an empty matrix is passed
%           is 0.5.
%   numSamp The algorithm is stochastic. This optional input is the number
%           of samples to use. The default if this parameter is omitted or
%           an empty matrix is passed is max(100*(N-1),1);
%
%OUTPUTS: PSCI The xDimXxDim fused covariance matrix.
%         xSCI The merged estimate. This requires xHyp to be given on
%              the input. The covariance matrix associated with the
%              merged estimate is PSCI.
%
%This function implements the sampling covariance intersection algorithm of
%[1]. If only a single estimate is passed, then it is just returned. This
%algorithm was developed to fuse tracks from multiple sensors in a
%networked tracking environment.
%
%EXAMPLE:
% This is similar to the example given in the paper for when standard
% covariance intersection might be bad. , except the means are
% not identical.
% x=zeros(2,2);
% P=zeros(2,2,2);
% x(:,1)=[1;-2];
% P(:,:,1)=[1,0;0,100];
% x(:,2)=[-2;-1];
% P(:,:,2)=[100,0;0,1];
%
% [PM,xM]=covarianceIntersectSCI(P,[],x);
% figure()
% clf
% hold on
% drawEllipse(x(:,1),inv(P(:,:,1)),[],'--r')
% drawEllipse(x(:,2),inv(P(:,:,2)),[],'--g')
% drawEllipse(xM,inv(PM),[],'-b')
%
%REFERENCES:
%[1] X. Tian, Y. Bar-Shalom, and G. Chen, "A no-loss covariance
%    intersection algorithm for track-to-track fusion," in Proceedings of
%    SPIE: Signal and Data Processing of Small targets, vol. 7698,
%    Orlando, FL, 5 Apr. 2010.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(CovHyp,1);
numHyp=size(CovHyp,3);

if(nargin<5||isempty(numSamp))
    numSamp=max(100*(numHyp-1),1);
end

if(nargin<4||isempty(u))
    u=0.5;
end

if(nargin<2||isempty(covType))
   covType=0;
end

if(covType==0)%If covariance matrices are passed.
    PInvHyp=zeros(numDim,numDim,numHyp);

    for curHyp=1:numHyp
        PInvHyp(:,:,curHyp)=inv(CovHyp(:,:,curHyp));
    end
else%If inverse covariance matrices are passed.
    PInvHyp=CovHyp;
end

if(numHyp>0)
    P0=zeros(numDim,numDim);
    for curHyp=1:numHyp
        P0=P0+PInvHyp(:,:,curHyp);
    end
    P0Inv=P0;
    P0=inv(P0);

    S0=chol(P0,'lower');

    x=S0*randn(numDim,numSamp);

    maxVals=zeros(numSamp,1);
    for curSamp=1:numSamp
        xCur=x(:,curSamp);

        maxVal=-1;
        for curHyp=1:numHyp
            val=xCur'*PInvHyp(:,:,curHyp)*xCur;
            maxVal=max(val,maxVal);
        end
        maxVals(curSamp)=maxVal;
    end

    rMax=-1;
    rMin=Inf;
    for curSamp=1:numSamp
        xCur=x(:,curSamp);
        curRat=xCur'*P0Inv*xCur/maxVals(curSamp);
        rMax=max(curRat,rMax);
        rMin=min(curRat,rMin);
    end

    PSCI=P0./(u*rMin+(1-u)*rMax);
else%Given 1 or 0, don't change it.
    if(numHyp==0)
        xSCI=[];
        PSCI=[];
        return;
    end
    return
end

if(nargin>2&&~isempty(xHyp))
    xSCI=zeros(numDim,1);

    for curHyp=1:numHyp
        xSCI=xSCI+PInvHyp(:,:,curHyp)*xHyp(:,curHyp);
    end
    xSCI=P0*xSCI;
else
    xSCI=[];
end

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:
%Voice Based Biometric System
%By Ambavi K. Patel.

function w1=nanclr(w1)  % to replace  'not a number' vlues
[r1,c1]=size(w1);
for i=1:r1
    temp=isnan(w1(i,1:c1));
    for j=1:c1
        if temp(j)==1
            w1(i,j)=0.0001;
        else
        end;
    end;
end;


text:
classdef(Abstract) AbstractOptimizationVariable < matlab.mixin.SetGet & matlab.mixin.Heterogeneous
    %AbstractOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        id(1,1) double
    end

    methods
        x = getXsForVariable(obj)

        [lb, ub] = getBndsForVariable(obj)

        [lb, ub] = getAllBndsForVariable(obj)

        setBndsForVariable(obj, lb, ub)

        useTf = getUseTfForVariable(obj)

        setUseTfForVariable(obj, useTf)

        updateObjWithVarValue(obj, x)

        nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)

        function [xS, lbS, ubS] = getScaledXsForVariable(obj)
            x = obj.getXsForVariable();
            [lb, ub] = obj.getBndsForVariable();

            xS = x;
            lbS = lb;
            ubS = ub;
            for(i=1:length(x))
                xi = x(i);
                lbi = lb(i);
                ubi = ub(i);

                bndDiff = ubi - lbi;
                bndCenter = (lbi + ubi)/2;
                if(bndDiff > 1E-10)
                    xS(i) = (xi - bndCenter)/(bndDiff/2);
                    lbS(i) = -1;
                    ubS(i) = 1;
                else
                    xS(i) = xi;
                    lbS(i) = lbi;
                    ubS(i) = ubi;
                end
            end
        end

        function updateObjWithScaledVarValue(obj, xS)
            [lb, ub] = obj.getBndsForVariable();

            x = xS;
            for(i=1:length(xS))
                xSi = xS(i);
                lbi = lb(i);
                ubi = ub(i);

                bndDiff = ubi - lbi;
                bndCenter = (lbi + ubi)/2;

                if(bndDiff > 1E-10)
                    x(i) = xSi * (bndDiff/2) + bndCenter;
                else
                    x(i) = xSi;
                end
            end

            obj.updateObjWithVarValue(x);
        end

        function perturbVar(obj, pPct)
            x = obj.getXsForVariable();
            [lb, ub] = obj.getBndsForVariable();

            if(isempty(x))
                return;
            end

            for(i=1:length(x)) %#ok<*NO4LP>
                xi = x(i);
                lbi = lb(i);
                ubi = ub(i);

                fact = xi;
                if(abs(fact) < 1E-10)
                    fact = 1;
                end

                p1 = xi - (pPct/100)*fact;
                p2 = xi + (pPct/100)*fact;

                lbRnd = max(min(p1,p2),lbi);
                ubRnd = min(max(p1,p2),ubi);

                xRnd = lbRnd + (ubRnd - lbRnd)*rand();
                x(i) = xRnd;
            end

            obj.updateObjWithVarValue(x);
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            useTf = obj.getUseTfForVariable();
            varsStoredInRad = false(size(useTf));
        end
    end

    methods(Sealed)
        function numVars = getNumOfVars(obj)
            [lb,~] = obj.getBndsForVariable();
            numVars = numel(lb);
        end

        function tf = eq(A,B)
            tf = [A.id] == [B.id];
        end

        function tf = ne(A,B)
            tf = [A.id] ~= [B.id];
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% std_detachplots() -  Given a figure with subplots and several lines per axis, will add a callback to
% each axis specified in the 'figtitles' input. The callback consist in a figure with all the detached
% individuals lines.
%
% Usage:
%   >>   std_detachplots('','','data',data 'figtitles', alltitlestmp,'sbtitles',sbtitles,'handles', handles);
%
% Inputs:
%      data        - Cell array containing the data matrices for each plot in the same order showed in the figure
%      figtitles   - Cell array of the titles of each individual axes in
%                    the  figure. The titles must correspond. The function
%                    use this value to find the right hanlde of the axis
%      sbtitles    - Cell array of cell arrays with the titles for each
%                    detached line per axis. i.e. {{'Axis1 line1' 'Axis1 line2'} {'Axis2 line1' 'Axis2 line2'}}
%      handles     - Handles of the main figure who contain all the
%                    subplots
%      flagstd     - Flag to plot the Standar Deviation  {default: 1} means 'on'
%
% See also:
%
% Author: Ramon Martinez-Cancino, SCCN, 2014
%
% Copyright (C) 2014  Ramon Martinez-Cancino,INC, SCCN
%
% 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 std_detachplots(hObject,eventdata,varargin)

% display help if not enough arguments
if nargin < 2
    help std_detachplots;
    return;
end
icadefs;
try
    options = varargin;
    if ~isempty( varargin ),
        for i = 1:2:numel(options)
            g.(options{i}) = options{i+1};
        end
    else g= []; end;
catch
    disp('std_detachplots() error: calling convention {''key'', value, ... } error'); return;
end;

try g.data;         catch, g.data          =   [];  end; % Name of plots
try g.figtitles;    catch, g.figtitles     =   [];  end; % Name of plots
try g.sbtitles;     catch, g.sbtitles      =   [];  end; % Name of plots
try g.handles;      catch, g.handles       =   [];  end; % Handles of figure
try g.flagstd;      catch, g.flagstd       =   1;   end; % Plot std band around mean
try g.xlabel;       catch, g.xlabel        =   '';  end; % xlabel
try g.ylabel;       catch, g.ylabel        =   '';  end; % ylabel
try g.timevec;      catch, g.timevec       =   '';  end; % Time or freq vector
try g.filter;       catch, g.filter        =   '';  end; % Low pass filter freq

% Checking data
if isempty(g.handles) && any([isempty(g.data) isempty(g.figtitles)])
    error('std_detachplots : Check entries. Options ''handles'', ''data'' and ''figtitles'' must be provided');
end

if iscell(g.data)
    nplots = numel(g.data(:));
else
    nplots = 1 ;
    g.data = {g.data};
end

%  Checking sbtitles
if isempty(g.sbtitles)
    for i = 1:numel(g.data)
        c = 1;
        if ~isempty(g.data{i})
            for j = 1:size(g.data{i},2)
                g.sbtitles{i}{j} = {['Line ' num2str(c)]};
                c = c+1;
            end
        end
    end
end

% Plot goes here
%--------------------------------------------------------------------------
if isempty(g.handles)
    for i_nplots = 1 : nplots
        idata = g.data{i_nplots};

        % Filtering data to be plotted
        if ~isempty(g.filter), idata = myfilt(idata, 1000/(g.timevec(2)-g.timevec(1)), 0, g.filter); end;

        len = size(idata,2);
        if len > 0 % A non-empty cluster
            % Getting the mean
            meandata = mean(idata,2);
            stddata  = std(idata,0,2);

            if license('checkout', 'statistics_toolbox')
                SEM = std(idata,0,2)/sqrt(size(idata,2));       % Standard Error
                ts = tinv([0.025  0.975],size(idata,2)-1);      % T-Score
                lower    = meandata + ts(1)*SEM;                % CI
                upper    = meandata + ts(2)*SEM;                % CI
            else
                lower    = meandata-2*stddata;
                upper    = meandata+2*stddata;
            end

            hplot = figure('name', g.figtitles, 'NumberTitle','off');
            rowcols(2) = ceil(sqrt(len + 4));
            rowcols(1) = ceil((len+4)/rowcols(2));

            for k = 1:len
                %--- first sbplot row ----
                if k <= rowcols(2) - 2
                    figure(hplot);
                    sbplot(rowcols(1),rowcols(2),k+2);
                    hold on;
                    plotlines(k,idata,meandata,lower, upper,g.sbtitles{k},g);
                else
                    figure(hplot)
                    sbplot(rowcols(1),rowcols(2),k+4);
                    hold on;
                    plotlines(k,idata,meandata,lower, upper,g.sbtitles{k},g);
                end
            end
            % Plot all figure
            figure(hplot)
            sbplot(rowcols(1),rowcols(2),[1 rowcols(2)+2 ]);
            hold on;
            plotlines(1:len,idata,meandata,lower,upper,g.figtitles,g);
            set(gcf,'Color', BACKCOLOR);
            orient tall;
        end
    end
else
    xlabelval = '';
    ylabelval = '';
    % Match Children handles based on titles provided
    c = 0;
    for i = 1: nplots
        htemp          = findall(g.handles,'String', g.figtitles{i});
        if all([~isempty(htemp) ~isempty(g.data(i))])
            handlestemp{i} = htemp(1);
            % Getting x label from handles
            if isempty(xlabelval)
                xlabelval = get(get(get(handlestemp{i},'Parent'),'Xlabel'),'String'); % xlabelval = handlestemp{i}.Parent.XLabel.String;
            end
            % Getting y label from handles
            if isempty(ylabelval)
                ylabelval = get(get(get(handlestemp{i},'Parent'),'Ylabel'),'String'); % ylabelval = handlestemp{i}.Parent.YLabel.String;
            end
            % Getting timevec from handles
            if isempty(g.timevec)
                tmp = get(get(get(handlestemp{i},'Parent'),'Children'));
                g.timevec = tmp(i).XData; % g.timevec = handlestemp{i}.Parent.Children(i).XData;
            end
        else
            handlestemp{i} = [];
        end

        % Callback setting
        if ~isempty(handlestemp{i})
            c = c + 1;
            % For Axis
            set( get( handlestemp{i}, 'Parent'), 'ButtonDownFcn',{@std_detachplots,...
                'data'     ,g.data{i},...
                'timevec'  ,g.timevec,...
                'handles'  ,'',...
                'figtitles',g.figtitles{i},...
                'sbtitles' ,g.sbtitles{c},...
                'xlabel'   ,xlabelval,...
                'ylabel'   ,ylabelval,...
                'filter'   , g.filter,...
                });
            % For lines
            set( get(get( handlestemp{i}, 'Parent'), 'Children'), 'ButtonDownFcn',{@std_detachplots,...
                'data'     ,g.data{i},...
                'timevec'  ,g.timevec,...
                'handles'  ,'',...
                'figtitles',g.figtitles{i},...
                'sbtitles' ,g.sbtitles{c},...
                'xlabel'   ,xlabelval,...
                'ylabel'   ,ylabelval,...
                'filter'   , g.filter,...
                });
        end
    end
end

function plotlines(kindx,idata,meandata,lower,upper,sbtitle,g)
for i = 1:length(kindx)
    plot(g.timevec,idata(:,kindx(i)),'b','LineWidth', 0.1); % plot g.data
end
plot(g.timevec,meandata,'r','LineWidth', 0.1);              % plot mean

if g.flagstd                                                % plot std band
    eeglabciplot(lower,upper,g.timevec, 'r', 0.2);
    axis tight;
end

if length(kindx)>1
    xlabel(g.xlabel);
    ylabel(g.ylabel);
end
box on;
grid on;
axis tight;
title(sbtitle, 'interpreter', 'none');

% rapid filtering for ERP (from std_plotcurve)
% -----------------------
function tmpdata2 = myfilt(tmpdata, srate, lowpass, highpass);
    bscorrect = 1;
    if bscorrect
        % Getting initial baseline
        bs_val1  =  mean(tmpdata,1);
        bs1      = repmat(bs_val1, size(tmpdata,1), 1);
    end

    % Filtering
    tmpdata2 = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3)*size(tmpdata,4));
    tmpdata2 = eegfiltfft(tmpdata2',srate, lowpass, highpass)';
    tmpdata2 = reshape(tmpdata2, size(tmpdata,1), size(tmpdata,2), size(tmpdata,3), size(tmpdata,4));

    if bscorrect
        % Getting after-filter baseline
        bs_val2  =  mean(tmpdata2,1);
        bs2      = repmat(bs_val2, size(tmpdata2,1), 1);

        % Correcting the baseline
        realbs = bs1-bs2;
        tmpdata2 = tmpdata2 + realbs;
    end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


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

% FT_NIRS_TRANSFORM_ODs computes the transformation from optical densities (OD)
% to chromophore concentrations such as (de-)oxygenated  hemoglobin, or the
% other way around.
%
% Use as either
%   [data]      = ft_nirs_transform_ODs(cfg, data)
%
% The configuration "cfg" is a structure containing information about
% target of the transformation. The configuration should contain
%   cfg.channel            = Nx1 cell-array with selection of channels
%                            (default = 'nirs'), see FT_CHANNELSELECTION for
%                            more details
%   cfg.target             = Mx1 cell-array, can be 'O2Hb' (oxygenated hemo-
%                            globin), 'HHb' de-oxygenated hemoglobin') or
%                            'tHb' (total hemoglobin), or a combination of
%                            those (default: {'O2Hb', 'HHb'})
%
% Optional configuration settings are
%   cfg.age                = scalar, age of the subject (necessary to
%                            automatically select the appropriate DPF, or
%   cfg.dpf                = scalar, differential path length factor
%   cfg.dpffile            = string, location to a lookup table for the
%                            relation between participant age and DPF
%
% Note that the DPF might be different across channels, and is usually
% stored in the optode structure contained in the data.
%
% The function returns data transformed to the specified chromophore
% concentrations that were requested.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure.
%
% See also FT_NIRS_PREPARE_ODTRANSFORMATION, FT_APPLY_MONTAGE

% Options to be implemented:
%   cfg.opto          = structure with optode definition or filename, see FT_READ_SENS
%   cfg.siunits       = yes/no, ensure that SI units are used consistently
%   cfg.logarithm     = string, can be 'natural' or 'base10' (default = 'base10')
%   cfg.inverse       = string, whether the multiplicative inverse should be take
%
% Note that HomER uses the inverse, natural logarithm. The inverse is taken
% after the ln transform.

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2015-2016 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com
%
% Main programmer: Jörn M. Horschig
% $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 loadvar     data
ft_preamble provenance  data

ft_preamble debug

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes');

% set the defaults
cfg.channel = ft_getopt(cfg, 'channel', 'nirs');
cfg.target  = ft_getopt(cfg, 'target', {'O2Hb', 'HHb'});

% determine the NIRS channels
cfg.channel = ft_channelselection(cfg.channel, data);

% make a copy of the data
origdata = data;

% keep only the NIRS channels
data = ft_selectdata(cfg, data);

% also keep the non-NIRS channels
tmpcfg = cfg;
tmpcfg.channel = setdiff(origdata.label, cfg.channel);
dataextra = ft_selectdata(tmpcfg, origdata);

target = cfg.target;
if ~iscell(target)
  target     = {target};
end

computetHb  = any(ismember(target, 'tHb'));
computeHHb  = any(ismember(target, 'HHb'));
computeO2Hb = any(ismember(target, 'O2Hb'));

% cfg-handling is done inside here
[montage, cfg] = ft_nirs_prepare_ODtransformation(cfg, data);

% save montage in the cfg
cfg.montage = montage;

% apply the (combined) montages
dataout = ft_apply_montage(data, montage, 'keepunused', 'yes');

if computetHb
  % total hemoglobin is the sum of oxygenated and deoxygenated hemoglobin
  tmpcfg = [];
  tmpcfg.channel = '*[O2Hb]';
  dataO2Hb = ft_selectdata(tmpcfg, dataout);
  dataO2Hb.label = strrep(dataO2Hb.label, 'O2Hb', 'tHb');

  tmpcfg = [];
  tmpcfg.channel = '*[HHb]';
  dataHHb = ft_selectdata(tmpcfg, dataout);
  dataHHb.label = strrep(dataHHb.label, 'HHb', 'tHb');

  tmpcfg = [];
  tmpcfg.operation = 'add';
  tmpcfg.parameter = 'trial';
  dataTotal = ft_math(tmpcfg, dataO2Hb, dataHHb);

  if ~computeHHb && ~computeO2Hb
    % replace dataout
    dataout.label = dataTotal.label;
    dataout.trial = dataTotal.trial;
  else
    % concat to dataout
    dataout.label = vertcat(dataout.label, dataTotal.label);
    dataout.trial = cellfun(@vertcat, dataout.trial, dataTotal.trial, 'UniformOutput', false);
  end
end

% remove O2 or H if not desired
tmpcfg = [];
tmpcfg.channel = [];
if computetHb
  tmpcfg.channel = horzcat(tmpcfg.channel, {'*tHb]'});
end
if computeHHb
  tmpcfg.channel =  horzcat(tmpcfg.channel, {'*HHb]'});
end
if computeO2Hb
  tmpcfg.channel =  horzcat(tmpcfg.channel, {'*O2Hb]'});
end

dataout = ft_selectdata(tmpcfg, dataout);

if length(dataextra.label)>0
  % append the extra data
  data = ft_appenddata([], dataout, dataextra);
else
  data = dataout;
end

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

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

ft_postamble previous    data
ft_postamble provenance  data
ft_postamble history     data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pot = convert_to_pot(CPD, pot_type, domain, evidence)
% CONVERT_TO_POT Convert a tabular CPD to one or more potentials
% pot = convert_to_pot(CPD, pot_type, domain, evidence)

% This is the same as discrete_CPD/convert_to_pot,
% except we didn't want to the kernel to inherit methods like sample_node etc.

sz = CPD.sz;
ns = zeros(1, max(domain));
ns(domain) = sz;

odom = domain(~isemptycell(evidence(domain)));
T = convert_to_table(CPD, domain, evidence);

switch pot_type
 case 'u',
  pot = upot(domain, sz, T, 0*myones(sz));
 case 'd',
  ns(odom) = 1;
  pot = dpot(domain, ns(domain), T);
 case 'c',
  % Since we want the output to be a Gaussian, the whole family must be observed.
  % In other words, the potential is really just a constant.
  p = T.p;
  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));
  ns(domain) = 0;
  pot = cpot(domain, ns(domain), log(p));
 case 'cg',
  T = T(:);
  ns(odom) = 1;
  can = cell(1, length(T));
  for i=1:length(T)
    can{i} = cpot([], [], log(T(i)));
  end
  pot = cgpot(domain, [], ns, can);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function x_new = r8sd_cg ( n, ndiag, offset, a, b, x )

%*****************************************************************************80
%
%% R8SD_CG uses the conjugate gradient method on a R8SD linear system.
%
%  Discussion:
%
%    The R8SD storage format is for symmetric matrices whose only nonzero entries
%    occur along a few diagonals, but for which these diagonals are not all
%    close enough to the main diagonal for band storage to be efficient.
%
%    In that case, we assign the main diagonal the offset value 0, and
%    each successive superdiagonal gets an offset value 1 higher, until
%    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
%
%    Assuming there are NDIAG nonzero diagonals (ignoring subdiagonals!),
%    we then create an array B that has N rows and NDIAG columns, and simply
%    "collapse" the matrix A to the left:
%
%    For the conjugate gradient method to be applicable, the matrix A must
%    be a positive definite symmetric matrix.
%
%    The method is designed to reach the solution to the linear system
%      A * x = b
%    after N computational steps.  However, roundoff may introduce
%    unacceptably large errors for some problems.  In such a case,
%    calling the routine a second time, using the current solution estimate
%    as the new starting guess, should result in improved results.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer NDIAG, the number of diagonals that are stored.
%    NDIAG must be at least 1 and no more than N.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal storage.
%
%    Input, real A(N,NDIAG), the R8SD matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, real X(N), an estimate for the solution, which may be 0.
%
%    Output, real X_NEW(N), the approximate solution vector.  Note that
%    repeated calls to this routine, with the output X_NEW from the
%    previous call used as the input value of X, MAY improve the solution.
%
  x_new(1:n) = x(1:n);
%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap(1:n) = r8sd_mxv ( n, ndiag, offset, a, x_new );

  r(1:n) = b(1:n) - ap(1:n);
  p(1:n) = b(1:n) - ap(1:n);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP = A*P.
%
    ap(1:n) = r8sd_mxv ( n, ndiag, offset, a, p );
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p(1:n) * ap(1:n)';
    pr =  p(1:n) * r(1:n)';

    if ( pap == 0.0E+00 )
      return;
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x_new(1:n) = x_new(1:n) + alpha * p(1:n);
    r(1:n) = r(1:n) - alpha * ap(1:n);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r(1:n) * ap(1:n)';

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n) = r(1:n) + beta * p(1:n);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% predicting numerical values using Discreminate Analysis

clear all;

disp('===REgress Tree===');
disp('Reading featur vector');

featurs = csvread('forWeka_featuresonly.csv');

num_data = 5000;
size_training = floor(.6*num_data);

trainingset = featurs(1:size_training,:);
testset = featurs((size_training+1):num_data,:);

disp('Splitting up data into training/test sets');
[num,txt,raw] = xlsread('C:\MatlabNLP\examples\gsa\data\final104.xls');

% reading the description of each shoe
descriptions = raw(2:size(raw,1),2);
style_ratings = num(1:size(num,1),1);
comfort_ratings = num(1:size(num,1),4);
overal_ratings = num(1:size(num,1),5);

% only take m data points
m=num_data;
descriptions = descriptions(1:m);
style_ratings = style_ratings(1:m);
comfort_ratings = comfort_ratings(1:m);
overal_ratings = overal_ratings(1:m);

responsevals = [style_ratings, comfort_ratings, overal_ratings];

responsevals_training = responsevals(1:size_training,:);
responsevals_test = responsevals((size_training+1):num_data,:);
disp('Multiple linear regression');
% http://www.mathworks.com/help/toolbox/stats/classregtree.html

tic;

a = responsevals_training(:,3)';
b = responsevals_test(:,3)';

class = classregtree(trainingset,a');
yfit = round(eval(class,testset));

cMat2 = confusionmat(yfit,b)
toc;


text:
function [f] = elec_fit_sphere_optim(r, X, Y, Z, xo, yo, zo)

% elec_fit_sphere_optim - Optimization for elec_fit_sphere.m
%
% Called from elec_fit_sphere.m
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no implied or express warranties
% History:  02/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% with center (Xo,Yo,Zo) and radius r, the equation of a sphere is:
%
% r^2 = (x-xo)^2  +  (y-yo)^2  +  (z-zo)^2
%
% This function below creates a scalar value to
% return to the fminsearch function in elec_fit_sphere.

S = (X-xo).^2  +  (Y-yo).^2  +  (Z-zo).^2  -  r^2;

f = sum( S.^2 );

author: PatternRecognition

repo: OpenBMI


text:
function y = acoth(x)
%ACOTH        Implements  acoth(x)  for intervals
%
%   y = acoth(x)
%
%interval standard function implementation
%

% written  10/16/98     S.M. Rump
% modified 12/30/98     S.M. Rump  improved speed and use atanh
% modified 08/31/99     S.M. Rump  complex allowed, sparse input,
%                                  major revision, improved accuracy near 1
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 12/04/05     S.M. Rump  'realstdfctsexcptnignore' added and
%                                     some improvements, extreme values
%                                     for approximate part
% modified 09/06/07     S.M. Rump  approximate std fcts removed, exceptional arguments
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/20/08     S.M. Rump  check for zero omitted
% modified 10/18/08     S.M. Rump  StdFctsException ignore/NaN
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  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
    if issparse(x.mid)
      x.mid = full(x.mid);
      x.rad = full(x.rad);
    end
    y = atanh( 1./x );
    if rndold
      setround(rndold)
    end
    return
  end

  if issparse(x.inf)
    x.inf = full(x.inf);
    x.sup = full(x.sup);
  end
  % input x real and full
  % real range of definition:  [-inf,-1] and [1,inf]
  % take care for intersection( x , (-1,1) ) nonempty
  INTLAB_STDFCTS_EXCPTN = getappdata(0,'INTLAB_STDFCTS_EXCPTN');
  Index1 = ( abs(x.inf)<1 );           % (partially) exceptional indices
  Index2 = ( abs(x.sup)<1 );
  if ~isempty(find(Index1)) | ~isempty(find(Index2))   % handle input out-of-range
    if INTLAB_STDFCTS_EXCPTN<=1 % out-of-range input handled as complex
      if INTLAB_STDFCTS_EXCPTN==1
        warning('ACOTH: Real interval input out of range changed to be complex')
      end
      y = x;
      index = Index1 | Index2;
      %VVVV  y(index) = atanh(cintval(1./x(index)));
      s.type = '()'; s.subs = {index}; y = subsasgn(y,s,atanh(cintval(1./subsref(x,s))));
      %AAAA  Matlab bug fix
      index = ~index;
      if any(index(:))
        %VVVV  y(index) = atanh(1./x(index));
        s.type = '()'; s.subs = {index}; y = subsasgn(y,s,atanh(1./subsref(x,s)));
        %AAAA  Matlab bug fix
      end
      index = ( x==0 );
      if any(index(:))
        INTLAB_STDFCTS_PI = getappdata(0,'INTLAB_STDFCTS_PI');
        iPI2 = intval(j*INTLAB_STDFCTS_PI.PI2MID,INTLAB_STDFCTS_PI.PI2RAD,'midrad');
        %VVVV  y(index) = iPI2;
        s.type = '()'; s.subs = {index}; y = subsasgn(y,s,iPI2);
        %AAAA  Matlab bug fix
      end
      if rndold
        setround(rndold)
      end
      return
    end
    setappdata(0,'INTLAB_STDFCTS_EXCPTN_',1);
    if INTLAB_STDFCTS_EXCPTN==3    % ignore input out of range (ignore-mode)
      x.inf(Index1) = 1;                  % completely exceptional indices treated below
      x.sup(Index2) = -1;
      Index = ( Index1 & Index2 );        % completely exceptional indices
    end
  else
    Index = [];                           % make sure Index is not undefined
  end

  % input x real and full
  y = x;
  wng = warning;                         % get current warning mode
  warning off

  % treat positive intervals
  index1 = ( x.inf>0 );
  if any(index1(:))
    y.inf(index1) = acoth_pos(x.sup(index1),-1);
    y.sup(index1) = acoth_pos(x.inf(index1),1);
  end

  % treat negative intervals
  index1 = ( x.sup<0 );
  if any(index1(:))
    y.inf(index1) = - acoth_pos(-x.sup(index1),1);
    y.sup(index1) = - acoth_pos(-x.inf(index1),-1);
  end

  if INTLAB_STDFCTS_EXCPTN==3      % ignore input out of range (ignore-mode)
    if ~isempty(find(Index))              % completely exceptional arguments to NaN
      y.inf(Index) = NaN;
      y.sup(Index) = NaN;
    end
  else                                    % any input out of range to NaN (NaN-mode)
    Index = Index1 | Index2;
    if ~isempty(find(Index))              % exceptional arguments to NaN
      y.inf(Index) = NaN;
      y.sup(Index) = NaN;
    end
  end

  setround(rndold)
  warning(wng)                           % restore warning mode

function y = acoth_pos(x,rnd)
% local acoth for double vector x>=1 with rounding corresponding to rnd
%

  y = x;

  index = ( x<4 );
  if any(index(:))             % 1 <= x <= 4
    setround(rnd)              % acoth(x) = log( 1 + 2/(x-1) )
    e = x(index) - 1;          % e w/o rounding error
    e = 1 + 2./e;
    y(index) = log_rnd( e , rnd ) / 2;
    y(x==1) = inf;
  end

  index = ~index;              % x >= 4, difference x-1 not necessarily exact
  if any(index(:))             % acoth(x) = atanh(1/x)
    setround(rnd)
    e = 1./x(index);
    y(index) = atanh_pos( e , rnd );
  end

repo: OpenMAS


text:
function mono = increasing_except_at(xL,xU,p)

if xU <= p || xL >= p
    mono = 'increasing';
else
    mono = 'none';
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef FixedPointInFrame < AbstractGeometricPoint
    %FixedPointInFrame Summary of this class goes here
    %   Detailed explanation goes here

    properties
        cartElem(1,1) CartesianElementSet

        name(1,:) char

        lvdData LvdData

        %marker
        markerColor(1,1) ColorSpecEnum = ColorSpecEnum.Red;
        markerShape(1,1) MarkerStyleEnum = MarkerStyleEnum.RightTriangle;

        %track line
        plotTrkLine(1,1) logical = true;
        trkLineColor(1,1) ColorSpecEnum = ColorSpecEnum.Black;
        trkLineSpec(1,1) LineSpecEnum = LineSpecEnum.DottedLine;
    end

    methods
        function obj = FixedPointInFrame(rVect, frame, name, lvdData)
            obj.cartElem = CartesianElementSet(0, rVect(:), [0;0;0], frame);
            obj.name = name;
            obj.lvdData = lvdData;
        end

        function newCartElem = getPositionAtTime(obj, time, ~, inFrame)
            newCartElem = repmat(CartesianElementSet.getDefaultElements(), [1 length(time)]);
            for(i=1:length(time))
                newCartElem(i) = copy(obj.cartElem);
                newCartElem(i).time = time(i);
            end
            newCartElem = convertToFrame(newCartElem, inFrame);
        end

        function name = getName(obj)
            name = obj.name;
        end

        function setName(obj, name)
            obj.name = name;
        end

        function listboxStr = getListboxStr(obj)
            listboxStr = sprintf('%s (Fixed in Frame: %s)', obj.getName(), obj.cartElem.frame.getNameStr());
        end

        function useTf = openEditDialog(obj, ~)
%             useTf = lvd_EditFixedInFramePointGUI(obj, obj.lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditFixedInFramePointGUI_App(obj, obj.lvdData, output);
            useTf = output.output{1};
        end

        function tf = isVehDependent(obj)
            tf = false;
        end

        function tf = canBePlotted(obj)
            tf = true;
        end

        function bodyInfo = getOriginBody(obj)
            bodyInfo = obj.cartElem.frame.getOriginBody();
        end

        function tf = usesGroundObj(~, ~)
            tf = false;
        end

        function tf = usesGeometricPoint(~, ~)
            tf = false;
        end

        function tf = usesGeometricVector(~, ~)
            tf = false;
        end

        function tf = usesGeometricCoordSys(~, ~)
            tf = false;
        end

        function tf = usesGeometricRefFrame(~, ~)
            tf = false;
        end

        function tf = usesGeometricAngle(~, ~)
            tf = false;
        end

        function tf = usesGeometricPlane(~, ~)
            tf = false;
        end

        function tf = isInUse(obj, lvdData)
            tf = lvdData.usesGeometricPoint(obj);
        end

        function rVect = getRVect(obj)
            rVect = obj.cartElem.rVect;
        end

        function setRVect(obj, rVect)
            obj.cartElem.rVect = rVect;
        end

        function frame = getFrame(obj)
            frame = obj.cartElem.frame;
        end
    end
end


text:
function matlab_random_parallel_test ( )

%*****************************************************************************80
%
%% MATLAB_RANDOM_PARALLEL_TEST looks at random number generation in parallel Matlab.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MATLAB_RANDOM_PARALLEL_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Look at random number generation in parallel Matlab.\n' );

  matlab_random_parallel_test01 ( );
  matlab_random_parallel_test02 ( );
  matlab_random_parallel_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MATLAB_RANDOM_PARALLEL_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function noise = noiseTest(noiseType);

% NOISETEST Run some tests on the specified noise model.

% NOISE

if iscell(noiseType)
  % compound noise type
  noise.type = 'cmpnd';
  for i = 1:length(noiseType)
    noise.comp{i}.type = noiseType{i};
  end
else
  noise.type = noiseType;
end
noise.numProcess = 3;
numData = 20;

noise.C = 10;
noise.numData = numData;
noise = noiseParamInit(noise);
% Set the parameters randomly.
params = noiseExtractParam(noise);
params = randn(size(params))./sqrt(randn(size(params)).^2);
noise = noiseExpandParam(noise, params);

mu = randn(numData, noise.numProcess).*sqrt(1./(randn(numData,noise.numProcess).^2));

varsigma1 = 1./(randn(numData, noise.numProcess).^2);
varsigma2 = (randn(numData, noise.numProcess).^2);
varsigmaSwitch = rand(numData, noise.numProcess)> 0.5;
varsigma = varsigmaSwitch.*varsigma1 + (1-varsigmaSwitch).*varsigma2;
y = noiseOut(noise, mu, varsigma);

mu = randn(numData, noise.numProcess).*sqrt(1./(randn(numData,noise.numProcess).^2));

varsigma1 = 1./(randn(numData, noise.numProcess).^2);
varsigma2 = (randn(numData, noise.numProcess).^2);
varsigmaSwitch = rand(numData, noise.numProcess)> 0.5;
varsigma = varsigmaSwitch.*varsigma1 + (1-varsigmaSwitch).*varsigma2;

% Test for missing variables
if noise.missing
  index = randperm(size(y, 1)*size(y, 2));
  index = index(1:3);
  y(index) = NaN;
end

fprintf('mu values\n');
disp(mu)
fprintf('varsigma values\n');
disp(varsigma)
fprintf('y values\n');
disp(y)
epsilon = 1e-6;
params = noiseExtractParam(noise);
origParams = params;
for i = 1:length(params);
  params = origParams;
  params(i) = origParams(i) + epsilon;
  noise = noiseExpandParam(noise, params);
  Lplus(i) = noiseLogLikelihood(noise, mu, varsigma, y);
  params(i) = origParams(i) - epsilon;
  noise = noiseExpandParam(noise, params);
  Lminus(i) = noiseLogLikelihood(noise, mu, varsigma, y);
end
params = origParams;
noise = noiseExpandParam(noise, params);
[void, names] = noiseExtractParam(noise);
gLDiff = .5*(Lplus - Lminus)/epsilon;
g = noiseGradientParam(noise, mu, varsigma, y);

paramMaxDiff = max(max(abs(gLDiff-g)));
if paramMaxDiff > 2*epsilon
  l = 0;
  for i = 1:length(names)
    if l < length(names{i})
      l = length(names{i});
    end
  end

  fprintf([char(repmat(32, 1, l)) '\tanalytic   diffs     delta\n']);
  for i = 1:length(names)
    spaceLen = l - length(names{i});
    space = char(repmat(32, 1, spaceLen));
    fprintf([space names{i} ':\t%4.6f\t%4.6f\t%4.6f\n'], ...
            g(i), gLDiff(i), gLDiff(i) - g(i));
  end
end
Lplus = zeros(size(mu));
Lminus = zeros(size(mu));
origMu = mu;
for i = 1:size(mu, 1)
  for j = 1:size(mu, 2)
    mu = origMu;
    mu(i, j) = origMu(i, j) + epsilon;
    Lplus(i, j) = noiseLogLikelihood(noise, mu, varsigma, y);
    mu(i, j) = origMu(i, j) - epsilon;
    Lminus(i, j) = noiseLogLikelihood(noise, mu, varsigma, y);
  end
end
mu = origMu;
gMuDiff = .5*(Lplus - Lminus)/epsilon;

Lplus = zeros(size(varsigma));
Lminus = zeros(size(varsigma));
origVarsigma = varsigma;
for i = 1:size(varsigma, 1)
  for j = 1:size(varsigma, 2)
    varsigma = origVarsigma;
    varsigma(i, j) = origVarsigma(i, j) + epsilon;
    Lplus(i, j) = noiseLogLikelihood(noise, mu, varsigma, y);
    varsigma(i, j) = origVarsigma(i, j) - epsilon;
    Lminus(i, j) = noiseLogLikelihood(noise, mu, varsigma, y);
  end
end
varsigma = origVarsigma;
gVarsigmaDiff = .5*(Lplus - Lminus)/epsilon;

[g, gvs] = noiseGradVals(noise, mu, varsigma, y);

vsMaxDiff = max(max(abs(gvs-gVarsigmaDiff)));
muMaxDiff = max(max(abs(g-gMuDiff)));

if vsMaxDiff > 1e-7
  fprintf('y\n')
  disp(y)
  fprintf('varsigma\n')
  disp(varsigma)
  fprintf('mu\n')
  disp(mu)
  fprintf('gvs\n')
  disp(gvs)
  fprintf('diffs vs\n')
  disp(gVarsigmaDiff)
  fprintf('gvs -diffvgs\n')
  disp(gvs -gVarsigmaDiff)
end
fprintf('Param max diff: %2.6f.\n', paramMaxDiff)
fprintf('Mu max diff: %2.6f.\n', muMaxDiff)
fprintf('Varsigma max diff: %2.6f.\n', vsMaxDiff)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function ind = wsddnVOChash_lookup(hash,s)
% From the PASCAL VOC 2011 devkit

hsize=numel(hash.key);
h=mod(str2double(s([4 6:end])),hsize)+1;
ind=hash.val{h}(strmatch(s,hash.key{h},'exact'));

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
%%%%%%%%%%%%%%%%%%%%%%%
% Compute features for a set of video files from datasets
%
close all;
clear;
warning('off','all');
addpath(genpath('./lib/iqm/BRISQUE_release'));

%==================================================
% parameters
algo_name = 'BRISQUE_feat_sel';
data_name = 'KONVID_1K';

%% *You need to customize here*
if strcmp(data_name, 'TEST_VIDEOS')
    data_path = 'videos'; % dataset video path
elseif strcmp(data_name, 'KONVID_1K')
    data_path = '/media/ztu/Seagate-ztu-ugc/KONVID_1K/KoNViD_1k_videos';
elseif strcmp(data_name, 'LIVE_VQC')
    data_path = '/media/ztu/Seagate-ztu-ugc/LIVE_VQC/VideoDatabase';
elseif strcmp(data_name, 'YOUTUBE_UGC')
    data_path = '/media/ztu/Seagate-ztu-ugc/YT_UGC/original_videos';
elseif strcmp(data_name, 'LIVE_VQA')
    data_path = '/media/ztu/Seagate-ztu/LIVE_VQA/videos';
end

video_tmp = 'video_tmp';
if ~exist(video_tmp, 'dir'), mkdir(video_tmp); end
feat_path = '../features';
filelist_csv = fullfile(feat_path, [data_name,'_metadata.csv']);
filelist = readtable(filelist_csv);
num_videos = size(filelist,1);
out_path = './feat_sel_mats';
if ~exist(out_path, 'dir'), mkdir(out_path); end
out_feat_name = fullfile(out_path, [data_name,'_',algo_name,'_feats.mat']);
out_feat_frames_name = fullfile(out_path, [data_name,'_',algo_name,'_frames_feats.mat']);
feats_mat = [];
feats_mat_frames = cell( num_videos, 1 );

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

tic
% parallel
for i = 1:num_videos
    % get video full path and decoded video name
	if strcmp(data_name, 'TEST_VIDEOS')
	    video_name = fullfile(data_path,  filelist.video_name{i});
	    yuv_name = fullfile(video_tmp, [filelist.video_name{i}, '.yuv']);
	elseif strcmp(data_name, 'KONVID_1K')
	    video_name = fullfile(data_path,  [num2str(filelist.flickr_id(i)),'.mp4']);
	    yuv_name = fullfile(video_tmp, [num2str(filelist.flickr_id(i)), '.yuv']);
	elseif strcmp(data_name, 'LIVE_VQC')
	    video_name = fullfile(data_path, filelist.File{i});
	    yuv_name = fullfile(video_tmp, [filelist.File{i}, '.yuv']);
	elseif strcmp(data_name, 'YOUTUBE_UGC')
	    video_name = fullfile(data_path, filelist.category{i},...
		[num2str(filelist.resolution(i)),'P'],[filelist.vid{i},'.mkv']);
	    yuv_name = fullfile(video_tmp, [filelist.vid{i}, '.yuv']);
	elseif strcmp(data_name, 'LIVE_VQA')
	    strs = strsplit(filelist.filename{i}, '_');
	    video_name = fullfile(data_path, [strs{1}(1:2), '_Folder'], filelist.filename{i});
	    yuv_name = video_name;
	end
	fprintf('\n---\nComputing features for %d-th sequence: %s\n', i, video_name);
    % decode video and store in video_tmp dir
    if ~strcmp(video_name, yuv_name)
    cmd = ['ffmpeg -loglevel error -y -i ', video_name, ' -pix_fmt yuv420p -vsync 0 ', yuv_name];
    system(cmd);
    end

    % get video meta data
    width = filelist.width(i);
    height = filelist.height(i);
    framerate = round(filelist.framerate(i));
    nb_frames = filelist.nb_frames(i);

    % read YUV frame (credit: Dae Yeol Lee)
    fp_input = fopen(yuv_name, 'r');
    uv_width = width/2;
    uv_height = height/2;
    feats_frames = [];

    tic
    % calculate every frame
    for fr = 1:2:nb_frames-5

        try
            %% Start a file pointer
            fseek(fp_input,(fr-1)*1.5*width*height, 'bof'); % Frame read for 8 bit
            %% Y component
            %1) read y stream
            y_stream = fread(fp_input, width * height, 'uchar'); % for 8 bit
            % 2) reshape into a plane
            y_plane= reshape(y_stream, width, height).';
%             imshow(y_plane,[]);
            feats_frames(end+1,:) = brisque_feature(y_plane);
        catch
            continue
        end
    end
    fclose(fp_input);
    delete(yuv_name)
    feats_mat_frames{i} = feats_frames;

    % compute brisque mean and consistency within each 1-sec chunk!!
    cons_feats = [];
    mean_feats = [];
    n_temp_vecs = length(feats_frames(:,1));
    half_blk_len = floor(framerate/2);

    fprintf('Pooling mean and consistency features\n');

    for j = 1:half_blk_len:n_temp_vecs - half_blk_len

        j_start = j;
        j_end = j + half_blk_len;

        % compute consistency features
        cons_feats = [cons_feats; std(feats_frames(j_start:j_end, :))];
        mean_feats = [mean_feats; mean(feats_frames(j_start:j_end, :))];
    end

    feats = [mean(mean_feats) ...
             mean(cons_feats) ];

    feats_mat(i,:) = feats;
    toc
end
toc
save(out_feat_name, 'feats_mat');
save(out_feat_frames_name, 'feats_mat_frames');

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function nedge=savegts(v,f,fname,edges)
%
% nedge=savegts(v,f,fname,edges)
%
% save a surface mesh to GNU Triangulated Surface Format (GTS)
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2011/02/23
%
% input:
%      v: input, surface node list, dimension (nn,3)
%      f: input, surface face element list, dimension (be,3)
%      fname: output file name
%      edges: edge list, if ignored, savegts will compute
%
% output:
%      nedge: the number of unique edges in the mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fid=fopen(fname,'wt');
if(fid==-1)
    error('You do not have permission to save mesh files.');
end
v=v(:,1:3);
f=f(:,1:3);
[edges,idx,edgemap]=uniqedges(f);
nedge=size(edges,1);

fprintf(fid,'%d %d %d\n',size(v,1),nedge,size(f,1));
fprintf(fid,'%.16f %.16f %.16f\n',v');
fprintf(fid,'%d %d\n',edges');
fprintf(fid,'%d %d %d\n',edgemap');
fclose(fid);


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rfAmp,rfPhase,rfFreq,rfCoil,GzAmp,GyAmp,GxAmp,ADC,Ext,uts,ts,flags]=PSD_GRE3DSpiral
global VCtl
global VVar
CV1=10e-3;
CV10=0;
CV11=0;
CV12=0;
CV13=0;
CV14=0;
CV2=20e-3;
CV3=0;
CV4=0;
CV5=0;
CV6=0;
CV7=0;
CV8=0;
CV9=0;
rfAmpAll=[];
rfPhaseAll=[];
rfFreqAll=[];
rfCoilAll=[];
GzAmpAll=[];
GyAmpAll=[];
GxAmpAll=[];
ADCAll=[];
ExtAll=[];
rfTimeAll=[];
GzTimeAll=[];
GyTimeAll=[];
GxTimeAll=[];
ADCTimeAll=[];
ExtTimeAll=[];
SEtAll=[];
uts=[];
ts=[];
flags=[];
if VCtl.PlotSeq == 1
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
ts = [ts tS tE];
end
ts = [0 max(ts)-min(ts)];
return;
end
%==============Pulses 1==============
rfAmp=[];
rfPhase=[];
rfFreq=[];
rfCoil=[];
GzAmp=[];
GyAmp=[];
GxAmp=[];
ADC=[];
Ext=[];
rfTime=[];
GzTime=[];
GyTime=[];
GxTime=[];
ADCTime=[];
ExtTime=[];
Freq=1;
Notes='regular TR section';
AttributeOpt={'on','off'};
Switch=AttributeOpt{1};
TREnd=Inf;
TRStart=1;
tE=VCtl.TR;
tS=0;
if isempty(tS) | isempty(tE) | (tS>=tE)
error('SE setting is incorrect for Pulses 1!');
end
if VVar.TRCount<TRStart | VVar.TRCount>TREnd | mod(VVar.TRCount-TRStart,Freq)~=0 | strcmp(Switch,'off')
% do nothing
else
%--------------------
AttributeOpt={'on','off'};
p.AnchorTE=AttributeOpt{1};
AttributeOpt={'Non','Hamming','Hanning'};
p.Apod=AttributeOpt{1};
p.CoilID=1;
p.DupSpacing=0;
p.Duplicates=1;
p.FA=VCtl.FlipAng;
p.Notes='excitation rf';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.TBP=4;
p.dt=20e-6;
p.rfFreq=0;
p.rfPhase=0;
p.tEnd=1e-3;
p.tStart=0;
if strcmp(p.Switch,'on')
if strcmp(p.AnchorTE,'on')
switch VCtl.TEAnchor
case 'Start'
VCtl.TEAnchorTime=p.tStart;
case 'Middle'
VCtl.TEAnchorTime=(p.tStart+p.tEnd)/2;
case 'End'
VCtl.TEAnchorTime=p.tEnd;
end
end
[rfAmp1,rfPhase1,rfFreq1,rfCoil1,rfTime1]=rfSinc(p);
if strcmp(VCtl.MultiTransmit,'off')
if VCtl.MasterTxCoil==rfCoil1(1)
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
else
rfAmp=[rfAmp rfAmp1];
rfPhase=[rfPhase rfPhase1];
rfFreq=[rfFreq rfFreq1];
rfCoil=[rfCoil rfCoil1];
rfTime=[rfTime rfTime1];
end
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Gz1Sign=1;
p.Gz2Sign=0;
p.Notes='cartesian phase';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.t1End=CV2;
p.t1Start=CV1;
p.t2End=VCtl.TE+CV2;
p.t2Start=VCtl.TE+CV1;
p.tRamp=100e-6;
if strcmp(p.Switch,'on')
[GzAmp1,GzTime1]=GzCartesian(p);
GzAmp=[GzAmp GzAmp1];
GzTime=[GzTime GzTime1];
end
p=[];
%--------------------
p.Notes='spiral encoding';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.dt=0.004e-3;
p.tStart=VCtl.TE;
if strcmp(p.Switch,'on')
[GyAmp1,GyTime1]=GySpiral(p);
GyAmp=[GyAmp GyAmp1];
GyTime=[GyTime GyTime1];
end
p=[];
%--------------------
p.Notes='spiral encoding';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.dt=0.004e-3;
p.tStart=VCtl.TE;
if strcmp(p.Switch,'on')
[GxAmp1,GxTime1]=GxSpiral(p);
GxAmp=[GxAmp GxAmp1];
GxTime=[GxTime GxTime1];
end
p=[];
%--------------------
p.Notes='spiral readout';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=VCtl.TE;
if strcmp(p.Switch,'on')
[ADC1,ADCTime1]=ADCSpiral(p);
ADC=[ADC ADC1];
ADCTime=[ADCTime ADCTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=1;
p.Notes='reset K space location';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=1e-3;
if strcmp(p.Switch,'on')
[Ext1,ExtTime1]=ExtBit(p);
Ext=[Ext Ext1];
ExtTime=[ExtTime ExtTime1];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=5;
p.Notes='calculate remaining scan time';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=0;
if strcmp(p.Switch,'on')
[Ext2,ExtTime2]=ExtBit(p);
Ext=[Ext Ext2];
ExtTime=[ExtTime ExtTime2];
end
p=[];
%--------------------
p.DupSpacing=0;
p.Duplicates=1;
p.Ext=9;
p.Notes='real time recon';
AttributeOpt={'on','off'};
p.Switch=AttributeOpt{1};
p.tStart=VCtl.TE+CV2;
if strcmp(p.Switch,'on')
[Ext3,ExtTime3]=ExtBit(p);
Ext=[Ext Ext3];
ExtTime=[ExtTime ExtTime3];
end
p=[];
%--------------------
SEt=[tS tE];
rfAmp(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfPhase(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfFreq(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
rfCoil(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzAmp(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyAmp(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxAmp(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADC(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
Ext(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfTime(rfTime<0 | rfTime>SEt(2)-SEt(1)) = [];
GzTime(GzTime<0 | GzTime>SEt(2)-SEt(1)) = [];
GyTime(GyTime<0 | GyTime>SEt(2)-SEt(1)) = [];
GxTime(GxTime<0 | GxTime>SEt(2)-SEt(1)) = [];
ADCTime(ADCTime<0 | ADCTime>SEt(2)-SEt(1)) = [];
ExtTime(ExtTime<0 | ExtTime>SEt(2)-SEt(1)) = [];
rfAmp(abs(rfAmp)<eps) = 0;
rfTime = rfTime + SEt(1);
GzTime = GzTime + SEt(1);
GyTime = GyTime + SEt(1);
GxTime = GxTime + SEt(1);
ADCTime = ADCTime + SEt(1);
ExtTime = ExtTime + SEt(1);
rfAmpAll=[rfAmpAll rfAmp];
rfPhaseAll=[rfPhaseAll rfPhase];
rfFreqAll=[rfFreqAll rfFreq];
rfCoilAll=[rfCoilAll rfCoil];
GzAmpAll=[GzAmpAll GzAmp];
GyAmpAll=[GyAmpAll GyAmp];
GxAmpAll=[GxAmpAll GxAmp];
ADCAll=[ADCAll ADC];
ExtAll=[ExtAll Ext];
rfTimeAll=[rfTimeAll rfTime];
GzTimeAll=[GzTimeAll GzTime];
GyTimeAll=[GyTimeAll GyTime];
GxTimeAll=[GxTimeAll GxTime];
ADCTimeAll=[ADCTimeAll ADCTime];
ExtTimeAll=[ExtTimeAll ExtTime];
SEtAll=[SEtAll SEt];
end
%====================================
if isempty(rfTimeAll)
error('rf sequence line can not be empty! Master Tx coil element must be used.');
end
if isempty(GzTimeAll)
error('GzSS sequence line can not be empty!');
end
if isempty(GyTimeAll)
error('GyPE sequence line can not be empty!');
end
if isempty(GxTimeAll)
error('GxR sequence line can not be empty!');
end
if isempty(ADCTimeAll)
error('ADC sequence line can not be empty!');
end
if isempty(ExtTimeAll)
error('Ext sequence line can not be empty!');
end
SEflag=repmat([0 0 0 0 0 0]',[1 2]);
rfflag=repmat([1 0 0 0 0 0]',[1 max(size(rfTimeAll))]);
Gzflag=repmat([0 1 0 0 0 0]',[1 max(size(GzTimeAll))]);
Gyflag=repmat([0 0 1 0 0 0]',[1 max(size(GyTimeAll))]);
Gxflag=repmat([0 0 0 1 0 0]',[1 max(size(GxTimeAll))]);
ADCflag=repmat([0 0 0 0 1 0]',[1 max(size(ADCTimeAll))]);
Extflag=repmat([0 0 0 0 0 1]',[1 max(size(ExtTimeAll))]);
ts=[[min(SEtAll) max(SEtAll)] rfTimeAll GzTimeAll GyTimeAll GxTimeAll ADCTimeAll ExtTimeAll]-min(SEtAll);
flags=[SEflag rfflag Gzflag Gyflag Gxflag ADCflag Extflag];
[ts,ind]=sort(ts);
uts=unique(ts);
flags=flags(:,ind);
[rfTime,ind]=sort(rfTimeAll-min(SEtAll));
rfAmp=rfAmpAll(:,ind);
rfPhase=rfPhaseAll(:,ind);
rfFreq=rfFreqAll(:,ind);
rfCoil=rfCoilAll(:,ind);
[GzTime,ind]=sort(GzTimeAll-min(SEtAll));
GzAmp=GzAmpAll(:,ind);
[GyTime,ind]=sort(GyTimeAll-min(SEtAll));
GyAmp=GyAmpAll(:,ind);
[GxTime,ind]=sort(GxTimeAll-min(SEtAll));
GxAmp=GxAmpAll(:,ind);
[ADCTime,ind]=sort(ADCTimeAll-min(SEtAll));
ADC=ADCAll(:,ind);
[ExtTime,ind]=sort(ExtTimeAll-min(SEtAll));
Ext=ExtAll(:,ind);
rfAmp(1) = 0;
rfPhase(1) = 0;
rfFreq(1) = 0;
GzAmp(1) = 0;
GyAmp(1) = 0;
GxAmp(1) = 0;
ADC(1) = 0;
Ext(1) = 0;
rfAmp(end) = 0;
rfPhase(end) = 0;
rfFreq(end) = 0;
GzAmp(end) = 0;
GyAmp(end) = 0;
GxAmp(end) = 0;
ADC(end) = 0;
Ext(end) = 0;
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%--- help for mcmc.reload_draws ---
%
%  subset is a 1 x 2 cell array in which
%  - the first cell contains the columns to retain: can be empty, defaults
%  to all
%  - the second column contains the chains to retain: can be empty, defaults
%  to all
%

author: jmaih


text:
function f = p05_f ( m, n, x )

%*****************************************************************************80
%
%% P05_F evaluates the objective function for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Marcin Molga, Czeslaw Smutnicki,
%    Test functions for optimization needs.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of arguments.
%
%    Input, real X(M,N), the arguments.
%
%    Output, real F(N), the function evaluated at the arguments.
%
  f = zeros ( n, 1 );

  for j = 1 : n

    f(j) = 10 * m;

    for i = 1 : m
      f(j) = f(j) + x(i,j)^2 - 10.0 * cos ( 2.0 * pi * x(i,j) );
    end

  end

  return
end


text:
%  INTERNAL FUNCTION: third-order multivariate chain rule
%
%  ::
%
%    res=third_order(dvvv,dvv,dv,vzzz,vzz,vz)
%    res=third_order(dvvv,dvv,dv,vzzz,vzz,vz,options)
%
%  Args:
%
%     - **dvvv** [nd x nv^3 matrix]: matrix of third derivatives of the d
%       function with respect to its locations. The derivatives are unfolded
%       columnwise
%     - **dvv** [nd x nv^2 matrix]: matrix of second derivatives of the d
%       function with respect to its locations. The derivatives are unfolded
%       columnwise
%     - **dv** [nd x nv matrix]: jacobian of function with respect to the
%       locations of its arguments
%     - **vzzz** [nv x nz^3 matrix]: third derivatives of the locations with
%       respect to the variables to differentiate. The derivatives are unfolded
%       columnwise
%     - **vzz** [nv x nz^2 matrix]: second derivatives (hessian) of the
%       locations with respect to the variables to differentiate. The derivatives
%       are unfolded columnwise
%     - **vz** [nv x nz matrix]: jacobian of the locations with respect to the
%       variables to differentiate
%     - **options** [empty|struct]: When not empty, options is a structure with
%       fields:
%
%       - **large** [true|{false}] if true, a computation explicitly using the
%         kronecker product is avoided.
%       - **multiply** [true|{false}]: if true, explicit omega matrices are
%         constructed and then multiplied to other matrices to sum the
%         permutations. Else, a functional form is used instead.
%
%  Returns:
%     :
%
%     - **res** [nd x nz^3]: output matrix
%
%


text:
% FEATURE NORMALIZE function.
% Normalizes the features in X. Returns a normalized version of X where the mean value of
% each feature is 0 and the standard deviation is 1.
function [X_normalized, mu, sigma] = feature_normalize(X)
    X_normalized = X;
    mu = zeros(1, size(X_normalized, 2));
    sigma = zeros(1, size(X_normalized, 2));

    % Get average values for each feature (column) in X.
    mu = mean(X_normalized);

    % Calculate the standard deviation for each feature.
    sigma = std(X_normalized);

    % Subtract mean values from each feature (column) of every example (row)
    % to make all features be spread around zero.
    X_normalized = X_normalized - mu;

    % Normalize each feature values for each example so that all features
    % are close to [-1:1] boundaries.
    X_normalized = X_normalized ./ sigma;
end

author: trekhleb

repo: machine-learning-octave


text:
function varargout=rotqr2qc(varargin)
%
% For calling details please see v_rotqr2qc.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_rotqr2qc.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: rotqr2qc.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_rotqr2qc(varargin{:});
else
    v_rotqr2qc(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function PlotAdaptiveContour3D(u, levels, tol)

% function PlotAdaptiveContour3D(u, levels, tol)
% Purpose: adaptively refine the mesh to approximately locate isocontours

Globals3D;

% build interpolation matrix (coarse->fine)
EToVi = [1 5 7 8; 5 2 6 9; 7 6 3 10; 8 9 10 4; 8 5 7 9; 7 5 6 9; 8 9 7 10; 9 6 7 10];
VXi   = [-1  1 -1 -1  0  0 -1 -1  0 -1];
VYi   = [-1 -1  1 -1 -1  0  0 -1 -1  0];
VZi   = [-1 -1 -1  1 -1 -1 -1  0  0  0];

v1 = EToVi(:,1); v2 = EToVi(:,2); v3 = EToVi(:,3); v4 = EToVi(:,4);
ri = 0.5*(-(r+s+t+1)*VXi(v1) + (1+r)*VXi(v2) + (1+s)*VXi(v3) + (1+t)*VXi(v4) );
si = 0.5*(-(r+s+t+1)*VYi(v1) + (1+r)*VYi(v2) + (1+s)*VYi(v3) + (1+t)*VYi(v4) );
ti = 0.5*(-(r+s+t+1)*VZi(v1) + (1+r)*VZi(v2) + (1+s)*VZi(v3) + (1+t)*VZi(v4) );

interp = Vandermonde3D(N, ri(:), si(:), ti(:))*invV;
ri = [-1;1;-1;-1]; si = [-1;-1;1;-1]; ti = [-1;-1;-1;1]; refNp = length(ri);
interp1 = Vandermonde3D(N, ri(:), si(:), ti(:))*invV;

Nlevels = length(levels);

sk = 1;
F = spalloc(Np,Np,1);
for i=0:N % old ordering
  for j=0:N - i
    for k=0:N - i - j
      if(i+j+k<=1), F(sk,sk) = 1.; end;
      sk = sk+1;
    end
  end
end

hold on;
for n=1:Nlevels

  xref = x; yref = y; zref = z; uref = u; Kref = K; Jref = J;

  err = 1;
  while(err > tol)

    level = levels(n);

    umin = min(uref, [], 1);
    umax = max(uref, [], 1);

    refineflag = (umin <= level & umax >= level);

    toref = find( refineflag);
    Nref = length(toref);

    uref = reshape(interp*uref(:,toref), Np, 8*Nref);
    xref = reshape(interp*xref(:,toref), Np, 8*Nref);
    yref = reshape(interp*yref(:,toref), Np, 8*Nref);
    zref = reshape(interp*zref(:,toref), Np, 8*Nref);
    Jref = reshape(interp*Jref(:,toref), Np, 8*Nref)/8;

    Kref = 8*Nref;

    ufilt = V*(F*(invV*uref));
    err = max(max(abs(ufilt-uref)));

  end

  xref = interp1*xref; yref = interp1*yref; zref = interp1*zref; uref = interp1*uref;

  tets = reshape(1:4*Kref, 4, Kref)';

  x1 = xref(1,:)'; y1 = yref(1,:)'; z1 = zref(1,:)'; u1 = uref(1,:)' + 1e-10*rand(Kref,1);
  x2 = xref(2,:)'; y2 = yref(2,:)'; z2 = zref(2,:)'; u2 = uref(2,:)' + 1e-10*rand(Kref,1);
  x3 = xref(3,:)'; y3 = yref(3,:)'; z3 = zref(3,:)'; u3 = uref(3,:)' + 1e-10*rand(Kref,1);
  x4 = xref(4,:)'; y4 = yref(4,:)'; z4 = zref(4,:)'; u4 = uref(4,:)' + 1e-10*rand(Kref,1);

  trix = []; triy = []; triz = []; triu = [];;

  umin = min(uref, [], 1);
  umax = max(uref, [], 1);

  lev = levels(n);

  % find candidate elements
  ks = find(umin<= lev & umax >= lev);

  Nks = length(ks);

  if(Nks>0)

    fc = zeros(Nks,6);

    u1ks = u1(ks); u2ks = u2(ks); u3ks = u3(ks); u4ks = u4(ks);

    % find local coordinate at each of 6 edges
    c1 = (lev-u1ks)./(u2ks-u1ks);  fc(:,1) = (c1>=0 & c1<=1);
    c2 = (lev-u2ks)./(u3ks-u2ks);  fc(:,2) = (c2>=0 & c2<=1);
    c3 = (lev-u3ks)./(u1ks-u3ks);  fc(:,3) = (c3>=0 & c3<=1);
    c4 = (lev-u1ks)./(u4ks-u1ks);  fc(:,4) = (c4>=0 & c4<=1);
    c5 = (lev-u2ks)./(u4ks-u2ks);  fc(:,5) = (c5>=0 & c5<=1);
    c6 = (lev-u3ks)./(u4ks-u3ks);  fc(:,6) = (c6>=0 & c6<=1);

    % find triangle intersection
    tris = find(sum(fc, 2)==3); Ntris = length(tris);  tfc = fc(tris,:); ids = find(tfc');
    ktris = ks(tris);

    % trim list
    tc1 = c1(tris); tc2 = c2(tris); tc3 = c3(tris); tc4 = c4(tris); tc5 = c5(tris); tc6 = c6(tris);

    tx1 = x1(ktris);   ty1 = y1(ktris);   tz1 = z1(ktris);
    tx2 = x2(ktris);   ty2 = y2(ktris);   tz2 = z2(ktris);
    tx3 = x3(ktris);   ty3 = y3(ktris);   tz3 = z3(ktris);
    tx4 = x4(ktris);   ty4 = y4(ktris);   tz4 = z4(ktris);

    xc = zeros(Ntris, 6); yc = zeros(Ntris, 6); zc = zeros(Ntris, 6);
    xc(:,1) = (1-tc1).*tx1 + tc1.*tx2; yc(:,1) = (1-tc1).*ty1 + tc1.*ty2; zc(:,1) = (1-tc1).*tz1 + tc1.*tz2;
    xc(:,2) = (1-tc2).*tx2 + tc2.*tx3; yc(:,2) = (1-tc2).*ty2 + tc2.*ty3; zc(:,2) = (1-tc2).*tz2 + tc2.*tz3;
    xc(:,3) = (1-tc3).*tx3 + tc3.*tx1; yc(:,3) = (1-tc3).*ty3 + tc3.*ty1; zc(:,3) = (1-tc3).*tz3 + tc3.*tz1;
    xc(:,4) = (1-tc4).*tx1 + tc4.*tx4; yc(:,4) = (1-tc4).*ty1 + tc4.*ty4; zc(:,4) = (1-tc4).*tz1 + tc4.*tz4;
    xc(:,5) = (1-tc5).*tx2 + tc5.*tx4; yc(:,5) = (1-tc5).*ty2 + tc5.*ty4; zc(:,5) = (1-tc5).*tz2 + tc5.*tz4;
    xc(:,6) = (1-tc6).*tx3 + tc6.*tx4; yc(:,6) = (1-tc6).*ty3 + tc6.*ty4; zc(:,6) = (1-tc6).*tz3 + tc6.*tz4;
    xc = xc'; yc = yc'; zc = zc';

    ids = reshape(ids, 3, Ntris);

    trix = [trix, xc(ids)];
    triy = [triy, yc(ids)];
    triz = [triz, zc(ids)];
    triu = [triu, lev*ones(size(ids))];

    % find quadrilateral intersection
    quads = find(sum(fc, 2)==4); Nquads = length(quads);  qfc = fc(quads,:); ids = find(qfc');
    kquads = ks(quads);

    % quad list
    qc1 = c1(quads); qc2 = c2(quads); qc3 = c3(quads); qc4 = c4(quads); qc5 = c5(quads); qc6 = c6(quads);

    qx1 = x1(kquads);   qy1 = y1(kquads);   qz1 = z1(kquads);
    qx2 = x2(kquads);   qy2 = y2(kquads);   qz2 = z2(kquads);
    qx3 = x3(kquads);   qy3 = y3(kquads);   qz3 = z3(kquads);
    qx4 = x4(kquads);   qy4 = y4(kquads);   qz4 = z4(kquads);

    xc = zeros(Nquads, 6); yc = zeros(Nquads, 6); zc = zeros(Nquads, 6);
    xc(:,1) = (1-qc1).*qx1 + qc1.*qx2; yc(:,1) = (1-qc1).*qy1 + qc1.*qy2; zc(:,1) = (1-qc1).*qz1 + qc1.*qz2;
    xc(:,2) = (1-qc2).*qx2 + qc2.*qx3; yc(:,2) = (1-qc2).*qy2 + qc2.*qy3; zc(:,2) = (1-qc2).*qz2 + qc2.*qz3;
    xc(:,3) = (1-qc3).*qx3 + qc3.*qx1; yc(:,3) = (1-qc3).*qy3 + qc3.*qy1; zc(:,3) = (1-qc3).*qz3 + qc3.*qz1;
    xc(:,4) = (1-qc4).*qx1 + qc4.*qx4; yc(:,4) = (1-qc4).*qy1 + qc4.*qy4; zc(:,4) = (1-qc4).*qz1 + qc4.*qz4;
    xc(:,5) = (1-qc5).*qx2 + qc5.*qx4; yc(:,5) = (1-qc5).*qy2 + qc5.*qy4; zc(:,5) = (1-qc5).*qz2 + qc5.*qz4;
    xc(:,6) = (1-qc6).*qx3 + qc6.*qx4; yc(:,6) = (1-qc6).*qy3 + qc6.*qy4; zc(:,6) = (1-qc6).*qz3 + qc6.*qz4;
    xc = xc'; yc = yc'; zc = zc';

    ids1 = reshape(ids, 4, Nquads);
    ids2 = ids1([1;4;2;3],:);
    ids3 = ids1([1;3;4;2],:);

    % find length perimeters bounded by 4 nodes
    perims1 = QuadPerimeter3D(xc(ids1),yc(ids1),zc(ids1));
    perims2 = QuadPerimeter3D(xc(ids2),yc(ids2),zc(ids2));
    perims3 = QuadPerimeter3D(xc(ids3),yc(ids3),zc(ids3));

    % locate minimum perimeter quad
    idsA = find(perims1<min(perims2, perims3));
    idsB = find(perims2<min(perims1, perims3));
    idsC = find(perims3<min(perims1, perims2));

    % use only minimum perimeter quads
    ids  = [ids1(:,idsA),ids2(:,idsB),ids3(:,idsC)];

    % divide into two triangles (ok since quad must be convex)
    ids1 = [ids([1;2;3],:),ids([1;3;4],:)];

    % add triangles to list
    trix = [trix, xc(ids1)];
    triy = [triy, yc(ids1)];
    triz = [triz, zc(ids1)];
    triu = [triu, lev*ones(size(ids1))];

  end

  size(trix)

  ha = patch(trix,triy,triz,triu); shading interp;

  set(ha, 'EdgeColor', 'none')
end
hold off;

camlight; material shiny; lighting gouraud; view(3); axis equal;
axis off;

return;

repo: Aero-matlab