text:
function mm_io_test05 ( )

%*****************************************************************************80
%
%% MM_IO_TEST05 tests MM_SIZE_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MM_IO_TEST05:\n' );
  fprintf ( 1, '  MM_SIZE_WRITE writes the size line \n' );
  fprintf ( 1, '  to a Matrix Market file.\n' );

  output_filename = 'mm_io_test05.mm';

  id = '%%MatrixMarket';
  type = 'matrix';
  rep = 'array';
  field = 'real';
  symm = 'general';

  nrow = 5;
  ncol = 4;
  nnz = 13;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the size information we intend to write\n' );
  fprintf ( 1, '  to a Matrix Market file:\n' );
  fprintf ( 1, '\n' );

  mm_size_print ( output_filename, rep, symm, nrow, ncol, nnz );

  output_unit = fopen ( output_filename, 'w' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MM_IO_TEST05 - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'MM_IO_TEST05 - Fatal error!' );
  end

  mm_header_write ( output_unit, id, type, rep, field, symm );

  mm_comment_write ( output_unit, ' ' );
  mm_comment_write ( output_unit, 'This is the file mm_io_test05.mm' );
  mm_comment_write ( output_unit, ' ' );
  mm_comment_write ( output_unit, 'It is a partial Matrix Market file,' );
  mm_comment_write ( output_unit, 'containing a header line, comments, and sizes.' );
  mm_comment_write ( output_unit, ' ' );

  mm_size_write ( output_unit, rep, nrow, ncol, nnz );

  fclose ( output_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'MM_IO_TEST05:\n' );
  fprintf ( 1, '  The partial MM file "%s" containing \n', output_filename );
  fprintf ( 1, '  the header, comments, and sizes has been created.\n' );
  fprintf ( 1, '\n' );

  return
end


text:
function krn = spm_smoothkern(fwhm,x,t)
% Generate a Gaussian smoothing kernel
% FORMAT krn = smoothing_kernel(fwhm,x,t)
%     fwhm - full width at half maximum
%     x    - position
%     t    - either 0 (nearest neighbour) or 1 (linear).
%            if only two arguments are passed, then a value
%            of one is assumed.
%     krn  - value of kernel at position x
%
% For smoothing images, one should really convolve a Gaussian
% with a sinc function.  For smoothing histograms, the
% kernel should be a Gaussian convolved with the histogram
% basis function used. This function returns a Gaussian
% convolved with a triangular (1st degree B-spline) basis
% function (by default). A Gaussian convolved with a hat
% function (0th degree B-spline) can also be returned.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

if nargin<3, t = 1; end;

% Variance from FWHM
s = (fwhm/sqrt(8*log(2)))^2+eps;

% The simple way to do it. Not good for small FWHM
% krn = (1/sqrt(2*pi*s))*exp(-(x.^2)/(2*s));

if t==0
    % Gaussian convolved with 0th degree B-spline
    % int(exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= -0.5..0.5)
    w1  = 1/sqrt(2*s);
    krn = 0.5*(erf(w1*(x+0.5))-erf(w1*(x-0.5)));
    krn(krn<0) = 0;

elseif t==1
    % Gaussian convolved with 1st degree B-spline
    %  int((1-t)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= 0..1)
    % +int((t+1)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t=-1..0)
    w1  =  0.5*sqrt(2/s);
    w2  = -0.5/s;
    w3  = sqrt(s/2/pi);
    krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) - 2*erf(w1*x   ).* x)...
          +w3*(exp(w2*(x+1).^2)     + exp(w2*(x-1).^2)     - 2*exp(w2*x.^2));
    krn(krn<0) = 0;

else
    error('Only defined for nearest neighbour and linear interpolation.');
    % If anyone knows a nice formula for a sinc function convolved with a
    % a Gaussian, then that could be quite useful.
end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% run_simMegaPressShapedEdit.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% This script is run simply by editing the input parameters and then
% clicking "Run".
%
% DESCRIPTION:
% This script simulates a MEGA-PRESS experiment with fully shaped editing
% pulses.  Phase cycling of editing pulses is performed, and summation
% across phase cycles is performed.
%
% INPUTS:
% To run this script, edit the parameters below as desired and then click
% "run":
% editWaveform          = name of editing pulse waveform.
% editOnFreq            = freqeucny of edit on pulse[ppm]
% editOffFreq           = frequency of edit off pulse[ppm]
% editTp                = duration of editing pulses[ms]
% Npts                  = number of spectral points
% sw                    = spectral width [Hz]
% Bfield                = magnetic field strength [Tesla]
% lw                    = linewidth of the output spectrum [Hz]
% taus                  = vector of pulse sequence timings  [ms]
% spinSys               = spin system to simulate
% editPhCyc1            = vector of phase cycling steps for 1st editing pulse [degrees]
% editPhCyc2            = vector of phase cycling steps for 2nd editing pulse [degrees]
%
% OUTPUTS:
% outON             = Simulated MEGA-PRESS edit-ON spectrum.
% outOFF            = Simulated ExTE-MEGA-PRESS edit-OFF spectrum.
% outDIFF           = Simulated ExTE-MEGA-PRESS difference spectrum.

% ************INPUT PARAMETERS**********************************
editWaveform='sampleEditPulse.pta'; %name of editing pulse waveform.
editOnFreq=1.88; %freqeucny of edit on pulse[ppm]
editOffFreq=7.4; %frequency of edit off pulse[ppm]
editTp=20; %duration of editing pulses[ms]
Npts=2048; %number of spectral points
sw=2000; %spectral width [Hz]
Bfield=3; %magnetic field strength [Tesla]
lw=2; %linewidth of the output spectrum [Hz]
taus=[5,... %Time from excitation to 1st refoc pulse [ms]
    17,...  %Time from 1st refoc pulse to 1st editing pulse [ms]
    17,...  %Time from 1st editing pulse to 2nd refoc pulse [ms]
    17,...  %Time from 2nd refoc pusle to 2nd editing pulse [ms]
    12];    %Time from 2nd editing pulse to ADC onset [ms]
spinSys='GABA'; %spin system to simulate
centreFreq=3.0; %Centre Frequency of MR spectrum [ppm];
editPhCyc1=[0 90]; %phase cycling steps for 1st editing pulse [degrees]
editPhCyc2=[0 90 180 270]; %phase cycling steps for 2nd editing pulse [degrees]
% ************END OF INPUT PARAMETERS**********************************

%Load RF waveforms
editRF=io_loadRFwaveform(editWaveform,'inv',0);

gamma=42577000; %gyromagnetic ratio

%Load spin systems
load spinSystems
sys=eval(['sys' spinSys]);

%This is the step where the editing pulse waveform (initially a pulse with
%zero-frequency) is frequency shifted to produce and edit-on and an
%edit-off pulse;
editRFon=rf_freqshift(editRF,editTp,(centreFreq-editOnFreq)*Bfield*gamma/1e6);
editRFoff=rf_freqshift(editRF,editTp,(centreFreq-editOffFreq)*Bfield*gamma/1e6);

%Initialize structures:
outON_epc=cell(length(editPhCyc1),length(editPhCyc2));
outOFF_epc=cell(length(editPhCyc1),length(editPhCyc2));
outON=struct([]);
outOFF=struct([]);

%loop through space: Don't forget to initialize the parallel processing
%toolbox workers using 'matlabpool open N' (for N workers, 12 max).

%for X=1:length(x);

for EP1=1:length(editPhCyc1)
    for EP2=1:length(editPhCyc2)
                disp(['Executing First Edit phase cycle ' num2str(EP1) ' of ' num2str(length(editPhCyc1)) ', '...
                    ' and Second Edit phase cycle ' num2str(EP2) ' of ' num2str(length(editPhCyc2)) '!!!']);
                outON_epc{EP1}{EP2}=sim_megapress_shapedEdit(Npts,sw,Bfield,lw,taus,sys,...
                    editRFon,editTp,editPhCyc1(EP1),editPhCyc2(EP2));
                outOFF_epc{EP1}{EP2}=sim_megapress_shapedEdit(Npts,sw,Bfield,lw,taus,sys,...
                    editRFoff,editTp,editPhCyc1(EP1),editPhCyc2(EP2));

        if EP1==1 && EP2==1
            outON=outON_epc{EP1}{EP2};
            outOFF=outOFF_epc{EP1}{EP2};
        else
            outON=op_addScans(outON,outON_epc{EP1}{EP2});
            outOFF=op_addScans(outOFF,outOFF_epc{EP1}{EP2});
        end

    end %end of 1st editing phase cycle loop.
end %end of 2nd editing phase cycle loop.

figure
plot(outON.ppm,outON.specs,outOFF.ppm,outOFF.specs);

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [x0,nss] = sofa_get_secondary_sources(sofa,idx,coordinate_system)
%SOFA_GET_SECONDARY_SOURCES returns x0 from the given SOFA data set
%
%   Usage: [x0,nss] = sofa_get_secondary_sources(sofa,[idx],[coordinate_system])
%
%   Input parameters:
%       sofa              - impulse response data set (SOFA struct/file)
%       idx               - index of secondary sources that should be returned
%                           (default: return all secondary sources).
%       coordinate_system - coordinate system the position and direction of the
%                           secondary sources should be specified in:
%                             'cartesian' (default)
%                             'spherical'
%
%   Output parameters:
%       x0      - secondary source matrix [n 7]
%       nss     - number of secondary sources
%
%   SOFA_GET_SECONDARY_SOURCES(sofa,idx,coordinate_system) returns secondary
%   sources as defined in the given SOFA file or struct, specified by idx.  If
%   no idx is specified all secondary sources are returned. The coordinate
%   system the position and direction of the secondary sources are specified in
%   can be given by the string 'coordinate_system', 'cartesian' is assumed as
%   default.
%
%   see also: sofa_get_header, secondary_source_positions

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 1;
nargmax = 3;
narginchk(nargmin,nargmax)
if nargin==nargmax-2
    idx = ':';
    coordinate_system = 'cartesian';
elseif nargin==nargmax-1
    if ischar(idx)
        coordinate_system = idx;
        idx = ':';
    else
        coordinate_system = 'cartesian';
    end
else
    isargvector(idx);
    isargchar(coordinate_system);
end

%% ===== Computation ====================================================
header = sofa_get_header(sofa);

if strcmp('SimpleFreeFieldHRIR',header.GLOBAL_SOFAConventions)
    % For free field HRTFs the source positions are equivalent to the apparent
    % positons of the sources
    apparent_directions = SOFAcalculateAPV(header);
    [apparent_positions(:,1) apparent_positions(:,2) apparent_positions(:,3)] = ...
        sph2cart(rad(apparent_directions(:,1)), ...
                 rad(apparent_directions(:,2)), ...
                 apparent_directions(:,3));
    x0 = [apparent_positions(idx,:) ...
          direction_vector(apparent_positions(idx,:), ...
                           zeros(size(apparent_positions,1),3)) ...
          ones(size(apparent_positions(idx,1)))];

