text:
function [rvts,rvtsSS,rvtsD,rvtsSSD,diagnostics] = realized_twoscale_variance(price,time,timeType,samplingType,samplingInterval,subsamples,options)
% Estimated quadratic variation using the Two-Scale estimator of Ait-Sahalia, Mykland and Zhang
%
% USAGE:
%   RVTS = realized_twoscale_variance(PRICE)
%   [RVTS,RVTSSS,RVTSD,RVTSSSD,DIAGNOSTICS] = realized_twoscale_variance(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,SUBSAMPLES,OPTIONS)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - [OPTIONAL] m by 1 vector of times where TIME(i) corresponds to PRICE(i).
%   TIMETYPE         - [OPTIONAL] String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS.mmm, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - [OPTIONAL] String describing the type of sampling to use when
%                        filtering PRICE for estimating the fast scale estimator.
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVAL ticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL  - [OPTIONAL] Scalar integer or n by 1 vector whose meaning depends on the
%                         selected SAMPLINGTYPE
%   OPTIONS           - [OPTIONAL] Realized Two-Scale option structure initialized by calling
%                         realized_options('Multiscale'). See help realized_options for a description of
%                         available options.
%
% OUTPUTS:
%   RVTS        - Realized two-scale variance estimate
%   RVTSSS      - Realized two-scale variance estimate constructed by averaging RVTS across multiple
%                 initial observations.  If SAMPLINGTYPE is 'BusinessTime' and SAMPLINGINTERAL is 1
%                 then subsampling is not possible.
%   RVTSD       - Debiased version of realized two-scale variance estimate
%   RVTSSSD     - Debiased version of realized two-scale variance estimate constructed by averaging RVTS across multiple
%                 initial observations.
%   DIAGNOSTICS - Structure of useful diagnostic information.  Fields are
%                   BANDWIDTH             - Ratio of slow scale to fast scale where the fast scale
%                                           sampling is determined by SAMPLINGTYPE and SAMPLINGINTERVAL.
%                   NOISEVARIANCE         - Bandi-Russell noise variance estimate. Empty if user
%                                           supplies bandwidth.
%                   DEBIASEDNOISEVARIANCE - Bias adjusted estimate of the noise variance.
%                   IQESTIMATE            - Estimate of IQ used. Empty if not needed. Empty if user
%                                           supplies bandwidth.
%
% COMMENTS:
%   For best results:
%     - Use prices sampled close to the highest frequency available, if not using 'BusinessTime' and
%     - If using calendar time sampling, sample close to the limit of the data availability (e.g. 5-15
%       seconds for a stock that trades a few thousand times per day)
%     - When sampling less frequently (e.g. 5 - 30 minutes) standard realized variance is
%       probably more appropriate
%     - The debiased version should be preferred to the standard version
%
% EXAMPLE:
%  % Default usage with 'BusinessTime' sampling with interval 1 and automatic bandwidth selection
%  [RVTS,RVTSSS,RVTSD,RVTSSSD] = realized_twoscale_variance(PRICE)
%
%  % 5-tick two-scale variance with automatic bandwidth selection
%  [RVTS,RVTSSS,RVTSD,RVTSSSD] = realized_twoscale_variance(PRICE,TIME,'wall','BusinessTime',5)
%
%  % Two scale variance with a bandwidth of 30 sampling at every tick
%  options = realized_options('Twoscale');
%  options.bandwidth = 30;
%  [RVTS,RVTSSS,RVTSD,RVTSSSD] = realized_twoscale_variance(PRICE,TIME,'wall','BusinessTime',1,0,options)
%
%  % 5-tick two-scale variance with automatic bandwidth selection and dense subsampling
%  [RVTS,RVTSSS,RVTSD,RVTSSSD] = realized_twoscale_variance(PRICE,TIME,'wall','BusinessTime',5,4)
%
%  See also REALIZED_OPTIONS, REALIZED_NOISE_ESTIMATE, REALIZED_VARIANCE,
%  REALIZED_VARIANCE_OPTIMAL_SAMPLING, REALIZED_RANGE, REALIZED_QUANTILE_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 1
        m = length(price);
        time = linspace(9.5*3600,16*3600,m)';
        timeType = 'seconds';
        samplingType = 'businesstime';
        samplingInterval = 1;
        subsamples = 1;
        options = realized_options('twoscale');
    case 5
        subsamples = 1;
        options = realized_options('twoscale');
    case 6
        options = realized_options('twoscale');
    case 7
        % Nothing
    otherwise
        error('One, five, six or seven inputs required.')
end
samplingType=lower(samplingType);
timeType=lower(timeType);
[errorMessage, options] = realized_twoscale_variance_parameter_check(price,time,timeType,samplingType,samplingInterval,subsamples,options);

if ~isempty(errorMessage )
    error(errorMessage)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

logPrice = log(price);
filteredLogPrice = realized_price_filter(logPrice,time,timeType,samplingType,samplingInterval);
%options.bandwidth= 60;
if isempty(options.bandwidth)
    % Estimate the optimal bandwidth if not provided
    [noiseVariance, debiasedNoiseVariance, IQEstimate] = realized_noise_estimate(price, time, timeType, options);
    diagnostics.noiseVariance=noiseVariance;
    diagnostics.debiasedNoiseVariance=debiasedNoiseVariance;
    diagnostics.IQEstimate=IQEstimate;
    % Select the correct noise variance estimate
    if options.useDebiasedNoise
        selectedNoiseVariance = debiasedNoiseVariance;
    else
        selectedNoiseVariance = noiseVariance;
    end
    % Compute the bandwidth
    n = length(filteredLogPrice) - 1;
    cOpt = ((12 * selectedNoiseVariance^2) / IQEstimate)^(1/3);
    options.bandwidth = ceil(cOpt * n^(2/3));
end

bandwidth = ceil(options.bandwidth);
if bandwidth<2
    warning('oxfordRealized:smallBandwidth','The selected bandwidth is less then 2, and RTVS is not well defined in this case.  Setting bandwidth to 2 and proceeding.')
    bandwidth = 2;
end
diagnostics.bandwidth = bandwidth;

n = length(filteredLogPrice) - 1;
K = bandwidth;
nbar = (n - K + 1) / K;
[RV,count,naturalCount] = overlap_realized_variance(filteredLogPrice,K);
overlapScale = count / naturalCount;
RVq = RV / overlapScale;
RV1 = overlap_realized_variance(filteredLogPrice,1);
% Eq. 55
rvts = RVq - nbar/n * RV1;
% Eq. 64
rvtsD = (1-nbar/n)^(-1)*rvts;

% Subsampled RVq and RV1
subsampledLogPrices = realized_subsample(logPrice,time,timeType,samplingType,samplingInterval,subsamples);
rvqs = zeros(subsamples,1);
rv1s = zeros(subsamples,1);
totalCount  = 0;
totalCount0 = 0;
for i=1:subsamples
    filteredLogPrice =  subsampledLogPrices{i};
    [RV,count] = overlap_realized_variance(filteredLogPrice,K);
    rvqs(i) = RV;
    totalCount = totalCount + count;
    if i==1
        baseCount = count;
    end
    [rv1s(i),count] = overlap_realized_variance(filteredLogPrice,1);
    totalCount0 = totalCount0 + count;
    if i==1
        baseCount0 = count;
    end
end
RVqSS = sum(rvqs) * (baseCount / totalCount);
% RVq is now based on all overlapping blcoks, so need to normalize to the
% non-overlapping case
RVqSS = RVqSS / overlapScale;
RV1SS = sum(rv1s) * (baseCount0 / totalCount0);
rvtsSS = RVqSS - nbar/n * RV1SS;
rvtsSSD = (1-nbar/n)^(-1)*rvtsSS;
end

function [rv,count,naturalCount] = overlap_realized_variance(price,skip)

naturalCount = (length(price)-1)/skip;
m = length(price);
returns = price(1+skip:m) - price(1:m-skip);
count = length(returns);
rv = returns' * returns;

end

function [errorMessage,options] = realized_twoscale_variance_parameter_check(price,time,timeType,samplingType,samplingInterval,subsamples,options)
% Support function for realized_twoscale_variance that does input validation
%
% USAGE:
%   [ERRORMESSAGE,OPTIONS] = realized_variance_optimal_sampling_parameter_check(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINVERVAL,OPTIONS)
%
% INPUTS:
%   See realized_variance_optimal_sampling
%
% OUTPUT:
%   ERRORMESSAGE - String containing a description of the error if one is detected.  Empty if no error.
%
% COMMENTS:
%   See realized_options for a description of the other fields in OPTIONS
%
%  See also REALIZED_VARIANCE_OPTIMAL_SAMPLING

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008
errorMessage = [];

if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    errorMessage = 'PRICE must be a m by 1 vector.';
    return
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    errorMessage = 'TIME must be sorted and increasing';
    return
end
if size(time,2)>1 || length(time)~=length(price)
    errorMessage = 'TIME must be a m by 1 vector.';
    return
end

timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    errorMessage = 'TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.';
    return;
end

samplingType=lower(samplingType);
if ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform'})
    % Must be a scalar integer
    if ~isscalar(samplingInterval) || floor(samplingInterval)~=samplingInterval || samplingInterval<1
        errorMessage = 'SAMPLINGINTERVAL must be a positive integer for the SAMPLINGTYPE selected.';
        return
    end
else
    if size(samplingInterval,2)>size(samplingInterval,1)
        samplingInterval=samplingInterval';
    end
    if ~(any(samplingInterval>=t0) && any(samplingInterval<=tT))
        errorMessage = 'At least one sampling interval must be between min(TIME) and max(TIME) when using ''Fixed'' as SAMPLINGTYPE.';
        return
    end
    if any(diff(samplingInterval)<=0)
        errorMessage = 'When using ''Fixed'' as SAMPLINGTYPE the vector of sampling times in SAMPLINGINTERVAL must be sorted and strictly increasing.';
        return
    end
end

% Options

% List of flat top kernels
flatTopKernelList = {'bartlett','twoscale','2ndorder','epanechnikov',...
    'cubic','multiscale','5thorder','6thorder','7thorder','8thorder','parzen',...
    'th1','th2','th5','th16'};

% List of non flat top kernels
nonFlatTopKernelList = {'nonflatparzen','qs','fejer','thinf','bnhls'};

% Combined kernel list
kernelList = [flatTopKernelList nonFlatTopKernelList];

% Check fields for valid values
optionsFieldNames = fieldnames(options);
% Insert any missing fiedls
defaultOptions = realized_options('Optimal Sampling');
defaultFieldNames = fieldnames(defaultOptions);
missingFields = setdiff(defaultFieldNames,optionsFieldNames);
if ~isempty(missingFields)
    for i = 1:length(missingFields)
        options.(missingFields{i}) = defaultOptions.(missingFields{i});
    end
end

