text:
function [yAvg, tHRF, nTrials] = hmrG_SubjAvg_Nirs(yAvgSubjs, tHRFSubjs, nTrialsSubjs)
% SYNTAX:
% [yAvg, tHRF, nTrials] = hmrG_SubjAvg_Nirs(yAvgSubjs, tHRFSubjs, nTrialsSubjs)
%
% UI NAME:
% Subj_Average
%
% DESCRIPTION:
% Calculate avearge HRF of all subjects in a group .
%
% INPUTS:
% yAvgSubjs:
% tHRFSubjs:
% nTrialsSubjs:
%
% OUTPUTS:
% yavg: the averaged results
% tHRF: the time vector
% nTrials: the number of trials averaged for each condition across all
%          subjects
%
% USAGE OPTIONS:
% Subj_Average_on_Concentration_Data: [dcAvg, tHRF, nTrials] = hmrG_SubjAvg_Nirs(dcAvgSubjs, tHRFSubjs, nTrialsSubjs)
% Subj_Average_on_Delta_OD_Data:      [dodAvg, tHRF, nTrials] = hmrG_SubjAvg_Nirs(dodAvgSubjs, tHRFSubjs, nTrialsSubjs)
%
%

yAvg       = [];
tHRF       = [];
nTrials    = [];

subjCh = [];
nStim = 0;
grp1=[];
nSubj = length(yAvgSubjs);

for iSubj = 1:nSubj

    if isempty(yAvgSubjs{iSubj}) || isempty(tHRFSubjs{iSubj}) || isempty(nTrialsSubjs{iSubj})
        continue;
    end

    yAvg      = yAvgSubjs{iSubj};
    tHRF      = tHRFSubjs{iSubj};
    nTrials   = nTrialsSubjs{iSubj};

    nCond = size(nTrials,2);

    if ndims(yAvg) == (4-(nCond<2))

        if iSubj==1
            grp1 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
        end

        nCh  = size(yAvg,3);
        if isempty(subjCh)
            subjCh = zeros(nCh, nCond);
        end

        for iC = 1:nCond
            if nTrials(:,iC)==0
                continue;
            end

            if iSubj==1 | iC>nStim
                for iCh = 1:size(yAvg,3)
                    for iHb=1:3
                        grp1(:,iHb,iCh,iC) = interp1(tHRF,yAvg(:,iHb,iCh,iC),tHRF(:));
                    end
                end
                nStim = iC;
            else
                for iCh = 1:size(yAvg,3)
                    for iHb=1:3
                        grp1(:,iHb,iCh,iC) = grp1(:,iHb,iCh,iC) + interp1(tHRF,yAvg(:,iHb,iCh,iC),tHRF(:));
                    end
                end
            end
            subjCh(:,iC) = subjCh(:,iC) + 1; %#ok<*AGROW>
        end

        yAvg = [];
        if ~isempty(grp1)
            for iC = 1:size(grp1,4)
                for iCh = 1:size(grp1,3)
                    yAvg(:,1,iCh,iC) = grp1(:,1,iCh,iC) / subjCh(iCh,iC);
                    yAvg(:,2,iCh,iC) = grp1(:,2,iCh,iC) / subjCh(iCh,iC);
                    yAvg(:,3,iCh,iC) = grp1(:,3,iCh,iC) / subjCh(iCh,iC);
                end
            end
        end

    elseif ndims(yAvg) == (3-(nCond<2))

        if iSubj==1
            grp1 = zeros(size(yAvg,1),size(yAvg,2),nCond);
        end

        nCh  = size(yAvg,2);
        if isempty(subjCh)
            subjCh = zeros(nCh, nCond);
        end

        for iC = 1:nCond
            if nTrials(:,iC)==0
                continue;
            end

            for iWl = 1:2
                if iSubj==1 | iC>nStim
                    for iCh = 1:size(yAvg,2)
                        grp1(:,iCh,iC) = interp1(tHRF,yAvg(:,iCh,iC),tHRF(:));
                    end
                    nStim = iC;
                else
                    for iCh = 1:size(yAvg,3)
                        grp1(:,iCh,iC) = grp1(:,iCh,iC) + interp1(tHRF,yAvg(:,iCh,iC),tHRF(:));
                    end
                end
                subjCh(:,iC) = subjCh(:,iC) + 1;
            end
        end

        yAvg = [];
        if ~isempty(grp1)
            for iC = 1:size(grp1,3)
                for iCh = 1:size(grp1,2)
                    yAvg(:,iCh,iC) = grp1(:,iCh,iC) / subjCh(iCh,iC);
                    yAvg(:,iCh,iC) = grp1(:,iCh,iC) / subjCh(iCh,iC);
                    yAvg(:,iCh,iC) = grp1(:,iCh,iC) / subjCh(iCh,iC);
                end
            end
        end

    end
end

author: BUNPC


% track a target
function tracker = MOT_track(frame_idx, frame_image, frame_size, bboxes, tracker, opt, seq_name)

if tracker.state == opt.STATE_TRACKED || tracker.state == opt.STATE_ACTIVATED
    tracker.num_occluded = 0;
    tracker.num_tracked = tracker.num_tracked + 1;
    tracker = track_frame(tracker, frame_idx, frame_image, bboxes, [], seq_name, opt);

    % check if target exits from the camera view
    [~, ov] = calc_overlap(tracker.bboxes, numel(tracker.bboxes.fr), frame_size, 1);
    if ov < opt.exit_thre
        if opt.is_print
            fprintf('The target exits from the camera view.\n');
        end
        tracker.state = opt.STATE_STOP;
    end
end


text:
function model3d_new = transform_3d_rigid(model3d,R,T)

if(model_has_normals(model3d))
    [x,y,z,nx,ny,nz] = model_to_components(model3d);
else
    [x,y,z] = model_to_components(model3d);
end

