text:
function out = spm_run_bms_map (job)
% Run Bayesian Model Selection Maps
% SPM job execution function
% takes a harvested job data structure and calls SPM functions to perform
% Bayesian Inference for Model Selection of Log. Evidence Maps
% Input:
% job    - harvested job data structure (see matlabbatch help)
% Output:
% out    - computation results, usually a struct variable.
%
%
% Bayesian Inference on Model Space:
%
% The Random-effects 'RFX' method is described in Stephan et al. [1]
% 'Bayesian Model Selection for Group Studies'.
% Output files (for each model):
%       BMS.mat
%       Exceedance Probability Maps (*epm.<ext>),
%       Posterior Probability Maps (*ppm.<ext>),
%       Dirichlet Paramters (alpha) Maps (*alpha.<ext>).
%
% The Fixed-effects 'FFX' method adds together the log-evidences over
% subjects/sessions for each group, then compares the group log-ev's.
% This is also known as the Group Bayes Factor (GBF) approach [2].
% Output files (for each model):
%       BMS.mat
%       Posterior Probability Maps (*ppm.<ext>).
%
% BMS contains:
%     BMS.fname
%     BMS.map.ffx(rfx).data
%     BMS.map.ffx(rfx).ppm
%     BMS.map.ffx(rfx).xppm     - only for RFX (this is the expected posterior
%                                 probability map ie. posterior mean)
%     BMS.map.ffx(rfx).epm      - only for RFX (optional) - this is the
%                                 exceedance probability map
%     BMS.map.ffx(rfx).alpha    - only for RFX
%
% [1] Rosa et al., 2009, Bayesian Model Selection Maps for Group Studies,
% NeuroImage.
% [2] Stephan et al., 2009, Bayesian Model Selection for Group Studies,
% NeuroImage.
% [3] Penny et al., 2004, Comparing Dynamic Causal Models, NeuroImage.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Maria Joao Rosa
% $Id: spm_run_bms_map.m 5740 2013-11-13 12:00:04Z guillaume $

% Input
% -------------------------------------------------------------------------
direct  = job.dir{1};
fname   = fullfile(direct,'BMS.mat'); % Output filename (with full path)
mask    = length(job.mask{1});        % Mask image
if mask
   mask_image = spm_vol(job.mask);    % Mask image Vol
end
nsamps    = str2num(job.nsamp);       % Number of samples (nmodels > 3)
do_maps   = job.out_file;
do_ecp    = do_maps > 0;              % Compute Exceedance Probability
do_alpha  = do_maps > 1;              % Compute Alpha Parameters

% Nb. of subjects and models
% -------------------------------------------------------------------------
nsubjs    = size(job.sess_map,2);
nmodels   = size(job.sess_map{1}(1).mod_map,1);
nsess     = size(job.sess_map{1},2);
nnames    = size(job.mod_name,2);
names     = job.mod_name;

% method
% -------------------------------------------------------------------------
if strcmp(job.method_maps,'FFX');
    method = 'FFX';
else
    method = 'RFX';
end

% Name models
% -------------------------------------------------------------------------
if nnames < nmodels
    for nn=1:nmodels-nnames
        names = [names, sprintf('m%d',nn)];
    end
end

if size(unique(names(1:nmodels)),2) < nmodels,
    id = 'Indentical names for different models!';  % Same name!
    error(id);
end

if nsubjs == 1
   method = 'FFX';                              % If only 1 subject do FFX
end

% Sort out log-evidence images dimensions
% -------------------------------------------------------------------------
Vol_models(1,1) = spm_vol(job.sess_map{1}(1).mod_map(1));

first_vol       = Vol_models(1,1);
M               = first_vol{1}.mat;
DIM             = first_vol{1}.dim(1:3)';

xdim            = DIM(1);
ydim            = DIM(2);
zdim            = DIM(3);
[xords,yords]   = ndgrid(1:xdim,1:ydim);
xords           = xords(:)';
yords           = yords(:)';
I               = 1:xdim*ydim;
zords_init      = ones(1,xdim*ydim);

% Setup images
% -------------------------------------------------------------------------
switch method

    case 'FFX',   % Fixed Effects

        % Check if BMS.mat exists
        if exist(fullfile(job.dir{1},'BMS.mat'),'file')
           load(fname);
           if  isfield(BMS,'map') && isfield(BMS.map,'ffx')
               str = { 'Warning: existing BMS.mat file has been over-written!'};
               msgbox(str)
           end
        end

        % Save BMS data
        out.files{1} = fname;

        % Create PPM image files for each model
        model_ppm(1:nmodels) = struct(...
        'fname',    '',...
        'dim',      DIM',...
        'dt',       [spm_type('float32') spm_platform('bigend')],...
        'mat',      M,...
        'pinfo',    [1 0 0]',...
        'n', [1 1], ...
        'descrip',  '');

        % Load Vols for all subjects/models
        for i = 1:nmodels
            model_ppm(i).fname   = fullfile(direct,[sprintf('%s_model_ppm',names{i}) spm_file_ext]);
            model_ppm(i).descrip = sprintf('PPM: %s model',names{i});
            BMS.map.ffx.ppm{i}   = model_ppm(i).fname;

            for s = 1:nsubjs,
                for se = 1:nsess,
                    nsessi      = size(job.sess_map{s},2);
                    nmodelsi    = size(job.sess_map{s}(se).mod_map,1);
                    if (nsess == nsessi && nmodels == nmodelsi)
                        Vol_models(s,i,se) = spm_vol(job.sess_map{s}(se).mod_map(i));
                        tmp = Vol_models(s,i,se);
                    else
                        msgbox('The number of sessions/models should be the same for all subjects!')
                        return
                    end
                    % Stop if log-ev images have different dimensions
                    if tmp{1}.dim(1)~=xdim || tmp{1}.dim(2)~=ydim || tmp{1}.dim(3)~=zdim
                       error('Log-evidence images must have the same dimensions!')
                    end
                end
            end
        end

        % Create files
        model_ppm = spm_create_vol(model_ppm);
        BMS.fname = fname;

        % Save data and BMS
        BMS.fname = fname;
        BMS.map.ffx.data = job.sess_map;
        save(out.files{1},'BMS', spm_get_defaults('mat.format'));

    case 'RFX',  % Random Effects

        % Check if BMS.mat exists
        if exist(fullfile(job.dir{1},'BMS.mat'),'file')
           load(fname);
           if  isfield(BMS,'map') && isfield(BMS.map,'rfx')
               str = { 'Warning: existing BMS.mat file has been over-written!'};
               msgbox(str)
           end
        end

        % BMS structure
        out.files{1}   = fname;

        % Create PPM image files for each model
        model_exp_r(1:nmodels) = struct(...
        'fname',    '',...
        'dim',      DIM',...
        'dt',       [spm_type('float32') spm_platform('bigend')],...
        'mat',      M,...
        'pinfo',    [1 0 0]',...
        'n', [1 1], ...
        'descrip',  '');

        if do_ecp
            % Create EPM image files for each model
            model_xp(1:nmodels) = struct(...
            'fname',    '',...
            'dim',      DIM',...
            'dt',       [spm_type('float32') spm_platform('bigend')],...
            'mat',      M,...
            'pinfo',    [1 0 0]',...
            'n', [1 1], ...
            'descrip',  '');
        end

        if do_alpha
            % Create alpha image files for each model
            model_alpha(1:nmodels) = struct(...
            'fname',    '',...
            'dim',      DIM',...
            'dt',       [spm_type('float32') spm_platform('bigend')],...
            'mat',      M,...
            'pinfo',    [1 0 0]',...
            'n', [1 1], ...
            'descrip',  '');
        end

        % Load Vols for all subjects/models
        for i = 1:nmodels
            model_exp_r(i).fname   = fullfile(direct,[sprintf('%s_model_xppm',names{i}) spm_file_ext]);
            model_exp_r(i).descrip = sprintf('Exp_r: %s model',names{i});
            BMS.map.rfx.ppm{i}     = model_exp_r(i).fname;
            if do_ecp
            model_xp(i).fname      = fullfile(direct,[sprintf('%s_model_epm',names{i}) spm_file_ext]);
            model_xp(i).descrip    = sprintf('XP: %s model',names{i});
            BMS.map.rfx.epm{i}     = model_xp(i).fname;
            end
            if do_alpha
            model_alpha(i).fname   = fullfile(direct,[sprintf('%s_model_alpha',names{i}) spm_file_ext]);
            model_alpha(i).descrip = sprintf('Alpha: %s model',names{i});
            BMS.map.rfx.alpha{i}   = model_alpha(i).fname;
            end
            for s = 1:nsubjs,
                for se = 1:nsess,
                    nsessi      = size(job.sess_map{s},2);
                    nmodelsi    = size(job.sess_map{s}(se).mod_map,1);
                    if (nsess == nsessi && nmodels == nmodelsi)
                        Vol_models(s,i,se) = spm_vol(job.sess_map{s}(se).mod_map(i));
                        tmp = Vol_models(s,i,se);
                    else
                        msgbox('The number of sessions/models should be the same for all subjects!')
                        return
                    end
                    % Stop if log-ev images have different dimensions
                    if tmp{1}.dim(1)~=xdim || tmp{1}.dim(2)~=ydim || tmp{1}.dim(3)~=zdim
                       error('Log-evidence images must have the same dimensions!')
                    end
                end
            end
        end

        % Create files
        model_exp_r              = spm_create_vol(model_exp_r);
        if do_ecp, model_xp      = spm_create_vol(model_xp); end
        if do_alpha, model_alpha = spm_create_vol(model_alpha); end

        % Save data and BMS
        BMS.fname = fname;
        BMS.map.rfx.data = job.sess_map;
        save(out.files{1},'BMS', spm_get_defaults('mat.format'));

end

% Progress bar
% -------------------------------------------------------------------------
spm_progress_bar('Init',zdim,'BMS Maps (Inference)','Slices complete');

% Loop through image slices
% -------------------------------------------------------------------------
for z = 1:zdim,

    spm_progress_bar('Set',z);                  % Update progress bar
    j = repmat(NaN,xdim,ydim);                  % Init. image values

    fprintf('%s%30s',repmat(sprintf('\b'),1,30),'Computing maps...')
    str   = sprintf('Slice %d out of %d',z,zdim); % Display slice nb.
    fprintf('\r%-40s: %30s',str,' ')

    zords   = z*zords_init;                     % Slice z
    xyz     = [xords(I); yords(I); zords(I)];   % Slice coordinates
    nVox    = size(xyz,2);                      % Nb. of voxels per slice

    if mask
        % Voxels inside mask
        mask_xyz  = mask_image{1}.mat\M*[xyz(:,1:nVox);ones(1,nVox)];
        gamma     = spm_get_data(mask_image{1},mask_xyz);
        b         = find(gamma>0.5);            % Voxels in the mask
    else
        b         = 1:nVox;                     % All voxels
    end

    z_models        = NaN(nsubjs,nmodels,nVox);       % Data
    z_models(1,1,:) = spm_get_data(first_vol{1},xyz); % Data: all subs/mods
    non_nan         = find(~isnan(z_models(1,1,:)));  % Voxels ~NaN

    % Find voxels ~NaN and sum sessions
    % ---------------------------------------------------------------------
    for s = 1:nsubjs,
        for k = 1:nmodels,
                sum_tmp_data    = [];
            for ns = 1:nsess,
                tmp_data        = Vol_models(s,k,ns);
                sum_tmp_data    = [sum_tmp_data; spm_get_data(tmp_data{1},xyz)];
            end
                z_models(s,k,:) = sum(sum_tmp_data,1);
                non_nani        = find(~isnan(z_models(s,k,:)));
                non_nan         = intersect(non_nan,non_nani);
        end
    end

    % Voxels to be analysed
    non_nan = intersect(non_nan,b);
    Nvoxels = length(non_nan);

    % Method
    % ---------------------------------------------------------------------
    switch method

          % Fixed Effects
          % ---------------------------------------------------------------
          case 'FFX',

              if Nvoxels > 0                % Slice with ~NaN voxels

                zz     = sum(z_models,1);   % Sum all subjects/sessions
                mz     = mean(zz,2);        % Get mean of all models
                zzmean = zeros(1,nmodels,length(mz));
                for jj = 1:nmodels
                    zzmean(1,jj,:) = mz;
                end

                if nmodels==1
                    % Process log Bayes factor image
                    zz  = exp(zz);              % Exponentiate log-bf values
                    % Calculate posterior probabiliy
                    pz  = zeros(1,1,length(zz));
                    pz  = zz./(1+zz);
                    j(non_nan)   = pz(1,1,non_nan);
                    model_ppm(1) = spm_write_plane(model_ppm(k),j,z);
                else
                    zz  = zz-zzmean;            % Subtract mean
                    zz  = exp(zz);              % Exponentiate log-ev values
                    tzz = sum(zz,2);            % Sum exp(log-ev.)

                    % Calculate posterior probabiliy
                    pz  = zeros(1,nmodels,length(zz));
                    for k = 1:nmodels,
                        pz(1,k,:)    = zz(1,k,:)./tzz;
                        j(non_nan)   = pz(1,k,non_nan);
                        model_ppm(k) = spm_write_plane(model_ppm(k),j,z);
                    end
                end

              else
                % Nvoxels = 0
                for k = 1:nmodels,
                    % Write NaN for slice z
                    model_ppm(k) = spm_write_plane(model_ppm(k),j,z);
                end
              end

          % Fixed Effects
          % ---------------------------------------------------------------
          case 'RFX',

                if Nvoxels > 0
                    % Initialise results
                    exp_r_total              = zeros(Nvoxels,nmodels);
                    if do_ecp, xp_total      = zeros(Nvoxels,nmodels); end
                    if do_alpha, alpha_total = zeros(Nvoxels,nmodels); end

                    % Do BMS in all voxels of slice z
                    for n = 1:Nvoxels,
                        lme = z_models(:,:,non_nan(n));

                        if nmodels==1
                            % Provide evidence for dummy null model
                            lme=0.5*[lme, -lme];
                            [alpha,exp_r,xp] = spm_BMS(lme,nsamps,0,0,do_ecp);

                            exp_r_total(n,:)              = exp_r(1);  % Cond. Expecta.
                            if do_ecp, xp_total(n,:)      = xp(1); end % Exceeda. Prob.
                            if do_alpha, alpha_total(n,:) = alpha(1); end % Dirichlet par.
                        else

                            % Group BMS
                            [alpha,exp_r,xp] = spm_BMS(lme,nsamps,0,0,do_ecp);

                            exp_r_total(n,:)              = exp_r;  % Cond. Expecta.
                            if do_ecp, xp_total(n,:)      = xp; end % Exceeda. Prob.
                            if do_alpha, alpha_total(n,:) = alpha; end % Dirichlet par.
                        end
                    end

                    % Write images
                    for i = 1:nmodels,
                        j(non_nan)     = exp_r_total(:,i);
                        model_exp_r(i) = spm_write_plane(model_exp_r(i),j,z);
                        if do_ecp
                        j(non_nan)     = xp_total(:,i);
                        model_xp(i)    = spm_write_plane(model_xp(i),j,z);
                        end
                        if do_alpha
                        j(non_nan)     = alpha_total(:,i);
                        model_alpha(i) = spm_write_plane(model_alpha(i),j,z);
                        end
                    end
                else
                    % Write images when Nvoxels = 0
                    for i = 1:nmodels,
                        model_exp_r(i) = spm_write_plane(model_exp_r(i),j,z);
                        if do_ecp, model_xp(i) = spm_write_plane(model_xp(i),j,z); end
                        if do_alpha, model_alpha(i) = spm_write_plane(model_alpha(i),j,z); end
                    end
                end

    end

