text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)
% example code for multi-scale predictions, fusion and accuracy evaluation

function demo_predict_mscale_nyud()

rng('shuffle');

addpath('./my_utils');
dir_matConvNet='../libs/matconvnet/matlab';
run(fullfile(dir_matConvNet, 'vl_setupnn.m'));

run_config=[];
ds_config=[];

run_config.use_gpu=true;
% run_config.use_gpu=false;
run_config.gpu_idx=1;

%-------------------------------------------------------------------------------------------------------------------------
% settings for using trained model:

ds_name_subfix='nyud';

% result dir:
result_name=['result_' datestr(now, 'YYYYmmDDHHMMSS') '_predict_custom_data'];
result_dir=fullfile('../cache_data', ['test_examples_' ds_name_subfix], result_name);

% using a trained model:
run_config.trained_model_path='../model_trained/refinenet_res101_nyud.mat'; % resnet101 based refinenet
% run_config.trained_model_path='../model_trained/refinenet_res152_nyud.mat'; % resnet152 based refinenet

% provide class_info for the trained model:
ds_config.class_info=gen_class_info_nyud();

% control the size of input images to avoid excessively small or large images
run_config.input_img_short_edge_min=450;
run_config.input_img_short_edge_max=800;
%-------------------------------------------------------------------------------------------------------------------------

% specify the folder that contains testing images:
img_data_dir='../datasets/example_imgs_nyud';

% providing ground-truth mask folder for evaluation, or set to empty if not available
% if the folder of ground-truth masks are provided (an example is shown below),
% this code will perform evaluation after multi-scale prediction and fusion.
gt_mask_dir=[];

%-------------------------------------------------------------------------------------------------------------------------
% settings for accuracy evaluation after multi-scale prediction and fusion:
% the following example specifies images for testing and the ground-truth masks for evaluation:

% Notes: for the evaluation purpose, the image folder should only contain testing images, excluding any training images:
% Notes: the ground-truth masks should follow the class label definition in the variable specified above: 'ds_config.class_info'

% uncomment the following lines to enable accuracy evaluation, should replace by your folders:

% img_data_dir='../datasets/nyudv2/images_testonly';
% gt_mask_dir='../datasets/nyudv2/gt_testonly';
%-------------------------------------------------------------------------------------------------------------------------

%-------------------------------------------------------------------------------------------------------------------------
% settings for multi-scale prediction, you can consider 3 scales or 5 scales:

prediction_scales=[0.4 0.6 0.8 1 1.2]; % 5 scales
% prediction_scales=[0.6 0.8 1]; % 3 scales
% prediction_scales=[0.8]; % or only use 1 scale
%-------------------------------------------------------------------------------------------------------------------------

ds_config.img_data_dir=img_data_dir;
result_evaluate_param=[];
result_evaluate_param.gt_mask_dir=gt_mask_dir;
ds_config.result_evaluate_param=result_evaluate_param;

run_config.prediction_scales=prediction_scales;
run_config.model_name=result_name;
run_config.root_cache_dir=result_dir;

do_predict_and_evaluate(run_config, ds_config);

end

function do_predict_and_evaluate(run_config, ds_config)

run_config.gen_net_opts_fn=@gen_net_opts_model_type1;

run_config.run_evaonly=true;
ds_config.use_custom_data=true;
ds_config.use_dummy_gt=true;
run_config.use_dummy_gt=ds_config.use_dummy_gt;

ds_config.ds_name='tmp_data';
mkdir_notexist(run_config.root_cache_dir);

diary_dir=run_config.root_cache_dir;
mkdir_notexist(diary_dir);
diary(fullfile(diary_dir, 'output.txt'));
diary on

run_dir_name=fileparts(mfilename('fullpath'));
[~, run_dir_name]=fileparts(run_dir_name);
run_config.run_dir_name=run_dir_name;
run_config.run_file_name=mfilename();

ds_info=gen_dataset_info(ds_config);
my_diary_flush();

if run_config.use_gpu
	gpu_num=gpuDeviceCount;
	if gpu_num>=1
		gpuDevice(run_config.gpu_idx);
    else
        error('no gpu found!');
	end
end

train_opts=run_config.gen_net_opts_fn(run_config, ds_info.class_info);

imdb=my_gen_imdb(train_opts, ds_info);
data_norm_info=[];
data_norm_info.image_mean=128;
imdb.ref.data_norm_info=data_norm_info;

[net_config, net_exp_info]=prepare_running_model(train_opts);

prediction_scales=run_config.prediction_scales;
scale_num=length(prediction_scales);
predict_result_dirs=cell(scale_num, 1);

for s_idx=1:scale_num

	one_scale=prediction_scales(s_idx);
	one_result_dir=fullfile(run_config.root_cache_dir, sprintf('predict_result_%d', s_idx));
	predict_result_dirs{s_idx}=one_result_dir;

	fprintf('\n\n--------------------------------------------------\n\n');
	fprintf('conduct prediction using image scale: %1.2f  (current scale / total scales: %d/%d) \n', one_scale, s_idx, scale_num);

	train_opts.root_cache_dir=one_result_dir;
	train_opts.input_img_scale=one_scale;
	train_opts.eva_param=update_eva_param_mscale(train_opts.eva_param, train_opts);

	my_net_tool(train_opts, imdb, net_config, net_exp_info);

	fprintf('\n\n--------------------------------------------------\n\n');
	disp('results are saved in:');
	disp(train_opts.root_cache_dir);

	my_diary_flush();

end

if length(predict_result_dirs)>1

    fprintf('\n\n--------------------------------------------------\n\n');
    disp('fusing multiscale predictions');

    fuse_param=[];
    fuse_param.predict_result_dirs=predict_result_dirs;
    fuse_param.fuse_result_dir=fullfile(run_config.root_cache_dir, 'predict_result_final_fused');
    fuse_param.cache_fused_score_map=true;
    fuse_multiscale_results(fuse_param, ds_config.class_info);

    fprintf('\n\n--------------------------------------------------\n\n');
    disp('final fused results are saved in:');
    disp(fuse_param.fuse_result_dir);

    my_diary_flush();

    final_prediction_dir=fuse_param.fuse_result_dir;

else

    final_prediction_dir=predict_result_dirs{1};
end

% evaluate the final fused result, if ground-truth is provided
result_evaluate_param=ds_config.result_evaluate_param;
if ~isempty(result_evaluate_param) && ~isempty(result_evaluate_param.gt_mask_dir)
	result_evaluate_param.predict_result_dir=fullfile(final_prediction_dir, 'predict_result_mask');
	fprintf('\n\n--------------------------------------------------\n\n');
	disp('perform evaluation of the predicted masks');
	result_info=evaluate_predict_results(result_evaluate_param, ds_config.class_info);
    disp_evaluate_result(result_info, ds_config.class_info);

    eva_result_cached_file=fullfile(final_prediction_dir, 'eva_result_info.mat');
    fprintf('saving evaluation result to: %s\n', eva_result_cached_file);
    save(eva_result_cached_file, 'result_info');
end

my_diary_flush();
diary off

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%% t_mrdFibersLoad
%
% OLD - probably should be deleted.
%
% Initialize a fiber group, potentially cull them using old ER code, and
% save them (culled).
%
% See: t_mrdArcuatePredictions
%
% (c) Stanford VISTA Team

doculling = 0;

%% Read the dt6 file.

% We load up the dt6 and open the mrDiffusion window so that we can easily
% transform between spaces.
dt6Name = fullfile(dataDir,'dti40','dt6.mat');
[dtiF, dtiH] = mrDiffusion('off',dt6Name);

% get the xForm from image space 3 acpc space
xForm.img2acpc = dtiGet(dtiH,'img 2 acpc xform');

% get the xForm from acpc space to image space
xForm.acpc2img = dtiGet(dtiH,'acpc2img xform');

% Set the colors of the three fiber groups.
cList = {[200,50 ,50],[50,200,50],[50,50,200]};

%% Load the first STT group

fgName{1} = fullfile(mrvDataRootPath,'diffusion','fiberPrediction','fibers','R_Arcuate_Box_STT.mat');
temp = dtiReadFibers(fgName{1}); % Fiber coordinates in acpc space

if doculling
 % Remove redundant fibers from the group:
 fgAcpc(1) = dtiCullFibers(temp, dt6Name);
else
 fgAcpc(1) = temp;
end

%% read ConTrac fiber group
fgName{2} = fullfile(mrvDataRootPath,'diffusion','fiberPrediction','fibers','R_Arcuate_ctr_clean.mat');
temp = dtiReadFibers(fgName{2}); % Fiber coordinates in acpc space

if doculling
 % Remove redundant fibers from the group:
 fgAcpc(2) = dtiCullFibers(temp, dt6Name);
else
  fgAcpc(2) = temp;
end

%% Merge them first fiber group
fgName{3} = 'R_Arcuate_Box_merged_stt_ctr-culled.mat';
fgAcpc(3) = dtiMergeFiberGroups(fgAcpc(1),fgAcpc(2),fgName{3}); % Fiber coordinates in acpc space

%% Spatial transforms into image space

% Fibers are stored in ACPC space.
% In this space, they are in mm
%    fLengths = fgGet(fg,'fiber lengths');
%     \mrvNewGraphWin; hist(fLengths,100)
fgImg = struct(fgAcpc);
for ii = 1:length(fgAcpc)
 dtiH = dtiSet(dtiH,'add fiber group',fgAcpc(ii));

 % Have a look at the fibers.  Visualizations are in ACPC space.
 % fgMesh(fgAcpc{i},dtiH);

 % Create a version of the fiber group in image space for computations
 % There is a problem:  See dtiXformFiberCoords.  But since this script
 % should probably go away, I am not fixing it now.
 fgImg(ii) = dtiXformFiberCoords(fgAcpc(ii),xForm.acpc2img);

 fgImg(ii).colorRgb = cList{ii};

end

%% Save

% Becasuse dtiCullFibers takes a while for the large group, we save these
% out.
culledFiberGroups = ...
    fullfile(mrvDataRootPath,'diffusion','fiberPrediction','fibers','culledFiberGroups.mat');
save(culledFiberGroups,'fgImg');

%% To view the saved fibers ...

% guidata(dtiF,dtiH);
% figure(dtiF)

% To load the fiber groups:
%   load(culledFiberGroups,'fgImg')

%% End

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [] = rotation(st);

if nargin < 1,
   st= 1;
end;

fig = gcf;

ax = gca;

vv = get(ax,'view');

az = vv(1);
el = vv(2);

for azi = az:-abs(st):az-360,

   set(ax,'view',[azi el]);
   figure(fig);
   drawnow;

end;

repo: ekfmonoslam


text:
classdef TestRectify3Collinear
    %TestRectify3Collinear

    properties (Constant)
        fields = {'R1', 'R2', 'R3', 'P1', 'P2', 'P3', 'Q', ...
            'roi1', 'roi2', 'ratio'};
    end

    methods (Static)
        function test_1
            imgSz = [640 480];
            cam1 = eye(3);
            cam2 = eye(3);
            cam3 = eye(3);
            distort1 = zeros(1,5);
            distort2 = zeros(1,5);
            distort3 = zeros(1,5);
            R12 = [cv.getRotationMatrix2D(imgSz./2, 0, 1); 0 0 1];
            R13 = [cv.getRotationMatrix2D(imgSz./2, 0, 1); 0 0 1];
            T12 = [1; 0; 0];
            T13 = [1; 0; 0];
            S = cv.rectify3Collinear(cam1, distort1, cam2, distort2, ...
                cam3, distort3, imgSz, R12, T12, R13, T13, ...
                'ZeroDisparity',false, 'Alpha',-1, 'NewImageSize',imgSz./2);
            validateattributes(S, {'struct'}, {'scalar'});
            assert(all(ismember(TestRectify3Collinear.fields, fieldnames(S))));
            validateattributes(S.R1, {'numeric'}, {'size',[3 3]});
            validateattributes(S.R2, {'numeric'}, {'size',[3 3]});
            validateattributes(S.R3, {'numeric'}, {'size',[3 3]});
            validateattributes(S.P1, {'numeric'}, {'size',[3 4]});
            validateattributes(S.P2, {'numeric'}, {'size',[3 4]});
            validateattributes(S.P3, {'numeric'}, {'size',[3 4]});
            validateattributes(S.Q, {'numeric'}, {'size',[4 4]});
            validateattributes(S.roi1, {'numeric'}, ...
                {'vector', 'numel',4, 'integer'});
            validateattributes(S.roi2, {'numeric'}, ...
                {'vector', 'numel',4, 'integer'});
            validateattributes(S.ratio, {'numeric'}, {'scalar'});
        end

        function test_error_argnum
            try
                cv.rectify3Collinear();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function outdata=specscope(indata)
% record and plot audio spectrogram
%
% Usage: outdata=specscope(indata)
%
%  Input: indata (optional)
%    Displays a recorded piece of data, if an argument is passed
%    Otherwise displays audio data from an attached microphone
%
%  Output: outdata (optional)
%    If present, will return up to 10 minutes
%    of captured audio data.
%
%  Note: Parameters such as sampling frequency, number of tapers
%    and display refresh rate may be set below if desired.
%    You can also acquire data from a national instruments card.
%

close all;

%=======================================================================
% Check toolboxes
%=======================================================================

% Check for toolboxes
if not(exist('analoginput','file'));
    fprintf('You need to install the DAQ toolbox first\n');
    return
end
if not(exist('mtspecgramc','file'));
    fprintf('You need to install the Chronux toolbox first from http://chronux.org/\n');
    return
end

%=======================================================================
% Set parameters
%=======================================================================

global acq;

% Set defaults
acq.params.Fs = 44100;
acq.pause = 0;
acq.skips = 0;
acq.stop = 0;
acq.restart = 0;
acq.plot_frequency = 10;
acq.samples_acquired = 0;
acq.spectra = [];
acq.times = [];
defaults
audio_instr;
fig=create_ui;

%=======================================================================
% Check arguments, start DAQ
%=======================================================================

if nargout
   % save up to ten minutes data, preallocated...
    fprintf('Pre-allocating memory for data save - please be patient.\n');
   outdata=zeros( (acq.params.Fs * 60 * 10), 1 );
end

if nargin == 1;
    acq.indata = indata;
    acq.live = 0;
else
    % Create and set up and start analog input daq
    input=1;
    if input==1;
        acq.ai = analoginput('winsound');
        addchannel( acq.ai, 1 );
    else
        acq.ai = analoginput('nidaq', 1);
        addchannel(acq.ai, 0);
        set(acq.ai,'InputType','SingleEnded')
        set(acq.ai,'TransferMode','Interrupts')
        set(acq.ai,'TriggerType','Manual');
    end
    set( acq.ai, 'SampleRate', acq.params.Fs )
    acq.params.Fs = get( acq.ai, 'SampleRate' );
    set( acq.ai, 'SamplesPerTrigger', inf )
    start(acq.ai)
    acq.live = 1;

    if input==2;
        trigger(acq.ai);
    end
end

acq.samples_per_frame = acq.params.Fs / acq.plot_frequency;

%=======================================================================
% The scope main loop
%=======================================================================

acq.t0=clock;
acq.tn=clock;

% Loop over frames to acquire and display
while 1;

    % Check for quit signal
    if acq.stop;
        break;
    end

    % Calculate times
    calctime = acq.samples_acquired / acq.params.Fs;
    acq.samples_acquired = acq.samples_acquired + acq.samples_per_frame;
    acq.t1 = clock;
    elapsed = etime(acq.t1,acq.t0);

    % Get a small snippet of data
    if ( acq.live )
       data = getdata( acq.ai, acq.samples_per_frame );
    else
      while elapsed < acq.samples_acquired / acq.params.Fs
        pause( acq.samples_acquired / (acq.params.Fs) - elapsed );
        acq.t1=clock;
        elapsed = etime(acq.t1,acq.t0);
      end
      if acq.samples_acquired + 2 * acq.samples_per_frame >= length( acq.indata )
        acq.stop=1;
      end
      data = acq.indata(floor(acq.samples_acquired+1):floor(acq.samples_per_frame+acq.samples_acquired));
    end

    if nargout
        outdata(floor(acq.samples_acquired+1):floor(acq.samples_acquired+length(data))) = data(:);
    end

    if acq.restart;
       acq.restart = 0;
       acq.spectra = [];
       acq.times = [];
    end

    % Calculate spectrogram of data snippet
    if acq.deriv
      [s, t, f] = mtspecgramc(diff(data), acq.moving_window, acq.params );
    else
      [s, t, f] = mtspecgramc(data, acq.moving_window, acq.params );
    end

    % Add new spectra to that already calculated
    acq.times = [acq.times t+calctime];
    if acq.log
        acq.spectra = [acq.spectra log(s')];
    else
        acq.spectra = [acq.spectra s'];
    end

    % Remove old spectra once window reaches desired size
    while acq.times(1,end) - acq.times(1,1) > acq.display_size;

        % Ring buffer!
        y = length(t);
        acq.times(:,1:y) = [];
        acq.spectra(:,1:y) = [];

    end

    % Only plot if display is keeping up with real time and not paused
    show_plot=1;
    if nargin==0
       if get(acq.ai, 'SamplesAvailable' ) > 10 * acq.samples_per_frame && acq.pause==0
	  show_plot=0;
       end
    else
      if elapsed > calctime + 0.5
	show_plot=0;
      end
    end

    if acq.pause
       show_plot=0;
    end
    if show_plot

        if  acq.bgsub
            acq.mean_spectra = mean( acq.spectra, 2 );
        end

        % Normalize until full screen passes by if requested
        if acq.normalize>=1;
            if acq.normalize==1
                acq.tn=clock;
                acq.normalize=2;
            end
            if etime(clock,acq.tn)>1.25*acq.display_size
                acq.normalize=0;
            end
            mins = min(min(acq.spectra));
            maxs = max(max(acq.spectra));
        end

        % Scale the spectra based upon current offset and scale
        if acq.bgsub
           scaled_spectra = acq.offset + ( acq.scale ) * ( acq.spectra - repmat( acq.mean_spectra, [1,size(acq.spectra,2)]) ) / ( maxs - mins );
        else
            scaled_spectra = acq.offset + acq.scale * ( acq.spectra - mins ) / ( maxs - mins );
        end

        % Draw the image to the display
        image( acq.times, f, scaled_spectra ); axis xy;
        drawnow;

    else
        % Keep track of skipped displays
        acq.skips = acq.skips + 1;
    end

end

%=======================================================================
% Clean up
%=======================================================================

acq.t1=clock;
elapsed = etime(acq.t1,acq.t0);
fprintf( 'Elapsed time %f seconds\n', elapsed );

% Warn if many skips were encountered
if acq.skips > 5;
    fprintf( '\nWARNING:\nThis program skipped plotting %d times to keep pace.\n', acq.skips )
    fprintf( 'Run again without keyboard interaction or changing the figure size.\n' )
    fprintf( 'If this message reappears you should reduce the plot frequency parameter.\n\n' );
end

% Clean up the analoginput object
if acq.live
  stop(acq.ai);delete( acq.ai );clear acq.ai;
end

% Clean up the figure
delete(fig);
delete(gcf);

if nargout
   % save up to ten minutes data, preallocated...
    fprintf('Saving output data\n');
   outdata=outdata(1:floor(acq.samples_acquired));
end

return;

%
%
%=======================================================================
% Functions called
%=======================================================================
%
%

%=======================================================================
% Handle Keypresses
%=======================================================================

% Handle figure window keypress events
function keypress(varargin)

global acq;
keypressed=get(gcf,'CurrentCharacter');

% ignore raw control, shift, alt keys
if keypressed;

    % Save current frame as gif
    if strcmp( keypressed, 'g');
       saveas( acq.fig, sprintf( 'frame%d.png',acq.times(length(acq.times)) ) )
    end

    % Offset changes
    increment=1;
    if strcmp( keypressed, 'l');
	acq.offset = acq.offset - increment;
    elseif strcmp( keypressed, 'o');
        acq.offset = acq.offset + increment;

    % Scale changes
    elseif strcmp( keypressed, 'x');
        acq.scale = acq.scale - increment;
    elseif strcmp( keypressed, 's');
        acq.scale = acq.scale + increment;

    % Reset defaults
    elseif strcmp( keypressed, 'd');
        defaults
	acq.restart=1;
    % Normalize spectra
    elseif strcmp( keypressed, 'n');
	   request_normalize

    % Quit
    elseif strcmp( keypressed, 'q');
	request_quit

    % Pause
    elseif strcmp( keypressed, 'p');
	   request_pause

    % Help
    elseif strcmp( keypressed, 'h');
        audio_instr

    % Change colormaps for 0-9,a-c
    elseif strcmp( keypressed, '0' );
        colormap( 'jet' );
    elseif strcmp( keypressed, '1' );
        colormap( 'bone' );
    elseif strcmp( keypressed, '2' );
        colormap( 'colorcube' );
    elseif strcmp( keypressed, '3' );
        colormap( 'cool' );
    elseif strcmp( keypressed, '4' );
        colormap( 'copper' );
    elseif strcmp( keypressed, '5' );
        colormap( 'gray' );
    elseif strcmp( keypressed, '6' );
        colormap( 'hot' );
    elseif strcmp( keypressed, '7' );
        colormap( 'hsv' );
    elseif strcmp( keypressed, '8' );
        colormap( 'autumn' );
    elseif strcmp( keypressed, '9' );
        colormap( 'pink' );
    elseif strcmp( keypressed, 'a' );
        colormap( 'spring' );
    elseif strcmp( keypressed, 'b' );
        colormap( 'summer' );
    elseif strcmp( keypressed, 'c' );
        colormap( 'winter' );
    end

    update_display

end
return

%=======================================================================
% Defaults
%=======================================================================

% Reset defaults
function defaults()
  global acq;
  acq.params.raw_tapers = [2 3];
  acq.moving_window = [0.02 0.02];
  acq.params.tapers=dpsschk(acq.params.raw_tapers,round(acq.params.Fs*acq.moving_window(1)),acq.params.Fs);
  acq.offset = 0;
  acq.scale = 64;
  acq.display_size = 3;
  acq.params.fpass = [50 8000];
  acq.deriv=1;
  acq.log=1;
  acq.bgsub = 1;
  acq.params.pad= 0;
  acq.normalize = 2;
return

function update_display()
global acq;
    set(acq.tapers_ui,'String',sprintf( '%.0f %.0f', acq.params.raw_tapers(1), acq.params.raw_tapers(2) ));
    set(acq.window_ui,'String',sprintf( '%.2f %.2f', acq.moving_window(1), acq.moving_window(2) ));
    set(acq.offset_ui,'String',sprintf( '%d', acq.offset ));
    set(acq.scale_ui,'String',sprintf( '%d', acq.scale ));
    set(acq.display_size_ui,'String',sprintf( '%.1f', acq.display_size ));
    set(acq.frequency_ui,'String',sprintf( '%.1f %.1f', acq.params.fpass(1), acq.params.fpass(2)  ))
    set(acq.derivative_ui,'Value',acq.deriv);
    set(acq.log_ui,'Value',acq.log);
    set(acq.bgsub_ui,'Value',acq.bgsub);
    return

%=======================================================================
% Update ui controls
%=======================================================================

function request_quit(varargin)
	 global acq;
	 acq.stop=1;
return

function request_pause(varargin)
	 global acq;
	 acq.pause = not( acq.pause );
return

function request_normalize(varargin)
	global acq;
        acq.normalize = 2;
return

function update_defaults(varargin)
  global acq;
  defaults
  update_display
  acq.restart=1;
return

function update_tapers(varargin)
	 global acq;
     acq.params.raw_tapers = sscanf(get( gco, 'string' ),'%f %d')';
     acq.params.tapers=dpsschk(acq.params.raw_tapers,round(acq.params.Fs*acq.moving_window(1)),acq.params.Fs); % check tapers
return

function update_window(varargin)
	 global acq;
	 acq.moving_window = sscanf(get( gco, 'string' ),'%f %f');
         acq.params.tapers=dpsschk(acq.params.raw_tapers,round(acq.params.Fs*acq.moving_window(1)),acq.params.Fs);
	 acq.restart = 1;
return

function update_offset(varargin)
	 global acq;
	 acq.offset = sscanf(get( gco, 'string' ),'%f');
	 return

function update_scale(varargin)
	 global acq;
	 acq.scale = sscanf(get( gco, 'string' ),'%f');
	 return

function update_display_size(varargin)
	 global acq;
	 acq.display_size = sscanf(get( gco, 'string' ),'%f');
	 return

function update_fpass(varargin)
	 global acq;
	 acq.params.fpass = sscanf(get( gco, 'string' ),'%f %f');
         acq.restart = 1;
	 return

function update_deriv(varargin)
	 global acq;
     acq.deriv=get( gco, 'Value' );
     acq.normalize=1;
     return

function update_log(varargin)
	 global acq;
     acq.log=get( gco, 'Value' );
     acq.normalize=1;
     return

function update_bgsub(varargin)
	 global acq;
     acq.bgsub=get( gco, 'Value' );
     return

%=======================================================================
% UI display
%=======================================================================

function fig=create_ui()
    global acq;

    bgcolor = [1 1 1]; % .7 .7 .7
    % ===Create main figure==========================
    fig = figure('Position',centerfig(800,600),...
        'NumberTitle','off',...
        'Name','Real-time spectrogram',...
        'doublebuffer','on',...
        'HandleVisibility','on',...
	'Renderer', 'openGL', ...
	'KeyPressFcn', @keypress, ...
        'Color',bgcolor);

    acq.fig = fig;
    offset = 80;
    % ===text==========
    uicontrol(gcf,'Style','text',...
        'String', 'tapers',...
        'Position',[offset+225 20 45 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'moving win',...
        'Position',[offset+300 20 70 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'offset',...
        'Position',[offset+375 20 30 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'scale',...
        'Position',[offset+410 20 30 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 't axis',...
        'Position',[offset+445 20 30 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'f axis',...
        'Position',[offset+480 20 40 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'deriv',...
        'Position',[offset+550 20 35 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'log',...
        'Position',[offset+580 20 35 20],...
        'BackgroundColor',bgcolor);
    uicontrol(gcf,'Style','text',...
        'String', 'bgsub',...
        'Position',[offset+610 20 35 20],...
        'BackgroundColor',bgcolor);

    % ===The quit button===============================
    uicontrol('Style','pushbutton',...
        'Position',[offset+5 5 45 20],...
        'String','Quit',...
        'Interruptible','off',...
        'BusyAction','cancel',...
        'Callback',@request_quit);

    % ===The pause button===============================
    uicontrol('Style','pushbutton',...
        'Position',[offset+55 5 45 20],...
        'String','Pause',...
        'Interruptible','off',...
        'BusyAction','cancel',...
        'Callback',@request_pause);

    % ===The defaults button===============================
    uicontrol('Style','pushbutton',...
        'Position',[offset+105 5 50 20],...
        'String','Defaults',...
        'Interruptible','off',...
        'BusyAction','cancel',...
        'Callback',@update_defaults);

    % ===The normalize button===============================
    uicontrol('Style','pushbutton',...
        'Position',[offset+160 5 60 20],...
        'String','Normalize',...
        'Interruptible','off',...
        'BusyAction','cancel',...
        'Callback',@request_normalize );

    % ===Tapers============================================
    acq.tapers_ui = uicontrol(gcf,'Style','edit',...
        'String', sprintf( '%.0f %.0f', acq.params.raw_tapers(1), acq.params.raw_tapers(2) ),...
        'Position',[offset+225 5 70 20],...
        'CallBack', @update_tapers);

    % ===Window============================================
    acq.window_ui=uicontrol(gcf,'Style','edit',...
        'String', sprintf( '%.2f %.2f', acq.moving_window(1), acq.moving_window(2) ),...
        'Position',[offset+300 5 70 20],...
        'CallBack', @update_window);

    % ===Offset============================================
    acq.offset_ui = uicontrol(gcf,'Style','edit',...
        'String', sprintf( '%d', acq.offset ),...
        'Position',[offset+375 5 30 20],...
        'CallBack', @update_offset);

    % ===Scale============================================
    acq.scale_ui = uicontrol(gcf,'Style','edit',...
        'String', sprintf( '%d', acq.scale ),...
        'Position',[offset+410 5 30 20],...
        'CallBack', @update_scale);

    % ===display size======================================
    acq.display_size_ui = uicontrol(gcf,'Style','edit',...
        'String', sprintf( '%.1f', acq.display_size ),...
        'Position',[offset+445 5 30 20],...
        'CallBack', @update_display_size);

    % ===frequency axis=====================================
    acq.frequency_ui = uicontrol(gcf,'Style','edit',...
        'String', sprintf( '%.1f %.1f', acq.params.fpass(1), acq.params.fpass(2)  ),...
        'Position',[offset+480 5 80 20],...
        'CallBack', @update_fpass);

    % ===derivative=====================================
    acq.derivative_ui = uicontrol(gcf,'Style','checkbox',...
        'Value',acq.deriv,...
        'Position',[offset+565 5 20 20],...
        'CallBack', @update_deriv);

    % ===log=====================================
    acq.log_ui = uicontrol(gcf,'Style','checkbox',...
        'Value',acq.log,...
        'Position',[offset+590 5 20 20],...
        'CallBack', @update_log);

    % ===bgsub=====================================
    acq.bgsub_ui = uicontrol(gcf,'Style','checkbox',...
        'Value',acq.bgsub,...
        'Position',[offset+615 5 20 20],...
        'CallBack', @update_bgsub);

return

%=======================================================================
% Assorted functions
%=======================================================================

function pos = centerfig(width,height)
% Find the screen size in pixels
screen_s = get(0,'ScreenSize');
pos = [screen_s(3)/2 - width/2, screen_s(4)/2 - height/2, width, height];
return

function audio_instr()
% Show instructions

  fprintf('INSTRUCTIONS:\n');
  fprintf('Click on figure window first to activate controls.\n')
  fprintf('Adjust tapers, windows, scales, offsets and axes using the gui\n');
  fprintf('The deriv checkbox toggles derivative of the data\n');
  fprintf('The log checkbox toggles a log of the spectrum\n');
  fprintf('Press d or use defaults button to reset most parameters to defaults.\n')
  fprintf('Press n or use normalize button to normalize spectra based upon values in current display.\n')
  fprintf('Press 0-9,a-c to choose a colormap (default 0).\n')
  fprintf('Press p to pause and unpause display.\n')
  fprintf('Press o and l to adjust offset, or use offset textbox on gui.\n');
  fprintf('Press s and x to adjust scale, or use scale textbox on gui.\n');
  fprintf('Press h for this message.\n')
  fprintf('Press q to quit, or use quit button on gui.\n\n')

return

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function geometry_test022 ( )

%*****************************************************************************80
%
%% TEST022 tests DIRECTION_UNIFORM_3D, GET_SEED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST022\n' );
  fprintf ( 1, '  DIRECTION_UNIFORM_3D picks a random direction vector.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ vran, seed ] = direction_uniform_3d ( seed );
    fprintf ( 1, '  %10f  %10f  %10f\n', vran(1:dim_num) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function loadlibltfat(varargin)
global libltfat_intptrstr;

definput.keyvals.lib='libltfat.so';
definput.flags.phase={'load','reload','recompile'};
definput.flags.comptarget={'fulloptim','release','debug'};
definput.flags.verbosity={'quiet','verbose'};
definput.flags.corcpp={'c','cpp'};
definput.keyvals.compiler = [];

[flags,kv,lib]=ltfatarghelper({'lib'},definput,varargin);

[~,libname]=fileparts(lib);
currdir = fileparts(mfilename('fullpath'));
libltfatpath = [currdir, filesep, '..', filesep,'..',filesep,'..',filesep,'..',filesep];
libpath = [libltfatpath, filesep,'build',filesep,lib];
iscompiled = exist(libpath,'file');

makecmd = ['make -C ',libltfatpath];
if ~iscompiled || flags.do_recompile
    [status,result] = system([makecmd, ' clean']);
    if status ~=0, error(result);  end

    makecmd = [makecmd, ' MODULE=libltfat'];
    makecmd = [makecmd, ' munit -j12'];
    makecmd = [makecmd, ' MATLABROOT=', matlabroot];
    makecmd = [makecmd, sprintf(' COMPTARGET=%s',flags.comptarget)];

    if flags.do_cpp
        makecmd = [makecmd, ' USECPP=1'];
    end

    if kv.compiler
        makecmd = [makecmd, sprintf(' CC=%s',kv.compiler)];
    end

    if flags.do_verbose
        disp(makecmd);
        system(makecmd);
    else
        [status,result] = system(makecmd);
        if status ~=0, error(result);  end
    end
end

if libisloaded(libname)
    if flags.do_reload || flags.do_recompile
        unloadlibrary(libname);
    else
        error('%s: libltfat is already loaded. Use ''reload'' to force reload.',upper(mfilename));
    end
end

warning('off');
headerpath = [libltfatpath,'build',filesep,'ltfat.h'];
loadlibrary(libpath,headerpath,'mfilename','libltfatprotofile.m');
warning('on');

intbitsize = 8*calllib('libltfat','ltfat_int_size');
libltfat_intptrstr = sprintf('int%dPtr',intbitsize);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [theta] =  tapas_sem_multiv_init_theta(data, model, inference)
%% Obtain an initial sample with positive likelihood.
%
% aponteeduardo@gmail.com
% copyright (C) 2017
%

ptheta = model.graph{1}.htheta.model;

nc = size(model.graph{1}.htheta.T, 2);
ns = size(data, 1);
mu = model.graph{4}.htheta.y.mu;
nb = numel(mu);
np = size(ptheta.jm, 2);

ty = ptheta.x * model.graph{4}.htheta.y.mu;
ty = reshape(ty', numel(ty), 1);
theta = repmat(mat2cell(ty, np * ones(ns, 1), 1), 1, 1);

llh = tapas_sem_multiv_llh(data, struct('y', {theta}), model.graph{1}.htheta);
failing = find(llh == -inf);

tolerance = 1000;
while numel(failing) && tolerance
    for i = failing'
        theta{i} = sample_gaussian(ptheta);
    end
    llh = tapas_sem_multiv_llh(data, struct('y', {theta}), ...
        model.graph{1}.htheta);
    failing = find(llh == -inf);
    tolerance = tolerance - 1;
end

if numel(failing)
    error('tapas:sem:multiv:init', ...
    'It was not possible to initialize a sample with positive likelihood');
end

theta = repmat(theta, 1, nc);

end

function [theta] = sample_gaussian(ptheta)
%% Sample from a Gaussian prior.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

np = size(ptheta.jm, 2);

sample_pars = logical(sum(ptheta.jm, 2));

if size(ptheta.pm, 2) ~= np
    pm = diag(ptheta.pm);
else
    pm = ptheta.pm;
end

lt = chol(pm);
theta = lt \ ptheta.jm * randn(np, 1);
theta = ptheta.sm' * theta;

end

author: translationalneuromodeling

repo: tapas


text:
function mtrPlotRankCorrSummary(summaryFilename)

summary = load(summaryFilename);
paramNames = {'kLength','kSmooth','kMidSD'};
midP = [-2 18 0.175];
ccGrid = mtrFilenames2Paramlist(summary.pdbIDFiles,paramNames);

for pp = 1:length(paramNames)
    strThreshVec = {};
    parVecs = [];
    corrVecs = [];
    ovrVecs = [];
    for tt = 1:length(summary.threshVec)
        ccGrid(:,4) = summary.rhoS(:,tt);
        ccGrid(:,5) = summary.overlapVec(:,tt);
        indGrid = ones(size(ccGrid,1),1);
        for ss = 1:length(paramNames)
            if ss ~= pp
                indGrid = indGrid & ccGrid(:,ss) == midP(ss);
            end
        end
        subGrid = ccGrid(indGrid,:);
        [foo, sortI] = sort(subGrid(:,pp));
        parVecs(:,tt) = subGrid(sortI(:),pp);
        corrVecs(:,tt) = subGrid(sortI(:),4);
        ovrVecs(:,tt) = subGrid(sortI(:),5);
        strThreshVec{tt} = ['Top ' num2str(summary.threshVec(tt))];
    end
    figure; plot(parVecs,corrVecs);
    legend(strThreshVec);
    xlabel([paramNames{pp} ' parameter']);
    ylabel('Spearman rank correlation');
    figure; plot(parVecs,ovrVecs);
    legend(strThreshVec);
    xlabel([paramNames{pp} ' parameter']);
    ylabel('Overlap');
end


text:
function [Nm] = J2Nm(J)
% Convert energy or work from joules to newton-meters.
% Chad A. Greene 2012
Nm = J;
% You say tomatermort, I say tomato.

author: Sable

repo: mcbench-benchmarks


text:
function varargout = lvd_editLaunchVehicle(varargin)
% LVD_EDITLAUNCHVEHICLE MATLAB code for lvd_editLaunchVehicle.fig
%      LVD_EDITLAUNCHVEHICLE, by itself, creates a new LVD_EDITLAUNCHVEHICLE or raises the existing
%      singleton*.
%
%      H = LVD_EDITLAUNCHVEHICLE returns the handle to a new LVD_EDITLAUNCHVEHICLE or the handle to
%      the existing singleton*.
%
%      LVD_EDITLAUNCHVEHICLE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITLAUNCHVEHICLE.M with the given input arguments.
%
%      LVD_EDITLAUNCHVEHICLE('Property','Value',...) creates a new LVD_EDITLAUNCHVEHICLE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_editLaunchVehicle_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_editLaunchVehicle_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_editLaunchVehicle

% Last Modified by GUIDE v2.5 31-Jul-2020 17:15:00

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_editLaunchVehicle_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_editLaunchVehicle_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_editLaunchVehicle is made visible.
function lvd_editLaunchVehicle_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to lvd_editLaunchVehicle (see VARARGIN)

% Choose default command line output for lvd_editLaunchVehicle
    handles.output = hObject;

    lvdData = varargin{1};
    setappdata(hObject, 'lvdData', lvdData);

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_editLaunchVehicle wait for user response (see UIRESUME)
    uiwait(handles.lvd_editLaunchVehicle);

function populateGUI(handles, lvdData)
    jScrollPane = findjobj(handles.lvSummaryText);
    jViewPort = jScrollPane.getViewport;
    jEditbox = jViewPort.getComponent(0);
    jEditbox.setWrapping(false);
    jScrollPane.setHorizontalScrollBarPolicy(32);
    jEditbox.setEditable(false);
    jEditbox.revalidate();
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

function setSummText(hSummText, lv)
    set(hSummText, 'String', lv.getLvSummaryStr());

% --- Outputs from this function are returned to the command line.
function varargout = lvd_editLaunchVehicle_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
        lvdData.initStateModel.clearAllTankStatesAndRegenerate();
        lvdData.initStateModel.clearDuplicateEngineStates();

        close(handles.lvd_editLaunchVehicle);
    end

function lvSummaryText_Callback(hObject, eventdata, handles)
% hObject    handle to lvSummaryText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of lvSummaryText as text
%        str2double(get(hObject,'String')) returns contents of lvSummaryText as a double

% --- Executes during object creation, after setting all properties.
function lvSummaryText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lvSummaryText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in editStagesButton.
function editStagesButton_Callback(hObject, eventdata, handles)
% hObject    handle to editStagesButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditStagesGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editEnginesButton.
function editEnginesButton_Callback(hObject, eventdata, handles)
% hObject    handle to editEnginesButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditEnginesGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editTanksButton.
function editTanksButton_Callback(hObject, eventdata, handles)
% hObject    handle to editTanksButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditTanksGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editEngineToTankConnsButton.
function editEngineToTankConnsButton_Callback(hObject, eventdata, handles)
% hObject    handle to editEngineToTankConnsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditEngineToTanksConnsGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in closeButton.
function closeButton_Callback(hObject, eventdata, handles)
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    uiresume(handles.lvd_editLaunchVehicle);

% --- Executes on key press with focus on lvd_editLaunchVehicle or any of its controls.
function lvd_editLaunchVehicle_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_editLaunchVehicle (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            uiresume(handles.lvd_editLaunchVehicle);
        case 'enter'
            uiresume(handles.lvd_editLaunchVehicle);
        case 'escape'
            uiresume(handles.lvd_editLaunchVehicle);
    end

% --- Executes on button press in editTankToTankConnsButton.
function editTankToTankConnsButton_Callback(hObject, eventdata, handles)
% hObject    handle to editTankToTankConnsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditTankToTankConnsGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editTankFluidTypesButton.
function editTankFluidTypesButton_Callback(hObject, eventdata, handles)
% hObject    handle to editTankFluidTypesButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditTankFluidTypesGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editEpsSrcsButton.
function editEpsSrcsButton_Callback(hObject, eventdata, handles)
% hObject    handle to editEpsSrcsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditElectricalPowerSrcsGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editEpsStorageButton.
function editEpsStorageButton_Callback(hObject, eventdata, handles)
% hObject    handle to editEpsStorageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditElectricalPowerStoragesGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

% --- Executes on button press in editEpsSinksButton.
function editEpsSinksButton_Callback(hObject, eventdata, handles)
% hObject    handle to editEpsSinksButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_editLaunchVehicle, 'lvdData');
    lvd_EditElectricalPowerSinksGUI(lvdData);
    setSummText(handles.lvSummaryText, lvdData.launchVehicle);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function X=pathSolutionLeast(A, y, z, opts)
%
%% Fuction pathSolution:
%      Solving the pathwise solutions
%
%% Input & Output parameters
%  See the description of the related functions
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified 2 August 2009.
%
% Related functions:
%  sll_opts, initFactor,
%  eppVector, eppMatrix, eplb,
%
%  LeastR, LeastC,
%  nnLeastR, nnLeastC
%  glLeastR,  mtLeastR, mcLeastR
%%

switch(opts.fName)
    case 'LeastR'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(-z);   % sort z in a decresing order
        z_value=-z_value;            % z_value in a decreasing order
        n=size(A,2);                 % the dimensionality of the data
        X=zeros(n,z_num);            % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=LeastR(A, y, z_value(1), opts);

        X(:,z_ind(1))=x;             % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function LeastR
            [x, funVal]=LeastR(A, y, z_value(i), opts);

            X(:,z_ind(i))=x;         % store the solution
        end

    case 'LeastC'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(z);    % sort z in an ascending order
        n=size(A,2);                 % the dimensionality of the data
        X=zeros(n,z_num);            % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=LeastC(A, y, z_value(1), opts);

        X(:,z_ind(1))=x;             % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function LeastC
            [x, funVal]=LeastC(A, y, z_value(i), opts);

            X(:,z_ind(i))=x;         % store the solution
        end

    case 'glLeastR'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(-z);   % sort z in a decresing order
        z_value=-z_value;            % z_value in a decreasing order
        n=size(A,2);                 % the dimensionality of the data
        X=zeros(n,z_num);            % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=glLeastR(A, y, z_value(1), opts);

        X(:,z_ind(1))=x;             % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function glLeastR
            [x, funVal]=glLeastR(A, y, z_value(i), opts);

            X(:,z_ind(i))=x;         % store the solution
        end

    case 'mtLeastR'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(-z);   % sort z in a decresing order
        z_value=-z_value;            % z_value in a decreasing order
        n=size(A,2);                 % the dimensionality of the data
        k=length(opts.ind)-1;        % the number of tasks
        X=zeros(n,k,z_num);          % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=mtLeastR(A, y, z_value(1), opts);

        X(:,:,z_ind(1))=x;           % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function mtLeastR
            [x, funVal]=mtLeastR(A, y, z_value(i), opts);

            X(:,:, z_ind(i))=x;      % store the solution
        end

    case 'mcLeastR'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(-z);   % sort z in a decresing order
        z_value=-z_value;            % z_value in a decreasing order
        n=size(A,2);                 % the dimensionality of the data
        k=size(y,2);                 % the number of tasks
        X=zeros(n,k,z_num);          % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=mcLeastR(A, y, z_value(1), opts);

        X(:,:,z_ind(1))=x;           % store the solution

        % set .init for warm start
        opts.init=1;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function mcLeastR
            [x, funVal]=mcLeastR(A, y, z_value(i), opts);

            X(:,:, z_ind(i))=x;      % store the solution
        end

    case 'nnLeastR'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(-z);   % sort z in a decresing order
        z_value=-z_value;            % z_value in a decreasing order
        n=size(A,2);                 % the dimensionality of the data
        X=zeros(n,z_num);            % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=nnLeastR(A, y, z_value(1), opts);

        X(:,z_ind(1))=x;             % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function LeastR
            [x, funVal]=nnLeastR(A, y, z_value(i), opts);

            X(:,z_ind(i))=x;         % store the solution
        end

    case 'nnLeastC'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(z);    % sort z in an ascending order
        n=size(A,2);                 % the dimensionality of the data
        X=zeros(n,z_num);            % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=nnLeastC(A, y, z_value(1), opts);

        X(:,z_ind(1))=x;             % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function LeastC
            [x, funVal]=nnLeastC(A, y, z_value(i), opts);

            X(:,z_ind(i))=x;         % store the solution
        end

    case 'mtLeastC'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(z);    % sort z in an ascending order
        n=size(A,2);                 % the dimensionality of the data
        k=length(opts.ind)-1;        % the number of tasks
        X=zeros(n,k,z_num);          % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=mtLeastC(A, y, z_value(1), opts);

        X(:,:,z_ind(1))=x;           % store the solution

        % set .init for warm start
        opts.init=0;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function mtLeastC
            [x, funVal]=mtLeastC(A, y, z_value(i), opts);

            X(:,:, z_ind(i))=x;      % store the solution
        end

    case 'mcLeastC'
        z_num=length(z);             % the number of parameters
        [z_value, z_ind]=sort(z);    % sort z in an ascending order
        n=size(A,2);                 % the dimensionality of the data
        k=size(y,2);                 % the number of tasks
        X=zeros(n,k,z_num);          % set the size of output X

        % run the code to compute the first solution
        [x, funVal]=mcLeastC(A, y, z_value(1), opts);

        X(:,:,z_ind(1))=x;           % store the solution

        % set .init for warm start
        opts.init=1;                 % using .initFactor

        for i=2:z_num
            opts.x0=x;               % warm-start

            % run the function mcLeastC
            [x, funVal]=mcLeastC(A, y, z_value(i), opts);

            X(:,:, z_ind(i))=x;      % store the solution
        end

    otherwise
        fprintf('\n The function value specified in opts.fName is not supported!');
end

author: thomaskuestner


text:
function centroid = polygonCentroid3d(varargin)
%POLYGONCENTROID3D Centroid (or center of mass) of a polygon
%
%   PTC = polygonCentroid3d(POLY)
%   Computes center of mass of a polygon defined by POLY. POLY is a N-by-3
%   array of double containing coordinates of polygon vertices.
%
%   PTC = polygonCentroid3d(VX, VY, VZ)
%   Specifies vertex coordinates as three separate arrays.
%
%   Example
%     % compute centroid of a basic polygon
%     poly = [0 0 0; 10 0 10;10 10 20;0 10 10];
%     centro = polygonCentroid3d(poly)
%     centro =
%         5.0000    5.0000    10.0000
%
%   See also
%   polygons3d, polygonArea3d, polygonCentroid
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-09-18
% Copyright 2007 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

if nargin == 1
    % polygon is given as a single argument
    pts = varargin{1};

elseif nargin == 2
    % polygon is given as 3 corodinate arrays
    px = varargin{1};
    py = varargin{2};
    pz = varargin{3};
    pts = [px py pz];
end

% create supporting plane (assuming first 3 points are not colinear...)
plane = createPlane(pts(1:3, :));

% project points onto the plane
pts = planePosition(pts, plane);

% compute centroid in 2D
centro2d = polygonCentroid(pts);

% project back in 3D
centroid = planePoint(plane, centro2d);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [pk, final_manip] = SCO_null_space(robot, pathT)
%global robot
global parameters
global hfigures
hfigures.hpaths = figure;
hfigures.hcosts = figure;
hfigures.hee = figure;
hfigures.htheta = figure;
hfigures.hbest_costs = figure;
hfigures.hdtheta = figure;

%N waypoints at each trajectory
%N = parameters.N;
%K particles
K = parameters.K;

%generate K different paths
%starting from arbitrary positions.
for k=1:K
  q0 = uniform(-pi, pi, 1, robot.DOF)';
  [pathq, pathT] = plan_path_moore_penrose(robot, q0);
  %generate particle
  G{k}.pathq = pathq;
  G{k}.pathT = pathT;
  animate_local(robot, pathq)
  figure, plot(pathq'), legend('q1', 'q2', 'q3', 'q4')
end
%pick arbitrarily one...
initial_manip = compute_manip(robot, pathq);
G_inicial = G;

%animate one of the particles (random)
%this may not be the best particle, of course
animate_local(robot, pathq)

i=0;
best_probs = [];
%Main loop
%till convergence of trajectory cost function
while i < 20
    i
    i=i+1;
    G = add_noise_to_particle_set(G, pathT);

    %best particle
    [pk, best_prob] = best_particle(G);
    best_probs = [best_probs, best_prob];
    plot_info(G, pk, best_probs)
end

%plot particle info
best_particle_info(pk)
figure,
plot(best_probs)
title('best PROB at each iteration (sum of probs-weights along trajectory)')
%animate best particle
animate_local(robot, pk.pathq)
final_manip = compute_manip(robot, pk.pathq);
figure, plot(initial_manip), hold
plot(final_manip)
legend('initial_manip', 'final_manip')
title('manip at each time step')

figure, hold
for k=1:K
    plot(compute_manip(robot, G{k}.pathq))
end

% function plot_particles_costs(G)
% K = size(G,2);
% figure, hold
% for k=1:K
%   pk = G{k};
%   plot(pk.P)
% end
% title('Probabilities')

function plot_info(G, best_pk, best_probs)
global hfigures parameters
K = size(G,2);
figure(hfigures.hpaths),
clf
hold on
title('robot paths')
pp=[];
for k=1:K
    pk = G{k};
    plot(pk.pathq')
    pp(k) = sum(pk.P);
end
figure(hfigures.hcosts)
plot(sort(pp))
title('Trajectory sum of probabilities for each particle')

figure(hfigures.hbest_costs), clf, hold on
plot(best_probs)
title('Best trajectory probability at each time step (sum of prob-weights)')

function best_particle_info(pk)
figure,
subplot(3,1,1)
plot(pk.pathq(1,:))
subplot(3,1,2)
plot(pk.pathq(2,:))
subplot(3,1,3)
plot(pk.pathq(3,:))
title('joint values')
figure,
plot(pk.P)
title('Prob at each time step of particle SHOULD BE >0')
figure,
plot(pk.pathq')

title('BEST PARTICLE PATH')

function animate_robot(pk)
global robot
N = size(pk.path,2);
for i=1:N
   qi = pk.path(:,i);
   drawrobot3d(robot, qi)
end

%generate a prior path for the particles
function path = generate_path(robot, q0, qN, N)
deltaq=(qN-q0)/(N-1);
path = [];
for i=0:N-1
    qi = q0 + i*deltaq;
    path = [path qi];
    %drawrobot3d(robot, qi)
end

function mcost = mean_cost(G)
pk = G{1};
K = size(G,2);
N = size(pk.path,2);
cc = [];
for k=1:K
    pk = G{k};
    cc(k) = sum(pk.costs);
end
mcost = mean(cc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OBTAIN the best particle according to the defined cost function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pk, best_prob] = best_particle(G)
%global robot
pk = G{1};
K = size(G,2);
N = size(pk.pathq,2);
probs = [];

for k=1:K
    pk = G{k};
    probs(k) = sum(pk.P);
end
%we must minimize costs
%[val, index]=min(costs);
%alternatively maximize our defined weigths
[val, index]=max(probs);
fprintf('Best particle is index: k=%d\n', index)
pk = G{index};
best_prob = val;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute cost function to encode manipulability.
% Return:
% P: likelihood function
% q: the cost function it self
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [P, cost] = cost_function_manipulability(q)
global robot parameters
%compute nil-space!!
% J = manipulator_jacobian(robot, q);
manip = compute_manip(robot, q);
%define cost: cost is lower as manipulability is higher
cost = 1/(manip+0.01);
lambda = parameters.lambda_manip;
%return likelihood
P = exp(-(1/lambda)*cost);

function V = compute_link_velocity_1(theta, thetad)
V = 1*thetad(1)^2;

function V = compute_link_velocity_2(theta, thetad)
th1 = theta(1);
th2 = theta(2);

J = [-sin(th1)  -sin(th1+th2);
    cos(th1)    cos(th1+th2)];
V = J*[thetad(1) thetad(1)+thetad(2)]';
V = norm(V);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate a noisy particle set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G=generate_initial_particle_set(K, path)
for k=1:K
    G{k} = generate_particle(path);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate a path for each particle without noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pk = generate_particle(path)
global robot
%each path is in each row for qi
DOF = robot.DOF;
for i = 1:DOF
    pk.path(i,:) = path(i,:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add noise to the
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G=add_noise_to_particle_set(G, pathT)
K = size(G,2);
%for every particle in the set
for k=1:K
    pk = G{k};
    pathT = pk.pathT;
    pk = add_noise_to_particle_null_space(pk, pathT);
    G{k} = pk;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% !! Add noise along the null space!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pk = add_noise_to_particle_null_space(pk, pathT)
global parameters robot
N = size(pk.pathq, 2);
%each path is in each row for qi
DOF = robot.DOF;

% at each particle add some noise!
%move along the path of the particle
for i = 1:N
    %current joint position
    q1 = pk.pathq(:,i);
    %generate some samples along the nulls space
    qq = sample_from_null_space(q1);
    %populate with the current pose!!!
    qq = [qq q1];
    p = [];
    costs = [];
    %for all the new samples!
    for j=1:size(qq,2)
        %weight the samples,
        %compute cost of obstacles and manipulability
        [Pm, costm] = cost_function_manipulability(qq(:,j));
        if robot.DOF == 4
            [Po] = cost_function_distance_4dof(qq(:,j));
        else
            [Po] = cost_function_distance_7dof(qq(:,j));
        end
        p(j)=Pm*Po;
        %p(j)=Pm;
    end
    %select the best one
    %from the sample performed based on probability
    %max probatility!
    [prob, index] = max(p);
    q = qq(:,index);
    pk.pathq(:,i) = q;
    %and current probability
    pk.P(i) = prob;
end

%obtain samples from the null space
%near to q
function qq=sample_from_null_space(q1)
global robot
global parameters
%sigma_noise = parameters.noise_sigma_null_space;
alpha = parameters.alpha;
time_step = parameters.time_step;

q = q1;
the = 0;
qq=[];
%sigma_noise = 0.3;
%alpha = mvnrnd(0, sigma_noise);
%alpha = 0.3;
ss = +1;
%iterate at each time step: forward movement
while the < abs(alpha)
    %compute null_space
    if robot.DOF==4
        qdnull = null_space_4dof(robot, q);
    else
        qdnull = null_space_7dof(robot, q);
    end
    if norm(qdnull) < 0.01
        break
    end
    qdnull = qdnull/norm(qdnull);
    dq = ss*time_step*qdnull;
    the = the + norm(dq);
    q = q  + dq;
    qq = [qq q];
end

ss = -1;
q = q1;
the = 0;
%backward movement
while the < abs(alpha)
    %compute null_space
        %compute null_space
    if robot.DOF==4
        qdnull = null_space_4dof(robot, q);
    else
        qdnull = null_space_7dof(robot, q);
    end
    if norm(qdnull) < 0.01
        break
    end
    qdnull = qdnull/norm(qdnull);
    dq = ss*time_step*qdnull;
    the = the + norm(dq);
    q = q  + dq;
    qq = [qq q];
end

%
% Uniform value betweenn min_val and max_val
%
function delta = uniform(min_val, max_val, n, m)
delta = (max_val-min_val)*rand(n, m) + min_val;

repo: ARTE


text:
% TSERRORPLOT(Y, E)
%
% TSERRORPLOT(X, Y, E)
%
% TSERRORPLOT(X, Y, L, U)
%
% Plots D signals with errors.
%
% X : N x 1 matrix of inputs
% Y : N x D matrix of function values
% E : N x D or N x 1 or 1 x D or 1 x 1 matrix of errors
% L : N x D or N x 1 or 1 x D or 1 x 1 matrix of lower errors
% U : N x D or N x 1 or 1 x D or 1 x 1 matrix of upper errors
%
% See ERRORPLOT for details on the optional arguments.

% Copyright (c) 2010 Jaakko Luttinen

function hax = tserrorplot(Y, E, varargin)

% Check whether the input was (Y,E) or (X,Y,E) or (X,Y,L,U)
if nargin >= 3 && isnumeric(varargin{1})
  x = Y;
  Y = E;
  if nargin >= 4 && isnumeric(varargin{2})
    % (X,Y,L,U)
    L = varargin{1};
    U = varargin{2};
    params = varargin(3:end);
  else
    % (X,Y,E)
    L = varargin{1};
    U = varargin{1};
    params = varargin(2:end);
  end
else
  % (Y,E)
  x = (1:(size(Y,1)))';
  L = E;
  U = E;
  params = varargin;
end

% Parse optional arguments
options = struct();
[options, errmsg, remopts] = argparse(options, params{:});
error(errmsg);

% Scale L and U to proper sizes
L = bsxfun(@times, ones(size(Y)), L);
U = bsxfun(@times, ones(size(Y)), U);

figure

D = size(Y,2);
hax = zeros(D,1);

for d=1:D
  % Plot with errors
  hax(d) = subplot(D,1,d);
  errorplot(x(:),Y(:,d),L(:,d),U(:,d), remopts{:});
end

if nargout < 1
  clear hax;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
classdef sincos < dml.method
% SINCOS circular regression by decomposing into sine and cosine.
%
%   DESCRIPTION
%   Angle is represented as sine and cosine on which regularized linear
%   regression is performed
%
%   EXAMPLE
%
%   DEVELOPER
%   Marcel van Gerven (m.vangerven@donders.ru.nl)

    properties

      regressor = dml.enet('family','gaussian','alpha',1);

      sinreg; % regressor applied to sine
      cosreg; % regressor applied to cosine

    end

    methods

        function obj = sincos(varargin)

            obj = obj@dml.method(varargin{:});
        end

        function obj = train(obj,X,Y)

          obj.sinreg = obj.regressor.train(X,sin(Y));
          obj.cosreg = obj.regressor.train(X,cos(Y));

        end

        function Y = test(obj,X)

          Y = atan2(obj.sinreg.test(X),obj.cosreg.test(X));

        end

    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [l] = qt2l(qt)
% Convert volume from US liquid quarts to liters.
% Chad Greene 2012
l = qt*0.946352946;


text:
classdef nme_gen_acc < mp.nme_gen_ac & mp.form_acc

%   MATPOWER
%   Copyright (c) 2019-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% This is small script that demonstrate the computation of extrinsic
% parameters using 3D structures.
% This test was build from data provided by Daniel Small (thank you Daniel!)

%-- Image points (in pixels):

x = [479.5200  236.0800
  608.4100  415.3700
  461.0000   40.0000
  451.4800  308.7000
  373.9900  314.8900
  299.3200  319.1300
  231.5500  321.3700
  443.7300  282.9200
  378.3600  288.3000
  314.6900  292.7400
  255.4700  296.2300]';

% 3D world coordinates:

X = [ 0         0         0
   54.0000         0         0
         0         0   40.5000
   27.0000   -8.4685   -2.3750
   27.0000  -18.4685   -2.3750
   27.0000  -28.4685   -2.3750
   27.0000  -38.4685   -2.3750
   17.0000   -8.4685   -2.3750
   17.0000  -18.4685   -2.3750
   17.0000  -28.4685   -2.3750
   17.0000  -38.4685   -2.3750]';

%------------ Intrinsic parameters:
%--- focal:
fc = [ 395.0669  357.1178 ]';
%--- Principal point:
cc = [ 380.5387  230.5278 ]';
%--- Distortion coefficients:
kc = [-0.2601    0.0702   -0.0019   -0.0003         0]';
%--- Skew coefficient:
alpha_c = 0;

%----- Computation of the pose of the object in space
%----- (or the rigid motion between world reference frame and camera ref. frame)
[om,T,R] = compute_extrinsic(x,X,fc,cc,kc,alpha_c);

%--- Try to reproject the structure to see if the computed pose makes sense:
x2 = project_points2(X_1,omckk,Tckk,fc,cc,kc,alpha_c);

% Graphical output:
figure(2);
plot(x(1,:),x(2,:),'r+');
hold on;
plot(x2(1,:),x2(2,:),'go');
hold off;
axis('equal');
axis('image');
title('red crosses: data, green circles: reprojected structure -- IT WORKS!!!');
xlabel('x');
ylabel('y');


text:
classdef mme_line3p < mp.mm_element

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function name = name(obj)
            name = 'line3p';
        end

        function obj = data_model_update(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            nme = obj.network_model_element(nm);

            pp = nm.get_idx('port');
            nn = nme.np/2;

            %% branch active power flow
            for p = 1:nn
                s_fr = nm.soln.gs_(pp.i1.line3p(p):pp.iN.line3p(p)) * dm.base_kva;
                s_to = nm.soln.gs_(pp.i1.line3p(p+nn):pp.iN.line3p(p+nn)) * dm.base_kva;

                %% update in the data model
                dme.tab.(sprintf('pl%d_fr', p))(dme.on) = real(s_fr);
                dme.tab.(sprintf('ql%d_fr', p))(dme.on) = imag(s_fr);
                dme.tab.(sprintf('pl%d_to', p))(dme.on) = real(s_to);
                dme.tab.(sprintf('ql%d_to', p))(dme.on) = imag(s_to);
            end
        end
    end     %% methods
end         %% classdef


text:
function err = test_factorize (A)
%TEST_FACTORIZE test the accuracy of the factorization object
%
% Example
%   test_factorize (A) ;    % where A is a square matrix (sparse or dense)
%
% See also test_all, factorize1, factorize, inverse, mldivide

% Copyright 2009, Timothy A. Davis, University of Florida

if (nargin < 1)
    A = rand (100) ;
end

[m n] = size (A) ;
err = 0 ;
if (min (m,n) > 0)
    anorm = norm (A,1) ;
else
    anorm = 1 ;
end

for nrhs = 1:3

    for bsparse = 0:1

        % fprintf ('n %4d nrhs %d bsparse %d : ', n, nrhs, bsparse) ;

        b = rand (m,nrhs) ;
        if (bsparse)
            b = sparse (b) ;
        end

        %-------------------------------------------------------------------
        % test backslash and related methods
        %-------------------------------------------------------------------

        % method 0:
        x = A\b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 1:
        S = inverse (A) ;
        x = S*b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 2:
        if (m == n)
            F = factorize1 (A) ;
            x = F\b ;       err = check_resid (err, anorm, A, x, b) ;
        end

        % method 3:
        if (m == n)
            S = inverse (F) ;
            x = S*b ;       err = check_resid (err, anorm, A, x, b) ;
        end

        % method 4:
        F = factorize (A) ;
        x = F\b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 5:
        S = inverse (F) ;
        x = S*b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 6:
        if (m == n)
            [L,U,p] = lu (A, 'vector') ;
            x = U \ (L \ (b (p,:))) ;
            err = check_resid (err, anorm, A, x, b) ;
        end

        % method 7 (ack!)
        if (m == n)
            S = inv (A) ;
        else
            S = pinv (full (A)) ;
            if (issparse (A))
                S = sparse (S) ;
            end
        end
        x = S*b ;           err = check_resid (err, anorm, A, x, b) ;

        %-------------------------------------------------------------------
        % test mtimes
        %-------------------------------------------------------------------

        S = inverse (F) ;
        d = rand (n,1) ;
        x = F*d ;
        y = A*d ;
        z = S\d ;
        e = max (norm (x-y,1), norm (x-z,1)) ;
        if (e > 0)
            error ('mtimes error') ;
        end

        if (m == n)
            F = factorize1 (A) ;
            S = inverse (F) ;
            d = rand (n,1) ;
            x = F*d ;
            y = A*d ;
            z = S\d ;
            e = max (norm (x-y,1), norm (x-z,1)) ;
            if (e > 0)
                error ('mtimes error') ;
            end
        end

        %-------------------------------------------------------------------
        % test slash and related methods
        %-------------------------------------------------------------------

        b = rand (nrhs,n) ;
        if (bsparse)
            b = sparse (b) ;
        end

        % method 0:
        x = b/A ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 1:
        S = inverse (A) ;
        x = b*S ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 2:
        if (m == n)
            F = factorize1 (A) ;
            x = b/F ;       err = check_resid (err, anorm, A, x, b, 1) ;
        end

        % method 3:
        if (m == n)
            S = inverse (F) ;
            x = b*S ;       err = check_resid (err, anorm, A, x, b, 1) ;
        end

        % method 4:
        F = factorize (A) ;
        x = b/F ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 5:
        S = inverse (F) ;
        x = b*S ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 6:
        if (m == n)
            [L,U,p] = lu (A, 'vector') ;
            x = (b / U) / L ; x (:,p) = x ;
            err = check_resid (err, anorm, A, x, b, 1) ;
        end

        % method 7 (ack!)
        if (m == n)
            S = inv (A) ;
        else
            S = pinv (full (A)) ;
            if (issparse (A))
                S = sparse (S) ;
            end
        end
        x = b*S ;           err = check_resid (err, anorm, A, x, b, 1) ;

        %-------------------------------------------------------------------
        % test double
        %-------------------------------------------------------------------

        Y = double (inverse (A)) ;
        if (m == n)
            Z = inv (A) ;
        else
            Z = pinv (full (A)) ;
        end
        e = norm (Y-Z,1) ;
        if (n > 0)
            e = e / norm (Z,1) ;
        end
        err = max (e, err) ;

        %-------------------------------------------------------------------
        % test subsref
        %-------------------------------------------------------------------

        F = factorize (A) ;
        Y = inverse (A) ;
        if (numel (A) > 1)
            if (F (end) ~= A (end))
                error ('factorization subsref error') ;
            end
            if (F.A (end) ~= A (end))
                error ('factorization subsref error') ;
            end
        end
        if (n > 0)
            if (F (1,1) ~= A (1,1))
                error ('factorization subsref error') ;
            end
            if (F.A (1,1) ~= A (1,1))
                error ('factorization subsref error') ;
            end
            e = abs (Y (1,1) - Z (1,1)) ;
            err = max (e,err) ;
            if (m > 1 && n > 1)
                e = norm (Y (1:2,1:2) - Z (1:2,1:2), 1) ;
                err = max (e,err) ;
            end
        end
        if (m > 3 && n > 1)
            if (any (F (2:end, 1:2) - A (2:end, 1:2)))
                error ('factorization subsref error') ;
            end
            if (any (F (2:4, :) - A (2:4, :)))
                error ('factorization subsref error') ;
            end
            if (any (F (:, 1:2) - A (:, 1:2)))
                error ('factorization subsref error') ;
            end
        end

        %-------------------------------------------------------------------
        % test update/downdate
        %-------------------------------------------------------------------

        if (F.kind == 6)
            w = rand (n,1) ;
            b = rand (n,1) ;
            % update
            G = F + w ;
            x = G\b ;       err = check_resid (err, anorm, A+w*w', x, b) ;
            % downdate
            G = G - w ;
            x = G\b ;       err = check_resid (err, anorm, A, x, b) ;
            clear G
        end

        %-------------------------------------------------------------------
        % test size
        %-------------------------------------------------------------------

        [m1 n1] = size (F) ;
        [m n] = size (A) ;
        if (m1 ~= m || n1 ~= n)
            error ('size error') ;
        end
        [m1 n1] = size (Y) ;
        if (m1 ~= n || n1 ~= m)
            error ('pinv size error') ;
        end
        if (size (Y,1) ~= n || size (Y,2) ~= m)
            error ('pinv size error') ;
        end
        if (size (F,1) ~= m || size (F,2) ~= n)
            error ('size error') ;
        end

        %-------------------------------------------------------------------
        % test mtimes
        %-------------------------------------------------------------------

        d = rand (1,m) ;
        x = d*F ;
        y = d*A ;
        z = d/Y ;
        e = max (norm (x-y,1), norm (x-z,1)) ;
        if (e > 0)
            error ('mtimes error') ;
        end

        if (m == n)
            F = factorize1 (A) ;
            Y = inverse (F) ;
            d = rand (1,m) ;
            x = d*F ;
            y = d*A ;
            z = d/Y ;
            e = max (norm (x-y,1), norm (x-z,1)) ;
            if (e > 0)
                error ('mtimes error') ;
            end
        end

        %-------------------------------------------------------------------
        % test inverse
        %-------------------------------------------------------------------

        Y = double (inverse (inverse (A))) ;
        e = norm (A-Y,1) ;
        if (e > 0)
            error ('inverse error') ;
        end

    end
end

fprintf ('.') ;

if (err > 1e-8)
    fprintf ('error: %8.3e\n', err) ;
    error ('error is too high!') ;
end

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

function err = check_resid (err, anorm, A, x, b, transposed)
if (nargin < 6)
    transposed = 0 ;
end
[m n] = size (A) ;

if (transposed)
    if (m >= n)
        e = norm (A'*x'-b',1) / (anorm + norm (x,1)) ;
    else
        e = norm (A*(A'*x')-A*b',1) / (anorm + norm (x,1)) ;
    end
else
    if (m <= n)
        e = norm (A*x-b,1) / (anorm + norm (x,1)) ;
    else
        e = norm (A'*(A*x)-A'*b,1) / (anorm + norm (x,1)) ;
    end
end

if (min (m,n) > 1)
    if (issparse (A) && issparse (b))
        if (~issparse (x))
            error ('x must be sparse') ;
        end
    else
        if (issparse (x))
            error ('x must be full') ;
        end
    end
end
err = max (err, e) ;


text:
function [Y,L,S,params] = loadSyntheticProblem(varargin)
% [Y,L,S,params] = loadSyntheticProblem(problemName, referenceSolutionDirectory)
%
% Returns parameters and solution for the following problems
%   (parameters are tuned so that each problem has the
%    the same true solution (L,S) )
% Or, call this with no input arguments to set the path...
%
% "Lag": min   lambdaL*||L||_* + lambdaS*||S||_1 + .5||L+S-Y||_F^2
%
% "Sum": min   ||L||_* + lambdaSum*||S||_1
%        subject to ||L+S-Y||_F <= epsilon
%
% "Max": min   max( ||L||_* , lambdaMax*||S||_1  )
%        subject to ||L+S-Y||_F <= epsilon
%
% "Dual-sum"  min  .5||L+S-Y||_F^2
%         subject to ||L||_* + lambdaSum*||S||_1  <= tauSum
%
% "Max-sum"  min  .5||L+S-Y||_F^2
%         subject to max( ||L||_*, lambdaMax*||S||_1)  <= tauMax
%
% so params has the following fields:
%   .lambdaL
%   .lambdaS
%   .lambdaSum
%   .lambdaMax
%   .tauSum
%   .tauMax
%   .epsilon
%
%   .objective  This is a function of (L,S) and is the objective
%               of the "Lagrangian" formulation
%
%   .L1L2   set to 0, 'rows', or 'sum' for using l1 or l1l2 block norms
%
% Stephen Becker, March 6 2014

% addpath for the proxes and solvers
Mfilename = mfilename('fullpath');
baseDirectory = fileparts(Mfilename);
addpath(genpath(baseDirectory));

if nargin==0
    % we just setup the path, nothing else...
    return;
end

if nargin >= 2
    baseDirectory = varargin{2};
else
    baseDirectory = fullfile(baseDirectory,'referenceSolutions');
end
norm_nuke = @(L) sum(svd(L,'econ'));
vec = @(X) X(:);
problemString = varargin{1};

%%
% -- Common problem strings --
% problemString = 'small_randn_v2';
% problemString = 'medium_randn_v1';
% problemString = 'medium_exponentialNoise_v1';
% problemString = 'large_exponentialNoise_square_v1';
% problemString = 'large_exponentialNoise_rectangular_v1';
% problemString = 'huge_exponentialNoise_rectangular_v1'; % not yet done!
% problemString = 'nsa500_SNR80_problem1_run1';

fprintf('%s\n** Problem %s **\n\n',char('*'*ones(60,1)), problemString );
fileName      = fullfile(baseDirectory,[problemString,'.mat']);
if exist(fileName,'file')
    fprintf('Loading reference solution from file %s\n', fileName);
    load(fileName); % should contain L, S and Y at least, and preferably lambdaL and lambdaS
else
    [L,S,Y] = deal( [] );
end
% Default:
L1L2  = 0;

% For problems from the NSA software package:
if length(problemString)>3 && strcmpi( problemString(1:3), 'nsa' )
    % This sscanf line is broken in versions of Matlab after R2010a
    tmp = sscanf(problemString,'nsa%d_SNR%d_problem%d_run%d');
    m       = tmp(1); % either 500, 1000 or 1500
    SNR     = tmp(2); % either 80 or 45
    problem_type = tmp(3);
    run_counter  = tmp(4);
    % e.g.,  'nsa500_SNR80_problem1_run1'
    n  = m;
%     % Variable parameters
%     m = 500; n = m;
%     problem_type = 1; % between 1 and 4. Controls rank/sparsity
%     SNR          = 80;
%     run_counter  = 1; % between 1 and 10. Only affects seed

%     % Possible sparsity amounts:
%     p_array = ceil( m*n*[ 0.05, 0.1, 0.05, 0.1 ] );
%     r_array = ceil(min(m,n)*[  0.05, 0.05, 0.1, 0.1 ] );

    % Changing the conventions, March 14. They had 4 cases, for p = {.05,.1}
    %   and r = {.05, .1 } (these are sparsity and rank percentages)
    % I don't want to do that many tests since we will display graphically
    %   and not in tables. I also want more values of p and r. So I will
    %   make p and r be linked.
    p_array = ceil( m*n*[ 0.05, 0.1, 0.2, 0.3 ] );
    r_array = ceil(min(m,n)*[  0.05, 0.1, 0.2, 0.3 ] );

    p = p_array(problem_type);
    r = r_array(problem_type);
    pwr = r+p/(m*n)*100^2/3;
    stdev = sqrt(pwr/10^(SNR/10));
    delta = sqrt(min(m,n)+sqrt(8*min(m,n)))*stdev;
    lambdaSum = 1/sqrt( max(m,n) ); % "xi" in NSA
else
    switch problemString
        case 'small_randn_v1'
            my_rng(101);
            m   = 40; n = 42; lambdaL = .5; lambdaS = 1e-1;
            if isempty(Y)
                Y   = randn(m,n);
                printEvery = 500; tol = 1e-12;  maxIts = 3e3;
            end

        case 'small_randn_v2'
            my_rng(102);
            m   = 40; n = 42; lambdaL = .5; lambdaS = 1e-1;
            if isempty(Y)
                Y   = randn(m,n);
                printEvery = 500; tol = 1e-12; maxIts = 3e3;
            end
        case 'medium_randn_v1'
            my_rng(101);
            m   = 400; n = 500; lambdaL = .25; lambdaS = 1e-2;
            if isempty(Y)
                Y   = randn(m,n);
                printEvery = 1; %tol = 1e-4; maxIts = 100;
                tol = 1e-10; maxIts = 1e3;
            end
        case 'medium_exponentialNoise_v1'
            % Converges much faster than 'medium_randn_v1'
            my_rng(101);
            m   = 400; n = 500; lambdaL = .25; lambdaS = 1e-2;
            rk  = round(.5*min(m,n) );
            Q1  = haar_rankR(m,rk,false);
            Q2  = haar_rankR(n,rk,false);
            Y   = Q1*diag(.1+rand(rk,1))*Q2';
            mdn = median(abs(Y(:)));
            Y   = Y + exprnd( .1*mdn, m, n );

            printEvery = 5; tol = 1e-10; maxIts = 1000;

        case 'medium_exponentialNoise_v2'
            % Converges much faster than 'medium_randn_v1'
            % our quasi-Newton scheme does well here
            my_rng(101);
            m   = 400; n = 500; lambdaL = .25; lambdaS = 1e-2;
            rk  = round(.05*min(m,n) );
            Q1  = haar_rankR(m,rk,false);
            Q2  = haar_rankR(n,rk,false);
            Y   = Q1*diag(.1+rand(rk,1))*Q2';
            mdn = median(abs(Y(:)));
            Y   = Y + exprnd( .1*mdn, m, n );

            printEvery = 5; tol = 1e-10; maxIts = 1000;

        case 'medium_exponentialNoise_v3'
            % Converges much faster than 'medium_randn_v1'
            my_rng(101);
            m   = 400; n = 500; lambdaL = .25; lambdaS = 1e-2;
            rk  = round(.1*min(m,n) );
            Q1  = haar_rankR(m,rk,false);
            Q2  = haar_rankR(n,rk,false);
            Y   = Q1*diag(.1+rand(rk,1))*Q2';
            mdn = median(abs(Y(:)));
            Y   = Y + exprnd( .1*mdn, m, n );

            printEvery = 5; tol = 1e-10; maxIts = 1000;

        case 'large_exponentialNoise_square_v1'
            my_rng(201);
            m   = 1e3; n = m; lambdaL = .25; lambdaS = 1e-2;
            rk  = round(.5*min(m,n) );
            Q1  = haar_rankR(m,rk,false);
            Q2  = haar_rankR(n,rk,false);
            Y   = Q1*diag(.1+rand(rk,1))*Q2';
            mdn = median(abs(Y(:)));
            Y   = Y + exprnd( .1*mdn, m, n );

            printEvery = 5; tol = 1e-8; maxIts = 250;

        case 'large_exponentialNoise_rectangular_v1'
            my_rng(201);
            m   = 5e3; n = 200; lambdaL = .25; lambdaS = 5e-3; % OK
            rk  = round(.5*min(m,n) );
            Q1  = haar_rankR(m,rk,false);
            Q2  = haar_rankR(n,rk,false);
            Y   = Q1*diag(.1+rand(rk,1))*Q2';
            mdn = median(abs(Y(:)));
            Y   = Y + exprnd( .1*mdn, m, n );

            printEvery = 5; tol = 1e-8; maxIts = 100;

        case 'huge_exponentialNoise_rectangular_v1';
            disp('Warning: lambda values not yet tuned');
            my_rng(301);
            m   = 5e4; n = 200; lambdaL = .25; lambdaS = 5e-2;
            rk  = round(.5*min(m,n) );
            Q1  = haar_rankR(m,rk,false);
            Q2  = haar_rankR(n,rk,false);
            Y   = Q1*diag(.1+rand(rk,1))*Q2';
            mdn = median(abs(Y(:)));
            Y   = Y + exprnd( .1*mdn, m, n );

            printEvery = 5; tol = 1e-10; maxIts = 100;

        case 'escalator_small';
            % added March 29 2015
            load escalator_data_2015_small % has Y, mat, m, n, nFrames
            m   = m*n; % 64 x 79, 200 frames
            n   = nFrames;
            lambdaL     = 3;
            lambdaS     = 0.06;
            printEvery  = 1;
            tol         = 1e-9;
            maxIts      = 100;
            L0          = Y;
            S0          = 0*Y;

        case 'verySmall_l1'
            load verySmall_l1
            % Used CVX to generate solution

        case 'verySmall_l1l2'
            load verySmall_l1l2
            % Used CVX to generate solution
            L1L2 = 'rows';

        otherwise
            error('bad problem name');
    end
end

% Generate new reference soln if not precomputed
if isempty(L)
    fprintf('Generating reference solution...\n');
    if length(problemString)>3 && strcmpi( problemString(1:3), 'nsa' )
        % We solve Sum problem, which gives us a lot of parameters but not
        % the parameters for the Lagrangian problem unfortunately.

        % I am changing the seeds format
        seed = 601;
        my_rng( seed + run_counter );
        [Y, optimal_L, optimal_S, deltabar]=create_data_noisy_L2(m,n,p,r,stdev);

        tol = 0.0001;
        denoise_flag=0;

        %         [L,S,out] = nsa(Y,stdev,tol,denoise_flag,optimal_L,optimal_S);
        %         [L,S,out] = nsa(Y,-delta,tol,denoise_flag,optimal_L,optimal_S);

        [L,S,out] = nsa(Y,-delta,tol,denoise_flag,[],[],lambdaSum);
        fprintf('Discrepancy with intitial signal L is %.2e\n', ...
            norm(L-optimal_L,'fro')/norm(L,'fro') );
        fprintf('Discrepancy with intitial signal S is %.2e\n', ...
            norm(S-optimal_S,'fro')/norm(S,'fro') );

        lambdaL     = [];
        lambdaS     = [];
        maxIts      = [];
    else
        tic
        opts    = struct('printEvery',printEvery,'FISTA',true,'tol',tol,'maxIts',maxIts);
        %[L,S,errHist] = solver_RPCA_Lagrangian_simple(Y,lambdaL,lambdaS,[],opts);
        opts.quasiNewton  = true; opts.FISTA=false; opts.BB = false;
        opts.SVDstyle   = 1; % for highest accuracy
        if exist('L0','var')
            opts.L0     = L0;
        end
        if exist('S0','var')
            opts.S0     = S0;
        end
        [L,S,errHist] = solver_RPCA_Lagrangian(Y,lambdaL,lambdaS,[],opts);
        toc
    end

    save(fileName,'L','S','Y','lambdaL','lambdaS','tol','maxIts');
end

sigma   = svd(L,'econ');
rk      = sum( sigma > 1e-8 );
normL   = sum(sigma);

params = [];
params.L1L2        = L1L2;
if L1L2
    if strcmpi(L1L2,'rows')
        params.objective    = @(L,S) lambdaL*norm_nuke(L)+lambdaS*sum( sqrt( sum(S.^2,2) ) ) + ...
            .5*norm(vec(L+S-Y))^2;
        normS   = sum( sqrt( sum(S.^2,2) ) );
    elseif strcmpi(L1L2,'cols')
        params.objective    = @(L,S) lambdaL*norm_nuke(L)+lambdaS*sum( sqrt( sum(S.^2,1) ) ) + ...
            .5*norm(vec(L+S-Y))^2;
        normS   = sum( sqrt( sum(S.^2,1) ) );
    else error('bad value for L1L2');
    end
else
    params.objective    = @(L,S) lambdaL*norm_nuke(L)+lambdaS*norm(S(:),1) + ...
        .5*norm(vec(L+S-Y))^2;
    normS   = norm(S(:),1);
end
params.epsilon     = norm( L + S - Y, 'fro' );
params.lambdaL     = lambdaL;
params.lambdaS     = lambdaS;
if exist('lambdaSum','var')
    params.lambdaSum   = lambdaSum;
else
    params.lambdaSum   = lambdaS/lambdaL;
end
params.lambdaMax   = normL/normS;
params.tauSum      = normL + params.lambdaSum*normS;
params.tauMax      = max(normL,params.lambdaMax*normS );

fprintf('\tSize %d x %d\n\tL has nuclear norm %.3f and rank %d of %d possible\n', ...
    m,n,normL, rk, length(sigma) );
fprintf('\tS has l1 (or l1l2) norm %.3f and %.2f%% of its elements are nonzero\n', ...
    normS, 100*sum( abs(S(:)) > 1e-8)/numel(S) );
fprintf('\t||L+S-Y||/||Y|| is %.2e\n', params.epsilon/norm(Y,'fro' ) );

end % end of main function

function Q = haar_rankR(n, r, cmplx)
% Q = HAAR_RANKR( N , R)
%   returns a N x N unitary matrix Q
%   from the Haar measure on the
%   Circular Unitary Ensemble (CUE)
%
% Q = HAAR_RANKR( N, R, COMPLEX )
%   if COMPLEX is false, then returns a real matrix from the
%   Circular Orthogonal Ensemble (COE). By default, COMPLEX = true
%
% For more info, see http://www.ams.org/notices/200705/fea-mezzadri-web.pdf
% "How to Generate Random Matrices fromthe Classical Compact Groups"
% by Francesco Mezzadri (also at http://arxiv.org/abs/math-ph/0609050 )
%
% Other references:
%   "How to generate a random unitary matrix" by Maris Ozols
%
% Stephen Becker, 2/24/11. updated 11/1/12 for rank r version.
%
% To test that the eigenvalues are evenly distributed on the unit
% circle, do this:   a = angle(eig(haar(1000))); hist(a,100);
%
% This calls "randn", so it will change the stream of the prng

if nargin < 3, cmplx = true; end
if nargin < 2 || isempty(r) , r=n; end

if cmplx
    z = (randn(n,r) + 1i*randn(n,r))/sqrt(2.0);
else
    z = randn(n,r);
end
[Q,R] = qr(z,0); % time consuming
d = diag(R);
ph = d./abs(d);
% Q = multiply(q,ph,q) % in python. this is q <-- multiply(q,ph)
%   where we multiply each column of q by an element in ph

%Q = Q.*repmat( ph', n, 1 );
%  use bsxfun for this to make it fast...
Q = bsxfun( @times, Q, ph' );

end % end of haar_rankR

function varargout = my_rng(varargin)
%MY_RNG Control the random number generator used by RAND, RANDI, and RANDN (SRB version)
%   MY_RNG(SD) seeds the random number generator using the non-negative
%   integer SD so that RAND, RANDI, and RANDN produce a predictable
%   sequence of numbers.
% MODIFIED BY STEPHEN BECKER
%   See also RAND, RANDI, RANDN, RandStream, NOW.

%   See Choosing a Random Number Generator for details on these generators.

if exist('rng','file')
    [varargout{1:nargout}] = rng(varargin{:});
    return;
end
% Otherwise, imitate the functionality...

% -- SRB adding this --
error(nargchk(1,1,nargin));
error(nargoutchk(0,0,nargout));
arg1 = varargin{1};
% For R2008a, this doesn't work... (not sure what earliest version is)
if verLessThan('matlab','7.7')
    randn('state',arg1);
    rand('state',arg1);
else
    if verLessThan('matlab','8.2')
        RandStream.setDefaultStream(RandStream('mt19937ar', 'seed', arg1) );
    else
        RandStream.setGlobalStream(RandStream('mt19937ar', 'seed', arg1) );
    end
end
end % end of my_rng

author: stephenbeckr

repo: fastRPCA