text:
function  test_linear_svm()

    clc;
    clear;
    close all;

    %% Set algorithms
    if 0
        algorithms = gd_solver_list('ALL');
    else
        %algorithms = gd_solver_list('LS');
        %algorithms = gd_solver_list('NCG');
        %algorithms = gd_solver_list('BFGS');
        algorithms = {'SD-STD', 'SD-BKT'};
    end

    % # of classes (must not change)
    l = 2;

    %% prepare dataset
    if 1     % generate synthetic data
        n = 100;    % # of samples per class
        d = 3;      % # of dimensions
        std = 0.15; % standard deviation

        data = multiclass_data_generator(n, d, l, std);
        d = d + 1; % adding '1' row for intersect

        % train data
        x_train = [data.x_train; ones(1,l*n)];
        % assign y (label) {1,-1}
        y_train(data.y_train<=1.5) = -1;
        y_train(data.y_train>1.5) = 1;

        % test data
        x_test = [data.x_test; ones(1,l*n)];
        % assign y (label) {1,-1}
        y_test(data.y_test<=1.5) = -1;
        y_test(data.y_test>1.5) = 1;

    else    % load real-world data
        data = importdata('../data/mushroom/mushroom.mat');
        n = size(data.X,1);
        d = size(data.X,2) + 1;
        x_in = [data.X ones(n,1)]';
        y_in = data.y';

        perm_idx = randperm(n);
        x = x_in(:,perm_idx);
        y = y_in(perm_idx);

        % split data into train and test data
        % train data
        n_train = floor(n/8);
        x_train = x(:,1:n_train);
        y_train = y(1:n_train);
        x_train_class1 = x_train(:,y_train>0);
        x_train_class2 = x_train(:,y_train<0);
        n_class1 = size(x_train_class1,2);
        n_class2 = size(x_train_class2,2);

        % test data
        x_test = x(:,n_train+1:end);
        y_test = y(n_train+1:end);
        x_test_class1 = x_test(:,y_test>0);
        x_test_class2 = x_test(:,y_test<0);
        n_test_class1 = size(x_test_class1,2);
        n_test_class2 = size(x_test_class2,2);
        n_test = n_test_class1 + n_test_class2;

    end
    lambda = 0.1;
    w_opt = zeros(d,1);

    % set plot_flag
    if d > 4
        plot_flag = false;  % too high dimension
    else
        plot_flag = true;
    end

    %% define problem definitions
    problem = linear_svm(x_train, y_train, x_test, y_test, lambda);

    %% initialize
    w_init = randn(d,1);

    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% calculate solution
    if norm(w_opt)
    else
        % calculate solution
        w_opt = problem.calc_solution(1000);
    end
    f_opt = problem.cost(w_opt);
    fprintf('f_opt: %.24e\n', f_opt);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_iter = 100;
        options.verbose = true;
        options.f_opt = f_opt;
        options.store_w = true;

        switch algorithms{alg_idx}
            case {'SD-STD'}

                options.step_alg = 'fix';
                options.step_init = 1;
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-SCALE-EXACT'}

                options.sub_mode = 'SCALING';
                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'Newton-STD'}

                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-DAMP'}

                options.sub_mode = 'DAMPED';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'CG-PRELIM'}

                options.sub_mode = 'PRELIM';
                options.step_alg = 'exact';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-BKT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-EXACT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'exact';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-PRECON-EXACT'}

                options.sub_mode = 'PRECON';
                % diagonal scaling
                options.M = diag(diag(A));
                options.step_alg = 'exact';
                options.beta_alg = 'PR';

                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'NCG-BTK'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'NCG-WOLFE'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'strong_wolfe';
                options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'BFGS-H-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-H-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-B-BKT'}

                options.step_alg = 'backtracking';
                options.update_mode = 'B';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-B-EXACT'}

                options.step_alg = 'exact';
                options.update_mode = 'B';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'DAMPED-BFGS-BKT'}

                options.step_alg = 'backtracking';
                options.update_mode = 'Damping';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'DAMPED-BFGS-EXACT'}

                options.step_alg = 'exact';
                options.update_mode = 'Damping';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'BB'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = bb(problem, options);

            case {'SGD'}

                options.batch_size = 1;
                options.step = 0.1 * options.batch_size;
                %options.step_alg = 'decay';
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = sgd(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','cost', algorithms, w_list, info_list);
    display_graph('iter','gnorm', algorithms, w_list, info_list);

    % draw convergence sequence
    w_history = cell(1);
    cost_history = cell(1);
    for alg_idx=1:length(algorithms)
        w_history{alg_idx} = info_list{alg_idx}.w;
        cost_history{alg_idx} = info_list{alg_idx}.cost;
    end
    draw_convergence_sequence(problem, w_opt, algorithms, w_history, cost_history);

    % display classification results
    y_pred_list = cell(length(algorithms),1);
    accuracy_list = cell(length(algorithms),1);
    for alg_idx=1:length(algorithms)
        p = problem.prediction(w_list{alg_idx});
        % calculate accuracy
        accuracy_list{alg_idx} = problem.accuracy(p);

        fprintf('Classificaiton accuracy: %s: %.4f\n', algorithms{alg_idx}, problem.accuracy(p));

        % convert from {1,-1} to {1,2}
        p(p==-1) = 2;
        p(p==1) = 1;
        % predict class
        y_pred_list{alg_idx} = p;
    end

    % convert from {1,-1} to {1,2}
    y_train(y_train==-1) = 2;
    y_train(y_train==1) = 1;
    y_test(y_test==-1) = 2;
    y_test(y_test==1) = 1;
    if plot_flag
        display_classification_result(problem, algorithms, w_list, y_pred_list, accuracy_list, x_train, y_train, x_test, y_test);
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function [x,y]=simulate(a,b,c,d,x0,T)
%Simulation of dynamical system
x(1)=x0;
for t=2:T
x(t)=random('normal',mx(a,x(t-1)),vx(b,x(t-1)));
y(t)=random('normal',my(c,x(t)),vy(d,x(t-1)));
end

author: Sable

repo: mcbench-benchmarks


text:
function geometry_test2072 ( )

%*****************************************************************************80
%
%% TEST2072 tests TRIANGLE_POINT_NEAR_2D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 7;

  ptest = [ ...
     0.25,   0.25; ...
     0.75,   0.25; ...
     1.00,   1.00; ...
    11.00,   0.50; ...
     0.00,   1.00; ...
     0.50, -10.00; ...
     0.60,   0.60 ]';
  t = [ ...
    0.0, 1.0; ...
    0.0, 0.0; ...
    1.0, 0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST2072\n' );
  fprintf ( 1, '  For a triangle in 2D,\n' );
  fprintf ( 1, '  TRIANGLE_POINT_NEAR_2D computes the nearest\n' );
  fprintf ( 1, '  point to a point.\n' );

  r8mat_transpose_print ( 2, 3, t, '  Triangle vertices:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '           P                PN\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    p(1:2,1) = ptest(1:2,i);

    [ pn, dist ] = triangle_point_near_2d ( t, p );

    fprintf ( 1, '  %10f  %10f    %10f  %10f\n', p(1:2,1), pn(1:2,1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_sum2( )
% Test diskfun sum2() command.

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

f = @(x,y) 1 + x + y;
g = diskfun(f);
exact_int = pi;
pass(1) = abs(sum2(g) - exact_int) < tol;

f = @(x,y) 10*exp(-x.*y) +12*x.^5.*cos(2*y)-10./(1+exp(x-.3));
g = diskfun(f);
exact_int = 14.16225689953454;
pass(2) = abs(sum2(g) - exact_int) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [vhdr] = read_brainvision_vhdr(filename)

% READ_BRAINVISION_VHDR reads the known items from the BrainVision EEG
% header file and returns them in a structure.
%
% Use as
%   vhdr = read_brainvision_vhdr(filename)
%
% See also READ_BRAINVISION_EEG, READ_BRAINVISION_VMRK

% Copyright (C) 2003-2019, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

vhdr.DataFile         = read_ini(filename, 'DataFile=', '%s');
vhdr.MarkerFile       = read_ini(filename, 'MarkerFile=', '%s');
vhdr.DataFormat       = read_ini(filename, 'DataFormat=', '%s');
vhdr.DataOrientation  = read_ini(filename, 'DataOrientation=', '%s');
vhdr.BinaryFormat     = read_ini(filename, 'BinaryFormat=', '%s');
vhdr.NumberOfChannels = read_ini(filename, 'NumberOfChannels=', '%d');
vhdr.SamplingInterval = read_ini(filename, 'SamplingInterval=', '%f'); % microseconds

if ~isempty(vhdr.NumberOfChannels)
  % only show the warning once
  ft_warning once FieldTrip:fileio:UnknownResolution

  for i=1:vhdr.NumberOfChannels
    chan_str  = sprintf('Ch%d=', i);
    chan_info = read_ini(filename, chan_str, '%s');
    t = tokenize(chan_info, ',');
    vhdr.label{i} = t{1};
    vhdr.reference{i} = t{2};
    resolution = str2num(t{3}); % in microvolt
    if ~isempty(resolution)
      vhdr.resolution(i) = resolution;
    else
      ft_warning('FieldTrip:fileio:UnknownResolution', 'unknown resolution (i.e. recording units) for channel %d in "%s"', i, filename);
      vhdr.resolution(i) = 1;
    end
  end

  % switch the warning state back to the default
  ft_warning on FieldTrip:fileio:UnknownResolution
end

% compute the sampling rate in Hz
vhdr.Fs = 1e6/(vhdr.SamplingInterval);

% the number of samples is unkown to start with
vhdr.nSamples = Inf;

% confirm the names of the .vmrk and .eeg files as specified in the .vhdr file
[p, f, x] = fileparts(filename);
dataFile = fullfile(p, vhdr.DataFile);
[dum, dum, dataExt] = fileparts(vhdr.DataFile);

if ~isempty(vhdr.MarkerFile)
  markerFile = fullfile(p, vhdr.MarkerFile);
  [dum, dum, markerExt] = fileparts(vhdr.MarkerFile);
else
  % it is allowed to have a BrainVision dataset with only header and data, but no marker file
  markerFile = '';
  markerExt = '';
end

% in case they cannot be found, we can try whether they exist with the same name as the vhdr file
% this can happen if the researcher renamed the three files but not the header information
dataFileSameName = fullfile(p,[f dataExt]);
markerFileSameName = fullfile(p,[f markerExt]);

info = dir(dataFile);
if isempty(info)
  info = dir(filename);
  if ~isempty(info)
    ft_notice('could not find "%s" as specified in the .vhdr file, using "%s" instead', dataFile, dataFileSameName);
    vhdr.DataFile = [f dataExt]; % without full path
    dataFile = dataFileSameName; % with full path
  else
    ft_error('cannot determine the location of the data file "%s"', dataFile);
  end
else
  if ~strcmp(dataFile, dataFileSameName)
    ft_notice('could not find "%s" as specified in the .vhdr file, using "%s" instead', dataFile, dataFileSameName);
  end
end

if ~isempty(markerFile)
  info = dir(markerFile);
  if isempty(info)
    info = dir(markerFileSameName);
    if ~isempty(info)
      ft_notice('could not find "%s" as specified in the .vhdr file, using "%s" instead', markerFile, markerFileSameName);
      vhdr.MarkerFile = [f markerExt]; % without full path
      markerFile = markerFileSameName; % with full path
    else
      ft_error('cannot determine the location of the marker file "%s"', markerFile);
    end
  else
    if ~strcmp(markerFile, markerFileSameName)
      ft_notice('could not find "%s" as specified in the .vhdr file, using "%s" instead', markerFile, markerFileSameName);
    end
  end
end

% determine the number of samples by looking at the binary file
if strcmpi(vhdr.DataFormat, 'binary')
  % the data file is supposed to be located in the same directory as the header file
  % but that might be on another location than the present working directory

  info = dir(dataFile);
  switch lower(vhdr.BinaryFormat)
    case 'int_16'
      vhdr.nSamples = info.bytes./(vhdr.NumberOfChannels*2);
    case 'int_32'
      vhdr.nSamples = info.bytes./(vhdr.NumberOfChannels*4);
    case 'ieee_float_32'
      vhdr.nSamples = info.bytes./(vhdr.NumberOfChannels*4);
  end

elseif strcmpi(vhdr.DataFormat, 'ascii')
  vhdr.skipLines = 0;
  vhdr.skipColumns = 0;

  % Read ascii info from header (if available).
  dataPoints    = read_ini(filename, 'DataPoints=', '%d');
  skipLines     = read_ini(filename, 'SkipLines=', '%d');
  skipColumns   = read_ini(filename, 'SkipColumns=', '%d');
  decimalSymbol = read_ini(filename, 'DecimalSymbol=', '%s'); % This is not used in reading dataset yet

  if ~isempty(dataPoints); vhdr.nSamples = dataPoints; end
  if ~isempty(skipLines); vhdr.skipLines = skipLines; end
  if ~isempty(skipColumns); vhdr.skipColumns = skipColumns; end
  if ~isempty(decimalSymbol); vhdr.decimalSymbol = decimalSymbol; end

  if isempty(dataPoints) && strcmpi(vhdr.DataOrientation, 'vectorized')
    % this is a very inefficient fileformat to read data from, it looks like this:
    % Fp1   -2.129 -2.404 -18.646 -15.319 -4.081 -14.702 -23.590 -8.650 -3.957
    % AF3   -24.023 -23.265 -30.677 -17.053 -24.889 -35.008 -21.444 -15.896 -12.050
    % F7    -10.553 -10.288 -19.467 -15.278 -21.123 -25.066 -14.363 -10.774 -15.396
    % F3    -28.696 -26.314 -35.005 -27.244 -31.401 -39.445 -30.411 -20.194 -16.488
    % FC1   -35.627 -29.906 -38.013 -33.426 -40.532 -49.079 -38.047 -26.693 -22.852
    % ...
    fid = fopen_or_error(dataFile, 'rt');
    tline = fgetl(fid); % read the complete first line
    fclose(fid);
    t = tokenize(tline, ' ', true); % cut the line into pieces
    vhdr.nSamples = length(t) - 1;  % the first element is the channel label
  end
end

if isinf(vhdr.nSamples)
  ft_warning('cannot determine number of samples for this sub-fileformat');
end

% the number of trials is unkown, assume continuous data
vhdr.nTrials = 1;
vhdr.nSamplesPre = 0;

% ensure that the labels are in a column
vhdr.label = vhdr.label(:);
vhdr.reference = vhdr.reference(:);
vhdr.resolution = vhdr.resolution(:);

% read in impedance values
vhdr.impedances.channels = [];
vhdr.impedances.reference = [];
vhdr.impedances.ground = NaN;
vhdr.impedances.refChan = [];

fid = fopen_or_error(filename, 'rt');
while ~feof(fid)
  tline = fgetl(fid);
  if tline == -1
    break; % the end-of-file was reached with an empty line
  end
  if startsWith(tline, 'Impedance [')
    chanCounter = 0;
    refCounter = 0;
    impCounter = 0;
    while chanCounter<vhdr.NumberOfChannels && ~feof(fid)
      chan_info = fgetl(fid);
      if ~isempty(chan_info)
        impCounter = impCounter+1;
        [chanName,impedances] = strtok(chan_info,':');
        spaceList = strfind(chanName,' ');
        if ~isempty(spaceList)
          chanName = chanName(spaceList(end)+1:end);
        end
        if strfind(chanName,'REF_') == 1 % for situation where there is more than one reference
          refCounter = refCounter+1;
          vhdr.impedances.refChan(refCounter) = impCounter;
          if ~isempty(impedances)
            vhdr.impedances.reference(refCounter) = str2double(impedances(2:end));
          else
            vhdr.impedances.reference(refCounter) = NaN;
          end
        elseif strcmpi(chanName,'ref') % single reference
          refCounter = refCounter+1;
          vhdr.impedances.refChan(refCounter) = impCounter;
          if ~isempty(impedances)
            vhdr.impedances.reference(refCounter) = str2double(impedances(2:end));
          else
            vhdr.impedances.reference(refCounter) = NaN;
          end
        else
          chanCounter = chanCounter+1;
          if ~isempty(impedances)
            vhdr.impedances.channels(chanCounter,1) = str2double(impedances(2:end));
          else
            vhdr.impedances.channels(chanCounter,1) = NaN;
          end
        end
      end
    end
    if ~feof(fid)
      tline='';
      while ~feof(fid) && isempty(tline)
        tline = fgetl(fid);
      end
      if ~isempty(tline)
        if strcmp(tline(1:4),'Ref:')
          refCounter = refCounter+1;
          [chanName,impedances] = strtok(tline,':');
          if ~isempty(impedances)
            vhdr.impedances.reference(refCounter) = str2double(impedances(2:end));
          else
            vhdr.impedances.reference(refCounter) = NaN;
          end
        end
        if strcmpi(tline(1:4),'gnd:')
          [chanName,impedances] = strtok(tline,':');
          vhdr.impedances.ground = str2double(impedances(2:end));
        end
      end
    end
    if ~feof(fid)
      tline='';
      while ~feof(fid) && isempty(tline)
        tline = fgetl(fid);
      end
      if ~isempty(tline)
        if strcmpi(tline(1:4),'gnd:')
          [chanName,impedances] = strtok(tline,':');
          vhdr.impedances.ground = str2double(impedances(2:end));
        end
      end
    end
  end
end
fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ImageTakerFor3DShape < handle

    properties (Access = private)
        uMesh
        outDir
    end

    methods (Access = public)

        function obj = ImageTakerFor3DShape(cParams)
            obj.init(cParams)
        end

        function takeImages(obj)
            obj.plotUnfittedMesh();
            obj.takeIsoFrontImage();
            obj.takeBackImage();
            obj.takeFrontImage();
            obj.takeIsoBackImage();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.uMesh  = cParams.unfittedMesh;
            obj.outDir = cParams.outDir;
        end

        function takeIsoFrontImage(obj)
            fName = 'IsoFrontImage';
            az = 120;
            el = 30;
            fName = [obj.outDir,fName];
            view(az,el)
            obj.takeImage(fName)
        end

        function takeIsoBackImage(obj)
            fName = 'IsoBackImage';
            az = 240;
            el = 30;
            fName = [obj.outDir,fName];
            view(az,el)
            obj.takeImage(fName)
        end

        function takeFrontImage(obj)
            fName = 'YZfrontImage';
            az = 90;
            el = 0;
            fName = [obj.outDir,fName];
            view(az,el)
            obj.takeImage(fName)
        end

        function takeBackImage(obj)
            fName = 'YZbackImage';
            az = 90;
            el = 180;
            fName = [obj.outDir,fName];
            view(az,el)
            obj.takeImage(fName)
        end

        function plotUnfittedMesh(obj)
            f = figure();
            f.WindowState = 'maximized';
            obj.uMesh.plot();
        end

    end

    methods (Access = private, Static)

        function takeImage(name)
            print(name,'-dpng')
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function prob_test1341 ( )

%*****************************************************************************80
%
%% TEST01341 checks RIBESL against BESSEL_IX_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1341:\n' );
  fprintf ( 1, '  RIBESL computes values of Bessel functions\n' );
  fprintf ( 1, '  of NONINTEGER order.\n' );
  fprintf ( 1, '  BESSEL_IX_VALUES returns selected values of the\n' );
  fprintf ( 1, '  Bessel function In for NONINTEGER order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '      ALPHA         X             FX                        FX2\n' );
  fprintf ( 1, ...
    '                                  (table)                   (RIBESL)\n' );
  fprintf ( 1, '\n' );
  n_data = 0;

  while ( 1 )

    [ n_data, alpha, x, fx ] = bessel_ix_values ( n_data );

    if ( n_data == 0 );
      break
    end

    ize = 1;
    nb = floor ( alpha ) + 1;
    alpha_frac = alpha - floor ( alpha );

    [ b, ncalc ] = ribesl ( x, alpha_frac, nb, ize );
    fx2 = b(nb);

    fprintf ( 1, '  %12f  %12f  %24.16e  %24.16e\n', alpha, x, fx, fx2 );

  end

  return
end


text:
%EXPORT_PLY Export geometry in PLY format.
%
%   EXPORT_PLY( FILE_NAME, GEOM ) Export of 3D CAD geometries in ASCII
%   Polygon File Format/Stanford Triangle Format. FILE_NAME is a
%   string which specifies the file name to process. GEOM is the
%   geometry to export.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function Fitness = CalFitness(C,Population)
% Calculate the fitness of each solution in terms of a single level

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    PopObj = Population.objs;
    PopCon = Population.cons;
    if any(isnan(PopObj(:,1)))  % Lower level
        PopObj = PopObj(:,2);
        PopCon = PopCon(:,C+1:end);
    else                        % Upper level
        PopObj = PopObj(:,1);
        PopCon = PopCon(:,1:C);
    end
    if isempty(PopCon)
        PopCon = zeros(size(PopObj,1),1);
    else
        PopCon = sum(max(0,PopCon),2);
    end
    Feasible = PopCon <= 0;
    Fitness  = Feasible.*PopObj + ~Feasible.*(PopCon+1e10);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function out = convert_caffe2mat( out )
    assert(length(size(out)) <= 4, 'Only support at most 4-D data for convert.');
    out = single(permute(out, [2 1 3 4]));
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function [imageCha, ssim_map, metrics] = algo_SB_proxA_2D_real( obj,input )
% 2D variant of SB
% based on Goldstein and Oshers paper on Split Bregman
%
% input:
% obj           CS reconstruction object (holding all parameters)
% input         struct containing recon parameters and image
%
% output:
% imageCha      reconstructed channel individual image
% ssim_map      structural similarity map
% metrics       evaluation metrics
%
% (c) Marc Fischer, Thomas Kuestner, May 2015
% -------------------------------------------------------------------------

%%
timer_proxA = tic;

%% variables:
% internal flags
flag_wavetree = true;
flag_fast = true;
flag_extendImage = true;

% internal variables:
L = 1;
t_old = 1;
NLTV_struct.kernelratio = 3;
NLTV_struct.windowratio = 6;
NLTV_struct.nThreads = 1;  % mind this option if used with -singleCompThread on BWCluster
itrNLTV = obj.maxitr - 20;
% chambolle tv:
parsin.MAXITER=100; parsin.tv='iso'; % 'iso' or 'l1'

% from obj:
% maxitr = obj.maxitr;
maxitr = obj.iNINNER;
% maxitrOUTER = obj.maxitrOUTER;
maxitrOUTER = obj.iNOUTER;
maxitrINNER = ceil( maxitr / maxitrOUTER );
% n1 = obj.measPara.dim(1);
% n2 = obj.measPara.dim(2);
n1 = input.n1;
n2 = input.n2;
% nSlices = obj.measPara.dim(3);
nCha = obj.measPara.dim(5);
lambdaWave = obj.lambda;
lambdaTV = obj.lambdaTV;
lambdaGroup = obj.lambdaGroup;
NLTV_struct.filterstrength = obj.lambdaNLTV_h; % 0.03 % converted from NLTV h = 0.01 %old: used: 3e-10
lambdaNLTV = obj.lambdaNLTV;
lambdaNLTV_h = obj.lambdaNLTV_h;
mue = obj.mue;
regularizerWeights = obj.regularizerWeights;
flagTV = obj.flagTV;
flagTV_iso = obj.flagTV_iso;
flagWave = obj.flagWave;
flagGroup = obj.flagGroup;
flagNLTV = obj.flagNLTV;
flagSBNLTV = obj.flagSBNLTV;
waveletStages = obj.trafo.waveletStages;
waveletFilterName_l1 = obj.trafo.waveletFilterName_l1;
waveletFilterName_l12 = obj.trafo.waveletFilterName_l12;

% from input:
b=input.b;
mask = input.mask;
G_prox = input.G_prox;
Gt_prox = input.Gt_prox;
groupnorm_index = input.groupnorm_index;
waveS_l1 = input.waveS_l1;
waveS_l12 = input.waveS_l12;
waveS_l12proxA = input.waveS_l12proxA;
proxA_extend_y = waveS_l12proxA(waveletStages+2,1) - waveS_l12(waveletStages+2,1);
proxA_extend_x = waveS_l12proxA(waveletStages+2,2) - waveS_l12(waveletStages+2,2);
x_proxA = cell(1,nCha);

% im_ref = input.im_ref;
% im_ref_full = zeros(n1,n2);
% for j = 1:nCha
%     im_ref_full = im_ref_full + abs(im_ref{1,j}).^2;
% end;
% im_ref_full = sqrt(im_ref_full);
clear input

% initialize cells/vectors
for j=1:nCha
    FTb{1,j} = real(iFFT2D(b{1,j}));
    % y{1,j} = real(FTb{1,j});
    % y{1,j+nCha} = imag(FTb{1,j});
end;
% starting point:
x = FTb; % (x0 = FTb)
rhs = FTb; % different start (x0 = F^-1 K^-1 F FTb); %used atm
b_outer = b;
% rhs: right hand side

x_wave =  cell(1,nCha);
x_g_proxA = x_wave;
x_g_helper = cell(2,nCha);
x_wave_helper = rhs;
v_wave = rhs;
d_wave = rhs;
x_tv = rhs;
v_tv = rhs;
d_tv = rhs;
x_group = rhs;
v_group = rhs;
d_group = rhs;
x_nltv = rhs;
v_nltv = rhs;
d_nltv = rhs;
K_inv = rhs;
for j = 1:nCha
    x_wave{1,j} = 0;
    v_wave{1,j} = 0;
    d_wave{1,j} = 0;
    x_tv{1,j} = zeros(n1,n2);
    v_tv{1,j} = 0;
    d_tv{1,j} = 0;
    x_group{1,j} = 0;
    v_group{1,j} = 0;
    d_group{1,j} = 0;
    x_nltv{1,j} = 0;
    v_nltv{1,j} = 0;
    d_nltv{1,j} = 0;
end;

%% MAD dependent lambdas:
flag_MAD = true;
if flag_MAD
    x_wavedec = cell(1,nCha);
    threshold = zeros(1:nCha);
    for j=1:nCha % 2*nCha
        x_wavedec{1,j} = wavedec2(x{1,j},waveletStages,waveletFilterName_l1); % atm only based on l1-daubechie
        x_wave_fine_scale = size(x_wavedec{1,j},2) - (3*waveS_l1(waveletStages+1,1)*waveS_l1(waveletStages+1,2));
        threshold(j) = mad(x_wavedec{1,j}(x_wave_fine_scale:end),1);
    end;
    clear x_wavedec
else
    threshold(1:nCha) = 1;
end;

threshold_wave(j) = lambdaWave * threshold(j) * 2/L;
threshold_TV(j) = lambdaTV * threshold(j) * 2/L;
threshold_group(j) = lambdaGroup * threshold(j) * 2/L;
threshold_NLTV(j) = lambdaNLTV; % * threshold(j) * 2/L;
threshold_NLTV_h(j) = lambdaNLTV_h; %*threshold(j); % adjust carefully or NLTV won't find a solution. lambdaNLTV_h should be < 0.01

%% initialize metrics:
itr = 0;
metrics.xtime(itr+1)= 0;
% [metrics, ssim_map{1,1}] = get_metrics_itr( im_ref, im_ref_full, x, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );
ssim_map = [];

%% recon
% inverse for x - step:
for j = 1:nCha
    K_inv{1,j} = 1./(mue*mask + (flagWave*regularizerWeights(1) + flagTV*regularizerWeights(2) +  flagGroup*regularizerWeights(3))*ones(n1,n2)); % flagNLTV changes K_inv -> in separate
end;

dispProgress('Proximal Average', 0, maxitrOUTER*maxitrINNER);
itr = 0;
for itrOUTER = 1:maxitrOUTER  % total iter counter
    for itrINNER = 1:maxitrINNER

        %% x - step
        if itr > 1 % for i = 1 already set (x{1,j} = FTb{1,j})
            for j = 1:nCha
                % rhs part:
                if flagWave
                    x_wave{1,j} = waverec2((v_wave{1,j} - d_wave{1,j}),waveS_l1,waveletFilterName_l1); end;
                if flagTV
                    x_tv{1,j} = v_tv{1,j} - d_tv{1,j}; end;
                if flagGroup
                    x_g_proxA{1,j} = waverec2((v_group{1,j} - d_group{1,j}),waveS_l12proxA,waveletFilterName_l12);
                    x_group{1,j} =  x_g_proxA{1,j}(1:end-proxA_extend_y,1:end-proxA_extend_x);
                end;
                if flagNLTV
                    if itr == itrNLTV
                        K_inv{1,j} = 1./(mue*mask + (flagWave*regularizerWeights(1) + flagTV*regularizerWeights(2) +  flagGroup*regularizerWeights(3) + flagNLTV*regularizerWeights(4))*ones(n1,n2));
                        x_nltv{1,j} = x{1,j};
                    elseif itr >= itrNLTV
                        x_nltv{1,j} = v_nltv{1,j} - d_nltv{1,j};
                    end;
                end;
                rhs{1,j} = mue*FTb{1,j} + ( flagWave*regularizerWeights(1).*(x_wave{1,j}) + flagTV*regularizerWeights(2).*(x_tv{1,j}) + flagGroup*regularizerWeights(3).*(x_group{1,j}) + flagNLTV*regularizerWeights(4).*(x_nltv{1,j}));
                x{1,j} = real(iFFT2D(K_inv{1,j}.*FFT2D(rhs{1,j})));
            end;
        end;

        %% auxiliary variables:
        %% l1-Wavelet
        if flagWave
            for j = 1:nCha
                x_wave_helper{1,j} = wavedec2(x{1,j},waveletStages,waveletFilterName_l1);
                v_wave{1,j} = softthresh_real((x_wave_helper{1,j} + d_wave{1,j}),threshold_wave(j));
                d_wave{1,j} = d_wave{1,j} + x_wave_helper{1,j} - v_wave{1,j};
            end;
        end;

        %% TV:
        if flagTV
            if ~flagTV_iso
                for j = 1:nCha
                    v_tv{1,j} = MTV_2D((x{1,j} + d_tv{1,j}), threshold_TV(j), n1, n2);
                end;
            else
                for j = 1:nCha
                    if (itr < 2)
                        [v_tv{1,j}, P]=denoise_TV_One((x{1,j} + d_tv{1,j}), threshold_TV(j),-inf,inf,[],parsin);
                    else
                        [v_tv{1,j}, P]=denoise_TV_One((x{1,j} + d_tv{1,j}), threshold_TV(j),-inf,inf,P,parsin);
                    end;
                end;
            end;
            d_tv{1,j} = d_tv{1,j} + x{1,j} - v_tv{1,j};
        end;

        %% NLTV
        if flagNLTV
            if itr >= itrNLTV
                if flagSBNLTV
                    for j = 1:nCha
                        v_nltv{1,j} = SB_NLTVfunc_slim_rescale( (x{1,j} + d_nltv{1,j}),n1,n2, threshold_NLTV(j), threshold_NLTV_h(j) );
                    end;
                else
                    for j = 1:nCha
                        % if mod(itr,5) == 0 || itr == 1
                        v_nltv{1,j} = NLMF((x{1,j} + d_nltv{1,j}),NLTV_struct);
                        v_nltv{1,j} = (L.*(x{1,j} + d_nltv{1,j}) + 2*threshold_NLTV(j)*v_nltv{1,j})./(L+2*threshold_NLTV(j));
                        % end;
                    end;
                end;
                d_nltv{1,j} = d_nltv{1,j} + x{1,j} - v_nltv{1,j};
            else
                v_nltv{1,j} = 0;
                d_nltv{1,j} = 0;
            end;

        end;

        %% l12-Wavelet
        if flagGroup
            for j = 1:nCha
                if flag_extendImage
                    x_proxA{1,j} = extend_image(x{1,j}, waveS_l12proxA, waveletStages, proxA_extend_y, proxA_extend_x);
                    x_g_helper{1,j} = wavedec2(x_proxA{1,j},waveletStages,waveletFilterName_l12) + d_group{1,j};
                else
                    x_g_helper{1,j} = wavedec2(x{1,j},waveletStages,waveletFilterName_l12) + d_group{1,j};
                end;

                if flag_wavetree
                    x_g_helper{2,j} = (G_prox*x_g_helper{1,j}')';
                else
                    x_g_helper{2,j} = zeros(1,size(x_g_helper{1,j},2));
                end;
            end;

            % softthresh_proxA_cha returns cell! hence (1,:) assignement
            v_group(1,:) = softthresh_proxA_cha(x_g_helper,threshold_group(j),nCha,waveS_l12proxA,groupnorm_index);

            for j = 1:nCha
                d_group{1,j} = d_group{1,j} + x_g_helper{1,j} - v_group{1,j};
            end;
        end;

        %% metrics of current itr:
        itr = itr + 1;
%         disp(itr);
        dispProgress('Proximal Average', itr/(maxitrOUTER*maxitrINNER));

        metrics.xtime(itr+1)= toc(timer_proxA);
%         [metrics, ssim_map{1,2}] = get_metrics_itr( im_ref, im_ref_full, x, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );

        if itr >= maxitr break; end;
    end;
    if itr >= maxitr break; end;
    for j=1:nCha
        b_outer{1,j} = b_outer{1,j} + b{1,j} - FFT2D_mask(x{1,j},mask);
        FTb{1,j} = real(iFFT2D(b_outer{1,j}));
    end;
end;
dispProgress('Proximal Average', 'Close');

imageCha = x;
for j = 1:nCha
    imageCha{1,j} = turn_image( imageCha{1,j} );
end;
% for j = 1:nCha+1
%     ssim_map{1,1}{1,j} = turn_image( ssim_map{1,1}{1,j} );
%     ssim_map{1,2}{1,j} = turn_image( ssim_map{1,2}{1,j} );
% end;

end

author: thomaskuestner


text:
function out=leadsphere_chans(xloc,sensorloc,sensorori)
% usage: out=leadsphere_chans(xloc,sensorloc,sensorori)

% Copyright (C) 2003, Guido Nolte
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: leadsphere_all.m 2885 2011-02-16 09:41:58Z roboos $

[n,nsens]=size(sensorloc); %n=3 m=?
[n,ndip]=size(xloc);

xlocrep=reshape(repmat(xloc,1,nsens),3,ndip,nsens);
sensorlocrep=reshape(repmat(sensorloc,ndip,1),3,ndip,nsens);
sensororirep=reshape(repmat(sensorori,ndip,1),3,ndip,nsens);

r2=norms(sensorlocrep);
veca=sensorlocrep-xlocrep;
a=norms(veca);
adotr2=dotproduct(veca,sensorlocrep);

gradf1=scal2vec(1./r2.*(a.^2)+adotr2./a+2*a+2*r2);
gradf2=scal2vec(a+2*r2+adotr2./a);
gradf=gradf1.*sensorlocrep-gradf2.*xlocrep;

F=a.*(r2.*a+adotr2);

A1=scal2vec(1./F);
A2=A1.^2;

A3=crossproduct(xlocrep,sensororirep);
A4=scal2vec(dotproduct(gradf,sensororirep));
A5=crossproduct(xlocrep,sensorlocrep);

out=1e-7*(A3.*A1-(A4.*A2).*A5); %%GRB change

return;

function out=crossproduct(x,y)
[n,m,k]=size(x);
out=zeros(3,m,k);
out(1,:,:)=x(2,:,:).*y(3,:,:)-x(3,:,:).*y(2,:,:);
out(2,:,:)=x(3,:,:).*y(1,:,:)-x(1,:,:).*y(3,:,:);
out(3,:,:)=x(1,:,:).*y(2,:,:)-x(2,:,:).*y(1,:,:);
return;

function out=dotproduct(x,y)
[n,m,k]=size(x);
outb=x(1,:,:).*y(1,:,:)+x(2,:,:).*y(2,:,:)+x(3,:,:).*y(3,:,:);
out=reshape(outb,m,k);
return;

function result=norms(x)
[n,m,k]=size(x);
resultb=sqrt(x(1,:,:).^2+x(2,:,:).^2+x(3,:,:).^2);
result=reshape(resultb,m,k);
return;

function result=scal2vec(x)
[m,k]=size(x);
% result=zeros(3,m,k);
% for i=1:3
%     result(i,:,:)=x;
% end
result=reshape(repmat(x(:)', [3 1]), [3 m k]);
return

author: PatternRecognition

repo: OpenBMI


text:
function percent = parfor_progress(N)
%PARFOR_PROGRESS Progress monitor (progress bar) that works with parfor.
%   PARFOR_PROGRESS works by creating a file called parfor_progress.txt in
%   your working directory, and then keeping track of the parfor loop's
%   progress within that file. This workaround is necessary because parfor
%   workers cannot communicate with one another so there is no simple way
%   to know which iterations have finished and which haven't.
%
%   PARFOR_PROGRESS(N) initializes the progress monitor for a set of N
%   upcoming calculations.
%
%   PARFOR_PROGRESS updates the progress inside your parfor loop and
%   displays an updated progress bar.
%
%   PARFOR_PROGRESS(0) deletes parfor_progress.txt and finalizes progress
%   bar.
%
%   To suppress output from any of these functions, just ask for a return
%   variable from the function calls, like PERCENT = PARFOR_PROGRESS which
%   returns the percentage of completion.
%
%   Example:
%
%      N = 100;
%      parfor_progress(N);
%      parfor i=1:N
%         pause(rand); % Replace with real code
%         parfor_progress;
%      end
%      parfor_progress(0);
%
%   See also PARFOR.

% By Jeremy Scheff - jdscheff@gmail.com - http://www.jeremyscheff.com/

%error(nargchk(0, 1, nargin, 'struct'));
narginchk(0, 1); %because nargchk is depreciated

if nargin < 1
    N = -1;
end

percent = 0;
w = 50; % Width of progress bar

if N > 0
    f = fopen('parfor_progress.txt', 'w');
    if f<0
        error('Do you have write permissions for %s?', pwd);
    end
    fprintf(f, '%d\n', N); % Save N at the top of progress.txt
    fclose(f);

    if nargout == 0
        disp(['  0%[>', repmat(' ', 1, w), ']']);
    end
elseif N == 0
    delete('parfor_progress.txt');
    percent = 100;

    if nargout == 0
        disp([repmat(char(8), 1, (w+9)), char(10), '100%[', repmat('=', 1, w+1), ']']);
    end
else
    if ~exist('parfor_progress.txt', 'file')
        error('parfor_progress.txt not found. Run PARFOR_PROGRESS(N) before PARFOR_PROGRESS to initialize parfor_progress.txt.');
    end

    f = fopen('parfor_progress.txt', 'a');
    fprintf(f, '1\n');
    fclose(f);

    f = fopen('parfor_progress.txt', 'r');
    progress = fscanf(f, '%d');
    fclose(f);
    percent = (length(progress)-1)/progress(1)*100;

    if nargout == 0
        perc = sprintf('%3.0f%%', percent); % 4 characters wide, percentage
        disp([repmat(char(8), 1, (w+9)), char(10), perc, '[', repmat('=', 1, round(percent*w/100)), '>', repmat(' ', 1, w - round(percent*w/100)), ']']);
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%  DemoTV_UP.m
%
%  This is a short script that performs TV minimization recovery
%   from partial Fourier measurements
%
%  min Lambda ||x||_TV +1/2||b - A x||_2
%
%  Here A*A is assumed to be an orthogonal projector
%
%  Parameters to set : 1) n : the size of the image to reconstruct (default: 128)
%                      2) Dyna : the dynamic range of the original image (in dB) (default : 40)
%                      3) Lambda : value of the Lagrange multiplier (default: 0.05)
%
%  Created entities : I : original image
%                     Xnesta : image recovered with Nesta
%                     NA_nesta = number of calls of A or A* for the method
%
% Written by: Jerome Bobin, Caltech
% Email: bobin@acm.caltech.edu
% Created: May 2009
%
% NESTA Version 1.1
%   See also NESTA_UP and Core_Nesterov_UP
clear; clc;
n = 128;   %--- The data are n*n images
Dyna = 40; %--- Dynamic range in dB
Lambda = 0.05;
Setup_Nesta % -- add the solvers to the path
%% Setting up the experiment

N = n*n;
I=MakeRDSquares(n,7,Dyna);
x = reshape(I,n*n,1);
L = floor(55*(n/256)); %--- so as to get a compression ratio ~10
[M,Mh,mi,mhi] = LineMask(L,n);
OMEGA = mhi;
OMEGA = [OMEGA];
K = length(OMEGA);
A = @(z) A_fhp(z,OMEGA);
At = @(z) At_fhp(z,OMEGA,n);
b0 = A(x);
sigma = 0.1;
noise = sigma*randn(size(b0));
b = b0+noise;

fprintf('##############################################\n\n');
fprintf('NESTA: Total Variation minimization experiment\n\n');
fprintf('##############################################\n\n');
fprintf('Image size = %g x %g / Dynamic range : %g / Noise level : %g\n\n',n,n,Dyna,sigma);

%% Setting up the experiment for NESTA

fprintf(' NESTA starts\n\n');
U = @(z) z;
Ut = @(z) z;
mu = 0.2; %--- can be chosen to be small
opts = [];
opts.maxintiter = 5;
opts.TOlVar = 1e-4;
opts.verbose = 25;
opts.maxiter = 5000;
opts.U = U;
opts.Ut = Ut;
opts.stoptest = 1;
opts.typemin = 'tv';
counter();
Ac = @(z) counter(A,z);
Atc = @(z) counter(At,z);
tic;
[x_nesta,niter,resid,err] = NESTA_UP(Ac,Atc,b,Lambda,1,mu,opts);
t.NESTA = toc;
NA_nesta = counter();
tvnesta = calctv(n,n,reshape(x_nesta,n,n));
Xnesta = reshape(x_nesta,n,n);

% -- Print results

fprintf('\nNESTA: TV: %g -- ||b-Ax|| = %g -- # calls : %g\n\n',tvnesta,norm(b - A(Xnesta(:))),NA_nesta);
fprintf('relative l2 error compared to original signal: %.2e\n',norm(Xnesta-I,'fro')/norm(I,'fro') );
%% Plot results
figure(1); clf;

subplot(1,2,1);
image( I );  title('Original image');
map = colormap(hsv);
map(1,:) = [1,1,1];
colormap(map);
% colorbar('westoutside');
axis square; axis off

subplot(1,2,2);
image(Xnesta); title(sprintf(...
    'Reconstruction via NESTA\n%.1f sec',t.NESTA));
axis square; axis off

repo: REST


text:
classdef pfSections < ODFSections

  properties
    h1        % the pole figure which is splitted up
    h2        %
    omega
    sR
    referenceField
  end

  properties (Hidden=true)
    maxOmega
  end

  methods

    function oS = pfSections(CS1,CS2,varargin)

      oS = oS@ODFSections(CS1,CS2);

      if isa(CS1,'crystalSymmetry')
        oS.h1 = CS1.cAxisRec; % c*
        oS.h2 = CS1.aAxis; % a
      else
        oS.h1 = vector3d.Z;
        oS.h2 = vector3d.X;
      end

      oS.maxOmega = get_option(varargin,'maxOmega',2*pi / CS1.nfold(oS.h1));
      if angle(oS.h1,-oS.h1) < 1e-2
        oS.sR = CS2.fundamentalSector('upper',varargin{:});
      else
        oS.sR = CS2.fundamentalSector(varargin{:});
      end

      % get sections
      oS.omega = linspace(0,oS.maxOmega,1+get_option(varargin,'sections',6));
      oS.omega(end) = [];
      oS.omega = get_option(varargin,'omega',oS.omega,'double');

      oS.updateTol(oS.omega);

      oS.referenceField = S2VectorField.sigma;
      %oS.referenceField = @(r) S2VectorField.oneSingularity;
      %oS.referenceField = @(r) pfSections.polarField(r);

    end

    function ori = makeGrid(oS,varargin)

      oS.plotGrid = plotS2Grid(oS.sR,varargin{:});
      oS.gridSize = (0:numel(oS.omega)) * length(oS.plotGrid);

      ori = orientation.nan(oS.plotGrid.size(1),oS.plotGrid.size(2),numel(oS.omega),oS.CS1,oS.CS2);
      for iOmega = 1:numel(oS.omega)

        r2 = oS.vectorField(oS.plotGrid,oS.omega(iOmega));
        ori(:,:,iOmega) = reshape(orientation.map(oS.h1,oS.plotGrid,oS.h2,r2),size(oS.plotGrid));

      end

    end

    function ori = quiverGrid(oS,varargin)

      S2G = equispacedS2Grid(oS.sR,varargin{:});

      ori = orientation.nan(S2G.size(1),S2G.size(2),numel(oS.omega),oS.CS1,oS.CS2);
      for iOmega = 1:numel(oS.omega)

        r2 = oS.vectorField(S2G,oS.omega(iOmega));
        ori(:,:,iOmega) = reshape(orientation.map(oS.h1,S2G,oS.h2,r2),size(S2G));

      end
    end

    function n = numSections(oS)
      n = numel(oS.omega);
    end

    function [r,secPos] = project(oS,ori,varargin)

      % maybe this can be done more efficiently
      if ~check_option(varargin,'noSymmetry')
        ori = ori.symmetrise('proper').';
      end

      % determine pole figure position
      r = ori * oS.h1;

      % determine omega angle
      rF = ori * oS.h2;
      vF = vectorField(oS,r);
      omegaSec = angle(vF,rF,r);

      % determine section number
      secPos = oS.secList(omegaSec,oS.omega);

    end

    function ori = iproject(oS,rho,theta,iOmega)
      r1 = vector3d.byPolar(theta,rho);
      r2 = oS.vectorField(r1,oS.omega(iOmega));

      ori = orientation.map(oS.h1,r1,oS.h2,r2);
    end

    function h = plotSection(oS,ax,sec,v,data,varargin)

      % plot data
      h = plot(v,data{:},oS.sR,'TR',[int2str(oS.omega(sec)./degree),'^\circ'],...
        'parent',ax,varargin{:},'doNotDraw');

      wasHold = ishold(ax);
      hold(ax,'on');
      r = equispacedS2Grid(oS.sR,'resolution',15*degree);
      vF = oS.vectorField(r,oS.omega(sec));
      quiver(r,vF,'parent',ax,'doNotDraw','color',0.7*[1 1 1],'HitTest','off');
      if ~wasHold, hold(ax,'off'); end

    end

    function h = quiverSection(oS,ax,sec,v,data,varargin)

      % translate rotational data into tangential data
      if iscell(data) && isa(data{1},'quaternion')
        [v2,sec2] = project(oS,data{1},'noSymmetry');
        data{1} = v2 - v;
        data{1}(sec2 ~= sec)=NaN;
      end
      if check_option(varargin,'normalize'), data{1} = normalize(data{1}); end

      % plot data
      h = quiver(v,data{:},oS.sR,'TR',[int2str(oS.omega(sec)./degree),'^\circ'],...
        'parent',ax,varargin{:},'doNotDraw');

    end

    function vF = vectorField(oS,r,omega)

      vF = oS.referenceField.eval(r);

      if nargin == 3
        vF = rotation.byAxisAngle(r,omega) .* vF;
      end

    end

  end

end

% testing code
% r = equispacedS2Grid('resolution',20*degree)
% quiver(r,vector3d(vF))

author: mtex-toolbox

repo: mtex

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


text:
function [vol,dist] = volume(psi,radius,dist)

if nargin < 3
  N = 60;
  dist = zeros(numel(radius),N);
  for i = 1:length(radius)

    if radius(i) + 5 * psi.halfwidth < pi
      dist(i,1:N-9) = linspace(0,radius(i) + 5 * psi.halfwidth,N-9);
      dist(i,N-9:end) = linspace(radius(i) + 5 * psi.halfwidth,pi,10);
    else
      dist(i,:) = linspace(0,pi,N);
    end

  end

end

% make radius square as well
radius = repmat(radius(:),1,size(dist,2));

% weight function
weight = @(rot_angle,quat_dist,volume_radius) (max(0,min(4*pi,2*pi*(1-(cos(volume_radius/2) - ...
  cos(quat_dist/2) * cos(rot_angle/2))./...
  (sin(quat_dist/2) * sin(rot_angle/2))))) ...
  + max(0,min(4*pi,2*pi*(1-(cos(volume_radius/2) + ...
  cos(quat_dist/2) * cos(rot_angle/2))./...
  (sin(quat_dist/2) * sin(rot_angle/2))))))./pi./4;

% integrant
KV = @(rot_angle,quat_dist,volume_radius) weight(rot_angle,quat_dist,volume_radius) ...
  .* sin(rot_angle ./2).^2 .* psi.eval(cos(rot_angle./2)) ./pi.*2 ;

% perform quadrature
%vol = zeros(size(dist));
%for j = 1:length(dist)
%  vol(j) = quad(@(rot_angle) KV(rot_angle,dist(j),radius),0,min(pi,5*k.hw),1e-6);
%end

%vol = quadv(@(rot_angle) KV(rot_angle,dist,radius),0,min(pi,5*psi.halfwidth),1e-7);
vol = integral(@(rot_angle) KV(rot_angle,dist,radius),0,min(pi,5*psi.halfwidth),'AbsTol',1e-7,'ArrayValued',true);


text:
function [eta,etaf]=finishat(frac,tol,fmt)
%FINISHAT print estimated finish time of a long computation (FRAC,TOL,FMT)
% Usage: (1)  finishat(0);
%              for i=1:many
%                  ... computation ...
%                  finishat(i/many);
%              end
%
%         (2)  finishat(0);
%              for i=1:NI
%                  for j=1:NJ
%                      for k=1:NK
%                          ... computation ...
%                          finishat((((i-1)*NJ+j-1)*NK+k)/(NI*NJ*NK));
%                      end
%                  end
%              end
%
% Inputs: FRAC = fraction of total comutation that has been completed
%                As a special case, 0 is used to initialize the routine
%         TOL  = Tolerance in minutes. If the estimated time has changed by less
%                than this, then nothing will be printed. [default 10% of remaining time]
%         FMT  = Format string which should include %s for estimated finish time, %d for remaining minutes and %f for fraction complete
%
% Output: ETA  = string containing the expected finish time
%                specifying this will suppress printing message to std err (fid=2)
%         ETAF = expected finish time as a daynumber
%
% Example:       finishat(0);
%                for i=1:many
%                    long computation;
%                    finishat(i/many);
%                end

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

persistent oldt oldnw
if nargin<3
    fmt='Estimated finish at %s (%.2f done, %d min remaining)\n';

end
if frac<=0
    oldt=0;
    eta='Unknown';
    tic;
else
    nw=now;                             % current time as serial number
    sectogo=(1/frac-1)*toc;      % seconds to go
    newt=nw+sectogo/86400;       % add estimated time in days
    if nargin<2 || ~numel(tol)
        tol=max(0.1*(newt-nw)*1440,1);
    end
    if ~exist('oldt','var') || oldt==0 || (abs(newt-oldt)>tol/1440 && (nw-oldnw)>10/86400) || (nw-oldnw)>10/1440 || nargout>0
        oldt=newt;
        if floor(oldt)==floor(nw)
            df='HH:MM';
        else
            df='HH:MM dd-mmm-yyyy';
        end
        eta=datestr(oldt,df);
        if ~nargout
            ix=find(fmt=='%',1);
            while ~isempty(ix)
                fprintf(2,fmt(1:ix-1));
                fmt=fmt(ix:end);
                ix=find(fmt>='a' & fmt<='z',1); % find letter
                switch fmt(ix)
                    case 's'
                        fprintf(2,fmt(1:ix),eta);
                    case 'd'
                        fprintf(2,fmt(1:ix),round(sectogo/60));
                    case 'f'
                        fprintf(2,fmt(1:ix),frac);
                end
                fmt=fmt(ix+1:end);
                ix=find(fmt=='%',1);
            end
            fprintf(2,fmt);
        end
        oldnw=nw;                           %
    end
end
etaf=oldt;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function Nstar = ctranspose(N)
%'   Compute the adjoint of a CHEBOP.
%   NSTAR = N' returns the adjoint of a CHEBOP that has either periodic or
%   endpoint boundary conditions. If N is nonlinear then N is first
%   linearized around U = 0 and NSTAR is the adjoint of the linearization.
%
%   N' is calculated by a call to ADJOINT(N).  To linearize about
%   a nonzero function U, use ADJOINT(N, U).
%
%  Example:
%   L = chebop(-1,1); L.op = @(x,u) diff(u,2) + diff(u,1) + u;
%   L.lbc = 0; L.rbc = 1, Ls = L'
%
% See also CHEBOP/ADJOINT.

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

% This is just a wrapper for adjoint().
Nstar = adjoint(N);

end


text:
classdef CoreTextFileWriter < CoreBaseClass
    % CoreTextFileWriter. A helper class to assist with writing a text file
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        FileText
        FileId
    end

    methods
        function obj = CoreTextFileWriter(file_path, file_name, reporting)
            full_file_name = fullfile(file_path, file_name);
            file_id = fopen(full_file_name, 'w');

            if (file_id == -1)
                reporting.Error('CoreTextFileWriter:CannotOpenFile', ['Unable to open file ' full_file_name]);
            else
                obj.FileId = file_id;
            end
        end

        function Close(obj)
            if ~isempty(obj.FileId)
                fclose(obj.FileId);
                obj.FileId = [];
            end
        end

        function delete(obj)
            obj.Close;
        end

        function WriteLine(obj, text)
            fprintf(obj.FileId, [strrep(text, '%', '%%') '\r\n']);
        end

        function is_eof = Eof(obj)
            is_eof = isempty(obj.FileId);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function outdata = jatmos70(indata)

% Jacchia 1970 atmosphere main driver

% input

%  indata(1)  = geodetic altitude (kilometers)
%  indata(2)  = geodetic latitude (radians)
%  indata(3)  = geographic longitude (radians)
%  indata(4)  = calendar year (all digits)
%  indata(5)  = calendar month
%  indata(6)  = calendar day
%  indata(7)  = utc hours
%  indata(8)  = utc minutes
%  indata(9)  = geomagnetic index type
%               (1 = indata(12) is Kp, 2 = indata(12) is Ap)
%  indata(10) = solar radio noise flux (jansky)
%  indata(11) = 162-day average F10 (jansky)
%  indata(12) = geomagnetic activity index

% output

%  outdata(1)  = exospheric temperature (deg K)
%  outdata(2)  = temperature at altitude (deg K)
%  outdata(3)  = N2 number density (per meter-cubed)
%  outdata(4)  = O2 number density (per meter-cubed)
%  outdata(5)  = O number density (per meter-cubed)
%  outdata(6)  = A number density (per meter-cubed)
%  outdata(7)  = He number density (per meter-cubed)
%  outdata(8)  = H number density (per meter-cubed)
%  outdata(9)  = average molecular weight
%  outdata(10) = total density (kilogram/meter-cubed)
%  outdata(11) = log10(total density)
%  outdata(12) = total pressure (pascals)

% Orbital Mechanics with MATLAB

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

bfh = 440;

rgas = 8314.32;

% unload input data

z = indata(1);
xlat = indata(2);
xlng = indata(3);
iyr = indata(4);
mn = indata(5);
ida = indata(6);
ihr = indata(7);
min = indata(8);
i1 = indata(9);
f10 = indata(10);
f10b = indata(11);
gi = indata(12);

% compute solar coordinates

[sda, sha, dd, dy] = jtme(mn, ida, iyr, ihr, min, xlng);

% calculate exospheric temperature

te = jtinf(f10, f10b, gi, xlat, sda, sha, dy, i1);

% evaluate jacchia model

[dens, dl, em, tz, a(1), a(2), a(3), a(4), a(5), a(6)] = jacchia(z, te);

denlg = 0;
dummy = dl;
den = dl;

if (z <= 170)
   dummy = jslv(z, xlat, dd);
   denlg = dummy;
end

if (z >= 500)
    den = jslvh(den, a(5), xlat, sda);
    dl = den;
elseif (z > bfh)
    dhel1 = a(5);
    dhel2 = a(5);
    dlg1 = dl;
    dlg2 = dl;
    dlg2 = jslvh(dlg2, dhel2, xlat, sda);
    ih = z;
    [fdhel, fdlg] = jfair(dhel1, dhel2, dlg1, dlg2, ih);
    dl = fdlg;
    a(5) = fdhel;
end

dl = dl + denlg;

dens = 10 ^ dl;

% load output data array

outdata(1) = te;

outdata(2) = tz;

for i = 1:1:6
    outdata(i + 2) = 1000000 * (10 ^ a(i));
end

outdata(9) = em;
outdata(10) = dens * 1000;
outdata(11) = dl;

p = outdata(10) * rgas * tz / em;

outdata(12) = p;

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

function [dens, dl, em, tz, an, ao2, ao, aa, ahe, ah] = jacchia (z, t)

% jacchia 1970 atmosphere

global alpha ei

av = 6.02257e+23;
qn = 0.7811;
qo2 = 0.20955;
qa = 9.343e-03;
qhe = 0.00001289;
rgas = 8.31432;
t0 = 183;

tx = 444.3807 + 0.02385 * t - 392.8292 * exp(-0.0021357 * t);
a2 = 2 * (t - tx) / pi;
txt0 = tx - t0;
t1 = 1.9 * txt0 / 35;
t3 = -1.7 * txt0 / (35 ^ 3);
t4 = -.8 * txt0 / (35 ^ 4);
tz = jtemp(z, tx, t1, t3, t4, a2);

a = 90;
d = min(z, 105);

r = jgauss(a, d, 1, tx, t1, t3, t4, a2);

em = jmweight(d);

td = jtemp(d, tx, t1, t3, t4, a2);

dens = 0.000000021926 * em * exp(-r / rgas) / td;

factor = av * dens;
par = factor / em;
factor = factor / 28.96;

if (z <= 105)
   dl = log10(dens);
   an = log10(qn * factor);
   aa = log10(qa * factor);
   ahe = log10(qhe * factor);
   ao = log10(2 * par * (1 - em / 28.96));
   ao2 = log10(par * (em * (1 + qo2) / 28.96 - 1));
   ah = 0;
   return;
end

di(1) = qn * factor;
di(2) = par * (em * (1 + qo2) / 28.96 - 1);
di(3) = 2 * par * (1 - em / 28.96);
di(4) = qa * factor;
di(5) = qhe * factor;

r = jgauss(d, z, 2, tx, t1, t3, t4, a2);

for i = 1:1:5
    dit(i) = di(i) * (td / tz) ^ (1 + alpha(i)) * exp(-ei(i) * r / rgas);

    if (dit(i) <= 0)
       dit(i) = 0.000001;
    end
end

if (z > 500)
   a1 = 500;
   s = jtemp(a1, tx, t1, t3, t4, a2);
   di(6) = 10 ^ (73.13 - 39.4 * log10(s) + 5.5 * log10(s) * log10(s));
   r = jgauss(a1, z, 7, tx, t1, t3, t4, a2);
   dit(6) = di(6) * (s / tz) * exp(-ei(6) * r / rgas);
else
   dit(6) = 1;
end

dens = 0;

for i = 1:1:6
    dens = dens + ei(i) * dit(i) / av;
end

em = dens * av / (dit(1) + dit(2) + dit(3) + dit(4) + dit(5) + dit(6));

dl = log10(dens);

an = log10(dit(1));
ao2 = log10(dit(2));
ao = log10(dit(3));
aa = log10(dit(4));
ahe = log10(dit(5));
ah = log10(dit(6));

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

function [sda, sha, dd, dy] = jtme (mn, ida, iyr, ihr, min, xlng)

% solar declination and hour angle

dtr = pi / 180;

atr = dtr / 3600;

% compute julian date on january 1

jdate = julian(1, 1, iyr);

% compute julian date on day of interest

xmjd = julian(mn, ida, iyr);

% compute day number of the year

dd = xmjd - jdate;

% compute fraction of tropical year

dy = dd / 365.2422;

% compute greenwich mean time

gmt = 60 * ihr + min;

% compute solar coordinates

a = xmjd + gmt / 1440 - 2451545;

b = a / 36525 + 1;

o = r2r(0.056531 + 0.00023080893 * a);
m = r2r(0.140023 + 0.00445036173 * a);
l = r2r(0.779072 + 0.00273790931 * a);
h = r2r(0.606434 + 0.0366011013 * a);
n = r2r(0.053856 + 0.00145561327 * a);
f = r2r(0.993126 + 0.0027377785 * a);
r = r2r(0.347343 - 0.00014709391 * a);

obliq = atr * (84428 - 47 * b + 9 * cos(r));

p = 6910 * sin(f) + 72 * sin(2 * f) - 17 * b * sin(f) - 7 * cos(f - o) ...
    + 6 * sin(h - l) + 5 * sin(4 * f - 8 * n + 3 * o) ...
    - 5 * cos(2 * (f - m)) - 4 * sin(f - m) + 4 * cos(4 * f - 8 * n ...
    + 3 * o) + 3 * sin(2 * (f - m)) - 3 * sin(o) - 3 * sin(2 * (f - o));

p = l + atr * (p - 17 * sin(r));

a = sin(p) * cos(obliq);

b = cos(p);

% compute solar right ascension and declination

ras = atan3(a, b);

sda = asin(sin(obliq) * sin(p));

% compute greenwich sidereal time

gst = gast1(xmjd + ihr / 24 + min / 1440);

% compute right ascension of point of interest

rap = mod(gst + xlng, 2.0 * pi);

% compute hour angle

sha = mod(rap - ras, 2.0 * pi);

if (sha > pi)
   sha = sha - 2.0 * pi;
end

if (sha < -pi)
   sha = sha + 2.0 * pi;
end

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

function te = jtinf (f10, f10b, gi, xlat, sda, sha, dy, i1)

% exospheric temperature

c1 = 383;
c2 = 3.32;
c3 = 1.8;

d1 = 28;
d2 = 0.03;
d3 = 1;
d4 = 100;
d5 = -0.08;

e1 = 2.41;
e2 = 0.349;
e3 = 0.206;
e4 = 6.2831853;
e5 = 3.9531708;
e6 = 12.5663706;
e7 = 4.3214352;
e8 = 0.1145;
e9 = 0.5;
e10 = 6.2831853;
e11 = 5.974262;
e12 = 2.16;

beta = -0.6457718;
gamma = 0.7504916;
p = 0.1047198;
re = 0.31;

tc = c1 + c2 * f10b + c3 * (f10 - f10b);

eta = 0.5 * abs(xlat - sda);
theta = 0.5 * abs(xlat + sda);
tau = sha + beta + p * sin(sha + gamma);

if (tau > pi)
   tau = tau - 2 * pi;
end

if (tau < -pi)
   tau = tau + 2 * pi;
end

a1 = (sin(theta)) ^ 2.5;
a2 = (cos(eta)) ^ 2.5;
a3 = (cos(tau / 2)) ^ 3;
b1 = 1 + re * a1;
b2 = (a2 - a1) / b1;
tv = b1 * (1 + re * b2 * a3);
tl = tc * tv;

if (i1 == 1)
   tg = d1 * gi + d2 * exp(gi);
else
   tg = d3 * gi + d4 * (1 - exp(d5 * gi));
end

g3 = 0.5 * (1 + sin(e10 * dy + e11));

g3 = g3 ^ e12;

tau1 = dy + e8 * (g3 - e9);

g1 = e2 + e3 * (sin(e4 * tau1 + e5));

g2 = sin(e6 * tau1 + e7);

ts = e1 + f10b * g1 * g2;

te = tl + tg + ts;

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

function den = jslv (alt, xlat, day)

den = 0;

if (alt > 170)
   return;
end

z = alt - 90;
x = -0.0013 * z * z;
y = 0.0172 * day + 1.72;
p = sin(y);
sp = (sin(xlat)) ^ 2;
s = 0.014 * z * exp(x);
d = s * p * sp;

if (xlat < 0)
   d = -d;
end

den = d;

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

function den = jslvh (den, denhe, xlat, sda)

ezero = 10 ^ denhe;

a = abs(0.65 * (sda / 0.40909079));

b = 0.5 * xlat;

if (sda < 0)
   b = -b;
end

x = 0.7854 - b;

y = sin(x) ^ 3;

dhe = a * (y - 0.35356);

denhe = denhe + dhe;

d1 = 10 ^ denhe;

del = d1 - ezero;

rho = 10 ^ den;

drho = 6.646e-24 * del;

rho = rho + drho;

den = log10(rho);

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

function [fdhel, fdlg] = jfair (dhel1, dhel2, dlg1, dlg2, ih)

% fair density

global cz

bfh = 440;

i = fix((ih - bfh) / 10) + 1;

czi = cz(i);
szi = 1 - czi;

fdlg = (dlg1 * czi) + (dlg2 * szi);

fdhel = (dhel1 * czi) + (dhel2 * szi);

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

function r = jgauss (z1, z2, nmin, tx, t1, t3, t4, a2)

% gaussian quadrature

global ng xgauss cgauss altmin

r = 0;

for k = nmin:1:8
    ngauss = ng(k);
    a = altmin(k);
    d = min(z2, altmin(k + 1));
    rr = 0;
    del = 0.5 * (d - a);
    j = ngauss - 2;

    for i = 1:1:ngauss
        z = del * (xgauss(i, j) + 1) + a;

        rr = rr + cgauss(i, j) * jmweight(z) * jgrav(z) ...
             / jtemp(z, tx, t1, t3, t4, a2);
    end

    rr = del * rr;
    r = r + rr;

    if (d == z2)
       break;
    end
end

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

function agrav = jgrav(altitude)

% acceleration of gravity

agrav = 9.80665 / ((1 + altitude / 6356.766) ^ 2);

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

function mweight = jmweight (a)

% molecular weight function

global bdata

if (a > 105)
   mweight = 1;
else
   u = a - 100;

   wttmp = bdata(1);

   for i = 2:1:7
       wttmp = wttmp + bdata(i) * u ^ (i - 1);
   end

   mweight = wttmp;
end

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

function temp = jtemp (alt, tx, t1, t3, t4, a2)

% temperature function

bb = 0.0000045;

u = alt - 125.0;

if (u > 0)
   temp = tx + a2 * atan(t1 * u * (1 + bb * (u ^ 2.5)) / a2);
else
   temp = tx + t1 * u + t3 * (u ^ 3) + t4 * (u ^ 4);
end


text:
function a = gk324 ( m, n, x )

%*****************************************************************************80
%
%% GK324 returns the GK324 matrix.
%
%  Discussion:
%
%    This is Gregory and Karney example matrix 3.24.
%
%  Example:
%
%    M = N = 5
%
%    X = ( 11, 12, 13, 14 )
%
%     1  1  1  1  1
%    11  1  1  1  1
%    11 12  1  1  1
%    11 12 13  1  1
%    11 12 13 14  1
%
%  Properties:
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 51,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, integer M, N, the order of the matrix.
%
%    Input,
%    * real X(N-1), the first N-1 entries of the
%      last row, if M <= N,
%    or
%    * real X(N), the N entries of the last row,
%      if N < M.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n
      if ( i <= j )
        a(i,j) = 1.0;
      else
        a(i,j) = x(j);
      end
    end
  end

  return
end


text:
function [f_I] = ToolInstFreq(X,iHop,fs)
    % get phase
    phi         = angle(X);

    % phase offset
    omega       = pi*iHop/size(X,2)*(0:size(X,2)-1);

    % unwrapped difference
    deltaphi    = omega + princarg_I(phi(2,:)-phi(1,:)-omega);

    % instantaneous frequency
    f_I         = deltaphi/iHop/(2*pi)*fs;
end

function phase = princarg_I(phi)
    phase = mod(phi + pi,-2*pi) + pi;
end

author: alexanderlerch

repo: ACA-Code


text:
%% PATH TRACING, MEASUREMENTS, FRAGMENTATION

%% Path Tracing and Computation
%
%       An application of path tracing is intended for
% detecting  an external borders  of 1-connected black
% pixels' sets (spots), computing of  spots' and whole
% image quantitative  characteristics,  extracting and
% removing spots or group of spots.
%
%       An output includes:
% - two matrices with coordinates of each spot border
%   pixels;
% - matrice with data for each spot in a proper string:
%   - coordinates of points of optimal location (centers
%     of mass) for a spot and for a spot border,
%   - perimeter (as length of line connecting  centers of
%     border pixels and as a number of border pixels'
%     external edges),
%   - number of  pixels in a spot and its border,
%   - coordinates of the most remote pixels,
%   - a distance between them,
%   - coordinates of vertices of minimal  rectangle
%     with a spot inside it.
% - The last string of the above mentioned matrice
%   contains data for all spots in the image:
%   - mass center coordinates of all curves,
%   - mass center coordinates of all spots,
%   - a number of pixels of all borders,
%   - a number of pixels of all spots,
%   - coordinates of the most remote black points in the image,
%   - a distance between them,
%   - perimeter (a number of border pixels' external edges),
%   - perimeter (the total length of lines connecting
%     border pixels' centers),
% - a series of image files with spots.
%       Output data might be calculated from both the
% image and a numeric matrix.
% The author - Eduard Polityko, PHD.
% E-mail     - Edpolit@gmail.com
% Edition 07-Feb-2007

alfa='examp16.bmp';
[A B XC]=digisn(alfa);

%% Displaying results

figure; C=[0 0 0;1 1 1];
image(imread(alfa));    colormap(C);
title('Spots'); axis xy;    grid on
format short g;
rp=repmat('_',1,26);
disp('     Coordinates');   disp(rp)
disp('Centers of mass of contours')
disp(XC(1:end-1,1:2));
disp('Centers of mass of spots');
disp([XC(1:end-1,3:4)]);
disp(['Rectangles to crop ';'(ends of diagonal) '])
disp(XC(1:end-1,10:13))
disp('The most remote points')
disp([A(XC(1:end-1,7)),B(XC(1:end-1,7)),A(XC(1:end-1,8)),...
  B(XC(1:end-1,8))])
disp(rp);   disp('Maximal distances between')
disp('points:');    disp(XC(1:end-1,9));
disp('A number of black pixels:')
disp(['       Border      ','Spot'])
disp([XC(1:end-1,5) XC(1:end-1,6)])
disp('         Perimeter:')
disp('  a number of external edges')
disp('      of border pixels ');
disp(XC(1:end-1,14))
disp('  lengths of lines connecting')
disp('  centers of border pixels')
disp(XC(1:end-1,15));   disp([rp;rp])
disp('SUMMARY DATA');
disp('     Coordinates');   disp(rp)
disp('Center of mass of contours'); disp(XC(end,1:2));
disp('Center of mass of spots   '); disp(XC(end,3:4));
disp('The most remote points'); disp(XC(end,7:10));
disp(rp);   disp('Maximal distance between')
disp(['points - ' num2str(XC(end,11))]);    disp(' ')
disp('A number of black pixels:')
disp(['       Border      ','Spot'])
disp([XC(end,5) XC(end,6)]);
disp('        Perimeter:')
disp('  a number of external edges')
disp(['  of border pixels - ' num2str(XC(end,14))])
disp('  length of lines connecting')
disp('  centers of border pixels -');   disp(XC(end,15))
figure; zx1=plot(XC(end,1),XC(end,2),'+k',...
  XC(end,3),XC(end,4),'*k',[XC(end,7),XC(end,9)],...
  [XC(end,8),XC(end,10)],'ok');grid on
legend(['center of  '; 'curves mass'],...
  ['center of '; 'spots mass'],...,
  ['the most     ';'remote points'],'location', 'best');
title('Measurement');   set(zx1,'markersize',6)
hold on;    zx=plot(A,B,'.'); set(zx,'markersize',1)

%% Series images with spots
%Extract all spots and place them as image files
%"Newpic1.bmp", "Newpic2.bmp"...

figure; newpic='newpic';
digisn(alfa,newpic);
sXC=size(XC,1)-1;
for i=1:sXC
  sX=ceil(sXC/3);
  subplot(sX,3,i)
image(imread([newpic num2str(i),'.bmp']));colormap(C)
end

%% Spots removing
%Remove the first three spots from the image and place
%the rest as image file 'RestAfterRemoving3.bmp',
%where 3 is a number of removed spots)

figure; k=3;    digisn(alfa,-k);
image(imread(['RestAfterRemoving' num2str(k) '.bmp']));
title(['The rest after removing first ' num2str(k) ' spots']);
colormap(C);    grid on;

%% Keep spots
%Extract spots from #3 to #5 and place them as
%image file "Keep_35.bmp"
figure; k=[3;5];    digisn(alfa,k);
image(imread(['Keep' num2str(k(1)) '_' num2str(k(2)) '.bmp']));
title(['Spots from #' num2str(k(1)) ' to #' num2str(k(2))]);
colormap(C);grid on


text:
function Ypred = continuous_prediction_2class(Y,Ypred)
% Which one is the most likely class, out of a continuous estimation?
[N,q] = size(Y); % N is time points
L = length(unique(Y(:)));
if q == 1 && L == 2 % class is -1 or +1
    Ypred(Ypred<=0) = -1;
    Ypred(Ypred>0) = +1;
elseif q > 1 && L == 2 % dummy variable representation for classification
    [~,m] = max(Ypred,[],2);
    Ypred = zeros(size(Ypred));
    ind = sub2ind(size(Ypred),1:N, m');
    Ypred(ind) = 1;
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
% Version 1.000
%
% Code provided by Ruslan Salakhutdinov and Geoff Hinton
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

function [f, df] = CG_CLASSIFY(VV,Dim,XX,target);

l1 = Dim(1);
l2 = Dim(2);
l3= Dim(3);
l4= Dim(4);
l5= Dim(5);
N = size(XX,1);

% Do decomversion.
 w1 = reshape(VV(1:(l1+1)*l2),l1+1,l2);
 xxx = (l1+1)*l2;
 w2 = reshape(VV(xxx+1:xxx+(l2+1)*l3),l2+1,l3);
 xxx = xxx+(l2+1)*l3;
 w3 = reshape(VV(xxx+1:xxx+(l3+1)*l4),l3+1,l4);
 xxx = xxx+(l3+1)*l4;
 w_class = reshape(VV(xxx+1:xxx+(l4+1)*l5),l4+1,l5);

  XX = [XX ones(N,1)];
  w1probs = 1./(1 + exp(-XX*w1)); w1probs = [w1probs  ones(N,1)];
  w2probs = 1./(1 + exp(-w1probs*w2)); w2probs = [w2probs ones(N,1)];
  w3probs = 1./(1 + exp(-w2probs*w3)); w3probs = [w3probs  ones(N,1)];

  targetout = exp(w3probs*w_class);
  targetout = targetout./repmat(sum(targetout,2),1,10);
  f = -sum(sum( target(:,1:end).*log(targetout))) ;

IO = (targetout-target(:,1:end));
Ix_class=IO;
dw_class =  w3probs'*Ix_class;

Ix3 = (Ix_class*w_class').*w3probs.*(1-w3probs);
Ix3 = Ix3(:,1:end-1);
dw3 =  w2probs'*Ix3;

Ix2 = (Ix3*w3').*w2probs.*(1-w2probs);
Ix2 = Ix2(:,1:end-1);
dw2 =  w1probs'*Ix2;

Ix1 = (Ix2*w2').*w1probs.*(1-w1probs);
Ix1 = Ix1(:,1:end-1);
dw1 =  XX'*Ix1;

df = [dw1(:)' dw2(:)' dw3(:)' dw_class(:)']';

author: qiuwch

repo: DeepLearning

save_path: github-repos/MATLAB/qiuwch-DeepLearning


text:
clc;
load('detector1.mat');

I1=imread('C:\Shivanshu laptop data\Desktop data\Ortho_HOG\video\frame\frame205.jpg');

I=rgb2gray(I1);
J=edge(I,'Canny',[0.15 0.3]);
s=size(J);
u=s(1,1);
v=s(1,2);
y=[2*u/5 u u 2*u/5];%h
x=[v/2 v/4 3*v/4 v/2];%w
mask=poly2mask(x,y,u,v);
f=J.*mask;
y1=[0 u u 0];%h
x1=[0 0 v/2 v/2];%w
mask1=poly2mask(x1,y1,u,v);
f1=f.*mask1;
[H,T,R]=hough(f1);
y2=[0 u u 0];%h
x2=[v/2 v/2 v v];%w
mask2=poly2mask(x2,y2,u,v);
f2=f.*mask2;
[H,T,R]=hough(f1,'RhoResolution',2,'ThetaResolution',5);
P=houghpeaks(H,5,'threshold',ceil(0.5*max(H(:))));
lines1=houghlines(f1,T,R,P,'FillGap',15,'MinLength',7);

[H2,T2,R2]=hough(f2,'RhoResolution',2,'ThetaResolution',5);
P2=houghpeaks(H2,5,'threshold',ceil(0.5*max(H2(:))));
lines2=houghlines(f2,T2,R2,P2,'FillGap',15,'MinLength',7);

figure(6)
imshow(I);
title('0.6');

max_len1=0;
max_len2=0;
%AA1=zeros(1,length(lines1));
%AA2=zeros(1,length(lines2));
for k=1:length(lines1)
len1=norm(lines1(k).point1-lines1(k).point2);
xy1=[lines1(k).point1;lines1(k).point2];
if(len1>max_len1)
    max_len1=len1;
    xy_long=xy1;
end
end

for k=1:length(lines2)
len2=norm(lines2(k).point1-lines2(k).point2);
xy2=[lines2(k).point1;lines2(k).point2];
if(len2>max_len2)
    max_len2=len2;
    xy_2=xy2;
end
end

%for k=1:length(lines1)
 %   disp(k);
  %  xy=[lines1(k).point1;lines1(k).point2];
   % len=norm(lines(k).point1-lines(k).point2);
    %disp(xy(1,2));
    %if(xy(1,2)<=(v/2))
     %       AA1(1,k)=len;
    %else
     %   AA2(1,k)=len;
    %end
%end
%[max_len1,idx]=max(AA1);
%for k=1:length(lines)
%if(AA1(1,k)==max_len1)
 %       xy_long=[lines(k).point1;lines(k).point2];
  %      end
%end
%[max_len2,idx]=max(AA2);
%for k=1:length(lines)
%if(AA2(1,k)==max_len2)
 %       xy2=[lines(k).point1;lines(k).point2];
  %      end
%end

    %[ee,ff]=max(AA);
    %xy2=[lines(ff).point1;lines(ff).point2];
     %plot(xy_long(:,1),xy_long(:,2)+[20; 150],'LineWidth',2,'Color','green');%
     %plot(xy_2(:,1),xy_2(:,2)+[-10; -20],'LineWidth',2,'Color','green');

h1=imread('C:\Shivanshu laptop data\Desktop data\Ortho_HOG\video\frame\frame739.jpg');
I1=imread('C:\Shivanshu laptop data\Desktop data\Ortho_HOG\video\frame\frame194.jpg');
I=rgb2hsv(I1);
h=rgb2hsv(h1);
%figure(1)
%imshow(I1);
II=I(:,:,3);
HH=h(:,:,3);
for i=1:360
    for j=1:450
        if((II(i,j)<=1) && (II(i,j)>0.85))
            II(i,j)=II(i,j);
        else
           II(i,j)=0;
        end
    end
end

for i=1:360
    for j=1:450
        if((HH(i,j)<=1) && (HH(i,j)>0.12))
            HH(i,j)=HH(i,j);
        else
           HH(i,j)=0;
        end
    end
end
%figure(2)
%imshow(II);

s=size(II);
u=s(1,1);
v=s(1,2);
y=[2*u/5 u u 1*u/3 2*u/5];%h
x=[v/2 v/3 v v v/2];%w
mask=poly2mask(x,y,u,v);
f=II.*mask;
%figure(2);
%imshow(mask);
hh=HH.*mask;
hh=1-hh;

for i=1:360
  for  j=1:450
if(hh(i,j)<1)
    hh(i,j)=0;
end
end
end
f=f.*hh;

[bbox, score, label] = detect(detector1,f)
%Display detection results.
%detectedImg = insertObjectAnnotation(I1,'Rectangle',bbox);
%detectedImg = insertShape(I1,'Rectangle',bbox);
%figure(3)
%imshow(detectedImg)

% Display detection results
label_str = cell(1,1);
%conf_val = [score];
conf_lab = [label];
%for ii=1:3
    label_str =[sprintf('%s',conf_lab)];
%end
position = [bbox];
outputImage = insertObjectAnnotation(I1,'rectangle',position,label_str,...
    'TextBoxOpacity',0.9,'FontSize',10);
figure(4)
imshow(outputImage)
hold on
plot(xy_long(:,1),xy_long(:,2)+[20; 150],'LineWidth',2,'Color','green');%
     plot(xy_2(:,1),xy_2(:,2)+[-10; -20],'LineWidth',2,'Color','green');

author: Aarchishya


text:
function percentile = prctile(this, percentile, varargin)
% Computes given percentile in data (e.g. 50 for median)
%
%   Y = MrImage()
%   percentile = ...
%       Y.prctile(percentile, 'ParameterName1', 'ParameterValue1', ...)
%
% This is a method of class MrImage.
%
% IN
%   percentile  percent for percentile computation; default: 50 (median)
%   varargin    parameterName/Value pairs for selection of volumes/slices
%
% OUT
%
% EXAMPLE
%   Y.prctile(50, 'z', 1, 't', 3:100, ...,
%           'x', 55:75)
%
%   See also MrImage MrImage.hist

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

if nargin < 2
    percentile = 50;
end

if nargin < 3
    imgSelect = this;
else
    imgSelect = this.select(varargin{:});
end

percentile = prctile(imgSelect.data(:), percentile);

author: translationalneuromodeling

repo: tapas


text:
function [HDR]=bni2hdr(arg1,arg3,arg4,arg5,arg6)
% BNI2HDR converts BNI header information into BioSig Header information
%	HDR = BNI2HDR(HDR)
%
% INPUT:
%   HDR.H1 contains ascii header
%
% OUTPUT:
%   HDR.Label
%   HDR.T0
%       ...
%
% see also: SOPEN

%	$Id$
%	Copyright (c) 2007,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 3
% of the  License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

if isstruct(arg1),
	HDR=arg1;
elseif exist(arg1,'file');
	HDR.FileName = arg1;
	[HDR.FILE.Path,HDR.FILE.Name,HDR.FILE.Ext]=fileparts(arg1);
end;
if ~isfield(HDR,'H1')
	fid = fopen(fullfile(HDR.FILE.Path,HDR.FILE.Name,'.bni'),'rt');
	HDR.H1 = char(fread(fid,[1,inf],'uint8'));
	fclose(fid);
end;

s = HDR.H1;
if ~strcmp(s(1:28),'FileFormat = BNI-1-BALTIMORE')
	fprintf(HDR.FILE.stderr,'WARNING BNI2HDR: Header information is not Nicolet BNI format.\n');
end;

while ~isempty(s),
	[t,s]=strtok(s,[10,13]);
	[t1,t2]=strtok(t,' =');
	[t2,t3]=strtok(t2,' =');
	if strcmp(t1,'PatientId')
		HDR.Patient.Id = t2;
	elseif strcmpi(t1,'Sex')
		HDR.Patient.Sex = t2; % strncmpi(t2,'m',1)+2*strncmpi(t2,'f',1);
	elseif strncmpi(t1,'medication',10)
		HDR.Patient.Medication = t2;
	elseif strncmpi(t1,'diagnosis',10)
		HDR.Patient.Diagnosis = t2;
	elseif strcmpi(t1,'MontageRaw')
		[tmp1,tmp2,HDR.Label] = str2double(t2,',');
	elseif strcmpi(t1,'Age')
		HDR.Patient.Age = str2double(t2);
	elseif strcmp(t1,'Date')
		if any(t2=='/')
			t2(t2=='/')=' ';
			HDR.T0([2,3,1])=str2double(t2);
		end;
	elseif strcmp(t1,'Time')
		t2(t2==':') = ' ';
		HDR.T0(4:6) = str2double(t2);
	elseif strcmp(t1,'Rate')
		HDR.SampleRate = str2double(t2);
	elseif strcmp(t1,'NchanFile')
		HDR.NS = str2double(t2);
	elseif strcmp(t1,'UvPerBit')
		HDR.Cal = str2double(t2);
	elseif strcmp(t1,'[Events]')
		% not supported yet
		s = [];
	end;
end;


text:
% This example shows a simple comparison of two different algorithm for tensor completion:
%
%   -- ALS completion
%   -- Riemannian tensor completion (RTTC)
%
% in a very similar comparison as Figure 5.2. in
%
%   Michael Steinlechner, Riemannian optimization for high-dimensional tensor completion,
%   Technical report, March 2015, revised December 2015.
%   To appear in SIAM J. Sci. Comput.
%
% See this report for more details about the algorithms and the setup.
% The different to the therein described setup is only a reduced problem size (d, n, r) so
% that it takes less time to compute the results.

%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

rng(13);
d = 10;

ranks = [4, 6, 8];

cost = cell(1,length(ranks));
test = cell(1,length(ranks));
stats = cell(1,length(ranks));

for j = 1:length(ranks)
    r = ranks(j);
    rr = [1, r*ones(1,d-1), 1];

    nn = 20;
    n = nn*ones(1,d);

    opts = struct('maxiter', 50, 'tol', 0, 'reltol',0, 'gradtol',0);
    opts_tt = struct('maxiter', 60, 'tol', 0, 'reltol',0, 'gradtol',0);

    dof = d*nn*r^2;
    sizeOmega = 10*dof;
    sizeGamma = sizeOmega;

    Omega = makeOmegaSet_mod(n, sizeOmega);
    Gamma = makeOmegaSet_mod(n, sizeGamma);

    A = TTeMPS_rand( rr, n );
    A = 1/norm(A) * A;

    A_Omega = A(Omega);
    A_Gamma = A(Gamma);

    X0 = TTeMPS_rand( rr, n );
    X0 = 1/norm(X0) * X0;
    X0 = orthogonalize( X0, X0.order );

    [X,cost_als{j},test_als{j},stats_als{j}] = completion_als( A_Omega, Omega, A_Gamma, Gamma, X0, opts );
    [X,cost_tt{j},test_tt{j},stats_tt{j}] = completion_orth( A_Omega, Omega, A_Gamma, Gamma, X0, opts_tt );
end

l = lines(7);
midred = l(end,:);
darkred = brighten(l(end,:),-0.7);
lightred = brighten(midred,0.7);

midblue = l(1,:)
darkblue = brighten(midblue,-0.7);
lightblue = brighten(midblue,0.7);

subplot(1,2,1)
semilogy( test_als{1}(1:end),'color',darkred,'linewidth',2)
hold on
semilogy( test_als{2}(1:end),'color',midred,'linewidth',2)
semilogy( test_als{3}(1:end),'color',lightred,'linewidth',2)
semilogy( test_tt{1},'--','color',darkblue,'linewidth',2)
semilogy( test_tt{2},'--','color',midblue,'linewidth',2)
semilogy( test_tt{3},'--','color',lightblue,'linewidth',2)

xlabel('Iterations')
ylabel('Error on test set')
legend({'ALS, rank 4','ALS, rank 6', 'ALS, rank 8','RTTC, rank 4', 'RTTC, rank 6', 'RTTC, rank 8'})

subplot(1,2,2)
loglog( stats_als{1}.time(1:end), test_als{1}(1:end),'color',darkred,'linewidth',2)
hold on
loglog( stats_als{2}.time(1:end), test_als{2}(1:end),'color',midred,'linewidth',2)
loglog( stats_als{3}.time(1:end), test_als{3}(1:end),'color',lightred,'linewidth',2)
loglog( stats_tt{1}.time, test_tt{1},'--','color',darkblue,'linewidth',2)
loglog( stats_tt{2}.time, test_tt{2},'--','color',midblue,'linewidth',2)
loglog( stats_tt{3}.time, test_tt{3},'--','color',lightblue,'linewidth',2)
xlim([1e-1,1e3])

xlabel('Time [s]')
ylabel('Error on test set')

author: NicolasBoumal

repo: manopt