text:
function [] = aps_wrf_SAR()
% [] = aps_wrf_SAR()
% Script to load wrf processed data, on 37 pressure levels. Load data for area of interest,
% and calculate delays maps for each specified date.
% The DEM file inputed should have an asociated
% ".rsc" file, with the same filename as the DEM. The ".rsc" files should
% contain a WIDTH, LENGTH, X_FIRST, Y_FIRST, X_STEP, Y_STEP and optional a
% FORMAT string. The meris data is assumed to be structured in date folders.
% The batchfile contains the full path to the meris files in these folders.
% Note that the first line of the batchfile should read "files".
%
%
% INPUTS:
% demfile               Full path to the DEM file. The DEM needs to me in meters.
% xlims                 Limits in the x-direction, either in degrees
% ylims                 Limits in the y-direction, either in degrees
% demnull               The value for no DEM data, default is -32768.
% smpres                The output resolution, either in degrees
%                       Units needs to be consistend with xlims and ylims.
%
% OUTPUTS:
% It will give the computed ZENITH dry and wet delay map in cm for the selected region.
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% Richard Walters/ Bekaert David - Oxford / Leeds University 2012-2013
% modifications:
% 19/11/2013    DB      Change the computation of the saturated vapor
%                       pressure over ice and water
% 19/11/2013    DB      Change the computation of the hydrostatic delay to
%                       be directly computed from pressure and temperature
% 20/11/2013    DB      Check if wrf data exist to avoid crashing
% 21/11/2013    DB      Interpolate H,Temp, Hum spatially to cope with
%                       coarser land mask compared to the DEM
% 02/01/2014	DB      fix for NaN value over sea
% 04/08/2014    DB      Redefine meris_lat(lon)_range to region_lat(lon)_range
% 24/02/2015    DB      Save support information for WRF
% 27/09/2015    DB      Separate the DEM loading, ask for number of domains.
% 28/09/2015    DB      Also save the delay computation for a grid or z to include Pirate compartibility.
% 28/09/2015    DB      Include multi-core from matlab
% 07/07/2016    DB      Redefine hydrostatic delay to be based on surface pressure.

save_complete=0;        % save support information when 0

%% Constants
% Parameters from Jolivet et al 2011, GRL
Rd = 287.05;            % [j/kg/K] specific gas constants dry air
Rv = 461.495;           % [j/kg/K] water vapour
k1 = 0.776;             % [K/Pa]
k2 = 0.716;             % [K/Pa]
k3 = 3.75e3;            % [K^2/Pa]

%% Defaults
zref = 15000;       % zref for integration- tropopause
zincr = 15;         % z increment spacing for vertical interpolation before integration
vertres = 5;        % vertical resolution for comparing dem to vert profiles of delay

%% get the number of domains
repeat=1;
while repeat==1
    action_flag= str2num(input('How many domains do you have including the parent? ','s'));
    if isnumeric(action_flag)
        n_domains = action_flag;
        repeat=0;
    end
end

%% getting the variables from the parms_aps file
stamps_processed = getparm_aps('stamps_processed',1);
ll_matfile = getparm_aps('ll_matfile',1);
wrf_datapath = getparm_aps('wrf_datapath',1);
datestructure = 'yyyymmdd';                               % assumed date structure for ERA-I
UTC_sat =  getparm_aps('UTC_sat',1);

% loading the data
if strcmp(stamps_processed,'y')
   ps = load(ll_matfile);
   dates = ps.day;
   load psver
   fprintf('Stamps processed structure \n')
else
    psver = 2;
    ifgday_matfile = getparm_aps('ifgday_matfile',1);
    ifgs_dates = load(ifgday_matfile);
    ifgs_dates = ifgs_dates.ifgday;
    dates = reshape(ifgs_dates,[],1);
    dates = unique(dates);
    dates = datenum(num2str(dates),'yyyymmdd');

end

% getting the dates
n_dates = length(dates);

%% Compute and resample DEM
[dem,xmin,xmax,ymin,ymax,smpres,nncols,nnrows] = get_DEM;

% the region which is cropped from the ERA data and used to make the interpolation.
% Should be  larger than the region to which the delay is computed
lonmin = floor(xmin)-1;
lonmax= ceil(xmax)+1;
latmin = floor(ymin)-1;
latmax = ceil(ymax)+1;

% setting the maximum height of the DEM to limit the range at whcih ERA-!
% needs to be interpolated to
maxdem = ceil(max(max(dem))/100)*100+50; % max height of dem in metres

fprintf(['Interpolate to a maximum dem height of ', num2str(maxdem) ,' m\n'])

%% performing the calucluation for each  date