for i=1:length(optionsFieldNames)
    fieldName = optionsFieldNames{i};
    fieldValue = options.(fieldName);
    if ischar(fieldValue)
        fieldValue = lower(fieldValue);
        options.(fieldName) = fieldValue;
    end

    switch fieldName
        case {'medFrequencyKernel'}
            % Member of kernelList
            if ~ismember(fieldValue,kernelList)
                errorMessage = ['OPTIONS.' fieldName ' must be one of the listed types.'];
                return
            end
        case {'medFrequencyBandwidth'}
            % Non-negative scalar
            if ~isempty(fieldValue) && ~isnonnegativescalar(fieldValue)
                errorMessage = ['OPTIONS.' fieldName ' must a non-negative scalar.'];
                return
            end
        case {'useDebiasedNoise','useAdjustedNoiseCount'}
            % Logical or scalar
            if ~islogical(fieldValue) && ~ismember(fieldValue,[0 1])
                errorMessage = 'OPTIONS.useDebiasedNoise must be a logical value.';
                return
            end
        case {'IQEstimationSamplingType','medFrequencySamplingType','noiseVarianceSamplingType'}
            % One of the sampling types
            if ~ismember(fieldValue,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
                errorMessage = ['OPTIONS.' fieldName ' must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.'];
                return
            end
        case {'medFrequencySamplingInterval','noiseVarianceSamplingInterval','IQEstimationSamplingInterval'}
            % Non-negative scalar, less that 1 if timeType is unit
            if isempty(fieldValue) || ~isnonnegativescalar(fieldValue)
                errorMessage = ['OPTIONS.' fieldName ' must be a non-negative scalar between 0 and 1.'];
                return
            end
            if strcmp(timeType,'unit') && fieldValue>1
                errorMessage = ['OPTIONS.' fieldName ' must be less than 1 if TIMETYPE when ''unit''.'];
                return
            end
    end
end

if ~isempty(subsamples)
    if ~isscalar(subsamples) || subsamples<0 || floor(subsamples)~=subsamples
        errorMessage = 'SUBSAMPLES must be a non-negative scalar.';
    end
end
end

function condition = isnonnegativescalar(x)
% Function that returns logical true if that input is a non-empty scalar >=0
condition = ~isempty(x) && isscalar(x) && x>=0;
end

author: bashtage

repo: mfe-toolbox


text:
function out=mrdivide(x,y)

precisionx=0;   precisiony=0;
if isa(x,'mp')
 precisionx=x(1).precision;
end
if isa(y,'mp')
 precisiony=y(1).precision;
end
precision=max(precisionx,precisiony);
ex=numel(x);
ey=numel(y);

if ey==1
 out=x./y;
else
 error('full mrdivide hasn''t been written for mp objects yet')
end

author: opencobra

repo: cobratoolbox


text:
%%%  ModeSwitch_HIL
clear
path( './icon/',path);
load vehicle_local_position.mat;
%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
%Attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_ANGLE = 3;
Kp_YAW_AngleRate = 0.3;
Ki_YAW_AngleRate = 0.001;
Kd_YAW_AngleRate = 0.00;
%Position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
MAX_CONTROL_ANGLE_NAV_ROLL = 15;
MAX_CONTROL_ANGLE_NAV_PITCH  = 15;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
ModeSwitch_HIL

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
%==============================================================================
% Copyright (C) 2005, Jan Modersitzki and Nils Papenberg, see copyright.m,
% this file is part of the FLIRT Package, all rights reserved;
% http://www.math.uni-luebeck.de/SAFIR/FLIRT-MATLAB.html
%==============================================================================
% function [B,Bstr] = getDiffusiveMatrix(Omega,m);
% generates the elastic matrix B for the domain Omega with resolution m,
% by default, mu = 1, lambda = 0;
%==============================================================================
function [B,Bstr] = getDiffusiveMatrixStg(Omega,m)

Bstr = 'diffusive-stg';
h = Omega./m;

d11 = spdiags(ones(m(1),1)*[-1,1],[0,1],m(1),m(1)+1)/h(1);
d12 = spdiags(ones(m(2)-1,1)*[-1,1],[0,1],m(2)-1,m(2))/h(2);
d21 = spdiags(ones(m(1)-1,1)*[-1,1],[0,1],m(1)-1,m(1))/h(1);
d22 = spdiags(ones(m(2),1)*[-1,1],[0,1],m(2),m(2)+1)/h(2);

D11 = sparse(kron(speye(m(2)),d11));
D12 = sparse(kron(d12,speye(m(1)+1)));
D21 = sparse(kron(speye(m(2)+1),d21));
D22 = sparse(kron(d22,speye(m(1))));

% build the diffusive operator
%
%      | \nabla 0      |
%  B = |               |
%      | 0      \nabla |
%      |               |
% B[U] = [\partial_1 U1,\partial_2 U1, \partial_1 U2,\partial_2 U2]

n1 = (m(1)+1)*m(2);
n2 = m(1)*(m(2)+1);
j1 = size(D11,1);
j2 = size(D12,1);
j3 = size(D21,1);
j4 = size(D22,1);

B = [  D11,sparse(j1,n2);
       D12,sparse(j2,n2);
       sparse(j3,n1),D21;
       sparse(j4,n1),D22];
return;
%==============================================================================

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [options] = getOptions(options,inputArgs,varargin)
	% Gets default options for a function and replaces them with inputArgs inputs if they are present in Name-Value pair input (e.g. varargin).
	% Biafra Ahanonu
	% Started: 2013.11.04.
	%
	% inputs
		% options - structure passed by parent function with each fieldname containing an option to be used by the parent function.
		% inputArgs - an even numbered cell array, with {'option','value'} as the ordering. Normally pass varargin.
	% Outputs
		% options - options structure passed back to parent function with modified Name-Value inputs to function added.
	% NOTE
		% Use the 'options' name-value pair to input an options structure that will overwrite default options in a function, example below.
		% options.Stargazer = 1;
		% options.SHH = 0;
		% getMutations(mutationList,'options',options);
		%
		% This is in contrast to using name-value pairs, both will produce the same result.
		% getMutations(mutationList,'Stargazer',1,'SHH',0);
		%
		% The 'passArgs' name-value pair will pass through the parent functions varargin to child functions.
	% USAGE
		% function [input1,input2] = exampleFxn(input1,input2,varargin)
		% 	%========================
		% 	% DESCRIPTION
		% 	options.Stargazer = '';
		% 	% DESCRIPTION
		% 	options.SHH = '';
		% 	% DESCRIPTION
		% 	options.Option3 = '';
		% 	% get options
		% 	options = getOptions(options,varargin); % ***HERE IS WHERE getOptions IS USED***
		% 	% display(options)
		% 	% unpack options into current workspace
		% 	% fn=fieldnames(options);
		% 	% for i=1:length(fn)
		% 	% 	eval([fn{i} '=options.' fn{i} ';']);
		% 	% end
		% 	%========================
		% 	try
		% 		% Do something.
		% 		% How to use the passArgs feature.
		% 		childFunction(arg1,arg2,'passArgs',varargin);
		% 	catch err
		% 		disp(repmat('@',1,7))
		% 		disp(getReport(err,'extended','hyperlinks','on'));
		% 		disp(repmat('@',1,7))
		% 	end
		% end

	[options] = ciapkg.io.getOptions(options,inputArgs,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [ compdata ] = fiff_read_ctf_comp(fid,node,chs)

%
% [ compdata ] = fiff_read_ctf_comp(fid,node,chs)
%
% Read the CTF software compensation data from the given node
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.8  2007/11/13 10:55:32  msh
%   Specify the second argument to all calls to the exist function
%
%   Revision 1.7  2006/09/25 19:48:16  msh
%   Added projection item kinds to fiff_define_constants
%   Changed some fields to int32 in surface structures
%
%   Revision 1.6  2006/08/09 15:22:51  msh
%   Removed debug printout.
%
%   Revision 1.5  2006/06/22 21:22:46  msh
%   Take into account the possibility of calibrated compensation matrices
%
%   Revision 1.4  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.3  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.2  2006/04/12 10:29:02  msh
%   Made evoked data writing compatible with the structures returned in reading.
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_read_ctf_comp';

if nargin ~= 3
    error(me,'Incorrect number of arguments');
end

compdata = struct('ctfkind',{},'kind',{},'save_calibrated',{},'rowcals',{},'colcals',{},'data',{});
comps = fiff_dir_tree_find(node,FIFF.FIFFB_MNE_CTF_COMP_DATA);

for k = 1:length(comps)
    node = comps(k);
    %
    %   Read the data we need
    %
    mat  = fiff_read_named_matrix(fid,node,FIFF.FIFF_MNE_CTF_COMP_DATA);
    for p = 1:node.nent
        kind = node.dir(p).kind;
        pos  = node.dir(p).pos;
        if kind == FIFF.FIFF_MNE_CTF_COMP_KIND
            tag = fiff_read_tag(fid,pos);
            break;
        end
    end
    if ~exist('tag','var')
        error(me,'Compensation type not found');
    end
    %
    %   Get the compensation kind and map it to a simple number
    %
    one.ctfkind = tag.data;
    clear('tag');
    one.kind    = int32(-1);
    if one.ctfkind     == hex2dec('47314252')
        one.kind = int32(1);
    elseif one.ctfkind == hex2dec('47324252')
        one.kind = int32(2);
    elseif one.ctfkind == hex2dec('47334252')
        one.kind = int32(3);
    else
        one.kind = one.ctfkind;
    end
    for p = 1:node.nent
        kind = node.dir(p).kind;
        pos  = node.dir(p).pos;
        if kind == FIFF.FIFF_MNE_CTF_COMP_CALIBRATED
            tag = fiff_read_tag(fid,pos);
            break;
        end
    end
    if ~exist('tag','var')
        calibrated = false;
    else
        calibrated = tag.data;
    end
    one.save_calibrated = calibrated;
    one.rowcals = ones(1,size(mat.data,1));
    one.colcals = ones(1,size(mat.data,2));
    if ~calibrated
        %
        %   Calibrate...
        %
        %
        %   Do the columns first
        %
        for p  = 1:length(chs)
            ch_names{p} = chs(p).ch_name;
        end
        for col = 1:size(mat.data,2)
            p = strmatch(mat.col_names{col},ch_names,'exact');
            if isempty(p)
                error(me,'Channel %s is not available in data',mat.col_names{col});
            elseif length(p) > 1
                error(me,'Ambiguous channel %s',mat.col_names{col});
            end
            col_cals(col) = 1.0/(chs(p).range*chs(p).cal);
        end
        %
        %    Then the rows
        %
        for row = 1:size(mat.data,1)
            p = strmatch(mat.row_names{row},ch_names,'exact');
            if isempty(p)
                error(me,'Channel %s is not available in data',mat.row_names{row});
            elseif length(p) > 1
                error(me,'Ambiguous channel %s',mat.row_names{row});
            end
            row_cals(row) = chs(p).range*chs(p).cal;
        end
        mat.data            = diag(row_cals)*mat.data*diag(col_cals);
        one.rowcals         = row_cals;
        one.colcals         = col_cals;
    end
    one.data       = mat;
    compdata(k)    = one;
    clear('row_cals');
    clear('col_cals');
end

if length(compdata) > 0
    fprintf(1,'\tRead %d compensation matrices\n',length(compdata));
end

return;

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef DisjointSet < handle
%%DISJOINTSET A disjoint set class that can be used for partitioning a
%             general set. It can be used to cluster targets given
%             measurements. If one also wishes to keep track of which
%             measurements are associated with which target clusters, then
%             the subclass DisjointSetM should be used.
%
%Disjoint sets are described in Chapter 8 of [1], in Chapter 21 of [2], and
%Chapter 12.4 of [3].
%
%The main point of this class is to bring together cluster data so that a
%ClusterSet class can be created. A ClusterSet class allows the clusters to
%be indexed by a cluster number and the elements of the clusters to be
%subsequently indexed with another number. The algorithms here use path
%compression in the find algorithm and a union by rank algorithm to process
%things quickly.
%
%REFERENCES:
%[1] M. A. Weiss, Data Structures and Algorithm Analysis in C++, 2nd ed.
%    Reading, MA: Addison-Wesley, 1999.
%[2] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein,
%    Introduction to Algorithms, 2nd ed. Cambridge, MA: The MIT Press,
%    2001.
%[3] B. R. Preiss, Data Structures and Algorithms with Object-Oriented
%    Design Patterns in C++. New York, NY: John Wiley & Sons, Inc., 1999.
%
%November 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    properties
        numberOfSets
        setArray
        numInTree
    end
    methods
        function newSet=DisjointSet(setSize)
            %The constructor method
            newSet.setArray=-1*ones(setSize,1);
            newSet.numInTree=ones(setSize,1);
            newSet.numberOfSets=setSize;
        end

        function unionFromList(DSObj,uList)
        %%UNIONFROMLIST Given a list of targets in the same cluster, such
        %               as a list of targets that gate with a common
        %               measurement, cluster them together.
        %
        %INPUTS: DSObj The implicitly passed calling object.
        %        uList A one-dimensional list of indices of targets to
        %              merge into one cluster

            numInList=length(uList);
            if(numInList==0)
                return;
            end

            rootIdx1=DSObj.find(uList(1));
            for curTar=2:numInList
                rootIdx2=DSObj.find(uList(curTar));
                rootIdx1=DSObj.unionRoots(rootIdx1,rootIdx2);
            end
        end

        function unionFromBinMat(DSObj,binMat)
        %%UNIONFROMBINMAT Given a binary matrix specifying which
        %                 measurements gate with which targets, cluster
        %                 together all targets that gate with each other.
        %
        %INPUTS: DSObj The implicitly passed calling object.
        %        binMat A numTarXnumMeas binary matrix where numTar is
        %               the same as the length of setArray in DSObj and
        %               numMeas is an arbitrary number. If an empty matrix
        %               is passed, nothing is done.
        %
        %This algorithm takes the union of all elements of setArray in each
        %column corresponding to ones in the column. Thus, if this is a
        %gating matrix for targets and measurements, this will create all
        %of the target clusters.
        %
            if(isempty(binMat))
                return;
            end

            numTar=size(binMat,1);
            numMeas=size(binMat,2);

            %Check for input validity
            assert(numTar==length(DSObj.setArray))

            for curMeas=1:numMeas
                rootIdx1=0;
                for curTar=1:numTar
                    if(binMat(curTar,curMeas)~=0)
                        if(rootIdx1==0)
                            rootIdx1=DSObj.find(curTar);
                        else
                            rootIdx2=DSObj.find(curTar);
                            rootIdx1=DSObj.unionRoots(rootIdx1,rootIdx2);
                        end
                    end
                end
            end
        end

        function newCSet=createClusterSet(DSObj)
        %%CREATECLUSTERSET Turn the DisjointSet object into a ClusterSet
        %                  so that targets in each cluster can be easily
        %                  addressed.
        %
        %OUTPUTS: newCSet A ClusterSet object representing the clustering
        %                 given by this DisjointSet object.

            numTar=length(DSObj.setArray);

            rootIdxList=zeros(DSObj.numberOfSets,1);
            curNumAdded=zeros(DSObj.numberOfSets,1);

            rootRankList=zeros(DSObj.numberOfSets,1);
            clusterEls=zeros(numTar,1);

            %First, record the indices and ranks of all of the root nodes.
            rootNodesFound=0;
            for curItem=1:length(DSObj.setArray)
                if(DSObj.setArray(curItem)<0)
                    rootNodesFound=rootNodesFound+1;
                    rootIdxList(rootNodesFound)=curItem;
                    rootRankList(rootNodesFound)=DSObj.numInTree(curItem);
                    if(rootNodesFound==DSObj.numberOfSets)
                        break;
                    end
                end
            end

            %Next, sort the root nodes by their index, so that children can
            %be efficiently assigned to them.
            [rootIdxList,idx] = sort(rootIdxList,1,'ascend');
            rootRankList=rootRankList(idx);

            %rootRankList currently holds all of the ranks of the root
            %indices in rootIdxList for the clusters. The elements of the
            %clusters are going to be put in one big array, clusterEls.
            %We are going to  create another array, offsetArray from
            %rootRankList such that each element is the offset needed to
            %access the first element of a given cluster in clusterEls. The
            %first element of offsetArray consequently is zero. The others
            %are a cumulative sum of the previous elements in rootRankList.
            offsetArray=zeros(DSObj.numberOfSets,1);
            offsetArray(2:end)=cumsum(rootRankList(1:(end-1)));

            %Now, we shall scan though all of the elements in
            %DSObj.setArray and add then to the appropriate cluster in
            %clusterEls depending on their rank.
            for curItem=1:numTar
                %Find the root node for the given node.
                rootNode=DSObj.find(curItem);
                %Find the index of the cluster corresponding to that root
                %node.
                [~, clusterIdx]=binSearch(rootIdxList,rootNode);
                %Add the node to the cluster.
                clusterEls(1+offsetArray(clusterIdx)+curNumAdded(clusterIdx))=curItem;
                curNumAdded(clusterIdx)=curNumAdded(clusterIdx)+1;
            end

            newCSet=ClusterSet(clusterEls,rootRankList,offsetArray);
        end

        function mergedRootIdx=unionRoots(DSObj,rootIdx1,rootIdx2)
        %%UNIONROOTS Merge two sets given their root nodes. If one wishes
        %            to merge sets given arbitrary (not necessarily root)
        %            nodes in the sets, then the function union should be
        %            used.
        %
        %INPUTS: DSObj      The implicitly passed calling object.
        %        rootIdx1   The root index of the first set to merge.
        %        rootIdx2   The root index of the second set to merge.
        %
        %OUTPUTS: mergedRootIdx The index of the root of the merged
        %                       cluster. It will either be rootIdx1 or
        %                       rootIdx2.
        %
        %This is a union-by-rank algorithm.
        %
            %If the roots are the same, then there is nothing to merge.
            if(rootIdx1==rootIdx2)
                %Do measurement bookeeping, if desired.
                mergedRootIdx=rootIdx1;
                return;
            end

            %If root 2 is deeper, than make root2 the new root.
            if(DSObj.setArray(rootIdx2)< DSObj.setArray(rootIdx1))
                DSObj.setArray(rootIdx1)=rootIdx2;
                DSObj.numInTree(rootIdx2)=DSObj.numInTree(rootIdx1)+DSObj.numInTree(rootIdx2);
                mergedRootIdx=rootIdx2;
            else
                %If they are the same height, then update the height of
                %rootIdx1.
                if(DSObj.setArray(rootIdx1)== DSObj.setArray(rootIdx2))
                    DSObj.setArray(rootIdx1)=DSObj.setArray(rootIdx1)-1;
                end
                DSObj.setArray(rootIdx2)=rootIdx1;
                DSObj.numInTree(rootIdx1)=DSObj.numInTree(rootIdx1)+DSObj.numInTree(rootIdx2);
                mergedRootIdx=rootIdx1;
            end

            %Merging reduces the number of sets by one.
            DSObj.numberOfSets=DSObj.numberOfSets-1;
        end

        function rootIdx=find(DSObj,idx)
        %FIND Perform a search with path compression. That is, all nodes
        %     traversed are set to point directly to the root of the tree.
        %     This means that future searches will be faster.
        %
        %INPUTS:  DSObj     The implicitly passed calling object.
        %         idx       The index of a node.
        %
        %OUTPUTS: rootIdx  The index of the root node of the cluster to
        %                  which the node with index idx belongs.
        %
        %The use of path compression would imply that the ranks associated
        %with the root nodes should be recalculated. However, it is noted
        %in Chapters 8.5 and 8.6 of
        %M. A. Weiss, Data Structures and Algorithms for C++, 2nd ed.
        %Reading, MA: Addison-Wesley, 1999.
        %that failing to update the rank does not worsen the theoretical
        %worst-case run time.
        %
            if(DSObj.setArray(idx)<0)
                rootIdx=idx;
            else
                rootIdx=DSObj.find(DSObj.setArray(idx));
                DSObj.setArray(idx)=rootIdx;
            end

        end
    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


function Rets = Simulate_Jump_Diffusion_Log_Returns_func( N_sim, dt, r, q, sigma, jumpModel, jumpParams )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Simulates Log Returns of Jump Diffusion Models with jumps (including simple Black-Scholes, no jumps)
% Returns: paths of dimension (N_sim, M+1), since they include S_0
%          ... Simulates N_sim paths, each row is a full path starting from S_0, ending with S_M (M+1 points in path)
% Author: Justin Lars Kirkby
%
% -----------------
% Params
% -----------------
% N_sim = # paths
% M = #time steps on [0,T], time step is dt=T/M, so each path has M+1 points
% T = time to maturity, ie path is on [0,T]
% S_0 = initial underlying value (e.g. S_0=100)
% r = interst rate (e.g. r = 0.05)
% q = dividend yield (e.g. q = 0.05)
% sigma = diffusion parameter (e.g. sigma = 0.2)
%
%===================================
% jumpModel: 0 = NoJumps, 1 = NormalJumps, 2 = DEJumps, 3 = MixedNormalJumps
%===================================
% jumpParams = paramters container containing all necessary params for models,
%            : if jumpModel = 0, no jump params are needed
%            : if jumpModel > 0, jumpParams must contain lambda, kappa, and any other model specific params (see below)

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

%==============================
% Initialize Jump Model Params and JumpFunc (function handle)
%==============================
%%% NOTE:  Jump Model is of the form in LOG space
%%% X(m+1) = X(m) + drift + Brownian Component + sum(Jumps on [m,m+1])
%%% By Jump we mean log(Y), e.g. in Merton Model, Jump ~ Normal (since we are in log space )

if jumpModel > 0 %ie if there are jumps in the model
    lambda = jumpParams.lambda;
    kappa  = jumpParams.kappa;

    Zeta = r - q - lambda*kappa;  %NOTE: we are redefining r to include compensation for jump component
    lamdt = lambda*dt;

    if jumpModel == 1 %Normal Jumps, e.g. Merton
        muJ  = jumpParams.muJ;
        sigJ = jumpParams.sigJ;
        JumpFunc = @(n) sum(muJ +sigJ*randn(n,1)); %Generates n independent jumps and sums them

    elseif jumpModel == 2 %Double Exponenial Jumps
        p_up = jumpParams.p_up;
        eta1 = jumpParams.eta1;
        eta2 = jumpParams.eta2;
        JumpFunc = @(n) sum(DoubleExpoRnd(n,p_up, eta1,eta2));

    elseif jumpModel == 3 %Mixed normal Jumps
        p_up = jumpParams.p_up;
        a1 = jumpParams.a1;  b1 = jumpParams.b1;
        a2 = jumpParams.a2;  b2 = jumpParams.b2;
        JumpFunc = @(n) sum(MixedNormalRnd(n, p_up, a1,b1,a2,b2));
    end
else
    Zeta = r - q;
end

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

Sigsqdt = sigma*sqrt(dt);
drift = (Zeta - .5*sigma^2)*dt;

if jumpModel == 0  % Just a drifted Brownian motion
    Rets = drift + Sigsqdt*randn(N_sim,1);
else
    Poi = PoissonRnd(N_sim, lamdt);  %Generate Poisson Column Vector of size N_Sim
    for n = 1:N_sim
        if Poi(n)>0
            Poi(n) = JumpFunc(Poi(n));  %JumpFunc(Poi(n)) sums up Poi(n) many jumps from the jump distribution
        end
    end

    Rets = drift + Poi + Sigsqdt*randn(N_sim,1);
end

end


text:
function w = w0_func(z,N)
% Based on:
% FADDEEVA   Faddeeva function
%   W = FADDEEVA(Z) is the Faddeeva function, aka the plasma dispersion
%   function, for each element of Z. The Faddeeva function is defined as:
%
%     w(z) = exp(-z^2) * erfc(-j*z)
%
%   where erfc(x) is the complex complementary error function.
%
%   W = FADDEEVA(Z,N) can be used to explicitly specify the number of terms
%   to truncate the expansion (see (13) in [1]). N = 16 is used as default.
%
%   Example:
%       x = linspace(-10,10,1001); [X,Y] = meshgrid(x,x);
%       W = faddeeva(complex(X,Y));
%       figure;
%       subplot(121); imagesc(x,x,real(W)); axis xy square; caxis([-1 1]);
%       title('re(faddeeva(z))'); xlabel('re(z)'); ylabel('im(z)');
%       subplot(122); imagesc(x,x,imag(W)); axis xy square; caxis([-1 1]);
%       title('im(faddeeva(z))'); xlabel('re(z)'); ylabel('im(z)');
%
%   Reference:
%   [1] J.A.C. Weideman, "Computation of the Complex Error Function," SIAM
%       J. Numerical Analysis, pp. 1497-1518, No. 5, Vol. 31, Oct., 1994
%       Available Online: http://www.jstor.org/stable/2158232

% Called by: crbs6.m, crbs7.m

if nargin<2, N = []; end
if isempty(N), N = 50; end

w = zeros(size(z)); % initialize output

%%%%%
% for purely imaginary-valued inputs, use erf as is if z is real
idx = real(z)==0; %
w(idx) = exp(-z(idx).^2).*erfc(imag(z(idx)));

if all(idx), return; end
idx = ~idx;

%%%%%
% for complex-valued inputs

% make sure all points are in the upper half-plane (positive imag. values)
idx1 = idx & imag(z)<0;
z(idx1) = conj(z(idx1));

M = 2*N;
M2 = 2*M;
k = (-M+1:1:M-1)'; % M2 = no. of sampling points.
L = sqrt(N/sqrt(2)); % Optimal choice of L.

theta = k*pi/M;
t = L*tan(theta/2); % Variables theta and t.
f = exp(-t.^2).*(L^2+t.^2);
f = [0; f]; % Function to be transformed.
a = real(fft(fftshift(f)))/M2; % Coefficients of transform.
a = flipud(a(2:N+1)); % Reorder coefficients.

Z = (L+1i*z(idx))./(L-1i*z(idx));
p = polyval(a,Z); % Polynomial evaluation.
w(idx) = 2*p./(L-1i*z(idx)).^2 + (1/sqrt(pi))./(L-1i*z(idx)); % Evaluate w(z).

% convert the upper half-plane results to the lower half-plane if necesary
w(idx1) = conj(2*exp(-z(idx1).^2) - w(idx1));

w=-w*sqrt(-1)*pi; % converts to 'w0_func' used in crbs_molecular.m

author: Sable

repo: mcbench-benchmarks


text:
function vect = convertxyz2Vect(xyz)

% CONVERTXYZ2VECT convert an xyz matrix into a vector.
%
%	Description:
%
%	channelsFinal = gradientDescent(skel, initialChannels, finalChannels,
%	error,step)
%
%	 Returns:
%	 VECT - a line vector
%	 Arguments:
%	 XYZ - the xyz matrix to be transformed

for i=1:size(xyz,1)
    for j=1:size(xyz,2)
        vect((i-1)*3 + j) = xyz(i,j);
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_ft_combineplanar

% WALLTIME 00:20:00
% MEM 2gb
% DEPENDENCY ft_combineplanar ft_preprocessing ft_timelockanalysis ft_prepare_neighbours ft_megplanar

% template loading should be modified once template MEG data are available
% (cf.issue #1834)
subjectfilename = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.ds');

cfg            = [];
cfg.dataset    = subjectfilename;
cfg.channel    = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
data = ft_preprocessing(cfg);

cfg = [];
data = ft_timelockanalysis(cfg,data);

cfg              = [];
cfg.feedback     = 'no';
cfg.method       = 'template';
cfg.neighbours   = ft_prepare_neighbours(cfg, data);
cfg.planarmethod = 'sincos';
data = ft_megplanar(cfg, data);

cfg = [];
cfg.updatesens = 'no';
dataout = ft_combineplanar(cfg,data);


text:
% Copyright (C) 2017-2018 Titus Cieslewski, RPG, University of Zurich,
%   Switzerland
%   You can contact the author at <titus at ifi dot uzh dot ch>
% Copyright (C) 2017-2018 Siddharth Choudhary, College of Computing,
%   Georgia Institute of Technology, Atlanta, GA, USA
% Copyright (C) 2017-2018 Davide Scaramuzza, RPG, University of Zurich,
%   Switzerland
%
% This file is part of dslam_open.
%
% dslam_open 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.
%
% dslam_open 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 dslam_open. If not, see <http://www.gnu.org/licenses/>.

% Visualize trajectory outputs
import gtsam.*
close all;
nr_robots = 10;
dataDir = 'data';
cmap = hsv(nr_robots);

    figure(1);
    cla; hold on;
for robot_i = 1:nr_robots
    robot_i
    initial_file = sprintf('%s/%d_optimized.g2o', dataDir, robot_i-1)
    [graph, initial] = gtsam.readG2o(initial_file, true);

    plot3DSubGraph(graph, initial, '-',  cmap(robot_i,:), '');
    hold on;
end

axis equal;

author: uzh-rpg


text:
% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

%GEOMETRICAL PARAMETERS

clear;

%Length from the bearing to the hinge
l0 = 27.5;

%Length of the vertical wall that supports the hinge
h0 = 6.5;

%Lenght of the vertical support of the bearing
hb = 4.0;

%Bearing radius
br = 1.1;

%Array of different max and min positions of the bearing
h1a = [6.5, 7.1, 7.7, 8.3, 8.9];
h2a = [20.5, 20.5, 20.5, 20.5, 20.5];

%Number of cycles per cam turn (1 to 3).
nc = 1;

%Minimum radius of the camshaft
rmin = 5;

%BREATHING CYCLE PARAMETERS

%Duration of the inhale cycle / duration of the whole cycle
lambda1 = 0.500;
lambda2 = 0.04;

%Soft transition between inhale and exhale cycles
dpsi21 = 0.05;
dpsi12 = 0.04;

%Adjust parameters of the inhale curve
ga1 = 3.1;
gb1 = .85;
ff1 = 50;

%Adjust parametes of the exhale curve
ga2 = 3.1;
gb2 = .80;
ff2 = 50;

%GENERATION OF THE CURVES AND THE CAMSHAFT

%Generation of the geometry

l = sqrt(l0^2+hb^2);

for i = 1:numel(h1a);

    ymin(i) = h1a(i);
    ymax(i) = h2a(i);

    alphamin(i) = acos((h0 - h1a(i)) / l);
    alphamax(i) = acos((h0 - h2a(i)) / l);

    xmin(i) = l*sin(alphamin(i));
    xmax(i) = l*sin(alphamax(i));

    d(i) = sqrt((xmin(i)-xmax(i))^2 + (ymin(i)-ymax(i))^2);

    alphatan(i) = (alphamin(i) + alphamax(i)) / 2;

    xtan(i) = l*sin(alphatan(i));
    ytan(i) = h0 - l*cos(alphatan(i));

    xsup(i) = xtan(i) + (d(i)/2)*cos(alphatan(i));
    xinf(i) = xtan(i) - (d(i)/2)*cos(alphatan(i));

    ysup(i) = ytan(i) + (d(i)/2)*sin(alphatan(i));
    yinf(i) = ytan(i) - (d(i)/2)*sin(alphatan(i));

    xcam(i) = xtan(i) + (d(i)/2 + rmin + br)*cos(alphatan(i));
    ycam(i) = ytan(i) + (d(i)/2 + rmin + br)*sin(alphatan(i));

end;

%Time increment
dt = 0.01;

%time coordinate during the whole cycle
theta = 0:dt:2*pi;

%Generation of the soft transition between inhale and exhale curves
psi1 = [];
psi2 = [];

for i = 1:numel(theta);

    if theta(i) < (lambda2-dpsi21/2)*2*pi;
        psi1(i) = 0;

    elseif theta(i) < (lambda2+dpsi21/2)*2*pi;
        psi1(i) = 0.5 + 0.5*cos((theta(i)-(lambda2+dpsi21/2)*2*pi)/(2*dpsi21));

    else theta(i) < (lambda1-dpsi12/2)*2*pi;
        psi1(i) = 1;

    end;

    psi2(i) = 1 - psi1(i);

end;

%Generation of the complete breathing cycle rho(theta)

rho = [];
rho1 = [];
rho2 = [];
rho2next = [];

rhomin = 1000;
rhomax = 0;

for i = 1:numel(theta);

    %Inhale curve
    rho1(i) = ff1*gampdf(theta(i), ga1, gb1);
    rho1next(i) = ff1*gampdf(theta(i)+2*pi, ga1, gb1);

    %Exhale curve
    rho2(i) = ff2*gampdf(theta(i), ga2, gb2);
    rho2next(i) = ff2*gampdf(theta(i)+2*pi, ga2, gb2);

    rho(i) = max(rho1(i), rho1next(i));

    %Capturing min and max in order to generate the normalized curve
    if rho(i) > rhomax
        rhomax = rho(i);
    end;

    if rho(i) < rhomin
        rhomin = rho(i);
    end;

end;

%Generation of a normalised curve and camshaft

rhonorm = [];
rhocam = [];

for i = 1:numel(theta);

    rhonorm(i) = (rho(i)-rhomin)/(rhomax-rhomin);

    for n = 1:numel(d)

        rhocam(n,i) = rmin + rhonorm(i)*d(n);

    end;

end;

%Generation of the first derivate of the camshaft geometry to analize and
%validate the design

drho = [];
drhonorm = [];
drhocam = [];

a = rmin;
b = 1/dt;

for i = 1 : numel(rho)-1;

    drho(i) = b*(rho(i+1)-rho(i));
    drhonorm(i) = b*(rhonorm(i+1)-rhonorm(i));
    drhocam(i) = b*(rhocam(i+1)-rhocam(i)) + a;

end;

drho(numel(rho)) = b*(rho(1)-rho(numel(rho)));
drhonorm(numel(rhonorm)) = b*(rhonorm(1)-rhonorm(numel(rhonorm)));
drhocam(numel(rhocam)) = b*(rhocam(1)-rhocam(numel(rhocam))) + a;

%X and Y coordinates during two cycles (for 2-cycle camshaft plot)

theta2 = [theta/2, (2*pi+theta)/2];

rhocam2 = [rhocam, rhocam];
drhocam2 = [drho, drho];

%X and Y coordinates during three cycles (for 3-cycle camshaft plot)

theta3 = [theta/3, (2*pi+theta)/3, (4*pi+theta)/3];

rhocam3 = [rhocam, rhocam, rhocam];
drhocam3 = [drho, drho, drho];

%GENERATION OF THE PLOTS

%Trying to print it out in real size (FAIL)
%set(gcf,'PaperUnits','centimeters');
%set(gcf,'PaperSize',[42 29.7]);

fpos = figure('Name', 'Dimensions', 'Units', 'centimeters', 'NumberTitle', 'off');
set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 28, 20]);
fpos = gcf;
hold on

xlim([-5 35]);
ylim([0 30]);

for i = 1:numel(xmin)

    plot([0, xmin(i)], [h0, ymin(i)], '-x')
    plot([0, xmax(i)], [h0, ymax(i)], '-x')

    plot([0, xtan(i)], [h0, ytan(i)], '--xb')

    plot([xsup(i), xinf(i)], [ysup(i), yinf(i)], '--xr')
    scatter(xcam(i), ycam(i))

end;

hold off

%Plot of the breathing cycle interpolation by curves

fcycle = figure('Name', 'Breath Cycle curves', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
plot(theta, psi1, '-k')
plot(theta, psi2, '-k')
plot(theta, rho1, '-g')
plot(theta, rho1next, '-g')
plot(theta, rho2, '-g')
plot(theta, rho2next, '-g')
%plot(theta, drho, '-r')
plot(theta, rho, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of the normalized breathing cycle and first derivate

fnorm = figure('Name', 'Normalized Breath Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
%plot(theta, drhonorm, '-r')
plot(theta, rhonorm, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of a 1-cycle camshaft

for n = 1:numel(d)

    fcam1 = figure('Name', '1-Cycle Camshaft', 'Units', 'centimeters', 'NumberTitle', 'off');
    %hold on

    %polarplot(theta, drhocam, '-r')
    polar(theta, rhocam(n,:), '-b')

    %rlim([0 25]);
    hold off

    set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
    set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
    fcam3 = gcf;

end

% %Plot of a 2-cycle camshaft
%
% fcam2 = figure('Name', '2-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');
%
% %hold on;
% %polar(theta, drhocam2, '-r')
% polar(theta2, rhocam2, '-b')
% %hold off;
%
% set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
% set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
% fcam3 = gcf;
%
%
% %Plot of a 3-cycle camshaft
%
% fcam3 = figure('Name', '3-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');
%
% %hold on;
% %polar(theta,drhocam3)
% polar(theta3, rhocam3)
% %hold off;
%
% set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
% set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
% fcam3 = gcf;

% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

author: ProtofyTeam

repo: OxyGEN

save_path: github-repos/MATLAB/ProtofyTeam-OxyGEN


%function [varargout] = date_(varargin)
function     [dt_]  = ...
    date_(x, time, strfreq, t0, t1, options)
%    date_(x, time, freq, ndfy, ndfm, ndly, ndlm, turnphase, phase, cycle, thresh, nrep, complete)
% st,trinary,notentp,dura_,ampl_,cumm_,excc_,durcv_,amplcv_,exccv_
%***********************************************************************************************/
%*                                                                                */
%*  Modified BBQ program                                                                       */
%*                                                                                */
%*  Computes turning points and imposes restrictions in one step                               */
%*                                                                                */
%*    NB: rule is peak greater than turnphase on either side,                                  */
%*             trough less than turnphase on either side.           (can be modified)          */
%*    Restriction: Phase quarters/months, Cycle  quarters/months Alternating Troughs and Peaks,*/
%*                  if two peaks(troughs) in a row chooses highest(lowest), Trough must be     */
%*                  lower then preceeding peak, and No turning point within phase length       */
%*                  of end points                                                              */
%*                                                                                */
%*                                                                                */
%*                                                                                */
%*    Date: 15th November 2005                                                                  */
%*     Author: James Engel - Code modified from Adrian Pagan and Don Harding's BBQ code     */
%*
%*     further adjustments  made  by  F. Canova
%***********************************************************************************************/
%

ndfy = t0(1); % year start
ndfm = t0(2); % month/quarter start,
ndly = t1(1); % year start,
ndlm = t1(2); % month/quarter start,

if length(x) ~= length(time)
    error('data (x) and time must have the same length');
end

if strcmp(strfreq,'q') == 1
    freq =1;
    % default settings cycle parameters q data
    turnphase   = 2;
    phase       = 2;        % censoring rules
    cycle       = 5;        % lenght  of cycle

    if ndfm>4 || ndlm>4
        error('Something went wrong. Frequency is quarterly; therefore, starting and ending quarters cannot be larger than 4!');
    end

elseif strcmp(strfreq,'m') == 1
    freq =2;
    % default settings cycle parameters m data
    turnphase   = 6;
    phase       = 6;        % censoring rules
    cycle       = 15;        % lenght  of cycle

else
    error('strfreq must be ''q'' (quarterly) or ''m'' (monthly)');
end

notentp=0;

if freq==1 % q
    nd= 4*(ndly-1-ndfy) + (5-ndfm) + (ndlm);     % Number of data points %
    stepsize = 1/4;
elseif freq==2 % m
    nd= 12*(ndly-1-ndfy) + (13-ndfm) + (ndlm);
    stepsize = 1/12;
end

thresh      = 10.4;     % bypasses phase and cycle restriction if peak to trough is > than thresh
nrep        = 1;        % 1 if analyze real data
complete    = 1;        % if= 1- use complete cycles,if =0 -use incomplete cycles (excess still on complete cycle)

if nargin > 5
    if isfield(options,'turnphase')==1
        turnphase = options.turnphase;
    end
    % censoring rules
    if isfield(options,'phase')==1
        phase = options.phase;
    end
    % lenght of cycle
    if isfield(options,'cycle')==1
        cycle = options.cycle;
    end
    % bypasses phase and cycle restriction if peak to trough is > than thresh
    if isfield(options,'thresh')==1
        thresh = options.thresh;
    end
    % 1 if analyze real data
    if isfield(options,'nrep')==1
        nrep = options.nrep;
    end
    % if= 1- use complete cycles,if =0 -use incomplete cycles (excess still on complete cycle)
    if isfield(options,'complete')==1
        complete = options.complete;
    end

end

pdcv=0;tdcv=0;pacv=0;tacv=0;pdm=0;pdma=0;tdm=0;tdma=0;pdcm=0;tdcm=0;
pdem=0;tdem=0;tdema=0;pdema=0;epcv=0;etcv=0;

[bcp5, bct5, nbp, nbt] = rawall(x(1:nd),turnphase,nd,phase,cycle,thresh);   % calculates turning points with restrictions %

if nbp+nbt<=2
    notentp=notentp+1;
else

    %ntr=nbt;npk=nbp;

    nr  = [nbt;nbp];
    nv  = max(nr);

    %               pdc=zeros(nv,1);tdc=zeros(nv,1);
    %               pda=zeros(nv,1);tda=zeros(nv,1);td=zeros(nv,1);pd=zeros(nv,1);
    pdc = zeros(nv,1);
    tdc = zeros(nv,1); % pde=zeros(nv,1);tde=zeros(nv,1);
    %               tdea=zeros(nv,1);pdea=zeros(nv,1);

    % calculate peak to trough durations & amps
    % p st&s for peaks,t for troughs,p gives cntractions,t expansions
    % code is pd,td is durations; pda,tda is amps; pdc,tdc is cum move pde, tde exces
    % excess is measured differently to avoid case that amps is close to zero in part cycle so that denom can become neg.so use cum movements as denom
    % now vs triangle in early paper

    if bcp5(1,1) < bct5(1,1)      % Peaks are first

        nr  = [nbt;nbp];
        r   = nbt;
        pd  = bct5(1:r,1)-bcp5(1:r,1);
        pda = x(bct5(1:r,1))-x(bcp5(1:r,1));
        k   = 1;

        while k<=r
            %pdc(k)=sumc(x(bcp5(k,1):bct5(k,1),1)-x(bcp5(k,1),1));
            pdc(k)=sum(x(bcp5(k,1):bct5(k,1),1)-x(bcp5(k,1),1));
            k=k+1;
        end
    else                      % troughs are First

        r   = nbt-1;
        pd  = bct5(2:r+1,1)-bcp5(1:r,1);
        pda = x(bct5(2:r+1,1))-x(bcp5(1:r,1));

        k=1;
        while k<=r
            %pdc(k)=sumc(x(bcp5(k,1):bct5(k+1,1),1)-x(bcp5(k,1),1));
            pdc(k)=sum(x(bcp5(k,1):bct5(k+1,1),1)-x(bcp5(k,1),1));
            k=k+1;
        end

        %r1 = r;

    end

    % calculate trough to peak durations & amplitudes

    if bct5(1,1) < bcp5(1,1)        %  Troughs are first
        r   = nbp;
        td  = bcp5(1:r,1)-bct5(1:r,1);
        tda = x(bcp5(1:r,1))-x(bct5(1:r,1));
        k   = 1;
        while k<=r
            %tdc(k)=sumc(x(bct5(k,1):bcp5(k,1),1)-x(bct5(k,1),1));
            tdc(k)=sum(x(bct5(k,1):bcp5(k,1),1)-x(bct5(k,1),1));
            k=k+1;
        end

    else                      % peaks are First
        r   = nbp-1;
        td  = bcp5(2:r+1,1)-bct5(1:r,1);
        tda = x(bcp5(2:r+1,1))-x(bct5(1:r,1));

        k=1;
        while k<=r

            %tdc(k)=sumc(x(bct5(k,1):bcp5(k+1,1),1)-x(bct5(k,1),1));
            tdc(k)=sum(x(bct5(k,1):bcp5(k+1,1),1)-x(bct5(k,1),1));
            k=k+1;
        end

    end
    pdc=pdc(1:size(pd,1));
    tdc=tdc(1:size(td,1));

    % compute excesses
    %excess is percentage of triangle area
    za      = (pd.*pda)/2;
    pde     = 100*(pdc-za-.5*pda)./za;

    pdea    = 100*(pdc-((pd.*pda)/2))./za;
    za      = (td.*tda)/2;
    tde     = 100*(tdc-za-.5*tda)./za;

    tdea    = 100*(tdc-((td.*tda)/2))./za;

    %***********************************************************************************************/

    if complete == 0    % switch: 1-only use complete cycles,0-use incomplete cycles (excess still on complete cycle)%

        bct5u   = bct5;
        bcp5u   = bcp5;

        if bcp5(1,1) < bct5(1,1)     % modifies code to include incomplete cycles %

            bct5u  = [1;bct5];

        elseif 	bct5(1,1) < bcp5(1,1)

            bcp5u = [1;bcp5];

        end

        nbtu = size(bct5u,1);
        nbpu = size(bcp5u,1);

        if bcp5u(nbpu,1) < bct5u(nbtu,1)    % modifies code to include incomplete cycles %

            bcp5u  = [bcp5u;nd];

        elseif 	bct5u(nbt,1) < bcp5u(nbp,1)

            bct5u = [bct5u;nd];

        end

        ntr     = size(bct5u,1);
        npk     = size(bcp5u,1);
        % nr      = [ntr;npk];
        % nv=max(nr);

        %        pdc=zeros(nv,1);tdc=zeros(nv,1);
        %        pda=zeros(nv,1);tda=zeros(nv,1);td=zeros(nv,1);pd=zeros(nv,1);
        %        pdc=zeros(nv,1);tdc=zeros(nv,1);

        if bcp5u(1,1) < bct5u(1,1)       % Peaks are first

            %nr  = [ntr;npk];
            r   = ntr;
            pd  = bct5u(1:r,1)-bcp5u(1:r,1);
            pda = x(bct5u(1:r,1))-x(bcp5u(1:r,1));
            k   = 1;
            while k<=r
                %pdc(k)=sumc(x(bcp5u(k,1):bct5u(k,1),1)-x(bcp5u(k,1),1));
                pdc(k)=sum(x(bcp5u(k,1):bct5u(k,1),1)-x(bcp5u(k,1),1));
                k=k+1;
            end
        else                      % troughs are First
            r   = ntr-1;
            pd  = bct5u(2:r+1,1)-bcp5u(1:r,1);
            pda = x(bct5u(2:r+1,1))-x(bcp5u(1:r,1));
            k   = 1;
            while k<=r
                pdc(k)=sum(x(bcp5u(k,1):bct5u(k+1,1),1)-x(bcp5u(k,1),1));
                %pdc(k)=sumc(x(bcp5u(k,1):bct5u(k+1,1),1)-x(bcp5u(k,1),1));
                k=k+1;
            end

%             r1=r;

        end

        % calculate trough to peak durations & amplitudes

        if bct5u(1,1) < bcp5u(1,1)       %  Troughs are first

            r   = npk;
            td  = bcp5u(1:r,1)-bct5u(1:r,1);
            tda = x(bcp5u(1:r,1))-x(bct5u(1:r,1));
            k   = 1;
            while k<=r
                %tdc(k)=sumc(x(bct5u(k,1):bcp5u(k,1),1)-x(bct5u(k,1),1));
                tdc(k)=sum(x(bct5u(k,1):bcp5u(k,1),1)-x(bct5u(k,1),1));
                k=k+1;
            end

        else                     % peaks are First
            r   = npk-1;
            td  = bcp5u(2:r+1,1)-bct5u(1:r,1);
            tda = x(bcp5u(2:r+1,1))-x(bct5u(1:r,1));

            k=1;
            while k<=r
                %tdc(k)=sumc(x(bct5u(k,1):bcp5u(k+1,1),1)-x(bct5u(k,1),1));
                tdc(k)=sum(x(bct5u(k,1):bcp5u(k+1,1),1)-x(bct5u(k,1),1));
                k=k+1;
            end

        end
        pdc=pdc(1:size(pd,1));
        tdc=tdc(1:size(td,1));

    end

    %***********************************************************************************************/

    % cumulate.when nrep=1 then gives raw data,otherwsise sums over monte carlo amounts

    pdcv    = pdcv + stdc(pd)/meanc(pd);       % compute cv's
    tdcv    = tdcv + stdc(td)/meanc(td);
    pacv    = pacv + stdc(pda)/meanc(pda);
    tacv    = tacv + stdc(tda)/meanc(tda);
    epcv    = epcv + stdc(pde)/meanc(pde);  %  cv of xss
    etcv    = etcv + stdc(tde)/meanc(tde);

    pdm     = pdm+meanc(pd);          % durations
    pdma    = pdma+meanc(pda);

    tdm     = tdm+meanc(td);          % durations
    tdma    = tdma+meanc(tda);

    pdcm    = pdcm+meanc(pdc);       % cumulative
    tdcm    = tdcm+meanc(tdc);

    % pdem=pdem+meanc(pde);
    tdem    = tdem+meanc(tde);       % xss
    tdema   = tdema+meanc(tdea);
    pdem    = pdem+meanc(pde);
    pdema   = pdema+meanc(pdea);

end

if nrep==1&&(nbp+nbt>2)

    %nbp and nbt=number of peaks and number of troughs %

    disp('peaks at')
    fprintf('%.2f\n', time(bcp5))
    %disp(bcp5(1:nbp))

    disp('troughs at')
    fprintf('%.2f\n', time(bct5))
    %disp(bct5(1:nbt))

    % total  number of  turning  points
    ttp=length(bcp5)+length(bct5);

    % collecting  turning  points in  trinary  indicator
    trinary=zeros(length(x),1);
    trinary(bcp5)=1;
    trinary(bct5)=-1;

end

if nrep==1&&(nbp+nbt>2)
    %remove below if want states ...
    %st=zeros(nd,1);

    [st] = states(bcp5,bct5,nbp,nbt,nd);
    dt_.st(:,nrep) = st;

    %determine obs in which states have been completed%...
    na  = min([bct5(1);bcp5(1)])';
    nb  = max([bct5(nbt);bcp5(nbp)])';
    %nb-na+1;
    z   = [x(na:nb) st(na:nb)];

    %disp('dated series')
    %z1= [time(time_start+1:time_end)' x st];
    %format  short  g
    %round(z1, 2)
    %%fprintf('%.2f\n', '%.2f\n', '%.2f\n', z1)

end
nrep1   = nrep-notentp;

% Compute final statistics
% duration contractions/duration expansions
dt_.dura = [pdm/nrep1 tdm/nrep1];
% amplitudes contractions/amplitude expansions
dt_.ampl  = 100*[pdma/nrep1 tdma/nrep1];
% cumulative contractions/cumulative expansions
dt_.cumm  = 100*[pdcm/nrep1 tdcm/nrep1];
% excess movements percent of triangle area
% contractions/expansions
dt_.excc  = [pdem/nrep1 tdem/nrep1];
% cv of durations contractions/expansions
dt_.durcv = [pdcv/nrep1 tdcv/nrep1];
% cv of amplitude contractions/expansions
dt_.amplcv = [pacv/nrep1 tacv/nrep1];
% cv of excess movements contractions/expansions
dt_.exccv = [epcv/nrep1 etcv/nrep1];
%
dt_.trinary = trinary;
%
dt_.notentp = notentp;

function m=meanc(x)
% get mean of each column
m=mean(x);
if size(m,2)>1
   m=m';
end

function m=stdc(x)
% get standard deviation of each column
m=std(x);
if size(m,2)>1
   m=m';
end


text:
function [modelFile, model] = trainSceneModel(imSet, paths, typ, classMapping, sceneMapping, objectDir),

	train = imSet{1};
	val = imSet{2};

	pt = getMetadata(sceneMapping);
	numScene = length(pt.sceneName);
	sceneName = pt.sceneName;

	pt = getMetadata(classMapping);
	numObjectClass = length(pt.className);

	imSet = {train, val}; useVal = 1;
	imList{1} = getImageSet(imSet{1});
	imList{2} = getImageSet(imSet{2});
	trainingParam.fileSuffix = sprintf('tr-%s_val-%s_useVal-%d', imSet{1}, imSet{2}, useVal);

	trainingParam.featureParam = getSceneParams(typ, struct('dirName', objectDir, 'numObjectClass', numObjectClass));

	trainingParam.classifierParam = getClassifierParam('svm-scene', struct('useVal', useVal, 'numClass', numScene));
	trainingParam.gtParam = struct('sceneMapping', sceneMapping);

	% Load the features.
	for i = 1:2,
		clear f;
		parfor j = 1:length(imList{i}),
			f{j} = getSceneFeatures(imList{i}{j}, paths, trainingParam.featureParam);
		end
		F{i} = cat(2, f{:});
		gt{i} = getGroundTruth(imList{i}, 'scene', '', sceneMapping);
	end

	for i = 1:2,
		X{i} = F{i};
		Y{i} = gt{i};
		W{i} = ones(size(Y{i}));
	end

	svmParam = trainingParam.classifierParam;

	model = svmMulticlassTrain(X, Y, W, svmParam);
	modelFile = fullfile(paths.modelDir, strcat(sprintf('scene-%s_%s', trainingParam.featureParam.featureStr, trainingParam.fileSuffix), '.mat'))

	% Write the scene model!
	save(modelFile, 'model', 'trainingParam', 'F', 'gt');
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function mtrPlotLengthPenalty(length_absorption)
%
% mtrPlotLengthPenalty([length_absorption=0.8])
%

if(~exist('length_absorption','var')|isempty(length_absorption))
    length_absorption = 0.800;
end

FA = [0.88, 0.7664, 0.50, 0.30, 0.15];
cor = {'k', 'g', 'r', 'm', 'c'};
stdev = 1.0978*exp(-1.9567*FA) - 0.1437;
%maxprior = 1/(2*pi)*normpdf(0,0,0.628)/(normcdf(pi,0,0.628)-0.5);
maxprior = normpdf(0,0,0.628)/(normcdf(pi,0,0.628)-0.5);
maxlike = normpdf(0,0,stdev)./(normcdf(pi/2,0,stdev)-0.5);
maxcortexlike = normpdf(0,0,0.15)./(normcdf(pi/2,0,0.15)-0.5);
%length_penalty = pi/5;
length_penalty = 1-length_absorption;

x = [50:2:250];
for i = 1:length(maxlike)
    log_segment_score = log(maxprior * maxlike(i) * length_penalty);
    y = (x-4)/2*log_segment_score + 2*log(0.999)*maxprior*maxcortexlike;
    plot(x,y,cor{i},'LineWidth',2)
    hold on;
end
hold off;

%legend('STT Path','FA = 0.88','FA = 0.77','FA = 0.50','FA = 0.30','FA = 0.15');
legend('FA = 0.88','FA = 0.77','FA = 0.50','FA = 0.30','FA = 0.15');
% Find where stddev will make the function not penalize or reward
stdsearch = [0:0.01:pi/2];
maxlike = normpdf(0,0,stdsearch)./(normcdf(pi/2,0,stdsearch)-0.5);
abs_post = abs(maxprior*maxlike*length_penalty-1);
[junk, mi] = min(abs_post);
stdev_zero = stdsearch(mi);
fa_zero = log((stdev_zero+0.1437)/1.0978)/(-1.9567)


text:
function [params,X,means,wX,wM,dwM] = pre_sphere(params, X, dim)
% Spheres the data with the eigenvalue decomposition of the covariance matrix
%   [params,X,means,wX,wM,dwM] = pre_sphere(params, X, dim)
%     X    Data to be sphered
%     dim  Requested result data dimension
%     wX   Sphered data
%     wM   Sphering matrix
%     dwM  De-sphering matrix

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<2
    params.name = 'PCA sphering';
    params.description = 'Spheres the data with the eigenvalue decomposition of the covariance matrix.';
    return;
end

if nargin<3; dim=0; end

[xdim,tdim] = size(X);
% removing the mean
means = mean(X,2);
X = X - repmat(means,1,tdim);

[wM, dwM] = dss_sphere(X, dim);

wX = wM * X;


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun  vpa (@var{x})
%% @defunx vpa (@var{x}, @var{n})
%% Create a variable-precision floating point number.
%%
%% @var{x} can be a string, a sym or a double.  Example:
%% @example
%% @group
%% x = vpa('1/3', 32)
%%   @result{} x = (sym) 0.33333333333333333333333333333333
%% a = sym(1)/3;
%% x = vpa(a, 32)
%%   @result{} x = (sym) 0.33333333333333333333333333333333
%% @end group
%% @end example
%%
%% If @var{n} is omitted it defaults to the current value of
%% @code{digits()}.
%%
%% Be careful when creating a high-precision float from a
%% double as you will generally only get 15 digits:
%% @example
%% @group
%% vpa(1/3)
%%   @result{} (sym) 0.3333333333333333148296162...
%% vpa(sqrt(2));
%% ans^2
%%   @result{} (sym) 2.0000000000000002734323463...
%% @end group
%% @end example
%%
%% For the same reason, passing numbers with decimal points
%% may produce undesirable results:
%% @example
%% @group
%% vpa(0.1)
%%   @result{} (sym) 0.1000000000000000055511151...
%% @end group
%% @end example
%%
%% Instead, enclose the decimal number in a string:
%% @example
%% @group
%% vpa('0.1')
%%   @result{} (sym) 0.10000000000000000000000000000000
%% @end group
%% @end example
%%
%% Very simple expressions can also be enclosed in quotes:
%% @example
%% @group
%% vpa('sqrt(2)')
%%   @result{} (sym) 1.4142135623730950488016887242097
%% @end group
%% @end example
%%
%% But be careful as this can lead to unexpected behaviour, such as
%% low-precision results if the string contains decimal points:
%% @example
%% @group
%% vpa('cos(0.1)')
%%   @print{} warning: string expression involving decimals is
%%   @print{}          dangerous, see "help vpa"...
%%   @result{} ans = (sym) 0.99500416527802...
%% @end group
%% @end example
%%
%% Instead, it is preferrable to use @code{sym} or @code{vpa} on the
%% inner-most parts of your expression:
%% @example
%% @group
%% cos(vpa('0.1'))
%%   @result{} (sym) 0.99500416527802576609556198780387
%% vpa(cos(sym(1)/10))
%%   @result{} (sym) 0.99500416527802576609556198780387
%% @end group
%% @end example
%%
%% @seealso{sym, vpasolve, digits}
%% @end defun

function r = vpa(x, n)

  if (nargin == 1)
    n = digits();
  elseif (nargin ~= 2)
    print_usage ();
  end

  if (isa(x, 'sym'))
    cmd = {
        'x, n = _ins'
        'return sympy.N(x, n),' };
    r = pycall_sympy__ (cmd, x, n);
  elseif (ischar (x))
    x = strtrim (x);
    isfpnum = ...
      ~isempty (regexp (x, '^[-+]*?[\d_]*\.[\d_]*(e[+-]?[\d_]+)?[ij]?$'));
    if (~isfpnum && ~isempty (strfind (x, '.')))
      warning ('OctSymPy:vpa:precisionloss', ...
               'string expression involving decimals is dangerous, see "help vpa"')
    end
    if (isfpnum && any (strcmp (x(end), {'i', 'j'})))
      r = sym (1i)*vpa (x(1:end-1), n);
      return
    end
    if (any (strcmp (x, {'inf', 'Inf', '+inf', '+Inf'})))
      x = 'S.Infinity';
    elseif (strcmp (x, '-inf') || strcmp (x, '-Inf'))
      x = '-S.Infinity';
    elseif (strcmp (x, 'I'))
      x = 'Symbol("I")';
    elseif (any (strcmp (x, {'1i', '1j'})))
      x = 'S.ImaginaryUnit';
    end
    % Want Float if its '2.3' but N if its 'sqrt(2)'
    cmd = {
        'x, n = _ins'
        'try:'
        '    return sympy.Float(x, n),'
        'except ValueError:'
        '    return sympy.N(x, n),' };
    r = pycall_sympy__ (cmd, x, n);
  elseif (isfloat (x) && ~isreal (x))
    r = vpa (real (x), n) + sym (1i)*vpa (imag (x), n);
  elseif (isfloat(x) && isscalar(x) == 1)
    if (isnan (x))
      x = 'S.NaN';
    elseif (isinf (x) && x < 0)
      x = '-S.Infinity';
    elseif (isinf (x))
      x = 'S.Infinity';
    elseif (isequal (x, pi))
      x = 'S.Pi';
    elseif (isequal (x, -pi))
      x = '-S.Pi';
    elseif (isequal (x, exp (1)))
      x = exp (sym (1));
    elseif (isequal (x, -exp (1)))
      x = -exp (sym (1));
    end
    cmd = {
      'x, n = _ins'
      'return sympy.N(x, n)' };
    r = pycall_sympy__ (cmd, x, n);
  elseif (isinteger(x) && isscalar(x) == 1)
    cmd = {
        'x, n = _ins'
        'return sympy.N(x, n),' };
    r = pycall_sympy__ (cmd, x, n);
  elseif (~isscalar(x))
    cmd = { sprintf('return sympy.ZeroMatrix(%d, %d).as_mutable(),', size(x)) };
    r = pycall_sympy__ (cmd);
    for i=1:numel(x)
      r(i) = vpa(x(i), n);
    end
  else
    x
    class(x)
    error('conversion to vpa with those arguments not (yet) supported');
  end

end

%!test
%! a = vpa(0, 4);
%! b = double(a);
%! assert(b == 0)

%!test
%! a = vpa(pi, 4);
%! b = sin(a);
%! assert(abs(double(b)) < 1e-4)

%!test
%! % vpa from double is ok, doesn't warn (c.f., sym(2.3))
%! a = vpa(2.3);
%! assert(true)

%!test
%! % vpa from double not more than 16 digits
%! a = vpa(sqrt(pi), 32);
%! b = sin(a^2);
%! assert(abs(double(b)) > 1e-20)
%! assert(abs(double(b)) < 1e-15)

%!test
%! a = vpa(sym(pi), 32);
%! b = sin(a);
%! assert(abs(double(b)) < 1e-30)

%!test
%! a = vpa(sym(pi), 256);
%! b = sin(a);
%! assert(abs(double(b)) < 1e-256)

%!test
%! % pi str
%! a = vpa('pi', 32);
%! b = sin(a);
%! assert(abs(double(b)) < 1e-32)

%!test
%! % pi str
%! a = vpa('pi', 32);
%! b = vpa(sym('pi'), 32);
%! assert (double (a - b) == 0)

%!test
%! spi = sym(pi);
%! a = vpa(spi, 10);
%! b = double(a);
%! assert(~isAlways(spi == a))

%!test
%! % matrix of sym
%! a = [sym(pi) 0; sym(1)/2 1];
%! b = [pi 0; 0.5 1];
%! c = vpa(a, 6);
%! assert(max(max(abs(double(c)-b))) < 1e-6)

%!test
%! % matrix of double
%! b = [pi 0; 0.5 1];
%! c = vpa(b, 6);
%! assert(max(max(abs(double(c)-b))) < 1e-6)

%!test
%! % integer type
%! a = vpa(int32(6), 64);
%! b = vpa(6, 64);
%! assert (isequal (a, b))

%!test
%! % matrix of int
%! b = int32([pi 0; 6.25 1]);
%! c = vpa(b, 6);
%! assert (isequal (double(c), [3 0; 6 1]))

%!test
%! % can pass pi directly to vpa
%! a = vpa(sym(pi), 128);
%! b = vpa(pi, 128);
%! assert (isequal (a, b))

%!test
%! % if sym does sth special for e so should vpa
%! a = vpa(sym(exp(1)), 64);
%! b = vpa(exp(1), 64);
%! assert (isequal (a, b))

%!test
%! % can pass pi directly to vpa, even in array
%! a = vpa(sym([2 pi]), 128);
%! b = vpa([2 pi], 128);
%! assert (isequal (a, b))

%!test
%! % can pass i directly to vpa
%! a = vpa(sym(i));
%! b = vpa(i);

%!test
%! % 'i' and 'I' just make vars
%! a = vpa(sym(1i));
%! b = vpa('i');
%! c = vpa('I');
%! assert (~isequal (a, b))
%! assert (~isequal (a, c))

%!test
%! % '1i' and '1j' strings
%! a = vpa(sym(1i));
%! b = vpa('1i');
%! c = vpa('1j');
%! assert (isequal (a, b))
%! assert (isequal (a, c))

%!test
%! % Issue #868, precision loss on '0.33j'
%! a = vpa('0.33j', 40);
%! b = vpa('0.33i', 40);
%! assert (double (abs (imag (a)*100/33) - 1) < 1e-39)
%! assert (isequal (a, b))

%!test
%! % inf/-inf do not become symbol('inf')
%! S = {'oo', '-oo', 'inf', 'Inf', '-inf', '+inf'};
%! for j = 1:length(S)
%!   a = vpa(S{j});
%!   b = vpa(sym(S{j}));
%!   assert (isequal (a, b))
%! end

%!test
%! a = vpa('2.3', 20);
%! s = strtrim(disp(a, 'flat'));
%! assert (strcmp (s, '2.3000000000000000000'))

%!test
%! % these should *not* be the same
%! a = vpa(2.3, 40);
%! b = vpa('2.3', 40);
%! sa = sympy (a);
%! sb = sympy (b);
%! assert (~isequal (a, b))
%! assert (abs(double(a - b)) > 1e-20)
%! assert (abs(double(a - b)) < 1e-15)
%! assert (~strcmp(sa, sb))

%!test
%! % these should *not* be the same
%! x = vpa('1/3', 32);
%! y = vpa(sym(1)/3, 32);
%! z = vpa(1/3, 32);
%! assert (isequal (x, y))
%! assert (~isequal (x, z))

%!test
%! % big integers
%! a = int64(12345678);
%! a = a*a;
%! b = vpa(a);
%! c = vpa('152415765279684');
%! assert (isequal (b, c))

%!test
%! % big integers (workaround poor num2str, works in 4.0?)
%! a = int64(1234567891);  a = a*a;
%! b = vpa(a);
%! c = vpa('1524157877488187881');
%! assert (isequal (b, c))

%!warning <dangerous> vpa ('sqrt(2.0)');

%!warning <dangerous>
%! if (pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%! a = vpa('2**0.5');
%! b = vpa(sqrt(sym(2)));
%! assert (isequal (a, b))
%! else
%! warning('dangerous')  % fake it until we drop 1.4
%! end

%!test
%! a = vpa('2.3e1');
%! b = vpa(' 2.3e+1 ');
%! assert (isequal (a, b))
%! a = vpa('21e-1');
%! b = vpa('2.1');
%! assert (isequal (a, b))

%!test
%! % Issue #859, operations on immutable matrices
%! x = vpa (sym ([1 2]));
%! % If vpa no longer makes an ImmutableDenseMatrix,
%! % may need to adjust or remove this test.
%! assert (~ isempty (strfind (sympy (x), 'Immutable')))
%! y = sin(x);
%! y2 = [sin(vpa(sym(1))) sin(vpa(sym(2)))];
%! assert (isequal (y, y2))

repo: octsympy


text:
function my_diary_flush()

tmp_status=get(0,'Diary');
if strcmp(tmp_status, 'on')
    diary off
    diary on
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function test_scalingfactor

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_scalingfactor ft_convert_units

% since the function to test is in a private directory, we explicitely have to cd into that directory
[ftver, ftpath] = ft_version;
cd(fullfile(ftpath, 'utilities', 'private'));

assert(ft_scalingfactor('m', 'mm') == 1000);
assert(ft_scalingfactor('mm', 'm') == 0.001);
assert(ft_scalingfactor('T/cm', 'fT/m') == 1e17);


text:
function [varargout]=subQuadCatmullClark(varargin)

% function [Fs,Vs,Cs,CV]=subQuadCatmullClark(F,V,n,fixBoundaryOpt)
% ------------------------------------------------------------------------
%
%
%
% Change log:
% 2022/04/28 Minor array (Vn) allocation improvement
% ------------------------------------------------------------------------

%% parse input

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        n=1;
        fixBoundaryOpt=0;
    case 3
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        fixBoundaryOpt=0;
    case 4
        F=varargin{1};
        V=varargin{2};
        n=varargin{3};
        fixBoundaryOpt=varargin{4};
end

C=(1:1:size(F,1))'; %Face colors or indices

%%

if n>0
    for qIter=1:1:n
        M=patchConnectivity(F,V,{'ev','ef','vf','ve','vv'});

        edgeVertexMat=M.edge.vertex;
        edgeFaceMat=M.edge.face;

        vertexFaceMat=M.vertex.face;
        vertexEdgeMat=M.vertex.edge;
        vertexVertexMat=M.vertex.vertex;

        numPoints = size(V,1);
        numEdges = size(edgeVertexMat,1);

        % Get indices of the three edges associated with each face
        A = sparse(edgeVertexMat(:,1),edgeVertexMat(:,2),(1:numEdges)+numPoints,numPoints,numPoints,numEdges);
        A = max(A,A'); %Copy symmetric

        %Indices for A matrix
        indA_12=F(:,1)+(F(:,2)-1)*numPoints;
        indA_23=F(:,2)+(F(:,3)-1)*numPoints;
        indA_34=F(:,3)+(F(:,4)-1)*numPoints;
        indA_41=F(:,4)+(F(:,1)-1)*numPoints;

        %Get indices for vertex array
        indV_12=full(A(indA_12));
        indV_23=full(A(indA_23));
        indV_34=full(A(indA_34));
        indV_41=full(A(indA_41));

        indV_mid=(1:1:size(F,1))'+numPoints+size(edgeVertexMat,1);

        %Create faces array
        Fs=[[F(:,1)  indV_12 indV_mid indV_41];...
            [indV_12 F(:,2)  indV_23 indV_mid];...
            [indV_mid indV_23 F(:,3)  indV_34 ];...
            [indV_41 indV_mid indV_34 F(:,4) ]];

        %Face centre points
        Vm=patchCentre(F,V);

        %Create vertex arrays

        Vne=(V(edgeVertexMat(:,1),:)+ V(edgeVertexMat(:,2),:))/2;

        if size(edgeFaceMat,2)==1
            logicBoundaryEdges=true(size(edgeFaceMat,1),1);
        else
            logicBoundaryEdges=edgeFaceMat(:,2)==0;
        end

        V_F_mean=zeros(size(V));
        logicValid=vertexFaceMat>0;
        X=nan(size(vertexFaceMat));
        for q=1:1:size(V,2)
            X(logicValid)=Vm(vertexFaceMat(logicValid),q);
            V_F_mean(:,q)=mean(X,2,'omitnan');
        end

        V_E_mean=zeros(size(V));
        logicValid=vertexEdgeMat>0;
        X=nan(size(vertexEdgeMat));
        for q=1:1:size(V,2)
            X(logicValid)=Vne(vertexEdgeMat(logicValid),q);
            V_E_mean(:,q)=mean(X,2,'omitnan');
        end

        N=sum(logicValid,2);
        Vv=(V_F_mean+2*V_E_mean+(N-3).*V)./N;

        Vn=zeros(numEdges,size(V,2));
        if nnz(logicBoundaryEdges)>0
            %Use normal mid-edge nodes for boundary edges
            Vn(logicBoundaryEdges,:)=Vne(logicBoundaryEdges,:);
            indBoundaryVertices=unique(edgeVertexMat(logicBoundaryEdges,:));
            if fixBoundaryOpt==1
                %Replace boundary nodes with original
                Vv(indBoundaryVertices,:)=V(indBoundaryVertices,:);
            else
                vvm=vertexVertexMat(indBoundaryVertices,:);
                logicCheck=~ismember(vvm,indBoundaryVertices);
                vvm(logicCheck)=0;
                vvm=sort(vvm,2,'descend');
                vvm=vvm(:,[1 2]);

                Vv(indBoundaryVertices,:)= 6/8*V(indBoundaryVertices,:) + 1/8*(V(vvm(:,1),:)+V(vvm(:,2),:));
            end

        end

        Vn(~logicBoundaryEdges,:)=(V(edgeVertexMat(~logicBoundaryEdges,1),:)+V(edgeVertexMat(~logicBoundaryEdges,2),:)+...
            Vm(edgeFaceMat(~logicBoundaryEdges,1),:) + Vm(edgeFaceMat(~logicBoundaryEdges,2),:) )/4;

        Vs = [Vv; Vn; Vm]; %Join point sets

        if qIter>1
            CV=[CV; 1*ones(size(Vn,1),1); 2*ones(size(Vm,1),1);];
        else
            CV=[zeros(size(V,1),1); 1*ones(size(Vn,1),1); 2*ones(size(Vm,1),1);];
        end

    %Override input for looping
    C=repmat(C,[size(Fs,1)/size(F,1),1]);
    F=Fs;
    V=Vs;
    end

else
    CV=zeros(size(V,1),1);
end

varargout{1}=F;
varargout{2}=V;
varargout{3}=C;
varargout{4}=CV;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


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

[lat,lon] = meshgrat(tmap,tmapleg);
gx = rex(:,1);
gy = rey(1,:)';

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

ren = interp2(rex,rey,re,lon,lat);

l = ren < -8.8; ren(l) = -8.8;
mi = min(min(ren));

l = isnan(ren);
ren(l) = mi-2;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[33.5 35 ],'MapLonLimit',[-117.3 -116])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',8);
tightmap
view([0 90])
% hl = lightangle(45,25);
hl = camlight
lighting phong
set(gca,'projection','perspective');

load worldlo
%h = displaym(POline); set(h(1),'color',[0.9 0.9 0.9],'Linewidth',1.7)
% h2 = displaym(PPpoint);
%   h = displaym(PPtext); trimcart(h);

plotm(faults(:,2), faults(:,1),'k','Linewidth',1.4);
% plotm(mainfault(:,2), mainfault(:,1),'m','Linewidth',3);

pl = plotm(b(:,2),b(:,1),'ok');
set(pl,'LineWidth',0.3,'MarkerSize',2,...
    'MarkerFaceColor','k','MarkerEdgeColor','k')
pl = plotm(main(:,2),main(:,1),'hw');
set(pl,'LineWidth',1,'MarkerSize',20,...
    'MarkerFaceColor','w','MarkerEdgeColor','k')

zdatam(handlem('allline'),10000) % keep line on surface
%zdatam(handlem('alltext'),10000) % keep line on surface
caxis([0.0 0.10])
j = jet;
%j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];

colormap(j); brighten(0.3);

axis off; set(gcf,'color','w')

setm(gca,'ffacecolor','w')
setm(gca,'fedgecolor','k','flinewidth',2);

setm(gca,'mlabellocation',0.5)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',0.5)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','k','Fontweight','bold','FontSize',14,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.7 0.15 0.01 0.3],'TickDir','out','Ycolor','k','Xcolor','k',...
    'Fontweight','bold','FontSize',14);
set(gcf,'Inverthardcopy','off');

author: CelsoReyes


text:
function imgOut = KuangNormalizedGamma(img, gamma_value)
%
%
%       img = KuangNormalizedGamma(img, gamma_value)
%
%
%       Image is clamped if its values are over [a,b]
%
%       Input:
%           -img: the input img to be clamped.
%           -gamma_value: a gamma value to be applied to img.
%
%       Output:
%           -imgOut: an image with normalized gamma encoding.
%
%     Copyright (C) 2015  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('gamma_value', 'var'))
    gamma_value = 1.0;
end

max_img = max(img(:));

imgOut = ((img / max_img).^gamma_value) * max_img;

end

author: banterle


text:
function zz=lpcpz2zz(pz)
%LPCPZ2ZZ LPC: Power spectrum roots to LPC poles ZZ=(PZ)
% pz are the roots of the power spectrum polynomial pp(cos(w))

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcpz2zz.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zs=sqrt(pz.^2-1);
zz=pz-sign(real(conj(pz).*zs)).*zs;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [dVVectNTW] = getNTWdvVect(dVVect, rVect, vVect)
%getNTWdvVect Summary of this function goes here
%   Detailed explanation goes here
    vVect = reshape(vVect, 3,1);
    dVVect = reshape(dVVect, 3,1);

    tHat = vVect/norm(vVect);
    wHat = cross(rVect,vVect)/norm(cross(rVect,vVect));
    nHat = cross(tHat,wHat)/norm(cross(tHat,wHat));
    ECI2TWNRotMat = [tHat,wHat,nHat];
    dVVectNTW = ECI2TWNRotMat \ dVVect;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [b,R,se]=drxlr_wfit(y,x,w,p)
%DREX subfunction
%Written by Issam El Naqa 2003-2005
%Extracted for generalized use 2005, AJH
%
% Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
% This file is part of the Dose Response Explorer System (DREES).
%
% DREES development has been led by:  Issam El Naqa, Aditya Apte, Gita Suneja, and Joseph O. Deasy.
%
% DREES has been financially supported by the US National Institutes of Health under multiple grants.
%
% DREES is distributed under the terms of the Lesser GNU Public License.
%
%     This version of DREES 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.
%
% DREES 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 DREES.  If not, see <http://www.gnu.org/licenses/>.

sw = sqrt(w);
[r c] = size(x);
yw = y .* sw;
xw = x .* sw(:,ones(1,c));
[Q,R]=qr(xw,0);
b = R\(Q'*yw);
RI = R\eye(p);
C = RI * RI';
se = sqrt(max(eps,diag(C)));
return

author: mvallieres

repo: radiomics


text:
function dfupdatelegend(dffig,reset)
%DFUPLDATELEGEND Update legend in dfittool window

%   $Revision: 1.1.6.5 $  $Date: 2004/01/24 09:36:01 $
%   Copyright 2001-2004 The MathWorks, Inc.

if nargin<2, reset=false; end

% If figure not passed in, find figure that contains this thing
while ~isequal(get(dffig,'parent'),0),
  dffig = get(dffig,'parent');
end

% Remember info about old legend, if any
ax = get(dffig,'CurrentAxes');
if isempty(get(ax,'Children'))
   legh = [];
else
   legh = legend('-find',ax);
end

% Try to put the legend where the curves are not likely to be.
% Position "0" is supposed to do this, but it is deprecated
% and can be slow.  Instead use a heurisitic.
ftype = dfgetset('ftype');
if ischar(ftype) && ~isempty(ftype) && ismember(ftype, {'pdf' 'survivor'})
   % The survivor function is decreasing.  Most pdf functions drop more
   % quickly to the right.  "top right" is probably good for these functions.
   legendpos = 'TR';
else
   % The remaining ones are increasing, so "top left" is probably good.
   legendpos = 'TL';
end
legend(ax, 'off');

% Maybe no legend has been requested
if isequal(dfgetset('showlegend'),'off')
   return
end

% Get data line handles and labels
hh = flipud(findobj(ax,'Type','line'));
hData = findobj(hh,'flat','Tag','dfdata');
n = length(hData);

textData = cell(n,1);
for j=1:length(hData)
   nm = '';
   ds = get(hData(j),'UserData');
   if ~isempty(ds) && ishandle(ds) && ~isempty(findprop(ds,'name'))
      nm = ds.name;
   end
   if isempty(nm)
      hData(j) = NaN;
   else
      textData{j} = nm;
   end
end
t = ~isnan(hData);
textData = textData(t);
hData = hData(t);
sortData = 1000*(1:length(hData));
if isempty(sortData)
   maxnum = 0;
else
   maxnum = max(sortData) + 1000;
end

% Indent bounds if there are two or more data set lines
if n>1
   pre = '  ';
else
   pre = '';
end

% Deal with confidence bounds, if any, around empirical cdf
n = length(hData);
textDataBounds = cell(n,1);
hDataBounds = repmat(NaN,n,1);
sortDataBounds = zeros(n,1);
for j=1:n
   ds = get(hData(j),'UserData');
   if ds.showbounds
      hbounds = ds.boundline;
      if ~isempty(hbounds) && ishandle(hbounds) ...
                           && ~isempty(get(hbounds,'YData'))
         textDataBounds{j} = [pre 'confidence bounds'];
         hDataBounds(j) = hbounds;
         sortDataBounds(j) = sortData(j) + .5;
      end
   end
end
if any(isnan(hDataBounds))
   t = isnan(hDataBounds);
   textDataBounds(t) = [];
   hDataBounds(t) = [];
   sortDataBounds(t) = [];
end

% Indent fits if there are two or more data lines
if (length(hData)>1)
   pre = '  ';
else
   pre = '';
end

% Get fit line handles and labels
hFit = findobj(hh,'flat','Tag','distfit');
sortFit = NaN*hFit;
n = length(hFit);
hFitConf = NaN*zeros(n,1);
textFit = cell(n,1);
nms = cell(n,1);
for j=1:length(hFit)
   try
      fit = get(hFit(j),'UserData');
      nm = fit.name;
   catch
      nm = '';
   end
   if isempty(nm)
      hFit(j) = NaN;
   else
      nms{j} = nm;
      textFit{j} = [pre nm];

      % Find the dataset for this fit
      ds = fit.dshandle;
      sortFitj = maxnum + j;
      for k=1:length(hData)
         if isequal(ds.name,textData{k})
            sortFitj = sortData(k) + j;
            break;
         end
      end
      sortFit(j) = sortFitj;

      % Look for bounds
      b = get(fit,'boundline');
      if ~isempty(b)
         hFitConf(j) = b(1);
      end
   end
end
t = ~isnan(hFit);
nms = nms(t);
textFit = textFit(t);
hFit = hFit(t);
sortFit = sortFit(t);
hFitConf = hFitConf(t);

% Indent bounds if there are two or more fits
if (length(hFit)>1)
   pre = [pre '  '];
end

% Get confidence bound line handles and labels
n = length(hFitConf);
textFitBounds = cell(n,1);
sortFitBounds = zeros(size(hFitConf));
for j=1:length(hFitConf)
   if ~isnan(hFitConf(j)) && ishandle(hFitConf(j)) ...
                          && ~isempty(get(hFitConf(j),'XData'))
      textFitBounds{j} = sprintf('%sconfidence bounds (%s)',pre,nms{j});
      sortFitBounds(j) = sortFit(j) + 0.5;
   else
      hFitConf(j) = NaN;
   end
end
t = ~isnan(hFitConf);
textFitBounds = textFitBounds(t);
hFitConf = hFitConf(t);
sortFitBounds = sortFitBounds(t);

% Combine everything together for the legend
h = [hData(:); hDataBounds(:); hFit(:); hFitConf(:)];
c = [textData; textDataBounds; textFit; textFitBounds];
s = [sortData(:); sortDataBounds(:); sortFit(:); sortFitBounds(:)];

% Sort so related things are together
[s,j] = sort(s);
c = c(j);
h = h(j);

% Create the legend
if (length(h)>0)
   ws = warning;
   lw = lastwarn;
   warning('off');
   oldu = get(ax,'Units');
   set(ax,'Units','normalized');
   try
      havecoords = isnumeric(legendpos) && length(legendpos)>1;
      if havecoords
          temppos = 0;
      else
          temppos = legendpos;
      end
      [legh,objh] = legend(ax,h,c,temppos);
      if havecoords
         set(legh,'Units','points');
         newpos = get(legh,'Position');
         newpos(1:2) = legendpos(1:2);
         set(legh,'Position',newpos);
      end
   catch
      objh = [];
   end
%   set(ax,'Units',oldu);
   warning(ws);
   lastwarn(lw);

   % Avoid treating ds/fit names as TeX strings
   objh = findobj(objh,'flat','Type','text');
   set(objh,'Interpreter','none','Hittest','off');
end

% Set a resize function that will handle legend and layout
set(dffig,'resizefcn','dfittool(''adjustlayout'');');

author: zouchuhang

repo: LayoutNet