text:
classdef SessClass < TreeNodeClass

    properties % (Access = private)
        runs;
    end

    methods

        % ----------------------------------------------------------------------------------
        function obj = SessClass(varargin)
            obj@TreeNodeClass(varargin);

            obj.type  = 'sess';
            obj.runs = RunClass().empty;
            if nargin==0
                obj.name  = '';
                return;
            end

            if nargin<4
                if isa(varargin{1}, 'SessClass')
                    if nargin==1
                        obj.Copy(varargin{1});
                    elseif nargin==2
                        obj.Copy(varargin{1}, varargin{2});
                    end
                    return;
                elseif isa(varargin{1}, 'FileClass')
                    [~, ~, obj.name] = varargin{1}.ExtractNames();
                else
                    obj.name = varargin{1};
                end
            elseif nargin==4
                if ~isa(varargin{1}, 'FileClass')
                    [~, obj.name] = varargin{1}.ExtractNames();
                else
                    obj.name = varargin{1};
                end
                obj.iGroup = varargin{2};
                obj.iSubj = varargin{3};
                obj.iSess = varargin{4};
            end
        end

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            for ii = 1:length(obj.runs)
                nbytes = nbytes + obj.runs(ii).MemoryRequired();
            end
            nbytes = nbytes + obj.procStream.MemoryRequired();
        end

        % ----------------------------------------------------------------------------------
        % Copy processing params (procInut and procResult) from
        % S to obj if obj and S are equivalent nodes
        % ----------------------------------------------------------------------------------
        function Copy(obj, obj2, conditional)
            % Copy SessClass object obj2 to SessClass object obj. Conditional option applies
            % only to all the runs under this group. If == 'conditional' ONLY derived data,
            % that is, only from procStream but NOT from acquired data is copied for all the runs.
            %
            % Conversly unconditional copy copies all properties in the runs under this session
            if nargin==3 && strcmp(conditional, 'conditional')
                if obj.Mismatch(obj2)
                    return
                end
                for i = 1:length(obj.runs)
                    j = obj.existRun(i,obj2);
                    if (j>0)
                        obj.runs(i).Copy(obj2.runs(j), 'conditional');
                    elseif i<=length(obj2.runs)
                        obj.runs(i).Copy(obj2.runs(i), 'conditional');
                    else
                        obj.Mismatch();
                    end
                end
                obj.Copy@TreeNodeClass(obj2, 'conditional');
            else
                if nargin<3
                    conditional = '';
                end
                for i=1:length(obj2.runs)
                    obj.runs(i) = RunClass(obj2.runs(i), conditional);
                end
                obj.Copy@TreeNodeClass(obj2);
            end
        end

        % --------------------------------------------------------------
        function CopyStims(obj, obj2)
            obj.CondNames = obj2.CondNames;
            for ii = 1:length(obj.runs)
                obj.runs(ii).CopyStims(obj2.runs(ii));
            end
        end

        % ----------------------------------------------------------------------------------
        % Check whether run R exists in this session and return
        % its index if it does exist. Else return 0.
        % ----------------------------------------------------------------------------------
        function j = existRun(obj, k, S)
            j=0;
            for i=1:length(S.runs)
                [~,rname1] = fileparts(obj.runs(k).name);
                [~,rname2] = fileparts(S.runs(i).name);
                if strcmp(rname1, rname2)
                    j=i;
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        % Subjects obj1 and obj2 are considered equivalent if their names
        % are equivalent and their sets of runs are equivalent.
        % ----------------------------------------------------------------------------------
        function B = equivalent(obj1, obj2)
            B=1;
            if ~strcmp(obj1.name, obj2.name)
                B=0;
                return;
            end
            for i=1:length(obj1.runs)
                j = existRun(obj1, i, obj2);
                if j==0
                    B=0;
                    return;
                end
            end
            for i=1:length(obj2.runs)
                j = existRun(obj2, i, obj1);
                if j==0
                    B=0;
                    return;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function b = CheckForNameConflict(obj, run)
            b = false;
            [~, runname] = fileparts(run.GetName());
            if strcmp(runname, obj.name)
                b = true;
            end
        end

        % ----------------------------------------------------------------------------------
        function Add(obj, run)
            [~,f,e] = fileparts(run.GetName());
            if strcmp(f, obj.name)
                msg{1} = sprintf('WARNING: The run being added (%s) has the same name as the session (%s) containing it. ', [f,e], obj.name);
                msg{2} = sprintf('The run names should not have the same name as the session, otherwise ');
                msg{3} = sprintf('it may cause incorrect results in processing. Do you want to change this run''s name?');
                obj.logger.Write('%s\n', [msg{:}]);

                % Deconflict name of session with name of run if there is no
                % session folder for this session
                if ~ispathvalid(['./',obj.name], 'dir')
                    obj.name = [obj.name, '_s'];
                    obj.logger.Write('Renaming session to %s\n', obj.name);
                end
            end

            % Add run to this session
            jj = 0;
            for ii = 1:length(obj.runs)
                if strcmp(obj.runs(ii).GetName, run.GetName())
                    jj = ii;
                    break;
                end
            end
            if jj==0
                jj = length(obj.runs)+1;
                run.SetIndexID(obj.iGroup, obj.iSubj, obj.iSess, jj);
                run.SetPath(obj.path);                      % Inherit root path from session
                obj.runs(jj) = run;
                obj.logger.Write('         Added run  "%s"  to session  "%s" .\n', obj.runs(jj).GetFileName, obj.GetFileName);
            end
            obj.children = obj.runs;
        end

        % ----------------------------------------------------------------------------------
        function list = DepthFirstTraversalList(obj)
            list{1} = obj;
            for ii = 1:length(obj.runs)
                list{ii+1,1} = obj.runs(ii);
            end
        end

        % ----------------------------------------------------------------------------------
        % Deletes derived data in procResult
        % ----------------------------------------------------------------------------------
        function Reset(obj, option)
            if ~exist('option','var')
                option = 'down';
            end
            if strcmp(option, 'down')
                for jj = 1:length(obj.runs)
                    obj.runs(jj).Reset();
                end
            end
            Reset@TreeNodeClass(obj);
        end

        % ----------------------------------------------------------------------------------
        function err = LoadSubBranch(obj)
            err = -1;
            if isempty(obj)
                return;
            end
            err1 = obj.procStream.Load([obj.path, obj.GetOutputFilename()]);
            err2 = obj.runs(1).Load();
            if err1==0 && err2==0
                err = 0;
            end
        end

        % ----------------------------------------------------------------------------------
        function FreeMemorySubBranch(obj)
            if isempty(obj)
                return;
            end
            obj.runs(1).FreeMemory()
        end

        % ----------------------------------------------------------------------------------
        function CreateOutputDir(obj)
            if ispathvalid([obj.pathOutputAlt, obj.outputDirname, obj.name])
                return;
            end
            if ~ispathvalid([obj.pathOutputAlt, obj.name])
                return;
            end
            mkdir([obj.pathOutputAlt, obj.outputDirname, obj.name]);
        end

        % ----------------------------------------------------------------------------------
        function LoadInputVars(obj, tHRF_common)
			obj.inputVars.nTrialsRuns = [];
            for iRun = 1:length(obj.runs)
                % Set common tHRF: make sure size of tHRF, dcAvg and dcAvg is same for
                % all runs. Use smallest tHRF as the common one.
                obj.runs(iRun).procStream.output.SettHRFCommon(tHRF_common, obj.runs(iRun).name, obj.runs(iRun).type);

                obj.inputVars.dodAvgRuns{obj.runs(iRun).iRun}    = obj.runs(iRun).procStream.output.GetVar('dodAvg');
                obj.inputVars.dodAvgStdRuns{obj.runs(iRun).iRun} = obj.runs(iRun).procStream.output.GetVar('dodAvgStd');
                obj.inputVars.dodSum2Runs{obj.runs(iRun).iRun}   = obj.runs(iRun).procStream.output.GetVar('dodSum2');
                obj.inputVars.dcAvgRuns{obj.runs(iRun).iRun}     = obj.runs(iRun).procStream.output.GetVar('dcAvg');
                obj.inputVars.dcAvgStdRuns{obj.runs(iRun).iRun}  = obj.runs(iRun).procStream.output.GetVar('dcAvgStd');
                obj.inputVars.dcSum2Runs{obj.runs(iRun).iRun}    = obj.runs(iRun).procStream.output.GetVar('dcSum2');
                obj.inputVars.tHRFRuns{obj.runs(iRun).iRun}      = obj.runs(iRun).procStream.output.GetTHRF();
                obj.inputVars.mlActRuns{obj.runs(iRun).iRun}     = obj.runs(iRun).procStream.output.GetVar('mlActAuto');
                obj.inputVars.nTrialsRuns{obj.runs(iRun).iRun}   = obj.runs(iRun).procStream.output.GetVar('nTrials');
                if ~isempty(obj.runs(iRun).procStream.output.GetVar('misc'))
                    if isfield(obj.runs(iRun).procStream.output.misc, 'stim') == 1
                        obj.inputVars.stimRuns{obj.runs(iRun).iRun}      = obj.runs(iRun).procStream.output.misc.stim;
                    else
                        obj.inputVars.stimRuns{obj.runs(iRun).iRun}      = obj.runs(iRun).GetVar('stim');
                    end
                else
                    obj.inputVars.stimRuns{obj.runs(iRun).iRun}      = obj.runs(iRun).GetVar('stim');
                end
                obj.inputVars.AauxRuns{obj.runs(iRun).iRun}      = obj.runs(iRun).procStream.output.GetVar('Aaux');
                obj.inputVars.tIncAutoRuns{obj.runs(iRun).iRun}  = obj.runs(iRun).procStream.output.GetVar('tIncAuto');
                obj.inputVars.rcMapRuns{obj.runs(iRun).iRun}     = obj.runs(iRun).procStream.output.GetVar('rcMap');

                % a) Find all variables needed by proc stream
                args = obj.procStream.GetInputArgs();

                % b) Find these variables in this run
                for ii = 1:length(args)
                    if ~eval( sprintf('isproperty(obj.inputVars, ''%s'')', args{ii}) )
                        eval( sprintf('obj.inputVars.%s = obj.GetVar(args{ii});', args{ii}) );
                    end
                end

                % Free run memory
                obj.runs(iRun).FreeMemory()
            end

        end

        % ----------------------------------------------------------------------------------
        function Calc(obj, options)
            if ~exist('options','var') || isempty(options)
                options = 'overwrite';
            end

            if strcmpi(options, 'overwrite')
                % Recalculating result means deleting old results, if
                % option == 'overwrite'
                obj.procStream.output.Flush();
            end

            if obj.DEBUG
                obj.logger.Write('Calculating processing stream for group %d, session %d\n', obj.iGroup, obj.iSubj);
            end

            % Calculate all runs in this session and generate common tHRF
            tHRF_common = {};
            for iRun = 1:length(obj.runs)
                obj.runs(iRun).Calc();

                % Find smallest tHRF among the sessions and make this the common one.
                tHRF_common = obj.runs(iRun).procStream.output.GeneratetHRFCommon(tHRF_common);
            end

            % Update call application GUI using it's generic Update function
            if ~isempty(obj.updateParentGui)
                obj.updateParentGui('DataTreeClass', [obj.iGroup, obj.iSubj, obj.iSess, obj.iRun]);
            end

            % Load all the variables that might be needed by procStream.Calc() to calculate proc stream for this session
            obj.LoadInputVars(tHRF_common);

            Calc@TreeNodeClass(obj);

            if obj.DEBUG
                fprintf('Completed processing stream for group %d, session %d\n', obj.iGroup, obj.iSubj);
                fprintf('\n');
            end

        end

        % ----------------------------------------------------------------------------------
        function Print(obj, indent)
            if ~exist('indent', 'var')
                indent = 4;
            else
                indent = indent+4;
            end
            Print@TreeNodeClass(obj, indent);
            for ii = 1:length(obj.runs)
                obj.runs(ii).Print(indent);
            end
        end

        % ---------------------------------------------------------------
        function PrintProcStream(obj)
            fcalls = obj.procStream.GetFuncCallChain();
            obj.logger.Write('Session processing stream:\n');
            for ii = 1:length(fcalls)
                obj.logger.Write('%s\n', fcalls{ii});
            end
            obj.logger.Write('\n');
            obj.runs(1).PrintProcStream();
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return;
            end
            for ii = 1:length(obj.runs)
                if ~obj.runs(ii).IsEmpty()
                    b = false;
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmptyOutput(obj)
            b = true;
            if isempty(obj)
                return;
            end
            for ii = 1:length(obj.runs)
                if ~obj.runs(ii).IsEmptyOutput()
                    b = false;
                    break;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function SaveAcquiredData(obj)
            for ii = 1:length(obj.runs)
                obj.runs(ii).SaveAcquiredData();
            end
        end

        % ----------------------------------------------------------------------------------
        function b = AcquiredDataModified(obj)
            b = false;
            for ii = 1:length(obj.runs)
                if obj.runs(ii).AcquiredDataModified()
                    b = true;
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function varval = GetVar(obj, varname)
            % First call the common code for all levels
            varval = obj.GetVar@TreeNodeClass(varname);

            % Now call the session specific part
            if isempty(varval)
                varval = obj.runs(1).GetVar(varname);
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Public Set/Get methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function probe = GetProbe(obj, option)
            probe = obj.runs(1).GetProbe();