elseif strcmp('SingleRoomDRIR',header.GLOBAL_SOFAConventions)
    to_be_implemented;

elseif strcmp('MultiSpeakerBRIR',header.GLOBAL_SOFAConventions)
    source_position = SOFAconvertCoordinates( ...
        header.SourcePosition,header.SourcePosition_Type,'cartesian');
    emitter_positions = SOFAconvertCoordinates( ...
        header.EmitterPosition,header.EmitterPosition_Type,'cartesian');
    emitter_directions = SOFAconvertCoordinates( ...
        header.EmitterView,header.EmitterView_Type,'cartesian');
    x0 = [bsxfun(@minus,emitter_positions(idx,:),source_position) ...
          emitter_directions(idx,:) ...
          ones(size(emitter_positions(idx,1)))];

else
    error('%s: %s convention currently not supported.', ...
        upper(mfilename),header.GLOBAL_SOFAConventions);
end

nss = size(x0,1);

if strcmp('cartesian',coordinate_system)
    return;
elseif strcmp('spherical',coordinate_system)
    [x0(:,1) x0(:,2) x0(:,3)] = cart2sph(x0(:,1),x0(:,2),x0(:,3));
else
    error('%s: %s is not a supported coordinate system.', ...
        upper(mfilename),coordinate_system);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [x, infos] = svrmu_nmf(V, rank, in_options)