xyz = [x y z];
xyz = xyz * R';
xyz  = xyz + repmat(T',length(x),1);

if(model_has_normals(model3d))
    nxyz = [nx ny nz] * R';
    model3d_new = components_to_model(xyz(:,1),xyz(:,2),xyz(:,3),nxyz(:,1),nxyz(:,2),nxyz(:,3));
else
    model3d_new = components_to_model(xyz(:,1),xyz(:,2),xyz(:,3),[],[],[]);
end

if(model_has_colors(model3d))
    model3d_new.r = model3d.r;
    model3d_new.g = model3d.g;
    model3d_new.b = model3d.b;
end

if(model_has_faces(model3d))
    model3d_new.face = model3d.face;
end

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function estimate = spm_cfg_dcm_est
% SPM Configuration file for DCM estimation
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin & Peter Zeidman
% $Id: spm_cfg_dcm_est.m 7479 2018-11-09 14:17:33Z peter $

% -------------------------------------------------------------------------
% dcmmat Select DCM_*.mat
% -------------------------------------------------------------------------
dcmmat         = cfg_files;
dcmmat.tag     = 'dcmmat';
dcmmat.name    = 'Select DCM_*.mat';
dcmmat.help    = {'Select DCM_*.mat files.'};
dcmmat.filter  = 'mat';
dcmmat.ufilter = '^DCM_.*\.mat$';
dcmmat.num     = [1 Inf];

% -------------------------------------------------------------------------
% dcmmat Select GCM_.*.mat
% -------------------------------------------------------------------------
gcmmat         = cfg_files;
gcmmat.tag     = 'gcmmat';
gcmmat.name    = 'Select GCM_*.mat';
gcmmat.help    = {'Select GCM_*.mat files.'};
gcmmat.filter  = 'mat';
gcmmat.ufilter = '^GCM_.*\.mat$';
gcmmat.num     = [1 1];

%--------------------------------------------------------------------------
% dir Directory
%--------------------------------------------------------------------------
dir         = cfg_files;
dir.tag     = 'dir';
dir.name    = 'Directory';
dir.help    = {'Select the directory where the output will be written.'};
dir.filter  = 'dir';
dir.ufilter = '.*';
dir.num     = [1 1];

% -------------------------------------------------------------------------
% name Model name
%--------------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {['Specify a name for the group DCM file. The prefix GCM_' ...
                'and suffix .mat are automatically added']};
name.strtype = 's';
name.num     = [0 Inf];

% -------------------------------------------------------------------------
% subj Create single subject
%--------------------------------------------------------------------------
subj      = cfg_branch;
subj.tag  = 'subj';
subj.name = 'Subject';
subj.val  = {dcmmat};
subj.help = {'Subject with one or more models.'};

% -------------------------------------------------------------------------
% multiple_models Create set of subjects
%--------------------------------------------------------------------------
subjects        = cfg_repeat;
subjects.tag    = 'subjects';
subjects.name   = 'Per subject';
subjects.values = {subj};
subjects.help   = {'Create the subjects and select the models for each'};
subjects.num    = [1 Inf];

% -------------------------------------------------------------------------
% subj Create single model
%--------------------------------------------------------------------------
model      = cfg_branch;
model.tag  = 'model';
model.name = 'Model';
model.val  = {dcmmat};
model.help = {'Corresponding model for each subject'};

% -------------------------------------------------------------------------
% subjects Create set of models
%--------------------------------------------------------------------------
models        = cfg_repeat;
models.tag    = 'models';
models.name   = 'Per model';
models.values = {model};
models.help   = {'Select DCM.mat files per model'};
models.num    = [1 Inf];

% -------------------------------------------------------------------------
% output_single Output one .mat file for the group
%--------------------------------------------------------------------------
output_single         = cfg_branch;
output_single.tag     = 'single';
output_single.name    = 'Create group GCM_*.mat file';
output_single.val     = { dir name };
output_single.help    = {['Creates a single group-level DCM file ' ...
                          'containing a subjects x models cell array.']};

% -------------------------------------------------------------------------
% output_overwrite_gcm Output a GCM file with existing name
%--------------------------------------------------------------------------
output_overwrite_gcm         = cfg_branch;
output_overwrite_gcm.tag     = 'overwrite_gcm';
output_overwrite_gcm.name    = 'Overwrite existing GCM/DCM files';
output_overwrite_gcm.val     = {};
output_overwrite_gcm.help    = {['Overwrites existing group-level DCM file ' ...
                                 'with estimated models.']};

% -------------------------------------------------------------------------
% output_separate Output one .mat file per model
%--------------------------------------------------------------------------
output_separate         = cfg_branch;
output_separate.tag     = 'separate';
output_separate.name    = 'Only save individual DCM files';
output_separate.val     = {};
output_separate.help    = {'Updated existing individual DCM.mat files'};

% -------------------------------------------------------------------------
% output Choice of how many DCM.mat files to output
%--------------------------------------------------------------------------
output         = cfg_choice;
output.tag     = 'output';
output.name    = 'Output';
output.values  = { output_single output_overwrite_gcm output_separate };
output.val     = { output_single };
output.help    = {'How to store the estimated DCMs. The options are: ' ...
                  ['1. Create group GCM_*.mat file - this will create a single '...
                   'file containing a cell array, with one row per subject ' ...
                   'and one column per DCM, containing the filenames of ' ...
                   'the DCMs. The original DCM files will be overwritten ' ...
                   'with the outcome of the estimation. ' ...
                   '(Alternatively, if the input is a GCM file containing ' ...
                   'DCM structures rather than filenames, then the output ' ...
                   'will also be an array containing DCM structures.)'] ...
                  ['2. Overwrite existing GCM/DCM files - as above but with ' ...
                   'with an existing GCM filename.'] ...
                  ['3. Only save individual DCM files - does not change the ' ...
                   'GCM file and simply overwrites each subject''s DCM ' ...
                   'with estimated values.']};

% -------------------------------------------------------------------------
% way Choice of ways to select DCMs (nested models)
%--------------------------------------------------------------------------
dcms        = cfg_choice;
dcms.tag    = 'dcms';
dcms.name   = 'Select DCMs';
dcms.values = {models subjects gcmmat};
dcms.val    = {gcmmat};
dcms.help   = {['Select one DCM per subject, multiple DCMs per subject ' ...
                  'or an existing group DCM file.'] ...
                 ['If multiple DCMs are selected per subject, then ' ...
                  'the first DCM for each subject should a ''full'' ' ...
                  'model containing all connections of interest. Subsequent ' ...
                  '(nested) DCMs will have certain connections switched ' ...
                  'off.']};

% -------------------------------------------------------------------------
% est_type Estimation type
%--------------------------------------------------------------------------
est_type        = cfg_menu;
est_type.tag    = 'est_type';
est_type.name   = 'Estimation type';
est_type.labels = {'Full + BMR (default)',...
                   'Full + BMR PEB (more accurate but slower)',...
                   'Full (not recommended)',...
                   'None (collate only)'};
est_type.values = {1,2,3,4};
est_type.val    = {1};
est_type.help   = {['Full + BMR: Estimates the full (first) model for ' ...
                    'each subject then uses Bayesian Model Reduction (BMR) '...
                    'to rapidly infer the evidence / parameters for any ' ...
                    'subsequent nested models.'] ...
                   ['Full + BMR PEB: Iteratively estimates the full (first) '...
                    'model for each subject, then sets the priors on each' ...
                    'each parameter to the group mean (from a PEB model)' ...
                    'then re-estimates. This improves estimation '...
                    'by overcoming local optima, but takes longer.' ] ...
                   ['Full: Estimates all models individually. Provided for '...
                    'backward compatibility.']...
                   ['None: Creates a group level DCM file without ' ...
                    'performing estimation']};

% -------------------------------------------------------------------------
% analysis Analysis
% -------------------------------------------------------------------------
fmri_analysis         = cfg_menu;
fmri_analysis.tag     = 'analysis';
fmri_analysis.name    = 'Analysis';
fmri_analysis.labels  = {'default (time series)','cross-spectral densities'};
fmri_analysis.values  = {'time','csd'};
fmri_analysis.val     = {'time'};
fmri_analysis.help    = {['Whether to analyse in the time domain (for task-' ...
                     'based studies or stochastic DCM) or in the frequency '...
                     'domain (for resting state analysis with DCM for CSD']};

fmri         = cfg_branch;
fmri.tag     = 'fmri';
fmri.name    = 'MRI specific options';
fmri.val     = {fmri_analysis};
fmri.help    = {'MRI specific options'};

% -------------------------------------------------------------------------
% estimate Estimate
% -------------------------------------------------------------------------
estimate      = cfg_exbranch;
estimate.tag  = 'estimate';
estimate.name = 'DCM estimation';
estimate.val  = { dcms output est_type fmri };
estimate.help = {['Estimate the parameters and free energy (log model ' ...
                  'evidence) of first level DCMs for fMRI. Models ' ...
                  'are assembled into a Subjects x Models array and ' ...
                  'saved in group GCM_*.mat file']};
estimate.prog = @spm_run_dcm_est;
estimate.vout = @vout_dcm;

% -------------------------------------------------------------------------
% fmri Dynamic Causal Model for fMRI
% -------------------------------------------------------------------------
est         = cfg_choice;
est.tag     = 'est';
est.name    = 'DCM estimation';
est.help    = {'Estimation of Dynamic Causal Models.'};
est.values  = { estimate };

%==========================================================================
function out = spm_run_dcm_est(job)
%==========================================================================

dcms = job.dcms;

% Get selected estimation option
EST_FULL_BMR     = 1;
EST_FULL_BMR_PEB = 2;
EST_FULL         = 3;
EST_NONE         = 4;

est_type = job.est_type;

% Get selected input option
INPUT_DCM_BY_MODEL   = 1;
INPUT_DCM_BY_SUBJECT = 2;
INPUT_GCM            = 3;

if isfield(dcms,'model')
    input_type = INPUT_DCM_BY_MODEL;
elseif isfield(dcms,'subj')
    input_type = INPUT_DCM_BY_SUBJECT;
elseif isfield(dcms,'gcmmat')
    input_type = INPUT_GCM;
else
    error('Unknown input type');
end

% Get selected output option
OUTPUT_GCM_NEW       = 1;
OUTPUT_GCM_OVERWRITE = 2;
OUTPUT_DCM           = 3;

if isfield(job.output,'single')
    output_type = OUTPUT_GCM_NEW;
elseif isfield(job.output,'overwrite_gcm')
    output_type = OUTPUT_GCM_OVERWRITE;
elseif isfield(job.output,'separate')
    output_type = OUTPUT_DCM;
else
    error('Unknown output type');
end

if (output_type == OUTPUT_GCM_OVERWRITE && input_type ~= INPUT_GCM)
    error('To overwrite an existing GCM file, the input must be a GCM');
end

% Build subjects x models filename matrix
switch input_type
    case INPUT_DCM_BY_MODEL
        ns = length(dcms.model(1).dcmmat);
        nm = length(dcms.model);
        P  = cell(ns,nm);

        for m = 1:nm
            if length(dcms.model(m).dcmmat) ~= ns
                error(['Please ensure all models have the same number of ' ...
                       'subjects']);
            end

            P(:,m) = dcms.model(m).dcmmat;
        end

    case INPUT_DCM_BY_SUBJECT
        ns  = length(dcms.subj);
        nm  = length(dcms.subj(1).dcmmat);
        P = cell(ns,nm);

        for s = 1:ns
            if length(dcms.subj(s).dcmmat) ~= nm
                error(['Please ensure all subjects have the same number of ' ...
                       'models']);
            end

            P(s,:) = dcms.subj(s).dcmmat';
        end

    case INPUT_GCM
        GCM = load(dcms.gcmmat{1});
        GCM = GCM.GCM;
        ns = size(GCM,1);
        nm = size(GCM,2);

        if ischar(GCM{1})
            P = GCM;
        else
            P = '';
        end
end

% Validate
if output_type == OUTPUT_DCM && isempty(P)
    error(['Cannot save individual DCM files when the input is a GCM ' ...
           'array of DCM structures']);
end

% Load all models into memory
if ~isempty(P)
    GCM = spm_dcm_load(P);
end

% Set timeseries or CSD estimation (fMRI)
for s = 1:ns
    for m = 1:nm
        if strcmpi(job.fmri.analysis,'CSD')
            GCM{s,m}.options.analysis   = 'CSD';
            GCM{s,m}.options.induced    = 1;
            GCM{s,m}.options.stochastic = 0;
        else
            if isfield(GCM{s,m},'options') && isfield(GCM{s,m},'analysis')
                GCM{s,m}.options = rmfield(GCM{s,m}.options,'analysis');
            end
            GCM{s,m}.options.induced = 0;
        end
    end
end

% Check that models 2-N are nested forms of the full model
if nm > 1 && (est_type == EST_FULL_BMR || est_type == EST_FULL_BMR_PEB)
    idx1 = spm_find_pC(GCM{1});
    for m = 2:nm
        idx = spm_find_pC(GCM{1,m});
        if any(setdiff(idx,idx1))
            error(['Model %d is not a nested model of model 1. This is ' ...
                  'required for Bayesian Model Reduction (BMR). Please ' ...
                  'introduce a full model as Model 1, or switch estimation ' ...
                  'type to ''full'''],m);
        end
    end
end

% Estimate models if requested
switch est_type
    case EST_FULL_BMR
        GCM(:,1) = spm_dcm_fit(GCM(:,1));

        if nm > 1
            GCM = spm_dcm_bmr(GCM);
        end
    case EST_FULL_BMR_PEB
        GCM = spm_dcm_peb_fit(GCM);
    case EST_FULL
        GCM = spm_dcm_fit(GCM);
    case EST_NONE
        % Do nothing
end

% Save individual DCM .mat files if requested
if ~isempty(P) && (est_type ~= EST_NONE)
    for s = 1:ns
        for m = 1:nm
            DCM = GCM{s,m};
            F   = DCM.F;
            Ep  = DCM.Ep;
            Cp  = DCM.Cp;
            save(P{s,m}, 'DCM', 'F', 'Ep', 'Cp', ...
                spm_get_defaults('mat.format'));
        end
    end
end

% If filenames were provided, set the GCM to contain the filenames
if ~isempty(P)
    GCM = P; %#ok<NASGU>
end

% Save GCM
if output_type == OUTPUT_GCM_NEW
    % Create single GCM mat file
    dir  = job.output.single.dir{1};
    name = ['GCM_' job.output.single.name '.mat'];
    filename = fullfile(dir,name);
    save(filename,'GCM', spm_get_defaults('mat.format'));

    out.gcmmat = {filename};
elseif output_type == OUTPUT_GCM_OVERWRITE
    % Update existing gcm file
    filename = dcms.gcmmat{1};
    save(filename,'GCM', spm_get_defaults('mat.format'));

    out.gcmmat = {filename};
end

if ~isempty(P)
    out.dcmmat = P(:,1);
end

%==========================================================================
function dep = vout_dcm(job)
%==========================================================================
dep(1)            = cfg_dep;
dep(1).sname      = 'DCM mat File(s) - full models';
dep(1).src_output = substruct('.','dcmmat');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});

if isfield(job.output,'single') || ...
        isfield(job.output,'overwrite_gcm')
    dep(2)            = cfg_dep;
    dep(2).sname      = 'GCM mat File(s)';
    dep(2).src_output = substruct('.','gcmmat');
    dep(2).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
end

author: spm


text:
function PDt(option, up)
%PDT detects PPG pulse peaks in PPG signals.
%	            PDt(option, up)
%
%	Inputs:
%		option          the option which has led to this function being used
%       up              universal parameters structure
%
%	Outputs:
%       ...
%

fprintf('\n--- Detecting Pulse Peaks ');
log_int_respSig = 1;             % Has value 1 unless this is a final respiratory signal

