text:
classdef edgefx
    %   EDGEFX: Edgefunction class
    %   Constructs edgefunctions that are based on numerous geomteric and
    %   topographic criteria that guide the spatially variability of elemental
    %   resolution when building a mesh.
    %   Copyright (C) 2018  Keith Roberts & William Pringle
    %
    %   The following properties are available for input to the edgefx method:
    %      Keyword     Default Val	Variable type
    %      'dis'        defval;     scalar
    %      'fs'         defval;     scalar
    %      'wl'         defval;     scalar/array
    %      'slp'        defval;     scalar/array
    %      'ch'         defval;     scalar
    %      'min_el_ch'  100;     	scalar
    %      'AngOfRe'    60;        	scalar
    %      'max_el'     inf;        scalar/array
    %      'max_el_ns'  inf;     	scalar
    %      'g'          0.20;       scalar/array
    %      'geodata'    defval;     geodata class
    %      'lmsl'       defval;     geodata class
    %      'dt'         -1;         scalar
    %      'fl'         defval;     scalar
    %      'Channels'   defval;   	cell array
    %      'h0'         defval;     scalar
    %      [defval = 0; % placeholder value if arg is not passed.]
    %
    %      The 'wl','slp','max_el','g' options inputs can be scalars to apply
    %      the parameter globally or they can be arrays to bound each parameter
    %      value within topographic elevation bounds, using the following stucture:
    %      [parameter_value1 z_min1 z_max1;
    %       parameter_value2 z_min2 z_max2;
    %       ...
    %       parameter_valuen z_minn z_maxn];
    %
    %       where, z_min is the minimum elevation bound and z_max is the maximum elevation bound.
    %        - similar to axis bounds on matlab plots such as caxis, xlim, ylim, ...
    %        - z is oriented same as the DEM, i.e., +ve value of z is above datum (overland)
    %
    %   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/>.
    properties
        nx % number of x coords
        ny % number of y coords
        x0y0 % 2-tuple of bot. left of domain

        lmsl % lmsl boundary as a geodata class

        bbox
        boubox
        used % edge function keywords

        dis  % percent the edgelength changes in space
        fd   % distance function
        hhd  % matrix containing distance fx values.
        fs   % number of elements to resolve a feature
        fsd  % matrix containing feature size fx values.
        wl   % wavelength
        wld  % matrix containing wavelength fx values.
        slp  % slope edge function
        slpd % matrix containing filtered slope fx values.
        fl   % slope filter parameters
        ch   % channels
        chd  % matrix containing channel fx values.
        AngOfRe % Angle of reslope
        g    % max allowable grade of mesh.
        dt   % theoretical simulateable timestep

        F    % edge function in gridded interpolant format.

        h0    % min. resolution of edgelength function in meters
        gridspace % edgelength function resolution in WGS84 degrees.
        max_el % max. resolution in domain
        max_el_ns % max. resolution +-<0.01 from shoreline
        boudist % distance in WGS84 degrees to boundary of mesh

        min_el_ch % min. element size along the channel
        Channels % cell array of streams.

    end

    methods

        %% Class constructor/pass it a shape class and edgefx params.
        function obj = edgefx(varargin)
            % Check for m_map dir
            M_MAP_EXISTS=0 ;
            if exist('m_proj','file')==2
              M_MAP_EXISTS=1 ;
            end
            if M_MAP_EXISTS~=1
              error('Where''s m_map? Chief, you need to read the user guide')
            end

            % Check for utilties dir
            UTIL_DIR_EXISTS=0 ;
            if exist('inpoly.m','file')
              UTIL_DIR_EXISTS=1 ;
            end
            if UTIL_DIR_EXISTS~=1
              error('Where''s the utilities directory? Chief, you need to read the user guide')
            end

            p = inputParser;

            defval = 0; % placeholder value if arg is not passed.
            % add name/value pairs
            addOptional(p,'dis',defval);
            addOptional(p,'fs',defval);
            addOptional(p,'wl',defval);
            addOptional(p,'slp',defval);
            addOptional(p,'ch',defval);
            addOptional(p,'min_el_ch',100);
            addOptional(p,'AngOfRe',60);
            addOptional(p,'max_el',inf);
            addOptional(p,'max_el_ns',inf);
            addOptional(p,'g',0.20);
            addOptional(p,'geodata',defval)
            addOptional(p,'lmsl',defval)
            addOptional(p,'dt',-1);
            addOptional(p,'fl',defval);
            addOptional(p,'Channels',defval);
            addOptional(p,'h0',defval);

            % parse the inputs
            parse(p,varargin{:});
            % store the inputs as a struct
            inp = p.Results;
            % get the fieldnames of the edge functions
            inp = orderfields(inp,{'max_el','min_el_ch','AngOfRe',...
                'geodata','lmsl','Channels',...
                'dis','fs','fl','g','max_el_ns',...
                'wl','slp','ch','dt','h0'});
            fields = fieldnames(inp);
            % loop through and determine which args were passed.
            % also, assign reasonable default values if some options were

            % not assigned.
            for i = 1 : numel(fields)
                type = fields{i};
                switch type
                    % parse aux options first
                    case('max_el')
                        obj.max_el = inp.(fields{i});
                        if ~isempty(obj.max_el)
                            obj.max_el = inp.(fields{i});
                        end
                    case('max_el_ns')
                        obj.max_el_ns= inp.(fields{i});
                        if obj.max_el_ns ~=0
                            obj.max_el_ns = inp.(fields{i});
                        end
                        % obj.max_el_ns = obj.max_el_ns/111e3;
                    case('g')
                        obj.g= inp.(fields{i});
                        if obj.g ~=0
                            obj.g = inp.(fields{i});
                        end
                    case('fl')
                        obj.fl = inp.(fields{i});
                        if obj.fl ~=0
                            obj.fl = inp.(fields{i});
                        end
                    case('geodata')
                        if isa(inp.(fields{i}),'geodata')
                            feat = inp.(fields{i});
                        end
                    case('lmsl')
                        if isa(inp.(fields{i}),'geodata')
                            obj.lmsl = inp.(fields{i});
                        end
                    case('dt')
                        obj.dt=inp.(fields{i});
                        if obj.dt >= 0
                            obj.dt = inp.(fields{i});
                        end
                    case('h0')
                        obj.h0=inp.(fields{i});
                        if obj.h0 ~= 0
                            obj.h0 = inp.(fields{i});
                        else
                            obj.h0 = feat.h0;
                        end
                    case('min_el_ch')
                        obj.min_el_ch=inp.(fields{i});
                        if obj.min_el_ch ~= 100
                            obj.min_el_ch = inp.(fields{i});
                        end
                    case('AngOfRe')
                        obj.AngOfRe=inp.(fields{i});
                        if obj.AngOfRe ==60
                            obj.AngOfRe = inp.(fields{i});
                        end
                    case('Channels')
                        obj.Channels=inp.(fields{i});
                        if ~isempty(obj.Channels) ~=0
                            obj.Channels = inp.(fields{i});
                        end
                end
            end

            % kjr april 28, 2018-form mesh size grid on-the-fly
            obj.fd       = @dpoly;
            obj.x0y0     = feat.x0y0 + sqrt(eps);
            obj.gridspace = obj.h0/111e3;
            obj.nx       = ceil(abs(feat.x0y0(1)-feat.bbox(1,2))/obj.gridspace);
            obj.ny       = ceil(abs(feat.x0y0(2)-feat.bbox(2,2))/obj.gridspace);
            obj.bbox     = feat.bbox;
            obj.boubox   = feat.boubox;

            % WJP: Do Mercator projection for calculating
            % distances from shoreline
            if feat.bbox(1,2) > 180
                m_proj('mercator','lat', [-88 90], 'long',[0 360]);
            else
                m_proj('mercator','lat', [-88 90],'long',[-180 180]);
            end
            % now turn on the edge functions
            for i = 1 : numel(fields)
                type = fields{i};
                switch type
                    case('dis')
                        obj.dis = inp.(fields{i});
                        if obj.dis ~= 0
                            disp('Building distance function...');
                            obj = distfx(obj,feat);
                            obj.used{end+1} = 'dis';
                        end
                    case('fs')
                        obj.fs  = inp.(fields{i});
                        if obj.fs ~= 0
                            disp('Building feature size function...');
                            try
                                obj = featfx(obj,feat);
                                obj.used{end+1} = 'fs';
                            catch
                                disp('man, no medial points!')
                            end
                        end
                    case('wl')
                        obj.wl  = inp.(fields{i});
                        if obj.wl(1)~=0 && ~isempty(feat.Fb)
                            disp('Building wavelength function...');
                            obj = wlfx(obj,feat);
                            obj.used{end+1} = 'wl';
                        end
                    case('slp')
                        obj.slp  = inp.(fields{i});
                        if obj.slp(1) ~= 0 && ~isempty(feat.Fb)
                            disp('Building slope function...');
                            obj = slpfx(obj,feat);
                            obj.used{end+1} = 'slp';
                        end
                    case('ch')
                        obj.ch  = inp.(fields{i});
                        if obj.ch ~= 0
                            disp('Building channel function...');
                            obj = chfx(obj,feat);
                            obj.used{end+1} = 'ch';
                        end
                    case{'g','geodata','lmsl','max_el','min_el_ch',...
                            'Channels','max_el_ns','h0','dt','fl','AngOfRe'}
                        % dummy to avoid warning
                    otherwise
                        warning('Unexpected edge function name/value pairs.')
                end
                if i==numel(fields)
                    obj = finalize(obj,feat);
                    disp('Finalized edge function!');
                end
            end
        end

        %% Traditional distance function, linear distance from coastline.
        function obj = distfx(obj,feat)
            [d,obj] = get_dis(obj,feat);
            obj.boudist = d;
            obj.hhd = obj.h0 + obj.dis*abs(d);
        end
        %% Feature size function, approximates width of nearshore geo.
        function obj = featfx(obj,feat)
            if isempty(feat.mainland) && isempty(feat.inner)
                % No medial points so use distance function using grade
                warning('No mainland or inner')
                xg = CreateStructGrid(obj);
                obj.fsd = xg*NaN;
                clearvars xg
                return
            end
            % Make sure we don't create a singularity on the coast in the
            % distance function!
            [d,obj] = get_dis(obj,feat);
            obj.boudist = d ;

            [xg,yg] = CreateStructGrid(obj);
            % use a harvestine assumption
            dx = obj.h0*cosd(yg(1,:)); % for gradient function
            dy = obj.h0;               % for gradient function
            % Calculate the gradient of the distance function.
            [ddy,ddx] = EarthGradient(d,dy,dx);
            d_fs = sqrt(ddx.^2 + ddy.^2);
            clearvars ddx ddy
            % WJP: This fix is to put a medial point in narrow channels
            rf = []; cf = [];
            for ii = 2:size(d,1)-1
                for jj = 2:size(d,2)-1
                    if d(ii,jj) < 0
                        if (d(ii-1,jj) >= 0 && ...
                                d(ii+1,jj) >= 0) || ...
                                (d(ii,jj-1) >= 0 && ...
                                d(ii,jj+1) >= 0)
                            rf(end+1) = ii;
                            cf(end+1) = jj;
                        end
                    end
                end
            end
            % Find singularties in the distance function that are
            % within the poly to get the medial axis.
            % Lets only create medial points in places that are
            % sufficiently far away from the coastline to capture the
            % feature and thus aren't spurious.
            d_fs = reshape(d_fs,[],1); d = reshape(d,[],1);
            lg = d_fs < 0.90  & d < -0.5*obj.h0 ;
            [r,c] = ind2sub([obj.nx obj.ny],find(lg));
            % create the coordinates using x0y0 + r*h0
            r = [r; rf']; c = [c; cf']; %WJP: add on "fixed points"
            x_kp = obj.x0y0(1) + (r-1)*obj.gridspace;
            y_kp = obj.x0y0(2) + (c-1)*obj.gridspace;
            clearvars lg r c cf rf
            % Ensure there are enough medial points nearby
            % to avoid surpurious medial points.
            % WJP: We want a line of points larger than around 7. This
            % corresponds to having three points back or forward up the
            % line. Let's check for those three closest points and ensure
            % they are within about co*h0 distance from each other.
            % co is the cutoff distance = 2*sqrt(2) (two diagonal dist)
            co = 2*sqrt(2);
            if length(x_kp) > 12
                [~, dmed] = WrapperForKsearch([x_kp,y_kp],[x_kp,y_kp],4);
                prune = dmed(:,2) > co*obj.h0 | dmed(:,3) > 2*co*obj.h0 | ...
                    dmed(:,4) > 3*co*obj.h0;
                x_kp( prune ) = [];
                y_kp( prune ) = [];
            end

            % Now get the feature size along the coastline
            if length(x_kp) > 12
                % Use KD-tree
                dPOS = 0*d;
                noblks = ceil(length(dPOS)*2*8*1e-9);
                blklen = floor(length(dPOS)/noblks);
                ns = 1;
                for blks = 1:noblks
                    if blks == noblks
                        ne = length(dPOS);
                    else
                        ne = ns + blklen - 1;
                    end
                    [~,dPOS(ns:ne)] = WrapperForKsearch([x_kp,y_kp],...
                        [xg(ns:ne)',yg(ns:ne)'],1);
                    ns = ne + 1;
                end
            else
                % No medial points so use distance function using grade
                warning('No medial points, resorting to distance function')
                d = reshape(d,obj.nx,[]);
                obj.fsd = obj.h0 + obj.g*abs(d);
                clear x_kp y_kp d d_fs dPOS xg yg
                return
            end
            % reshape back
            d = reshape(d,obj.nx,[]); dPOS = reshape(dPOS,obj.nx,[]);
            % Feature_size is distance from medial axis plus distance to
            % coastline. min_el is then feature_size*2/R where R is
            % number of elements to model the feature
            % WJP: d needs to be the absolute value because when we
            % include the floodplain since dPOS is not signed but d is,
            % they will cancel out producing really fine resolution
            W = dPOS+abs(d);
            if obj.fs < 0
                % automatic feature size
                fs_auto = min(-obj.fs,ceil(W/obj.h0));
                obj.fsd = 2*W./fs_auto;
            else
                obj.fsd = 2*W/obj.fs;
            end
            clear x_kp y_kp d d_fs dPOS xg yg
        end

        function [d,obj] = get_dis(obj,feat)
            % Function used by distfx and featfx to return distance
            % and make the Fdis interpolant
            d = feval(obj.fd,obj,feat);
            d = reshape(d,obj.nx,[]);
        end

        %% Wavelength edge function.
        function obj = wlfx(obj,feat)

            % interpolate DEM's bathy linearly onto our edgefunction grid.
            [xg,yg] = CreateStructGrid(obj);
            tmpz    = feat.Fb(xg,yg);
            if ~isempty(feat.Fb2)
               tmpz2 = feat.Fb2(xg,yg);
               tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
               clear tmpz2
            end
            % Initialise wld
            obj.wld = NaN([obj.nx,obj.ny]);
            grav    = 9.807;
            period  = 12.42*3600; % M2 period in seconds

            for param = obj.wl'
                if numel(param)==1
                    % no bounds specified.
                    wlp = param(1);
                    z_min = -inf;
                    z_max = +inf;
                else
                    wlp = param(1);
                    z_min = param(2);
                    z_max = param(3);
                end
                % limit min. depth to 1 m
                twld = period*sqrt(grav*max(abs(tmpz),1))/wlp;
                limidx = (tmpz >= z_min & tmpz < z_max);
                % Set wld with depth mask applied
                obj.wld(limidx) = twld(limidx);
                clearvars twld
            end
            clearvars tmpz xg yg;
        end

      %% Topographic length scale/slope edge function.
        function obj = slpfx(obj,feat)

            [xg,yg] = CreateStructGrid(obj);

            tmpz    = feat.Fb(xg,yg);
            if ~isempty(feat.Fb2)
               tmpz2 = feat.Fb2(xg,yg);
               tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
               clear tmpz2
            end
            tmpz(tmpz > 50) = 50; % ensure no larger than 50 m above land
            % use a Harvestine assumption
            dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function
            dy = obj.h0;               % for gradient function
            % lets filter the bathy to get only relevant features
            % loop over each set of bandpass filter lengths
            tmpz_f = zeros(size(tmpz));
            if obj.fl(1) < 0 && obj.fl(1) ~= -999
                disp('INFO: Rossby radius of deformation filter on.') ;
                rbfilt = abs(obj.fl(1));
                barot  = 1;
                if length(obj.fl) > 1 && obj.fl(2) < 0
                    disp('INFO: Using 1st-mode baroclinic Rossby radius.');
                    barot = 0;
                else
                    disp('INFO: Using barotropic Rossby radius.') ;
                end
                obj.fl = [];
                filtit = 1;
            elseif obj.fl(1) == 0
                disp('INFO: Slope filter is off.');
                obj.fl = [];
                tmpz_f = tmpz;
                filtit = 0 ;
            elseif obj.fl(1) == -999
                disp('INFO: Slope filter is LEGACY.');
                obj.fl = [];
                tmpz_f = tmpz;
                filtit = -999;
            else
                filtit = 0;
                for lambda = obj.fl'
                    if length(lambda) == 1 || lambda(2) == 0
                        % do a low pass filter
                        tmpz_ft  = filt2(tmpz,dy,lambda(1),'lp') ;
                    elseif all(lambda ~= 0)
                        % do a bandpass filter
                        tmpz_ft  = filt2(tmpz,dy,lambda,'bp') ;
                    else
                        % highpass filter not recommended
                        warning(['Highpass filter on bathymetry in slope' ...
                            'edgelength function is not recommended'])
                        tmpz_ft  = filt2(tmpz,dy,lambda(2),'hp') ;
                    end
                    tmpz_f = tmpz_f + tmpz_ft;
                end
            end

            % Rossby radius of deformation filter
            if filtit == 1
                tic
                bs = NaN([obj.nx,obj.ny]);
                % break into 10 deg latitude chunks, or less if higher res
                div = ceil(min(1e7/obj.nx,10*obj.ny/(max(yg(:))-min(yg(:)))));
                grav = 9.807;
                nb = ceil(obj.ny/div);
                n2s = 1;
                for jj = 1:nb
                    n2e = min(obj.ny,n2s + div - 1);
                    % Rossby radius of deformation filter
                    % See Shelton, D. B., et al. (1998): Geographical variability of the first-baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.
                    ygg = yg(:,n2s:n2e);
                    dxx = mean(dx(n2s:n2e));
                    f = 2*7.29e-5*abs(sind(ygg));
                    if barot
                        % barotropic case
                        c = sqrt(grav*max(1,-tmpz(:,n2s:n2e)));
                    else
                        % baroclinic case (estimate Nm to be 2.5e-3)
                        Nm = 2.5e-3;
                        c = Nm*max(1,-tmpz(:,n2s:n2e))/pi;
                    end
                    rosb = c./f;
                    clear f;
                    % update for equatorial regions
                    I = abs(ygg) < 5; Re = 6371e3;
                    twobeta = 4*7.29e-5*cosd(ygg(I))/Re;
                    rosb(I) = sqrt(c(I)./twobeta);
                    clear twobeta
                    % limit rossby radius to 10,000 km for practical purposes
                    rosb(rosb > 1e7) = 1e7;
                    % Keep lengthscales rbfilt * barotropic
                    % radius of deformation
                    rosb = min(10,max(0,floor(log2(rosb/dy/rbfilt))));
                    edges = double(unique(rosb(:)));
                    bst = rosb*0;
                    for i = 1:length(edges)
                        if edges(i) > 0
                            mult = 2^edges(i);
                            xl = 1; xu = obj.nx;
                            if (max(xg(:)) > 179 && min(xg(:)) < -179) || ...
                                    (max(xg(:)) > 359 && min(xg(:)) < 1)
                                % wraps around
                                xr = [obj.nx-mult/2+1:1:obj.nx xl:xu 1:mult/2];
                            else
                                xr = xl:xu;
                            end
                            yl = max(1,n2s-mult/2);
                            yu = min(obj.ny,n2e+mult/2);
                            if max(yg(:)) > 89 && yu == obj.ny
                                % create mirror around pole
                                yr = [yl:yu yu-1:-1:2*obj.ny-n2e-mult/2];
                            else
                                yr = yl:yu;
                            end
                            if mult == 2
                                tmpz_ft = filt2(tmpz(xr,yr),...
                                           [dxx dy],dy*2.01,'lp');
                            else
                                tmpz_ft = filt2(tmpz(xr,yr),...
                                           [dxx dy],dy*mult,'lp');
                            end
                            % delete the padded region
                            tmpz_ft(1:find(xr == 1,1,'first')-1,:) = [];
                            tmpz_ft(obj.nx+1:end,:) = [];
                            tmpz_ft(:,1:find(yr == n2s,1,'first')-1) = [];
                            tmpz_ft(:,n2e-n2s+2:end) = [];
                        else
                            tmpz_ft = tmpz(:,n2s:n2e);
                        end
                        [by,bx] = EarthGradient(tmpz_ft,dy,dx(n2s:n2e)); % get slope in x and y directions
                        tempbs  = sqrt(bx.^2 + by.^2);          % get overall slope
                        bst(rosb == edges(i)) = tempbs(rosb == edges(i));
                    end
                    bs(:,n2s:n2e) = bst;
                    n2s = n2e + 1;
                end
                toc
                % legacy filter
            elseif filtit==-999
                bs = NaN([obj.nx,obj.ny]);
                % Rossby radius of deformation filter
                f = 2*7.29e-5*abs(sind(yg));
                % limit to 1000 km
                rosb = min(1000e3,sqrt(9.81*abs(tmpz))./f);
                % autmatically divide into discrete bins
                [~,edges] = histcounts(rosb);
                tmpz_ft  = tmpz; dyb = dy;
                % get slope from filtered bathy for the segment only
                [by,bx] = EarthGradient(tmpz_ft,dy,dx); % get slope in x and y directions
                tempbs  = sqrt(bx.^2 + by.^2); % get overall slope
                for i = 1:length(edges)-1
                    sel = rosb >= edges(i) & rosb <= edges(i+1);
                    rosbylb = mean(edges(i:i+1));
                    if rosbylb > 2*dyb
                        disp(['i = ',num2str(i), ' rl/dx = ',num2str(rosbylb/dyb)])
                        tmpz_ft  = filt2(tmpz_ft,dyb,rosbylb,'lp');
                        dyb = rosbylb;
                        % get slope from filtered bathy for the segment only
                        [by,bx] = EarthGradient(tmpz_ft,dy,dx); % get slope in x and y directions
                        tempbs  = sqrt(bx.^2 + by.^2); % get overall slope
                    else
                        % otherwise just use the same tempbs from before
                    end
                    % put in the full one
                    bs(sel) = tempbs(sel);
                end
            else
                % get slope from (possibly filtered) bathy
                [by,bx] = EarthGradient(tmpz_f,dy,dx); % get slope in x and y directions
                bs      = sqrt(bx.^2 + by.^2); % get overall slope
            end
            clear bx by tmpz_f tmpz_ft

            % Allow user to specify depth ranges for slope parameter.
            obj.slpd = NaN([obj.nx,obj.ny]);
            for param = obj.slp'
                if numel(param)==1
                    % no bounds specified. valid in this range.
                    slpp = param(1);
                    z_min = -inf;
                    z_max = +inf;
                else
                    slpp = param(1);
                    z_min = param(2);
                    z_max = param(3);
                end
                % Calculating the slope function
                dp = max(1,-tmpz);
                tslpd = (2*pi/slpp)*dp./(bs+eps);
                % apply slope with mask
                limidx = (tmpz >= z_min & tmpz < z_max);
                obj.slpd(limidx) = tslpd(limidx);
                clearvars tslpd
            end
            clearvars tmpz xg yg
        end

        %% Channel edgefunction
        function obj = chfx(obj,feat)

            ang_of_reslope=obj.AngOfRe ;      % Estimate width of channel using tangent of this angle times depth.

            % STEP 1: Calculate the width of each channel using a v-shape approx to
            % channel's cross sectional area.
            for jj = 1:length(obj.Channels)
                if (isempty(obj.Channels{jj})); continue ;end
                pts{jj} = obj.Channels{jj};
                dp{jj} = feat.Fb(pts{jj});                               % depth at x,y channel locations
                radii{jj}=(tand(ang_of_reslope)*max(1,-dp{jj}));         % estimate of channel's width in degrees at x,y locations ASSUMING angle of reslope
                tempbb{jj} = feat.boubox;
            end

            [xg,yg] = CreateStructGrid(obj);

            % STEP 2: For each channel point, set the resolution around a
            % neighborhood of points as |h|/ch where ch is non-dim # between 0.5-1.5
            obj.chd = NaN([obj.nx,obj.ny]);
            % prune the potentially large channel network to only ones
            % partially enclosed in boubox
            isin = cellfun(@inpolysum,pts,tempbb);
            pts(isin==0) = []; radii(isin==0)= [];

            %figure; plot(xg,yg,'k.');
            tic
            tidx = [];
            for jj =1:length(pts)                                           % For each channel jj
                sel = pts{jj};
                in  = inpoly(sel,feat.boubox(1:end-1,:));
                sel(~in,:) = [];
                for jjj = 1 : size(sel,1)                                   % For each channel point jjj on channel jj.
                    % will the stencil be too large (nidx^2)?
                    nidx = ceil(radii{jj}(jjj)/obj.h0);                       % stencil of points around channel point.
                    % find linear index of channel point in 2-d matrix.
                    [lidx] = FindLinearIdx(sel(jjj,1),sel(jjj,2),xg,yg);
                    % convert to r,c or grid indices.
                    [r,c]     = ind2sub(size(xg),lidx);
                    [cid,rid] = ndgrid(c-nidx:c+nidx,r-nidx:r+nidx);    % grid of nearby points
                    % ensure that all these points are within the domain.
                    rid = min(max(rid(:),1),size(xg,1));
                    cid = min(max(cid(:),1),size(xg,2));
                    % convert back to linear indices.
                    temp = sub2ind(size(xg),rid,cid);
                    tidx = [tidx; temp];
                end
            end
            toc
            tmpz = feat.Fb(xg(tidx),yg(tidx));
            if ~isempty(feat.Fb2)
               tmpz2 = feat.Fb2(xg(tidx),yg(tidx));
               tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
               clear tmpz2
            end
            dp = max(1,-tmpz);
            obj.chd(tidx) = dp/obj.ch;
            obj.chd(obj.chd < obj.min_el_ch) = obj.min_el_ch;
            clearvars Fb pts dp radii tempbb xg yg
        end

        %% General function to plot
        function plot(obj,type)
            % form grid here.
            % working out plottable size
            [xgrid,ygrid] = CreateStructGrid(obj);

            if nargin == 2
                switch type
                    case('dis')
                        val = obj.hhd;
                        tt = 'Distance EdgeLength Function';
                    case('fs')
                        val = obj.fsd;
                        tt = 'Feature Size EdgeLength Function';
                    case('wl')
                        val = obj.wld;
                        tt = 'Wavelength EdgeLength Function';
                    case('ch')
                        val = obj.chd;
                        tt = 'Channel EdgeLength Function';
                    case('slp')
                        val = obj.slpd;
                        tt = 'Slope EdgeLength Function';
                    otherwise
                        warning('Unexpected plot type. No plot created.')
                end
            else
                val = obj.F.Values;
                tt = 'Total EdgeLength Function';
            end

            % working out stride
            mem = inf; stride = obj.gridspace;
            while mem > 1
                xx = obj.x0y0(1):stride:obj.bbox(1,2);
                yy = obj.x0y0(2):stride:obj.bbox(2,2);
                xs = whos('xx'); ys = whos('yy');
                mem = xs.bytes*ys.bytes/1e9;
                stride = stride*2;
            end
            stride = ceil(stride/obj.gridspace);
            I = [1:stride:size(xgrid,1)];
            J = [1:stride:size(xgrid,2)];

            figure;
            m_proj('merc',...
                'long',[obj.bbox(1,1) obj.bbox(1,2)],...
                'lat',[obj.bbox(2,1) obj.bbox(2,2)])
            hold on; m_fastscatter(xgrid(I,J),ygrid(I,J),val(I,J));
            cb = colorbar; ylabel(cb,'edgelength in meters');
            colormap(lansey)
            m_grid() %'xtick',10,'tickdir','out',...
            % 'yaxislocation','left','fontsize',7);
            title(tt);
        end

        %% Finalize edge function
        % Add grading, cfl limiting and bound enforcement.
        function obj = finalize(obj,feat)
            % package the edge functions into a known order.
            counter = 0;
            hh = zeros([obj.nx,obj.ny]);
            for i = 1 : numel(obj.used)
                type = obj.used{i};
                switch type
                    case('dis')
                        counter = counter + 1;
                        hh(:,:,counter) = obj.hhd;
                        obj.hhd = single(obj.hhd);
                    case('fs')
                        counter = counter + 1;
                        hh(:,:,counter) = obj.fsd;
                        obj.fsd = single(obj.fsd);
                    case('wl')
                        counter = counter + 1;
                        hh(:,:,counter) = obj.wld;
                        obj.wld = single(obj.wld);
                    case('slp')
                        counter = counter + 1;
                        hh(:,:,counter) = obj.slpd;
                        obj.slpd = single(obj.slpd);
                    case('ch')
                        counter = counter + 1;
                        hh(:,:,counter) = obj.chd;
                        obj.ch = single(obj.chd);
                    otherwise
                        error('FATAL:  Could not finalize edge function');
                end
            end

            [hh_m] = min(hh,[],3);
            clearvars hh

            [xg,yg] = CreateStructGrid(obj);

            % Before grading, ensure edgefx reflects the spacing in
            % proximity to constraints.
            if ~isempty(feat.weirs)
                % form temporary mesh size function interpolat
                % if spacing is negative -999 then it will use the mesh
                % size function to determine the spacing.
                %Ftmp = griddedInterpolant(xg,yg,hh_m,'linear','nearest');

                for i =1:length(feat.weirs)
                    if iscell(feat.weirs)
                        % get spacing in meters along the face of the weir
                        weir_spacing = feat.weirs{i}(1,4);
                        % get crestline points in wgs84
                        txy = feat.weirs{i}(:,1:2) ;
                    else
                        % get spacing in meters along the face of the weir
                        weir_spacing = feat.weirs(i).min_ele;
                        % get crestline points in wgs84
                        txy = [feat.weirs(i).X feat.weirs(i).Y];
                    end
                    xy = [];
                    [xy(:,2),xy(:,1)] = my_interpm(txy(:,2),txy(:,1),...
                        (obj.h0/2)/111e3);
                    % edit the mesh size function in proximity to the weir
                    for j = 1 : length(xy)
                        nidx = 10 ;
                        % return the indices into the mesh size function
                        [lidx] = FindLinearIdx(xy(j,1),xy(j,2),xg,yg);
                        [r,c]  = ind2sub(size(xg),lidx);
                        [cid,rid]=ndgrid((c-nidx:c+nidx)',(r-nidx:r+nidx)');    % grid of nearby points
                        rid = [rid(:);r]; cid=[cid(:);c];
                        hh_m(rid,cid) = weir_spacing ;
                    end
                end
            end

            % enforce all mesh resolution bounds,grade and enforce the CFL in planar metres
            if ~isinf(obj.max_el_ns) && ~isempty(obj.boudist)
                nearshore = abs(obj.boudist) < 2/3*obj.h0;
                hh_m(nearshore & hh_m > obj.max_el_ns) = obj.max_el_ns;
            end
            hh_m(hh_m < obj.h0 ) = obj.h0;

            % Compute tmpz if necessary
            if length(obj.max_el) > 1 || length(obj.g) > 1 || obj.dt(1) >= 0
                tmpz    = feat.Fb(xg,yg);
                if ~isempty(feat.Fb2)
                   tmpz2 = feat.Fb2(xg,yg);
                   tmpz(isnan(tmpz)) = tmpz2(isnan(tmpz));
                   clear tmpz2
                end
            end

            for param = obj.max_el'
                if numel(param)==1 && param~=0
                    mx   = obj.max_el(1);
                    limidx = hh_m > mx | isnan(hh_m);

                    hh_m( limidx ) = mx;
                else
                    mx  = param(1);
                    z_min = param(2);
                    z_max = param(3);

                    limidx = (tmpz >= z_min & tmpz < z_max) & ...
                             (hh_m > mx | isnan(hh_m));

                    hh_m( limidx ) = mx;
                end
            end

            % Make sure this is called before releasing memory...
            if obj.dt(1) == 0
                % Find min allowable dt based on dis or fs function
                if any(~cellfun('isempty',strfind(obj.used,'dis')))
                    hh_d = obj.hhd;
                elseif any(~cellfun('isempty',strfind(obj.used,'fs')))
                    hh_d = obj.fsd;
                else
                    error(['FATAL: cannot use automatic timestep ' ...
                           'limiter without specifying dis or fs option']);
                end
            end

            obj = release_memory(obj) ;

            disp('Relaxing the mesh size gradient');
            if all(abs(obj.bbox(1,:)) == 180)
                % for global mesh make it cyclical
                hh_m = [hh_m(end,:); hh_m; hh_m(1,:)];
            end
            hfun = reshape(hh_m',[numel(hh_m),1]);

            dx = obj.h0*cosd(min(yg(1,:),85)); % for gradient function
            dy = obj.h0;                       % for gradient function

            % make g a function of space
            dmy     = xg*0 ;
            for param = obj.g'
                if numel(param)==1 && param~=0
                    lim   = obj.g(1);
                    dmy  = dmy + lim ;
                else
                    lim  = param(1);
                    z_min  = param(2);
                    z_max  = param(3);

                    limidx = (tmpz >= z_min & tmpz < z_max) ;

                    dmy( limidx ) = lim;
                end
            end
            if all(abs(obj.bbox(1,:)) == 180)
                % for global mesh make it cyclical
                dmy = [dmy(end,:); dmy; dmy(1,:)];
            end
            fdfdx = reshape(dmy',[numel(dmy),1]);
            clearvars dmy;

            [hfun,flag] = limgradStruct(obj.ny,dx,dy,hfun,...
                fdfdx,sqrt(length(hfun)));
            if flag == 1
                disp('Gradient relaxing converged!');
            else
                error(['FATAL: Gradient relaxing did not converge, '
                    'please check your edge functions']);
            end
            % reshape it back
            hh_m = reshape(hfun,obj.ny,[])';
            clearvars hfun fdfdx
            if all(abs(obj.bbox(1,:)) == 180)
                % for global mesh make it cyclical
                hh_m = hh_m(2:end-1,:);
                hh_m(end,:) = hh_m(1,:);
            end

            % Limit CFL if dt >= 0, dt = 0 finds dt automatically.
            if obj.dt(1) >= 0
                if isempty(feat.Fb)
                    error('No DEM supplied Can''t CFL limit.')
                end
                % Estimate wavespeed in ocean (second term represents orbital
                % velocity at 0 degree phase for 1-m amp. wave).
                grav = 9.807;
                % Limit the minimum depth to 1 m (ignore overland)
                tmpz(tmpz > - 1) = -1;
                u = sqrt(grav*abs(tmpz)) + sqrt(grav./abs(tmpz));

                % Desired timestep to bound edgefx
                desDt = obj.dt(1);
                if isnan(desDt)
                    disp('No timestep enforcement due to no distance function')
                end

                % Determine Courant min. and max. bounds
                if length(obj.dt) < 2
                    % default behavior if no bounds are passed
                    % for explicit type schemes we ensure the maxCr is
                    % bounded
                    minCr = eps; % Cannot be zero!!!
                    maxCr = 0.5;
                elseif length(obj.dt) == 2
                    if obj.dt(1)==0.0
                        error('dt = 0 is only intended for restricting the upper bound of Cr')
                    end
                    % minimum Courant number bound
                    minCr = obj.dt(2);
                    maxCr = inf;
                elseif length(obj.dt) == 3
                    % minimum and maximum Courant bound
                    minCr = obj.dt(2);
                    maxCr = obj.dt(3);
                else
                    error('Length of dt name-value pair should be <=3')
                end

                if minCr > maxCr
                    error('MinCr is > MaxCr, please switch order in dt name value call')
                end

                % Automatic timestep selection option (for ADCIRC models)
                if desDt == 0
                    hh_d(hh_d < obj.h0) = obj.h0;
                    obj.dt = min(min(maxCr*hh_d./u));
                    desDt  = obj.dt;
                end

                % Enforcement of Courant bounds in mesh size function
                fprintf(1, [ ...
                    'Enforcing timestep of ',num2str(desDt),' seconds ', ...
                    'with Courant number bounds of ',num2str(minCr),' to ',num2str(maxCr),'\n', ...
                    ] ) ;

                Cr = (desDt*u)./hh_m; % Courant number for given desDt for mesh size variations

                % resolution that meets max. Cr. criteria for desDt
                dxn_min = u*desDt/maxCr;
                % resolution that meets min. Cr. criteria for desDt
                dxn_max = u*desDt/minCr;

                % resolve min. Cr violations
                hh_m( Cr <= minCr) = dxn_max( Cr <= minCr); %--in planar metres

                % resolve max. Cr violations
                hh_m( Cr > maxCr) = dxn_min( Cr > maxCr);   %--in planar metres

                clear dxn_min dxn_max Cr
                clear u hh_d tmpz
            end

            obj.F = griddedInterpolant(xg,yg,hh_m,'linear','nearest');

            clearvars xg yg
        end

        function [xg,yg] = CreateStructGrid(obj)
            [xg,yg] = ndgrid(obj.x0y0(1) + (0:obj.nx-1)'*obj.gridspace, ...
                obj.x0y0(2) + (0:obj.ny-1)'*obj.gridspace);
        end

        function obj = release_memory(obj)
            % releases heavy components from data structures
            if isa(obj,'edgefx')
                disp('--------------------------------------------------');
                disp('Releasing individual edge functions from memory...');
                disp('--------------------------------------------------');
                if ~isempty(obj.fsd)
                    obj.fsd = [];
                end

                if ~isempty(obj.wld)
                    obj.wld = [];
                end

                if ~isempty(obj.slpd)
                    obj.slpd = [];
                end

                if ~isempty(obj.chd)
                    obj.chd = [];
                end

                if ~isempty(obj.ch)
                    obj.ch = [];
                end

                if ~isempty(obj.hhd)
                    obj.hhd = [];
                end

                if ~isempty(obj.boudist)
                    obj.boudist = [];
                end

                if ~isempty(obj.Channels)
                    obj.Channels = [];
                end

                if ~isempty(obj.ch)
                   obj.ch = [];
                end

                if ~isempty(obj.boudist)
                   obj.boudist = [];
                end
            end
        end

    end

end

author: CHLNDDEV


text:
classdef CellsClassifier < handle

    methods (Access = public, Static)

        function [full,empty,cut] = classifyCells(levelSet,connec)
            nodes = connec;
            indexes = (1:size(connec,1))';
            lsNodes = levelSet(nodes);
            phiCase = sum((sign(lsNodes)<0),2);

            full  = phiCase == size(nodes,2);
            empty = phiCase == 0;
            cut = indexes(~(full | empty));
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

function EU = SimpleCalcExpectedUtility(I)

  % Inputs: An influence diagram, I (as described in the writeup).
  %         I.RandomFactors = list of factors for each random variable.  These are CPDs, with
  %              the child variable = D.var(1)
  %         I.DecisionFactors = factor for the decision node.
  %         I.UtilityFactors = list of factors representing conditional utilities.
  % Return Value: the expected utility of I
  % Given a fully instantiated influence diagram with a single utility node and decision node,
  % calculate and return the expected utility.  Note - assumes that the decision rule for the
  % decision node is fully assigned.

  % In this function, we assume there is only one utility node.
  F = [I.RandomFactors I.DecisionFactors];
  U = I.UtilityFactors(1);
  EU = [];
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % YOUR CODE HERE
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Fnew = VariableElimination(F,setdiff(unique([F.var]),U.var));
  %Fnew = F;
  out = Fnew(1);
  for i = 2:length(Fnew)
	  out = FactorProduct(out,Fnew(i));
  end

  hha = U;
  hha.val = ones(1,length(U.val));
  out = FactorProduct(hha,out);
  U = FactorProduct(hha,U);
  %out = FactorMarginalization(out,setdiff(unique([F.var]),U.var));

EU = [sum(out.val.*U.val)];

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function val=tanDerivVal(z,n)
%%TANDERIVVAL Return the value of the nth derivative of tan(x) with respect
%             to x evaluated at x=z.
%
%INPUTS: z A matrix of real or complex values at which the nth derivative
%          of the tangent function is desired.
%        n The number of derivatives to take. n>=0.
%
%OUTPUTS: val The value of the nth derivative of the tangent function taken
%             at all of the points in z. val has the same dimensions as z.
%
%This function implements the algorithm of [1].
%
%REFERENCES:
%[1] M. E. Hoffman, "Derivative polynomials for tangent and secant," The
%    American Mathematical Monthly, vol. 102, no. 1, pp. 23-30, Jan. 1995.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numEls=numel(z);

val=zeros(size(z));
for curEl=1:numEls
    x=z(curEl);
    u=tan(x);

    P=zeros(n+1,1);
    P(1)=u;

    %The loop implements Equation 5.
    for k=1:n
        sumVal=0;
        for i=0:(k-1)
            sumVal=sumVal+binomial(k-1,i)*P(i+1)*P(k-1-i+1);
        end
        if(k==1)
            sumVal=sumVal+1;
        end
        P(k+1)=sumVal;
    end
    val(curEl)=P(end);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function s = s_blank_delete ( s )

%*****************************************************************************80
%
%% S_BLANK_DELETE removes blanks from a string, left justifying the remainder.
%
%  Discussion:
%
%    All TAB characters are also removed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be transformed.
%
%    Output, string S, the transformed string.
%
  TAB = 9;

  s_len = length ( s );
  s2 = '';

  put = 0;

  for get = 1 : s_len

    if ( s(get) ~= ' ' && s(get) ~= TAB )
      put = put + 1;
      s2(put) = s(get);
    end

  end

  s = s2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function theta = lineAngle(varargin)
%LINEANGLE Computes angle between two straight lines.
%
%   A = lineAngle(LINE);
%   Returns the angle between horizontal, right-axis and the given line.
%   Angle is fiven in radians, between 0 and 2*pi, in counter-clockwise
%   direction.
%
%   A = lineAngle(LINE1, LINE2);
%   Returns the directed angle between the two lines. Angle is given in
%   radians between 0 and 2*pi, in counter-clockwise direction.
%
%   See also
%   lines2d, angles2d, createLine, normalizeAngle
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   2004-02-19 added support for multiple lines.
%   2011-01-20 use bsxfun

nargs = length(varargin);
if nargs == 1
    % angle of one line with horizontal
    line = varargin{1};
    theta = mod(atan2(line(:,4), line(:,3)) + 2*pi, 2*pi);

elseif nargs==2
    % angle between two lines
    theta1 = lineAngle(varargin{1});
    theta2 = lineAngle(varargin{2});
    theta = mod(bsxfun(@minus, theta2, theta1)+2*pi, 2*pi);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function test_ft_preproc_dftfilter

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_preproc_dftfilter

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  fn = localfunctions;
  for i=1:numel(fn)
    feval(fn{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testIssue1770(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% this is an issue related to numerical precision which is too strictly
% checked in ft_preproc_dftfilter (with spectral interpolation)

fs = 500;
data = [];
data.time{1} = -1:1/fs:2.4980; %-> observation: using linspace instead works fine

lnoise(1,:) = (1+hanning(1750))'.*sin((2.*pi.*data.time{1}).*50);
lnoise(2,:) = (1+hanning(1750))'.*sin((2.*pi.*data.time{1}).*100);

data.trial{1} = randn(2,1750)+lnoise;
data.label = {'a';'b'};
data.fsample = fs;

lineFreq = 50;

cfg              = [];
cfg.dftfilter    = 'yes'; % apply line noise filter with spectrum interpolation
cfg.dftfreq      = [lineFreq lineFreq*2]; % line noise and harmonic
cfg.dftreplace   = 'neighbour'; % spectral interpolation
cfg.dftbandwidth = [1 2]; % width of window to be interpolated
cfg.dftneighbourwidth = [2 2]; % width of window from which to interpolate
datafilt1 = ft_preprocessing(cfg, data);

cfg.dftreplace   = 'neighbour_fft';
datafilt2 = ft_preprocessing(cfg, data); % this should now work thanks to some eps leniency

cfg           = [];
cfg.dftfilter = 'yes';
cfg.dftfreq   = [lineFreq lineFreq*2];
datafilt3     = ft_preprocessing(cfg, data);

figure

subplot(2,2,1);plot(datafilt1.time{1}, data.trial{1}-datafilt1.trial{1}); ylim([-2.1 2.1]);xlabel('estimated linenoise neighbour');
subplot(2,2,2);plot(datafilt2.time{1}, data.trial{1}-datafilt2.trial{1}); ylim([-2.1 2.1]);xlabel('estimated linenoise neighbour_fft','interpreter','none');
subplot(2,2,3);plot(datafilt2.time{1}, data.trial{1}-datafilt3.trial{1}); ylim([-2.1 2.1]);xlabel('estimated linenoise static dft');
subplot(2,2,4);plot(data.time{1}, lnoise); ylim([-2.1 2.1]);xlabel('simulated linenoise');

figure; plot(data.time{1}, datafilt1.trial{1}-datafilt2.trial{1});
% the difference here can be explained by the fact that neighbour_fft takes
% an asymmetric band around 50/100 Hz, due to rounding in nearest I guess

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some other instances, testing ft_preproc_dftfilter directly
tim = (0:1999)./1000;
dat = randn(1, 2000) + (1+hanning(2000))'.*sin(2.*pi.*(tim).*50);
filt  = ft_preproc_dftfilter(dat, 1000, 50, 'dftreplace', 'neighbour');
filt2 = ft_preproc_dftfilter(dat, 1000, 50, 'dftreplace', 'neighbour_fft');

figure;plot(tim, dat-filt); hold on;plot(tim, (1+hanning(2000))'.*sin(2.*pi.*(tim).*50));
figure;plot(tim, dat-filt2); hold on;plot(tim, (1+hanning(2000))'.*sin(2.*pi.*(tim).*50));
figure;plot(filt, filt2, 'o');

tim = (0:2000)./1000;
dat = randn(1, 2001) + (1+hanning(2001))'.*sin(2.*pi.*(tim).*50);
filt = ft_preproc_dftfilter(dat, 1000, 50, 'dftreplace', 'neighbour');
try
  filt2 = ft_preproc_dftfilter(dat, 1000, 50, 'dftreplace', 'neighbour_fft');
  ft_error('if the code ends up here, then something suddenly started working');
catch
  % this is supposed to happen
end

figure;plot(tim, dat-filt); hold on;plot(tim, (1+hanning(2001))'.*sin(2.*pi.*(tim).*50));

dat = randn(3, 2001) + [(1+hanning(2001))'.*sin(2.*pi.*(tim).*53) ; (1+hanning(2001))'.*sin(2.*pi.*(tim).*79 - 0.025) ; (1+hanning(2001))'.*sin(2.*pi.*(tim).*127 + 0.002)];
filt = ft_preproc_dftfilter(dat, 1000, [53, 79, 127], 'dftreplace', 'neighbour', 'dftneighbourwidth', [1 1 1]);
figure;
subplot(2,2,1); plot(tim, dat(1,:)-filt(1,:)); hold on;plot(tim, (1+hanning(2001))'.*sin(2.*pi.*(tim).*53));
subplot(2,2,2); plot(tim, dat(2,:)-filt(2,:)); hold on;plot(tim, (1+hanning(2001))'.*sin(2.*pi.*(tim).*79 - 0.025));
subplot(2,2,3); plot(tim, dat(3,:)-filt(3,:)); hold on;plot(tim, (1+hanning(2001))'.*sin(2.*pi.*(tim).*127 + 0.002));

tim = (0:5000)./678.253;
krn = (1+hanning(5000)')./2;
dat = randn(1, 5001) + ([krn(1:2500) ones(1,2501)]).*sin(2.*pi.*(tim).*50);
filt = ft_preproc_dftfilter(dat, 678.253, 50, 'dftreplace', 'neighbour', 'dftneighbourwidth', 2, 'dftbandwidth', 2);
figure;hold on;plot(tim, ([krn(1:2500) ones(1,2501)]).*sin(2.*pi.*(tim).*50));plot(tim, dat-filt)

tim = (0:4999)./1000;
krn = (1+hanning(5000)')./2;
dat = randn(1, 5000) + ([krn(1:2500) ones(1,2500)]).*sin(2.*pi.*(tim).*50);
filt = ft_preproc_dftfilter(dat, 1000, 50, 'dftreplace', 'neighbour', 'dftneighbourwidth', 2, 'dftbandwidth', 2);
figure;hold on;plot(tim, ([krn(1:2500) ones(1,2500)]).*sin(2.*pi.*(tim).*50));plot(tim, dat-filt)

%%%%%%%%%
% Code chunk from Sabine
samples_off = 2; % samples to add, so a full 50 Hz cycle doesn't fit

lengthsec= 4;  % data length in seconds
fs = 1000; % sampling rate
acfreq = 50; % set the powerline frequency, e.g., 50 or 60 Hz

datlength = (lengthsec*fs);
dat = gausswin(datlength, round(datlength/100) )';
t= (0:(length(dat)-1))/fs;
noise = cos( 2 * pi * acfreq * t );

% spectrum interpolation
filt_spec = ft_preproc_dftfilter(dat + noise, fs, acfreq, 'dftreplace', 'neighbour', ...
  'dftbandwidth', 2, 'dftneighbourwidth', 2 );

% add samples to the data to introduce leakage
datlength2 = (lengthsec*fs) + samples_off;  % LEAK, add samples, so a full cycle doesn't fit
dat2 = gausswin(datlength2, round(datlength2/100) )';
t2 = (0:(length(dat2)-1))/fs;
noise2 = cos( 2 * pi * acfreq * t2 );

% comment out error message in the dftfilter function, before running this
filt_spec_leak = ft_preproc_dftfilter(dat2 + noise2, fs, acfreq, 'dftreplace', 'neighbour', ...
  'dftbandwidth', 2, 'dftneighbourwidth', 2 );

figure;
hold all
plot( t , dat + noise,'k')
ylabel('Amplitude (a.u.)')
xlabel('Time (s)')
ylim( [ -1.25 2.05 ] )
xlim( [ 0 4.005 ] )
legend('Gaussian with 50 Hz line noise')
title('Gaussian with added 50 Hz sinusoid of constant amplitude - 4.003 s')

figure;
hold all
plot( t , filt_spec, 'b')
plot( t2 , filt_spec_leak,'r')
plot(t, dat,'k')
xlabel('Time (s)')
ylabel('Amplitude (a.u.)')
legend('Spectrum Interpolation Full Cycle', 'Spectrum Interpolation Not Full Cycle', ...
  'Original Clean Gaussian')

ylim( [ -0.5 1.5] )
xlim( [ 0 4 ] )
title(['DFT neighbour function used for mixed signal (gaussian + line noise)  - Data length = ' ...
  num2str(datlength/fs) ' s and ' num2str(datlength2/fs) ' s']);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [parent_pop] = srn(parent_pop)
%   This procedure implements srn (Srinivas and Deb's) function.
%   The canonical srn function is defined as below --
%       f1 = 2.0 + ((x_1 - 2.0)^2.0) + ((x_2 - 1.0)^2.0)
%       f2 = (9.0 * x_1) - ((x_2 - 1.0) ^ 2.0)
%       s.t.
%           c_1(x) = 1.0 - ((x_1^2.0) + (x_2^2.0))/225.0
%           c_2(x) = (3.0 * x_2/10.0) - (x_1/10.0) - 1.0
%       where
%           -20.0 <= x_i <= 20.0    (i = 1,2)

global nreal ;
global ncon ;
global nobj ;

cindex = nreal + nobj + 1 : nreal + nobj + ncon ;

x = parent_pop(:,1:nreal);
f1 = 2.0 + ((x(:,1) - 2.0) .^ 2.0) + ((x(:,2) - 1.0) .^ 2.0);
f2 = (9.0 .* x(:,1)) - ((x(:,2) - 1.0) .^ 2.0);

c = parent_pop(:,cindex) ;
c(:,1) = 1.0 - (((x(:,1) .^ 2.0) + (x(:,2) .^ 2.0)) ./ 225.0);
c(:,2) = (3.0 .* (x(:,2) ./ 10.0)) - (x(:,1) ./ 10.0) - 1.0 ;

parent_pop(:, (nreal+1)) = f1 ;
parent_pop(:, (nreal+2)) = f2 ;
parent_pop(:, cindex) = c ;
end

author: chudur-budur


text:
function [Qx,Qy]=CubicBezier1(Px,Py,n)
% Input:
% Px contains x-coordinates of control points [Px0,Px1,Px2,Px3]
% Py contains y-coordinates of control points [Py0,Py1,Py2,Py3]
% n is number of intervals (uniform parameterization)

%Output:
% Qx contains parameteric evaluted x-value
% Qy contains parameteric evaluted y-value

% Equation of Bezier Curve, utilizes Horner's rule for efficient computation.
% Q(t)=(-P0 + 3*(P1-P2) + P3)*t^3 + 3*(P0-2*P1+P2)*t^2 + 3*(P1-P0)*t + Px0

Px0=Px(1);
Py0=Py(1);
Px1=Px(2);
Py1=Py(2);
Px2=Px(3);
Py2=Py(3);
Px3=Px(4);
Py3=Py(4);

cx3=-Px0 + 3*(Px1-Px2) + Px3;
cy3=-Py0 + 3*(Py1-Py2) + Py3;
cx2=3*(Px0-2*Px1+Px2);
cy2=3*(Py0-2*Py1+Py2);
cx1=3*(Px1-Px0);
cy1=3*(Py1-Py0);
cx0=Px0;
cy0=Py0;

dt=1/n;
Qx(1)=Px0; % Qx at t=0
Qy(1)=Py0; % Qy at t=0
for i=1:n
    t=i*dt;
    Qx(i+1)=((cx3*t+cx2)*t+cx1)*t + cx0;
    Qy(i+1)=((cy3*t+cy2)*t+cy1)*t + cy0;
end

% % % --------------------------------
% % % Author: Dr. Murtaza Khan
% % % Email : drkhanmurtaza@gmail.com
% % % --------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function dataNdConcat = concat(this, dataNdArray, concatDims, tolerance)
%Concatenate an array of MrDataNd along one specified dimension
%
%   Y = MrDataNd()
%   dataNdConcat = Y.concat(dataNdArray, concatDims, tolerance)
%
% This is a method of class MrDataNd. It is very similar to combine with
% the difference that non-singleton dimensions can be concatenated (e.g.,
% volume 1-10 and 11-20).
% Note: Multi-dimensional concatenation is supported,
% but it will leave parts of the array empty
%   - e.g., volume 1-10 of slice 1-9 concatenated with volume 11-20
%     of slice 10-20 would leave the "cross-terms" of the dim Matrix like
%     slice 1-9 of volume 11-10 undefined (we set them to 0).
%
% NOTE: The order of data in the dataNdArray will *not* define the
%       concatenation order. Ratherthe actual values of the samplingPoints
%       in dimInfo along the concat dimensions specify the position of the
%       data.
%
% IN
%   dataNdArray     cell(nDatasets,1) of MrDataNd to be concatenated
%                       OR
%                   single MrDataNd object. In this case, the input object
%                   and the calling object will be concatenated
%   concatDims       index or string (label) of dimension along which data
%                   should be concatenated
%   tolerance                   dimInfos are only combined, if their
%                               information is equal for all but the
%                               concatDims (because only one
%                               representation is retained for those,
%                               usually from the first of the dimInfos).
%                               However, sometimes numerical precision,
%                               e.g., rounding errors, preclude the
%                               combination. Then you can increase this
%                               tolerance;
%                               default: single precision (eps('single')
%                               ~1.2e-7)
% OUT
%
% EXAMPLE
%   concat
%
%   See also MrDataNd MrDataNd.combine MrDimInfo.combine

% Author:   Lars Kasper
% Created:  2019-04-15
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

if nargin < 4
    tolerance = eps('single');
end

if ~iscell(dataNdArray)
    dataNdArray = {this, dataNdArray};
end

nImages = numel(dataNdArray);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Split all images in array along specified dimensions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dataNdArraySplit = {};
for iImage = 1:nImages
    dataNdArraySplit = [dataNdArraySplit; ...
        reshape(dataNdArray{iImage}.split('splitDims', concatDims), ...
        [], 1)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Combine all image arrays and reconcatenate!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Note that not the order here will be important, but the actual value of
% the slice position in dimInfo

dataNdConcat = dataNdArraySplit{1}.combine(dataNdArraySplit, concatDims, ...
    tolerance);

author: translationalneuromodeling

repo: tapas


text:
function x = At_fft2( b, N, mask, P, m, n, bComplex )
% 2D IFFT
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

if( nargin < 4 ), P = 1:N; end

K = length(b);
fx = complex(zeros(m,n),zeros(m,n));
fx(mask) = 1/sqrt(2)*b(1:K/2) + 1i *1/sqrt(2)* b(K/2+1:K);
x = zeros(N,1);

if(bComplex)
    fct = sqrt(N/2); % complex input image
else
    fct = sqrt(N); % real input image
end

tmp = fct*fftshift(ifft2(ifftshift(fx)));

if(~bComplex)
    % real input image
    x(P) = abs(tmp(:));
else
    % complex input image
    x(P) = [1/sqrt(2)*real(tmp(:)); 1/sqrt(2)*imag(tmp(:))];
end

% x(P) = real(tmp(:)); % phantom

end

author: thomaskuestner


text:
function [U,S,output] = lmlra_nls(T,U0,S0,options)
%LMLRA_NLS LMLRA by nonlinear least squares.
%   [U,S,output] = lmlra_nls(T,U0,S0) computes the factor matrices U{1},
%   ..., U{N} and core tensor S belonging to a low multilinear rank
%   approximation of the N-th order tensor T by minimizing
%   0.5*frob(T-lmlragen(U,S))^2. Each term U{r} is a cell array of N factor
%   matrices U{r}{n}, followed by a core tensor U{r}{N+1}. The algorithm is
%   initialized with the factor matrices U0{n} and core tensor S0. The
%   structure output returns additional information:
%
%      output.Name  - The name of the selected algorithm.
%      output.<...> - The output of the selected algorithm.
%
%   lmlra_nls(T,U0,S0,options) may be used to set the following options:
%
%      options.Algorithm =   - The desired optimization method.
%      [@nls_gncgs| ...
%       {@nls_gndl}|@nls_lm]
%      options.M =           - The preconditioner to use when
%      [{'block-Jacobi'}|...   options.LargeScale is true.
%       false]
%      options.<...>         - Parameters passed to the selected method,
%                              e.g., options.TolFun, options.TolX and
%                              options.PlaneSearchOptions. See also help
%                              [options.Algorithm].
%
%   See also lmlra_minf.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.
%   [2] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables," SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.

% Forward problem to btd_nls as a one-term BTD.
if nargin < 4, options = struct; end
[U,output] = btd_nls(T,{[U0(:).',S0]},options);
S = U{1}{end};
U = U{1}(1:end-1);

author: andrewssobral

repo: lrslibrary


text:
function cd2EditedFile()
% cd2EditedFile: goes to the the file opened in editor
%
% references at:
% <<http://blogs.mathworks.com/community/2011/05/16/matlab-editor-api-examples/>>

% Copyright 2012, Clemens Ager
%%
a = fileparts(matlab.desktop.editor.getActiveFilename);
cd(a);


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio]...
    = olmar2(fid, data, varargins, opts)
% This program starts the OLMAR-2 algorithm
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ....
%           = olmar2_start(fid, data, varargins, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = olmar2_start(fid, data, {10, 0.5, 0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
epsilon = varargins{1}; % reversion parameter \epsilon
alpha = varargins{2};       % Alpha size
tc = varargins{3};      % transaction cost fee rate

% Run the OLMAR-2 algorithm
[cum_ret, cumprod_ret, daily_ret, daily_portfolio]...
    = olmar2_run(fid, data, epsilon, alpha, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
% PROPAGATION_KERNEL_GRID calculates the propagation kernel between data
% organized in grids.
%
% This function contains an implementation of the propagation kernel
% for grid graphs. Propagation kernels for general graphs are described in:
%
% Neumann, M., Patricia, N., Garnett, R., and Kersting, K. Efficient
% Graph Kernels by Randomization. (2012). Machine Learning And
% Knowledge Discovery in Databases: European Conference, (ECML/PKDD
% 2012), pp. 378-392.
%
% This implementation supports (approximately) preserving any of the
% following distances between the feature vectors:
%
% - \ell^1 or \ell^2 for arbitrary feature vectors
% - the total variation or Hellinger distance for
% distribution-valued features
%
% Depending on the chosen distance, the input features will be
% transformed appropriately.
%
% This implementation also supports arbitrary transformations to be
% used. See the transformations/ directory for example implementations
% (isotropic diffusion). Transformations must satisfy the following very
% general interface:
%
% A = transformation(A);
%
% Where the input is a cell array of N feature matrices of size (n x m x d).
% One feature matrix per grid graph. (n x m) are the grid dimensions and d
% is the feature dimension. The transformation function computes a new set
% of features given an old set.
%
% Usage:
%
% K = propagation_kernel_grid(A, transformation, num_iterations, varargin)
%
% Inputs:
%
% A: a cell array of N feature matrices of sizes (n x m x d) to hash. n and
% m can be different for each grid graph. If either the total variation or
% Hellinger distance is chosen, the feature matrices should sum to 1 in
% their 3rd dimension.
%
% transformation: a function handle to a feature transformation
% function satisfying the above-described API.
%
% num_iterations: the number of iterations to use for the kernel
% computation.
%
% Optional inputs (specified as name/value pairs):
%
% 'w': the bin width to use during the hashing
% computation (default: 1e-4)
%
% 'distance': a string indicating the distance to approximately
% preserve; the following values are supported:
% 'l1': \ell^1
% 'l2': \ell^2
% 'tv': total variation distance (equivalent to \ell^1)
% 'hellinger': Hellinger distance
%
% The input is not case sensitive. See
% calculate_hashes for more information.
% (default: 'l1').
%
% 'base_kernel': a function handle to the base kernel to use. The
% kernel will be called as:
%
% K = base_kernel(counts),
%
% where counts is a (N x k) matrix, where m is the number of graphs and k
% is the number of unique hashes during this step of the computation.
% counts(i, j) contains the number of times hash j occurs in graph i.
% The default base kernel is the linear one:
%
% @(counts) (counts * counts');
%
% Outputs:
%
% K: an (N x N) matrix containing the computed propagation kernel.
%
% See also CALCULATE_HASHES_GRID, LABEL_DIFFUSION_CONVOLUTION.
%
%
% Based on the implementation of propagation_kernel by
% (c) Roman Garnett, 2012--2014.
%
% Copyright (c) Marion Neumann, 2014.

function K = propagation_kernel_grid(A, transformation, ...
                                                num_iterations, varargin)
  verbose = false;

  % parse optional inputs
  options = inputParser;

  % which distance to use
  options.addOptional('distance', 'l1', ...
                      @(x) ismember(lower(x), {'l1', 'l2', 'tv', 'hellinger'}));
  % width for hashing
  options.addOptional('w', 1e-4, ...
                      @(x) (isscalar(x) && (x > 0)));
  % base kernel for counts vector, default is linear kernel
  options.addOptional('base_kernel', ...
                      @(counts) (counts * counts'), ...
                      @(x) (isa(x, 'function_handle')));
  options.addOptional('num_steps', 1, ...
                      @(x) (isscalar(x) && (x > 0)));   % do num_stps feature updates per kernel contribution

  options.parse(varargin{:});
  options = options.Results;

  num_graphs = size(A,1);
  K = zeros(num_graphs);    % initialize output
  iteration = 0;
  while (true)

    % hashing
    if (verbose); fprintf('...computing hashvalues\n'); end
    hash_labels = calculate_hashes_grid(A, options.distance, options.w);

    % GET set of existing hash labels
    sth = cellfun(@(im) (unique(im(:))), hash_labels,'UniformOutput',false);
    label_set = unique(cell2mat(sth));
    min_label = min(label_set);
    num_labels = numel(label_set);
    label_set = label_set-min_label+1;

    if num_labels > 200000
            fprintf('skipping kernel contribution at iteration %d (num hashlabels too large = %d)\n',...
                    iteration, num_labels);
    else

        % aggregate counts on graphs
        if (verbose); fprintf('...computing kernel contribution\n'); end

        add_contrib = true;
        counts = zeros(num_graphs,num_labels);
        for i=1:num_graphs
            hash_labels{i} = hash_labels{i}-min_label+1;
            try
                curr_counts = accumarray(hash_labels{i}(:), 1);
            catch
                add_contrib=false;
                warning('skipping kernel contribution at iteration %d (num hashlabels = %d)\n', iteration, num_labels);
                break
            end
            idx = ismember(label_set,find(curr_counts));
            counts(i,idx) = curr_counts(find(curr_counts));
        end

        % contribution specified by base kernel on count vectors
        if add_contrib
            K = K + options.base_kernel(counts);
        end
    end
    % avoid unnecessary transformation on last iteration
    if (iteration == num_iterations)
      break;
    end

    % apply transformation to features for next step
    if (verbose); fprintf('...label update\n'); end
    for step=1:options.num_steps
        A = transformation(A);
    end
    iteration = iteration + 1;
  end
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function varargout = test_smooth( f, N, DOM )
% TEST_SMOOTH Runs diagnostic checks on a TFOCS smooth function object.
% TEST_SMOOTH( F )
%   tests whether the function handle F works as a TFOCS smooth
%   function object.
%
%   Requirements: f = F(X) must return the value of the function at X.
%       [f,g] = F(X) must return the value, f, and the gradient, g.
%
%   For an example, see SMOOTH_QUAD.M
%
% TEST_SMOOTH( F, N )
%    specifies the size of the domain space.  If N is a scalar,
%    then the domain space is the set of  N x 1 vectors.
%    If N is a matrix of the form [n1, n2], the the domain
%    space is the set of n1 x n2 matrices.
%
% TEST_SMOOTH( ..., DOM )
%    specifies the domain of F.  DOM can be of the form [a,b]
%    which signifies the 1D set (a,b).
%
% OK = TEST_SMOOTH(...)
%    returns "true" if all tests are passed, and "false" otherwise.
%
% See also private/tfocs_smooth, smooth_huber

OK = false;
PLOT = true;

error(nargchk(1,3,nargin));

if ~isa(f,'function_handle')
    fail('TFOCS smooth function must be a FUNCTION HANDLE');
    return;
end
fprintf('== Testing the smooth function %s ==\n', func2str(f) );

FORCE_N = false;  % always require a vector input
if nargin < 2 || isempty(N),   N = 10;
else FORCE_N = true;
end
if nargin < 3 || isempty(DOM), DOM = [-1,1];    end

a = DOM(1); b = DOM(2);

x = (a+b)/2;
fprintf('Testing scalar inputs... \n');
try
    vi = f(x);
catch
    fail('TFOCS smooth function failed to return a function value');
    return;
end
if isinf(vi)
    fail('TFOCS smooth function tester: default domain is invalid. Please specify valid domain');
    return;
end
fprintf('\t\t\t\t...passed. \n');

% Now, also ask for the gradient
fprintf('Testing gradient output... \n');
try
    [vi,gi] = f(x);
catch
    fail('TFOCS smooth function failed to return a derivative value');
    return;
end
fprintf('\t\t\t\t...passed. \n');

% Now, try a vector
if isscalar(N)
    x = repmat(x,N,1);
else
    x = repmat(x,N(1),N(2));
    % if a > 0, assume we want a PSD matrix:
    if a >= 0
        x = ones(N(1),N(2)) + eye(N(1),N(2) );
    end
end
fprintf('Testing vector inputs... \n');
try
    [vi,gi] = f(x);
catch
    fail('TFOCS smooth function failed when supplied a vector input');
    return;
end
fprintf('\t\t\t\t...passed. \n');

% 1D example.  Does not require function to be vectorized
n   = 100;  % number of grid points
h   = (b-a)/n;
grid = (a+h/2):h:(b-h/2);

v    = zeros(size(grid));
g    = zeros(size(grid));
first = @(x) x(1);
for i = 1:length(grid)
    v(i) = f(grid(i));
    if isinf(v(i))
        g(i) = v(i);
    else
        [v(i),gi] = f(grid(i));
        g(i) = first(gi);
    end
end

if PLOT
    figure;
    clf;
    plot(grid,v,'.-');
    hold all
    plot(grid,g,'.-');
    legend('function','derivative');
%     title(func2str(f),'interpreter','none');
    line([a,b], 0*[1,1],'color','k' )
    if a < 0 && b > 0
        line( 0*[1,1], get(gca,'ylim'),'color','k' );
    end
end

OK = true;

if nargout > 0
    varargout{1} = OK;
end

disp('Test passed succesfully.');

function fail(str)
    disp('Test failed.  Reason:');
    disp(str);

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
%% Clear, close
clear; close; clc;

%% Data load
size.subject = 9;
size.trial = 288;
size.class = 4;
size.chan = 22;
size.sr = 250;

for i = 1:size.subject
    [train(i).s, train(i).h] = sload(char("C:\Users\KHJ-work\Documents\MATLAB\BCICIV_2a_gdf\A0" + int2str(i) + "T.gdf"));
    [eval(i).s, eval(i).h] = sload(char("C:\Users\KHJ-work\Documents\MATLAB\BCICIV_2a_gdf\A0" + int2str(i) + "E.gdf"));
    train(i).event_typ = find(train(i).h.EVENT.TYP >= 769 & train(i).h.EVENT.TYP <= 772 );
    eval(i).event_typ = find(eval(i).h.EVENT.TYP == 783);
    train(i).h.Classlabel = load(char("C:\Users\KHJ-work\Documents\MATLAB\BCICIV_2a_gdf\A0" + int2str(i) + "T.mat"));
    eval(i).h.Classlabel = load(char("C:\Users\KHJ-work\Documents\MATLAB\BCICIV_2a_gdf\A0" + int2str(i) + "E.mat"));
end

%% Make Segment
for i = 1:size.subject
    for j = 1:size.trial
        train(i).samp_seg(:, 1:size.chan, j) = train(i).s(train(i).h.TRIG(j) + 626:train(i).h.TRIG(j) + 1125, 1:size.chan);
        eval(i).samp_seg(:, 1:size.chan, j) = eval(i).s(eval(i).h.TRIG(j) + 626:eval(i).h.TRIG(j) + 1125, 1:size.chan);
        if j < size.trial
            train(i).ref_seg(:, 1:size.chan, j) = train(i).s(train(i).h.TRIG(j) + 1626:train(i).h.TRIG(j) + 2000, 1:size.chan);
            eval(i).ref_seg(:, 1:size.chan, j) = eval(i).s(eval(i).h.TRIG(j) + 1626:eval(i).h.TRIG(j) + 2000, 1:size.chan);
        end
    end
end

%% Remove NaN value trial
for i = 1:size.subject
    for j = 1:size.trial
        if find(isnan(train(i).samp_seg(:, :, j)))
            train(i).samp_seg(:, :, j) = zeros;
            train(i).samp_remain(j) = 0;
        else
            train(i).samp_remain(j) = 1;
        end

        if find(isnan(eval(i).samp_seg(:, :, j)))
            eval(i).samp_seg(:, :, j) = zeros;
            eval(i).samp_remain(j) = 0;
        else
            eval(i).samp_remain(j) = 1;
        end

        if j < size.trial
            if find(isnan(train(i).ref_seg(:, :, j)))
                train(i).ref_seg(:, :, j) = zeros;
                train(i).ref_remain(j) = 0;
            else
                train(i).ref_remain(j) = 1;
            end

            if find(isnan(eval(i).ref_seg(:, :, j)))
                eval(i).ref_seg(:, :, j) = zeros;
                eval(i).ref_remain(j) = 0;
            else
                eval(i).ref_remain(j) = 1;
            end
        end
    end
end

%% Bandpass filtering
for i = 1:size.subject
    train(i).samp_segF = reshape(bandpass(reshape(train(i).samp_seg, [500 * size.trial size.chan]), [8 30], size.sr), [500, size.chan, size.trial]);
    eval(i).samp_segF = reshape(bandpass(reshape(eval(i).samp_seg, [500 * size.trial size.chan]), [8 30], size.sr), [500, size.chan, size.trial]);
    train(i).ref_segF = reshape(bandpass(reshape(train(i).ref_seg, [375 * (size.trial - 1) size.chan]), [8 30], size.sr), [375, size.chan, (size.trial - 1)]);
    eval(i).ref_segF = reshape(bandpass(reshape(eval(i).ref_seg, [375 * (size.trial - 1) size.chan]), [8 30], size.sr), [375, size.chan, (size.trial - 1)]);
end

%% Covariance matrix
for i = 1:size.subject
    for j = 1:size.trial
        train(i).samp_cov(:, :, j) = (1 / (length(find(train(i).samp_remain)) - 1)) * (train(i).samp_segF(:, :, j)' * train(i).samp_segF(:, :, j));
        eval(i).samp_cov(:, :, j) = (1 / (length(find(eval(i).samp_remain)) - 1)) * (eval(i).samp_segF(:, :, j)' * eval(i).samp_segF(:, :, j));
        if j < size.trial
            train(i).ref_cov(:, :, j) = (1 / (length(find(train(i).ref_remain)) - 1)) * (train(i).ref_segF(:, :, j)' * train(i).ref_segF(:, :, j));
            eval(i).ref_cov(:, :, j) = (1 / (length(find(eval(i).ref_remain)) - 1)) * (eval(i).ref_segF(:, :, j)' * eval(i).ref_segF(:, :, j));
        end
    end
end

%% Riemannian mean of whole sample and refference
for i = 1:size.subject
    train_samp_tmp = train(i).samp_cov;
    train_samp_tmp(:, :, train(i).samp_remain == 0) = [];
    eval_samp_tmp = eval(i).samp_cov;
    eval_samp_tmp(:, :, eval(i).samp_remain == 0) = [];
    train(i).samp_Rmean = mean(train_samp_tmp, 3);
    eval(i).samp_Rmean = mean(eval_samp_tmp, 3);

    train_ref_tmp = train(i).ref_cov;
    train_ref_tmp(:, :, train(i).ref_remain == 0) = [];
    eval_ref_tmp = eval(i).ref_cov;
    eval_ref_tmp(:, :, eval(i).ref_remain == 0) = [];
    train(i).ref_Rmean = mean(train_ref_tmp, 3);
    eval(i).ref_Rmean = mean(eval_ref_tmp, 3);

    for j = 1:5
        Tsamp = tangentspace(train_samp_tmp, train(i).samp_Rmean);
        Tmean = mean(Tsamp, 3);
        train(i).samp_Rmean = untangentspace(Tmean, train(i).samp_Rmean);

        Tsamp = tangentspace(eval_samp_tmp, eval(i).samp_Rmean);
        Tmean = mean(Tsamp, 3);
        eval(i).samp_Rmean = untangentspace(Tmean, eval(i).samp_Rmean);

        Tsamp = tangentspace(train_samp_tmp, train(i).ref_Rmean);
        Tmean = mean(Tsamp, 3);
        train(i).ref_Rmean = untangentspace(Tmean, train(i).ref_Rmean);

        Tsamp = tangentspace(eval_samp_tmp, eval(i).ref_Rmean);
        Tmean = mean(Tsamp, 3);
        eval(i).ref_Rmean = untangentspace(Tmean, eval(i).ref_Rmean);
    end
end

%% Riemannian mean of each class
for i = 1:size.subject
    for j = 1:size.class
        train_samp_tmp = train(i).samp_cov;
        train_samp_tmp = train_samp_tmp(:, :, train(i).h.Classlabel.classlabel == j & train(i).samp_remain');
        eval_samp_tmp = eval(i).samp_cov;
        eval_samp_tmp = eval_samp_tmp(:, :, eval(i).h.Classlabel.classlabel == j & eval(i).samp_remain');

        train(i).class_Rmean(:, :, j) = mean(train_samp_tmp, 3);
        eval(i).class_Rmean(:, :, j) = mean(eval_samp_tmp, 3);

        for k = 1:5
            Tsamp = tangentspace(train_samp_tmp, train(i).class_Rmean(:, :, j));
            Tmean = mean(Tsamp, 3);
            train(i).class_Rmean(:, :, j) = untangentspace(Tmean, train(i).class_Rmean(:, :, j));

            Tsamp = tangentspace(eval_samp_tmp, eval(i).class_Rmean(:, :, j));
            Tmean = mean(Tsamp, 3);
            eval(i).class_Rmean(:, :, j) = untangentspace(Tmean, eval(i).class_Rmean(:, :, j));
        end
    end
end

%% Rmdm classification
for i = 1:size.subject
    tmp_res = Rmdm(eval(i).samp_cov(:, :, find(eval(i).samp_remain)), train(i).class_Rmean);
    acc1(i) = length(find(tmp_res == eval(i).h.Classlabel.classlabel(find(eval(i).samp_remain))')) / length(eval(i).samp_remain);
    tmp_res = Rmdm(train(i).samp_cov(:, :, find(train(i).samp_remain)), eval(i).class_Rmean);
    acc2(i) = length(find(tmp_res == train(i).h.Classlabel.classlabel(find(train(i).samp_remain))')) / length(train(i).samp_remain);
    def_acc(i) = (acc1(i) + acc2(i)) / 2;
end

%% Covariance matrix using affine transform
for i = 1:size.subject
    for j = 1:size.trial
        train(i).aff_cov(:, :, j) = (train(i).ref_Rmean ^ (-1 / 2)) * train(i).samp_cov(:, :, j) * ((train(i).ref_Rmean) ^ (-1 / 2));
        eval(i).aff_cov(:, :, j) = (eval(i).ref_Rmean ^ (-1 / 2)) * eval(i).samp_cov(:, :, j) * ((eval(i).ref_Rmean) ^ (-1 / 2));
    end
end

%% Riemannian mean of each class of affine transformed samples
for i = 1:size.subject
    for j = 1:size.class
        train_aff_tmp = train(i).aff_cov;
        train_aff_tmp = train_aff_tmp(:, :, train(i).h.Classlabel.classlabel == j & train(i).samp_remain');
        eval_aff_tmp = eval(i).aff_cov;
        eval_aff_tmp = eval_aff_tmp(:, :, eval(i).h.Classlabel.classlabel == j & eval(i).samp_remain');

        train(i).aff_class_Rmean(:, :, j) = mean(train_aff_tmp, 3);
        eval(i).aff_class_Rmean(:, :, j) = mean(eval_aff_tmp, 3);

        for k = 1:5
            Tsamp = tangentspace(train_aff_tmp, train(i).aff_class_Rmean(:, :, j));
            Tmean = mean(Tsamp, 3);
            train(i).aff_class_Rmean(:, :, j) = untangentspace(Tmean, train(i).aff_class_Rmean(:, :, j));

            Tsamp = tangentspace(eval_aff_tmp, eval(i).aff_class_Rmean(:, :, j));
            Tmean = mean(Tsamp, 3);
            eval(i).aff_class_Rmean(:, :, j) = untangentspace(Tmean, eval(i).aff_class_Rmean(:, :, j));
        end
    end
end

%% Rmdm classification of affine transformed samples
for i = 1:size.subject
    tmp_res = Rmdm(eval(i).aff_cov(:, :, find(eval(i).samp_remain)), train(i).aff_class_Rmean);
    acc1(i) = length(find(tmp_res == eval(i).h.Classlabel.classlabel(find(eval(i).samp_remain))')) / length(eval(i).samp_remain);
    tmp_res = Rmdm(train(i).aff_cov(:, :, find(train(i).samp_remain)), eval(i).aff_class_Rmean);
    acc2(i) = length(find(tmp_res == train(i).h.Classlabel.classlabel(find(train(i).samp_remain))')) / length(train(i).samp_remain);
    aff_acc(i) = (acc1(i) + acc2(i)) / 2;
end

%% MINE
for i = 1:size.subject
    for j = 1:size.subject
        for k = 1:size.class
            for l = find(train(j).h.Classlabel.classlabel' == k & train(j).samp_remain)
                tmp_samp(:, :, l) = (train(i).aff_class_Rmean(:, :, k) ^ (-1 / 2)) * train(j).aff_cov(:, :, l) * (train(i).aff_class_Rmean(:, :, k) ^ (-1 / 2));
            end
        end

        if i == 1 && j == 1
            train(i).transfer_sample = tmp_samp;
            continue
        end

        train(i).transfer_sample = cat(3, train(i).transfer_sample, tmp_samp);
    end
end

%%
for i = 1:size.subject
    if i == 1
        train_remain = train(i).samp_remain;
        train_class = train(i).h.Classlabel.classlabel;

    else
        train_remain = cat(2, train_remain, train(i).samp_remain);
        train_class = cat(1, train_class, train(i).h.Classlabel.classlabel);
    end
end

%%
for i = 1:size.subject
    for j = 1:size.class
        train_mine_tmp = train(i).transfer_sample;
        train_mine_tmp = train_mine_tmp(:, :, train_class == j & train_remain');

        train(i).mine_class_Rmean(:, :, j) = mean(train_mine_tmp, 3);

        for k = 1:5
            Tsamp = tangentspace(train_mine_tmp, train(i).mine_class_Rmean(:, :, j));
            Tmean = mean(Tsamp, 3);
            train(i).mine_class_Rmean(:, :, j) = untangentspace(Tmean, train(i).mine_class_Rmean(:, :, j));
        end
    end
end

%%
for i = 1:size.subject
    for j = 1:size.class
        for k = find(eval(j).h.Classlabel.classlabel' == j & eval(j).samp_remain)
            tmp_samp(:, :, k) = (train(i).aff_class_Rmean(:, :, j) ^ (-1 / 2)) * eval(i).aff_cov(:, :, k) * (train(i).aff_class_Rmean(:, :, j) ^ (-1 / 2));
        end
    end

    tmp_res = Rmdm(eval(i).aff_cov(:, :, find(eval(i).samp_remain)), train(i).mine_class_Rmean);
    acc1(i) = length(find(tmp_res == eval(i).h.Classlabel.classlabel(find(eval(i).samp_remain))')) / length(eval(i).samp_remain);
end

author: PatternRecognition

repo: OpenBMI


text:
classdef OnlineNaiveBayes < handle

  properties

    means;
    s_temp;
    std2s;
    counts;
    forward_bias;
    first;
    normpdf = @(x, mu, sigma)exp(-0.5 * ((x - mu)./sigma).^2) ./ (sqrt(2*pi) .* sigma);

  end

  methods

    function a = OnlineNaiveBayes(forward_bias)
      %Naive Bayes with a forward bias for online learning.
      %forward_bias = 0 (default) is standard naive bayes.

      if (nargin < 1)
        forward_bias = 0;
      end
      a.forward_bias = forward_bias;
      a.counts = ones(1, 3) * 2;
      a.first = true;

    end

    function train(a, inputs, outputs)

      arrayfun(@(x)a.online_train(inputs(x, :), outputs(x)), 1:size(inputs, 1));

    end

    function online_train(a, input, output)

      if (a.first)
        features = size(input, 2);
        a.means = zeros(3, features);
        a.s_temp = ones(3, features);
        a.std2s = ones(3, features);
        a.first = false;
      end
      output = (output + 3) / 2;
      indices = [output; 3];
      input = repmat(input, 2, 1);
      a.counts(indices) = a.counts(indices) + 1;
      r_counts = repmat(a.counts(indices), size(input, 2), 1)';
      new_means = a.means(indices, :) + ...
        (1 + (a.forward_bias * (a.counts(3) > 5))) .* ...
        (input - a.means(indices, :)) ./ r_counts;
      a.s_temp(indices, :) = a.s_temp(indices, :) + abs((input - a.means(indices, :)) .* (input - new_means));
      a.std2s(indices, :) = sqrt(a.s_temp(indices, :) ./ (r_counts - 1));
      a.means(indices, :) = new_means;

    end

    function outputs = test(a, inputs)

      probs = cell2mat(arrayfun(@(x)prod(a.normpdf(inputs, ...
        repmat(a.means(x, :), size(inputs, 1), 1), ...
        repmat(a.std2s(x, :), size(inputs, 1), 1)), 2) ...
        * (a.counts(x) / a.counts(3)), 1:2, 'UniformOutput', false));
      [v i] = max(probs(:, 1:2), [], 2);
      outputs = (i * 2) - 3;

    end

    function margins = margins(a, inputs)

      probs = cell2mat(arrayfun(@(x)prod(a.normpdf(inputs, ...
        repmat(a.means(x, :), size(inputs, 1), 1), ...
        repmat(a.std2s(x, :), size(inputs, 1), 1)), 2) ...
        * (a.counts(x) / a.counts(3)), 1:2, 'UniformOutput', false));
      [v i] = max(probs(:, 1:2), [], 2);
      margins = ((i * 2) - 3) .* (v ./ sum(probs(:, 1:2), 2));

    end

    function inputs = generate(a, outputs)

      index = (outputs == -1) + 1;
      mean = a.means(index, :);
      std = a.std2s(index, :);
      inputs = normrnd(mean, std);

    end

  end

end

author: epfl-lasa


text:
function [fPerMc, res2] = calc_GoodnessFit(mCatalog, fMc, fBinning)
% function [fPerMc, res2] = calc_FitComp(mCatalog, fMc, fBinning);
% -----------------------------------------------------------------------------
%
% Function to calculate  Goodness of fit percentage of a fitted to observed
% magnitude frequency distribution
%
% Incoming variables:
% mCatalog
% fMc
% fBinning
%
% Outgoing variables:
% fPerMc : Percentage of fit
% res2   : Result
%
% Author: J. Woessner
% woessner@seismo.ifg.ethz.ch
% last update: 23.01.03

vSel = mCatalog(:,6) > (fMc - (fBinning/2));
nNumberEvents = length(mCatalog(vSel,6));
if nNumberEvents >= 25
    [fDummy fBValue fDummy,  fDummy] =  bmemag(mCatalog(vSel,:));

    fStartMag = fMc; % Starting magnitude (hypothetical Mc)

    % log10(N)=A-B*M
    vMag = [fStartMag:fBinning:10]; % Ending magnitude must be sufficiently high
    vNumber = 10.^(log10(nNumberEvents)-fBValue*(vMag - fStartMag));
    vNumber = round(vNumber);

    % Find the last bin with an event
    nLastEventBin = min(find(vNumber == 0)) - 1;
    if isempty(nLastEventBin)
        nLastEventBin = length(vNumber);
    end
    % Determine set of all magnitude bins with number of events > 0
    ct = round((vMag(nLastEventBin)-fStartMag)*(1/fBinning) + 1);

    PM=vMag(1:ct);
    vNumber = vNumber(1:ct);
    [bval, vDummy] = hist(mCatalog(vSel,6),PM);
    b3 = fliplr(cumsum(fliplr(bval)));    % N for M >= (counted backwards)
    res2 = sum(abs(b3 - vNumber))/sum(b3)*100;
    fPerMc = 100-res2;
end

author: CelsoReyes


text:
function point = pointOnLine(line, pos)
%POINTONLINE Create a point on a line at a given position on the line
%
%   P = pointOnLine(LINE, POS);
%   Creates the point belonging to the line LINE, and located at the
%   distance D from the line origin.
%   LINE has the form [x0 y0 dx dy].
%   LINE and D should have the same number N of rows. The result will have
%   N rows ans 2 column (x and y positions).
%
%   See also:
%   lines2d, points2d, onLine, onLine, linePosition
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 07/04/2004.
%

angle = lineAngle(line);
point = [line(:,1) + pos .* cos(angle), line(:,2) + pos .* sin(angle)];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function n = p34_n ( )

%*****************************************************************************80
%
%% P34_N returns the number of variables for problem 34.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 2;

  return
end


text:
function signal = flt_delayembed(varargin)
% Apply delay embedding to epoched data.
% Signal = flt_delayembed(Signal,NumLags)
%
% Delay embedding is essentially appending to each multi-channel samples the subsequent k
% multi-channel samples, and thereby multiplies the number of channels by k. Delay embedding is a
% practical tool to extend linear spatial models (e.g., spatial filters or independent components)
% to linear spatio-temporal (and therefore implicitly spatio-spectral) models, just by applying
% those models to delay-embedded data. As a result, approaches that can learn optimal spatial
% filters (finding sources of interest) can be repurposed to learning optimal spatio-spectral
% filters (jointly finding sources and frequencies of interest).
%
% The tradeoff associated with delay-embedding is that the complexity of the models increases and
% they become harder to estimate, which might require more data or better
% constraints, priors, or regularization.
%
% In:
%   Signal : Epoched data set to be processed
%
%   NumLags : the number of lags that shall be used for delay-embedding (default: 1)
%
%   IncludeIntermediates : Include intermediate lags. If this is set to false, only the 0''th and
%                          the N''th lag will be embedded. (default: true)
%
% Out:
%   Signal : the processed signal; will have more channels
%
% Notes:
%   The temporal filters that can be designed for a small number of lags are often limited to
%   high-frequency responses; the frequency range can often be extended without increasing model
%   complexity by first resampling the data to the lowest acceptable sampling rate (e.g., 60 Hz).
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-11-17

% flt_delayembed_version<1.00> -- for the cache

if ~exp_beginfun('filter') return; end

% requires epoched data, works best on spatially filtered data
declare_properties('name','DelayEmbedding', 'depends','set_makepos', 'follows',{'flt_project','flt_window'}, 'independent_channels',false, 'independent_trials',true);

% declare arguments
arg_define(varargin,...
    arg_norep({'signal','Signal'}), ...
    arg({'numlags','NumLags'}, 1, uint32([1 1 20 1000]), 'Number of lags. For delay-embedding.'), ...
    arg({'includeIntermediates','IncludeIntermediates'}, true, [], 'Include intermediate lags. If this is set to false, only the 0''th and the N''th lag will be embedded.'));

for k=quickif(includeIntermediates,0:numlags,[0 numlags])
    tmp{k+1} = signal.data(:,k+(1:end-numlags),:); end
signal.data = cat(1,tmp{:});
[signal.nbchan,signal.pnts,signal.trials] = size(signal.data);
signal.chanlocs = struct('labels', cellfun(@num2str,num2cell(1:signal.nbchan),'UniformOutput',false));

exp_endfun;

repo: REST


text:
function varargout=sde_correlate(c,r1)
%SDE_CORRELATE  Correlated values.
%   R2 = SDE_CORRELATE(C, R1) returns the matrix R2 of M correlated values given
%   the N-by-N correlation matrix C and the M-by-N matrix of (uncorrelated)
%   values R1. The first column of R2 is always SQRT(C(1,1)) times the first
%   column of R1.
%
%   [R2, D] = SDE_CORRELATE(C, R1) also returns the N-by-N diffusion matrix, D,
%   created from the correlation matrix, C.
%
%   D = SDE_CORRELATE(C) without a second input argument returns only the N-by-N
%   diffusion matrix, D, created from the correlation matrix, C.
%
%   Note:
%       C may either be a symmetric positive semidefinite correlation matrix or
%       covariance matrix.
%
%   Example:
%       % Plot uncorrelated and correlated normally-distributed points
%       r1 = randn(1e4,2); corr(r1)
%       r2 = sde_correlate([1 -0.8;-0.8 1],r1); corr(r2)
%       figure; plot(r1(:,1),r1(:,2),'b.',r2(:,1),r2(:,2),'r*'); axis equal;
%
%   See also: SDE_DECORRELATE, RAND, RANDN, RANDSTREAM, RANDSTREAM/RANDN

%   Andrew D. Horchler, horchler @ gmail . com, Created 5-20-13
%   Revision: 1.2, 7-17-13

if nargout > min(nargin,2)
    error('SDETools:sde_correlate:TooManyOutputs',...
          'Too many output arguments for number of supplied inputs.');
end

if isempty(c) || isempty(r1)
    if nargin == 1
        varargout{1} = [];
    else
        varargout{2} = [];
    end
else
    [m,n] = size(c);
    if ndims(c) ~= 2 || m ~= n              %#ok<ISMAT>
        error('SDETools:sde_correlate:NonSquareMatrix',...
              'The correlation matrix must be square.');
    end

    if sde_isdiag(c)
        num = rank(c);
        c = diag(c);
        if any(c) < 0
            error('SDETools:sde_correlate:NegativeDiagonal',...
                 ['All elements of the diagonal of the correlation matrix '...
                  'must be non-negative.']);
        end

        d = sqrt(c);
        isDiag = isscalar(c);
    else
        [d,num] = cholcov(c);
        if isempty(d)
            error('SDETools:sde_correlate:NonSymmetricSemiDefinite',...
                 ['The correlation matrix must be a symmetric positive '...
                  'semidefinite matrix.']);
        end
        isDiag = false;
    end

    if nargin == 1
        if num ~= m
            warning('SDETools:sde_correlate:NonFullRank1',...
                    'The correlation matrix is not full rank.');
        end

        if isDiag
            varargout{1} = diag(d);
        else
            varargout{1} = d;
        end
    else
        if ndims(r1) ~= 2 || size(r1,2) ~= m	%#ok<ISMAT>
            error('SDETools:sde_correlate:DimensionMismatch',...
                 ['The number of columns in the matrix of normally '...
                  'distributed values must equal the dimension of the '...
                  'correlation matrix.']);
        end
        if num ~= m
            warning('SDETools:sde_correlate:NonFullRank2',...
                   ['The correlation matrix is not full rank. Only the '...
                    'first %d columns of the input values will be '...
                    'correlated.'],num);
        end

        if isDiag
            varargout{1} = bsxfun(@mtimes,r1,d);
            if nargout == 2
                varargout{2} = diag(d);
            end
        else
            varargout{1} = r1*d;
            if nargout == 2
                varargout{2} = d;
            end
        end
    end
end

author: horchler

repo: SDETools

save_path: github-repos/MATLAB/horchler-SDETools