text:
%
% Created by Saurabh Tavildar on 3/23/16.
%
classdef PolarCode < handle
    % Polar code
    properties

        block_length;
        info_length;
        frozen_bits;
        n;
        design_epsilon;

        bit_reversed_order;
        info_bits;

        %sc_decoder
        u_sc;
        x_sc;

        %scl decoder
        list_size;
        p_scl;
        c_scl;
        i_scl;

        inactivePathIndices;
        inactivePathIndicesSize;
        activePathArray;
        pathIndexToArrayIndex;
        inactiveArrayIndices;
        inactiveArrayIndicesSize;
        arrayReferenceCount;

        lambda_offset;
        list_offset;

        llr_based_computation;
        llr_scl;
        llr_path_metric;

        % CRC
        crc_matrix;
        crc_size;

        %code construction
        cc_method;
        cc_parameter;
        cc_misc;

        %
        info_bit_order
    end

    methods

        %% code construction
        function obj = PolarCode(block_length, info_length, design_epsilon, crc_size)

            obj.block_length = block_length;
            obj.info_length = info_length;
            obj.n = log2(block_length);
            obj.design_epsilon = design_epsilon;

            if nargin < 4
                crc_size = 0;
            end

            obj.crc_size = crc_size;

            if obj.crc_size ~= 0
                obj.crc_matrix = floor(2*rand(obj.crc_size, info_length));
            end

            obj.bit_reversed_order = bitrevorder((1:obj.block_length)');

            channels = PolarCode.calculate_channel_polarization(design_epsilon, obj.n );
            channels = channels(obj.bit_reversed_order);
            [~, info_bits_sorted] = sort(channels, 'ascend');
            obj.info_bit_order = info_bits_sorted(obj.block_length:-1:1);
            obj.frozen_bits = ones(1,obj.block_length);
            obj.frozen_bits(info_bits_sorted(1:obj.info_length + obj.crc_size)) = 0;
            obj.info_bits = info_bits_sorted(1:obj.info_length + obj.crc_size);
            obj.llr_based_computation = 0;

            obj.cc_method = 'bhattacharya';
            obj.cc_parameter = design_epsilon;
            obj.cc_misc = '';

            disp(['Bhattacharya code construction done. BLER estimate:', num2str(sum(channels(obj.info_bits)))]);

        end

        function monte_carlo_code_construction(obj, design_snr_db, num_runs, constellation_name, receiver_algo)
            if (nargin < 3) || isempty(num_runs)
                num_runs = 100e3;
            end
            if (nargin < 4) || isempty(constellation_name)
                constellation_name = 'bpsk';
            end
            if (nargin < 5) || isempty(receiver_algo)
                receiver_algo = 'bicm';
            end

            obj.cc_method = 'monte-carlo';
            obj.cc_parameter = design_snr_db;
            obj.cc_misc = [constellation_name, '_', receiver_algo, '_', num2str(num_runs)];

            txt_file_name = ['CodeConstructionData/MC_block_length_', obj.get_unique_string(), '.txt'];

            if exist(txt_file_name, 'file')
                channels = load(txt_file_name);
                disp(['Found monte carlo BER pattern: ', txt_file_name]);
            else
                disp(['Didnt find monte carlo BER here: ', txt_file_name]);
                disp('Hence, determining the  monte carlo ber.');
                channels = obj.monte_carlo(design_snr_db, num_runs, constellation_name, receiver_algo);
                fileID = fopen(txt_file_name,'w');
                for c  = channels
                    fprintf(fileID,'%d \n',c);
                end
                fclose(fileID);
            end
            [~, channel_order] = sort(channels, 'ascend');
            if channels(channel_order(obj.info_length + obj.crc_size)) < 100
                disp('Warning: not enough runs to get a reliable code.');
                disp(['Worst # of errors = ', num2str(channels(channel_order(obj.info_length + obj.crc_size)))]);
            end

            obj.info_bits = channel_order(1:obj.info_length + obj.crc_size);

            obj.frozen_bits = ones(1,obj.block_length);
            obj.frozen_bits(obj.info_bits) = 0;
            bler_estimate = sum(channels(obj.info_bits))/num_runs;

            disp(['Monte carlo code construction done. Bler estimate = ', num2str(bler_estimate)]);

        end

        function num_err = monte_carlo(obj, design_snr_db, num_runs, constellation_name, receiver)

            num_err = zeros(obj.block_length, 1);
            modulation = Constellation(constellation_name);

            for i_run  = 1 : num_runs
                if mod(i_run, ceil(num_runs/10)) == 1
                    disp(['Code construction iteration running = ', num2str(i_run)]);
                end
                if strcmp(receiver, 'bicm')

                    dummy_info = (rand(1, obj.block_length) < 0.5);
                    dummy_coded = PolarCode.polar_encode(dummy_info);

                elseif strcmp(receiver, 'mlc')

                    num_codes = modulation.n_bits;
                    dummy_info = (rand(num_codes, obj.block_length/num_codes) < 0.5);
                    dummy_coded = zeros(obj.block_length, 1);
                    for layer = 1 : num_codes
                        dummy_coded(layer:num_codes:obj.block_length) = PolarCode.polar_encode(dummy_info(layer,:))';
                    end

                end

                mod_sym = modulation.modulate(dummy_coded);
                sigma = sqrt(1/2) *  10^(-design_snr_db/20);
                noise = sigma * randn(length(mod_sym), 1);
                y = mod_sym + noise;

                if strcmp(receiver, 'bicm')

                    p1 = 0.5 * ones(obj.block_length, 1);
                    eff_block_length = floor(obj.block_length/modulation.n_bits)*modulation.n_bits;
                    [p1(1:eff_block_length),~] = modulation.compute_llr_bicm(y, sigma^2);
                    [~, ber_tmp] = PolarCode.polar_decode_monte(p1, dummy_info);
                    num_err = num_err + ber_tmp';

                elseif strcmp(receiver, 'mlc')

                    decoded_coded = zeros(obj.block_length/num_codes, num_codes);
                    ber = zeros(obj.block_length/num_codes, num_codes);

                    for layer = 1 : num_codes
                        u = decoded_coded(:, 1:layer-1);
                        [p1,~] = modulation.compute_llr_mlc(y, sigma^2, u);
                        [decoded_coded(:, layer), ber(:, layer)] = PolarCode.polar_decode_monte(p1, dummy_info(layer, :));
                    end
                    ber = ber(:);
                    num_err = num_err + ber;

                end
            end
        end

        function [bler_estimate] = ga_code_construction(obj, design_snr_db, constellation_name, receiver_algo)

            modulation = Constellation(constellation_name);

            if strcmp(receiver_algo, 'bicm')

                if modulation.n_bits == 1
                    capacity = modulation.get_bicm_capacity(design_snr_db);
                else
                    capacity = modulation.get_polarized_capacity(design_snr_db);
                end

            elseif strcmp(receiver_algo, 'mlc')

                capacity = modulation.get_mlc_capacity(design_snr_db);

            end

            num_codes = modulation.n_bits;

            mean_llrs = get_bpsk_llr_for_capacity(capacity);

            llr_vec = zeros(obj.block_length, 1);

            if modulation.n_bits == 3  && strcmp(receiver_algo, 'mlc') %non power of 2 not supported
                disp('MLC code construction not supported');
            else
                for i_code = 1 : num_codes
                    bit_loc = (i_code-1) * obj.block_length/num_codes + 1 : i_code * obj.block_length/num_codes;
                    llr_vec(bit_loc) = mean_llrs(i_code);
                end
            end

            bit_rev_order_subcode = bitrevorder(1:obj.block_length/ num_codes);
            channels = zeros(obj.block_length, 1);

            for i_code = 1 : num_codes
                bit_loc = (i_code-1) * obj.block_length/num_codes + 1 : i_code * obj.block_length/num_codes;
                tmp_c = calculate_awgn_polarization(llr_vec((bit_loc)).', obj.n - log2(num_codes));
                channels(bit_loc) = tmp_c(bit_rev_order_subcode);
            end

            [~, channel_order] = sort(channels, 'descend');
            obj.info_bit_order = channel_order(obj.block_length:-1:1);

            obj.info_bits = channel_order(1:obj.info_length + obj.crc_size);
            obj.frozen_bits = ones(1,obj.block_length);
            obj.frozen_bits(obj.info_bits) = 0;
            bler_estimate =  sum(qfunc(  sqrt(channels(obj.info_bits))/sqrt(2)));

            obj.cc_method = 'gauss-approx';
            obj.cc_parameter = design_snr_db;
            obj.cc_misc = [constellation_name, '_', receiver_algo];

            disp(['GA code construction done. BLER estimate:', num2str(bler_estimate)]);

        end

        function [unique_string] = get_unique_string(obj)
            unique_string = [num2str(obj.block_length), '_', num2str(length(obj.info_bits)),  ...
                '_cc_method_', obj.cc_method, '_cc_param_', num2str(obj.cc_parameter), '_', obj.cc_misc];
        end

        %% encoder

        function coded_bits = encode(obj, info_bits)

            info_bits_padded = zeros(obj.block_length, 1);
            if obj.crc_size ~= 0
                crc = mod(obj.crc_matrix*info_bits', 2);
                info_bits = [info_bits, crc'];
            end
            info_bits_padded(obj.info_bits) = info_bits;
            coded_bits = PolarCode.polar_encode(info_bits_padded);

        end

        %%  Decoder helper

        function [b] = crc_check(obj, info_bits)
            info = info_bits(1:obj.info_length);
            crc  = info_bits(obj.info_length + 1:obj.info_length + obj.crc_size);
            crc_check =  mod(obj.crc_matrix*info', 2)';
            b = 1 - any(crc ~= crc_check);
        end

        %%  SC decoder

        function decoded_bits = decode_sc_p1(obj, p1)

            [decoded_bits, ~] = PolarCode.polar_decode(p1, obj.frozen_bits);
            decoded_bits = decoded_bits(obj.info_bits(1:obj.info_length));

        end

        %%  SCL decoder

        function [u] = decode_scl_p1(obj, p1, p0, list_size)

            obj.list_size = list_size;
            obj.llr_based_computation  = 0;
            obj.initializeDataStructures();
            l_index = obj.assignInitialPath();
            s_index = obj.getArrayPointer_P(0, l_index);
            obj.p_scl(obj.get_i_scl(0, 0, s_index) + 1:obj.get_i_scl(0, obj.block_length - 1, s_index) + 1 , 1) = p0;
            obj.p_scl(obj.get_i_scl(0, 0, s_index) + 1:obj.get_i_scl(0, obj.block_length - 1, s_index) + 1 , 2) = p1;
            u = obj.polar_decode_scl();

        end

        function [u] = decode_scl_llr(obj, llr, list_size)

            obj.list_size = list_size;
            obj.llr_based_computation = 1;
            obj.initializeDataStructures();
            l_index = obj.assignInitialPath();
            s_index = obj.getArrayPointer_P(0, l_index);
            obj.llr_scl( obj.get_i_scl(0, 0, s_index) + 1:obj.get_i_scl(0, obj.block_length - 1, s_index) + 1) = llr;
            u = obj.polar_decode_scl();

        end

        function [u] = polar_decode_scl(obj)

            for phi = 0 : obj.block_length - 1

                obj.recursivelyCalcP_scl(obj.n, phi);

                if obj.frozen_bits(phi + 1) == 1
                    obj.continuePaths_FrozenBit(phi);
                else
                    obj.continuePaths_UnfrozenBit(phi);
                end

                if mod(phi, 2) == 1
                    obj.recursivelyUpdateC_scl(obj.n, phi);
                end

            end

            l_index = obj.findMostProbablePath(1);
            c_m = obj.getArrayPointer_C(obj.n, l_index);
            info = obj.i_scl(c_m+1,:);
            u = info(obj.info_bits(1:obj.info_length));

        end

        function index = get_i_scl(obj, lambda, beta, list_index)
            index = beta + obj.lambda_offset(lambda + 1) + obj.list_offset(list_index + 1);
        end

        function init_i_scl(obj)
            obj.lambda_offset = (2.^(obj.n - (0: obj.n)) - 1);
            obj.list_offset = (0:obj.list_size)*(2 * obj.block_length - 1);

        end

        function initializeDataStructures(obj)

            obj.inactivePathIndices = zeros(obj.list_size,1);
            obj.inactivePathIndicesSize = 0;
            % the above two variables are used to define a stack

            obj.activePathArray =  zeros(obj.list_size,1);
            obj.pathIndexToArrayIndex = zeros(obj.n  + 1, obj.list_size);

            obj.inactiveArrayIndices = zeros(obj.n  + 1, obj.list_size);
            obj.inactiveArrayIndicesSize = zeros(obj.n + 1, 1);
            % the above two variables are used to define a vector of stacks

            obj.arrayReferenceCount = zeros(obj.n  + 1, obj.list_size);

            if obj.llr_based_computation
                obj.llr_scl = zeros(obj.list_size * (2 * obj.block_length - 1), 1);
                obj.llr_path_metric =  zeros(obj.list_size, 1);
            else
                obj.p_scl = zeros(obj.list_size * (2 * obj.block_length - 1), 2);
            end

            obj.c_scl = zeros(obj.list_size * (2 * obj.block_length - 1), 2);
            obj.i_scl = zeros(obj.list_size, obj.block_length);
            obj.init_i_scl();

            for lambda = 0 : obj.n
                for i_list = 0 : obj.list_size - 1
                    obj.inactiveArrayIndices(lambda + 1, i_list + 1) = i_list;

                end
                obj.inactiveArrayIndicesSize(lambda + 1) = obj.list_size;
            end

            for i_list = 0 : obj.list_size - 1
                obj.activePathArray(i_list + 1) = 0;
                obj.inactivePathIndices(i_list + 1) = i_list;
            end

            obj.inactivePathIndicesSize  = obj.list_size;

        end

        function l_index = assignInitialPath(obj)
            l_index = obj.inactivePathIndices(obj.inactivePathIndicesSize);
            obj.inactivePathIndicesSize = obj.inactivePathIndicesSize - 1;
            obj.activePathArray(l_index + 1) = 1;

            for lambda = 0 : obj.n
                s = obj.inactiveArrayIndices(lambda + 1, obj.inactiveArrayIndicesSize(lambda + 1));
                obj.inactiveArrayIndicesSize(lambda + 1) = obj.inactiveArrayIndicesSize(lambda + 1) - 1;
                obj.pathIndexToArrayIndex(lambda + 1, l_index + 1) = s;
                obj.arrayReferenceCount(lambda + 1, l_index + 1) = 1;
            end
        end

        function l_p_index = clonePath(obj, l_index)
            l_p_index = obj.inactivePathIndices(obj.inactivePathIndicesSize);
            obj.inactivePathIndicesSize = obj.inactivePathIndicesSize - 1;
            obj.activePathArray(l_p_index + 1) = 1;

            if obj.llr_based_computation
                obj.llr_path_metric(l_p_index + 1) = obj.llr_path_metric(l_index + 1);
            end

            for lambda = 0 : obj.n
                s = obj.pathIndexToArrayIndex(lambda + 1, l_index + 1);
                obj.pathIndexToArrayIndex(lambda + 1, l_p_index + 1) = s;
                obj.arrayReferenceCount(lambda + 1, s + 1) = obj.arrayReferenceCount(lambda + 1, s + 1) + 1;
            end
        end

        function killPath(obj, l_index)
            obj.activePathArray(l_index + 1) = 0;
            obj.inactivePathIndices(obj.inactivePathIndicesSize + 1) = l_index;
            obj.inactivePathIndicesSize = obj.inactivePathIndicesSize  + 1;
            if obj.llr_based_computation
                obj.llr_path_metric(l_index + 1) = 0;
            end
            for lambda = 0 : obj.n
                s = obj.pathIndexToArrayIndex(lambda + 1, l_index + 1);
                obj.arrayReferenceCount(lambda + 1, s + 1) = obj.arrayReferenceCount(lambda + 1, s + 1) - 1;
                if obj.arrayReferenceCount(lambda + 1, s + 1) == 0
                    obj.inactiveArrayIndices(lambda + 1, obj.inactiveArrayIndicesSize(lambda + 1) + 1) = s;
                    obj.inactiveArrayIndicesSize(lambda + 1)  = obj.inactiveArrayIndicesSize(lambda + 1)  + 1;
                end
            end
        end

        function [b] = pathIndexInactive(obj, l_index)
            b = 1 - obj.activePathArray(l_index + 1);
        end

        function [s_p] = getArrayPointer_P(obj, lambda, l_index)
            s = obj.pathIndexToArrayIndex(lambda + 1, l_index + 1);
            m = obj.n;
            if obj.arrayReferenceCount(lambda + 1, s + 1) == 1
                s_p = s;
            else
                s_p = obj.inactiveArrayIndices(lambda + 1, obj.inactiveArrayIndicesSize(lambda + 1));
                i_s_p = obj.lambda_offset(lambda + 1) + obj.list_offset(s_p + 1) + 1: ...
                    obj.lambda_offset(lambda + 1) + obj.list_offset(s_p + 1) + 2^(m - lambda);
                i_s = obj.lambda_offset(lambda + 1) + obj.list_offset(s + 1) + 1: ...
                    obj.lambda_offset(lambda + 1) + obj.list_offset(s + 1) + 2^(m - lambda);
                obj.c_scl(i_s_p, :) = obj.c_scl(i_s, :);
                if obj.llr_based_computation
                    obj.llr_scl(i_s_p) = obj.llr_scl(i_s);
                else
                    obj.p_scl(i_s_p, :) = obj.p_scl(i_s, :);
                end
                obj.inactiveArrayIndicesSize(lambda + 1) = obj.inactiveArrayIndicesSize(lambda + 1) - 1;
                obj.arrayReferenceCount(lambda + 1, s + 1) = obj.arrayReferenceCount(lambda + 1, s + 1) - 1;
                obj.arrayReferenceCount(lambda + 1, s_p + 1) = 1;
                obj.pathIndexToArrayIndex(lambda + 1, l_index + 1) = s_p;
            end
        end

        function [s_p] = getArrayPointer_C(obj, lambda, l_index)
            s = obj.pathIndexToArrayIndex(lambda + 1, l_index + 1);
            m = obj.n;
            if obj.arrayReferenceCount(lambda + 1, s + 1) == 1
                s_p = s;
            else
                s_p = obj.inactiveArrayIndices(lambda + 1, obj.inactiveArrayIndicesSize(lambda + 1));
                i_s_p = obj.lambda_offset(lambda + 1) + obj.list_offset(s_p + 1) + 1: ...
                    obj.lambda_offset(lambda + 1) + obj.list_offset(s_p + 1) + 2^(m - lambda);
                i_s = obj.lambda_offset(lambda + 1) + obj.list_offset(s + 1) + 1: ...
                    obj.lambda_offset(lambda + 1) + obj.list_offset(s + 1) + 2^(m - lambda);
                if obj.llr_based_computation
                    obj.llr_scl(i_s_p) = obj.llr_scl(i_s);
                else
                    obj.p_scl(i_s_p, :) = obj.p_scl(i_s, :);
                end
                obj.c_scl(i_s_p, :) = obj.c_scl(i_s, :);
                obj.inactiveArrayIndicesSize(lambda + 1) = obj.inactiveArrayIndicesSize(lambda + 1) - 1;
                obj.arrayReferenceCount(lambda + 1, s + 1) = obj.arrayReferenceCount(lambda + 1, s + 1) - 1;
                obj.arrayReferenceCount(lambda + 1, s_p + 1) = 1;
                obj.pathIndexToArrayIndex(lambda + 1, l_index + 1) = s_p;
            end
        end

        function recursivelyCalcP_scl(obj, lambda, phi)

            if lambda == 0
                return;
            end

            psi = floor(phi/2);
            if mod(phi, 2) == 0
                obj.recursivelyCalcP_scl(lambda - 1, psi);
            end

            sigma  = 0;
            p_index_3_base_list = zeros(obj.list_size, 1);
            l_index_1_list = zeros(obj.list_size, 1);
            for l_index = 0 : obj.list_size - 1
                if obj.activePathArray(l_index + 1) == 0
                    continue;
                end
                l_index_1_list(l_index+1) = obj.getArrayPointer_P(lambda, l_index);
                l_index_2 = obj.getArrayPointer_P(lambda - 1, l_index);
                l_index_3 = obj.getArrayPointer_C(lambda, l_index);

                p_index_1_base = obj.lambda_offset(lambda) + obj.list_offset(l_index_2 + 1) + 1;
                p_index_3_base_list(l_index+1) = obj.lambda_offset(lambda + 1) + obj.list_offset(l_index_1_list(l_index+1) + 1) + 1;
                c_index_3_base = obj.lambda_offset(lambda + 1) + obj.list_offset(l_index_3 + 1) + 1;
                for beta = 0: 2^(obj.n - lambda) - 1
                    p_index_1 = p_index_1_base + 2 * beta;
                    p_index_2 = p_index_1_base + 2 * beta + 1;
                    p_index_3 = p_index_3_base_list(l_index+1) + beta;
                    if mod(phi, 2) == 0
                        if obj.llr_based_computation
                            if max( abs(obj.llr_scl ( p_index_1)), abs(obj.llr_scl ( p_index_2)) ) < 40
                                obj.llr_scl(p_index_3) = PolarCode.cnop_llr( obj.llr_scl ( p_index_1), obj.llr_scl ( p_index_2));
                                % log( (exp( obj.llr_scl ( p_index_1) + obj.llr_scl ( p_index_2)) + 1) ...
                                %    /(exp( obj.llr_scl ( p_index_1))  + exp( obj.llr_scl ( p_index_2))) );
                            else
                                obj.llr_scl(p_index_3) = sign( obj.llr_scl ( p_index_1)) * sign(obj.llr_scl ( p_index_2)) * min(abs(obj.llr_scl ( p_index_2)), abs(obj.llr_scl ( p_index_1)));
                            end
                        else
                            obj.p_scl(p_index_3, 1) =  0.5 * ( obj.p_scl ( p_index_1 , 1) *  obj.p_scl(p_index_2, 1)  +  obj.p_scl ( p_index_1 , 2) *  obj.p_scl(p_index_2, 2));
                            obj.p_scl(p_index_3, 2) =  0.5 * ( obj.p_scl ( p_index_1 , 2) *  obj.p_scl(p_index_2, 1)  +  obj.p_scl ( p_index_1 , 1) *  obj.p_scl(p_index_2, 2));
                            sigma = max(obj.p_scl(p_index_3, 1), sigma);
                            sigma = max(obj.p_scl(p_index_3, 2), sigma);
                        end
                    else
                        u_p = obj.c_scl( c_index_3_base + beta, 1);
                        if obj.llr_based_computation
                            obj.llr_scl(p_index_3) = (-1)^u_p * obj.llr_scl(p_index_1) +  obj.llr_scl(p_index_2);
                        else
                            obj.p_scl(p_index_3, 1)  = 0.5 * obj.p_scl (p_index_1, mod(u_p, 2) + 1) * ...
                                obj.p_scl(p_index_2, 1);
                            obj.p_scl(p_index_3, 2)  = 0.5 * obj.p_scl (p_index_1, mod(u_p + 1, 2) + 1) * ...
                                obj.p_scl(p_index_2, 2);
                            sigma = max(obj.p_scl(p_index_3, 1), sigma);
                            sigma = max(obj.p_scl(p_index_3, 2), sigma);
                        end
                    end
                end
            end

            for l_index = 0 : obj.list_size - 1
                if sigma == 0 %typically happens because of underflow
                    break;
                end
                if obj.activePathArray(l_index + 1) == 0
                    continue;
                end
                p_range = p_index_3_base_list(l_index+1): p_index_3_base_list(l_index+1) +  2^(obj.n - lambda) - 1;
                obj.p_scl(p_range, :) = obj.p_scl(p_range, :)/sigma;
            end

        end

        function recursivelyUpdateC_scl(obj, lambda, phi)
            if mod(phi, 2) == 0
                disp('Error: phi should always be odd in this function call');
            end
            psi = floor(phi/2);

            for l_index = 0 : obj.list_size - 1
                if obj.activePathArray(l_index + 1) == 0
                    continue;
                end

                l_index_1 = obj.getArrayPointer_C(lambda, l_index);
                l_index_2 = obj.getArrayPointer_C(lambda - 1, l_index);

                p_index_1 = obj.lambda_offset(lambda + 1) + obj.list_offset(l_index_1 + 1)  + 1;
                p_index_2 = obj.lambda_offset(lambda) + obj.list_offset(l_index_2 + 1)  +  1;

                for beta = 0: 2^(obj.n - lambda) - 1
                    obj.c_scl(p_index_2 + 2*beta, mod(psi, 2) + 1) = ...
                        mod( obj.c_scl(p_index_1 + beta, 1) + ...
                        obj.c_scl(p_index_1 + beta, 2),  2);
                    obj.c_scl(p_index_2 + 2 * beta + 1, mod(psi, 2) + 1) = ...
                        obj.c_scl(p_index_1 + beta, 2);
                end

            end

            if mod(psi, 2) == 1
                obj.recursivelyUpdateC_scl(lambda - 1, psi);
            end

        end

        function continuePaths_FrozenBit(obj, phi)

            for l_index = 0 : obj.list_size -1

                if obj.activePathArray(l_index + 1) == 0
                    continue;
                end

                l_index_1 = obj.getArrayPointer_C(obj.n, l_index);
                obj.c_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1, mod(phi, 2) + 1) = 0;
                if obj.llr_based_computation
                    obj.llr_path_metric(l_index_1 + 1) = obj.llr_path_metric(l_index_1 + 1) ...
                        + log(1 + exp(-obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1)));
                end
            end

        end

        function    continuePaths_UnfrozenBit(obj, phi)

            probForks = -realmax * ones(obj.list_size, 2);
            index = 0;
            for l_index = 0 : obj.list_size - 1

                if obj.activePathArray(l_index + 1)
                    l_index_1 = obj.getArrayPointer_P(obj.n, l_index);
                    if obj.llr_based_computation
                        % computing negative of path metric so that an
                        % ascending order can be used for sorting
                        probForks(l_index + 1, 1) =  - (obj.llr_path_metric(l_index_1 + 1) ...
                            + log(1 + exp(-obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1))));
                        probForks(l_index + 1, 2) =  - ( obj.llr_path_metric(l_index_1 + 1) ...
                            + log(1 + exp(obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1))));
                    else
                        probForks(l_index + 1, 1) = obj.p_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1, 1);
                        probForks(l_index + 1, 2) = obj.p_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1, 2);
                    end
                    index = index + 1;

                end
            end

            rho = min(2*index, obj.list_size);
            contForks = zeros(obj.list_size, 2);
            prob = sort(probForks(:), 'descend');

            threshold = prob(rho);
            num_populated = 0;
            for l_index = 0 : obj.list_size - 1
                for j_index = 1 : 2
                    if num_populated == rho
                        break;
                    end
                    if  probForks(l_index + 1, j_index) > threshold
                        contForks(l_index + 1, j_index) = 1;
                        num_populated = num_populated + 1;
                    end
                end
            end

            if num_populated < rho
                for l_index = 0 : obj.list_size - 1
                    for j_index = 1 : 2
                        if num_populated == rho
                            break;
                        end
                        if  probForks(l_index + 1, j_index) == threshold
                            contForks(l_index + 1, j_index) = 1;
                            num_populated = num_populated + 1;
                        end
                    end
                end
            end

            for l_index = 0 : obj.list_size - 1
                if obj.activePathArray(l_index + 1) == 0
                    continue;
                end

                if (contForks(l_index + 1, 1) == 0) && ( contForks(l_index + 1, 2) == 0)
                    obj.killPath(l_index);
                end
            end

            for l_index = 0 : obj.list_size - 1
                if (contForks(l_index + 1, 1) == 0) && ( contForks(l_index + 1, 2) == 0)
                    continue;
                end
                l_index_1 = obj.getArrayPointer_C(obj.n, l_index);
                if (contForks(l_index + 1, 1) == 1) && ( contForks(l_index + 1, 2) == 1)
                    obj.c_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1, mod(phi,2) + 1) = 0;
                    obj.i_scl(l_index_1 + 1, phi + 1) = 0;

                    l_p = obj.clonePath(l_index);

                    l_index_2 = obj.getArrayPointer_C(obj.n, l_p);
                    obj.i_scl(l_index_2 + 1, 1 : phi) = obj.i_scl(l_index_1 + 1, 1 : phi);
                    obj.c_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_2 + 1) + 1, mod(phi,2) + 1) = 1;
                    obj.i_scl(l_index_2 + 1, phi + 1) = 1;
                    if obj.llr_based_computation
                        obj.llr_path_metric(l_index + 1) = obj.llr_path_metric(l_index + 1) ...
                            + log(1 + exp(-obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1)));
                        obj.llr_path_metric(l_p + 1) = obj.llr_path_metric(l_p + 1) ...
                            + log(1 + exp(obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_2 + 1) + 1)));
                    end

                else
                    if contForks(l_index + 1, 1) == 1
                        obj.c_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1, mod(phi,2) + 1) = 0;
                        obj.i_scl(l_index_1 + 1, phi + 1) = 0;
                        if obj.llr_based_computation
                            obj.llr_path_metric(l_index + 1) = obj.llr_path_metric(l_index + 1) ...
                                + log(1 + exp(-obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1)));
                        end
                    else
                        obj.c_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1, mod(phi,2) + 1) = 1;
                        obj.i_scl(l_index_1 + 1, phi + 1) = 1;
                        if obj.llr_based_computation
                            obj.llr_path_metric(l_index + 1) = obj.llr_path_metric(l_index + 1) ...
                                + log(1 + exp(obj.llr_scl(obj.lambda_offset(obj.n + 1) + obj.list_offset(l_index_1 + 1) + 1)));
                        end
                    end
                end

            end

        end

        function [l_p_index] = findMostProbablePath(obj, crc_check)
            l_p_index = 0;

            p_max = 0;
            if obj.llr_based_computation
                p_max = realmax;
            end
            path_with_crc = 0;
            for l_index = 0 : obj.list_size -1

                if obj.activePathArray(l_index + 1) == 0
                    continue;
                end

                c_index = obj.getArrayPointer_C( obj.n, l_index);
                if (crc_check) && (obj.crc_size ~= 0)
                    a = obj.i_scl(c_index+1,:);
                    u = a(obj.info_bits);
                    if obj.crc_check(u) == 0
                        continue;
                    end
                end
                path_with_crc = 1;
                if obj.llr_based_computation
                    if p_max > obj.llr_path_metric(l_index + 1)
                        p_max = obj.llr_path_metric(l_index + 1);
                        l_p_index = l_index;
                    end
                else
                    p_index = obj.getArrayPointer_P( obj.n, l_index);
                    if p_max < obj.p_scl(obj.get_i_scl( obj.n, 0, p_index) + 1, obj.c_scl(obj.get_i_scl( obj.n, 0, c_index) + 1, 2) + 1)
                        l_p_index = l_index;
                        p_max  = obj.p_scl(obj.get_i_scl( obj.n, 0, p_index) + 1, obj.c_scl(obj.get_i_scl( obj.n, 0, c_index) + 1, 2) + 1);
                    end
                end
            end

            if (crc_check) && (path_with_crc == 0) % no path with crc check found
                l_p_index = obj.findMostProbablePath(0);
            end

        end

        %% helper function
        function [ bler, ber ] = get_bler_quick(obj, ebno_vec, list_size_vec)

            snr_db_vec = ebno_vec + 10*log10(obj.info_length/obj.block_length);

            num_block_err = zeros(length(ebno_vec), length(list_size_vec));
            num_bit_err = zeros(length(ebno_vec), length(list_size_vec));
            num_runs = zeros(length(ebno_vec), length(list_size_vec));
            max_err = 50;
            max_runs = 500;

            for i_run = 1 : max_runs

                if mod(i_run, ceil(max_runs/10)) == 1
                    disp(['Sim iteration running = ', num2str(i_run)]);
                end
                info = rand(1 , obj.info_length) < 0.5;
                coded_bits = obj.encode(info);
                bpsk = 2 * coded_bits - 1;
                sigma = sqrt(1/2);
                noise = sigma * randn(1, obj.block_length);
                prev_decoded = zeros(length(list_size_vec), length(ebno_vec));

                for i_ebno = 1 : length(ebno_vec)
                    snr_db = snr_db_vec(i_ebno);
                    received_bits = 10^(snr_db/20) * bpsk + noise;
                    p1 = exp(-(received_bits - 10^(snr_db/20)).^2/(2 * sigma^2))/sigma/sqrt(2*pi);
                    p0 = exp(-(received_bits + 10^(snr_db/20)).^2/(2 * sigma^2))/sigma/sqrt(2*pi);

                    for i_list = 1 : length(list_size_vec)
                        if num_block_err(i_ebno, i_list) > max_err
                            continue;
                        end

                        num_runs(i_ebno,  i_list) =  num_runs(i_ebno, i_list) + 1;

                        run_sim = 1;

                        for i_ebno2 = 1 : i_ebno
                            if prev_decoded(i_list, i_ebno2)
                                run_sim = 0;
                            end
                        end

                        if (run_sim == 0)
                            % This is a hack to speed up simulations --
                            % it assumes that this run will be decoded correctly since it was
                            % decoded correctly for a lower EbNo
                            continue;
                        end
                        if list_size_vec(i_list) == 1
                            decoded_bits = obj.decode_sc_p1(p1./(p1+p0));
                        else
                            %  decoded_bits = obj.decode_scl_p1(p1, p0, list_size_vec(i_list));
                            decoded_bits = obj.decode_scl_llr(log(p0./p1), list_size_vec(i_list));
                        end
                        err = any(info ~= decoded_bits);
                        if err
                            num_block_err(i_ebno, i_list) =  num_block_err(i_ebno, i_list) + 1;
                            num_bit_err(i_ebno, i_list) = num_bit_err(i_ebno, i_list) + sum(info ~= decoded_bits);
                        else
                            prev_decoded(i_list, i_ebno) = 1;
                        end
                    end
                end

            end
            bler = num_block_err./num_runs;
            ber = num_bit_err./num_runs;

        end

    end

    methods (Static)
        % the next four functions are modified version of code in:
        % http://pfister.ee.duke.edu/courses/ecen655/polar.pdf
        function x = polar_encode(info_bits_padded)

            N = length(info_bits_padded);
            if (N == 1)
                x = info_bits_padded;
            else
                u1u2 = mod(info_bits_padded(1:2:end) + info_bits_padded(2:2:end) , 2);
                u2 = info_bits_padded(2:2:end);
                x = [PolarCode.polar_encode(u1u2) PolarCode.polar_encode(u2)];
            end

        end

        function [u, x] = polar_decode(y,f)
            N = length(y);
            if (N==1)
                if (f == 0)  % info bit
                    x = (1-sign(1-2*y))/2;
                else  % frozen bit -- assumed to be zero
                    x = 0;
                end
                u = x;
            else
                u1est = PolarCode.cnop(y(1:2:end),y(2:2:end));
                [uhat1,u1hardprev] = PolarCode.polar_decode(u1est,f(1:N/2));
                u2est = PolarCode.vnop(PolarCode.cnop(u1hardprev,y(1:2:end)),y(2:2:end));
                [uhat2,u2hardprev] = PolarCode.polar_decode(u2est,f(N/2+1:end));
                u = [uhat1 uhat2];
                x = reshape([PolarCode.cnop(u1hardprev,u2hardprev); u2hardprev],1,[]);
            end
        end

        function z = cnop(w1,w2) % notation = probability of bit = 1
            z = w1.*(1-w2) + w2.*(1-w1);
        end

        function z = vnop(w1,w2) % notation = probability of bit = 1
            z = w1.*w2 ./ (w1.*w2 + (1-w1).*(1-w2));
        end

        function [x, ber] = polar_decode_monte(y, dummy_info)
            N = length(y);
            if (N==1)
                if (y > 0.5 && dummy_info == 1) || (y <= 0.5 && dummy_info == 0)
                    ber = 0;
                else
                    ber = 1;
                end
                x = dummy_info;
            else
                u1est = PolarCode.cnop(y(1:2:end),y(2:2:end));
                [u1hardprev, ber1] = PolarCode.polar_decode_monte(u1est,dummy_info(1:N/2));
                u2est = PolarCode.vnop(PolarCode.cnop(u1hardprev',y(1:2:end)),y(2:2:end));
                [u2hardprev, ber2] = PolarCode.polar_decode_monte(u2est,dummy_info(N/2+1:end));
                ber = [ber1 ber2];
                x = reshape([PolarCode.cnop(u1hardprev,u2hardprev); u2hardprev],1,[]);
            end
        end

        function [x, u] = polar_decode_capacity(y, dummy_info)
            N = length(y);
            if (N==1)
                u = y;
                x = dummy_info;
            else
                u1est = PolarCode.cnop(y(1:2:end),y(2:2:end));
                [u1hardprev, u1] = PolarCode.polar_decode_capacity(u1est,dummy_info(1:N/2));
                u2est = PolarCode.vnop(PolarCode.cnop(u1hardprev',y(1:2:end)),y(2:2:end));
                [u2hardprev, u2] = PolarCode.polar_decode_capacity(u2est,dummy_info(N/2+1:end));
                u = [u1 u2];
                x = reshape([PolarCode.cnop(u1hardprev,u2hardprev); u2hardprev],1,[]);
            end
        end

        function [x, u] = polar_decode_capacity_llr(y, dummy_info)
            N = length(y);
            if (N==1)
                u = y;
                x = dummy_info;
            else
                u1est = PolarCode.cnop_llr(y(1:2:end),y(2:2:end));
                [u1hardprev, u1] = PolarCode.polar_decode_capacity_llr(u1est,dummy_info(1:N/2));
                u2est = PolarCode.vnop_llr((1 - 2*u1hardprev').*y(1:2:end),y(2:2:end));
                [u2hardprev, u2] = PolarCode.polar_decode_capacity_llr(u2est,dummy_info(N/2+1:end));
                u = [u1 u2];
                x = reshape([PolarCode.cnop(u1hardprev,u2hardprev); u2hardprev],1,[]);
            end
        end

        function l = cnop_llr(l1,l2)
            l = 2 * atanh(tanh(l1/2).*tanh(l2/2));
        end

        function l = vnop_llr(l1,l2)
            l = l1 + l2;
        end

        function channels  = calculate_channel_polarization( epsilon, n)

            N  = 2^n;
            log_domain_enabled = 1;
            if log_domain_enabled
                epsilon = log(epsilon);
            end
            if length(epsilon) == 1
                channels = epsilon*ones(1, N);
            elseif length(epsilon) == N
                channels = epsilon;
            end

            for i = 1:n
                c1 = channels(1:2:N);
                c2 = channels(2:2:N);
                if log_domain_enabled
                    channels = [log(exp(c1) + exp(c2) - exp(c1+c2)), c1 + c2];
                else
                    channels = [c1 + c2 - c1.*c2, c1.*c2];
                end
            end

            if log_domain_enabled
                channels = exp(channels);
            end

        end

    end

end

author: tavildar

repo: Polar

save_path: github-repos/MATLAB/tavildar-Polar


text:
function [rmain,r1]= ReasInteract(var1,newcat,rfact,xmeff)

	%ORIGINAL FROM ZMAP: NEEDS WORK FOR MAPSEIS
	%SUBFUNCTION: ReasenbergDeclus.m MIGHT NOT BE NEEDED
	%---------------------------------------------------

	%--------->NOT REALLY NEEDED

	%interact.m                                        A.Allmann
	% calculates the interaction zones of the earthquakes
	% in [km]
	%Last modification 6/95

	%global newcat rfact xmeff

	if var1==1
		rmain = 0.011*10.^(0.4*newcat(:,6)); %interaction zone for mainshock
		%rmain = 10.^(-2.44+(.59*newcat(:,6)));

		%tm1=find(rmain==0.011);             %these eqs got no magnitude in the catalog
		%tm2= 0.011*10^(0.4*xmeff);          %assume that for eqs with magnitude 0
		%rmain(tm1)=tm2*ones(1,length(tm1)); %the real magnitude is around xmeff

		r1    =rfact*rmain;                  %interaction zone if included in a cluster

	end

	%That function can be directly integrated it consist out of two functions anyway (DE)

end

author: CelsoReyes


text:
function result=edge_association(A,B,F)
% Author:  Qu Xiao-Bo    <quxiaobo [at] xmu.edu.cn>  and Hu Chang-wei Hu , June 26, 2009
%          Postal address:
% Rom 509, Scientific Research Building # 2,Haiyun Campus, Xiamen University,Xiamen,Fujian, P. R. China, 361005
% Website: http://quxiaobo.go.8866.org

%reference paper:Objective image fusion performance measure
A=double(A);
B=double(B);
F=double(F);
[row,column]=size(A);
[gA,aA]=get_g_a(A);
[gB,aB]=get_g_a(B);
[gF,aF]=get_g_a(F);
GAF=get2G(gA,gF);
GBF=get2G(gB,gF);
AAF=get2A(aA,aF);
ABF=get2A(aB,aF);
QgAF=getQg(GAF,0.9994,-15,0.5);
QgBF=getQg(GBF,0.9994,-15,0.5);
QaAF=getQa(AAF,0.9879,-22,0.8);
QaBF=getQa(ABF,0.9879,-22,0.8);
QAF=getQ(QgAF,QaAF);
QBF=getQ(QgBF,QaBF);
a=sum(sum(QAF.*gA+QBF.*gB));
b=sum(sum(gA+gB));
result=a/(b+eps);
%%
function [g,a]=get_g_a(im)

s1=[1 2 1; 0 0 0; -1 -2 -1];
s2=[-1 0 1;-2 0 2;-1 0 1];
sx=conv2(im,s1,'same');
sy=conv2(im,s2,'same');
g=sqrt(sx.^2+sy.^2);
a=atan(sy./(sx+eps));

end
%%
function Q=getQ(Qg,Qa)
Q=Qg.*Qa;
end
%%
function Qa=getQa(A,Ta,Ka,Oa)
Qa=Ta./(1+exp(Ka.*(A-Oa)));
end
%%
function Qg=getQg(G,Tg,Kg,Og)
Qg=Tg./(1+exp(Kg.*(G-Og)));
end
%%
function A=get2A(aIm1,aIm2)
A=1-abs(((aIm1-aIm2).*2)./pi);
end
%%
function G=get2G(gIm1,gIm2)
G=min(gIm1,gIm2)./(max(gIm1,gIm2)+eps);
end
%%
end

author: thfylsty


text:
%SIMILARRECTS  Class for grouping object candidates, detected by Cascade Classifier, HOG etc.
%
%     b = cv.SimilarRects(r1, r2)
%     b = cv.SimilarRects(r1, r2, 'OptionName',optionValue, ...)
%
% ## Input
% * __r1__ First input rectangle, 4-element vector `[x,y,w,h]`.
% * __r2__ Second input rectangle.
%
% ## Output
% * __b__ logical output indicate whether rectangles are similar.
%
% ## Options
% * __EPS__ Relative difference between sides of the rectangles to consider
%   similar. default 0.2
%
% See also: cv.groupRectangles, cv.groupRectangles_meanshift
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [seriesC, typeC] = extract_all_series(dcmdir_PATIENT)
%"extract_all_series"
%   Returns a cell array of all series contained within a dcmdir.PATIENT
%   datastructure, regardless of parent study.  Also returns an optional
%   type cell, containing strings describing the type of each series in
%   seriesC.
%
%JRA 06/15/06
%YWU modified 03/03/08
%
%Usage:
%   seriesC = extract_all_scan_series(dcmdir_PATIENT);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

seriesC     = {};
typeC       = {};

for studyNum = 1:length(dcmdir_PATIENT.STUDY)

    study = dcmdir_PATIENT.STUDY(studyNum); %wy {} --> ()

    for seriesNum = 1:length(study.SERIES)

        seriesC{end+1} = study.SERIES(seriesNum); %wy {} --> ()

        %Remove the info field, leaving only the field describing what type
        %of series this is.
       % disp(study.SERIES(seriesNum));
        tmp = rmfield(study.SERIES(seriesNum), 'info'); %wy {} --> ()
        tmpFields = fieldnames(tmp);
        if length(tmpFields) == 2; %wy 1-->2
            imgModality = tmp.Modality;
            imgModality = parseModality(imgModality);
            typeC{end+1} = imgModality;
            %typeC{end+1} = tmp.Modality; % wy tmpFields{:};
        else
           error('Improperly formed dcmdir_PATIENT structure passed to extract_all_series.');
        end

    end

end

return;

function imgModality = parseModality(imgModality)
%function imgModality = parseModality(imgModality)
%
% This function parses image modality and returns valid DICOM modality. For
% example, "CT SCAN", "M_CT" will be considered as belonging to modality "CT"

imgModality = upper(imgModality);

indStruct = strfind(imgModality,'STRUCT');
indDose = strfind(imgModality,'DOSE');
indPlan = strfind(imgModality,'PLAN');

if ~isempty([indStruct indDose indPlan])
    return;
end
indCT = strfind(imgModality,'CT');
indPT = strfind(imgModality,'PT');
indPET = strfind(imgModality,'PET');

if ~isempty(indCT) && isempty([indPT indPET])
    imgModality = 'CT';
elseif ~isempty(indPT) && isempty([indCT indPET])
    imgModality = 'PT';
elseif ~isempty(indPET) && isempty([indCT indPT])
    imgModality = 'PT';
end

return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [C_df,Df] = extract_DF_F(Y,A,C,P,options,AY)

% extract DF/F signals after performing NMF
% inputs:  Y raw data (d X T matrix, d # number of pixels, T # of timesteps)
%          A matrix of spatial components (d x K matrix, K # of components)
%          C matrix of temporal components (K x T matrix)
%          P neuron structure, used to read the baseline activity for each
%                    component of C
%          options structure used for specifying method for determining DF
%           default method is the median of the trace. By changing
%           options.df_prctile an arbitray percentile can be used (between 0 and 100).
%           a moving window can also be established by specifying options.df_window
%          AY pass is the product A'*Y otherwise compute it

% outputs:  C_df temporal components in the DF/F domain
%           Df   background for each component to normalize the filtered raw data

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2016

%memmaped = isobject(Y);
defoptions = CNMFSetParms;
if nargin < 5 || isempty(options)
    options = defoptions;
end
if ~isfield(options,'df_prctile') || isempty(options.df_prctile)
    options.df_prctile = defoptions.df_prctile;
end
if ~isfield(options,'df_window') || isempty(options.df_window)
    options.df_window = defoptions.df_window;
end
if ~isfield(options,'full_A') || isempty(options.full_A); full_A = defoptions.full_A; else full_A = options.full_A; end

[K,T] = size(C);

if ~exist('AY','var')
    AY = mm_fun(A,Y);
end

Bas = zeros(K,T);

if ~(nargin < 5 || isempty(P))
    bas_val = cell2mat(P.b);
    Ntr = size(bas_val,2);
    if Ntr > 1
        ln = diff(P.cs_frtrs);
        for i = 1:Ntr
            Bas(:,:,P.cs_frtrs(i)+1:P.cs_frtrs(i+1)) = repmat(bas_val(:,i),1,ln(i));
        end
    else
        Bas = repmat(bas_val,1,T);
    end
end

nA = sqrt(sum(A.^2))';
AA = A'*A;
AA(1:K+1:end) = 0;

Cf = bsxfun(@times,C - Bas,nA(:).^2);

C2 = AY - AA*C;

if isempty(options.df_window) || (options.df_window > size(C,2))
    if options.df_prctile == 50
        Df = median(C2,2);
    else
        Df = prctile(C2,options.df_prctile,2);
    end
    C_df = bsxfun(@times,Cf,1./Df(:));
else
    if options.df_prctile == 50
        if verLessThan('matlab','2015b')
            warning('Median filtering at the boundaries might be inaccurate due to zero padding.')
            Df = medfilt1(C2,options.df_window,[],2);
        else
            Df = medfilt1(C2,options.df_window,[],2,'truncate');
        end
    else
        Df = zeros(size(C2));
        for i = 1:size(Df,1);
            df_temp = running_percentile(C2(i,:), options.df_window, options.df_prctile);
            Df(i,:) = df_temp(:)';
        end
    end
    C_df = Cf./Df;
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function WriteLDRStack(stack, name, format)
%
%       WriteLDRStack(dir_name, format)
%
%       This function writes an LDR stack into a set of LDR images
%
%        Input:
%           -stack: the stack to be written
%           -name: the folder name where the stack is stored.
%           -format: an LDR format for reading LDR images.
%
%        Output:
%           -exposure: a stack of exposure values from images in dir_name
%
%     Copyright (C) 2019  Francesco Banterle
%
%     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/>.
%

n = size(stack, 4);

for i=1:n
    imwrite(stack(:,:,:,i), [name, '_', num2str(100000 + i), '.', format]);
end

end

author: banterle


text:
clear
close all

% Set defaults
font_size = 11;
set(0, 'DefaultAxesFontSize',font_size)

line_width = 1;
set(0, 'DefaultLineLineWidth',line_width)

font_size = 15;
set(0, 'DefaultTextFontSize',font_size)

subplot(2,1,1)

dualtree_eg1

box off

% remove defaults
set(0, 'DefaultAxesFontSize','remove')
set(0, 'DefaultLineLineWidth','remove')
set(0, 'DefaultTextFontSize','remove')

if 1
    % paperposition
    paperpos = [2 2 7 5];    % [left, bottom, width, height]

    set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperPosition', paperpos);
end

print -depsc figures/dualtree_eg1_figure

% then open eps file and save as png file.

author: thomaskuestner


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtMshRefine.m                                |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Triangular mesh of a planar square            |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Create mesh
Nvtx = 1e3;
mesh = mshSquare(Nvtx,[2 2]);

% Colours
mesh.col(1:20)                 = -1;
mesh.col(length(mesh)+(-20:0)) = 1;

% Graphical representation
figure
plot(mesh)
hold on
plotNrm(mesh,'r')
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')

% Refine all element with midpoint algorithm
tic
meshr = midpoint(mesh);
toc

% Graphical representation
figure
plot(meshr)
hold on
plotNrm(meshr,'r')
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')

% Midpoint algorithm for selected indices
I = find(mesh.col~=0);
tic
meshr = midpoint(mesh,I);
toc

% Graphical representation
figure
plot(meshr)
hold on
plotNrm(meshr,'r')
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')

% Refinement with recursive midpoint algorithm and fixed order
ord = (1 + mesh.vtx(mesh.elt(:,1),1))/2;
ord = floor(4*ord);
tic
meshr = refine(mesh,ord);
toc

% Graphical representation
figure
plot(meshr)
hold on
plotNrm(meshr,'r')
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')

% Refinement with recursive midpoint algorithm and function order
fct = @(X) floor(5*(1+X(:,1))/2);
tic
meshr = refine(mesh,fct);
toc

% Graphical representation
figure
plot(meshr)
hold on
plotNrm(meshr,'r')
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')

% Refinement with recursive midpoint algorithm and fixed edge length
stp    = meshr.stp;
lambda = stp(2)/4;
tic
meshr = refine(meshr,lambda);
toc

% Graphical representation
figure
plot(meshr)
hold on
plotNrm(meshr,'r')
axis equal
xlabel('X'); ylabel('Y'); zlabel('Z')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan