text:
function int_val = spline_linear_int ( ndata, tdata, ydata, a, b )

%*****************************************************************************80
%
%% SPLINE_LINEAR_INT evaluates the integral of a piecewise linear spline.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NDATA, the number of data points defining the spline.
%
%    Input, real TDATA(NDATA), YDATA(NDATA), the values of the independent
%    and dependent variables at the data points.  The values of TDATA should
%    be distinct and increasing.
%
%    Input, real A, B, the interval over which the integral is desired.
%
%    Output, real INT_VAL, the value of the integral.
%
  int_val = 0.0;

  if ( a == b )
    return;
  end

  a_copy = min ( a, b );
  b_copy = max ( a, b );
%
%  Find the interval [ TDATA(A_LEFT), TDATA(A_RIGHT) ] that contains, or is
%  nearest to, A.
%
  [ a_left, a_right ] = r8vec_bracket ( ndata, tdata, a_copy );
%
%  Find the interval [ TDATA(B_LEFT), TDATA(B_RIGHT) ] that contains, or is
%  nearest to, B.
%
  [ b_left, b_right ] = r8vec_bracket ( ndata, tdata, b_copy );
%
%  If A and B are in the same interval...
%
  if ( a_left == b_left )

    tval = ( a_copy + b_copy ) / 2.0;

    yp = ( ydata(a_right) - ydata(a_left) ) / ...
         ( tdata(a_right) - tdata(a_left) );

    yval = ydata(a_left) + ( tval - tdata(a_left) ) * yp;

    int_val = yval * ( b_copy - a_copy );

    return;
  end
%
%  Otherwise, integrate from:
%
%  A               to TDATA(A_RIGHT),
%  TDATA(A_RIGHT)  to TDATA(A_RIGHT+1),...
%  TDATA(B_LEFT-1) to TDATA(B_LEFT),
%  TDATA(B_LEFT)   to B.
%
%  Use the fact that the integral of a linear function is the
%  value of the function at the midpoint times the width of the interval.
%
  tval = ( a_copy + tdata(a_right) ) / 2.0E+00;

  yp = ( ydata(a_right) - ydata(a_left) ) / ...
       ( tdata(a_right) - tdata(a_left) );

  yval = ydata(a_left) + ( tval - tdata(a_left) ) * yp;

  int_val = int_val + yval * ( tdata(a_right) - a_copy );

  for i_left = a_right : b_left - 1

    tval = ( tdata(i_left+1) + tdata(i_left) ) / 2.0E+00;

    yp = ( ydata(i_left+1) - ydata(i_left) ) / ...
         ( tdata(i_left+1) - tdata(i_left) );

    yval = ydata(i_left) + ( tval - tdata(i_left) ) * yp;

    int_val = int_val + yval * ( tdata(i_left + 1) - tdata(i_left) );

  end

  tval = ( tdata(b_left) + b_copy ) / 2.0E+00;

  yp = ( ydata(b_right) - ydata(b_left) ) / ...
       ( tdata(b_right) - tdata(b_left) );

  yval = ydata(b_left) + ( tval - tdata(b_left) ) * yp;

  int_val = int_val + yval * ( b_copy - tdata(b_left) );

  if ( b < a )
    int_val = -int_val;
  end

  return
end


text:
classdef nme_branch_dc < mp.nme_branch & mp.form_dc

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

%     properties
%         name = 'branch';
%     end

    methods
        function obj = build_params(obj, nm, dm)
            build_params@mp.nme_branch(obj, nm, dm);    %% call parent

            dme = obj.data_model_element(dm);
            nl = obj.nk;

            tm = ones(nl, 1);           %% default tap ratio = 1
            i = find(dme.tm);           %% indices of non-zero tap ratios
            tm(i) = dme.tm(i);          %% assign non-zero tap ratios

            b = 1 ./ dme.x;             %% series susceptance
            b = b ./ tm;
            Pfinj = b .* (-dme.ta);
            obj.B = sparse( ...
                [1:nl 1:nl nl+1:2*nl nl+1:2*nl]', ...
                [1:nl nl+1:2*nl 1:nl nl+1:2*nl]', ...
                [b; -b; -b; b], ...
                2*nl, 2*nl );
            obj.p = [Pfinj + dme.g_fr; -Pfinj + dme.g_to];
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function nperm = perm_enum ( n )

%*****************************************************************************80
%
%% PERM_ENUM enumerates the permutations on N digits.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of values being permuted.
%    N must be nonnegative.
%
%    Output, integer NPERM, the number of distinct elements.
%
  nperm = i4_factorial ( n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function kbrDisplay(model, spacing)

% KBRDISPLAY Display parameters of the KBR model.
% FORMAT
% DESC displays the parameters of the kernel based regression
% model and the model type to the console.
% ARG model : the model to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG model : the model to display.
% ARG spacing : how many spaces to indent the display of the model by.
%
% SEEALSO : kbrCreate, modelDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2007

% MLTOOLS

if nargin > 1
  spacing = repmat(32, 1, spacing);
else
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('Kernel based regression model:\n')
fprintf(spacing);
fprintf('Kernel type:\n')
kernDisplay(model.kern, length(spacing)+2)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ ar ] = ctransform(a)
% Copula-transform array - rank and scale to [0, 1]
    [as ai] = sort(a, 2);
    [aa ar] = sort(ai, 2);
    ar = (ar - 1) / (size(ar, 2) - 1);
end

author: Sable

repo: mcbench-benchmarks


text:
% run_pressproc_GEauto.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,out_w,out_noproc,out_w_noproc]=run_pressproc_GEauto(filestring,aaDomain,tmaxin,iterin);
%
% DESCRIPTION:
% Fully automated processing script for GE PRESS MRS data in p-file
% format (GE raw data).  Includes combination of reciever channels,
% removal of bad averages, freqeuncy drift correction, and leftshifting.
% This pipeline requires no user interaction. This function automatically
% generates an html report to describe the results of each processing step.
%
% INPUTS:
% filestring    = String variable for the name of the p-file
%                   the water suppressed data;
% aaDomain      = (Optional) Perform the spectral registration (drift correction) using
%                   the full spectrum ('t'), or only a limited frequency
%                   range ('f').  Default is 'f'.
% tmaxin        = (Optional).  Duration (in sec.) of the time domain signal
%                   used in the spectral registration (drift correction).
%                   Default is 0.2 sec.
% iterin        = (Optional).  Maximum number of allowed iterations for the spectral
%                   registration to converge. Default is 20.
%
% OUTPUTS:
% out           = Fully processed, water suppressed output spectrum.
% out_w         = Fully processed, water unsuppressed output spectrum.
% out_noproc    = Water suppressed output spectrum without pre-
%                   processing (No bad-averages removal, no frequency drift
%                   correction).
% out_w_noproc  = Water unsuppressed output spectrum without pre-
%                   processing.

function [out,outw,out_noproc,outw_noproc]=run_pressproc_GEauto(filestring,aaDomain,tmaxin,iterin);

if nargin<4
    iterin=20;
    if nargin<3
        tmaxin=0.2;
        if nargin<2
            aaDomain='f';
        end
    end
end

%make a new directory for the output report and figures:
mkdir(['./report']);
mkdir(['./report/figs']);

% %read in both datasets:
[raw,raww]=io_loadspec_GE(filestring,1);
raw=op_complexConj(raw);

if isstruct(raww)
    water=true;
    raww=op_complexConj(raww);
end

%first step should be to combine coil channels.  To do this find the coil
%phases from the water unsuppressed data.
if water
    coilcombos=op_getcoilcombos(raww,1);
    [outw_cc,fidw_pre,specw_pre,phw,sigw]=op_addrcvrs(raww,1,'w',coilcombos);
else
    coilcombos=op_getcoilcombos(op_averaging(raw),1);
end
[out_cc,fid_pre,spec_pre,ph,sig]=op_addrcvrs(raw,1,'w',coilcombos);
[out_av_cc,fid_av_pre,spec_av_pre]=op_addrcvrs(op_averaging(raw),1,'w',coilcombos);
raw_av=op_averaging(raw);

%generate unprocessed spectrum:
out_noproc=op_averaging(out_cc);
if water
    outw_noproc=op_averaging(outw_cc);
end

%Generate plots showing coil channels before and after phase alignment
%figure('position',[0 50 560 420]);
h=figure('visible','off');
subplot(1,2,1);
plot(raw_av.ppm,real(raw_av.specs(:,:,1)));xlim([1 5]);
set(gca,'FontSize',8);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude (a.u.)','FontSize',10);
title('Before correction','FontSize',12);
box off;
subplot(1,2,2);
plot(raw_av.ppm,real(spec_av_pre(:,:,1)));xlim([1 5]);
set(gca,'FontSize',12);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude(a.u.)','FontSize',10);
title('After correction','FontSize',12);
box off;
set(h,'PaperUnits','centimeters');
set(h,'PaperPosition',[0 0 20 10]);
saveas(h,['./report/figs/coilReconFig'],'jpg');
saveas(h,['./report/figs/coilReconFig'],'fig');
close(h);

%%%%%%%%OPTIONAL REMOVAL OF BAD AVERAGES FROM FIRST DATASET%%%%%%%%%%%%%%%%%%%%
close all;
out_cc2=out_cc;
nBadAvgTotal=0;
nbadAverages=1;
rmbadav='y';
close all;
if rmbadav=='n' || rmbadav=='N'
    out_rm=out_cc;
    nsd='N/A';
else
    sat='n'
    while sat=='n' || sat=='N'
        nsd=4; %Setting the number of standard deviations;
        iter=1;
        nbadAverages=1;
        nBadAvgTotal=0;
        out_cc2=out_cc;
        while nbadAverages>0
            [out_rm,metric{iter},badAverages]=op_rmbadaverages(out_cc2,nsd,'t');
            badAverages;
            nbadAverages=length(badAverages);
            nBadAvgTotal=nBadAvgTotal+nbadAverages;
            out_cc2=out_rm;
            iter=iter+1;
            disp([num2str(nbadAverages) ' bad averages removed on this iteration.']);
            disp([num2str(nBadAvgTotal) ' bad averages removed in total.']);
            close all;
        end
        %figure('position',[0 50 560 420]);
        %Make figure to show pre-post removal of averages
        h=figure('visible','off');
        subplot(1,2,1);
        plot(out_cc.ppm,real(out_cc.specs(:,:)));xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Before','FontSize',12);
        box off;
        subplot(1,2,2);
        plot(out_rm.ppm,real(out_rm.specs(:,:)));xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('After','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 15]);
        saveas(h,['./report/figs/rmBadAvg_prePostFig'],'jpg');
        saveas(h,['./report/figs/rmBadAvg_prePostFig'],'fig');
        close(h);

        %figure('position',[0 550 560 420]);
        h=figure('visible','off');
        plot([1:length(metric{1})],metric{1},'.r',[1:length(metric{iter-1})],metric{iter-1},'x','MarkerSize',16);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Deviation Metric','FontSize',10);
        legend('Before rmBadAv','After rmBadAv');
        legend boxoff;
        title('Deviation Metric','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 10]);
        saveas(h,['./report/figs/rmBadAvg_scatterFig'],'jpg');
        saveas(h,['./report/figs/rmBadAvg_scatterFig'],'fig');
        close(h);

        %sat1=input('are you satisfied with the removal of bad averages? ','s');
        sat='y';

    end
end

%NOW ALIGN AVERAGES:  A.K.A. Frequency Drift Correction.
driftCorr='y';
if driftCorr=='n' || driftCorr=='N'
    out_av=op_averaging(out_rm);
    if water
        outw_av=op_averaging(outw_cc);
    end
    fs=0;
    phs=0;
else
    if water
        outw_aa=op_alignAverages(outw_cc,0.2,'n');
    end
    sat='n';
    out_rm2=out_rm;
    while sat=='n' || sat=='N'
        fsPoly=100;
        phsPoly=1000;
        fscum=zeros(out_rm2.sz(out_rm2.dims.averages),1);
        phscum=zeros(out_rm2.sz(out_rm2.dims.averages),1);
        iter=1;
        while (abs(fsPoly(1))>0.001 || abs(phsPoly(1))>0.01) && iter<iterin
            iter=iter+1
            close all
            tmax=0.25+0.03*randn(1);
            ppmmin=1.6+0.1*randn(1);
            ppmmaxarray=[3.5+0.1*randn(1,2),4+0.1*randn(1,3),5.5+0.1*randn(1,1)];
            ppmmax=ppmmaxarray(randi(6,1));
            switch aaDomain
                case 't'
                    [out_aa,fs,phs]=op_alignAverages(out_rm2,tmax,'y');
                case 'f'
                    [out_aa,fs,phs]=op_alignAverages_fd(out_rm2,ppmmin,ppmmax,tmax,'y');
                otherwise
                    error('ERROR: avgAlignDomain not recognized!');
            end

            fsPoly=polyfit([1:out_aa.sz(out_aa.dims.averages)]',fs,1)
            phsPoly=polyfit([1:out_aa.sz(out_aa.dims.averages)]',phs,1)
            iter

            fscum=fscum+fs;
            phscum=phscum+phs;

            if driftCorr=='y' || driftCorr=='Y'
                out_rm2=out_aa;
            end
        end
        h=figure('visible','off');
        subplot(1,2,1);
        plot(out_rm.ppm,real(out_rm.specs(:,:)));xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('Before','FontSize',12);
        box off;
        subplot(1,2,2);
        plot(out_aa.ppm,real(out_aa.specs(:,:)));xlim([1 5]);
        set(gca,'FontSize',8);
        set(gca,'XDir','reverse');
        xlabel('Frequency (ppm)','FontSize',10);
        ylabel('Amplitude(a.u.)','FontSize',10);
        title('After','FontSize',12);
        box off;
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 20 15]);
        saveas(h,['./report/figs/alignAvgs_prePostFig'],'jpg');
        saveas(h,['./report/figs/alignAvgs_prePostFig'],'fig');
        close(h);

        h=figure('visible','off');
        plot([1:out_aa.sz(out_aa.dims.averages)],fscum,'.-','LineWidth',2);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Frequency Drift [Hz]','FontSize',10);
        box off;
        legend('Frequency Drift','Location','SouthEast');
        legend boxoff;
        title('Estimated Freqeuncy Drift','FontSize',12);
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 10 10]);
        saveas(h,['./report/figs/freqDriftFig'],'jpg');
        saveas(h,['./report/figs/freqDriftFig'],'fig');
        close(h);

        h=figure('visible','off');
        plot([1:out_aa.sz(out_aa.dims.averages)],phscum,'.-','LineWidth',2);
        set(gca,'FontSize',8);
        xlabel('Scan Number','FontSize',10);
        ylabel('Phase Drift [Deg.]','FontSize',10);
        box off;
        legend('Phase Drift','Location','SouthEast');
        legend boxoff;
        title('Estimated Phase Drift','FontSize',12);
        set(h,'PaperUnits','centimeters');
        set(h,'PaperPosition',[0 0 10 10]);
        saveas(h,['./report/figs/phaseDriftFig'],'jpg');
        saveas(h,['./report/figs/phaseDriftFig'],'fig');
        close(h);

        sat='y';
        if sat=='n'
            iter=0;
            p1=100;
            fscum=zeros(out_rm.sz(2:end));
            phscum=zeros(out_rm.sz(2:end));
            fs2cum=zeros(out_cc.sz(2:end));
            phs2cum=zeros(out_cc.sz(2:end));
            out_rm2=out_rm;
            out_cc2=out_cc;
        end
        totalFreqDrift=mean(max(fscum)-min(fscum));
        totalPhaseDrift=mean(max(phscum)-min(phscum));
        close all
    end
    %now combine the averages averages
    out_av=op_averaging(out_aa);
    if water
        outw_av=op_averaging(outw_aa);
    end
end

%now leftshift
out_ls=op_leftshift(out_av,out_av.pointsToLeftshift);
if water
    outw_ls=op_leftshift(outw_av,outw_av.pointsToLeftshift);
end

%now do automatic zero-order phase correction (Use water Peak):
out_ls_zp=op_zeropad(out_ls,16);
%index=find(abs(out_ls_zp.specs)==max(abs(out_ls_zp.specs(out_ls_zp.ppm>2.9 & out_ls_zp.ppm<3.1))));
%ph0=-phase(out_ls_zp.specs(index))*180/pi;
[out_ls_zp_ph,ph0]=op_autophase(out_ls_zp,4,5.5);
out_ph=op_addphase(out_ls,ph0);
%And now for water unsuppressed data:
if water
    outw_ls_zp=op_zeropad(outw_ls,16);
    %indexw=find(abs(outw_ls_zp.specs)==max(abs(outw_ls_zp.specs(outw_ls_zp.ppm>4 & outw_ls_zp.ppm<5.5))));
    %ph0w=-phase(outw_ls_zp.specs(indexw))*180/pi;
    outw_ph=op_addphase(outw_ls,ph0);
    outw_ls_zp_ph=op_addphase(outw_ls_zp,ph0);
end

%do same phase corection on unprocessed data
out_noproc=op_addphase(op_leftshift(out_noproc,out_noproc.pointsToLeftshift),ph0);
if water
    outw_noproc=op_addphase(op_leftshift(outw_noproc,outw_noproc.pointsToLeftshift),ph0);
end

%Frequency shift all spectra so that Creatine appears at 3.027 ppm:
[~,frqShift]=op_ppmref(out_ls_zp_ph,2.9,3.1,3.027);
out=op_freqshift(out_ph,frqShift);
out_noproc=op_freqshift(out_noproc,frqShift);
%And now for water unsuppressed data (user water peak and set to 4.65 ppm):
if water
    [~,frqShiftw]=op_ppmref(outw_ls_zp_ph,4,5.5,4.65);
    outw=op_freqshift(outw_ph,frqShiftw);
    outw_noproc=op_freqshift(outw_noproc,frqShiftw);
end

%Make figure to show the final spectrum:
h=figure('visible','off');
plot(out.ppm,out.specs,'linewidth',2);xlim([0.2 5.2]);
set(gca,'FontSize',8);
set(gca,'XDir','reverse');
xlabel('Frequency (ppm)','FontSize',10);
ylabel('Amplitude(a.u.)','FontSize',10);
legend('diff');
legend boxoff;
box off;
title('Result: Final Spectrum','FontSize',12);
set(h,'PaperUnits','centimeters');
set(h,'PaperPosition',[0 0 20 10]);
saveas(h,['./report/figs/finalSpecFig'],'jpg');
saveas(h,['./report/figs/finalSpecFig'],'fig');

% writ=input('Write? (y or n):  ','s');
writ='y';
if writ=='y' || writ=='Y'
     RF=io_writelcm(out,['./' filestring '_lcm'],out.te);
    if water
        RF=io_writelcm(outw,['./' filestring '_w_lcm'],outw.te);
    end
end

close all;

%write an html report:
fid=fopen(['./report/report.html'],'w+');
fprintf(fid,'<!DOCTYPE html>');
fprintf(fid,'\n');
logoPath=which('FID-A_LOGO.jpg');
fprintf(fid,'\n',logoPath);
fprintf(fid,'\n
FID-A Processing Report
');
fprintf(fid,'\n
Processing pipeline applied to PRESS data using run_pressproc_auto.m
');
fprintf(fid,'\n
FILENAME: %s/%s
',pwd,filestring);
fprintf(fid,'\n
DATE: %s
',date);
fprintf(fid,'\n\n

');
fprintf(fid,'\n\n
Results of multi-coil combination:
');
fprintf(fid,'\n',pwd);
fprintf(fid,'\n\n

');
fprintf(fid,'\n\n
Results of removal of bad averages:
');
fprintf(fid,'\n
Original number of averages: \t%5.6f
',raw.sz(raw.dims.averages));
fprintf(fid,'\n
Number of bad Averages removed:  \t%5.6f
',nBadAvgTotal);
fprintf(fid,'\n
Number of remaining averages in processed dataset:  \t%5.6f
',out_rm.sz(out_rm.dims.averages));
fprintf(fid,'\n
Bad Averages Removal Threshold was:  \t%2.2f
',nsd);
fprintf(fid,'\n',pwd,pwd);
fprintf(fid,'\n\n

');
fprintf(fid,'\n\n

');
fprintf(fid,'\n\n
Results of spectral registration:
');
fprintf(fid,'\n
Total frequency drift was: \t%5.6f
',max(totalFreqDrift));
fprintf(fid,'\n
Total phase drift was: \t%5.6f
',max(totalPhaseDrift));
fprintf(fid,'\n',pwd);
fprintf(fid,'\n\n

');
fprintf(fid,'\n',pwd,pwd);
fprintf(fid,'\n\n

');
fprintf(fid,'\n\n
Final Result:
');
fprintf(fid,'\n',pwd);
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [tm,D2]=momftfr(tfr,tmin,tmax,time);
%MOMFTFR Frequency moments of a time-frequency representation.
%	[TM,D2]=MOMFTFR(TFR,TMIN,TMAX,TIME) computes the frequeny
%	moments of a time-frequency representation.
%
%	TFR    : time-frequency representation ([Nrow,Ncol]size(TFR)).
%	TMIN   : smallest column element of TFR taken into account
%	                            (default : 1)
%	TMAX   : highest column element of TFR taken into account
%	                            (default : Ncol)
%	TIME   : true time instants (default : 1:Ncol)
%	TM     : averaged time          (first order moment)
%	D2     : squared time duration  (second order moment)
%
%	Example :
%	 sig=fmlin(128,0.1,0.4);
%	 [tfr,t,f]=tfrwv(sig); [tm,D2]=momftfr(tfr);
%	 subplot(211); plot(f,tm); subplot(212); plot(f,D2);
%
%	See also MOMTTFR, MARGTFR.

%	F. Auger, August 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

[tfrrow,tfrcol]=size(tfr);
if (nargin==1),
 tmin=1; tmax=tfrcol; time=tmin:tmax;
elseif (nargin==2),
 tmax=tfrcol; time=tmin:tmax;
elseif (nargin==3),
 time=tmin:tmax;
end;

if (tmin>tmax)|(tmin<=0)|(tmax>tfrcol),
 error('1<=TMIN<=TMAX<=Ncol');
end;

E  = sum(tfr(:,tmin:tmax).');
tm = (time    * tfr(:,tmin:tmax).' ./E).';
D2 = (time.^2 * tfr(:,tmin:tmax).' ./E).' - tm.^2;

author: HeLiangHIT


text:
function r8plu_mul_test ( )

%*****************************************************************************80
%
%% R8PLU_MUL_TEST tests R8PLU_MUL;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8PLU_MUL_TEST\n' );
  fprintf ( 1, '  R8PLU_MUL computes the product A*x=b\n' );
  fprintf ( 1, '  using the compressed PLU factors of A.\n' );
  fprintf ( 1, '  Using initial random number seed = %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8mat_uniform_01 ( n, n, seed );

  r8mat_print ( n, n, a, '  The matrix A:' );
%
%  Set the right hand side B1.
%
  for i = 1 : n
    x(i) = i;
  end

  b(1:n) = a(1:n,1:n) * x(1:n)';

  r8vec_print ( n, b, '  The right hand side B (computed from A):' );
%
%  Factor the matrix.
%
  [ pivot, lu, info ] = r8mat_to_r8plu ( n, a );
%
%  Compute the matrix-vector product.
%
  b = r8plu_mul ( n, pivot, lu, x );

  r8vec_print ( n, b, '  The right hand side B (computed from PLU):' );

  return
end


text:
function [d,fp,dt,tc,t]=readhtk(file)
%READHTK  read an HTK parameter file [D,FP,DT,TC,T]=(FILE)
%
% Input:
%    FILE = name of HTX file
% Outputs:
%       D = data: column vector for waveforms, one row per frame for other types
%      FP = frame period in seconds
%      DT = data type (also includes Voicebox code for generating data)
%             0  WAVEFORM     Acoustic waveform
%             1  LPC          Linear prediction coefficients
%             2  LPREFC       LPC Reflection coefficients:  -lpcar2rf([1 LPC]);LPREFC(1)=[];
%             3  LPCEPSTRA    LPC Cepstral coefficients
%             4  LPDELCEP     LPC cepstral+delta coefficients (obsolete)
%             5  IREFC        LPC Reflection coefficients (16 bit fixed point)
%             6  MFCC         Mel frequency cepstral coefficients
%             7  FBANK        Log Fliter bank energies
%             8  MELSPEC      linear Mel-scaled spectrum
%             9  USER         User defined features
%            10  DISCRETE     Vector quantised codebook
%            11  PLP          Perceptual Linear prediction
%            12  ANON
%      TC = full type code = DT plus (optionally) one or more of the following modifiers
%               64  _E  Includes energy terms
%              128  _N  Suppress absolute energy
%              256  _D  Include delta coefs
%              512  _A  Include acceleration coefs
%             1024  _C  Compressed
%             2048  _Z  Zero mean static coefs
%             4096  _K  CRC checksum (not implemented yet)
%             8192  _0  Include 0'th cepstral coef
%            16384  _V  Attach VQ index
%            32768  _T  Attach delta-delta-delta index
%       T = text version of type code e.g. LPC_C_K

%   Thanks to Dan Ellis (ee.columbia.edu) for sorting out decompression.
%   Thanks to Stuart Anderson (whispersys.com) for making it work on 64 bit machines.

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

fid=fopen(file,'r','b');
if fid < 0
    error(sprintf('Cannot read from file %s',file));
end
nf=fread(fid,1,'int32');             % number of frames
fp=fread(fid,1,'int32')*1.E-7;       % frame interval (converted to seconds)
by=fread(fid,1,'int16');            % bytes per frame
tc=fread(fid,1,'int16');            % type code (see comments above for interpretation)
tc=tc+65536*(tc<0);
cc='ENDACZK0VT';                    % list of suffix codes
nhb=length(cc);                     % number of suffix codes
ndt=6;                              % number of bits for base type
hb=floor(tc*pow2(-(ndt+nhb):-ndt));
hd=hb(nhb+1:-1:2)-2*hb(nhb:-1:1);   % extract bits from type code
dt=tc-pow2(hb(end),ndt);            % low six bits of tc represent data type

% hd(7)=1 CRC check
% hd(5)=1 compressed data
if (dt==5)  % hack to fix error in IREFC files which are sometimes stored as compressed LPREFC
    fseek(fid,0,'eof');
    flen=ftell(fid);        % find length of file
    fseek(fid,12,'bof');
    if flen>14+by*nf        % if file is too long (including possible CRCC) then assume compression constants exist
        dt=2;               % change type to LPREFC
        hd(5)=1;            % set compressed flag
        nf=nf+4;            % frame count doesn't include compression constants in this case
    end
end

if any(dt==[0,5,10])        % 16 bit data for waveforms, IREFC and DISCRETE
    d=fread(fid,[by/2,nf],'int16').';
    if ( dt == 5),
        d=d/32767;                    % scale IREFC
    end
else
    if hd(5)                            % compressed data - first read scales
        nf = nf - 4;                    % frame count includes compression constants
        ncol = by / 2;
        scales = fread(fid, ncol, 'float');
        biases = fread(fid, ncol, 'float');
        d = ((fread(fid,[ncol, nf], 'int16')+repmat(biases,1,nf)).*repmat(1./scales,1,nf)).';
    else                              % uncompressed data
        d=fread(fid,[by/4,nf],'float').';
    end
end;
fclose(fid);
if nargout > 4
    ns=sum(hd);                 % number of suffixes
    kinds={'WAVEFORM' 'LPC' 'LPREFC' 'LPCEPSTRA' 'LPDELCEP' 'IREFC' 'MFCC' 'FBANK' 'MELSPEC' 'USER' 'DISCRETE' 'PLP' 'ANON' '???'};
    t=[kinds{min(dt+1,length(kinds))} reshape(['_'*ones(1,ns);cc(hd>0)],1,2*ns)];
end

author: jtkim-kaist

repo: VAD

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


text:
function d=quadratic_form_distance(XI,XJ,A)
  % Implementation of the Quadratic form distance
  % (cf. "The Earth Movers' Distance as a Metric for Image Retrieval",
  %      Y. Rubner, C. Tomasi, L.J. Guibas, 2000)
  %
  % @author: B. Schauerte
  % @date:   2009
  % @url:    http://cvhci.anthropomatik.kit.edu/~bschauer/

  % Copyright 2009 B. Schauerte. 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 B. SCHAUERTE ''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 B. SCHAUERTE 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 B. Schauerte.

	if nargin < 3
		error('Not enough arguments, i.e. similarity matrix is missing')
	end

  m=size(XJ,1); % number of samples of p
  p=size(XI,2); % dimension of samples

  assert(p == size(XJ,2)); % equal dimensions
  assert(size(XI,1) == 1); % pdist requires XI to be a single sample

  d=zeros(m,1); % initialize output array

	diff=zeros(1,p);
  for i=1:m
		diff=(XI - XJ(i,:));
		d(i) = sqrt( diff * A * diff' );
  end

author: opencobra

repo: cobratoolbox


text:
function [ fxx, fxy, fyy ] = f03_f2 ( n, x, y )

%*****************************************************************************80
%
%% F03_F2 returns second derivatives of function 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real FXX(N,1), FXY(N,1), FYY(N,1), second derivatives.
%
  t1(1:n,1) = 5.4 * y(1:n,1);
  t2(1:n,1) = 1.0 + ( 3.0 * x(1:n,1) - 1.0 ).^2;

  fxx(1:n,1) = 3.0 * ( 1.25 + cos ( t1(1:n,1) ) ) .* ( 3.0 * t2(1:n,1) - 4.0 ) ...
    ./ ( t2(1:n,1).^3 );
  fxy(1:n,1) = 5.4 * ( 3.0 * x(1:n,1) - 1.0 ) .* sin ( t1(1:n,1) ) ...
    ./ ( t2(1:n,1) .* t2(1:n,1) );
  fyy(1:n,1) = - 4.86 * cos ( t1(1:n,1) ) ./ t2(1:n,1);

  return
end


text:
function labels_plot(expdirs,jabs,T0,T1,varargin)
% Labelgram
% labels_plot(expdirs,jabs,T0,T1,varargin)

if ~exist('expdirs','var')
  expdirs = {};
end
assert(iscellstr(expdirs));
if ~exist('jabs','var')
  jabs = [];
end
if isempty(jabs)
  [tfsuccess,jabs] = ExpPP.uiGetJabFiles;
  if ~tfsuccess
    return;
  end
end
validateattributes(T0,{'numeric'},{'scalar' 'integer'});
validateattributes(T1,{'numeric'},{'scalar' 'integer'});
assert(T1>=T0,'T1 must be greater than or equal to T0.');
% if ~exist('nframesplot','var')
%   nframesplot = [];
% elseif ~isempty(nframesplot)
%   assert(isnumeric(nframesplot));
%   assert(isscalar(nframesplot) || numel(nframesplot)==2,'Invalid nframesplot');
% end

assert(isscalar(jabs),'Currently expect single jab/macguffin.');
if iscellstr(jabs)
  jabfiles = jabs;
  jabs = loadAnonymous(jabs{1});
elseif isa(jabs,'Macguffin');
  jabfiles = {[]};
else
  assert(false,'''jabs'' argument must be either a jab filename or object.');
end

if isempty(expdirs)
  expdirs = jabs.expDirNames;
end
Nexp = numel(expdirs);

jab = jabs;
%jab.modernize();

% get the labels for expdirs
assert(numel(jab.expDirNames)==numel(jab.labels));
[tf,loc] = ismember(expdirs,jab.expDirNames);
lbls = Labels.labels(Nexp);
lbls(tf) = jab.labels(loc(tf));

% AL: loading Trx takes forever
% % Trx/T0/T1/nframesplot
% trxfname = jab.file.trxfilename;
% trxstarts = nan(Nexp,1);
% trxends = nan(Nexp,1);
% for i = 1:Nexp
%   trxfile = fullfile(expdirs{i},trxfname);
%   tmp = load(trxfile,'trx');
%   trxstarts(i) = min([tmp.trx.firstframe]);
%   trxends(i) = max([tmp.trx.endframe]);
% end
% T0 = min(trxstarts);
% T1 = max(trxends)-1;
% if ~isempty(nframesplot)
%   if numel(nframesplot)==1,
%     T1 = min(T1,T0+nframesplot-1);
%   else
%     assert(numel(nframesplot)==2);
%     T0 = nframesplot(1);
%     T1 = min(T1,nframesplot(2));
%   end
% end

% go
behnames = Labels.verifyBehaviorNames(jab.behaviors.names);
m = Labels.labelMatrix(lbls,T0,T1,behnames);
ethogram_plot_core(figure,m,T0,T1,expdirs,behnames,repmat(jabfiles,numel(behnames),1));

author: kristinbranson

repo: JAABA


text:
function blend_test05 ( )

%*****************************************************************************80
%
%% BLEND_TEST05 checks out BLEND_IJK_0D1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
  m1 = 4;
  m2 = 3;
  m3 = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLEND_TEST05\n' );
  fprintf ( 1, '  BLEND_IJK_0D1 interpolates data in a table,\n' );
  fprintf ( 1, '  from corner data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The table is %d rows by %d columns by %d layers.\n', ...
    m1, m2, m3 );
%
%  Load data on the faces.
%
  for i = 1 : m1
    r = ( i - 1 ) / ( m1 - 1 );
    for j = 1 : m2
      s = ( j - 1 ) / ( m2 - 1 );
      for k = 1 : m3
        t = ( k - 1 ) / ( m3 - 1 );

        num_extreme = 0;
        if ( i == 1 | i == m1 )
          num_extreme = num_extreme + 1;
        end
        if ( j == 1 | j == m2 )
          num_extreme = num_extreme + 1;
        end
        if ( k == 1 | k == m3 )
          num_extreme = num_extreme + 1;
        end

        if ( num_extreme == 3 )
          x(i,j,k) = quad_rst ( r, s, t, 1 );
        else
          x(i,j,k) = 0.0;
        end

      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data given to BLEND_IJK_0D1:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : m3
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Layer K = %d\n', k );
    fprintf ( 1, '\n' );
    for i = 1 : m1
      fprintf ( 1, '  %10f  %10f  %10f\n', x(i,1:m2,k) );
    end
  end

  x = blend_ijk_0d1 ( x, m1, m2, m3 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Values interpolated by BLEND_IJK_0D1:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : m3
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Layer K = %d\n', k );
    fprintf ( 1, '\n' );
    for i = 1 : m1
      fprintf ( 1, '  %10f  %10f  %10f\n', x(i,1:m2,k) );
    end
  end
%
%  Load all data.
%
  for i = 1 : m1
    r = ( i - 1 ) / ( m1 - 1 );
    for j = 1 : m2
      s = ( j - 1 ) / ( m2 - 1 );
      for k = 1 : m3
        t = ( k - 1 ) / ( m3 - 1 );
        x(i,j,k) = quad_rst ( r, s, t, 1 );
      end
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact data:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : m3
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Layer K = %d\n', k );
    fprintf ( 1, '\n' );
    for i = 1 : m1
      fprintf ( 1, '  %10f  %10f  %10f\n', x(i,1:m2,k) );
    end
  end

  return
end


text:
function matlab2hypreParVectors(matlab_mvector,num_procs,vector_filename,varargin)
%MATLAB2HYPREPARVECTORS Converts MATLAB Vectors to HYPRE IJ Format.
%
%   This program was tested on Matlab 7.1.
%
%   Description:
%     matlab2hypreParVectors(matlab_mvector,num_procs) converts MATLAB
%     formatted vectors to HYPRE formatted ones, and writes the files to
%     the disk. This function takes as an input an NxM matlab matrix
%     matlab_mvector, where M represents the number of vectors and N stands
%     for the size of each vector and writes the matrix to the number of
%     files specified by giving num_procs argument
%
%     matlab2hypreParVectors(matlab_mvector,num_procs,vector_filename)
%     writes the HYPRE formatted vectors to the output file named
%     vector_filename e.g. if vector_filename is 'vectors' the program
%     will generate first 4 files and the output files will be:
%       vectors.0.0, vectors.0.1, vectors.0.INFO.0, vectors.0.INFO.1.
%     the Mth output will be, consequently,
%       vectors.m.0
%       vectors.m.1
%       vectors.m.INFO.0
%       vectors.m.INFO.1.
%
%     matlab2hypreParVectors(matlab_mvector,num_procs,vector_filename,
%       precision) converts from MATLAB to HYPRE with specified precision
%       given by floating point precision.
%
%     Examples:
%       matlab2hypreParVectors(A,2,'vectors')
%       1) Create a matrix A in matlab, consisting of M vectors of size
%         N to be converted.
%       2)Set num_procs to 2. (Specifies the number of processors)
%       3)'vectors' specifies the hypre output filenames:
%         vectors.0.0
%         vectors.0.1
%         vectors.0.INFO.0
%         vectors.0.INFO.1
%
%       matlab2hypreParVectors(A,2,'vectors','15.10e')
%       changes format of floating point numbers in file
%
%   See also matlabIJ2hypre, testIJmatlabhypre.m,
%      hypreIJ2matlab.m, hypreParVectors2matlab.m
%
%   Author: Diana Zakaryan, Dept of Mathematics, University of Colorado,
%      Denver, 15-Mar-2005.
%

for j=1:size(matlab_mvector,2)
    % create filename
    filename = strcat(vector_filename, '.', num2str(j-1));
    B=matlab_mvector(:,j);

    % call single vector convertor
    matlab2hypreParVector(B,num_procs,filename,varargin);
end

function matlab2hypreParVector(matlab_mvector,num_procs,filename,varargin)

% check
if isvector(matlab_mvector)~=1
    error('The argument must be a vector.')
end

% default print format
my_format='20.19e';
if (nargin>4)
    my_format=varargin(1);
end
prt_format=strcat('%',char(my_format));

% intialization
n=size(matlab_mvector,1);

% generate partitioning
part=generate_part(n,num_procs);

% generate Hypre input files
for i=1:num_procs
    s=int2str(i-1);
    filename2=strcat(filename,'.',s);
    fprintf('Generating file: %s\n',filename2);
    X=matlab_mvector(part(i)+1:part(i+1));
    nrows=part(i+1)-part(i);
    dlmwrite(filename2, nrows, 'precision', '%d');
    dlmwrite(filename2, X, '-append','precision',prt_format);

    % writing INFO file
    filename2=strcat(filename,'.','INFO','.',s);
    fprintf('Generating INFO file: %s\n',filename2);
    Y = [n;0;nrows];
    dlmwrite(filename2, Y, 'precision', '%d');
    clear ('X', 'Y');
end
return

function [partitioning]=generate_part(length,num_procs)
% [partitioning]=generate_part(length,num_procs)
% generate partitioning across processes
% See Hypre getpart.c
partitioning=zeros(num_procs+1,1);
size = floor(length/num_procs);
rest= length - size*num_procs;
for i=1:num_procs-1
    partitioning(i+1) = partitioning(i)+size;
    if (i==1)
        partitioning(i+1)=partitioning(i+1)+rest;
    end
end
partitioning(num_procs+1)=length;
return


text:
% Authors: A. Iscen, G. Tolias, Y. Avrithis, T. Furon, O. Chum. 2017.
% create truncated Laplacian
function L = trunclap(A, sub, alpha)

    T = transition_matrix(A(sub, sub));
    L = speye(size(T)) - alpha * T;

author: ahmetius

repo: diffusion-retrieval


text:
classdef IMMOEA_F1 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing IM-MOEA

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = (1+5*repmat(2:obj.D,size(X,1),1)/obj.D).*X(:,2:obj.D) - repmat(X(:,1),1,obj.D-1);
            g = 1 + 9*mean(t.^2,2);
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-sqrt(PopObj(:,1)./g));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function value = r8_round_i4 ( x )

%*****************************************************************************80
%
%% R8_ROUND_I4 rounds an R8 to the nearest integral value, returning an I4.
%
%  Discussion:
%
%    In MATLAB, it is essentially true that there is little difference between
%    this function and R8_ROUND, because we store our integers in what amounts
%    to a real variable.
%
%  Example:
%
%        X        R8_ROUND_I4
%
%      1.3         1
%      1.4         1
%      1.5         1 or 2
%      1.6         2
%      0.0         0
%     -0.7        -1
%     -1.1        -1
%     -1.6        -2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value.
%
%    Output, integer R8_ROUND_I4, the rounded value.
%
  if ( x < 0.0 )
    value = - floor ( - x + 0.5 );
  else
    value =   floor ( + x + 0.5 );
  end

  return
end


text:
function rule_num = keast_rule_num ( DUMMY )

%*****************************************************************************80
%
%% KEAST_RULE_NUM returns the number of Keast rules for the tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Keast,
%    Moderate Degree Tetrahedral Quadrature Formulas,
%    Computer Methods in Applied Mechanics and Engineering,
%    Volume 55, Number 3, May 1986, pages 339-348.
%
%  Parameters:
%
%    Output, integer RULE_NUM, the number of rules.
%
  rule_num = 10;

  return
end


text:
function [map,lo,hi] = cubehelix(N,start,rots,sat,gamma,irange,domain)
% Generate an RGB colormap of Dave Green's Cubehelix colorscheme. With range and domain control.
%
% (c) 2017 Stephen Cobeldick
%
% Returns a colormap with colors defined by Dave Green's Cubehelix colorscheme.
% The colormap nodes are selected along a tapered helix in the RGB color cube,
% with a continuous increase in perceived intensity. Black-and-white printing
% using postscript results in a monotonically increasing grayscale colorscheme.
%
% This function offers two extra controls over the Cubehelix colorscheme:
%  <irange> specifies the intensity levels of the colormap's endnodes (lightness).
%  <domain> subsamples a part of the helix, so the endnodes are color (not gray).
% These options are both explained in the section below 'Range and Domain'.
%
%%% Syntax:
%  map = cubehelix;
%  map = cubehelix(N);
%  map = cubehelix(N,start,rots,sat,gamma);
%  map = cubehelix(N,start,rots,sat,gamma,irange);
%  map = cubehelix(N,start,rots,sat,gamma,irange,domain);
%  map = cubehelix(N,[start,rots,sat,gamma],...)
%  map = cubehelix([],...)
% [map,lo,hi] = cubehelix(...)
%
% Cubehelix is defined here: http://astron-soc.in/bulletin/11June/289392011.pdf
% For more information and examples: http://www.mrao.cam.ac.uk/~dag/CUBEHELIX/
%
% Note: The original specification (the links above) misnamed the saturation
% option as "hue". In this function the saturation option is named "sat".
%
% See also BREWERMAP RGBPLOT COLORMAP COLORBAR SURF CONTOURF IMAGE CONTOURCMAP JET LBMAP
%
%% Range and Domain %%
%
% Using the default <irange> and <domain> vectors ([0,1]) creates colormaps
% exactly the same as Dave Green's original algorithm: from black to white.
%
% The option <irange> sets the intensity level of the colormap's endnodes:
%  cubehelix(3, [0.5,-1.5,1,1], [0.2,0.8]) % irange=[0.2,0.8]
%  ans = 0.2          0.2          0.2     % <- gray, not black
%        0.62751      0.47498      0.28642
%        0.8          0.8          0.8     % <- gray, not white
%
% The option <domain> sets the sampling window for the Cubehelix, such
% that the tapered-helix does not taper all the way to unsaturated (gray).
% This allows the colormap to end with colors rather than gray shades:
%  cubehelix(3, [0.5,-1.5,1,1], [0.2,0.8], [0.3,0.7]) % domain=[0.3,0.7]
%  ans = 0.020144     0.29948      0.15693 % <- color, not gray shade
%        0.62751      0.47498      0.28642
%        0.91366      0.71351      0.95395 % <- color, not gray shade
%
% The function "colormap_view" demonstrates the effects of these options.
%
%% Examples %%
%
%%% New colors for the COLORMAP example:
% load spine
% image(X)
% colormap(cubehelix)
%
%%% New colors for the SURF example:
% [X,Y,Z] = peaks(30);
% surfc(X,Y,Z)
% colormap(cubehelix([],0.7,-0.7,2,1,[0.2,0.8],[0.4,0.8]))
% axis([-3,3,-3,3,-10,5])
%
%% Input and Output Arguments %%
%
%%% Inputs (*=default):
%  N     = NumericScalar, an integer to define the colormap length.
%        = *[], uses the length of the current figure's colormap.
%  start = NumericScalar, *0.5, the helix's start color (modulus 3): R=1, G=2, B=3.
%  rots  = NumericScalar, *-1.5, the number of R->G->B rotations over the scheme length.
%  sat   = NumericScalar, *1,    controls how saturated the colors are.
%  gamma = NumericScalar, *1,    can be used to emphasize low or high intensity values.
%  irange = NumericVector, *[0,1], range of brightness levels of the scheme's endnodes. Size 1x2.
%  domain = NumericVector, *[0,1], domain of the Cubehelix calculation (endnode positions). Size 1x2.
%
%%% Outputs:
%  map = NumericMatrix, a colormap of RGB values between 0 and 1. Size Nx3
%  lo  = LogicalMatrix, true where <map> values<0 were clipped to 0. Size Nx3
%  hi  = LogicalMatrix, true where <map> values>1 were clipped to 1. Size Nx3
%
% [map,lo,hi] = cubehelix(N, start,rots,sat,gamma, irange, domain)
% OR
% [map,lo,hi] = cubehelix(N, [start,rots,sat,gamma], irange, domain)

%% Input Wrangling %%
%
if nargin==0 || (isnumeric(N)&&isempty(N))
	N = size(get(gcf,'colormap'),1);
else
	assert(isnumeric(N)&&isscalar(N),'First input <N> must be a scalar numeric.')
	assert(isreal(N)&&isfinite(N)&&fix(N)==N,'First input <N> must be real and whole: %g+%gi',N,imag(N))
	N = double(N);
end
%
if N==0
	map = ones(0,3);
	lo = false(0,3);
	hi = false(0,3);
	return
end
%
iss = @(x)isnumeric(x)&&isreal(x)&&isscalar(x)&&isfinite(x);
isn = @(x,n)isnumeric(x)&&isreal(x)&&numel(x)==n&&all(isfinite(x(:)));
%
% Parameters:
if nargin<2
	% Default parameter values.
	start = 0.5;
	rots  = -1.5;
	sat   = 1.5;
	gamma = .65;
elseif nargin<5
	% Parameters are in a vector.
	if nargin>2
		irange = rots;
	end
	if nargin>3
		domain = sat;
	end
	assert(isn(start,4)&&isvector(start),'Second input can be a 1x4 real numeric of parameter values.')
	start = double(start);
	gamma = start(4); sat = start(3); rots = start(2); start = start(1);
else
	% Parameters as individual scalar values.
	assert(iss(start),'Input <start> must be a real scalar numeric.')
	assert(iss(rots), 'Input <rots> must be a real scalar numeric.')
	assert(iss(sat),  'Input <sat> must be a real scalar numeric.')
	assert(iss(gamma),'Input <gamma> must be a real scalar numeric.')
	start=double(start); rots=double(rots); sat=double(sat); gamma=double(gamma);
end
%
% Range:
if any(nargin==[0,1,2,5])
	irange = [0,1];
else
	assert(isn(irange,2),'Input <irange> must be a 1x2 real numeric.')
	irange = double(irange);
end
%
% Domain:
if any(nargin==[0,1,2,3,5,6])
	domain = [0,1];
else
	assert(isn(domain,2),'Input <domain> must be a 1x2 real numeric.')
	domain = double(domain);
end
%
%% Core Function %%
%
vec = linspace(domain(1),domain(2),abs(N)).';
ang = 2*pi * (start/3+1+rots*vec);
csm = [cos(ang),sin(ang)].';
fra = vec.^gamma;
amp = sat .* fra .* (1-fra)/2;
%
tmp = linspace(0,1,abs(N)).'.^gamma;
tmp = irange(1)*(1-tmp) + irange(2)*(tmp);
%
cof = [-0.14861,1.78277;-0.29227,-0.90649;1.97294,0];
%
vec = sign(N)*(1:abs(N)) - min(0,N-1);
for m = abs(N):-1:1
	n = vec(m);
	map(m,:) = tmp(n) + amp(n) * (cof*csm(:,n));
end
%
lo = map<0;
hi = map>1;
map = max(0,min(1,map));
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cubehelix
% Copyright (c) 2017 Stephen Cobeldick
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%license

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak