adaptiveQ.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - adaptiveQ
%    Copyright (C) 2009-2015 Peter Rakyta, Ph.D.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.
%
%> @addtogroup basic Basic Functionalities
%> @{
%> @file adaptiveQ.m
%> @brief A class providing adaptive distribution of the transverse momentum points.
%> @} 
%> @brief A class providing adaptive distribution of the transverse momentum points.
%> A minimal list of transverse momentums (#qvec) and the corresponding function values (#y) of fhandle are dynamicaly created in order to obtain the values of fhandle over the elements of #interpq via interpolation. (Typically #qvec contains significantly less elements than #interpq.)
%> @Available
%> EQuUs v4.6 or later
%%
classdef adaptiveQ < handle & Messages

properties ( Access = protected )
%>  The number of elements in qvec in the first iteration.
    start_qnum
%>  The list of transverse momentums where the function fhandle is desired to be calculated.
    interpq
%>  The interpolated values of fhandle over the elements of interpq.
    interpy 
%>  The dynamically created vector of the transverse momentums.
    qvec
%>  The function values of fhandle calculated over the elements of qvec.
    y 
%>  Maximum of the iterations.
    maxiter
%>  Logical value. If true, the iteration cintinues with the loaded data, or false (default) otherwise.
    resume 
%>  The filename to export the calculated data.
    outFileName
%>  Method of the interpolation (see http://www.mathworks.com/help/matlab/ref/interp1.html )
    interp_method
end



methods ( Access = public )
    
%% adaptiveQ
%> @brief Constructor of the class.
%> @param Opt An instance of the structure Opt.
%> @param interpq An array of transverse momentum points at which fhandle is calculated via interpolation.
%> @param varargin Cell array of optional parameters. For details see @InputParsing.
%> @return An instance of the class
function obj = adaptiveQ( Opt, interpq, varargin )
    obj = obj@Messages( Opt );
    obj.start_qnum = 4;
    obj.interpq = interpq; % the transverse momentums of our interest (usually equidistant).
    obj.interpy = []; % the interpolated function values at momentums given in interpq.    
    obj.qvec    = []; % the adaptive 1D grid of transverse momentums used to iteroplate over interpq
    obj.y       = []; % the function values calculated on the adaptive grid qvec    
    obj.maxiter = 10e2;
    obj.resume        = [];
    obj.outFileName   = [];
    obj.interp_method = [];
    
    obj.InputParsing( varargin{:} )
end

%% runAdaptiveIterations
%> @brief Runs the adaptive iterations to obtain the elemets of qvec where the function fhandle changes the most.
%> @param fhandle The function handle to be calculated over the transverse momentum points: #y = fhandle( #qvec ), where #y is an MxN matrix, and #qvec is a vector of transverse momentums containing of M elements.
    function runAdaptiveIterations( obj, fhandle )
        
        obj.display( '*************************************')
        obj.display( 'Starting adaptive iterations for transverse momentum computations.', 1 );
        obj.display( '*************************************')
        
        if obj.resume && ~isempty(obj.outFileName)
            
            try
                HandleForLoad = LoadFromFile( obj.outFileName );
                obj.qvec = HandleForLoad.LoadVariable('qvec', 'NoEmpty', 'On');
                obj.y = HandleForLoad.LoadVariable('y', 'NoEmpty', 'On');
                obj.interpq = HandleForLoad.LoadVariable('interpq', 'NoEmpty', 'On');
                obj.interpy = HandleForLoad.LoadVariable('interpy', 'NoEmpty', 'On');
                HandleForLoad.ClearLoadedVariables();
            catch err
               obj.display( 'Failed to load data. Restarting calculations.') 
            end
        end
        
        if isempty( obj.interpy )
            obj.initialize( fhandle )
        end
        
        for idx = 1:obj.maxiter
            
            % break the iteration if convergence is reched
            unaccurateq = obj.getUnaccurateq( );
            if isempty(unaccurateq)
                obj.display('convergence reached!!!!!!!!!!');
                break
            end
            
            newq = obj.defNewq( unaccurateq );
            newy = feval( fhandle, newq );
            obj.joinq( newq, newy );
            
            obj.display( ' ');
            obj.display( '*************************************', 1)
            obj.display( ['Iteration ', num2str(idx), ' finished.'], 1 );
            obj.display( '*************************************', 1)
               obj.display( ' ', 1);
            
            if ~isempty( obj.outFileName )
                    qvec = obj.qvec;
                    y = obj.y;
                    interpy = obj.interpy;
                    interpq = obj.interpq;
                    maxiter = obj.maxiter;
                    resume = obj.resume;
                    start_qnum = obj.start_qnum;
                save( obj.outFileName, 'qvec', 'y', 'interpy', 'interpq', 'maxiter', 'resume', 'start_qnum');
            end
        end
        
        
    end

%% initialize
%> @brief  Initializes #y = fhandle( #qvec ) for the first iteration.
%> @param fhandle The function handle to be calculated over the transverse momentum points: #y = fhandle( #qvec ), where #y is an MxN matrix, and #qvec is a vector of transverse momentums containing of M elements.
    function initialize( obj, fhandle )
        
        qmax = max( obj.interpq );
        qmin = min( obj.interpq );
        deltaq = (qmax-qmin)/(obj.start_qnum-1);
        obj.qvec = transpose(qmin:deltaq:qmax);   
        

        obj.y = feval( fhandle, obj.qvec );  
        obj.interpy = NaN(length( obj.interpq), size(obj.y,2) );
        
        % the first and last elements coincide with the iterpolation array
        obj.interpy(1,:) = obj.y(1,:);
        obj.interpy(end,:) = obj.y(end,:);
    end
    
    
%% Write
%> @brief Sets the value of an attribute in the interface.
%> @param MemberName The name of the attribute to be set.
%> @param input The value to be set.
    function Write(obj, MemberName, input)  
       
        try
          obj.(MemberName) = input;
        catch
               error(['EQuUs:', class(obj), ':Read'], ['No property to write with name: ',  MemberName]);
        end
                
    end
%% Read
%> @brief Query for the value of an attribute in the interface.
%> @param MemberName The name of the attribute to be set.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        
        try
          ret = obj.(MemberName);
        catch
               error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ',  MemberName]);
        end            
        
    end
%% Clear
%> @brief Clears the value of an attribute in the interface.
%> @param MemberName The name of the attribute to be cleared.
    function Clear(obj, MemberName)
          try
          obj.(MemberName) = [];
        catch
               error(['EQuUs:', class(obj), ':Clear'], ['No property to clear with name: ',  MemberName]);
        end  
        
    end    



    
end % methods public


methods ( Access = protected )

%% defNewq
%> @brief Defines the new transverse momentums to performs the next iteration.
%> @param unaccurateq The list of transverse moments at which the interpolated function fhandle is unaccurate. 
%> @return Returns with the list of new transverse moments at which the function fhandle should be calculated to improve the accuracy of the interpolation. 
    function newq = defNewq( obj, unaccurateq )
        % removing qs higher than max(oldq)
        indexes = logical( unaccurateq > obj.qvec(end) );
        unaccurateq(indexes) = [];
        
        % identifying the intervals in oldq
        indexes = fix( interp1( obj.qvec, 1:length(obj.qvec), unaccurateq ) );
        indexes = unique(indexes);
        
        % determining the new transverse momentums in the middle of the
        % previously determined intervals
        increments = diff( obj.qvec );
        newq = obj.qvec(indexes) + increments(indexes)/2;
        
    end

%% getUnaccurateq
%> @brief Determine the transverse momentum points at which the interpolated function fhandle is unaccurate.
%> @return Returns with the list of unaccurate transverse moments at which the interpolated function fhandle is unaccurate.
    function unaccurateq = getUnaccurateq( obj )
        %the raltive tolareance of the accuracy
        tolerance = 7e-3;
        
        unaccurateq = [];
        maxvalue = max(max(abs(obj.y)));
        
        for idx = 1:size(obj.y,2)
            interpy_new = interp1( obj.qvec, obj.y(:,idx), obj.interpq, obj.interp_method);
        
            if isempty( maxvalue )
                unaccurateq = obj.interpq;
                return
            end
        
            NaNindexes = isnan( obj.interpy(:,idx) );
        
            relative_differences = abs( (interpy_new - obj.interpy(:,idx))/maxvalue );
            indexes = logical( relative_differences > tolerance );
            unaccurateq = [unaccurateq;obj.interpq( indexes | NaNindexes )];
            obj.interpy(:,idx) = interpy_new;
        end
        
        unaccurateq = unique( unaccurateq );
        obj.display(['Ratio of unaccurate transverse momentums: ', num2str(length(unaccurateq)/length(obj.interpq)*100, '%2.2f' ), '%'], 1);
    end

%% joinq
%> @brief Joins the new transverse moment points with the old ones
%> @param newq The new transverse momentum points.
%> @param newy The values of the function fhandle calculated at points newq.
    function joinq( obj, newq, newy )
        obj.qvec = [obj.qvec; newq];
        obj.y       = [obj.y; newy];
        
        [obj.qvec,IX] = sort( obj.qvec );
        obj.y         = obj.y(IX,:);
        
    end

end % methods protected

methods ( Access = private )

%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'resume' Logical value. If true, the iteration continues with the loaded data, or false (default) otherwise.
%> @param 'outFileName' The filename to save the calculated data. If it is empty (default) no data are exported.
%> @param 'interp_method' The interpolation method to be used to calculate #interpy. The default method is <a href="http://www.mathworks.com/help/matlab/ref/spline.html">spline</a>.
    function InputParsing( obj, varargin )
        
        p = inputParser;
        p.addParameter('resume', false);
        p.addParameter('outFileName', []);
        p.addParameter('interp_method', 'spline', @ischar);
        
        p.parse(varargin{:});

       
        obj.resume        = p.Results.resume;
        obj.outFileName       = p.Results.outFileName;
        obj.interp_method  = p.Results.interp_method;
        
        
    end



end % methods private

end