text:
% =========================================================================
% (c) 2016 Ronald Nissel, ronald.nissel@gmail.com
% =========================================================================
% This script simulates an FBMC and OFDM transmission over a doubly-flat
% channel, including channel estimation. The pilot symbol aided channel
% estimation in FBMC is based on R. Nissel, M. Rupp, "On Pilot-Symbol Aided
% Channel Estimation in FBMC-OQAM", IEEE ICASSP, 2016.

clear; close all;
addpath('./Theory');

M_SNR_OFDM_dB = [0:5:30];           % Signal-to-Noise Ratio in dB
NrRepetitions = 1000;               % Number of Monte Carlo repetition (different channel realizations)
QAM_ModulationOrder = 16;           % QAM signal constellation order, 4, 16, 64, 256, 1024,...

%% FBMC Object
FBMC = Modulation.FBMC(...
    12,...                          % Number subcarriers
    30,...                          % Number FBMC symbols
    15e3,...                        % Subcarrier spacing (Hz)
    15e3*14*12,...                  % Sampling rate (Samples/s)
    15e3*20,...                     % Intermediate frequency first subcarrier (Hz)
    false,...                       % Transmit real valued signal
    'Hermite-OQAM',...              % Prototype filter (Hermite, PHYDYAS, RRC) and OQAM or QAM,
    8, ...                          % Overlapping factor (corresponding to the prototype filter length)
    0, ...                          % Initial phase shift
    true ...                        % Polyphase implementation
    );

%% OFDM Object
OFDM = Modulation.OFDM(...
    12,...                          % Number subcarriers
    15,...                          % Number OFDM Symbols
    15e3,...                        % Subcarrier spacing (Hz)
    15e3*14*12,...                  % Sampling rate (Samples/s)
    15e3*20,...                     % Intermediate frequency first subcarrier (Hz)
    false,...                       % Transmit real valued signal
    0, ...                          % Cyclic prefix length (s), LTE: 1/15e3/14
    (8-1/2)*1/15e3*1/2 ...          % Zero guard length (s)
    );

%% PAM and QAM Object
PAM = Modulation.SignalConstellation(sqrt(QAM_ModulationOrder),'PAM');
QAM = Modulation.SignalConstellation(QAM_ModulationOrder,'QAM');

%% Channel Estimation Objects
ChannelEstimation_OFDM = ChannelEstimation.PilotSymbolAidedChannelEstimation(...
    'Diamond',...                           % Pilot pattern
    [...                                    % Matrix that represents the pilot pattern parameters
    OFDM.Nr.Subcarriers,...                 % Number of subcarriers
    6; ...                                  % Pilot spacing in the frequency domain
    OFDM.Nr.MCSymbols,...                   % Number of FBMC/OFDM Symbols
    4 ...                                   % Pilot spacing in the time domain
    ],...
    'linear'...                             % Interpolation(Extrapolation) method 'linear','spline','FullAverage,'MovingBlockAverage'
    );
ChannelEstimation_FBMC = ChannelEstimation.PilotSymbolAidedChannelEstimation(...
    'Diamond',...                           % Pilot pattern
    [...                                    % Matrix that represents the pilot pattern parameters
    FBMC.Nr.Subcarriers,...                 % Number of subcarriers
    6; ...                                  % Pilot spacing in the frequency domain
    FBMC.Nr.MCSymbols,...                   % Number of FBMC/OFDM Symbols
    8 ...                                   % Pilot spacing in the time domain
    ],...
    'linear'...                             % Interpolation(Extrapolation) method 'linear','spline','FullAverage,'MovingBlockAverage',...
    );

%% Imaginary Interference Cancellation Objects
AuxiliaryMethod = ChannelEstimation.ImaginaryInterferenceCancellationAtPilotPosition(...
    'Auxiliary', ...                                    % Cancellation method
    ChannelEstimation_FBMC.GetAuxiliaryMatrix(1), ...   % PilotMatrix
    FBMC.GetFBMCMatrix, ...                             % Imaginary interference matrix
    16, ...                                             % Cancel 16 closest interferers
    2 ...                                               % Pilot to data power offset
    );
CodingMethod = ChannelEstimation.ImaginaryInterferenceCancellationAtPilotPosition(...
    'Coding', ...                                       % Cancellation method
    ChannelEstimation_FBMC.PilotMatrix, ...             % PilotMatrix
    FBMC.GetFBMCMatrix, ...                             % Imaginary interference matrix
    16, ...                                             % Cancel 16 closest interferers
    2 ...                                               % Pilot to data power offset
    );

BER_FBMC_Aux = nan(length(M_SNR_OFDM_dB),NrRepetitions);
BER_FBMC_Cod = nan(length(M_SNR_OFDM_dB),NrRepetitions);
BER_FBMC_perfect = nan(length(M_SNR_OFDM_dB),NrRepetitions);
BER_OFDM = nan(length(M_SNR_OFDM_dB),NrRepetitions);
BER_OFDM_perfect = nan(length(M_SNR_OFDM_dB),NrRepetitions);
for i_rep = 1:NrRepetitions
    for i_SNR = 1:length(M_SNR_OFDM_dB)
        SNR_OFDM_dB = M_SNR_OFDM_dB(i_SNR);
        Pn_time = OFDM.PHY.SamplingRate/(OFDM.PHY.SubcarrierSpacing*OFDM.Nr.Subcarriers)*10^(-SNR_OFDM_dB/10);

        %% Generate Random BitStream
        BinaryDataStream_FBMC_Aux = randi([0 1],AuxiliaryMethod.NrDataSymbols*log2(PAM.ModulationOrder),1);
        BinaryDataStream_FBMC_Cod = randi([0 1],CodingMethod.NrDataSymbols*log2(PAM.ModulationOrder),1);
        BinaryDataStream_OFDM     = randi([0 1],(OFDM.Nr.Subcarriers*OFDM.Nr.MCSymbols-ChannelEstimation_OFDM.NrPilotSymbols)*log2(QAM.ModulationOrder),1);

        %% Transmitted Data Symbols
        xD_FBMC_Aux = PAM.Bit2Symbol(BinaryDataStream_FBMC_Aux);
        xD_FBMC_Cod = PAM.Bit2Symbol(BinaryDataStream_FBMC_Cod);
        xD_OFDM     = QAM.Bit2Symbol(BinaryDataStream_OFDM);

        %% Transmitted Pilot Symbols
        xP_FBMC = PAM.SymbolMapping(randi(PAM.ModulationOrder,[ChannelEstimation_FBMC.NrPilotSymbols 1]));
        xP_FBMC = xP_FBMC./abs(xP_FBMC);
        xP_OFDM = QAM.SymbolMapping(randi(QAM.ModulationOrder,[ChannelEstimation_OFDM.NrPilotSymbols 1]));
        xP_OFDM = xP_OFDM./abs(xP_OFDM);

        %% Transmitted Symbols
        x_FBMC_Aux = reshape(AuxiliaryMethod.PrecodingMatrix*[xP_FBMC;xD_FBMC_Aux],[FBMC.Nr.Subcarriers FBMC.Nr.MCSymbols]);
        x_FBMC_Cod = reshape(CodingMethod.PrecodingMatrix*[xP_FBMC;xD_FBMC_Cod],[FBMC.Nr.Subcarriers FBMC.Nr.MCSymbols]);
        x_OFDM = nan(OFDM.Nr.Subcarriers,OFDM.Nr.MCSymbols);
        x_OFDM(ChannelEstimation_OFDM.PilotMatrix==1) = xP_OFDM;
        x_OFDM(ChannelEstimation_OFDM.PilotMatrix==0) = xD_OFDM;

        %% Transmitted FBMC Signal (time domain)
        s_FBMC_Aux = FBMC.Modulation(x_FBMC_Aux);
        s_FBMC_Cod = FBMC.Modulation(x_FBMC_Cod);
        s_OFDM     = OFDM.Modulation(x_OFDM);

        %% Channel (doubly flat fading and AWGN)
        h = sqrt(1/2)*(randn+1j*randn);
    %     h = 1; % Pure AWGN
        n_FBMC = sqrt(Pn_time/2)*(randn(size(s_FBMC_Cod))+1j*randn(size(s_FBMC_Cod)));
        n_OFDM = sqrt(Pn_time/2)*(randn(size(s_OFDM))+1j*randn(size(s_OFDM)));

        r_FBMC_Aux = h*s_FBMC_Aux + n_FBMC;
        r_FBMC_Cod = h*s_FBMC_Cod + n_FBMC;
        r_OFDM     = h*s_OFDM + n_OFDM;

        %% Demodulate OFDM and FBMC signal
        y_FBMC_Aux = FBMC.Demodulation(r_FBMC_Aux);
        y_FBMC_Cod = FBMC.Demodulation(r_FBMC_Cod);
        y_OFDM     = OFDM.Demodulation(r_OFDM);

        %% LS channel estimates at pilot positions
        hP_LS_FBMC_Aux = y_FBMC_Aux(ChannelEstimation_FBMC.PilotMatrix==1)./xP_FBMC/sqrt(AuxiliaryMethod.PilotToDataPowerOffset*AuxiliaryMethod.DataPowerReduction);
        hP_LS_FBMC_Cod = y_FBMC_Cod(ChannelEstimation_FBMC.PilotMatrix==1)./xP_FBMC/sqrt(CodingMethod.PilotToDataPowerOffset);
        hP_LS_OFDM     = y_OFDM(ChannelEstimation_OFDM.PilotMatrix==1)./xP_OFDM;

        %% Channel Estimation using Interpolation
        h_FBMC_Aux = ChannelEstimation_FBMC.ChannelInterpolation(hP_LS_FBMC_Aux);
        h_FBMC_Cod = ChannelEstimation_FBMC.ChannelInterpolation(hP_LS_FBMC_Cod);
        h_OFDM     = ChannelEstimation_OFDM.ChannelInterpolation(hP_LS_OFDM);

        %% Equalized received symbols at data position
        y_EQ_FBMC_Aux = real(y_FBMC_Aux(AuxiliaryMethod.PilotMatrix==0)./h_FBMC_Aux(AuxiliaryMethod.PilotMatrix==0)/sqrt(AuxiliaryMethod.DataPowerReduction));
        y_EQ_FBMC_Cod = real(CodingMethod.PrecodingMatrix(:,CodingMethod.NrPilotSymbols+1:end)'*(y_FBMC_Cod(:)./h_FBMC_Cod(:)));
        y_EQ_FBMC_perfect = real(CodingMethod.PrecodingMatrix(:,CodingMethod.NrPilotSymbols+1:end)'*(y_FBMC_Cod(:)./h));

        y_EQ_OFDM = y_OFDM(ChannelEstimation_OFDM.PilotMatrix==0)./h_OFDM(ChannelEstimation_OFDM.PilotMatrix==0);
        y_EQ_OFDM_perfect = y_OFDM(ChannelEstimation_OFDM.PilotMatrix==0)./h;

        %% Detect BitStream
        DetectedBitStream_FBMC_Aux = PAM.Symbol2Bit(real(y_EQ_FBMC_Aux(:)));
        DetectedBitStream_FBMC_Cod = PAM.Symbol2Bit(real(y_EQ_FBMC_Cod(:)));
        DetectedBitStream_FBMC_perfect = PAM.Symbol2Bit(real(y_EQ_FBMC_perfect(:)));

        DetectedBitStream_OFDM = QAM.Symbol2Bit(y_EQ_OFDM(:));
        DetectedBitStream_OFDM_perfect = QAM.Symbol2Bit(y_EQ_OFDM_perfect(:));

        %% Calculate BER
        BER_FBMC_Aux(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Aux~=DetectedBitStream_FBMC_Aux);
        BER_FBMC_Cod(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Cod~=DetectedBitStream_FBMC_Cod);
        BER_FBMC_perfect(i_SNR,i_rep) = mean(BinaryDataStream_FBMC_Cod~=DetectedBitStream_FBMC_perfect);

        BER_OFDM(i_SNR,i_rep) = mean(BinaryDataStream_OFDM~=DetectedBitStream_OFDM);
        BER_OFDM_perfect(i_SNR,i_rep) = mean(BinaryDataStream_OFDM~=DetectedBitStream_OFDM_perfect);

    end

    if mod(i_rep,100)==0
       disp([int2str(i_rep/NrRepetitions*100) '%']);
    end
end

%% Theoretical BEP for perfect channel knowledge
% BEP_4QAM = 1/2-1./(2*sqrt(2*(1+10.^(-M_SNR_OFDM_dB/10))-1));
M_SNR_OFDM_dB_morePoints = min(M_SNR_OFDM_dB):0.5:max(M_SNR_OFDM_dB);
BEP_perfect = BitErrorProbabilityDoublyFlatRayleigh(M_SNR_OFDM_dB_morePoints,QAM.SymbolMapping,QAM.BitMapping);

%% Plot BER and BEP
figure();
semilogy(M_SNR_OFDM_dB,mean(BER_FBMC_Aux,2),'red -o');
hold on;
semilogy(M_SNR_OFDM_dB,mean(BER_FBMC_Cod,2),'blue -o');
semilogy(M_SNR_OFDM_dB,mean(BER_OFDM,2),'black -o');
semilogy(M_SNR_OFDM_dB,mean(BER_FBMC_perfect,2),'blue -x');
semilogy(M_SNR_OFDM_dB,mean(BER_OFDM_perfect,2),'black -x');
semilogy(M_SNR_OFDM_dB_morePoints,BEP_perfect','black');
xlabel('SNR for OFDM (dB)');
ylabel('BER, BEP');
legend('Simulation: FBMC Auxiliary','Simulation: FBMC Coding','Simulation: OFDM', 'Simulation FBMC perfect CSI', 'Simulation OFDM perfect CSI','Theory perfect CSI','Location','SouthWest');

%% Plot Pilot Pattern
figure();
ChannelEstimation_OFDM.PlotPilotPattern;
title('OFDM');
figure();
ChannelEstimation_FBMC.PlotPilotPattern(AuxiliaryMethod.PilotMatrix)
title('FBMC Auxiliary');
figure();
ChannelEstimation_FBMC.PlotPilotPattern(-(CodingMethod.ConsideredInterferenceMatrix<0)+(CodingMethod.ConsideredInterferenceMatrix>0))
title('FBMC Coding');

%% Calculate and Plot Expected Transmit Power Over Time
[Power_FBMC_Aux,t_FBMC] = FBMC.PlotTransmitPower(AuxiliaryMethod.PrecodingMatrix*AuxiliaryMethod.PrecodingMatrix');
[Power_FBMC_Cod,~] = FBMC.PlotTransmitPower(CodingMethod.PrecodingMatrix*CodingMethod.PrecodingMatrix');
[Power_OFDM,t_OFDM] = OFDM.PlotTransmitPower;
figure();
plot(t_FBMC,Power_FBMC_Aux,'red');
hold on;
plot(t_FBMC,Power_FBMC_Cod,'blue');
plot(t_OFDM,Power_OFDM,'black ');
legend({'FBMC Auxiliary','FBMC Coding','OFDM'});
ylabel('Transmit Power');
xlabel('Time(s)');

%% Calculate Power Spectral Density
[PSD_FBMC_Aux,t_FBMC] = FBMC.PlotPowerSpectralDensity(AuxiliaryMethod.PrecodingMatrix*AuxiliaryMethod.PrecodingMatrix');
[PSD_FBMC_Cod,~] = FBMC.PlotPowerSpectralDensity(CodingMethod.PrecodingMatrix*CodingMethod.PrecodingMatrix');
[PSD_OFDM,t_OFDM] = OFDM.PlotPowerSpectralDensity;
figure();
plot(t_FBMC,10*log10(PSD_FBMC_Aux),'red');
hold on;
plot(t_FBMC,10*log10(PSD_FBMC_Cod),'blue');
plot(t_OFDM,10*log10(PSD_OFDM),'black ');
legend({'FBMC Auxiliary','FBMC Coding','OFDM'});
ylabel('Power Spectral Density (dB)');
xlabel('Frequency (Hz)');

author: rnissel

repo: Channel-Estimation


text:
function [output_signal,output_time] = mvstat_final(Inputdata_signal,Inputdata_time,windowrange,reductionrange,fn)
%------------ Introduction ----------------------
%This file is created to calculate the mathematical operation like 'max'
%on a pre-defined window range and with a pre-defined reduction range of an vector data.
%This function will need the following inputs in order to calculate the
%wished results:
%   1- Inputdata_signal: A vector containg your Y-Values that need to be
%   filtered.
%   2- Inputdata_time: A vector containg your x-Values that need to be
%   filtered (Normally is a time signal where the sampling rate can be
%   extracted).
%   3- windowrange: in unit of time (s) is the range to calculate the needed
%   mathematical function (The moving range ex. mvmean(y,x,0.02,0.02) or mvrms(y,x,0.02,0.02))
%   4- reductionrange: in unit of time (s) is the scalling reference for
%   the new data
%   5- fn : is used to define the wished mathmatical operation (ex. @max, @min, @mean, @rms, @int)
%The relation between the windowrange input and the reductionrange input
%must be an integer value of 1:1.2:1..10:1. If not the reductionrange would
%be corrected to get one of those relation. this specification give us the
%possibility to calculate the 'mvstat_final' without using a for-Loop which will save a whole amount of time
% for example:
% fn = @mean
% [output_signal,out_put_time] = mvstat_final (y,x,10,10)
% the moving average will be calculated as a result of this function. every
% 10 s the mean value for the last 10 s will be calculated.
% fn = @max
% [output_signal,out_put_time] = mvstat_final (y,x,20,10)
% the moving max will be calculated as a result of this function. every
% 10 s the maximum value for the last 20 s will be calculated.
% fn = @min
% [output_signal,out_put_time] = mvstat_final (y,x,5,10)
% the moving min will be calculated as a result of this function. every
% 10 s the minmum value for the last 5 s will be calculated.
%Notice: when the input values of windowrange and reductionrange are not
%the same then at the beginning of the calculation the fn value of the a
%available data will be calculated.
% Is m.file will also allows you to calculate the intgration of your input
% data in a spesific windowrange. in this case the reductionrange input
% should be defined as in the example:
%fn = @int;
%[Ua_sin,Ua_sin_time] = mvstat(Ua_sin,U1_time,(1/f_grid),delta_t,fn); Ua_sin = Ua_sin*2*f_grid;
% in this example the integration of your input data (in this case U1_sin)
% will be calculated ever one period (normally this mean 0.02 s or 1/50 or
% 1/f_grid in our case). and then the reductionrange is in this case equal
% to the sampling time of the inputdata (delta = U1_time(2)-U1_time(1)).
% this work was done in order to bring FAMOS (IMC) software function into
% matlab.
%Auther: Msc.Eng. Aubai Alkhatib
%Datum: 10.09.2013
%-------- End of introduction ------------------
delta = Inputdata_time(2)-Inputdata_time(1);
length_window = windowrange/delta;
length_reduction = reductionrange/delta;
time = Inputdata_time(end);
if windowrange > reductionrange
    number = windowrange/reductionrange;
    integer = floor(number);
    fract = number-integer;
    if fract >= 0.5
        if 0 == 1

        else
            A = windowrange/(integer+1);
        end
    else
        A = windowrange/(integer);
    end
    reductionrange_new = A;
     integer_A = floor(A);
    fract_A = A-integer_A;
    if 0 == 1
        if windowrange/(reductionrange+1) < windowrange/reductionrange
            reductionrange_new = reductionrange;
            while rem(windowrange,reductionrange_new) > 0
                reductionrange_new = reductionrange_new - 1;
            end
        else
            reductionrange_new = reductionrange;
            while rem(windowrange,reductionrange_new) == 0
                reductionrange_new = reductionrange_new + 1;
            end
        end
    end
    [~,b] = rat(A);
    if b > 1
        length_reduction_new = round(reductionrange_new/delta);
    else
        length_reduction_new = reductionrange_new/delta;
    end
    i = length_reduction_new;
    ii = 1;
    output = zeros();
    while i < length(Inputdata_signal)
        if i >= length_window
            if fract_A == 0
                switch func2str(fn)
                    case 'mean'
                        output(ii) = mean(Inputdata_signal((i-length_window)+1:i));
                    case 'max'
                        output(ii) = max(Inputdata_signal((i-length_window)+1:i));
                    case 'sum'
                        output(ii) = sum(Inputdata_signal((i-length_window)+1:i));
                    case 'min'
                        output(ii) = min(Inputdata_signal((i-length_window)+1:i));
                    case 'int'
                        output(ii) = sum(Inputdata_signal((i-length_window)+1:i))*reductionrange_new;
                    otherwise
                        output(ii) = rms(Inputdata_signal((i-length_window)+1:i));
                end
            else
                switch func2str(fn)
                    case 'mean'
                        output(ii) = mean(Inputdata_signal((i-length_window):i));
                    case 'max'
                        output(ii) = max(Inputdata_signal((i-length_window):i));
                    case 'sum'
                        output(ii) = sum(Inputdata_signal((i-length_window):i));
                    case 'min'
                        output(ii) = min(Inputdata_signal((i-length_window):i));
                    case 'int'
                        output(ii) = sum(Inputdata_signal((i-length_window)+1:i))*reductionrange_new;
                    otherwise
                        output(ii) = rms(Inputdata_signal((i-length_window):i));
                end
            end
        else
            switch func2str(fn)
                case 'mean'
                    output(ii) = mean(Inputdata_signal(1:i));
                case 'max'
                    output(ii) = max(Inputdata_signal(1:i));
                case 'sum'
                    output(ii) = sum(Inputdata_signal(1:i));
                case 'int'
                    output(ii) = sum(Inputdata_signal(1:i))*reductionrange_new;
                case 'min'
                    output(ii) = min(Inputdata_signal(1:i));
                otherwise
                    output(ii) = rms(Inputdata_signal(1:i));
            end
        end
        i = i + length_reduction_new;
        ii = ii + 1;
    end

    if strcmp(func2str(fn),'int')
        output_signal = output(length_window:length(output));
        output_time = windowrange:reductionrange_new:(length(output)*reductionrange_new);
    elseif isinteger(reductionrange_new)
        output_signal = output;
        output_time = 0:reductionrange_new:(length(output)*reductionrange_new)-reductionrange_new;
    else
        output_signal = output;
        output_time = 0:(length_reduction_new*delta):(length(output)*(length_reduction_new*delta))-(length_reduction_new*delta);
    end
elseif reductionrange > windowrange
    Length_Output = fix(time/windowrange);
    Length_Input_New = Length_Output*windowrange/delta;
    delta_new = length_window;
    delta_time = windowrange;
    Input_reshaped = reshape(Inputdata_signal(1:Length_Input_New),delta_new,Length_Input_New/delta_new);
    switch func2str(fn)
        case 'mean'
            output_signal = mean(Input_reshaped);
        case 'max'
            output_signal = max(Input_reshaped);
        case 'sum'
            output_signal = sum(Input_reshaped);
        case 'int'
            output = sum(Input_reshaped)*delta_time;
        case 'min'
            output_signal = min(Input_reshaped);
        otherwise
            output_signal = rms(Input_reshaped);
    end
    if strcmp(func2str(fn),'int')
        output_signal = output(length_window:length(output));
        output_time = windowrange:delta_time:(Length_Output*delta_time)-delta_time;
    else
        output_time = 0:delta_time:(Length_Output*delta_time)-delta_time;
    end
elseif windowrange == reductionrange
    Length_Output = fix(time/windowrange);
    Length_Input_New = Length_Output*reductionrange/delta;
    delta_new = length_reduction;
    delta_time = windowrange;
    Input_reshaped = reshape(Inputdata_signal(1:Length_Input_New),delta_new,Length_Input_New/delta_new);
    switch func2str(fn)
        case 'mean'
            output_signal = mean(Input_reshaped);
        case 'max'
            output_signal = max(Input_reshaped);
        case 'sum'
            output_signal = sum(Input_reshaped);
        case 'int'
            output = sum(Input_reshaped)*delta_time;
        case 'min'
            output_signal = min(Input_reshaped);
        otherwise
            output_signal = rms(Input_reshaped);
    end
    if strcmp(func2str(fn),'int')
        output_signal = output(length_window:length(output));
        output_time = windowrange:delta_time:(Length_Output*delta_time)-delta_time;
    else
        output_time = 0:delta_time:(Length_Output*delta_time)-delta_time;
    end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function [scanArray3M,mask3M] = perturbImageAndSeg(scanArray3M,mask3M,planC,...
         scanNum,perturbString,angle2DStdDeg,volScaleStdFraction,superPixVol)
% [scanArray3M,mask3M] = perturbImageAndSeg(scanArray3M,mask3M,planC,scanNum,...
%                        perturbString,angle2DStdDeg,volScaleStdFraction,superPixVol)
%
% perturbString is a string containing perturbation operations. It can be a
% combination of any letters RNTVC.
%
% APA, 2/28/2019
% AI , 3/22/2019 Added inputs angle2DStdDeg,volScaleStdFraction,superPixVol

% Noise
if ismember('N',perturbString)
    scanArray3M = perturbImageNoise(scanArray3M);
end

% Rotation
if ismember('R',perturbString)
    angl = normrnd(0,angle2DStdDeg);
    [scanArray3M,mask3M] = perturbImageRotation(scanArray3M,mask3M,angl);
end

% Volume adaptation
if ismember('V',perturbString)
    scale = normrnd(1,volScaleStdFraction);
    mask3M = perturbImageVolume(mask3M,scale);
end

% Contour randomization
if ismember('C',perturbString)
    mask3M = perturbImageContourBySuperpix(mask3M,scanArray3M,superPixVol,planC,scanNum);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%kkgenman 'kgenman: generate man3 pages from source code headers'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgenman.pane file
%
% Parameters:
% String: tb 'Toolbox Name', default: ' ': 'toolbox the library belongs to'
% String: lib 'Library Name', default: ' ': 'library to generate man pages for'
% String: section 'Man Section', default: ' ': 'Man section to put the man page into'
% String: dir 'Generate man pages and install', default: '.': 'generate and install in the directory specified'
% Toggle: install 'Generate and install the man pages in $TOOLBOX/man/man3', default: 0: 'Install the man pages in $TOOLBOX/man/man3'
% Toggle: list 'List Man Pages Only', default: 0: 'List the man pages rather than generate them'
% Toggle: all 'All src files in library object', default: 0: 'create man pages for all src in the library object'
% String: i 'Input Source Code File', default: ' ': 'input source code file'
%
% Example: kkgenman( {'tb',' ';'lib',' ';'section',' ';'dir','.';'install',0;'list',0;'all',0;'i',' '})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgenman - generate man page(s) for a software object
%
%  DESCRIPTION
%
% This program parses a khoros public routine header for information, and
% generates a man 3 page for the public library routine.  This routine takes a
% as input the library object the user wants man 3 pages for, and creates all
% the man 3 files for every public routine in a single file, or in the whole
% library.  It then installs the man page in the directory of choice, or
% the final installation directory in $TOOLBOX/man/object_name/man.
%
% Note, all the information stored between the lines ".SH ADDITIONAL INFO" and
% ".SH SIDE EFFECTS" is saved and propagated to the new man page.  Thus,
% the user can put extra information in this section of the man page, such as an
% example section, if desired.
%
% This program will detect when headers have missing fields, and recover
% gracefully.  If a field is missing, it will indicate which field, routine,
% and filename the error occurred in, and continue on with that field left
% empty.  Therefore, it is the programmers responsibility to fix the header,
% and re-run kgenman on that file, so that the man page will be complete.
%
% The example field below has an example of the public header along with an
% explaination of each field that goes in it.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkgenman(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kkgenman(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'tb', ' ';'lib', ' ';'section', ' ';'dir', '.';'install', 0;'list', 0;'all', 0;'i', ' '};
maxval={0,0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0,0};
istoggle=[0,0,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'String','String','String','String','Toggle','Toggle','Toggle','String'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'kgenman"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function write_images_to_webrtDB(patient_id)

global stateS planC
indexS = planC{end};

% absolutePathForImageFiles = fullfile('/Users/aptea/Documents/MSKCC/Projects/webRT',['patient_',num2str(patient_id)]);

absolutePathForImageFiles = fullfile('/Volumes/DeasyLab3/Repository/RIVIEW',['patient_',num2str(patient_id)]);

% Open planC in CERR viewer
if isfield(stateS,'handle')
    hCSV = stateS.handle.CERRSliceViewer;
else
    hCSV = [];
end
if isempty(hCSV) || ~exist('hCSV') || ~ishandle(hCSV)
    CERR('CERRSLICEVIEWER')
end

stateS.CTToggle = 1;
stateS.CTDisplayChanged = 1;
sliceCallBack('OPENWORKSPACEPLANC')

% Capture Scan Images and store pointers in DB

%Find the scan with scanUID matching this scan's CERR scan UID

for scanNum = 1:length(planC{indexS.scan})

    % Write scan images to file
    write_scan_images_to_disk(absolutePathForImageFiles,scanNum,'transverse')
    write_scan_images_to_disk(absolutePathForImageFiles,scanNum,'coronal')
    write_scan_images_to_disk(absolutePathForImageFiles,scanNum,'sagittal')

    % Write dose images to file
    dosesV = getScanAssociatedDose(scanNum,'all');
    for doseNum = dosesV
        write_dose_images_to_disk(absolutePathForImageFiles,doseNum,scanNum,'transverse')
        write_dose_images_to_disk(absolutePathForImageFiles,doseNum,scanNum,'coronal')
        write_dose_images_to_disk(absolutePathForImageFiles,doseNum,scanNum,'sagittal')
    end

%     % Toggle plane locators off
%
%     % Display Transverse View in larger window
%
%     % Toggle Dose off
%
%     % Toggle Structures off
%
%     % Capture Axis and save scan as png file
%
%     % Toggle Dose on
%
%     % Toggle Scan off
%
%     % Capture Axis and save dose as png file
%
%
%     % Find structures associated to this scan
%     assocScanV = getStructureAssociatedScan(1:numStructs);
%     structsInScanV = find(assocScanV == scanNum);
%
%     % Loop over structures to store contours in the database
%
%     for structNum = structsInScanV
%
%         px = axes2pix(ncols,xLim,planC{indexS.structures}(12).contour(74).segments.points(:,1));
%         py = axes2pix(nrows,-yLim,planC{indexS.structures}(12).contour(74).segments.points(:,2));
%         ptsM = [px py];
%         aa = ptsM';
%         mat2str(aa(:))
%
%     end

end


text:
function p = basicHKL(cs,varargin)
% plot symmetry
%
% Input
%  cs - symmetry
%
% Output
%
% Options
%  antipodal      - include <VectorsAxes.html antipodal symmetry>

[h,k,l] = meshgrid(-1:1);

p = Miller(h(:),k(:),l(:),cs);

p(isinf(p.dspacing)) = [];

p = unique(p.symmetrise,'noSymmetry');

author: mtex-toolbox

repo: mtex

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


text:
classdef MaxFixedStepsReachedValidator < AbstractLaunchVehicleDataValidator
    %MaxSimTimeReachedValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = MaxFixedStepsReachedValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            warnEvtNums = [];
            evts = obj.lvdData.script.evts;
            for(i=1:length(evts)) %#ok<*NO4LP>
                evt = evts(i);

                integratorOptions = evt.integratorObj.getOptions();
                integrationStep = integratorOptions.getIntegratorStepSize();
                integrationMaxFixedSteps = integratorOptions.getIntegratorMaxNumFixedSteps();

                if(integrationStep > 0) %fixed step requested
                    stateLogEntries = obj.lvdData.stateLog.getAllStateLogEntriesForEvent(evt);
                    numEntries = numel(stateLogEntries);

                    if(numEntries >= integrationMaxFixedSteps)
                        warnEvtNums(end+1) = evt.getEventNum(); %#ok<AGROW>
                    end
                end

                if(not(isempty(warnEvtNums)))
                    eventStr = makeEventsStr(unique(warnEvtNums));
                    str = sprintf('Maximum number of fixed step sizes reached on some events. (Events: %s)\n%s', eventStr);
                    warnings(end+1) = LaunchVehicleDataValidationWarning(str); %#ok<AGROW>
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%
% References:
%
% C. Lu. A Library of ADMM for Sparse and Low-rank Optimization. National University of Singapore, June 2016.
% https://github.com/canyilu/LibADMM.
% C. Lu, J. Feng, S. Yan, Z. Lin. A Unified Alternating Direction Method of Multipliers by Majorization
% Minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 40, pp. 527-541, 2018
%

addpath(genpath(cd))
clear

%% Examples for testing the low-rank matrix based models
% For detailed description of the sparse models, please refer to the Manual.

%% generate toy data
d = 10;
na = 200;
nb = 100;

A = randn(d,na);
X = randn(na,nb);
B = A*X;
b = B(:,1);

opts.tol = 1e-6;
opts.max_iter = 1000;
opts.rho = 1.2;
opts.mu = 1e-3;
opts.max_mu = 1e10;
opts.DEBUG = 0;

%% RPCA
n1 = 100;
n2 = 200;
r = 10;
L = rand(n1,r)*rand(r,n2); % low-rank part

p = 0.1;
m = p*n1*n2;
temp = rand(n1*n2,1);
[~,I] = sort(temp);
I = I(1:m);
Omega = zeros(n1,n2);
Omega(I) = 1;
E = sign(rand(n1,n2)-0.5);
S = Omega.*E; % sparse part, S = P_Omega(E)

Xn = L+S;

lambda = 1/sqrt(max(n1,n2));
opts.loss = 'l1';
opts.DEBUG = 1;
tic
[Lhat,Shat,obj,err,iter] = rpca(Xn,lambda,opts);
toc
rel_err_L = norm(L-Lhat,'fro')/norm(L,'fro')
rel_err_S = norm(S-Shat,'fro')/norm(S,'fro')

err
iter

%% low rank matrix completion (lrmc) and regularized lrmc

n1 = 100;
n2 = 200;
r = 5;
X = rand(n1,r)*rand(r,n2);

p = 0.6;
omega = find(rand(n1,n2)<p);
M = zeros(n1,n2);
M(omega) = X(omega);
[Xhat,obj,err,iter] = lrmc(M, omega, opts);
rel_err_X = norm(Xhat-X,'fro')/norm(X,'fro')

E = randn(n1,n2)/100;
M = X+E;
lambda = 0.1;
[Xhat,obj,err,iter] = lrmcR(M, omega, lambda, opts);

%% low rank representation (lrr)
lambda = 0.001;
opts.loss = 'l21';
tic
[X,E,obj,err,iter] = lrr(A,A,lambda,opts);
toc
obj
err
iter

%% latent LRR (latlrr)
lambda = 0.1;
opts.loss = 'l1';
tic
[Z,L,obj,err,iter] = latlrr(A,lambda,opts);
toc
obj
err
iter

%% low rank and sparse representation (lrsr)
lambda1 = 0.1;
lambda2 = 4;
opts.loss = 'l21';
tic
[X,E,obj,err,iter] = lrsr(A,B,lambda1,lambda2,opts);
toc
obj
err
iter

%% improved graph clustering (igc)
n = 100;
r = 5;
X = rand(n,r)*rand(r,n);
C = rand(size(X));
lambda = 1/sqrt(n);
opts.loss = 'l1';
opts.DEBUG = 1;
tic
[L,S,obj,err,iter] = igc(X,C,lambda,opts);
toc
err
iter

%% multi-task low-rank affinity pursuit (mlap)
n1 = 100;
n2 = 200;
K = 10;
X = rand(n1,n2,K);
lambda = 0.1;
alpha = 0.2;
opts.loss = 'l1';
tic
[Z,E,obj,err,iter] = mlap(X,lambda,alpha,opts);
toc
err
iter

%% robust multi-view spectral clustering (rmsc)
n = 100;
r = 5;
m = 10;
X = rand(n,n,m);
lambda = 1/sqrt(n);
opts.loss = 'l1';
opts.DEBUG = 1;
tic
[L,S,obj,err,iter] = rmsc(X,lambda,opts);
toc
err
iter

%% sparse spectral clustering (sparsesc)
lambda = 0.001;
n = 100;
X = rand(n,n);
W = abs(X'*X);
I = eye(n);
D = diag(sum(W,1));
L = I - sqrt(inv(D))*W*sqrt(inv(D));
k = 5;
[P,obj,err,iter] = sparsesc(L,lambda,k,opts);
obj
err
iter

author: canyilu

repo: LibADMM-toolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Y, W, Thres, sigma] = TwoLayerAGH_Train(X, Anchor, r, s, sigma, options)
%
% Two-Layer Anchor Graph Hashing Training
% Written by Wei Liu (wliu@ee.columbia.edu)
% X = nXdim input data
% Anchor = mXdim anchor points (m<<n)
% r = number of hash bits
% s = number of nearest anchors
% sigma: Gaussian RBF kernel width parameter; if no input (sigma=0) this function will
%        self-estimate and return a value.
% Y = nXr binary codes (Y_ij in {1,0})
% W = mX(r/2) projection matrix in spectral space
% Thres = 2X(r/2) threshold matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,dim] = size(X);
m = size(Anchor,1);

%% get Z
Z = zeros(n,m);
Dis = sqdist(X',Anchor');
% clear X;
% clear Anchor;

val = zeros(n,s);
pos = val;
for i = 1:s
    [val(:,i),pos(:,i)] = min(Dis,[],2);
    tep = (pos(:,i)-1)*n+[1:n]';
    Dis(tep) = 1e60;
end
clear Dis;
clear tep;

if sigma == 0
   sigma = mean(val(:,s).^0.5);
end

switch lower(options.CodingMethod)
    case {lower('Cosine')}
        if isfield(options,'bNormalized') && options.bNormalized
        else
            X = NormalizeFea(X);
        end
        Anchor = NormalizeFea(Anchor);
        S = full(X*Anchor');
        linidx = [1:size(pos,1)]';
        val = S(sub2ind(size(S),linidx(:,ones(1,size(pos,2))),pos));
    case {lower('Gaussian')}
        val = exp(-val/(1/1*sigma^2));
    otherwise
        error('method does not exist!');
end

val = repmat(sum(val,2).^-1,1,s).*val; %% normalize
tep = (pos-1)*n+repmat([1:n]',1,s);
Z([tep]) = [val];
Z = sparse(Z);
clear tep;
clear val;
clear pos;

%% compute eigensystem
lamda = sum(Z);
M = Z'*Z;
M = diag(lamda.^-0.5)*M*diag(lamda.^-0.5);
[W,V] = eig(full(M));
clear M;
eigenvalue = diag(V)';
clear V;
[eigenvalue,order] = sort(eigenvalue,'descend');
W = W(:,order);
clear order;
ind = find(eigenvalue>0 & eigenvalue<1-1e-3);
eigenvalue = eigenvalue(ind);
W = W(:,ind);
W = diag(lamda.^-0.5)*W(:,1:r/2)*diag(eigenvalue(1:r/2).^-0.5)*sqrt(n);
clear ind;
clear eigenvalue;
clear lambda;

%% get binary codes
Y = Z*W;
Y1 = Y;
Thres = zeros(2,r/2);
for k = 1:r/2
    pos = find(Y(:,k)>0);
    neg = setdiff([1:n],pos);
    n1 = length(pos);
    n2 = n-n1;
    get = Y(:,k);
    get(neg) = -1*get(neg);
    right = full(get'*Z);
    clear get;
    right = right./lamda;
    left = sum(Z(pos,:));
    s1 = sum(left.*right,2);
    left1 = left./lamda;
    s2 = n1-sum(left.*left1,2);
    s3 = sum(Y(pos,k),1);
    clear left
    clear left1;
    clear right;

    Thres(1,k) = (2*s3+n2*(s3-s1)/s2)/n;
    Thres(2,k) = (n1*(s3-s1)/s2-2*s3)/n;
    Y1(pos,k) = Y1(pos,k)-Thres(1,k);
    Y1(neg,k) = -Y1(neg,k)+Thres(2,k);
    clear pos;
    clear neg;
end
clear Z;
Y = [Y,Y1];
clear Y1;

Y = (Y>0); %% logical format

author: ZJULearning

repo: MatlabFunc


text:
% LTFAT - Block processing
%
%  Zdenek Prusa, 2013 - 2018.
%
%  Basic methods
%    BLOCK          - Setup a new block-stream
%    BLOCKDEVICES   - List available audio I/O devices
%    BLOCKREAD      - Read samples from file/device
%    BLOCKPLAY      - Play block (sound output)
%    BLOCKPANEL     - Block-stream control GUI
%    BLOCKPANELGET  - Obtain parameter(s) from GUI
%    BLOCKDONE      - Closes block-stream and frees resources
%    BLOCKWRITE     - Appends data to a wav file
%
%  Block-adapted transforms
%    BLOCKFRAMEACCEL     - Prepare a frame for a block-stream processing
%    BLOCKFRAMEPAIRACCEL - Prepare a pair of frames for a block-stream processing
%    BLOCKANA            - Block analysis
%    BLOCKSYN            - Block synthesis
%
%  Running visualisation
%    BLOCKFIGURE   - Initialize figure for redrawing
%    BLOCKPLOT     - Append coefficients to the running plot
%
%  Other
%    LTFATPLAY     - Replacement for the sound command allowing selecting an output device
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function p = prior_gamma(varargin)
%PRIOR_GAMMA  Gamma prior structure
%
%  Description
%    P = PRIOR_GAMMA('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates Gamma prior structure in which the named parameters
%    have the specified values. Any unspecified parameters are set
%    to default values.
%
%    P = PRIOR_GAMMA(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parametrisation is done by Bayesian Data Analysis,
%    second edition, Gelman et.al. 2004.
%
%    Parameters for Gamma prior [default]
%      sh       - shape [4]
%      is       - inverse scale [1]
%      sh_prior - prior for sh [prior_fixed]
%      is_prior - prior for is [prior_fixed]
%
%  See also
%    PRIOR_*

% Copyright (c) 2000-2001,2010 Aki Vehtari
% Copyright (c) 2010 Jaakko Riihim�ki

% 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.FunctionName = 'PRIOR_GAMMA';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('sh',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('sh_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('is',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('is_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Gamma';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Gamma')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('sh',ip.UsingDefaults)
    p.sh = ip.Results.sh;
  end
  if init || ~ismember('is',ip.UsingDefaults)
    p.is = ip.Results.is;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('sh_prior',ip.UsingDefaults)
    p.p.sh=ip.Results.sh_prior;
  end
  if init || ~ismember('is_prior',ip.UsingDefaults)
    p.p.is=ip.Results.is_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_gamma_pak;
    p.fh.unpak = @prior_gamma_unpak;
    p.fh.lp = @prior_gamma_lp;
    p.fh.lpg = @prior_gamma_lpg;
    p.fh.recappend = @prior_gamma_recappend;
  end

end

function [w, s] = prior_gamma_pak(p)

  w=[];
  s={};
  if ~isempty(p.p.sh)
    w = log(p.sh);
    s=[s; 'log(Gamma.sh)'];
  end
  if ~isempty(p.p.is)
    w = [w log(p.is)];
    s=[s; 'log(Gamma.is)'];
  end
end

function [p, w] = prior_gamma_unpak(p, w)

  if ~isempty(p.p.sh)
    i1=1;
    p.sh = exp(w(i1));
    w = w(i1+1:end);
  end
  if ~isempty(p.p.is)
    i1=1;
    p.is = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_gamma_lp(x, p)

  lp = sum(-p.is.*x + (p.sh-1).*log(x) +p.sh.*log(p.is)  -gammaln(p.sh));

  if ~isempty(p.p.sh)
    lp = lp + p.p.sh.fh.lp(p.sh, p.p.sh) + log(p.sh);
  end
  if ~isempty(p.p.is)
    lp = lp + p.p.is.fh.lp(p.is, p.p.is) + log(p.is);
  end
end

function lpg = prior_gamma_lpg(x, p)

  lpg = (p.sh-1)./x - p.is;

  if ~isempty(p.p.sh)
    lpgsh = (sum(-digamma1(p.sh) + log(p.is) + log(x)) + p.p.sh.fh.lpg(p.sh, p.p.sh)).*p.sh + 1;
    lpg = [lpg lpgsh];
  end
  if ~isempty(p.p.is)
    lpgis = (sum(p.sh./p.is+x) + p.p.is.fh.lpg(p.is, p.p.is)).*p.is + 1;
    lpg = [lpg lpgis];
  end

end

function rec = prior_gamma_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.sh)
    rec.sh(ri,:) = p.sh;
  end
  if ~isempty(p.p.is)
    rec.is(ri,:) = p.is;
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [newFluenceMap rowLeafPositions] = mergeFluenceMap(inflMap, rowLeafPositions, resolution);
% JC Sept 18, 2007
% ALG: For each pixel, search the 8 neighbors, merge them, if the change of
% the intensity is less than the threshold.

% General flow:
% Step, group the inter-leaf leakage part, i.e.
%type A: interleaf leakage part, with width 1mm.
%type B: the full leaf width, eg. 0.5cm or 1cm

% How to distinguish them?
% At the first column, check the gradient.
%

% The first colume of inflMap
% Do not use, since the
% s = diff(inflMap(:,1));

% Why s = diff(inflMap(:,2)) is always zero?
% The second colume of inflMap
% s = diff(inflMap(:,2));
% Use first column. Worked for both "bar pattern" and the "complimental bar
% pattern"
s = diff(inflMap(:,1));

newFluenceMap = struct('width', [], 'intensity', []);

newRow = 1;
y=0;

widthLeaf = 5/resolution;       % The width of the MLC leaf is 5mm. (Should work for 10mm=1cm leaf width)
% Havent' tried on resolution = 0.5mm.
iStart = mod(size(inflMap,1), widthLeaf)+1;
newFluenceMap(newRow).width = iStart-1;    %mm
newFluenceMap(newRow).intensity = inflMap(1+1,:);  % second row.

for i = iStart: length(s),
    if (s(i) ~= 0)
        newRow = newRow + 1;
        newFluenceMap(newRow).width = resolution;
        newFluenceMap(newRow).intensity = inflMap(i+1,:);
    else
        if (i == iStart)
            newRow = newRow + 1;
            newFluenceMap(newRow).width = resolution;
        end
        newFluenceMap(newRow).width = newFluenceMap(newRow).width + resolution;
        newFluenceMap(newRow).intensity = inflMap(i+1,:);
    end
end

y = 0;
rowLeafPositions = zeros(size(newFluenceMap,2)+1,1);
rowLeafPositions(1) = y+1;
for i=1:length(newFluenceMap)
    % Following works for resolution == 1.
    % y =y+newFluenceMap(i).width;
    y =y+newFluenceMap(i).width/resolution;
    rowLeafPositions(i+1) = y+1;
end

[Ncol Nrow] = size(inflMap);

% for i = 1 : Nrow,
figure;hAxis1 = axes;imagesc(inflMap);
set(gcf, 'renderer', 'zbuffer')
figure;hAxis2 = axes;hold on;
set(gcf, 'renderer', 'zbuffer')
xL = get(hAxis1, 'xlim');
yL = get(hAxis1, 'ylim');
set(hAxis2, 'xlim', xL);
set(hAxis2, 'ylim', yL);
axis(hAxis2, 'manual');
%     w_colors = floor((w_field ./ max(w_field))*255)+1;
%set(gcf, 'doublebuffer', 'on');

%y=0.05; %(?)
y=0;
for i=1:length(newFluenceMap)
    % Following works for resolution == 1.
    % y =y+newFluenceMap(i).width;
    y =y+newFluenceMap(i).width/resolution;
    for j=1:length(newFluenceMap(i).intensity)
% The display scheme here is problematic. To be worked on.
        patch([xL(1) xL(1) xL(2) xL(2) xL(1)], [y-newFluenceMap(i).width/resolution y y y-newFluenceMap(i).width/resolution y-newFluenceMap(i).width/resolution], newFluenceMap(i).intensity(j));
    end
end
axis([hAxis1 hAxis2], 'ij');
kids = get(hAxis2, 'children');
%     set(kids, 'edgecolor', 'none');
cMap = colormap('jet');
set(hAxis2, 'color', cMap(1,:));

return;


text:
function [CEImage] = runCLAHE(Image,XRes,YRes,Min,Max,NrX,NrY,NrBins,Cliplimit)

%  "Contrast Limited Adaptive Histogram Equalization"
%  by Karel Zuiderveld, karel@cv.ruu.nl
%  in "Graphics Gems IV", Academic Press, 1994
%  (Ported to Matlab by Leslie Smith)
%
%   These functions implement Contrast Limited Adaptive Histogram Equalization.
%   The main routine (CLAHE) expects an input image that is stored contiguously in
%   memory;  the CLAHE output image overwrites the original input image and has the
%   same minimum and maximum values (which must be provided by the user).
%   This implementation assumes that the X- and Y image resolutions are an integer
%   multiple of the X- and Y sizes of the contextual regions. A check on various other
%   error conditions is performed.
%
%
%    Image - The input/output image
%    XRes - Image resolution in the X direction
%    YRes - Image resolution in the Y direction
%    Min - Minimum greyvalue of input image (also becomes minimum of output image)
%    Max - Maximum greyvalue of input image (also becomes maximum of output image)
%    NrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
%    NrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
%    NrBins - Number of greybins for histogram ("dynamic range")
%    Cliplimit - Normalized cliplimit (higher values give more contrast)
%  The number of "effective" greylevels in the output image is set by uiNrBins; selecting
%  a small value (eg. 128) speeds up processing and still produce an output image of
%  good quality. The output image will have the same minimum and maximum value as the input
%  image. A clip limit smaller than 1 results in standard (non-contrast limited) AHE.

% [XRes,YRes]=size(Image);
% CEimage = Image;
CEImage = zeros(XRes,YRes);
if Cliplimit == 1
    return
end

NrBins=max(NrBins,128);
XSize = round(XRes/NrX);
YSize = round(YRes/NrY);
NrPixels = XSize*YSize;
XSize2 = round(XSize/2);
YSize2 = round(YSize/2);

if Cliplimit > 0
    ClipLimit = max(1,Cliplimit*XSize*YSize/NrBins);
else
    ClipLimit = 1E8;
end

LUT=makeLUT(Min,Max,NrBins);
% avgBin = NrPixels/NrBins;
Bin=1+LUT(round(Image));

Hist = makeHistogram(Bin,XSize,YSize,NrX,NrY,NrBins);
if Cliplimit > 0
    Hist = clipHistogram(Hist,NrBins,ClipLimit,NrX,NrY);
end
Map=mapHistogram(Hist,Min,Max,NrBins,NrPixels,NrX,NrY);

% Interpolate
xI = 1;
for i = 1:NrX+1
    if i == 1
        subX = XSize/2;
        xU = 1;
        xB = 1;
    elseif i == NrX+1
        subX = XSize/2;
        xU = NrX;
        xB = NrX;
    else
        subX = XSize;
        xU = i - 1;
        xB = i;
    end
    yI = 1;
    for j = 1:NrY+1
        if j == 1
            subY = YSize/2;
            yL = 1;
            yR = 1;
        elseif j == NrY+1
            subY = YSize/2;
            yL = NrY;
            yR = NrY;
        else
            subY = YSize;
            yL = j - 1;
            yR = j;
        end
        UL = Map(xU,yL,:);
        UR = Map(xU,yR,:);
        BL = Map(xB,yL,:);
        BR = Map(xB,yR,:);
        subImage = Bin(xI:xI+subX-1,yI:yI+subY-1);
        subImage = interpolate(subImage,UL,UR,BL,BR,subX,subY);
        CEImage(xI:xI+subX-1,yI:yI+subY-1) = subImage;
        yI = yI + subY;
    end
    xI = xI + subX;
end


text:
function r = mpower(a,b)
%MPOWER       Hessian power  a ^ n
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  % for scalars use .^ with improved diameter for even exponent
  if ( prod(size(a))==1 ) & ( prod(size(b))==1 )
    r = a .^ b ;
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  if isa(b,'double') & isreal(b) & prod(size(b))==1 & b==round(b)
    if b<0
      error('negative exponent in mpower')
    end
    [m n] = size(a);
    if m~=n
      error('hessian mpower of non-square matrix')
    end
    if b==0
      if issparse(a)
        r = typeadj( hessian(speye(size(a))) , typeof(a) );
      else
        r = typeadj( hessian(eye(size(a))) , typeof(a) );
      end
    else                        % b is integer, at least 1
      b = b - 1;
      r = a;
      while b>0
        if mod(b,2)==1
          r = r*a;
        end
        b = floor(b/2);
        if b~=0
          a = a*a;
        end
      end
    end
  else
    error('invalid call of hessian mpower ^')
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
classdef SOP_F8 < PROBLEM
% <single> <real> <expensive/none>
% Generalized Schwefel's function 2.26

%------------------------------- Reference --------------------------------
% X. Yao, Y. Liu, and G. Lin, Evolutionary programming made faster, IEEE
% Transactions on Evolutionary Computation, 1999, 3(2): 82-102.
%------------------------------- 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 = 1;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D) - 500;
            obj.upper    = zeros(1,obj.D) + 500;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = -sum(PopDec.*sin(sqrt(abs(PopDec))),2);
        end
        %% Generate the minimum objective value
        function R = GetOptimum(obj,N)
            R = -418.9829*obj.D;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [AngleR] = AngleSensorRoomba(serPort);
%[AngleR] = AngleSensorRoomba(serPort)
% Displays the angle in radians and degrees that Create has turned since the angle was last requested.
% Counter-clockwise angles are positive and Clockwise angles are negative.

% By; Joel Esposito, US Naval Academy, 2011

%Initialize preliminary return values
AngleR = nan;

try

%Flush Buffer
N = serPort.BytesAvailable();
while(N~=0)
fread(serPort,N);
N = serPort.BytesAvailable();
end

warning off
global td

fwrite(serPort, [142]);  fwrite(serPort,20);

AngleR = fread(serPort, 1, 'int16')*pi/180;
pause(td)
catch
    disp('WARNING:  function did not terminate correctly.  Output may be unreliable.')
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%

randidx = randperm(size(X, 1));
centroids = X(randidx(1:K), :);

% =============================================================

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function navego_nav2csv(nav)
% navego_nav2csv: exports navigation data to a .csv file.
%
% INPUT
%   nav, INS/GNSS estimations.
%
% OUTPUT
%	.csv file, with the following columns: time (seconds), roll (deg),
%       pitch (deg), yaw (deg), vel N (m/s), vel E (m/s), vel D (m/s),
%        latitude (deg), longitud (deg), altitude (m).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%
% Version: 002
% Date:    2021/12/09
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

R2D = (180/pi);     % Radians to degrees

% File name
if(is_octave)
    dtime = datestr(date, 'yymmdd_HHMMSS' );
else
    dtime = datetime('now','TimeZone','local','Format','yyMMdd-HHmm');
end

file_name = sprintf('navego_%s.csv', dtime);

% Headline
ht = 'time (seconds), ';
ha = 'roll (deg), pitch (deg), yaw (deg), ';
hv = 'vel N (m/s), vel E (m/s), vel D (m/s), ';
hp = 'latitude (deg), longitud (deg), altitude (m),';

% Navigation data
nav_data = [nav.t nav.roll*R2D nav.pitch*R2D nav.yaw*R2D ...
              nav.vel(:,1) nav.vel(:,2) nav.vel(:,3) ...
              nav.lat*R2D nav.lon*R2D nav.h  ]' ;

% Navigation data precision to be saved
pt = '%.6f, ';
pa = repmat('%.3f, ', 1, 3);
pv = repmat('%.3f, ', 1, 3);
pp = repmat('%.7f, ', 1, 2);
ph = '%.3f, ';
format_pattern = [pt pa pv pp ph '\n'];

% Saving the .cvs file
fprintf('navego_nav2cvs: saving navigation data to file %s ... \n', file_name)

fd = fopen(file_name, 'w');

% Saving headline
fprintf(fd, '%s%s%s%s \n', ht, ha, hv, hp);

% Saving data
fprintf(fd, format_pattern , nav_data );

fclose(fd);

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function test_suite=test_map_pca
% tests for cosmo_map_pca
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_map_pca_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_map_pca(varargin{:}),'');

    ds=cosmo_synthetic_dataset();
    pca_params=struct();
    pca_params.mu=[0.0864 0.9248 -0.2001 1.2522 1.0349 0.0568];
    pca_params.coef=[ 0.0014    0.7569    0.3369   -0.0540    0.3234;...
                      0.4052   -0.5164    0.4105   -0.3355   -0.0290;...
                      0.3133    0.0275    0.6618    0.0901    0.1745;...
                     -0.4294   -0.2190   -0.0420   -0.5306    0.6797;...
                      0.7090    0.0436   -0.5070    0.0744    0.4815;...
                     -0.2251   -0.3313    0.1456    0.7677    0.4127];
    pca_params.retain=true(1,5);

    % test mutually exclusive parameters
    aet(ds,'pca_params',pca_params,'pca_explained_count',3);
    aet(ds,'pca_explained_ratio',.5,'pca_explained_count',3);
    aet(ds,'pca_explained_ratio',.5,'pca_params',pca_params);

    % wrong size for coef
    bad_pca_params=pca_params;
    bad_pca_params.coef=1;
    aet(ds,'pca_params',bad_pca_params);

    % missing field
    bad_pca_params=pca_params;
    bad_pca_params=rmfield(bad_pca_params,'mu');
    aet(ds,'pca_params',bad_pca_params);

    % bad explained_count
    aet(ds,'pca_explained_count',-1);
    aet(ds,'pca_explained_count',1.5);
    aet(ds,'pca_explained_count',7);

    % bad explained ratio
    aet(ds,'pca_explained_ratio',-1);
    aet(ds,'pca_explained_ratio',-.001);
    aet(ds,'pca_explained_ratio',1.001);

function test_map_pca_max_feature_count_exceptions()
    default_max_feature_count=1000;

    nsamples=20;
    for max_feature_count=[NaN, ...
                            round(rand()*10+10), ...
                            default_max_feature_count]
        for delta=(-1:1)
            opt=struct();

            if isnan(max_feature_count)
                nfeatures=default_max_feature_count+delta;
            else
                nfeatures=max_feature_count+delta;
                opt.max_feature_count=max_feature_count;
            end

            data=randn(nsamples,nfeatures);

            expect_exception=delta>0;

            for use_struct=[false,true]
                if use_struct
                    ds=struct();
                    ds.samples=data;
                else
                    ds=data;
                end

                func_handle=@()cosmo_map_pca(ds,opt);

                if expect_exception
                    assertExceptionThrown(func_handle);
                else
                    % should not raise an exception
                    func_handle();
                end
            end
        end
    end

function test_map_pca_basics
    nfeatures=ceil(rand()*10+10);
    nsamples_train=ceil(rand()*10+10)+nfeatures;
    nsamples_test=ceil(rand()*10+10)+nfeatures;

    train_samples=randn(nsamples_train,nfeatures);
    test_samples=randn(nsamples_test,nfeatures);

    for count=[1 ceil(nfeatures/2) nfeatures ceil(rand()*nfeatures)]
        helper_test_map_pca_with_count(train_samples,...
                                                test_samples,count);
    end

    for ratio=[.1 .5 1 rand()]
        helper_test_map_pca_with_ratio(train_samples,...
                                                test_samples,ratio);
    end

    helper_test_map_pca_with_no_arguments(train_samples,...
                                                test_samples);

function helper_test_map_pca_with_ratio(train_samples,...
                                                test_samples,ratio)
    ratio_opt=struct();
    ratio_opt.pca_explained_ratio=ratio;

    % check that ratio and count give the same output
    [unused,param]=cosmo_pca(train_samples);
    count=find(cumsum(param.explained)>=ratio*100,1);
    if isempty(count)
        count=numel(param.explained);
    end

    count_opt=struct();
    count_opt.pca_explained_count=count;

    did_throw=true;
    try
        [samples_ratio,params_ratio]=cosmo_map_pca(train_samples,ratio_opt);
        did_throw=false;
    end

    if did_throw
        cosmo_map_pca(train_samples,count_opt)
        assertExceptionThrown(@()...
                cosmo_map_pca(train_samples,count_opt));
        assertExceptionThrown(@()...
                helper_test_dataset_samples_correspondence(train_samples,...
                                            ratio_opt));

    else
        [samples_count,params_count]=cosmo_map_pca(train_samples,...
                                            count_opt);

        assert_almost_equal(samples_ratio,samples_count);
        assert_almost_equal(params_ratio,params_count);

        % verify that raw samples and dataset structures give the same
        % result
        helper_test_dataset_samples_correspondence(train_samples,...
                                                    ratio_opt);
        helper_test_dataset_samples_correspondence(test_samples,...
                                                'pca_params',params_ratio);
    end

    % verify correct result by delegating to count
    if ~did_throw
        helper_test_map_pca_with_count(train_samples,test_samples,count);
    else
        assertExceptionThrown(@()...
            helper_test_map_pca_with_count(train_samples,test_samples,...
                                            count));
    end

function helper_test_map_pca_with_count(train_samples,test_samples,count)
    assert(numel(count)==1);
    opt=struct();
    opt.pca_explained_count=count;

    % do PCA in two ways; we assume that cosmo_pca works correctly, and
    % compare its result to cosmo_map_pca
    [ys,params]=cosmo_map_pca(train_samples,opt);
    [pca_samples,pca_params]=cosmo_pca(train_samples,count);

    assert_almost_equal(ys,pca_samples);

    % verify when using samples
    helper_test_dataset_samples_correspondence(train_samples, ...
                                    'pca_explained_count',count);

    % verify params
    expected_params=struct();
    expected_params.mu=pca_params.mu;
    expected_params.coef=pca_params.coef;
    nfeatures=size(train_samples,2);
    expected_params.retain=(1:nfeatures)<=count;
    assert_almost_equal(expected_params,params);

    [zs,params2]=cosmo_map_pca(test_samples,'pca_params',params);
    assertEqual(params2,params);

    expected_samples=bsxfun(@minus,test_samples,params.mu)*params.coef;
    assert_almost_equal(zs,expected_samples);

    helper_test_dataset_samples_correspondence(test_samples,...
                                'pca_params',params);

function helper_test_map_pca_with_no_arguments(train_samples,test_samples)
    [ys,params]=cosmo_map_pca(train_samples);

    count=size(train_samples,2);
    [ys_count,params_count]=cosmo_map_pca(train_samples,...
                                    'pca_explained_count',count);
    assert_almost_equal(ys,ys_count);
    assert_almost_equal(params,params_count);

    helper_test_dataset_samples_correspondence(train_samples);
    helper_test_dataset_samples_correspondence(test_samples,...
                                            'pca_params',params);

function helper_test_dataset_samples_correspondence(samples, varargin)
    ds=struct();
    ds.samples=samples;

    [result_ds,param_ds]=cosmo_map_pca(ds,varargin{:});
    [result,param]=cosmo_map_pca(samples,varargin{:});

    % verify that output is the same
    expected_ds=struct();
    expected_ds.samples=result;

    expected_ds.fa.comp=1:size(result,2);
    expected_ds.a.fdim.labels={'comp'};
    expected_ds.a.fdim.values={1:size(samples,2)};

    assert_almost_equal(expected_ds,result_ds);

    % parameters must be identical
    assert_almost_equal(param_ds,param);

function assert_almost_equal(x,y,msg,tol)
    % helper that supports structs recursively
    if nargin<4
        tol=1e-5;
    end
    if nargin<3
        msg='';
    end

    if isstruct(x)
        assertTrue(isstruct(y),[msg ' - x is a struct but y is not']);

        fns=fieldnames(x);
        assertEqual(sort(fns),sort(fieldnames(y)),...
                                    [msg ' - fieldname mismatch'])

        n=numel(fns);
        for k=1:n
            fn=fns{k};

            assert_almost_equal(x.(fn),y.(fn),...
                            sprintf('%s - mismatch for field %s',msg,fn))
        end

    elseif iscell(x) && iscell(y)
        assertEqual(size(x),size(y),[msg ' - cell size mismatch']);
        for k=1:numel(x)
            assert_almost_equal(x{k},y{k},...
                        sprintf(' - cell element %d different',k));
        end
    elseif ischar(x) && ischar(y)
        assertEqual(x,y,sprintf('%s - %s ~= %s',msg,x,y));
    elseif islogical(x) && islogical(y)
        assertEqual(x,y,sprintf('%s - boolean array mismatch',msg));
    elseif isnumeric(x) && isnumeric(y)
        if isempty(x) && isempty(y)
            return;
        end
        assertElementsAlmostEqual(x,y,'absolute',tol,[msg ' - different']);
    else
        error('unsupported data type');
    end

function test_pca_map_ratio_unity()
    sizes=[4,10;...
           10,4;...
           4,4];

    for k=size(sizes,1);
        sz=sizes(k,:);
        samples=randn(sz);
        [xs,params]=cosmo_map_pca(samples,'pca_explained_ratio',1);

        assertEqual(params.retain,true(1,sz(end)-1));
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function predData = predictData(rez, samps)

W = gather_try(rez.W);
U = gather_try(rez.U);

samps = samps(1):samps(end);
buff = size(W,1); % width of a spike in samples

predData = zeros(size(U,1),numel(samps)+buff*4);

spikeTimes = rez.st3(:,1); % use the original spikes st2, not the merged and split

inclSpikes = spikeTimes>samps(1)-buff/2 & spikeTimes<samps(end)+buff/2;
st = spikeTimes(inclSpikes); % in samples

inclTemps = uint32(rez.st3(inclSpikes,2)); % use the original spikes st2, not the merged and split
amplitudes = rez.st3(inclSpikes,3);

for s = 1:sum(inclSpikes)
    ampi = amplitudes(s);
    Wi = sq(W(:, inclTemps(s), :));
    Ui = sq(U(:, inclTemps(s), :));

    % this is the reconstruction of the spatial part

    theseSamps = st(s)+(1:buff)-19-samps(1)+buff*2;
    predData(:,theseSamps) = predData(:,theseSamps) + Ui * Wi' * ampi;

    %     predData(:,theseSamps) = predData(:,theseSamps) + ...
    %         squeeze(U(:,inclTemps(s),:)) * squeeze(W(:,inclTemps(s),:))' * amplitudes(s);
end

predData = predData(:,buff*2+1:end-buff*2).*rez.ops.scaleproc;

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function fun = calcS1Density(x,varargin)

if ~check_option(varargin,'sigma')
  % automatic bandwidth selection
  k = ((280*pi^0.5)/729)^0.2/2/pi;
  s = min(std(x), iqr(x)/1.349);
  sigma = k*s*length(x)^-0.2;
else
  sigma = get_option(varargin,'sigma');
end

% compute bandwidth dependent from sigma
N = round(4/sigma);

y = get_option(varargin,'weights',ones(size(x)));
fun = S1FunHarmonic.quadrature(x,y,'bandwidth',N,varargin{:});

% convolution
fun.fhat(1) = 0;
fun.fhat(2:end) = fun.fhat(2:end) .* exp(- 0.5*sigma.^2 * (-N:N).^2).';

% normalize
fun.fhat = fun.fhat ./ fun.fhat(N+2);

end

function test

x = gB.direction.rho;
norm(gB.direction)

out = linspace(-pi,pi,1000);

plot(out,abs(plan2.f))

%plan.fhat
figure(2)
histogram(omega)

% make Gaussian kernel
ti = 1:length(Tp);
mu = length(ti)/2 + 1;
G = (exp( -((ti - mu).^2)./ (2*sigma^2) ))';
% convolve kernel with data
pdf = ifft(fft(Tp).*fft(G));
% remove padding
pdf = pdf(tn-t2+1:tn-t2+tn);
% normalize values
pdf = pdf./sum(pdf);

end


text:
function [hm,HMpix] = euc2hmg(pix)

% EUC2HMG Euclidean to Homogeneous point transform.
%   EUC2HMG(E) transforms the Euclidean point E onto homogeneous space by
%   appending 1 at the last coordinate.
%
%   [h,H_e] = EUC2HMG(E) returns the Jacobian of the transformation.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

hm = [pix;ones(1,size(pix,2))];

if nargout > 1 % Jac -- OK

    if size(pix,2) == 1

        HMpix = [eye(numel(pix));zeros(1,numel(pix))];

    else
        error('??? Jacobians not available for multipla points.')
    end

end

return

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function p = circle_imp_points_3d ( r, pc, nc, n )

%*****************************************************************************80
%
%% CIRCLE_IMP_POINTS_3D returns points on an implicit circle in 3D.
%
%  Discussion:
%
%    Points P on an implicit circle in 3D satisfy the equations:
%
%      ( P(1) - PC(1) )^2
%    + ( P(2) - PC(2) )^2
%    + ( P(3) - PC(3) )^2 = R^2
%
%    and
%
%      ( P(1) - PC(1) ) * NC(1)
%    + ( P(2) - PC(2) ) * NC(2)
%    + ( P(3) - PC(3) ) * NC(3) = 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real PC(3,1), the center of the circle.
%
%    Input, real NC(3,1), a nonzero vector that is normal to
%    the plane of the circle.  It is customary, but not necessary,
%    that this vector have unit norm.
%
%    Input, integer N, the number of points desired.
%    N must be at least 1.
%
%    Output, real P(3,N), the coordinates of points
%    on the circle.
%

%
%  Get two unit vectors N1 and N2 which are orthogonal to each other,
%  and to NC.
%
  [ n1, n2 ] = plane_normal_basis_3d ( pc, nc );
%
%  Rotate R units away from PC in the plane of N1 and N2.
%
  for i = 1 : n

    theta = ( 2.0 * pi * ( i - 1 ) ) / n;

    p(1:3,i) = ( pc(1:3,1) ...
      + r * ( cos ( theta ) * n1(1:3,1) ...
            + sin ( theta ) * n2(1:3,1) ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ADEM_lorenz
% Action-DEM demo specifying an attractor (in terms of the parameters of
% desired equations of motion) This demo first exposes an agent to a Lorenz
% attractor world such that it can learn the three parameters of the Lorenz
% system. It is then placed in a test world with a fixed point attractor to
% see if it has remembered the chaotic dynamics it learnt in the training
% environment
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_lorenz.m 4804 2012-07-26 13:14:18Z karl $

% generative model
%==========================================================================
G(1).E.s = 1/2;                        % smoothness
G(1).E.n = 4;                          % smoothness
G(1).E.d = 2;                          % smoothness
G(1).E.linear = 3;                     % nonlinear

% dynamics
%--------------------------------------------------------------------------
f0      = '[v; 0; 0] + [-1 1 1; -1/2 -1 0; 0 1 -1]*x/P';
fA      = '[v; 0; 0] + [-1 1 1; -1/2 -1 0; 0 1 -1]*x/P + a';
fL      = '[v; 0; 0] + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x/64';

% parameters
%--------------------------------------------------------------------------
P0      = 16;
PL      = [10; -8/3; 32];
PM      = [0;     0;  0];

% P(1): Prandtl number
% P(2): 8/3
% P(3): Rayleigh number

% level 1
%--------------------------------------------------------------------------
G(1).x  = [1; 1; 24];
G(1).f  = inline(fL ,'x','v','a','P');
G(1).g  = inline('x','x','v','a','P');
G(1).pE = PL;
G(1).V  = exp(8);                           % error precision
G(1).W  = exp(8);                           % error precision

% level 2
%--------------------------------------------------------------------------
G(2).a  = [0;0;0];                          % action variables
G(2).v  = 0;                                % inputs
G(2).V  = exp(16);
G       = spm_ADEM_M_set(G);

% plot flow fields and equilibrium densities
%==========================================================================
spm_figure('GetWin','Figure 1');

x{1}    = linspace(-20,20,32);
x{2}    = linspace(-32,32,32);
x{3}    = linspace(  10,40,8);

% controlled flow (P0)
%--------------------------------------------------------------------------
subplot(3,2,1)
spm_fp_display_density(G,x);
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('controlled','Fontsize',16)
a     = axis;

% exemplar trajectory
%--------------------------------------------------------------------------
U.u   = sparse(1024,G(1).m);
t     = spm_int_J(G(1).pE,G,U);
subplot(3,2,2)
plot(t(:,1),t(:,2),'r')
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('controlled','Fontsize',16)
axis(a);

% uncontrolled flow (P0)
%--------------------------------------------------------------------------
G0       = G;
G0(1).f  = inline(f0,'x','v','P');
G0(1).pE = P0;

subplot(3,2,3)
spm_fp_display_density(G0,x);
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('before','Fontsize',16)
a = axis;

% exemplar trajectory
%--------------------------------------------------------------------------
t     = spm_int_J(G0(1).pE,G0,U);
subplot(3,2,4)
plot(t(:,1),t(:,2),'r')
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('before','Fontsize',16)
axis(a);
drawnow

% recognition model: learn the controlled environmental dynamics
%==========================================================================

% make a naive model (M)
%--------------------------------------------------------------------------
M       = G;
M(1).f  = inline(fL ,'x','v','P');
M(1).g  = inline('x','x','v','P');
M(1).pE = PM;
M(1).pC = 32;

% teach naive model by exposing it to a controlled environment (G)
%--------------------------------------------------------------------------
n     = 256;
C     = sparse(n,1);

DEM.M = M;
DEM.G = G;
DEM.C = C;
DEM.U = C;

% optimise recognition model and show results of learning
%--------------------------------------------------------------------------
DEM    = spm_ADEM(DEM);

spm_figure('GetWin','Figure 2');
spm_DEM_qP(DEM.qP,DEM.pP)

% replace priors with learned conditional expectation and plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');

M       = DEM.M;
M(1).pE = DEM.qP.P{1};

subplot(3,2,5)
spm_fp_display_density(M,x);
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('after','Fontsize',16)
a   = axis;

% exemplar trajectory
%--------------------------------------------------------------------------
t   = spm_int_J(M(1).pE,M,U);
subplot(3,2,6)
plot(t(:,1),t(:,2),'r')
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('after','Fontsize',16)
axis(a);
drawnow

% evaluate performance under active inference
%==========================================================================
N       = 256;
C       = spm_conv(randn(1,N)*4,8);

% create uncontrolled environment: A (with action)
%--------------------------------------------------------------------------
G(1).f  = inline(f0 ,'x','v','a','P');       % no action
G(1).f  = inline(fA ,'x','v','a','P');       % with action
G(1).pE = P0;

% make the recognition model confident about its predictions
%--------------------------------------------------------------------------
M(1).x  = G(1).x;
M(1).V  = exp(8);
M(1).W  = exp(8);
M(1).pC = [];

G(1).W  = exp(16);

% create DEM structure
%--------------------------------------------------------------------------
DEM.G   = G;
DEM.M   = M;
DEM.C   = sparse(1,N);                       % no perturbation
DEM.C   = C;                                 % with perturbation
DEM.U   = sparse(1,N);
DEM     = spm_ADEM(DEM);

spm_DEM_qU(DEM.qU,DEM.pU)

% plot behaviour
%--------------------------------------------------------------------------
subplot(2,2,3)
spm_fp_display_density(M,x); hold on
plot(DEM.pU.v{1}(1,:),DEM.pU.v{1}(2,:),'b'), hold off
xlabel('position','Fontsize',12)
ylabel('velocity','Fontsize',12)
title('learnt','Fontsize',16)

return

% Notes: alternative specification of Lorenz system
%==========================================================================

% dynamics
%--------------------------------------------------------------------------
fL      = '[v; 0; 0] + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x/32';

% parameters
%--------------------------------------------------------------------------
PL      = [10; -8/3; 32];
P.A     = [-PL(1) PL(1) 0; PL(3) -1 0; 0 0 PL(2)]/32;
P.B     = [0 0  0;
    0 0 -1;
    0 1  0]/32;

% PL(1): Prandtl number
% PL(2): 8/3
% PL(3): Rayleigh number

author: spm


text:
function a = le( x, y )

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '<=' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function obj=set_limit_extra(obj,x_extra,y_extra,z_extra)
%set_limit_extra Add some breathing room around data in plots
%
% Example syntax gramm_object.set_limit_extra([0.05 0.2],0.05)
% first argument is XLim extra, second is YLim extra, extra
% room is expressed as ratio to original limits. It is possible to set
% extra room separately for upper and lower limits by giving a 2x1 array.
% If not, the lower and uupper limit extra will be the same

if nargin<4
    z_extra=0;
end

%If only one element is given we use it for upper and lower extra
if numel(x_extra==1)
    x_extra=[x_extra x_extra];
end
if numel(y_extra==1)
    y_extra=[y_extra y_extra];
end
if numel(z_extra==1)
    z_extra=[z_extra z_extra];
end

for obj_ind=1:numel(obj)
    obj(obj_ind).xlim_extra=x_extra;
    obj(obj_ind).ylim_extra=y_extra;
    obj(obj_ind).zlim_extra=z_extra;
end

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
%   FUZZY RULE BASE COMPRESSION SOFTWARE
%
%   This software implements advanced inference in fuzzy rule based systems
%   by filtration of non-monotonic rules. This software is for mandani type
%   system only.
%
%
%   Neelamugilan Gobalakrishnan
%   2006
%   $Revision:1.0 $ $Date: 15/03/2006 $

clc
clear all

disp('***********************************************************');
disp('*                                                         *');
disp('*         FUZZY RULE BASE COMPRESSION SOFTWARE            *');
disp('*                                                         *');
disp('*     This software implements advanced inference in      *');
disp('*     fuzzy rule based systems by filtration of non-      *');
disp('*     monotonic rules.                                    *');
disp('*                                                         *');
disp('***********************************************************');

fprintf('\n\n');

condStatus = 'start'; %Condition status for choice while loop

while strcmp(condStatus,'start') | strcmp(condStatus,'cont')
    clear choise
    fprintf('Enter 1 to create a new system\n');
    fprintf('Enter 2 to load a system (*.fis)\n\n');
    fprintf('Please enter your choice : ');
    choise = input('');
    if choise == 1 % Create a new system
        fprintf('\n');
        systemName = input('Enter the name of this system : ','s');
        % Call createinttable function to create integer table
        [sys intTable] = createinttable(systemName);
        % Arranged the rules into groups
        [groupedTable outNum] = arrangerules(sys,intTable);

        % Get the initial values for each inputs
        for i=1:size(sys.input, 2)
            fprintf('\nEnter the initial %s : ', sys.input(i).name);
            initialVal(i) = input('');
        end

        displayOpt = 'on'; % Option to display the compressed rules
        % Find dominant rules for provided initial values
        compressRules = finddomrules(sys,intTable,outNum,initialVal,displayOpt);

        newSysName = input('Please enter the name for compressed system : ','s');
        fprintf('\n');
        % Create the system with compressed rules
        createfis(newSysName,sys,compressRules);

        % Simulates the fuzzy inference system
        a = readfis(newSysName);
        outResult = evalfis([initialVal],a);

        fprintf('\n');
        fprintf('The %s value for \n', sys.output(outNum).name);

        for i=1:size(sys.input, 2)
            fprintf('%s(%.2f) \n',sys.input(i).name ,initialVal(i));
        end

        fprintf('= %.2f', outResult);
        fprintf('\n\n');

        view_surface = input('Do u want to view the output surface (yes/no) : ','s');

        % Display the output surface
        while strcmp(view_surface,'yes')
            fprintf('\nEnter 1 to view the original system, output surface\n');
            fprintf('Enter 2 to view the compressed system, output surface\n\n');
            fprintf('Please enter your option : ');
            option = input('');
            if option == 1
                % Create output surface for the original system
                createoutsurf(sys,systemName,intTable,outNum,option);
                fprintf('\nDo u want to view the output surface with\n');
                view_surface = input('new points or compressed system? (yes/no) : ','s');
            elseif option == 2
                % Create output surface for the compressed system
                createoutsurf(sys,newSysName,intTable,outNum,option);
                fprintf('\nDo u want to view the output surface with\n');
                view_surface = input('new points or original system? (yes/no) : ','s');
            else
                disp('Invalid value has been entered! ');
                fprintf('\n');
                view_surface = 'yes';
            end
        end
        condStatus = 'stop'; % choice while loop stopping condition
        % Load a system
    elseif choise == 2
        [fileName,pathName]=uigetfile('*.fis','Read FIS');
        if isequal(fileName,0) || isequal(pathName,0)
          % If fileName is zero, "cancel" was hit, or there was an error.
          errorStr='No file was loaded';
          return
        end
        systemName = [pathName fileName];

        sys = readfis(systemName);

        % Create the integer table
        for i=1:size(sys.rule,2)
            intTable(i,:) = [sys.rule(i).antecedent sys.rule(i).consequent];
        end

        fprintf('\n');
        disp('-------------------------------------------------------------');
        fprintf(' ');
        disp('                    Complete integer table                   ');
        fprintf(' ');
        disp('-------------------------------------------------------------');

        fprintf('\n');
        disp('-------------------------------------------------------------');
        fprintf(' ');

        for i=1:size(sys.input,2)
            fprintf('Input%d ', i);
        end

        for i=1:size(sys.output,2)
            fprintf('Output%d ', i);
        end

        fprintf('\n');
        disp('-------------------------------------------------------------');
        fprintf('\n');
        disp(intTable);
        disp('-------------------------------------------------------------');
        fprintf('\n');

        % Arranged the rules into groups
        [groupedTable outNum] = arrangerules(sys,intTable);

        % Get the initial values for each inputs
        for i=1:size(sys.input, 2)
            fprintf('\nEnter the initial %s : ', sys.input(i).name);
            initialVal(i) = input('');
        end

        displayOpt = 'on'; % Option to display the compressed rules
        % Find dominant rules for provided initial values
        compressRules = finddomrules(sys,intTable,outNum,initialVal,displayOpt);

        fprintf('\n');
        newSysName = input('Please enter the name for compressed system : ','s');
        % Create the system with compressed rules
        createfis(newSysName,sys,compressRules);

        % Simulates the fuzzy inference system
        a = readfis(newSysName);
        outResult = evalfis([initialVal],a);

        fprintf('\n');
        fprintf('The %s value for \n', sys.output(outNum).name);

        for i=1:size(sys.input, 2)
            fprintf('%s(%.2f) \n',sys.input(i).name ,initialVal(i));
        end

        fprintf('= %.2f', outResult);
        fprintf('\n\n');

        view_surface = input('Do u want to view the output surface (yes/no) : ','s');

        % Display the output surface
        while strcmp(view_surface,'yes')
            fprintf('\nEnter 1 to view the original system, output surface\n');
            fprintf('Enter 2 to view the compressed system, output surface\n\n');
            fprintf('Please enter your option : ');
            option = input('');
            if option == 1
                % Create output surface for the original system
                createoutsurf(sys,systemName,groupedTable,outNum,option);
                fprintf('\nDo u want to view the output surface with\n');
                view_surface = input('new points or compressed system? (yes/no) : ','s');
            elseif option == 2
                % Create output surface for the compressed system
                createoutsurf(sys,newSysName,groupedTable,outNum,option);
                fprintf('\nDo u want to view the output surface with\n');
                view_surface = input('new points or original system? (yes/no) : ','s');
            else
                fprintf('\n');
                disp('Invalid value has been entered! ');
                view_surface = 'yes';
            end
        end
        condStatus = 'stop'; % choice while loop stopping condition
    else
        disp('Invalid value has been entered! ');
        condStatus = 'cont'; % Choice while loop continue condition
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function p_data = p02_data ( dim_num, data_num )

%*****************************************************************************80
%
%% P02_DATA returns the data for problem p02.
%
%  Discussion:
%
%    Two pairs of identical X values have now been slightly separated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the dependent
%    variables.
%
%    Input, integer DATA_NUM, the number of data points.
%
%    Output, real P_DATA(DIM_NUM,DATA_NUM), the data.
%
  p_data = [ ...
     0.00,   0.00; ...
     1.34,   5.00; ...
     5.00,   8.66; ...
    10.00,  10.00; ...
    10.60,  10.40; ...
    10.70,  12.00; ...
    10.705, 28.60; ...
    10.80,  30.20; ...
    11.40,  30.60; ...
    19.60,  30.60; ...
    20.20,  30.20; ...
    20.295, 28.60; ...
    20.30,  12.00; ...
    20.40,  10.40; ...
    21.00,  10.00; ...
    26.00,   8.66; ...
    29.66,   5.00; ...
    31.00,   0.00  ]';

  return
end