%             for run = obj.runs
%                if ~(probe == run.GetProbe())
%                   warning(['Probe ', run.name, 'differs from ', obj.runs(1).name])
%                end
%             end
        end

        % ----------------------------------------------------------------------------------
        function SD = GetSDG(obj,option)
            if exist('option','var')
                SD = obj.runs(1).GetSDG(option);
            else
                SD = obj.runs(1).GetSDG();
            end
        end

        % ----------------------------------------------------------------------------------
        function srcpos = GetSrcPos(obj,option)
            if exist('option','var')
                srcpos = obj.runs(1).GetSrcPos(option);
            else
                srcpos = obj.runs(1).GetSrcPos();
            end
        end

        % ----------------------------------------------------------------------------------
        function detpos = GetDetPos(obj,option)
            if exist('option','var')
                detpos = obj.runs(1).GetDetPos(option);
            else
                detpos = obj.runs(1).GetDetPos();
            end
        end

        % ----------------------------------------------------------------------------------
        function bbox = GetSdgBbox(obj)
            bbox = obj.runs(1).GetSdgBbox();
        end

        % ----------------------------------------------------------------------------------
        function wls = GetWls(obj)
            wls = obj.runs(1).GetWls();
        end

        % ----------------------------------------------------------------------------------
        function [iDataBlks, iCh] = GetDataBlocksIdxs(obj, iCh)
            if nargin<2
                iCh = [];
            end
            [iDataBlks, iCh] = obj.runs(1).GetDataBlocksIdxs(iCh);
        end

        % ----------------------------------------------------------------------------------
        function n = GetDataBlocksNum(obj)
            n = obj.runs(1).GetDataBlocksNum();
        end

        % ----------------------------------------------------------------------------------
        function SetConditions(obj, CondNames)
            if nargin==1
                CondNames = {};
                for ii = 1:length(obj.runs)
                    obj.runs(ii).SetConditions();
                    CondNames = [CondNames, obj.runs(ii).GetConditions()];
                end
            elseif nargin==2
                for ii = 1:length(obj.runs)
                    obj.runs(ii).SetConditions(CondNames);
                end
            end
            obj.CondNames = unique(CondNames);
        end

        % ----------------------------------------------------------------------------------
        function CondNames = GetConditionsActive(obj)
            CondNames = obj.CondNames;
            for ii = 1:length(obj.runs)
                CondNamesRun = obj.runs(ii).GetConditionsActive();
                for jj = 1:length(CondNames)
                    k = find(strcmp(['-- ', CondNames{jj}], CondNamesRun)); %#ok<EFIND>
                    if ~isempty(k)
                        CondNames{jj} = ['-- ', CondNames{jj}];
                    end
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function aux = GetAuxiliary(obj)
            aux = [];
        end

        % ----------------------------------------------------------------------------------
        function [fn_error, missing_args, prereqs] = CheckProcStreamOrder(obj)
            missing_args = {};
            fn_error = 0;
            prereqs = '';
            for i = 1:length(obj.runs)
                [fn_error, missing_args, prereqs] = obj.runs(i).CheckProcStreamOrder;
                if ~isempty(missing_args)
                    return
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function r = ListOutputFilenames(obj, options)
            if ~exist('options','var')
                options = '';
            end
            r = obj.GetOutputFilename(options);
            fprintf('  %s %s\n', obj.path, r);
            for ii = 1:length(obj.runs)
                obj.runs(ii).ListOutputFilenames(options);
            end
        end

        % --------------------------------------------------------------------------
        function ApplyParamEditsToAllRuns(obj, iFcall, iParam, val)
            for jj = 1:length(obj.runs)
                obj.runs(jj).procStream.EditParam(iFcall, iParam, val);
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Private methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods (Access = public)

        % ----------------------------------------------------------------------------------
        function BackwardCompatability(obj)
            obj.BackwardCompatability@TreeNodeClass();
            for ii = 1:length(obj.runs)
                obj.runs(ii).BackwardCompatability();
            end
        end

    end  % Private methods

end

author: BUNPC


text:
function [s,c,r,t]=v_atan2sc(y,x)
%V_ATAN2SC    sin and cosine of atan(y/x) [S,C,R,T]=(Y,X)
%
% Outputs:
%    s    sin(t) where tan(t) = y/x
%    c    cos(t) where tan(t) = y/x
%    r    sqrt(x^2 + y^2)
%    t    arctan of y/x
%
% Note y and x can be arrays but must be the same size. The outputs will
% all be the same size as y.

%      Copyright (C) Mike Brookes 2007
%      Version: $Id: v_atan2sc.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sz=size(y);
s=zeros(sz);
c=NaN(sz);
r=zeros(sz);
t=NaN(sz);
m=y==0;
if any(m(:)) % handle case when y=0 and possibly x=0 also
    t(m)=(x(m)<0);
    c(m)=1-2*t(m);
    r(m)=abs(x(m));
    t(m)=t(m)*pi;
end
m=abs(y)>abs(x) & isnan(c);
if any(m(:))
    q = x(m)./y(m);
    u = sqrt(1+q.^2).*sign(y(m)); % avoids underflow even if x and y are very small
    s(m) = 1./u;
    c(m) = s(m).*q;
    r(m) = y(m).*u;
end
m=isnan(c);
if any(m(:))
    q = y(m)./x(m);
    u = sqrt(1+q.^2).*sign(x(m));
    c(m) = 1./u;
    s(m) = c(m).*q;
    r(m) = x(m).*u;
end
m=isnan(t);
if nargout>3 && any(m(:))
    t(m)=atan2(s(m),c(m));
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtFfmHelmholtzBWdirSPMD.m                    |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Solve dirichlet scatering problem with        |
%|  `---'  |                Brackage-Werner formulation                   |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N   = 1e3
tol = 1e-3
typ = 'P1'
gss = 3
X0  = [0 0 -1]

% Spherical mesh
sphere = mshSphere(N,1);
sigma  = dom(sphere,gss);
figure
plot(sphere)
axis equal

% Radiative mesh
square     = mshSquare(5*N,[5 5]);
square.vtx = [square.vtx(:,1) zeros(size(square.vtx,1),1) square.vtx(:,2)];
hold on
plot(square)

% Frequency adjusted to maximum esge size
stp = sphere.stp;
k   = 1/stp(2)
f   = (k*340)/(2*pi)

% Incident wave
PW = @(X) exp(1i*k*X*X0');

% Incident wave representation
plot(sphere,real(PW(sphere.vtx)))
plot(square,real(PW(square.vtx)))
title('Incident wave')
xlabel('X');   ylabel('Y');   zlabel('Z');
hold off
view(0,10)

%%% PREPARE OPERATORS
disp('~~~~~~~~~~~~~ PREPARE OPERATORS ~~~~~~~~~~~~~')

% Green kernel function --> G(x,y) = exp(ik|x-y|)/|x-y|
Gxy         = '[exp(ikr)/r]';
gradyGxy{1} = 'grady[exp(ikr)/r]1';
gradyGxy{2} = 'grady[exp(ikr)/r]2';
gradyGxy{3} = 'grady[exp(ikr)/r]3';

% Finite elements
u = fem(sphere,typ);
v = fem(sphere,typ);

% Coupling coeff
beta = 1i*k*0.5;

% Number of pool
Nlab = length(Composite);

% Domain decomposition for u
[Ilab,sigmaLab,uLab] = femSubdivide(sigma,u,Nlab,10);
drawnow

% Parallel loop for Full matrix
tic
spmd
    % Initialize composite
    M  = cell(1,numlabs);

    % Normal loop
    for j = 1:Nlab
        % Mass matrix
        Id = integral(sigmaLab{labindex},uLab{labindex},uLab{j});

        % Single layer
        Sr = 1/(4*pi) .* regularize(sigmaLab{labindex},sigmaLab{j},...
            uLab{labindex},'[1/r]',uLab{j});
        S  = 1/(4*pi) .* integral(sigmaLab{labindex},sigmaLab{j},...
            uLab{labindex},Gxy,k,uLab{j},tol);

        % Double layer
        Dr = 1/(4*pi) .* regularize(sigmaLab{labindex},sigmaLab{j},...
            uLab{labindex},'grady[1/r]',ntimes(uLab{j}));
        D  = 1/(4*pi) .* integral(sigmaLab{labindex},sigmaLab{j},...
            uLab{labindex},gradyGxy,k,ntimes(uLab{j}),tol);

        % Dirichlet Brackage-Werner :  [1i*k*beta*S - (Id/2 + D)]mu
        M{j} = beta.*(S+Sr) - (0.5*Id + (D+Dr));
    end
end
toc

% Define LHS
LHS = @(V) spmdProduct(Ilab,M,V);

% Finite element incident wave trace --> \int_Sx psi(x)' pw(x) dx
RHS = - integral(sigma,u,PW);

tic
LHS(RHS);
toc

%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% Solve linear system : [1i*k*beta*S - (Id/2 + D)] = P0
tic
mu = mgcr(LHS,RHS,[],tol,100);
toc

% Jump for derivative
lambda = beta * mu;

%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
theta = 2*pi/1e3 .* (1:1e3)';
nu    = [sin(theta),zeros(size(theta)),cos(theta)];

% Green kernel function
Ginf      = '[exp(-ikxy)]';
gradxGinf = {'gradx[exp(-ikxy)]1','gradx[exp(-ikxy)]2','gradx[exp(-ikxy)]3'};

% Finite element infinite operators
Sinf = 1/(4*pi) .* integral(nu,sigma,Ginf,k,v,tol);
Dinf = 1/(4*pi) .* integral(nu,sigma,gradxGinf,k,ntimes(v),tol);

% Finite element radiation
sol = Sinf*lambda - Dinf*mu;

% Analytical solution
ref = sphereHelmholtz('inf','dir',1,k,nu);
norm(ref-sol,2)/norm(ref,2)
norm(ref-sol,'inf')/norm(ref,'inf')

% Graphical representation
figure
plot(theta,log(abs(sol)),'b',theta,log(abs(ref)),'--r')

%%% DOMAIN SOLUTION
disp('~~~~~~~~~~~~~ RADIATION ~~~~~~~~~~~~~')

% Green kernel function --> G(x,y) = exp(ik|x-y|)/|x-y|
Gxy      = '[exp(ikr)/r]';
gradyGxy = {'grady[exp(ikr)/r]1','grady[exp(ikr)/r]2','grady[exp(ikr)/r]3'};

% Finite element mass matrix --> \int_Sx psi(x)' psi(x) dx
Id = integral(sigma,u,v);

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' G(x,y) psi(y) dx dy
tic
Sbnd = 1/(4*pi) .* (integral(sigma,sigma,u,Gxy,k,v,tol) + ...
    regularize(sigma,sigma,u,'[1/r]',v));
toc

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' dny G(x,y) psi(y) dx dy
tic
Dbnd = 1/(4*pi) .* (integral(sigma,sigma,u,gradyGxy,k,ntimes(v),tol) + ...
    regularize(sigma,sigma,u,'grady[1/r]',ntimes(v)));
toc

% Boundary solution
Psca = Id\(Sbnd*lambda - (0.5*Id*mu + Dbnd*mu));
Pinc = PW(u.dof);
Pbnd = Pinc + Psca;

% Finite element radiative operator --> \int_Sy G(x,y) psi(y) dy
tic
Sdom = 1/(4*pi) .* (integral(square.vtx,sigma,Gxy,k,v,tol) + ...
    regularize(square.vtx,sigma,'[1/r]',v));
toc

% Finite element radiative operator --> \int_Sx \int_Sy psi(x)' grady(G(x,y)) ny.psi(y) dx dy
tic
Ddom = 1/(4*pi) .* ( integral(square.vtx,sigma,gradyGxy,k,ntimes(v),tol) + ...
    regularize(square.vtx,sigma,'grady[1/r]',ntimes(v)) );
toc

% Domain solution
Psca = Sdom*lambda - Ddom*mu;
Pinc = PW(square.vtx);
Pdom = Pinc + Psca;

% Annulation sphere interieure
r             = sqrt(sum(square.vtx.^2,2));
Pdom(r<=1.01) = Pinc(r<=1.01);

% Graphical representation
figure
plot(sphere,abs(Pbnd))
axis equal;
hold on
plot(square,abs(Pdom))
title('Total field solution')
colorbar
hold off
view(0,10)

%%% ANAYTICAL SOLUTIONS FOR COMPARISONS
% Analytical solution
Pbnd = sphereHelmholtz('dom','dir',1,k,1.001*sphere.vtx) + PW(sphere.vtx);
Pdom = sphereHelmholtz('dom','dir',1,k,square.vtx) + PW(square.vtx);

% Solution representation
figure
plot(sphere,abs(Pbnd))
axis equal;
hold on
plot(square,abs(Pdom))
title('Analytical solution')
colorbar
hold off
view(0,10)

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function r = cosh(a)
%COSH         Taylor hyperbolic cosine  cosh(a)
%

% written  05/21/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                   % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = getappdata(0,'INTLAB_TAYLOR_ORDER');

  st = a.t;
  r = a;
  N = size(a.t,2);
  st(1,:) = sinh(a.t(1,:));
  r.t(1,:) = cosh(a.t(1,:));
  for j=2:K
    at_ = a.t(2:j,:);           % some 3 % faster
    st(j,:) = sum( repmat((1:j-1)',1,N).*r.t(j-1:-1:1,:).*at_ , 1 ) ./ (j-1);
    r.t(j,:) = sum( repmat((1:j-1)',1,N).*st(j-1:-1:1,:).*at_ , 1 ) ./ (j-1);
  end
  r.t(K+1,:) = sum( repmat((1:K)',1,N).*st(K:-1:1,:).*a.t(2:K+1,:) , 1 ) ./ K;

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function SNR = Convert_SNR(SNR_dB)

% Convert back from dBs to SNR value
%SNR_dB = 10 * Log10(SNR)

SNR = 10^(SNR_dB/10);

author: Sable

repo: mcbench-benchmarks


text:
function value = dmach ( job )

%*****************************************************************************80
%
%% DMACH computes machine parameters of floating point arithmetic.
%
%  Discussion:
%
%    This routine is for testing only.  It is not required by LINPACK.
%
%    If there is trouble with the automatic computation of these quantities,
%    they can be set by direct assignment statements.
%    We assume the computer has
%
%      B = base of arithmetic
%      T = number of base B digits
%      L = smallest possible exponent
%      u = largest possible exponent
%
%    then
%
%      EPS = B**(1-T)
%      TINY = 100.0D+00 *B**(-L+T)
%      HUGE = 0.01D+00 *B**(U-T)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer JOB:
%    1: requests EPSILON;
%    2: requests TINY;
%    3: requests HUGE.
%
%    Output, real VALUE, the requested value.
%
  s = 1.0;

  while ( 1 )

    tiny = s;
    s = s / 2.0;

    if ( s * 1.0 == 0.0 | s == 0.0 )
      break
    end

  end

  tiny = ( tiny / eps ) * 100.0;
  huge = 1.0 / tiny;

  if ( job == 1 )
    value = eps;
  elseif ( job == 2 )
    value = tiny;
  elseif ( job == 3 )
    value = huge;
  else
    xerbla ( 'DMACH', 1 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [F,M,Cq,Cp,QE,Qp] = spm_eeg_invert_EBoptimise(AY,UL,opttype,Qp,Qe,Qe0)
%% function [F,M,Cq,Cp,QE,Qp] = spm_eeg_invert_EBoptimise(AY,UL,opttype,Qp,Qe,Qe0)
%% Empirical Bayes optimization of priors Qp and Qe to fit data AY based on lead fields UL
% AY concatenated dimension reduced trials of M/EEG data
% UL dimension reduced lead field
% Qp source level priors- where Qp{i}.q holds an eigenmode. So source covariance
%                        component is Qp{i}.q*Qp{i}.q'.
%                         Alternately Qp{i} could be full source covariance component
% Qe sensor noise prior
% Qe0 floor of noise power to signal power (posteiror estimate of sensor noise will always be at
% least this big)
% opttype- how to optimize 'ARD','GS' or 'REML'

%% QE sensor noise posterior
%% Cp source level posterior (source by source variance)
%% F free energy
%% M MAP operator
%% Cq conditional variance
%% F free energy
%% Qp contains the posterior in same form as prior

% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes
% $Id: spm_eeg_invert_EBoptimise.m 6494 2015-07-06 10:23:04Z gareth $

if ~iscell(Qe),
    Qe={Qe};
end;

AYYA=AY*AY'; %% covariance over trials
Q0          = Qe0*trace(AYYA)*Qe{1}; %% fixed (min) level of sensor space variance
% Get source-level priors (using all subjects)
%--------------------------------------------------------------------------
j=1;
ploton=0;
if ploton,
    figure;
end;

[LQpL,Q,sumLQpL,QE,Cy,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp,Qe,ploton);

sparseind=[];
for k=1:length(Qp), %% split up priors into those with sparse (eigenmode) and full matrix form
    if isfield(Qp{k},'q'),
        sparseind=[sparseind k];
    end;
end;
fullind=setxor([1:length(Qp)],sparseind);
Qp_sp=Qp(sparseind);
Qp_full=Qp(fullind);

if isfield(opttype{j},'GSopt'),
    % Greedy search over MSPs
    %% needs to work with sparse covariance matrices Qp{i}.q
    %------------------------------------------------------------------

    if isempty(AY),
        error('NEED AY for greedy search');
    end;

    if ~isempty(sparseind),
        Qp=Qp_sp;
        [LQpL,Q,sumLQpL,QE,Cy,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp,Qe,ploton);

        MVB   = spm_mvb(AY,UL,[],Q,Qe,16);
        %%  THE VERSION BELOW IS MORE PEDESTRIAN BUT EASIER TO FOLLOW
        %  MVB_grb = spm_mvb_slow_grb( AY,UL,Q,QE,16,Q0 );

        Ne=length(Qe);
        QE=zeros(size(Qe{1}));
        for j=1:Ne,
            QE=QE+MVB.h(j)*Qe{j};
        end;

        Qcp           = Q*MVB.cp; %% mvb works with unity dipole moments so here they get scaled back up by priors

        qful=Qcp*Q';

        Qp=compact_form({qful});

        Qp=[Qp_full Qp]; %% keep full components too
        Qe={QE};
        F=max(MVB.F);

    else
        disp('Skipping GS as no eigenmode priors');
        F=-Inf;
    end;

end; %%GS

if  isfield(opttype{j},'ARDopt'),
    %% needs to work with sparse (svd decomposed) source covariance matrices
    Nn=size(AY,2); %% number of data samples used to make up covariance matrix
%    [LQpL,Q,sumLQpL,QE,Cy,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp,Qe,ploton);
    if ~isempty(sparseind),
        Qp=Qp_sp;

        [LQpL,Q,sumLQpL,QE,Cy,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp,Qe,ploton);

        ardthresh=opttype{j}.ARDopt; %%
        fprintf('ARD removing components less than 1/%3.2f max\n',ardthresh);

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

        %% SPM_SP_REML STARTS WITH EIGEN MODES Lq RATHER THAN FULL COV MATRICES
        fprintf('entering REML for ARD')

        %[Cy,h,Ph,F0] = spm_sp_reml(AYYA,[],[Qe Lq],Nn);
        [Cy,h,Ph,F0] = spm_sp_reml(AYYA,[],[Qe Lq],Nn);

        if isnan(F0),
            warning('NAN in ARD stage, dropping out');
            return;
        end;
        % Spatial priors (QP)
        %------------------------------------------------------------------
        % h provides the final weights of the hyperparameters
        Ne    = length(Qe);
        Np    = length(Qp);

        hp    = h(Ne + (1:Np));

        dropind=find(hp<max(hp)/ardthresh);
        fprintf('Removing %d of %d components\n',length(dropind),length(hp));
        hp(dropind)=0;
        keepind=find(hp>max(hp)/ardthresh);

        h=[h(1:Ne) hp(keepind)'];

        [LQpL,Q,sumLQpL,QE,Csensor,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp(keepind),Qe,ploton,h);

        [Cy,h2,Ph,F] = spm_sp_reml(AYYA,[],[Qe Lq],Nn);
        [LQpL,Q,sumLQpL,QE,Csensor,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp(keepind),Qe,ploton,h2);
        % Accumulate empirical priors (New set of patches for the second inversion)
        fprintf('ARD improvement %3.2f\n',F-F0);
        Qp=compact_form({Cp});
        Qp=[Qp_full Qp]; %% keep full components too
    else
        disp('Skipping ARD as no eigenmode priors');
        F=-Inf;
    end; % if LQ empty

end; % ARD

if  isfield(opttype{j},'REMLopt'),
    %  ReML: can work with both full and sparse source covariances
    %------------------------------------------------------------------
    Nn=size(AY,2); %% number of data samples used to make up covariance matrix

    Qp=[Qp_full Qp_sp];

    LQpL=[];

    if numel(Qp)>5,
        warning('Compressing %d priors into one full and one sparse for REML stage',length(Qp));
        count=0;
        Qp=[];
        if ~isempty(Qp_full),
            [dumLQpL,Qdum,sumLQpL,QE,Cy,M,Cp]=spm_eeg_assemble_priors(UL,Qp_full,Qe,ploton);
            count=count+1;
            LQpL{count}=sumLQpL;
            Qp{count}=Cp;
        end;
        if ~isempty(Qp_sp),
            [dumLQpL,Qdum,sumLQpL,QE,Cy,M,Cp]=spm_eeg_assemble_priors(UL,Qp_sp,Qe,ploton);
            count=count+1;
            LQpL{count}=sumLQpL;
            Qp{count}=Cp;
        end;
    else
        [LQpL,Q,sumLQpL,QE]=spm_eeg_assemble_priors(UL,Qp,Qe,ploton);
    end; % if

    %%% NOW OPTMIZE MIXTURE OF PRIOR COVARIANCE COMPS WITH REML

    [Cy,h,Ph,F] = spm_reml_sc(AYYA,[],[Qe LQpL],Nn,-4,16,Q0);

    %% Now convert the original priors to scaled versions of themselves

    [LCpL,Q,sumLCpL,QE,Cy,M,Cp,Cq]=spm_eeg_assemble_priors(UL,Qp,Qe,ploton,h);

    %% THIS NEXT LINE IS JUST A CHECK THAT THE POSTERIORS WORK AS PRIORS F2 should be greater or equal to F

    [Cy2,h2,Ph2,F2] = spm_reml_sc(AYYA,[],[{QE} {sumLCpL}],Nn,-4,16,Q0);

    if F2+1<F,
        F2-F
        error('something wrong here');
    end;

    Qp{1}=Cp;
    %Qp=compact_form({Cp});

end; %% REML opt

% re-do ReML (with informative hyperpriors)
%----------------------------------------------------------------------

function [Qp,Q]=compact_form(QP)

if isfield(QP{1},'q'),
    for j=1:length(QP),
        v=1;
        if isfield(QP{j},'v'),
            v=QP{j}.v;
        end;
        Q(:,j)=QP{j}.q*sqrt(v);
        Qp{j}.q=Q(:,j);
    end;
    disp('already sparse, returning');
    return;
end;

if length(QP)>1,
    error('only for single cells');
end;

[q,v] = spm_svd(QP{1});

if size(v,1)==size(QP{1},1),
    disp('Cannot sparsify returning full cov matrix');
    Qp=QP;
    return
end;

Q=q*sqrt(v); %% eigen vector

for i=1:size(Q,2),
    Qp{i}.q=Q(:,i);
end;
return;

author: spm


text:
function [year,month,day,hour,minute,second]=UTC2Cal(Jul1,Jul2,giveDayFrac)
%%UTC2Cal Convert universal coordinated times (UTC) given as a
%         two-part pseudo-Julian dates to a date in terms of the Gregorian
%         calendar in years, months, days, hours, minutes and seconds or
%         in terms of the Gregorian calendar in years, months, days and
%         the fraction of a day.
%
%INPUTS: Jul1, Jul2 Vectors or matrices of the two parts of Julian dates
%                   given in UTC. The units of the date are days. The full
%                   date is the sum of both terms. The date is broken into
%                   two parts to provide more bits of precision. It does
%                   not matter how the date is split.
%       giveDayFrac An optional boolean variable specifying whether the
%                   output should be as years months, days and a fraction
%                   of a day. If this parameter is omitted, the default is
%                   false. That means that the output will be given as
%                   years, months, days, hours, minutes, and seconds.
%
%OUTPUTS: Regardless of the value of giveDayFrac, the first three outputs
%         are the same. They are:
%           year A vector or matrix of years in the Gregorian calendar
%                under UTC time, one for for each Julian date.
%          month A vector or matrix of months in the Gregorian calendar
%                under UTC time, one for each Julian date. 1<=month<=12
%            day A vector or matrix of days in the Gregorian calendar
%                under UTC time, one for each Julian date. Days count from
%                1.
%          If giveDayFrac is omitted or is false, then three additional
%          outputs are
%           hour  A vector or matrix of hours in the Gregorian calendar
%                 under UTC time, one for each Julian date. 0<=hour<=23.
%          minute A vector or matrix of minutes in the Gregorian calendar
%                 under UTC time, one for each Julian date.
%                 0<=minute<=59.
%          second A vector or matrix of seconds in the Gregorian calendar
%                 under UTC time, one for each Julian date. This includes
%                 the possibility of a leap second on the day in question.
%          On the other hand, if giveDayFrac is true, then the one
%          additional output is
%          dayFrac A vector or matrix of values >=0 and <1 indicating the
%                  fraction of the day elapsed, one for each Julian date.
%
%The UTC date is only pseudo-Julian, because there is not a fixed number
%of seconds in a Julian day. The convention used in the IAU standard is
%that the Julian day matches the UTC day regardless of whether the UTC day
%is 86399, 86400 or 86401 SI seconds (depending on the presence of leap
%seconds).
%
%UTC began at 1960 January 1.0 (JD 2436934.5) and this function should not
%be called with an earlier date.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[year,month,day,hour,minute,second]=UTC2Cal(Jul1,Jul2);
%or using
%[year,month,day,dayFrac]=UTC2Cal(Jul1,Jul2,true);
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function L=sqrtpsdinv(X)

   [u,d]=schur(X);
   d=diag(d);
   d(d<=0)=inf;
   L=u*diag(1./sqrt(d))*u';
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [lstmStates, trainData, attnInfos] = rnnLayerForward(W_rnn, W_emb, prevState, input, masks, params, rnnFlags, trainData, model, charData)
% Running Multi-layer RNN for one time step.
% Input:
%   W_rnn: recurrent connections of multiple layers, e.g., W_rnn{ll}.
%   prevState: previous hidden state, e.g., for LSTM, prevState.c{ll}, prevState.h{ll}.
%   input: indices for the current batch
%   isTest: 1 -- don't store intermediate results in each state
%   isAttn: for attention, require attnData to be non-empty, has
%     attnData.srcHidVecsOrig and attnData.srcLens.
%   isDecoder: 0 -- encoder, 1 -- decoder
% Output:
%   nextState
%
% Thang Luong @ 2015, <lmthang@stanford.edu>

T = size(input, 2);
lstmStates = cell(T, 1);

% attention, encoder
attnInfos = cell(T, 1);
if rnnFlags.attn && rnnFlags.decode == 0
  assert(T == params.numSrcHidVecs);
  trainData.srcHidVecsOrig = zeroMatrix([params.lstmSize, params.curBatchSize, T], params.isGPU, params.dataType); % params.numSrcHidVecs
end

for tt=1:T % time
  if rnnFlags.charSrcRep && rnnFlags.decode == 0 % char representation, encoder
    inputEmb = zeroMatrix([params.lstmSize, params.curBatchSize], params.isGPU, params.dataType);

    % charData.rareFlags: to know which words are rare
    % charData.rareWordReps: the actual rare word representations
    % rareWordMap: to find out indices in rareWordReps
    rareIds = find(charData.rareFlags(:, tt));
    freqIds = find(~charData.rareFlags(:, tt));

    if params.assert
      assert(all(ismember(rareIds, find(masks(:, tt)==1))));
    end

    % embeddings for rare words
    if ~isempty(rareIds)
      inputEmb(:, rareIds) = charData.rareWordReps(:, charData.rareWordMap(input(rareIds, tt)));
    end

    % embeddings for frequent words
    inputEmb(:, freqIds) = W_emb(:, input(freqIds, tt));
  else
    inputEmb = W_emb(:, input(:, tt));
  end

  % multi-layer RNN
  [prevState, attnInfos{tt}] = rnnStepLayerForward(W_rnn, inputEmb, prevState, masks(:, tt), params, rnnFlags, trainData, model);

  % encoder, attention
  if rnnFlags.attn && rnnFlags.decode == 0
    trainData.srcHidVecsOrig(:, :, tt) = prevState{end}.h_t;
  end

  % store all states
  lstmStates{tt} = prevState;
end

author: lmthang


text:
function submit(part)
%SUBMIT Submit your code and output to the ml-class servers
%   SUBMIT() will connect to the ml-class server and submit your solution

  fprintf('==\n== [ml-class] Submitting Solutions | Programming Exercise %s\n==\n', ...
          homework_id());
  if ~exist('part', 'var') || isempty(part)
    partId = promptPart();
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid homework part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  [login password] = loginPrompt();
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to ml-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    % Get Challenge
    [login, ch, signature] = getChallenge(login);
    if isempty(login) || isempty(ch) || isempty(signature)
      % Some error occured, error string in first return element.
      fprintf('\n!! Error: %s\n\n', login);
      return
    end

    % Attempt Submission with Challenge
    ch_resp = challengeResponse(login, password, ch);
    [result, str] = submitSolution(login, ch_resp, partId, output(partId), ...
                                 source(partId), signature);

    fprintf('\n== [ml-class] Submitted Homework %s - Part %d - %s\n', ...
            homework_id(), partId, partNames{partId});
    fprintf('== %s\n', strtrim(str));
    if exist('OCTAVE_VERSION')
      fflush(stdout);
    end
  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '1';
end

function [partNames] = validParts()
  partNames = { 'Warm up exercise ', ...
                'Computing Cost (for one variable)', ...
                'Gradient Descent (for one variable)', ...
                'Feature Normalization', ...
                'Computing Cost (for multiple variables)', ...
                'Gradient Descent (for multiple variables)', ...
                'Normal Equations'};
end

function srcs = sources()
  % Separated by part
  srcs = { { 'warmUpExercise.m' }, ...
           { 'computeCost.m' }, ...
           { 'gradientDescent.m' }, ...
           { 'featureNormalize.m' }, ...
           { 'computeCostMulti.m' }, ...
           { 'gradientDescentMulti.m' }, ...
           { 'normalEqn.m' }, ...
         };
end

function out = output(partId)
  % Random Test Cases
  X1 = [ones(20,1) (exp(1) + exp(2) * (0.1:0.1:2))'];
  Y1 = X1(:,2) + sin(X1(:,1)) + cos(X1(:,2));
  X2 = [X1 X1(:,2).^0.5 X1(:,2).^0.25];
  Y2 = Y1.^0.5 + Y1;
  if partId == 1
    out = sprintf('%0.5f ', warmUpExercise());
  elseif partId == 2
    out = sprintf('%0.5f ', computeCost(X1, Y1, [0.5 -0.5]'));
  elseif partId == 3
    out = sprintf('%0.5f ', gradientDescent(X1, Y1, [0.5 -0.5]', 0.01, 10));
  elseif partId == 4
    out = sprintf('%0.5f ', featureNormalize(X2(:,2:4)));
  elseif partId == 5
    out = sprintf('%0.5f ', computeCostMulti(X2, Y2, [0.1 0.2 0.3 0.4]'));
  elseif partId == 6
    out = sprintf('%0.5f ', gradientDescentMulti(X2, Y2, [-0.1 -0.2 -0.3 -0.4]', 0.01, 10));
  elseif partId == 7
    out = sprintf('%0.5f ', normalEqn(X2, Y2));
  end
end

function url = challenge_url()
  url = 'http://www.ml-class.org/course/homework/challenge';
end

function url = submit_url()
  url = 'http://www.ml-class.org/course/homework/submit';
end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          while ~feof(fid)
            line = fgets(fid);
            src = [src line];
          end
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature] = getChallenge(email)
  str = urlread(challenge_url(), 'post', {'email_address', email});

  str = strtrim(str);
  [email, str] = strtok (str, '|');
  [ch, str] = strtok (str, '|');
  [signature, str] = strtok (str, '|');
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'homework', homework_id(), ...
            'part', num2str(part), ...
            'email', email, ...
            'output', output, ...
            'source', source, ...
            'challenge_response', ch_resp, ...
            'signature', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Password: ', 's');
end

function [str] = challengeResponse(email, passwd, challenge)
  salt = ')~/|]QMB3[!W`?OVt7qC"@+}';
  str = sha1([challenge sha1([salt email passwd])]);
  sel = randperm(numel(str));
  sel = sort(sel(1:16));
  str = str(sel);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function [ x, know ] = p06_sol ( m, know )