for d = 1:n_dates
    date = datestr(dates(d),datestructure);
    date_str_file = datestr(dates(d),'yyyy-mm-dd');

    % the netcdf files
    infile = [wrf_datapath filesep date filesep 'wrfplev_d0' num2str(n_domains) '_' date_str_file '_' UTC_sat ':00'];
    infile2 = [wrf_datapath filesep date filesep 'wrfout_d0' num2str(n_domains) '_' date_str_file '_' UTC_sat ':00'];

    if exist(infile,'file')==2
        % the output file names
        outfile = [wrf_datapath filesep date filesep date '_ZWD.xyz'];
        hydroutfile = [wrf_datapath filesep date filesep date '_ZHD.xyz'];
        outfile2 = [wrf_datapath filesep date filesep date '_ZWD.z'];
        hydroutfile2 = [wrf_datapath filesep date filesep date '_ZHD.z'];

        % getting the information from the output file
        ncid = netcdf.open(infile,'NC_NOWRITE');
        [numdims,numvars,numglobalatts,unlimdimid] = netcdf.inq(ncid);
        [dimname, dimlen] = netcdf.inqDim(ncid,0); % dimensions
        % 3 = pressure*37
        % 6 = temperature

        % variables at each node 20 (0-19)
        [varname, vartype, dimids, natts] = netcdf.inqVar(ncid,0);

        % load variables
        Temp = netcdf.getVar(ncid,6);       % temp in K
        Hum = netcdf.getVar(ncid,7);        % relative humidity in percent
        H = netcdf.getVar(ncid,8);          % Geopotential Height in m
        Plevs = netcdf.getVar(ncid,3);      % 37 pressure levels in Pa
        netcdf.close(ncid)

        % flip the N-S axis
        Temp = flipdim(Temp,2);
        Temp = permute(Temp,[2,1,3]);
        Hum = flipdim(Hum,2);               % flip the N-S axis
        Hum = permute(Hum,[2,1,3]);
        H = flipdim(H,2);                   % flip the N-S axis
        H = permute(H,[2,1,3]);

        % replacing the non-usefull values by NaN
        Temp(Temp==-999)=NaN;
        Hum(Hum==-999)=NaN;
        H(H==-999)=NaN;

        % getting the longitude and latitude coordinates
        numlat = size(Temp,1);
        numlon = size(Temp,2);
        ncid = netcdf.open(infile2,'NC_NOWRITE');
        lats = netcdf.getVar(ncid,1);
        latgrid = repmat(lats,[1,1,37]);
        latgrid = permute(latgrid,[2,1,3]);
        latgrid = flipdim(latgrid,1);
        lons = netcdf.getVar(ncid,2);
        longrid = repmat(lons,[1,1,37]);
        longrid = permute(longrid,[2,1,3]);
        longrid = flipdim(longrid,1);
        % close netcdf file
        netcdf.close(ncid)

        % Load pressure levels and replicate vector to make 3d matrix same size as variables
        Plevs = Plevs./100; % convert into hPa (or millibars) 1000 hPa at surface
        Pressure = repmat(Plevs,[1,numlat,numlon]);
        Pressure = permute(Pressure,[2,3,1]);

        % Get list of points to look at in analysis
        [xx,yy] = meshgrid(1:numlon,1:numlat);
        latlist = reshape(latgrid(:,:,1),[],1);
        lonlist = reshape(longrid(:,:,1),[],1);
        xlist = reshape(xx,[],1);
        ylist = reshape(yy,[],1);

        % When the landmask is coarser it can introduce extreme values at
        % sharp transitions. To cope with this I perform an interpoaltion
        % over the land points, such the land mask is removed.
        % interpolate the land
        n_grid_points = size(Pressure,1)*size(Pressure,2);
        ix_vert =find(sum(sum(isnan(Temp)))==n_grid_points);
        if isempty(ix_vert)
            ix_vert = size(Temp,3);
        end
        longrid_vector = reshape(longrid(:,:,1),[],1);
        latgrid_vector= reshape(latgrid(:,:,1),[],1);

        % check for land mask only those layers that have land within them.
        parfor kkk=1:ix_vert(1)-1
             ix = find(isnan(Temp(:,:,kkk))==1);
             if ~isempty(ix)
                 Temp_temp = reshape(Temp(:,:,kkk),[],1);
                 Hum_temp = reshape(Hum(:,:,kkk),[],1);
                 H_temp = reshape(H(:,:,kkk),[],1);

                 longrid_vector_temp = longrid_vector;
                 latgrid_vector_temp = latgrid_vector;

                 % the land points that need to be interpolated
                 longrid_vector_needed = longrid_vector(ix);
                 latgrid_vector_needed = latgrid_vector(ix);

                 % removing the land poitns with currently NaN values
                 Temp_temp(ix)=[];
                 Hum_temp(ix)=[];
                 H_temp(ix)=[];

                 longrid_vector_temp(ix)=[];
                 latgrid_vector_temp(ix)=[];

                 % do the interpolation for the landpoints
                 Temp_vector_needed = griddata(double(longrid_vector_temp),double(latgrid_vector_temp),double(Temp_temp),double(longrid_vector_needed),double(latgrid_vector_needed),'linear');
                 Hum_vector_needed = griddata(double(longrid_vector_temp),double(latgrid_vector_temp),double(Hum_temp),double(longrid_vector_needed),double(latgrid_vector_needed),'linear');
                 H_vector_needed = griddata(double(longrid_vector_temp),double(latgrid_vector_temp),double(H_temp),double(longrid_vector_needed),double(latgrid_vector_needed),'linear');
                 Temp_new = Temp(:,:,kkk);
                 Hum_new = Hum(:,:,kkk);
                 H_new = H(:,:,kkk);

                 % replacing the values
                 Temp_new(ix)=Temp_vector_needed;
                 Hum_new(ix)=Hum_vector_needed;
                 H_new(ix)=H_vector_needed;

                 % storing them back in the array
                 Temp(:,:,kkk) = Temp_new;
                 Hum(:,:,kkk) = Hum_new;
                 H(:,:,kkk) =  H_new;
             end
        end
        clear Temp_new Hum_new H_new ix Temp_temp Hum_temp H_temp longrid_vector_temp latgrid_vector_temp

        % select the grid points within the area of interest
        % the projection of the grid is not uniform for the longitude as the center is defined for a specific longitude coordinate
        latlist_temp = sort(abs(diff(latlist)));
        lat_res = abs(latlist_temp(ceil(length(latlist_temp)*0.95)))*1.5;

        lonlist_temp = diff(lonlist);
        lon_res = max(abs(lonlist_temp))*1.5;

        % make sure the grid selected is covered and can be interpolated by WRF. Add a resolution cell to make sure it is the case.
        ix = find(ymin-lat_res<= latlist & latlist <= ymax+lat_res  & xmin-lon_res<= lonlist & lonlist<= xmax+lon_res) ;
        xlist = xlist(ix);
        ylist = ylist(ix);
        latlist = latlist(ix);
        lonlist = lonlist(ix);
        latlist=double(latlist);
        lonlist=double(lonlist);

        numy = length(unique(latlist));
        numx = length(unique(lonlist));
        ulatlist = unique(latlist);
        ulonlist = unique(lonlist);
        uxlist = unique(xlist);
        uylist = unique(ylist);

        if save_complete==0
            % saving the information for support plotting
            wrf.wrf_lonlat = [lonlist latlist];
            wrf.region = [[xmin xmin xmax xmax xmin]' [ymin ymax ymax ymin ymin]'];
            deminfo.xmin = xmin;
            deminfo.xmax = xmax;
            deminfo.ymax = ymax;
            deminfo.ymin = ymin;
            deminfo.dem = dem;
            wrf.deminfo = deminfo;
            clear deminfo
            % checking if the file already exist. Yes append otherwiuze create it
            if exist('tca_support.mat','file')==2
              save('tca_support.mat','-append','wrf')
            else
              save('tca_support.mat','wrf')
            end
            save_complete=1;
        end

        % Could not find the wrf used equation as they appear to be mixed
        % with latent heat etc. Istead I used the equations used in ERA-I
        % (see IFS documentation part 2: Data assimilation (CY25R1)). Calculate
        % saturated water vapour pressure (svp) for water (svpw) using Buck 1881
        % and for ice (swpi) from Alduchow and Eskridge (1996) euation AERKi
        svpw = 6.1121.*exp((17.502.*(Temp-273.16))./(240.97+Temp-273.16));
        svpi = 6.1121.*exp((22.587.*(Temp-273.16))./(273.86+Temp-273.16));
        tempbound1 = 273.16; %0
        tempbound2 = 250.16; %-23

        svp = svpw;
        % Faster expression
        wgt = (Temp - tempbound2)/(tempbound1 - tempbound2);
        svp = svpi+ (svpw - svpi).*wgt.^2; % not sure why have wgt^2, Romain does it in his pyaps code.
        ix_bound1 = find(Temp > tempbound1);
        svp(ix_bound1) = svpw(ix_bound1);
        ix_bound2 = find(Temp < tempbound2);
        svp(ix_bound2) = svpi(ix_bound2);

        e = Hum./100.*svp;
        % NB, P = Pressure in hydrostatic equation, NOT dry partial pressure as
        % stated (incorrectly) in Jolivet et al. 2011.
        P = Pressure;

        %% Convert Geopotential to Geopotential Height and then to Geometric Height
        g0 = 9.80665;

        % map of g with latitude
        g = 9.80616.*(1 - 0.002637.*cosd(2.*latgrid) + 0.0000059.*(cosd(2.*latgrid)).^2);
        % map of Re with latitude
        Rmax = 6378137;
        Rmin = 6356752;
        Re = sqrt(1./(((cosd(latgrid).^2)./Rmax^2) + ((sind(latgrid).^2)./Rmin^2)));

        %% Calculate Geometric Height, Z
        Z = (H.*Re)./(g/g0.*Re - H);

        % Find middle of scene to work out glocal and Rlocal for use later
        midx = round(mean(uxlist));
        midy = round(mean(uylist));
        glocal = g(midy,midx,1);
        Rlocal = Re(midy,midx,1);

        %% Find middle of scene to work out glocal and Rlocal for use later
        cdslices = maxdem/vertres +1;
        cdstack = zeros(size(lonlist,1),cdslices);
        cdstack_dry = zeros(size(lonlist,1),cdslices);
        cdstack_wet = zeros(size(lonlist,1),cdslices);

        XI=(0:zincr:zref)';
        gh = glocal.*(Rlocal./(Rlocal+XI)).^2; %gravity with height for XI height range

        % Interpolate Temp P and e from 0:20:15000 m
        % then integrate using trapz to estimate delay as function of height
        parfor m = 1:size(lonlist,1);
            xn = xlist(m);
            yn = ylist(m);

            %interpolate at res zincr before doing integration
            X=double(squeeze(Z(yn,xn,:)));
            YP=double(squeeze(P(yn,xn,:)));
            Ye=double(squeeze(e(yn,xn,:)));
            YT=double(squeeze(Temp(yn,xn,:)));

            YT(isnan(Ye))=[];
            YP(isnan(Ye))=[];
            X(isnan(Ye))=[];
            Ye(isnan(Ye))=[];

            YeI = interp1(X,Ye,XI,'spline')*100; %change from hPa to Pa
            YPI = interp1(X,YP,XI,'spline')*100; %change from hPa to Pa
            YTI = interp1(X,YT,XI,'spline');

            tmp1 = ((k2-(Rd*k1/Rv)).*YeI./YTI + k3.*YeI./(YTI.^2));
            Lw = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp1)));
            % L is zenith delay one way in metres
            tmp2 = k1.*YPI./YTI;
            %Ld = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp2)));             % This is using P/T expression (Hanssen, 2001)
            gm = glocal;
            Ld = (10^-6).*((k1*Rd/gm).*(YPI - YPI(zref/zincr +1)));                 % This is P0 expression (Hanssen, 2001)

            % Interpolate important part (i.e. total delay at elevations
            % less than maxdem) at high res i.e. vertres, and put in cdstack.
            cdI=(0:vertres:maxdem)';
            LdI=interp1(XI,Ld,cdI,'spline');
            LwI=interp1(XI,Lw,cdI,'spline');

            %cdstack(j,i,:) = LsI;
            cdstack_dry(m,:) = LdI;
            cdstack_wet(m,:) = LwI;
        end

        % Interpolate each cdstack layer onto a grid given by the DEM extents
        % in UTM m.
        xsmpres = (xmax-xmin)/nncols;
        ysmpres = (ymax-ymin)/nnrows;
        [xi,yi] = meshgrid(xmin+0.5*xsmpres:xsmpres:xmax-0.5*xsmpres,ymin+0.5*ysmpres:ysmpres:ymax-0.5*ysmpres);

        cdstack_interp_dry = zeros(nnrows,nncols,cdslices);
        parfor n = 1:cdslices
           slicelist = reshape(cdstack_dry(:,n),[],1);
           newslice = griddata(lonlist,latlist,slicelist,xi,yi,'linear');
           %F = TriScatteredInterp(lonlist,latlist,slicelist,'linear');
           %newslicelist = F(xi,yi);
           %newslice = reshape(newslicelist,nnrows,nncols);
           cdstack_interp_dry(:,:,n)= flipud(newslice); % flipud due to ypositive downpage for matlab figs, and ypositive uppage for utmy
        end

        cdstack_interp_wet = zeros(nnrows,nncols,cdslices);
        parfor n = 1:cdslices
           slicelist = reshape(cdstack_wet(:,n),[],1);
           newslice = griddata(lonlist,latlist,slicelist,xi,yi,'linear');
           %F = TriScatteredInterp(lonlist,latlist,slicelist,'linear');
           %newslicelist = F(xi,yi);
           %newslice = reshape(newslicelist,nnrows,nncols);
           cdstack_interp_wet(:,:,n)= flipud(newslice); % flipud due to ypositive downpage for matlab figs, and ypositive uppage for utmy
        end
        % keeping the coordinates in the same grid as the data
        xi = flipud(xi);
        yi = flipud(yi);

        % Pull out delays from cdstack layers that match dem heights
        wetcorrection = ones(nnrows,nncols);
        hydrcorrection = ones(nnrows,nncols);
        rounddem = round(dem/vertres);
        rounddem(dem < 0)=0;

        % cope with the case that NaN are present. This can be sea-level
        rounddem(isnan(dem))=0;

        for i=1:nnrows
            for j=1:nncols
                wetcorrection(i,j) = cdstack_interp_wet(i,j,rounddem(i,j)+1);
            end
        end

        for i=1:nnrows
            for j=1:nncols
                hydrcorrection(i,j) = cdstack_interp_dry(i,j,rounddem(i,j)+1);
            end
        end

        % note that this is a one way Zenith delay and not a slant delay.
        % Units are in cm
        % Output wet correction
        wetcorrection = wetcorrection*100;                 % delay in [cm]
        fid = fopen(outfile,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(wetcorrection,[],1)]';
        tally = fwrite(fid,data_write,'double');
        fclose(fid);
        fid = fopen(outfile2,'w');
        fwrite(fid,wetcorrection','real*4');
        fclose(fid);

        % Output dry correction
        hydrcorrection = hydrcorrection*100;                 % delay in [cm]
        fid = fopen(hydroutfile,'w');
        data_write =  [reshape(xi,[],1) reshape(yi,[],1) reshape(hydrcorrection,[],1)]';
        tally = fwrite(fid,data_write,'double');
        fclose(fid);
        fid = fopen(hydroutfile2,'w');
        fwrite(fid,hydrcorrection','real*4');
        fclose(fid);

        fprintf([num2str(d) ' completed out of ' num2str(n_dates) '\n' ])

    else
        fprintf(['No wrf data for ' date ' \n'])
    end
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function [idcs, nSmpls] = class_samples( lbls, nActive )

% CLASS_SAMPLES
% COPYRIGHT: Patrick Sauer, 2012
% VARGPLVM

    if size(lbls,1) ~= 1
        lbls = lbls';
    end

    clids  = unique(lbls);
    nEx    = zeros(1, numel(clids));
    nSmpls = zeros(1, numel(clids));
    idcs   = cell( 1, numel(clids));

    count = 1;
    for c = clids
        tmp = find(lbls==c);
        idcs{count} = tmp;
        nEx(count)  = numel(tmp);
        nSmpTmp =  ceil(nEx(count)/length(lbls)*nActive);
        if nSmpTmp > nEx(count)
            nSmpTmp = nEx(count);
        end
        nSmpls(count) = nSmpTmp;

        count = count + 1;
    end

    if sum(nEx) < nActive
       error('There must be at least as many training samples as inducing variables.');
    end

    nRes = nActive - sum(nSmpls);
    if nRes > 0
        while nRes > 0
            [nSmplsSort, ind] = sort(nSmpls,'ascend');
            for count = ind
                if nSmpls(count) < nEx(count)
                    nSmpls(count) = nSmpls(count)+1;
                    nRes = nRes-1;
                end

                if nRes == 0
                    break;
                end
            end
        end
    else
       while nRes < 0
          [mx, ind] = max(nSmpls);
          nSmpls(ind) = nSmpls(ind) - 1;
          nRes = nRes + 1;
       end
    end
end

author: SheffieldML

repo: vargplvm


text:
EbN0=0:5:30;
MSE=zeros(1,length(EbN0));
for i=1:length(EbN0)
    SNR=EbN0(i);
    sim('v_blast_t4_r4_8PSK_crrect_MSE');
    MSE(i)=mse;
    save MSE;
end
semilogy(EbN0,MSE)
xlabel('EbN0(dB)')
ylabel('MSE')
grid on

author: Sable

repo: mcbench-benchmarks


text:
function blas1_c_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests CABS2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  CABS2 returns the L2 norm of a complex number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Real      Imaginary\n' );
  fprintf ( 1, '      Part      Part           CABS2(Z)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ c, seed ] = c4_uniform_01 ( seed );
    c = 5.0 * c;
    c_norm = cabs2 ( c );

    fprintf ( 1, '  %10f  %10f     %10f\n', real ( c ), imag ( c ), c_norm );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Sparse covariance estimation for Gaussian variables
% Joëlle Skaf - 04/24/08
% (a figure is generated)
%
% Suppose y \in\reals^n is a Gaussian random variable with zero mean and
% covariance matrix R = \Expect(yy^T), with sparse inverse S = R^{-1}
% (S_ij = 0 means that y_i and y_j are conditionally independent).
% We want to estimate the covariance matrix R based on N independent
% samples y1,...,yN drawn from the distribution, and using prior knowledge
% that S is sparse
% A good heuristic for estimating R is to solve the problem
%           maximize    logdet(S) - tr(SY) - lambda*sum(sum(abs(S)))
%           subject to  S >= 0
% where Y is the sample covariance of y1,...,yN, and lambda is a sparsity
% parameter to be chosen or tuned.
% A figure showing the sparsity (number of nonzeros) of S versus lambda
% is generated.

% Input data
randn('state',0);
n = 10;
N = 100;
Strue = sprandsym(n,0.5,0.01,1);
nnz_true = sum(Strue(:)>1e-4);
R = inv(full(Strue));
y_sample = sqrtm(R)*randn(n,N);
Y = cov(y_sample');
Nlambda = 20;
lambda = logspace(-2, 3, Nlambda);
nnz = zeros(1,Nlambda);

for i=1:Nlambda
    disp(['i = ' num2str(i) ', lambda(i) = ' num2str(lambda(i))]);
    % Maximum likelihood estimate of R^{-1}
    cvx_begin sdp quiet
        variable S(n,n) symmetric
        maximize log_det(S) - trace(S*Y) - lambda(i)*sum(sum(abs(S)))
        S >= 0
    cvx_end
    nnz(i) = sum(S(:)>1e-4);
end

figure;
semilogx(lambda, nnz);
hold on;
semilogx(lambda, nnz_true*ones(1,Nlambda),'r');
xlabel('\lambda');
legend('nonzeros in S', 'nonzeros in R^{-1}');

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function hedge = acqhedge_vbmc(action,hedge,stats,options)
%ACQPORTFOLIO Evaluate and update portfolio of acquisition functions.

switch lower(action(1:3))
    case 'acq'
        % Choose acquisition function based on hedge strategy

        if isempty(hedge)
            % Initialize hedge struct
            hedge.g = zeros(1,numel(options.SearchAcqFcn));
            hedge.n = numel(options.SearchAcqFcn);
            hedge.count = 0;
            hedge.lambda = 0.2;     % Lapse rate - random choice
            hedge.beta = 1;
            hedge.decay = options.AcqHedgeDecay^(options.FunEvalsPerIter);
        end

        hedge.count = hedge.count + 1;
        hedge.p = exp(hedge.beta*(hedge.g - max(hedge.g)))./sum(exp(hedge.beta*(hedge.g - max(hedge.g))));
        hedge.p = hedge.p*(1-hedge.lambda) + hedge.lambda/hedge.n;

        hedge.chosen = find(rand() < cumsum(hedge.p),1);
        hedge.phat = Inf(size(hedge.p));
        hedge.phat(hedge.chosen) = hedge.p(hedge.chosen);

    case 'upd'
        % Update value of hedge portfolio based on uncertainty reduction

        HedgeCutoff = 5;

        if ~isempty(hedge)
            iter = stats.iter(end);
            min_iter = max(1,iter-options.AcqHedgeIterWindow);

            min_sd = min(stats.elbo_sd(min_iter:iter-1));
            er_sd = max(0, log(min_sd / stats.elbo_sd(iter)));

            elcbo = stats.elbo - options.ELCBOImproWeight*stats.elbo_sd;
            max_elcbo = max(elcbo(min_iter:iter-1));
            er_elcbo = max(0,elcbo(iter) - max_elcbo)/options.TolImprovement;
            if er_elcbo > 1; er_elcbo = 1 + log(er_elcbo); end

            min_r = min(stats.rindex(min_iter:iter-1));
            er_r = max(0, log(min_r / stats.rindex(iter)));

            % er = 0.5*er_sd + 0.5*er_elcbo;  % Reward
            er = er_r;

            for iHedge = 1:hedge.n
                hedge.g(iHedge) = hedge.decay*hedge.g(iHedge) + er/hedge.phat(iHedge);
            end

            % Apply cutoff value on hedge
            hedge.g = min(hedge.g,HedgeCutoff);
            hedge.g
        end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
clear;clc;close all;
path(path,'toolbox') ;

sk_filename ='../result/simplejoint_v4770_contract_t(4)_nn(30)_WL(9.621282)_WH(1.000000)_sl(3.000000)_skeleton.mat';

%%
SHOW_JOINTS = true;
SHOW_ROOT_JOINT = true;
SHOW_CYCLES = true;
SHOW_IRRELEVANT_EXTRAMA = true;

t1 = 0.1; % for inner branch nodes
a1 = pi*5.0/7.0; % for inner branch nodes,
t2 = 0.1; % for irrelevant extrama;
t3 = 6; % for small cycles;
sprintf('thresholds for remove 1)inner nodes: %f angle, %f length \n 2)irrelevant extrama: %f\n3)small cycles: %d \n',t1, a1, t2, t3)

load(sk_filename,'P');
figure;set(gcf,'color','white');movegui('northwest');set(gcf,'Renderer','OpenGL');view3d rot;
scatter3(P.pts(:,1),P.pts(:,2),P.pts(:,3),20,'.','MarkerEdgeColor', GS.PC_COLOR);hold on;
showoptions.sizep=600;showoptions.sizee=3;
plot_skeleton(P.spls,P.spls_adj,showoptions);
axis off; axis equal; camorbit(0,0,'camera'); axis vis3d; view(-90,0);
%% 0: find joints and NaN, first
[joints, segments] = find_joints(P.pts, P.spls, P.corresp, P.spls_adj, SHOW_JOINTS);
%% 1: removing small cycles measured by topological length
[P.spls, P.corresp, P.spls_adj, joints, segments] = remove_small_cycles(P.spls, P.corresp, P.spls_adj, joints, segments,t3, SHOW_CYCLES);
%% 2: find root joint, global distance relative to root_id, "size of skeleton"
[root_id, global_dist] = find_root_node(P.spls, P.spls_adj, joints, SHOW_ROOT_JOINT);
%% 3: unify joints if there are no branch node inbetween of them. The effect is bad for some cases.
if ~isempty(joints)
    [P.spls, P.corresp, P.spls_adj, joints, root_id] = merge_nearby_joints(P.spls, P.corresp, P.spls_adj, joints, root_id);
end
%% 4: removing irrelevant extrama
[P.spls, P.corresp, P.spls_adj, joints, segments] = remove_irrelevant_extrama(P.spls, P.corresp, P.spls_adj, joints, segments,global_dist, t2, SHOW_IRRELEVANT_EXTRAMA);
%% 5: removing inner branch nodes
[P.spls, P.corresp, P.spls_adj, joints, segments] = remove_inner_nodes(P.spls, P.corresp, P.spls_adj, joints, segments,global_dist,t1,a1);
%% 6: form hierarchy, and rebuild skel_adj remove null node.
P.root_id = root_id;% The first node is root node if there is no joints.
% [P, graph] = build_hierarchy(P,global_dist); % do not contain cycles
[P.spls,P.corresp,P.spls_adj, graph] = build_graph(P.spls,P.corresp,P.spls_adj);% may contain cycles
%% 7: move skeleton nodes to center of its correspondence point cloud
figure('Name','Build hierarchy','NumberTitle','off');set(gcf,'color','white');hold on; movegui('south');
plot_skeleton(P.spls,P.spls_adj);
%scatter3(P.pts(:,1),P.pts(:,2),P.pts(:,3),60,'.','MarkerEdgeColor', GS.PC_COLOR);hold on;
for i = 1:size(P.spls,1)
    verts = P.pts(P.corresp==i,:);
    if size(verts,1) == 1
        P.spls(i,:) = verts;
    else
        P.spls(i,:) = mean(verts);
    end
end
for i = 1:size(P.spls_adj,1)
    links = find( P.spls_adj(i,:)==1 );
    if length(links) == 2
        verts = P.spls(links,:);
        P.spls(i,:) = mean(verts);
    end
end
plot_embedded_graph(P.spls, graph, 'b', 'linewidth', 2);
axis off; axis equal; camorbit(0,0,'camera'); axis vis3d; view(-90,0);view3d rot;
%% finally
figure;set(gcf,'color','white');movegui('northwest');set(gcf,'Renderer','OpenGL');view3d rot;
scatter3(P.pts(:,1),P.pts(:,2),P.pts(:,3),20,'.','MarkerEdgeColor', GS.PC_COLOR);hold on;
plot_skeleton(P.spls, P.spls_adj);
axis off; axis equal; camorbit(0,0,'camera'); axis vis3d; view(0,90);
%%
% P.skel_dist = global_dist;
%P = rmfield(P, 'knn_idx');
%P = rmfield(P, 'knn_dist');

% write_cg();
fid = fopen([sk_filename(1:end-4), '.cg'],'w');
if( fid~=-1 )
    fprintf(fid, '# D:3 NV:%d NE:%d\n', size(P.spls,1), size(graph,1));
    fprintf(fid, 'v %f %f %f\n', P.spls');
    fprintf(fid, 'e %d %d\n', graph');
    fclose(fid);
end

% write_mapping();
fid = fopen([sk_filename(1:end-4), '.map'],'w');
if( fid~=-1 )
    fprintf(fid, '%d\n', P.corresp);
    fclose(fid);
end
save([sk_filename(1:end-4), '_r.mat'], 'P');

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function pass = test_linearSystem2(pref)
% Test solution of a 2x2 system
% Toby Driscoll (extended by AB)

% Note: This test is taken from chebop_systemsolve1 in V4.

%% Setup:
if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

%% Smooth domain:
d = [-1 1];
A = chebop(@(x, u) [diff(u{1}) + u{1} + 2*u{2} ; ...
    diff(u{1}) - u{1} + diff(u{2})], d);
A.lbc = @(u) u{1}+diff(u{1});
A.rbc = @(u) diff(u{2});
x = chebfun('x',d);

f = [ exp(x) ; chebfun(1,d) ];
u = mldivide(A, f, pref);

% Residual of differential equation:
err1 = norm(A(u) - f);

% Want to check BCs as well.
bcFunLeft = A.lbc(u);
bcFunRight = A.rbc(u);
err2 = [norm(bcFunLeft(d(1))), norm(bcFunRight(d(end)))];

% Happy?
pass(1) = norm( err1 ) < tol;
pass(2) = all(err2 < tol);
%% Piecewise:
A.domain = [-1 0 1];
u = mldivide(A, f, pref);

% Norm of the residual:
err3 = norm( A(u) - f );

% Want to check BCs as well:
bcFunLeft = A.lbc(u);
bcFunRight = A.rbc(u);
err4 = [norm(bcFunLeft(d(1))), norm(bcFunRight(d(end)))];

% And check that we're continuous over the breakpoint:
u1jump = jump(u{1}, 0);
u2jump = jump(u{2}, 0);
err5 = [norm(u1jump), norm(u2jump)];

% Happy?
pass(3) = err3 < tol;
pass(4) = all(err4 < tol);
pass(5) = all(err5 < tol);

%% Now use a function handle, rathen than anonymous function.

% Let's change the discretization as well, just because we can:
pref.discretization = @ultraS;
A = chebop(@myop, d);
A.lbc = @(u) u{1}+diff(u{1});
A.rbc = @(u) diff(u{2});

% Here we need to pass information about the number of variables that A
% operates on:
A.numVars = 2;
u = mldivide(A, f, pref);

% Residual of differential equation:
err6 = norm(A(u) - f);

% Want to check BCs as well.
bcFunLeft = A.lbc(u);
bcFunRight = A.rbc(u);
err7 = [norm(bcFunLeft(d(1))), norm(bcFunRight(d(end)))];

% Happy?
pass(6) = norm( err6 ) < tol;
pass(7) = all(err7 < tol);

%% Now try piecewise again!
A.domain = [-1 0 1];
u = mldivide(A, f, pref);

% Norm of the residual:
err8 = norm( A(u) - f );

% Want to check BCs as well:
bcFunLeft = A.lbc(u);
bcFunRight = A.rbc(u);
err9 = [norm(bcFunLeft(d(1))), norm(bcFunRight(d(end)))];

% And check that we're continuous over the breakpoint:
u1jump = jump(u{1}, 0);
u2jump = jump(u{2}, 0);
err10 = [norm(u1jump), norm(u2jump)];

% Happy?
pass(8)  = err8 < tol;
pass(9)  = all(err9 < tol);
pass(10) = all(err10 < tol);

%% Finally, let's include x in the anonmoys function argument list:
d = [-1 1];
A = chebop(@(u) [diff(u{1}) + u{1} + 2*u{2} ; ...
    diff(u{1}) - u{1} + diff(u{2})], d);
A.lbc = @(u) u{1}+diff(u{1});
A.rbc = @(u) diff(u{2});

% And for fun, let's try another discretization
pref.discretization = @ultraS;
u = mldivide(A, f, pref);

% Residual of differential equation:
err11 = norm(A(u) - f);

% Want to check BCs as well.
bcFunLeft = A.lbc(u);
bcFunRight = A.rbc(u);
err12 = [norm(bcFunLeft(d(1))), norm(bcFunRight(d(end)))];

% Happy?
pass(11) = norm( err11 ) < tol;
pass(12) = all(err12 < tol);

end

function out = myop(x, u)

out = [diff(u{1}) + u{1} + 2*u{2} ; ...
    diff(u{1}) - u{1} + diff(u{2})];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Texture Evolution
%
%%

author: mtex-toolbox

repo: mtex

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


text:
function pass = test_compose(pref)
% Test COMPOSE().

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1e3*pref.cheb2Prefs.chebfun2eps;

F = spherefunv(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z);
g = chebfun3(@(x,y,z) x + y + z);
h_true = spherefun(@(x,y,z) x + y + z);
h = compose(F, g);
pass(1) = ( norm(h - h_true) < tol );

G = chebfun3v(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z);
H_true = spherefunv(@(x,y,z) x, @(x,y,z) y, @(x,y,z) z);
H = compose(F, G);
pass(2) = ( norm(H - H_true) < tol );

end


text:
function DataSet = prtDataGenXor(N)
%prtDataGenXor  Generates XOR example data
%
%   DATASET = prtDataGenXor returns a prtDataSetClass with randomly
%   generated data according to the following distribution.
%
%       H0: N([-1 -1],eye(2))
%       H1: N([2 2],[1 .5; .5 1])
%
%   % Example
%
%   ds = prtDataGenXor;
%   plot(ds)
%
%   See also: prtDataSetClass, prtDataGenBiModal, prtDataGenIris,
%   prtDataGenMary, prtDataGenNoisySinc, prtDataGenOldFaithful,
%   prtDataGenSpiral, prtDataGenUnimodal, prtDataGenUnimodal, prtDataGenXor

if nargin == 0
    nSamples = 200;
else
    nSamples = N;
end
if nargin < 5
    mu01 = [0 3];
    sigma01 = eye(2)/3;
    mu02 = [3 0];
    sigma02 = eye(2)/3;

    mu11 = [3 3];
    sigma11 = eye(2)/3;
    mu12 = [0 0];
    sigma12 = eye(2)/3;
end
rv(1) = prtRvMvn('mu',mu01,'sigma',sigma01);
rv(2) = prtRvMvn('mu',mu02,'sigma',sigma02);

rv(3) = prtRvMvn('mu',mu11,'sigma',sigma11);
rv(4) = prtRvMvn('mu',mu12,'sigma',sigma12);

X = cat(1,draw(rv(1),nSamples/2),draw(rv(2),nSamples/2),draw(rv(3),nSamples/2),draw(rv(4),nSamples/2));
Y = prtUtilY(nSamples,nSamples);

DataSet = prtDataSetClass(X,Y,'name','XOR Data');

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% Script: run_loma_yearly_stresstensor.m

% Starting ZMAP
sPath = pwd
cd /home/jowoe/zmap
startup
cd(sPath)

% Load the parameter file
load Params_Loma_ConstRad3km_Nmin50_0.01deg_T365d.mat

% Do loop over different radii to select events
for fRadius = 3:1:5
    sString = ['Radius: ', num2str(fRadius) 'km'];
    params.fRadius = fRadius;
    vResults = [];
    disp(sString)
    tstart1= cputime
    !date
    % Perform the calculation
    [vResults] = gui_CalcStressInv(params);
    !date
    tend1=cputime-tstart1
    !rm tmp*
end

author: CelsoReyes


text:
classdef FminconOptimizer < AbstractGradientOptimizer
    %FminconOptimizer Summary of this class goes here
    %   Detailed explanation goes here

    properties(Access = private)
        options(1,1) FminconOptions = FminconOptions();
    end

    methods
        function obj = FminconOptimizer()
            obj.options = FminconOptions();
        end

        function [exitflag, message] = optimize(obj, lvdOpt, writeOutput, callOutputFcn, hLvdMainGUI)
            [x0All, actVars, varNameStrs] = lvdOpt.vars.getTotalScaledXVector();
            [lbAll, ubAll, lbUsAll, ubUsAll] = lvdOpt.vars.getTotalScaledBndsVector();
            typicalX = lvdOpt.vars.getTypicalScaledXVector();

            if(isempty(x0All) && isempty(actVars))
                exitflag = 0;
                message = 'No variables enabled on script.  Aborting optimization.';

                return;
            end

            bool = x0All < lbAll;
            x0All(bool) = lbAll(bool);

            bool = x0All > ubAll;
            x0All(bool) = ubAll(bool);

            evtNumToStartScriptExecAt = obj.getEvtNumToStartScriptExecAt(lvdOpt, actVars);
            evtToStartScriptExecAt = lvdOpt.lvdData.script.getEventForInd(evtNumToStartScriptExecAt);

            opts = obj.options.getOptionsForOptimizer(typicalX);

            objFuncWrapper = @(x) lvdOpt.objFcn.evalObjFcn(x, evtToStartScriptExecAt);

            if(opts.SpecifyConstraintGradient)
                nonlcon = @(x) lvdOpt.constraints.evalConstraintsWithGradients(x, true, evtToStartScriptExecAt, true, []);
            else
                nonlcon = @(x) lvdOpt.constraints.evalConstraints(x, true, evtToStartScriptExecAt, true, []);
            end

            if(obj.options.computeOptimalStepSizes == true)
                optimalStepSizes = CustomFiniteDiffsCalculationMethod.determineOptimalStepSizes(@(x) CustomFiniteDiffsCalculationMethod.combinedConstrFun(x, lvdOpt.lvdData), x0All, hLvdMainGUI);

                if(not(isempty(optimalStepSizes)))
                    opts.FiniteDifferenceStepSize = optimalStepSizes;
                    opts.TypicalX = ones(length(optimalStepSizes),1);
                end
            end

            if(lvdOpt.gradAlgo == LvdOptimizerGradientCalculationAlgoEnum.BuiltIn)
                objFunToPass = objFuncWrapper;
                opts = optimoptions(opts, 'SpecifyObjectiveGradient',false);
            elseif(lvdOpt.gradAlgo == LvdOptimizerGradientCalculationAlgoEnum.FiniteDifferences)
                gradCalcMethod = lvdOpt.customFiniteDiffsCalcMethod;
                objFunToPass = @(x) obj.objFuncWithGradient(objFuncWrapper, x, gradCalcMethod, obj.usesParallel());
                opts = optimoptions(opts, 'SpecifyObjectiveGradient',true);

                sparsityTF = gradCalcMethod.shouldComputeSparsity();
                if(sparsityTF && not(isempty(hLvdMainGUI)))
%                     hMsgBox = msgbox('Computing sparsity.  Please wait...');
                    hMsgBox = uiprogressdlg(hLvdMainGUI, 'Message','Computing sparsity.  Please wait...', 'Title','Computing Sparsity', 'Indeterminate',true, 'Icon','info');
                else
                    hMsgBox = NaN;
                end

                fAtX0 = objFunToPass(x0All);
                gradCalcMethod.computeGradientSparsity(objFuncWrapper, x0All, fAtX0, obj.usesParallel());

                if(sparsityTF && isgraphics(hMsgBox))
                    close(hMsgBox); drawnow;
                end
            end

            problem = createOptimProblem('fmincon', 'objective',objFunToPass, 'x0', x0All, 'lb', lbAll, 'ub', ubAll, 'nonlcon', nonlcon, 'options', opts);
            problem.lvdData = lvdOpt.lvdData; %need to get lvdData in somehow

            %%% Run optimizer
            recorder = ma_OptimRecorder();
            celBodyData = lvdOpt.lvdData.celBodyData;

            if(callOutputFcn)
                propNames = lvdOpt.lvdData.launchVehicle.tankTypes.getFirstThreeTypesCellArr();
%                 handlesObsOptimGui = ma_ObserveOptimGUI(celBodyData, problem, true, writeOutput, [], varNameStrs, lbUsAll, ubUsAll);

                out = AppDesignerGUIOutput();
                ma_ObserveOptimGUI_App(out);
                handlesObsOptimGui = out.output{1};

%                 outputFnc = @(x, optimValues, state) ma_OptimOutputFunc(x, optimValues, state, handlesObsOptimGui, problem.objective, problem.lb, problem.ub, celBodyData, recorder, propNames, writeOutput, varNameStrs, lbUsAll, ubUsAll);
                hOptimStatusLabel = handlesObsOptimGui.optimStatusLabel;
                hFinalStateOptimLabel = handlesObsOptimGui.finalStateOptimLabel;
                hDispAxes = handlesObsOptimGui.dispAxesPanel;
                hCancelButton = handlesObsOptimGui.cancelButton;
                optimStartTic = tic();

                outputFnc = @(x, optimValues, state) FminconOptimizer.getOutputFunction(x, optimValues, state, hOptimStatusLabel, hFinalStateOptimLabel, hDispAxes, hCancelButton, ...
                                                                                objFuncWrapper, problem.lb, problem.ub, celBodyData, recorder, propNames, writeOutput, varNameStrs, lbUsAll, ubUsAll, optimStartTic);
                problem.options.OutputFcn = outputFnc;
            end

            [exitflag, message] = lvd_executeOptimProblem(celBodyData, writeOutput, problem, recorder, callOutputFcn);

            if(callOutputFcn)
                close(handlesObsOptimGui.ma_ObserveOptimGUI);
            end
        end

        function options = getOptions(obj)
            options = obj.options;
        end

        function setGradientCalculationMethod(obj, newGradCalcMethod)
            obj.gradCalcMethod = newGradCalcMethod;
        end

        function openOptionsDialog(obj)
%             lvd_editFminconOptionsGUI(obj);

            output = AppDesignerGUIOutput({false});
            lvd_editFminconOptionsGUI_App(obj, output);
        end

        function tf = usesParallel(obj)
            tf = obj.options.useParallel.optionVal;
        end

        function numWorkers = getNumParaWorkers(obj)
            numWorkers = obj.options.getNumParaWorkers();
        end
    end

    methods(Access=private)
        function [f, g, stateLog] = objFuncWithGradient(~, objFun, x, gradCalcMethod, useParallel)
            [f, stateLog] = objFun(x);

            if(nargout >= 2)
                g = gradCalcMethod.computeGrad(objFun, x, f, useParallel);
            end
        end
    end

    methods(Static, Access=private)
        function stop = getOutputFunction(x, optimValues, state, hOptimStatusLabel, hFinalStateOptimLabel, hDispAxes, hCancelButton, ...
                                          objFcn, lb, ub, celBodyData, recorder, propNames, writeOutput, varLabels, lbUsAll, ubUsAll, optimStartTic)
            switch state
                case 'iter'
                    stop = get(hCancelButton,'Value');

                    recorder.iterNums(end+1) = optimValues.iteration;
                    recorder.xVals(end+1) = {x};
                    recorder.fVals(end+1) = optimValues.fval;
                    recorder.maxCVal(end+1) = optimValues.constrviolation;
                case {'init','interrupt','done'}
                    stop = get(hCancelButton,'Value');
            end

            if(stop == true)
                return;
            end

            [~, stateLog] = objFcn(x);

%             finalStateLogEntry = stateLog.getFinalStateLogEntry();
%             finalStateLogEntryMA = finalStateLogEntry.getMAFormattedStateLogMatrix(true);

            stateLogMA = stateLog.getMAFormattedStateLogMatrix(true);

            if(strcmpi(state,'init') || strcmpi(state,'iter'))
                try
                    FminconOptimizer.writeOptimStatus(hOptimStatusLabel, optimValues, state, writeOutput, optimStartTic);
                    ma_UpdateStateReadout(hFinalStateOptimLabel, 'final', propNames, stateLogMA, celBodyData);
                    FminconOptimizer.generatePlots(x, optimValues, state, hDispAxes, lb, ub, varLabels, lbUsAll, ubUsAll);
                    drawnow;
                catch ME
                    warning(ME.message);
                end
            end
        end

        function writeOptimStatus(hOptimStatusLabel, optimValues, state, writeOutput, timer)
            elapTime = toc(timer);

            outStr = {};
            outStr{end+1} = ['State                = ', state];
            outStr{end+1} = '                        ';
            outStr{end+1} = ['Iterations           = ', num2str(optimValues.iteration)];
            outStr{end+1} = ['Function Evals       = ', num2str(optimValues.funccount)];
            outStr{end+1} = ['Objective Value      = ', num2str(optimValues.fval)];
            outStr{end+1} = ['Constraint Violation = ', num2str(optimValues.constrviolation)];
            outStr{end+1} = ['Optimality           = ', num2str(optimValues.firstorderopt)];
            outStr{end+1} = ['Step Size            = ', num2str(optimValues.stepsize)];
            outStr{end+1} = '                       ';
            outStr{end+1} = ['Elapsed Time         = ', num2str(elapTime), ' sec'];

            set(hOptimStatusLabel, 'String', outStr);

            switch state
                case 'iter'
                    formatstr = ' %- 12.1i %- 12.0i %- 12.6g %- 12.3g %- 12.3g %- 12.3g';

                    iter = optimValues.iteration;
                    fcnt = optimValues.funccount;
                    val  = optimValues.fval;
                    feas = optimValues.constrviolation;
                    optm = optimValues.firstorderopt;
                    step = optimValues.stepsize;

                    hRow = sprintf(formatstr,iter,fcnt,val,feas,optm,step);
                    writeOutput(hRow,'append');
                case 'init'
                    hdrStr = sprintf('%- 13s%- 13s%- 13s%- 13s%- 13s%- 13s', 'Iteration','Fcn-Count','f(x)-Value', 'Feasibility', 'Optimality', 'Norm. Step');
                    writeOutput(hdrStr,'append');
            end
        end

        function generatePlots(x, optimValues, state, hDispAxes, lb, ub, varLabels, lbUsAll, ubUsAll)
            persistent fValPlotIsLog tLayout hPlot1 hPlot2 hPlot3

            if(isempty(fValPlotIsLog))
                fValPlotIsLog = true;
            end

            switch state
                case 'init'
                    if(isvalid(hDispAxes))
%                         set(hDispAxes,'Visible','on');
%                         subplot(hDispAxes);
                        tLayout = tiledlayout(hDispAxes, 3,1);
%                         axes(hDispAxes);
                    end
                    fValPlotIsLog = true;
            end

            hPlot1 = nexttile(tLayout, 1);
            if(strcmpi(state,'init'))
%                 hPlot1 = subplot(hDispAxes, 3,1,1);

                hPlot1.XTickLabel= [];
                hPlot1.YTickLabel= [];
                hPlot1.ZTickLabel= [];
%                 axes(hPlot1);
            end
            optimplotxKsptot(x, optimValues, state, lb, ub, varLabels, lbUsAll, ubUsAll);

            hPlot2 = nexttile(tLayout, 2);
            if(strcmpi(state,'init'))
%                 hPlot2 = subplot(3,1,2);

                hPlot2.XTickLabel= [];
                hPlot2.YTickLabel= [];
                hPlot2.ZTickLabel= [];
                h = hPlot2;
            else
                h = hPlot2;
%                 axes(hPlot2);
            end
            if(optimValues.fval<=0)
                fValPlotIsLog = false;
                set(h,'yscale','linear');
            end
            optimplotfvalKsptot(x, optimValues, state);
            if(fValPlotIsLog)
                set(h,'yscale','log');
            else
                set(h,'yscale','linear');
            end
            grid on;
            grid minor;

            hPlot3 = nexttile(tLayout, 3);
            if(strcmpi(state,'init'))
%                 hPlot3 = subplot(3,1,3);

                hPlot3.XTickLabel= [];
                hPlot3.YTickLabel= [];
                hPlot3.ZTickLabel= [];
                h = hPlot3;
            else
                h = hPlot3;
%                 axes(hPlot3);
            end
            optimplotconstrviolationKsptot(x, optimValues, state);

            if(not(isempty(h.Children)))
                hLine = h.Children(1);
                if(isa(hLine,'matlab.graphics.chart.primitive.Line'))
                    yDataLine = hLine.YData;
                    if(abs(max(yDataLine) / min(yDataLine)) >= 10 && all(yDataLine > 0))
                        set(h,'yscale','log');
                    else
                        set(h,'yscale','linear');
                    end
                else
                    set(h,'yscale','linear');
                end
            end

            grid on;
            grid minor;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function cost = rcrossvalidatelssvm(model,Y, L, omega, estfct,combinefct)
%%%%%%%%%%%%%%%%%%%%%
% INTERNAL FUNCTION %
%%%%%%%%%%%%%%%%%%%%%
% Estimate the model performance of a model with l-fold robust crossvalidation

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @
% http://www.esat.kuleuven.be/sista/lssvmlab

%
% initialisation and defaults
%
%if size(X,1)~=size(Y,1), error('X and Y have different number of datapoints'); end
[nb_data,y_dim] = size(Y);
d = size(model.xtrain,2);
% LS-SVMlab
eval('model = initlssvm(model{:});',' ');
model.status = 'changed';

eval('L;','L=min(round(sqrt(size(model.xfull,1))),10);');
eval('estfct;','estfct=''mse'';');
eval('combinefct;','combinefct=''mean'';');

py = Y;
[~,Y] = postlssvm(model,[],Y);

gams = model.gamcsa;
eval('sig2s = model.kernel_parscsa;','sig2s=[];')
eval('deltas = model.deltacsa;','deltas=[];')
%
%initialize: no incremental  memory allocation
%
costs = zeros(L,length(gams));
block_size = floor(nb_data/L);

% check whether there are more than one gamma or sigma
for j =1:numel(gams)
    if strcmp(model.kernel_type,'RBF_kernel') || strcmp(model.kernel_type,'RBF4_kernel')
        model = changelssvm(changelssvm(model,'gam',gams(j)),'kernel_pars',sig2s(j));
        eval('model.delta=deltas(j);','model.delta=[];')
    elseif strcmp(model.kernel_type,'lin_kernel')
        model = changelssvm(model,'gam',gams(j));
        eval('model.delta=deltas(j);','model.delta=[];')
    elseif strcmp(model.kernel_type,'poly_kernel')
        model = changelssvm(changelssvm(model,'gam',gams(j)),'kernel_pars',[sig2s(1,j);sig2s(2,j)]);
        eval('model.delta=deltas(j);','model.delta=[];')
    else
        model = changelssvm(changelssvm(model,'gam',gams(j)),'kernel_pars',[sig2s(1,j);sig2s(2,j);sig2s(3,j)]);
        eval('model.delta=deltas(j);','model.delta=[];')
    end

    S = ones(nb_data,1);
    Atot = kernel_matrix2(omega,model.kernel_type,model.kernel_pars,d)+eye(nb_data)./model.gam;
    %
    %
    % start loop over l validations
    %
    for l = 1:L,

        % divide in data and validation set, trainings data set is a copy
        % of permutated_data, validation set is just a logical index
        if l==L,
            train = 1:block_size*(l-1);
            validation = block_size*(l-1)+1:nb_data;
        else
            train = [1:block_size*(l-1) block_size*l+1:nb_data];
            validation = block_size*(l-1)+1:block_size*l;
        end

        A = [0 S(train)';S(train) Atot(train,train)];
        b = [0;py(train)];

        % Solve linear system
        sol = linsolve(A,b,struct('SYM',true));

        % Determine residuals ek
        ek = sol(2:end)./model.gam;
        g = model.gam;
        %for i=2:size(A,1), A(i,i) = A(i,i) - 1/g; end
        A = A-eye(size(train,2)+1)./g; A(1,1)=0;
        Ah = A;
        %
        % robust estimation of the variance
        %
        for k = 1:30
            vare = 1.483*median(abs((ek)-median(ek)));
            alphaold = sol(2:end);
            %
            % robust re-estimation of the alpha's and the b
            %
            cases = reshape((ek./vare),1,size(ek,1));
            W = g*weightingscheme(cases,model.weights,model.delta);

            for t=1:size(train,2), A(t+1,t+1) = A(t+1,t+1)+1./W(t); end

            sol = linsolve(A,b,struct('SYM',true));

            ek = sol(2:end)./W';
            A = Ah;
            if norm(abs(alphaold-sol(2:end)),'fro')<=1e-3,
                %fprintf('\n Converged after %.0f iteration(s)', k);
                k = inf;
            end
            model.status = 'changed';
        end

        % regression
        % Simulate system on validation data
        yh = Atot(train,validation)'*sol(2:end) + ones(numel(validation),1)*sol(1);
        [~,yh] = postlssvm(model,[],yh);
        z = yh - Y(validation,:);
        eval('costs(l,j) = feval(estfct,z);')
    end
end

cost = feval(combinefct, costs);

author: peterhcharlton

repo: RRest


text:
function [ok,tp,tq] = linenear(pa,pb,pc,pd)
%LINENEAR calc. the nearest points on line segments embedded
%in d-dimensions. Line paramters are bounded on [-1,+1].

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 10/10/2017
%-----------------------------------------------------------

    m1 = (pa+pb) * +.5 ;
    D1 = (pb-pa) * +.5 ;

    m2 = (pc+pd) * +.5 ;
    D2 = (pd-pc) * +.5 ;

    r1 = sum(m2.*D1,2) ...
       - sum(m1.*D1,2) ;
    r2 = sum(m1.*D2,2) ...
       - sum(m2.*D2,2) ;

    A1 = sum(D1.*D1,2) ;
    A2 =-sum(D1.*D2,2) ;
    A3 =-sum(D1.*D2,2) ;
    A4 = sum(D2.*D2,2) ;

    dd = A1.*A4 - A2.*A3 ;

    tp = A4.*r1 - A2.*r2 ;
    tq =-A3.*r1 + A1.*r2 ;

    rt = max(abs([A1,A2,A3,A4]),[],2);
    rt = rt * eps ^ .8 ;

    ok = abs(dd) > +rt ;

    tp(~ok) = +0. ;
    tq(~ok) = +0. ;

    tp(ok) = tp(ok) ./ dd(ok) ;
    tq(ok) = tq(ok) ./ dd(ok) ;

end

author: CHLNDDEV


text:
function [A,varargout] = selectMinbyRow(A,varargin)
% find minimum in each row

s = size(A);
[A,ind] = min(A,[],2);
ind = sub2ind(s,1:s(1),ind.');

for i = 1:nargout-1
  varargout{i} = varargin{i}(ind.');
end


text:
% clear classes
% WaitObj = prtUtilProgressBar(0,'Wait some seconds ...','reset',true,'allowStop',true);
% for p=0:0.01:1
%     pause(0.01);
%     WaitObj.update(p);
%     if WaitObj.isCanceled
%         WaitObj.close();
%         break
%     end
% end
%
% %%
% clear classes
% WaitObj = prtUtilProgressBar(0,'Wait some seconds ...','autoClose',false,'reset',true);
% WaitObj2 = prtUtilProgressBar(0,'Wait some more seconds ...','autoClose',false);
% for p=0:0.1:1
%     WaitObj.update(p);
%     for q=0:0.05:1
%         WaitObj2.update(q);
%         pause(0.001);
%     end
% end
% WaitObj.closeAll();
%
% %%
%
% clear classes
% WaitObj = prtUtilProgressBar(0,'Wait some seconds ...','autoClose',false,'reset',true,'allowStop',true);
% WaitObj2 = prtUtilProgressBar(0,'Wait some more seconds ...','autoClose',false);
% for p=0:0.1:1
%     WaitObj.update(p);
%     for q=0:0.051:1
%         WaitObj2.update(q);
%         pause(0.001);
%     end
%     if WaitObj.isCanceled
%         WaitObj.closeAll();
%         break
%     end
% end
% WaitObj.closeAll();
%

classdef prtUtilProgressBar
    properties
        titleStr = 'Progress...';
        autoClose = true;
        reset = false;
    end

    properties (SetAccess = protected)
        figureHandle = [];
    end

    properties (Dependent = true)
        isCanceled
    end

    properties (SetAccess = private, GetAccess = private, Hidden = true)
        isCanceledLocal = false;
    end

    properties (Hidden = true)
        allowStop = false;
        isDead = false;
        globalTitleStr = '-auto';

        barIndex  % Used internally

        % Restricted update parameters
        minTimeBetweenUpdates = 0.2; % in seconds
        minPercentBetweenUpdates = 0; % in percent [0 100] (0 ignores this)
    end

    properties (Constant = true, GetAccess = protected);
        % Graphics parameters
        windowWidth = 512;
        windowHorzMargin = 16; % On both sides
        windowTopMargin = 24;  % Distance between top of window and top of top bar;
        windowBottomMargin = 16; % Distance between bottom of window and bottom cancel button

        barBackgroundColor = [1 1 1];
        barForegroundColor = [0.2157    0.4941    0.7216];

        barWidthMin = 0.1;
        barHeight = 22;
        barVertMargin = 24;
        titleTextMargin = 8;

        cancelButtonHeight = 36;
        cancelButtonVertMargin = 16; % Distance between top of cancel button and bottom of first bar
        cancelButtonWidth = 176;

        textHorzMargin = 4;

        % Time Estimate Smoothing
        nSamplesForTimeEstSmoothing = 25;
    end

    methods
        function Obj = set.isCanceled(Obj,val)
            Obj.isCanceledLocal = val;
        end
        function val = get.isCanceled(Obj)

            if ishandle(Obj.figureHandle)
                myFig = Obj.figureHandle;
                ud = get(myFig,'userdata');

                if ~isempty(ud) && isfield(ud,'prtUtilProgressBarIsCanceled')
                    val = (Obj.allowStop && ud.prtUtilProgressBarIsCanceled) || Obj.isCanceledLocal;
                else
                    val = false;
                end

                prtUtilProgressBarData = guidata(Obj.figureHandle);
                if ~isempty(prtUtilProgressBarData)

                else
                    val = false;
                end
            else
                val = false;
            end

%             if ishandle(Obj.figureHandle)
%                 prtUtilProgressBarData = guidata(Obj.figureHandle);
%                 if ~isempty(prtUtilProgressBarData)
%                     val = (Obj.allowStop && prtUtilProgressBarData.isCanceled) || Obj.isCanceledLocal;
%                 else
%                     val = false;
%                 end
%             else
%                 val = false;
%             end
        end

        function Obj = prtUtilProgressBar(percentage,titleStr,varargin)
            %Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
            if nargin < 1
                percentage = 0;
            end
            if nargin > 1
                Obj.titleStr = titleStr;
            end

            Obj = prtUtilAssignStringValuePairs(Obj, varargin{:});

            % Attempt to locate parent
            oldHandleVisilibity = get(0,'ShowHiddenHandles');
            set(0,'ShowHiddenHandles','on');
            oldFigureHandle = findobj('tag','PrtProgressBar','HandleVisibility','off');
            set(0,'ShowHiddenHandles',oldHandleVisilibity);

            if Obj.reset && ~isempty(oldFigureHandle)
                close(oldFigureHandle)
                oldFigureHandle = [];
            end

            if isempty(oldFigureHandle)
                % No existing prtUtilProgressBar windows open
                Obj.barIndex = 1;

                Obj = Obj.createWindow();
                Obj = Obj.addBar();

            else
                % There is a prtUtilProgressBar window open
                % We want to tack into that

                Obj.figureHandle = oldFigureHandle;

                prtProgressBarData = guidata(Obj.figureHandle);
                Obj.barIndex = prtProgressBarData.nBars + 1;

                Obj = Obj.addBar();
            end

            Obj.update(percentage); % Actually update the bars so that we can initialize partially done
        end
        function Obj = update(Obj, percentage)
            assert(numel(percentage)==1 &&  ~(percentage > 1+eps || percentage < -eps),'prt:prtUtilProgressBar','progress must be a scalar percentage in the range [0 1]');

            if  Obj.isCanceled || Obj.isDead || ~ishandle(Obj.figureHandle) || strcmpi(get(Obj.figureHandle,'BeingDeleted'),'on')
                Obj.isCanceled = true;
                Obj.isDead = false;
                return
            end

            % Must get bar structure from figure's guidata
            prtUtilProgressBarData = guidata(Obj.figureHandle);

            if Obj.barIndex > prtUtilProgressBarData.nBars
                % We deleted at least one bar (possibly through auto close
                % We gotta create a new bar(s)

                nBarsNeeded = Obj.barIndex-prtUtilProgressBarData.nBars;
                for iNewBar = 1:nBarsNeeded
                    Obj = addBar(Obj);
                end
                Obj.update(percentage);
                return;
            end

            bar = prtUtilProgressBarData.bars(Obj.barIndex);

            % Get Current time with seconds = 1 instead of days = 1 like
            % now() provides
            currentTime = 86400*now();

            bar.percentage = percentage;

            % We only want to store times and percentages if the
            % percentages are different than the last time we saved.
            % Store iteration times and percentages
            if bar.lastIterationsPercentages(end) ~= percentage
                % An actual change occured for this bar
                if bar.percentage == 0
                    % Much like Donnie we have no frame of reference
                    bar.lastIterationsTimes(:) = nan;
                    bar.lastIterationsPercentages(:) = 0;
                    bar.startTime = currentTime;
                    bar.estimatedRemainingTime = Inf;
                else
                    bar.lastIterationsTimes = circshift(bar.lastIterationsTimes,[0 -1]);
                    bar.lastIterationsTimes(end) = currentTime;

                    bar.lastIterationsPercentages = circshift(bar.lastIterationsPercentages,[0 -1]);
                    bar.lastIterationsPercentages(end) = percentage;
                end
            end

            % Update Patch position
            cPatchPosition = get(bar.foregroundPatchHandle,'Position');
            totalPatchArea = Obj.windowWidth-cPatchPosition(1)*2; % assumes equal x margin on both sides, which is currently true

            % Update Bar
            set(bar.foregroundPatchHandle,'Position',[cPatchPosition(1:2), min(max(0.1,totalPatchArea*percentage),totalPatchArea), cPatchPosition(4)]);

            % Update title string (sub title)
            set(bar.topTextHandle,'string',Obj.titleStr);

            % Update texts
            if sum(isfinite(bar.lastIterationsTimes)) == 1
                % First update with percetage > 0
                % Must reference start time instead of previously saved
                % iteration times
                timeChanges = currentTime - bar.startTime;
                percentChanges = percentage;
            else
                timeChanges = diff(bar.lastIterationsTimes);
                percentChanges = diff(bar.lastIterationsPercentages);
            end

            % Remove negative steps...
            percentChanges = percentChanges(percentChanges>=0);
            timeChanges = timeChanges(percentChanges>=0);

            % Based on the current percentages and the back log of
            % percentages estime the time remaining

            estimatedTimesPerPercent = timeChanges./percentChanges;
            estimatedSecondsPerPercent = prtUtilNanMean(estimatedTimesPerPercent(~isnan(estimatedTimesPerPercent)));

            if isempty(estimatedSecondsPerPercent)
                estimatedRemainingTimeInSeconds = inf;
            else
                estimatedRemainingTimeInSeconds = round(estimatedSecondsPerPercent*(1-percentage));
            end

            if abs(estimatedRemainingTimeInSeconds-bar.estimatedRemainingTime) > 0.75 % Only update if we differ significantly. This eliminates jitter.
                bar.estimatedRemainingTime = estimatedRemainingTimeInSeconds;

                timeDifferentials = currentTime - cat(1,prtUtilProgressBarData.bars.timeLastGraphicsUpdate);
                remainingTimes = cat(1,prtUtilProgressBarData.bars.estimatedRemainingTime) - timeDifferentials;
                remainingTimes(remainingTimes<0) = 0;
                elapsedTimes = currentTime - cat(1,prtUtilProgressBarData.bars.startTime);

                for iBar = 1:prtUtilProgressBarData.nBars

                    if isfinite(remainingTimes(iBar))
                        estimatedRemainingTimeInSeconds = fix(remainingTimes(iBar));
                        estimatedRemainingTimeInMinutes = estimatedRemainingTimeInSeconds/60;
                        estimatedRemainingHours = floor(estimatedRemainingTimeInMinutes/60);
                        estimatedRemainingMinutes = floor(mod(estimatedRemainingTimeInMinutes,60));
                        estimatedRemainingSeconds = mod(estimatedRemainingTimeInSeconds,60);

                        timeRemainingStr = sprintf('%u:%02u:%02u',estimatedRemainingHours,estimatedRemainingMinutes,estimatedRemainingSeconds);

                        set(prtUtilProgressBarData.bars(iBar).rightTextHandle,'string',timeRemainingStr);
                    end

                    if isfinite(elapsedTimes(iBar))
                        elapsedRemainingTimeInSeconds = fix(elapsedTimes(iBar));
                        elapsedRemainingTimeInMinutes = elapsedRemainingTimeInSeconds/60;
                        elapsedRemainingHours = floor(elapsedRemainingTimeInMinutes/60);
                        elapsedRemainingMinutes = floor(mod(elapsedRemainingTimeInMinutes,60));
                        elapsedRemainingSeconds = mod(elapsedRemainingTimeInSeconds,60);

                        timeElapsedStr = sprintf('%u:%02u:%02u',elapsedRemainingHours,elapsedRemainingMinutes,elapsedRemainingSeconds);

                        set(prtUtilProgressBarData.bars(iBar).leftTextHandle,'string',timeElapsedStr);
                    end
                end
            end

            % Always update the center
            set(bar.centerTextHandle,'string', sprintf('%d%%',round(percentage*100)));

            if ((currentTime-bar.timeLastGraphicsUpdate) > Obj.minTimeBetweenUpdates) || ((percentage-bar.percentLastGraphicsUpdate) > Obj.minPercentBetweenUpdates)
                bar.timeLastGraphicsUpdate = currentTime;
                bar.percentLastGraphicsUpdate = percentage;

                Obj.updateFigureName(prtUtilProgressBarData);

                drawnow;
            end

            % Save bar data in figure
            prtUtilProgressBarData.bars(Obj.barIndex) = bar;
			if ishandle(Obj.figureHandle)
				guidata(Obj.figureHandle,prtUtilProgressBarData);
			end

            if all(abs(percentage-1)<1e-9)

                if Obj.autoClose
                    % close the object
                    Obj.close();
                end
            end

        end

        function Obj = set.globalTitleStr(Obj,str)
            assert(ischar(str),'globalTitleStr must be a character array');

            oldGlobalTitleStr = Obj.globalTitleStr;
            Obj.globalTitleStr = str;

            if ishandle(Obj.figureHandle) %#ok<MCSUP>
                prtUtilProgressBarData = guidata(Obj.figureHandle); %#ok<MCSUP>
                isAuto = strcmpi(str,'-auto');
                if (isAuto && strcmpi(oldGlobalTitleStr,prtUtilProgressBarData.globalTitleStr)) || ~isAuto
                    prtUtilProgressBarData.globalTitleStr = Obj.globalTitleStr;
                    guidata(Obj.figureHandle, prtUtilProgressBarData); %#ok<MCSUP>
                end
            end
        end

        function updateFigureName(Obj,gData)

            if ~ishandle(Obj.figureHandle)
                return
            end

            if nargin < 2 || isempty(gData)
                gData = guidata(Obj.figureHandle);
            end

            if ~isfield(gData,'globalTitleStr')
                guidata(Obj.figureHandle,gData);
                Obj.globalTitleStr = Obj.globalTitleStr;
                gData = guidata(Obj.figureHandle);
            end

            if strcmpi(gData.globalTitleStr,'-auto')
                % Automatically set to the active sub-titles title
                set(Obj.figureHandle,'Name',Obj.titleStr);
            else
                % Set to the global one.
                set(Obj.figureHandle,'Name',gData.globalTitleStr);
            end
        end

        function Obj = close(Obj)

            Obj = removeBar(Obj);

            prtUtilProgressBarData = guidata(Obj.figureHandle);

            if prtUtilProgressBarData.nBars == 0 && ishandle(Obj.figureHandle) && strcmpi(get(Obj.figureHandle,'BeingDeleted'),'off')
                delete(Obj.figureHandle);
            end
        end

        function Obj = closeAll(Obj)
            if ishandle(Obj.figureHandle) && strcmpi(get(Obj.figureHandle,'BeingDeleted'),'off')
                delete(Obj.figureHandle);
            end
        end

        function Obj = set.titleStr(Obj,titleStr)
            assert(ischar(titleStr),'titleStr must be a character array');
            Obj.titleStr = titleStr;
            if ishandle(Obj.figureHandle) %#ok<MCSUP>
                set(Obj.figureHandle,'Name',titleStr); %#ok<MCSUP>
                prtUtilProgressBarData = guidata(Obj.figureHandle); %#ok<MCSUP>
                set(prtUtilProgressBarData.bars(Obj.barIndex).topTextHandle,'string',Obj.titleStr); %#ok<MCSUP>
            end

        end
    end

    methods (Access = private)
        function Obj = createWindow(Obj)

            nBars = 1;
            if Obj.allowStop
                windowHeight = Obj.windowBottomMargin + Obj.windowTopMargin + Obj.cancelButtonHeight + Obj.cancelButtonVertMargin + (nBars-1)*Obj.barVertMargin + nBars*Obj.barHeight;
            else
                windowHeight = Obj.windowBottomMargin + Obj.windowTopMargin + Obj.cancelButtonVertMargin + (nBars-1)*Obj.barVertMargin + nBars*Obj.barHeight;
            end

            windowSize = [Obj.windowWidth windowHeight];

            cancelButtonXStart = floor((Obj.windowWidth - Obj.cancelButtonWidth)/2);
            cancelButtonYStart = Obj.windowBottomMargin;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % Create Figure
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            screenSize = get(0,'ScreenSize');
            Obj.figureHandle = figure('DoubleBuffer','on',...
                                      'HandleVisibility','off',...
                                      'MenuBar','none',...
                                      'NumberTitle','off',...
                                      'Resize','off',...
                                      'Tag','PrtProgressBar',...
                                      'ToolBar','none',...
                                      'DockControls','Off',...%'WindowStyle','modal',...
                                      'Name',Obj.titleStr,...
                                      'Position',[floor((screenSize(3:4)-windowSize)/2) windowSize]);...% Center window
                                      % 'CloseRequestFcn','');

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

            % Create Axes
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            axesHandle = axes('Parent',Obj.figureHandle,...
                                  'Position',[0 0 1 1],...
                                  'Visible','off',...
                                  'XLim',[0 windowSize(1)],...
                                  'YLim',[0 windowSize(2)]);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % Cancel Button
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if Obj.allowStop
                cancelButtonHandle = uicontrol('Parent',Obj.figureHandle,...
                                                   'Units','pixels',...
                                                   'position',[cancelButtonXStart cancelButtonYStart Obj.cancelButtonWidth Obj.cancelButtonHeight],...
                                                   'style','pushbutton',...
                                                   'callback',@(myHandle,evenData)prtUtilProgressBar.cancelCallBack(myHandle,evenData),...
                                                   'String','Stop','BusyAction','Queue');
            else
                cancelButtonHandle = [];
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            prtUtilProgressBarData.nBars = 0;
            prtUtilProgressBarData.axesHandle = axesHandle;
            prtUtilProgressBarData.cancelButtonHandle = cancelButtonHandle;
            prtUtilProgressBarData.bars = struct([]);
            prtUtilProgressBarData.isCanceled = false;
            prtUtilProgressBarData.nBarsSpace = 0;
            prtUtilProgressBarData.globalTitleStr = Obj.globalTitleStr;

            guidata(Obj.figureHandle,prtUtilProgressBarData);

        end

        function Obj = removeBar(Obj)

            % Get Figure user data
            prtUtilProgressBarData = guidata(Obj.figureHandle);
            prtUtilProgressBar.deleteBarHandles(prtUtilProgressBarData.bars(Obj.barIndex));
            prtUtilProgressBarData.bars(Obj.barIndex) = [];

            prtUtilProgressBarData.nBars = prtUtilProgressBarData.nBars - 1;
            nBars = prtUtilProgressBarData.nBars; % Shorthand
            nBarsSpace = prtUtilProgressBarData.nBarsSpace;
            % Derived Graphics parameters
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            hasStop = ~isempty(prtUtilProgressBarData.cancelButtonHandle);

            if hasStop
                barYStarts = Obj.windowBottomMargin + Obj.cancelButtonHeight + Obj.cancelButtonVertMargin + ((1:nBarsSpace)-1)*(Obj.barHeight + Obj.barVertMargin);
            else
                barYStarts = Obj.windowBottomMargin + ((1:nBarsSpace)-1)*(Obj.barHeight + Obj.barVertMargin);
            end

            barYStarts = flipud(barYStarts(:));
            barYStarts = barYStarts(1:nBars);

            for iBar = 1:nBars
                cBarYStart = barYStarts(iBar);

                % An old bar, need to move the y positions up
                BarStruct = prtUtilProgressBarData.bars(iBar);

                handleFieldNamesToChange = {'backgroundPatchHandle';
                    'foregroundPatchHandle';};

                for iField = 1:length(handleFieldNamesToChange)
                    cHandle = BarStruct.(handleFieldNamesToChange{iField});
                    cPos = get(cHandle,'position');
                    cPos(2) = cBarYStart;
                    set(cHandle,'position', cPos);
                end

                handleFieldNamesToChange = {'leftTextHandle';
                    'rightTextHandle';
                    'centerTextHandle';};

                for iField = 1:length(handleFieldNamesToChange)
                    cHandle = BarStruct.(handleFieldNamesToChange{iField});
                    cPos = get(cHandle,'position');
                    cPos(2) = cBarYStart + Obj.barHeight/2;
                    set(cHandle,'position', cPos);
                end

                cHandle = BarStruct.topTextHandle;
                cPos = get(cHandle,'position');
                cPos(2) = cBarYStart + Obj.barHeight + Obj.titleTextMargin;
                set(cHandle,'position', cPos);
            end

            guidata(Obj.figureHandle,prtUtilProgressBarData);

            Obj.isDead = true;

        end

        function Obj = addBar(Obj)

            % Get Figure user data
            prtUtilProgressBarData = guidata(Obj.figureHandle);

            prtUtilProgressBarData.nBars = prtUtilProgressBarData.nBars + 1;
            nBars = prtUtilProgressBarData.nBars; % Shorthand
            nBarsSpace = prtUtilProgressBarData.nBarsSpace; % Shorthand
            cBar = nBars;

            % Derived Graphics parameters
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            xStart = Obj.windowHorzMargin;
            xStop = Obj.windowWidth - Obj.windowHorzMargin;

            barWidth = xStop - xStart;

            textLeftX = xStart + Obj.textHorzMargin;
            textCenterX = mean([xStop xStart]);
            textRightX = xStop - Obj.textHorzMargin;

            hasStop = ~isempty(prtUtilProgressBarData.cancelButtonHandle);

            if nBars > nBarsSpace

                % Expand window
                if hasStop
                    windowHeight = Obj.windowBottomMargin + Obj.windowTopMargin + Obj.cancelButtonHeight + Obj.cancelButtonVertMargin + (nBars-1)*Obj.barVertMargin + nBars*Obj.barHeight;
                else
                    windowHeight = Obj.windowBottomMargin + Obj.windowTopMargin + Obj.cancelButtonVertMargin + (nBars-1)*Obj.barVertMargin + nBars*Obj.barHeight;
                end
                windowSize = [Obj.windowWidth windowHeight];

                screenSize = get(0,'ScreenSize');
                set(Obj.figureHandle,'Position',[floor((screenSize(3:4)-windowSize)/2) windowSize]);% Center window

                set(prtUtilProgressBarData.axesHandle,'Position',[0 0 1 1],...
                                                   'XLim',[0 windowSize(1)],...
                                                   'YLim',[0 windowSize(2)]);

                prtUtilProgressBarData.nBarsSpace = nBars;
                nBarsSpace = nBars;
            end

            % Force Set function call; Be sure to juggle guidata properly
            guidata(Obj.figureHandle,prtUtilProgressBarData);
            Obj.globalTitleStr = Obj.globalTitleStr;
            prtUtilProgressBarData = guidata(Obj.figureHandle);

            if hasStop
                %barYStarts = cancelButtonYStart + Obj.cancelButtonHeight + Obj.cancelButtonVertMargin + ((1:nBars)-1)*(Obj.barHeight + Obj.barVertMargin);
                barYStarts = Obj.windowBottomMargin + Obj.cancelButtonHeight + Obj.cancelButtonVertMargin + ((1:nBarsSpace)-1)*(Obj.barHeight + Obj.barVertMargin);
            else
                %barYStarts = cancelButtonYStart + ((1:nBars)-1)*(Obj.barHeight + Obj.barVertMargin);
                barYStarts = Obj.windowBottomMargin + ((1:nBarsSpace)-1)*(Obj.barHeight + Obj.barVertMargin);
            end
            barYStarts = flipud(barYStarts(:));
            barYStarts = barYStarts(1:nBars);

            for iBar = 1:nBars
                cBarYStart = barYStarts(iBar);

                if iBar == cBar
                    % The new bar

                    % Create Bars
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    BarStruct.backgroundPatchHandle = rectangle('Position',[xStart cBarYStart barWidth Obj.barHeight],...
                        'EdgeColor','none',...
                        'FaceColor',Obj.barBackgroundColor,...
                        'parent',prtUtilProgressBarData.axesHandle);
                    BarStruct.foregroundPatchHandle = rectangle('Position',[xStart cBarYStart Obj.barWidthMin Obj.barHeight],...
                        'EdgeColor','none',...
                        'FaceColor',Obj.barForegroundColor,...
                        'parent',prtUtilProgressBarData.axesHandle);

                    % Create Texts
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    BarStruct.leftTextHandle = text(textLeftX,cBarYStart + Obj.barHeight/2,'','HorizontalAlignment','Left','VerticalAlignment','middle','parent',prtUtilProgressBarData.axesHandle);
                    BarStruct.rightTextHandle = text(textRightX,cBarYStart + Obj.barHeight/2,'','HorizontalAlignment','Right','VerticalAlignment','middle','parent',prtUtilProgressBarData.axesHandle);
                    BarStruct.centerTextHandle = text(textCenterX,cBarYStart + Obj.barHeight/2,'','HorizontalAlignment','Center','VerticalAlignment','middle','parent',prtUtilProgressBarData.axesHandle);
                    BarStruct.topTextHandle = text(xStart, cBarYStart + Obj.barHeight + Obj.titleTextMargin, Obj.titleStr,'HorizontalAlignment','Left','VerticalAlignment','middle','parent',prtUtilProgressBarData.axesHandle);

                    % Initialize other variables
                    BarStruct.lastIterationsTimes = nan(1,Obj.nSamplesForTimeEstSmoothing);
                    BarStruct.lastIterationsPercentages = nan(1,Obj.nSamplesForTimeEstSmoothing);
                    BarStruct.startTime = 86400*now;
                    BarStruct.estimatedRemainingTime = inf;
                    BarStruct.percentage = 0;

                    BarStruct.timeLastGraphicsUpdate = -Inf;
                    BarStruct.percentLastGraphicsUpdate = -Inf;

                    if isempty(prtUtilProgressBarData.bars)
                        prtUtilProgressBarData.bars = BarStruct;
                    else
                        prtUtilProgressBarData.bars = cat(1,prtUtilProgressBarData.bars,BarStruct);
                    end
                else
                    % An old bar, need to move the y positions up
                    BarStruct = prtUtilProgressBarData.bars(iBar);

                    handleFieldNamesToChange = {'backgroundPatchHandle';
                                                'foregroundPatchHandle';};

                    for iField = 1:length(handleFieldNamesToChange)
                        cHandle = BarStruct.(handleFieldNamesToChange{iField});
                        cPos = get(cHandle,'position');
                        cPos(2) = cBarYStart;
                        set(cHandle,'position', cPos);
                    end

                    handleFieldNamesToChange = {'leftTextHandle';
                                                'rightTextHandle';
                                                'centerTextHandle';};

                    for iField = 1:length(handleFieldNamesToChange)
                        cHandle = BarStruct.(handleFieldNamesToChange{iField});
                        cPos = get(cHandle,'position');
                        cPos(2) = cBarYStart + Obj.barHeight/2;
                        set(cHandle,'position', cPos);
                    end

                    cHandle = BarStruct.topTextHandle;
                    cPos = get(cHandle,'position');
                    cPos(2) = cBarYStart + Obj.barHeight + Obj.titleTextMargin;
                    set(cHandle,'position', cPos);

                end
            end

            guidata(Obj.figureHandle,prtUtilProgressBarData);

        end
    end

    methods (Static)
        function cancelCallBack(myHandle,evenData) %#ok<INUSD>
            set(myHandle,'String','Stopping...');
            %prtUtilProgressBarData = guidata(get(myHandle,'parent'));

            myFig = get(myHandle,'parent');
            ud = get(myFig,'userdata');
            ud.prtUtilProgressBarIsCanceled = true;
            set(myFig,'userData',ud);
            %prtUtilProgressBarData.isCanceled = true;
            %guidata(get(myHandle,'parent'), prtUtilProgressBarData);
            drawnow;
        end
        function deleteBarHandles(barStruct)
            delete(barStruct.leftTextHandle);
            delete(barStruct.rightTextHandle);
            delete(barStruct.centerTextHandle);
            delete(barStruct.backgroundPatchHandle);
            delete(barStruct.foregroundPatchHandle);
            delete(barStruct.topTextHandle);
        end

        function forceClose()

            oldHandleVisilibity = get(0,'ShowHiddenHandles');
            set(0,'ShowHiddenHandles','on');
            figHandle = findobj('tag','PrtProgressBar','HandleVisibility','off');
            set(0,'ShowHiddenHandles',oldHandleVisilibity);

            if ~isempty(figHandle)
                close(figHandle)
            end
        end
    end
end


text:
function [y,feature_names,cache] = ComputeDiffNeighborMeanWindowFeatures(x,varargin)

x = x(:)';
N = numel(x);
y = nan(0,N);
feature_names = {};

%% default parameters

[windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii] = ...
  SetDefaultWindowParameters();

% use all transformation types by default
%trans_types = 'all';
trans_types = uint8(15);

% for debugging purposes
SANITY_CHECK = true;

% whether to cache results
DOCACHE = true;

% initialize empty cache
cache = InitializeCache();

% initialize feature_types already computed to empty
feature_types = {};

relativeParams = [];

%% parse parameters

[...
  windows,...
  window_radii,window_offsets,...
  min_window_radius,max_window_radius,nwindow_radii,...
  feature_types,...
  trans_types,...
  SANITY_CHECK,...
  DOCACHE,...
  cache,...
  relativeParams,...
  ] = myparse(varargin,...
  'windows',windows,...
  'window_radii',window_radii,'window_offsets',window_offsets,...
  'min_window_radius',min_window_radius,'max_window_radius',max_window_radius,'nwindow_radii',nwindow_radii,...
  'feature_types',feature_types,...
  'trans_types',trans_types,...
  'sanitycheck',SANITY_CHECK,...
  'docache',DOCACHE,...
  'cache',cache,...
  'relativeParams',relativeParams); %#ok<ASGLU>

%% whether we've specified to use all trans types by default
%if ischar(trans_types) && strcmpi(trans_types,'all'),
%  trans_types = {'none','abs','flip','relative'};
%end

%% select default windows from various ways of specifying windows

[windows,window_radii,windowi2radiusi,nradii] = ...
  SetWindowParameters(...
  windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii);

%% main computation

%if ismember('relative',trans_types)
if bitand(8,trans_types)
  if DOCACHE && ~isempty(cache.relX)
    modX = cache.relX;
  else
    modX = convertToRelative(x,relativeParams);
    cache.relX = modX;
  end
end

for radiusi = 1:nradii,
  r = window_radii(radiusi);
  w = 2*r+1;

  % doesn't make sense for r == 0
  if r == 0,
    continue;
  end

  if DOCACHE && ismember(r,cache.mean.radii),
    cache_i = find(r == cache.mean.radii,1);
    res = cache.mean.data{cache_i};
  else
    res = MeanWindowCore(x,w);
    if DOCACHE,
      cache.mean.radii(end+1) = r;
      cache.mean.data{end+1} = res;
    end
  end

%  if ismember('relative',trans_types)
  if bitand(8,trans_types)
    if DOCACHE && ismember(r,cache.meanRel.radii),
      cache_i = find(r == cache.meanRel.radii,1);
      resRel = cache.meanRel.data{cache_i};
    else
      resRel = MeanWindowCore(modX,w);
      % store for future computations
      if DOCACHE,
        cache.meanRel.radii(end+1) = r;
        cache.meanRel.data{end+1} = resRel;
      end
    end
  end

  % all offsets for this radius
  windowis = find(windowi2radiusi == radiusi);
  for windowi = windowis',
    off = windows(windowi,2);
    % frame t, radius r, offset off:
    % [t-r+off, t+r+off]
    % so for r = 0, off = 1, we want [t+1,t+1]
    % which corresponds to res(t+r+off)
    % so we want to grab for 1+r+off through N+r+off
    res1 = x - padgrab2(res,nan,1,1,1+r+off,N+r+off);

%    if ismember('none',trans_types),
    if bitand(1,trans_types),
      y(end+1,:) = res1; %#ok<*AGROW>
      feature_names{end+1} = {'stat','diff_neighbor_mean','trans','none','radius',r,'offset',off};
    end

%    if ismember('abs',trans_types),
    if bitand(2,trans_types),
      y(end+1,:) = abs(res1);
      feature_names{end+1} = {'stat','diff_neighbor_mean','trans','abs','radius',r,'offset',off};
    end

%    if ismember('flip',trans_types),
    if bitand(4,trans_types),
      res2 = res1.*sign(x);
      y(end+1,:) = res2;
      feature_names{end+1} = {'stat','diff_neighbor_mean','trans','flip','radius',r,'offset',off};
    end

%    if ismember('relative',trans_types),
    if bitand(8,trans_types),
      resRel1 = modX - padgrab2(resRel,nan,1,1,1+r+off,N+r+off);
      y(end+1,:) = resRel1;
      feature_names{end+1} = {'stat','diff_neighbor_mean','trans','relative','radius',r,'offset',off};
    end

    if SANITY_CHECK,
      funcType = 'DiffNeighborMean';
%      if ismember('none',trans_types),
      if bitand(1,trans_types),
        fastY = res1; %#ok<*AGROW>
        res_dumb = nan(1,N);
        for n_dumb = 1:N,
          res_dumb(n_dumb) = x(n_dumb) - nanmean(padgrab2(x,nan,1,1,n_dumb-r+off,n_dumb+r+off));
        end
        checkSanity(fastY,res_dumb,r,off,funcType,'none');
      end

%      if ismember('abs',trans_types),
      if bitand(2,trans_types),
        fastY = abs(res1);
        res_dumb = nan(1,N);
        for n_dumb = 1:N,
          res_dumb(n_dumb) = abs(x(n_dumb) - nanmean(padgrab2(x,nan,1,1,n_dumb-r+off,n_dumb+r+off)));
        end
        checkSanity(fastY,res_dumb,r,off,funcType,'abs');
      end

%      if ismember('flip',trans_types),
      if bitand(4,trans_types),
        res2 = res1; res2(x<0) = -res2(x<0); fastY = res2;
        res_dumb = nan(1,N);
        for n_dumb = 1:N,
          res_dumb(n_dumb) = x(n_dumb) - nanmean(padgrab2(x,nan,1,1,n_dumb-r+off,n_dumb+r+off));
          if x(n_dumb) < 0,
            res_dumb(n_dumb) = -res_dumb(n_dumb);
          end
        end
        checkSanity(fastY,res_dumb,r,off,funcType,'flip');
      end

    end

  end
end

author: kristinbranson

repo: JAABA


text:
function im = clouds(dim)
    fname = '/tmp/zz.ppm';

    if nargin < 1
        dim = [256 256];
    end

    system( sprintf('ppmforge -clouds -width %d -height %d 1> %s 2> /dev/null', dim(1), dim(2), fname) );

    [im,n] = imread(fname, 'PPM');

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [vDrate] = calc_InterpBgrate(mGrid,fMmin)
% function [vDrate] = calc_InterpBgrate(mGrid,fMmin)
% --------------------------------------------------
% Get background rates from a Frankel model on the grid points of a rate
% change grid. Sums the rates for magnitudes larger than fMmin
%
% Incoming:
% mGrid : X-Y grid to interpolate on ([Lon Lat])
% fMmin : Minimum magnitude (generally the same as Mc)
%
% Outgoing:
% vDrate : Vector of DAILY seismicity rates
%
% last update: 09.07.2004
% jochen.woessner@sed.ethz.ch

% Load Frankel model with daily background rate
% Contains variable nullRates :[Lonmin Lonmax Latmin Latmax Dpethmin Depthmax Magmin Magmax DailyRate
% Weight]

% Load CFS file
[sFilename, sPathname] = uigetfile('*.mat', 'Pick Frankel DailyRate MAT-file');
sHelp = [sPathname sFilename];
sFile = [sFilename(1:length(sFilename)-4)];
load(sHelp)

% Loop over grid node positions
for nCnt = 1:length(mGrid(:,1))
    fXcoord = mGrid(nCnt,1);
    fYcoord = mGrid(nCnt,2);
    vSel = (fXcoord >= nullRates(:,1) & fXcoord < nullRates(:,2) &...
        fYcoord >= nullRates(:,3) & fYcoord < nullRates(:,4) & fMmin > nullRates(:,7));
    if ~isempty(nullRates(vSel,9))
        fAllSumMagRate = max(cumsum(nullRates(vSel,9)));
        vDrate(nCnt) = fAllSumMagRate;
    else
        vDrate(nCnt) = nan;
    end
end