text:
classdef Translation_Test < edu.stanford.covert.cell.sim.ProcessTestCase
    methods
        function this = Translation_Test(name)
            this = this@edu.stanford.covert.cell.sim.ProcessTestCase(name);
        end

        function testSimpleFixture(this)
            m = this.process;

            %% initial conditions: one mRNA, components of one ribosome
            m.stimuli(:) = 0;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_gtp) = 1e6;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.boundEnzymes(:) = 0;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_translationFactors) = 1;
            m.enzymes(m.enzymeIndexs_initiationFactor3) = 0;
            m.enzymes(m.enzymeIndexs_elongationGFactor) = 2;
            m.enzymes(m.enzymeIndexs_ribosome30SIF3) = 1;
            m.enzymes(m.enzymeIndexs_ribosome50S) = 1;
            m.ribosome.states = [];
            m.ribosome.boundMRNAs = [];
            m.ribosome.mRNAPositions = [];
            m.ribosome.tmRNAPositions = [];
            m.polypeptide.boundMRNAs = [];
            m.polypeptide.nascentMonomerLengths = [];
            m.polypeptide.proteolysisTagLengths = [];

            m.polypeptide.monomerAASequences = {'MKKKKKKKKKKKKKKKKK'};
            m.polypeptide.monomerTRNASequences = {[10 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20]'};
            m.polypeptide.monomerLengths = 18;
            m.monomers = 0;

            m.mRNAs = 1;
            m.aminoacylatedTRNAs = repmat(20, size(m.rna.matureTRNAIndexs));
            m.aminoacylatedTMRNA = 0;

            %% initiation
            m.evolveState();
            assertEqual(m.ribosome.activeValue, m.ribosome.states(1));
            assertEqual(0, m.polypeptide.nascentMonomerLengths(1));
            assertEqual(0, m.polypeptide.proteolysisTagLengths(1));
            assertEqual(1, m.polypeptide.boundMRNAs(1));
            assertEqual(1, m.ribosome.boundMRNAs(1));
            assertEqual(0, m.monomers(1));
            assertEqual(zeros(20,1), m.substrates(m.substrateIndexs_aminoAcids));
            assertEqual(1e6-1, m.substrates(m.substrateIndexs_gtp));
            assertEqual(1, m.substrates(m.substrateIndexs_gdp));
            assertEqual(1, m.substrates(m.substrateIndexs_phosphate));
            assertEqual(1e6-1, m.substrates(m.substrateIndexs_water));
            assertEqual(1, m.substrates(m.substrateIndexs_hydrogen));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor1));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor2));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor3));
            assertEqual([2;1;1;1], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual(1, m.enzymes(m.enzymeIndexs_terminationFactor));
            assertEqual(1, m.enzymes(m.enzymeIndexs_recyclingFactor));
            assertEqual(repmat(20, size(m.rna.matureTRNAIndexs)), m.aminoacylatedTRNAs);
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30SIF3));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome50S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome70S));
            assertEqual(1, m.boundEnzymes(m.enzymeIndexs_ribosome70S));

            %% elongation
            m.evolveState();
            assertEqual(m.ribosome.activeValue, m.ribosome.states(1));
            assertEqual(16, m.polypeptide.nascentMonomerLengths(1));
            assertEqual(0, m.polypeptide.proteolysisTagLengths(1));
            assertEqual(1, m.polypeptide.boundMRNAs(1));
            assertEqual(1, m.ribosome.boundMRNAs(1));
            assertEqual(0, m.monomers(1));
            assertEqual(zeros(20,1), m.substrates(m.substrateIndexs_aminoAcids));
            assertEqual(1e6-33, m.substrates(m.substrateIndexs_gtp));
            assertEqual(33, m.substrates(m.substrateIndexs_gdp));
            assertEqual(33, m.substrates(m.substrateIndexs_phosphate));
            assertEqual(1e6-34, m.substrates(m.substrateIndexs_water));
            assertEqual(49, m.substrates(m.substrateIndexs_hydrogen));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor1));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor2));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor3));
            assertEqual([1;0;0;0], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual([1;1;1;1], m.boundEnzymes(m.enzymeIndexs_elongationFactors));
            assertEqual(1, m.enzymes(m.enzymeIndexs_terminationFactor));
            assertEqual(1, m.enzymes(m.enzymeIndexs_recyclingFactor));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30SIF3));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome50S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome70S));
            assertEqual(1, m.boundEnzymes(m.enzymeIndexs_ribosome70S));
            assertEqual([19 20 20 20 20 20 20 20 20 20 5]', m.aminoacylatedTRNAs(10:20));

            %% elongation + termination
            m.evolveState();
            assertEqual(m.ribosome.notExistValue, m.ribosome.states(1));
            assertEqual(0, m.polypeptide.nascentMonomerLengths(1));
            assertEqual(0, m.polypeptide.proteolysisTagLengths(1));
            assertEqual(0, m.polypeptide.boundMRNAs(1));
            assertEqual(0, m.ribosome.boundMRNAs(1));
            assertEqual(1, m.monomers(1));
            assertEqual(zeros(20,1), m.substrates(m.substrateIndexs_aminoAcids));
            assertEqual(1e6-39, m.substrates(m.substrateIndexs_gtp));
            assertEqual(39, m.substrates(m.substrateIndexs_gdp));
            assertEqual(39, m.substrates(m.substrateIndexs_phosphate));
            assertEqual(1e6-40, m.substrates(m.substrateIndexs_water));
            assertEqual(57, m.substrates(m.substrateIndexs_hydrogen));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor1));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor2));
            assertEqual(1, m.enzymes(m.enzymeIndexs_initiationFactor3));
            assertEqual([1;0;0;0], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual([1;1;1;1], m.boundEnzymes(m.enzymeIndexs_elongationFactors));
            assertEqual(1, m.enzymes(m.enzymeIndexs_ribosome30S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30SIF3));
            assertEqual(1, m.enzymes(m.enzymeIndexs_ribosome50S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome70S));
            assertEqual(0, m.boundEnzymes(m.enzymeIndexs_ribosome70S));
            assertEqual([19 20 20 20 20 20 20 20 20 20 3]', m.aminoacylatedTRNAs(10:20));

            %% freeing the elongation factors
            m.evolveState();
            assertEqual([2;1;1;1], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual([0;0;0;0], m.boundEnzymes(m.enzymeIndexs_elongationFactors));
        end

        function testInitiationFromUncomplexedRibosome30SAndIF3(this)
            m = this.process;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_gtp) = 1;
            m.substrates(m.substrateIndexs_water) = 1;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_initiationFactors) = 1;
            m.enzymes(m.enzymeIndexs_ribosome30S) = 1;
            m.enzymes(m.enzymeIndexs_ribosome50S) = 1;
            m.boundEnzymes(:) = 0;
            m.mRNAs(:) = 1;
            m.ribosome.states(:) = m.ribosome.notExistValue;
            m.evolveState();
            assertTrue(any(m.ribosome.states==m.ribosome.activeValue),...
                'initiation failed');
        end

        function testInitiation(this)
            this.initializeProcessForSuccessfulInitiation();
            m = this.process;

            m.evolveState();
            assertTrue(any(m.ribosome.states==m.ribosome.activeValue),...
                'initiation failed');
        end

        function testNoInitiationWithoutAvailableEnergy(this)
            this.initializeProcessForSuccessfulInitiation();
            m = this.process;
            m.substrates(m.substrateIndexs_gtp) = 0;

            m.evolveState();
            assertTrue(~any(m.ribosome.states == m.ribosome.activeValue),...
                'initiation happened');
        end

        function testNoInitiationWithoutMRNAs(this)
            this.initializeProcessForSuccessfulInitiation();
            m = this.process;
            m.mRNAs(:) = 0;

            m.evolveState();
            assertTrue(~any(m.ribosome.states==m.ribosome.activeValue),...
                'initiation happened');
        end

        function testNoInitiationWithOneEnzymeMissing(this)
            m = this.process;
            essential_enzyme_indexs = [...
                m.enzymeIndexs_initiationFactor1,...
                m.enzymeIndexs_initiationFactor2,...
                m.enzymeIndexs_ribosome50S];
            for i=1:numel(essential_enzyme_indexs)
                this.initializeProcessForSuccessfulInitiation();
                m.enzymes(essential_enzyme_indexs(i)) = 0;

                m.evolveState();
                assertTrue(...
                    ~any(m.ribosome.states==m.ribosome.activeValue),...
                    ['initiation happened, i=' num2str(i)]);
            end

            % Ribosome 30S has two forms
            this.initializeProcessForSuccessfulInitiation();
            m.enzymes(m.enzymeIndexs_ribosome30S) = 0;
            m.enzymes(m.enzymeIndexs_ribosome30SIF3) = 0;

            m.evolveState();
            assertTrue(...
                ~any(m.ribosome.states == m.ribosome.activeValue),...
                'initiation happened');
        end

        function testElongation(this)
            this.initializeProcessForSuccessfulElongation(3);
            m = this.process;

            m.evolveState();
            assertEqual(3, m.polypeptide.nascentMonomerLengths);
            assertEqual(1, m.polypeptide.boundMRNAs);
            assertEqual(1, m.ribosome.boundMRNAs);
            assertEqual(0, m.monomers(1));
            assertTrue(~any(m.enzymes(m.enzymeIndexs_elongationFactors)));
            assertTrue(all(m.boundEnzymes(m.enzymeIndexs_elongationFactors)));
        end

        function testNoElongationWithoutEnoughEnergy(this)
            this.initializeProcessForSuccessfulElongation(3);
            m = this.process;
            m.substrates(m.substrateIndexs_gtp) = 1;  % 2 GTP required per AA

            m.evolveState();
            assertEqual(0, m.polypeptide.nascentMonomerLengths);
            assertEqual(1, m.polypeptide.boundMRNAs);
            assertEqual(1, m.ribosome.boundMRNAs);
            assertEqual(0, m.monomers(1));
            assertTrue(all(m.enzymes(m.enzymeIndexs_elongationFactors)));
            assertTrue(~any(m.boundEnzymes(m.enzymeIndexs_elongationFactors)));
        end

        function testNoElongationWithoutWater(this)
            this.initializeProcessForSuccessfulElongation(3);
            m = this.process;
            m.substrates(:) = 1e6;
            m.substrates(m.substrateIndexs_water) = 0;

            m.evolveState();
            assertEqual(0, m.polypeptide.nascentMonomerLengths);
            assertEqual(1, m.polypeptide.boundMRNAs);
            assertEqual(1, m.ribosome.boundMRNAs);
            assertEqual(0, m.monomers(1));
            assertTrue(all(m.enzymes(m.enzymeIndexs_elongationFactors)));
            assertTrue(~any(m.boundEnzymes(m.enzymeIndexs_elongationFactors)));
            assertEqual(0, m.substrates(m.substrateIndexs_water));
        end

        function testNoSubstrates(this)
            this.initializeProcessForSuccessfulElongation(3);

            m = this.process;
            m.substrates(:) = 0;

            initial_substrates = m.substrates;
            initial_enzymes = m.enzymes;
            initial_boundEnzymes = m.boundEnzymes;
            initial_aminoacylatedTRNAs = m.aminoacylatedTRNAs;
            initial_ribosomeStates = m.ribosome.states;
            initial_boundMRNAs = m.ribosome.boundMRNAs;
            initial_monomers = m.monomers;
            initial_monomerLengths = m.polypeptide.monomerLengths;

            m.evolveState();

            %no substrate used
            assertEqual(initial_substrates, m.substrates);
            assertEqual(initial_aminoacylatedTRNAs, m.aminoacylatedTRNAs);

            %all active ribosomes remain active and don't advance
            assertEqual(...
                initial_ribosomeStates(initial_ribosomeStates == m.ribosome.activeValue),...
                m.ribosome.states(initial_ribosomeStates == m.ribosome.activeValue));
            assertEqual(...
                initial_boundMRNAs(initial_ribosomeStates == m.ribosome.activeValue),...
                m.ribosome.boundMRNAs(initial_ribosomeStates == m.ribosome.activeValue));
            assertEqual(...
                initial_boundMRNAs(initial_ribosomeStates == m.ribosome.activeValue),...
                m.polypeptide.boundMRNAs(initial_ribosomeStates == m.ribosome.activeValue));

            %no enzyme gain or loss or binding
            initial_enzymes(m.enzymeIndexs_ribosome30SIF3) = 0;
            initial_enzymes(m.enzymeIndexs_ribosome30S) = 0;
            initial_enzymes(m.enzymeIndexs_initiationFactor3) = 0;

            assertEqual(initial_enzymes', m.enzymes');
            assertEqual(initial_boundEnzymes', m.boundEnzymes');

            %no additional monomers or nascent monomer elongation
            assertEqual(initial_monomers, m.monomers);
            assertEqual(initial_monomerLengths, m.polypeptide.monomerLengths);
        end

        function testManyTranscripts(this)
            m = this.process;

            %% initial conditions: three mRNAs, components of three ribosomes
            m.stimuli(:) = 0;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_gtp) = 1e6;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.boundEnzymes(:) = 0;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_translationFactors) = 6;
            m.enzymes(m.enzymeIndexs_elongationGFactor) = 12;
            m.enzymes(m.enzymeIndexs_ribosome30S) = 6;
            m.enzymes(m.enzymeIndexs_ribosome50S) = 6;
            m.ribosome.states = zeros(0, 1);
            m.ribosome.boundMRNAs = zeros(0, 1);
            m.ribosome.mRNAPositions = zeros(0, 1);
            m.ribosome.tmRNAPositions = zeros(0, 1);
            m.polypeptide.boundMRNAs = zeros(0, 1);
            m.polypeptide.nascentMonomerLengths = zeros(0, 1);
            m.polypeptide.proteolysisTagLengths = zeros(0, 1);

            m.polypeptide.monomerAASequences = {
                'MMKRTYQPSKLKRAKTH';
                'MGFMARMATAQGRKVLR';
                'MQRRFKNRAQLTVSSER'};
            m.polypeptide.monomerTRNASequences = {
                [10 7 20 17 28 22 21 6 3 20 19 27 13 1 20 4 33]';
                [10 14 12 7 1 15 7 1 4 1 21 18 13 20 29 19 15]';
                [10 21 36 13 12 20 32 13 1 21 34 30 29 24 3 31 13]'};

            m.polypeptide.monomerLengths = [17 17 17]';
            m.monomers = [0 0 0]';

            m.mRNAs = [1 1 1]';
            m.aminoacylatedTRNAs = repmat(100, size(m.rna.matureTRNAIndexs));
            m.aminoacylatedTMRNA = 0;

            % 3 initiated, 3 idle
            m.evolveState();
            assertEqual(3, numel(find(m.ribosome.states==m.ribosome.activeValue)));
            assertTrue(~any(m.polypeptide.nascentMonomerLengths));
            assertEqual(0, sum(m.monomers));

            % 3 elongating, 3 initiated
            m.evolveState();
            assertEqual(6, numel(find(m.ribosome.states==m.ribosome.activeValue)));
            assertEqual([16 16 16 0 0 0]', m.polypeptide.nascentMonomerLengths);
            assertEqual(0, sum(m.monomers));

            % 3 terminated, 3 elongating
            m.evolveState();
            assertEqual(3, numel(find(m.ribosome.states==m.ribosome.activeValue)));
            assertEqual([0 0 0 16 16 16]', m.polypeptide.nascentMonomerLengths);
            assertEqual(3, sum(m.monomers));

            % 3 initiated, 3 terminated
            m.evolveState();
            assertEqual(3, numel(find(m.ribosome.states==m.ribosome.activeValue)));
            assertEqual([0 0 0 0 0 0]', m.polypeptide.nascentMonomerLengths);
            assertEqual(6, sum(m.monomers));

            % 3 elongating, 3 initiated
            m.evolveState();
            assertEqual(6, numel(find(m.ribosome.states==m.ribosome.activeValue)));
            assertEqual([16 16 16 0 0 0]', m.polypeptide.nascentMonomerLengths);
            assertEqual(6, sum(m.monomers));

            assertEqual(6, numel(find(m.ribosome.boundMRNAs)));
            assertEqual(6, numel(find(m.polypeptide.boundMRNAs)));
            assertEqual(100 * numel(m.rna.matureTRNAIndexs) - 150, sum(m.aminoacylatedTRNAs));

            % substrate usage = initiation + elongation + termination
            assertEqual(1e6-(12 + 2 * 150 + 2 * 6), m.substrates(m.substrateIndexs_gtp));
            assertEqual(12 + 2 * 150 + 2 * 6, m.substrates(m.substrateIndexs_gdp));
            assertEqual(12 + 2 * 150 + 2 * 6, m.substrates(m.substrateIndexs_phosphate));
            assertEqual(1e6-(12 + 2*150 + 3 * 6 + 3), m.substrates(m.substrateIndexs_water));
            assertEqual(12 + 3 * 150 + 2 * 6, m.substrates(m.substrateIndexs_hydrogen));

            assertEqual([6;6;6], m.enzymes(m.enzymeIndexs_initiationFactors));
            assertEqual([9;3;3;3], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual([3;3;3;3], m.boundEnzymes(m.enzymeIndexs_elongationFactors));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30SIF3));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome50S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome70S));
            assertEqual(6, m.boundEnzymes(m.enzymeIndexs_ribosome70S));
        end

        function testBindingProbabilities(this)
            m = this.process;

            %plentiful substrate and enzyme and a small number of ribosomes
            numRibosomes = 50;
            m.substrates(:) = 1e6;
            m.enzymes(:) = 1e6;
            m.enzymes(m.enzymeIndexs_ribosome30S) = 0;
            m.enzymes(m.enzymeIndexs_ribosome30SIF3) = numRibosomes;
            m.enzymes(m.enzymeIndexs_ribosome50S) = numRibosomes;
            m.enzymes(m.enzymeIndexs_ribosome70S) = 0;

            %no bound ribosomes
            m.boundEnzymes(:) = 0;
            m.ribosome.states = repmat(m.ribosome.notExistValue, numRibosomes, 1);
            m.ribosome.boundMRNAs = zeros(numRibosomes, 1);
            m.polypeptide.boundMRNAs = zeros(numRibosomes, 1);
            m.polypeptide.nascentMonomerLengths = zeros(numRibosomes, 1);
            m.polypeptide.proteolysisTagLengths = zeros(numRibosomes, 1);
            m.ribosome.mRNAPositions = m.polypeptide.nascentMonomerLengths;
            m.ribosome.tmRNAPositions = m.polypeptide.proteolysisTagLengths;

            %binding probabilities
            m.mRNAs(1:numRibosomes) = 1;
            m.mRNAs(numRibosomes+1:end) = 0;

            m.evolveState();
            assertTrue(all(m.ribosome.states == m.ribosome.activeValue), 'not all ribosomes bound');
            assertEqual(1:numRibosomes, unique(m.ribosome.boundMRNAs)', 'not all mRNAs bound');
            assertEqual(1:numRibosomes, unique(m.polypeptide.boundMRNAs)', 'not all mRNAs bound');
        end

        function testRibosomeStalling(this)
            m = this.process;

            %% initial conditions: one ribosome bound to one mRNA
            m.stimuli(:) = 0;
            m.substrates(:) = 0;
            m.substrates(m.substrateIndexs_gtp) = 1e6;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_translationFactors) = 1;
            m.enzymes(m.enzymeIndexs_elongationGFactor) = 2;
            m.enzymes(m.enzymeIndexs_tmRNABindingProtein) = 1;
            m.enzymes(m.enzymeIndexs_peptidylTRNAHydrolase) = 1;
            m.boundEnzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ribosome70S) = 1;
            m.boundTMRNA = 0;

            m.polypeptide.monomerAASequences = {'MKKKKKVVVVV'};
            m.polypeptide.monomerTRNASequences = {[10 20 20 20 20 20 29 29 29 29 29]'};
            m.polypeptide.monomerLengths = 11;
            m.monomers = 0;

            m.mRNAs = 1;
            m.ribosome.boundMRNAs = 1;
            m.polypeptide.boundMRNAs = 1;
            m.ribosome.states = m.ribosome.activeValue;
            m.ribosome.mRNAPositions = 0;
            m.ribosome.tmRNAPositions = 0;
            m.polypeptide.nascentMonomerLengths = 0;
            m.polypeptide.proteolysisTagLengths = 0;
            m.aminoacylatedTRNAs = repmat(20, size(m.rna.matureTRNAIndexs));
            m.aminoacylatedTRNAs(29) = 0;  % withhold valine
            m.aminoacylatedTMRNA = 1;
            m.tmRNABindingProbability = 1;

            %% elongation
            m.evolveState();
            assertEqual(m.ribosome.activeValue, m.ribosome.states(1));
            assertEqual(6, m.polypeptide.nascentMonomerLengths(1));
            assertEqual(0, m.polypeptide.proteolysisTagLengths(1));
            assertEqual(1, m.polypeptide.boundMRNAs(1));
            assertEqual(1, m.ribosome.boundMRNAs(1));
            assertEqual(0, m.monomers(1));

            %% stalling
            m.evolveState();
            assertEqual(m.ribosome.stalledValue, m.ribosome.states(1));
            assertEqual(6, m.polypeptide.nascentMonomerLengths(1));
            assertEqual(1, m.polypeptide.proteolysisTagLengths(1));
            assertEqual(1, m.polypeptide.boundMRNAs(1));
            assertEqual(1, m.ribosome.boundMRNAs(1));
            assertEqual(0, m.monomers(1));
            assertEqual(0, m.aminoacylatedTMRNA);
            assertEqual(0, m.enzymes(m.enzymeIndexs_tmRNA));
            assertEqual(1, m.boundTMRNA);
            assertEqual(zeros(20,1), m.substrates(m.substrateIndexs_aminoAcids));
            assertEqual(1e6 - 12, m.substrates(m.substrateIndexs_gtp));
            assertEqual(12, m.substrates(m.substrateIndexs_gdp));
            assertEqual(12, m.substrates(m.substrateIndexs_phosphate));
            assertEqual(1e6 - 12 + 6 -6 -1, m.substrates(m.substrateIndexs_water));
            assertEqual(12+6, m.substrates(m.substrateIndexs_hydrogen));

            %% elongation and termination of proteolysis tag
            m.aminoacylatedTRNAs(29) = 20;  % provide valine again
            m.evolveState();
            m.evolveState();
            assertEqual(1, numel(m.polypeptide.abortedSequences));
            assertEqual([1;0;0;0], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual([1;1;1;1], m.boundEnzymes(m.enzymeIndexs_elongationFactors));
            assertEqual(1, m.enzymes(m.enzymeIndexs_ribosome30S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome30SIF3));
            assertEqual(1, m.enzymes(m.enzymeIndexs_ribosome50S));
            assertEqual(0, m.enzymes(m.enzymeIndexs_ribosome70S));
            assertEqual(0, m.boundEnzymes(m.enzymeIndexs_ribosome70S));
            assertEqual([19 18 18 20 20 20 20 20 20 20 14], m.aminoacylatedTRNAs(10:20)');

            %% freeing the elongation factors
            m.evolveState();
            assertEqual([2;1;1;1], m.enzymes(m.enzymeIndexs_elongationFactors));
            assertEqual([0;0;0;0], m.boundEnzymes(m.enzymeIndexs_elongationFactors));
        end

        function testDetailedMassBalance(this)
            m = this.process;

            m.mRNAs(:) = 0;
            m.mRNAs(1) = 1;

            m.enzymes = m.enzymes + m.boundEnzymes;
            m.boundEnzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_ribosome30S) = 1;
            m.enzymes(m.enzymeIndexs_ribosome30SIF3) = 0;
            m.enzymes(m.enzymeIndexs_ribosome50S) = 1;
            m.enzymes(m.enzymeIndexs_ribosome70S) = 0;

            m.ribosome.boundMRNAs = 0;
            m.polypeptide.boundMRNAs = 0;
            m.ribosome.states = m.ribosome.notExistValue;
            m.ribosome.mRNAPositions = 0;
            m.ribosome.tmRNAPositions = 0;
            m.polypeptide.nascentMonomerLengths = 0;
            m.polypeptide.proteolysisTagLengths = 0;

            m.substrates(:) = 1e6;
            m.aminoacylatedTRNAs(:) = 1e6;

            initSubstrates = m.substrates;

            for i = 1:200;
                m.evolveState();
                if m.monomers(1) == 1
                    break;
                end
            end

            finalSubstrates = initSubstrates;
            gtpCost = 3 + 2 * m.monomer.lengths(1);
            finalSubstrates(m.substrateIndexs_gtp) = finalSubstrates(m.substrateIndexs_gtp) - gtpCost;
            finalSubstrates(m.substrateIndexs_gdp) = finalSubstrates(m.substrateIndexs_gdp) + gtpCost;
            finalSubstrates(m.substrateIndexs_phosphate) = finalSubstrates(m.substrateIndexs_phosphate) + gtpCost;
            finalSubstrates(m.substrateIndexs_hydrogen) = finalSubstrates(m.substrateIndexs_hydrogen) ...
                + gtpCost ...           %energy
                + m.monomer.lengths(1); %peptidyl-tRNA bond
            finalSubstrates(m.substrateIndexs_water) = finalSubstrates(m.substrateIndexs_water) ...
                - gtpCost ...               %energy
                - m.monomer.lengths(1) ...  %peptidyl-tRNA bond
                + (m.monomer.lengths(1)-1); %peptide bond

            assertEqual(1, m.monomers(1));
            assertEqual(0, m.ribosome.boundMRNAs);
            assertEqual(finalSubstrates(setdiff(1:end, m.substrateIndexs_aminoAcids))', ...
                m.substrates(setdiff(1:end, m.substrateIndexs_aminoAcids))');
            assertEqual(finalSubstrates, m.substrates);
        end

        function testNoBiasedPausing(this)
            m = this.process;
            rib = m.ribosome;
            pol = m.polypeptide;

            %replace tRNA sequences with codons corresponding to limiting
            %tRNA
            pol.monomerTRNASequences = cellfun(@(seq) ones(1, 200), pol.monomerTRNASequences, 'UniformOutput', false);

            %all ribosomes
            %- active
            %- translating first mRNA at first codon which we've set to
            %  correspond to the limiting tRNA
            nRibs = 10;
            rib.states = repmat(rib.activeValue, [nRibs 1]);
            rib.boundMRNAs = ones(nRibs, 1);
            rib.mRNAPositions = zeros(nRibs, 1);
            rib.tmRNAPositions = zeros(nRibs, 1);

            %synchronize polypepide and ribosome states
            pol.boundMRNAs = rib.boundMRNAs;
            pol.nascentMonomerLengths = rib.mRNAPositions;
            pol.proteolysisTagLengths = rib.tmRNAPositions;
            pol.abortedPolypeptides = zeros(0, 3);

            %synchronize ribosome and protein complex states
            m.enzymes(:) = 0;
            m.boundEnzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ribosome70S) = nRibs;
            m.enzymes(m.enzymeIndexs_translationFactors) = 100;
            m.boundEnzymes(m.enzymeIndexs_translationFactors) = 100;

            %simulate
            iterMax = 100;
            for i = 1:iterMax
                %lot of tRNA, but 1 tRNA limiting
                m.aminoacylatedTRNAs(:) = 100;
                m.aminoacylatedTRNAs(1) = 1;
                m.freeTRNAs(:) = 0;

                %simulate
                m.evolveState();
            end

            %assert that all ribosomes move
            assertEqual(nRibs, nnz(rib.mRNAPositions));
            assertIn(range(rib.mRNAPositions), [0 20]);
        end

        function testGeneEssentiality(this)
            m = this.process;

            m.substrates(:) = 1e6;
            m.enzymes(:) = 1e6;
            numRibosomes = 10;
            m.enzymes(m.enzymeIndexs_ribosome30SIF3) = 0;
            m.enzymes(m.enzymeIndexs_ribosome30S) = numRibosomes;
            m.enzymes(m.enzymeIndexs_ribosome50S) = numRibosomes;
            m.enzymes(m.enzymeIndexs_ribosome70S) = 0;
            m.boundEnzymes(m.enzymeIndexs_ribosome70S) = 0;
            m.aminoacylatedTRNAs = repmat(1e6, size(m.rna.matureTRNAIndexs));
            m.aminoacylatedTMRNA = 1e6;
            m.mRNAs(:) = 1;
            m.monomers(:) = 0;

            m.ribosome.states = zeros(0, 1);
            m.ribosome.boundMRNAs = zeros(0, 1);
            m.ribosome.mRNAPositions = zeros(0, 1);
            m.ribosome.tmRNAPositions = zeros(0, 1);
            m.polypeptide.boundMRNAs = zeros(0, 1);
            m.polypeptide.nascentMonomerLengths = zeros(0, 1);
            m.polypeptide.proteolysisTagLengths = zeros(0, 1);

            this.helpTestGeneEssentiality({
                'MG_173';   %translation initiation factor IF-1
                'MG_142';   %translation initiation factor IF-2
                'MG_196';   %translation initiation factor IF-3
                'MG_089';   %translation elongation factor G
                'MG_026';   %translation elongation factor P
                'MG_451';   %translation elongation factor Tu
                'MG_433';   %translation elongation factor Ts
                'MG_435';   %ribosome recycling factor
                'MG_473';   %ribosomal protein L33 type 2
                'MG_070';   %ribosomal protein S2
                'MG_081';   %ribosomal protein L11
                'MG_082';   %ribosomal protein L1
                'MG_087';   %ribosomal protein S12
                'MG_088';   %ribosomal protein S7
                'MG_090';   %ribosomal protein S6
                'MG_092';   %ribosomal protein S18
                'MG_093';   %ribosomal protein L9
                'MGrrnA5S';   %5S ribosomal rRNA
                'MGrrnA16S';  %16S ribosomal rRNA
                'MGrrnA23S';  %23S ribosomal rRNA
                'MG_150';   %ribosomal protein S10
                'MG_151';   %ribosomal protein L3
                'MG_152';   %ribosomal protein L4/L1 family
                'MG_153';   %ribosomal protein L23
                'MG_154';   %ribosomal protein L2
                'MG_155';   %ribosomal protein S19
                'MG_156';   %ribosomal protein L22
                'MG_157';   %ribosomal protein S3
                'MG_158';   %ribosomal protein L16
                'MG_159';   %ribosomal protein L29
                'MG_160';   %ribosomal protein S17
                'MG_161';   %ribosomal protein L14
                'MG_162';   %ribosomal protein L24
                'MG_163';   %ribosomal protein L5
                'MG_164';   %ribosomal protein S14
                'MG_165';   %ribosomal protein S8
                'MG_166';   %ribosomal protein L6
                'MG_167';   %ribosomal protein L18
                'MG_168';   %ribosomal protein S5
                'MG_169';   %ribosomal protein L15
                'MG_174';   %ribosomal protein L36
                'MG_175';   %ribosomal protein S13
                'MG_176';   %ribosomal protein S11
                'MG_178';   %ribosomal protein L17
                'MG_197';   %ribosomal protein L35
                'MG_198';   %ribosomal protein L20
                'MG_481';   %30S ribosomal protein S21
                'MG_232';   %ribosomal protein L21
                'MG_234';   %ribosomal protein L27
                'MG_257';   %ribosomal protein L31
                'MG_258';   %peptide chain release factor 1
                'MG_311';   %ribosomal protein S4
                'MG_325';   %ribosomal protein L33
                'MG_361';   %ribosomal protein L10
                'MG_362';   %ribosomal protein L7/L12
                'MG_363';   %ribosomal protein L32
                'MG_522';   %ribosomal protein S20
                'MG_417';   %ribosomal protein S9
                'MG_418';   %ribosomal protein L13
                'MG_424';   %ribosomal protein S15
                'MG_426';   %ribosomal protein L28
                'MG_444';   %ribosomal protein L19
                'MG_446';   %ribosomal protein S16
                'MG_466';   %ribosomal protein L34
                },...
                @(m,i) any(m.monomers > i.monomers),...
                struct('lengthSec', 75));
        end
    end

    methods
        function initializeProcessForSuccessfulInitiation(this)
            m = this.process;
            m.substrates(m.substrateIndexs_gtp) = 1e6;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.enzymes(:) = 10;
            m.boundEnzymes(:) = 0;
            m.aminoacylatedTRNAs = repmat(20, size(m.rna.matureTRNAIndexs));
            m.mRNAs(:) = 1;
            m.ribosome.states(:) = m.ribosome.notExistValue;
        end

        function initializeProcessForSuccessfulElongation(this, numAAs)
            m = this.process;
            m.substrates(m.substrateIndexs_gtp) = 2 * numAAs;
            m.substrates(m.substrateIndexs_water) = 1e6;
            m.mRNAs(:) = 1;
            m.enzymes(:) = 0;
            m.enzymes(m.enzymeIndexs_elongationFactors) = 1;
            m.boundEnzymes(:) = 0;
            m.boundEnzymes(m.enzymeIndexs_ribosome70S) = 1;
            m.ribosome.boundMRNAs = 1;
            m.polypeptide.boundMRNAs = 1;
            m.ribosome.states = m.ribosome.activeValue;
            m.ribosome.mRNAPositions = 0;
            m.ribosome.tmRNAPositions = 0;
            m.polypeptide.nascentMonomerLengths = 0;
            m.polypeptide.proteolysisTagLengths = 0;
            m.aminoacylatedTRNAs(:) = numAAs;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%{
Authors:
Jonatas Lopes de Paiva
Claudio Fabiano Motta Toledo
Helio Pedrini

%}

function pop = createPop(sizePop, img, denoisedImages, beta, lambda  )

population(sizePop).fitness = 0;
population(sizePop).cromo = [];
population(sizePop).id = '';

mut1 = denoisedImages(1).img;
mut2 = denoisedImages(2).img;
mut3 = denoisedImages(3).img;

fit1 = calcFitness(mut1, img, beta, lambda);
fit2 = calcFitness(mut2, img, beta, lambda);
fit3 = calcFitness(mut3, img, beta, lambda);

population(1).fitness = fit1;
population(1).cromo = mut1;
population(1).id = char(java.util.UUID.randomUUID.toString);

population(2).fitness = fit2;
population(2).cromo = mut2;
population(2).id = char(java.util.UUID.randomUUID.toString);

population(3).fitness = fit3;
population(3).cromo = mut3;
population(3).id = char(java.util.UUID.randomUUID.toString);

for i = 4:sizePop
    tmp = denoisedImages(randi(3)).img;
    population(i).cromo = mutation(tmp);
    population(i).fitness = calcFitness(population(i).cromo, img, beta, lambda);
    population(i).id = char(java.util.UUID.randomUUID.toString);
end

pop = arrangePop(population);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function pass = test_univariate(~)
%TEST_UNIVARIATE   Test various univariate TREEVAR methods.
%   This methods generates a number of anonymous functions that we convert to
%   first order format. We then evaluate them, and check that they match the
%   expected results.

%% Setup
% Initialize a TREEVAR variable.
dom = [0 2];

% Used for setting up anonymous functions we convert to first order form.
rhs = 3.2;
alp = 2;

% Arguments we eventually evaluate the functions at:
tArg = 1;
uArg = [.4 .2];

% List of methods that we want to test:
testMethods = {@abs, @acos, @acosd, @acot, @acoth, @acsc, @acscd, @acsch, ...
    @airy, @asec, @asecd, @asech, @asin, @asind, @asinh, @atan, @atand, ...
    @atanh, @conj, @cos, @cosd, @cosh, @cot, @cotd, @coth, @csc, @cscd, ...
    @csch, @exp, @expm1, @imag, @log, @log10, @log1p, @log2, @pow2, @real, ...
    @sec, @secd, @sech, @sin, @sind, @sinh, @sqrt, @tan, @tand, @tanh, ...
    @uminus, @uplus};

% Store comparison errors:
errors = zeros(length(testMethods), 1);
for mCounter = 1:length(testMethods)
    % Current method we're testing:
    method = testMethods{mCounter};

    % Construct an anonymous function that calls the current method of interest:
    myFun = @(u) diff(u, 2) + diff(u) +  alp*method(u);

    % Convert MYFUN to first order system:
    anonFun = treeVar.toFirstOrder(myFun, rhs, dom);

    % The correct first order reformulation of MYFUN:
    correctFun = @(t,u) [u(2); rhs - u(2) - alp*method(u(1))];

    % Compare the result of evaluating the automatically converted ANONFUN and
    % the manually constructed CORRECTFUN. The difference in the outputs should
    % be very small (if any):
    errors(mCounter) = norm(anonFun(tArg, uArg) - correctFun(tArg, uArg));
end

pass = errors < 10*eps;

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pass = test_volt(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

K = @(u, v) exp(-(u-v).^2);     % A simple kernel.
f = chebfun(@sin, pref);        % A simple chebfun.
F = volt(K, f);                 % Call volt().

% Test against some V4 results:
pass(1) = abs(F(.5) - (-0.013808570536509)) < 10*vscale(F)*eps;
pass(2) = abs(norm(F) - 0.334612395278957) < 10*vscale(F)*eps;

% Test 3rd input argument. (Simply make sure we don't crash!)
try
    F2 = volt(K, f, 1);
    pass(3) = 1;
catch
    pass(3) = 0;
end

end


text:
%GENDATW Sample dataset by given weigths
%
%   B = GENDATW(A,V,N)
%
% INPUT
%   A    Dataset
%   V    Vector with weigths for each object in A
%   N    Number of objects to be generated (default size A);
%
% OUTPUT
%   B    Dataset
%
% DESCRIPTION
% The dataset A is sampled using the weigths in V as a prior distribution.

function b = gendatw(a,v,n)

isdataset(a);
if nargin < 3, n  = size(a,1); end

v = v./sum(v);
if any(v<0)
    error('Weights should be positive');
end

mins = 0;
nn = 0;
while(mins < 2)
	N = genclass(n,v);
	L = [];
	while any(N > 0)
	    L = [L find(N > 0)];
	    N = N-1;
	end
	b = a(L,:);
	mins = min(classsizes(b));
	nn = nn + 1;
	if nn > 100
		error('Problems with weighted subsampling: classes have disappeared. Please enlarge training set.')
	end
end

return

author: marianux

repo: ecg-kit

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


text:
function out = comp_loss(E,loss)

switch loss
    case 'l1'
        out = norm(E(:),1);
    case 'l21'
        out = 0;
        for i = 1 : size(E,2)
            out = out + norm(E(:,i));
        end
    case 'l2'
        out = 0.5*norm(E,'fro')^2;
end

author: canyilu

repo: LibADMM-toolbox


text:
function [Fn,Ln] = fibonacci(n,modulus)
% fibonacci: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
% usage: [Fn,Ln] = fibonacci(n)
% usage: [Fn,Ln] = fibonacci(n,modulus)
%
% Compute the nth Fibonacci number as well as the nth Lucas
% Lucas number. In the event that all members of these
% sequences from 1 to n are desired, then a simple,
% direct loop would be more efficient, and more direct.
%
% Both the Fibonacci numbers and the Lucas numbers
% are defined by the same basic recursion:
%
%  F(n) = F(n-1) + F(n-2)
%  L(n) = L(n-1) + L(n-2)
%
% The difference is the starting point. The Fibonacci
% numbers start with F(1) = F(2) = 1, whereas the Lucas
% sequence starts with L(1) = 1, and L(2) = 3. The first
% few members of these sequences are:
%
%  Fibonacci: [1 1 2 3 5 8 13 21 ... ]
%  Lucas:     [1 3 4 7 11 18 29 ... ]
%
% These sequences are also defined for n = 0 and for
% negative values of n.
%
% For efficiency, fibonacci uses a variety of tricks to
% maximize speed. While computation of fibonacci numbers
% is commonly done recursively, fibonacci does so using a
% direct iterative scheme given the binary representation
% of n. In addition, several Fibonacci and Lucas number
% identities are employed to maximize throughput.
%
% The methods employed by fibonacci will be O(log2(n)) in time.
%
%
% Arguments: (input)
%  n - any non-negative integer, vpi or numeric.
%      n must be less than 2^53 (in theory. Even that
%      number would be impossibly huge to compute.
%      Don't bother to try it.) In practice, recognize
%      that F(1e6) is a number with 208995 digits.
%
%      When n is a vector or array, fibonacci will
%      generate each of the indicated Fibonacci and
%      Lucas numbers.
%
%  modulus - (OPTIONAL) - allows the computation of the
%      indicated Fibonacci/Lucas numbers modulo a given
%      modulus. This enables the computation of such
%      numbers for truly immense index.
%
% Arguments: (output)
%  Fn, Ln - scalar vpi number, containing the nth
%      Fibonacci number and nth Lucas numbers in
%      their respective sequences.
%
% Example:
%  [Fn,Ln] = fibonacci(150)
%
% Fn =
%    9969216677189303386214405760200
% Ln =
%    22291846172619859445381409012498
%
%  See also:
%
%  Author: John D'Errico
%  e-mail: woodchips@rochester.rr.com
%  Release: 1.0
%  Release date: 4/30/09

if (nargin < 1) || (nargin > 2)
  error('fibonacci accepts only 1 or 2 arguments')
elseif any(n(:)~=round(n(:))) || any(abs(n) > 2^53)
  error('n must be an integer, <= 2^53 in absolute value')
end

% was a modulus provided?
if (nargin < 2)
  modulus = [];
end

% The first 15 Fibonacci and Lucas numbers to
% start things off efficiently.
Fseq = [0 1 1 2 3  5  8 13 21 34  55  89 144 233 377  610];
Lseq = [2 1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364];

if ~isempty(modulus)
  Fseq = mod(Fseq,modulus);
  Lseq = mod(Lseq,modulus);
end

% intialize Fn and Ln to the proper size,
% in case n is a vector or array
Fn = repmat(vpi(0),size(n));
Ln = Fn;

% much faster if n is not a vpi but a double.
% also, we don't need to worry about n being
% larger than 2^53, as this would have a vast
% number of digits.
n = double(n(:));

% catch any zeros in n first.
k = (n(:) == 0);
% Fn(k) is already zero
if any(k)
  Ln(k) = vpi(2);
end

% Negative values for n will be inconvenient
% in a loop, so make them all positive. deal
% with any signs later.
nsign = 2*(n>=0) - 1;
if any(n<0)
  neven = iseven(n);
  n = abs(n);
end

if numel(n) > 1
  % sort them and work upwards
  [n,ntags] = sort(n);

  flag = false;
  for i = 1:numel(n)
    if n(i) <= 15
      % small values of n are pre-computed
      Fn(i) = Fseq(n(i) + 1);
      Ln(i) = Lseq(n(i) + 1);
    elseif (i == 1) || ((n(i) - n(i-1)) > 15)
      % the smallest value of n is > 15,
      % or the difference from the last value
      % computed was too large
      [Fn(i),Ln(i)] = fibonacci(n(i),modulus);

      flag = false;
    elseif n(i) == n(i-1)
      % will happen if there were positive and
      % negative values in the list
      Fn(i) = Fn(i-1);
      Ln(i) = Ln(i-1);

      flag = false;
    elseif (n(i) == (n(i-1)+1))
      % we can use an addition formula to
      % get Fn & Ln. Which one?
      if flag
        % the last two values of n were
        % separated by 1, so just use the
        % two term recurrence
        Fn(i) = (Fn(i-1) + Fn(i-2));
        Ln(i) = (Ln(i-1) + Ln(i-2));
        if ~isempty(modulus)
          Fn(i) = mod(Fn(i),modulus);
          Ln(i) = mod(Ln(i),modulus);
        end
      else
        % we must use the addition formula
        Fn(i) = (Fn(i-1) + Ln(i-1))./2;
        Ln(i) = (5 .*Fn(i-1) + Ln(i-1))./2;
        if ~isempty(modulus)
          Fn(i) = mod(Fn(i),modulus);
          Ln(i) = mod(Ln(i),modulus);
        end

        % flag indicates whether the last two
        % numbers were consecutive
        flag = true;
      end

    else
      % 1 < n(i) - n(i-1) <= 15
      m = n(i) - n(i-1);
      Fm = Fseq(m+1);
      Lm = Lseq(m+1);

      % use an addition formula
      Fn(i) = (Fm.*Ln(i-1) + Lm.*Fn(i-1))./2;
      Ln(i) = (5 .*Fm.*Fn(i-1) + Lm.*Ln(i-1))./2;
      if ~isempty(modulus)
        Fn(i) = mod(Fn(i),modulus);
        Ln(i) = mod(Ln(i),modulus);
      end

      flag = false;
    end

  end % for i = 1:numel(n)

  % shuffle Fn and Ln to reflect the sort
  Fn(ntags) = Fn;
  Ln(ntags) = Ln;

else
  % For only one value of n, compute it individually.
  % Uses an efficient iterative (recursive, but not
  % really so) scheme.

  % get the binary representation of n. Thus
  % nbin is a character vector, of length
  % ceil(log2(n)).
  nbin = dec2bin(n);

  % get the 4 highest order bits from nbin
  k = min(numel(nbin),4);

  % start the sequence from the top
  % few bits of n.
  nhigh = bin2dec(nbin(1:k));

  Fn = vpi(Fseq(nhigh+1));
  Ln = vpi(Lseq(nhigh+1));

  % We need to loop forwards. Essentially, we started
  % with the highest order bit(s) of the binary representation
  % for n. Look at each successively lower order bit.
  % If the next bit is 0, then we are essentially doubling
  % the index at this step. If the next bit is odd, then
  % we are moving to 2*n+1.
  for k = 5:numel(nbin)
    bit = (nbin(k) == '1');
    if bit
      % the next bit was odd. Use
      % a 2*n+1 rule to step up.
      F2n = Fn.*Ln;
      % we want to do this...
      %   L2n = (5 .*Fn.*Fn + Ln.*Ln)./2;
      % Instead use the identity that
      %   5F(n)^2 + L(n)^2 = 2*L(n)^2 + 4*(-1)^(n+1)
      % to make that expression more efficiently
      % computed. See that the form used below
      % has only a single multiplication between a
      % pair of large integers, whereas the prior
      % form for L2n would have had several multiples
      % as well as a divide.
      L2n = Ln.*Ln + 2*(-1)^(nhigh+1);

      Fn = (L2n + F2n)./2;
      Ln = (5 .*F2n + L2n)./2;
      if ~isempty(modulus)
        Fn = mod(Fn,modulus);
        Ln = mod(Ln,modulus);
      end
    else
      % the next bit was even. Use the 2*n
      % rule to step up.
      F2n = Fn.*Ln;
      Ln = Ln.*Ln + 2*(-1)^(nhigh+1);
      Fn = F2n;
      if ~isempty(modulus)
        Fn = mod(Fn,modulus);
        Ln = mod(Ln,modulus);
      end
    end

    % update the top bits of n
    nhigh = 2*nhigh + bit;
  end % for k = 5:numel(nbin)

end % if numel(n) > 1

% if n was negative, then we may need to apply a
% sign change to Fn and Ln.
if any(nsign < 0)
  k = neven & (nsign < 0);
  Fn(k) = -Fn(k);

  k = (~neven) & (nsign < 0);
  Ln(k) = -Ln(k);

  if ~isempty(modulus)
    Fn(k) = mod(Fn(k),modulus);
    Ln(k) = mod(Ln(k),modulus);
  end
end

% ==================================
% End mainline, begin subfunctions.
% ==================================

function result = iseven(n)
% tests if a scalar value is an even integer, works
% for either numeric or vpi inputs
if isnumeric(n)
  result = (mod(n,2) == 0);
elseif isa(n,'vpi')
  % must have been a vpi
  result = (mod(trailingdigit(n,1),2) == 0);
else
  error('n must be either numeric or vpi')
end

author: Sable

repo: mcbench-benchmarks


text:
function X=JumpDiffusionKou(m,s,l,p,e1,e2,ts,J)

T=ts(end);
% simulate number of jumps
N=poissrnd(l*T,J,1);

Jumps=[];
L=length(ts);
for j=1:J
    % simulate jump arrival time
    t=T*rand(N(j),1);
    t=sort(t);

    % simulate jump size
    ww=binornd(1,p,N(j),1);
    S=ww.*exprnd(e1,N(j),1)-(1-ww).*exprnd(e2,N(j),1);

    % put things together
    CumS=cumsum(S);
    for n=1:L
        Events=sum(t<=ts(n));
        Jumps_ts(n)=0;
        if Events
            Jumps_ts(n)=CumS(Events);
        end
    end

    Jumps=[Jumps
        Jumps_ts];
end
for l=1:L
    Dt=ts(l);
    if l>1
        Dt=ts(l)-ts(l-1);
    end
    D_Diff(:,l)=m*Dt + s*sqrt(Dt)*randn(J,1);
end

X=[zeros(J,1) cumsum(D_Diff,2)+Jumps];


text:
function test_issue1564

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_artifact_clip ft_artifact_threshold ft_artifact_jump ft_rejectartifact

global ft_default
ft_default.representation = 'table';

%%

ntrials = 10;
nchans = 3;
fsample = 1000;

data = [];
data.label = arrayfun(@num2str, 1:nchans, 'UniformOutput', false);
for i=1:ntrials
  data.time{i} = (1:fsample)/fsample;
  data.trial{i} = randn(nchans, fsample);
  data.sampleinfo(i,1) = (i-1)*fsample + 1 + (i*100); % create some gaps between the trials
  data.sampleinfo(i,2) = data.sampleinfo(i,1) + fsample - 1;
end

if strcmp(ft_default.representation, 'numeric')
  data.trialinfo = [1:ntrials]';
else
  trialnum = arrayfun(@num2str, 1:ntrials, 'UniformOutput', false)';
  data.trialinfo = table(trialnum);
end

%%

cfg = [];
cfg.continuous = 'no';
ft_databrowser(cfg, data);

%%

cfg = [];
cfg.continuous = 'yes';
ft_databrowser(cfg, data);

%%

cfg = [];
cfg.continuous = [];
ft_databrowser(cfg, data);

%%
% close all

data_clip = data;
data_clip.trial{1}(1,201:700) =  2;
data_clip.trial{2}(2,301:400) = -2;
data_clip.trial{2}(2,501:600) = -2;

cfg = [];
cfg.artfctdef.clip.amplthreshold = 0;
[cfg, artifact] = ft_artifact_clip(cfg, data_clip);

% cfg = [];
% cfg.artfctdef.clip.amplthreshold = '0%';
% [cfg, artifact] = ft_artifact_clip(cfg, data_clip)

% cfg = [];
% cfg.artfctdef.clip.amplthreshold = '100%';
% [cfg, artifact] = ft_artifact_clip(cfg, data_clip)

% explore the detected artifacts in the databrowser
ft_databrowser(cfg, data_clip);

%%

cfg.artfctdef.reject = 'partial';
cfg.artfctdef.minaccepttim = 0;
cfg.artfctdef.feedback = 'yes';
data_rej = ft_rejectartifact(cfg, data_clip);
if ~istable(data.trialinfo)
  assert(isequal(data_rej.trialinfo, [1 1 2 2 2 3 4 5 6 7 8 9 10]'));
end

cfg.artfctdef.reject = 'complete';
data_rej = ft_rejectartifact(cfg, data_clip);
if ~istable(data.trialinfo)
  assert(isequal(data_rej.trialinfo, [3 4 5 6 7 8 9 10]'));
end

cfg.artfctdef.reject = 'nan';
data_rej = ft_rejectartifact(cfg, data_clip);
if ~istable(data.trialinfo)
  assert(isequal(data_rej.trialinfo, [1 2 3 4 5 6 7 8 9 10]'));
end

%%
close all

data_thresh = data;
data_thresh.trial{1}(1,201:700) = data_thresh.trial{1}(1,201:700) + 11;
data_thresh.trial{2}(2,301:400) = data_thresh.trial{2}(2,301:400) + 12;
data_thresh.trial{2}(3,501:600) = data_thresh.trial{2}(3,501:600) - 13;

cfg = [];
cfg.artfctdef.threshold.channel = 'all';
cfg.artfctdef.threshold.bpfilter = 'no';
cfg.artfctdef.threshold.min = -8;
cfg.artfctdef.threshold.max =  8;

[cfg, artifact] = ft_artifact_threshold(cfg, data_thresh);

% explore the detected artifacts in the databrowser
ft_databrowser(cfg, data_thresh);

cfg.artfctdef.reject = 'partial';
cfg.artfctdef.minaccepttim = 0;
cfg.artfctdef.feedback = 'yes';
data_rej = ft_rejectartifact(cfg, data_thresh);
if ~istable(data.trialinfo)
  assert(isequal(data_rej.trialinfo, [1 1 2 2 2 3 4 5 6 7 8 9 10]'));
end

cfg.artfctdef.reject = 'complete';
data_rej = ft_rejectartifact(cfg, data_thresh);
if ~istable(data.trialinfo)
  assert(isequal(data_rej.trialinfo, [3 4 5 6 7 8 9 10]'));
end

cfg.artfctdef.reject = 'nan';
data_rej = ft_rejectartifact(cfg, data_thresh);
if ~istable(data.trialinfo)
  assert(isequal(data_rej.trialinfo, [1 2 3 4 5 6 7 8 9 10]'));
end

%%
% Although I identified in issue 1564 that the detected jump artifacts could also
% have their channel in the outpout table, in the end I did not implement it. It
% would be a lot of work and not really worth it.

close all

data_jump = data;
data_jump.trial{1}(1,401:end) = data_jump.trial{1}(1,401:end) + 15;
data_jump.trial{2}(2,601:end) = data_jump.trial{2}(2,601:end) - 15;

cfg = [];
cfg.artfctdef.jump.channel = 'all';
cfg.artfctdef.jump.trlpadding = 0;
cfg.artfctdef.jump.fltpadding = 0;
cfg.artfctdef.jump.artpadding = 0.1;
cfg.artfctdef.jump.cutoff = 20;
cfg.artfctdef.jump.interactive = 'no';

[cfg, artifact] = ft_artifact_jump(cfg, data_jump);

% explore the detected artifacts in the databrowser
ft_databrowser(cfg, data_jump);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Display Start message
disp('Running Fading Channel BER Simulation set')
disp('To produce smoother simulation curves, increase')
disp('maxNumErr inside Fading_Simulate.m to collect more')
disp('errors.')
disp('By default, this script does not run the Monte Carlo simulations.')
disp('It simply loads results from a previous run.')
sim_mode = input('Do you want to (R)un simulations or (L)oad results from disk? [R/L]','s');
if isempty(sim_mode)
    sim_mode = 'l';
end
sim_mode = lower(sim_mode);

%% Plot theoretical curves
SNRs = 0:25;
[h_fig, h_lines] = Fading_BER_Curves(SNRs);
%% Run Monte Carlo simulations
total_time = tic;
% Pre-allocate
fBER = zeros(6,length(SNRs));
figure(h_fig)
drawnow
if strcmp(sim_mode,'r')
    hold on
    % Create place-holder plots
    simLines = semilogy(SNRs, fBER,'*');
    fBER = Fading_Simulate(SNRs, simLines);
else
    % Load from Disk (NO SIMULATIONS ARE RUN)
    load fBER  % Loads simulation results from disk instead of running them
    hold on
    simLines = semilogy(SNRs, fBER,'*');
end
toc(total_time)

%% Add Legend
% Do MATLAB graphics magic to create concise legend
% See "Controling Legends" in MATLAB doc

% group lines together
simGrp = hggroup('DisplayName','Simulation');
theoGrp = hggroup('DisplayName','Theoretical');
set(simLines,'Parent',simGrp)
set(h_lines,'Parent',theoGrp)
set(get(get(simGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
set(get(get(theoGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
legend show

figure1 = h_fig;
% Create textbox
annotation(figure1,'textbox',[0.6046 0.6266 0.1248 0.04538],...
    'String',{'No Diversity'},...
    'EdgeColor','none',...
    'BackgroundColor',[1 1 1]);

% Create textbox
annotation(figure1,'textbox',[0.63 0.4612 0.07359 0.04538],...
    'String',{'L = 2'},...
    'EdgeColor','none',...
    'BackgroundColor',[1 1 1]);

% Create textbox
annotation(figure1,'textbox',[0.6018 0.3718 0.07359 0.04538],...
    'String',{'L = 3'},...
    'EdgeColor','none',...
    'BackgroundColor',[1 1 1]);

% Create textbox
annotation(figure1,'textbox',[0.321 0.2733 0.08804 0.04538],...
    'String',{'AWGN'},...
    'EdgeColor','none',...
    'BackgroundColor',[1 1 1]);

% Create textarrow
annotation(figure1,'textarrow',[0.2853 0.4502],[0.5005 0.5],...
    'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'L = 4'});

% Create textarrow
annotation(figure1,'textarrow',[0.2912 0.4502],[0.4445 0.4455],...
    'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'L = 6'});

% Create textarrow
annotation(figure1,'textarrow',[0.3091 0.4517],[0.3945 0.3955],...
    'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'L = 8'});

figure(h_fig)


text:
function esvm_show_M(models, M)
% Show the M calibration matrix which, as the co-activation matrix, is
% a type of visual memex matrix.
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

for i = 1:length(models)

  [aa,bb] = sort(M.w{i},'descend');

  figure(1)
  clf
  subplot(4,4,1)
  imagesc(esvm_get_exemplar_icon(models,i))
  axis image
  axis off
  title(sprintf('Source exemplar %d',i))

  for q = 1:15
    %if aa(q) <= 1
    %  break
    %end
    subplot(4,4,q+1)

    flipper = 0;
    if bb(q) > length(models)
      bb(q) = bb(q) - length(models);
      flipper = 1;
      I=flip_image(esvm_get_exemplar_icon(models,bb(q)));
    else
      I=esvm_get_exemplar_icon(models,bb(q));
    end
    imagesc(I)
    axis image
    axis off
    title(sprintf('%d F=%d %.2f',bb(q),flipper,aa(q)))
  end

  pause
end

author: quantombone

repo: exemplarsvm


text:
function Dist = bst_surfdist(Points, Vertices, Faces)
% BST_SURFDIST: Compute the distances between points and a surface.
%
% USAGE:  Dist = bst_surfdist(Points, Vertices, Faces)
%
% DESCRIPTION:
%     Exact distance computation, which checks all 3 sets of distances: points
%     to vertices, points to edges, and points to faces, keeping the smallest
%     for each point.
%
% INPUTS:
%    - Points   : [Qx3] double matrix, points to compare to the mesh defined by Vertices/Faces
%    - Vertices : [Mx3] double matrix
%    - Faces    : [Nx3] double matrix
%
% OUTPUTS:
%    - Dist     : [Qx1] final distance between points and mesh

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Marc Lalancette, 2022

% TODO: A bit slow, look for alternatives
% This seems similar: https://www.mathworks.com/matlabcentral/fileexchange/52882-point2trimesh-distance-between-point-and-triangulated-surface
Epsilon = 1e-9; % nanometer
nP = size(Points, 1);
nF = size(Faces, 1);

% Prepare surface quantities, independent of points
% (In bst_meshfit, this can be done only once before iterative fitting.)
% Edges as indices
Edges = unique(sort([Faces(:,[1,2]); Faces(:,[2,3]); Faces(:,[3,1])], 2), 'rows');
% Edge direction "doubly normalized" so that later projection should be between 0 and 1.
EdgeDir = Vertices(Edges(:,2),:) - Vertices(Edges(:,1),:);
EdgeL = sqrt(sum(EdgeDir.^2, 2));
EdgeDir = bsxfun(@rdivide, EdgeDir, EdgeL);
% Edges as vectors
EdgesV = zeros(nF, 3, 3);
EdgesV(:,:,1) = Vertices(Faces(:,2),:) - Vertices(Faces(:,1),:);
EdgesV(:,:,2) = Vertices(Faces(:,3),:) - Vertices(Faces(:,2),:);
EdgesV(:,:,3) = Vertices(Faces(:,1),:) - Vertices(Faces(:,3),:);
% First edge to second edge: counter clockwise = up
FaceNormals = cross(EdgesV(:,:,1), EdgesV(:,:,2));
%FaceArea = sqrt(sum(FaceNormals.^2, 2));
FaceNormals = bsxfun(@rdivide, FaceNormals, sqrt(sum(FaceNormals.^2, 2)));
% Perpendicular vectors to edges, pointing inside triangular face.
for e = 3:-1:1
    EdgeTriNormals(:,:,e) = cross(FaceNormals, EdgesV(:,:,e));
end
FaceVertices = zeros(nF, 3, 3);
FaceVertices(:,:,1) = Vertices(Faces(:,1),:);
FaceVertices(:,:,2) = Vertices(Faces(:,2),:);
FaceVertices(:,:,3) = Vertices(Faces(:,3),:);

% Check distance to vertices
if license('test','statistics_toolbox')
    DistVert = pdist2(Vertices, Points, 'euclidean', 'Smallest', 1)';
else
    DistVert = zeros(nP, 1);
    for iP = 1:nP
        % Find closest surface vertex.
        DistVert(iP) = sqrt(min(sum(bsxfun(@minus, Points(iP, :), Vertices).^2, 2)));
    end
end
% Check distance to faces
DistFace = inf(nP, 1);
for iP = 1:nP
    % Considered Möller and Trumbore 1997, Ray-Triangle Intersection (https://stackoverflow.com/questions/42740765/intersection-between-line-and-triangle-in-3d), but this is simpler still.
    % Vectors from triangle vertices to point.
    Pyramid = bsxfun(@minus, Points(iP, :), FaceVertices);
    % Does the point project inside each face?
    InFace = all(sum(Pyramid .* EdgeTriNormals, 2) > -Epsilon, 3);
    if any(InFace)
        DistFace(iP) = min(abs(sum(Pyramid(InFace,:,1) .* FaceNormals(InFace,:), 2)));
    end
end
% Check distance to edges
DistEdge = inf(nP, 1);
for iP = 1:nP
    % Vector from first edge vertex to point.
    Pyramid = bsxfun(@minus, Points(iP, :), Vertices(Edges(:, 1), :));
    Projection = sum(Pyramid .* EdgeDir, 2);
    InEdge = Projection > -Epsilon & Projection < (EdgeL + Epsilon);
    if any(InEdge)
        DistEdge(iP) = sqrt(min(sum((Pyramid(InEdge,:) - bsxfun(@times, Projection(InEdge), EdgeDir(InEdge,:))).^2, 2)));
    end
end

Dist = min([DistVert, DistEdge, DistFace], [], 2);
end

author: brainstorm-tools


text:
function sFileOut = out_fopen_edf(OutputFile, sFileIn, ChannelMat, EpochSize, iChannels)
% OUT_FOPEN_EDF: Saves the header of a new empty EDF file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Martin Cousineau & Francois Tadel, 2017-2019

%% ===== PARSE INPUTS =====
% iChannels is set only if exporting only a subset of the input channels
if (nargin < 5) || isempty(iChannels)
    iChannels = [];
end
% Reject files with epochs
if (length(sFileIn.epochs) > 1)
    error('Cannot export epoched files to continuous EDF files.');
end
% Is the input file a native EDF file
fileSamples = round(sFileIn.prop.times .* sFileIn.prop.sfreq);
isRawEdf = strcmpi(sFileIn.format, 'EEG-EDF') && ~isempty(sFileIn.header) && isfield(sFileIn.header, 'patient_id') && isfield(sFileIn.header, 'signal');
nSamples = fileSamples(2) - fileSamples(1) + 1;
% Modify input headers (EDF export reuses the input header info directly)
if isRawEdf && ~isempty(iChannels)
    sFileIn.header.nsignal = length(iChannels);
    sFileIn.header.signal = sFileIn.header.signal(iChannels);
end

%% ===== GET MAXIMUM VALUES =====
if ~isRawEdf
    % Extracts the minimum and maximum values for each sensor over all the file (by blocks of 1s).
    % This helps optimizing the conversion of the recordings to int16 values.
    BlockSize = sFileIn.prop.sfreq; % 1-second blocks
    nBlocks = ceil(nSamples ./ BlockSize);
    % Initialize min/max matrices
    Fmax = 0 * ones(length(ChannelMat.Channel), 1);
    % Loop on all the blocks
    for iBlock = 1:nBlocks
        bst_progress('text', sprintf('Finding maximum values [%d%%]', round(iBlock/nBlocks*100)));
        % Get sample indices for a block of 1s
        SamplesBounds = round([(iBlock - 1) * BlockSize + fileSamples(1), min(fileSamples(2), fileSamples(1) + iBlock * BlockSize)]);
        % Read the block from the file
        Fblock = in_fread(sFileIn, ChannelMat, 1, SamplesBounds);
        % Keep only the files to be saved in the output file
        if ~isempty(iChannels)
            Fblock = Fblock(iChannels, :);
        end
        % Extract absolute max
        Fmax = max(Fmax, max(abs(Fblock),[],2));
    end
    % Make sure we don't have cases where the maximum is zero
    Fmax(Fmax == 0) = 1;
else
    Fmax =  1 * ones(length(ChannelMat.Channel), 1);
end

%% ===== WRITE EDF HEADER =====
bst_progress('text', 'Writing EDF+ header...');
% Get file comment
[fPath, fBase, fExt] = bst_fileparts(OutputFile);

% Initialize output file
sFileOut = sFileIn;
sFileOut.filename  = OutputFile;
sFileOut.condition = '';
sFileOut.format    = 'EEG-EDF';
sFileOut.byteorder = 'l';
sFileOut.comment   = fBase;
date = datetime;

% Create a new header structure
header            = struct();
header.version    = 0;
header.startdate  = datestr(date, 'dd.mm.yy');
header.starttime  = datestr(date, 'HH.MM.SS');
header.nsignal    = length(ChannelMat.Channel);
header.nrec       = nSamples / sFileIn.prop.sfreq;

% We need to choose a record length that produces a  record size less than
% 61440 bytes (at 2 bytes per sample).
% We'll enforce this only if the epoch size is divisible to avoid problems.
header.reclen = EpochSize / sFileIn.prop.sfreq;
factors = full_factors(EpochSize);

for iFactor = 1:length(factors)
    header.reclen = EpochSize / sFileIn.prop.sfreq / factors(iFactor);
    recordSize = EpochSize / factors(iFactor) * header.nsignal * 2;

    if recordSize <= 61440
        break;
    end
end
header.nrec = ceil(header.nrec / header.reclen);

% Add an additional channel at the end to save events if necessary.
% if ~isempty(sFileIn.events)
    % If the last channel is not already an Annotation channel: add it
    if ~strcmpi(ChannelMat.Channel(end).Name, 'Annotations') || ~strcmpi(ChannelMat.Channel(end).Type, 'EDF')
        header.nsignal = header.nsignal + 1;
        Fmax(end+1) = 1;
    end
    header.annotchan = header.nsignal;

    % Some EDF+ fields are required by strict viewers such as EDFbrowser
    header.unknown1    = 'EDF+C';
    header.patient_id  = 'UNKNOWN M 01-JAN-1900 Unknown_Patient';
    header.rec_id      = ['Startdate ', upper(datestr(date, 'dd-mmm-yyyy')), ...
                          ' Unknown_Hospital Unknown_Technician Unknown_Equipment'];

    % Compute annotations
    header.annotations = {};
    maxAnnotLength     = 0;

    for iEvt = 1:numel(sFileIn.events)
        event = sFileIn.events(iEvt);
        % EDF file start at 0s: removed the start file time
        event.times = event.times - sFileIn.prop.times(1);

        for iEpc = 1:length(event.epochs)
            startTime = event.times(1,iEpc);
            annot = sprintf('+%f', startTime);

            if (size(event.times,1) == 2)
                duration = event.times(2,iEpc) - startTime;
                annot    = [annot, sprintf('%c%f', char(21), duration)];
            end

            annot = [annot, sprintf('%c%s%c%c', char(20), event.label, char(20), char(0))];
            header.annotations{end + 1} = annot;

            if length(annot) > maxAnnotLength
                maxAnnotLength = length(annot);
            end
        end
    end
% else
%     header.annotchan   = -1;
%     header.unknown1    = '';
%     header.patient_id  = '';
%     header.rec_id      = '';
% end
header.hdrlen = 256 + 256 * header.nsignal;

% Channel information
for i = 1:header.nsignal
    if ~isRawEdf
        % Choose unit based on maximal value
        if Fmax(i) > 0.1
            header.signal(i).unit = 'V';
            header.signal(i).unit_gain = 1;
%         elseif Fmax(i) * 1e3 > 1
%             header.signal(i).unit = 'mV';
%             header.signal(i).unit_gain = 1e3;
        else
            header.signal(i).unit = 'uV';
            header.signal(i).unit_gain = 1e6;
        end

        chScale = 2^15 / Fmax(i);
        header.signal(i).digital_min  = -2^15;
        header.signal(i).digital_max  = 2^15 - 1;
        header.signal(i).physical_min = header.signal(i).digital_min / chScale .* header.signal(i).unit_gain;
        header.signal(i).physical_max = header.signal(i).digital_max / chScale .* header.signal(i).unit_gain;
        header.signal(i).filters      = '';
        header.signal(i).unknown2     = '';
    end

    % Approximate size of annotation channel
    header.signal(i).type  = '';
    if i == header.annotchan
        header.signal(i).label    = 'EDF Annotations';
        eventsPerRecord           = ceil(numel(header.annotations) / header.nrec);
        header.signal(i).nsamples = eventsPerRecord * maxAnnotLength + 15; % For first annotation of each record
        % Convert chars (1-byte) to 2-byte integers, the size of a sample
        header.signal(i).nsamples = int64((header.signal(i).nsamples + 1) / 2);
    else
        if strcmpi(ChannelMat.Channel(i).Type, 'EEG_NO_LOC')
            header.signal(i).label = ['EEG ', ChannelMat.Channel(i).Name];
        else
            header.signal(i).label = [ChannelMat.Channel(i).Type, ' ', ChannelMat.Channel(i).Name];
        end
        header.signal(i).nsamples = header.reclen * sFileIn.prop.sfreq;
    end
end

% Copy some values from the original header if possible
if isRawEdf
    header.patient_id = sFileIn.header.patient_id;
    header.rec_id     = sFileIn.header.rec_id;
    header.startdate  = sFileIn.header.startdate;
    header.starttime  = sFileIn.header.starttime;
    header.unknown1   = sFileIn.header.unknown1;
    for i = 1:sFileIn.header.nsignal
        header.signal(i).label        = sFileIn.header.signal(i).label;
        header.signal(i).type         = sFileIn.header.signal(i).type;
        header.signal(i).filters      = sFileIn.header.signal(i).filters;
        header.signal(i).unit         = sFileIn.header.signal(i).unit;
        header.signal(i).unknown2     = sFileIn.header.signal(i).unknown2;
        header.signal(i).physical_min = sFileIn.header.signal(i).physical_min;
        header.signal(i).physical_max = sFileIn.header.signal(i).physical_max;
        header.signal(i).digital_min  = sFileIn.header.signal(i).digital_min;
        header.signal(i).digital_max  = sFileIn.header.signal(i).digital_max;

        switch (header.signal(i).unit)
            case 'mV',                        header.signal(i).unit_gain = 1e3;
            case {'uV', char([166 204 86])},  header.signal(i).unit_gain = 1e6;
            otherwise,                        header.signal(i).unit_gain = 1;
        end
    end
end

% Open file
fid = fopen(OutputFile, 'w+', sFileOut.byteorder);
if (fid == -1)
    error('Could not open output file.');
end

% Write header
fwrite(fid, str_zeros(header.version, 8));
fwrite(fid, str_zeros(header.patient_id, 80));
fwrite(fid, str_zeros(header.rec_id, 80));
fwrite(fid, header.startdate);
fwrite(fid, header.starttime);
fwrite(fid, str_zeros(header.hdrlen, 8));
fwrite(fid, str_zeros(header.unknown1, 44));
fwrite(fid, str_zeros(header.nrec, 8));
fwrite(fid, str_zeros(header.reclen, 8));
fwrite(fid, str_zeros(header.nsignal, 4));

% Channel information
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).label, 16));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).type, 80));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).unit, 8));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).physical_min, 8));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).physical_max, 8));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).digital_min, 8));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).digital_max, 8));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).filters, 80));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).nsamples, 8));
end
for i = 1:header.nsignal
    fwrite(fid, str_zeros(header.signal(i).unknown2, 32));
end

% Sanity check on the written header
assert(ftell(fid) == header.hdrlen);

% Close file
fclose(fid);
% Copy header to the sFile structure
sFileOut.header = header;

end

%% ===== HELPER FUNCTIONS =====
function sout = str_zeros(sin, N)
    if (isnumeric(sin))
        sin_str = num2str(sin);
        if (length(sin_str) > N)
            if ((sin > 10^-(N-3)) && (sin < 10^N-1)) || (N < 7)
                sin_str = sprintf('%f', sin);
            else
                sin_str = sprintf(sprintf('%%0.%de',N-7), sin);
            end
        end
        sin = sin_str;
    end

    sout = char(double(' ') * ones(1,N));
    if (length(sin) <= N)
        sout(1:length(sin)) = sin;
    else
        sout = sin(1:N);
    end
end

% This extracts a subset of the full factors of a number, not just the
% prime factors like factors() returns.
% E.g.: factor(100)      = [2,2,5,5]
%       full_factors(100) = [2,4,5,10,20,25,50]
function factors = full_factors(n)
    factors = factor(n);
    if length(factors) == 1
        factors = [];
        return;
    end
    [occurences, factors] = hist(factors, unique(factors));
    factors = factors(occurences > 0);
    occurences = occurences(occurences > 0);
    nFactors = length(factors);
    possibleFactors = ones(sum(occurences) * nFactors * nFactors);
    iPos = 1;
    for iFactor = 1:nFactors
        for iOcc = 1:occurences(iFactor)
            fact = factors(iFactor) ^ iOcc;
            possibleFactors(iPos) = fact;
            iPos = iPos + 1;
            for iFactor2 = 1:nFactors
                if iFactor ~= iFactor2
                    possibleFactors(iPos) = fact * factors(iFactor2);
                    iPos = iPos + 1;
                end
            end
        end
    end
    factors = unique(possibleFactors(1:iPos-1));
end


text:
% GSPBOX - Graphs
%
%  Specific graphs
%    gsp_swiss_roll              -  Create swiss roll graph
%    gsp_david_sensor_network    -  Create the sensor newtwork from david
%    gsp_ring                    -  Create the ring graph
%    gsp_path                    -  Create the path graph
%    gsp_airfoil                 -  Create the airfoil graph
%    gsp_comet                   -  Create the comet graph
%    gsp_erdos_renyi             -  Create a erdos renyi graph
%    gsp_minnesota               -  Create Minnesota road graph
%    gsp_low_stretch_tree        -  Create a low stretch tree graph
%    gsp_sensor                  -  Create a random sensor graph
%    gsp_random_regular          -  Create a random regular graph
%    gsp_random_ring             -  Create a random ring graph
%    gsp_full_connected          -  Create a fully connected graph
%    gsp_nn_graph                -  Create a nearest neighbors graph
%    gsp_rmse_mv_graph           -  Create a nearest neighbors graph with missing values
%    gsp_sphere                  -  Create a spherical-shaped graph
%    gsp_cube                    -  Create a cubical-shaped graph
%    gsp_2dgrid                  -  Create a 2d-grid graph
%    gsp_torus                   -  Create a torus graph
%    gsp_logo                    -  Create a GSP logo graph
%    gsp_community               -  Create a community graph
%    gsp_bunny                   -  Create a bunny graph
%    gsp_spiral                  -  Create a spiral graph
%    gsp_stochastic_block_graph  -  Create a graph with the stochastic block model
%
%  Hypergraphs
%    gsp_nn_hypergraph           -  Create an hyper nearest neighbor graph
%    gsp_hypergraph              -  Create an hypergraph
%
%  Utils
%    gsp_graph_default_parameters-  Initialise all parameters for a graph
%    gsp_graph_default_plotting_parameters-  Initialise all plotting parameters for a graph
%    gsp_graph                   -  Create a graph from a weight matrix
%    gsp_update_weights          -  Update the weights of a graph
%    gsp_update_coordinates      -  Update the coordinate of a graph
%    gsp_components              -  Cuts non connected graph into several connected ones
%    gsp_subgraph                -  Create a subgraph
%    gsp_graph_product           -  Compute graph product between two graphs
%    gsp_line_graph              -  Create the Line Graph (or edge-to-vertex dual graph) of a graph
%    gsp_jtv_graph               -  Add time information to the graph structure
%
%
%  For help, bug reports, suggestions etc. please send email to
%  gspbox 'dash' support 'at' groupes 'dot' epfl 'dot' ch
%

repo: gspbox


text:
function opt= defopt_erps(varargin)
%DEFOPT_ERPS - Properties for plot_channel* functions, optimized for
%plotting ERPs

props= {'LineWidth',                    0.7
        'AxisType',                     'cross'
        'AxisTitleVerticalAlignment',   'top'
        'AxisTitleFontWeight',          'demi'
        'ScaleHPos',                    'left'
        'ColorOrder',             [1 .47 .02; .3 .3 .3; 0 .6 .8; .9 0 .9],
        };

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);

author: bbci


text:
function y = abs( x )

%Disciplined convex/geometric programming information for ABS:
%   ABS(X) is convex and nonmonotonic in X. Therefore, according to
%   the strict rules of DCP, X must be affine. However, because of
%   its special structure, CVX considers the sign of X as well. So,
%   for instance, if X is known to be nonnegative, then ABS(X)=X.

persistent P
if isempty( P ),
    P.map = cvx_remap( { 'constant' }, { 'p_nonconst' }, { 'n_nonconst' }, ...
        { 'r_affine' }, { 'c_affine' } );
    P.funcs = { @abs_cnst, @abs_posn, @abs_negn, @abs_affn, @abs_affn };
end
y = cvx_unary_op( P, x );

function y = abs_cnst( x )
% Constant
y = builtin( 'abs', x );

function y = abs_posn( x )
% Positive any
y = x;

function y = abs_negn( x )
% Negative any
y = -x;

function y = abs_affn( x ) %#ok
% Affine
cvx_begin
    epigraph variable y( size(x) ) nonnegative_
    { x, y } == lorentz( size(x), 0, ~isreal(x) ); %#ok
cvx_end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% PARSE A GENERIC INPUT SET AGAINST A DEFAULT CONFIG STRUCTURE
function [config] = GetParameterOverrides(defaultConfig,inputParameters)
% This function is designed to compare a vector of {'input',parameter} pairs
% against a structure(or object) with labelled fields.

% INPUTS:
% defaultConfig     - The reference structure/class with available fields.
% inputParameters   - The cell array of string, parameter pairs or
%                     comparative structure.
% OUTPUT:
% config            - The structure/class with the fields overitten.

% Input sanity check #1
if nargin < 2 || numel(inputParameters) == 0
    config = defaultConfig;
    return
end

% Input sanity check #2 - Check nested cell arrays
if iscell(inputParameters)
    t = 1;
    while numel(inputParameters) == 1 && t < 10
        inputParameters = inputParameters{:};
        t = t + 1;
    end
    assert(t ~= 10,'Failed to unload parameters.');
end

% Input-dependant behaviour
switch true
    case iscell(inputParameters)
        assert(mod(numel(inputParameters),2) == 0,'Entries are not of list of parameter:value pairs');
        config = parseCellArray(defaultConfig,inputParameters);
        return
    case isstruct(inputParameters)
        config = parseStructure(defaultConfig,inputParameters);
        return
    otherwise
        error('Please provide a list of parameter:value pairs or a structure.');
end
end

% Parse the input parameters against a list of parameter:value cell array
function [config] = parseCellArray(config,inputParameters)

% Move through the parameter pairs (i.e. ,'radius',1,...)
for parameterIndex = 1:numel(inputParameters)
    % For each user-defined parameter
    givenParameter = inputParameters{parameterIndex};
    if ~ischar(givenParameter)
        continue
    end

    % IF THE OBJECT HAS A PROPERTY BY THAT NAME
    if isstruct(config) && isfield(config,givenParameter)
        config.(givenParameter) = inputParameters{parameterIndex + 1};  % Make a substitution
    end
    if ~isstruct(config) && isprop(config,givenParameter)
        config.(givenParameter) = inputParameters{parameterIndex + 1};   % Make a substitution
    end
end
end
% Parse the input parameters against a structure
function [config] = parseStructure(config,inputStructure)

% Move through the fields of the second structure and override associated
% fields.
fields = fieldnames(inputStructure);
for f = 1:numel(fields)
    % Config is a struct
    if isstruct(config) && isfield(config,fields{f})
        config.(fields{f}) = inputStructure.(fields{f});        % Make a substitution
        continue
    end
    % Config is a class
    if ~isstruct(config) && isprop(config,fields{f})
        config.(fields{f}) = inputStructure.(fields{f});        % Make a substitution
        continue
    end
end
end

repo: OpenMAS


text:
function mrInit_updateSessiontSeries()
%
% USAGE: Takes a session that has already been initialized with an older
% version of mrInit and update it to the newest version. This version
% is of 2013-05-05 and it no longer uses tSeries data in the form of
% matrices, instead saving tSeries data to a nifti and using that.
%
% INPUT: N/A, no input is used. As long as the user is in the directory of
% the session to be migrated, this should work correctly.
%
% OUTPUT: N/A, no output is given. The error handling has been upgraded to
% use built in matlab try-catch blocks.
%
% This migration tool takes a series of folders of tSeries*.mat files and
% makes some assumptions about their orientation. Specifically, it is assumed
% that these are already in the normal display format.
% All of the tSeries data is loaded and the migration tool creates
% a nifti structure around this data, saves it to the filesystem and writes
% the information about its location in the session and datatype global
% variables, before saving these to mrSESSION.mat as well.

try
    loadSession;
    mrGlobals;

    %Before we reset mrSESSION, let's save a backup
    copyfile('mrSESSION.mat','mrSESSION_tSeriesMigrationBackup.mat');

    %Now that we have the number of scans, we know how many tSeries nifti
    %files we will need to create

    % Use local paths, not absolute paths
    % inplaneBasePath = fullfile(pwd,'Inplane');
    inplaneBasePath = 'Inplane';

    for dtNum = 1:numel(dataTYPES)
        fprintf('Starting dataTYPE number %d\n', dtNum);
        tSeriesOutPath = fullfile(inplaneBasePath, dtGet(dataTYPES(dtNum),'Name'));
        tSeriesInBasePath = fullfile(tSeriesOutPath,'TSeries');

        if exist(fullfile(tSeriesInBasePath,'Scan1','tSeries1.mat'),'file')
            %Checks to see if there is any inplane tSeries data in this
            %dataTYPE if there is, then does all of the processing

            numScans = dtGet(dataTYPES(dtNum),'N Scans');

            keepFrames = zeros(numScans,2);

            for scan = 1:numScans
                fprintf('Starting scan number %d\n', scan);
                % For each scan, go through each scan directory, read in all of the
                % matrix files and then build the data for a nifti from them.
                numSlices = dtGet(dataTYPES(dtNum),'N Slices', scan);
                tSeriesInFolder = fullfile(tSeriesInBasePath,['Scan' num2str(scan)]);
                dimSize = [dtGet(dataTYPES(dtNum),'N Frames', scan) dtGet(dataTYPES(dtNum),'Func Size', scan) dtGet(dataTYPES(dtNum),'N Slices', scan)];
                tSeries = zeros(dimSize);
                for slice = 1:numSlices
                    tSeriesInFile = fullfile(HOMEDIR, tSeriesInFolder,['tSeries' num2str(slice) '.mat']);
                    tSeriesIn = load(tSeriesInFile);
                    tSeriesIn = reshape(tSeriesIn.tSeries, dimSize(1:3));
                    tSeries(:,:,:,slice) = tSeriesIn;
                end %for

                % The following transform will preserve the orientation of the data
                % across the migration, such that inplane to vAnatomy xforms do not
                % need to change, nor do parameter maps, ROIs, etc. The permute([2 1])
                % followed by flipdim(2) effectively converts from x/y coordinates to
                % row/column coordinates. The reason we do this is that the old
                % vistasoft (prior to the move to NIFTIs) did this transform when
                % showing the slices as images, whereas the current code simply pulls
                % the data array from a nifti (after applying a standard xform to the
                % nifti), and uses an image tool such as imagesc slice by slice without
                % re-orienting the array. Since we no longer do this xform each time we
                % show the data, we must do it here as part of the migration if want
                % the migrated data to appear the same way as the pre-migrated data.
                % See also mrInit_updateInplaneSession.m.
                tSeries = permute(tSeries,[3 2 4 1]); %Standard format: freq phase slice time
                tSeries = flip(tSeries, 2);
                %Note: this needed to be changed to reflect the fact that
                %MATLAB stores values row, column, etc. and not column,
                %row,

                %Create the freq, phase and slice dimensions, assuming that we are in
                %standard format
                freqPhaseSliceDims = [1 2 3];

                %Create the slice information
                % We are assuming that the voxels are the same size in all
                % of the dataTYPES
                funcVoxel = sessionGet(mrSESSION,'Functional Voxel Size',scan);

                xform = [[diag(1./funcVoxel); 0 0 0], size(tSeries)'/2];

                funcVoxel(4) = dtGet(dataTYPES(dtNum),'Frame Period',scan);

                sliceInfo = [3 0 dtGet(dataTYPES(dtNum),'N Slices',scan)-1 funcVoxel(4)];

                %Build the nifti from the components above
                nii = niftiCreate('data',tSeries,'qto_xyz',xform,'freq_dim',freqPhaseSliceDims,'slice_code',sliceInfo);

                %However, this does not create the proper pix dims, so let's fix that
                nii = niftiSet(nii,'Pix dim',funcVoxel);

                keepFrames(scan, 1) = 0;
                keepFrames(scan, 2) = -1;

                mrSESSION = sessionSet(mrSESSION,'Keep Frames',keepFrames, scan);

                tSeriesOut = fullfile(tSeriesInBasePath,['tSeriesScan' num2str(scan) '.nii.gz']);

                nii = niftiSet(nii,'File Path',tSeriesOut);

                dataTYPES(dtNum) = dtSet(dataTYPES(dtNum),'Inplane Path',tSeriesOut,scan);
                dataTYPES(dtNum) = dtSet(dataTYPES(dtNum),'Keep Frames', keepFrames, scan);

                niftiWrite(nii,tSeriesOut);

                mrSESSION = sessionSet(mrSESSION,'Version','2.1');

                %Update the session variables
                save('mrSESSION.mat', 'mrSESSION','-append');
                save('mrSESSION.mat', 'dataTYPES','-append');

                % we already wrote the nifti 4 lines above. this is
                % redundant.
                %   writeFileNifti(nii);

                fprintf('Finished scan number %d\n', scan)
            end %for

        end %if
        fprintf('Finished dataTYPE number %d\n', dtNum)
    end %for

catch err
    warning(['There was an error when attempting to update your session.',...
        'No changes have been made to your system. Please run the update code again.']);
    rethrow(err);
end %try

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [hPlot] = plot_FMDboth(mCatalog, fMc, fAvalue, fBvalue, fBinning, nMinNum)
% function [hPlot, fBValue, fAValue, fStdDev, fMc, fMeanMag]
%   = plot_FMD(mCatalog, bCumulative, hAxes, sSymbol, sColor, bPlotB, nCalculateMC, fBinning)
% -------------------------------------------------------------------------------------------
% Plots a cumulative and non-cumulative frequency magnitude distribution including the b-value
%
% plot_FMD(mCatalog) opens a figure and plots the frequency magnitude distribution
%   with standard parameters
%
% Input parameters:
%   mCatalog        Earthquake catalog
%   fMc             Magnitude of completeness
%   fAValue         Calculateda-value
%   fBValue         Calculated b-value
%   fBinning        Magnitude binning of the catalog (default 0.1)
%   nMinNum         Minimum number of events
%
% Output parameters:
%   hPlot           Handle of the plot
%
% jowoe@gps.caltech.edu
% 25.01.2006

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Catalog size
[nRow,nCol] = size(mCatalog);

% Create the frequency magnitude distribution vector
[vFMD, vNonCFMD] = calc_FMD(mCatalog);

if ~exist('fBinning', 'var')
  fBinning = 0.1;
end

figure
% Plot cumulative frequency magnitude distribution
hPlot = semilogy(vFMD(1,:), vFMD(2,:));
set(hPlot,'Marker','^','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0.75 0.75 0.75],'Markersize',8,'Linewidth',2,'Linestyle','none');
hold on
hPlot2 = semilogy(vNonCFMD(1,:), vNonCFMD(2,:));
set(hPlot2,'Marker','d','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0.4 0.4 0.4],'Markersize',8,'Linewidth',2,'Linestyle','none');

% Add further plots to the axes
set(gca, 'NextPlot', 'add');

if (~isnan(fMc) & ~isnan(fBvalue) & nRow >= 2*nMinNum)
    % Determine the positions of 'x'-markers
    nIndexLo = find((vFMD(1,:) < fMc + 0.05) & (vFMD(1,:) > fMc - 0.05));
    fMagHi = vFMD(1,1);
    vSel = vFMD(1,:) <= fMagHi & vFMD(1,:) >= fMc-.0001;
    vMagnitudes = vFMD(1,vSel);

    % Plot the 'x'-marker
    hPlot = semilogy(vFMD(1,nIndexLo), vFMD(2,nIndexLo), 'x','Color',[0.8 0 0]);
    set(hPlot, 'LineWidth', [2.5], 'MarkerSize', 12);
    hPlot = semilogy(vFMD(1,1), vFMD(2,1), 'x','Color',[0.8 0 0]);
    set(hPlot, 'LineWidth', [2.5], 'MarkerSize', 12)

    % Plot the line representing the b-value
    vPoly = [-1*fBvalue fAvalue];
    fBFunc = 10.^(polyval(vPoly, vMagnitudes));
    hPlot = semilogy(vMagnitudes, fBFunc);
    set(hPlot, 'LineWidth', [2.0],'Linestyle','--','Color',[0.8 0 0]);
end
% Y-Limits
vYLim = get(gca,'YLim');
set(gca,'FontSize',12','Fontweight','bold','box','on','Ylim',[1 vYLim(2)],'Linewidth',1.5)
xlabel('Magnitude','FontSize',14','Fontweight','bold')
ylabel('Number of events','FontSize',14','Fontweight','bold')

author: CelsoReyes


text:
function f = logical(f)
%LOGICAL   CHEBFUN logical.
%   LOGICAL(F) returns a CHEBFUN which evaluates to one at all points where F is
%   non-zero and zero otherwise.

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

% Add breaks at the roots (since these will take the value 0 in the output).
f = addBreaksAtRoots(f);

% Loop over the FUNs:
for k = 1:numel(f.funs)
    f.funs{k} = logical(f.funs{k});
end

% pointValues:
tol = vscale(f)*eps;
f.pointValues = abs(f.pointValues) > tol;

end


text:
function [IterationMap,Metrics,err] = bst_connectivity_mra(F, Surface, Options)
% BST_CONNECTIVITY_MRA: Multiresolution analysis of connectivity
%
% USAGE:  [ItMap,Metrics] = bst_connectivity_mra(F, Surface, HeadModel.Gain(Results.GoodChannel,:))
%
% INPUTS:
%    - F           : N-source x N-sample time-series for analysis
%    - Surface     : tesselation structure (Faces,Vertices,VertConn)
%    - Options     : N-source gain matrix from the head model
%
% OUTPUTS:
%    - IterationMap  : N-source vector with each number corresponding to
%                      the iteration at which they left the process (indicative
%                      of their probability of high connectivity)
%    - Metrics       : N-iteration cells containing connectivity metric for
%                      each iteration
%    - err           : Error message, if any
%
% NOTE:
%    - Selection model that doesn't work:
%       Using the upper half of subdivision - Weird results
%       Using a high prctile value - Too discriminatory
%
%    - Selection model that work:
%       Hard threshold works well, but arbitrary threshold
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Sebastien Dery, 2013

    Verbose = 0;
    ShuffleVertices = 1;

    nSample = size(F,2);
    nVertices = length(Surface.Vertices);
    % Init varout
    IterationMap = zeros(nVertices,1);
    Metrics = [];
    err = [];

    %S1 = extend(3664, sSurf.VertConn, 4);
    %S2 = extend(10588, sSurf.VertConn, 4);
    %Source = [find(S1(3664,:)) find(S2(10588,:))];

    isAbsolute = 1;
    nChannel = 272;
    SplitFactor = 2;

    % Unconstrained function
    isNorm = 1;
    if isNorm
        XyzFunction = 'norm';
    else
        XyzFunction = 'none';
    end
    ScoutFunction = 'mean';
    nComponents = 1;

    % Equal to the number of recorded channel
    nTessScouts = nChannel / SplitFactor;

    % Remaining vertices
    RemVert = 1:nVertices;
    % Initial tesselation
    Vert = 1:length(RemVert);
    VertConn = Surface.VertConn(RemVert,RemVert);
    Scouts = tess_cluster_seed(Vert(randperm(length(Vert))), VertConn, nTessScouts);

    iter = 1;
    nDipole = round(length(RemVert) / nTessScouts);
    while (nDipole > 2)
        % Switch variable to prevent mixing
        OldScouts = Scouts;
        % Tesselate the scout vertices
        uScouts = unique(OldScouts(OldScouts ~= 0));
        Label = 0;
        for i=1:length(uScouts)
            iRows = find(OldScouts == uScouts(i));
            Scouts(iRows) = tess_cluster(Surface.VertConn(iRows,iRows), SplitFactor, ShuffleVertices, Verbose) + Label;
            Label = Label + 2;
        end

        %
        uScouts = unique(Scouts(Scouts ~= 0));
        nScouts = length(uScouts);
        %
        Fs = zeros(nScouts, nSample);
        %
        for i=1:nScouts
            % Get the scout vertices
            iRows = find(Scouts == uScouts(i));
            % Get the scout orientation
            ScoutOrient = Surface.VertNormals(iRows,:);
            % Get scout value
            Fs(i,:) = bst_scout_value(F(iRows,:), ScoutFunction, ScoutOrient, nComponents, XyzFunction);
        end

        % Compute normalised distance
        % D = dist(Center');
        % D = (D - min(D(:))) ./ (max(D(:)) - min(D(:)));
        % D = 1 ./ D;
        % Compute metric
        switch (Options.method)
            case 'corr'
                % R = real(bst_corrn(Fs, Fs, inputs.removemean));
                R = corrcoef([Fs Fs]');
            case 'cohere'
                % R = bst_coherence(Fs, Fs, Options);
                error('Not supported yet.');
            case 'granger'
                R = bst_granger(Fs, Fs, Options.GrangerOrder, Options);
            case 'spgranger'
                error('Not supported yet.');
            case 'plv'

            otherwise
        end
        % Remove diagonal
        R(eye(length(R)) == 1) = 0;
        %
        if isAbsolute
            R = abs(R);
        end

        % Higher than mean - Rows
        tLowCorr = max(R,[],2);
        minCorr = mean(tLowCorr);
        Keep = tLowCorr > minCorr;
        % If directional data, check the other way around
        if (strcmpi(Options.method,'granger') || strcmpi(Options.method,'spgranger'))
            % Higher than mean - Columns
            tLowCorr = max(R);
            minCorr = mean(tLowCorr);
            Keep = Keep | (tLowCorr > minCorr)';
        end
        % Remove everything else
        lRow = find(~Keep);
        if isempty(lRow)
            break;
        end
        %
        RowId = ismember(Scouts,lRow);
        % Keep track of when these dipoles left the iteration process
        IterationMap(RemVert(RowId == 1)) = iter;
        % Leave these vertices out
        Scouts(RowId == 1) = 0;
        % Keep metric matrix
        Metrics{iter} = R;
        % Update iterative variables
        nDipole = floor(sum(Scouts ~= 0) / (nScouts - length(lRow)));
        iter = iter + 1;
    end
    % Assign victorious last
    IterationMap(IterationMap == 0) = iter;
end


text:
%
% example file for generating a movie from a sequence of viewpoint
%
%% Olivier Salvado, Case Western Reserve University, 16-Sep-04

% construct a 3D scene
[X,Y] = meshgrid(-3:.125:3);
Z = peaks(X,Y);
surfl(X,Y,Z);
axis vis3d
axis off

% create the structure with current view
cam = CamSeqUpdate([]);
% change to view #2
% normally the new view is changed manually using the camera toolbar
campos([20 20 0])
camtarget([1 0 0])
% add the new view to be 5s after view #1
cam = CamSeqUpdate(cam,5);

% change to view #3
% normally the new view is changed manually using the camera toolbar
campos([5 5 50])
camtarget([0 1 1])
% add the new view to be 10s after view #2
cam = CamSeqUpdate(cam,10);

% change to view #4
% normally the new view is changed manually using the camera toolbar
zoom(0.8)
% add the new view to be 5s after view #3
cam = CamSeqUpdate(cam,5);

% generate the sequence
fps = 20;   % 20 frames per second
seq = CamSeqGen(cam,fps);

% play the movie and save the file
M = CamSeqPlay(seq,fps,'GoodLookingMovie.avi');