%*****************************************************************************80
%
%% P06_SOL returns known solutions for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = [ ...
      0.509282516910744, ...
      0.509282516910744, ...
      0.509282516910746, ...
      0.509282516910744 ]';
  else
    know = 0;
    x = zeros ( m, 1 );
  end

  return
end


text:
function S = Polygons_intersection_Polygon_cleanup(S,accuracy)
if nargin == 1
    accuracy = 1e-9;
end
if accuracy<eps
    error('Polygons_intersection_Polygon_cleanup:e0',...
          'Accuracy must be a positive real');
end

if isfield(S,'P')
    S_area = struct('A', {});
    for i=1:numel(S)
        S(i).area = 0;
        for j=1:numel(S(i).P)
            S_area(i).A(j) = polyarea(S(i).P(j).x,S(i).P(j).y);
            S(i).area      = S(i).area + (1-2*S(i).P(j).hole) * S_area(i).A(j);
        end
        To_delete = S_area(i).A < accuracy;
        S(i).P(To_delete) = [];
        S_area(i).A(To_delete) = [];
    end
elseif isfield(S,'x')
    S_area = zeros(numel(S),1);
    for i=1:numel(S)
        S_area(i) = polyarea(S(i).x,S(i).y);
    end
    To_delete = S_area < accuracy;
    S(To_delete) = [];
    S_area(To_delete) = [];
end


text:
function appdata = sliceomaticdnorm(d,norm)
% SLICEOMATICDNORM(rawdata,norm) - Normalize relevant data.

% Simplify the isonormals
d.reduce=d.reduce./norm.*100;
d.reducesmooth=d.reducesmooth./norm.*100;
d.smooth=d.smooth./norm.*100;
d.data=d.data./norm.*100;
appdata = d;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function sig = calcDoseSignature(doseIndex, planC)
%"calcDoseSignature"
%   Returns a signature to identify a dose distribution, for the purpose of
%   flagging DVHs as stale once the stored signature does not match an
%   existing dose distribution.
%
%   Currently uses the dose sum and the linear position of the
%   maxdose.  Some CRC signature may be possible in the future, if fast
%   code is available.
%
%JRA 12/3/04
%
%Usage:
%    function dS = calcDoseSignature(doseIndex, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

try
    dA = getDoseArray((doseIndex), planC);

    [maxDose, loc]  = max(dA(:));
    doseSum         = sum(dA(:));
    sig             = [doseSum, loc];
catch
    error('Could not calculate dose signature for requested dose distribution.');
end


text:
% Compare basic MatConvNet and Caffe blocks numerically
rootpath = fileparts(fileparts(mfilename('fullpath')));
run(fullfile(rootpath, 'matlab', 'vl_setupnn.m'));

inputScale = 100;

caffemodel = fullfile('data', 'tmp_caffemodels', 'test_model');
[~,~,~] = mkdir('data');
[~,~,~] = mkdir(fileparts(caffemodel));

%%
layers = {};
layers{end+1} = struct(...
  'name', 'conv', ...
  'type', 'conv', ...
  'stride', [2, 2], ...
  'pad', [1, 1, 1, 1], ...
  'weights', {{rand(3, 3, 10, 5, 'single'), rand(5, 1, 'single')}});

layers{end+1} = struct(...
  'name', 'relu', ...
  'type', 'relu');

layers{end+1} = struct(...
  'name', 'norm', ...
  'type', 'normalize', ...
  'param', [5, 1, 2e-5, 0.75]);

layers{end+1} = struct(...
  'name', 'softmax', ...
  'type', 'softmax');

%%
net_ = struct();
net_.meta.normalization.imageSize = [20, 20, 10];
net_.meta.normalization.averageImage = rand(1, 1, 10);

diffStats = zeros(numel(layers), 3);
for li = 1:numel(layers)
  net_.layers = layers(li);
  layerName = layers{li}.name;
  simplenn_caffe_deploy(net_, caffemodel, 'doTest', false, ...
    'outputBlobName', layerName, 'silent', true);
  res = simplenn_caffe_compare(net_, caffemodel, [], ...
    'randScale', inputScale, 'silent', true);
  diffStats(li, :) = res.(layerName);
end

fprintf('Results: \n');
layerNames = cellfun(@(l) l.name, layers, 'UniformOutput', false);
fprintf('Layer   %s\n', sprintf('% 10s', layerNames{:}));
fprintf('MeanErr %s\n', sprintf('% 10.2e', diffStats(:, 2)));
fprintf('MaxErr  %s\n', sprintf('% 10.2e', diffStats(:, 3)));

rmdir(fileparts(caffemodel), 's');

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [ m, d ] = easter_egr2 ( y )

%*****************************************************************************80
%
%% EASTER_EGR2 computes the month and day of Easter for a Common year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm P,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, page 376.
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, integer M, D, the month and day of Easter.
%
  if ( y <= 0 )
    m = -1;
    d = -1;
    return
  end

  a = floor ( y / 100 );
  b = a - floor ( a / 4 );
  c = mod ( y, 19 );
  d = mod ( 15 + 19 * c + b - floor ( ( a - floor ( ( a - 17 ) / 25 ) ) / 3 ), 30 );
  e = d - floor ( ( c + 11 * d ) / 319 );
  s = 22 + e + mod ( 140004 - y - floor ( y / 4 ) + b - e, 7 );

  m = 3 + floor ( s / 32 );
  d = i4_wrap ( s, 1, 31 );

  return
end


text:
function h = text(m,varargin)
% plot Miller indece
%
% Input
%  m  - Miller
%
% Options
%  symmetrised - plot symmetrically equivalent directions
%  antipodal   - include antipodal symmetry
%  labeled     - plot Miller indice as label
%  label       - plot user label
%
% See also
% vector3d/text

% symmetrise
if check_option(varargin,'symmetrised') && ~check_option(varargin,'skipSymmetrise')

  [m,l] = symmetrise(m,varargin{:},'unique','noAntipodal');

  % symmetrise labels
  if ~check_option(varargin,'labeled') && ~isempty(varargin)
    strings = ensurecell(varargin{1});
    if iscellstr(varargin{1}) && ~isempty(strings)

      if numel(strings)==1
        strings = repcell(strings{1},length(m),1);
      else
        strings = strings(repelem(1:numel(strings),l));
      end
      varargin{1} = strings;
    end
  end

  varargin = [varargin,{'skipSymmetrise','noAntipodal'}];
end

% ensure specific plot options
varargin = [varargin(1),m.CS.plotOptions,varargin(2:end)];

h = text@vector3d(m,varargin{:});

if nargout == 0, clear h; end

author: mtex-toolbox

repo: mtex

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


text:
function varargout = eig(f)
%EIG   Eigenvalues and eigenfunctions of a CHEBFUN2.
%   EIG(F) returns the eigenvalues of F. The number of nonzero eigenvalues
%   returned is at most the rank of the CHEBFUN2. F.domain needs to be
%   square: [a b a b], just as matrices need to be square for eigenvalues
%   to be defined.
%
%   S = EIG(F) returns the eigenvalues of F. S is a vector of eigenvalues.
%
%   [V, D] = EIG(F) returns the eigenvalues and eigenfunctions of F. V is a
%   quasimatrix of CHEBFUN objects and D is a diagonal matrix with the
%   eigenvalues on the diagonal.
%
%   If d is an eigenvalue of F with eigenfunction w, then this holds:
%   \int f(x,y) w(x) dx = d w(y).
%   This can be checked as follows:
%   [V, D] = eig(f);
%   and check that norm(f*V-V*D) is O(eps).

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

dom = f.domain;
if ( (dom(1) ~= dom(3)) || (dom(2) ~= dom(4)) ) % domain check
    error('CHEBFUN:CHEBFUN2:eig:domainerr',...
          'Domain of chebfun2 needs to be in form [a b a b]');
end

[u, S, v] = svd(f); % SVD of chebfun2
[V, D] = eig(v' * u * S); % eig(AB) = eig(BA)
V = u * S * V; % eigenfunctions

if ( nargout > 1 )
    for ii = 1:size(V,2)
    V(:,ii) = V(:,ii) / norm(V(:,ii)); % normalize eigenfunctions
    end
    varargout = { V, D };
else
    varargout = { diag(D) };
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function J = amedfilt2_calc(I) %#codegen
% 2-D Adaptive Median Filter
% This filter ignores edge effects and boundary conditions, as such, the
% output is a cropped version of the original image, where the amount
% cropped is equal to the maximum window size vertically and horizontally.

% Define smax as a constant
smax = 9;

% Initialize Output Image (J)
J = I;

% Calculate valid region limits for filter
[nrows ncols] = size(I);
ll = ceil(smax/2);
ul = floor(smax/2);

% Loop over the entire image ignoring edge effects
for rows = ll:nrows-ul
    for cols = ll:ncols-ul

        window_ind = -ul:ul;
        region = I(rows+window_ind,cols+window_ind);
        centerpixel = region(ll,ll);

        for s = 3:2:smax

            % We can collapse the ROI calculations into a single function
            [rmin,rmax,rmed] = roi_stats(region,smax,s);

            % adapt region size
            if rmed > rmin && rmed < rmax
                if centerpixel <= rmin || centerpixel >= rmax
                    J(rows,cols) = rmed;
                end

                % stop adapting
                break;
            end
        end
    end
end

function [rmin,rmax,rmed] = roi_stats(region,smax,s)
% Limits for ROI
ll = ceil(smax/2)-floor(s/2);
ul = ceil(smax/2)+floor(s/2);

v = ones(smax*smax,1);
count = 1;

for i = ll:ul
    for j = ll:ul
        v(count) = region(i,j);
        count = count+1;
    end
end

v = visort(v,s*s);
rmed = v(ceil(s*s/2));
rmin = v(1);
rmax = v(s*s);


text:
function [s,f] = ltas(x,fs,varargin)
%LTAS calculate the long-term average spectrum of a signal
%
%   S = IOSR.DSP.LTAS(X,FS) calculates the long-term average spectrum
%   (LTAS) of signal X, sampled at FS Hz. The spectrum is calculated from
%   the average power spectral density (PSD) obtained from a series of
%   overlapping FFTs; the FFT length is 4096, and the hop size is 2048. The
%   segments of X are Hann-windowed. The average PSD is then
%   Gaussian-smoothed to 1/3-octave resolution.
%
%   X can be a vector, matrix, or multidimensional array; LTAS will operate
%   along the first non-signleton dimension, and return the LTAS for each
%   corresponding row/column/etc.
%
%   S = IOSR.DSP.LTAS(X,FS,'PARAMETER','VALUE') allows numerous parameters
%   to be specified. These parameters are:-
%       'dim'     : {find(size(X)>1,1,'first')} | scalar
%           Specifies the dimension of operation (defaults to the first
%           non-singleton dimension).
%       'graph'   : {false} | true
%           Choose whether to plot a graph of the LTAS.
%       'hop'     : {NFFT/2} | scalar
%           Specifies the step size through X used to calculate each
%           segment. NFFT is determined by the 'win' parameter.
%       'noct'    : {3} | scalar
%           Apply 1/noct-octave smoothing to the frequency spectrum.
%           Setting 'noct' to 0 results in no smoothing.
%       'scaling' : {'none'} | 'max0'
%           Specifies any scaling to apply to S. By default, no scaling is
%           applied. If scaling is set to 'max0', S will be scaled to have
%           a maximum value of 0dB.
%       'units'   : {dB} | 'none'
%           Specifies the output units. By default the PSD is calculated in
%           dB. Otherwise the PSD is returned directly.
%       'win'     : {4096} | scalar | vector
%           Specifies the window or FFT length NFFT used to calculate the
%           spectrum. If 'win' is a scalar, it specifies the FFT length,
%           and a Hann window is applied to each segment. If 'win' is a
%           vector, NFFT is the length of the vector, and the vector is
%           multiplied with each segment.
%
%   [S,F] = IOSR.DSP.LTAS(...) returns the frequencies F for each
%   corresponding bin of S.
%
%   Example
%
%       % Plot the 1/6th-octave-smoothed LTAS of the Handel example
%       load handel.mat
%       figure
%       iosr.dsp.ltas(y,Fs,'noct',6,'graph',true);
%
%   See also IOSR.DSP.STFT, IOSR.DSP.SMOOTHSPECTRUM.

%   Copyright 2016 University of Surrey.

    %% parse input

    if nargin < 2
        error('iosr:ltas:nargin''Not enough input arguments')
    end

    options = struct(...
        'win',4096,...
        'hop',[],...
        'noct',3,...
        'dim',[],...
        'units','db',...
        'scaling','none',...
        'graph',false);

    % read parameter/value inputs
    if nargin>2 % if parameters are specified
        % read the acceptable names
        optionNames = fieldnames(options);
        % count arguments
        nArgs = length(varargin);
        if round(nArgs/2)~=nArgs/2
           error('iosr:ltas:nameValuePair','LTAS needs propertyName/propertyValue pairs')
        end
        % overwrite defults
        for pair = reshape(varargin,2,[]) % pair is {propName;propValue}
           IX = strcmpi(pair{1},optionNames); % find match parameter names
           if any(IX)
              % do the overwrite
              options.(optionNames{IX}) = pair{2};
           else
              error('iosr:ltas:unknownOption','%s is not a recognized parameter name',pair{1})
           end
        end
    end

    %% check and assign input

    % required inputs
    assert(isnumeric(x), 'iosr:ltas:invalidX', 'X must be numeric.');
    assert(isscalar(fs), 'iosr:ltas:invalidFs', 'FS must be a scalar.');
    assert(fs>0, 'iosr:ltas:invalidFs', 'FS must be greater than 0.');
    assert(isint(fs), 'iosr:ltas:invalidFs', 'FS must be an integer.');

    % determine fft parameters
    if numel(options.win)>1 && isvector(options.win)
        NFFT = length(options.win);
        win = options.win;
    elseif isscalar(options.win)
        NFFT = options.win;
        win = NFFT;
    else
        error('iosr:ltas:invalidWin','''WIN'' must be a vector or a scalar');
    end
    assert(isint(NFFT) && NFFT>0, 'iosr:ltas:invalidNfft', '''WIN'' must be a positive integer');

    % determine hop
    hop = options.hop;
    if isempty(hop)
        hop = fix(NFFT/2);
    end

    % smoothing
    Noct = options.noct;

    % dimension to operate along
    dim = options.dim;
    dims = size(x);
    if isempty(dim)
        dim = find(dims>1,1,'first');
    else
        assert(isnumeric(dim),'iosr:ltas:invalidDim', 'DIM must be an integer');
        assert(isint(dim), 'iosr:ltas:invalidDim', 'DIM must be an integer or empty');
        assert(dim>0, 'iosr:ltas:invalidDim', 'DIM must be greater than 0')
    end

    %% permute and rehape x to operate down columns

    % number of useful coefficients
    if mod(NFFT,2)==0
        Nout = (NFFT/2)+1;
    else
        Nout = (NFFT+1)/2;
    end

    % reorder and permute
    order = mod(dim-1:dim+length(dims)-2,length(dims))+1;
    dims_shift = dims(order);
    x = rearrange(x,order,[dims_shift(1) prod(dims_shift(2:end))]);
    dims_out_shift = dims_shift;
    dims_out_shift(1) = Nout;

    %% calculate spectra

    % choose units function
    switch lower(options.units)
        case 'db'
            units = @(x) 10*log10(x);
            max0scale = @(s) s-max(s(:));
            labelY = 'Power spectral density [dBFS]';
            yscale = 'linear';
        case 'none'
            units = @(x) x;
            max0scale = @(s) s./max(s(:));
            labelY = 'Power spectral density';
            yscale = 'log';
        otherwise
            error('iosr:ltas:unknownUnits','Unknown units option ''%s''',options.units);
    end

    % do calculations
    s = zeros(dims_out_shift);
    for c = 1:size(x,2)
        [S,f] = iosr.dsp.stft(x(:,c),win,hop,fs); % short-time ft
        s(:,c) = mean(abs(S/NFFT).^2,2); % mean PSD
        s(:,c) = units(s(:,c)); % put into dB
        s(:,c) = iosr.dsp.smoothSpectrum(s(:,c),f,Noct); % smooth
    end

    % invert permutation
    s = irearrange(s,order,dims_out_shift);

    % scale output
    switch lower(options.scaling)
        case 'max0'
            s = max0scale(s);
        case 'none'
            % do nothing
        otherwise
            error('iosr:matchEQ:unknownScaling','Unknown scaling option ''%s''',options.scaling);
    end

    %% plot

    assert(islogical(options.graph) && numel(options.graph)==1, 'iosr:ltas:invalidGraph', '''graph'' option must be logical.')
    if options.graph
        semilogx(f,rearrange(s,order,[dims_out_shift(1) prod(dims_out_shift(2:end))]));
        xlabel('Frequency [Hz]');
        ylabel(labelY);
        set(gca,'yscale',yscale);
        grid on
        if prod(dims_out_shift(2:end)) > 1
            legend(num2str((1:prod(dims_out_shift(2:end)))'));
        end
    end

end

function y = rearrange(x,order,shape)
%REARRANGE reshape and permute to make target dim column
    y = permute(x,order);
    y = reshape(y,shape);
end

function y = irearrange(x,order,shape)
%IREARRANGE reshape and permute to original size
    y = reshape(x,shape);
    y = ipermute(y,order);
end

function y = isint(x)
%ISINT check if input is whole number
    y = x==round(x);
end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function [DCM,BMR,BMA] = spm_dcm_bmr_all(DCM,field,OPT)
% Bayesian model reduction of all permutations of model parameters
% FORMAT [RCM,BMR,BMA] = spm_dcm_bmr_all(DCM,field,OPT)
%
% DCM      - A single estimated DCM (or PEB) structure:
%
%  DCM.M.pE  - prior expectation
%  DCM.M.pC  - prior covariance
%  DCM.Ep    - posterior expectation
%  DCM.Cp    - posterior covariances
%  DCM.beta  - prior expectation of reduced parameters (default: 0)
%  DCM.gamma - prior variance    of reduced parameters (default: 0)
%              NB: beta = 'pE' uses full priors
%
% field      - parameter fields in DCM{i}.Ep to optimise [default: {'A','B'}]
%             'All' will invoke all fields (i.e. random effects)
%             If Ep is not a structure, all parameters will be considered
%
% OPT        - Bayesian model selection or averaging: 'BMS' or 'BMA'
%              [default: 'BMA']
%
% Returns:
%
% DCM - Bayesian Model Average (BMA) over models in the final iteration of
%       the search:
%
%       DCM.Ep    - (BMA) posterior expectation
%       DCM.Cp    - (BMA) posterior covariance
%       DCM.Pp    - Model posterior over parameters (with and without)
%
% BMR -  (Nsub) summary structure reporting the model space from the last
%        iteration of the search:
%
%        BMR.name - character/cell array of parameter names
%        BMR.F    - free energies (relative to full model)
%        BMR.P    - and posterior (model) probabilities
%        BMR.K    - [models x parameters] model space (1 = off, 0 = on)
%
% BMA - Baysian model average (over reduced models; see spm_dcm_bma)
%
%--------------------------------------------------------------------------
% This routine searches over reduced (nested) models of a full model (DCM)
% using Bayesian model reduction and performs Bayesian Model Averaging.
% 'Reduced' means some free parameters (parameters with a non-
% zero prior covariance) are switched off by fixing their prior variance
% to zero.
%
% If there are fewer than nmax = 8 free parameters, all permutations of
% switching off parameters will be tested. Otherwise, this routine
% implements the following greedy search procedure. The nmax parameters
% are identified which, when switched off individually, produce the least
% reduction (greatest increase) in model evidence. All permutations of
% switching off these parameters are then evaluated and the best
% permutation is retained. This procedure is repeated until all nmax
% parameters are retained or there are no more parameters to consider.
% Finally, BMA is performed on the models from the last iteration.
%
% NB: The full model should be estimated prior to running this function.
%
% See also: spm_dcm_post_hoc - this routine is essentially a simplified
% version of spm_dcm_post_hoc
%__________________________________________________________________________
% Copyright (C) 2010-2014 Wellcome Trust Centre for Neuroimaging

% Karl Friston, Peter Zeidman
% $Id: spm_dcm_bmr_all.m 7717 2019-11-27 11:10:36Z peter $

%-specification of null prior covariance
%--------------------------------------------------------------------------
if isfield(DCM,'beta'),  beta  = DCM.beta;  else, beta  = 0; end
if isfield(DCM,'gamma'), gamma = DCM.gamma; else, gamma = 0; end

%-Check fields of parameter stucture (and options)
%--------------------------------------------------------------------------
if nargin < 3
    OPT   = 'BMA';
end
if nargin < 2 || isempty(field)
    field = {'A','B'};
end
if ischar(field)
    field = {field};
end

%-deal with filenames structure
%--------------------------------------------------------------------------
if ischar(DCM)
    DCM = load(DCM,'DCM');
    DCM = DCM.DCM;
end

% Get prior covariances
%--------------------------------------------------------------------------
if isstruct(DCM.M.pC), DCM.M.pC = diag(spm_vec(DCM.M.pC)); end
if spm_length(DCM.M.pE) ~= size(DCM.M.pC,1)
    DCM.M.pC = diag(spm_vec(DCM.M.pC));
end

% Get priors and posteriors
%--------------------------------------------------------------------------
qE  = DCM.Ep;
qC  = DCM.Cp;
pE  = DCM.M.pE;
pC  = DCM.M.pC;

% Remove (a priori) null space
%--------------------------------------------------------------------------
U   = spm_svd(pC);
qE  = U'*spm_vec(qE);
pE  = U'*spm_vec(pE);
qC  = U'*qC*U;
pC  = U'*pC*U;

%-Greedy search (GS) - eliminating parameters in a top down fashion
%==========================================================================

% Accumulated reduction vector (C)
%--------------------------------------------------------------------------
q   = diag(DCM.M.pC);
if sum(q < 1024)
    C   = double(q > mean(q(q < 1024))/1024);
else
    C   = double(q > 0);
end
GS  = 1;
while GS

    %-Find free coupling parameters
    %----------------------------------------------------------------------
    if isstruct(DCM.Ep)
        k = spm_fieldindices(DCM.Ep,field{:});
    else
        k = 1:spm_length(DCM.Ep);
    end
    k = k(find(C(k))); %#ok<FNDSB>

    % If there are too many parameters find those with the least evidence
    %----------------------------------------------------------------------
    nparam = length(k);
    nmax   = fix(max(nparam/4,8));
    if nparam > nmax

        % Model search over new prior without the i-th parameter
        %------------------------------------------------------------------
        Z     = zeros(1,nparam);
        for i = 1:nparam

            % Identify parameters to retain r and to remove s
            %--------------------------------------------------------------
            r   = C; r(k(i)) = 0; s = 1 - r;

            % Create reduced prior covariance matrix
            %--------------------------------------------------------------
            R   = U'*diag(r + s*gamma)*U;
            rC  = R*pC*R;

            % Create reduced prior means
            %--------------------------------------------------------------
            if isnumeric(beta)
                S  = U'*diag(r)*U;
                rE = S*pE + U'*s*beta;
            else
                rE = pE;
            end

            Z(i) = spm_log_evidence(qE,qC,pE,pC,rE,rC);
        end

        % Find parameters with the least evidence
        %------------------------------------------------------------------
        [z,i] = sort(-Z);
        k     = k(i(1:nmax));

        % Flag a greedy search
        %------------------------------------------------------------------
        GS = 1;

    elseif isempty(k)
        fprintf('\nThere are no free parameters in this model.\n')
        return
    else
        GS = 0;
    end

    % compare models
    %======================================================================
    for j = 1:2

        if j == 1
            % compare models with and without nmax parameters first
            %--------------------------------------------------------------
            K = repmat(logical([1;0]),1,numel(k));
        else
            % compare all combinations
            %--------------------------------------------------------------
            k = k(1:min(8,end));
            K = spm_perm_mtx(numel(k));
        end

        % Model search over new prior (covariance)
        %------------------------------------------------------------------
        nK    = size(K,1);
        G     = zeros(1,nK);
        for i = 1:nK

            % Identify parameters to retain (r) and to remove (s)
            %--------------------------------------------------------------
            r    = C; r(k(K(i,:))) = 0; s = 1 - r;

            % Create reduced prior covariance matrix
            %--------------------------------------------------------------
            R    = U'*diag(r + s*gamma)*U;
            rC   = R*pC*R;

            % Create reduced prior means
            %--------------------------------------------------------------
            if isnumeric(beta)
                S  = U'*diag(r)*U;
                rE = S*pE + U'*s*beta;
            else
                rE = pE;
            end

            G(i) = spm_log_evidence(qE,qC,pE,pC,rE,rC);
        end

        % if sufficient complexity reduction then omit combinations
        %------------------------------------------------------------------
        if G(1) - G(end) > nmax && nparam > nmax
            break;
        else
            nmax = 8;
        end
    end

    % posterior probability
    %----------------------------------------------------------------------
    p            = spm_softmax(G(:));

    %-Get selected model and prune redundant parameters
    %======================================================================
    [z,i]        = max(p);
    C(k(K(i,:))) = 0;

    % Continue greedy search if any parameters have been eliminated
    %----------------------------------------------------------------------
    nelim  = full(sum(K(i,:)));
    GS     = GS & nelim;

    % Show results
    % ---------------------------------------------------------------------
    fprintf('%i out of %i free parameters removed \n',nelim,nparam)

    if nmax <= 8
        spm_figure('Getwin','BMR - all'); clf
        subplot(3,2,1)
        if numel(G) > 32, plot(G,'k'), else, bar(G,'c'), end
        title('log-posterior','FontSize',16)
        xlabel('model','FontSize',12)
        ylabel('log-probability','FontSize',12)
        axis square

        subplot(3,2,2)
        if numel(G) > 32, plot(p,'k'), else, bar(p,'r'), end
        title('model posterior','FontSize',16)
        xlabel('model','FontSize',12)
        ylabel('probability','FontSize',12)
        axis square
        drawnow
    end

end

%-Inference over families (one family per coupling parameter)
%==========================================================================
for i = 1:length(k)
    Pk(1,i) = mean(p(~K(:,i)));
    Pk(2,i) = mean(p( K(:,i)));
end
Pk    = Pk(1,:)./sum(Pk);
Pp    = C;
Pp(k) = Pk;

%-Bayesian model selection or average
%==========================================================================
qE    = DCM.Ep;
qC    = DCM.Cp;
pE    = DCM.M.pE;
pC    = DCM.M.pC;
pE    = spm_vec(pE);

switch OPT

    case('BMA')
        % Bayesian model averaging
        %------------------------------------------------------------------
        Gmax     = max(G);

    case('BMS')
        % Bayesian model selection (place winning G outside Occam's window
        %------------------------------------------------------------------
        [Gmax,i] = max(G);
        G(i)     = G(i) + 16;
        Gmax     = Gmax + 16;

    otherwise
end

BMA   = {};
for i = 1:length(K)

    % if this mdel is in Occam's window, inlcude in BMA
    %----------------------------------------------------------------------
    if G(i) > (Gmax - 8)

        % reduced model
        %------------------------------------------------------------------
        r            = C;
        r(k(K(i,:))) = 0;
        s            = 1 - r;
        R            = diag(r + s*gamma);
        rC           = R*pC*R;
        S            = diag(r);
        if isnumeric(beta)
            rE       = S*spm_vec(pE) + s*beta;
        else
            rE       = pE;
        end

        % BMR
        %------------------------------------------------------------------
        [F,Ep,Cp]    = spm_log_evidence_reduce(qE,qC,pE,pC,rE,rC);
        BMA{end + 1} = struct('Ep',Ep,'Cp',Cp,'F',F);
    end
end

switch OPT
    case('BMA')

        % Bayesian model averaging
        %------------------------------------------------------------------
        BMA   = spm_dcm_bma(BMA);
        Ep    = BMA.Ep;
        Cp    = BMA.Cp;

    case('BMS')

        % Bayesian model selection
        %------------------------------------------------------------------
        BMA   = BMA{1};
        Ep    = BMA.Ep;
        Cp    = BMA.Cp;

    otherwise
end

if isstruct(Cp) || (spm_length(Cp) == spm_length(Ep))
    Cp = diag(spm_vec(Cp));
end

% Show full and reduced conditional estimates (for Bayesian average)
%--------------------------------------------------------------------------
spm_figure('Getwin','BMR - all');

if isstruct(DCM.Ep)
    i  = spm_find_pC(pC,DCM.Ep,field);
else
    i  = 1:spm_length(DCM.Ep);
end
qE     = spm_vec(qE);
Ep     = spm_vec(Ep);

j = i(ismember(i,1:length(spm_vec(Ep))));

% BMR summary and plotting
%--------------------------------------------------------------------------
try
    Pnames     = spm_fieldindices(DCM.Ep,k);
catch
    try
        Np     = numel(DCM.Pnames);
        Pnames = DCM.Pnames(rem(k - 1,Np) + 1);
    catch
        Pnames = 'all parameters';
    end
end
BMR.name = Pnames;
BMR.F    = G;
BMR.P    = p;
BMR.K    = K;
BMR.k    = k;

subplot(3,2,3), spm_plot_ci(qE(i),qC(i,i))
title('MAP (full)','FontSize',16)
axis square, a = axis;

subplot(3,2,4), spm_plot_ci(Ep(j),abs(Cp(j,j)))
title('MAP (reduced)','FontSize',16), axis square, axis(a)

subplot(3,2,5), imagesc(1 - K')
xlabel('model'), ylabel('parameter'), title('model space','FontSize',16)
set(gca,'YTickLabel',BMR.name);
axis tight, axis square

subplot(3,2,6)
Np = length(i);
bar(1:Np,diag(Pp(i)))
xlabel('parameter'), title(' posterior','FontSize',16)
axis square, drawnow, axis([0 (Np + 1) 0 1])

%-Save Bayesian parameter average (Ep,Cp) and family-wise inference (Pp)
%==========================================================================
if isstruct(DCM.Ep)
    if Np < 32
        legend(spm_fieldindices(DCM.Ep,i))
    end
    Pp    = spm_unvec(Pp,DCM.Ep);
    Ep    = spm_unvec(Ep,DCM.Ep);
end

DCM.Pp    = Pp;        % Model posterior over parameters (with and without)
DCM.Ep    = Ep;        % Bayesian model averages
DCM.Cp    = Cp;        % Bayesian model variance

% Clear free energy if supplied (which is no longer meaningful)
if isfield(DCM,'F')
    DCM = rmfield(DCM,'F');
end


text:
% test for 2D levy flight generation

save_image = 0;
n = 1024*4;
% exponent
alpha = 1;
% variance
sigma = 1;
% mean/median
beta = 0;
delta = 0;

type = 'axis';
type = 'isotropic';

alpha_list = linspace(0.6, 2, 9);
m = length(alpha_list);

clf;
for i=1:m
    alpha = alpha_list(i);
    x = gen_levy_flight(n,alpha,sigma,beta,delta,type);
    subplot(sqrt(m),sqrt(m), i);
    plot_curve(x);
    axis off;
    title( sprintf('\\alpha=%.1f', alpha) );
end

if save_image
    saveas(gcf, 'levy_flights.png');
end

author: gpeyre

repo: matlab-toolboxes


text:
function oev = eci2orb1 (mu, r, v)

% convert eci state vector to six classical orbital
% elements via equinoctial elements

% input

%  mu = central body gravitational constant (km**3/sec**2)
%  r  = eci position vector (kilometers)
%  v  = eci velocity vector (kilometers/second)

% output

%  oev(1) = semimajor axis (kilometers)
%  oev(2) = orbital eccentricity (non-dimensional)
%           (0 <= eccentricity < 1)
%  oev(3) = orbital inclination (radians)
%           (0 <= inclination <= pi)
%  oev(4) = argument of perigee (radians)
%           (0 <= argument of perigee <= 2 pi)
%  oev(5) = right ascension of ascending node (radians)
%           (0 <= raan <= 2 pi)
%  oev(6) = true anomaly (radians)
%           (0 <= true anomaly <= 2 pi)

% Orbital Mechanics with MATLAB

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

pi2 = 2.0 * pi;

% position and velocity magnitude

rmag = norm(r);

vmag = norm(v);

% position unit vector

rhat = r / rmag;

% angular momentum vectors

hv = cross(r, v);

hhat = hv / norm(hv);

% eccentricity vector

vtmp = v / mu;

ecc = cross(vtmp, hv);

ecc = ecc - rhat;

% semimajor axis

sma = 1.0 / (2.0 / rmag - vmag * vmag / mu);

p = hhat(1) / (1.0 + hhat(3));

q = -hhat(2) / (1.0 + hhat(3));

const1 = 1.0 / (1.0 + p * p + q * q);

fhat(1) = const1 * (1.0 - p * p + q * q);
fhat(2) = const1 * 2.0 * p * q;
fhat(3) = -const1 * 2.0 * p;

ghat(1) = const1 * 2.0 * p * q;
ghat(2) = const1 * (1.0 + p * p - q * q);
ghat(3) = const1 * 2.0 * q;

h = dot(ecc, ghat);

xk = dot(ecc, fhat);

x1 = dot(r, fhat);

y1 = dot(r, ghat);

% orbital eccentricity

eccm = sqrt(h * h + xk * xk);

% orbital inclination

inc = 2.0 * atan(sqrt(p * p + q * q));

% true longitude

xlambdat = atan3(y1, x1);

% check for equatorial orbit

if (inc > 0.00000001)
    raan = atan3(p, q);
else
    raan = 0.0;
end

% check for circular orbit

if (eccm > 0.00000001)
    argper = mod(atan3(h, xk) - raan, pi2);
else
    argper = 0.0;
end

% true anomaly

tanom = mod(xlambdat - raan - argper, pi2);

% load orbital element vector

oev(1) = sma;
oev(2) = eccm;
oev(3) = inc;
oev(4) = argper;
oev(5) = raan;
oev(6) = tanom;


text:
%%test network
cam1_size=size(test_data_cam1);
cam2_size=size(test_data_cam2);
cam1_feature1=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam1_size(1)/param.test_batch_size)-1
    for n=1:param.test_batch_size
        im_data=imresize((reshape(single(test_data_cam1(m*param.test_batch_size+n,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam1_feature1=[cam1_feature1;(squeeze(net.blobs('pool5/7x7_s1').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size+1:cam1_size(1)
        count=m-(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imresize((reshape(single(test_data_cam1(m,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam1_size(1)-(floor(cam1_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('pool5/7x7_s1').get_data))';
cam1_feature1=[cam1_feature1;result1(1:index,:)];
prob_feature=[];
name_cam1={};
for m=1:param.train_person_num
    test_label_index=find(label_test_cam1==(m+param.train_person_num));
    a=sum(cam1_feature1(test_label_index,:))/size(test_label_index,1);
    prob_feature=[prob_feature;a];
    name_cam1=[name_cam1,test_image_name_cam1{test_label_index(1)}];
end
%extract feature cam2
cam2_feature1=[];
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=0:floor(cam2_size(1)/param.test_batch_size)-1
    for n=1:param.test_batch_size
        im_data=imresize((reshape(single(test_data_cam2(m*param.test_batch_size+n,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,n)=im_data;
        test_batch_data(:,:,1,n)=test_batch_data(:,:,1,n)-104;
        test_batch_data(:,:,2,n)=test_batch_data(:,:,2,n)-117;
        test_batch_data(:,:,3,n)=test_batch_data(:,:,3,n)-123;
    end
    net.blobs('data').set_data(test_batch_data);
    net.forward_prefilled;
    cam2_feature1=[cam2_feature1;(squeeze(net.blobs('pool5/7x7_s1').get_data))'];
end
test_batch_data=zeros(224,224,3,param.test_batch_size,'single');
for m=(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size+1:cam2_size(1)
        count=m-(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size;
        im_data=imresize((reshape(single(test_data_cam2(m,:)),128,64,3)),[224 224]);
        im_data=im_data(:,:,[3,2,1]);
        im_data=permute(im_data,[2,1,3]);
        test_batch_data(:,:,:,count)=im_data;
        test_batch_data(:,:,1,count)=test_batch_data(:,:,1,count)-104;
        test_batch_data(:,:,2,count)=test_batch_data(:,:,2,count)-117;
        test_batch_data(:,:,3,count)=test_batch_data(:,:,3,count)-123;
end
net.blobs('data').set_data(test_batch_data);
net.forward_prefilled;
index=cam2_size(1)-(floor(cam2_size(1)/param.test_batch_size))*param.test_batch_size;
result1=(squeeze(net.blobs('pool5/7x7_s1').get_data))';
cam2_feature1=[cam2_feature1;result1(1:index,:)];
gallery_feature=[];
%%
name_cam2={};
for m=1:param.train_person_num
    test_label_index=find(label_test_cam2==(m+param.train_person_num));
    a=sum(cam2_feature1(test_label_index,:))/size(test_label_index,1);
    gallery_feature=[gallery_feature;a];
    name_cam2=[name_cam2,train_image_name_cam2{test_label_index(1)}];
end
% cal cmc
prob_norm=bsxfun(@rdivide,prob_feature,sum(abs(prob_feature).^2,2).^(1/2));
gallery_norm=bsxfun(@rdivide,gallery_feature,sum(abs(gallery_feature).^2,2).^(1/2));
% [~,it]=sort(prob_norm*gallery_norm',2,'descend');
score_matrix=prob_norm*gallery_norm';
rank1_hit=0;
rank5_hit=0;
rank10_hit=0;
rank20_hit=0;
problen=size(prob_feature,1);
for m=1:problen
    [~,location]=sort(score_matrix(m,:),2,'descend');
    if find(location(1)==m)
        rank1_hit=rank1_hit+1;
    end
    if find(location(1:5)==m)
        rank5_hit=rank5_hit+1;
    end
    if find(location(1:10)==m)
        rank10_hit=rank10_hit+1;
    end
    if find(location(1:20)==m)
        rank20_hit=rank20_hit+1;
    end
end
rank_acc=[rank1_hit/problen,rank5_hit/problen,rank10_hit/problen,rank20_hit/problen];
fin=fopen(param.result_save_file,'a');
fprintf(fin,'split_index:%d,iter:%d, rank1:%f,rank5:%f,rank10:%f,rank20:%f\n',split_index,iter,rank_acc(1),rank_acc(2),rank_acc(3),rank_acc(4));
fprintf('split_index:%d,iter:%d, rank1:%f,rank5:%f,rank10:%f,rank20:%f\n',split_index,iter,rank_acc(1),rank_acc(2),rank_acc(3),rank_acc(4));
fclose(fin);

author: liuyuisanai

repo: Quality-Aware-Network


text:
function im_o   =  Add_noise( im, v )

seed  =  0;
randn( 'state', seed );
noise     =  randn( size(im) );
% noise     =  noise/sqrt(mean2(noise.^2));
im_o      =  double(im) + v*noise;

author: lbasek

repo: image-denoising-benchmark