%    Copyright (C) 2012      Waves Audio LTD
%    Copyright (C) 2003-2008 Jean-Marc Valin
%
%    File: speex_mdf.m
%    Echo canceller based on the MDF algorithm (see below)
%
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions are
%    met:
%
%    1. Redistributions of source code must retain the above copyright notice,
%    this list of conditions and the following disclaimer.
%
%    2. Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%
%    3. The name of the author may not be used to endorse or promote products
%    derived from this software without specific prior written permission.
%
%    THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
%    IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
%    OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
%    DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
%    INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
%    (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
%    SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
%    HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
%    STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
%    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%
%    Notes from original mdf.c:
%
%    The echo canceller is based on the MDF algorithm described in:
%
%    J. S. Soo, K. K. Pang Multidelay block frequency adaptive filter,
%    IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-38, No. 2,
%    February 1990.
%
%    We use the Alternatively Updated MDF (AUMDF) variant. Robustness to
%    double-talk is achieved using a variable learning rate as described in:
%
%    Valin, J.-M., On Adjusting the Learning Rate in Frequency Domain Echo
%    Cancellation With Double-Talk. IEEE Transactions on Audio,
%    Speech and Language Processing, Vol. 15, No. 3, pp. 1030-1034, 2007.
%    http://people.xiph.org/~jm/papers/valin_taslp2006.pdf
%
%    There is no explicit double-talk detection, but a continuous variation
%    in the learning rate based on residual echo, double-talk and background
%    noise.
%
%    Another kludge that seems to work good: when performing the weight
%    update, we only move half the way toward the "goal" this seems to
%    reduce the effect of quantization noise in the update phase. This
%    can be seen as applying a gradient descent on a "soft constraint"
%    instead of having a hard constraint.
%
%    Notes for this file:
%
%    Usage:
%
%       speex_mdf_out = speex_mdf(Fs, u, d, filter_length, frame_size, dbg_var_name);
%
%       Fs                  sample rate
%       u                   speaker signal, column vector in range [-1; 1]
%       d                   microphone signal, column vector in range [-1; 1]
%       filter_length       typically 250ms, i.e. 4096 @ 16k FS
%                           must be a power of 2
%       frame_size          typically 8ms, i.e. 128 @ 16k Fs
%                           must be a power of 2
%       dbg_var_name        internal state variable name to trace.
%                           Default: 'st.leak_estimate'.
%
%    Jonathan Rouach <jonr@waves.com>
%

function  speex_mdf_out = speex_mdf(Fs, u, d, filter_length, frame_size, dbg_var_name)

fprintf('Starting Speex MDF (PBFDAF) algorithm.\n');

st = speex_echo_state_init_mc_mdf(frame_size, filter_length, 1, 1, Fs);

% which variable to trace
if nargin<6
    dbg_var_name = 'st.leak_estimate';
end
dbg = init_dbg(st, length(u));

[e, dbg] = main_loop(st, float_to_short(u), float_to_short(d), dbg);

speex_mdf_out.e = e/32768.0;
speex_mdf_out.var1 = dbg.var1;

    function x = float_to_short(x)
        x = x*32768.0;
        x(x< -32767.5) = -32768;
        x(x>  32766.5) =  32767;
        x = floor(0.5+x);
    end

    function [e, dbg] = main_loop(st, u, d, dbg)

        e = zeros(size(u));
        y = zeros(size(u));

        % prepare waitbar
%         try h_wb = waitbar(0, 'Processing...'); catch; end
        end_point = length(u);

        for n = 1:st.frame_size:end_point
            nStep = floor(n/st.frame_size)+1;

%             if mod(nStep, 128)==0 && update_waitbar_check_wasclosed(h_wb, n, end_point, st.sampling_rate)
%                 break;
%             end

            u_frame = u(n:n+st.frame_size-1);
            d_frame = d(n:n+st.frame_size-1);

            [out, st] = speex_echo_cancellation_mdf(st, d_frame, u_frame);

            e(n:n+st.frame_size-1) = out;
            y(n:n+st.frame_size-1) = d_frame - out;
            dbg.var1(:, nStep) = reshape( eval(dbg_var_name),  numel(eval(dbg_var_name)), 1);

        end

        try close(h_wb); catch; end

    end
    function st = speex_echo_state_init_mc_mdf(frame_size, filter_length, nb_mic, nb_speakers, sample_rate)

        st.K = nb_speakers;
        st.C = nb_mic;
        C=st.C;
        K=st.K;

        st.frame_size = frame_size;
        st.window_size = 2*frame_size;
        N = st.window_size;
        st.M = fix((filter_length+st.frame_size-1)/frame_size);
        M = st.M;
        st.cancel_count=0;
        st.sum_adapt = 0;
        st.saturated = 0;
        st.screwed_up = 0;

        %    /* This is the default sampling rate */
        st.sampling_rate = sample_rate;
        st.spec_average = (st.frame_size)/( st.sampling_rate);
        st.beta0 = (2.0*st.frame_size)/st.sampling_rate;
        st.beta_max = (.5*st.frame_size)/st.sampling_rate;
        st.leak_estimate = 0;

        st.e = zeros(N, C);
        st.x = zeros(N, K);
        st.input = zeros(st.frame_size, C);
        st.y = zeros(N, C);
        st.last_y = zeros(N, C);
        st.Yf = zeros(st.frame_size+1, 1);
        st.Rf = zeros(st.frame_size+1, 1);
        st.Xf = zeros(st.frame_size+1, 1);
        st.Yh = zeros(st.frame_size+1, 1);
        st.Eh = zeros(st.frame_size+1, 1);

        st.X = zeros(N, K, M+1);
        st.Y = zeros(N, C);
        st.E = zeros(N, C);
        st.W = zeros(N, K, M, C);
        st.foreground = zeros(N, K, M, C);
        st.PHI = zeros(frame_size+1, 1);
        st.power = zeros(frame_size+1, 1);
        st.power_1 = ones((frame_size+1), 1);
        st.window = zeros(N, 1);
        st.prop = zeros(M, 1);
        st.wtmp = zeros(N, 1);

        st.window = .5-.5*cos(2*pi*((1:N)'-1)/N);

        %        /* Ratio of ~10 between adaptation rate of first and last block */
        decay = exp(-2.4/M);
        st.prop(1, 1) = .7;
        for i=2:M
            st.prop(i, 1) = st.prop(i-1, 1) * decay;
        end

        st.prop = (.8 * st.prop)./sum(st.prop);

        st.memX = zeros(K, 1);
        st.memD = zeros(C, 1);
        st.memE = zeros(C, 1);
        st.preemph = .9;
        if (st.sampling_rate<12000)
            st.notch_radius = .9;
        elseif (st.sampling_rate<24000)
            st.notch_radius = .982;
        else
            st.notch_radius = .992;
        end

        st.notch_mem = zeros(2*C, 1);
        st.adapted = 0;
        st.Pey = 1;
        st.Pyy = 1;

        st.Davg1 = 0; st.Davg2 = 0;
        st.Dvar1 = 0; st.Dvar2 = 0;
    end

    function dbg = init_dbg(st, len)
        dbg.var1 = zeros(numel(eval(dbg_var_name)), fix(len/st.frame_size));
    end

    function [out, st] = speex_echo_cancellation_mdf(st, in, far_end)

        N = st.window_size;
        M = st.M;
        C = st.C;
        K = st.K;

        Pey_cur = 1;
        Pyy_cur = 1;

        out = zeros(st.frame_size, C);

        st.cancel_count = st.cancel_count + 1;

        ss=.35/M;
        ss_1 = 1-ss;

        for chan = 1:C
            % Apply a notch filter to make sure DC doesn't end up causing problems
            [st.input(:, chan), st.notch_mem(:, chan)] = filter_dc_notch16(in(:, chan), st.notch_radius, st.frame_size, st.notch_mem(:, chan));
            % Copy input data to buffer and apply pre-emphasis
            for i=1:st.frame_size
                tmp32 = st.input(i, chan)- (st.preemph* st.memD(chan));
                st.memD(chan) = st.input(i, chan);
                st.input(i, chan) = tmp32;
            end
        end

        for speak = 1:K
            for i =1:st.frame_size
                st.x(i, speak) = st.x(i+st.frame_size, speak);
                tmp32 = far_end(i, speak) - st.preemph * st.memX(speak);
                st.x(i+st.frame_size, speak) = tmp32;
                st.memX(speak) = far_end(i, speak);
            end
        end

        % Shift memory
        st.X = circshift(st.X, [0, 0, 1]);

        for speak = 1:K
            %  Convert x (echo input) to frequency domain
            % MATLAB_MATCH: we divide by N to get values as in speex
            st.X(:, speak, 1) = fft(st.x(:, speak)) /N;
        end

        Sxx = 0;
        for speak = 1:K
            Sxx = Sxx + sum(st.x(st.frame_size+1:end, speak).^2);
            st.Xf = abs(st.X(1:st.frame_size+1, speak, 1)).^2;
        end

        Sff = 0;
        for chan = 1:C

            %  Compute foreground filter
            st.Y(:, chan) = 0;
            for speak=1:K
                for j=1:M
                    st.Y(:, chan) = st.Y(:, chan) + st.X(:, speak, j) .* st.foreground(:, speak, j, chan);
                end
            end
            % MATLAB_MATCH: we multiply by N to get values as in speex
            st.e(:, chan) = ifft(st.Y(:, chan)) * N;
            st.e(1:st.frame_size, chan) = st.input(:, chan) - st.e(st.frame_size+1:end, chan);
            % st.e : [out foreground | leak foreground ]
            Sff = Sff + sum(abs(st.e(1:st.frame_size, chan)).^2);

        end

        % Adjust proportional adaption rate */
        if (st.adapted)
            st.prop = mdf_adjust_prop (st.W, N, M, C, K);
        end

        % Compute weight gradient */
        if (st.saturated == 0)
            for chan = 1:C
                for speak = 1:K
                    for j=M:-1:1
                        st.PHI = [st.power_1; st.power_1(end-1:-1:2)] .* st.prop(j) .* conj(st.X(:, speak, (j+1))) .* st.E(:, chan);
                        st.W(:, j) = st.W(:, j) + st.PHI;
                    end
                end
            end
        else
            st.saturated = st.saturated -1;
        end

        %FIXME: MC conversion required */
        % Update weight to prevent circular convolution (MDF / AUMDF)
        for chan = 1:C
            for speak = 1:K
                for j = 1:M
                    % This is a variant of the Alternatively Updated MDF (AUMDF) */
                    % Remove the "if" to make this an MDF filter */
                    if (j==1 || mod(2+st.cancel_count,(M-1)) == j)
                        st.wtmp = ifft(st.W(:, speak, j, chan));
                        st.wtmp(st.frame_size+1:N) = 0;
                        st.W(:, speak, j, chan) = fft(st.wtmp);
                    end
                end
            end
        end

        % So we can use power_spectrum_accum */
        st.Yf = zeros(st.frame_size+1, 1);
        st.Rf = zeros(st.frame_size+1, 1);
        st.Xf = zeros(st.frame_size+1, 1);

        Dbf = 0;

        for chan = 1:C
            st.Y(:, chan) = 0;
            for speak=1:K
                for j=1:M
                    st.Y(:, chan) = st.Y(:, chan) + st.X(:, speak, j) .* st.W(:, speak, j, chan);
                end
            end
            % MATLAB_MATCH: we multiply by N to get values as in speex
            st.y(:,chan) = ifft(st.Y(:,chan)) * N;
            % st.y : [ ~ | leak background ]
        end

        See = 0;

        % Difference in response, this is used to estimate the variance of our residual power estimate */
        for chan = 1:C
            st.e(1:st.frame_size, chan) = st.e(st.frame_size+1:N, chan) - st.y(st.frame_size+1:N, chan);
            Dbf = Dbf + 10 + sum(abs(st.e(1:st.frame_size, chan)).^2);
            st.e(1:st.frame_size, chan) = st.input(:, chan) - st.y(st.frame_size+1:N, chan);
            % st.e : [ out background | leak foreground ]
           See = See + sum(abs(st.e(1:st.frame_size, chan)).^2);
        end

        % Logic for updating the foreground filter */

        % For two time windows, compute the mean of the energy difference, as well as the variance */
        VAR1_UPDATE = .5;
        VAR2_UPDATE = .25;
        VAR_BACKTRACK = 4;
        MIN_LEAK = .005;

        st.Davg1 = .6*st.Davg1 + .4*(Sff-See);
        st.Davg2 = .85*st.Davg2 + .15*(Sff-See);
        st.Dvar1 = .36*st.Dvar1 + .16*Sff*Dbf;
        st.Dvar2 = .7225*st.Dvar2 + .0225*Sff*Dbf;

        update_foreground = 0;

        % Check if we have a statistically significant reduction in the residual echo */
        % Note that this is *not* Gaussian, so we need to be careful about the longer tail */
        if (Sff-See)*abs(Sff-See) > (Sff*Dbf)
            update_foreground = 1;
        elseif (st.Davg1* abs(st.Davg1) > (VAR1_UPDATE*st.Dvar1))
            update_foreground = 1;
        elseif (st.Davg2* abs(st.Davg2) > (VAR2_UPDATE*(st.Dvar2)))
            update_foreground = 1;
        end

        % Do we update? */
        if (update_foreground)

            st.Davg1 = 0;
            st.Davg2 = 0;
            st.Dvar1 = 0;
            st.Dvar2 = 0;
            st.foreground = st.W;
            % Apply a smooth transition so as to not introduce blocking artifacts */
            for chan = 1:C
                st.e(st.frame_size+1:N, chan) = (st.window(st.frame_size+1:N) .* st.e(st.frame_size+1:N, chan)) + (st.window(1:st.frame_size) .* st.y(st.frame_size+1:N, chan));
            end
        else
            reset_background=0;
            % Otherwise, check if the background filter is significantly worse */

            if (-(Sff-See)*abs(Sff-See)> VAR_BACKTRACK*(Sff*Dbf))
                reset_background = 1;
            end
            if ((-st.Davg1 * abs(st.Davg1))> (VAR_BACKTRACK*st.Dvar1))
                reset_background = 1;
            end
            if ((-st.Davg2* abs(st.Davg2))> (VAR_BACKTRACK*st.Dvar2))
                reset_background = 1;
            end

            if (reset_background)

                % Copy foreground filter to background filter */
                st.W = st.foreground;

                % We also need to copy the output so as to get correct adaptation */
                for chan = 1:C
                    st.y(st.frame_size+1:N, chan) = st.e(st.frame_size+1:N, chan);
                    st.e(1:st.frame_size, chan) = st.input(:, chan) - st.y(st.frame_size+1:N, chan);
                end

                See = Sff;
                st.Davg1 = 0;
                st.Davg2 = 0;
                st.Dvar1 = 0;
                st.Dvar2 = 0;
            end
        end

        Sey = 0;
        Syy = 0;
        Sdd = 0;

        for chan = 1:C

            % Compute error signal (for the output with de-emphasis) */
            for i=1:st.frame_size
                tmp_out = st.input(i, chan)- st.e(i+st.frame_size, chan);
                tmp_out = tmp_out + st.preemph * st.memE(chan);
                %  This is an arbitrary test for saturation in the microphone signal */
                if (in(i,chan) <= -32000 || in(i,chan) >= 32000)
                    if (st.saturated == 0)
                        st.saturated = 1;
                    end
                end
                out(i, chan) = tmp_out;
                st.memE(chan) = tmp_out;
            end

            % Compute error signal (filter update version) */
            st.e(st.frame_size+1:N, chan) = st.e(1:st.frame_size, chan);
            st.e(1:st.frame_size, chan) = 0;
            % st.e : [ zeros | out background ]

            % Compute a bunch of correlations */
            % FIXME: bad merge */
            Sey = Sey + sum(st.e(st.frame_size+1:N, chan) .* st.y(st.frame_size+1:N, chan));
            Syy = Syy + sum(st.y(st.frame_size+1:N, chan).^2);
            Sdd = Sdd + sum(st.input.^2);

            % Convert error to frequency domain */
            % MATLAB_MATCH: we divide by N to get values as in speex
            st.E = fft(st.e) / N;

            st.y(1:st.frame_size, chan) = 0;
            % MATLAB_MATCH: we divide by N to get values as in speex
            st.Y = fft(st.y) / N;

            % Compute power spectrum of echo (X), error (E) and filter response (Y) */
            st.Rf = abs(st.E(1:st.frame_size+1,chan)).^2;
            st.Yf = abs(st.Y(1:st.frame_size+1,chan)).^2;
        end

        % Do some sanity check */
        if (~(Syy>=0 && Sxx>=0 && See >= 0))
            % Things have gone really bad */
            st.screwed_up = st.screwed_up + 50;
            out = out*0;
        elseif Sff > Sdd+ N*10000
            % AEC seems to add lots of echo instead of removing it, let's see if it will improve */
            st.screwed_up = st.screwed_up + 1;
        else
            % Everything's fine */
            st.screwed_up=0;
        end

        if (st.screwed_up>=50)
            disp('Screwed up, full reset');
            st = speex_echo_state_reset_mdf(st);
        end

        % Add a small noise floor to make sure not to have problems when dividing */
        See = max(See, N* 100);

        for speak = 1:K
            Sxx = Sxx + sum(st.x(st.frame_size+1:end, speak).^2);
            st.Xf = abs(st.X(1:st.frame_size+1, speak, 1)).^2;
        end

        % Smooth far end energy estimate over time */
        st.power = ss_1*st.power+ 1 + ss*st.Xf;

        % Compute filtered spectra and (cross-)correlations */

        Eh_cur = st.Rf - st.Eh;
        Yh_cur = st.Yf - st.Yh;
        Pey_cur = Pey_cur + sum(Eh_cur.*Yh_cur) ;
        Pyy_cur = Pyy_cur + sum(Yh_cur.^2);
        st.Eh = (1-st.spec_average)*st.Eh + st.spec_average*st.Rf;
        st.Yh = (1-st.spec_average)*st.Yh + st.spec_average*st.Yf;

        Pyy = sqrt(Pyy_cur);
        Pey = Pey_cur/Pyy;

        % Compute correlation updatete rate */
        tmp32 = st.beta0*Syy;
        if (tmp32 > st.beta_max*See)
            tmp32 = st.beta_max*See;
        end
        alpha = tmp32/ See;
        alpha_1 = 1- alpha;

        % Update correlations (recursive average) */
        st.Pey = alpha_1*st.Pey + alpha*Pey;
        st.Pyy = alpha_1*st.Pyy + alpha*Pyy;

        if st.Pyy<1
            st.Pyy =1;
        end

        % We don't really hope to get better than 33 dB (MIN_LEAK-3dB) attenuation anyway */
        if st.Pey< MIN_LEAK * st.Pyy
            st.Pey = MIN_LEAK * st.Pyy;
        end

        if (st.Pey> st.Pyy)
            st.Pey = st.Pyy;
        end

        % leak_estimate is the linear regression result */
        st.leak_estimate = st.Pey/st.Pyy;

        % This looks like a stupid bug, but it's right (because we convert from Q14 to Q15) */
        if (st.leak_estimate > 16383)
            st.leak_estimate = 32767;
        end

        % Compute Residual to Error Ratio */
        RER = (.0001*Sxx + 3.*st.leak_estimate*Syy) / See;
        % Check for y in e (lower bound on RER) */
        if (RER < Sey*Sey/(1+See*Syy))
            RER = Sey*Sey/(1+See*Syy);
        end
        if (RER > .5)
            RER = .5;
        end

        % We consider that the filter has had minimal adaptation if the following is true*/
        if (~st.adapted && st.sum_adapt > M && st.leak_estimate*Syy > .03*Syy)
            st.adapted = 1;
        end

        if (st.adapted)
            % Normal learning rate calculation once we're past the minimal adaptation phase */
            for i=1:st.frame_size+1

                % Compute frequency-domain adaptation mask */
                r = st.leak_estimate*st.Yf(i);
                e = st.Rf(i)+1;
                if (r>.5*e)
                    r = .5*e;
                end
                r = 0.7*r + 0.3*(RER*e);
                %st.power_1[i] = adapt_rate*r/(e*(1+st.power[i]));*/
                st.power_1(i) = (r/(e*st.power(i)+10));
            end
        else
            % Temporary adaption rate if filter is not yet adapted enough */
            adapt_rate=0;

            if (Sxx > N* 1000)

                tmp32 = 0.25* Sxx;
                if (tmp32 > .25*See)
                    tmp32 = .25*See;
                end
                adapt_rate = tmp32/ See;
            end
            st.power_1 = adapt_rate./(st.power+10);

            % How much have we adapted so far? */
            st.sum_adapt = st.sum_adapt+adapt_rate;
        end

        % FIXME: MC conversion required */
        st.last_y(1:st.frame_size) = st.last_y(st.frame_size+1:N);
        if (st.adapted)
            % If the filter is adapted, take the filtered echo */
            st.last_y(st.frame_size+1:N) = in-out;
        end

    end

    function [out,mem] = filter_dc_notch16(in, radius, len, mem)
        out = zeros(size(in));
        den2 = radius*radius + .7*(1-radius)*(1-radius);
        for i=1:len
            vin = in(i);
            vout = mem(1) + vin;
            mem(1) = mem(2) + 2*(-vin + radius*vout);
            mem(2) = vin - (den2*vout);
            out(i) = radius*vout;
        end

    end

    function prop = mdf_adjust_prop(W, N, M, C, K)
        prop = zeros(M,1);
        for i=1:M
            tmp = 1;
            for chan=1:C
                for speak=1:K
                    tmp = tmp + sum(abs(W(1:N/2+1, K, i, C)).^2);
                end
            end
            prop(i) = sqrt(tmp);
        end
        max_sum = max(prop, 1);
        prop = prop + .1*max_sum;
        prop_sum = 1+sum(prop);
        prop = .99*prop / prop_sum;
    end

    % Resets echo canceller state */
    function st = speex_echo_state_reset_mdf(st)

        st.cancel_count=0;
        st.screwed_up = 0;
        N = st.window_size;
        M = st.M;
        C=st.C;
        K=st.K;

        st.e = zeros(N, C);
        st.x = zeros(N, K);
        st.input = zeros(st.frame_size, C);
        st.y = zeros(N, C);
        st.last_y = zeros(N, C);
        st.Yf = zeros(st.frame_size+1, 1);
        st.Rf = zeros(st.frame_size+1, 1);
        st.Xf = zeros(st.frame_size+1, 1);
        st.Yh = zeros(st.frame_size+1, 1);
        st.Eh = zeros(st.frame_size+1, 1);

        st.X = zeros(N, K, M+1);
        st.Y = zeros(N, C);
        st.E = zeros(N, C);
        st.W = zeros(N, K, M, C);
        st.foreground = zeros(N, K, M, C);
        st.PHI = zeros(N, 1);
        st.power = zeros(st.frame_size+1, 1);
        st.power_1 = ones((st.frame_size+1), 1);
        st.window = zeros(N, 1);
        st.prop = zeros(M, 1);
        st.wtmp = zeros(N, 1);

        st.memX = zeros(K, 1);
        st.memD = zeros(C, 1);
        st.memE = zeros(C, 1);

        st.saturated = 0;
        st.adapted = 0;
        st.sum_adapt = 0;
        st.Pey = 1;
        st.Pyy = 1;
        st.Davg1 = 0;
        st.Davg2 = 0;
        st.Dvar1 = 0;
        st.Dvar2 = 0;

    end

    function was_closed = update_waitbar_check_wasclosed(h, n, end_point, Fs)
        was_closed = 0;

        % update waitbar
        try
            waitbar(n/end_point, h, ['Processing... ', num2str(n/Fs, '%.2f'), 's / ', num2str(end_point/Fs, '%.2f'), 's' ]);
        catch ME
            % if it's no longer there (closed by user)
            if (strcmp(ME.identifier(1:length('MATLAB:waitbar:')), 'MATLAB:waitbar:'))
                was_closed = 1; % then get out of the loop
            end
        end

    end

end


text:
function pass = test_constructor2( pref )
% This tests the chebfun3 constructor.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e2 * pref.cheb3Prefs.chebfun3eps;

% Test building Chebfun3 objects from sample data:
exactCoeffs = rand(3,3,3);
constructorCoeffs = chebcoeffs3(chebfun3(exactCoeffs, 'coeffs'));
pass(1) = norm(exactCoeffs(:) - constructorCoeffs(:)) < 10*tol;

exactCoeffs = rand(4,4,4);
constructorCoeffs = chebcoeffs3(chebfun3(exactCoeffs, 'coeffs'));
pass(2) = norm(exactCoeffs(:) - constructorCoeffs(:)) < 10*tol;

% Make a chebfun3 based on TRIGTECH by calling the 'trig' flag:
f1 = chebfun3(@(x,y,z) cos(pi*x).*sin(pi*y).*cos(pi*z),'trig');
f2 = chebfun3(@(x,y,z) cos(pi*x).*sin(pi*y).*cos(pi*z), ...
    [-1 1 -1 1 -1 1], 'trig');
pass(3) = norm( f1 - f2 ) < tol;

f1 = chebfun3(@(x,y,z) cos(pi*cos(pi*x) + pi*sin(pi*y) + ...
    pi*cos(pi*z)),'trig');
f2 = chebfun3(@(x,y,z) cos(pi*cos(pi*x) + pi*sin(pi*y) + ...
    pi*cos(pi*z)), [-1 1 -1 1 -1 1], 'trig');
pass(4) = ( norm(f1 -f2) < 10*tol );

% Check underlying tech is a TRIGTECH:
colTech = get(f1.cols.funs{1}, 'tech');
rowTech = get(f1.rows.funs{1}, 'tech');
tubeTech = get(f1.tubes.funs{1}, 'tech');
pass(5) = isa(colTech(), 'trigtech');
pass(6) = isa(rowTech(), 'trigtech');
pass(7) = isa(tubeTech(), 'trigtech');

% Make sure the 'periodic' flag works as well:
f1 = chebfun3(@(x,y,z) cos(pi*cos(pi*x) + pi*sin(pi*y) + ...
    pi*cos(pi*z)),'periodic');
f2 = chebfun3(@(x,y,z) cos(pi*cos(pi*x) + pi*sin(pi*y) + ...
    pi*cos(pi*z)), [-1 1 -1 1 -1 1], 'periodic');
pass(8) = ( norm(f1 -f2) < 10*tol );

% Check underlying tech is a TRIGTECH:
colTech = get(f1.cols.funs{1}, 'tech');
rowTech = get(f1.rows.funs{1}, 'tech');
tubeTech = get(f1.tubes.funs{1}, 'tech');
pass(9) = isa(colTech(), 'trigtech');
pass(10) = isa(rowTech(), 'trigtech');
pass(11) = isa(tubeTech(), 'trigtech') ;

% Test making a chebfun3 from a scalar coefficient:
f = chebfun3(1, 'coeffs');
pass(12) = norm(f - 1) < tol ;

% Test passing an 'eps' value. Make sure it affects both rank and lengths.
ff = @(x,y,z) -x.*sin(sqrt(x+15+2*y+3*z));
f = chebfun3(ff);
fEps = chebfun3(ff, 'eps', 1e-8);
[m, n, p] = length(f);
[mEps, nEps, pEps] = length(fEps);
pass(13) = mEps < m && nEps < n && pEps < p ;

[r1, r2, r3] = rank(f);
[r1Eps, r2Eps, r3Eps] = rank(fEps);
pass(14) = r1Eps < r1 && r2Eps < r2 && r3Eps < r3 ;

% Construct from a string of constant type:
f = chebfun3('pi');
pass(15) = f(0,0,0) - pi < tol;

% Construct from a string of a univariate function
f = chebfun3('cos(alpha)');
pass(16) = f(0,0,0) - cos(0) < tol;

% Construct from a string of a bivariate function
f = chebfun3('x+y');
pass(17) = f(0.25,0.5,0) - 0.75 < tol;

% Construct from a string of a trivariate function
f = chebfun3('cos(x+y+z)');
pass(18) = f(0.25,0.5,-0.3) - cos(0.25+0.5-0.3) < tol;

% Test that the 'equi' flag outputs an error message if used for
% adaptive construction of a chebfun3.
try
    f = chebfun3(@(x, y, z) z.*x.^2.*exp(sin(x + y)), 'equi');
    pass(19) = false;
catch ME
    pass(19) = strcmp(ME.identifier, 'CHEBFUN:CHEBFUN3:constructor:equi');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Q=QTaylor(deltaT,xCur,curT,D,dadx,method)
%%QTAYLOR Get the process noise covariance under a nonlinear continuous-
%         time random process specified by the Langevin equation forward in
%         time by a step-size of deltaT using a Taylor scheme with additive
%         noise.
%
%INPUTS: deltaT The size of the single step over which to generate the
%               process noise covariance matrix.
%          xCur The initial target state at time curT.
%          curT The time of the initial state xCur.
%             D The diffusion function in the continuous-time stochastic
%               dynamic model. It takes the state and a time variable as
%               its arguments, D(x,t). Since the process noise must be
%               additive, D(x,t) should not depend on x.
%          dadx A xDimXxDim matrix of the derivative of the drift function
%               with respect to the state at state xCur and time curT.
%               This can also be a function that takes the state and a time
%               variable as its arguments, dadx(x,t).
%        method Set to 0 for the shorter Euler-Maruyama expansion, 1 for
%               the order 1.5 strong Taylor (default), and 2 for the order
%               2.0 weak Taylor.
%
%OUTPUTS: Q The process noise covariance matrix under a nonlinear
%           continuous-time random process specified by the Langevin
%           equation forward in time by a step-size of deltaT using a
%           strong Taylor scheme with additive noise.
%
%The process noise covariance matrix is derived from the stochastic order
%1.5 Taylor scheme with additive noise described in 10.4 of [1].
%
%REFERENCES:
%[1] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic
%    Differential Equations. Berlin: Springer, 1999.
%
%April 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isa(dadx,'function_handle'))
    dadx=dadx(xCur,curT);
end
if(nargin<6||isempty(method))
    method=1;
end

DCur=D(xCur,curT);

if(method==0)
    Q=deltaT*(DCur*DCur.');
else
    term1=deltaT*(DCur*DCur.');
    term2=(deltaT^2/2)*(DCur*(dadx*DCur).'+(DCur*(dadx*DCur).').');
    if(method==1)
        term3=(deltaT^3/3)*(dadx*DCur)*((dadx*DCur).');
    elseif(method==2)
        term3=(deltaT^3/4)*(dadx*DCur)*((dadx*DCur).');
    end

    Q=term1+term2+term3;
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:
function out = degree
% x * degree =  x in standard arc measure

out = pi/180;

author: mtex-toolbox

repo: mtex

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


text:
function pstruct = tapas_kf_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

pstruct.g_0       = pvec(1);
pstruct.mu_0      = pvec(2);
pstruct.om        = pvec(3);
pstruct.pi_u      = pvec(4);

return;

author: translationalneuromodeling

repo: tapas


text:
%{
clc;
lrs_load_conf;
load(fullfile(lrs_conf.lrs_dir,'dataset','trafficdb','traffic_patches.mat'));
V = im2double(imgdb{100});
[M,m,n,p] = convert_video3d_to_2d(V);
%}

D = M;
L = zeros(size(D));
S = zeros(size(D));
E = zeros(size(D));

%CM = ones(size(D));
%CM = zeros(size(D));
[numr,numc] = size(D);
CM = randi([0 1],numr,numc); % simulated confidence map (binary matrix)

Omega = find(CM ~= 0);
[I, J] = ind2sub([numr numc],Omega);

% Motion-Assisted Matrix Restoration (Ye et al. 2015)
if(strcmp(algorithm_id,'MAMR'))
  lambda = 1;
  [L,S,iter1] = core_MAMR(D,lambda,I,J);
end

% Robust Motion-Assisted Matrix Restoration (Ye et al. 2015)
if(strcmp(algorithm_id,'RMAMR'))
  lambda = 10;
  [L,S,E,iter1] = core_RMAMR(D,lambda,I,J);
end

M_hat = L + S + E;
error = norm(M_hat(:)-M(:))/norm(M(:));
disp(['Error: ' num2str(error)]);

%{
show_2dvideo(M,m,n);
show_2dvideo(M.*CM,m,n);
show_2dvideo(L,m,n);
show_2dvideo(S,m,n);
show_2dvideo(Z,m,n);
%}

author: andrewssobral


text:
function ok = StartSIFT(interactive,nosplash)
% This function initializes the Source Information Flow Toolbox (SIFT)
% Author: Tim Mullen, 2011, SCCN/INC/UCSD

ok = false;

if nargin<1
    interactive = false;
end
if nargin<2
    nosplash = false;
end

fprintf('Initializing SIFT...\n');

% temporary hack to remove measure projection toolbox from path
% (MPT interferes with SIFT)
% if exist('measure_projection')
%     hlp_removeMeasureProjection;
% end

if ispc, SLASH = '\'; else SLASH = '/'; end
if ~isdeployed
    siftroot = fileparts(which('StartSIFT'));
else
    try
        siftroot = utl_whichfile('StartSIFT');
    catch
        siftroot = [ctfroot filesep 'dependency' filesep 'SIFT-private'];
        if ~exist(siftroot,'dir'),siftroot = pwd;end
    end
end
disp(siftroot)
ArfitURL    = 'http://www.gps.caltech.edu/~tapio/arfit/arfit.zip';
ARfitTargetPath  = [siftroot SLASH 'external' SLASH 'arfit'];

% add subfolders to path (if not already present)
if ~exist('vis_TimeFreqGrid.m','file') && ~isdeployed
    disp('SIFT: Adding subfolders.')
    addpath(genpath(siftroot));
end

% call up the splash screen
if interactive && ~nosplash
    gui_splashscreen;
end

% optionally download arfit (if not already present)
if ~exist(ARfitTargetPath,'file') && interactive
    res = input('SIFT: Would you like to download and install the ARFIT toolbox as a SIFT plugin (recommended)? ''y''/''n'': ','s');
    if strcmpi(res,'y')

        fprintf('SIFT: Downloading and installing ARFIT from %s ...\n',ArfitURL);

        try
            outdir = unzip(ArfitURL,ARfitTargetPath);
        catch e
            switch e.identifier
                case 'MATLAB:unzip:urlwriteError'
                    fprintf('%s Aborting download.\n',e.message);
                case 'MATLAB:unzip:invalidZipFile'
                    fprintf('Unable to unpack file %s Aborting installation.\n', e.message);
                    res=rmdir(ARfitTargetPath,'s'); % clean up
            end

            return;
        end
        if ~isdeployed
            addpath(genpath(ARfitTargetPath));
            fprintf('SIFT: ARFIT installed to %s%c\n',fileparts(outdir{1}),SLASH);
        end
    end
end

ok = true;

fprintf('Start SIFTing!\n');

repo: REST


text:
function Animate_System(States,Params,Trajectory,Goal_Idx)

Fig_Num = 2;

N_Frames = Params.Sim.Frames_Per_Second * diff(Params.Sim.Tspan) / Params.Sim.Speed_Up_Video;
Frames = linspace(1,Params.Sim.Nsteps,N_Frames);
Frames = round(Frames);

Time_Target = linspace(1,diff(Params.Sim.Tspan)/Params.Sim.Speed_Up_Video,N_Frames);
TimeStep=0;

Npts = Params.Traj.Npts;

x = States(1,:);
y = States(2,:);
th = States(3,:);
phi = States(4,:);

Lt = Params.Dyn.Lt;
Lc = Params.Dyn.Lc;

%[dStates, A, B, Positions] = Derive_EoM()

Pa = [...   %Position of the back of the trailer
    x;
    y];

Pb = [...    %Position of the back of the cab
 x - Lt*sin(th);
 y + Lt*cos(th)];

Pc = [...   %Position of the front of the cab
 x - Lc*(cos(phi).*sin(th) + cos(th).*sin(phi)) - Lt*sin(th);
 y + Lc*(cos(phi).*cos(th) - sin(phi).*sin(th)) + Lt*cos(th)];

if strcmp(Params.Sim.Direction,'reverse')
    Dir = 2;
else
    Dir = 1;
end
Road = zeros(2,Npts);
for i=1:Npts
    Road(:,i) = Trajectory.States{Dir,i}(1:2);
end

All_Points = [Pa,Pb,Pc,Road];

Xmin = min(All_Points(1,:));
Xmax = max(All_Points(1,:));
Ymin = min(All_Points(2,:));
Ymax = max(All_Points(2,:));

figure(Fig_Num)
tic;
for i=Frames
    clf(Fig_Num);
    hold on
    axis([Xmin,Xmax,Ymin,Ymax])
    axis equal
    Trailer_X = [Pa(1,i); Pb(1,i)];
    Trailer_Y = [Pa(2,i); Pb(2,i)];
    Cab_X = [Pb(1,i); Pc(1,i)];
    Cab_Y = [Pb(2,i); Pc(2,i)];

    Idx = Goal_Idx(i);
    if Idx < 1
        Idx = Params.Traj.Npts;
        %At the end of the road
    end

    Target = Road(:,Idx);

    %Plot road
        plot(Road(1,:),Road(2,:),'k','LineWidth',Params.Sim.RoadWidth)
    %Plot the traces:
        plot(Pa(1,1:i),Pa(2,1:i),'b','LineWidth',4)
        plot(Pc(1,1:i),Pc(2,1:i),'g','LineWidth',4)
    %Plot truck
        plot(Trailer_X,Trailer_Y,'b','LineWidth',20)
        plot(Cab_X, Cab_Y,'g','LineWidth',17)
    %Plot the control target
        plot(Target(1),Target(2),'ro','MarkerSize',18,'LineWidth',3)
    title('Simulation of the Tractor Trailer Truck')
    axis equal

    TimeStep = TimeStep + 1;
    Time_Error = Time_Target(TimeStep)-toc;
    if Time_Error > 0
        pause(Time_Error)
    else
        pause(0.001)
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function jsubs = cfg2jobsubs(item, subs)

% function jsubs = cfg2jobsubs(item, subs)
% Return the subscript into the job tree for a given subscript vector into
% the val part of the cfg tree.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg2jobsubs.m 1862 2008-06-30 14:12:49Z volkmar $

rev = '$Rev: 1862 $'; %#ok

% Only de-reference subscripts into .val{x} field
if isequal(subs(1), substruct('.','val')) && strcmp(subs(2).type, '{}')
    % return possible indices into .val{1} field
    citem = subsref(item, subs(1:2));
    if numel(subs) > 2
        jsubs1 = cfg2jobsubs(citem, subs(3:end));
    else
        jsubs1 = [];
    end;
    if numel(item.values)==1 && isa(item.values{1},'cfg_branch') && ~item.forcestruct,
        cjsubs = substruct('.', gettag(citem), '()',subs(2).subs);
    else
        if numel(item.values)==1 && ~item.forcestruct,
            cjsubs = substruct('{}', subs(2).subs);
        else
            cjsubs = substruct('{}', subs(2).subs, '.', gettag(citem));
        end;
    end;
    jsubs = [cjsubs jsubs1];
else
    cfg_message('matlabbatch:cfg2jobsubs:wrongsubs', 'Inappropriate subscript reference in item ''%s''.', item.tag);
    jsubs = struct('type',{},'subs',{});
end;

author: spm


text:
classdef prtPreProcElm < prtPreProc
    %prtPreProcElm - Extreme Learning Machine
    %   A 1 Hidden Layer Neural Network with Random Weights
    %
    % Example:
    %    dsTrain = prtDataGenXor;
    %    dsTest = prtDataGenXor;
    %    algo = prtPreProcElm('nNeurons',1000) + prtClassLr;
    %    algo = algo.train(dsTrain);
    %    yOut = algo.run(dsTest);
    %    close all;
    %    prtScoreRoc(yOut);

    properties (SetAccess=private)
        name = 'Extreme Learning Machine'
        nameAbbreviation = 'ELM'
    end

    properties (SetAccess = protected)

    end

    properties
        nNeurons = 100;
        activationFunction = @(x)1./(1 + exp(-x));

        weights = [];
        bias = []
    end

    methods

               % Allow for string, value pairs
        function self = prtPreProcElm(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,dataSet)
            self.weights = rand(dataSet.nFeatures, self.nNeurons)*2 - 1;
            self.bias = rand(1,self.nNeurons);
        end

        function dataSet = runAction(self,dataSet)
           dataSet.X = self.activationFunction(bsxfun(@plus, dataSet.X*self.weights, self.bias));
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% DEMBARENCORANK Do ranking experiments on data from Barenco et al in Genome Biology.

% SHEFFIELDML

colordef white
[y, yvar, gene, times, scale, rawExp, rawVar] = gpsimLoadBarencoData;

load demBarenco1.mat

% Get the default options structure.
options = gpsimOptions;
options.includeNoise = 1;
% Fix one decay (from the fourth gene --- p21) to 0.8 hr^-1, and
% the corresponding sensitivity (see just after eqn 2 in the
% mathematical methods of Barenco et al.)
options.fix(1).index = 8;
options.fix(1).value = expTransform(0.8, 'xtoa');
options.fix(2).index = 9;
options.fix(2).value = expTransform(1, 'xtoa');

% initialise the model.
model.type = 'cgpsim';  % This new model type is a hack to run
                        % the model in a hierarchical manner.
                        % need to do this more elegantly later.
for i =1:3
  model.comp{i} = gpsimCreate(5, 1, times, y{i}, yvar{i}, options);
end

% Learn the model.
model = modelOptimise(model, [], [], 1, 3000);

% Each component of the model gives us a prediction given one set
% of replicates. The replicates are from an independent cell line,
% so the probabilistic assumption here is that the samples of f(t)
% are independent for each prediction (but inverse width, decay,
% basal rate and sensitivity parameters are all shared).
for j = 1:length(model.comp)

  % Generate predictions of the functions.
  % to do this we need to compute the K_xf portions of the kernel
  % (simXrbfKernCompute does this for us).
  predt = [linspace(-2, 14, 100) 0:2:12]';
  proteinKern = kernCreate(model.comp{1}.t, 'rbf');
  proteinKern.inverseWidth = ...
      model.comp{j}.kern.comp{1}.inverseWidth;
  K = [];
  for i=1:model.comp{j}.kern.numBlocks
    K = [K; simXrbfKernCompute(model.comp{j}.kern.comp{i}, proteinKern, ...
                               model.comp{j}.t, predt)];

  end
  predF = K'*model.comp{j}.invK*model.comp{j}.y;
  varF = kernDiagCompute(proteinKern, predt) - sum(K.*(model.comp{j}.invK*K), 1)';

  % Take out predictions at data points.
  % Use them to get the scale for the other data.
  dataF = predF(end-6:end);
  dataVarF = varF(end-6:end);
  predF(end-6:end) = [];
  varF(end-6:end) = [];
  predt(end-6:end) = [];
  scalePred = sqrt(var(dataF));

  % Info from Martino Paper:
  % 'True f' from Figure 3.
  % Got these figures with a ruler ...
  % Don't actually plot these below, but they are stored for reference
  truef = [0 1.6 2.6 2.5 2.6 1.6 0.9];
  truef = truef/sqrt(var(truef))*scalePred;

  % Figure 2(a) histograms;
  B = [2.6 1.5 0.5 0.2 1.35]; % From Martino paper ... but don't know the scale
  B = B/mean(B)*mean(model.comp{1}.B); % do a rough rescaling so
                                       % that the scales match.
  S = [3 0.8 0.7 1.8 0.7]/1.8; % From Martino paper ... but here we
                               % know the scale, because p21 is
                               % fixed to 1.
  D = [1.2 1.6 1.75 3.2 2.3]*0.8/3.2; % From Martino paper, again
                                      % we know the scale because
                                      % p21 is fixed to 0.8.

  % Martino f from Figure 2(b), again measured with a ruler.
  barencof = [0 2.7 3.9 2.3 1.5 1.6 1.4]/(1.8*mean(S))*mean(model.comp{1}.S);
  barencof = barencof/sqrt(var(barencof))*scalePred;

  figure, lin = plot(predt, predF, '-');
  hold on,
  bh = plot(predt, predF + 2*sqrt(varF), '--');
  bh = [bh plot(predt, predF - 2*sqrt(varF), '--')];
  lin = [lin plot(0:2:12, barencof, 'rx')];
  set(bh, 'lineWidth', 3);
  set(lin, 'lineWidth', 4);
  set(lin, 'markersize', 20);
  set(gca, 'fontname', 'arial', 'fontsize', 24, 'xlim', xlim)
  fileName = ['demBarenco1_profile' num2str(j)];
  print('-deps', ['../tex/diagrams/' fileName]);
  pos = get(gcf, 'paperposition')
  origpos = pos;
  pos(3) = pos(3)/2;
  pos(4) = pos(4)/2;
  set(gcf, 'paperposition', pos);
  lineWidth = get(gca, 'lineWidth');
  set(gca, 'lineWidth', lineWidth*2);
  %print('-dpng', ['../html/' fileName])
  set(gca, 'lineWidth', lineWidth);
  set(gcf, 'paperposition', origpos)

end

order = [1 5 3 4 2];

counter = 0;

% Plot first basal transcription rates.
figure
bar([model.comp{1}.B(order); B]', 0.6); colormap([0 0 0; 1 1 1]);
set(gca, 'xticklabel', {'DDB2', 'hPA26', 'TNFRSF20b', 'p21', 'BIK'})
fileName = ['demBarenco1_basal'];
print('-deps', ['../tex/diagrams/' fileName]);
pos = get(gcf, 'paperposition')
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
print('-dpng', ['../html/' fileName])
set(gcf, 'paperposition', origpos)
set(gca, 'lineWidth', lineWidth);

% Plot the sensitivities.
figure
bar([model.comp{1}.S(order); S]', 0.6); colormap([0 0 0; 1 1 1]);
set(gca, 'xticklabel', {'DDB2', 'hPA26', 'TNFRSF20b', 'p21', ...
                    'BIK'})
fileName = ['demBarenco1_sensitivity'];
print('-deps', ['../tex/diagrams/' fileName]);
pos = get(gcf, 'paperposition')
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
print('-dpng', ['../html/' fileName])
set(gcf, 'paperposition', origpos)
set(gca, 'lineWidth', lineWidth);

% Finally plot degradation rates.
figure
bar([model.comp{1}.D(order); D]', 0.6); colormap([0 0 0; 1 1 1]);
set(gca, 'xticklabel', {'DDB2', 'hPA26', 'TNFRSF20b', 'p21', ...
                    'BIK'})
fileName = ['demBarenco1_decay'];
print('-deps', ['../tex/diagrams/' fileName]);
pos = get(gcf, 'paperposition')
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
print('-dpng', ['../html/' fileName])
set(gcf, 'paperposition', origpos)
set(gca, 'lineWidth', lineWidth);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function check = r8sp_check ( m, n, nz_num, row, col )

%*****************************************************************************80
%
%% R8SP_CHECK checks that a R8SP matrix data structure is properly sorted.
%
%  Discussion:
%
%    This routine assumes that the data structure has been sorted,
%    so that the entries of ROW are ascending sorted, and that the
%    entries of COL are ascending sorted, within the group of entries
%    that have a common value of ROW.
%
%    The R8SP storage format stores the row, column and value of each nonzero
%    entry of a sparse matrix.
%
%    The R8SP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
%    ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of
%    the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in
%    the matrix.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and
%    column indices of the nonzero elements.
%
%    Output, logical CHECK, is TRUE if the matrix is properly defined.
%
  check = 1;
%
%  Check 1 <= ROW(*) <= M.
%
  for k = 1 : nz_num

    if ( row(k) < 1 | m < row(k) )
      check = 0;
      return
    end

  end
%
%  Check 1 <= COL(*) <= N.
%
  for k = 1 : nz_num

    if ( col(k) < 1 | n < col(k) )
      check = 0;
      return
    end

  end
%
%  Check that ROW(K) <= ROW(K+1).
%
  for k = 1 : nz_num - 1

    if ( row(k+1) < row(k) )
      check = 0;
      return
    end

  end
%
%  Check that, if ROW(K) == ROW(K+1), that COL(K) < COL(K+1).
%
  for k = 1 : nz_num - 1

    if ( row(k) == row(k+1) )
      if ( col(k+1) <= col(k) )
        check = 0;
        return
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [A,C] = update_spatial_components_nb(Y,C,A_,P,options)

% update spatial footprints and background through Basis Pursuit Denoising
% for each pixel i solve the problem
%   [A(i,:),b(i)] = argmin sum(A(i,:))
%       subject to || Y(i,:) - A(i,:)*C + b(i)*f || <= sn(i)*sqrt(T);
% for each pixel the search is limited to a few spatial components

% INPUTS:
% Y:    raw data
% C:    temporal components
% f:    temporal background
% A_:   current estimate of spatial footprints (used for determining search locations only)
% P:    dataset parameters (used for noise values and interpolated entries)

% options    parameter struct (for noise values and other parameters)

% OUTPUTS:
% A:    new estimate of spatial footprints
% C:    temporal components (updated only when spatial components are completely removed)

% Written by:
% Pengcheng Zhou, Columbia University, 2015,
% with modifications from update_spatial_components.m
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

warning('off', 'MATLAB:maxNumCompThreads:Deprecated');
memmaped = isobject(Y);
if memmaped
    sizY = size(Y,'Y');
    d = prod(sizY(1:end-1));
    T = sizY(end);
else
    [d,T] = size(Y);
end
if nargin < 5 || isempty(options); options = []; end
if ~isfield(options,'d1') || isempty(options.d1); d1 = input('What is the total number of rows? \n'); options.d1 = d1; else d1 = options.d1; end          % # of rows
if ~isfield(options,'d2') || isempty(options.d2); d2 = input('What is the total number of columns? \n'); options.d2 = d2; else d2 = options.d2; end          % # of columns
if ~isfield(options,'show_sum'); show_sum = 0; else show_sum = options.show_sum; end            % do some plotting while calculating footprints
if ~isfield(options,'interp'); Y_interp = sparse(d,T); else Y_interp = options.interp; end      % identify missing data
if ~isfield(options,'use_parallel'); use_parallel = ~isempty(which('parpool')); else use_parallel = options.use_parallel; end % use parallel toolbox if present
if ~isfield(options,'search_method'); method = []; else method = options.search_method; end     % search method for determining footprint of spatial components

if nargin < 2 || (isempty(A_) && isempty(C))  % at least either spatial or temporal components should be provided
    error('Not enough input arguments')
else
    if ~isempty(C); K = size(C,1); elseif islogical(A_); K = size(A_,2); else K = size(A_2,2) - options.nb; end
end

if nargin < 5 || isempty(P); P = preprocess_data(Y,1); end  % etsimate noise values if not present
if nargin < 4 || isempty(A_);
    IND = ones(d,size(C,1));
else
    if islogical(A_)     % check if search locations have been provided, otherwise estimate them
        IND = A_;
        if isempty(C)
            INDav = double(IND)/diag(sum(double(IND)));
            C = max(INDav'*Y,0);
        end
    else
        IND = determine_search_location(A_(:,1:K),method,options);
    end
end

options.sn = P.sn;
if ~memmaped
    Y(P.mis_entries) = NaN; % remove interpolated values
end

Cf = C;

if use_parallel         % solve BPDN problem for each pixel
    Nthr = max(2*maxNumCompThreads,round(d*T/2^24));
    siz_row = [floor(d/Nthr)*ones(Nthr-1,1);d-floor(d/Nthr)*(Nthr-1)];
    indeces = [0;cumsum(siz_row)];
    if ~memmaped
        Ycell = mat2cell(Y,siz_row,T);
    else
        Ycell = cell(Nthr,1);
    end
    INDc =  mat2cell(IND,siz_row,K);
    Acell = cell(Nthr,1);
    Psnc = mat2cell(options.sn(:),siz_row,1);
    for nthr = 1:Nthr
        Acell{nthr} = zeros(siz_row(nthr),size(Cf,1));
        if memmaped
            Ytemp = double(Y.Yr(indeces(nthr)+1:indeces(nthr+1),:));
        else
            Ytemp = Ycell{nthr};
        end
        for px = 1:siz_row(nthr)
            fn = ~isnan(Ytemp(px,:));       % identify missing data
            ind = find(INDc{nthr}(px,:));
            if ~isempty(ind);
                %[~, ~, a, ~] = lars_regression_noise(Ycell{nthr}(px,fn)', Cf(ind2,fn)', 1, Psnc{nthr}(px)^2*T);
                [~, ~, a, ~] = lars_regression_noise(Ytemp(px,fn)', Cf(ind,fn)', 1, Psnc{nthr}(px)^2*T);
                Acell{nthr}(px,ind) = a';
            end
        end
    end
    A = cell2mat(Acell);
else
    A = zeros(d,K);
    sA = zeros(d1,d2);
    for px = 1:d   % estimate spatial components
        fn = ~isnan(Y(px,:));       % identify missing data
        ind = find(IND(px,:));
        if ~isempty(ind);
            [~, ~, a, ~] = lars_regression_noise(Y(px,fn)', Cf(ind,fn)', 1, options.sn(px)^2*T);
            a(isnan(a)) = 0;
            A(px,ind) = a';
            sA(px) = sum(a);
        end
        if show_sum
            if mod(px,d1) == 0;
               figure(20); imagesc(sA); axis square;
               title(sprintf('Sum of spatial components (%i out of %i columns done)',round(px/d1),d2)); drawnow;
            end
        end
    end
end

A(isnan(A))=0;
A = sparse(A);
A = threshold_components(A,options);  % post-processing of components

fprintf('Updated spatial components \n');

ff = find(sum(A(:,1:K))==0);           % remove empty components
if ~isempty(ff)
    K = K - length(ff);
    A(:,ff) = [];
    C(ff,:) = [];
end

author: zhoupc


text:
function line = edgeToLine(edge)
%EDGETOLINE Convert an edge to a straight line
%
%   LINE = edgeToLine(EDGE);
%   Returns the line containing the edge EDGE.
%
%   Example
%       edge = [2 3 4 5];
%       line = edgeToLine(edge);
%       figure(1); hold on; axis([0 10 0 10]);
%       drawLine(line, 'color', 'g')
%       drawEdge(edge, 'linewidth', 2)
%
%   See also
%   edges2d, lines2d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-07-23,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

line = [edge(:, 1:2) edge(:, 3:4)-edge(:, 1:2)];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [a,d]=dualdiag(w,b)
%DUALDIAG Simultaneous diagonalisation of two hermitian matrices [A,D]=(W,B)
% Given two hermitian matrices W and B with W positive definite, this routine
% calculates A such that A'*W*A=I and A'*B*A=diag(D). The D will be in descending order.
%
% Suppose we have several N-dimensional data row-vectors arising from each of C different classes of data.
% for each class, c, we can form the mean data vector m(c) and the within-class covariance matrix W(c)
% We can then form the between class covariance matrix B by taking the covariance of the mean vectors m(1), m(2), ...
% and also the averaged within-class covariance matrix W by averaging W(1), W(2), ...
% If we then take A=dualdiag(W,B) and postmultiply all our original data vectors by A, we obtain new
% data vectors for which the average within-class covariance matrix is the identity and for which
% the first few components contain most of the information that is useful in discriminating between classes.

% An alternative algorithm that is 20% faster but slightly less accurate is:
% n=size(w,1);
% [v,l]=eig(w\b);
% [s,i]=sort(-diag(l));
% s=-s;
% d=l(i*(n+1)-n);
% q=sqrt(diag(v'*w*v))'.^(-1);
% a=v(:,i).*q(ones(n,1),i);

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

[y,l]=eig(w+w');
z=y*diag(sqrt(diag(l*0.5)).^(-1));
[u,s,v]=svd(z'*(b+b')*z);
d=diag(s)*0.5;
a=z*u;

author: decouples


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (C) 1995-2020 The Octave Project Developers
%
% See the file COPYRIGHT.md in the top-level directory of this
% distribution or <https://octave.org/copyright/>.
%
% This file is part of Octave.
%
% Octave 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.
%
% Octave 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 Octave; see the file COPYING.  If not, see
% <https://www.gnu.org/licenses/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% -*- texinfo -*-
% @deftypefn  {} {} cov (@var{x})
% @deftypefnx {} {} cov (@var{x}, @var{opt})
% @deftypefnx {} {} cov (@var{x}, @var{y})
% @deftypefnx {} {} cov (@var{x}, @var{y}, @var{opt})
% Compute the covariance matrix.
%
% If each row of @var{x} and @var{y} is an observation, and each column is
% a variable, then the @w{(@var{i}, @var{j})-th} entry of
% @code{cov (@var{x}, @var{y})} is the covariance between the @var{i}-th
% variable in @var{x} and the @var{j}-th variable in @var{y}.
% @tex
% $$
% \sigma_{ij} = {1 \over N-1} \sum_{i=1}^N (x_i - \bar{x})(y_i - \bar{y})
% $$
% where $\bar{x}$ and $\bar{y}$ are the mean values of @var{x} and @var{y}.
% @end tex
% @ifnottex
%
% @example
% cov (@var{x}) = 1/(N-1) * SUM_i (@var{x}(i) - mean(@var{x})) * (@var{y}(i) - mean(@var{y}))
% @end example
%
% @noindent
% where @math{N} is the length of the @var{x} and @var{y} vectors.
%
% @end ifnottex
%
% If called with one argument, compute @code{cov (@var{x}, @var{x})}, the
% covariance between the columns of @var{x}.
%
% The argument @var{opt} determines the type of normalization to use.
% Valid values are
%
% @table @asis
% @item 0:
%   normalize with @math{N-1}, provides the best unbiased estimator of the
% covariance [default]
%
% @item 1:
%   normalize with @math{N}, this provides the second moment around the mean
% @end table
%
% Compatibility Note:: Octave always treats rows of @var{x} and @var{y}
% as multivariate random variables.
% For two inputs, however, @sc{matlab} treats @var{x} and @var{y} as two
% univariate distributions regardless of their shapes, and will calculate
% @code{cov ([@var{x}(:), @var{y}(:)])} whenever the number of elements in
% @var{x} and @var{y} are equal.  This will result in a 2x2 matrix.
% Code relying on @sc{matlab}'s definition will need to be changed when
% running in Octave.
% @seealso{corr}
% @end deftypefn

function c = cov (x, y, opt)

if (nargin < 2)
    y = [];
end
if (nargin < 3)
    opt = 0;
end
if (nargin < 1 || nargin > 3)
    error('Usage: c = cov (x, y = [], opt = 0)');
end

if (~(isnumeric (x) || islogical (x)) || ~(isnumeric (y) || islogical (y)))
    error ('cov: X and Y must be numeric matrices or vectors');
end

if (ndims (x) ~= 2 || ndims (y) ~= 2)
    error ('cov: X and Y must be 2-D matrices or vectors');
end

if (nargin == 2 && isscalar (y))
    opt = y;
end

if (opt ~= 0 && opt ~= 1)
    error ('cov: normalization OPT must be 0 or 1');
end

% Special case, scalar has zero covariance
if (isscalar (x))
    if (isa (x, 'single'))
        c = single (0);
    else
        c = 0;
    end
    return;
end

if (isrow (x))
    x = x.';
end
n = size(x,1);

if (nargin == 1 || isscalar (y))
    %x = center (x, 1);
    x = bsxfun(@minus, x, mean (x, 1));
    c = x' * x / (n - 1 + opt);
else
    if (isrow (y))
        y = y.';
    end
    if (size(y,1) ~= n)
        error ('cov: X and Y must have the same number of observations');
    end
    %x = center (x, 1);
    %y = center (y, 1);
    x = bsxfun(@minus, x, mean (x, 1));
    y = bsxfun(@minus, y, mean (y, 1));
    c = x' * y / (n - 1 + opt);
end

end

author: brainstorm-tools


text:
function r8bto_print ( m, l, a, title )

%*****************************************************************************80
%
%% R8BTO_PRINT prints a R8BTO matrix.
%
%  Discussion:
%
%    The R8BTO storage format is for a block Toeplitz matrix. The matrix
%    can be regarded as an L by L array of blocks, each of size M by M.
%    The full matrix has order N = M * L.  The L by L matrix is Toeplitz,
%    that is, along its diagonal, the blocks repeat.
%
%    Storage for the matrix consists of the L blocks of the first row,
%    followed by the L-1 blocks of the first column (skipping the first row).
%    These items are stored in the natural way in an (M,M,2*L-1) array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the order of the blocks of the matrix A.
%
%    Input, integer L, the number of blocks in a row or column of A.
%
%    Input, real A(M,M,2*L-1), the R8BTO matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8bto_print_some ( m, l, a, 1, 1, m*l, m*l, title );

  return
end


text:
function test_ft_channelnormalise

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_channelnormalise

fs = 500;
nchan = 32;
start_time = -1; % seconds
end_time = 2.5; % seconds
nsamples = (end_time - start_time) * fs + 1;

data = [];
data.time{1} = linspace(start_time, end_time, nsamples);
data.trial{1} = randn(nchan,nsamples);
data.label = cellstr(num2str((1:nchan).'));

cfg = [];
dataout = ft_channelnormalise(cfg, data);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [gSig, gSiz, ring_radius, min_corr, min_pnr] = set_parameters(obj)
%% choose parameters used in CNMF-E

%% Author: Pengcheng Zhou, Columbia University, 2017
%zhoupc1988@gmail.com

gSiz = obj.options.gSiz;
gSig = obj.options.gSig;
ring_radius = obj.options.ring_radius;
min_corr = obj.options.min_corr;
min_pnr = obj.options.min_pnr;

%% choose gSiz, which is the neuron diameter
fprintf('\n*** Determine the neuron diameter (gSiz) and the gaussian width (gSig) of the filtering kernel.\n');
fprintf('\t-------------------------- GUIDE --------------------------\n');
fprintf('\tgSiz is usually slightly larger than a neuron;\n');
fprintf('\tgSig is usually selected as 1/4 of gSiz when the data is 1p\n');
fprintf('\tgSig is usually selected as 1 when the data is 2p.\n');
fprintf('\tIntegers are preferred for both values, but not necessary.\n');
fprintf('\tOdd number is preferred for gSiz. \n');
fprintf('\tWhen you want to make a change, an image will be shown.\n');
fprintf('\tYou can zoon in to see the size of a typical neuron.\n');
fprintf('\t--------------------------  END  --------------------------\n\n');

fprintf('You current selection of parameters (gSiz, gSig) are (%.1f, %.1f).\n', gSiz, gSig);
temp = input('Do you want to make a change? (y/n)    ', 's');
if strcmpi(temp, 'n')
    fprintf('Your values for (gSiz, gSig) will stay the same\n');
else
    Y = get_patch_data(obj.P.mat_data, [], [1,100]);
    [d1,d2, T] = size(Y);
    Y = reshape(Y, d1*d2, T);
    fprintf('Wait for a few second ....\n');
    if obj.options.center_psf
    [u, v] = nnmf(double(Y),1);
        img = reshape(max(double(Y)-u*v, [], 2), d1,d2);
    else
            img = reshape(max(Y,[],2), d1,d2);
    end
    fig_1 = figure;
    imagesc(img); axis equal;
    while true
        temp = input('type new values for (gSiz, gSig):  ', 's');
        try
            new_values = str2num(temp); %#ok<*ST2NM>
            gSiz = new_values(1);
            gSig = new_values(2);
            obj.options.gSig = gSig;
            obj.options.gSiz = gSiz;
            fprintf('Good! You current selection of parameters (gSiz, gSig) are (%.1f, %.1f).\n', gSiz, gSig);
            break;
        catch
            warning('values are bad. try again');
            continue;
        end
    end
    try
        close(fig_1);
    end
end

%%
%% choose the radius of the ring
if strcmpi('obj.options.background_model', 'ring')
    fprintf('\n*** Determine the ring radius for estimating background fluctuations.\n');
    fprintf('\t-------------------------- GUIDE --------------------------\n');
    fprintf('\tThe ring radius is selected to be larger than neuron diameters. \n');
    fprintf('\tThus pixels in the center only shares the common background fluctuations with pixels on the ring.\n');
    fprintf('\tA usual ratio between ring radius and neuron diameter is between 1 to 3. 1.5 is the default value.\n');
    fprintf('\tSmall values have a change of being smaller than neuron size.\n');
    fprintf('\tLarge values increase the computation cost.\n');
    fprintf('\t--------------------------  END  --------------------------\n\n');

    fprintf('You current selection of ring_radius is %d pixels.\n', ring_radius);
    temp = input('Do you want to make a change? (y/n)    ', 's');
    if strcmpi(temp, 'n')
        fprintf('Your values for ring_radius will stay the same\n');
    else
        while true
            temp = input('type new values for ring_radius:  ', 's');
            try
                ring_radius = str2double(temp);
                obj.options.ring_radius = ring_radius;
                fprintf('Good! You current selection of parameters ring_radius is %2d.\n', ring_radius);
                break;
            catch
                warning('values are bad. try again');
                continue;
            end
        end
    end
end
%% choose parameters for min_corr and min_pnr
fprintf('\n*** Determine min_corr and min_pnr for screening seed pixels.\n');
fprintf('\t-------------------------- GUIDE --------------------------\n');
fprintf('\tA seed pixel is used for initializing one neuron.\n');
fprintf('\tSeed pixel usually locates at the center of each neuron.\n');
fprintf('\tSeed pixel is determined using two values: min_corr and min_pnr.\n');
fprintf('\tmin_pnr stands for minimum local correlation. \n');
fprintf('\tPNR stands for minimum peak-to-noise ratio.\n');
fprintf('\tWe want to pick thresholds that all neuron centers are above the thresholds,\n');
fprintf('\tand include the least number of non-ROI pixels. \n');
fprintf('\tThere is a trade-of between missing neurons and acquiring more false positives.\n');
fprintf('\tWhen it is hard to determine, we tend to choose missing neurons, i.e., larger threshlds.\n');
fprintf('\tThose missed neurons can be picked up from the residual video later.\n');
fprintf('\tTip: use Data Cursor to help you pick thresholds.\n');
fprintf('\t--------------------------  END  --------------------------\n\n');

fprintf('You current selection of parameters (min_corr, min_pnr) are (%.3f, %.3f).\n', min_corr, min_pnr);
temp = input('Do you want to make a change? (y/n)    ', 's');
if strcmpi(temp, 'n')
    fprintf('Your values for (min_corr, min_pnr) will stay the same\n');
else

    if isempty(obj.Cn) || isempty(obj.PNR)
        fprintf('computing the correlation image and the peak-to-noise ratio image....\n');
        [cn, pnr] = obj.correlation_pnr_parallel([1, 5000]);
        obj.Cn = cn;
        obj.PNR = pnr;
        fprintf('Done\n');
    else
        cn = obj.Cn;
        pnr = obj.PNR;
    end

    %find all local maximum as initialization point
    tmp_d = max(1,round(gSiz/4));
    v_max = ordfilt2(cn.*pnr, tmp_d^2, true(tmp_d));
    ind = (v_max==cn.*pnr);

    figure('papersize', [12, 3]);
    init_fig;
    subplot(131);
    imagesc(cn, [0, 1]);
    title('local corr. image');
    axis equal off tight;
    subplot(132);
    pnr_vmax = max(pnr(:))*0.8;
    imagesc(pnr, [3, pnr_vmax]);
    axis equal off tight;
    title('PNR image');
    subplot(133);
    imagesc(cn.*pnr, [3, pnr_vmax]);
    hold on;

    tmp_ind = ind & (cn>=min_corr) & (pnr>=min_pnr);
    [r, c] = find(tmp_ind);
    ax_seeds = plot(c, r, '.m', 'markersize', 10);
    axis equal off tight;
    title('candidate seed pixels');
    ylabel('PNR');
    xlabel('Cn');

    while true
        temp = input('type new values for (min_corr, min_pnr):  ', 's');
        try
            new_values = str2num(temp); %#ok<*ST2NM>
            min_corr = new_values(1);
            min_pnr = new_values(2);
            obj.options.min_corr = min_corr;
            obj.options.min_pnr = min_pnr;
            tmp_ind = ind & (cn>=min_corr) & (pnr>=min_pnr);
            [r, c] = find(tmp_ind);
            delete(ax_seeds);
            subplot(133);
            ax_seeds = plot(c, r, '.m', 'markersize', 10);

            temp = input('Is this good? (y/n)   ', 's');
            if strcmpi(temp, 'y')
                fprintf('Good! You current selection of parameters (min_corr, min_pnr) are (%2d, %2d).\n', min_corr, min_pnr);
                break;

            end
        catch
            warning('values are bad. try again');
            continue;
        end
    end
end

fprintf('\nThe parameters used in CNMF-E are \ngSig:\t\t%d\ngSiz:\t\t%d\nring_radius\t%d\nmin_corr:\t%.3f\nmin_pnr:\t%.3f\n', ...
    gSig, gSiz, ring_radius, min_corr, min_pnr);
end


text:
%% FUNCTION bsa_ihb
%    Singular Projection
%
%% OBJECTIVE
%   min 1/2*||x - a||_2^2
%    s.t. b'*x = r, 0<= x <= u,  b > 0
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Pinghua Gong and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% Related papers
%
% [1] KC. Kiwiel. On linear-time algorithms for the continuous
%     quadratic knapsack problem, Journal of Optimization Theory
%     and Applications, 2007
%

function [x_star,t_star,iter] = bsa_ihb(a,b,r,u)

% initilization
break_flag = 0;
t_l = a./b; t_u = (a - u)./b;
T = [t_l;t_u];
t_L = -inf; t_U = inf;
g_tL = 0; g_tU = 0;

iter = 0;
while ~isempty(T)
    iter = iter + 1;
    g_t = 0;
    t_hat = median(T);

    U = t_hat < t_u;
    M = (t_u <= t_hat) & (t_hat <= t_l);

    if sum(U)
       g_t = g_t + b(U)'*u(U);
    end
    if sum(M)
        g_t = g_t + sum(b(M).*(a(M) - t_hat*b(M)));
    end

    if g_t > r
        t_L = t_hat;
        T = T(T > t_hat);
        g_tL = g_t;
    elseif g_t < r
        t_U = t_hat;
        T = T(T < t_hat);
        g_tU = g_t;
    else
        t_star = t_hat;
        break_flag = 1;
        break;
    end
end
if ~break_flag
     t_star = t_L - (g_tL -r)*(t_U - t_L)/(g_tU - g_tL);
end
x_star = min(max(0,a - t_star*b),u);
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function xslt6p1(in,xsl,out);
%XSLT6p1 - XSLT implementation for MATLAB 6.1
%
% xslt6p1(xml_in,xsl_in,html_out)
%
% All three arguments are file names.
%
% The versions of xerces.jar and saxon.jar which ship with MATLAB 6.1
% must be on your classpath

% Copyright 2005-2010 The MathWorks, Inc.

in = i_url(in);
xsl = i_url(xsl);

tf = javax.xml.transform.TransformerFactory.newInstance;
streamstyle = javax.xml.transform.stream.StreamSource(xsl);
transformer = tf.newTransformer(streamstyle);

streamin = javax.xml.transform.stream.StreamSource(in);
streamout = javax.xml.transform.stream.StreamResult(out);

transformer.transform(streamin,streamout);

%%%%%%%%%%%%%%%%%%%
function url = i_url(filename)

filename = which(filename);
% This bit might be Windows-specific.
url = [ 'file:///' strrep(filename,'\','/') ];

author: Sable

repo: mcbench-benchmarks


text:
function [rU, WU, CU, HU] = uncovered_engine(X, T, ...
	initial_equity, distribution, risk_free_rate, stock_cost)
%uncovered_engine - Generate scenarios for an uncovered and a covered buy-write strategy.
%
%	[rU, WU, CU, HU] = uncovered_engine(X, T, ...
%	 	initial_equity, distribution, risk_free_rate, stock_cost);
%
% Inputs:
%	Stock Information
%		X - Stock total return prices including initial price [scalar].
%		T - Duration of investment period in years (terminal time - initial time) [scalar].
%	Fund Details
%		initial_equity - Initial (uncovered) total value of stock and cash held in asset [scalar].
%		distribution - Annualized fund distribution [scalar].
%	Other Details
%		risk_free_rate - Annualized risk-free rate [scalar].
%	Costs/Frictions
%		stock_cost - Proportional cost to buy or sell stock [scalar].
%
% Outputs:
%	rU - Total return for uncovered strategy [scalar].
%	WU - Sequence of total wealth for uncovered strategy [vector].
%	CU - Sequence of cash for uncovered strategy [vector].
%	HU - Sequence of holdings for uncovered strategy [vector].
%
% Comments:
%	1) This function generates scenarios for an uncovered portfolio strategy. To introduce some
%		degree of "realism" into the model, several inputs can be specified to control the
%		simulations. The next few comments provide additional details on these inputs.
%	2) X and T are stock total return prices and "times" in years that are assumed to be generated
%		by a geometric Brownian motion process with stochastic differential equation in the form
%			dX(t) = mu*X(t)*dt + volatility*X(t)*dB(t)
%		for t > 0.

% Copyright (C) 2012 The MathWorks, Inc.

% initialization

N = numel(X) - 1;				% N is the number of samples in X excluding the initial price
tau = T/N;						% tau is the time interval between samples in "years"

periods_per_day = floor(N/(252*T) + 0.5);	% periods_per_day is number of periods in a "day"

% initial position

current_price = X(1);									% initial price
current_time = 0;										% initial time

current_shares = floor(initial_equity/current_price);	% initial number of shares
current_cash = max(0, (initial_equity - current_shares*current_price));	% initial cash

initial_wealth = current_cash + current_shares*current_price;

% generate scenarios

% track shares, strike, cash, and wealth over time for testing

if nargout > 2
	HU = zeros(N+1,1);		% (H)oldings in stocks
	CU = zeros(N+1,1);		% (C)ash
	WU = zeros(N+1,1);		% (W)ealth

	HU(1) = current_shares;
	CU(1) = current_cash;
	WU(1) = current_cash + current_shares*current_price;
end

% loop over investment period for uncovered strategy

for iter = 2:N+1

	current_price = X(iter);
	current_time = tau*(iter - 1);

	% accrue interest on cash account at end of each day
	if mod(iter, periods_per_day) == 1		% note that period 1 happens outside loop
		current_cash = current_cash*(1 + risk_free_rate*tau*periods_per_day);
	end

	% buy more stock if enough cash available
	if current_cash > (current_price + stock_cost)
		adjusted_stock_price = current_price + stock_cost;
		shares_purchased = floor(current_cash/adjusted_stock_price);
		current_shares = current_shares + shares_purchased;
		current_cash = current_cash - shares_purchased*adjusted_stock_price;
	end

	% update test variables
	if nargout > 2
		HU(iter) = current_shares;
		CU(iter) = current_cash;
		WU(iter) = current_cash + current_shares*current_price;
	end
end

% at the end of the period, pay a distribution and sell shares if not enough cash

distribution = distribution*T;

needed_cash = distribution*(current_cash + current_shares*X(end)) - current_cash;

if needed_cash > 0
	adjusted_stock_price = X(end) + stock_cost;
	shares_sold = ceil(needed_cash/adjusted_stock_price);

	current_shares = current_shares - shares_sold;
	current_cash = current_cash + shares_sold*(X(end) - stock_cost);

	if nargout > 2
		HU(end) = current_shares;
		CU(end) = current_cash;
		WU(end) = current_cash + current_shares*X(end);
	end
end

% final scenario returns

rU = (current_cash + current_shares*X(end))/initial_wealth - 1;


text:
function e = line_cvt_energy ( n, a, b, x )

%*****************************************************************************80
%
%% LINE_CVT_ENERGY computes the CVT energy for a given set of generators.
%
%  Discussion:
%
%    Given a set of generators G over the line [A,B], then the energy
%    is defined as
%      E = integral ( a <= x <= b ) ( x - g(x) )^2 dx
%    where g(x) is the nearest generator to the point x.
%
%    For the 1D case, this integral can be evaluated exactly as the
%    sum of integrals over each subinterval:
%
%      E(i) = integral ( xl <= x <= xr ) ( x - x(i) )^2 dx
%           = ( ( x(i) - xl )^3 + ( xr - x(i) )^3 ) / 3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%
%    Input, real A, B, the left and right endpoints.
%
%    Input, real X(N), the generator locations.
%
%    Output, real E, the energy of the generator distribution.
%
  e = 0.0;

  for j = 1 : n

    if ( j == 1 )
      xl = a;
    else
      xl = ( x(j-1) + x(j) ) / 2.0;
    end

    if ( j == n )
      xr = b;
    else
      xr = ( x(j) + x(j+1) ) / 2.0;
    end

    e = e + ( ( x(j) - xl )^3 + ( xr - x(j) )^3  ) / 3.0;

  end

  return
end


text:
function string = format_interval(time)
% format_interval Format a time interval
%
% Formats a given time interval in seconds to a string containing
% days, hours and minutes.
%
% Input:
% - time (dobule): A numeric representation of time interval in seconds.
%
% Output:
% - string (string): A string representation of the interval
%

days = floor(time / (24*60*60));

time = time - days * (24*60*60);

hours = floor(time / (60*60));

time = time - hours * (60*60);

minutes = floor(time / (60));

%time = time - minutes * (60);

string = sprintf('%d days %d hours %d minutes', days, hours, minutes);

author: votchallenge

repo: toolkit-legacy


text:
function [x,tspan] = eulermaruyama(f,L,tspan,x0,Q)
%% EULERMARUYAMA - Numerical SDE solver: The Euler-Maruyama method
%
% Syntax:
%   [x,tspan] = eulermaruyama(f,L,tspan,x0,Q)
%
% In:
%   f      - Drift function, f(x,t)
%   L      - Diffusion function, L(x,t)
%   tspan  - Time steps to simulate, [t0,...,tend]
%   x0     - Initial condition
%   Q      - Spectral density (default: standard Brownian motion)
%
% Out:
%   x      - Solved values
%   tspan  - Time steps
%
% Description:
%   Integrates the system of stochatic differential equations
%     dx = f(x,t) dt + L(x,t) dbeta,  for x(0) = x0
%   over the time interval defined in tspan.
%
% Copyright:
%   2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%%

  % Check if Q given
  if nargin<5 || isempty(Q), Q = eye(size(L(x0,tspan(1)),2)); end

  % Cholesky factor of Q
  cQ = chol(Q,'lower');

  % Number of steps
  steps = numel(tspan);

  % Allocate space
  x = zeros(size(x0,1),steps);

  % Initial state
  x(:,1) = x0;

  % Iterate
  for k=2:steps

    % Time discretization
    dt = tspan(k)-tspan(k-1);

    % Increment
    db = sqrt(dt)*cQ*randn(size(Q,1),1);

    % Step
    x(:,k) = x(:,k-1) + ...
        f(x(:,k-1),tspan(k-1))*dt + ...
        L(x(:,k-1),tspan(k-1))*db;

  end

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
%CREATETREFOILKNOT Create a 3D mesh around a trefoild curve
%
%   [X, Y, Z] = createTrefoilKnot;
%   [X, Y, Z] = createTrefoilKnot(NPTS);
%
%   Example
%   createTrefoilKnot
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2015-01-07,    using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015 INRA - Cepia Software Platform.

%% Constants

% number of vertices of trefoil curve
nPoints = 200;

% thickness of the 3D mesh
thickness = .5;

% number of corners around each curve vertex
nCorners = 16;

%% Create trefoil curve

% parameterisation variable
t = linspace(0, 2*pi, nPoints + 1);
t(end) = [];

% trefoil curve coordinates
curve(:,1) = sin(t) + 2 * sin(2 * t);
curve(:,2) = cos(t) - 2 * cos(2 * t);
curve(:,3) = -sin(3 * t);

% display curve
figure;
drawPolyline3d(curve, 'LineWidth', 4, 'color', 'b');
axis equal; view(3);
axis([-4 4 -4 4 -2 2]);

%% Create surrounding mesh

% compute mesh
[v2, f2] = curveToMesh(curve, thickness, nCorners);

% display mesh
drawMesh(v2, f2, 'FaceAlpha', 0.5);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [S,U] = scat_joint_timefreq(X,cascade)
	Y = concatenate_freq(X);

	S = {};
	U = {};

	for m = 0:length(X)-1
		r = 1;

		S{m+1} = {};
		U{m+1} = {};

		for k = 1:length(Y{m+1}.signal)
			signal = Y{m+1}.signal{k};

			sz_orig = size(signal);
			if numel(sz_orig) == 2
				sz_orig(3) = 1;
			end

			ind = r:r+size(signal,1)-1;

			signal = reshape(signal,[sz_orig(1) sz_orig(2) sz_orig(3)]);

			if m > 0
				% note that here scat is 2d
				[S_fr,U_fr] = scat(signal,cascade);

				% needed for the case of U, are not init by scat
				if ~isfield(U_fr{1}.meta,'bandwidth')
					U_fr{1}.meta.bandwidth = 2*pi*ones(2,1);
				end
				if ~isfield(U_fr{1}.meta,'resolution')
					U_fr{1}.meta.resolution = 0*ones(2,1);
				end
				if ~isfield(U_fr{1}.meta,'j1')
					U_fr{1}.meta.j1 = -1*ones(1,1);
				end
				if ~isfield(U_fr{1}.meta,'j2')
					U_fr{1}.meta.j2 = -1*ones(1,1);
				end
			else
				S_fr = {};

				S_fr{1}.signal = {signal};
				S_fr{1}.meta.bandwidth = 2*pi*ones(2,1);
				S_fr{1}.meta.resolution = 0*ones(2,1);
				S_fr{1}.meta.j1 = -1*ones(1,1);
				S_fr{1}.meta.j2 = -1*ones(1,1);

				U_fr = {};

				U_fr{1}.signal = {signal};
				U_fr{1}.meta.bandwidth = 2*pi*ones(2,1);
				U_fr{1}.meta.resolution = 0*ones(2,1);
				U_fr{1}.meta.j1 = -1*ones(1,1);
				U_fr{1}.meta.j2 = -1*ones(1,1);
			end

			if isempty(S{m+1})
				for mp = 0:length(S_fr)-1
					S{m+1}{mp+1}.signal = {};
					U{m+1}{mp+1}.signal = {};
					rp(mp+1) = 1;
					rb(mp+1) = 1;
				end
			end

			for mp = 0:length(S_fr)-1
				% what do these do?
				%S_fr{mp+1} = unpad_layer_1d(S_fr{mp+1},[size(signal,1), size(signal,2)]);
				%U_fr{mp+1} = unpad_layer_1d(U_fr{mp+1},[size(signal,1), size(signal,2)]);

				for kp = 1:length(S_fr{mp+1}.signal)
					for t = 0:1
						if t == 0
							X_fr = S_fr;
							X = S;
							rc = rp;
						else
							X_fr = U_fr;
							X = U;
							rc = rb;
						end

						nsignal = X_fr{mp+1}.signal{kp};

						j1_count = size(nsignal,1);
						t_count = size(nsignal,2);

						ds = X_fr{mp+1}.meta.resolution(:,kp);

						j1_inds = ind(1:2^ds(1):end);

						nsignal = reshape(nsignal,[j1_count t_count sz_orig(3)]);

						for j1 = 1:j1_count
							X{m+1}{mp+1}.signal{rc(mp+1)} = ...
								reshape(nsignal(j1,:,:),t_count,sz_orig(3));
							X{m+1}{mp+1}.meta.bandwidth(1,rc(mp+1)) = ...
								X_fr{mp+1}.meta.bandwidth(2,kp)/2*pi* ...
								Y{m+1}.meta.bandwidth(j1_inds(j1));
							X{m+1}{mp+1}.meta.resolution(1,rc(mp+1)) = ...
								X_fr{mp+1}.meta.resolution(2,kp)+ ...
								Y{m+1}.meta.resolution(j1_inds(j1));
							X{m+1}{mp+1}.meta.j(:,rc(mp+1)) = ...
								Y{m+1}.meta.j(:,j1_inds(j1));
							X{m+1}{mp+1}.meta.fr_bandwidth(:,rc(mp+1)) = ...
								X_fr{mp+1}.meta.bandwidth(1,kp);
							X{m+1}{mp+1}.meta.fr_resolution(:,rc(mp+1)) = ...
								X_fr{mp+1}.meta.resolution(1,kp);
							X{m+1}{mp+1}.meta.tf_j1(:,rc(mp+1)) = ...
								X_fr{mp+1}.meta.j1(:,kp);
							X{m+1}{mp+1}.meta.tf_j2(:,rc(mp+1)) = ...
								X_fr{mp+1}.meta.j2(:,kp);
							rc(mp+1) = rc(mp+1)+1;
						end

						if t == 0
							S_fr = X_fr;
							S = X;
							rp = rc;
						else
							U_fr = X_fr;
							U = X;
							rb = rc;
						end
					end
				end
			end

			r = r+size(signal,1);
		end
	end

	for m = 0:length(S)-1
		temp = flatten_scat(S{m+1});
		temp = temp{1};
		temp.meta.tf_order = temp.meta.order;
		temp.meta = rmfield(temp.meta,'order');
		S{m+1} = temp;
		temp = flatten_scat(U{m+1});
		temp = temp{1};
		temp.meta.tf_order = temp.meta.order;
		temp.meta = rmfield(temp.meta,'order');
		U{m+1} = temp;
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function f = not(f)
%~   CLASSICFUN logical NOT.
%   NOT(F) returns a CLASSICFUN which evaluates to one at all points where F is
%   zero and one otherwise.  F cannot have any roots in its domain.  If F
%   does have roots, then LOGICAL(F) will return garbage with no warning. F may
%   be complex.

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

f.onefun = not(f.onefun);

end

author: chebfun

repo: chebfun


text:
classdef RawBatteryAndThermometerDataClass < BatteryAndThermometerDataBaseClass

    %% Public 'read-only' properties
    properties (SetAccess = private)
        FileNameAppendage = '_RawBattAndTherm.csv';
    end

    %% Public methods
    methods (Access = public)
        function obj = RawBatteryAndThermometerDataClass(varargin)
            fileNamePrefix = varargin{1};
            for i = 2:2:nargin
                if  strcmp(varargin{i}, 'SampleRate'), obj.SampleRate = varargin{i+1};
                else error('Invalid argument.');
                end
            end
            obj.Import(fileNamePrefix);

            % Set protected parent class variables
            obj.ThermometerUnits = 'lsb';
            obj.BatteryUnits = 'lsb';
        end
    end
end

author: xioTechnologies

repo: Gait-Tracking-With-x-IMU


text:
function tests = tapas_physio_findpeaks_test()
% Tests whether current findpeaks function of Matlab's signal processing
% toolbox delivers same results as previous version used in reference data
%
%    tests = tapas_physio_findpeaks_test()
%
% This relies on the GE/PPU3T dataset which is particularly challenging in
% terms of peak detection because of low SNR and drift of the PPU signal
%
% IN
%
% OUT
%
% EXAMPLE
%   tapas_physio_findpeaks_test
%
%   See also

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

tests = functiontests(localfunctions);
end

function test_ge_ppu3t_peaks(testCase)
%% Compares previously saved cpulse (detected cardiac pulses) from
% physio-structure to same output when re-running current version of
% GE PPU3T example
% both SPM or matlab-script based execution is possible (check parameter
% doUseSpm below!)

doUseSpm = true;

% run GE example and extract physio
pathPhysioPublic = fullfile(fileparts(mfilename('fullpath')), '..', '..', '..');
% TODO: Make generic!
pathExamples =  fullfile(pathPhysioPublic, '..', 'examples');

if doUseSpm
    pathCurrentExample = fullfile(pathExamples, 'GE/PPU3T');
    pathNow = pwd;
    cd(pathCurrentExample); % for prepending absolute paths correctly
    fileExample = fullfile(pathCurrentExample, 'ge_ppu3t_spm_job.mat');
    load(fileExample, 'matlabbatch');

    % remove unnecessary verbosity and processing of resp data
    matlabbatch{1}.spm.tools.physio.verbose.level = 0;
    matlabbatch{1}.spm.tools.physio.log_files.respiration = {''};

    spm_jobman('run', matlabbatch);
    cd(pathNow);

    % retrieve physio struct from saved file
    matlabbatch{1}.spm.tools.physio.model.output_physio = fullfile(pathCurrentExample, ...
        matlabbatch{1}.spm.tools.physio.save_dir{1}, ...
        matlabbatch{1}.spm.tools.physio.model.output_physio);
    load(matlabbatch{1}.spm.tools.physio.model.output_physio, 'physio');
    actPhysio = physio;
else % has verbosity...cannot switch it off
    fileExample = fullfile(pathExamples, 'GE/PPU3T/ge_ppu3t_matlab_script.m');
    run(fileExample); % will output a PhysIO=struct
    actPhysio = physio;
end

% load physio from reference data
fileReferenceData = fullfile(pathExamples, 'TestReferenceResults', 'utils', ...
    'physio_findpeaks_ge_ppu3t.mat');
load(fileReferenceData, 'physio');
expPhysio = physio;

% extract cpulse from actual and expected solution and compare
actSolution = actPhysio.ons_secs.cpulse;
expSolution = expPhysio.ons_secs.cpulse;

verifyEqual(testCase, actSolution, expSolution);

end

function test_ge_ppu3t_findpeaks_compatible(testCase)
%% Removes signal processing findpeaks from path and checks whether
% deprecated modified version from older Matlab release delivers same
% results
% Is also true, when no tapas_physio_findpeaks_compatible function exists

pathSignalToolbox = fileparts(which('findpeaks'));

hasCompatibleFindPeaks = ~isempty(which('tapas_physio_findpeaks_compatible'));

if ~isempty(pathSignalToolbox) && hasCompatibleFindPeaks
    rmpath(pathSignalToolbox);
    fprintf('\n\t!!! Temporarily removed %s from path to check performance vs legacy findpeaks\n', ...
        pathSignalToolbox);
    test_ge_ppu3t_peaks(testCase);
    addpath(pathSignalToolbox);
else
    assumeEqual(testCase, 1, 1);
end
end

author: translationalneuromodeling


text:
function test07()
% function test07()
%
% Computes a robust mean of angles
%

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    % Generate the data
    N = 250;
    theta_true = pi*(2*rand(1)-1);
    thetas = zeros(N,1);
    p_true = .25;
    for i = 1 : N
        if rand(1) < p_true
            thetas(i) = theta_true + .3*randn(1);
        else
            thetas(i) = pi*(2*rand(1)-1);
        end
    end
    X = [cos(thetas') ; sin(thetas)'];

    % Pick the manifold
    problem.M = spherefactory(2);

    % Parameters
    p = p_true;
    kappa = 3;
    c2 = besseli(0, 2*kappa);

    % Define the problem cost function
    problem.cost = @cost;
    function [val store] = cost(x, store)

        if ~isfield(store, 'fi')
            store.fi = (p/c2)*exp(2*kappa*X'*x);
        end
        fi = store.fi;

        val = -sum(log(fi + (1-p)));

    end

    % And its gradient
    problem.grad = @grad;
    function [g store] = grad(x, store)

        if ~isfield(store, 'fi')
            store.fi = (p/c2)*exp(2*kappa*X'*x);
        end
        fi = store.fi;

        g = -X*(2*kappa*(fi./(fi+(1-p))));
        g = g - (x'*g)*x;

    end

    % Check differentials consistency.
    % checkgradient(problem);

    % Solve with trust-regions and FD approximation of the Hessian
    warning('off', 'manopt:getHessian:approx');

    % Test many random initial guess
%     best_cost = inf;
%     best_x = [];
%     for i = 1 : 5
%         [x cst] = trustregions(problem);
%         if cst < best_cost
%             best_cost = cst;
%             best_x = x;
%         end
%     end
%     theta_found = angle(best_x(1)+1i*best_x(2));

    % Do a histogram of the data and select the initial guess as the peak
    [freqs, bins] = hist(thetas, 30);
    [~, ind] = max(freqs);
    x0 = [cos(bins(ind)) ; sin(bins(ind))];
    x = trustregions(problem, x0);
    theta_found = angle(x(1)+1i*x(2));

    fprintf('True theta: %g\nHist theta: %g\nFound theta: %g\n', theta_true, angle(x0(1)+1i*x0(2)), theta_found);

    subplot(2,1,1);
    hist(thetas, 30);
    xlim([-pi pi]);

    th = linspace(-pi, pi);
    fth = zeros(size(th));
    for i = 1 : length(th)
        fth(i) = cost([cos(th(i));sin(th(i))], struct());
    end
    subplot(2,1,2);
    plot(th, -fth);
    xlim([-pi pi]);

%     keyboard;

end

author: NicolasBoumal

repo: manopt