end % Loop over slices

% Clear progress bar
% -------------------------------------------------------------------------
spm_progress_bar('Clear');
disp('Done.');

author: spm


text:
function out = gpuBenchReport( varargin )
%gpuBenchReport  create an HTML report based on some GPU performance data
%
%   gpuBenchReport(data) creates a new HTML report based on the supplied
%   PerformanceData and opens it in the system browser.
%
%   gpuBenchReport() creates an HTML report based only on the pre-stored
%   performance data.
%
%   filename = gpuBenchReport(...) returns the location of the main page
%   for the report instead of viewing it immediately.
%
%   Examples:
%   >> gpuBenchReport
%   >> data = gpuBench;
%   >> gpuBenchReport( data )
%
%   See also: gpuBench

%   Copyright 2011 The MathWorks, Inc.

narginchk( 0, 2 );
nargoutchk( 0, 1 );

% Try to get some comparison data
if nargin>0
    assert( all( cellfun( 'isclass', varargin, 'gpubench.PerformanceData' ) ) );
    [allData,allDataSummary] = getComparisonData( [varargin{:}]' );
else
    [allData,allDataSummary] = getComparisonData();
end
if isempty(allData)
    error( 'gpuBenchReport:NoData', 'No data to report' );
end

N = numel( allData );
gpubench.multiWaitbar( 'Creating GPUBench report...', 0 );

% Ensure the output folder exists
reportDir = fullfile( tempdir(), 'GPUBenchReport' );
if ~exist( reportDir, 'dir' )
    mkdir( reportDir );
end
copyFiles( reportDir );

% Create the summary page for this device
makeSummaryBarChart( reportDir, allDataSummary );
reportname = makeSummaryPage( reportDir, allDataSummary );
gpubench.multiWaitbar( 'Creating GPUBench report...', 'Increment', 1/(N+1) );

% Now create the detail pages for all devices
for ii=1:numel( allData )
    makePerformancePlots( reportDir, allData, ii );
    makeDetailPage( reportDir, allData, allDataSummary, ii );
    gpubench.multiWaitbar( 'Creating GPUBench report...', 'Increment', 1/(N+1) );
end

gpubench.multiWaitbar( 'Creating GPUBench report...', 'close' );

if nargout
    out = reportname;
else
    web( reportname );
end

end % gpuBenchReport

%-------------------------------------------------------------------------%
function [allData,allDataSummary] = getComparisonData( data )

if nargin<1 || isempty(data)
    data = [];
    % No user data, so use the current release
    thisRelease = regexp( version, 'R\d*[ab]', 'match' );
    thisRelease = thisRelease{1};
else
    % Work out which data-file to use from the release
    thisRelease = data(1).MATLABRelease;
    % Flag the user's data so that we can highlight it
    for ii=1:numel(data)
        data(ii).IsSelected = true;
    end
end

% Try to load the data for this release
otherData = loadDataForRelease( thisRelease, data );

% If no data was found, warn and try an older release
if isempty(otherData)
    allDataFiles = dir( fullfile( 'data', 'R*.mat' ) );
    allDataFiles = {allDataFiles.name};

    % Try each file in turn, starting from the most recent
    for ii=numel(allDataFiles):-1:1
        release = allDataFiles{ii}(1:6);
        otherData = loadDataForRelease( release, data );
        if ~isempty(otherData)
            break;
        end
    end
    % If we reach here without any data then no file had what we wanted.
    % Give up!
    if isempty( otherData )
        error( 'GPUBenchReport:NoData', ...
        ['No pre-stored data was found. Please re-download and install ', ...
        'gpuBench from the File Exchange.'] );
    else
        % Indicate that we're using an older release
        warning( 'GPUBenchReport:NoDataForRelease', ...
            ['There is no pre-stored data for this MATLAB release (%s). ' ...
            'Using data for %s instead. ', ...
            'Check the File Exchange for an updated version of gpuBench ', ...
            'that has data for this release.'], ...
            thisRelease, release );
    end
end

% Construct the summary statistics from all the results and then sort the
% original data using the summary score.
if nargin>0
    allData = [data;otherData];
else
    allData = otherData;
end

if isempty( allData )
    % No data at all, so warn and show something older
    error( 'GPUBenchReport:NoData', ...
        ['There is no pre-stored data for this MATLAB release (%s). ' ...
        'Check the File Exchange for an updated version or use gpuBench ', ...
        'to view the results for your system.'], release );
elseif isempty( otherData )
end

allDataSummary = gpubench.SummaryData( allData );
allData = allData(allDataSummary.SortOrder);
end % getComparisonData

%-------------------------------------------------------------------------%
function otherData = loadDataForRelease( release, existingData )
datafilename = fullfile( gpubench.getDataDir(), [release, '.mat'] );
if exist( datafilename, 'file' ) == 2
    otherData = load( datafilename );
    if ~isfield( otherData, 'results' ) || ~isa( otherData.results, 'gpubench.PerformanceData' )
        error( 'GPUBench:BadDataFile', 'Data file was corrupted, please delete it: %s', datafilename );
    end
    otherData = otherData.results;

    % We don't want to see a result for the exact same card
    if ~isempty(existingData)
        keep = true( size( otherData ) );
        for ii=1:numel( otherData )
            for jj=1:numel(existingData)
                if ~existingData(jj).IsHostData
                    keep(ii) = keep(ii) ...
                        && ~strcmpi( otherData(ii).GPUInfo.Name, existingData(jj).GPUInfo.Name );
                end
            end
            otherData(ii).IsSelected = false;
        end
        otherData = otherData(keep);
    end
else
    otherData = gpubench.PerformanceData.empty(0,1);
end
end % loadDataForRelease

%-------------------------------------------------------------------------%
function makeSummaryBarChart( outDir, summaryData )

assert( isa( summaryData, 'gpubench.SummaryData' ) );

N = numel( summaryData.DeviceName );
deviceNames = summaryData.DeviceName;
functionNames = summaryData.LongName;
flopsResults = summaryData.PeakFLOPS;
isSelected = summaryData.IsSelectedDevice;

% Sort by data-type
types = unique( summaryData.Datatype );
numOfType = zeros(size(types));
colsForType = cell(size(types));
for ii=1:numel(types)
    match = strcmp( summaryData.Datatype, types{ii} );
    numOfType(ii) = sum( match );
    colsForType{ii} = find( match );
end
colOrder = [colsForType{:}];
functionNames = functionNames(colOrder);
flopsResults = flopsResults(:,colOrder);

figh = gpubench.makeFigure( 'Summary' );
set(figh,'Position',[10 10 650 650])
bars = barh( flopsResults'/1e9 );
set( gca, 'YTickLabel', functionNames, 'YDir', 'Reverse' )
xlabel( sprintf('GFLOPS\n(higher is better)') )

% Highlight the selected entry
for ii=1:N
    if isSelected(ii)
        deviceNames{ii} = ['{\bf',deviceNames{ii},' (yours)}'];
    end
end

gpubench.legend( deviceNames{:}, 'Location', 'SouthEast' );
grid on
title( 'Performance Summary' )

% Set colors to fade from blue to red
colors = [
    linspace(0,0.6,N)
    zeros(1,N)
    linspace(0.75,0,N)
    ]';
for ii=1:N
    set( bars(ii), 'FaceColor', colors(ii,:), 'EdgeColor', 0.5*colors(ii,:) );
end

% highlight the user's result
highlightColor = [0 1 0];
selectedResults = find(summaryData.IsSelectedDevice);
N = numel(selectedResults);
for ii=1:N
    ratio = (N-ii+2)/(N+1);
    set( bars(selectedResults(ii)), ...
        'FaceColor', highlightColor*ratio, ...
        'EdgeColor', [0 0 0], ...
        'Linewidth', 2 );
end

% Add dividers between types
hold on
x = get(gca,'XLim');
for ii=1:numel(numOfType)-1
    plot( x, (numOfType(ii)+0.5)*[1 1], 'k-' )
end

gpubench.addAreaShadows( gca() );
gpubench.addGradientToAxes( gca() );

% Save the image to file for the HTML to pick up
filename = 'summarychart.png';
gpubench.captureFigure( figh, fullfile( outDir, filename ), false );
close( figh );
end % makeSummaryBarChart

%-------------------------------------------------------------------------%
function makePerformancePlots( outDir, data, thisIdx )
%Create a FLOPS plot for each results in data with the "thisIdx" result
%highlighted.

fileBase = sprintf('device%d',thisIdx);

plotNames = cell( size(data) );
for ii=1:numel(data)
    plotNames{ii} = data(ii).getDeviceName();
end
plotNames{thisIdx} = ['{\bf',plotNames{thisIdx}, ' (selected)}'];

results = data(thisIdx).Results;

for rr=1:numel( results )
    name = [results(rr).FunctionName, ' (', results(rr).DataType, ')'];
    figh = gpubench.makeFigure( name );
    color = get( gca(), 'ColorOrder' );

    % Plot each device's curve, if they have the relevant data
    for ii=1:numel(data)
        colorIdx = mod(ii-1,size(color,1))+1;
        % Find the corresponding result
        if hasResult( data(ii), name )
            linewidth = 1+2*(ii==thisIdx);
            gpubench.plotFLOPS( getResult( data(ii), name ), color(colorIdx,:), linewidth )
        else
            % We still need to plot something for the legend to be correct
            plot( nan, nan, 'Color', color(colorIdx,:), 'LineWidth', 1 );
        end
    end
    % Add a highlight around the peak-flops point
    colorIdx = mod(thisIdx-1,size(color,1))+1;
    thisResult = getResult( data(thisIdx), name );
    [maxVal,maxIdx] = max( 1e-9 * thisResult.NumOps ./ thisResult.Times);
    plot( thisResult.Sizes(maxIdx), maxVal, ...
        'Color', color(colorIdx,:), ...
        'Marker', 'o', ...
        'MarkerSize', 16, ...
        'Linewidth', 2 );

    title( name );
    gpubench.legend( plotNames{:}, 'Location', 'NorthWest' );
    gpubench.addGradientToAxes( gca() );

    % Due to a bug in ZBuffer renderer, the log-lines disappear when we add
    % a patch to the plot. Stick them back in now.
    xtick = str2num( get( gca, 'XTickLabel' ) ); %#ok<ST2NM>
    ylim = get( gca, 'YLim' );
    for ii=1:numel(xtick)-1
        xminortick = (2:9) * power(10,xtick(ii));
        xdata = [xminortick;xminortick;nan(1,numel(xminortick))];
        ydata = [repmat(ylim',1,numel(xminortick)); nan(1,numel(xminortick))];
        line( 'XData', xdata(:), 'YData', ydata(:), ...
            'LineStyle', ':', ...
            'Color', 0.6*[1 1 1] );
    end

    % Save the image to file for the HTML to pick up
    filename = [fileBase,'-',results(rr).FunctionName,'-',results(rr).DataType,'.png'];
    gpubench.captureFigure( figh, fullfile( outDir, filename ) );
    close( figh );
end
end % makePerformancePlots

%-------------------------------------------------------------------------%
function copyFiles( outDir )
% Copy the required stylesheet and other files into the report folder

dataDir = gpubench.getDataDir();

files = {
    'gpubench.css'
    'title.png'
    'background.png'
    'warning.png'
    };
for ii=1:numel(files)
    copyfile( fullfile( dataDir, files{ii} ), outDir, 'f' );
end
end % copyFiles

%-------------------------------------------------------------------------%
function reportname = makeSummaryPage( outDir, summaryData )
% Create the summary page for this device

assert( isa( summaryData, 'gpubench.SummaryData' ) );

% Find the user's data
userIdx = find( summaryData.IsSelectedDevice, 1, 'first' );
if isempty( userIdx )
    % No user data, so use the first one
    pageName = '';
else
    pageName = [': ', summaryData.DeviceName{userIdx}];
end
reportname = fullfile( outDir, 'index.html' );

fid = fopen( reportname, 'wt' );
fprintf( fid, '\n' );
fprintf( fid, '  <title>GPU Comparison Report%s</title>\n', pageName );
fprintf( fid, '  <link rel="stylesheet" type="text/css" href="gpubench.css"/>\n' );
fprintf( fid, '  <meta author="Generated by GPUBench."/>\n' );
fprintf( fid, '\n' );
fprintf( fid, '\n' );

% All of the content goes in a giant table to control the width
fprintf( fid, '

\n' );
fprintf( fid, '  \n' );

% Main title
fprintf( fid, '
GPU Comparison Report%s
\n', pageName );

% Summary section
fprintf( fid, '
Summary of results
\n' );
% Add the description
writeBoilerPlate( fid, 'summary_intro.html' )

names = summaryData.FunctionName;
numCols = numel( names );

% Print the table header
fprintf( fid, '

\n' );
fprintf( fid, '
\n' );

% Titles for types
types = getDatatypes( summaryData );
colsForType = getColsForType( summaryData, types );
numOfType = cellfun( 'length', colsForType );
colOrder = [colsForType{:}];

fprintf( fid, '
	' );
for ii=1:numel(types)
    fprintf( fid, '      \n', ...
        numOfType(ii) );
    fprintf( fid, '        Results for data-type ''%s''
(In GFLOPS)	\n', ...
        types{ii} );
end

% Titles for each result
fprintf( fid, '
	' );
for cc=1:numCols
    fprintf( fid, '%s	', names{colOrder(cc)} );
end
fprintf( fid, '
\n' );
% Now the body
for rr=1:numel(summaryData.DeviceName)
    fprintf( fid, '
' );
    if (summaryData.IsSelectedDevice(rr))
        nameStr = ['',summaryData.DeviceName{rr},''];
        cellformat = '%1.2f	';
    else
        nameStr = summaryData.DeviceName{rr};
        cellformat = '%1.2f	';
    end
    fprintf( fid, '%s	', rr, nameStr );
    for cc=1:numCols
        colIdx = colOrder(cc);
        fprintf( fid, cellformat, rr, colIdx, summaryData.PeakFLOPS(rr,colIdx) / 1e9 );
    end
    fprintf( fid, '
\n' );
end
fprintf( fid, '
\n\n' );
fprintf( fid, '(click any device name or result to see the detailed data)
\n\n' );

% Add the summary image
fprintf( fid, '\n\n' );
fprintf( fid, '

\n' );

% Footer
fprintf( fid, '
\n' );
writeGeneratedBy( fid );

% Close the main layout table
fprintf( fid, '

\n' );
fprintf( fid, '\n' );
fprintf( fid, '\n' );
fclose( fid );
end % makeSummaryPage

%-------------------------------------------------------------------------%
function makeDetailPage( outDir, data, summaryData, thisIdx )
% Create page of detailed information for one device
name = summaryData.DeviceName{thisIdx};
fileBase = sprintf('device%d',thisIdx);

reportname = fullfile( outDir, [fileBase,'.html'] );

fid = fopen( reportname, 'wt' );
fprintf( fid, '\n' );
fprintf( fid, '  <title>GPU Performance Details: %s</title>\n', name );
fprintf( fid, '  <link rel="stylesheet" type="text/css" href="gpubench.css"/>\n' );
fprintf( fid, '  <meta author="Generated by GPUBench."/>\n' );
fprintf( fid, '\n' );
fprintf( fid, '\n' );

% All of the content goes in a giant table to control the width
fprintf( fid, '

\n' );
fprintf( fid, '  \n' );

% Main title
fprintf( fid, '
GPU Performance Details: %s
\n', name );

% Contents section
fprintf( fid, '

\n' );
fprintf( fid, '    Contents:\n' );
fprintf( fid, '  	\n' );
names = summaryData.LongName;
fprintf( fid, '
\n' );
fprintf( fid, '
System Configuration
\n' );

% Sort the contents by type
types = getDatatypes( summaryData );

for tt=1:numel( types )
    fprintf( fid, '
Results for datatype %s
\n', types{tt} );
    colsForType = getColsForType( summaryData, types{tt} );
    for nn=1:numel( colsForType )
        myCol = colsForType(nn);
        fprintf( fid, '
%s
\n', myCol, names{myCol} );
    end
    fprintf( fid, '

\n' );
end
fprintf( fid, '
\n' );
fprintf( fid, '

\n' );
fprintf( fid, '
\n' );

% Add a section showing the operating environment
fprintf( fid, '
System Configuration
\n' );

% If this isn't the user's system, make that clear
if ~data(thisIdx).IsSelected
    fprintf( fid, '

\n' );
    fprintf( fid, '    	\n' );
    fprintf( fid, '    Note that this is previously stored data and does not reflect your system configuration.	\n' );
    fprintf( fid, '

\n' );
end

fprintf( fid, '
MATLAB Release: %s
\n', data(thisIdx).MATLABRelease );
fprintf( fid, '

\n' );
fprintf( fid, '
Host
\n' );
writeStructTable( fid, data(thisIdx).CPUInfo );
% - only add the GPU if we ran on it
if ~data(thisIdx).IsHostData
    fprintf( fid, '  	\n' );
    fprintf( fid, '
GPU
\n' );
    writeStructTable( fid, data(thisIdx).GPUInfo );
end
fprintf( fid, '

\n' );
fprintf( fid, '
\n' );

% Add one section per result
names = summaryData.LongName;
for nn=1:numel( names )
    fprintf( fid, '
Results for %s
\n', nn, names{nn} );
    if ~hasResult( data(thisIdx), names{nn} )
        % No results for this function
        fprintf( fid, '
No results found for %s.
\n', names{nn} );
        continue;
    end
    myResult = getResult( data(thisIdx), names{nn} );
    % See if there's a description for this function
    writeBoilerPlate( fid, [myResult.FunctionName,'.html'] )

    fprintf( fid, '

\n' );

    fprintf( fid, '

\n');
    fprintf( fid, '      Raw data for %s - %s\n', name, names{nn} );
    fprintf( fid, '

\n');

    % Print the table header
    fprintf( fid, '
\n' );
    fprintf( fid, '
' );
    fprintf( fid, 'Array size
(elements)	' );
    fprintf( fid, 'Num
Operations	' );
    fprintf( fid, 'Time
(ms)	' );
    fprintf( fid, 'GigaFLOPS	' );
    fprintf( fid, '
\n' );
    % Now one row per size
    sizes = myResult.Sizes;
    flops = myResult.NumOps;
    times = myResult.Times;
    [~,peakIdx] = max( flops ./ times );
    baseFormatStr1 = '';
    baseFormatStr2 = '	';
    for ss=1:numel(sizes)
        % Highlight the peak FLOPS row
        if ss==peakIdx
            formatStr1 = [baseFormatStr1,''];
            formatStr2 = ['',baseFormatStr2];
        else
            formatStr1 = baseFormatStr1;
            formatStr2 = baseFormatStr2;
        end
        fprintf( fid, '
' );
        fprintf( fid, [formatStr1,'%s',formatStr2], num2strWithCommas(sizes(ss)) );
        fprintf( fid, [formatStr1,'%s',formatStr2], num2strWithCommas(flops(ss)) );
        fprintf( fid, [formatStr1,'%2.2f',formatStr2], times(ss)*1000 );
        fprintf( fid, [formatStr1,'%2.2f',formatStr2], flops(ss)/times(ss)/1e9 );
        fprintf( fid, '
\n' );
    end
    fprintf( fid, '
\n' );
    fprintf( fid, '(N gigaflops = Nx109 operations per second)
\n\n' );

    fprintf( fid, '

\n');

    % Add the image
    fprintf( fid, '  	\n' );
    fprintf( fid, '    \n', ...
        fileBase, myResult.FunctionName, myResult.DataType );
    fprintf( fid, '

\n' );
end

% Footer
fprintf( fid, '
\n' );

fprintf( fid, '

\n' );
writeGeneratedBy( fid );
fprintf( fid, '  	\n' );
fprintf( fid, '    Back to summary\n' );
fprintf( fid, '

\n' );

% Close the main layout table
fprintf( fid, '

\n' );
fprintf( fid, '\n' );
fprintf( fid, '\n' );
fclose( fid );
end % makeDetailPage

%-------------------------------------------------------------------------%
function writeStructTable( fid, data )
assert( isstruct( data ) && isscalar( data ) );
fprintf( fid, '
\n' );
fields = fieldnames( data );
for ff=1:numel( fields )
    fprintf( fid, '
%s	', fields{ff} );
    fprintf( fid, '' );
    x = data.(fields{ff});
    if ischar( x )
        fprintf( fid, '%s', x );
    elseif isinteger( x )
        fprintf( fid, '%d', x );
    else
        % Try to let MATLAB do the right thing
        fprintf( fid, '%g', x );
    end
    fprintf( fid, '
\n' );
end
fprintf( fid, '
\n' );
end % writeStructTable

%-------------------------------------------------------------------------%
function writeBoilerPlate( outFid, filename )
%Read some boiler-plate HTML and paste it into the supplied output file
filename = fullfile( gpubench.getDataDir(), filename );
inFid = fopen( filename, 'rt' );
if inFid<=0
    warning( 'gpuBenchReport:MissingBoilerPlateFile', ...
        'Input file could not be opened: %s', filename );
    return;
end
txt = fread( inFid );
fwrite( outFid, txt );
fclose( inFid );
end % writeBoilerPlate

%-------------------------------------------------------------------------%
function writeGeneratedBy( outFid )
%Write the "generated by" string into the footer

fprintf( outFid, 'Generated by gpuBench v%s: %s\n', ...
    gpubench.version(), datestr( now(), 'yyyy-mm-dd HH:MM:SS' ) );
end % writeGeneratedBy

%-------------------------------------------------------------------------%
function str = num2strWithCommas( num )
%Convert an integer into a string with commas separating sets of 3 digits
%
%  e.g. num2StrWithCommas(12345678) = '12,345,678'

% First convert using the standard method
baseStr = num2str( abs(num) );
% now insert some commas.
% pad to a multiple of 3
padding = 3 - (mod(length(baseStr)-1,3)+1);
str = [repmat(' ',1,padding), baseStr];
numCols = length(str)/3;
str = [reshape(str,3,numCols);repmat(',',1,numCols)];
str = strtrim( str(1:end-1) );
% Finally, re-insert the sign
if num<0
    str = ['-',str];
end
end % num2StrWithCommas

repo: Aero-matlab


text:
function [a_vec, a_vec_flt] = KiserTMOv(hdrv, filenameOutput, tmo_alpha, tmo_dn_clamping, tmo_white, tmo_gamma, tmo_quality, tmo_video_profile)
%
%
%      KiserTMOv(hdrv, filenameOutput, tmo_alpha, tmo_dn_clamping, tmo_white, tmo_gamma, tmo_quality, tmo_video_profile)
%
%
%       Input:
%           -hdrv: a HDR video structure; use hdrvread to create a hdrv
%           structure
%           -filenameOutput: output filename (if it has an image extension,
%           single files will be generated)
%           -tmo_alpha: \alpha_A, \alpha_B, \alpha_C coefficients
%           costants in the paper (Equation 3a, 3b, and 3c)
%           -tmo_dn_clamping: a boolean value (0 or 1) for setting black
%           and white levels clamping
%           -tmo_white: white point;
%           -tmo_gamma: gamma for encoding the frame
%           -tmo_quality: the output quality in [1,100]. 100 is the best quality
%           1 is the lowest quality.%
%           -tmo_video_profile: the compression profile (encoder) for compressing the stream.
%           Please have a look to the profile of VideoWriter from the MATLAB
%           help. Depending on the version of MATLAB some profiles may be not
%           be present.
%
%     Copyright (C) 2013-17 Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%     The paper describing this technique is:
%     "Real-time Automated Tone Mapping System for HDR Video"
% 	  by Chris Kiser, Erik Reinhard, Mike Tocci and Nora Tocci
%     in IEEE International Conference on Image Processing, 2012
%
%

if(~exist('tmo_alpha', 'var'))
    tmo_alpha = 0.98;
end

if(~exist('tmo_dn_clamping', 'var'))
    tmo_dn_clamping = 0;
end

if(~exist('tmo_white', 'var'))
    tmo_white = 1e6;
end

if(tmo_white <= 0.0)
    tmo_white = 1e6;
end

if(~exist('tmo_gamma', 'var'))
    tmo_gamma = -2.2;
end

if(tmo_gamma < 0)
    bsRGB = 1;
else
    bsRGB = 0;
end

if(~exist('tmo_quality', 'var'))
    tmo_quality = 95;
end

if(~exist('tmo_video_profile', 'var'))
    tmo_video_profile = 'MPEG-4';
end

name = RemoveExt(filenameOutput);
ext  = fileExtension(filenameOutput);

bVideo = 0;
writerObj = 0;

if(strcmp(ext, 'avi') == 1 | strcmp(ext, 'mp4') == 1)
    bVideo = 1;
    writerObj = VideoWriter(filenameOutput, tmo_video_profile);
    writerObj.FrameRate = hdrv.FrameRate;
    writerObj.Quality = tmo_quality;
    open(writerObj);
end

hdrv = hdrvopen(hdrv);

disp('Tone Mapping...');

beta_clamping   = 0.999;
beta_clamping_c = (1.0 - beta_clamping);

a_vec = [];
a_vec_flt = [];

maxLprev = 0;

mkdir(name);

for i=1:hdrv.totalFrames
    disp(['Processing frame ', num2str(i)]);
    [frame, hdrv] = hdrvGetFrame(hdrv, i);

    %only physical values
    frame = RemoveSpecials(frame);
    frame(frame < 0) = 0;

    if(tmo_dn_clamping)
        %clamp black and white levels
        L = RemoveSpecials(lum(frame));
        %compute CDF's histogram
        [histo, bound, ~] = HistogramHDR(L, 256, 'log10', [], 1);
        histo_cdf = cumsum(histo);
        histo_cdf = histo_cdf/max(histo_cdf(:));
        [~, ind] = min(abs(histo_cdf - beta_clamping));
        maxL = 10^(ind * (bound(2) - bound(1)) / 256 + bound(1));

        [~, ind] = min(abs(histo_cdf-beta_clamping_c));
        minL = 10^(ind * (bound(2) - bound(1)) / 256 + bound(1));

        frame(frame > maxL) = maxL;
        frame(frame < minL) = minL;
        frame = frame - minL;
    end

    %compute statistics for the current frame
    L = lum(frame);
    Lav = logMean(L);
    maxL = max(L(:));
    A = maxL - Lav;
    B = Lav - min(L(:));

    if(i == 1)
        maxLprev = maxL;
        Aprev = A;
        Bprev = B;
        aprev = 0.18 * 2^(2 * (B - A) / (A + B));
    end

    %leaky integration
    tmo_alpha_c = 1.0 - tmo_alpha;
    An = tmo_alpha_c * Aprev + tmo_alpha * A;
    Aprev = An;

    Bn = tmo_alpha_c * Bprev + tmo_alpha * B;
    Bprev = Bn;

    a = SceneKey(An, Bn);
    an = tmo_alpha_c * aprev + tmo_alpha * a;
    aprev = an;

    %tone mapping
    [frameOut, ~, ~] = ReinhardTMO(frame, an, tmo_white, 'global');

    %example
    a_vec = [a_vec, maxL];
    maxLn =  0.5 * maxLprev + 0.5 *maxL;
    maxLprev = maxLn;
    a_vec_flt = [a_vec_flt, maxLn ];

    %gamma/sRGB encoding
    if(bsRGB)
        frameOut = ClampImg(ConvertRGBtosRGB(frameOut, 0), 0, 1);
    else
        frameOut = ClampImg(GammaTMO(frameOut, tmo_gamma, 0, 0), 0, 1);
    end

    if(bVideo)
        writeVideo(writerObj, frameOut);
    else
        imwrite(frameOut, [name, '/frame_', sprintf('%.10d',i), '.', ext]);
    end

    %update statistics for the next frame
    Aprev = A;
    Bprev = B;
    aprev = a;
end
disp('OK');

if(bVideo)
    close(writerObj);
end

hdrvclose(hdrv);

end

function a = SceneKey(An, Bn)
    a = 0.18 * 2^(2 * (Bn - An) / (An + Bn));
end

author: banterle


text:
function [E0,V0,tau0,kaf,kas,eps,alpha] = BOLD_parameters()
% prior values for the Balloon HRF model

E0 = 0.34;
V0 = 4;%0.02;
tau0 = 2;%0.98;
kaf = 0.41;
kas = 0.65;
eps = 1;
alpha = 0.32;

author: MBB-team

repo: VBA-toolbox


text:
function chrpak_test065 ( )

%*****************************************************************************80
%
%% TEST065 tests I4_TO_UNARY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 June 2007
%
%  Author:
%
%    John Burkardt
%
  test_num = 3;

  i4_test = [ -5, 0, 7 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST065\n' );
  fprintf ( 1, '  I4_TO_UNARY converts an integer to unary.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        I4     UNARY\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    i4 = i4_test(test);

    s = i4_to_unary ( i4 );

    fprintf ( 1, '  %8d  %s\n', i4, s );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Quan Wang, Kim L. Boyer,
% The active geometric shape model: A new robust deformable shape model and its applications,
% Computer Vision and Image Understanding, Volume 116, Issue 12, December 2012, Pages 1178-1194,
% ISSN 1077-3142, 10.1016/j.cviu.2012.08.004.
%
% For commercial use, please contact the authors.

function [xc yc a b phi]=fit_arb_ellipse_force(init,increment,threshold,bound,field_x,field_y,iter,show_fitness)

x0=init(1);
y0=init(2);
a0=init(3);
b0=init(4);
phi0=init(5);
d_xc=increment(1);
d_yc=increment(2);
d_a=increment(3);
d_b=increment(4);
d_phi=increment(5);
t_xc=threshold(1);
t_yc=threshold(2);
t_a=threshold(3);
t_b=threshold(4);
t_phi=threshold(5);
a_low=bound(1);
a_up=bound(2);
b_low=bound(3);
b_up=bound(4);

%%  Fit an arbitrary ellipse in the force filed. xc, yc, a, b and phi will be returned.
%   The parametric equations of the arbitrary ellipse:
%       x=xc+a*cos(theta)*cos(phi)-b*sin(theta)*sin(phi)
%       y=yc+a*cos(theta)*sin(phi)+b*sin(theta)*cos(phi)
%   x0, y0: initial position of center
%   a0, b0: initial semi-major and semi-minor axes
%   phi0: initial orientation
%   d_xc, d_yc: increment of xc and yc in each loop
%   d_a, d_b: increment of a and b in each loop
%   d_phi: increment of phi in each loop
%   t_xc, t_yc: threshold needed to update xc and yc
%   t_a, t_b: threshold needed to update a and b
%   t_phi: threshold needed to update phi
%   a_low, a_up: lower bound and upper bound of a
%   b_low, b_up: lower bound and upper bound of b
%   field_x: the x component of force field
%   field_y: the y component of force field
%   iter: number of iterations
%   show_fitness: a flag of whether to show the fitness function

% note: m rows and n columns, but x is column and y is row here
[m n]=size(field_x);

xc=x0;
yc=y0;
a=a0;
b=b0;
phi=phi0;

for it=1:iter
    [x,y,theta]=arb_ellipse_in_image(m,n,xc,yc,a,b,phi);
    N=max(size(x));

    %% torque along the ellpise about center
    torque_v=[0 0 0]; % the torque as a vector
    for i=1:max(size(theta))
        torque_v=torque_v+cross([x(i)-xc,y(i)-yc,0],...
            [field_x(y(i),x(i)),field_y(y(i),x(i)),0]);
    end
    torque=torque_v(3);
    torque=torque/N^2;
    if torque>t_phi
        phi=phi+d_phi;
    elseif torque<-t_phi
        phi=phi-d_phi;
    end

    %% F_around
    F_round=[0,0];
    for i=1:max(size(theta))
        F_round=F_round+[field_x(y(i),x(i)),field_y(y(i),x(i))];
    end
    F_round=F_round/max(size(theta));

    %% F_left, which is on the left quarter of the ellipse, and inward
    index = find( theta>pi*3/4 & theta<pi*5/4 );
    F_left=0;
    for i=index
        F_left=F_left+dot([field_x(y(i),x(i)),field_y(y(i),x(i))],[cos(phi),sin(phi)]);
    end
    F_left=F_left/max(size(index));

    %% F_right, which is on the right quarter of the ellipse, and inward
    index = find( theta<pi/4 | theta>pi*7/4 );
    F_right=0;
    for i=index
        F_right=F_right+dot([field_x(y(i),x(i)),field_y(y(i),x(i))],[-cos(phi),-sin(phi)]);
    end
    F_right=F_right/max(size(index));

    %% F_up, which is on the up quarter of the ellipse, and inward
    index = find( theta>pi/4 & theta<pi*3/4 );
    F_up=0;
    for i=index
        F_up=F_up+dot([field_x(y(i),x(i)),field_y(y(i),x(i))],[sin(phi),-cos(phi)]);
    end
    F_up=F_up/max(size(index));

    %% F_down, which is on the low quarter of the ellipse, and inward
    index = find( theta>pi*5/4 & theta<pi*7/4 );
    F_down=0;
    for i=index
        F_down=F_down+dot([field_x(y(i),x(i)),field_y(y(i),x(i))],[-sin(phi),cos(phi)]);
    end
    F_down=F_down/max(size(index));

    %% update xc and yc
    F_left_right=dot(F_round,[1,0]);
    if F_left_right>t_xc
        xc=xc+d_xc;
    elseif F_left_right<-t_xc
        xc=xc-d_xc;
    end

    F_down_up=dot(F_round,[0,1]);
    if F_down_up>t_yc
        yc=yc+d_yc;
    elseif F_down_up<-t_yc
        yc=yc-d_yc;
    end

    %% update xc and yc again according to diagonal force
    F_diag1=dot(F_round,[0.7071,0.7071]);
    if F_diag1>t_xc+t_yc
        xc=xc+d_xc;
        yc=yc+d_yc;
    elseif F_diag1<-t_xc-t_yc
        xc=xc-d_xc;
        yc=yc-d_yc;
    end

    F_diag2=dot(F_round,[-0.7071,0.7071]);
    if F_diag2>t_xc+t_yc
        xc=xc-d_xc;
        yc=yc+d_yc;
    elseif F_diag2<-t_xc-t_yc
        xc=xc+d_xc;
        yc=yc-d_yc;
    end

    %% update a and b

    if F_left+F_right>t_a
        a=a-d_a;
    elseif F_left+F_right<-t_a
        a=a+d_a;
    end

    if F_up+F_down>t_b
        b=b-d_b;
    elseif F_up+F_down<-t_b
        b=b+d_b;
    end

    if b>a
        temp=a;a=b;b=temp;
        phi=mod(phi+pi/2,pi);
    end

    %% restrict a and b using lower and upper bounds
    if a>a_up
        a=a_up;
    end
    if a<a_low
        a=a_low;
    end
    if b>b_up
        b=b_up;
    end
    if b<b_low
        b=b_low;
    end

    %% fitness function
    if show_fitness==1
        beta=0.9;
        fit0=0;
        fit1=0;
        fit2=0;
        [x,y,theta]=arb_ellipse_in_image(m,n,xc,yc,a,b,phi);
        [x1,y1,theta1]=arb_ellipse_in_image(m,n,xc,yc,a*beta,b*beta,phi);
        [x2,y2,theta2]=arb_ellipse_in_image(m,n,xc,yc,a/beta,b/beta,phi);
        for i=1:max(size(theta))
            fit0=fit0+norm([field_x(y(i),x(i)),field_y(y(i),x(i))]);
        end
        for i=1:max(size(theta1))
            fit1=fit1+norm([field_x(y1(i),x1(i)),field_y(y1(i),x1(i))]);
        end
        for i=1:max(size(theta2))
            fit2=fit2+norm([field_x(y2(i),x2(i)),field_y(y2(i),x2(i))]);
        end
        fit0=fit0/max(size(theta));
        fit1=fit1/max(size(theta1));
        fit2=fit2/max(size(theta2));
        fit_save(it)=fit0-fit1/2-fit2/2;
    end

end

if show_fitness==1
    figure;hold on;
    plot(1:iter,fit_save,'b','LineWidth',2);
    legend('fitness function');
    grid on;
    xlabel('iteration');
    ylabel('fitness function');
    title('fitness function in each iteration');
end

author: Sable

repo: mcbench-benchmarks


text:
function sys = optimizer(Constraints,Objective,options,x,u)
%OPTIMIZER  Container for optimization problem
%
%   OPT = OPTIMIZER(Constraints,Objective,options,x,u) exports an object that
%   contains precompiled numerical data to be solved for varying arguments
%   x, returning the optimal value of the expression u.
%
%   OPTIMIZER works most efficiently if the varying data x enters the
%   optmization problem affinely. For the general case, much more logic has
%   to be applied when instantiating the numerical data for a parametric
%   value, and when compiling the model, it is harder fpr YALMIP to
%   understand what kind of model it will be once the parameters are fixed.
%
%   By default, display is turned off (since optimizer is used in
%   situations where many problems are solved repeatedly. To turn on
%   display, set the verbose option in sdpsetting to 2.
%
%   Example
%
%    The following problem creates an LP with varying upper and lower
%    bounds on the decision variable.
%
%     A = randn(10,3);
%     b = rand(10,1)*19;
%     c = randn(3,1);
%
%     z = sdpvar(3,1);
%     sdpvar UB LB
%
%     Constraints = [A*z <= b, LB <= z <= UB];
%     Objective = c'*z
%     % We want the optimal z as a function of [LB;UB]
%     optZ = optimizer(Constraints,Objective,[],[LB; UB],z);
%
%     % Compute the optimal z when LB=1, UB = 3;
%     zopt = optZ([1; 3])
%
%     % Compute two solutions, one for (LB,UB) [1;3] and one for (LB,UB) [2;6]
%     zopt = optZ([[1; 3], [2;6]])
%
%     % A second output argument can be used to catch infeasibility
%     [zopt,infeasible] = optZ([1; 3])
%
%     % To avoid the need to vectorize in order to handle multiple
%     parameters, a cell-based definition can be used
%
%     optZ = optimizer(Constraints,Objective,[],{LB,UB},{z,sum(z)})
%     [zopt,infeasible] = optZ({1,3});
%     zopt{1}
%     zopt{2}

if nargin < 5
    error('OPTIMIZER requires 5 inputs');
end

% With the new optional cell-based format, the internal format is always a
% vector with all information stacked, both in and out. Hence, we need to
% save original sizes before stacking things up
xoriginal = x;
if ~isa(x,'cell')
    xoriginal = {x};
end
if isa(x,'cell')
    xvec = [];
  for i = 1:length(x)
      if ~(isa(x{i},'sdpvar') | isa(x{i},'ndsdpvar'))
          error(['The parameters must be SDPVAR objects. Parameter #' num2str(i) ' is a ' upper(class(x{i}))]);
      end
      if is(x{i},'complex')
          x{i} = [real(x{i});imag(x{i})];
          complexInput(i) = 1;
      else
          complexInput(i) = 0;
      end
      sizeOrigIn{i} = size(x{i});
      z = x{i}(:);
      mask{i} = uniqueNonZeroRows(z);
      xvec = [xvec;z(mask{i})];
  end
  x = xvec;
else
    if ~isreal(x)%,'complex')
        complexInput(1) = 1;
        x = [real(x);imag(x)];
    else
        complexInput(1) = 0;
    end
    sizeOrigIn{1} = size(x);
    x = x(:);
    mask{1} = uniqueNonZeroRows(x);
    x = x(mask{1});
end
nIn = length(x);
mIn = 1;

if isa(x,'sdpvar')
    if ~is(x,'lpcone')
    error('All parameter arguments have to be simple variables (i.e., not expressions such a+b or 1+a)');
    end
end

if isa(u,'cell')
    uvec = [];
    for i = 1:length(u)
        if is(u{i},'complex')
            complexOutput(i) = 1;
            u{i} = [real(u{i});imag(u{i})];
        else
            complexOutput(i) = 0;
        end
        sizeOrigOut{i} = size(u{i});
        uvec = [uvec;u{i}(:)];
    end
    u = uvec;
else
    if (isa(u,'sdpvar') || isa(u,'ndsdpvar')) && is(u,'complex')
        complexOutput(1) = 1;
        u = [real(u);imag(u)];
    else
        complexOutput(1) = 0;
    end
    sizeOrigOut{1} = size(u);
    u = u(:);
end
nOut = length(u);
mOut = 1;

if isempty(options)
    options = sdpsettings;
end

if ~isa(options,'struct')
    error('Third argument in OPTIMIZER should be an options structure.');
end

% Silent by default. If we want displays, set to 2
options.verbose = max(options.verbose-1,0);

% Since code is based on a fake equality, we must avoid bound propagation
% based on equalities
options.avoidequalitybounds=1;

% Normalize...
if isa(Constraints,'constraint')
    Constraints = lmi(Constraints);
end

if ~isempty(Constraints)
    if ~isa(Constraints,'constraint') &  ~isa(Constraints,'lmi')
        error('The first argument in OPTIMIZER should be a set of constraints');
    end
end

if ~isempty(Constraints)
    if any(is(Constraints,'sos'))
        tempOps = options;
        tempOps.sos.model = 2;
        tempOps.verbose = max(0,tempOps.verbose-1);
        parameter_sos = [x;u;recover(depends(Objective))];
        parameter_sos = depends(parameter_sos);
        for i = 1:length(Constraints)
            if ~is(Constraints,'sos')
                parameter_sos = [parameter_sos depends(Constraints(i))];
            end
        end
        parameter_sos = recover(parameter_sos);
        [Constraints,Objective] = compilesos(Constraints,Objective,tempOps,parameter_sos);
    end
end

if ~isequal(options.solver,'')
    % User has specified solver. Let us impose this solver forcefully to
    % the compilation code, in order to handle nonlinear parameterizations
    if ~strcmp(options.solver(1),'+')
        options.solver = ['+' options.solver];
    end
end

if options.removeequalities
    error('''removeequalities'' in optimizer objects not allowed.');
end

if ~isempty(Constraints) & any(is(Constraints,'uncertain'))
    [Constraints,Objective,failure] = robustify(Constraints,Objective,options);
    [aux1,aux2,aux3,model] = export((x == repmat(pi,nIn*mIn,1))+Constraints,Objective,options,[],[],0);
else
    [aux1,aux2,aux3,model] = export((x == repmat(pi,nIn*mIn,1))+Constraints,Objective,options,[],[],0);
end

if ~isempty(aux3)
    if isstruct(aux3)
        if ismember(aux3.problem, [-9 -5 -4 -3 -2 -1 14])
            error(['Failed exporting the model: ' aux3.info])
        end
    end
end

if norm(model.F_struc(1:nIn*mIn,1)-repmat(pi,length(x),1),inf) > 1e-10
    error('Failed exporting the model (try to specify another solver)')
end

% Try to set up an optimal way to compute the output
base = getbase(u);
if isempty(u) || (is(u,'linear') & all(sum(base | base,2) == 1) & all(sum(base,2)==1) & all(base(:,1)==0))
    % This is just a vector of variables
    z = [];
    map = [];
    uvec = u(:);
    % Setup to do fast getvariables(uvec(i))
    U = getbase(uvec);
    Uvar = getvariables(uvec);
    U = U(:,2:end);
    [ii,jj,ss] = find(U');
    for i = 1:length(uvec)
        %var = getvariables(uvec(i)); Slow
        %var = Uvar(find(U(i,:)));
        var = Uvar(ii(i));
        mapIndex = find(var == model.used_variables);
        if ~isempty(mapIndex)
            map = [map;mapIndex];
        else
            map = [map;0];
        end
    end
else
    % Some expression which we will use assign and double to evaluate
    vars = depends(u);
    z = recover(vars);
    map = [];
    for i = 1:length(z)
        var = vars(i);
        mapIndex = find(var == model.used_variables);
        if ~isempty(mapIndex)
            map = [map;mapIndex];
        else
            map = [map;0];
        end
    end
end

if ~isempty(u) && (isempty(map) | min(size(map))==0)
    error('The requested decision variable (argument 4) is not in model');
end

model.getsolvertime = 0;

model.solver.callhandle = str2func(model.solver.call);

model.options = pruneOptions(model.options);
model.hashCache = gen_rand_hash(0,size(model.monomtable,2),1);

sys.recover = aux2;
sys.model = model;
sys.dimin = [nIn mIn];
sys.dimout = [nOut mOut];
sys.diminOrig = sizeOrigIn;
sys.dimoutOrig = sizeOrigOut;
sys.complexInput = complexInput;
sys.complexOutput = complexOutput;
sys.mask = mask;
sys.map = map;
sys.input.xoriginal = xoriginal;
sys.input.expression = x;
sys.output.expression = u;
sys.output.z = z;
sys.lastsolution = [];
sys.ParametricSolution = [];
sys.model.infeasible = 0;
% This is not guaranteed to give the index in the order the variables where
% given (tested in test_optimizer2
% [a,b,c] = find(sys.model.F_struc(1:prod(sys.dimin),2:end));
% Could be done using
[b,a,c] = find(sys.model.F_struc(1:prod(sys.dimin),2:end)');
% but let us be safe
%b = [];
%for i = 1:prod(sys.dimin)
%    b = [b;find(sys.model.F_struc(i,2:end))];
%end
sys.model.parameterIndex = b;
used_in = find(any(sys.model.monomtable(:,b),2));
Q = sys.model.Q;
Qa = Q;Qa(:,b)=[];Qa(b,:)=[];
Qb = Q(:,b);Qb(b,:)=[];
if nnz(Q)>0
    zeroRow = find(~any(Q,1));
    Qtest = Q;Q(zeroRow,:)=[];Q(:,zeroRow)=[];
    problematicQP = nonconvexQuadratic(Qtest);%min(eig(full(Qtest)))<-1e-14;
else
    problematicQP = 0;
end
if  any(sum(sys.model.monomtable(used_in,:),2)>1) | any(sum(sys.model.monomtable(used_in,:) | sys.model.monomtable(used_in,:),2) > 1) | problematicQP | ~isempty(sys.model.evalMap) | any(any(sys.model.monomtable<0))
    sys.nonlinear = 1;
else
    sys.nonlinear = 0;
end
sys.F = Constraints;
sys.h = Objective;
sys.ops = options;

sys.complicatedEvalMap = 0;
% Are all nonlinear operators acting on simple parameters? Elimination
% strategy will only be applied on simple problems such as x<=exp(par)
if ~model.solver.evaluation
    for i = 1:length(sys.model.evalMap)
        if ~all(ismember(sys.model.evalMap{i}.variableIndex,sys.model.parameterIndex))
            sys.complicatedEvalMap = 1;
        end
        if length(sys.model.evalMap{i}.arg)>2
            sys.complicatedEvalMap = 1;
        end
    end
    if sys.complicatedEvalMap
        error('Parameters are currently only allowed to enter function such as exp, sin etc as exp(a), sin(b) etc.')
    end
end

sys.model.evalParameters = [];
if sys.nonlinear
    % These artificial equalities are removed if we will use eliminate variables

    %   sys.model.F_struc(1:length(sys.parameters),:) = [];
    %   sys.model.K.f = sys.model.K.f - length(sys.parameters);

    % Which variables are simple nonlinear operators acting on parameters
    evalParameters = [];
    for i = 1:length(sys.model.evalMap)
        if all(ismember(sys.model.evalMap{i}.variableIndex,sys.model.parameterIndex))
            evalParameters = [evalParameters;sys.model.evalMap{i}.computes(:)];
        end
    end
    sys.model.evalParameters = evalParameters;
end

% In case we perform partial instantiation, we have to remember where we
% came from originallty when finally solving problems
sys.instatiatedvalues = zeros(length(model.used_variables),1);
sys.orginal_usedvariables = sys.model.used_variables;
sys.orginal_parameterIndex = sys.model.parameterIndex;

sys.input.stochastics = cell(1,length(sys.diminOrig));
if ~isempty(Constraints)
    randDefinitions = find(is(Constraints,'random'));
    if ~isempty(randDefinitions)

        for i = 1:length(randDefinitions)
            Fi = Constraints(randDefinitions(i));
            randDef{i}.distribution = struct(struct(Fi).clauses{1}.data).extra.distribution;
            randDef{i}.variables = sdpvar(Fi);

            for j = 1:length(sys.diminOrig)
                if isequal(getbase(sys.input.xoriginal{j}),getbase(randDef{i}.variables)) && isequal(getvariables(sys.input.xoriginal{j}),getvariables(randDef{i}.variables))
                    sys.input.stochastics{j} = randDef{i}.distribution;
                end
            end
        end
    end
end

% Remove place holder constraints. No longer used
sys.model.F_struc(1:prod(sys.dimin),:)=[];
sys.model.K.f = sys.model.K.f-prod(sys.dimin);

sys = class(sys,'optimizer');
sys = optimizer_precalc(sys);

function i = uniqueNonZeroRows(x);
B = getbase(x);
% Quick check for trivially unique rows, typical 99% case
[n,m] = size(B);
if n == m-1 && nnz(B)==n
    if isequal(B,[spalloc(n,1,0) speye(n)])
        i = 1:n;
        return
    end
end
if  length(unique(B*randn(size(B,2),1))) == n && nnz(B*randn(size(B,2),1)) == size(B,1)
    i = 1:n;
    return
end
[temp,i] = unique(B*randn(size(B,2),1));
z = find(~any(B,2));
i = setdiff(i,z);
i = i(:);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


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

%Empty workspace and close figures
close all;
clear;

%Define the range of BS antennas
Mvalues = [1 10 100];

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

%Nominal angle of the desired UE
thetaDesired = pi/6;

%Range of nominal angles of the interfering UE
varphiInterfererDegrees = -180:1:180;
varphiInterfererRadians = varphiInterfererDegrees*(pi/180);

%Define the antenna spacing (in number of wavelengths)
antennaSpacing = 1/2; %Half wavelength distance

%Define the range of the effective SNR in (3.13) for the desired UE
SNR1dB = 10;
SNR1 = 10.^(SNR1dB/10);

%Define the range of the effective SNR in (3.13) for the interfering UE
SNR2dB = SNR1dB-10;
SNR2 = 10.^(SNR2dB/10);

%Preallocate matrices for storing the simulation results
relativeInterferenceCoherent = zeros(length(varphiInterfererRadians),length(thetaDesired),length(Mvalues));

%Compute correlation matrix of the desired UE
R1 = functionRlocalscattering(max(Mvalues),thetaDesired,ASDdeg,antennaSpacing);

%% Go through all angles of the interfering UE
for n = 1:length(varphiInterfererRadians)

    %Compute correlation matrix of the interfering UE
    R2 = functionRlocalscattering(max(Mvalues),varphiInterfererRadians(n),ASDdeg,antennaSpacing);

    %Go through all numbers of antennas
    for m = 1:length(Mvalues)

        R1m = R1(1:Mvalues(m),1:Mvalues(m));
        R2m = R2(1:Mvalues(m),1:Mvalues(m));

        %Compute numerator and denominator of the ratio between desired
        %signal term and coherent interference term in (4.17)
        numerator = SNR2^2*abs(trace(R1m*((SNR1*R1m+SNR2*R2m+eye(Mvalues(m)))\R2m)))^2;
        denominator = SNR1^2*abs(trace(R1m*((SNR1*R1m+SNR2*R2m+eye(Mvalues(m)))\R1m)))^2;

        relativeInterferenceCoherent(n,m) = numerator/denominator;

    end

end

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

plot(varphiInterfererDegrees,10*log10(relativeInterferenceCoherent(:,1)),'k-','LineWidth',1);
plot(varphiInterfererDegrees,10*log10(relativeInterferenceCoherent(:,2)),'r--','LineWidth',1);
plot(varphiInterfererDegrees,10*log10(relativeInterferenceCoherent(:,3)),'b-.','LineWidth',1);

xlabel('Angle of interfering UE [degree]');
ylabel('Coherent interf. power over signal power [dB]');
xlim([-180 180]);
ylim([-40 -10]);

legend('M=1','M=10','M=100','Location','NorthWest');

author: emilbjornson

repo: massivemimobook


text:
function title = p32_title ( )

%*****************************************************************************80
%
%% P32_TITLE returns the title of problem P32
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 32, The Lorenz System';

  return
end


text:
function [finer_approximation]=gsp_pyramid_synthesis_single(G,coarse_approximation,prediction_error,keep_inds,param)
%GSP_PYRAMID_SYNTHESIS_SINGLE Synthesize a single level of the graph pyramid transform
%   Usage:  [finer_approximation]=gsp_pyramid_synthesis_single(G,coarse_approximation,prediction_error,keep_inds,param);
%
%   Input parameters:
%         G                    : Graph structure on which the signal resides.
%         coarse_approximation : Coarse approximation on a reduced graph.
%         prediction_error     : Prediction error that was made when forming the coarse approximation.
%         keep_inds            : Indices of the vertices to keep when downsampling the graph and signal.
%         param                : Additional parameters.
%   Output parameters:
%         finer_approximation  : Coarse approximation of the signal on a higher resolution graph.
%
%   'gsp_pyramid_synthesis_single(G,coarse_approximation,prediction_error,keep_inds,param)'
%   synthesizes a single level of the graph pyramid transform.
%
%   *param* is a structure containing optional arguments including
%
%   * *param.regularize_epsilon* : Interpolation parameter.
%   * *param.least_squares* : Set to 1 to use the least squares synthesis
%     (default=0)
%   * *param.use_landweber* : Set to 1 to use the Landweber iteration
%     approximation in the least squares synthesis.
%   * *param.landweber_its* : Number of iterations in the Landweber
%     approximation for least squares synthesis.
%   * *param.landweber_tau* : Parameter for the Landweber iteration.
%   * *param.h_filter* : A graph spectral filter. This filter is
%     required for least squares synthesis, but not for the direct
%     synthesis method
%
%   Please read the documentation of |gsp_filter_analysis| for other
%   optional arguments.
%
%   See also: gsp_graph_multiresolution gsp_pyramid_synthesis
%             gsp_pyramid_analysis gsp_pyramid_analysis_single
%             gsp_pyramid_cell2coeff
%
%   Demo: gsp_demo_pyramid
%
%   References: shuman2013framework

%   Author : David I Shuman, Nathanael Perraudin.
%   Date   : 26 November 2015
%   Testing: test_pyramid

if ~isfield(param, 'least_squares')
    param.least_squares=0;
end

if ~isfield(G,'lmax')
    G=gsp_estimate_lmax(G);
end

% Compute the next coarse approximation at a higher resolution level
if ~param.least_squares % i.e., use the direct synthesis method
%     upsampled_coarse_approximation=zeros(G.N,size());
%     upsampled_coarse_approximation(keep_inds,:)=coarse_approximation;

    prediction=gsp_interpolate(G,coarse_approximation,keep_inds,param);
    finer_approximation=prediction+prediction_error;
else
    if ~isfield(param, 'h_filter')
        error('h-filter not provided');
    end
    if ~isfield(param,'use_landweber')
        if (G.N>3000 || ~isfield(G,'e') || ~isfield(G,'U') )
            param.use_landweber=1;
        else
            param.use_landweber=0;
        end
    end
    Nbar=length(keep_inds);
    if ~isfield(param, 'regularize_epsilon')
        param.regularize_epsilon=.005;
    end
    if param.use_landweber
        if ~isfield(param, 'landweber_its')
            param.landweber_its=50;
        end
        if ~isfield(param, 'landweber_tau')
            param.landweber_tau=1;
        end
        % This is the Landweber iteration to apply the pseudoinverse
        % (the efficiency of this set of operations can probably be improved)
        x_old=zeros(G.N,1);
        z=[coarse_approximation;prediction_error];
        PhiV1t=zeros(Nbar,G.N);
        green_kernel=@(x) 1./(x+param.regularize_epsilon);
        for i=1:Nbar
            PhiV1t(i,:)=(gsp_filter(G,green_kernel,gsp_delta(G.N,keep_inds(i)),param))';
        end
        for k=1:param.landweber_its
            [x_bar,y_bar]=gsp_pyramid_analysis_single(G,x_old,keep_inds,param.h_filter,param);
            z_prime=[x_bar;y_bar];
            z_delt=z-z_prime;
            alpha_new=PhiV1t*z_delt((Nbar+1):end);
            x_up=zeros(G.N,1);
            x_up(keep_inds)=z_delt(1:Nbar);
            regularized_L= G.L+param.regularize_epsilon*eye(G.N);
            elim_inds=setdiff(1:G.N,keep_inds);
            next_term=regularized_L(keep_inds,keep_inds)*alpha_new-regularized_L(keep_inds,elim_inds)*(regularized_L(elim_inds,elim_inds)\(regularized_L(elim_inds,keep_inds)*alpha_new));
            next_up=zeros(G.N,1);
            next_up(keep_inds)=next_term;
            x_new=x_old+param.landweber_tau*(gsp_filter(G,param.h_filter,x_up-next_up,param)+z_delt((Nbar+1):end));
            x_old=x_new;
        end
        finer_approximation=x_new;
    else % When the graph is small enough, we can do a full eigendecomposition and compute the full analysis operator T_a
        if ( ~isfield(G,'e') || ~isfield(G,'U') )
            G=gsp_compute_fourier_basis(G);
        end
        H=G.U*diag(param.h_filter(G.e))*G.U';
        Phi=G.U*diag(1./(param.regularize_epsilon+G.e))*G.U';
        S=sparse(1:Nbar,keep_inds,ones(Nbar,1),Nbar,G.N);
        Ta=[S*H;eye(G.N)-Phi(:,keep_inds)*(Phi(keep_inds,keep_inds)\(S*H))];
        finer_approximation=(Ta'*Ta)\(Ta'*[coarse_approximation;prediction_error]);
    end
end

end

repo: gspbox


text:
function [tree,asgn] = vl_hikmeans(data,K,nleaves)
% VL_HIKMEANS  Hierachical integer K-means
%   [TREE,ASGN] = VL_HIKMEANS(DATA,K,NLEAVES) applies integer K-menas
%   recursively to cluster the data DATA, returing a structure TREE
%   representing the clusters and a vector ASGN with the data to
%   cluster assignments. The depth of the recursive partition is
%   computed so that at least NLEAVES are generated.
%
%   VL_HIKMEANS() is built on top of VL_IKMEANS() and requires the
%   data to be of class UINT8.
%
%   TREE is a structure representing the hierarchical clusters.  Each
%   node of the tree is also a structure with fields:
%
%   DEPTH::
%     Depth of the tree (only at the root node)
%
%   CENTERS::
%     K cluster centers
%
%   SUB::
%     Array of K node structures representing subtrees
%     (this field is missing at leaves).
%
%   ASGN is a matrix with one column per datum and height equal to the
%   depth of the tree. Each column encodes the branch of the tree that
%   correspond to each datum.
%
%   Example::
%     ASGN(:,7) = [1 5 3] means that the tree as depth equal to 3 and
%     that the datum X(:,7) corresponds to the branch
%     ROOT->SUB(1)->SUB(5)->SUB(3).
%
%   See also: VL_HIKMEANSPUSH(), VL_HIKMEANSHIST(), VL_IKMEANS(), VL_HELP().

% Copyright (C) 2014 Andrea Vedaldi.
% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function bvhModify(handle, channels, skel, padding)

% BVHMODIFY Helper code for visualisation of bvh data.
% FORMAT
% DESC updates a bvh representation in a 3-D plot.
% ARG handle : a vector of handles to the structure to be updated.
% ARG channels : the channels to update the skeleton with.
% ARG skel : the skeleton structure.
%
% SEEALSO : skelModify, bvhVisualise
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% MOCAP

if nargin<4
  padding = 0;
end

skelModify(handle, channels, skel, padding);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef TestXYFunctions < matlab.uitest.TestCase
    properties(TestParameter)
        classToTest = getclasses('XYfun');
    end

    properties
        zap         =  getZmapAnalysisPkg();
    end

    methods(Test)
        function testInitializable(testCase, classToTest)
            disp(classToTest);
            myclass = str2func(classToTest);
            obj = myclass(testCase.zap,'DelayProcessing',true,'InteractiveMode',false,'AutoShowPlots',false);
            testCase.assumeClass(obj,myclass);
            testCase.assertEqual(obj.RawCatalog,    zap.Catalog);
            testCase.assertEqual(obj.EventSelector, zap.EventSel);
            testCase.assertEqual(obj.Grid,          zap.Grid);
            testCase.assertEqual(obj.Shape,         zap.Shape);
            testCase.assertTrue(any( ReturnDetails.Properties.VariableNames == obj.active_col ));

            % XYfun.objectHandle(abc)
        end

        function testCalculatable(testCase, classToTest)
            disp(classToTest);
            myclass = str2func(classToTest);
            obj = myclass(testCase.zap,'DelayProcessing',true,'InteractiveMode',false,'AutoShowPlots',false);
            results = obj.Calculate();
            testCase.assumeNotEmpty(results);

        end
    end

end

function zap = getZmapAnalysisPkg()

    % get catalog
    catalogFile = 'resrc/sample/SED_fdsn_2000_on.mat';
    c=load(catalogFile,'catalog');

    % set up sampling parameters
    num_events  = 200;
    max_radius = 50;
    dist_units = 'km';
    evsel = EventSelectionParameters('NumClosestEventsUpToRadius', num_events, max_radius,'DistanceUnits',dist_units);

    % set up grid
    g_opts.dx = 10;
    g_opts.dy = 10;
    g_opts.dz = 10;
    g_opts.xyunits = 'kilometer';
    g_opts.FollowMeridians = false;
    g_opts.GridEntireArea = false;

    %center on ETHZ
    fixedOpts.UseFixedAnchorPoint = true;
    fixedOpts.XAnchor = 47.3763;
    fixedOpts.YAnchor = 8.5477;
    fixedOpts.ZAnchor = 0;

    gopt= GridOptions(g_opts, fixedOpts);
    g = ZmapGrid('testgrid', gopt);

    % set up shape
    sh = load_shape('eq_data/switzerland_shape.csv');

    zap = ZmapAnalysisPkg([],c,evsel,g,sh);

    %zap = ZmapAnalysisPkg.fromGlobal('primeCatalog');

end

function c = getfuns(packageName)
    a=what(packageName);
    a.m(a.m == "Contents.m") = [];
    funNames=replace(a.m,'.m','');
    funNames = strcat(packageName, '.', funNames);
    c = cellfun(@str2func,funNames,'UniformOutput',false);
end

function c = getclasses(packageName)
    a=what(packageName);
    a.m(a.m == "Contents.m") = [];
    funNames=replace(a.m,'.m','');
    c = strcat(packageName, '.', funNames);
end

author: CelsoReyes


text:
% Residue Curve Map for Reactive Systems / Methyl Acetate Chemistry
% Author's Data: Housam BINOUS
% Department of Chemical Engineering
% National Institute of Applied Sciences and Technology
% Tunis, TUNISIA
% Email: binoushousam@yahoo.com

c=rand(50,3);

for i=0:9,

xO1=0.00867394;
x02=0.608674;
x03=0.191326;
x04=0.191326;

T0=330;
y01=0.00025;

X01=0.1*i;
X02=0.1*(9-i);

tf=20;

x0 = [xO1  x02  x03  x04 y01 T0 X01 X02];

opts = odeset('Mass','M','MassSingular','yes');

[t,x] = ode15s('RCM_MethylAcetate',[0 tf],x0,opts);

x1=x(:,7);
x2=x(:,8);

figure(1);
x1=x(:,7);
x2=x(:,8);

AXIS([0 1 0 1])
hold on

plot(x1,x2,'color',[c(i+1,1) c(i+1,2) c(i+1,3)],'LineWidth',2);

end

for i=0:9,

xO1=0.00867394;
x02=0.608674;
x03=0.191326;
x04=0.191326;

T0=330;
y01=0.00025;

X01=0.1*i;
X02=0.1*(9-i);

tf=20;

x0 = [xO1  x02  x03  x04 y01 T0 X01 X02];

opts = odeset('Mass','M','MassSingular','yes');

[t,x] = ode15s('RCM_MethylAcetate2',[0 tf],x0,opts);

x1=x(:,7);
x2=x(:,8);

figure(1);
x1=x(:,7);
x2=x(:,8);

AXIS([0 1 0 1])
hold on

plot(x1,x2,'color',[c(i+1,1) c(i+1,2) c(i+1,3)],'LineWidth',2);

end


text:
classdef Kernel
    % Kernel
    % Represents the kernel G(x)
    properties (SetAccess = protected, GetAccess = public)
        func, % anonymous function representing G(x)
        der, % anonymous function representing G'(x)
        scalFunc, % anonymous function,
        % scalFunc(a,b,rho) =2 \pi C(rho) \int_{a}^b rG'(r) J_0(\rho r)
        % where C(rho) is the normalization constant, si Cp.m in folder
        % radialQuad.
        normFunc,  % anonymous function,
        % normFunc(a,b) = \sqrt{2\pi \int_{a}^b r G'(r)^2 dr}
        startFreq = 0, % real number
        % frequency where the greatest Bessel coefficients are located.
        gamma_est, % Estimator used to predict the number of coefficients in
        % the Bessel series
        customRadialQuad = [],
        % If non-empty, contains the method to compute the radial
        % quadrature optimized for this kernel.
        singular = true;
        % Set to false if x G'(x)^2 is integrable near 0.
        lim0;
    end

    methods
        % Class constructor
        function[kernel] = Kernel(func,der)
            if nargin == 0
                func = @(x)(0*x);
                der = @(x)(0*x);
            end
            kernel.func = func;
            kernel.der = der;
            fun1 = @(rho)(@(x)((x.*der(x)).*...
                -rho(:)'.*Cp(rho(:)').*besselj(1,rho(:)'*x)));
            % H10 scalar product
            fun2 = @(x)(x.*der(x).^2);
            % H10 norm
            kernel.scalFunc = @(a,b,rho)(2*pi*integral(fun1(rho),...
                a,b,'ArrayValued',true)'); % -2\pi \int_{a}^b rf'(r)ep'(r)dr
            kernel.normFunc = @(a,b)(sqrt(2*pi*integral(fun2,a,b)));
            kernel.startFreq = 0; % We don't know a priori where the energy
            % is located in the spectrum.
            kernel.gamma_est = @(tol)deal(0,7); % No fine tuning of gamma
            kernel.customRadialQuad = [];
            if ~ or(isinf(func(0)),isnan(func(0)))
                kernel.singular = false;
                kernel.gamma_est = @(tol)deal(0,3); % No fine tuning of gamma
                kernel.lim0 = func(0);
            else
                kernel.lim0 = 0;
            end
        end
    end

    methods (Access = public)

        % Efficient Bessel Decomposition for the kernel.
        function[onlineEBD,rq,loc] = offlineEBD(this,X,Y,a,tol)

            % Rescaling
            rMax = rMaxCalc(X,Y); % diameter
            x = X/rMax; % in B(0,1)
            y = Y/rMax; % in B(0,1)
            k1 = this.dilatation(rMax); % k1 must be approximated on {a < r < 1}

            % Radial quadrature (Bessel decomposition in 1D)
            rq = k1.radialQuadKernel(a,tol);

            % Circular quadrature
            q2d = Quad2D(rq);

            % Local correction
            loc = localCorrections(x,y,a,k1,rq,tol,false);
            % We are ready to compute fast convolutions.
            onlineEBD = @(v)(q2d.conv(x,y,v) + loc*v);
        end
        function[MVx,MVy,rq,locx,locy] = offline_dEBD(this,X,Y,a,tol)

            % Rescaling
            rMax = rMaxCalc(X,Y); % diameter
            x = X/rMax; % in B(0,1)
            y = Y/rMax; % in B(0,1)
            k1 = this.dilatation(rMax); % k1 must be approximated on {a < r < 1}
            rq = k1.radialQuadKernel(a,tol*rMax,'grad',true);
            % the tolerance is set to tol*rMax because G'(rMax*u) = 1/rMax
            % k1'(u).
            q2d = Quad2D(rq);
            [q2dx,q2dy] = grad(q2d);
            loc = localCorrections(x,y,a,k1,rq,tol,true);
            locx = loc{1}; locy = loc{2};
            MVx = @(v)(-(q2dx.conv(x,y,v))/rMax + locx*v/rMax);
            MVy = @(v)(-(q2dy.conv(x,y,v))/rMax + locy*v/rMax);
        end

        % Display
        function[] = disp(this)
            fprintf('Kernel : function %s \n',func2str(this.func));
        end

        function[out] = eval(this,x)
            fun = this.func;
            out = fun(x);
            out(abs(x) < 1e-12) = this.lim0;
        end
        function[out] = evalDer(this,x)
            fun = this.der;
            out = fun(x);
            out(abs(x) < 1e-12) = 0;
        end
        % Setters
        function[this] = setScalFunc(this,f)
            this.scalFunc = f;
        end
        function[this] = setNormFunc(this,g)
            this.normFunc = g;
        end
        function[this] = setStartFreq(this,k)
            this.startFreq = k;
        end
        % Addition of two kernels
        function[C] = plus(k1,k2)
            C = Kernel(@(x)(k1.func(x) + k2.func(x)),@(x)(k1.der(x) + k2.der(x)));
            sf1 = k1.scalFunc;
            sf2 = k2.scalFunc;
            C.scalFunc =  @(a,b,rho)(sf1(a,b,rho) + sf2(a,b,rho));
            C.customRadialQuad = @(a,tol,varargin)(k1.radialQuadKernel(a,tol/2,varargin{:}) ...
                + k2.radialQuadKernel(a,tol/2,varargin{:}));
            C.singular = or(k1.singular,k2.singular);
            C.lim0 = k1.lim0 + k2.lim0;
        end
        % Multiplication by a constant
        function[c] = mtimes(lambda,this)
            if and(isa(lambda,'double'),isscalar(lambda))
                assert(isa(this,'Kernel'));
                c = this;
                c.func = @(x)(lambda*this.func(x));
                c.der = @(x)(lambda*this.der(x));
                c.scalFunc = @(a,b,rho)(lambda*this.scalFunc(a,b,rho));
                c.normFunc = @(a,b)(abs(lambda)*this.normFunc(a,b));
                c.lim0 = lambda*this.lim0;
                c.customRadialQuad = @(a,tol,varargin)(lambda*this.radialQuadKernel(a,tol/lambda,varargin{:}));
            else
                assert(and(isa(this,'double'),isscalar(this)))
                c = times(this,lambda);
            end
        end
    end
    methods (Access = public)
        % Wrapper
        function[rq] = radialQuadKernel(this,a,tol,varargin)
            if ~ isempty(this.customRadialQuad)
                rq = this.customRadialQuad(a,tol,varargin{:});
            else
                rq = RadialQuadrature(a,this,tol,varargin{:});
            end
        end
        % Dilatation of the kernel : G(x) -> G(R*x).
        function[this] = dilatation(old,lambda)
            oldFunc = old.func;
            oldDer = old.der;
            ffunc = @(x)(oldFunc(lambda*x));
            dder = @(x)(lambda*oldDer(lambda*x));
            this = Kernel(ffunc,dder);
            this = this.setStartFreq(lambda*old.startFreq);
            this.lim0 = old.lim0;
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function bario = calcBario(DSO3,ori,tetra)
% compute bariocentric coordinates for an orientation
%
% Input
%  DSO3  -
%  ori   - @orientation
%  tetra - indices to tetrahegons
%
% Output
%  bario - bariocentric coordinates
%

% compute vertices
vertices = DSO3.subSet(DSO3.tetra(tetra,:));

% translate everything such that ori becomes the identity
vertices = repmat(inv(ori(:)),1,4) .* reshape(vertices,[],4);

% and project to fundamental region
vertices = reshape(project2FundamentalRegion(vertices),[],4);

% project to three dimensional space
xyz = vertices.Rodrigues;

% compute bariocentic coordinates
bario = calcZeroBario(xyz(:,1),xyz(:,2),xyz(:,3),xyz(:,4));

end

author: mtex-toolbox

repo: mtex

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


text:
function data = h5load(filename, path)
%
% data = H5LOAD(filename)
% data = H5LOAD(filename, path_in_file)
%
% Load data in a HDF5 file to a Matlab structure.
%
% Parameters
% ----------
%
% filename
%     Name of the file to load data from
% path_in_file : optional
%     Path to the part of the HDF5 file to load
%

% Author: Pauli Virtanen <pav@iki.fi>
% This script is in the Public Domain. No warranty.

if nargin > 1
    path_parts = regexp(path, '/', 'split');
else
    path = '';
    path_parts = [];
end

loc = H5F.open(filename, 'H5F_ACC_RDONLY', 'H5P_DEFAULT');
try
    data = load_one(loc, path_parts, path);
    H5F.close(loc);
catch exc
    H5F.close(loc);
    rethrow(exc);
end

% ---------------------------------------------------------------------
function data = load_one(loc, path_parts, full_path)
% Load a record recursively.

while ~isempty(path_parts) & strcmp(path_parts{1}, '')
    path_parts = path_parts(2:end);
end

data = struct();

num_objs = H5G.get_num_objs(loc);

%
% Load groups and datasets
%
for j_item = 0:num_objs-1
    objtype = H5G.get_objtype_by_idx(loc, j_item);
    objname = H5G.get_objname_by_idx(loc, j_item);

    if objtype == 0
        % Group
        name = regexprep(objname, '.*/', '');

        if isempty(path_parts) | strcmp(path_parts{1}, name)
            if ~isempty(regexp(name,'^[a-zA-Z].*'))
                try
                    group_loc = H5G.open(loc, name);
                catch
                    d=1;
                end
                try
                    sub_data = load_one(group_loc, path_parts(2:end), full_path);
                    H5G.close(group_loc);
                catch exc
                    H5G.close(group_loc);
                    rethrow(exc);
                end
                if isempty(path_parts)
                    data = setfield(data, name, sub_data);
                else
                    data = sub_data;
                    return
                end
            end
        end

    elseif objtype == 1
        % Dataset
        name = regexprep(objname, '.*/', '');

        if isempty(path_parts) | strcmp(path_parts{1}, name)
            if ~isempty(regexp(name,'^[a-zA-Z].*'))
                dataset_loc = H5D.open(loc, name);
                try
                    sub_data = H5D.read(dataset_loc, 'H5ML_DEFAULT', 'H5S_ALL','H5S_ALL','H5P_DEFAULT');
                    H5D.close(dataset_loc);
                catch exc
                    H5D.close(dataset_loc);
                    rethrow(exc);
                end

                sub_data = fix_data(sub_data);

                if isempty(path_parts)
                    data = setfield(data, name, sub_data);
                else
                    data = sub_data;
                    return
                end
            end
        end
    end
end

% Check that we managed to load something if path walking is in progress
if ~isempty(path_parts)
    error(sprintf('Path "%s" not found in the HDF5 file', full_path));
end

% ---------------------------------------------------------------------
function data = fix_data(data)
% Fix some common types of data to more friendly form.

if isstruct(data)
    fields = fieldnames(data);
    if length(fields) == 2 & strcmp(fields{1}, 'r') & strcmp(fields{2}, 'i')
        if isnumeric(data.r) & isnumeric(data.i)
            data = data.r + 1j*data.i;
        end
    end
end

if isnumeric(data) & ndims(data) > 1
    % permute dimensions
    data = permute(data, fliplr(1:ndims(data)));
end

author: BUNPC


text:
function ng = ball_grid_count ( n, r, c )

%*****************************************************************************80
%
%% BALL_GRID computes grid points inside a ball.
%
%  Discussion:
%
%    The grid is defined by specifying the radius and center of the ball,
%    and the number of subintervals N into which the horizontal radius
%    should be divided.  Thus, a value of N = 2 will result in 5 points
%    along that horizontal line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Input, real R, the radius of the ball.
%
%    Input, real C(3), the coordinates of the center of the ball.
%
%    Output, integer NG, the number of grid points inside the ball.
%
  ng = 0;

  for i = 0 : n
    x = c(1) + r * 2 * i / ( 2 * n + 1 );

    for j = 0 : n
      y = c(2) + r * 2 * j / ( 2 * n + 1 );

      for k = 0 : n

        z = c(3) + r * 2 * k / ( 2 * n + 1 );

        if ( r * r < ( x - c(1) ).^2 ...
                   + ( y - c(2) ).^2 ...
                   + ( z - c(3) ).^2 )
          break
        end
        ng = ng + 1;
        if ( 0 < i )
          ng = ng + 1;
        end
        if ( 0 < j )
          ng = ng + 1;
        end
        if ( 0 < k )
          ng = ng + 1;
        end
        if ( 0 < i && 0 < j )
          ng = ng + 1;
        end
        if ( 0 < i && 0 < k )
          ng = ng + 1;
        end
        if ( 0 < j && 0 < k )
          ng = ng + 1;
        end
        if ( 0 < i && 0 < j && 0 < k )
          ng = ng + 1;
        end
      end
    end
  end

  return
end


text:
function printClosedloopData(mpciter, u, x, t_Elapsed,output_data)
    global fst_output_data ;
    global snd_output_data ;
    if size(x,2)==2
        fprintf(' %3d  |%+11.4f %+11.4f         %+11.4f  %+11.4f %+8.3f', ...
            mpciter, u(1), u(2), x(1), x(2), t_Elapsed);
        fst_output_data = [output_data; [mpciter, u(1), u(2), x(1), x(2), t_Elapsed]];
    else
        fprintf(' %3d  |%+11.4f %+11.4f %+11.4f %+11.4f  %+11.4f  %+11.4f', ...
            mpciter, u(1), u(2), u(3), x(1), x(2), x(3));
        snd_output_data = [output_data; [mpciter, u(1), u(2), u(3), x(1), x(2), x(3)]];
    end

end

author: juchengquan


text:
function R = create_2d_image_resize_matrix(insz, outsz)
    t1 = (insz(1) - outsz(1)) * repmat(linspace(0, 1, outsz(1))', [1, outsz(2)]);
    t2 = (insz(2) - outsz(2)) * repmat(linspace(0, 1, outsz(2)), [outsz(1), 1]);
%     R = create_2d_bicubic_interp_matrix(cat(3, t1, t2), insz);
%      R = create_2d_bicubic_interp_matrix(cat(3, t1, t2), insz);
     R = create_2d_bilinear_interp_matrix(cat(3, t1, t2), insz);
end

repo: pTVreg


text:
function reply = default(query,value)
%DEFAULT gets response to IFISS prompt
%   reply = default(query,value);
%   input
%          query   character string: asks a question
%          value   integer: the default response
%
%   IFISS function: AR 31 August 2005, HCE 2 January 2010.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage (see readme.m)
global BATCH FID
if exist('BATCH') & BATCH==1,
   replycell=textscan(FID,'%f%*[^\n]',1);
   reply=deal(replycell{:});
   disp(query)
   disp(reply)
else
   if length(query)>=26 && strcmp(query(10:26),'timestep sequence'),
      reply=input([query,' : '],'s');
   else
      reply=input([query,' : ']);
   end
   if isempty(reply), reply=value; end
end
return


text:
function [net,opts]=train_net(net,opts)

    if ~isfield(opts,'datatype')
        opts.datatype='single';
    end

    opts.training=1;

    if ~isfield(opts.parameters,'learning_method')
        opts.parameters.learning_method='sgd';
    end

    if ~isfield(opts,'display_msg')
        opts.display_msg=1;
    end
    opts.TrainMiniBatchError=[];
    opts.TrainMiniBatchError_Top5=[];
    opts.TrainMiniBatchLoss=[];

    tic

    opts.order=randperm(opts.n_train);

    if opts.parameters.selective_sgd==1
        [ net,opts ] = selective_sgd( net,opts );
    end

    batch_dim=length(size(opts.train));
    idx_nd=repmat({':'},[1,batch_dim]);

    for mini_b=1:opts.n_train_batch

        idx=opts.order(1+(mini_b-1)*opts.parameters.batch_size:mini_b*opts.parameters.batch_size);
        idx_nd{end}=idx;
        opts.idx_nd=idx_nd;
        opts.idx=idx;
        res(1).x=opts.train(idx_nd{:});

        %classification
        if strcmp(net.layers{end}.type,'softmaxloss')
            res(1).class=opts.train_labels(idx);
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%forward%%%%%%%%%%%%%%%%%%%
        [ net,res,opts ] = net_ff( net,res,opts );
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%backward%%%%%%%%%%%%%%%%
        opts.dzdy=1.0;

        [ net,res,opts ] = net_bp( net,res,opts );
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %%summarize the current batch
        loss=double(gather(mean(res(end).x(:))));

        if strcmp(net.layers{end}.type,'softmaxloss')
            err=error_multiclass(res(1).class,res);
            opts.TrainMiniBatchError=[opts.TrainMiniBatchError;err(1)/opts.parameters.batch_size];
            opts.TrainMiniBatchError_Top5=[opts.TrainMiniBatchError_Top5;err(2)/opts.parameters.batch_size];
            if opts.display_msg==1
                disp(['Minibatch loss: ', num2str(loss),...
                    ', top 1 err: ', num2str(opts.TrainMiniBatchError(end)),...
                    ',top 5 err:,',num2str(opts.TrainMiniBatchError_Top5(end))])
            end
        end

        opts.TrainMiniBatchLoss=[opts.TrainMiniBatchLoss;loss];
        if (~isfield(opts.parameters,'iterations'))
            opts.parameters.iterations=0;
        end
        opts.parameters.iterations=opts.parameters.iterations+1;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%stochastic gradients descent%%%%%%%%%%%%%%%%%%%%%%%%
        [ net,res,opts ] = opts.parameters.learning_method( net,res,opts );
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    end

    %%summarize the current epoch
     if ~isfield(opts,'results')||~isfield(opts.results,'TrainEpochLoss')
        opts.results.TrainEpochLoss=[];
        opts.results.TrainEpochError=[];
        opts.results.TrainEpochError_Top5=[];
     end

    opts.results.TrainEpochLoss=[opts.results.TrainEpochLoss;mean(opts.TrainMiniBatchLoss(:))];

    if strcmp(net.layers{end}.type,'softmaxloss')
        opts.results.TrainEpochError=[opts.results.TrainEpochError;mean(opts.TrainMiniBatchError(:))];
        opts.results.TrainEpochError_Top5=[opts.results.TrainEpochError_Top5;mean(opts.TrainMiniBatchError_Top5(:))];
        disp(['Epoch ',num2str(opts.parameters.current_ep),...
         ', training loss: ', num2str(opts.results.TrainEpochLoss(end)),...
                ', top 1 err: ', num2str(opts.results.TrainEpochError(end)),...
                ',top 5 err:,',num2str(opts.results.TrainEpochError_Top5(end))])

    end

    if opts.RecordStats==1
        if ~isfield(opts,'results')||~isfield(opts.results,'TrainMiniBatchLoss')
            opts.results.TrainMiniBatchLoss=[];
            opts.results.TrainMiniBatchError=[];
            opts.results.TrainMiniBatchError_Top5=[];
        end
        opts.results.TrainMiniBatchLoss=[opts.results.TrainLoss;opts.TrainMiniBatchLoss];
        opts.results.TrainMiniBatchError=[opts.results.TrainError;opts.TrainMiniBatchError];
        opts.results.TrainMiniBatchError_Top5=[opts.results.TrainError_Top5;opts.TrainMiniBatchError_Top5];

    end

    toc;

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
clear;
clc
close all;
addpath('../../library/nav_lib');

%% inital state
SV_pos_ecef = [26580000, 0, 0]';
SV_vel_evef = [0 1755 3170]';
lat = 0.785398;
lon =  0.174533;
ght = 1000;
R0 = 6378137;
e = 0.0818191908425;
omaga = 7.29E-05;
c = 299792458;

[~, RE ] = ch_earth(lat, ght);

usr_ecef = ch_lla2ecef(lat, lon, ght);

 [az, el] = satellite_az_el(SV_pos_ecef, usr_ecef);
 rad2deg(el)
 rad2deg(az)

author: yandld


text:
function RegDemons(handles)
%
% 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/>.

global planC stateS;
indexS = planC{end};

    originF = [str2double(get(handles.originFx, 'string')) ...
               str2double(get(handles.originFy, 'string')) ...
               str2double(get(handles.originFz, 'string'))];
    spacingF = [str2double(get(handles.spacingFx, 'string')) ...
               str2double(get(handles.spacingFy, 'string')) ...
               str2double(get(handles.spacingFz, 'string'))];

    originM = [str2double(get(handles.originMx, 'string')) ...
               str2double(get(handles.originMy, 'string')) ...
               str2double(get(handles.originMz, 'string'))];
    spacingM = [str2double(get(handles.spacingMx, 'string')) ...
               str2double(get(handles.spacingMy, 'string')) ...
               str2double(get(handles.spacingMz, 'string'))];

    Hist_Level = str2double(get(handles.Demons_HistLevel, 'string'));
    MatchPoints = str2double(get(handles.Demons_MatchPoints, 'string'));
    iternum = str2double(get(handles.Demons_iternum, 'string'));
    sd = str2double(get(handles.Demons_sd, 'string'));

    output = cell(1, 8);

    aIndex = stateS.currentAxis;
    hAxis = stateS.handle.CERRAxis(aIndex);
    axisInfo = get(hAxis, 'userdata');

    FImg = axisInfo.scanObj(stateS.imageRegistrationBaseDataset).data2M;
    MImg = axisInfo.scanObj(stateS.imageRegistrationMovDataset).data2M;

%     FImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
%     MImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;
%
%     %     downsample the input datasets
%     if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))
%
%         tic;
%         FImg = imdesample3d(FImg,2,2);
%         MImg = imdesample3d(MImg,2,2);
%         toc;
%         DesampleTime = num2str(floor(toc/60));
%
%     end

    %     call registration method
    tic;
    [im, DF_X, DF_Y] = Demon(int16(FImg), originF, spacingF, ...
                             int16(MImg), originM, spacingM, ...
                             Hist_Level, MatchPoints, iternum, sd);

    toc;

    RegTime = num2str(floor(toc/60));

    imtool(im, [min(im(:)) max(im(:))]);
    imtool(DF_X.^2+DF_Y.^2); colormap colorcube;

% %update the transM;
%
%     Tv = [-Offset(4) Offset(5) 0];
%     rot = eye(3);
%     rot(1) = cos(Offset(1)*pi/180); rot(2,2) = rot(1);
%     rot(1,2) = -sin(Offset(1)*pi/180);
%     rot(2,1) = -rot(1,2);
%
%     TM = eye(4);
%     TM(:,4) = [Tv 1];
%
%     RM = eye(4);
%     RM(1:3, 1:3) = rot;
%
%     newTransform = TM*RM;
%
%     scanSetM = stateS.imageRegistrationMovDataset;
%     oldTransM = getTransM(stateS.imageRegistrationMovDatasetType, scanSetM, planC);
%     if isempty(oldTransM), oldTransM = eye(4); end;
%     planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = (newTransform * oldTransM);
% %     planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = newTransform;
%
%     %planC{indexS.dose}(scanSetM).transM = (newTransform * oldTransM);
%     sliceCallBack('refresh');
%

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR