text:
function SieveAnalysis5(Granulometria)
% _____________________________________________________
% See Contents.m
%  With this routine you can get the sand sieve analysis. Routine
% displays the computation of the particle-size distribution in different
% windows. Add you can get the main parameters of statistics analysis:
% percentiles, mean, standard deviation, kurtosis, etc.
% Syntax:
%              >> Granulometria[ %your data];
%              >> SieveAnalysis5(Granulometria)
% Update: 1. I've received a couples of emails, in which  some
%                people had a lot of problems with how to give the inputs;
%                for this reason, I decided to delete the GUI. Now, the
%                user only has to do, it's to create a variable with the
%                name "Granulometria" and inside this variable, you need to
%                write your data. First column, the sieve mesh (in mm) and
%                the second, the weights (in g). For example:
%                Granulometria = [4.0 0.88; 2.0 1.08; 1.0 1.49; 0.50 3.58; 0.25 11.50;
%                                                0.125 21.50; 0.0625 9.81; 0.0313 2.70];
%                2. I deleted an enormous bugs.
%
% Last modification : 03/21/09
% Author: PhD(c) Gabriel Ruiz
% This routine is provided "as is" without warranty of any kind.
% Please, you don't attribute it.
% If you'll detect any mistake or bug, may you communicate to me,
% please?
% v1.1.5
% Copyright 2006.
% _____________________________________________________
     clc;

    screen = get(0, 'screensize');
    if screen(1,3) == 1152 && screen(1,4) == 864
                 screen(1,3) = screen(1,3) - 128; screen(1,4) = screen(1,4) - 96;
        elseif screen(1,3) == 1280 && screen(1,4) == 768
                 screen(1,3) = screen(1,3) - 256;
        elseif screen(1,3) == 1280 && screen(1,4) == 800
                 screen(1,3) = screen(1,3) - 256; screen(1,4) = screen(1,4) - 32;
        elseif screen(1,3) == 1280 && screen(1,4) == 960
                 screen(1,3) = screen(1,3) - 256; screen(1,4) = screen(1,4) - 192;
        elseif screen(1,3) == 1280 && screen(1,4) == 1024
                 screen(1,3) = screen(1,3) - 256; screen(1,4) = screen(1,4) - 256;
        elseif screen(1,3) == 1400 && screen(1,4) == 1050
                 screen(1,3) = screen(1,3) - 376; screen(1,4) = screen(1,4) - 282;
        elseif screen(1,3) == 1600 && screen(1,4) == 900
                 screen(1,3) = screen(1,3) - 576; screen(1,4) = screen(1,4) - 132;
        elseif screen(1,3) == 1600 && screen(1,4) == 1200
                 screen(1,3) = screen(1,3) - 576; screen(1,4) = screen(1,4) - 432;
    end
    vpan = screen(1,4)-35; hpan = screen(1,3);
    colores = [1 1 0.85 ]; colortexto = [0.50 0.25 0.25];
    WB = Granulometria(:,2);
    milimeters = Granulometria(:,1);
    r = length(milimeters);

        for i = 1 : r
                Worksheet_sieve_data{1,1}(i) = WB(i);
        end
        Sum_Col_1 = sum( Worksheet_sieve_data{1,1} );
        for i = 1: r
                Worksheet_sieve_data{1,2}(i) = ( Worksheet_sieve_data{1,1}(i) / Sum_Col_1) * 100;
                if i == 1
                        Worksheet_sieve_data{1,3}(i) = Worksheet_sieve_data{1,2}(i);
                        Worksheet_sieve_data{1,4}(i) = 100;
                else
                        Worksheet_sieve_data{1,3}(i) = Worksheet_sieve_data{1,2}(i) + Worksheet_sieve_data{1,3}(i-1);
                        Worksheet_sieve_data{1,4}(i) = Worksheet_sieve_data{1,4}(1) - Worksheet_sieve_data{1,3}(i);
                end
        end

        diamX = Worksheet_sieve_data{1,4} ;
        diamX = diamX';
        diamY = milimeters;
        D_5 = interp1(diamX, diamY, 5, 'pchip');
        D_10 = interp1(diamX, diamY, 10, 'pchip');
        D_16 = interp1(diamX, diamY, 16, 'pchip');
        D_25 = interp1(diamX, diamY, 25, 'pchip');
        D_30 = interp1(diamX, diamY, 30, 'pchip');
        D_50 = interp1(diamX, diamY, 50, 'pchip');
        D_60 = interp1(diamX, diamY, 60, 'pchip');
        D_75 = interp1(diamX, diamY, 75, 'pchip');
        D_84 = interp1(diamX, diamY, 84, 'pchip');
        D_95 = interp1(diamX, diamY, 95, 'pchip');

        resulchar.d_5 = num2str(D_5);
        resulchar.d_10 = num2str(D_10);
        resulchar.d_16 = num2str(D_16);
        resulchar.d_25 = num2str(D_25);
        resulchar.d_30 = num2str(D_30);
        resulchar.d_50 = num2str(D_50);
        resulchar.d_60 = num2str(D_60);
        resulchar.d_75 = num2str(D_75);
        resulchar.d_84 = num2str(D_84);
        resulchar.d_95 = num2str(D_95);

        D_5phi = -log(D_5) / log(2);      %Matlab has the function log2, but I prefer to use this form.
        D_10phi = -log(D_10) / log(2);
        D_16phi = -log(D_16) / log(2);
        D_25phi = -log(D_25) / log(2);
        D_30phi = -log(D_30) / log(2);
        D_50phi = -log(D_50) / log(2);
        D_60phi = -log(D_60) / log(2);
        D_75phi = -log(D_75) / log(2);
        D_84phi = -log(D_84) / log(2);
        D_95phi = -log(D_95) / log(2);

        Mean_grain_size = ( D_16phi + D_50phi + D_84phi ) / 3;
        Mean_grain_size_mm = 2 ^ -(Mean_grain_size);
        resulchar.MeanGS = num2str(Mean_grain_size_mm);

        if Mean_grain_size_mm <= 0.075
                    resulchar.ASTM = 'Fine Soil';
            elseif Mean_grain_size_mm >= 0.076 && Mean_grain_size_mm <= 0.425
                    resulchar.ASTM = 'Fine sand';
            elseif Mean_grain_size_mm >= 0.426 && Mean_grain_size_mm <= 2
                    resulchar.ASTM = 'Medium sand';
            elseif Mean_grain_size_mm >= 2.1 && Mean_grain_size_mm <= 4.75
                    resulchar.ASTM = 'Coarse sand';
            elseif Mean_grain_size_mm >= 4.76 && Mean_grain_size_mm <= 19
                    resulchar.ASTM = 'Fine gravel';
            elseif Mean_grain_size_mm >= 19.1 && Mean_grain_size_mm <= 75
                    resulchar.ASTM = 'Coarse gravel';
        end

        if Mean_grain_size_mm >= 0.0625 && Mean_grain_size_mm <= 0.125
                    resulchar.Wentworth = 'Very fine sand';
            elseif Mean_grain_size_mm >= 0.126 && Mean_grain_size_mm <= 0.250
                    resulchar.Wentworth = 'Fine sand';
            elseif Mean_grain_size_mm >= 0.251 && Mean_grain_size_mm <= 0.50
                    resulchar.Wentworth = 'Medium sand';
            elseif Mean_grain_size_mm >= 0.51 && Mean_grain_size_mm <= 1
                    resulchar.Wentworth = 'Coarse sand';
            elseif Mean_grain_size_mm >= 1.01 && Mean_grain_size_mm <= 2
                    resulchar.Wentworth = 'Very coarse sand';
            elseif Mean_grain_size_mm >= 2.01 && Mean_grain_size_mm <= 4.76
                    resulchar.Wentworth = 'Granule';
            elseif Mean_grain_size_mm >= 4.77 && Mean_grain_size_mm <= 8
                    resulchar.Wentworth = 'Small pebble';
            elseif Mean_grain_size_mm >= 8.01 && Mean_grain_size_mm <= 16
                    resulchar.Wentworth = 'Medium pebble';
            elseif Mean_grain_size_mm >= 16.01 && Mean_grain_size_mm <= 19.03
                    resulchar.Wentworth = 'Large pebble';
        end

        gravel = zeros(r,1);
        for i = 1 : r
            if milimeters(i,1) <= 75 && milimeters(i,1) >= 4.75
                gravel(i,1) = Worksheet_sieve_data{1,2}(i);
            end
        end
        Gravel = sum(gravel(:,1));
        resulchar.Grava = num2str(Gravel);

        sand = zeros(r,1);
        for i = 1 : r
            if milimeters(i,1) <= 4.74 && milimeters(i,1) >= 0.075
                    sand(i,1) = Worksheet_sieve_data{1,2}(i);
            end
        end
        Sand = sum(sand(:,1));
        resulchar.arena = num2str(Sand);

         fine = zeros(r,1);
        for i = 1 : r
            if milimeters(i,1) <= 0.074
                    fine(i,1) = Worksheet_sieve_data{1,2}(i);
            end
        end
        Fine = sum(fine(:,1));
        resulchar.fino =num2str(Fine);
        Muestrapor = [ Gravel Sand Fine ];

        Standard_Deviation = ( ( D_84phi - D_16phi ) / 4 ) + ( ( D_95phi - D_5phi ) / 6 );   % E.3.
        Standard_Deviation_mm = 2 ^ -(Standard_Deviation);
        resulchar.SDe = num2str(Standard_Deviation_mm);
        if Standard_Deviation < 0.34
                    resulchar.SD = 'Very well sorted';
            elseif Standard_Deviation >= 0.35 && Standard_Deviation <= 0.49
                    resulchar.SD = 'Well sorted';
            elseif Standard_Deviation >= 0.50 && Standard_Deviation <= 0.71
                    resulchar.SD = 'Moderately well sorted';
            elseif Standard_Deviation >= 0.72 && Standard_Deviation <= 0.99
                    resulchar.SD = 'Moderately sorted';
            elseif Standard_Deviation >= 1.00 && Standard_Deviation <= 1.99
                    resulchar.SD = 'Poorly sorted';
            elseif Standard_Deviation >= 2.00 && Standard_Deviation <= 3.99
                    resulchar.SD = 'Very poorly sorted';
            elseif Standard_Deviation >= 4.00
                    resulchar.SD = 'Extremely poorly sorted';
        end

        Skewness =  ( ( D_84phi + D_16phi - ( 2 * D_50phi ) ) / ( 2 * ( D_84phi - D_16phi ) ) ) + ...
                              ( ( D_95phi + D_5phi - ( 2 * D_50phi ) ) / ( 2 * ( D_95phi - D_5phi ) ) );
        Skewness_mm = 2 ^ -(Skewness);
        resulchar.S_ke = num2str(Skewness_mm);
        if Skewness < -0.29
                    resulchar.Sk = 'Very coarse skewed';
            elseif Skewness >= -0.30 && Skewness <= -0.09
                    resulchar.Sk = 'Coarse skewed';
            elseif Skewness >= -0.10 && Skewness <= 0.09
                    resulchar.Sk = 'Near symmetrical';
            elseif Skewness >= 0.10 && Skewness <= 0.29
                    resulchar.Sk = 'Fine skewed';
            elseif Skewness >= 0.30
                    resulchar.Sk = 'Very fine skewed';
        end

        Kurtosis = ( D_95phi - D_5phi ) / ( 2.44 * ( D_75phi - D_25phi ) );
        Kurtosis_mm = 2 ^ -(Kurtosis);
        resulchar.K_ur = num2str(Kurtosis_mm);
         if Kurtosis < 0.64
                    resulchar.kurt = 'Very platykurtic (flat)';
            elseif Kurtosis >= 0.65 && Kurtosis <= 0.89
                    resulchar.kurt = 'Platykurtic';
            elseif Kurtosis >= 0.90 && Kurtosis <= 1.10
                    resulchar.kurt = 'Mesokurtic ';
            elseif Kurtosis >= 1.11 && Kurtosis <= 1.49
                    resulchar.kurt = 'Leptokurtic (peaked)';
            elseif Kurtosis >= 1.50 && Kurtosis <= 2.99
                    resulchar.kurt = 'Very leptokurtic';
            elseif Kurtosis >= 3.00
                    resulchar.kurt = 'Extremely leptokurtic';
        end

        Cc_Dispersal = D_30 ^ 2 / ( D_10 * D_60 ) ;
        resulchar.CC = num2str(Cc_Dispersal);
        Cu_Hazen_uniformity = D_60 / D_10;
        resulchar.Cu =num2str(Cu_Hazen_uniformity);
        if Cu_Hazen_uniformity >= 6 && ( ( Cc_Dispersal >= 1 ) && ( Cc_Dispersal <= 3 ) )
                    resulchar.SUCS = 'SW';
        else
                    resulchar.SUCS = 'SP';
        end

        figure(3);
        set(figure(3), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , ...
                'Name' , '1.Sieve Analysis - Histogram' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [  hpan-880   vpan-440   hpan-664   vpan-464 ] );
        histo_h =bar(milimeters,Worksheet_sieve_data{1,2}, 'group');
        colormap hsv;
        xlabel('Grain Size in mm' , 'Fontsize' , 7);
        ylabel('Percent Finer by Weight (%)' , 'Fontsize' , 7);
        title('Histogram' , 'Color' , 'r' , 'FontWeight' , 'Bold' , 'Fontsize' , 8);
        set(histo_h,  'FaceColor', 'r');
        set(gca,  'Xgrid' , 'on' , 'Ygrid' , 'on' , 'Xcolor' , [0 0 0.37] , 'Ycolor' , [0 0 0.37] , 'Box' , 'off' , ...
                   'FontWeight' , 'Bold' , 'FontSize' , 7 , 'FontAngle' , 'Oblique' , 'XMinorTick' , 'on' , 'YMinorTick' , 'on');

        figure(4);
        set(figure(4), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [  hpan-859   vpan-469  hpan-664  vpan-464 ] , ...
                'name' , '2.Sieve Analysis - Frecuency Curve' );
        FCruve_h = plot(milimeters, Worksheet_sieve_data{1,2}, '-', 'Color' , 'b' , 'Linewidth' , 1.5);
        xlabel('Grain Size in mm' , 'Fontsize' , 7);
        ylabel('Percent Finer by Weight (%)' , 'Fontsize' , 7);
        title('Frecuency Curve' , 'Color' , 'b' , 'FontWeight' , 'Bold' , 'Fontsize' , 8);
        set(gca,  'Xgrid' , 'on' , 'Ygrid' , 'on' , 'Xcolor' , [0 0 0.37] , 'Ycolor' , [0 0 0.37] , 'Box' , 'off' , ...
                   'FontWeight' , 'Bold' , 'FontSize' , 7 , 'FontAngle' , 'Oblique' , 'XMinorTick' , 'on' , 'YMinorTick' , 'on');

        figure(5);
        set(figure(5), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [hpan-838  vpan-499  hpan-664  vpan-464] , ...
                'name' , '3. Sieve Analysis - Cumulative Arithmetic Curve' );
        ACurve_h = plot(milimeters, Worksheet_sieve_data{1,3}, '-', 'Color' , [0 0.5 0] , 'Linewidth' , 1.5);
        xlabel('Grain Size in mm' , 'Fontsize' , 7);
        ylabel('Percent Finer by Weight (%)' , 'Fontsize' , 7);
        title('Frecuency Cumulative  Curve' , 'Color' , [0 0.5 0] , 'FontWeight' , 'Bold' , 'Fontsize' , 8);
        set(gca,  'Xgrid' , 'on' , 'Ygrid' , 'on' , 'Xcolor' , [0 0 0.37] , 'Ycolor' , [0 0 0.37] , 'Box' , 'off' , ...
                   'FontWeight' , 'Bold' , 'FontSize' , 7 , 'FontAngle' , 'Oblique' , 'XMinorTick' , 'on' , 'YMinorTick' , 'on');

        figure(6);
        set(figure(6), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [ hpan-816   vpan-529   hpan-664  vpan-464] , ...
                'name' , '4. Sieve Analysis - Cumulative Probability Arithmetic Curve' );
        CCurve_h = plot(milimeters, Worksheet_sieve_data{1,4}, '-','Color' , [1 0.5 0] , 'Linewidth' , 1.5);
        xlabel('Grain Size in mm', 'Fontsize' , 7);
        ylabel('Probability Finer by Weight (%)' , 'Fontsize' , 7);
        title('Cumulative Probability Curve ' , 'Color' , [1 0.5 0] , 'FontWeight' , 'Bold' , 'Fontsize' , 8);
        set(gca,  'Xgrid' , 'on' , 'Ygrid' , 'on' , 'Xcolor' , [0 0 0.37] , 'Ycolor' , [0 0 0.37] , 'Box' , 'off' , ...
                   'FontWeight' , 'Bold' , 'FontSize' , 7 , 'FontAngle' , 'Oblique' , 'XMinorTick' , 'on' , 'YMinorTick' , 'on');

        figure(7);
        set(figure(7), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [   hpan-778   vpan-587    hpan-348    vpan-464  ] , ...
                 'name' , '6. Sieve Analysis - Cumulative Probability Semi-Log Curve'  );
        CCurvesl_h = semilogx(milimeters, Worksheet_sieve_data{1,4}, '-o','Color' , [0 0 0.5] , ...
                                      'Linewidth' , 2.5, 'MarkerFaceColor' , [ 0 0 0.5 ] , 'MarkerSize' , 2, ...
                                      'MarkerEdgeColor' , [ 0 0.5 0.5] );
        axis tight;
        xlabel('Grain Size in mm', 'Fontsize' , 7);
        ylabel('Probability Finer by Weight (%)' , 'Fontsize' , 7);
        title('Cumulative Probability Curve ' , 'Color' , [0 0 0.5] , 'FontWeight' , 'Bold' , 'Fontsize' , 8);
        set(gca,  'Xgrid' , 'on' , 'Ygrid' , 'on' , 'Xcolor' , [0 0 0.37] , 'Ycolor' , [0 0 0.37] , 'Box' , 'off' , 'XDir' , 'normal' , ...
                   'FontWeight' , 'Bold' , 'FontSize' , 7 , 'FontAngle' , 'Oblique' , 'XMinorTick' , 'on' , 'YMinorTick' , 'on', ...
                   'Xlim' , [0.01 10]  , 'XTickLabel', { '0.01' ; '0.1' ; '1' ;  '10' } );
        dimark = [0.074 4.76];
        hold on;
        escalY = get(gca,'YLim');
        escalY = [ 0 escalY(1,2) ];
        lenmark = length(dimark);
        xmark = reshape([dimark;dimark;ones(1,lenmark)*nan;],lenmark*3,1);
        ymark = repmat([escalY nan],1,lenmark);
        plot(xmark,ymark, 'r', 'LineWidth', 2 );
        punto07 = patch(D_50 , 50 , [1 0 0] );
        punto08 = patch(0.01, 50 ,  [1 0 0] );
        punto09 = patch(D_50 , 0.01 , [1 0 0]);
         l4 = line('Xdata' , [D_50 D_50] , 'Ydata' , [0 50] , 'Color' ,  [0.75 0.75 0.75] , ...
                        'Linestyle' , '-' , 'Linewidth' , 2 ) ;
         l5 = line('Xdata' , [0.01 D_50] , 'Ydata' , [50 50] , 'Color' ,  [0.75 0.75 0.75] , ...
                        'Linestyle' , '-' , 'Linewidth' , 2 );
        set(punto07, 'marker' , 'o' , 'MarkerEdgeColor', 'r' , 'MarkerFaceColor', 'r' , ...
                 'MarkerSize' , 2);
        set(punto08, 'marker' , 'o' , 'MarkerEdgeColor', 'r' , 'MarkerFaceColor', 'r' , ...
                 'MarkerSize' , 2);
        set(punto09, 'marker' , 'o' , 'MarkerEdgeColor', 'r' , 'MarkerFaceColor', 'r', ...
                 'MarkerSize' , 2);
        td50 = text(D_50+.005 , 3, num2str(D_50)); set(td50, 'FontSize' , 8 , ...
                            'FontAngle' , 'Normal' , 'FontWeight' , 'Bold' );
        tfinos = text(0.011 , 95, 'Fine'); set(tfinos, 'FontSize' , 8 , 'FontWeight' , 'Bold');
        tsand = text(0.085 , 95, 'Sand'); set(tsand, 'FontSize' , 8 , 'FontWeight' , 'Bold' );
        tgrav = text(5 , 95, 'Gravel'); set(tgrav, 'FontSize' , 8 , 'FontWeight' , 'Bold' );

        figure(8);
        set(figure(8), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [   hpan-796   vpan-558    hpan-664  vpan-464  ] , ...
                'name' , '5. Sieve Analysis - Cumulative Probability Logarithmic Curve' );
        CCurvel_h = loglog(milimeters, Worksheet_sieve_data{1,4}, '-','Color' , [0.5 0 0.5] , 'Linewidth' , 1.5);
        xlabel('Grain Size in mm' ,'Fontsize' , 7);
        ylabel('Probability Finer by Weight (%)', 'Fontsize' , 7);
        title('Cumulative Probability Curve ' , 'Color' , [0.5 0 0.5] , 'FontWeight' , 'Bold' , 'Fontsize' , 8);
        set(gca,  'Xgrid' , 'on' , 'Ygrid' , 'on' , 'Xcolor' , [0 0 0.37] , 'Ycolor' , [0 0 0.37] , 'Box' , 'off' , ...
                   'FontWeight' , 'Bold' , 'FontSize' , 7 , 'FontAngle' , 'Oblique' , 'XMinorTick' , 'on' , 'YMinorTick' , 'on'); clc;

        figure(9);
        set(figure(9), 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [  hpan-757   vpan-616    hpan-664  vpan-464  ] , ...
                'name' , '7. Sieve Analysis - % Percent of Sand'  );
        if Muestrapor(1,1) == 0 && Muestrapor(1,2) ~= 0 && Muestrapor(1,3) ~= 0
                    hotcake(1,1) =Muestrapor(1,2);
                    hotcake(1,2) =Muestrapor(1,3);
                    explode = [ 0 1 ];
                    piecg = pie(hotcake, explode);
                    hdlleg = legend('Sand', 'Fine', -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
                    set(piecg(1,4), 'FontUnits', 'pixels', 'FontSize' , 10);
         elseif Muestrapor(1,1) ~= 0 && Muestrapor(1,2) ~= 0 && Muestrapor(1,3) == 0
                    hotcake(1,1) =Muestrapor(1,1);
                    hotcake(1,2) =Muestrapor(1,2);
                    explode = [ 0 1 ];
                    piecg = pie(hotcake, explode);
                    hdlleg = legend('Gravel' , 'Sand',  -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
                    set(piecg(1,4), 'FontUnits', 'pixels', 'FontSize' , 10);
        elseif Muestrapor(1,1) ~= 0 && Muestrapor(1,2) ~= 0 && Muestrapor(1,3) ~= 0
                    explode = [ 0 1 0 ];
                    piecg = pie(Muestrapor, explode);
                    hdlleg = legend('Gravel' , 'Sand', 'Fine', -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
                    set(piecg(1,4), 'FontUnits', 'pixels', 'FontSize' , 10);
                    set(piecg(1,6), 'FontUnits', 'pixels', 'FontSize' , 10);
        elseif Muestrapor(1,1) ~= 0 && Muestrapor(1,2) == 0 && Muestrapor(1,3) ~= 0
                    hotcake(1,1) = Muestrapor(1,1);
                    hotcake(1,2) = Muestrapor(1,3);
                    explode = [ 0 1 ];
                    piecg = pie(hotcake, explode);
                    hdlleg = legend('Gravel' , 'Fine',  -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
                    set(piecg(1,4), 'FontUnits', 'pixels', 'FontSize' , 10);
        elseif Muestrapor(1,1) ~= 0 && Muestrapor(1,2) == 0 && Muestrapor(1,3) == 0
                    hotcake = Muestrapor(1,1) ;
                    piecg = pie(hotcake);
                    hdlleg = legend('Gravel' , -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
        elseif Muestrapor(1,1) == 0 && Muestrapor(1,2) ~= 0  && Muestrapor(1,3) == 0
                    hotcake = Muestrapor(1,2) ;
                    piecg = pie(hotcake);
                    hdlleg = legend('Sand' , -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
         elseif Muestrapor(1,1) == 0 && Muestrapor(1,2) == 0 && Muestrapor(1,3) ~= 0
                    hotcake = Muestrapor(1,3) ;
                    piecg = pie(hotcake);
                    hdlleg = legend('Fine' , -1);
                    set(piecg(1,2), 'FontUnits', 'pixels', 'FontSize' , 10);
        end
        set(hdlleg, 'Color' , [ 0.93 0.91 0.85 ], 'Box', 'on'); clc;

        result = figure(11);
        set(result, 'Units' , 'Pixels' , 'NumberTitle' , 'Off' , 'Resize' , 'on' , 'Color' , [0.93 0.91 0.85] , ...
                'Position' , [   10+hpan/3+hpan-686  35   hpan-690 (vpan/3)*2-50] , ...
                'name' , '8. Sieve Analysis - Results' );
         P4 = 12;
        uicontrol(result, 'Style' , 'Text' , 'String', 'The main results of sieve analysis were:' , 'Unit' , 'Pixels' , 'Position' , [0 418 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left' );
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_5 (mm) = ' , resulchar.d_5  ] , 'Unit' , 'Pixels' , 'Position' , [0 405 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_10 (mm) = ' , resulchar.d_10  ] , 'Unit' , 'Pixels' , 'Position' , [0 392 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_16 (mm) = ' , resulchar.d_16  ] , 'Unit' , 'Pixels' , 'Position' , [0 379 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_30 (mm) = ' , resulchar.d_30  ] , 'Unit' , 'Pixels' , 'Position' , [0 366 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_50 (mm) = ' , resulchar.d_50  ] , 'Unit' , 'Pixels' , 'Position' , [0 353 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_60 (mm) = ' , resulchar.d_60  ] , 'Unit' , 'Pixels' , 'Position' , [0 340 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_75 (mm) = ' , resulchar.d_75  ] , 'Unit' , 'Pixels' , 'Position' , [0 327 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_84 (mm) = ' , resulchar.d_84  ] , 'Unit' , 'Pixels' , 'Position' , [0 314 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85], 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['D_95 (mm) = ' , resulchar.d_95  ] , 'Unit' , 'Pixels' , 'Position' , [0 301 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Mean Grain Size  (mm) = ' , resulchar.MeanGS  ] , 'Unit' , 'Pixels' , 'Position' , [0 288 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Standard Deviation ' , resulchar.SDe  ] , 'Unit' , 'Pixels' , 'Position' , [0 273 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Skewness = ' , resulchar.S_ke  ] , 'Unit' , 'Pixels' , 'Position' , [0 260 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Kurtosis = ' , resulchar.K_ur  ] , 'Unit' , 'Pixels' , 'Position' , [0 247 200 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Coeficient of uniformity (Cc) = ' , resulchar.CC  ] , 'Unit' , 'Pixels' , 'Position' , [0 236 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Coeficient of curvature (Cu) = ' , resulchar.Cu  ] , 'Unit' , 'Pixels' , 'Position' , [0 223 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Standard Deviation Classification = ' , resulchar.SD  ] , 'Unit' , 'Pixels' , 'Position' , [0 210 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Skewness Classification = ' , resulchar.Sk  ] , 'Unit' , 'Pixels' , 'Position' , [0 197 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Kurtosis Classification = ' , resulchar.kurt  ] , 'Unit' , 'Pixels' , 'Position' , [0 184 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['ASTM Classification = ' , resulchar.ASTM  ] , 'Unit' , 'Pixels' , 'Position' , [0 169 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['% of Gravel = ' , resulchar.Grava  ] , 'Unit' , 'Pixels' , 'Position' , [0 158 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['% of Sand = ' , resulchar.arena  ] , 'Unit' , 'Pixels' , 'Position' , [0 145 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['% of Fine = ' , resulchar.fino  ] , 'Unit' , 'Pixels' , 'Position' , [0 132 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['Wentworth classification = ' , resulchar.Wentworth  ] , 'Unit' , 'Pixels' , 'Position' , [0 117 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');
        uicontrol(result, 'Style' , 'Text' , 'String', ['USCS classification = ' , resulchar.SUCS  ] , 'Unit' , 'Pixels' , 'Position' , [0 103 300 P4] , ...
                                   'BackGroundColor' , [0.93 0.91 0.85] , 'FontWeight' , 'Bold', 'ForegroundColor' , [0.5 0 0 ] , 'HorizontalAlignment', 'left');


text:
function [s1 s2 distances] = findTwoNearest(Input,nodes) %#codegen

NumOfNodes = size(nodes,2);

distances = zeros(1,NumOfNodes);
for i=1:NumOfNodes
    distances(i) = norm(Input - nodes(:,i));
end

[sdistances indices] = sort(distances);

s1 = indices(1);
s2 = indices(2);

% MEX Code Generation:

% mexcfg = coder.config('mex');
% mexcfg.DynamicMemoryAllocation = 'AllVariableSizeArrays';
% codegen -config mexcfg findTwoNearest -args {coder.typeof(In(:,n),[Inf 1]),coder.typeof(nodes,[Inf Inf])}

author: Sable

repo: mcbench-benchmarks


text:
function y = negbinztr_lpdf(x,l,r)
%NEGBINZTR_LPDF Zero trunc. negative binomial log probability density function
%
%  Description
%    Y = NEGBINZTR_LPDF(X,L,R) returns the log of zero truncated
%    Negative binomial probability density function with location
%    parameter L and dispersion parameter R (0<R<infty).
%
%    Negative binomial has different parameterizations and we use the form
%      p(x|l,r) = (r/(r+l))^r * gamma(r+y)/(gamma(r)*gamma(y+1))
%                             * (l/(r+l))^y
%    which approaches Poisson distribution when R goes to infinity.
%    Zero truncated Negative Binomial has p(x==0|l,r)=0.
%
%    The size of Y is the common size of X, L and R. A scalar input
%    functions as a constant matrix of the same size as the other
%    input.
%
%     Note that the density function is zero unless X is an
%     integer.
%
%  See also
%    NEGBINZTR_PDF, NEGBIN_LPDF, NEGBIN_PDF
%
% Copyright (c) 2010 Jarno Vanhatalo
% Copyright (c) 2010-2011  Aki Vehtari

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

if isscalar(r) && isscalar(l)
  if l<0
    y = repmat(NaN,size(x));
  else
    y = repmat(-Inf,size(x));
  end
  if l==0
    y(x==1) = 0;
  elseif l>0
    k = (x >= 1 & x == round(x));
    if (any(k))
      lp0=r.*(log(r) - log(r+l));
      y(k)=lp0 + gammaln(r+x(k)) - gammaln(r) - gammaln(x(k)+1) + x(k).*(log(l) - log(r+l)) -log(1-exp(lp0));
    end
  end
else
  y = repmat(-Inf,size(x));
  y(l < 0) = NaN;
  y(x==1 & l==0) = 0;

  k = (x >= 1 & x == round(x) & l > 0);
  if (any(k))
    lp0=r(k).*(log(r(k)) - log(r(k)+l(k)));
    y(k)=lp0 + gammaln(r(k)+x(k)) - gammaln(r(k)) - gammaln(x(k)+1) + x(k).*(log(l(k)) - log(r(k)+l(k))) -log(1-exp(lp0));
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
classdef filterer < dml.method
% FILTERER filtering approach to feature selection.
%
%   DESCRIPTION
%   The filterer produces an ordering of the features based on a univariate measure
%   and uses that ordering to greedily determine the 'optimal' feature set
%   or to select the best m features.
%
%   If filter is a string, e.g. 'ttest' then it calls rankfeatures in the bioinformatics
%   toolbox. If it is a function it will just evaluate that function on
%   each feature and the output.
%
%   EXAMPLE
%   dml.filterer('validator',dml.crossvalidator('mva',dml.naive,'stat','accuracy'))
%
%   DEVELOPER
%   Marcel van Gerven (m.vangerven@donders.ru.nl)

  properties

    maxfeatures = Inf; % maximum number of used features

    filter = @(x,y)(abs(corr(x,y))); % the used filter

    value % the value  of the filter for each feature

    order % the ordering of the features according to value

    criterion % the value of the validator for feature 1:N

    validator % the used crossvalidator

    subset % subset of selected features

  end

  methods

    function obj = filterer(varargin)

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

    end

    function obj = train(obj,X,Y)

      % multiple datasets
      if iscell(X) || iscell(Y)
        obj = ft_mv_ndata('mvmethod',obj);
        obj = obj.train(X,Y);
        return;
      end

      maxf = min(obj.maxfeatures,size(X,2));

      if obj.verbose
        if isempty(obj.validator)
          fprintf('selecting %d features based on %s filter\n',maxf,obj.filter);
        else
          fprintf('selecting optimal features based on %s filter\n',obj.filter);
        end
      end

      if isempty(obj.validator)
      % just take the best N features based on a univariate criterion

        % compute function values
        nf = size(X,2);
        obj.value = zeros(1,nf);
        if ischar(obj.filter)

          if obj.verbose
            fprintf('computing %s filter using rankfeatures\n');
          end
          Z = rankfeatures(X',Y','Criterion',obj.filter);
          for j=1:nf
            obj.value(j) = Z(j);
          end

        else

          for j=1:nf

            if obj.verbose
              fprintf('computing %s filter for feature %d of %d\n',obj.filter,j,nf);
            end

            obj.value(j) = obj.filter(X(:,j),Y);
          end

        end

        % get ordering of the features
        [tmp,obj.order] = sort(obj.value,'descend');

        % select the best n features
        obj.subset = obj.order(1:maxf);

      else
        % currently we use a simple approach. A better approach would be
        % the following:
        %
        % use procedure to determine the best feature subset.
        % first we use the procedure to compute the univariate filter
        % values per fold. Then we determine the average best number of features
        % over folds using the crossvalidation procedure. Then we recompute the
        % univariate filter values for all data and return the optimal
        % subset

        % compute function values
        nf = size(X,2);
        obj.value = zeros(1,nf);
        if ischar(obj.filter)

          if obj.verbose
            fprintf('computing %s filter using rankfeatures\n');
          end
          Z = rankfeatures(X',Y','Criterion',obj.filter);
          for j=1:nf
            obj.value(j) = Z(j);
          end

        else
          for j=1:nf

            if obj.verbose
              fprintf('computing %s filter for feature %d of %d\n',obj.filter,j,nf);
            end

            obj.value(j) = obj.filter(X(:,j),Y);

          end

        end

        % get ordering of the features
        [tmp,obj.order] = sort(obj.value,'descend');

        obj.subset = [];
        metric = -Inf;

        obj.criterion = zeros(1,maxf);
        for f=1:maxf

          if obj.verbose
            fprintf('evaluating %d out of %d features; ',f,maxf);
          end

          cv = obj.validator.train(X(:,obj.order(1:f)),Y);

          m = cv.statistic;

          if obj.verbose, fprintf('criterion: %.2g\n',m); end

          if m > metric
            obj.subset = obj.order(1:f);
            metric = m;
          end

          obj.criterion(f) = m;

        end
      end

    end

    function Y = test(obj,X)

      Y = X(:,obj.subset);

    end

    function m = model(obj)
      % return logical array with ones indicating selected features
      % m.selected logical array with ones indicating selected features
      % m.importance importance according to the used filter

      indim = obj.indims;
      if length(indim) == 1, indim = [1 indim]; end

      m1 = zeros(indim);
      m1(obj.subset) = 1;
      m.selected = m1;

      m1 = zeros(indim);
      m1(:) = obj.value;
      m.importance = m1;

    end

  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function roidb_new = do_proposal_test(conf, model_stage, imdb, roidb)
    aboxes                      = proposal_test(conf, imdb, ...
                                        'net_def_file',     model_stage.test_net_def_file, ...
                                        'net_file',         model_stage.output_model_file, ...
                                        'cache_name',       model_stage.cache_name);

    aboxes                      = boxes_filter(aboxes, model_stage.nms.per_nms_topN, model_stage.nms.nms_overlap_thres, model_stage.nms.after_nms_topN, conf.use_gpu);

    roidb_regions               = make_roidb_regions(aboxes, imdb.image_ids);

    roidb_new                   = roidb_from_proposal(imdb, roidb, roidb_regions, ...
                                        'keep_raw_proposal', false);
end

function aboxes = boxes_filter(aboxes, per_nms_topN, nms_overlap_thres, after_nms_topN, use_gpu)
    % to speed up nms
    if per_nms_topN > 0
        aboxes = cellfun(@(x) x(1:min(length(x), per_nms_topN), :), aboxes, 'UniformOutput', false);
    end
    % do nms
    if nms_overlap_thres > 0 && nms_overlap_thres < 1
        if use_gpu
            for i = 1:length(aboxes)
                aboxes{i} = aboxes{i}(nms(aboxes{i}, nms_overlap_thres, use_gpu), :);
            end
        else
            parfor i = 1:length(aboxes)
                aboxes{i} = aboxes{i}(nms(aboxes{i}, nms_overlap_thres), :);
            end
        end
    end
    aver_boxes_num = mean(cellfun(@(x) size(x, 1), aboxes, 'UniformOutput', true));
    fprintf('aver_boxes_num = %d, select top %d\n', round(aver_boxes_num), after_nms_topN);
    if after_nms_topN > 0
        aboxes = cellfun(@(x) x(1:min(length(x), after_nms_topN), :), aboxes, 'UniformOutput', false);
    end
end

function regions = make_roidb_regions(aboxes, images)
    regions.boxes = aboxes;
    regions.images = images;
end

author: ShaoqingRen


text:
function value = index10 ( i_min, i, i_max, j_min, j, j_max )

%*****************************************************************************80
%
%% INDEX10 indexes a 2D array by rows, with a zero base.
%
%  Discussion:
%
%    Entries of the array are indexed starting at entry (I_MIN,J_MIN),
%    and increasing the column index first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I_MIN, I, I_MAX, for row indices,
%    the minimum, the index, and the maximum.
%
%    Input, integer J_MIN, J, J_MAX, for column indices,
%    the minimum, the index, and the maximum.
%
%    Output, integer VALUE, the index of element (I,J).
%
  index_min = 0;

  value = index_min ...
             +                         ( j - j_min ) ...
             + ( i - i_min ) * ( j_max + 1 - j_min );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function k=kurtosis(x, flag, dim)
% Estimate the 4th standardized moment from a vector or matix of data.
%
%   K=kurtosis(X,FLAG,DIM)
%
% INPUTS
%   X     - Data to be used in skewness calculation
%   FLAG  - Not implemented
%   DIM   - Dimension to compute the kurtosis on if x is not a vector

% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 2/25/2006

if nargin==1
    dim=1;
end

ncm1 = mean(x,dim);
ncm2 = mean(x.^2,dim);
ncm3 = mean(x.^3,dim);
ncm4 = mean(x.^4,dim);

cm2 = ncm2-ncm1.^2;
cm4 = ncm4 - 4*ncm3.*ncm1 + 6*ncm1.^2.*ncm2 - 4.*ncm1.^4 + ncm1.^4;

k = cm4./cm2.^2;

author: bashtage

repo: mfe-toolbox


text:
function b = isany(list, func)
% BOOL = ISANY(ARRAY)
% BOOL = ISANY(ARRAY, FUNC)
%   Without FUNC, ISANY returns true if any element of ARRAY is logically
%   true.  With FUNC, every element first has FUNC applied to it, and this
%   result is used instead - thus if any application of FUNC across ARRAY
%   is true, then ISANY(ARRAY, FUNC) is true.
%

  if nargin < 2, func = @(x) x; end
  % We don't use foldl so as to incur better performance via short circuit.
  b = 0;
  if iscell(list)
    for i = 1:numel(list)
      if func(list{i})
        b = 1;
        return;
      end
    end
  else
    for i = 1:numel(list)
      if func(list(i))
        b = 1;
        return;
      end
    end
  end
end


text:
function [imageCha, ssim_map, metrics] = algo_FCSA_SLEP_2D_real( obj,input )
% 2D real-valued variant of FCSA
% based on Huang et al. paper on FCSA
% with Liu and Ye's idea to use the Moreau-Yosida regularization to solve the group lasso problem, using the SLEP package.
%
% input:
% obj           CS reconstruction object (holding all parameters)
% input         struct containing recon parameters and image
%
% output:
% imageCha      reconstructed channel individual image
% ssim_map      structural similarity map
% metrics       evaluation metrics
%
% (c) Marc Fischer, Thomas Kuestner, May 2015
% -------------------------------------------------------------------------

timer_SLEP = tic;

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

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

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

% from input:
b=input.b;
mask = input.mask;
G_slep = input.G_slep;
Gt_slep = input.Gt_slep;
ind = input.ind;
nodesize = input.nodesize;
waveS_l1 = input.waveS_l1;
waveS_l12 = input.waveS_l12;
waveS_SLEP = input.waveS_SLEP; % also change name in CS_2D
SLEP_extend_y = waveS_SLEP(waveletStages+2,1) - waveS_l12(waveletStages+2,1);
SLEP_extend_x = waveS_SLEP(waveletStages+2,2) - waveS_l12(waveletStages+2,2);
z_proxA = cell(1,nCha);

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

% initialize cells/vectors
for j=1:nCha
    FTb{1,j} = real(iFFT2D(b{1,j}));
    % y{1,j} = real(FTb{1,j});
    % y{1,j+nCha} = imag(FTb{1,j});
end;
z = FTb; % starting point

x_wave = cell(1,nCha);
x_helper = x_wave;
x_wave_helper = x_wave;
x_tv = x_wave;
x_nltv = x_wave;
for j=1:nCha
    x_nltv{1,j} = 0;
end;
x_g = x_wave;
x_g_helper = x_wave;
x_g_proxA = x_wave;
y = x_wave;
slep_thresh = x_wave;

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

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

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

%% recon
dispProgress('Proximal Average', 0, maxitr);
for itr = 1:maxitr  % total iter counter

    t_new = (1+sqrt(1+4*t_old^2))/2;
    t_old = t_new;

    y_old = z;   % y_old = y for complex case

%% landweber step
    for j = 1:nCha
        x_helper{1,j} = real(iFFT2D(FFT2D_mask(z{1,j},mask))) -FTb{1,j} ;
        z{1,j} = z{1,j} - real(x_helper{1,j})/L;
    end;

%% l1-Wavelet
    if flagWave
            for j = 1:nCha % 2*nCha
                    x_wave_helper{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l1);
                    x_wave_helper{1,j} = softthresh_real(x_wave_helper{1,j},threshold_wave(j));
                    x_wave{1,j} = waverec2(x_wave_helper{1,j},waveS_l1,waveletFilterName_l1);
            end;
    end;

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

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

%% l12-Wavelet - via SLEP
    if flagGroup
        for j = 1:nCha
            if flag_extendImage
                z_proxA{1,j} = extend_image(z{1,j}, waveS_SLEP, waveletStages, SLEP_extend_y, SLEP_extend_x);
                x_g_helper{1,j} = wavedec2(z_proxA{1,j},waveletStages,waveletFilterName_l12);
            else
                x_g_helper{1,j} = wavedec2(z{1,j},waveletStages,waveletFilterName_l12);
            end;

            slep_thresh{1,j} = x_g_helper{1,j}(G_slep);
            slep_thresh{1,j} = altra( slep_thresh{1,j}, size(slep_thresh{1,j},2), ind, nodesize);
            slep_thresh{1,j} = slep_thresh{1,j}';
            slep_thresh{1,j} = slep_thresh{1,j}(Gt_slep);

            x_g_proxA{1,j} = waverec2(slep_thresh{1,j},waveS_SLEP,waveletFilterName_l12);
            x_g{1,j} =  x_g_proxA{1,j}(1:end-SLEP_extend_y,1:end-SLEP_extend_x);
        end;
    end;

%% add prox(.)
    for j = 1:nCha
        y{1,j} = zeros(n1,n2);
        if flagWave y{1,j} = y{1,j} + x_wave{1,j}.*regularizerWeights(1); end;
        if flagTV y{1,j} = y{1,j} + x_tv{1,j}.*regularizerWeights(2); end;
        if flagGroup y{1,j} = y{1,j} + x_g{1,j}.*regularizerWeights(3); end;
        if flagNLTV y{1,j} = y{1,j} + x_nltv{1,j}.*regularizerWeights(4); end;

        if itr < itrNLTV
            y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagWave.*regularizerWeights(1) + flagGroup.*regularizerWeights(3));
        else
            y{1,j} = y{1,j}/(flagTV.*regularizerWeights(2) + flagNLTV.*regularizerWeights(4) + flagWave.*regularizerWeights(1) + flagGroup.*regularizerWeights(3));
        end;

        if flag_fast
        y{1,j}=y{1,j}+((t_old-1)/t_new).*(y{1,j}-y_old{1,j});
        end;
    end;

%% metrics of current itr:
%     disp(itr);
    dispProgress('Proximal Average', itr/maxitr);

   	metrics.xtime(itr+1)= toc(timer_SLEP);
    for j = 1:nCha
        z{1,j} = y{1,j};
    end;
%     [metrics, ssim_map{1,2}] = get_metrics_itr( im_ref, im_ref_full, z, itr, maxitr, nCha, n1, n2, metrics, obj.K_1, obj.K_2, obj.W_size, obj.W_sigma );

end;
dispProgress('Proximal Average', 'Close');

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

end

author: thomaskuestner


text:
function tests = test_wl_transform
  tests = functiontests(localfunctions);
end

function test_fwd_po_wavelet_transform(testCase)
    x = [0 5 -3 3 -4 -6 -1 6 1 -7 -2 -7 1 -8 1 -3];
    h = spx.wavelet.daubechies.quad_mirror_filter(20);
    L = 1;
    w = spx.wavelet.transform.forward_periodized_orthogonal(h, x, L);
    w_expected = [-0.886174910726689  -7.599106463564170  -4.126070413696197  -3.831372392026132...
    -1.964972811707912  -6.298466192390180  -5.868195699184612   3.365505053418349...
    -3.291527529820000   3.514441633073000   5.705304137811999  -1.015576057620000...
    1.086658992536000  -1.158475051501000  -6.297575180086999  -5.614318756272999];
    verifyEqual(testCase, w, w_expected, 'AbsTol', 1e-6);
    verifyEqual(testCase, norm(x), norm(w), 'AbsTol', 1e-6);
    x = x';
    w = spx.wavelet.transform.forward_periodized_orthogonal(h, x, L);
    verifyEqual(testCase, w, w_expected', 'AbsTol', 1e-6);
    verifyEqual(testCase, norm(x), norm(w), 'AbsTol', 1e-6);

    % We now restrict L to 2
    L = 2;
    x = x';
    w = spx.wavelet.transform.forward_periodized_orthogonal(h, x, L);
    w_expected = [0.004217817988112  -9.167655098625628  -0.377456049193954  -2.459106670217828...
    -1.964972811707912  -6.298466192390180  -5.868195699184612   3.365505053418349...
    -3.291527529820000   3.514441633073000   5.705304137811999  -1.015576057620000...
    1.086658992536000  -1.158475051501000  -6.297575180086999  -5.614318756272999];
    verifyEqual(testCase, w, w_expected, 'AbsTol', 1e-6);
    verifyEqual(testCase, norm(x), norm(w), 'AbsTol', 1e-6);
    x = x';
    w = spx.wavelet.transform.forward_periodized_orthogonal(h, x, L);
    verifyEqual(testCase, w, w_expected', 'AbsTol', 1e-6);
    verifyEqual(testCase, norm(x), norm(w), 'AbsTol', 1e-6);
end

function test_inv_po_wavelet_transform(testCase)
    x = [0 5 -3 3 -4 -6 -1 6 1 -7 -2 -7 1 -8 1 -3];
    h = spx.wavelet.daubechies.quad_mirror_filter(20);
    L = 1;
    w = spx.wavelet.transform.forward_periodized_orthogonal(h, x, L);
    y = spx.wavelet.transform.inverse_periodized_orthogonal(h, w, L);
    verifyEqual(testCase, x, y, 'AbsTol', 1e-6);
end

author: indigits

repo: sparse-plex


text:
function [h,h1,h2] = maxflatI(K,M)
% [h,h1,h2] = maxflatI(K,M)
% Maxflat Type-I FIR filter
%   2K zeros at z=-1
%   2M zeros away from z=-1
%   h = conv(h1,h2);
%   h1 : all zeros at z=-1
%   h2 : all other zeros
%
% Note: if K = M+1, then h is halfband.
%
% Reference:
% O. Herrmann, "On the approximation problem in Nonrecursive
% Digital Filter Design", IEEE Trans. on Circuit Theory,
% Vol. 18, No. 3, May 1971, pp. 411-413
%
% % Example
% [h,h1,h2] = maxflatI(4,6);

% Ivan Selesnick
% selesi@nyu.edu
% NYU - School of Engineering

h2 = 1;
hi = 1;
c  = 1;
for k = 1:M
   hi = conv(hi,[-1 2 -1]/4);
   c  = c*(K-1+k)/k;
   h2 = [0 h2 0] + c*hi;
end

h1 = 1;
for k = 1:2*K
   h1 = conv(h1,[1 1]/2);
end

h = conv(h1,h2);


text:
function ellipse_t = fit_ellipse( x,y,axis_handle )
%
% fit_ellipse - finds the best fit to an ellipse for the given set of points.
%
% Format:   ellipse_t = fit_ellipse( x,y,axis_handle )
%
% Input:    x,y         - a set of points in 2 column vectors. AT LEAST 5 points are needed !
%           axis_handle - optional. a handle to an axis, at which the estimated ellipse
%                         will be drawn along with it's axes
%
% Output:   ellipse_t - structure that defines the best fit to an ellipse
%                       a           - sub axis (radius) of the X axis of the non-tilt ellipse
%                       b           - sub axis (radius) of the Y axis of the non-tilt ellipse
%                       phi         - orientation in radians of the ellipse (tilt)
%                       X0          - center at the X axis of the non-tilt ellipse
%                       Y0          - center at the Y axis of the non-tilt ellipse
%                       X0_in       - center at the X axis of the tilted ellipse
%                       Y0_in       - center at the Y axis of the tilted ellipse
%                       long_axis   - size of the long axis of the ellipse
%                       short_axis  - size of the short axis of the ellipse
%                       status      - status of detection of an ellipse
%
% Note:     if an ellipse was not detected (but a parabola or hyperbola), then
%           an empty structure is returned

% =====================================================================================
%                  Ellipse Fit using Least Squares criterion
% =====================================================================================
% We will try to fit the best ellipse to the given measurements. the mathematical
% representation of use will be the CONIC Equation of the Ellipse which is:
%
%    Ellipse = a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0
%
% The fit-estimation method of use is the Least Squares method (without any weights)
% The estimator is extracted from the following equations:
%
%    g(x,y;A) := a*x^2 + b*x*y + c*y^2 + d*x + e*y = f
%
%    where:
%       A   - is the vector of parameters to be estimated (a,b,c,d,e)
%       x,y - is a single measurement
%
% We will define the cost function to be:
%
%   Cost(A) := (g_c(x_c,y_c;A)-f_c)'*(g_c(x_c,y_c;A)-f_c)
%            = (X*A+f_c)'*(X*A+f_c)
%            = A'*X'*X*A + 2*f_c'*X*A + N*f^2
%
%   where:
%       g_c(x_c,y_c;A) - vector function of ALL the measurements
%                        each element of g_c() is g(x,y;A)
%       X              - a matrix of the form: [x_c.^2, x_c.*y_c, y_c.^2, x_c, y_c ]
%       f_c            - is actually defined as ones(length(f),1)*f
%
% Derivation of the Cost function with respect to the vector of parameters "A" yields:
%
%   A'*X'*X = -f_c'*X = -f*ones(1,length(f_c))*X = -f*sum(X)
%
% Which yields the estimator:
%
%       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%       |  A_least_squares = -f*sum(X)/(X'*X) ->(normalize by -f) = sum(X)/(X'*X)  |
%       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% (We will normalize the variables by (-f) since "f" is unknown and can be accounted for later on)
%
% NOW, all that is left to do is to extract the parameters from the Conic Equation.
% We will deal the vector A into the variables: (A,B,C,D,E) and assume F = -1;
%
%    Recall the conic representation of an ellipse:
%
%       A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
%
% We will check if the ellipse has a tilt (=orientation). The orientation is present
% if the coefficient of the term "x*y" is not zero. If so, we first need to remove the
% tilt of the ellipse.
%
% If the parameter "B" is not equal to zero, then we have an orientation (tilt) to the ellipse.
% we will remove the tilt of the ellipse so as to remain with a conic representation of an
% ellipse without a tilt, for which the math is more simple:
%
% Non tilt conic rep.:  A`*x^2 + C`*y^2 + D`*x + E`*y + F` = 0
%
% We will remove the orientation using the following substitution:
%
%   Replace x with cx+sy and y with -sx+cy such that the conic representation is:
%
%   A(cx+sy)^2 + B(cx+sy)(-sx+cy) + C(-sx+cy)^2 + D(cx+sy) + E(-sx+cy) + F = 0
%
%   where:      c = cos(phi)    ,   s = sin(phi)
%
%   and simplify...
%
%       x^2(A*c^2 - Bcs + Cs^2) + xy(2A*cs +(c^2-s^2)B -2Ccs) + ...
%           y^2(As^2 + Bcs + Cc^2) + x(Dc-Es) + y(Ds+Ec) + F = 0
%
%   The orientation is easily found by the condition of (B_new=0) which results in:
%
%   2A*cs +(c^2-s^2)B -2Ccs = 0  ==> phi = 1/2 * atan( b/(c-a) )
%
%   Now the constants   c=cos(phi)  and  s=sin(phi)  can be found, and from them
%   all the other constants A`,C`,D`,E` can be found.
%
%   A` = A*c^2 - B*c*s + C*s^2                  D` = D*c-E*s
%   B` = 2*A*c*s +(c^2-s^2)*B -2*C*c*s = 0      E` = D*s+E*c
%   C` = A*s^2 + B*c*s + C*c^2
%
% Next, we want the representation of the non-tilted ellipse to be as:
%
%       Ellipse = ( (X-X0)/a )^2 + ( (Y-Y0)/b )^2 = 1
%
%       where:  (X0,Y0) is the center of the ellipse
%               a,b     are the ellipse "radiuses" (or sub-axis)
%
% Using a square completion method we will define:
%
%       F`` = -F` + (D`^2)/(4*A`) + (E`^2)/(4*C`)
%
%       Such that:    a`*(X-X0)^2 = A`(X^2 + X*D`/A` + (D`/(2*A`))^2 )
%                     c`*(Y-Y0)^2 = C`(Y^2 + Y*E`/C` + (E`/(2*C`))^2 )
%
%       which yields the transformations:
%
%           X0  =   -D`/(2*A`)
%           Y0  =   -E`/(2*C`)
%           a   =   sqrt( abs( F``/A` ) )
%           b   =   sqrt( abs( F``/C` ) )
%
% And finally we can define the remaining parameters:
%
%   long_axis   = 2 * max( a,b )
%   short_axis  = 2 * min( a,b )
%   Orientation = phi
%
%

% initialize
orientation_tolerance = 1e-3;

% empty warning stack
warning( '' );

% prepare vectors, must be column vectors
x = x(:);
y = y(:);

% remove bias of the ellipse - to make matrix inversion more accurate. (will be added later on).
mean_x = mean(x);
mean_y = mean(y);
x = x-mean_x;
y = y-mean_y;

% the estimation for the conic equation of the ellipse
X = [x.^2, x.*y, y.^2, x, y ];
a = sum(X)/(X'*X);

% check for warnings
if ~isempty( lastwarn )
    disp( 'stopped because of a warning regarding matrix inversion' );
    ellipse_t = [];
    return
end

% extract parameters from the conic equation
[a,b,c,d,e] = deal( a(1),a(2),a(3),a(4),a(5) );

% remove the orientation from the ellipse
if ( min(abs(b/a),abs(b/c)) > orientation_tolerance )

    orientation_rad = 1/2 * atan( b/(c-a) );
    cos_phi = cos( orientation_rad );
    sin_phi = sin( orientation_rad );
    [a,b,c,d,e] = deal(...
        a*cos_phi^2 - b*cos_phi*sin_phi + c*sin_phi^2,...
        0,...
        a*sin_phi^2 + b*cos_phi*sin_phi + c*cos_phi^2,...
        d*cos_phi - e*sin_phi,...
        d*sin_phi + e*cos_phi );
    [mean_x,mean_y] = deal( ...
        cos_phi*mean_x - sin_phi*mean_y,...
        sin_phi*mean_x + cos_phi*mean_y );
else
    orientation_rad = 0;
    cos_phi = cos( orientation_rad );
    sin_phi = sin( orientation_rad );
end

% check if conic equation represents an ellipse
test = a*c;
switch (1)
case (test>0),  status = '';
case (test==0), status = 'Parabola found';  warning( 'fit_ellipse: Did not locate an ellipse' );
case (test<0),  status = 'Hyperbola found'; warning( 'fit_ellipse: Did not locate an ellipse' );
end

% if we found an ellipse return it's data
if (test>0)

    % make sure coefficients are positive as required
    if (a<0), [a,c,d,e] = deal( -a,-c,-d,-e ); end

    % final ellipse parameters
    X0          = mean_x - d/2/a;
    Y0          = mean_y - e/2/c;
    F           = 1 + (d^2)/(4*a) + (e^2)/(4*c);
    [a,b]       = deal( sqrt( F/a ),sqrt( F/c ) );
    long_axis   = 2*max(a,b);
    short_axis  = 2*min(a,b);

    % rotate the axes backwards to find the center point of the original TILTED ellipse
    R           = [ cos_phi sin_phi; -sin_phi cos_phi ];
    P_in        = R * [X0;Y0];
    X0_in       = P_in(1);
    Y0_in       = P_in(2);

    % pack ellipse into a structure
    ellipse_t = struct( ...
        'a',a,...
        'b',b,...
        'phi',orientation_rad,...
        'X0',X0,...
        'Y0',Y0,...
        'X0_in',X0_in,...
        'Y0_in',Y0_in,...
        'long_axis',long_axis,...
        'short_axis',short_axis,...
        'status','' );
else
    % report an empty structure
    ellipse_t = struct( ...
        'a',[],...
        'b',[],...
        'phi',[],...
        'X0',[],...
        'Y0',[],...
        'X0_in',[],...
        'Y0_in',[],...
        'long_axis',[],...
        'short_axis',[],...
        'status',status );
end

% check if we need to plot an ellipse with it's axes.
if (nargin>2) & ~isempty( axis_handle ) & (test>0)

    % rotation matrix to rotate the axes with respect to an angle phi
    R = [ cos_phi sin_phi; -sin_phi cos_phi ];

    % the axes
    ver_line        = [ [X0 X0]; Y0+b*[-1 1] ];
    horz_line       = [ X0+a*[-1 1]; [Y0 Y0] ];
    new_ver_line    = R*ver_line;
    new_horz_line   = R*horz_line;

    % the ellipse
    theta_r         = linspace(0,2*pi);
    ellipse_x_r     = X0 + a*cos( theta_r );
    ellipse_y_r     = Y0 + b*sin( theta_r );
    rotated_ellipse = R * [ellipse_x_r;ellipse_y_r];

    % draw
    hold_state = get( axis_handle,'NextPlot' );
    set( axis_handle,'NextPlot','add' );
    plot( new_ver_line(1,:),new_ver_line(2,:),'r' );
    plot( new_horz_line(1,:),new_horz_line(2,:),'r' );
    plot( rotated_ellipse(1,:),rotated_ellipse(2,:),'r' );
    set( axis_handle,'NextPlot',hold_state );
end


text:
function [param, names] = gpnddisimExtractParam(model)

% NDDISIMEXTRACTPARAM Extract the parameters of a NDDISIM model.
% FORMAT
% DESC extracts the model parameters from a structure containing
% the information about a Gaussian process for single input motif
% modelling.
% ARG model : the model structure containing the information about
% the model.
% RETURN params : a vector of parameters from the model.
%
% SEEALSO : gpdisimCreate, gpdisimExpandParam, modelExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007
%
% COPYRIGHT : Jaakko Peltonen, 2011

% GPSIM

if nargout>1
  [param, names] = kernExtractParam(model.kern);

%  fprintf(1,'gpdisimExtractParam: parameters from kernels\n');
%  param
%  names

  %  model.mu
  %  length(model.mu)
  %for i=1:length(model.mu),
  for i=1:model.numGenes,
    names{end+1}=['Basal transcription ' num2str(i)];
  end
  if model.numGenes>0,
    if (model.use_disimstartmean==1),
      for i=1:model.numGenes,
	names{end+1}=['NDDISIM startmean ' num2str(i)];
      end;
    end;
  end;

  names{end+1}=['NDSIM mean ' num2str(i)];
else
  param = kernExtractParam(model.kern);
end

%length(param)
%model.numGenes
%length(model.B)
if model.numGenes>0,
  for k=1:model.numGenes,
    param = [param doTransform(model.B(k), 'xtoa',model.bTransform(k))];
  end;
end;
%length(param)

if model.numGenes>0,
  if (model.use_disimstartmean==1),
    for i=1:model.numGenes,
      param = [param doTransform(model.disimStartMean(k),'xtoa',model.disimStartMeanTransform(k))];
    end;
  end;
end;

param = [param doTransform(model.simMean,'xtoa',model.simMeanTransform)];

if isfield(model, 'fix')
%model.fix
  for i = 1:length(model.fix)
    param(model.fix(i).index) = model.fix(i).value;
  end
end
param = real(param);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [lo, hi] = afb2D_A(x, af, d)

% 2D Analysis Filter Bank
% (along one dimension only)
%
% [lo, hi] = afb2D_A(x, af, d);
% INPUT:
%    x - NxM matrix, where min(N,M) > 2*length(filter)
%           (N, M are even)
%    af - analysis filter for the columns
%    af(:, 1) - lowpass filter
%    af(:, 2) - highpass filter
%    d - dimension of filtering (d = 1 or 2)
% OUTPUT:
%     lo, hi - lowpass, highpass subbands
%
% % Example
% x = rand(32,64);
% [af, sf] = farras;
% [lo, hi] = afb2D_A(x, af, 1);
% y = sfb2D_A(lo, hi, sf, 1);
% err = x - y;
% max(max(abs(err)))

lpf = af(:, 1);     % lowpass filter
hpf = af(:, 2);     % highpass filter

if d == 2
   x = x';
end

N = size(x,1);
L = size(af,1)/2;
x = cshift2D(x,-L);

lo = upfirdn(x, lpf, 1, 2);
lo(1:L, :) = lo(1:L, :) + lo([1:L]+N/2, :);
lo = lo(1:N/2, :);

hi = upfirdn(x, hpf, 1, 2);
hi(1:L, :) = hi(1:L, :) + hi([1:L]+N/2, :);
hi = hi(1:N/2, :);

if d == 2
   lo = lo';
   hi = hi';
end


text:
function test_ft_inside_headmodel

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_inside_headmodel

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  func = localfunctions;
  for i=1:numel(func)
    fprintf('evaluating %s\n', func2str(func{i}));
    feval(func{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testSphere(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

headmodel = [];
headmodel.r = 120;
headmodel.o = [0 0 40];
headmodel.unit = 'mm';

assert( ft_inside_headmodel([  0   0  40], headmodel));
assert( ft_inside_headmodel([119   0  40], headmodel));
assert(~ft_inside_headmodel([121   0  40], headmodel)); % outside
assert( ft_inside_headmodel([  0   0 -79], headmodel));
assert(~ft_inside_headmodel([  0   0 -81], headmodel)); % outside

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testMesh(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[pos, tri] = mesh_sphere(162);

headmodel = [];
headmodel.bnd.pos = pos*120;
headmodel.bnd.pos(:,3) = headmodel.bnd.pos(:,3) + 40;
headmodel.bnd.tri = tri;
headmodel.unit = 'mm';

assert( ft_inside_headmodel([  0   0  40], headmodel));
assert( ft_inside_headmodel([119   0  40], headmodel));
assert(~ft_inside_headmodel([121   0  40], headmodel)); % outside
assert( ft_inside_headmodel([  0   0 -79], headmodel));
assert(~ft_inside_headmodel([  0   0 -81], headmodel)); % outside


text:
function [xVals,tVals,dxdtVals,exitCode,nextStepSize]=RKAdaptiveOverRange(xStart,tSpan,f,initStepSize,numRefinedSteps,order,solutionChoice,RelTol,AbsTol,maxSteps)
%%RKADAPTIVEOVERRANGE Perform explicit Runge-Kutta integration over a
%                     given range of values using an adaptive step size.
%                     Runge-Kutta methods are derivative-free techniques
%                     for solving ordinary differential equations. That is,
%                     integrating dx/dt=f(x,t) given initial conditions
%                     (xStart,tSpan(1)). More information on available
%                     algorithms for the steps is given in the comments to
%                     the function RungeKStep.
%
%INPUTS: xStart The NX1 state vector at time tSpan(1).
%         tSpan A 2X1 or 1X2 vector where tSpan(1) is the starting time and
%               tSpan2 is the desired stopping time for the integration.
%             f The function handle for f(x,t)=dxdt over which integration
%               is to be performed. The output is NX1-dimensional.
%  initStepSize An optional initial step size (in t) to use for the
%               integration. If omitted or an empty matrix is passed, an
%               ad-hoc method described below is used to find an initial
%               step size.
% numRefinedSteps An optional parameter specifying the number of steps to
%               interpolate between each normal step. With high-order
%               formulas, the steps produced might not be very nice to
%               use for plotting, so setting this to a value above 0 adds
%               more intermediate values (using the polynomial returned by
%               the RKInterpPolys function) so that plotted data can look
%               nicer. If omitted or an empty matrix is passed, the default
%               value of 0 (no extra/ interpolated steps) is used.
%               Interpolation is done at the main order of the Runge-Kutta
%               method.
% order,solutionChoice  A pair of optional parameters that specify the
%               highest order of the embedded Runge-Kutta pair to use as
%               well as the specific algorithm to use. Details are given in
%               the comments to the RungeKStep function. If omitted or
%               empty matrices are passed, the default order of 5 is used
%               and the default solutionChoice of 0 is used.
%        RelTol The maximum relative error tolerance allowed, a
%               positive scalar (its use is explained in more detail
%               below). If omitted or an empty matrix is passed, the
%               default value of 1e-3 is used.
%        AbsTol The absolute error tolerance allowed, a positive scalar, of
%               the same for all components of x, or a positive NX1 vector.
%               If omitted or an empty matrix is passed, the default value
%               of 1e-6 is used.
%      maxSteps The maximum allowable number of steps to perform the
%               integration. If omitted, the default of 1024 is used.
%
%OUTPUTS: xVals The NXnumSteps set of values of x along the path.
%               xVals(:,1) is xStart and xVals(:,end) is the value at final
%               time tSpan(2). If the Runge-Kutta integration failed, e.g.
%               due to encountering a NaN or being unable to get a
%               sufficient step size, an empty matrix is returned.
%         tVals A numStepsX1 vector of values of t corresponding to the
%               values of x in xVals. tVals(1) is equal to tSpan(1) and
%               tVals(end) is equal to tSpan(2). If integration fails, an
%               empty matrix is returned.
%      dxdtVals A numStepsXN array containing derivatives of x evaluated at
%               the times in tVals. The derivatives are just the result of
%               evaluating f(x,t) at each point in (xVals,tVals). This
%               combined with xVals and tVals could be used in functions to
%               perform Hermite interpolation. If integration fails, an
%               empty matrix is returned.
%      exitCode A code indicating how the algorithm terminated. Possible
%               values are
%               0: Integration was successful.
%               1: Unable to get a small enough step size.
%               2: Maximum number of steps reached without completion.
%               3: Non-finite number encountered.
%  nextStepSize The step size (in t) that would be used for the next
%               integration step. If empty, then the integration ended
%               because the step size was too small.
%
%The basic idea behind adaptive stepsize control can be ascertained from
%Chapter 5.2 of [3]. However, the implementation here has a number of
%changes. For example, unlike the basic algorithm of [3], depending on the
%user's choice of order and solutionChoice, the main Runge-Kutta formula,
%which is the one whose result is saved, might have a higher or lower order
%than the subsidiary formula, which is only used to check the error to
%determine how the stepsize should change. The smallest order of the two is
%the order that is used in formulae for adaptively changing the step size.
%
%The algorithm needs an initial stepsize to start. If the user does not
%provide an intial stepsize, then the initial stepsize is just set to the
%width of the integration region in t divided by the maximum number of
%steps. That is, initStepSize = (tSpan(2)-tSpan(1))/maxSteps;
%
%Once an initial stepsize has been determined, the function RungeKStep is
%used to perform a single step of the Runge-Kutta algorithm. Then one must
%determine whether the stepsize was small enough. If the stepsize is small
%enough, then the stepsize should still be adaptively changed for the nest
%step. The procedure for updating the stepsize is essentially the procedure
%for non-stiff equations from [2], Equation 2.3). The method has been
%slightly modified in that the error being adjusted is not just the
%absolute error, as explained below. Also, the maximum stepsize is limited
%to 1/5 the integration region and increases in the stepsize are limited
%so as not to be too large.
%
%We would like the error in each step to be less than the absolute and the
%relative error tolerances. For an absolute error tolerance, this means
%that
%abs(xPredMain-xPredSubsid)<AbsTol
%for all components of the predicted value xPred in a particular step of
%size deltaT. For a relative error tolerance, this means that
%abs(xPredMain-xPredSubsid)./max(max(abs(xPred),abs(xCur)),RelTol)<RelTol
%for all components. The extra RelTol in the denominator on the left helps
%deal with issues when the true value is close to zero.
%
%The requirements for absolute and relative tolerances can be combined such
%that the absolute tolerance only comes into play when the numbers are
%sufficiently small that the relative tolerance becomes problematic. The
%two lines
% normFactor=max(max(abs(xPred),abs(xHatPred)),AbsTol/RelTol);
% theError=max(abs((xPred-xHatPred)./normFactor));
%Produce a single quantity, theError, which, if less than RelTol, means
%that the stepsize was too large. If the denominator of the normFactor is
%too small, then it is just AbsTol/RelTol, in which case the RelTol values
%on both sides of the equal sign cancel. And the comparison is just the
%comparison for the absolute tolerance (looking only at the largest
%error component). Thus theError term uses a threshold factor of
%AbsTol/RelTol. This is used as a similar threshold ratio when creating an
%ad-hoc initial estimate.
%
%REFERENCES:
%[1] I. Gladwell, L. F. Shampine, and R. W. Brankin, "Automatic selection
%    of the initial step size for an ODE solver," Journal of Computational
%    and Applied Mathematics, vol. 18, no. 2, pp. 175-192, May 1987.
%[2] G. Hall, "A new stepsize strategy for explicit Runge-Kutta codes,"
%    Advances in Computational Mathematics, vol. 3, no. 4, pp. 343-352, May
%    1995.
%[3] R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Boston, MA:
%    Brooks/ Cole, 2011.
%
%February 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<10)
    maxSteps=1024;
end

if(nargin<9||isempty(AbsTol))
    AbsTol=1e-6;
end

if(nargin<8||isempty(RelTol))
   RelTol=1e-3;
end

if(nargin<7||isempty(solutionChoice))
   solutionChoice=0;
end

if(nargin<6||isempty(order))
    order=5;
end

if(nargin<5||isempty(numRefinedSteps))
    numRefinedSteps=0;
end

xDim=size(xStart,1);
tStart=tSpan(1);
tEnd=tSpan(2);
tDiff=diff(tSpan);
tDiffMag=abs(tDiff);
deltaTSign=sign(tDiff);
%The sign is separated out so that one can use the algorithm running
%backwards as well as forwards in time.

%The maximum step size is arbitrarily set to 1/3 the total distance.
deltaTMaxMag=tDiffMag/3;

if(nargin<4||isempty(initStepSize))
    %If no initial step size is given, then just use the smallest uniform
    %step size for the given value of maxSteps.
    initStepSize=tDiffMag/maxSteps;
elseif(initStepSize<eps(tStart))
   error('The initial step size provided is too small.')
end

%Allocate the maximum possible space needed for the return parameters.
numSpaces=1+maxSteps*(numRefinedSteps+1);
xVals=zeros(xDim,numSpaces);
dxdtVals=zeros(xDim,numSpaces);
tVals=zeros(numSpaces,1);

xVals(:,1)=xStart;
tVals(1)=tStart;
dxdtVals(:,1)=f(xStart,tStart);
lastStep=1;

deltaT=deltaTSign*initStepSize;
for curStep=2:(maxSteps+1)
    xCur=xVals(:,lastStep);
    tCur=tVals(lastStep);
    dxdtCur=dxdtVals(:,lastStep);

    %The absolute value of the minimum allowable step size. It is set so
    %that the step must makes something of a difference compared to the
    %numerical precision.
    deltaTMinMag=2^4*eps(tCur);

    %If we would have overstepped the end, change deltaT to be the end.
    %Then, we can terminate at the end.
    if(deltaTSign*(tCur+deltaT)>tEnd*deltaTSign)
        deltaTMag=abs(tEnd-tCur);
        deltaT=deltaTMag*deltaTSign;
        %Allow the last step to be very small.
        deltaTMinMag=min(deltaTMinMag,deltaTMag);
    end

    [deltaTNew,xNew,tNew,k,dxdtCur,exitCode]=performOneAdaptiveRKStep(xCur,tCur,f,deltaT,deltaTMinMag,deltaTMaxMag,dxdtCur,order,solutionChoice,AbsTol,RelTol);

    %If an error occurred, then just return.
    if(exitCode~=0)
        xVals=[];
        tVals=[];
        dxdtVals=[];
        nextStepSize=abs(deltaTNew);
        return
    end

    dxdtVals(:,lastStep)=dxdtCur;

    if(numRefinedSteps~=0)
        %If extra interpolation between the steps is to be performed.
        %Get a Hermite interpolating polynomial over the step.
        [interpPolyA,interpPolyC]=RKInterpPolys(xCur,tCur,xNew,tNew,f,order,solutionChoice,k);
        %The interpolating polynomials take stapes parameterized by
        %a value between zero and 1 representing the distance from
        %tCur to tNew.
        deltaTTaken=tNew-tCur;

        deltaFraction=1/(numRefinedSteps+1);
        curStepFrac=deltaFraction;
        for curRStep=1:numRefinedSteps
            lastStep=lastStep+1;
            tVals(lastStep)=tVals(lastStep-1)+deltaFraction*deltaTTaken;

            xVals(:,lastStep)=polyValNewton(curStepFrac,interpPolyA,interpPolyC);
            dxdtVals(:,lastStep)=f(xVals(:,lastStep),tVals(lastStep));
            curStepFrac=curStepFrac+deltaFraction;
        end
    end

    lastStep=lastStep+1;
    xVals(:,lastStep)=xNew;
    tVals(lastStep)=tNew;
    dxdtVals(:,lastStep)=dxdtCur;
    deltaT=deltaTNew;

    if(tVals(lastStep)==tEnd)
        %If integration completed, reshape the outputs to match the actual
        %number of steps.
        xVals=xVals(:,1:lastStep);
        tVals=tVals(1:lastStep);
        dxdtVals=dxdtVals(:,1:lastStep);
        exitCode=0;
        nextStepSize=deltaTMag;
        return
    end
end

%If we get here, then the maximum number of iterations was reached without
%reaching the end.
xVals=[];
dxdtVals=[];
tVals=[];
exitCode=2;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text: %MRDIVIDE Mapping overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ n_data, a, x, fx ] = jacobi_dn_values ( n_data )

%*****************************************************************************80
%
%% JACOBI_DN_VALUES returns some values of the Jacobi elliptic function DN(A,X).
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      JacobiDN[ x, a ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, the parameter of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  a_vec = [ ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.0E+00, ...
     0.5E+00, ...
     0.5E+00, ...
     0.5E+00, ...
     0.5E+00, ...
     0.5E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00, ...
     1.0E+00 ];

  fx_vec = [ ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.9975093485144243E+00, ...
     0.9901483195224800E+00, ...
     0.9429724257773857E+00, ...
     0.8231610016315963E+00, ...
     0.7108610477840873E+00, ...
     0.9950207489532265E+00, ...
     0.9803279976447253E+00, ...
     0.8868188839700739E+00, ...
     0.6480542736638854E+00, ...
     0.2658022288340797E+00, ...
     0.3661899347368653E-01, ...
     0.9803279976447253E+00, ...
     0.8868188839700739E+00, ...
     0.6480542736638854E+00, ...
     0.2658022288340797E+00  ];

  x_vec = [ ...
      0.1E+00, ...
      0.2E+00, ...
      0.5E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      0.1E+00, ...
      0.2E+00, ...
      0.5E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      0.1E+00, ...
      0.2E+00, ...
      0.5E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      4.0E+00, ...
     -0.2E+00, ...
     -0.5E+00, ...
     -1.0E+00, ...
     -2.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
% RRT for Redundant manipulator example.
% by Olzhas Adiyatov
% 08/28/2013

map = struct('name', 'bench_redundant_3.mat', 'start_point', [0 0], 'goal_point', [35 35]);
max_iter = 10e3;
is_benchmark = false;
rand_seed = 40;
variant = 'FNRedundantManipulator';
result = rrt_star(map, max_iter, is_benchmark, rand_seed, variant);

author: olzhas