% Stochastic variance reduced multiplicative update for non-negative matrix factorization (SVRMU-NMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
%
% Reference:
%       H. Kasai,
%       "Stochastic variance reduced multiplicative update for nonnegative matrix factorization,"
%       IEEE ICASSP, 2018.
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Mar. 22, 2017
%
% Change log:
%
%       Mar. 15, 2018 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.repeat_inneriter  = 1;
    local_options.W_sub_mode        = 'STD';
    local_options.H_sub_mode        = 'STD';
    local_options.accel             = false;
    local_options.ls                = false;
    local_options.precon            = false;
    local_options.h_repeat          = 1;
    local_options.rep_mode          = 'fix';
    local_options.stepsize_ratio    = 1; % stepsize ratio
    local_options.robust            = false;
    local_options.lambda            = 1;
    local_options.tol_optgap        = 1e-2;
    local_options.x_init_robust     = false;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    Wt = init_factors.W;
    H = init_factors.H;
    R = init_factors.R;

    % determine sub_mode
    if options.accel
        options.H_sub_mode = 'ACC';
    else
        if options.ls
            options.H_sub_mode = 'LS';
        else
            options.H_sub_mode = 'STD';
        end
        options.h_repeat = 1;
    end

    if options.precon
        options.W_sub_mode = 'Precon';
        fprintf('Unfortunately, Precon variant does not work well.\n');
    else
        options.W_sub_mode = 'STD';
    end

    if options.robust
        mode = 'R-SVRMU-NMF';
    else
        mode = 'SVRMU-NMF';
        R = zeros(m, n);
    end

    if strcmp(options.rep_mode, 'adaptive')
        rhoh = 1+(m+m*rank)/(1*(rank+1));
        alpha = 2;
        delta = 0.01;
    end

    % initialize
    method_name = sprintf('%s (%s,%s)', mode, options.W_sub_mode, options.H_sub_mode);
    epoch = 0;
    grad_calc_count = 0;
    l = zeros(m, options.batch_size) + options.lambda;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % permute samples
    if options.permute_on
        perm_idx = randperm(n);
    else
        perm_idx = 1:n;
    end
    V = V(:, perm_idx);
    H = H(:, perm_idx);

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, R, options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s (%s,%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', mode, options.W_sub_mode, options.H_sub_mode, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % store W and H, and calculate full grad
        W_0 = Wt;
        H_0 = H;
        if ~options.robust
            W0_H0_H0T =  W_0  * (H_0 * H_0')/n;
        else
            R_0 = R;
            W0_H0_H0T =  (W_0  * H_0 + R_0 ) * H_0'/n;
        end
        V_H0T = V * H_0'/n;

        if strcmp(options.W_sub_mode, 'Precon')
            invH0H0t = inv(H_0 * H_0');
            invH0H0t = max(invH0H0t,0);
        end

        grad_calc_count = grad_calc_count + m * n;

        % main inner loop
        for s = 1 : options.repeat_inneriter
            for t = 1 : options.batch_size : n - 1

                % Retrieve vt and ht
                start_idx = t;
                end_idx = t+options.batch_size-1;
                vt = V(:, start_idx:end_idx);
                ht = H(:, start_idx:end_idx);

                ht_0 = H_0(:, start_idx:end_idx);

                if ~options.robust

                    % uddate ht
                    Wtv = Wt' * vt;
                    WtW = Wt' * Wt;
                    if strcmp(options.H_sub_mode, 'ACC')
                        if strcmp(options.rep_mode, 'adaptive')
                            gamma = 1;
                            eps0 = 1;
                            j = 1;
                            rhoh_alpha = rhoh*alpha;
                            %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                            ht0 = ht;
                            while j <= rhoh_alpha && gamma >= delta*eps0
                                ht = ht .* (Wtv) ./ (WtW * ht);
                                ht = ht + (ht<eps) * eps;
                                if j == 1
                                    eps0 = norm(ht0-ht);
                                end
                                gamma = norm(ht0-ht);
                                j = j+1;
                            end
                        else
                            for iii=1:options.h_repeat
                                ht = ht .* (Wtv) ./ (WtW * ht);
                                ht = ht + (ht<eps) * eps;
                            end
                        end
                    elseif strcmp(options.H_sub_mode, 'LS')
                        ht = calc_nls_nmf(vt, Wt, 1e-16);
                        ht = ht + (ht<eps) * eps;
                    else
                        ht = ht .* (Wtv) ./ (WtW * ht);
                        ht = ht + (ht<eps) * eps;
                    end

                    % update W

                    if strcmp(options.W_sub_mode, 'Precon')
                        invhht = inv(ht * ht');
                        invhht = max(invhht,0);
                        invh0h0t = inv(ht_0 * ht_0');
                        invh0h0t = max(invh0h0t,0);
                        %Delta_minus = (W_0 * (ht_0 * ht_0') * invh0h0t + vt * ht' * invhht)/options.batch_size + V_H0T * invH0H0t;
                        %Delta_plus  = (Wt * (ht * ht') * invhht + vt * ht_0' * invh0h0t)/options.batch_size + W0_H0_H0T * invH0H0t;
                        Delta_minus = (W_0 * (ht_0 * ht_0') * invhht + vt * ht' * invhht)/options.batch_size + V_H0T * invhht;
                        Delta_plus  = (Wt * (ht * ht') * invhht + vt * ht_0' * invhht)/options.batch_size + W0_H0_H0T * invhht;
                    else
                        Delta_minus = (W_0 * (ht_0 * ht_0') + vt * ht')/options.batch_size + V_H0T;
                        Delta_plus  = (Wt * (ht * ht') + vt * ht_0')/options.batch_size + W0_H0_H0T;
                    end
                else

                    rt = R(:, start_idx:end_idx);
                    rt_0 = R_0(:, start_idx:end_idx);

                    % uddate ht
                    Wtv = Wt' * vt;
                    if strcmp(options.H_sub_mode, 'ACC')
                        if strcmp(options.rep_mode, 'adaptive')
                            gamma = 1;
                            eps0 = 1;
                            j = 1;
                            rhoh_alpha = rhoh*alpha;
                            %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                            ht0 = ht;
                            while j <= rhoh_alpha && gamma >= delta*eps0
                                Wh_r = Wt * ht + rt;
                                ht = ht .* (Wtv) ./ (Wt' * Wh_r);
                                ht = ht + (ht<eps) * eps;
                                rt = rt .* vt ./ (Wh_r + l);
                                if j == 1
                                    eps0 = norm(ht0-ht);
                                end
                                gamma = norm(ht0-ht);
                                j = j+1;
                            end
                        else
                            for iii=1:options.h_repeat
                                Wh_r = Wt * ht + rt;
                                ht = ht .* (Wtv) ./ (Wt' * Wh_r);
                                ht = ht + (ht<eps) * eps;
                                rt = rt .* vt ./ (Wh_r + l);
                            end
                        end
                    else
                        Wh_r = Wt * ht + rt;
                        ht = ht .* (Wtv) ./ (Wt' * Wh_r);
                        ht = ht + (ht<eps) * eps;
                        rt = rt .* vt ./ (Wh_r + l);
                    end

                    % update W
                    Delta_minus = ((W_0 * ht_0 + rt_0) * ht_0' + vt * ht')/options.batch_size + V_H0T;
                    Delta_plus  = ((Wt * ht + rt) * ht' + vt * ht_0')/options.batch_size + W0_H0_H0T;
                end

                if options.stepsize_ratio == 1
                    Wt = Wt .* (Delta_minus ./ Delta_plus);
                else
                    Wt = (1-options.stepsize_ratio) * Wt + options.stepsize_ratio * Wt .* (Delta_minus ./ Delta_plus);
                end

                Wt = Wt + (Wt<eps) * eps;

                % store new h
                H(:, start_idx:end_idx) = ht;

                % Update R
                if options.robust
                    R(:, start_idx:end_idx) = rt;
                end

                grad_calc_count = grad_calc_count + m * options.batch_size;
            end
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, R, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = Wt;
    x.H(:,perm_idx) = H;
    x.R = R;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function y = isreal( x )
y = isreal( cvxaff( x ) );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
classdef prtFeatSelLlnn < prtFeatSel
    % prtFeatSelLLnn  Local Learning based feature selection
    %
    %   FEATSEL = prtFeatSelLlnn returns a local learning based feature selection
    %   object.
    %
    %    FEATSEL = prtFeatSelLlnn(PROPERTY1, VALUE1, ...) constructs a
    %    prttFeatSelExhaustive object FEATSEL with properties as specified by
    %    PROPERTY/VALUE pair
    %
    %    A prtFeatSelExhaustive object has the following properties:
    %
    %    selectedFeatures       - The indices of the features selected,
    %                             a read-only parameter, found by training.
    %    verbosePlot            - Toggles plotting on/off during training
    %    nMaxIterations         - The maximum number of iterations
    %    normalizedWeightCutOff - The weight threshold for keeping a feature
    %
    %    The following features are settable, and are related to the
    %    training algorithm. Please see reference for further information.
    %
    %    kernelSigma
    %    sparsnessLambda
    %    vGradNMaxSteps
    %    vGradInitStepSize
    %    vGradChangeThreshold
    %    vGradNMaxStepSizeChanges
    %    weightChangeThreshold
    %
    %   Reference:
    %   http://www.computer.org/portal/web/csdl/doi/10.1109/TPAMI.2009.190
    %
    %   A prtFeatSelExhaustive object inherits the TRAIN and RUN methods
    %   from prtClass.
    %
    %   Example:
    %
    %   dataSet = prtDataGenSpiral;   % Create a 2 dimensional data set
    %   nNoiseFeatures = 100;      % Append 100 irrelevant features
    %   dataSet = prtDataSetClass(cat(2,dataSet.getObservations,randn([dataSet.nObservations, nNoiseFeatures])), dataSet.getTargets);
    %   featSel = prtFeatSelLlnn('verbosePlot',true);  % Create the feature
    %                                                  % selection object.
    %   featSel.nMaxIterations = 10;                   % Set the max # of
    %                                                  % iterations.
    %   featSel = featSel.train(dataSet);              % Train
    %
    %   See Also:  prtFeatSelStatic, prtFeatSelSfs, prtFeatSelExhaustive

    properties (SetAccess=private)
        % Required by prtAction
        name = 'Local Learning Nearest Neighbor'
        nameAbbreviation = 'LLNN'
    end

    properties

        kernelSigma = 1;
        sparsnessLambda = 2;
        vGradNMaxSteps = 100;
        vGradInitStepSize = 1;
        vGradChangeThreshold = 1e-4;
        vGradNMaxStepSizeChanges = 12;
        weightChangeThreshold = 0.01;
        nMaxIterations = 25;         % The maximum number of iterations
        normalizedWeightCutOff = 0.05;  % The weight threshold to include a feature.

        % Learned
        weights = [] ;           % The weights of each feature
        selectedFeatures = [];   % The selected features

        verbosePlot = false;     % Toggles plotting on/off during training
    end

    methods

        % Constructor %%
        function Obj = prtFeatSelLlnn(varargin)

            % Allow for string, value pairs
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
    end

    methods (Access=protected,Hidden=true)

        % Train %%
        function Obj = trainAction(Obj,DS)

            % Remove stuff from data set, too many calls to getObservations
            X = DS.getObservations();
            Y = DS.getTargetsClassInd();

            vExpDistanceFunction = @(v,zBar)exp(-zBar*(v.^2));

            vFitnessFunction = @(v,cExpDistance)sum(log(1+cExpDistance)) + Obj.sparsnessLambda*sum(v.^2);

            w = ones(DS.nFeatures,1);
            wOld = w;
            wChanges = nan(Obj.nMaxIterations,1);
            for iter = 1:Obj.nMaxIterations

                % Find zBar
                % Eq: 3-5
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                cSigma = Obj.kernelSigma*sqrt(sum((w./max(w))>Obj.normalizedWeightCutOff)); % Modify kernel to the "dimensionality" imposed by w

                zBar = zeros(DS.nObservations, DS.nFeatures);
                for iSamp = 1:DS.nObservations

                    % L1 Distance from this point in each dim
                    cD = abs(bsxfun(@minus,X(iSamp,:),X));
                    cD(iSamp,:) = inf;

                    % Label of this point
                    cY = Y(iSamp);

                    % Kernel mapped distance to this point
                    cK = exp(-sum(bsxfun(@times,cD,w(:)'),2)/cSigma); % Hard coded L1 Distance

                    % Calculate the prob of nearest miss and nearst hit
                    % Eqs: 4, 5
                    cKMiss = cK;
                    cKMiss(Y==cY) = 0; % Cant be a miss if you are the same type
                    cKMiss(iSamp) = 0; % Cant be your own miss
                    cKMiss = cKMiss./sum(cKMiss);

                    cKHit = cK;
                    cKHit(Y~=cY) = 0; % Cant be a hit if you are a different type
                    cKHit(iSamp) = 0; % Cant be your own hit
                    cKHit = cKHit./sum(cKHit);

                    % Eq: 3
                    cD(iSamp,:) = 0; % inf*0 = nan
                    zBar(iSamp,:) = sum(bsxfun(@times,cD, cKMiss),1) - sum(bsxfun(@times,cD, cKHit),1);

                end
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                % Optimize v
                % Eqs: 9, 10
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                v = sqrt(w);

                cExpDistance = vExpDistanceFunction(v,zBar);
                vScore = vFitnessFunction(v,cExpDistance); % Eq: 9

                vScoreOld = vScore;
                for iStep = 1:Obj.vGradNMaxSteps
                    cStepSize = Obj.vGradInitStepSize*2; % we multiply by 2 since we divide by 2 at every iteration (even the first)

                    % Transform cExpDistance into necessary term for the
                    % gradient, we use the same name... just because
                    cExpDistance = cExpDistance./(1 + cExpDistance);
                    cExpDistance(isnan(cExpDistance)) = 0; % inf/(inf + 1)

                    % Eq: 10
                    cChange = Obj.sparsnessLambda*ones(size(w)) - sum(bsxfun(@times,cExpDistance, zBar),1)';

                    for iStepSize = 1:Obj.vGradNMaxStepSizeChanges
                        cStepSize = cStepSize / 2; % Lower the step size, we get here if 1) first iteration (we multiplied the step size by 2) or 2) step size was too big so we need to shrink it.

                        vNew = v - cStepSize*cChange.*v; % Eq: 10

                        cExpDistance = vExpDistanceFunction(vNew,zBar);

                        vNewScore = vFitnessFunction(vNew,cExpDistance);

                        if vNewScore < vScore
                            vScore = vNewScore;
                            v = vNew;
                            % Note that cExpDistance will be reused but
                            % immediately over written with its alter ego
                            % this saves a little mem I guess
                            break
                        else
                            % We increased the fitness (bad) so the step
                            % size must be too big.
                            % Continue in the loop decrease the step size
                            % by 50%.
                        end
                    end

                    if abs(vScore-vScoreOld)/mean([vScore vScoreOld]) < Obj.vGradChangeThreshold
                        break
                    else
                        vScoreOld = vScore;
                    end
                end
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                w = v.^2;

                cWChange = norm(abs(w-wOld));

                %%% Other ideas for exit criterion %%%
                % cWChange = norm(abs(w-wOld))/sqrt(length(w));
                % cWChange = mean(abs(w-wOld)./mean(cat(2,w,wOld),2)); % Average percent change
                % wNorm = w ./ max(w);
                % wOldNorm = wOld ./ max(wOld);
                % cWChange = mean(abs(wNorm-wOldNorm)./mean(cat(2,wNorm,wOldNorm),2)); % Normalized Average percent change
                % cWChange = mean(abs(wNorm-wOldNorm)); % Normalized Average change

                if Obj.verbosePlot
                    subplot(2,1,1)
                    stem(w./max(w));
                    title('Feature Importance','FontSize',14)
                    xlabel('Feature')
                    ylabel('Normalized Feature Weight');
                    xlim([1 length(w)]);

                    subplot(2,1,2)
                    wChanges(iter) = cWChange;
                    plot(1:iter,wChanges(1:iter),'b-',1:iter,wChanges(1:iter),'rx')
                    title('Optimization Exit Criterion','FontSize',14)
                    xlabel('Iteration')
                    ylabel('Weight Change From Last Iteration')
                    xlim([0.5 iter+0.5]);
                    set(gca,'XTick',1:iter);

                    drawnow;
                end

                if cWChange < Obj.weightChangeThreshold
                    break
                end
                wOld =w;

            end

            Obj.weights = w;
            Obj.selectedFeatures = find(w > Obj.normalizedWeightCutOff);
        end

        % Run %
        function DataSet = runAction(Obj,DataSet)
            DataSet = DataSet.retainFeatures(Obj.selectedFeatures);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%% (Internal) Generate a PIDs work list
%
%
%   [pid_starts, pid_ends] = TaskPartition( task_size, cant_pid)
%
% Arguments:
%
%      + task_size: A positive integer with the size of the job
%
%      + cant_pid: The amount of PIDs to work
%
% Output:
%
%      + pid_starts: An index array of size cant_pid x 1 with the starting
%      indexes for each PID
%
%      + pid_ends: An index array of size cant_pid x 1 with the ending
%      indexes for each PID
%
% Example:
%
%       [pid_starts, pid_ends] = TaskPartition( 10, 2)
%
%       pid_starts = [ 1 6 ]
%       pid_ends =   [ 5 10 ]
%
% See also ECGwrapper
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function [pid_starts, pid_ends] = TaskPartition( task_size, cant_pid)

%no es recomendable hacerlo mas grande de cant_recs la particion.
cant_pid = min(cant_pid, task_size);

things2do = fix(task_size / cant_pid);

remainder = rem(task_size, cant_pid);

cantThingsXpid = repmat(things2do, cant_pid,1);

cantThingsXpid(1:remainder) = cantThingsXpid(1:remainder) + 1;

pid_ends = cumsum(cantThingsXpid);
pid_starts = [1;pid_ends(1:end-1)+1];

author: marianux

repo: ecg-kit

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


text:
function handles = enableDataClustered(hObject, handles)

set(handles.btnSaveData, 'Enable', 'on');
set(handles.btnClusterSave, 'Enable', 'on');

handles = togglePlotSilhouette(hObject, handles, 'on');
handles = togglePlotClusteredData(hObject, handles, 'on');

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function ffmsh_2d_data_print ( title, v_num, e_num, t_num, v_xy, v_l, e_v, ...
  e_l, t_v, t_l )

%*****************************************************************************80
%
%% FFMSH_2D_DATA_PRINT prints FFMSH data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 December 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string TITLE, a title.
%
%    Input, integer V_NUM, the number of vertices.
%
%    Input, integer E_NUM, the number of boundary edges.
%
%    Input, integer T_NUM, the number of triangles.
%
%    Input, real V_XY(2,V_NUM), vertex coordinates.
%
%    Input, integer V_L(V_NUM), vertex labels.
%
%    Input, integer E_V(2,E_NUM), edge vertices.
%
%    Input, integer E_L(E_NUM), vertex labels.
%
%    Input, integer T_V(3,T_NUM), triangle vertices.
%
%    Input, integer T_L(T_NUM), triangle labels.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  i4vec_print (              v_num, v_l,  '  Vertex labels:' );
  r8mat_transpose_print ( 2, v_num, v_xy, '  Vertex coordinates:' );
  i4vec_print (              e_num, e_l,  '  Edge labels:' );
  i4mat_transpose_print ( 2, e_num, e_v,  '  Edge vertices:' );
  i4vec_print (              t_num, t_l,  '  Triangle labels:' );
  i4mat_transpose_print ( 3, t_num, t_v,  '  Triangle vertices:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function res = vl_simplenn(net, x, dzdy, res, varargin)
%VL_SIMPLENN  Evaluate a SimpleNN network.
%   RES = VL_SIMPLENN(NET, X) evaluates the convnet NET on data X.
%   RES = VL_SIMPLENN(NET, X, DZDY) evaluates the convnent NET and its
%   derivative on data X and output derivative DZDY.
%
%   This function process networks using the SimpleNN wrapper
%   format. Such networks are 'simple' in the sense that they consist
%   of a linear sequence of computational layers. You can use the
%   `dagnn.DagNN` wrapper for more complex topologies, or write your
%   own wrapper around MatConvNet computational blocks for even
%   greater flexibility.
%
%   The format of the network structure NET and of the result
%   structure RES are described in some detail below. Most networks
%   expect the input data X to be standardized, for example by
%   rescaling the input image(s) and subtracting a mean. Doing so is
%   left to the user, but information on how to do this is usually
%   contained in the `net.meta` field of the NET structure (see
%   below).
%
%   The NET structure needs to be updated as new features are
%   introduced in MatConvNet; use the `VL_SIMPLENN_TIDY()` function
%   to make an old network current, as well as to cleanup and check
%   the structure of an existing network.
%
%   Networks can run either on the CPU or GPU. Use VL_SIMPLENN_MOVE()
%   to move the network parameters between these devices.
%
%   To print or obtain summary of the network structure, use the
%   VL_SIMPLENN_DISPLAY() function.
%
%   VL_SIMPLENN_DISPLAY(..., 'OPT', VAL, ...) takes the following
%   options:
%
%   `Mode`:: `normal`
%      Specifies the mode of operation. It can be either `'normal'` or
%      `'test'`. In test mode, dropout and batch-normalization are
%      bypassed. Note that, when a network is deployed, it may be
%      preferable to *remove* such blocks altogether.
%
%   `ConserveMemory`:: `true`
%      Aggressively delete intermediate results. This in practice has
%      a very small performance hit and allows training much larger
%      models. However, it can be useful to disable it for
%      debugging.
%
%   `CuDNN`:: `true`
%      Use CuDNN when available.
%
%   `Accumulate`:: `false`
%      Accumulate gradients in back-propagation instead of rewriting
%      them. This is useful to break the computation in sub-batches.
%
%   ## The result format
%
%   SimpleNN returns the result of its calculations in the RES
%   structure array. RES(1) contains the input to the network, while
%   RES(2), RES(3), ... contain the output of each layer, from first
%   to last. Each entry has the following fields:
%
%   - `res(i+1).x`: the output of layer `i`. Hence `res(1).x` is the
%     network input.
%
%   - `res(i+1).aux`: any auxiliary output data of layer i. For example,
%     dropout uses this field to store the dropout mask.
%
%   - `res(i+1).dzdx`: the derivative of the network output relative
%     to the output of layer `i`. In particular `res(1).dzdx` is the
%     derivative of the network output with respect to the network
%     input.
%
%   - `res(i+1).dzdw`: a cell array containing the derivatives of the
%     network output relative to the parameters of layer `i`. It can
%     be a cell array for multiple parameters.
%
%   ## The network format
%
%   The network is represented by the NET structure, which contains
%   two fields:
%
%   - `net.layers` is a cell array with the CNN layers.
%
%   - `net.meta` is a grab-bag of auxiliary application-dependent
%     information, including for example details on how to normalize
%     input data, the class names for a classifiers, or details of
%     the learning algorithm. The content of this field is ignored by
%     VL_SIMPLENN().
%
%   SimpleNN is aware of the following layers:
%
%   Convolution layer::
%     The convolution layer wraps VL_NNCONV(). It has fields:
%
%     - `layer.type` contains the string `'conv'`.
%     - `layer.weights` is a cell array with filters and biases.
%     - `layer.stride` is the sampling stride (e.g. 1).
%     - `layer.pad` is the padding (e.g. 0).
%
%   Convolution transpose layer::
%     The convolution transpose layer wraps VL_NNCONVT(). It has fields:
%
%     - `layer.type` contains the string `'convt'`.
%     - `layer.weights` is a cell array with filters and biases.
%     - `layer.upsample` is the upsampling factor (e.g. 1).
%     - `layer.crop` is the amount of output cropping (e.g. 0).
%
%   Max pooling layer::
%     The max pooling layer wraps VL_NNPOOL(). It has fields:
%
%     - `layer.type` contains the string `'pool'`.
%     - `layer.method` is the pooling method (either 'max' or 'avg').
%     - `layer.pool` is the pooling size (e.g. 3).
%     - `layer.stride` is the sampling stride (usually 1).
%     - `layer.pad` is the padding (usually 0).
%
%   Normalization (LRN) layer::
%     The normalization layer wraps VL_NNNORMALIZE(). It has fields:
%
%     - `layer.type` contains the string `'normalize'` or `'lrn'`.
%     - `layer.param` contains the normalization parameters (see VL_NNNORMALIZE()).
%
%   Spatial normalization layer::
%     The spatial normalization layer wraps VL_NNSPNORM(). It has fields:
%
%     - `layer.type` contains the string `'spnorm'`.
%     - `layer.param` contains the normalization parameters (see VL_NNSPNORM()).
%
%   Batch normalization layer::
%     This layer wraps VL_NNBNORM(). It has fields:
%
%     - `layer.type` contains the string `'bnorm'`.
%     - `layer.weights` contains is a cell-array with, multiplier and
%       biases, and moments parameters
%
%     Note that moments are used only in `'test'` mode to bypass batch
%     normalization.
%
%   ReLU and Sigmoid layers::
%     The ReLU layer wraps VL_NNRELU(). It has fields:
%
%     - `layer.type` contains the string `'relu'`.
%     - `layer.leak` is the leak factor (e.g. 0).
%
%     The sigmoid layer is the same, but for the sigmoid function,
%     with `relu` replaced by `sigmoid` and no leak factor.
%
%   Dropout layer::
%     The dropout layer wraps VL_NNDROPOUT(). It has fields:
%
%     - `layer.type` contains the string `'dropout'`.
%     - `layer.rate` is the dropout rate (e.g. 0.5).
%
%     Note that the block is bypassed in `test` mode.
%
%   Softmax layer::
%     The softmax layer wraps VL_NNSOFTMAX(). It has fields
%
%     - `layer.type` contains the string`'softmax'`.
%
%   Log-loss layer and softmax-log-loss::
%     The log-loss layer wraps VL_NNLOSS(). It has fields:
%
%     - `layer.type` contains `'loss'`.
%     - `layer.class` contains the ground-truth class labels.
%
%     The softmax-log-loss layer wraps VL_NNSOFTMAXLOSS() instead. it
%     has the same parameters, but `type` contains the `'softmaxloss'`
%     string.
%
%   P-dist layer::
%     The p-dist layer wraps VL_NNPDIST(). It has fields:
%
%     - `layer.type` contains the string  `'pdist'`.
%     - `layer.p` is the P parameter of the P-distance (e.g. 2).
%     - `layer.noRoot` it tells whether to raise the distance to
%     the P-th power (e.g. `false`).
%     - `layer.epsilon` is the regularization parameter for the derivatives.
%
%   Custom layer::
%     This can be used to specify custom layers.
%
%     - `layer.type` contains the string `'custom'`.
%     - `layer.forward` is  a function handle computing the block.
%     - `layer.backward` is a function handle computing the block derivative.
%
%     The first function is called as
%
%          res(i+1) = layer.forward(layer, res(i), res(i+1))
%
%     where RES is the structure array specified before. The second function is
%     called as
%
%          res(i) = layer.backward(layer, res(i), res(i+1))
%
%     Note that the `layer` structure can contain additional custom
%     fields if needed.
%
%   See also: dagnn.DagNN, VL_SIMPLENN_TIDY(),
%   VL_SIMPLENN_DISPLAY(), VL_SIMPLENN_MOVE().

% Copyright (C) 2014-15 Andrea Vedaldi.
% 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).

opts.res = [] ;
opts.conserveMemory = false ;
opts.sync = false ;
opts.mode = 'normal' ;
opts.accumulate = false ;
opts.cudnn = true ;
opts.backPropDepth = +inf ;
opts = vl_argparse(opts, varargin);

n = numel(net.layers) ;

if (nargin <= 2) || isempty(dzdy)
  doder = false ;
else
  doder = true ;
end

if opts.cudnn
  cudnn = {'CuDNN'} ;
else
  cudnn = {'NoCuDNN'} ;
end

switch lower(opts.mode)
  case 'normal'
    testMode = false ;
  case 'test'
    testMode = true ;
  otherwise
    error('Unknown mode ''%s''.', opts. mode) ;
end

gpuMode = isa(x, 'gpuArray') ;

if nargin <= 3 || isempty(res)
  res = struct(...
    'x', cell(1,n+1), ...
    'dzdx', cell(1,n+1), ...
    'dzdw', cell(1,n+1), ...
    'aux', cell(1,n+1), ...
    'time', num2cell(zeros(1,n+1)), ...
    'backwardTime', num2cell(zeros(1,n+1))) ;
end
res(1).x = x ;

% -------------------------------------------------------------------------
%                                                              Forward pass
% -------------------------------------------------------------------------

for i=1:n
  l = net.layers{i} ;
  res(i).time = tic ;
  switch l.type
    case 'conv'
      res(i+1).x = vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, ...
        'pad', l.pad, ...
        'stride', l.stride, ...
        l.opts{:}, ...
        cudnn{:}) ;

    case 'convt'
      res(i+1).x = vl_nnconvt(res(i).x, l.weights{1}, l.weights{2}, ...
        'crop', l.crop, ...
        'upsample', l.upsample, ...
        'numGroups', l.numGroups, ...
        l.opts{:}, ...
        cudnn{:}) ;

    case 'pool'
      res(i+1).x = vl_nnpool(res(i).x, l.pool, ...
        'pad', l.pad, 'stride', l.stride, ...
        'method', l.method, ...
        l.opts{:}, ...
        cudnn{:}) ;

    case {'normalize', 'lrn'}
      res(i+1).x = vl_nnnormalize(res(i).x, l.param) ;

    case 'softmax'
      res(i+1).x = vl_nnsoftmax(res(i).x) ;

    case 'loss'
      res(i+1).x = vl_nnloss(res(i).x, l.class) ;

    case 'softmaxloss'
      res(i+1).x = vl_nnsoftmaxloss(res(i).x, l.class) ;

    case 'relu'
      if l.leak > 0, leak = {'leak', l.leak} ; else leak = {} ; end
      res(i+1).x = vl_nnrelu(res(i).x,[],leak{:}) ;

    case 'sigmoid'
      res(i+1).x = vl_nnsigmoid(res(i).x) ;

    case 'noffset'
      res(i+1).x = vl_nnnoffset(res(i).x, l.param) ;

    case 'spnorm'
      res(i+1).x = vl_nnspnorm(res(i).x, l.param) ;

    case 'dropout'
      if testMode
        res(i+1).x = res(i).x ;
      else
        [res(i+1).x, res(i+1).aux] = vl_nndropout(res(i).x, 'rate', l.rate) ;
      end

    case 'bnorm'
      if testMode
        res(i+1).x = vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, 'moments', l.weights{3}) ;
      else
        res(i+1).x = vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}) ;
      end

    case 'pdist'
      res(i+1) = vl_nnpdist(res(i).x, l.p, 'noRoot', l.noRoot, 'epsilon', l.epsilon) ;

    case 'custom'
      res(i+1) = l.forward(l, res(i), res(i+1)) ;

    otherwise
      error('Unknown layer type ''%s''.', l.type) ;
  end

  % optionally forget intermediate results
  forget = opts.conserveMemory ;
  forget = forget & (~doder || strcmp(l.type, 'relu')) ;
  forget = forget & ~(strcmp(l.type, 'loss') || strcmp(l.type, 'softmaxloss')) ;
  forget = forget & (~isfield(l, 'rememberOutput') || ~l.rememberOutput) ;
  if forget
    res(i).x = [] ;
  end

  if gpuMode & opts.sync
    wait(gpuDevice) ;
  end
  res(i).time = toc(res(i).time) ;
end

% -------------------------------------------------------------------------
%                                                             Backward pass
% -------------------------------------------------------------------------

if doder
  res(n+1).dzdx = dzdy ;
  for i=n:-1:max(1, n-opts.backPropDepth+1)
    l = net.layers{i} ;
    res(i).backwardTime = tic ;
    switch l.type

      case 'conv'
        [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
          vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, res(i+1).dzdx, ...
          'pad', l.pad, ...
          'stride', l.stride, ...
          l.opts{:}, ...
          cudnn{:}) ;

      case 'convt'
        [res(i).dzdx, dzdw{1}, dzdw{2}] = ...
          vl_nnconv(res(i).x, l.weights{1}, l.weights{2}, res(i+1).dzdx, ...
          'crop', l.crop, ...
          'upsample', l.upsample, ...
          'numGroups', l.numGroups, ...
          l.opts{:}, ...
          cudnn{:}) ;

      case 'pool'
        res(i).dzdx = vl_nnpool(res(i).x, l.pool, res(i+1).dzdx, ...
                                'pad', l.pad, 'stride', l.stride, ...
                                'method', l.method, ...
                                l.opts{:}, ...
                                cudnn{:}) ;

      case {'normalize', 'lrn'}
        res(i).dzdx = vl_nnnormalize(res(i).x, l.param, res(i+1).dzdx) ;

      case 'softmax'
        res(i).dzdx = vl_nnsoftmax(res(i).x, res(i+1).dzdx) ;

      case 'loss'
        res(i).dzdx = vl_nnloss(res(i).x, l.class, res(i+1).dzdx) ;

      case 'softmaxloss'
        res(i).dzdx = vl_nnsoftmaxloss(res(i).x, l.class, res(i+1).dzdx) ;

      case 'relu'
        if l.leak > 0, leak = {'leak', l.leak} ; else leak = {} ; end
        if ~isempty(res(i).x)
          res(i).dzdx = vl_nnrelu(res(i).x, res(i+1).dzdx, leak{:}) ;
        else
          % if res(i).x is empty, it has been optimized away, so we use this
          % hack (which works only for ReLU):
          res(i).dzdx = vl_nnrelu(res(i+1).x, res(i+1).dzdx, leak{:}) ;
        end

      case 'sigmoid'
        res(i).dzdx = vl_nnsigmoid(res(i).x, res(i+1).dzdx) ;

      case 'noffset'
        res(i).dzdx = vl_nnnoffset(res(i).x, l.param, res(i+1).dzdx) ;

      case 'spnorm'
        res(i).dzdx = vl_nnspnorm(res(i).x, l.param, res(i+1).dzdx) ;

      case 'dropout'
        if testMode
          res(i).dzdx = res(i+1).dzdx ;
        else
          res(i).dzdx = vl_nndropout(res(i).x, res(i+1).dzdx, ...
                                     'mask', res(i+1).aux) ;
        end

      case 'bnorm'
        [res(i).dzdx, dzdw{1}, dzdw{2}, dzdw{3}] = ...
          vl_nnbnorm(res(i).x, l.weights{1}, l.weights{2}, res(i+1).dzdx) ;
        % multiply the moments update by the number of images in the batch
        % this is required to make the update additive for subbatches
        % and will eventually be normalized away
        dzdw{3} = dzdw{3} * size(res(i).x,4) ;

      case 'pdist'
        res(i).dzdx = vl_nnpdist(res(i).x, l.p, res(i+1).dzdx, ...
                                 'noRoot', l.noRoot, 'epsilon', l.epsilon) ;
      case 'custom'
        res(i) = l.backward(l, res(i), res(i+1)) ;

    end % layers

    switch l.type
      case {'conv', 'convt', 'bnorm'}
        if ~opts.accumulate
          res(i).dzdw = dzdw ;
        else
          for j=1:numel(dzdw)
            res(i).dzdw{j} = res(i).dzdw{j} + dzdw{j} ;
          end
        end
        dzdw = [] ;
    end
    if opts.conserveMemory
      res(i+1).dzdx = [] ;
    end
    if gpuMode & opts.sync
      wait(gpuDevice) ;
    end
    res(i).backwardTime = toc(res(i).backwardTime) ;
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ar_omp_phase_transition_gaussian_dict_gaussian_data.mat';
N = 1024;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.NumTrials = 200;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K).gaussian;
%recovery_solver = @(Phi, K, y) atom_ranked_omp(Phi, y, K).z;
pta.run(dict_model, data_model, @ar_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
function varargout = process_plv1( varargin )
% PROCESS_PLV1: Compute the Phase locking value between one signal and all the others, in one file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2020

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Phase locking value 1xN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 670;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 0);
    % === FREQ BANDS
    sProcess.options.freqbands.Comment = 'Frequency bands for the Hilbert transform:';
    sProcess.options.freqbands.Type    = 'groupbands';
    sProcess.options.freqbands.Value   = bst_get('DefaultFreqBands');
    % === PLV METHOD
    sProcess.options.plvmethod.Comment = {'PLV: Phase locking value', 'ciPLV:  Corrected imaginary phase locking value', 'wPLI: Weighted phase lag index'; 'plv', 'ciplv', 'wpli'};
    sProcess.options.plvmethod.Type    = 'radio_label';
    sProcess.options.plvmethod.Value   = 'plv';
    % === KEEP TIME
    sProcess.options.keeptime.Comment = 'Keep time information, and estimate the PLV across trials
(requires the average of many trials)';
    sProcess.options.keeptime.Type    = 'checkbox';
    sProcess.options.keeptime.Value   = 0;
    % === PLV MEASURE
    sProcess.options.plvmeasure.Comment = {'None (complex)', 'Magnitude', 'Measure:'};
    sProcess.options.plvmeasure.Type    = 'radio_line';
    sProcess.options.plvmeasure.Value   = 2;
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    if ~isempty(sProcess.options.plvmethod.Value)
        iMethod = find(strcmpi(sProcess.options.plvmethod.Comment(2,:), sProcess.options.plvmethod.Value));
        if ~isempty(iMethod)
            Comment = str_striptag(sProcess.options.plvmethod.Comment{1,iMethod});
        else
            Comment = sProcess.Comment;
        end
    else
        Comment = sProcess.Comment;
    end
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end
    % Keep time or not: different methods
    OPTIONS.Method = sProcess.options.plvmethod.Value;
    if sProcess.options.keeptime.Value
        OPTIONS.Method = [OPTIONS.Method 't'];
    end
    % Hilbert and frequency bands options
    OPTIONS.Freqs = sProcess.options.freqbands.Value;
    OPTIONS.isMirror = 0;
    % PLV measure
    if isfield(sProcess.options, 'plvmeasure') && isfield(sProcess.options.plvmeasure, 'Value') && ~isempty(sProcess.options.plvmeasure.Value)
        switch (sProcess.options.plvmeasure.Value)
            case 1,  OPTIONS.PlvMeasure = 'none';
            case 2,  OPTIONS.PlvMeasure = 'magnitude';
        end
    else
        OPTIONS.PlvMeasure = 'magnitude';
    end

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS);
end

author: brainstorm-tools


text:
function x = solve_tril(T,b)
% SOLVE_TRIL      Left division by lower triangular matrix.
% SOLVE_TRIL(T,b) is the same as T\b but requires T to be lower triangular
% and runs faster.

error('You must run install_lightspeed in order to use this function.')

author: andrewssobral

repo: lrslibrary


text:
function [PSNR, y_est] = BM3D(y, z, sigma, profile, print_to_screen)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  BM3D is an algorithm for attenuation of additive white Gaussian noise from
%  grayscale images. This algorithm reproduces the results from the article:
%
%  [1] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image Denoising
%      by Sparse 3D Transform-Domain Collaborative Filtering,"
%      IEEE Transactions on Image Processing, vol. 16, no. 8, August, 2007.
%      preprint at http://www.cs.tut.fi/~foi/GCF-BM3D.
%
%
%  FUNCTION INTERFACE:
%
%  [PSNR, y_est] = BM3D(y, z, sigma, profile, print_to_screen)
%
%  ! The function can work without any of the input arguments,
%   in which case, the internal default ones are used !
%
%  BASIC USAGE EXAMPLES:
%
%     Case 1) Using the default parameters (i.e., image name, sigma, etc.)
%
%      [PSNR, y_est] = BM3D;
%
%     Case 2) Using an external noisy image:
%
%      % Read a grayscale image and scale its intensities in range [0,1]
%      y = im2double(imread('Cameraman256.png'));
%      % Generate the same seed used in the experimental results of [1]
%      randn('seed', 0);
%      % Standard deviation of the noise --- corresponding to intensity
%      %  range [0,255], despite that the input was scaled in [0,1]
%      sigma = 25;
%      % Add the AWGN with zero mean and standard deviation 'sigma'
%      z = y + (sigma/255)*randn(size(y));
%      % Denoise 'z'. The denoised image is 'y_est', and 'NA = 1' because
%      %  the true image was not provided
%      [NA, y_est] = BM3D(1, z, sigma);
%      % Compute the putput PSNR
%      PSNR = 10*log10(1/mean((y(:)-y_est(:)).^2))
%      % show the noisy image 'z' and the denoised 'y_est'
%      figure; imshow(z);
%      figure; imshow(y_est);
%
%     Case 3) If the original image y is provided as the first input
%      argument, then some additional information is printed (PSNRs,
%      figures, etc.). That is, "[NA, y_est] = BM3D(1, z, sigma);" in the
%      above code should be replaced with:
%
%      [PSNR, y_est] = BM3D(y, z, sigma);
%
%
%  INPUT ARGUMENTS (OPTIONAL):
%
%     1) y (matrix M x N): Noise-free image (needed for computing PSNR),
%                           replace with the scalar 1 if not available.
%     2) z (matrix M x N): Noisy image (intensities in range [0,1] or [0,255])
%     3) sigma (double)  : Std. dev. of the noise (corresponding to intensities
%                          in range [0,255] even if the range of z is [0,1])
%     4) profile (char)  : 'np' --> Normal Profile
%                          'lc' --> Fast Profile
%     5) print_to_screen : 0 --> do not print output information (and do
%                                not plot figures)
%                          1 --> print information and plot figures
%
%  OUTPUTS:
%     1) PSNR (double)          : Output PSNR (dB), only if the original
%                                 image is available, otherwise PSNR = 0
%     2) y_est (matrix M x N): Final estimate (in the range [0,1])
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2006-2011 Tampere University of Technology.
% All rights reserved.
% This work should only be used for nonprofit purposes.
%
% AUTHORS:
%     Kostadin Dabov, email: dabov _at_ cs.tut.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% In case, a noisy image z is not provided, then use the filename
%%%%  below to read an original image (might contain path also). Later,
%%%%  artificial AWGN noise is added and this noisy image is processed
%%%%  by the BM3D.
%%%%
image_name = [
%     'montage.png'
     'Cameraman256.png'
%     'boat.png'
%     'Lena512.png'
%     'house.png'
%     'barbara.png'
%     'peppers256.png'
%     'fingerprint.png'
%     'couple.png'
%     'hill.png'
%     'man.png'
    ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Quality/complexity trade-off profile selection
%%%%
%%%%  'np' --> Normal Profile (balanced quality)
%%%%  'lc' --> Low Complexity Profile (fast, lower quality)
%%%%
%%%%  'high' --> High Profile (high quality, not documented in [1])
%%%%
%%%%  'vn' --> This profile is automatically enabled for high noise
%%%%           when sigma > 40
%%%%
%%%%  'vn_old' --> This is the old 'vn' profile that was used in [1].
%%%%           It gives inferior results than 'vn' in most cases.
%%%%
if (exist('profile') ~= 1)
    profile         = 'np'; %% default profile
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Specify the std. dev. of the corrupting noise
%%%%
if (exist('sigma') ~= 1),
    sigma               = 25; %% default standard deviation of the AWGN
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Following are the parameters for the Normal Profile.
%%%%

%%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
transform_2D_Wiener_name = 'dct';     %% transform used for the Wiener filt. of size N1_wiener x N1_wiener
transform_3rd_dim_name   = 'haar';    %% transform used in the 3-rd dim, the same for HT and Wiener filt.

%%%% Hard-thresholding (HT) parameters:
N1                  = 8;   %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
Nstep               = 3;   %% sliding step to process every next reference block
N2                  = 16;  %% maximum number of similar blocks (maximum size of the 3rd dimension of a 3D array)
Ns                  = 39;  %% length of the side of the search neighborhood for full-search block-matching (BM), must be odd
tau_match           = 3000;%% threshold for the block-distance (d-distance)
lambda_thr2D        = 0;   %% threshold parameter for the coarse initial denoising used in the d-distance measure
lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D transform domain
beta                = 2.0; %% parameter of the 2D Kaiser window used in the reconstruction

%%%% Wiener filtering parameters:
N1_wiener           = 8;
Nstep_wiener        = 3;
N2_wiener           = 32;
Ns_wiener           = 39;
tau_match_wiener    = 400;
beta_wiener         = 2.0;

%%%% Block-matching parameters:
stepFS              = 1;  %% step that forces to switch to full-search BM, "1" implies always full-search
smallLN             = 'not used in np'; %% if stepFS > 1, then this specifies the size of the small local search neighb.
stepFSW             = 1;
smallLNW            = 'not used in np';
thrToIncStep        = 8;  % if the number of non-zero coefficients after HT is less than thrToIncStep,
                          % than the sliding step to the next reference block is incresed to (nm1-1)

if strcmp(profile, 'lc') == 1,

    Nstep               = 6;
    Ns                  = 25;
    Nstep_wiener        = 5;
    N2_wiener           = 16;
    Ns_wiener           = 25;

    thrToIncStep        = 3;
    smallLN             = 3;
    stepFS              = 6*Nstep;
    smallLNW            = 2;
    stepFSW             = 5*Nstep_wiener;

end

% Profile 'vn' was proposed in
%  Y. Hou, C. Zhao, D. Yang, and Y. Cheng, 'Comment on "Image Denoising by Sparse 3D Transform-Domain
%  Collaborative Filtering"', accepted for publication, IEEE Trans. on Image Processing, July, 2010.
% as a better alternative to that initially proposed in [1] (which is currently in profile 'vn_old')
if (strcmp(profile, 'vn') == 1) | (sigma > 40),

    N2                  = 32;
    Nstep               = 4;

    N1_wiener           = 11;
    Nstep_wiener        = 6;

    lambda_thr3D        = 2.8;
    thrToIncStep        = 3;
    tau_match_wiener    = 3500;
    tau_match           = 25000;

    Ns_wiener           = 39;

end

% The 'vn_old' profile corresponds to the original parameters for strong noise proposed in [1].
if (strcmp(profile, 'vn_old') == 1) & (sigma > 40),

    transform_2D_HT_name = 'dct';

    N1                  = 12;
    Nstep               = 4;

    N1_wiener           = 11;
    Nstep_wiener        = 6;

    lambda_thr3D        = 2.8;
    lambda_thr2D        = 2.0;
    thrToIncStep        = 3;
    tau_match_wiener    = 3500;
    tau_match           = 5000;

    Ns_wiener           = 39;

end

decLevel = 0;        %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
thr_mask = ones(N1); %% N1xN1 mask of threshold scaling coeff. --- by default there is no scaling, however the use of different thresholds for different wavelet decompoistion subbands can be done with this matrix

if strcmp(profile, 'high') == 1, %% this profile is not documented in [1]

    decLevel     = 1;
    Nstep        = 2;
    Nstep_wiener = 2;
    lambda_thr3D = 2.5;
    vMask = ones(N1,1); vMask((end/4+1):end/2)= 1.01; vMask((end/2+1):end) = 1.07; %% this allows to have different threhsolds for the finest and next-to-the-finest subbands
    thr_mask = vMask * vMask';
    beta         = 2.5;
    beta_wiener  = 1.5;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Note: touch below this point only if you know what you are doing!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Check whether to dump information to the screen or remain silent
dump_output_information = 1;
if (exist('print_to_screen') == 1) & (print_to_screen == 0),
    dump_output_information = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Create transform matrices, etc.
%%%%
[Tfor, Tinv]   = getTransfMatrix(N1, transform_2D_HT_name, decLevel);     %% get (normalized) forward and inverse transform matrices
[TforW, TinvW] = getTransfMatrix(N1_wiener, transform_2D_Wiener_name, 0); %% get (normalized) forward and inverse transform matrices

if (strcmp(transform_3rd_dim_name, 'haar') == 1) | (strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
    %%% If Haar is used in the 3-rd dimension, then a fast internal transform is used, thus no need to generate transform
    %%% matrices.
    hadper_trans_single_den         = {};
    inverse_hadper_trans_single_den = {};
else
    %%% Create transform matrices. The transforms are later applied by
    %%% matrix-vector multiplication for the 1D case.
    for hpow = 0:ceil(log2(max(N2,N2_wiener))),
        h = 2^hpow;
        [Tfor3rd, Tinv3rd]   = getTransfMatrix(h, transform_3rd_dim_name, 0);
        hadper_trans_single_den{h}         = single(Tfor3rd);
        inverse_hadper_trans_single_den{h} = single(Tinv3rd');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2D Kaiser windows used in the aggregation of block-wise estimates
%%%%
if beta_wiener==2 & beta==2 & N1_wiener==8 & N1==8 % hardcode the window function so that the signal processing toolbox is not needed by default
    Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
        0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
        0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
        0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
        0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
        0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
        0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
        0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924];
    Wwin2D_wiener = Wwin2D;
else
    Wwin2D           = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
    Wwin2D_wiener    = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the aggregation of the Wiener filt. part
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% If needed, read images, generate noise, or scale the images to the
%%%% [0,1] interval
%%%%
if (exist('y') ~= 1) | (exist('z') ~= 1)
    y        = im2double(imread(image_name));  %% read a noise-free image and put in intensity range [0,1]
    randn('seed', 0);                          %% generate seed
    z        = y + (sigma/255)*randn(size(y)); %% create a noisy image
else  % external images

    image_name = 'External image';

    % convert z to double precision if needed
    z = double(z);

    % convert y to double precision if needed
    y = double(y);

    % if z's range is [0, 255], then convert to [0, 1]
    if (max(z(:)) > 10), % a naive check for intensity range
        z = z / 255;
    end

    % if y's range is [0, 255], then convert to [0, 1]
    if (max(y(:)) > 10), % a naive check for intensity range
        y = y / 255;
    end
end

if (size(z,3) ~= 1) | (size(y,3) ~= 1),
    error('BM3D accepts only grayscale 2D images.');
end

% Check if the true image y is a valid one; if not, then we cannot compute PSNR, etc.
y_is_invalid_image = (length(size(z)) ~= length(size(y))) | (size(z,1) ~= size(y,1)) | (size(z,2) ~= size(y,2));
if (y_is_invalid_image),
    dump_output_information = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Print image information to the screen
%%%%
if dump_output_information == 1,
    fprintf('Image: %s (%dx%d), sigma: %.1f\n', image_name, size(z,1), size(z,2), sigma);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Step 1. Produce the basic estimate by HT filtering
%%%%
tic;
y_hat = bm3d_thr(z, hadper_trans_single_den, Nstep, N1, N2, lambda_thr2D,...
	lambda_thr3D, tau_match*N1*N1/(255*255), (Ns-1)/2, (sigma/255), thrToIncStep, single(Tfor), single(Tinv)', inverse_hadper_trans_single_den, single(thr_mask), Wwin2D, smallLN, stepFS );
estimate_elapsed_time = toc;

if dump_output_information == 1,
    PSNR_INITIAL_ESTIMATE = 10*log10(1/mean((y(:)-double(y_hat(:))).^2));
    fprintf('BASIC ESTIMATE, PSNR: %.2f dB\n', PSNR_INITIAL_ESTIMATE);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Step 2. Produce the final estimate by Wiener filtering (using the
%%%%  hard-thresholding initial estimate)
%%%%
tic;
y_est = bm3d_wiener(z, y_hat, hadper_trans_single_den, Nstep_wiener, N1_wiener, N2_wiener, ...
    'unused arg', tau_match_wiener*N1_wiener*N1_wiener/(255*255), (Ns_wiener-1)/2, (sigma/255), 'unused arg', single(TforW), single(TinvW)', inverse_hadper_trans_single_den, Wwin2D_wiener, smallLNW, stepFSW, single(ones(N1_wiener)) );
wiener_elapsed_time = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Calculate the final estimate's PSNR, print it, and show the
%%%% denoised image next to the noisy one
%%%%
y_est = double(y_est);

PSNR = 0; %% Remains 0 if the true image y is not available
if (~y_is_invalid_image), % checks if y is a valid image
    PSNR = 10*log10(1/mean((y(:)-y_est(:)).^2)); % y is valid
end

if dump_output_information == 1,
    fprintf('FINAL ESTIMATE (total time: %.1f sec), PSNR: %.2f dB\n', ...
        wiener_elapsed_time + estimate_elapsed_time, PSNR);

    figure, imshow(z); title(sprintf('Noisy %s, PSNR: %.3f dB (sigma: %d)', ...
        image_name(1:end-4), 10*log10(1/mean((y(:)-z(:)).^2)), sigma));

    figure, imshow(y_est); title(sprintf('Denoised %s, PSNR: %.3f dB', ...
        image_name(1:end-4), PSNR));

end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some auxiliary functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
%
% Create forward and inverse transform matrices, which allow for perfect
% reconstruction. The forward transform matrix is normalized so that the
% l2-norm of each basis element is 1.
%
% [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
%
%  INPUTS:
%
%   N               --> Size of the transform (for wavelets, must be 2^K)
%
%   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
%                       listed by 'help wfilters' (bi-orthogonal wavelets)
%                       'DCrand' -- an orthonormal transform with a DC and all
%                       the other basis elements of random nature
%
%   dec_levels      --> If a wavelet transform is generated, this is the
%                       desired decomposition level. Must be in the
%                       range [0, log2(N)-1], where "0" implies
%                       full decomposition.
%
%  OUTPUTS:
%
%   Tforward        --> (N x N) Forward transform matrix
%
%   Tinverse        --> (N x N) Inverse transform matrix
%

if exist('dec_levels') ~= 1,
    dec_levels = 0;
end

if N == 1,
    Tforward = 1;
elseif strcmp(transform_type, 'hadamard') == 1,
    Tforward    = hadamard(N);
elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
    Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
       0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
       0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
      -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
       0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
                       0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
                       0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
                       0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
    Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
       0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
       0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
       0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
       0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
       0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
       0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
       0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
    Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
       0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
       0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
       0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
       0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
       0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
       0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
       0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
elseif strcmp(transform_type, 'dct') == 1,
    Tforward    = dct(eye(N));
elseif strcmp(transform_type, 'dst') == 1,
    Tforward    = dst(eye(N));
elseif strcmp(transform_type, 'DCrand') == 1,
    x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
    if (Q(1) < 0),
        Q = -Q;
    end;
    Tforward = Q';
else %% a wavelet decomposition supported by 'wavedec'
    %%% Set periodic boundary conditions, to preserve bi-orthogonality
    dwtmode('per','nodisp');

    Tforward = zeros(N,N);
    for i = 1:N
        Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
    end
end

%%% Normalize the basis elements
Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';

%%% Compute the inverse transform matrix
Tinverse = inv(Tforward);

return;

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


text:
%--------------------------------------------------------------------------------------------------------
% The system is created based on the principles described in the following paper
% Jimmy SJ. Ren, Li Xu, Qiong Yan, Wenxiu Sun, "Shepard Convolutional Neural Networks",
% Advances in Neural Information Processing Systems (NIPS 2015)
% email: jimmy.sj.ren@gmail.com
%--------------------------------------------------------------------------------------------------------
addpath applications/Shepard_CNN/Shepard_super_res/
addpath applications/Shepard_CNN/utility/
addpath applications/image_denoise/utility/
addpath utils/
addpath cuda/
addpath mem/
addpath layers/
addpath layers_adapters/
addpath pipeline/

clearvars -global config;
clearvars -global mem;
clear gen_mask_patch_cat_idx_for_super_res;
clear;
global config;

% set5
% -- baby bird butterfly head woman
I = im2double(imread('applications/Shepard_CNN/Shepard_super_res/images/Set5/butterfly_GT.bmp'));
% set 14
% -- baboon barbara bridge coastguard comic face flowers foreman lenna man monarch pepper ppt3 zebra
%I = im2double(imread('applications/Shepard_CNN/Shepard_super_res/images/Set14/pepper.bmp'));
if(size(I,3) == 3)
    I = rgb2ycbcr(I);
end

I = modcrop(I, 3);
I = padarray(I, [3 3], 'replicate');

% downsample the image
D3chs = imresize(I, [size(I,1)/3 size(I,2)/3], 'bicubic');
D = D3chs(:,:,1);

% super-resolution x3
U = imresize(D, [size(I,1) size(I,2)], 'nearest');
prepare_sr_net(size(U, 1), size(U, 2), 'applications/Shepard_CNN/Shepard_super_res/results/shepard_layer_x3/sr_x3.mat');
mask = config.NEW_MEM([1 0 0; ...
                       0 0 0; ...
                       0 0 0]);
mask = repmat(mask, size(U, 1)/3, size(U, 2)/3, 1);

final_output = apply_net_with_mask(config.NEW_MEM(U) .* mask, mask);
final_output = gather(final_output);
D3chs = imresize(D3chs, [size(final_output,1),size(final_output,2)]);
final_output = shave(final_output, [3 3]);
D3chs = shave(D3chs, [3 3]);

final_output = cat(3, final_output, D3chs(:,:,2:3));
final_output = ycbcr2rgb(double(final_output));
I_bicubic = ycbcr2rgb(D3chs);

figure('name', 'bicubic'), imshow(I_bicubic);
figure('name', 'Sheperd CNN'), imshow(final_output);
drawnow();

author: jimmy-ren


text:
function map = tab20b(N)
% Qualitative colormap from MatPlotLib, for plots using the line ColorOrder.
% In MatPlotLib 2 it is named VEGA20b, for MatPlotLib 3 was renamed TAB20b.
%
% Copyright (c) 2017-2019 Stephen Cobeldick
%
%%% Syntax:
%  map = tab20b
%  map = tab20b(N)
%
% For MatPlotLib 2.0 improved colormaps were created for plot lines of
% categorical data. The new colormaps are introduced here:
% <http://matplotlib.org/2.0.0rc2/users/dflt_style_changes.html>
% VEGA10/TAB10 is the default Line Color Order for MatPlotLib 2 and 3.
%
% MATLAB axes ColorOrder (note that this is NOT the axes COLORMAP):
% <https://www.mathworks.com/help/matlab/creating_plots/defining-the-color-of-lines-for-plotting.html>
% <https://www.mathworks.com/help/matlab/graphics_transition/why-are-plot-lines-different-colors.html>
%
%% Examples %%
%
%%% PLOT using matrices:
% N = 20;
% axes('ColorOrder',tab20b(N),'NextPlot','replacechildren')
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*sin(x+2*n*pi/N), X(:), 1:N);
% plot(X,Y, 'linewidth',4)
%
%%% PLOT in a loop:
% N = 20;
% set(0,'DefaultAxesColorOrder',tab20b(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*sin(x+2*n*pi/N), X(:), 1:N);
% for n = 1:N
%     plot(X(:),Y(:,n), 'linewidth',4);
%     hold all
% end
%
%%% LINE using matrices:
% N = 20;
% set(0,'DefaultAxesColorOrder',tab20b(N))
% X = linspace(0,pi*3,1000);
% Y = bsxfun(@(x,n)sqrt(n)*cos(x+2*n*pi/N), X(:), 1:N);
% line(X(:),Y)
%
%% Input and Output Arguments %%
%
%%% Inputs (*=default):
% N = NumericScalar, N>=0, an integer to define the colormap length.
%   = *[], use the length of the current figure's colormap (see COLORMAP).
%
%%% Outputs:
% map = NumericMatrix, size Nx3, a colormap of RGB values between 0 and 1.
%
% See also TAB10 TAB20 TAB20C SET TWILIGHT VIRIDIS LINES COLORMAP PARULA

if nargin<1
	N = size(get(gcf,'colormap'),1);
else
	assert(isscalar(N)&&isreal(N),'First argument must be a real numeric scalar.')
	assert(fix(N)==N&&N>=0,'First argument must be a positive integer.')
end
%
hex = ['#393b79';'#5254a3';'#6b6ecf';'#9c9ede';'#637939';'#8ca252';'#b5cf6b';'#cedb9c';'#8c6d31';'#bd9e39';'#e7ba52';'#e7cb94';'#843c39';'#ad494a';'#d6616b';'#e7969c';'#7b4173';'#a55194';'#ce6dbd';'#de9ed6'];
raw = sscanf(hex.','#%2x%2x%2x',[3,Inf]).';
%
map = raw(1+mod(0:N-1,size(raw,1)),:) / 255;
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tab20b


text:
%% Efficient subpixel image registration by cross-correlation.
% Registers two images (2-D rigid translation) within a  fraction
% of a pixel specified by the user. Instead of computing a zero-padded FFT
% (fast Fourier transform), this code uses selective upsampling by a
% matrix-multiply DFT (discrete FT) to dramatically reduce computation time and memory
% without sacrificing accuracy. With this procedure all the image points are used to
% compute the upsampled cross-correlation in a very small neighborhood around its peak. This
% algorithm is referred to as the single-step DFT algorithm in [1].
%
% [1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
% "Efficient subpixel image registration algorithms," Opt. Lett. 33,
% 156-158 (2008).
%
% -----------------------------------------------------------------------
%
% Copyright (c) 2016, Manuel Guizar Sicairos, James R. Fienup, University of Rochester
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of the University of Rochester nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
% --------------------------------------------------------------------------

%% Syntax
% The code receives the FFT of the reference and the shifted images, and an
% (integer) upsampling factor. The code expects FFTs with DC in (1,1) so do not use
% fftshift.
%
%    output = dftregistration(fft2(f),fft2(g),usfac);
%
% The images are registered to within 1/usfac of a pixel.
%
% output(1) is the normalized root-mean-squared error (NRMSE) [1] between f and
% g.
%
% output(2) is the global phase difference between the two images (should be
% zero if images are real-valued and non-negative).
%
% output(3) and output(4) are the row and column shifts between f and g respectively.
%
%    [output Greg] = dftregistration(fft2(f),fft2(g),usfac);
%
% Greg is an optional output, it returns the Fourier transform of the registered version of g,
% where the global phase difference [output(2)] is also compensated.

%% Obtain a reference and shifted images
% To illustrate the use of the algorithm, lets obtain a reference and a
% shifted image. First we read the reference image f(x,y)
f = im2double(imread('cameraman.tif'));

%%
% Define g(x,y) as a version of f(x,y) shifted by fractional values of a
% pixel and multiplied by a global phase.
deltar = -3.48574;
deltac = 8.73837;
phase = 2;
[nr,nc]=size(f);
Nr = ifftshift((-fix(nr/2):ceil(nr/2)-1));
Nc = ifftshift((-fix(nc/2):ceil(nc/2)-1));
[Nc,Nr] = meshgrid(Nc,Nr);
g = ifft2(fft2(f).*exp(1i*2*pi*(deltar*Nr/nr+deltac*Nc/nc))).*exp(-1i*phase);
figure(1);
subplot(1,2,1);
imshow(abs(f));
title('Reference image, f(x,y)')
subplot(1,2,2);
imshow(abs(g));
title('Shifted image, g(x,y)')
%%
% We have shifted the image by 8.73837 and -3.48574 pixels in the x and
% y direction, respectively, and added a phase of 2 radians to g(x,y). The shift
% was implemented by applying a linear phase on its
% FT, thus we have assumed that the images wrap around (features leaving one
% side of the window reappear on the opposite side) and that the image is band-limited
% (interpolated by a sinc function). Cross-correlation image registration by DFTs
% (both the matrix-multiply DFT and the zero-padded FFT) share these assumptions.
%
% This registration technique is well suited to compare images that are captured
% in Fourier domain (i.e. to evaluate an image reconstruction by holography
% or phase retrieval) which are strictly band-limited and exhibit the
% wrap-around effect.
%
% Even though the registration code assumes band-limited images that wrap around, we
% have obtained very good results when applying it to
% band-limited microscope images, and aliased imagery. That is when shifting
% the image brings in new content instead of wrapping it around or when the
% images are not band-limited.

%% Sample Image Registration
% dftregistration.m receives the FT of f and g and the upsampling factor.
% The code expects DC of the FTs at (1,1) so don't use fftshift.
%
% We now use the image registration code to register f and g within 0.01
% pixels by specifying an upsampling parameter of 100
usfac = 100;
[output, Greg] = dftregistration(fft2(f),fft2(g),usfac);
display(output),

%%
% The pixel shift error (difference between the true and obtained shifts)
% is 0.0016 and 0.0043 in the x and y directions respectively. Well within
% the expected accuracy of 0.01. Notice that using the conventional zero-padded
% FFT approach with the same accuracy, would
% require computation of a 25,600x25,600 FFT, which would require more than
% 19 Gbytes of RAM and a very comfortable chair.
%
% The following plot shows the reference image and the registered image.
figure(1);
subplot(1,2,1);
imshow(abs(f));
title('Reference image, f(x,y)')
subplot(1,2,2);
imshow(abs(ifft2(Greg)));
title('Registered image, gr(x,y)')
%% Disclaimer
% I have made every effort to evaluate the proper working of this code
% under many different conditions. However, it is the responsibility of
% the user to ensure that this registration code is adequate and working
% correcntly for their application.
%
% Feel free to e-mail me with questions or comments.

author: AbdoKamel

repo: sidd-ground-truth-image-estimation


text:
function option_num = p18_option_num ( )

%*****************************************************************************80
%
%% P18_OPTION_NUM returns the number of options for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer OPTION_NUM, the number of options.
%
  option_num = 1;

  return
end


text:
function g = grad(sF, v) % gradient

if nargin == 2 % direct evaluation
  v = v(:);
  bario = sF.tri.calcBario(v);
  g = vector3d(zeros(3, length(v)));
  for i = 1:length(v)
    I = find(bario(i, :));
    if length(I) == 3
      f = sF.values(I);
      v = sF.vertices(I);
      g(i) = vector3d(v.xyz \ (f(:)-(f(1)+f(2)+f(3))/3*ones(3, 1)));
    end
  end
else % return S2VectorField
  mp = sF.tri.midPoints;
  g = S2VectorFieldTri(mp, sF.grad(mp));
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [out,Xt,str,ts] = rfnn_mimo_scatter(t,Xt,u,flag,itaVector,alphaVector, NumInVars,NumInTerms,NumOutVars,x0,T)

% This program is an implementation of the on-line RFNN (MIMO) system.
% The structure of the network is determined by the user.
% The input space is partitiond using the scatter-type method.
% All parameters of the network are estimated by Gradient Descent (GD)
% through error backpropagation.

    ninp = NumInVars;
    nout = NumOutVars;
   ninps = ninp+nout+1;  % number of inputs to sfunction [ x y LE ]
   NumRules = NumInTerms;  % Scatter-Type Input Space Partitioning.
       ns = 4*NumInVars*NumInTerms + NumOutVars*NumRules;
     nds = 3*NumInVars*NumInTerms + NumOutVars*NumRules;
     % Learning Rates
     ita1 = itaVector(1); ita2 = itaVector(2);
     ita3 = itaVector(3); ita4 = itaVector(4);
     % Momentum Constants.
     alpha1 = alphaVector(1); alpha2 = alphaVector(2);
     alpha3 = alphaVector(3); alpha4 = alphaVector(4);
%  ----------------------- % initial informations --------------
if abs(flag)==0

    out = [0,ns+nds,nout+ns+nds,ninps,0,1,1];    % states, outputs, inputs, ?, df, #ts
    str = [];                                 % API block consistency
     ts = T;                                  % sample time
    Xt = x0;
%  ----------------------- % state derivatives -----------------
elseif abs(flag) == 2

          x = u(1:ninp);
          e = u(ninp+1:ninp+nout);
   learning = u(ninp+nout+1);

if learning == 1

   % Unroll the states:
   off=1;
   off_end=NumInVars*NumInTerms;
   mean2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   sigma2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Theta2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumOutVars*NumRules-1;
   W = reshape(Xt(off:off_end),NumOutVars,NumRules);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Out2 = reshape(Xt(off:off_end),NumInVars,NumInTerms);

   % Unroll the differential states:
   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   dmean2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   dsigma2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   dTheta2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumOutVars*NumRules-1;
   dW = reshape(Xt(off:off_end),NumOutVars,NumRules);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                    FEEDFORWARD OPERATION                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LAYER 2 - INPUT TERM NODES
  Out2_pr = Out2;
  In2 = x*ones(1,NumInTerms) + Out2_pr.*Theta2;
  Out2 = exp(-((In2-mean2)./sigma2).^2);

% LAYER 3 - RULE (PRODUCT) NODES
 precond = Out2.';
 Out3 = prod(precond,2);

%%%%%%%%%%%% END OF NETWORK FUNCTIONALITY SECTION %%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%			 					                 PARAMETER LEARNING SECTION	                       			                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BACKWARD PASS. Error Backpropagation
% LAYER 3
delta3 = e.'*W;

% LAYER 2
ThetaE = zeros(NumInVars,NumInTerms);
for i=1:NumInVars
     for j=1:NumInTerms
   		  ThetaE(i,j) = (Out3(j)/Out2(i,j))*delta3(j);
     end
 end

% LAYER 2 PARAMETER ADJUSTMENT BY GRADIENT DESCENT.
deltamean2  =  2*ThetaE.*Out2.*(In2-mean2)./((sigma2).^2);
deltasigma2 =  2*ThetaE.*Out2.*((In2-mean2).^2)./(sigma2.^3);
deltaTheta2 = -2*ThetaE.*Out2.*(In2-mean2).*Out2_pr./sigma2.^2;

dmean2 = ita2*deltamean2 + alpha2*dmean2;
 mean2  = mean2 + dmean2;

dsigma2 = ita3*deltasigma2 + alpha3*dsigma2;
  sigma2 = sigma2 + dsigma2;

dTheta2 = ita4*deltaTheta2 + alpha4*dTheta2;
 Theta2  = Theta2 + dTheta2;

% LAYER 4 PARAMETER ADJUSTMENT
 deltaW = e*Out3.';
 dW = ita1*deltaW + alpha1*dW;
  W = W + dW;

%%%%%%%%%%%   END OF PARAMETER LEARNING PROCESS %%%%%%%%%
% State Vector Storage.
% Xt = [mean2 sigma2 Theta2 W Out2 dmean2 dsigma2 dTheta2 dW];

Xt = [ reshape(mean2,NumInVars*NumInTerms,1);
          reshape(sigma2,NumInVars*NumInTerms,1);
          reshape(Theta2,NumInVars*NumInTerms,1);
          reshape(W,NumOutVars*NumRules,1);
          reshape(Out2,NumInVars*NumInTerms,1);
          reshape(dmean2,NumInVars*NumInTerms,1);
          reshape(dsigma2,NumInVars*NumInTerms,1);
          reshape(dTheta2,NumInVars*NumInTerms,1);
          reshape(dW,NumOutVars*NumRules,1);];
end

out=Xt;

%  ----------------------- % outputs -------------------------
elseif flag == 3

  % Unpack the network's parameters first...
   off=1;
   off_end=NumInVars*NumInTerms;
   mean2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   sigma2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Theta2 =reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumOutVars*NumRules - 1;
   W = reshape(Xt(off:off_end),NumOutVars,NumRules);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Out2 = reshape(Xt(off:off_end),NumInVars,NumInTerms);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                                                       FEEDFORWARD OPERATION                                                       %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % LAYER 2 - INPUT TERM NODES
   x = u(1:ninp);

 Out2_pr = Out2;
 In2 = x*ones(1,NumInTerms) + Out2_pr.*Theta2;
 Out2 = exp(-((In2-mean2)./sigma2).^2);

% LAYER 3 - RULE (PRODUCT) NODES
precond = Out2.';
 Out3 = prod(precond,2);

 % LAYER 4
  outact = W*Out3;

  % Block Outputs Vector Formation.
   out=[outact;Xt];

else
   out=[];
end

author: Sable

repo: mcbench-benchmarks