for subj = up.paramSet.subj_list

    %% Skip if this processing has been done previously
    iden_resp_sig_file_ending
    savepath = [up.paths.data_save_folder, num2str(subj), ending];
    filecontents = whos('-file', savepath);
    var_names = extractfield(filecontents, 'name');
    temp = strfind(var_names, [option(1:3), up.paths.filenames.elim_vhf]); rel_els = find(~cellfun(@isempty,temp));
    rel_var_names = var_names(rel_els);
    for rel_var_name_no = 1 : length(rel_var_names)
        for current_opt_no = 1 : length(up.al.options.PDt)
            eval(['save_name = ''' option(1:3), up.paths.filenames.pulse_peaks, up.al.options.PDt{current_opt_no}, rel_var_names{rel_var_name_no}(4:end) ''';']);
            exist_log = check_exists(savepath, save_name);
            if exist_log
                continue
            end

            %% Load relevant data
            loadpath = [up.paths.data_save_folder, num2str(subj), up.paths.filenames.int_respSigs];
            rel_name = rel_var_names{rel_var_name_no};
            load(loadpath, rel_name);
            eval(['s = ' rel_name ';']);
            eval(['fs = ' rel_name '.fs;']);

            %% High-pass filter data for peak detection
            s_filt = elim_sub_cardiac(s, up);

            %% Detect Pulse Peaks in HPF'd version
            [peaks,onsets] = feval(up.al.options.PDt{current_opt_no}, s_filt, fs, up);

            eval([save_name '.fs = fs;']);
            % Identify Pulse Peaks in original version
            eval([save_name '.p.v = s.v(peaks);']);
            eval([save_name '.p.t = s.t(peaks);']);
            eval([save_name '.tr.v = s.v(onsets);']);
            eval([save_name '.tr.t = s.t(onsets);']);
            % Identify start and end times of raw data (for resampling)
            eval([save_name '.timings.t_start = s.t(1);']);
            eval([save_name '.timings.t_end = s.t(end);']);

            %% Save processed data
            save_or_append_data

            clear peaks onsets s_filt s fs temp

        end
    end

end

end

function [peaks,onsets] = IMS(s_filt,fs,up)

[peaks,onsets,artifs] = adaptPulseSegment(s_filt.v,fs);

end

author: peterhcharlton

repo: RRest


text:
function val=sincFun(x)
%%SINCFUN Evaluate the sinc function. This is defined to be
%         sin(pi*x)/(pi*x). Note that some authors use the definition
%         sin(x)/x. This function accounts for the singularity at zero.
%
%INPUTS: x A matrix of the values at which to evaluate the sinc function.
%
%OUTPUTS: vals A matrix the same size as x where the sinc function has been
%              evaluted at each of the elements.
%
%The sinc function is described in [1]. It arises in many areas of signal
%processing.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/SincFunction.html
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    y=pi*x;
    sel=(y==0);
    val=sin(y)./y;
    val(sel)=1;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% This script estimates the vowel space ratio of a speaker with respect to
% his/her reference population.
%
%
% Description
% Based on the tracked formants F1 and F2 for the voiced regions of speech
% we compute the vowel space for each recorded subject individually [1].
% We define the vowel space as the frequency region covered by the triangle
% in the two dimensional frequency space spanned by F1 and F2 for the vowels
% /i/ (as in heed), /a/ (as in hod), and /u/ (as in who'd) following [2].
% These three vowels represent the vowels with the most extreme positions
% of the tongue and are therefore located in the extremes of this triangularly
% shaped two-dimensional frequency space [3]. The method is developed
% following the recommendations in [4]. Reference frequencies of american
% english vowels are found in [5].
%
%
% [ratio, centroids] = getVowelSpace(formant_data, gender, mode, scaling, plot_formants)
%
% Inputs
%  formant_dat      : [samples] [Nx2] input formant observations F1 and F2
%  gender           : [nominal] [1x1] female = 0; male = 1; child = 2
%  mode             : [code]    [string] 'triangle' (default); 'polygon'
%                               (add /ae/ (as in had) to the vowel space)
%  scaling          : [boolean] [1x1] bark scaling on/off (default);
%                               mel_scaling optional in code.
%  plot_formants    : [boolean] [1x1] visualize vowel space on/off
%                               (default)
%
% Outputs
%  ratio            : [1x1] vowel space ratio
%  centroids        : [1x12] vector containing prototype vowel locations
%
% References
% [1] Scherer, S., Morency, L.-P., Gratch, J., and Pestian, J., EDUCED VOWEL SPACE
% IS A ROBUST INDICATOR OF PSYCHOLOGICAL DISTRESS: A CROSS-CORPUS ANALYSIS,
% Proceedings of ICASSP 2015.
% [2] H.-M. Liu, F.-M. Tsao, and P. K. Kuhl, ?The effect of reduced vowel
% working space on speech intelligibility in mandarin-speaking young
% adults with cerebral palsy,? Journal of the Acoustical Society of America,
% vol. 117, no. 6, pp. 3879?3889, 2005.
% [3] B. Lindblom, ?Explaining phonetic variation: A sketch of the h&h
% theory,? Speech Production and Speech Modeling, pp. 403?439, 1990.
% [4] S. Sandoval, V. Berisha, R. L. Utianski, J. M. Liss, and A. Spanias,
% ?Automatic assessment of vowel space area,? Journal of the Acoustical
% Society of America, vol. 134, no. 5, pp. 477?483, 2013.
% [5] J. Hillenbrand, L. A. Getty, M. J. Clark, and K. Wheeler, ?Acoustic
% characteristics of american english vowels,? Journal of the Acoustical
% Society of America, vol. 97, no. 5, pp. 3099?3111, 1995.
%
% Copyright (c) 2014 University of Southern California, Institute for
% Creative Technologies
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Stefan Scherer scherer@ict.usc.edu
%
function [ratio, centroids] = getVowelSpace(formant_data, gender, mode, scaling, plot_formants)

%% initialization
mel_scaling = @(x)(2595.*log10(ones(size(x))+x./700));
bark_scaling = @(x)(13.*atan(0.00076.*x)+3.5.*atan((x./7500).^2));
min_observations = 1000;

switch nargin
    case 1
        gender = 1;
        warning('getVowelSpace.m : no gender specified. assumes male gender.');
        plot_formants = 0;
        scaling = false;
        mode = 'triangle';
    case 2
        if(~ismember(gender, [0,1,2]))
            error('getVowelSpace.m : wrong use of gender variable; permitted values: 1 .. male; 0 .. female; 2 .. child');
        end
        plot_formants = 0;
        scaling = false;
        mode = 'triangle';
    case 3
        if(~ismember(gender, [0,1,2]))
            error('getVowelSpace.m : wrong use of gender variable; permitted values: 1 .. male; 0 .. female; 2 .. child');
        end
        if(~ismember(mode, {'triangle','polygon'}))
            error('getVowelSpace.m : wrong use of gender variable; permitted values: 1 .. male; 0 .. female; 2 .. child');
        end
        scaling = false;
        plot_formants = 0;
    case 4
        if(~ismember(gender, [0,1,2]))
            error('getVowelSpace.m : wrong use of gender variable; permitted values: 1 .. male; 0 .. female; 2 .. child');
        end
        plot_formants = 0;
    case 5
        if(~ismember(gender, [0,1,2]))
            error('getVowelSpace.m : wrong use of gender variable; permitted values: 1 .. male; 0 .. female; 2 .. child');
        end
    otherwise
        error('getVowelSpace.m : wrong use of getVowelSpace.');
end

if(size(formant_data,2) < 2)
    error('getVowelSpace.m : wrong use of getVowelSpace. you must provide at least the first two formants!');
end

%% Formants for vowels
% list of formants; based on [4]
% this list is quite flexible.. you can adjust that to your dialect/language if you like!
men_f1 = [342, 427, 476, 580, 588, 768, 652, 497, 469, 378, 623, 474];
men_f2 = [2322, 2034, 2089, 1799, 1952, 1333, 997, 910, 1122, 997, 1200, 1379];

fem_f1 = [437, 483, 536, 731, 669, 936, 781, 555, 519, 459, 753, 523];
fem_f2 = [2761, 2365, 2530, 2058, 2349, 1551, 1136, 1035, 1225, 1105, 1426, 1588];

child_f1 = [452, 511, 564, 749, 717, 1002, 803, 597, 568, 494, 794, 586];
child_f2 = [3081, 2552, 2656, 2267, 2501, 1688, 1210, 1137, 1490, 1345, 1546, 1719];

if strcmp(mode, 'polygon')
    target_vowels = [1, 5, 6, 10];
else
    target_vowels = [1, 6, 10];
end

%% initialize data for processing
if(gender == 1)
    init_formants = [men_f1; men_f2]';
    formant_data = formant_data(formant_data(:,1) < 900 & formant_data(:,2) > 800 & formant_data(:,1) > 250 & formant_data(:,2) < 2450, :); % remove some outliers of formants
elseif(gender == 0)
    init_formants = [fem_f1; fem_f2]';
    formant_data = formant_data(formant_data(:,1) < 1000 & formant_data(:,2) > 1000 & formant_data(:,1) > 350 & formant_data(:,2) < 2800, :); % remove some outliers of formants
elseif(gender == 2)
    init_formants = [child_f1; child_f2]';
    formant_data = formant_data(formant_data(:,1) < 1000 & formant_data(:,2) > 1000 & formant_data(:,1) > 350 & formant_data(:,2) < 2800, :); % remove some outliers of formants
end

if(size(formant_data,1) < min_observations)
    warning('getVowelSpace.m : More formant observations are required for this method to work robustly.');
    ratio = 0;
    centroids = [];
    return
end

%% bark scaling
if scaling
    formant_data = bark_scaling(formant_data);
    init_formants = bark_scaling(init_formants);
    % optional:
    % formant_data = mel_scaling(formant_data);
    % init_formants = mel_scaling(init_formants);
end

%% vector quantization
warning off;
[~,centroids] = kmeans(formant_data(:,1:2),[],'start',init_formants, 'emptyaction', 'singleton');
warning on;

%% temp data for clustering
% distance measure can be changed!
pairwise_distances = pdist2(init_formants, centroids, 'seuclidean');
[~, idx] = min(pairwise_distances(target_vowels, :), [], 2);

%% triangle ratio
convHullPoints=centroids(idx,:);
vowelSpacePoints = init_formants(target_vowels,:);
[convOrder , cur_space] = convhull(convHullPoints);
[spaceOrder , vowel_space] = convhull(vowelSpacePoints);
ratio = cur_space./vowel_space;

%% plot vowel space
if(plot_formants)
    figure
    scatter(formant_data(:,1),formant_data(:,2), '.', 'LineWidth', 0.05);
    hold on;
    plot(centroids(:,1), centroids(:,2), 'ro')
    plot(centroids(idx,1), centroids(idx,2), 'go')
    xlabel('Formant 1');
    ylabel('Formant 2');
    plot(init_formants(:,1), init_formants(:,2), 'kx')
    plot(init_formants(target_vowels,1), init_formants(target_vowels,2), 'mx')
    plot(convHullPoints(convOrder,1),convHullPoints(convOrder,2),'g-')
    plot(vowelSpacePoints(spaceOrder,1),vowelSpacePoints(spaceOrder,2),'m-')
    set(gca,'XDir','reverse');
    set(gca,'YDir','reverse');
    set(gca,'XAxisLocation','top','YAxisLocation','right');
    title(sprintf('vowel space ratio %2.3f (#obs = %d)', ratio, size(formant_data, 1)));
    drawnow;
    hold off;
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function SaveParamsAsNIfTI(paramsfile, roifile, targetfile, outputpref)
%
% function SaveParamsAsNIfTI(paramsfile, roifile, targetfile, outputpref)
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

% load the fitting results
fprintf('loading the fitted parameters from : %s\n', paramsfile);
load(paramsfile);

% load the roi file
fprintf('loading the roi from : %s\n', roifile);
load(roifile);

% load the target volume
fprintf('loading the target volume : %s\n', targetfile);
target = nifti(targetfile);
xsize = target.dat.dim(1);
ysize = target.dat.dim(2);
if length(target.dat.dim) == 2
    zsize = 1;
else
    zsize = target.dat.dim(3);
end
total = xsize*ysize*zsize;

% determine the volumes to be saved
% the number of fitted parameters minus the two fibre orientation
% parameters.
idxOfFitted = find(model.GD.fixed(1:end-2)==0);
noOfVols = length(idxOfFitted);
vols = zeros(total,noOfVols);
% the volume for the objective function values
fobj_ml_vol = zeros(total,1);
% the volume for the error code
error_code_vol = zeros(total,1);
% some special cases
if (strfind(model.name, 'Watson'))
    idxOfKappa = find(ismember(model.paramsStr, 'kappa')==1);
    odi_vol = zeros(total,1);
end

% the volumes for the fibre orientations
fibredirs_x_vol = zeros(total,1);
fibredirs_y_vol = zeros(total,1);
fibredirs_z_vol = zeros(total,1);
% compute the fibre orientations from the estimated theta and phi.
fibredirs = GetFibreOrientation(model.name, mlps);

% determine the volumes with MCMC fitting to be saved
if (model.noOfStages==3)
    idxOfFittedMCMC = find(model.MCMC.fixed(1:end-2)==0);
    noOfVolsMCMC = length(idxOfFittedMCMC);
    volsMCMC = zeros(total,noOfVolsMCMC);

    if (strfind(model.name, 'Watson'))
        idxOfKappa = find(ismember(model.paramsStr, 'kappa')==1);
	     if (model.MCMC.fixed(idxOfKappa)==0)
            odi_volMCMC = zeros(total,1);
        end
    end
end

% convert to volumetric maps
fprintf('converting the fitted parameters into volumetric maps ...\n');
for i=1:size(mlps,1)
    % compute the index to 3D
    volume_index = (idx(i,3)-1)*ysize*xsize + (idx(i,2)-1)*xsize + idx(i,1);

    % fitted parameters other than the fiber orientations
    for j=1:length(idxOfFitted)
        vols(volume_index,j) = mlps(i, idxOfFitted(j));
    end

    % objective function values
    fobj_ml_vol(volume_index) = fobj_ml(i);

    % error codes
    error_code_vol(volume_index) = error_code(i);

    % fiber orientations
    fibredirs_x_vol(volume_index) = fibredirs(1,i);
    fibredirs_y_vol(volume_index) = fibredirs(2,i);
    fibredirs_z_vol(volume_index) = fibredirs(3,i);

    % special cases
    if (strfind(model.name, 'Watson'))
        odi_vol(volume_index) = atan2(1, mlps(i,idxOfKappa)*10)*2/pi;
    end

	 % MCMC fitted parameters
    if (model.noOfStages==3)
        for j=1:length(idxOfFittedMCMC)
            volsMCMC(volume_index,j) = mean(squeeze(mcmcps(i,:,j)));
        end
        if (strfind(model.name, 'Watson'))
            % Warning: Here we hard-coded the index to kappa!!!
            odi_volMCMC(volume_index) = atan2(1, mean(squeeze(mcmcps(i,:,3)))*10)*2/pi;
        end
    end

end

% save as NIfTI
fprintf('Saving the volumetric maps of the fitted parameters ...\n');

niftiSpecs.dim = [xsize ysize zsize];
niftiSpecs.mat = target.mat;
niftiSpecs.mat_intent = target.mat_intent;
niftiSpecs.mat0 = target.mat0;
niftiSpecs.mat0_intent = target.mat0_intent;

% the fitted parameters other than the fiber orientations
for i=1:length(idxOfFitted)
    output = [outputpref '_' cell2mat(model.paramsStr(idxOfFitted(i))) '.nii'];
    SaveAsNIfTI(reshape(squeeze(vols(:,i)), [xsize ysize zsize]), niftiSpecs, output);
end

% the special cases
if (strfind(model.name, 'Watson'))
    output = [outputpref '_' 'odi.nii'];
    SaveAsNIfTI(reshape(odi_vol, [xsize ysize zsize]), niftiSpecs, output);
end

% the objective function values
output = [outputpref '_' 'fmin.nii'];
SaveAsNIfTI(reshape(fobj_ml_vol, [xsize ysize zsize]), niftiSpecs, output);

% the error codes
output = [outputpref '_' 'error_code.nii'];
SaveAsNIfTI(reshape(error_code_vol, [xsize ysize zsize]), niftiSpecs, output);

% the fibre orientations
output = [outputpref '_' 'fibredirs_xvec.nii'];
SaveAsNIfTI(reshape(fibredirs_x_vol, [xsize ysize zsize]), niftiSpecs, output);
output = [outputpref '_' 'fibredirs_yvec.nii'];
SaveAsNIfTI(reshape(fibredirs_y_vol, [xsize ysize zsize]), niftiSpecs, output);
output = [outputpref '_' 'fibredirs_zvec.nii'];
SaveAsNIfTI(reshape(fibredirs_z_vol, [xsize ysize zsize]), niftiSpecs, output);

% the MCMC fitted parameters
if (model.noOfStages==3)
    for i=1:length(idxOfFittedMCMC)
        output = [outputpref '_' cell2mat(model.paramsStr(idxOfFittedMCMC(i))) '_MCMC.nii'];
        SaveAsNIfTI(reshape(squeeze(volsMCMC(:,i)), [xsize ysize zsize]), niftiSpecs, output);
    end
    if (strfind(model.name, 'Watson'))
        output = [outputpref '_' 'odi_MCMC.nii'];
        SaveAsNIfTI(reshape(odi_volMCMC, [xsize ysize zsize]), niftiSpecs, output);
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function fX=SchurMatFunEval(f,X)
%%SCHURMATFUNEVAL Evaluate the matrix function f(X), where the function f
%         is a function that has a Taylor series representation using the
%         Schur-Parlett method. In this approach, f(x) only needs to be
%         evaluated on _scalar_ quantities. Thus, this function can be used
%         to turn scalar functions into matrix functions. Matrix functions
%         with a Taylor series expansions satisfy the identities
%         X*f(X)=f(X)*X and f(inv(B)*X*B)=B*f(X)*inv(B) and include
%         functions such as the matrix exponential and square root.
%
%INPUTS: f A function handle that takes a scalar input and returns a scalar
%          output.
%        X The real or complex nXn matrix that should be evaluated via the
%          function f.
%
%OUTPUTS: fX The value f(X), where the function f applied to scalar terms
%            is applied to the matrix X. This is not the same as applying f
%            to every scalar element in X.
%
%This function implements Algorithm 9.1.1 of Chapter 9.1.4 of [1].
%
%EXAMPLE:
%Here, we show that the result agrees with the mpow function to take a
%matrix square root.
% X=randn(6,6);
% f=@(x)sqrt(x);
% fX=SchurMatFunEval(f,X);
% diffVal=fX-mpower(X,1/2)
%One will see that the entries in diffVal are all on the order of 1e-14 or
%less, indicating that the two are probably equal within finite precision
%bounds.
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%January 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=length(X);

%We must use the 'complex' option or else T is not guaranteed to be
%upper triangular.
[Q,T]=schur(X,'complex');

%The diagonal terms.
F=zeros(n,n);
for i=1:n
    F(i,i)=f(T(i,i));
end

for p=1:(n-1)
    for i=1:(n-p)
        j=i+p;
        s=T(i,j)*(F(j,j)-F(i,i));
        for k=(i+1):(j-1)
            s=s+T(i,k)*F(k,j)-F(i,k)*T(k,j);
        end
        F(i,j)=s/(T(j,j)-T(i,i));
    end
end

if(isreal(X))
    fX=real(Q*F*Q');
else
    fX=Q*F*Q';
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.


text:
%% set kalman parameters
function k = klmi(ocn, ncn, nzn, k, frame, st)
if frame  == st.st.st                                    % all objects in the first frame
k(size(ocn, 2)).s = [];                                  % initialize kalman filter
for      i = 1 : size(ocn, 2)                            % for every detected objects
k(i).s     = kalmani(ocn(:, i));                         % initialize every kalman filter
end
elseif frame ~= st.st.st                                 % not associated objects
n          = size(k, 2);
for      i = n + 1 : n + size(ncn, 2)                    % for every new objects
k(i).s     = kalmani(ncn(:, i - n));                     % initialize every new kalman filter
k(i).sz    = nzn(:, i - n);
end
end
end
%% initialize kalman
function s = kalmani(int)
dt         = 0.1;                                        % time step
s.A        = [1 dt 0 0; 0 1 0 0; 0 0 1 dt; 0 0 0 1];     % state transition matrix
s.H        = [1 0 0 0; 0 0 1 0];                         % observation matrix
s.Q        = 5 * eye(4);                                 % process noise covariance
s.R        = [0.5  0; 0 0.5];                            % measurement error covariance
s.x        = [int(1), 0, int(2), 0]';                    % a priori 'state vector' estimate
s.P        = 5 * eye(4);                                 % a priori estimate 'error covariance'
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ out ] = Distance1D( hist1, hist2 )
%hist1 and hist2 are sets of three historgrams.  One for each color (R,G,B)

dR = EMD1D(hist1(:,1),hist2(:,1));
dG = EMD1D(hist1(:,2),hist2(:,2));
dB = EMD1D(hist1(:,3),hist2(:,3));

out = dR + dG + dB;

end


text:
function [A,sumV] = quadcof(N,NW,order)
% Helper function to calculate the nonstationary quadratic inverse matrix
% Usage: [A,sumV] = quadcof(N,NW,order)
% N             (number of samples)
% NW: Time bandwidth product
% order: order (number of coefficients, upto 4NW)
%
% Outputs:
%
% A: quadratic inverse coefficient matrix
% sumV: sum of the quadratic inverse eigenvectors

A = zeros(2*NW,2*NW,order);
V = quadinv(N,NW);
[P,alpha] = dpss(N,NW,'calc');

for ii = 1:order

  for jj = 1:2*NW
    for kk = 1:2*NW
	A(jj,kk,ii) = sqrt(alpha(jj)*alpha(kk))*...
                      sum(P(:,jj).*P(:,kk).*V(:,ii));
    end;
  end;
end;
sumV=sum(V)/N;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function ellipse_t = ellFit(x,y,axis_handle,colorStr)
% Finds the best fit to an ellipse for the given set of points.
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=3215&objectType=File
%
% Format:   ellipse_t = ellFit( x,y, [axis_handle])
%
% Input:    x,y         - a set of points in 2 column vectors. AT LEAST 5 points are needed !
%           axis_handle - optional. a handle to an axis, at which the estimated ellipse
%                         will be drawn along with it's axes
%
% Output:   ellipse_t - structure that defines the best fit to an ellipse
%                       a           - sub axis (radius) of the X axis of the non-tilt ellipse
%                       b           - sub axis (radius) of the Y axis of the non-tilt ellipse
%                       phi         - orientation in radians of the ellipse (tilt)
%                       X0          - center at the X axis of the non-tilt ellipse
%                       Y0          - center at the Y axis of the non-tilt ellipse
%                       X0_in       - center at the X axis of the tilted ellipse
%                       Y0_in       - center at the Y axis of the tilted ellipse
%                       long_axis   - size of the long axis of the ellipse
%                       short_axis  - size of the short axis of the ellipse
%                       status      - status of detection of an ellipse
%
% Note:     if an ellipse was not detected (but a parabola or hyperbola), then
%           an empty structure is returned
%
% Example:  (Verify recovery)
%     a = 1.45; b = 3.2; theta = 0.15; radSpacing = 0.05;
%     [x,y] = ellipsePoints(a,b,theta,radSpacing);
%     figure(1), plot(x,y,'-'); axis equal
%     ellS = ellFit(x,y);
%     [ellS.a, ellS.b, ellS.phi]
%     [a, b, theta]
%
%     a = 1.45; b = 3.2; theta = pi/5; radSpacing = 0.05;
%     [x,y] = ellipsePoints(a,b,theta,radSpacing);
%     figure(1), plot(x,y,'-'); axis equal
%     for ii=1:10
%       ellS = ellFit(x + randn(size(x))/5, y + randn(size(y))/5);
%       aHat(ii) = ellS.a; bHat(ii) = ellS.b; thetaHat(ii) =  ellS.phi;
%     end
%     [mean(aHat), mean(bHat), mean(thetaHat)]
%
%     ellS = ellFit(x,y,figure(1));

%
% =====================================================================================
%                  Ellipse Fit using Least Squares criterion
% =====================================================================================
% We will try to fit the best ellipse to the given measurements. the mathematical
% representation of use will be the CONIC Equation of the Ellipse which is:
%
%    Ellipse = a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
%
% The fit-estimation method of use is the Least Squares method (without any weights)
% The estimator is extracted from the following equations:
%
%    g(x,y;A) := a*x^2 + b*x*y + c*y^2 + d*x + e*y = f
%
%    where:
%       A   - is the vector of parameters to be estimated (a,b,c,d,e)
%       x,y - is a single measurement
%
% We will define the cost function to be:
%
%   Cost(A) := (g_c(x_c,y_c;A)-f_c)'*(g_c(x_c,y_c;A)-f_c)
%            = (X*A+f_c)'*(X*A+f_c)
%            = A'*X'*X*A + 2*f_c'*X*A + N*f^2
%
%   where:
%       g_c(x_c,y_c;A) - vector function of ALL the measurements
%                        each element of g_c() is g(x,y;A)
%       X              - a matrix of the form: [x_c.^2, x_c.*y_c, y_c.^2, x_c, y_c ]
%       f_c            - is actually defined as ones(length(f),1)*f
%
% Derivation of the Cost function with respect to the vector of parameters "A" yields:
%
%   A'*X'*X = -f_c'*X = -f*ones(1,length(f_c))*X = -f*sum(X)
%
% Which yields the estimator:
%
%       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%       |  A_least_squares = -f*sum(X)/(X'*X) ->(normalize by -f) = sum(X)/(X'*X)  |
%       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% (We will normalize the variables by (-f) since "f" is unknown and can be accounted for later on)
%
% NOW, all that is left to do is to extract the parameters from the Conic Equation.
% We will deal the vector A into the variables: (A,B,C,D,E) and assume F = -1;
%
%    Recall the conic representation of an ellipse:
%
%       A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
%
% We will check if the ellipse has a tilt (=orientation). The orientation is present
% if the coefficient of the term "x*y" is not zero. If so, we first need to remove the
% tilt of the ellipse.
%
% If the parameter "B" is not equal to zero, then we have an orientation (tilt) to the ellipse.
% we will remove the tilt of the ellipse so as to remain with a conic representation of an
% ellipse without a tilt, for which the math is more simple:
%
% Non tilt conic rep.:  A`*x^2 + C`*y^2 + D`*x + E`*y + F` = 0
%
% We will remove the orientation using the following substitution:
%
%   Replace x with cx+sy and y with -sx+cy such that the conic representation is:
%
%   A(cx+sy)^2 + B(cx+sy)(-sx+cy) + C(-sx+cy)^2 + D(cx+sy) + E(-sx+cy) + F = 0
%
%   where:      c = cos(phi)    ,   s = sin(phi)
%
%   and simplify...
%
%       x^2(A*c^2 - Bcs + Cs^2) + xy(2A*cs +(c^2-s^2)B -2Ccs) + ...
%           y^2(As^2 + Bcs + Cc^2) + x(Dc-Es) + y(Ds+Ec) + F = 0
%
%   The orientation is easily found by the condition of (B_new=0) which results in:
%
%   2A*cs +(c^2-s^2)B -2Ccs = 0  ==> phi = 1/2 * atan( b/(c-a) )
%
%   Now the constants   c=cos(phi)  and  s=sin(phi)  can be found, and from them
%   all the other constants A`,C`,D`,E` can be found.
%
%   A` = A*c^2 - B*c*s + C*s^2                  D` = D*c-E*s
%   B` = 2*A*c*s +(c^2-s^2)*B -2*C*c*s = 0      E` = D*s+E*c
%   C` = A*s^2 + B*c*s + C*c^2
%
% Next, we want the representation of the non-tilted ellipse to be as:
%
%       Ellipse = ( (X-X0)/a )^2 + ( (Y-Y0)/b )^2 = 1
%
%       where:  (X0,Y0) is the center of the ellipse
%               a,b     are the ellipse "radiuses" (or sub-axis)
%
% Using a square completion method we will define:
%
%       F`` = -F` + (D`^2)/(4*A`) + (E`^2)/(4*C`)
%
%       Such that:    a`*(X-X0)^2 = A`(X^2 + X*D`/A` + (D`/(2*A`))^2 )
%                     c`*(Y-Y0)^2 = C`(Y^2 + Y*E`/C` + (E`/(2*C`))^2 )
%
%       which yields the transformations:
%
%           X0  =   -D`/(2*A`)
%           Y0  =   -E`/(2*C`)
%           a   =   sqrt( abs( F``/A` ) )
%           b   =   sqrt( abs( F``/C` ) )
%
% And finally we can define the remaining parameters:
%
%   long_axis   = 2 * max( a,b )
%   short_axis  = 2 * min( a,b )
%   Orientation = phi
%
%

if notDefined('colorStr');
    colorStr = 'r';
end
% initialize
orientation_tolerance = 1e-3;

% empty warning stack
warning( '' );

% prepare vectors, must be column vectors
x = x(:);
y = y(:);

% remove bias of the ellipse - to make matrix inversion more accurate. (will be added later on).
mean_x = mean(x);
mean_y = mean(y);
x = x-mean_x;
y = y-mean_y;

% the estimation for the conic equation of the ellipse
X = [x.^2, x.*y, y.^2, x, y ];
a = sum(X)/(X'*X);

% check for warnings
if ~isempty( lastwarn )
    disp( 'stopped because of a warning regarding matrix inversion' );
    ellipse_t = [];
    return
end

% extract parameters from the conic equation
[a,b,c,d,e] = deal( a(1),a(2),a(3),a(4),a(5) );

% remove the orientation from the ellipse
if ( min(abs(b/a),abs(b/c)) > orientation_tolerance )

    orientation_rad = 1/2 * atan( b/(c-a) );
    cos_phi = cos( orientation_rad );
    sin_phi = sin( orientation_rad );
    [a,b,c,d,e] = deal(...
        a*cos_phi^2 - b*cos_phi*sin_phi + c*sin_phi^2,...
        0,...
        a*sin_phi^2 + b*cos_phi*sin_phi + c*cos_phi^2,...
        d*cos_phi - e*sin_phi,...
        d*sin_phi + e*cos_phi );
    [mean_x,mean_y] = deal( ...
        cos_phi*mean_x - sin_phi*mean_y,...
        sin_phi*mean_x + cos_phi*mean_y );
else
    orientation_rad = 0;
    cos_phi = cos( orientation_rad );
    sin_phi = sin( orientation_rad );
end

% check if conic equation represents an ellipse
test = a*c;
switch (1)
case (test>0),  status = '';
case (test==0), status = 'Parabola found';  warning( 'fit_ellipse: Did not locate an ellipse' );
case (test<0),  status = 'Hyperbola found'; warning( 'fit_ellipse: Did not locate an ellipse' );
end

% if we found an ellipse return it's data
if (test>0)

    % make sure coefficients are positive as required
    if (a<0), [a,c,d,e] = deal( -a,-c,-d,-e ); end

    % final ellipse parameters
    X0          = mean_x - d/2/a;
    Y0          = mean_y - e/2/c;
    F           = 1 + (d^2)/(4*a) + (e^2)/(4*c);
    [a,b]       = deal( sqrt( F/a ),sqrt( F/c ) );
    long_axis   = 2*max(a,b);
    short_axis  = 2*min(a,b);

    % rotate the axes backwards to find the center point of the original TILTED ellipse
    R           = [ cos_phi sin_phi; -sin_phi cos_phi ];
    P_in        = R * [X0;Y0];
    X0_in       = P_in(1);
    Y0_in       = P_in(2);

    % pack ellipse into a structure
    ellipse_t = struct( ...
        'a',a,...
        'b',b,...
        'phi',orientation_rad,...
        'X0',X0,...
        'Y0',Y0,...
        'X0_in',X0_in,...
        'Y0_in',Y0_in,...
        'long_axis',long_axis,...
        'short_axis',short_axis,...
        'status','' );
else
    % report an empty structure
    ellipse_t = struct( ...
        'a',[],...
        'b',[],...
        'phi',[],...
        'X0',[],...
        'Y0',[],...
        'X0_in',[],...
        'Y0_in',[],...
        'long_axis',[],...
        'short_axis',[],...
        'status',status );
end

% check if we need to plot an ellipse with it's axes.
if (nargin>2) && ~isempty( axis_handle ) && (test>0)

    % rotation matrix to rotate the axes with respect to an angle phi
    R = [ cos_phi sin_phi; -sin_phi cos_phi ];

    % the axes
    ver_line        = [ [X0 X0]; Y0+b*[-1 1] ];
    horz_line       = [ X0+a*[-1 1]; [Y0 Y0] ];
    new_ver_line    = R*ver_line;
    new_horz_line   = R*horz_line;

    % the ellipse
    theta_r         = linspace(0,2*pi);
    ellipse_x_r     = X0 + a*cos( theta_r );
    ellipse_y_r     = Y0 + b*sin( theta_r );
    rotated_ellipse = R * [ellipse_x_r;ellipse_y_r];

    % draw
    hold_state = get( axis_handle,'NextPlot' );
    set( axis_handle,'NextPlot','add' );
%     plot( new_ver_line(1,:),new_ver_line(2,:),colorStr );
%     plot( new_horz_line(1,:),new_horz_line(2,:),colorStr );
    plot( rotated_ellipse(1,:),rotated_ellipse(2,:),colorStr );
    set( axis_handle,'NextPlot',hold_state );
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ROIInterpretedType = initROIInterpretedType
% ROIInterpretedType
%       Written DK Sept-2009
%
%
% This function lists in alphabetical order ROIInterpretedType as per DICOM PS 3.3 - 2008
% Page 519.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

ROIInterpretedType = struct(...
    'AVOIDANCE'                                  ,  'region in which dose is to be minimized',...
    'BOLUS'                                          ,  'patient bolus to be used for external beam therapy',...
    'BRACHY_ACCESSORY'                    ,  'brachytherapy accessory device',...
    'BRACHY_CHANNEL'                       ,  'brachytherapy channel',...
    'BRACHY_CHNL_SHLD'                   ,  'brachytherapy channel shield',...
    'BRACHY_SRC_APP'                        ,  'brachytherapy source applicator',...
    'CAVITY'                                         ,  'patient anatomical cavity',...
    'CONTRAST_AGENT'                       ,  'volume into which a contrast agent has been injected',...
    'CONTROL'                                     ,  'ROI to be used in control of dose optimization and calculation',...
    'CTV'                                              ,  'Clinical Target Volume (as defined in ICRU50)',...
    'DOSE_REGION'                             ,  'ROI to be used as a dose reference',...
    'EXTERNAL'                                    ,  'external patient contour',...
    'FIXATION'                                     ,  'external patient fixation or immobilisation device',...
    'GTV'                                              ,  'Gross Tumor Volume (as defined in ICRU50)',...
    'IRRAD_VOLUME'                            ,  'Irradiated Volume (as defined in ICRU50)',...
    'ISOCENTER'                                  ,  'treatment isocenter to be used for external beam therapy',...
    'MARKER'                                       ,  'patient marker or marker on a localizer',...
    'ORGAN'                                         ,  'patient organ',...
    'PTV'                                               ,  'Planning Target Volume (as defined in ICRU50)',...
    'REGISTRATION'                             ,  'registration ROI',...
    'SUPPORT'                                      ,  'external patient support device',...
    'TREATED_VOLUME'                        ,  'Treated Volume (as defined in ICRU50)'...
);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function obj = Q_pTop_func_new(in1,in2,in3)
coefs_tq1_11 = in2(31);
coefs_tq2_11 = in2(32);
coefs_tq3_11 = in2(33);
coef_Jacob1_qt_syms11 = in3(41);
coef_Jacob1_qt_syms12 = in3(45);
coef_Jacob1_qt_syms13 = in3(49);
coef_Jacob1_qt_syms14 = in3(53);
coef_Jacob1_qt_syms21 = in3(81);
coef_Jacob1_qt_syms22 = in3(85);
coef_Jacob1_qt_syms23 = in3(89);
coef_Jacob1_qt_syms24 = in3(93);
coef_Jacob1_qt_syms28 = in3(109);
coef_Jacob1_qt_syms29 = in3(113);
coef_Jacob1_qt_syms30 = in3(117);
coef_Jacob1_qt_syms31 = in3(121);
coef_Jacob1_qt_syms33 = in3(129);
coef_Jacob1_qt_syms34 = in3(133);
coef_Jacob1_qt_syms35 = in3(137);
coef_Jacob1_qt_syms36 = in3(141);
coef_Jacob2_qt_syms11 = in3(42);
coef_Jacob2_qt_syms12 = in3(46);
coef_Jacob2_qt_syms13 = in3(50);
coef_Jacob2_qt_syms14 = in3(54);
coef_Jacob2_qt_syms21 = in3(82);
coef_Jacob2_qt_syms22 = in3(86);
coef_Jacob2_qt_syms23 = in3(90);
coef_Jacob2_qt_syms24 = in3(94);
coef_Jacob2_qt_syms28 = in3(110);
coef_Jacob2_qt_syms29 = in3(114);
coef_Jacob2_qt_syms30 = in3(118);
coef_Jacob2_qt_syms31 = in3(122);
coef_Jacob2_qt_syms33 = in3(130);
coef_Jacob2_qt_syms34 = in3(134);
coef_Jacob2_qt_syms35 = in3(138);
coef_Jacob2_qt_syms36 = in3(142);
coef_Jacob3_qt_syms11 = in3(43);
coef_Jacob3_qt_syms12 = in3(47);
coef_Jacob3_qt_syms13 = in3(51);
coef_Jacob3_qt_syms14 = in3(55);
coef_Jacob3_qt_syms21 = in3(83);
coef_Jacob3_qt_syms22 = in3(87);
coef_Jacob3_qt_syms23 = in3(91);
coef_Jacob3_qt_syms24 = in3(95);
coef_Jacob3_qt_syms28 = in3(111);
coef_Jacob3_qt_syms29 = in3(115);
coef_Jacob3_qt_syms30 = in3(119);
coef_Jacob3_qt_syms31 = in3(123);
coef_Jacob3_qt_syms33 = in3(131);
coef_Jacob3_qt_syms34 = in3(135);
coef_Jacob3_qt_syms35 = in3(139);
coef_Jacob3_qt_syms36 = in3(143);
coef_Jacob4_qt_syms11 = in3(44);
coef_Jacob4_qt_syms12 = in3(48);
coef_Jacob4_qt_syms13 = in3(52);
coef_Jacob4_qt_syms14 = in3(56);
coef_Jacob4_qt_syms21 = in3(84);
coef_Jacob4_qt_syms22 = in3(88);
coef_Jacob4_qt_syms23 = in3(92);
coef_Jacob4_qt_syms24 = in3(96);
coef_Jacob4_qt_syms28 = in3(112);
coef_Jacob4_qt_syms29 = in3(116);
coef_Jacob4_qt_syms30 = in3(120);
coef_Jacob4_qt_syms31 = in3(124);
coef_Jacob4_qt_syms33 = in3(132);
coef_Jacob4_qt_syms34 = in3(136);
coef_Jacob4_qt_syms35 = in3(140);
coef_Jacob4_qt_syms36 = in3(144);
pinvG1_1 = in1(1);
pinvG1_2 = in1(4);
pinvG1_3 = in1(7);
pinvG2_1 = in1(2);
pinvG2_2 = in1(5);
pinvG2_3 = in1(8);
pinvG3_1 = in1(3);
pinvG3_2 = in1(6);
pinvG3_3 = in1(9);
obj = reshape([coef_Jacob1_qt_syms14+coefs_tq1_11.*coef_Jacob1_qt_syms11.*pinvG1_1+coefs_tq1_11.*coef_Jacob1_qt_syms12.*pinvG2_1+coefs_tq1_11.*coef_Jacob1_qt_syms13.*pinvG3_1+coefs_tq2_11.*coef_Jacob1_qt_syms11.*pinvG1_2+coefs_tq2_11.*coef_Jacob1_qt_syms12.*pinvG2_2+coefs_tq2_11.*coef_Jacob1_qt_syms13.*pinvG3_2+coefs_tq3_11.*coef_Jacob1_qt_syms11.*pinvG1_3+coefs_tq3_11.*coef_Jacob1_qt_syms12.*pinvG2_3+coefs_tq3_11.*coef_Jacob1_qt_syms13.*pinvG3_3,coef_Jacob2_qt_syms14+coefs_tq1_11.*coef_Jacob2_qt_syms11.*pinvG1_1+coefs_tq1_11.*coef_Jacob2_qt_syms12.*pinvG2_1+coefs_tq1_11.*coef_Jacob2_qt_syms13.*pinvG3_1+coefs_tq2_11.*coef_Jacob2_qt_syms11.*pinvG1_2+coefs_tq2_11.*coef_Jacob2_qt_syms12.*pinvG2_2+coefs_tq2_11.*coef_Jacob2_qt_syms13.*pinvG3_2+coefs_tq3_11.*coef_Jacob2_qt_syms11.*pinvG1_3+coefs_tq3_11.*coef_Jacob2_qt_syms12.*pinvG2_3+coefs_tq3_11.*coef_Jacob2_qt_syms13.*pinvG3_3,coef_Jacob3_qt_syms14+coefs_tq1_11.*coef_Jacob3_qt_syms11.*pinvG1_1+coefs_tq1_11.*coef_Jacob3_qt_syms12.*pinvG2_1+coefs_tq1_11.*coef_Jacob3_qt_syms13.*pinvG3_1+coefs_tq2_11.*coef_Jacob3_qt_syms11.*pinvG1_2+coefs_tq2_11.*coef_Jacob3_qt_syms12.*pinvG2_2+coefs_tq2_11.*coef_Jacob3_qt_syms13.*pinvG3_2+coefs_tq3_11.*coef_Jacob3_qt_syms11.*pinvG1_3+coefs_tq3_11.*coef_Jacob3_qt_syms12.*pinvG2_3+coefs_tq3_11.*coef_Jacob3_qt_syms13.*pinvG3_3,coef_Jacob4_qt_syms14+coefs_tq1_11.*coef_Jacob4_qt_syms11.*pinvG1_1+coefs_tq1_11.*coef_Jacob4_qt_syms12.*pinvG2_1+coefs_tq1_11.*coef_Jacob4_qt_syms13.*pinvG3_1+coefs_tq2_11.*coef_Jacob4_qt_syms11.*pinvG1_2+coefs_tq2_11.*coef_Jacob4_qt_syms12.*pinvG2_2+coefs_tq2_11.*coef_Jacob4_qt_syms13.*pinvG3_2+coefs_tq3_11.*coef_Jacob4_qt_syms11.*pinvG1_3+coefs_tq3_11.*coef_Jacob4_qt_syms12.*pinvG2_3+coefs_tq3_11.*coef_Jacob4_qt_syms13.*pinvG3_3,coef_Jacob1_qt_syms24+coefs_tq1_11.*coef_Jacob1_qt_syms21.*pinvG1_1+coefs_tq1_11.*coef_Jacob1_qt_syms22.*pinvG2_1+coefs_tq1_11.*coef_Jacob1_qt_syms23.*pinvG3_1+coefs_tq2_11.*coef_Jacob1_qt_syms21.*pinvG1_2+coefs_tq2_11.*coef_Jacob1_qt_syms22.*pinvG2_2+coefs_tq2_11.*coef_Jacob1_qt_syms23.*pinvG3_2+coefs_tq3_11.*coef_Jacob1_qt_syms21.*pinvG1_3+coefs_tq3_11.*coef_Jacob1_qt_syms22.*pinvG2_3+coefs_tq3_11.*coef_Jacob1_qt_syms23.*pinvG3_3,coef_Jacob2_qt_syms24+coefs_tq1_11.*coef_Jacob2_qt_syms21.*pinvG1_1+coefs_tq1_11.*coef_Jacob2_qt_syms22.*pinvG2_1+coefs_tq1_11.*coef_Jacob2_qt_syms23.*pinvG3_1+coefs_tq2_11.*coef_Jacob2_qt_syms21.*pinvG1_2+coefs_tq2_11.*coef_Jacob2_qt_syms22.*pinvG2_2+coefs_tq2_11.*coef_Jacob2_qt_syms23.*pinvG3_2+coefs_tq3_11.*coef_Jacob2_qt_syms21.*pinvG1_3+coefs_tq3_11.*coef_Jacob2_qt_syms22.*pinvG2_3+coefs_tq3_11.*coef_Jacob2_qt_syms23.*pinvG3_3,coef_Jacob3_qt_syms24+coefs_tq1_11.*coef_Jacob3_qt_syms21.*pinvG1_1+coefs_tq1_11.*coef_Jacob3_qt_syms22.*pinvG2_1+coefs_tq1_11.*coef_Jacob3_qt_syms23.*pinvG3_1+coefs_tq2_11.*coef_Jacob3_qt_syms21.*pinvG1_2+coefs_tq2_11.*coef_Jacob3_qt_syms22.*pinvG2_2+coefs_tq2_11.*coef_Jacob3_qt_syms23.*pinvG3_2+coefs_tq3_11.*coef_Jacob3_qt_syms21.*pinvG1_3+coefs_tq3_11.*coef_Jacob3_qt_syms22.*pinvG2_3+coefs_tq3_11.*coef_Jacob3_qt_syms23.*pinvG3_3,coef_Jacob4_qt_syms24+coefs_tq1_11.*coef_Jacob4_qt_syms21.*pinvG1_1+coefs_tq1_11.*coef_Jacob4_qt_syms22.*pinvG2_1+coefs_tq1_11.*coef_Jacob4_qt_syms23.*pinvG3_1+coefs_tq2_11.*coef_Jacob4_qt_syms21.*pinvG1_2+coefs_tq2_11.*coef_Jacob4_qt_syms22.*pinvG2_2+coefs_tq2_11.*coef_Jacob4_qt_syms23.*pinvG3_2+coefs_tq3_11.*coef_Jacob4_qt_syms21.*pinvG1_3+coefs_tq3_11.*coef_Jacob4_qt_syms22.*pinvG2_3+coefs_tq3_11.*coef_Jacob4_qt_syms23.*pinvG3_3,coef_Jacob1_qt_syms31+coefs_tq1_11.*coef_Jacob1_qt_syms28.*pinvG1_1+coefs_tq1_11.*coef_Jacob1_qt_syms29.*pinvG2_1+coefs_tq1_11.*coef_Jacob1_qt_syms30.*pinvG3_1+coefs_tq2_11.*coef_Jacob1_qt_syms28.*pinvG1_2+coefs_tq2_11.*coef_Jacob1_qt_syms29.*pinvG2_2+coefs_tq2_11.*coef_Jacob1_qt_syms30.*pinvG3_2+coefs_tq3_11.*coef_Jacob1_qt_syms28.*pinvG1_3+coefs_tq3_11.*coef_Jacob1_qt_syms29.*pinvG2_3+coefs_tq3_11.*coef_Jacob1_qt_syms30.*pinvG3_3,coef_Jacob2_qt_syms31+coefs_tq1_11.*coef_Jacob2_qt_syms28.*pinvG1_1+coefs_tq1_11.*coef_Jacob2_qt_syms29.*pinvG2_1+coefs_tq1_11.*coef_Jacob2_qt_syms30.*pinvG3_1+coefs_tq2_11.*coef_Jacob2_qt_syms28.*pinvG1_2+coefs_tq2_11.*coef_Jacob2_qt_syms29.*pinvG2_2+coefs_tq2_11.*coef_Jacob2_qt_syms30.*pinvG3_2+coefs_tq3_11.*coef_Jacob2_qt_syms28.*pinvG1_3+coefs_tq3_11.*coef_Jacob2_qt_syms29.*pinvG2_3+coefs_tq3_11.*coef_Jacob2_qt_syms30.*pinvG3_3,coef_Jacob3_qt_syms31+coefs_tq1_11.*coef_Jacob3_qt_syms28.*pinvG1_1+coefs_tq1_11.*coef_Jacob3_qt_syms29.*pinvG2_1+coefs_tq1_11.*coef_Jacob3_qt_syms30.*pinvG3_1+coefs_tq2_11.*coef_Jacob3_qt_syms28.*pinvG1_2+coefs_tq2_11.*coef_Jacob3_qt_syms29.*pinvG2_2+coefs_tq2_11.*coef_Jacob3_qt_syms30.*pinvG3_2+coefs_tq3_11.*coef_Jacob3_qt_syms28.*pinvG1_3+coefs_tq3_11.*coef_Jacob3_qt_syms29.*pinvG2_3+coefs_tq3_11.*coef_Jacob3_qt_syms30.*pinvG3_3,coef_Jacob4_qt_syms31+coefs_tq1_11.*coef_Jacob4_qt_syms28.*pinvG1_1+coefs_tq1_11.*coef_Jacob4_qt_syms29.*pinvG2_1+coefs_tq1_11.*coef_Jacob4_qt_syms30.*pinvG3_1+coefs_tq2_11.*coef_Jacob4_qt_syms28.*pinvG1_2+coefs_tq2_11.*coef_Jacob4_qt_syms29.*pinvG2_2+coefs_tq2_11.*coef_Jacob4_qt_syms30.*pinvG3_2+coefs_tq3_11.*coef_Jacob4_qt_syms28.*pinvG1_3+coefs_tq3_11.*coef_Jacob4_qt_syms29.*pinvG2_3+coefs_tq3_11.*coef_Jacob4_qt_syms30.*pinvG3_3,coef_Jacob1_qt_syms36+coefs_tq1_11.*coef_Jacob1_qt_syms33.*pinvG1_1+coefs_tq1_11.*coef_Jacob1_qt_syms34.*pinvG2_1+coefs_tq1_11.*coef_Jacob1_qt_syms35.*pinvG3_1+coefs_tq2_11.*coef_Jacob1_qt_syms33.*pinvG1_2+coefs_tq2_11.*coef_Jacob1_qt_syms34.*pinvG2_2+coefs_tq2_11.*coef_Jacob1_qt_syms35.*pinvG3_2+coefs_tq3_11.*coef_Jacob1_qt_syms33.*pinvG1_3+coefs_tq3_11.*coef_Jacob1_qt_syms34.*pinvG2_3+coefs_tq3_11.*coef_Jacob1_qt_syms35.*pinvG3_3,coef_Jacob2_qt_syms36+coefs_tq1_11.*coef_Jacob2_qt_syms33.*pinvG1_1+coefs_tq1_11.*coef_Jacob2_qt_syms34.*pinvG2_1+coefs_tq1_11.*coef_Jacob2_qt_syms35.*pinvG3_1+coefs_tq2_11.*coef_Jacob2_qt_syms33.*pinvG1_2+coefs_tq2_11.*coef_Jacob2_qt_syms34.*pinvG2_2+coefs_tq2_11.*coef_Jacob2_qt_syms35.*pinvG3_2+coefs_tq3_11.*coef_Jacob2_qt_syms33.*pinvG1_3+coefs_tq3_11.*coef_Jacob2_qt_syms34.*pinvG2_3+coefs_tq3_11.*coef_Jacob2_qt_syms35.*pinvG3_3,coef_Jacob3_qt_syms36+coefs_tq1_11.*coef_Jacob3_qt_syms33.*pinvG1_1+coefs_tq1_11.*coef_Jacob3_qt_syms34.*pinvG2_1+coefs_tq1_11.*coef_Jacob3_qt_syms35.*pinvG3_1+coefs_tq2_11.*coef_Jacob3_qt_syms33.*pinvG1_2+coefs_tq2_11.*coef_Jacob3_qt_syms34.*pinvG2_2+coefs_tq2_11.*coef_Jacob3_qt_syms35.*pinvG3_2+coefs_tq3_11.*coef_Jacob3_qt_syms33.*pinvG1_3+coefs_tq3_11.*coef_Jacob3_qt_syms34.*pinvG2_3+coefs_tq3_11.*coef_Jacob3_qt_syms35.*pinvG3_3,coef_Jacob4_qt_syms36+coefs_tq1_11.*coef_Jacob4_qt_syms33.*pinvG1_1+coefs_tq1_11.*coef_Jacob4_qt_syms34.*pinvG2_1+coefs_tq1_11.*coef_Jacob4_qt_syms35.*pinvG3_1+coefs_tq2_11.*coef_Jacob4_qt_syms33.*pinvG1_2+coefs_tq2_11.*coef_Jacob4_qt_syms34.*pinvG2_2+coefs_tq2_11.*coef_Jacob4_qt_syms35.*pinvG3_2+coefs_tq3_11.*coef_Jacob4_qt_syms33.*pinvG1_3+coefs_tq3_11.*coef_Jacob4_qt_syms34.*pinvG2_3+coefs_tq3_11.*coef_Jacob4_qt_syms35.*pinvG3_3],[4,4]);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Compares the speed of CG implementations
clc;
clearvars;
close all;
%rng default;

N = 32;

T = 1000;
t1 = 0;
t2 = 0;
for i=1:T
    if mod(i, 10) == 0
        fprintf('.');
    end
    A = spx.dict.simple.gaussian_mtx(N, N);
    A = A' * A;
    x = randn(N, 1);
    b  = A * x;
    tolerance = 1e-2;
    max_iterations = N * 5;

    tic;
    [x, res, iter ] = cgsolve(A, b, tolerance, max_iterations, 0);
    t1 = t1 + toc;

    options.tolerance = tolerance;
    options.max_iterations = max_iterations;
    tic;
    result = spx.fast.cg(A, b, options);
    t2 = t2 + toc;
end

fprintf('\n');
fprintf('Time taken by MATLAB implementation: %0.3f\n', t1);
fprintf('Time taken by C++ implementation: %0.3f\n', t2);
fprintf('Gain: %0.2f\n', t1/t2);

author: indigits

repo: sparse-plex


text:
function [ y , ysup ] = exp_rnd(x,rnd)
%EXP_RND      Rigorous bounds for exp(x) and input vector x according to rnd
%
%   [ y , ysup ] = exp_rnd(x,rnd)
%
%If specified with two output arguments, lower and upper bound is computed
%independent of rnd, otherwise y is rounded according to rnd.
%Rounding need not be to nearest after leaving exp_rnd
%

% written  12/30/98     S.M. Rump
% modified 08/31/99     S.M. Rump  extreme input, improved accuracy,
%                                  major revision
% modified 02/14/01     S.M. Rump  improved accuracy
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  INTLAB_STDFCTS_EXP = getappdata(0,'INTLAB_STDFCTS_EXP');

  infsup = ( nargout==2 );
  if infsup
    rnd = -1;
  end

  setround(0)
  xint = fix(x);                            % integer part of x, chopped
  xfrac = x - xint;                         % fractional part, -1 < xfrac < 1

  indexlarge = ( xint>709 );                % care for exceptions
  indexsmall = ( xint<-744 );
  index = indexlarge | indexsmall | isnan(x);
  xint(index) = 0;
  xfrac(index) = 0;

  xs = pow2( floor(pow2(xfrac,14)) , -14 );  % max. 14 bits of mantissa of xfrac,
                                             % no bit below 2^-14
  d = xfrac - xs;                            % 0 <= d < 2^-14, exactly repr.
  expxs = exp(xs);
  factor = INTLAB_STDFCTS_EXP.EPS;
  Exint = INTLAB_STDFCTS_EXP.POW(745+xint);
  Exinteps = INTLAB_STDFCTS_EXP.POWSUP(745+xint);

  % general case, exp(xfrac) = exp(xs)*exp(d)
  if infsup          % 0 <= err <= exp(d)*d^4/4! < 0.2501*d*d^3/3!

    % calculate upper bound
    setround(1)
    % exp(d)*expxs <= corr + expxs
    corr = (( ( 1+0.2501*d ) .* expxs.*d/3 + expxs ).*d/2 + expxs ).*d;

    % exp(xs) <= expxs*(1+EPS)
    % exp(xfrac) = exp(d)*exp(xs) <= (corr+expxs)*(1+EPS)
    ysup = ( corr + corr*factor ) + expxs*factor ;

    % exp(xfrac) <= expxs+ysup,  exp(x) = exp(xfrac)*exp(xint)
    ysup = expxs.*Exint + ( (expxs+ysup).*Exinteps + ysup.*Exint );
    ysup(x==0) = 1;

  end

  setround(rnd)
  if rnd==-1
    Exinteps = INTLAB_STDFCTS_EXP.POWINF(745+xint);
    corr1 = 1;
  else
    corr1 = 1 + 0.2501*d;
  end

  % exp(d)*expxs ~ corr + expxs    subject to rounding
  corr = (( corr1 .* expxs.*d/3 + expxs ).*d/2 + expxs ).*d;

  % exp(xs)  in  expxs*(1+/-EPS)
  % exp(xfrac) = exp(d)*exp(xs)  in  (corr+expxs)*(1+rnd*EPS)
  y = ( corr + (corr*rnd)*factor ) + (expxs*rnd)*factor;

  % exp(xfrac) ~ expxs+y,  exp(x) = exp(xfrac)*exp(xint)
  y = expxs.*Exint + ( (expxs+y).*Exinteps + y.*Exint );
  y(x==0) = 1;

  % large or small input, exceptions
  if any(index)

    if any(indexlarge)
      if infsup
        y(indexlarge) = realmax;
        ysup(indexlarge) = inf;
      else
        if rnd==-1
          y(indexlarge) = realmax;
        else
          y(indexlarge) = inf;
        end
      end
    end

    if any(indexsmall)
      INTLAB_INTVAL_ETA = realmin*eps;  % smallest positive denormalized fl-pt
      if infsup
        y(indexsmall) = 0;
        ysup(indexsmall) = INTLAB_INTVAL_ETA;
      else
        if rnd==-1
          y(indexsmall) = 0;
        else
          y(indexsmall) = INTLAB_INTVAL_ETA;
        end
      end
    end

    index = isnan(x);
    if any(index)
      y(index) = NaN;
      if infsup
        ysup(index) = NaN;
      end
    end

  end

  index = ( xint<-720 );
  if any(index)
    y(index) = max( y(index) , 0 );
  end

repo: OpenMAS


text:
function [ know, x ] = p03_sol ( n )

%*****************************************************************************80
%
%% P03_SOL returns the solution for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, real X(N), the solution, if known.
%
  know = 0;

  x = zeros ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%==========================================================================
% Recursively walks through a Matlab hierarchy and converts it to the
% hierarchy of java.util.ArrayListS and java.util.MapS. Then calls
% Snakeyaml to write it to a file.
%=========================================================================
function result = WriteYaml(filename, data, flowstyle)
    if ~exist('flowstyle','var')
        flowstyle = 0;
    end;
    if ~ismember(flowstyle, [0,1])
        error('Flowstyle must be 0,1 or empty.');
    end;
    result = [];
    [pth,~,~] = fileparts(mfilename('fullpath'));
    try
        import('org.yaml.snakeyaml.*');
        javaObject('Yaml');
    catch
        dp = [pth filesep 'external' filesep 'snakeyaml-1.9.jar'];
        if not(ismember(dp, javaclasspath ('-dynamic')))
        	javaaddpath(dp); % javaaddpath clears global variables...!?
        end
        import('org.yaml.snakeyaml.*');
    end;
    javastruct = scan(data);
    dumperopts = DumperOptions();
    dumperopts.setLineBreak(...
        javaMethod('getPlatformLineBreak',...
        'org.yaml.snakeyaml.DumperOptions$LineBreak'));
    if flowstyle
        classes = dumperopts.getClass.getClasses;
        flds = classes(3).getDeclaredFields();
        fsfld = flds(1);
        if ~strcmp(char(fsfld.getName), 'FLOW')
            error(['Accessed another field instead of FLOW. Please correct',...
            'class/field indices (this error maybe caused by new snakeyaml version).']);
        end;
        dumperopts.setDefaultFlowStyle(fsfld.get([]));
    end;
    yaml = Yaml(dumperopts);

    output = yaml.dump(javastruct);
    if ~isempty(filename)
        fid = fopen(filename,'w');
        fprintf(fid,'%s',char(output) );
        fclose(fid);
    else
        result = output;
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan(r)
    if ischar(r)
        result = scan_char(r);
    elseif iscell(r)
        result = scan_cell(r);
    elseif isord(r)
        result = scan_ord(r);
    elseif isstruct(r)
        result = scan_struct(r);
    elseif isnumeric(r)
        result = scan_numeric(r);
    elseif islogical(r)
        result = scan_logical(r);
    elseif isa(r,'DateTime')
        result = scan_datetime(r);
    else
        error(['Cannot handle type: ' class(r)]);
    end
end

%--------------------------------------------------------------------------
%
%
function result = scan_numeric(r)
    if isempty(r)
        result = java.util.ArrayList();
    else
        result = java.lang.Double(r);
    end
end

%--------------------------------------------------------------------------
%
%

function result = scan_logical(r)
    if isempty(r)
        result = java.util.ArrayList();
    else
        result = java.lang.Boolean(r);
    end
end

%--------------------------------------------------------------------------
%
%
function result = scan_char(r)
    if isempty(r)
        result = java.util.ArrayList();
    else
        result = java.lang.String(r);
    end
end

%--------------------------------------------------------------------------
%
%
function result = scan_datetime(r)
    % datestr 30..in ISO8601 format
    %java.text.SimpleDateFormat('yyyymmdd'T'HH:mm:ssz" );

    [Y, M, D, H, MN,S] = datevec(double(r));
	result = java.util.GregorianCalendar(Y, M-1, D, H, MN,S);
	result.setTimeZone(java.util.TimeZone.getTimeZone('UTC'));

    %tz = java.util.TimeZone.getTimeZone('UTC');
    %cal = java.util.GregorianCalendar(tz);
    %cal.set
    %result = java.util.Date(datestr(r));
end

%--------------------------------------------------------------------------
%
%
function result = scan_cell(r)
    if(isrowvector(r))
        result = scan_cell_row(r);
    elseif(iscolumnvector(r))
        result = scan_cell_column(r);
    elseif(ismymatrix(r))
        result = scan_cell_matrix(r);
    elseif(issingle(r));
        result = scan_cell_single(r);
    elseif(isempty(r))
        result = java.util.ArrayList();
    else
        error('Unknown cell content.');
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_ord(r)
    if(isrowvector(r))
        result = scan_ord_row(r);
    elseif(iscolumnvector(r))
        result = scan_ord_column(r);
    elseif(ismymatrix(r))
        result = scan_ord_matrix(r);
    elseif(issingle(r))
        result = scan_ord_single(r);
    elseif(isempty(r))
        result = java.util.ArrayList();
    else
        error('Unknown ordinary array content.');
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_cell_row(r)
    result = java.util.ArrayList();
    for ii = 1:size(r,2)
        result.add(scan(r{ii}));
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_cell_column(r)
    result = java.util.ArrayList();
    for ii = 1:size(r,1)
        tmp = r{ii};
        if ~iscell(tmp)
            tmp = {tmp};
        end;
        result.add(scan(tmp));
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_cell_matrix(r)
    result = java.util.ArrayList();
    for ii = 1:size(r,1)
        i = r(ii,:);
        result.add(scan_cell_row(i));
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_cell_single(r)
    result = java.util.ArrayList();
    result.add(scan(r{1}));
end

%--------------------------------------------------------------------------
%
%
function result = scan_ord_row(r)
    result = java.util.ArrayList();
    for i = r
        result.add(scan(i));
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_ord_column(r)
    result = java.util.ArrayList();
    for i = 1:size(r,1)
        result.add(scan_ord_row(r(i)));
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_ord_matrix(r)
    result = java.util.ArrayList();
    for i = r'
        result.add(scan_ord_row(i'));
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_ord_single(r)
    result = java.util.ArrayList();
    for i = r'
        result.add(r);
    end;
end

%--------------------------------------------------------------------------
%
%
function result = scan_struct(r)
    result = java.util.LinkedHashMap();
    for i = fields(r)'
        key = i{1};
        val = r.(key);
        result.put(key,scan(val));
    end;
end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function [baseMVA, bus, gen, branch, areas, gencost] = c30
% CASE_IEEE30    Power flow data for IEEE 30 bus test case.
%   Please see 'help caseformat' for details on the case file format.
%   This data was converted from IEEE Common Data Format
%   (ieee30cdf.txt) on 20-Sep-2004 by cdf2matp, rev. 1.11
%   See end of file for warnings generated during conversion.
%
%   Converted from IEEE CDF file from:
%       http://www.ee.washington.edu/research/pstca/
%
%  08/20/93 UW ARCHIVE           100.0  1961 W IEEE 30 Bus Test Case

%   MATPOWER
%   $Id: case_ieee30.m,v 1.2 2004/09/21 01:45:05 ray Exp $

%%-----  Power Flow Data  -----%%
%% system MVA base
baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
bus = [
	1	3	0	0	0	0	1	1.06	0	132	1	1.06	0.94;
	2	2	21.7	12.7	0	0	1	1.043	-5.48	132	1	1.06	0.94;
	3	1	2.4	1.2	0	0	1	1.021	-7.96	132	1	1.06	0.94;
	4	1	7.6	1.6	0	0	1	1.012	-9.62	132	1	1.06	0.94;
	5	2	94.2	19	0	0	1	1.01	-14.37	132	1	1.06	0.94;
	6	1	0	0	0	0	1	1.01	-11.34	132	1	1.06	0.94;
	7	1	22.8	10.9	0	0	1	1.002	-13.12	132	1	1.06	0.94;
	8	2	30	30	0	0	1	1.01	-12.1	132	1	1.06	0.94;
	9	1	0	0	0	0	1	1.051	-14.38	1	1	1.06	0.94;
	10	1	5.8	2	0	19	1	1.045	-15.97	33	1	1.06	0.94;
	11	2	0	0	0	0	1	1.082	-14.39	11	1	1.06	0.94;
	12	1	11.2	7.5	0	0	1	1.057	-15.24	33	1	1.06	0.94;
	13	2	0	0	0	0	1	1.071	-15.24	11	1	1.06	0.94;
	14	1	6.2	1.6	0	0	1	1.042	-16.13	33	1	1.06	0.94;
	15	1	8.2	2.5	0	0	1	1.038	-16.22	33	1	1.06	0.94;
	16	1	3.5	1.8	0	0	1	1.045	-15.83	33	1	1.06	0.94;
	17	1	9	5.8	0	0	1	1.04	-16.14	33	1	1.06	0.94;
	18	1	3.2	0.9	0	0	1	1.028	-16.82	33	1	1.06	0.94;
	19	1	9.5	3.4	0	0	1	1.026	-17	33	1	1.06	0.94;
	20	1	2.2	0.7	0	0	1	1.03	-16.8	33	1	1.06	0.94;
	21	1	17.5	11.2	0	0	1	1.033	-16.42	33	1	1.06	0.94;
	22	1	0	0	0	0	1	1.033	-16.41	33	1	1.06	0.94;
	23	1	3.2	1.6	0	0	1	1.027	-16.61	33	1	1.06	0.94;
	24	1	8.7	6.7	0	4.3	1	1.021	-16.78	33	1	1.06	0.94;
	25	1	0	0	0	0	1	1.017	-16.35	33	1	1.06	0.94;
	26	1	3.5	2.3	0	0	1	1	-16.77	33	1	1.06	0.94;
	27	1	0	0	0	0	1	1.023	-15.82	33	1	1.06	0.94;
	28	1	0	0	0	0	1	1.007	-11.97	132	1	1.06	0.94;
	29	1	2.4	0.9	0	0	1	1.003	-17.06	33	1	1.06	0.94;
	30	1	10.6	1.9	0	0	1	0.992	-17.94	33	1	1.06	0.94;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin
gen = [
	1	260.2	-16.1	10	0	1.06	100	1	360.2	0;
	2	40	50	50	-40	1.045	100	1	140	0;
	5	0	37	40	-40	1.01	100	1	100	0;
	8	0	37.3	40	-10	1.01	100	1	100	0;
	11	0	16.2	24	-6	1.082	100	1	100	0;
	13	0	10.6	24	-6	1.071	100	1	100	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status
branch = [
	1	2	0.0192	0.0575	0.0528	9900	0	0	0	0	1;
	1	3	0.0452	0.1652	0.0408	9900	0	0	0	0	1;
	2	4	0.057	0.1737	0.0368	9900	0	0	0	0	1;
	3	4	0.0132	0.0379	0.0084	9900	0	0	0	0	1;
	2	5	0.0472	0.1983	0.0418	9900	0	0	0	0	1;
	2	6	0.0581	0.1763	0.0374	9900	0	0	0	0	1;
	4	6	0.0119	0.0414	0.009	9900	0	0	0	0	1;
	5	7	0.046	0.116	0.0204	9900	0	0	0	0	1;
	6	7	0.0267	0.082	0.017	9900	0	0	0	0	1;
	6	8	0.012	0.042	0.009	9900	0	0	0	0	1;
	6	9	0	0.208	0	9900	0	0	0.978	0	1;
	6	10	0	0.556	0	9900	0	0	0.969	0	1;
	9	11	0	0.208	0	9900	0	0	0	0	1;
	9	10	0	0.11	0	9900	0	0	0	0	1;
	4	12	0	0.256	0	9900	0	0	0.932	0	1;
	12	13	0	0.14	0	9900	0	0	0	0	1;
	12	14	0.1231	0.2559	0	9900	0	0	0	0	1;
	12	15	0.0662	0.1304	0	9900	0	0	0	0	1;
	12	16	0.0945	0.1987	0	9900	0	0	0	0	1;
	14	15	0.221	0.1997	0	9900	0	0	0	0	1;
	16	17	0.0524	0.1923	0	9900	0	0	0	0	1;
	15	18	0.1073	0.2185	0	9900	0	0	0	0	1;
	18	19	0.0639	0.1292	0	9900	0	0	0	0	1;
	19	20	0.034	0.068	0	9900	0	0	0	0	1;
	10	20	0.0936	0.209	0	9900	0	0	0	0	1;
	10	17	0.0324	0.0845	0	9900	0	0	0	0	1;
	10	21	0.0348	0.0749	0	9900	0	0	0	0	1;
	10	22	0.0727	0.1499	0	9900	0	0	0	0	1;
	21	22	0.0116	0.0236	0	9900	0	0	0	0	1;
	15	23	0.1	0.202	0	9900	0	0	0	0	1;
	22	24	0.115	0.179	0	9900	0	0	0	0	1;
	23	24	0.132	0.27	0	9900	0	0	0	0	1;
	24	25	0.1885	0.3292	0	9900	0	0	0	0	1;
	25	26	0.2544	0.38	0	9900	0	0	0	0	1;
	25	27	0.1093	0.2087	0	9900	0	0	0	0	1;
	28	27	0	0.396	0	9900	0	0	0.968	0	1;
	27	29	0.2198	0.4153	0	9900	0	0	0	0	1;
	27	30	0.3202	0.6027	0	9900	0	0	0	0	1;
	29	30	0.2399	0.4533	0	9900	0	0	0	0	1;
	8	28	0.0636	0.2	0.0428	9900	0	0	0	0	1;
	6	28	0.0169	0.0599	0.013	9900	0	0	0	0	1;
];
mm=length(branch(:,1));
for i=1:mm
    if branch(i,9)==0;
        branch(i,9)=1;
    else
    end
end
%%-----  OPF Data  -----%%
%% area data
areas = [
	1	1;
];

%% generator cost data
%	1	startup	shutdown	n	x0	y0	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
gencost = [
	2	0	0	3	.5*0.02	2	0;
	2	0	0	3	.5*0.0175	1.75	0;
	2	0	0	3	.5*0.0625	1	0;
	2	0	0	3	.5*0.00834	3.25	0;
	2	0	0	3	.5*0.025	3	0;
	2	0	0	3	.5*0.025	3	0;
];

return;

% Warnings from cdf2matp conversion:
%
% ***** Qmax = Qmin at generator at bus    1 (Qmax set to Qmin + 10)
% ***** area data conversion not yet implemented (creating dummy area data)
% ***** MVA limit of branch 1 - 2 not given, set to 9900
% ***** MVA limit of branch 1 - 3 not given, set to 9900
% ***** MVA limit of branch 2 - 4 not given, set to 9900
% ***** MVA limit of branch 3 - 4 not given, set to 9900
% ***** MVA limit of branch 2 - 5 not given, set to 9900
% ***** MVA limit of branch 2 - 6 not given, set to 9900
% ***** MVA limit of branch 4 - 6 not given, set to 9900
% ***** MVA limit of branch 5 - 7 not given, set to 9900
% ***** MVA limit of branch 6 - 7 not given, set to 9900
% ***** MVA limit of branch 6 - 8 not given, set to 9900
% ***** MVA limit of branch 6 - 9 not given, set to 9900
% ***** MVA limit of branch 6 - 10 not given, set to 9900
% ***** MVA limit of branch 9 - 11 not given, set to 9900
% ***** MVA limit of branch 9 - 10 not given, set to 9900
% ***** MVA limit of branch 4 - 12 not given, set to 9900
% ***** MVA limit of branch 12 - 13 not given, set to 9900
% ***** MVA limit of branch 12 - 14 not given, set to 9900
% ***** MVA limit of branch 12 - 15 not given, set to 9900
% ***** MVA limit of branch 12 - 16 not given, set to 9900
% ***** MVA limit of branch 14 - 15 not given, set to 9900
% ***** MVA limit of branch 16 - 17 not given, set to 9900
% ***** MVA limit of branch 15 - 18 not given, set to 9900
% ***** MVA limit of branch 18 - 19 not given, set to 9900
% ***** MVA limit of branch 19 - 20 not given, set to 9900
% ***** MVA limit of branch 10 - 20 not given, set to 9900
% ***** MVA limit of branch 10 - 17 not given, set to 9900
% ***** MVA limit of branch 10 - 21 not given, set to 9900
% ***** MVA limit of branch 10 - 22 not given, set to 9900
% ***** MVA limit of branch 21 - 22 not given, set to 9900
% ***** MVA limit of branch 15 - 23 not given, set to 9900
% ***** MVA limit of branch 22 - 24 not given, set to 9900
% ***** MVA limit of branch 23 - 24 not given, set to 9900
% ***** MVA limit of branch 24 - 25 not given, set to 9900
% ***** MVA limit of branch 25 - 26 not given, set to 9900
% ***** MVA limit of branch 25 - 27 not given, set to 9900
% ***** MVA limit of branch 28 - 27 not given, set to 9900
% ***** MVA limit of branch 27 - 29 not given, set to 9900
% ***** MVA limit of branch 27 - 30 not given, set to 9900
% ***** MVA limit of branch 29 - 30 not given, set to 9900
% ***** MVA limit of branch 8 - 28 not given, set to 9900
% ***** MVA limit of branch 6 - 28 not given, set to 9900

author: Sable

repo: mcbench-benchmarks


text:
function [y, dzdg, dzdb] = vl_nnnonorm(x, g, b, varargin)
%VL_NNNONORM applies weights and biases, but does no normalization
%   Y = VL_NNNONORM(X,G,B) applies a set of gains and biases to
%   the input X with shape HxWxCxN. "No normalization" is defined as:
%
%      Y(i,j,k,t) = G(k') * X(i,j,k,t) + B(k')
%
%   where
%      k' = group_idx(k,C,G), where N_G is the number of groups and
%        group_idx(k,C,G) := floor(k / (C/N_G)).
%
%   VL_NNGNORM(..., 'option', value, ...) takes the following option:
%
%   This layer was largely inspired by this blog post:
%       http://www.offconvex.org/2018/03/02/acceleration-overparameterization/
%
% Copyright (C) 2018 Samuel Albanie
% All rights reserved.

  [~,dzdy] = vl_argparsepos(struct(), varargin) ;

  expectedSz = [1 1 size(x,3) 1] ;
  sg = size(g) ; sb = size(b) ;
  assert(all(expectedSz(1:numel(sg)) == sg), 'GAINS have unexpected size') ;
  assert(all(expectedSz(1:numel(sb)) == sb), 'BIASES have unexpected size') ;

  if isempty(dzdy)
    y = bsxfun(@times, g, x) ; % apply gain
    y = bsxfun(@plus, y, b) ; % add bias
  else
    dzdy = dzdy{1} ;
    dzdb = chanSum(dzdy) ;
    dzdg = chanSum(x .* dzdy) ;
    dzdx = bsxfun(@times, dzdy, g) ;
    y = dzdx ;
  end

% -----------------------
function res = chanSum(x)
% -----------------------
  res = sum(sum(sum(x, 1), 2), 4) ;

repo: MFT


text:
% OP_GRADSYMU_V_OTIMES_N: assemble the matrix A = [A(i,j)], A(i,j) = (epsilon (gradsym u)_j, (v \otimes n)_i ).
%
%   A = op_gradsymu_v_otimes_n (spu, spv, msh, coeff);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_vector/sp_eval_boundary_side)
%   spv:   structure representing the space of test functions (see sp_vector/sp_eval_boundary_side)
%   msh:   structure containing the domain partition and the quadrature rule for the boundary,
%           since it must contain the normal vector (see msh_cartesian/msh_eval_boundary_side)
%   coeff: vector-valued function f, evaluated at the quadrature points
%
% OUTPUT:
%
%   A: assembled matrix
%
% Copyright (C) 2015, 2017, 2020 Rafael Vazquez
%
%    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/>.

function varargout = op_gradsymu_v_otimes_n (spu, spv, msh, coeff)

  gradu = reshape (spu.shape_function_gradients, spu.ncomp, [], ...
		   msh.nqn, spu.nsh_max, msh.nel);

  shpv = reshape (spv.shape_functions, spv.ncomp, msh.nqn, spv.nsh_max, msh.nel);

  ncomp = size (gradu, 1);
  ndir = size (gradu, 2);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      gradu_iel = reshape (gradu(:,:,:,:,iel), spu.ncomp, ndir, msh.nqn, spu.nsh_max);
      gradu_iel = 0.5 * (gradu_iel + permute(gradu_iel, [2 1 3 4]));
      gradu_iel = reshape (gradu_iel, spu.ncomp*ndir, msh.nqn, 1, spu.nsh_max);
      shpv_iel = reshape (shpv(:, :, :, iel), spv.ncomp, msh.nqn, spv.nsh_max);

      v_otimes_n_iel = zeros (ncomp, ndir, msh.nqn, spv.nsh_max);
      normal_iel = msh.normal (:, :, iel);
      for ii = 1:spv.ncomp
        for jj = 1:ndir
          v_otimes_n_iel(ii,jj,:,:) = bsxfun (@times, shpv_iel(ii,:,:), normal_iel(jj,:));
        end
      end
      v_otimes_n_iel = reshape (v_otimes_n_iel, ncomp*ndir, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1, msh.nqn, 1]);
      v_oxn_times_jw = bsxfun (@times, jacdet_iel, v_otimes_n_iel);
      tmp1 = sum (bsxfun (@times, v_oxn_times_jw, gradu_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_gradsymu_v_otimes_n: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_gradsymu_v_otimes_n: wrong number of output arguments')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function issingle = getSinglePersonImages(annolist_all)

load('/BS/leonid-pose/work/data/new_dataset/dataset_release_candidate1/test_annolist_all_annotated','annolist');
load('/home/andriluk/IMAGES/human_pose_dataset/dataset/dataset_release_candidate1/dataset_info_v12');
imgidxs = find(DATASET.img_train == 0);

names = cell(length(annolist),1);
idxs = [];
for imgidx = 1:length(annolist)
    if (isfield(annolist(imgidx),'annorect') && length(annolist(imgidx).annorect) == 1)
        names{imgidx} = ['im' padZeros(num2str(imgidxs(imgidx)),5)];
        idxs = [idxs; imgidx];
    end
end

names = names(idxs);

issingle = zeros(length(annolist_all),1);
for imgidx = 1:length(annolist_all)
    [~,n] = fileparts(annolist_all(imgidx).image.name);
    issingle(imgidx) = ~isempty(find(strcmp(names,n(1:7))));
end

% assert(length(names) == sum(idxs));

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut