text:
% predicting numerical values using CCA
% large values and small values are corrected

clear all;
format long
disp('===== Linear CCA for Regression ====');
disp('Reading featur vector');

output_mean_mse_array = zeros(7,9);
colinoutput = 0;
for helel_data = [100,200,400,800,1200,1500, 2000,5000,10000]
    colinoutput = colinoutput+1;
    Indices = crossvalind('Kfold',helel_data , 10);%crossvalind('Kfold', 26548, 10);
    figure;

    subplot(7,2,1);

    for feat = 1:7
        MSEarray =[];
        elapsedarray =[];
        for crossvalidateIter = 1:10
            (fprintf('%d, %.5f',crossvalidateIter, mean(MSEarray)));
            possiblefeaturizations =  {'tfidfbucket','all','logmultinomial', 'logmultinomial2', 'logmultinomial3','bernouli', 'tfidf','multinomial'};
            %featurization = 'bernouli'%'tfidf'%'tfidf'%'multinomial'%'tfidf' %'multinomial'; % 'bernouli', 'tfidf'
            featurization  = possiblefeaturizations{feat};
            featurs = csvread('data\forWeka_featuresonly.csv');
            featurs = featurs(:,2:size(featurs,2));

            num_data = helel_data%size(featurs,1); %5000;

            %disp(sprintf('Number of datapoints %d',num_data))
            if strcmp(featurization,'multinomial')
                %just pass
            elseif strcmp(featurization,'bernouli')
                featurs = bernoulli(featurs);
            elseif strcmp(featurization,'tfidf')
                featurs = tfidf(featurs);
            elseif strcmp(featurization,'logmultinomial')
                featurs = log(featurs+1)./log(2);
            elseif strcmp(featurization,'logmultinomial2')
                featurs = round(log(featurs+1)./log(2));
            elseif strcmp(featurization,'logmultinomial3')
                featurs = (log(featurs+1));
            elseif strcmp(featurization,'all')
                featurs = [bernoulli(featurs), tfidf(featurs), (log(featurs+1))] ;
            elseif strcmp(featurization,'tfidfbucket')
                featurs = [log(1000*tfidf(featurs)+1)] ;

            end

            size_training = floor(.9*num_data);

            trainingset = featurs(Indices~=crossvalidateIter,:);
            testset = featurs(Indices==crossvalidateIter,:);

            %disp('Splitting up data into training/test sets');
            [num,txt,raw] = xlsread('data\final104.xls');

            % reading the description of each shoe
            descriptions = raw(2:size(raw,1),2);
            style_ratings = num(1:size(num,1),1);
            comfort_ratings = num(1:size(num,1),4);
            overal_ratings = num(1:size(num,1),5);

            % only take m data points
            m=num_data;
            descriptions = descriptions(1:m);
            style_ratings = style_ratings(1:m);
            comfort_ratings = comfort_ratings(1:m);
            overal_ratings = overal_ratings(1:m);

            responsevals = [style_ratings, comfort_ratings, overal_ratings];

            responsevals_training = responsevals(Indices~=crossvalidateIter,:);
            responsevals_test = responsevals(Indices==crossvalidateIter,:);

            %disp('Adjusted Linear CCA ');
            % http://www.mathworks.com/help/toolbox/stats/classregtree.html

            tic;

            %predictions = [];
            %actual = [];

            [Wx, Wy, r, U, V]  = canoncorr(trainingset,responsevals_training);
            % in [A,B,r,U,V],  U and V are cannonical scores
            % U = (X-repmat(mean(X),N,1))*A
            % V = (Y-repmat(mean(Y),N,1))*B

            %recheck this part
            N =  size(testset,1);
            predictions = ((testset-repmat(mean(testset),N,1))*Wx*pinv(Wy))+repmat(mean(responsevals_test),N,1);
            actual = responsevals_test;
            predictions( predictions>5)=5;
            predictions( predictions<1)=1;
            MSE = mean(sum(((predictions-actual).^2)'));
            elapsed = toc;
            MSEarray = [MSEarray MSE];
            elapsedarray = [elapsedarray elapsed];

        end
        featurization  = possiblefeaturizations{feat}
        fprintf('avergae MSE for adjusted CCA is %0.10f\n', mean(MSEarray));
        subplot(7,2,feat*2-1)
        plot(MSEarray,'-o');
        title(strcat('MSE adjusted CCA (',featurization,') ', sprintf(' avergae MSE = %0.10f\n', mean(MSEarray))));
        xlabel('10 fold cross-validation (iteration no)')
        ylabel('MSE')
        subplot(7,2,feat*2)
        plot(elapsedarray,'-o');
        title(strcat('Elapsed time for adjusted CCA (',featurization, ') ' , sprintf(' avergae elapsed time = %0.10f\n', mean(elapsedarray))));
        xlabel('10 fold cross-validation (iteration no)')
        ylabel('Elapsed Time')
        mtit(sprintf('%d datapoints',helel_data))
        drawnow;
        output_mean_mse_array(feat,colinoutput) = mean(MSEarray);
    end

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [cum] = lgpdens_cum(bb,x1,x2)
% [CUM] = LGPDENS_CUM(BB)
%
% Description
%   Given Bayesian Bootstrap estimated density, integrates it from point x1
%   to point x2
%
% Copyright (c) 2012 Ernesto Ulloa
% Copyright (c) 2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.addRequired('bb',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('x1',@(x) ~isempty(x) && isreal(x))
ip.addRequired('x2', @(x) ~isempty(x) && isreal(x))
ip.parse(bb,x1,x2)

 [p,pq,xt]=lgpdens(bb);
 I1=min(find(xt>x1));
 I2=max(find(xt<x2));
 sd=xt(2)-xt(1);
 cum=sd*trapz(p(I1:I2));

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [Xpareto,Ypareto,X,Y,XParetoGP,YParetoGP,YParetoGPstd,hypf] = TSEMO_V4(f,X,Y,lb,ub,Opt)
%  Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2020-23-05.
%  Version 4 by Eric Bradford (eric.bradford@ntnu.no), Artur M.
%  Schweidtmann (Artur.Schweidtmann@avt.rwth-aachen.de) and Alexei Lapkin
%  (aal35@cam.ac.uk)

% ALGORITHM  Executes the optimization algorithm.

%   OUTPUTS
%       Xpareto      Pareto front approximation from X inputs       [np,D]
%       Ypareto      Pareto front approximation from Y outputs      [np,O]
%       X            Data input vector                              [n+meval,D]
%       Y            Data output vector                             [n+meval,O]
%       XParetoGP    Pareto front approximation of GP model inputs  [npGP,D]
%       YParetoGP    Pareto front approximation of GP model outputs [npGP,D]
%       YParetoGPstd Pareto front standard deviation of GP model    [npGP,D]
%       hypf         Final hyperparameter values []                 [D+2,O]

%   INPUTS
%       f       objective function
%       X       input values of given data set  [n,D]
%       Y       output values of given data set [n,O]
%       lb      lower bounds of input variables [1,D]
%       ub      upper bounds of input variables [1,D]
%       Opt     options struct of algorithm

%   with
%       n       number of given data points
%       meval   maximum number of function evaluations
%       np      number of pareto points of the final dataset
%       npGP    number of pareto points of GP model dataset
%       D       dimension of input space
%       O       dimension of output space

%% Initialize option structure
it = 1;
Opt = set_option_structure(Opt,X,Y);

%% Write initial text file for TSEMO_log
create_log_file(X,Y,Opt,f,lb,ub)

for i = 1:ceil(Opt.maxeval/Opt.NoOfBachSequential)
    tic;
    %% Scale Variables
    [Xnew,Ynew] = ScaleVariables(X,Y,lb,ub,Opt) ;

    %% Training of GP
    for j = 1:Opt.Gen.NoOfGPs
        Opt.GP(j).hyp = TrainingOfGP(Xnew,Ynew(:,j),Opt.GP(j));
    end

    %% Draw samples of GPs
    for j = 1:Opt.Gen.NoOfGPs
        Opt.Sample(j).f = posterior_sample(Xnew,Ynew(:,j),Opt.GP(j));
    end

    %% Determine pareto front of function samples
    [Sample_pareto,Sample_xpareto,Sample_nadir] = Find_sample_pareto(Opt,i);
    Opt.warmstart_pareto = Sample_xpareto;

    %% Sampling point at maximum hypervolume improvement
    [index,hv_imp] = hypervolume_improvement_index(Ynew,Sample_nadir,Sample_pareto,Opt);

    xNew = Sample_xpareto(index,:);
    Xnew = [Xnew;xNew];

    for j = 1:Opt.Gen.NoOfInputDim
        xnewtrue(:,j) = xNew(:,j)*(ub(j)-lb(j)) + lb(j);
    end

    for l = 1 : size(xnewtrue,1)
        ytrue(l,:) = f(xnewtrue(l,:));
    end

    X = [X;xnewtrue];
    Y = [Y;ytrue];
    Ynew = zeros(size(Y,1),Opt.Gen.NoOfGPs);
    for j = 1:Opt.Gen.NoOfGPs
        Ynew(:,j) = (Y(:,j) - mean(Y(:,j)))/std(Y(:,j));
    end

    front = paretofront(Y);
    Xpareto = X(front,:);
    Ypareto = Y(front,:);

    %% Update log each iterations
    update_log_file(it,hv_imp,toc,xnewtrue,ytrue,Opt,Y,ub,lb)

    %% Display
    if it == 1
        fprintf('%10s %10s %10s \n','Iteration', 'HypImp', 'Time(s)');
    end
    fprintf('%10d %10.4g %10.3g \n', it, hv_imp, toc);

    it = it+1;

    %% Determine final hyperparameter values for analysis
    if i == ceil(Opt.maxeval/Opt.NoOfBachSequential)
        for j = 1:Opt.Gen.NoOfGPs
            Opt.GP(j).hyp = TrainingOfGP(Xnew,Ynew(:,j),Opt.GP(j));
        end

        hypf = zeros(Opt.Gen.NoOfInputDim+2,Opt.Gen.NoOfGPs);
        for j = 1:Opt.Gen.NoOfGPs
            covhyp = exp(Opt.GP(j).hyp.cov);
            hypf(:,j) = [covhyp(1:Opt.Gen.NoOfInputDim).*(ub-lb)';covhyp(end)*std(Y(:,j));exp(Opt.GP(j).hyp.lik)*std(Y(:,j))];
        end

        %% Obtain Pareto front from spectral Gaussian process model
        for j = 1:Opt.Gen.NoOfGPs
            [Opt.Mean(j).f,Opt.Mean(j).varf] = mean_sample(Xnew,Ynew(:,j),Opt.GP(j));
        end

        [Mean_pareto,Mean_xpareto] = Find_mean_pareto(Opt);

        for j = 1:Opt.Gen.NoOfInputDim
            XParetoGP(:,j) = Mean_xpareto(:,j)*(ub(j)-lb(j)) + lb(j);
        end

        for j = 1:Opt.Gen.NoOfGPs
            YParetoGP(:,j) = Mean_pareto(:,j)*std(Y(:,j)) + mean(Y(:,j));
        end

        for j = 1:Opt.Gen.NoOfGPs
            for k = 1:Opt.pop
                YParetoGPstd(k,j) = sqrt(Opt.Mean(j).varf(Mean_xpareto(k,:)))*std(Y(:,j));
            end
        end

        %% Update log with final results
        final_log_update(Xpareto,Ypareto,X,Y,XParetoGP,YParetoGP,hypf,Opt)
    end
end
return

function create_log_file(X,Y,Opt,f,lb,ub)
try
    function_name = func2str(f);
    string1 = '';
    string2 = {};
    string3 = '';
    string4 = '';
    string5 = {};
    string6 = '';
    string7 = '';
    for i = 1:size(X,2)
        string1 = strcat(string1,'%8.4f ');
        string2 = {string2{:},strcat('x',num2str(i))};
        string3 = strcat(string3,'%+8s ');
    end
    for i = 1:size(Y,2)
        string4 = strcat(string4,'%8.4f ');
        string5 = {string5{:},strcat('f',num2str(i))};
        string6 = strcat(string6,'%+8s ');
        string7 = strcat(string7,'%8d ');
    end

    TSEMO_log = fopen( 'TSEMO_log.txt', 'w');
    fprintf(TSEMO_log,'\n %s %s \n', 'TSEMO log file created on',date);
    fprintf(TSEMO_log,'\n %s','This file shows the initial specifications of TSEMO and logs the output.');
    fprintf(TSEMO_log,'\n %s', '¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯');

    fprintf(TSEMO_log,'\n %s \n', 'License information');
    fprintf(TSEMO_log,'\n %s \n', 'BSD 2-Clause License');
    fprintf(TSEMO_log,'\n %s', 'Copyright (c) 2017, Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin');
    fprintf(TSEMO_log,'\n %s \n', 'All rights reserved.');
    fprintf(TSEMO_log,'\n %s', 'Redistribution and use in source and binary forms, with or without');
    fprintf(TSEMO_log,'\n %s \n', 'modification, are permitted provided that the following conditions are met:');
    fprintf(TSEMO_log,'\n %s   ', '*Redistributions of source code must retain the above copyright notice, this');
    fprintf(TSEMO_log,'\n %s \n', ' list of conditions and the following disclaimer.');
    fprintf(TSEMO_log,'\n %s   ', '*Redistributions in binary form must reproduce the above copyright notice,');
    fprintf(TSEMO_log,'\n %s   ', ' this list of conditions and the following disclaimer in the documentation');
    fprintf(TSEMO_log,'\n %s \n', ' and/or other materials provided with the distribution.');
    fprintf(TSEMO_log,'\n %s', '¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯');

    fprintf(TSEMO_log,'\n %s \n', 'Problem specifications');
    fprintf(TSEMO_log,'\n %s %s \n', 'Function used:  ',function_name);
    fprintf(TSEMO_log,'\n %s %d', 'Number of inputs:  ',size(X,2));
    fprintf(TSEMO_log,'\n %s %d \n', 'Number of outputs: ',size(Y,2));
    fprintf(TSEMO_log,'\n %s', 'Lower bounds of decision variables:');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,string1,lb);
    fprintf(TSEMO_log,'\n \n %s', 'Upper bounds of decision variables:');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,string1,ub);
    fprintf(TSEMO_log,'\n %s', '¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯');

    fprintf(TSEMO_log,'\n %s \n', 'Algorithm options');
    fprintf(TSEMO_log,'\n %s %d', 'Maximum number of function evaluations: ',Opt.maxeval);
    fprintf(TSEMO_log,'\n %s %d', 'Sample batch size:                      ',Opt.NoOfBachSequential);
    fprintf(TSEMO_log,'\n %s %d \n', 'Number of algorithm iterations:         ',ceil(Opt.maxeval/Opt.NoOfBachSequential));
    fprintf(TSEMO_log,'\n %s %d', 'Genetic algorithm population size:       ',Opt.pop);
    fprintf(TSEMO_log,'\n %s %d \n', 'Genetic algorithm number of generations: ',Opt.Generation);
    fprintf(TSEMO_log,strcat('\n','%s',string6,'\n'),'                                         ',string5{:});
    fprintf(TSEMO_log,strcat('%s',string7,'\n'), ' Number of spectral sampling points:     ',Opt.GP(1:size(Y,2)).nSpectralpoints);
    fprintf(TSEMO_log,strcat('%s',string7), ' Type of matern function:                ',Opt.GP(1:size(Y,2)).matern);
    fprintf(TSEMO_log,strcat('\n','%s',string7), ' Direct evaluations per input dimension: ',Opt.GP(1:size(Y,2)).fun_eval);
    fprintf(TSEMO_log,'\n %s', '¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯');

    fprintf(TSEMO_log,'\n %s \n', 'Initial data set');
    fprintf(TSEMO_log,'\n %s %d \n', 'Number of initial data points: ',size(X,1));
    fprintf(TSEMO_log,'\n %s', 'Initial input data matrix:');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,strcat(string1,'\n'),X');
    fprintf(TSEMO_log,'\n %s', 'Initial output data matrix:');
    fprintf(TSEMO_log,strcat('\n',string6,'\n'),string5{:});
    fprintf(TSEMO_log,strcat(string4,'\n'),Y');
    fprintf(TSEMO_log,'\n %s', '¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯');
    fclose(TSEMO_log) ;
catch
    warning('There was an error when writing the log file. Maybe the file was currently open in another program. The algorithm continues to run but some parts of the log file are maybe not correct.') ;
    fclose('all') ;
end
return

function update_log_file(it,hv_imp,toc,xnewtrue,ytrue,Opt,Y,ub,lb)
try
    string1 = '';
    string2 = {};
    string3 = '';
    string4 = '';
    string5 = {};
    string6 = '';
    string7 = '';
    string8 = {};
    for i = 1:Opt.Gen.NoOfInputDim
        string1 = strcat(string1,'%8.4f ');
        string2 = {string2{:},strcat('x',num2str(i))};
        string3 = strcat(string3,'%+8s ');
        string8 = {string8{:},strcat('lambda',num2str(i))};
    end
    string8 = {string8{:},'sigmaf'};
    string8 = {string8{:},'sigman'};
    for i = 1:Opt.Gen.NoOfGPs
        string4     = strcat(string4,'%8.4f ');
        string5     = {string5{:},strcat('f',num2str(i))};
        string6     = strcat(string6,'%+8s ');
        string7     = strcat(string7,'%8d ');
        hypcov      = exp(Opt.GP(i).hyp.cov);
        hypmat(:,i) = [hypcov(1:end-1).*(1./(ub-lb))';hypcov(end)*std(Y(:,i));exp(Opt.GP(i).hyp.lik)*std(Y(:,i))];
    end

    TSEMO_log = fopen( 'TSEMO_log.txt', 'a');
    fprintf(TSEMO_log,'\n %s %d \n', 'Algorithm iteration',it);
    fprintf(TSEMO_log,'\n %s %8.4f', 'Predicted hypervolume improvement: ',hv_imp);
    fprintf(TSEMO_log,'\n %s %8.4f \n', 'Time taken: ',toc);
    fprintf(TSEMO_log,'\n %s', 'Proposed evaluation point(s): ');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,strcat(string1,'\n'),xnewtrue');
    fprintf(TSEMO_log,'\n %s', 'Corresponding observation(s): ');
    fprintf(TSEMO_log,strcat('\n',string6,'\n'),string5{:});
    fprintf(TSEMO_log,strcat(string4,'\n'),ytrue);
    fprintf(TSEMO_log,'\n %s', 'Current hyperparameter values: ');
    fprintf(TSEMO_log,strcat('\n','%+16s',string6,'\n'),'Hyperparameter',string5{:});
    for i = 1:Opt.Gen.NoOfInputDim+2
        fprintf(TSEMO_log,strcat('%+16s',string4,'\n'),string8{i},hypmat(i,:));
    end
    fprintf(TSEMO_log,'\n %s', '¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯');
    fclose(TSEMO_log) ;
catch
    warning('There was an error when writing the log file. Maybe the file was currently open in another program. The algorithm continues to run but some parts of the log file are maybe not correct.') ;
    fclose('all') ;
end
return

function final_log_update(Xpareto,Ypareto,X,Y,XParetoGP,YParetoGP,hypf,Opt)
try
    string1 = '';
    string2 = {};
    string3 = '';
    string4 = '';
    string5 = {};
    string6 = '';
    string7 = '';
    string8 = {};
    for i = 1:Opt.Gen.NoOfInputDim
        string1 = strcat(string1,'%8.4f ');
        string2 = {string2{:},strcat('x',num2str(i))};
        string3 = strcat(string3,'%+8s ');
        string8 = {string8{:},strcat('lambda',num2str(i))};
    end
    string8 = {string8{:},'sigmaf'};
    string8 = {string8{:},'sigman'};
    for i = 1:Opt.Gen.NoOfGPs
        string4     = strcat(string4,'%8.4f ');
        string5     = {string5{:},strcat('f',num2str(i))};
        string6     = strcat(string6,'%+8s ');
        string7     = strcat(string7,'%8d ');
    end

    TSEMO_log = fopen( 'TSEMO_log.txt', 'a');
    fprintf(TSEMO_log,'\n %s \n', 'Final algorithm output');
    fprintf(TSEMO_log,'\n %s', 'Final input data matrix:');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,strcat(string1,'\n'),X');
    fprintf(TSEMO_log,'\n %s', 'Final output data matrix:');
    fprintf(TSEMO_log,strcat('\n',string6,'\n'),string5{:});
    fprintf(TSEMO_log,strcat(string4,'\n'),Y');
    fprintf(TSEMO_log,'\n %s', 'Input data matrix of corresponding Pareto front:');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,strcat(string1,'\n'),Xpareto');
    fprintf(TSEMO_log,'\n %s', 'Output data matrix of corresponding Pareto front:');
    fprintf(TSEMO_log,strcat('\n',string6,'\n'),string5{:});
    fprintf(TSEMO_log,strcat(string4,'\n'),Ypareto');
    fprintf(TSEMO_log,'\n %s', 'Input data matrix of Pareto front of final Gaussian process model:');
    fprintf(TSEMO_log,strcat('\n',string3,'\n'),string2{:});
    fprintf(TSEMO_log,strcat(string1,'\n'),XParetoGP');
    fprintf(TSEMO_log,'\n %s', 'Output data matrix of Pareto front of final Gaussian process model:');
    fprintf(TSEMO_log,strcat('\n',string6,'\n'),string5{:});
    fprintf(TSEMO_log,strcat(string4,'\n'),YParetoGP');
    fprintf(TSEMO_log,'\n %s', 'Final hyperparameter values: ');
    fprintf(TSEMO_log,strcat('\n','%+16s',string6,'\n'),'Hyperparameter',string5{:});
    for i = 1:Opt.Gen.NoOfInputDim+2
        fprintf(TSEMO_log,strcat('%+16s',string4,'\n'),string8{i},hypf(i,:));
    end
    fclose(TSEMO_log) ;
catch
    warning('There was an error when writing the log file. Maybe the file was currently open in another program. The algorithm continues to run but some parts of the log file are maybe not correct.') ;
    fclose('all') ;
end
return

function [Xnew,Ynew] = ScaleVariables(X,Y,lb,ub,Opt)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

%% Scales input and output variabels
Xnew = zeros(size(X)) ; % scaled inputs
Ynew = zeros(size(Y)) ; % scaled outputs

%% Scale input variables to [0,1]
for i = 1 : size(X,2)
    Xnew(:,i) = (X(:,i)-lb(i)) / (ub(i)-lb(i)) ;
end

%% Scale output variables to zero mean and unit variance
MeanOfOutputs = zeros(Opt.Gen.NoOfGPs,1) ;
stdOfOutputs = zeros(Opt.Gen.NoOfGPs,1) ;
for i = 1 : size(Y,2)
    MeanOfOutputs(i) = mean(Y(:,i)); % calculate mean
    stdOfOutputs(i) = std(Y(:,i)) ; % calculate standard deviation
    Ynew(:,i) = (Y(:,i) - MeanOfOutputs(i)) / stdOfOutputs(i) ; % scale outputs
end
return

function [OptGPhyp] = TrainingOfGP(Xnew,Ynew,OptGP)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

% Function which minimizes the neg-loglikelihood to find hyperparameters
%% Initialize variables
Opt.GP = OptGP ;
[n,D] = size(Xnew) ;

%% Set initial hyperparameters
h1 = Opt.GP.h1; % number of hyperparameters from covariance
h2 = Opt.GP.h2; % number of hyperparameters from likelihood

%% Calculation of squared-distance matrix
a = Xnew' ;
K_M = zeros(n,n*D) ;
for i = 1:D
    K_M(:,(i-1)*n+1:i*n) = sqdist(a(i,:),a(i,:)) ;
end

%% Minimize log-negative likeliehood
% Objective Function
obj_fun.f = @(hypVar) NLikelihood(hypVar,Xnew,Ynew,K_M,Opt.GP);

% Define bounds
lb              = ones(h1+h2,1) *  log(sqrt(10^(-3))) ;  % see Jones paper
ub              = ones(h1+h2,1) *  log(sqrt(10^(3)))  ;  % see Jones paper
lb(h1+h2)    = -6;
ub(h1+h2)    = Opt.GP.noiselimit;
bounds = [lb,ub];
opts.maxevals = Opt.GP.fun_eval*(h1+h2);
opts.maxits =  100000*(h1+h2);
opts.maxdeep = 100000*(h1+h2);
opts.showits = 0;

% Defintion of options for global search
[~,x0] = Direct(obj_fun,bounds,opts);

% Defintion of options for fmincon solver
LSoptions.Algorithm = 'interior-point';
LSoptions.DerivativeCheck = 'off';
LSoptions.TolCon = 1e-12;
LSoptions.Display = 'off';
LSoptions.Hessian = 'bfgs';
LSoptions.TolFun = 1e-12;
LSoptions.PlotFcns = [];
LSoptions.GradConstr = 'off';
LSoptions.GradObj = 'on';
LSoptions.TolX = 1e-14;
LSoptions.UseParallel = 0;

% Solve optimization problem
hypResult = fmincon(obj_fun.f,x0,[],[],[],[],lb,ub,[],LSoptions);

%% Return optimal hyperparameters
Opt.GP.hyp.cov  = hypResult(1:h1);
Opt.GP.hyp.lik  = hypResult(h1+1:h1+h2);
OptGPhyp = Opt.GP.hyp ;
return

function [NLL,dNLL] = NLikelihood(hypVar, Xnew, Ynew, K_M, OptGP)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

% Calculates the log-negative likelihood
%% Initialize variables
[n,D]       = size(Xnew) ;
Opt.GP      = OptGP ;
if Opt.GP.cov ~= inf
    d           = Opt.GP.cov;    % type of Martern
else
    d           = 1;
end
h1          = Opt.GP.h1 ;    % number of hyperparameters from covariance
h2          = Opt.GP.h2;     % number of hyperparameters from likelihood
hyp.cov     = hypVar(1:h1);
hyp.lik     = hypVar(h1+1:h1+h2);
ell         = exp(hyp.cov(1:D));
sf2         = exp(2*hypVar(D+1));
K           = zeros(n,n) ;

%% Calculate covariance matrix
for i = 1:D
    K = K_M(:,(i-1)*n+1:i*n) * d/ell(i)^2 + K;
end

if Opt.GP.cov ~= inf
    sqrtK = sqrt(K) ;
    expnK = exp(-sqrtK) ;
else
    expnK = exp(-1/2*K);
    sqrtK = [];
end

if      Opt.GP.cov == 3, t = sqrtK ; m =  (1 + t).*expnK;
elseif  Opt.GP.cov == 1,             m =  expnK;
elseif  Opt.GP.cov == 5, t = sqrtK ; m =  (1 + t.*(1+t/3)).*expnK;
elseif  Opt.GP.cov == inf,           m  = expnK;
end
K = sf2*m;
K =  K + eye(n)*exp(hyp.lik*2) ;
K = (K+K')/2 ; % This guarantees a symmetric matrix

%% Calculate inverse of covariance matrix
try
    CH = chol(K) ;
    invK = CH\(CH'\eye(n));
catch
    CH = chol(K+eye(n)*1e-4);
    invK = CH\(CH'\eye(n));
    warning('Covariance matrix in Nlikelihood is not positive semi-definite')
end

%% Calculate determinant of covariance matrix
logDetK = 2*sum(log(abs(diag(CH)))) ;

%% Calculate hyperperpriors
logprior = 0 ;
dlogpriorcov = zeros(1,h1) ;
for i = 1 : h1
    [A, dlogpriorcov(i)] =  priorGauss(Opt.GP.priorcov(1),Opt.GP.priorcov(2), hyp.cov(i) ) ;
    logprior = logprior + A ;
end

dlogpriorlik = zeros(1,h2);
for i = 1 : h2
    [A, dlogpriorlik(i)] =  priorGauss(Opt.GP.priorlik(1),Opt.GP.priorlik(2), hyp.lik(i) ) ;
    logprior = logprior + A ;
end

%% Calculate negative log-likeliehood
NLL = n/2*log(2*pi) + 1/2 * logDetK + 1/2 * Ynew'*invK*Ynew - logprior ;

%% Gradient calculation
if nargout == 2 % do only if No of output variable is 2 (if necessary)
    dsq_M = zeros(n,n*D) ;
    for i = 1 : D
        dsq_M(:,(i-1)*n+1:i*n)  = K_M(:,(i-1)*n+1:i*n) * (d)/ell(i)^2 ;
    end

    c = invK*Ynew;
    for i = 1 : h1
        dK = covMaternanisotropic(Opt.GP.cov,hyp.cov, sqrtK, expnK, dsq_M, Xnew, [], i);
        b = invK* dK ;
        dNLL_f.cov(i) = 1/2*trace(b) - 1/2*Ynew'*b*c ;
    end

    for i = 1 : h2
        dK = 2 * exp(hyp.lik(i)) * eye(n) * exp(hyp.lik(i));
        b = invK* dK ;
        dNLL_f.lik(i) = 1/2*trace(b) - 1/2*Ynew'*b*c ;
    end

    dNLL = [dNLL_f.cov';dNLL_f.lik'] - [dlogpriorcov';dlogpriorlik'];
end

return

function f = posterior_sample(Xnew,Ynew,Opt)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2018-08-07.

% extration of variables from problem structure
nSpectralpoints = Opt.nSpectralpoints;
[n,D] = size(Xnew);
ell = exp(Opt.hyp.cov(1:D));
sf2 = exp(2*Opt.hyp.cov(D+1));
sn2 = exp(2*Opt.hyp.lik);

% Sampling of W and b
sW1  = lhsdesign(nSpectralpoints,D,'criterion','none');
sW2  = lhsdesign(nSpectralpoints,D,'criterion','none');
if Opt.cov ~= inf
    W = repmat(1./(ell)', nSpectralpoints, 1).*norminv(sW1).*sqrt(Opt.cov./chi2inv(sW2,Opt.cov));
else
    W = randn(nSpectralpoints,D) .* repmat(1./ell', nSpectralpoints, 1);
end

b = 2*pi*lhsdesign(nSpectralpoints,1,'criterion','none');

% Calculation of phi
phi = sqrt(2 * sf2 / nSpectralpoints) * cos(W * Xnew' + repmat(b, 1, n));

% Sampling of theta according to phi
A = phi * phi' + sn2 * eye(nSpectralpoints);

invA      = invChol(A);
mu_theta  = invA*phi*Ynew;
cov_theta = sn2*invA;
cov_theta = (cov_theta+cov_theta')/2;
theta     = mvnrnd(mu_theta,cov_theta)';

% Posterior sample (function) according to theta
f = @(x) (theta' * sqrt(2 * sf2 / nSpectralpoints) * cos(W * x' + repmat(b,1,size(x,1))))';
return

function [f,varf] = mean_sample(Xnew,Ynew,Opt)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2020-23-05.

% extration of variables from problem structure
nSpectralpoints = Opt.nSpectralpoints;
[n,D] = size(Xnew);
ell = exp(Opt.hyp.cov(1:D));
sf2 = exp(2*Opt.hyp.cov(D+1));
sn2 = exp(2*Opt.hyp.lik);

% Sampling of W and b
sW1  = lhsdesign(nSpectralpoints,D,'criterion','none');
sW2  = lhsdesign(nSpectralpoints,D,'criterion','none');
if Opt.cov ~= inf
    W = repmat(1./(ell)', nSpectralpoints, 1).*norminv(sW1).*sqrt(Opt.cov./chi2inv(sW2,Opt.cov));
else
    W = randn(nSpectralpoints,D) .* repmat(1./ell', nSpectralpoints, 1);
end

b = 2*pi*lhsdesign(nSpectralpoints,1,'criterion','none');

% Calculation of phi
phi = sqrt(2 * sf2 / nSpectralpoints) * cos(W * Xnew' + repmat(b, 1, n));

% Sampling of theta according to phi
A = phi * phi' + sn2 * eye(nSpectralpoints);

invA      = invChol(A);
mu_theta  = invA*phi*Ynew;

% Mean approximation (function) according to theta
phi_x = @(x) sqrt(2 * sf2 / nSpectralpoints) * cos(W * x' + repmat(b,1,size(x,1)));
f     = @(x) (mu_theta' * phi_x(x))';
varf  = @(x) sn2 + sn2*phi_x(x)'*invA*phi_x(x);
return

function v = hypervolumemonte(P,r,N)
% Copyright (c) 2009, Yi Cao
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% HYPERVOUME    Hypervolume indicator as a measure of Pareto front estimate.
%   V = HYPERVOLUME(P,R,N) returns an estimation of the hypervoulme (in
%   percentage) dominated by the approximated Pareto front set P (n by d)
%   and bounded by the reference point R (1 by d). The estimation is doen
%   through N (default is 1000) uniformly distributed random points within
%   the bounded hyper-cuboid.
%
%   V = HYPERVOLUMN(P,R,C) uses the test points specified in C (N by d).
%
% See also: paretofront, paretoGroup

% Version 1.0 by Yi Cao at Cranfield University on 20 April 2008

% Example
%{
% an random exmaple
F=(randn(100,3)+5).^2;
% upper bound of the data set
r=max(F);
% Approximation of Pareto set
P=paretofront(F);
% Hypervolume
v=hypervolume(F(P,:),r,100000);
%}
% https://se.mathworks.com/matlabcentral/fileexchange/19651-hypervolume-indicator

% Check input and output
error(nargchk(2,3,nargin));
error(nargoutchk(0,1,nargout));

P=P*diag(1./r);
[n,d]=size(P);
if nargin<3
    N=1000;
end
if ~isscalar(N)
    C=N;
    N=size(C,1);
else
    C=rand(N,d);
end

fDominated=false(N,1);
lB=min(P);
fcheck=all(bsxfun(@gt, C, lB),2);

for k=1:n
    if any(fcheck)
        f=all(bsxfun(@gt, C(fcheck,:), P(k,:)),2);
        fDominated(fcheck)=f;
        fcheck(fcheck)=~f;
    end
end

v=sum(fDominated)/N;
return

function K = covMaternanisotropic(d, hyp, sqrtK,expnK, dsq_M, x, z, i)
% Copyright (c) 2005-2017 Carl Edward Rasmussen & Hannes Nickisch. All rights reserved.
%
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    1. Redistributions of source code must retain the above copyright notice,
%       this list of conditions and the following disclaimer.
%    2. Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY CARL EDWARD RASMUSSEN & HANNES NICKISCH ``AS IS''
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL CARL EDWARD RASMUSSEN & HANNES NICKISCH OR CONTRIBUTORS BE LIABLE
% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
% ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
% EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% The views and conclusions contained in the software and documentation
% are those of the authors and should not be interpreted as representing official policies,
% either expressed or implied, of Carl Edward Rasmussen & Hannes Nickisch.
%
% The code and associated documentation is available from http://gaussianprocess.org/gpml/code.</pre>

[n,D] = size(x);
sf2 = exp(2*hyp(D+1));

if nargin<7                                                        % covariances
    if      d == 3, t = sqrtK ; m =  (1 + t).*expnK;
    elseif  d == 1,             m =  expnK;
    elseif  d == 5, t = sqrtK ; m =  (1 + t.*(1+t/3)).*expnK;
    elseif  d == inf, m = expnK;
    end
    K = sf2*m;
else                                                               % derivatives
    if i<=D                                               % length scale parameter
        Ki = dsq_M(:,(i-1)*n+1:i*n) ;
        if     d == 3,             dm = expnK;
        elseif d == 1, t = sqrtK ; dm = (1./t).*expnK;
        elseif d == 5, t = sqrtK ; dm = ((1+t)/3).*expnK;
        elseif d == inf; dm = -1/2*expnK;
        end

        K = sf2*dm.*Ki;
        K(Ki<1e-12) = 0;                                    % fix limit case for d=1
    elseif i==D+1                                            % magnitude parameter
        if      d == 3, t = sqrtK ; m =  (1 + t).*expnK;
        elseif  d == 1,             m =  expnK;
        elseif  d == 5, t = sqrtK ; m =  (1 + t.*(1+t/3)).*expnK;
        elseif  d == inf,           m = expnK;
        end
        K = 2*sf2*m;
    end
end
return

function D = sqdist(X1, X2)
% Copyright (c) 2016, Mo Chen
% All rights reserved.
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% https://se.mathworks.com/matlabcentral/fileexchange/24599-pairwise-distance-matrix
% Pairwise square Euclidean distance between two sample sets
% Input:
%   X1, X2: dxn1 dxn2 sample matrices
% Output:
%   D: n1 x n2 square Euclidean distance matrix
% Written by Mo Chen (sth4nth@gmail.com).

D = bsxfun(@plus,dot(X2,X2,1),dot(X1,X1,1)')-2*(X1'*X2);
D(D<0) = 0 ; % check due to numerical errors
return

function [lp,dlp] = priorGauss(mu,s2,x)
% Copyright (c) 2005-2017 Carl Edward Rasmussen & Hannes Nickisch. All rights reserved.
%
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    1. Redistributions of source code must retain the above copyright notice,
%       this list of conditions and the following disclaimer.
%    2. Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY CARL EDWARD RASMUSSEN & HANNES NICKISCH ``AS IS''
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL CARL EDWARD RASMUSSEN & HANNES NICKISCH OR CONTRIBUTORS BE LIABLE
% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
% ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
% EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% The views and conclusions contained in the software and documentation
% are those of the authors and should not be interpreted as representing official policies,
% either expressed or implied, of Carl Edward Rasmussen & Hannes Nickisch.
%
% The code and associated documentation is available from http://gaussianprocess.org/gpml/code.</pre>

lp  = -(x-mu).^2/(2*s2) - log(2*pi*s2)/2;
dlp = -(x-mu)/s2;
return

function [f,const] = pareto_objective(x,OptSample)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

Opt.Sample = OptSample;
f = zeros(size(x,1),size(Opt.Sample,2));
for i = 1:size(Opt.Sample,2)
    f(:,i) = Opt.Sample(i).f(x);
end
const = [];
return

function [Sample_pareto,Sample_xpareto,Sample_nadir] = Find_sample_pareto(Opt,it)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

D = Opt.Gen.NoOfInputDim;               % Number of input dimensions
options = nsgaopt();                    % create default options structure
options.popsize = Opt.pop;              % populaion size
options.maxGen  = Opt.Generation;       % max generation
options.numObj = Opt.Gen.NoOfGPs;       % number of objectives
options.numVar = D;                     % number of design variables
options.numCons = 0;                    % number of constraints
options.outputfuns = [];                % saving pop
options.lb = zeros(1,D);                % lower bound of x
options.ub = ones(1,D);                 % upper bound of x
options.objfun = @(x) pareto_objective(x,Opt.Sample);        % objective function handle
options.useParallel = 'no';             % parallel computation is non-essential here
if it > 1
options.initfun = {@(opt,pop) initpop_new(opt,pop,Opt)};
end
[~,result] = evalc('nsga2(options);');  % begin the optimization!

Sample_xpareto = zeros(Opt.pop,D);
Sample_pareto = zeros(Opt.pop,Opt.Gen.NoOfGPs);
result = result.pops(Opt.Generation,:);

for k = 1:Opt.pop
    Sample_xpareto(k,:) = result(k).var;
    Sample_pareto(k,:) = result(k).obj;
end

for k = 1:Opt.Gen.NoOfGPs
    Sample_nadir(k) = max(Sample_pareto(:,k));
end
return

function pop = initpop_new(opt,pop,Opt)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2020-22-05.

for p = 1:opt.popsize
    pop(p).var = Opt.warmstart_pareto(p,:);
end

return

function [Mean_pareto,Mean_xpareto] = Find_mean_pareto(Opt)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

D = Opt.Gen.NoOfInputDim;               % Number of input dimensions
options = nsgaopt();                    % create default options structure
options.popsize = Opt.pop;              % populaion size
options.maxGen  = Opt.Generation;       % max generation
options.numObj = Opt.Gen.NoOfGPs;       % number of objectives
options.numVar = D;                     % number of design variables
options.numCons = 0;                    % number of constraints
options.outputfuns = [];                % saving pop
options.lb = zeros(1,D);                % lower bound of x
options.ub = ones(1,D);                 % upper bound of x
options.objfun = @(x) pareto_objective(x,Opt.Mean); % objective function handle
options.useParallel = 'no';             % parallel computation is non-essential here
[~,result] = evalc('nsga2(options);');  % begin the optimization!

Mean_xpareto = zeros(Opt.pop,D);
Mean_pareto  = zeros(Opt.pop,Opt.Gen.NoOfGPs);
result       = result.pops(Opt.Generation,:);

for k = 1:Opt.pop
    Mean_xpareto(k,:) = result(k).var;
    Mean_pareto(k,:) = result(k).obj;
end

return

function hv = hypervolume_2D(Yfront,r)
% The mex-file used follows the description of:
% 'M. Emmerich, K. Yang, A. Deutz, H. Wang and C. M. Fonseca. A
% Multicriteria Generalization of Bayesian Global Optimization.'
% Group website: http://liacs.leidenuniv.nl/~csmoda/index.php?page=code

AYfront = remove_points_above_reference(Yfront,r);
if isempty(AYfront)
    hv = 0;
else
    normvec = min(AYfront,[],1);
    A =  AYfront-repmat(normvec,size(AYfront,1),1);
    A =  A * diag(1./(r-normvec));
    A = -A + ones(size(A));
    A = sortrows(A,2);
    hyp_percentage = hypervolume2D(A,[0,0]);
    hv = prod(r-normvec)*hyp_percentage;
end
return

function hv = hypervolume_3D(Yfront,r)
% The mex-file used follows the description of:
% 'K. Yang, M. Emmerich, A. Deutz and C. M. Fonseca. A
% Computing 3-D Expected Hypervolume Improvement and Related Integrals in
% Asymptotically Optimal Time.'
% Group website: http://liacs.leidenuniv.nl/~csmoda/index.php?page=code

AYfront = remove_points_above_reference(Yfront,r);
if isempty(AYfront)
    hv = 0;
else
    normvec = min(AYfront,[],1);
    A = AYfront-repmat(normvec,size(AYfront,1),1) ;
    A = A *diag(1./(r-normvec));
    A = -A + ones(size(A)) ;
    A = sortrows(A,3);
    hyp_percentage = hypervolume3D(A,[0,0,0],[1,1,1]);
    hv = prod(r-normvec)*hyp_percentage;
end
return

function [index,hv_imp] = hypervolume_improvement_index(Ynew,Sample_nadir,Sample_pareto,Opt)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

r = Sample_nadir + 0.01*(max(Sample_pareto)-min(Sample_pareto));
index = [];

for i = 1 : Opt.NoOfBachSequential

    Yfront = Ynew(paretofront(Ynew),:);

    if size(Ynew,2) == 2
        hvY = hypervolume_2D(Yfront,r);

        for k = 1:size(Sample_pareto,1)
            A = [Ynew;Sample_pareto(k,:)];
            Afront = A(paretofront(A),:);
            hv = hypervolume_2D(Afront,r);
            hv_improvement(k) = hv-hvY;
        end

    elseif size(Ynew,2) == 3
        hvY = hypervolume_3D(Yfront,r);

        for k = 1:size(Sample_pareto,1)
            A = [Ynew;Sample_pareto(k,:)];
            Afront = A(paretofront(A),:);
            hv = hypervolume_3D(Afront,r);
            hv_improvement(k) = hv-hvY;
        end

    else
        AYfront = remove_points_above_reference(Yfront,r);
        normvec = min(AYfront,[],1);
        hyp_percentage = hypervolumemonte(AYfront-repmat(normvec,size(AYfront,1),1),r-normvec,3000);
        hvY = prod(r-normvec)*hyp_percentage;

        hv_improvement = zeros(size(Sample_pareto,1),1);
        for k = 1:size(Sample_pareto,1)
            B = [Ynew;Sample_pareto(k,:)];
            Bfront = B(paretofront(B),:);
            ABfront = remove_points_above_reference(Bfront,r);
            if isempty(ABfront)
                hv_improvement(k) = 0;
            else
                normvec = min(ABfront,[],1);
                hyp_percentage = hypervolumemonte(ABfront-repmat(normvec,size(ABfront,1),1),r-normvec,10000);
                hv = prod(r-normvec)*hyp_percentage;
                hv_improvement(k) = hv-hvY;
            end
        end
    end

    if i == 1
        hvY0 = hvY;
    end

    [~,Currentindex] = max(hv_improvement);
    Ynew = [Ynew;Sample_pareto(Currentindex,:)];
    index = [index;Currentindex];
end
hv_imp = hv_improvement(index(end))+hvY-hvY0;
return

function A = remove_points_above_reference(Afront,r)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

[A,~] = sortrows(Afront);
for p = 1:size(Afront,2)
    A = A(A(:,p)<=r(p),:);
end
return

function Opt = set_option_structure(Opt,X,Y)
% Copyright (c) by Eric Bradford, Artur M. Schweidtmann and Alexei Lapkin, 2017-13-12.

%% Extraction of input and output dimensions from X and Y
Opt.Gen.NoOfGPs = size(Y,2);           % Number of GPs to be trained
Opt.Gen.NoOfInputDim = size(X,2);      % Number of inputs dimensions
for i = 1 : Opt.Gen.NoOfGPs
    %% Set up GP options
    Opt.GP(i).cov = Opt.GP(i).matern;      % Matern type 1 / 3 / 5 / inf

    %% Set hyperpriors (MAP)
    Opt.GP(i).noiselimit = 0;                      % Upper bound on noise
    Opt.GP(i).var        = 10;                     % Upper bound on signal variance
    Opt.GP(i).h1         = Opt.Gen.NoOfInputDim+1; % Number of hyperparameters from covariance
    Opt.GP(i).h2         = 1;                      % Number of hyperparameters from likelihood

    %% priorGauss (mean, var)
    Opt.GP(i).priorlik  = [-6 ,Opt.GP(i).var];
    Opt.GP(i).priorcov  = [ 0 ,Opt.GP(i).var];

    %% Initial values for hyperparameters
    Opt.GP(i).hyp.cov       = zeros(1, Opt.GP(i).h1);
    Opt.GP(i).hyp.lik       = log(1e-2);
end
return

author: Eric-Bradford

repo: TS-EMO


text:
function [fMc, mcCalculator] = calc_Mc(mCatalog, calcMethod, binInterval, mcCorrectionFactor)
    % CALC_MC Calculates the magnitude of completeness for a given catalog
    %
    % fMc = CALC_MC(mCatalog, nMethod, binInterval, mcCorrectionFactor)
    % [fMc, mc_calculator] = CALC_MC(...) returns a function handle to the calculation
    % method, so it can be reused in heavy loops.  MC_CALCULATOR has the form:
    %    fMc =  MY_CALCULATOR(catalog, bins, correction);
    %
    % --------------------------------------------------------------------
    %
    %
    % Input parameters:
    %   mCatalog       Earthquake catalog for determing the magnitude of completeness
    %   calcMethod        Method to determine the magnitude of completeness
    %                     see: McMethods for a list of valid values
    %   binInterval       Binning of catalog's magnitudes (default 0.1)
    %   mcCorrectionFactor  Correction term to be added to fMc (default 0)
    %
    % Output parameters:
    %   fMc            Magnitude of completeness
    %   mc_calculator
    %
    %
    % Copyright (C) 2004 by Danijel Schorlemmer, Jochen Woessner

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

    % Magnitude binning
    if ~exist('binInterval', 'var') || isempty(binInterval)
        binInterval = 0.1;
    end
    if ~isa(calcMethod,'McMethods')
        error('Expected an actual method (McMethods) but received something else.\n See McMethods');
    end
    % Correction
    if ~exist('mcCorrectionFactor', 'var') || isempty(mcCorrectionFactor)
        mcCorrectionFactor = 0;
    end
    useNumeric = isnumeric(mCatalog) || (ischarlike(mCatalog) && mCatalog == "AsMagnitudes");
    if useNumeric
        switch calcMethod
            case McMethods.MaxCurvature
                methodFun = @calc_McMaxCurvature;

            case  McMethods.FixedMc
                methodFun = @min;

            case  McMethods.Mc90
                methodFun = @(C)method_mc90(C,binInterval);

            case  McMethods.Mc95
                methodFun = @(C)method_mc95(C,binInterval);

            case McMethods.McBestCombo
                methodFun = @(C)method_bestcombo(C,binInterval);

            case  McMethods.McEMR
                error('this only works with full catalogs')
                %methodFun = @(C)calc_McEMR(C,binInterval); %requires catalog Magnitude AND Date

            case  McMethods.McDueB_ShiBolt
                methodFun = @(C)calc_Mcdueb(C, binInterval);

            case McMethods.McDueB_Bootstrap
                nSample     = 500;
                nWindowSize = 5;
                nMinEvents  = 50;
                methodFun   = @(C) calc_McduebBst(C, binInterval, nWindowSize, nMinEvents,nSample);

            case McMethods.McDueB_Cao
                methodFun = @(C, ~)calc_McduebCao(C);

            otherwise
                error('unknown Mc method');
        end
    else % catalog object
        switch calcMethod
            case McMethods.MaxCurvature
                methodFun = @(C)calc_McMaxCurvature(C.Magnitude);

            case  McMethods.FixedMc
                methodFun = @(C) min(C.Magnitude);

            case  McMethods.Mc90
                methodFun = @(C)method_mc90(C.Magnitude,binInterval);

            case  McMethods.Mc95
                methodFun = @(C)method_mc95(C.Magnitude,binInterval);

            case McMethods.McBestCombo
                methodFun = @(C)method_bestcombo(C.Magnitude,binInterval);

            case  McMethods.McEMR
                methodFun = @(C)calc_McEMR(C,binInterval); %requires catalog Magnitude AND Date

            case  McMethods.McDueB_ShiBolt
                methodFun = @(C)calc_Mcdueb(C.Magnitude, binInterval);

            case McMethods.McDueB_Bootstrap
                nSample     = 500;
                nWindowSize = 5;
                nMinEvents  = 50;
                methodFun   = @(C) calc_McduebBst(C.Magnitude, binInterval, nWindowSize, nMinEvents,nSample);

            case McMethods.McDueB_Cao
                methodFun = @(C, ~)calc_McduebCao(C.Magnitude);

            otherwise
                error('unknown Mc method');
        end
    end

    % lock the method into this calculation
    mcCalculator = @(C) do_calculation(methodFun, C, mcCorrectionFactor);

    if ~isempty(mCatalog)
        % do the calculation
        fMc = mcCalculator(mCatalog);
    else
        fMc = [];
    end

end

function fMc = do_calculation(methodFun, mCatalog, mcCorrectionFactor)
    if isempty(mCatalog) || ischarlike(mCatalog)
        fMc = nan;
        return
    end

    fMc = methodFun(mCatalog);

    % Check fMc
    if isempty(fMc)
        fMc = nan;
    end

    % Apply correction
    fMc = fMc + mcCorrectionFactor;
end

% functions to return that all have same signature
function fMc = method_mc90(mags, binInterval)
    [~, ~, fMc] = calc_McBest(mags, binInterval);
end

function fMc = method_mc95(mags, binInterval)
    [~, fMc] = calc_McBest(mags, binInterval);
end

function fMc = method_bestcombo(mags, binInterval)
        [~, Mc95, Mc90] = calc_McBest(mags, binInterval);
        if ~isnan(Mc95)
            fMc = Mc95;
        elseif ~isnan(Mc90)
            fMc = Mc90;
        else
            fMc = calc_McMaxCurvature(mags);
        end
end

author: CelsoReyes


text:
% FITELLIPSE  Least-squares fit of ellipse to 2D points.
%        A = FITELLIPSE(X, Y) returns the parameters of the best-fit
%        ellipse to 2D points (X, Y).
%        The returned vector A contains the center, radii, and orientation
%        of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)
%

% Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999
%
%  @Article{Fitzgibbon99,
%   author = "Fitzgibbon, A.~W.and Pilu, M. and Fisher, R.~B.",
%   title = "Direct least-squares fitting of ellipses",
%   journal = pami,
%   year = 1999,
%   volume = 21,
%   number = 5,
%   month = may,
%   pages = "476--480"
%  }
%
% This is a more bulletproof version than that in the paper, incorporating
% scaling to reduce roundoff error, correction of behaviour when the input
% data are on a perfect hyperbola, and returns the geometric parameters
% of the ellipse, rather than the coefficients of the quadratic form.
%
%  See also http://research.microsoft.com/en-us/um/people/awf/ellipse/
%
function a = fitEllipse(X, Y)
    % normalize data
    mx = mean(X);
    my = mean(Y);
    sx = (max(X)-min(X))/2;
    sy = (max(Y)-min(Y))/2;

    x = (X-mx)/sx;
    y = (Y-my)/sy;

    % Force to column vectors
    x = x(:);
    y = y(:);

    % Build design matrix
    D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

    % Build scatter matrix
    S = D'*D;

    % Build 6x6 constraint matrix
    C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

    % Solve eigensystem
    if 0
        % Old way, numerically unstable if not implemented in matlab
        [gevec, geval] = eig(S,C);

        % Find the negative eigenvalue
        I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));

        % Extract eigenvector corresponding to negative eigenvalue
        A = real(gevec(:,I));
    else
        % New way, numerically stabler in C [gevec, geval] = eig(S,C);

        % Break into blocks
        tmpA = S(1:3,1:3);
        tmpB = S(1:3,4:6);
        tmpC = S(4:6,4:6);
        tmpD = C(1:3,1:3);
        tmpE = inv(tmpC)*tmpB';
        [evec_x, eval_x] = eig(inv(tmpD) * (tmpA - tmpB*tmpE));

        % Find the positive (as det(tmpD) < 0) eigenvalue
        I = find(real(diag(eval_x)) < 1e-8 & ~isinf(diag(eval_x)));

        % Extract eigenvector corresponding to negative eigenvalue
        A = real(evec_x(:,I));

        % Recover the bottom half...
        evec_y = -tmpE * A;
        A = [A; evec_y];
    end

    % unnormalize
    par = [
    A(1)*sy*sy,   ...
        A(2)*sx*sy,   ...
        A(3)*sx*sx,   ...
        -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
        -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
        A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
        - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
        + A(6)*sx*sx*sy*sy   ...
        ]';

    % Convert to geometric radii, and centers

    thetarad = 0.5*atan2(par(2),par(1) - par(3));
    cost = cos(thetarad);
    sint = sin(thetarad);
    sin_squared = sint.*sint;
    cos_squared = cost.*cost;
    cos_sin = sint .* cost;

    Ao = par(6);
    Au =   par(4) .* cost + par(5) .* sint;
    Av = - par(4) .* sint + par(5) .* cost;
    Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;
    Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;

    % ROTATED = [Ao Au Av Auu Avv]

    tuCentre = - Au./(2.*Auu);
    tvCentre = - Av./(2.*Avv);
    wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;

    uCentre = tuCentre .* cost - tvCentre .* sint;
    vCentre = tuCentre .* sint + tvCentre .* cost;

    Ru = -wCentre./Auu;
    Rv = -wCentre./Avv;

    Ru = sqrt(abs(Ru)).*sign(Ru);
    Rv = sqrt(abs(Rv)).*sign(Rv);

    a = [uCentre, vCentre, Ru, Rv, thetarad];
end

author: kavli-ntnu


text:
%% Colour Filter
%% This function can be used for separating or filter out Red components,
%% Green components and Blue Components of colors from the color images.
%% Function C = colorfilter(img,color)
%% input    img = Color image (The input image should be a color image)
%%          color = which color want to be filtered out from the color image
%%                  it may be 'R' / 'r' / 'G' / 'g' / 'B' / 'b'
%%  Example: C = colorfilter(img,'r');
%%      Posted date   : 26 - 06 - 2008
%%      Modified date : 08 - 07 - 2008
%%
%% Developed By : K.Kannan & Jeny Rajan
%%                  Medical Imaging Research Group (MIRG), NeST, Trivandrum.
%%
function C = colorfilter(img,color)

[row col plane] = size(img);
img = double(img);
C = zeros(row,col,plane);
if plane ~= 3
    disp('Input should be a color image');
    return;
end
factor = max(img(:)) * 0.2;
switch color
    case {'R','r'}
        for i = 1:row
            for j = 1:col
                if (img(i,j,1) > factor && img(i,j,1) == max([img(i,j,1) img(i,j,2) img(i,j,3)]))
                    C(i,j,1:3) = img(i,j,1:3);
                else
                    C(i,j,1:3) = (img(i,j,1) * 0.3) + (img(i,j,2) * 0.59) + (img(i,j,3) * 0.11);
                end
            end
        end
    case {'G','g'}
        for i = 1:row
            for j = 1:col
                if (img(i,j,2) > factor && img(i,j,2) == max([img(i,j,1) img(i,j,2) img(i,j,3)]))
                    C(i,j,1:3) = img(i,j,1:3);
                else
                    C(i,j,1:3) = (img(i,j,1) * 0.3) + (img(i,j,2) * 0.59) + (img(i,j,3) * 0.11);
                end
            end
        end
    case {'B','b'}
        for i = 1:row
            for j = 1:col
                if (img(i,j,3) > factor && img(i,j,3) == max([img(i,j,1) img(i,j,2) img(i,j,3)]))
                    C(i,j,1:3) = img(i,j,1:3);
                else
                    C(i,j,1:3) = (img(i,j,1) * 0.3) + (img(i,j,2) * 0.59) + (img(i,j,3) * 0.11);
                end
            end
        end
    otherwise
        disp('unknown method');
end
C = uint8(C);
figure,imshow(uint8(img),[]);
figure,imshow(uint8(C),[]);

author: Sable

repo: mcbench-benchmarks


text:
function ppwd = pwd_imp_circexp(pm,Npw)
%PWD_IMP_CIRCEXP converts a circular basis expansion of a sound field to its
%two-dimensional plane wave decomposition
%
%   Usage: ppwd = pwd_imp_circexp(pm,[Npw])
%
%   Input parameters:
%       pm      - circular basis expansion [N x (M+1)]
%       Npw     - number of equi-angular distributed plane waves, optional,
%                 default: 2*M+1
%
%   Output parameters:
%       ppwd    - plane wave decomposition [N x Npw]
%
%   See also: driving_function_imp_localwfs_sbl_ps,
%   driving_function_imp_localwfs_sbl_pw

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ==================================
nargmin = 1;
nargmax = 2;
narginchk(nargmin,nargmax);
isargmatrix(pm);
M = size(pm,2)-1;
if nargin == nargmin
  Npw = 2*M+1;
else
  isargpositivescalar(Npw);
end

%% ===== Computation ====================================================
% Implementation of
%                 ___
% _               \
% p(phipw, t) =   /__     p (t) i^m  e^(-i m phipw)
%               m=-M..M    m
% with
%
% phipw = n * 2*pi/Npw

pm = [conj(pm(:,end:-1:2)), pm];  % append coefficients for negative m
ppwd = inverse_cht(bsxfun(@times,pm,1i.^(-M:M)),Npw);

author: sfstoolbox

repo: sfs-matlab


text:
function h=m_streamline(long,lat,u,v,varargin)
% M_STREAMLINE Makes a quiverplot on a map (QUIVER-style)
%    M_STREAMLINE(LONG,LAT,U,V) draws well-spaced streamlines (with direction
%    arrows) with components (U,V) at the points (LONG,LAT) on the currently
%    defined map.The arrays LONG and LAT, which define the coordinates for U and
%    V, must be monotonic, but do not need to be uniformly spaced. The
%    matrices LONG,LAT,U,V must all be the same size. U and V contain the
%    eastward and northward components of velocity (in m/s or equivalent, NOT
%    degrees lat/long  per sec or equivalent). Arrow scaling is automatic.
%
%   Note - this is basically a call to STREAMSLICE, which has some limitations and
%   cannot be used for ALL projections. In particular, STREAMSLICE requires that
%
%            [X,Y]=m_ll2xy(LONG,LAT,'clip','point')
%
%   returns X/Y matrices that themselves must be monotonic and plaid (as if produced
%   by MESHGRID). This means that a) none of the points are from outside the map
%   boundaries (since the M_LL2XY call will turn them into NaN), and b) you are
%   probably limited to cylindrical projections (miller, mercator, equidistant
%   cylindrical).
%
%   M_STREAMLINE(...,density) modifies the automatic spacing of the streamlines.
%   Density must be greater than 0. The default value is 1; higher values
%   produce  more streamlines on each plane. For example, 2 produces
%   approximately twice as many streamlines, while 0.5 produces approximately
%   half as many.
%
%   M_STREAMLINE(...,'arrowsmode') determines if direction arrows are present
%   or not. arrowmode can be:
%
%   arrows   - Draw direction arrows on the streamlines (default).
%   noarrows - Do not draw direction arrows.
%
%   M_STREAMLINE(...,'method') specifies the interpolation method to use.
%   method can be
%
%   linear   - Linear interpolation (default)
%   cubic    - Cubic interpolation
%   nearest  - Nearest-neighbor interpolation
%   See interp3 for more information on interpolation methods.
%
%   M_STREAMLINE(axes_handle,...) plots into the axes object with the handle
%   axes_handle instead of into the current axes object (gca).
%
%   h = M_STREAMLINE(...) returns a vector of handles to the line objects
%   created.
%
%   [vertices arrowvertices] = M_STREAMLINE(...) returns two cell arrays of
%   vertices  for drawing the streamlines and the arrows. You can pass these
%   values to any of the streamline drawing functions (streamline, streamribbon,
%   streamtube).
%

% Shi Weiheng (tfoterye@gmail.com) 29/Jun/17
%
% Based on m_quiver written by Prof Rich Pawlowicz.
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
%  Nov/2017 - cleaned up comments.

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

[X,Y]=m_ll2xy(long,lat,'clip','point');

if any(isnan(X(:)))
  error(['M_Map : ' mfilename ' : InvalidInputs - input data includes points outside the map area']);
end

[XN ,YN ]=m_ll2xy([long(:) long(:)]',[lat(:) lat(:)+.001]','clip','off');
[XE ,YE ]=m_ll2xy([long(:) long(:)+(.001)./cos(lat(:)*pi/180)]',[lat(:) lat(:)]','clip','off');
mU=u.*reshape(diff(XE),size(lat))*1000 + v.*reshape(diff(XN),size(lat))*1000;
mV=u.*reshape(diff(YE),size(lat))*1000 + v.*reshape(diff(YN),size(lat))*1000;

h=streamslice(X,Y,mU,mV,varargin{:});
set(h,'tag','m_streamline');

if nargout==0
 clear h
end

author: goGPS-Project


text:
clear, clc;

% This is an example for running the function eplb
%
%  Problem:
%
%  min  1/2 || x - y||^2
%  s.t. \|x\|_1 <= z
%
% For detailed description of the function, please refer to the Manual.
%
%% ------------   History --------------------
% First version on August 10, 2008.
%
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1;

n=1000;
v=randn(n, 1);
z=20;
lambda0=0;

% run the function epp
[x, lambda, iter_step]=eplb(v, n, z, lambda0);

author: zzstefan

repo: BrainNetClass


text:
function sphere_triangle_quad_test03 ( )

%*****************************************************************************80
%
%% SPHERE_TRIANGLE_TEST03 tests SPHERE01_TRIANGLE_QUAD_ICOS1C.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0, 0; ...
    1, 0, 0; ...
    0, 1, 0; ...
    0, 0, 1; ...
    2, 0, 0; ...
    0, 2, 2; ...
    2, 2, 2; ...
    0, 2, 4; ...
    0, 0, 6; ...
    1, 2, 4; ...
    2, 4, 2; ...
    6, 2, 0; ...
    0, 0, 8; ...
    6, 0, 4; ...
    4, 6, 2; ...
    2, 4, 8; ...
   16, 0, 0 ]';
  n_mc1 = 1000;
  n_mc2 = 10000;
  n_mc3 = 100000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_TRIANGLE_QUAD_TEST03\n' );
  fprintf ( 1, '  SPHERE01_TRIANGLE_QUAD_ICOS1C approximates the\n' );
  fprintf ( 1, '  integral of a function over a spherical triangle on\n' );
  fprintf ( 1, '  the surface of the unit sphere using a centroid rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We do not have an exact result, so we compare each\n' );
  fprintf ( 1, '  estimate to the final one.\n' );
%
%  Choose three points at random to define a spherical triangle.
%
  [ v1, seed ] = sphere01_sample ( 1, seed );
  [ v2, seed ] = sphere01_sample ( 1, seed );
  [ v3, seed ] = sphere01_sample ( 1, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Vertices of random spherical triangle:\n' );
  fprintf ( 1, '\n' );
  r8vec_transpose_print ( 3, v1, '  V1:' );
  r8vec_transpose_print ( 3, v2, '  V2:' );
  r8vec_transpose_print ( 3, v3, '  V3:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  FACTOR   N   RESULT\n' );

  for j = 1 : 17

    e(1:3) = e_test(1:3,j);

    e = polyterm_exponent ( 'SET', e );

    polyterm_exponent ( 'PRINT', e );

%   factor = 2 ^ 11;
    factor = 2 ^ 7;
    [ best, node_num ] = sphere01_triangle_quad_icos1c ( v1, v2, v3, factor, ...
      @polyterm_value_3d );

    factor = 1;
    for factor_log = 0 : 5

      [ result, node_num ] = sphere01_triangle_quad_icos1c ( v1, v2, v3, ...
        factor, @polyterm_value_3d );

      error = abs ( result - best );

      fprintf ( 1, '  %4d  %8d  %16.8g  %10.2e\n', ...
        factor, node_num, result, error );

      factor = factor * 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = ma_InsertDocking(varargin)
% MA_INSERTDOCKING MATLAB code for ma_InsertDocking.fig
%      MA_INSERTDOCKING, by itself, creates a new MA_INSERTDOCKING or raises the existing
%      singleton*.
%
%      H = MA_INSERTDOCKING returns the handle to a new MA_INSERTDOCKING or the handle to
%      the existing singleton*.
%
%      MA_INSERTDOCKING('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MA_INSERTDOCKING.M with the given input arguments.
%
%      MA_INSERTDOCKING('Property','Value',...) creates a new MA_INSERTDOCKING or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ma_InsertDocking_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ma_InsertDocking_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 ma_InsertDocking

% Last Modified by GUIDE v2.5 20-Jan-2019 16:28:47

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ma_InsertDocking_OpeningFcn, ...
                   'gui_OutputFcn',  @ma_InsertDocking_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 ma_InsertDocking is made visible.
function ma_InsertDocking_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 ma_InsertDocking (see VARARGIN)

% Choose default command line output for ma_InsertDocking
handles.output = hObject;
handles.ma_MainGUI = varargin{1};

% Update handles structure
guidata(hObject, handles);

numOsc = populateOtherSCCombo(handles, handles.dockingOtherScCombo);
if(numOsc == 0)
    warndlg({'There are no other spacecraft available to dock to.', ...
             'To dock, please add an "Other Spacecraft" via the ', ...
             'Spacecraft -> Other Spacecraft menu.'},'No Other S/C','modal');
    close(handles.ma_InsertDocking);
    return;
end

if(length(varargin)>1)
    event = varargin{2};
    setappdata(hObject,'lossConverts',event.massloss.lossConvert);
    populateGUIWithEvent(handles, event);
    set(hObject,'UserData',event);
else
    set(hObject,'UserData',[]);
    setappdata(hObject,'lossConverts',getDefaultLossConvert(handles));
end

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

function populateGUIWithEvent(handles, event)
    set(handles.titleLabel, 'String', 'Edit Docking');
    set(handles.ma_InsertDocking, 'Name', 'Edit Docking');
    set(handles.dockingNameText, 'String', event.name);

    colorStr = getStringFromLineSpecColor(event.lineColor);
    colorValue = findValueFromComboBox(colorStr, handles.coastLineColorCombo);
    set(handles.coastLineColorCombo,'value',colorValue);

    styleStr = getLineStyleFromString(event.lineStyle);
    styleValue = findValueFromComboBox(styleStr, handles.dockingLineStyleCombo);
 	set(handles.dockingLineStyleCombo,'Value',styleValue);

    contents = handles.lineWidthCombo.String;
    contentsDouble = str2double(contents);
    ind = find(contentsDouble == event.lineWidth);
    set(handles.lineWidthCombo,'Value',ind);

    undockTime = event.undockTime;
    set(handles.undockAfterText, 'String', fullAccNum2Str(undockTime));

    set(handles.massLossCheckbox,'Value',event.massloss.use);
    massLossCheckbox_Callback(handles.massLossCheckbox, [], handles);

% --- Outputs from this function are returned to the command line.
function varargout = ma_InsertDocking_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
        maData = getappdata(handles.ma_MainGUI,'ma_data');

        name = get(handles.dockingNameText, 'String');

        contents = cellstr(get(handles.coastLineColorCombo,'String'));
        colorStr = contents{get(handles.coastLineColorCombo,'Value')};
        lineSpecColor = getLineSpecColorFromString(colorStr);

        contents = cellstr(get(handles.dockingLineStyleCombo,'String'));
        lineStyleStr = contents{get(handles.dockingLineStyleCombo,'Value')};
        lineStyle = getLineStyleStrFromText(lineStyleStr);

        contents = handles.lineWidthCombo.String;
        contentsDouble = str2double(contents);
        contensInd = get(handles.lineWidthCombo,'Value');
        lineWidth = contentsDouble(contensInd);

        oScIndex = get(handles.dockingOtherScCombo,'Value');
        otherSC = maData.spacecraft.otherSC{oScIndex};
        oScId = otherSC.id;

        undockTime = str2double(get(handles.undockAfterText,'String'));

        massLoss = struct('use',logical(get(handles.massLossCheckbox,'Value')), 'lossConvert',getappdata(handles.dockingOtherScCombo,'lossConverts'));

        varargout{1} = ma_createDocking(name, oScId, undockTime, lineSpecColor, lineStyle, lineWidth, massLoss);
        close(handles.ma_InsertDocking);
    end

function dockingNameText_Callback(hObject, eventdata, handles)
% hObject    handle to dockingNameText (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 dockingNameText as text
%        str2double(get(hObject,'String')) returns contents of dockingNameText as a double

% --- Executes during object creation, after setting all properties.
function dockingNameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to dockingNameText (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 saveCloseButton.
function saveCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.ma_InsertDocking);
    else
        msgbox(errMsg,'Errors were found while inserting a docking event.','error');
    end

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

function undockAfterText_Callback(hObject, eventdata, handles)
% hObject    handle to undockAfterText (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 undockAfterText as text
%        str2double(get(hObject,'String')) returns contents of undockAfterText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function undockAfterText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to undockAfterText (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 selection change in dockingOtherScCombo.
function dockingOtherScCombo_Callback(hObject, eventdata, handles)
% hObject    handle to dockingOtherScCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns dockingOtherScCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from dockingOtherScCombo

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

% Hint: popupmenu 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 massLossCheckbox.
function massLossCheckbox_Callback(hObject, eventdata, handles)
% hObject    handle to massLossCheckbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of massLossCheckbox
    if(get(hObject,'Value')==1)
        set(handles.editLossConvertButton,'Enable','on');
    else
        set(handles.editLossConvertButton,'Enable','off');
    end

% --- Executes on button press in editLossConvertButton.
function editLossConvertButton_Callback(hObject, eventdata, handles)
% hObject    handle to editLossConvertButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    event = get(handles.ma_InsertDocking,'UserData');
    eventLossConvert = getappdata(handles.ma_InsertDocking,'lossConverts');
    if(~isempty(event) && isstruct(eventLossConvert) && ~isequal(eventLossConvert, getDefaultLossConvert(handles)))
        lossConverts = ma_MassLossesConversionsGUI(handles.ma_MainGUI, eventLossConvert);
    else
        lossConverts = ma_MassLossesConversionsGUI(handles.ma_MainGUI);
    end

    if(~isempty(lossConverts))
        setappdata(handles.ma_InsertDocking,'lossConverts',lossConverts);
    end

    setEditLossConvertButtonTooltipString(handles);

function setEditLossConvertButtonTooltipString(handles)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    names = maData.spacecraft.propellant.names;

	lossConverts = getappdata(handles.ma_InsertDocking,'lossConverts');
    if(~isempty(lossConverts))
        maxNumRes = max(horzcat(lossConverts.resLost, lossConverts.resConvert));
        resRates = zeros(1,maxNumRes);

        for(i=1:length(lossConverts)) %#ok<*NO4LP>
            lossConvert = lossConverts(i);
            resRates(lossConvert.resLost) = resRates(lossConvert.resLost) - lossConvert.resLostRates(lossConvert.resLost);

            sumRates = sum(lossConvert.resLostRates);
            resRates(lossConvert.resConvert) = resRates(lossConvert.resConvert) + sumRates*lossConvert.resConvertPercent(lossConvert.resConvert);
        end

        str = 'Right Click for Copy/Paste\n\nTotal Resource Losses & Conversions\n-----------------------\n';
        for(i=1:length(resRates))
            str = [str, ...
                   sprintf('%s: %1.6f mT/Hr\n', names{i}, resRates(i)*3600)];
        end
    else
        str = 'Right Click for Copy/Paste\n\nTotal Resource Losses & Conversions\n-----------------------\nNone!';
    end

    str = sprintf(str);
    set(handles.editLossConvertButton,'TooltipStr',str);

function errMsg = validateInputs(handles)
    errMsg = {};

    revs = str2double(get(handles.undockAfterText,'String'));
    enteredStr = get(handles.undockAfterText,'String');
    numberName = 'Undock After Time';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(revs, numberName, lb, ub, isInt, errMsg, enteredStr);

% --- Executes on selection change in coastLineColorCombo.
function coastLineColorCombo_Callback(hObject, eventdata, handles)
% hObject    handle to coastLineColorCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns coastLineColorCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from coastLineColorCombo

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

% Hint: popupmenu 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 key release with focus on ma_InsertDocking or any of its controls.
function ma_InsertDocking_WindowKeyReleaseFcn(hObject, eventdata, handles)
% hObject    handle to ma_InsertDocking (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was released, in lower case
%	Character: character interpretation of the key(s) that was released
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) released
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'escape'
            close(handles.ma_InsertDocking);
    end

% --- Executes on selection change in dockingLineStyleCombo.
function dockingLineStyleCombo_Callback(hObject, eventdata, handles)
% hObject    handle to dockingLineStyleCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns dockingLineStyleCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from dockingLineStyleCombo

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

% Hint: popupmenu 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 selection change in lineWidthCombo.
function lineWidthCombo_Callback(hObject, eventdata, handles)
% hObject    handle to lineWidthCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns lineWidthCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from lineWidthCombo

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

% Hint: popupmenu 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

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function order = gl_level_to_order ( dim_num, level )

%*****************************************************************************80
%
%% GL_LEVEL_TO_ORDER converts a GL level to a GL order.
%
%  Discussion:
%
%    Gauss-Legendre rules cannot be nested, but in order to facilitate
%    comparison with Clenshaw Curtis grids, we need to be able to
%    convert a "level" into an order.  Except for the
%    first case of LEVEL = 0, the relationship is
%
%      ORDER = 2**LEVEL + 1
%
%    Level      Order
%
%    0          1
%    1          3 =  2 + 1
%    2          5 =  4 + 1
%    3          9 =  8 + 1
%    4         17 = 16 + 1
%    5         33 = 32 + 1
%
%    In this routine, we assume that a vector of levels is given,!
%    and the corresponding orders are desired.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL(DIM_NUM), the level.
%
%    Output, integer ORDER(DIM_NUM), the order (number of points) of the
%    corresponding rule.
%
  for dim = 1 : dim_num

    if ( level(dim) < 0 )
      order(dim) = -1;
    elseif ( level(dim) == 0 )
      order(dim) = 1;
    else
      order(dim) = ( 2^level(dim) ) + 1;
    end

  end

  return
end


text:
function lerch_values_test ( )

%*****************************************************************************80
%
%% LERCH_VALUES_TEST demonstrates the use of LERCH_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2004
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LERCH_VALUES_TEST:\n' );
  fprintf ( 1, '  LERCH_VALUES stores values of\n' );
  fprintf ( 1, '  the Lerch function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '           Z        S             A            FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, z, s, a, fx ] = lerch_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %4d  %12f  %24.16f\n', z, s, a, fx );

  end

  return
end


text:
function submit(partId)
%SUBMIT Submit your code and output to the ml-class servers
%   SUBMIT() will connect to the ml-class server and submit your solution

  fprintf('==\n== [ml-class] Submitting Solutions | Programming Exercise %s\n==\n', ...
          homework_id());
  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid homework part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  [login password] = loginPrompt();
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to ml-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    % Get Challenge
    [login, ch, signature] = getChallenge(login);
    if isempty(login) || isempty(ch) || isempty(signature)
      % Some error occured, error string in first return element.
      fprintf('\n!! Error: %s\n\n', login);
      return
    end

    % Attempt Submission with Challenge
    ch_resp = challengeResponse(login, password, ch);
    [result, str] = submitSolution(login, ch_resp, partId, output(partId), ...
                                 source(partId), signature);

    fprintf('\n== [ml-class] Submitted Homework %s - Part %d - %s\n', ...
            homework_id(), partId, partNames{partId});
    fprintf('== %s\n', strtrim(str));
    if exist('OCTAVE_VERSION')
      fflush(stdout);
    end
  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '6';
end

function [partNames] = validParts()
  partNames = { 'Gaussian Kernel', ...
                'Parameters (C, sigma) for Dataset 3', ...
                'Email Preprocessing' ...
                'Email Feature Extraction' ...
                };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'gaussianKernel.m' }, ...
           { 'dataset3Params.m' }, ...
           { 'processEmail.m' }, ...
           { 'emailFeatures.m' } };
end

function out = output(partId)
  % Random Test Cases
  x1 = sin(1:10)';
  x2 = cos(1:10)';
  ec = 'the quick brown fox jumped over the lazy dog';
  wi = 1 + abs(round(x1 * 1863));
  wi = [wi ; wi];
  if partId == 1
    sim = gaussianKernel(x1, x2, 2);
    out = sprintf('%0.5f ', sim);
  elseif partId == 2
    load('ex6data3.mat');
    [C, sigma] = dataset3Params(X, y, Xval, yval);
    out = sprintf('%0.5f ', C);
    out = [out sprintf('%0.5f ', sigma)];
  elseif partId == 3
    word_indices = processEmail(ec);
    out = sprintf('%d ', word_indices);
  elseif partId == 4
    x = emailFeatures(wi);
    out = sprintf('%d ', x);
  end
end

function url = challenge_url()
  url = 'http://www.ml-class.org/course/homework/challenge';
end

function url = submit_url()
  url = 'http://www.ml-class.org/course/homework/submit';
end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          while ~feof(fid)
            line = fgets(fid);
            src = [src line];
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature] = getChallenge(email)
  str = urlread(challenge_url(), 'post', {'email_address', email});

  str = strtrim(str);
  [email, str] = strtok (str, '|');
  [ch, str] = strtok (str, '|');
  [signature, str] = strtok (str, '|');
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'homework', homework_id(), ...
            'part', num2str(part), ...
            'email', email, ...
            'output', output, ...
            'source', source, ...
            'challenge_response', ch_resp, ...
            'signature', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Password: ', 's');
end

function [str] = challengeResponse(email, passwd, challenge)
  salt = ')~/|]QMB3[!W`?OVt7qC"@+}';
  str = sha1([challenge sha1([salt email passwd])]);
  sel = randperm(numel(str));
  sel = sort(sel(1:16));
  str = str(sel);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave