text:
function [datsmooth] = ft_preproc_smooth(dat, n, tol)

% FT_PREPROC_SMOOTH performs boxcar smoothing with specified length.
% Edge behavior is improved by implicit padding with the mean over
% half the boxcar length at the edges of the data segment.
%
% Use as
%   [dat] = ft_preproc_smooth(dat, n)
%
% Where dat is an Nchan x Ntime data matrix, and n is the length
% of the boxcar smoothing kernel
%
% If the data contains NaNs, these are ignored for the computation, but
% retained in the output.
%
% See also PREPROC

% Undocumented options:
%  n can also be a vector containing a custom smoothing kernel
%  n can also be 'regsmooth', using a regularised estimate of the first temporal derivative
%
% Copyright (C) 2010, Stefan Klanke
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

% create smoothing kernel
regflag = false;
if isequal(n, 'regsmooth')
  regflag = true;
  n = 0;
elseif isscalar(n)
  krn = ones(1,n)/n;
else
  krn = n(:)'./sum(n(:));
  n   = numel(krn);
end

% deal with padding
dat = ft_preproc_padding(dat, 'localmean', ceil(n/2));

% do the smoothing
if regflag
  if nargin<3
    tol = 1e-9;
  end
  datsmooth = smooth_regularised(dat, tol);
elseif n<100
  % heuristic: for large kernel the convolution is faster when done along
  % the columns, weighing against the costs of doing the transposition.
  % the threshold of 100 is a bit ad hoc.
  datsmooth = convn(dat,   krn,   'same');
else
  datsmooth = convn(dat.', krn.', 'same').';
end

% cut the eges
datsmooth = ft_preproc_padding(datsmooth, 'remove', ceil(n/2));

function out = smooth_regularised(dat, tol)

tol = tol.*std(dat(:));
n   = size(dat, 2);
B    = eye(n);

% create Toeplitz matrix F
r = [1 zeros(1,n-1)];
d = [1 zeros(1,n-1)];
m = 2;
for k=1:m
 d = filter([1 -1],1,d);
end
F = toeplitz(d,r);

% create Toeplitz matrix G
c = ones(n,1);
G = toeplitz(c,r);

% compute gamma parameter

% SVD
H       = B*(G/F);
[U,D,V] = svd(H);
diagD2  = diag(D).^2;
diagD   = diag(D);
epsi    = U'*B*dat.';

gammarange = 10.^(-10:0.2:10);
ngamma     = numel(gammarange);
onevec_n   = ones(n,1);
onevec_ngamma = ones(1,ngamma);

nchan = size(dat,1);
ro    = zeros(n,ngamma,nchan);
for k = 1:nchan
  ro(:,:,k) = (gammarange(onevec_n,:).*epsi(:,k.*onevec_ngamma))./(diagD2(:,onevec_ngamma) + gammarange(onevec_n,:));
end
wrss = squeeze(sum(ro.^2));

gamma = zeros(nchan,1);
ni = zeros(n,nchan);
for k = 1:nchan
  gamma(k,1) = gammarange(find(wrss(:,k)<tol,1,'last'));
  ni(:,k) = (diagD.*epsi(:,k))./(diagD2 + gamma(k));
end

ddata_smooth = transpose(F\V*ni);
out          = ddata_smooth*G';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [dM1m dM2m dM3m Edges dEdges A1m A2m M1m M2m M3m W1m W2m COm dW1m dW2m dCOm dAreasImgm dA1m dA2m dA3m] = showResults2(yOpt, yOpt_inv, corners, filename, images)

load(filename);

m=[2^maxLevel 2^maxLevel];
Omega = [1, 1];
X=getGrid(Omega,m);
h=Omega./m;

yCentered=stg2center(yOpt,m,'Py');

% % Check whether yOpt + yOpt_inv generaets original image
% dA1i1=interpolation(TD1,Omega,yCentered);
% dA1m1=flipud(reshape(dA1i1,m)');
%
% dA1i2=interpolation(dA1m1,Omega,yOpt_inv);
% dA1m2=flipud(reshape(dA1i2,m)');
%
% figure
% subplot(311)
% imagesc(TD1)
% subplot(312)
% imagesc(dA1m1)
% subplot(313)
% imagesc(dA1m2)

% interpolate raw data
interpolation('set','MODE','linear-periodic','period',2*pi);
M2i=interpolation(images.M2,Omega,X);
dA2i=interpolation(TD2,Omega,yCentered);
interpolation('set','MODE','linear');
% M1i=interpolation(RD1,Omega,X);
M1i=interpolation(images.M1,Omega,X);
M3i=interpolation(images.M3,Omega,X);
dA1i=interpolation(TD1,Omega,yCentered);
if isfield(images,'A3')
    dA3i=interpolation(images.A3,Omega,yCentered);
    dA3m=ip2mat(dA3i,m);
end

if max(areasImg(:))==6 % LO/TO maps
    areasImg = areaNumTransform(areasImg);
end
dAreasImgi=(interpolation(areasImg,Omega,yCentered));
W1i=interpolation(WD1,Omega,X);
W2i=interpolation(WD2,Omega,X);
COi=interpolation(images.CO,Omega,X);

M2m=ip2mat(M2i,m);
M1m=ip2mat(M1i,m);
M3m=ip2mat(M3i,m);
dA1m=ip2mat(dA1i,m);
dA2m=ip2mat(dA2i,m);
W1m=ip2mat(W1i,m);
W2m=ip2mat(W2i,m);
COm=ip2mat(COi,m);
dAreasImgm=ip2mat(dAreasImgi,m);

% transform raw data into atlas space
interpolation('set','MODE','linear-periodic','period',2*pi);
dM2i=interpolation(images.M2,Omega,yOpt_inv);
A2i=interpolation(TD2,Omega,X);
interpolation('set','MODE','linear');
dM1i=interpolation(images.M1,Omega,yOpt_inv);
dM3i=interpolation(images.M3,Omega,yOpt_inv);
A1i=interpolation(TD1,Omega,X);
dW1i=interpolation(WD1,Omega,yOpt_inv);
dW2i=interpolation(WD2,Omega,yOpt_inv);
dCOi=interpolation(COm,Omega,yOpt_inv);
AreasImgi=(interpolation(areasImg,Omega,X));

dM2m=ip2mat(dM2i,m);
dM1m=ip2mat(dM1i,m);
dM3m=ip2mat(dM3i,m);
A1m=ip2mat(A1i,m);
A2m=ip2mat(A2i,m);
AreasImgm=ip2mat(AreasImgi,m);
dW1m=ip2mat(dW1i,m);
dW2m=ip2mat(dW2i,m);
dCOm=ip2mat(dCOi,m);

% remove walls to draw edges
% areasImg(find(areasImg==5))=NaN;
% areasImg(find(areasImg==6))=NaN;

% get interporated areas of interest
OV=ceil(interpolation(areasImg,Omega,X));
OVm=ip2mat(OV,m);
OVm(isnan(OVm))=0;
if isfield(images,'A3') % for LO-1/2, TO-1/2
    Edges = detectAreaBorders(OVm);
else % for V1-3
    Edges=double(edge(OVm, 'Canny'));
end

% get deformed areas of interest
OV=round(interpolation(areasImg,Omega,yCentered));
OVm=ip2mat(OV,m);
OVm(isnan(OVm))=0;
if isfield(images,'A3') % for LO-1/2, TO-1/2
    dEdges = detectAreaBorders(OVm);
else % for V1-3
    dEdges=double(edge(OVm, 'Canny'));
end

% images_size = size(images.A1,1);
% figure
% subplot(231)
% hold on
% imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(M1m,dEdges,hsv,15))
% % I=flipdim(mergedImage(dM2m,Edges,hsv,2*pi),1);
% % image(h(1)/2:h(1):1,h(2)/2:h(2):1,I);
% axis image; axis xy;axis off;
% set(gca,'YDir','reverse');
%
% subplot(234)
% hold on
% imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(dM1m,Edges,hsv,15))
% plot(corners{1}(:,1),corners{1}(:,2),'k+','LineWidth',2)
% plot(corners{2}(:,1),corners{2}(:,2),'k+','LineWidth',2)
% plot(corners{3}(:,1),corners{3}(:,2),'k+','LineWidth',2)
% plot(corners{4}(:,1),corners{4}(:,2),'k+','LineWidth',2)
% % I=flipdim(mergedImage(dM2m,Edges,hsv,2*pi),1);
% % image(h(1)/2:h(1):1,h(2)/2:h(2):1,I);
% axis image; axis xy;axis off;
% set(gca,'YDir','reverse');
%
% subplot(232)
% hold on
% imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(M2m,dEdges,hsv,2*pi))
% % I=flipdim(mergedImage(dM2m,Edges,hsv,2*pi),1);
% % image(h(1)/2:h(1):1,h(2)/2:h(2):1,I);
% axis image; axis xy;axis off;
% set(gca,'YDir','reverse');
%
% subplot(235)
% hold on
% imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(dM2m,Edges,hsv,2*pi))
% plot(corners{1}(:,1),corners{1}(:,2),'k+','LineWidth',2)
% plot(corners{2}(:,1),corners{2}(:,2),'k+','LineWidth',2)
% plot(corners{3}(:,1),corners{3}(:,2),'k+','LineWidth',2)
% plot(corners{4}(:,1),corners{4}(:,2),'k+','LineWidth',2)
% % I=flipdim(mergedImage(dM2m,Edges,hsv,2*pi),1);
% % image(h(1)/2:h(1):1,h(2)/2:h(2):1,I);
% axis image; axis xy;axis off;
% set(gca,'YDir','reverse');
%
% tmp=cool_springCmap;
% subplot(233)
% hold on
% imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(M3m,dEdges,tmp(129:224,:),15))
% % I=flipdim(mergedImage(dM2m,Edges,hsv,2*pi),1);
% % image(h(1)/2:h(1):1,h(2)/2:h(2):1,I);
% axis image; axis xy;axis off;
% set(gca,'YDir','reverse');
%
% subplot(236)
% hold on
% imagesc(1:(images_size-1)/(2^maxLevel-1):images_size,1:(images_size-1)/(2^maxLevel-1):images_size,mergedImage(dM3m,Edges,tmp(129:224,:),15))
% plot(corners{1}(:,1),corners{1}(:,2),'k+','LineWidth',2)
% plot(corners{2}(:,1),corners{2}(:,2),'k+','LineWidth',2)
% plot(corners{3}(:,1),corners{3}(:,2),'k+','LineWidth',2)
% plot(corners{4}(:,1),corners{4}(:,2),'k+','LineWidth',2)
% % I=flipdim(mergedImage(dM2m,Edges,hsv,2*pi),1);
% % image(h(1)/2:h(1):1,h(2)/2:h(2):1,I);
% axis image; axis xy;axis off;
% set(gca,'YDir','reverse');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
hydro = struct();

hydro = readWAMIT(hydro,'cubes.out',[]);
hydro = radiationIRF(hydro,20,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,200,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [xi,w]=thirdOrderCrossPolyCubPoints(numDim,algorithm)
%%THIRDORDERCROSSPOLYCUBPOINTS Generate third order cubature points for
%               integration over the region sum(abs(x))<=1 with weighting
%               function w(x)=1. In two dimensions, this is a square, in
%               three an octohedron, in n a cross polytope (or
%               n-dimensional octahedron).
%
%INPUTS: numDim An integer specifying the dimensionality of the points to
%               be generated.
%     algorithm A value indicating which algorithm should be used. Possible
%               values are:
%               0 (The default if omitted or an empty matrix is passed)
%                 Formula G_n 3-1 in [1], pg. 304, 2*numDim points.
%               1 Formula G_n 3-2 in [1], pg. 304, 2^numDim points,
%                 numDim<4.
%
%OUTPUTS: xi A numDim X numCubaturePoints matrix containing the cubature
%            points. (Each "point" is a vector)
%          w A numCubaturePoints X 1 vector of the weights associated with
%            the cubature points.
%
%REFERENCES:
%[1] A.H. Stroud, Approximate Calculation of Multiple Integrals. Cliffs,
%    NJ: Prentice-Hall, Inc., 1971.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

if(nargin<2||isempty(algorithm))
   algorithm=0;
end

switch(algorithm)
    case 0%G_n 3-1 in [1], pg. 304, 2*numDim points.
        V=exp(numDim*log(2)-gammaln(numDim+1));
        n=numDim;

        r=sqrt(2*n/((n+1)*(n+2)));

        w=V/(2*n)*ones(2*n,1);
        xi=fullSymPerms([r;zeros(numDim-1,1)]);
    case 1%G_n 3-2 in [1], pg. 304, 2^numDim points, numDim<4.
        if(numDim>=4)
           error('This algorithm requires numDim<4')
        end

        V=exp(numDim*log(2)-gammaln(numDim+1));
        n=numDim;

        r=sqrt(2/((n+1)*(n+2)));

        w=2^(-n)*V*ones(2^n,1);
        xi=PMCombos(r*ones(numDim,1));
    otherwise
        error('Unknown algorithm specified');
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef testEnergyEquivalenceVoigtAndTensorNotationForIAniTensor ...
        < testEnergyEquivalenceVoigtAndTensorNotation

    methods (Access = public)

        function obj = testEnergyEquivalenceVoigtAndTensorNotationForIAniTensor()
            obj@testEnergyEquivalenceVoigtAndTensorNotation();
        end

    end

    methods

        function generateFourthOrderTensor(obj)
            obj.Ch = Stiffness3DTensor();
            obj.Ch.createRandomTensor();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function st_io_test ( )

%*****************************************************************************80
%
%% ST_IO_TEST tests the ST_IO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_IO_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ST_IO library.\n' );

  st_io_test01 ( );
  st_io_test02 ( );
  st_io_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_IO_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = maxk(varargin)
% function res = MAXK(list, k)
%
% If LIST is a vector, MAXK returns in RES the K largest elements of LIST
%   RES is sorted in descending order
% [res loc] = MAXK(...)
%   Location of the largest elements: RES=LIST(LOC)
% If list is a matrix, MAXK operates along the first dimension
% Use MAXK(..., dim) to operate along the dimension dim
%     MAXK(..., dim, 'sorting', false) to disable the post-sorting step
%                                      (true by default)
%
% Author Bruno Luong <brunoluong@yahoo.com>
% Contributor: Matt Fig
% Last update: 07/April/2009
%              10/Jan/2010: possibility to disable post-sorting step

nout=cell(1,max(1,nargout));
[nout{:}] = minmaxk(@maxkmex, varargin{:});
varargout=nout;

author: andrewssobral

repo: lrslibrary


text:
%RNA production medium test
%
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 9/10/2010
classdef RNA_Test < TestCase
    properties
        simulation
    end

    methods
        function this = RNA_Test(name)
            this = this@TestCase(name);
        end

        function setUp(this)
            sim = edu.stanford.covert.cell.sim.SimulationFixture.load([], {
                'Metabolism'
                'Transcription'
                'RNAProcessing'
                'RNAModification'
                'tRNAAminoacylation'
                'RNADecay'
                });

            this.simulation = sim;
        end

        function testRnaProduction(this)
            import edu.stanford.covert.util.ConstantUtil;

            %% references
            sim = this.simulation;
            %this.seedSimulation(sim, 10);
            %sim.applyOptions('verbosity', 1);

            g = sim.gene;
            mass = sim.state('Mass');
            m = sim.state('Metabolite');
            rna = sim.state('Rna');
            time = sim.state('Time');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');
            transcription = sim.process('Transcription');

            %% constants
            rnaMWs = rna.molecularWeights(rna.aminoacylatedIndexs);

            %% turn off decay
            rna.decayRates(setdiff(1:end, rna.intergenicIndexs)) = 0;

            %% keep track of initial state
            initRNAs = ...
                + rna.counts(rna.matureIndexs, sim.compartment.cytosolIndexs) ...
                + rna.counts(rna.aminoacylatedIndexs, sim.compartment.cytosolIndexs);
            initRNAs(setdiff(1:end, rna.matureMRNAIndexs)) = ...
                initRNAs(setdiff(1:end, rna.matureMRNAIndexs)) + ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * sum(pc.counts(pc.matureIndexs, :) + pc.counts(pc.boundIndexs, :), 2);

            initMonomers = ...
                + pm.counts(pm.matureIndexs, :) ...
                + pm.counts(pm.boundIndexs, :);
            initComplexs = ...
                + pc.counts(pc.matureIndexs, :) ...
                + pc.counts(pc.boundIndexs, :);

            %% simulate
            iterMax = 1000;

            ntps = zeros(4, iterMax);
            if sim.verbosity > 0
                fprintf('%5s %6s %6s %6s %6s\n', 'Iter ', ' ATP  ', ' CTP  ', ' GTP  ', ' UTP  ');
                fprintf('%5s %6s %6s %6s %6s\n', '=====', '======', '======', '======', '======');
            end
            for i = 1:iterMax
                if sim.verbosity > 0 && mod(i, 100) == 1
                    fprintf('%5d %6d %6d %6d %6d\n', i, ...
                        transcription.substrates(transcription.substrateIndexs_ntp(1)), ...
                        transcription.substrates(transcription.substrateIndexs_ntp(2)), ...
                        transcription.substrates(transcription.substrateIndexs_ntp(3)), ...
                        transcription.substrates(transcription.substrateIndexs_ntp(4)));
                end

                %mock protein synthesis
                pm.counts(pm.matureIndexs, :) = ...
                    + pm.counts(pm.matureIndexs, :) ...
                    + pm.randStream.stochasticRound(initMonomers * log(2) / time.cellCycleLength * exp(i * log(2) / time.cellCycleLength));
                pc.counts(pc.matureIndexs, :) = ...
                    + pc.counts(pc.matureIndexs, :) ....
                    + pc.randStream.stochasticRound(initComplexs * log(2) / time.cellCycleLength * exp(i * log(2) / time.cellCycleLength));

                %simulate
                sim.evolveState();

                ntps(:, i) = transcription.substrates(transcription.substrateIndexs_ntp);
            end

            newRNAs = ...
                - initRNAs ...
                + rna.counts(rna.matureIndexs, sim.compartment.cytosolIndexs) ...
                + rna.counts(rna.aminoacylatedIndexs, sim.compartment.cytosolIndexs);
            newRNAs(setdiff(1:end, rna.matureMRNAIndexs)) = ...
                newRNAs(setdiff(1:end, rna.matureMRNAIndexs)) + ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * sum(initComplexs, 2);

            %% assertions

            %NTPs are limiting
            assertIn(min(min(ntps(:, end-100+1:end))), [0 max(2000, min(ntps(:, 1)))]);

            %NTPs are mostly used -- not too many are free
            assertIn(...
                (m.counts(m.ntpIndexs([2 4]), sim.compartment.cytosolIndexs)' * m.molecularWeights(m.ntpIndexs([2 4]))) / ...
                (rna.counts(:, sim.compartment.cytosolIndexs)' * rna.molecularWeights), ...
                [0 0.03]);

            %mature/aminoacylated RNA counts increased
            assertAllEqual(true, newRNAs >= 0);
            assertElementsAlmostEqual(newRNAs' * rnaMWs / ConstantUtil.nAvogadro, ...
                (exp(iterMax * log(2) / time.cellCycleLength) - 1) * ...
                mass.initialFractionNTPsInRNAs * mass.cellInitialDryWeight * mass.dryWeightFractionRNA, ...
                'relative', 2.5e-1, 0);

            assertIn(sum(sum(rna.counts(rna.nascentIndexs, :))), [0 20]);
            assertIn(sum(sum(rna.counts(rna.processedIndexs, :))), [0 40]);
            assertIn(sum(sum(rna.counts(rna.intergenicIndexs, :))), [0 10]);
            assertAllEqual(0, rna.counts(rna.boundIndexs, :));
            assertAllEqual(0, rna.counts(rna.damagedIndexs, :));
            assertAllEqual(0, rna.counts(rna.misfoldedIndexs, :));

            %RNA production matches expectations
            rates = rna.nascentRNAMatureRNAComposition * transcription.computeRNAPolymeraseTUBindingProbabilities();
            rates = rates(:, 1);
            assertIn(corr(rates, newRNAs), [0.50 1]);
            assertIn(180 / pi * acos(rates' * newRNAs / ...
                (sqrt(rates' * rates) * sqrt(newRNAs' * newRNAs))), ...
                [0 60]);
        end

        function testRnaProductionAndDecay(this)
            import edu.stanford.covert.util.ConstantUtil;

            %% references
            sim = this.simulation;
            %sim.applyOptions('verbosity', 1);
            this.seedSimulation(sim, 1000);

            g = sim.gene;
            comp = sim.compartment;
            mass = sim.state('Mass');
            m = sim.state('Metabolite');
            rna = sim.state('Rna');
            time = sim.state('Time');
            pm = sim.state('ProteinMonomer');
            pc = sim.state('ProteinComplex');
            rnaPol = sim.state('RNAPolymerase');
            transcript = sim.state('Transcript');
            transcription = sim.process('Transcription');

            %% keep track of initial state
            initRNAs = ...
                + rna.counts(rna.matureIndexs, comp.cytosolIndexs) ...
                + rna.counts(rna.aminoacylatedIndexs, comp.cytosolIndexs);
            initRNAs(setdiff(1:end, rna.matureMRNAIndexs)) = ...
                initRNAs(setdiff(1:end, rna.matureMRNAIndexs)) + ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * ...
                sum(pc.counts(pc.matureIndexs, :) + pc.counts(pc.boundIndexs, :), 2);
            initMonomers = ...
                + pm.counts(pm.matureIndexs, :) ...
                + pm.counts(pm.boundIndexs, :);
            initComplexs = ...
                + pc.counts(pc.matureIndexs, :) ...
                + pc.counts(pc.boundIndexs, :);

            %% simulate
            iterMax = 30000;

            rnaExp = zeros(size(initRNAs));
            rnaProd = zeros(size(rna.nascentIndexs));
            oldNascentRNAs = transcription.RNAs;
            rnaPolStateOcc = zeros(4, iterMax);
            ntps = zeros(4, iterMax);
            ndps = zeros(4, iterMax);
            nmps = zeros(4, iterMax);
            phosphates = zeros(3, iterMax);
            boundTUs = zeros(numel(rna.nascentIndexs), iterMax);
            rnaDryWt = zeros(1, iterMax);
            transcriptDryWt = zeros(1, iterMax);
            nucDryWt = zeros(1, iterMax);
            if sim.verbosity > 0
                fprintf('%5s %6s %6s %6s %6s\n', 'Iter ', ' ATP  ', ' CTP  ', ' GTP  ', ' UTP  ');
                fprintf('%5s %6s %6s %6s %6s\n', '=====', '======', '======', '======', '======');
            end
            for i = 1:iterMax
                if sim.verbosity > 0 && mod(i, 100) == 1
                    fprintf('%5d %6d %6d %6d %6d\n', i, ...
                        transcription.substrates(transcription.substrateIndexs_ntp(1)), ...
                        transcription.substrates(transcription.substrateIndexs_ntp(2)), ...
                        transcription.substrates(transcription.substrateIndexs_ntp(3)), ...
                        transcription.substrates(transcription.substrateIndexs_ntp(4)));
                end

                %mock protein synthesis
                pm.counts(pm.matureIndexs, :) = ...
                    + pm.counts(pm.matureIndexs, :) ...
                    + pm.randStream.stochasticRound(initMonomers * log(2) / time.cellCycleLength * exp(i * log(2) / time.cellCycleLength));
                pc.counts(pc.matureIndexs, :) = ...
                    + pc.counts(pc.matureIndexs, :) + ...
                    + pc.randStream.stochasticRound(initComplexs * log(2) / time.cellCycleLength * exp(i * log(2) / time.cellCycleLength));

                %simulate
                sim.evolveState();

                %store data
                rnaExp = rnaExp ...
                    + rna.counts(rna.matureIndexs, comp.cytosolIndexs) ...
                    + rna.counts(rna.aminoacylatedIndexs, comp.cytosolIndexs);
                rnaExp(setdiff(1:end, rna.matureMRNAIndexs)) = ...
                    rnaExp(setdiff(1:end, rna.matureMRNAIndexs)) + ...
                    sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * sum(initComplexs, 2);
                rnaProd = rnaProd + max(0, transcription.RNAs - oldNascentRNAs);
                oldNascentRNAs = transcription.RNAs;
                rnaPolStateOcc(:, i) = [
                    rnaPol.nActive
                    rnaPol.nSpecificallyBound
                    rnaPol.nNonSpecificallyBound
                    rnaPol.nFree
                    ];
                ntps(:, i) = transcription.substrates(transcription.substrateIndexs_ntp);
                ndps(:, i) = m.counts(m.ndpIndexs, comp.cytosolIndexs);
                nmps(:, i) = m.counts(m.nmpIndexs, comp.cytosolIndexs);
                phosphates(:, i) = m.counts([m.diphosphateIndexs; m.phosphateIndexs; m.hydrogenIndexs], comp.cytosolIndexs);
                boundTUs(:, i) = histc(transcript.boundTranscriptionUnits(rnaPol.states >= rnaPol.activelyTranscribingValue), (1:numel(rna.nascentIndexs))');

                rnaDryWt(:, i) = rna.dryWeight(1);
                transcriptDryWt(:, i) = transcript.dryWeight;
                nucDryWt(:, i) = (...
                    + ntps(:, i)' * m.molecularWeights(m.ntpIndexs) ...
                    + ndps(:, i)' * m.molecularWeights(m.ndpIndexs) ...
                    + nmps(:, i)' * m.molecularWeights(m.nmpIndexs)) / ConstantUtil.nAvogadro;
            end

            %% assertions

            %NTPs are limiting
            assertIn(min(min(ntps(:, end-100+1:end))), [0 exp(log(2) * iterMax / time.cellCycleLength) * min(ntps(:, 1))]);
            assertIn(min(min(ndps(:, end-100+1:end))), [0 exp(log(2) * iterMax / time.cellCycleLength) * min(ndps(:, 1))]);
            assertIn(min(min(nmps(:, end-100+1:end))), [0 exp(log(2) * iterMax / time.cellCycleLength) * min(nmps(:, 1))]);
            assertIn(min(transcriptDryWt(:, end-100+1:end)), [0 2.5 * exp(log(2) * iterMax / time.cellCycleLength) * transcriptDryWt(1)]);

            %NTPs are mostly used -- not too many are free
            assertIn(...
                (m.counts(m.ntpIndexs([2 4]), comp.cytosolIndexs)' * m.molecularWeights(m.ntpIndexs([2 4]))) / ...
                ((rna.dryWeight(1) + transcript.dryWeight) * edu.stanford.covert.util.ConstantUtil.nAvogadro), ...
                [0 0.05]);

            rnas = ...
                + rna.counts(rna.matureIndexs, comp.cytosolIndexs) ...
                + rna.counts(rna.aminoacylatedIndexs, comp.cytosolIndexs);
            rnas(setdiff(1:end, rna.matureMRNAIndexs)) = ...
                rnas(setdiff(1:end, rna.matureMRNAIndexs)) + ...
                sum(pc.proteinComplexComposition(setdiff(1:end, g.mRNAIndexs), :, :), 3) * sum(initComplexs, 2);
            assertElementsAlmostEqual(rnas' * rna.molecularWeights(rna.aminoacylatedIndexs) / edu.stanford.covert.util.ConstantUtil.nAvogadro, ...
                exp(iterMax * log(2) / time.cellCycleLength) * mass.cellInitialDryWeight * mass.dryWeightFractionRNA, ...
                'relative', 2e-1, 0);

            %total NTP incorporation similar to production by metabolism
            assertElementsAlmostEqual(...
                (edu.stanford.covert.util.ComputationUtil.invertCompositionMatrix(rna.nascentRNAMatureRNAComposition) * ...
                (mass.cellInitialDryWeight * mass.dryWeightFractionRNA * edu.stanford.covert.util.ConstantUtil.nAvogadro / ...
                (rna.expression(rna.matureIndexs)' * rna.molecularWeights(rna.matureIndexs)) * rna.expression(rna.matureIndexs) .* ...
                (log(2) / time.cellCycleLength + rna.decayRates(rna.matureIndexs))))' * rna.lengths(rna.nascentIndexs) * ...
                time.cellCycleLength / log(2) * (exp(log(2) / time.cellCycleLength * iterMax)-1), ...
                rnaProd' * rna.lengths(rna.nascentIndexs), ...
                'relative', 0.40);

            %RNA expression matches expectations
            assertIn(corr(rna.expression(rna.matureIndexs), rnaExp), [0.9 1]);
            assertIn(180 / pi * acos(rna.expression(rna.matureIndexs)' * rnaExp / ...
                (sqrt(rna.expression(rna.matureIndexs)' * rna.expression(rna.matureIndexs)) * sqrt(rnaExp' * rnaExp))), ...
                [0 20]);

            %RNA polymerase occupancy matches expectations
            assertElementsAlmostEqual(rnaPol.stateExpectations, ...
                sum(rnaPolStateOcc, 2) / sum(rnaPolStateOcc(:)), ...
                'relative', 0.75, 0.20); %todo: tighten

            %RNAs matured
            assertIn(sum(sum(rna.counts(rna.nascentIndexs, :))), [0 20]);
            assertIn(sum(sum(rna.counts(rna.processedIndexs, :))), [0 30]);
            assertIn(sum(sum(rna.counts(rna.intergenicIndexs, :))), [0 10]);
            assertAllEqual(0, rna.counts(rna.boundIndexs, :));
            assertAllEqual(0, rna.counts(rna.damagedIndexs, :));
            assertAllEqual(0, rna.counts(rna.misfoldedIndexs, :));
        end

        function sim = seedSimulation(~, sim, seed)
            if sim.verbosity > 0
                fprintf('Seed: %d\n', seed)
            end

            sim.applyOptions('seed', seed);
            sim.seedRandStream();

            for i = 1:numel(sim.states)
                o = sim.states{i};
                o.seed = seed;
                o.seedRandStream();
            end

            for i = 1:numel(sim.processes)
                o = sim.processes{i};
                o.seed = seed;
                o.seedRandStream();
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function l=v_gmmlpdf(varargin)
%V_GMMLPDF obsolete function - please use GAUSSMIXP instead
l=v_gaussmixp(varargin{1:nargin});

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% Copyright (C) 2014-2016, 2018-2019, 2022 Colin B. Macdonald
%% Copyright (C) 2022 Chris Gorman
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym heaviside (@var{x})
%% @defmethodx @@sym heaviside (@var{x}, @var{zero_value})
%% Symbolic Heaviside step function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.4")'))
%% y = heaviside (x)
%%   @result{} y = (sym) θ(x)
%% @end group
%% @end example
%%
%% There are various conventions for 'heaviside(sym(0))'; this function
%% returns the midpoint by default.
%% @example
%% @group
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.8")'))
%% heaviside(sym(0))
%%   @result{} (sym) 1/2
%% @end group
%% @end example
%%
%% The optional second argument overrides the default.  For example
%% to make the function right-continuous,
%% @example
%% @group
%% heaviside(0, sym(1))
%%   @result{} (sym) 1
%% @end group
%% @end example
%% or left-continuous.
%% @example
%% @group
%% heaviside(0, sym(0))
%%   @result{} (sym) 0
%% @end group
%% @end example
%%
%% If passed a matrix, the heaviside function is computed for each
%% element.
%% @example
%% @group
%% @c doctest: +SKIP_UNLESS(pycall_sympy__ ('return Version(spver) > Version("1.8")'))
%% heaviside([sym(-1) sym(0) sym(1) sym(x)])
%%   @result{} (sym) [0  1/2  1  θ(x)]  (1×4 matrix)
%% @end group
%% @end example
%% @seealso{heaviside, @@sym/dirac}
%% @end defmethod

function y = heaviside(x, h0)
  if (nargin == 1)
    y = elementwise_op ('Heaviside', x);
  elseif (nargin == 2)
    y = elementwise_op ('Heaviside', sym(x), sym(h0));
  else
    print_usage ();
  end
end

%!error heaviside (sym(1), 2, 3)

%!assert (isequal (heaviside (sym(1)), sym(1)))
%!assert (isequal (heaviside (-sym(1)), sym(0)))

%!assert (double (heaviside (1)), heaviside (1))

%!test
%! D = [1 -1; -10 20];
%! A = sym(D);
%! assert (double (heaviside (A)), heaviside (D))

%!test
%! H0 = sym([1 -2 0; 3 0 pi]);
%! A = heaviside (sym(0), H0);
%! assert (isequal (A, H0))

%!test
%! A = heaviside ([-1 0 1], sym(1)/2);
%! assert (isequal (A, [0 sym(1)/2 1]))

%!test
%! A = heaviside ([-1 0 1], sym(1)/2);
%! assert (isequal (A, [0 sym(1)/2 1]))

%!assert (isequaln (heaviside (sym(nan)), sym(nan)))

%!test
%! assert (isequaln (heaviside (sym(nan), sym(nan)), sym(nan)))
%! assert (isequaln (heaviside (0, sym(nan)), sym(nan)))
%! assert (isequaln (heaviside (2, sym(nan)), sym(1)))
%! assert (isequaln (heaviside (-2, sym(nan)), sym(0)))

%!test
%! % round trip
%! syms x
%! A = heaviside (1);
%! f = heaviside (x);
%! h = function_handle (f);
%! B = h (1);
%! assert (A, B, -eps)

%!test
%! % round trip
%! syms x h0
%! f = heaviside (x, h0);
%! h = function_handle (f, 'vars', {x h0});
%! A = heaviside (1, 1/2);
%! B = h (1, 1/2);
%! assert (A, B, -eps)
%! A = heaviside (0, 1/2);
%! B = h (0, 1/2);
%! assert (A, B, -eps)

repo: octsympy


text:
function Inds = sllabelinds(labels, labelset)
%SLLABELINDS Extract indices corresponding to specified labels
%
% $ Syntax $
%   - Inds = sllabelinds(labels, labelset)
%
% $ Arguments $
%   - labels:       The labels of samples
%   - labelset:     The set of labels whose indices to be extracted
%   - Inds:         The cell array of indices extracted for labelset
%
% $ Description $
%   - Inds = sllabelinds(labels, labelset) extracts the indices
%     corresponding to the labels specified in labelset. Suppose the
%     labelset is given by [l1, l2, ...], then Inds would be like
%     {[i11, i12, ...], [i21, i22, ...], ...}, where [i11, i12, ...] is
%     a row vector of indices corresponding to l1, so that
%     labels(i11) = labels(i12) = ... = l1.
%
% $ History $
%   - Created by Dahua Lin, on Aug 31, 2006
%

%% parse and verify input

if ~isvector(labels) || ~isnumeric(labels)
    error('sltoolbox:invalidarg', ...
        'labels should be a numeric vector');
end

if size(labels, 1) ~= 1
    labels = labels(:)';
end

%% re-arrange

[labels, si] = sort(labels, 2, 'ascend');
[nums, labelfound] = slcount(labels);
[sinds, einds] = slnums2bounds(nums);
[sfound, smap] = ismember(labelset, labelfound);

%% extract

c = length(labelset);
Inds = cell(1, c);
for i = 1 : c
    if sfound(i)
        mi = smap(i);
        curinds = si(sinds(mi):einds(mi));
        Inds{i} = curinds;
    else
        Inds{i} = [];
    end
end

author: lmthang


text:
function out = ind2char(i,s)
% calculates coordinate from one index
% i - index
% s - size of matrix

out = '(';
for d = 1:length(s)
	out = [out,int2str(mod(ceil(i/prod(s(1:d-1)))-1,s(d))+1),',']; %#ok<AGROW>
end
out(end) = ')';

author: mtex-toolbox

repo: mtex

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


text:
function digit = r8_digit ( x, idigit )

%*****************************************************************************80
%
%% R8_DIGIT returns a particular decimal digit of an R8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose NDIG-th decimal digit
%    is desired.  If X is zero, all digits will be returned as 0.
%
%    Input, integer IDIGIT, the position of the desired decimal digit.
%    A value of 1 means the leading digit, a value of 2 the second digit
%    and so on.
%
%    Output, integer DIGIT, the value of the IDIGIT-th decimal digit of X.
%
  if ( x == 0.0 )
    digit = 0;
    return
  end

  if ( idigit <= 0 )
    digit = 0;
    return
  end
%
%  Force X to lie between 1 and 10.
%
  x = abs ( x );

  while ( x < 1.0 )
    x = x * 10.0;
  end

  while ( 10.0 <= x )
    x = x / 10.0;
  end

  for i = 1 : idigit
    ival = floor ( x );
    x = ( x - ival ) * 10.0;
  end

  digit = ival;

  return
end


text:
% This is an approximation of random numbers wrt to wrapped Normal Distrib
%
% Any neater, faster, more precise and accurate generator is very welcome.
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  you can
%  redistribute it and/or modify it under the terms of the GNU Lesser General
%  Public License as published by the Free Software Foundation, either version 3
%  of the License, or (at your option) any later version. This file is
%  distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%  without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
%  PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
%  details.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function wgr = wrappednormrnd(mu, sigma, n, pregen)

    if nargin<4; pregen=false; end

    % The CDF is symetric around the mean. Thus, just need zero-mean random values
    % of unitary std and scale them wrt sigma.
    wgr = wrappednormrndunitystd(n, pregen);

    wgr = sigma.*wgr;

    wgr = wrap(wgr + mu);

return

function wgr = wrappednormrndunitystd(n, pregen)

    if nargin<2; pregen=0; end

    if pregen==1
        global wrappednormrnd;

        % If run for the first time, initialize the cumulative distrib fns
        if isempty(wrappednormrnd) || ~isfield(wrappednormrnd, 'icdf')
            fprintf('wrappednormrnd.m: Pre-compute the cumulative distribution function once for all ... ');

            xs = -pi:(2*pi)/1e6:pi;
            xs = xs(1:end-1);
            cdf = wrappednormcdf(xs, 0, 1);
            % Resample the inverse CDF with 1 million uniform steps
            % This also assumes that we don't need a thinner resolution of the distribution
            wrappednormrnd.icdf_ys = (0:1e-4:1);
            wrappednormrnd.icdf = interp1(cdf, xs, wrappednormrnd.icdf_ys, 'linear', 'extrap').';

            fprintf('done.\n');
            if nargin==0; return; end
        end

        % wgr = interp1(wrappednormrnd.icdf_ys, sigma*wrappednormrnd.icdf, rand(n,1), 'linear', 'extrap');
        % The following is a speed-up version of the above line.
        ri = (length(wrappednormrnd.icdf)-1)*rand(n,1)+1;
        rfi = floor(ri);
        rci = ceil(ri);
        wgr = wrappednormrnd.icdf(rfi);
        wgr = wgr + (ri-rfi).*(wrappednormrnd.icdf(rci)-wgr);

        % % Or a faster sampling (implies only 1e4 different random values can be generated).
        % idx = round(length(wrappednormrnd.icdf)*rand(n,1))+1;
        % idx = max(1,min(length(wrappednormrnd.icdf),idx));
        % wgr = wrappednormrnd.icdf(idx);
    else
        % The non-optimized random generator
        xs = -pi:(2*pi)/1e6:pi;
        xs = xs(1:end-1);
        cdf = wrappednormcdf(xs, 0, 1);

        wgr = interp1(cdf, xs, rand(n,1), 'linear', 'extrap');
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [areas, count] = calc_triangle_areas(verts, faces, space, focus)

if strcmp(space, 'object')
    tri_areas = cal_obj_area(verts, faces);
elseif strcmp(space, 'parameter')
    tri_areas = cal_par_area(verts, faces);
end

if strcmp(focus, 'vertex')
    areas = zeros(length(verts),1);
    count = zeros(length(verts),1);

% find all incident triangles upon each vertex and accumulate the areas of the triangles
    for i=1:size(faces,2)
        [u1, m1, n1] = unique(faces(:,i));
        areas(u1) = areas(u1) + accumarray(n1, tri_areas);
        count(u1) = count(u1) + accumarray(n1, ones(length(n1),1));
    end

elseif strcmp(focus, 'triangle')
    areas = tri_areas;
    count = ones(length(tri_areas),1);
end

return;

%
% calculate relative areas of triangles on object surface net
%

function obj_area = cal_obj_area(vertices,faces)

A = faces(:,1); B = faces(:,2); C = faces(:,3);
a = sqrt(sum(((vertices(A,:)-vertices(B,:)).^(2))'))';
b = sqrt(sum(((vertices(B,:)-vertices(C,:)).^(2))'))';
c = sqrt(sum(((vertices(C,:)-vertices(A,:)).^(2))'))';
s = (a+b+c)/2;
obj_area = sqrt(s.*(s-a).*(s-b).*(s-c));
obj_area = obj_area/sum(obj_area);

return;

%
% calculate relative areas of spherical triangles in parameter space
%

function par_area = cal_par_area(vs,faces)

angles = [];
for j = 1:3
    % note that the order of A B C is clockwise (see 08-22-02.htm notes)
    A = vs(faces(:,j),:);
    B = vs(faces(:,mod(j,3)+1),:);
    C = vs(faces(:,mod(j-2,3)+1),:);
    y = A(:,1).*B(:,2).*C(:,3) - A(:,1).*B(:,3).*C(:,2) + ...
        A(:,2).*B(:,3).*C(:,1) - A(:,2).*B(:,1).*C(:,3) + ...
        A(:,3).*B(:,1).*C(:,2) - A(:,3).*B(:,2).*C(:,1);
    x = B(:,1).*C(:,1) + B(:,2).*C(:,2) + B(:,3).*C(:,3) - ...
       (A(:,1).*C(:,1) + A(:,2).*C(:,2) + A(:,3).*C(:,3)).* ...
       (A(:,1).*B(:,1) + A(:,2).*B(:,2) + A(:,3).*B(:,3));
    angles(:,j) = atan2(y,x);
end
ind = find(angles<0);
angles(ind) = angles(ind) + 2*pi;
par_area = sum(angles')' - pi;
par_area = par_area/(4*pi);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [stat, cfg] = ft_statistics_crossvalidate(cfg, dat, design)

% FT_STATISTICS_CROSSVALIDATE performs cross-validation using a prespecified
% multivariate analysis given by cfg.mva
%
% Use as
%   stat = ft_timelockstatistics(cfg, data1, data2, data3, ...)
%   stat = ft_freqstatistics    (cfg, data1, data2, data3, ...)
%   stat = ft_sourcestatistics  (cfg, data1, data2, data3, ...)
%
% Options:
%   cfg.mva           = a multivariate analysis (default = {dml.standardizer dml.svm})
%   cfg.statistic     = a cell-array of statistics to report (default = {'accuracy' 'binomial'})
%   cfg.nfolds        = number of cross-validation folds (default = 5)
%   cfg.resample      = true/false; upsample less occurring classes during
%                       training and downsample often occurring classes
%                       during testing (default = false)
%
% Returns:
%   stat.statistic    = the statistics to report
%   stat.model        = the models associated with this multivariate analysis
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS, FT_SOURCESTATISTICS

% Copyright (c) 2007-2011, F.C. Donders Centre, Marcel van Gerven
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% do a sanity check on the input data
assert(isnumeric(dat),    'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');
assert(isnumeric(design), 'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');

cfg.mva       = ft_getopt(cfg, 'mva');
cfg.statistic = ft_getopt(cfg, 'statistic', {'accuracy', 'binomial'});
cfg.nfolds    = ft_getopt(cfg, 'nfolds',   5);
cfg.resample  = ft_getopt(cfg, 'resample', false);
cfg.cv        = ft_getopt(cfg, 'cv', []);
cfg.cv.type   = ft_getopt(cfg.cv, 'type', 'nfold');

% specify classification procedure or ensure it's the correct object
if isempty(cfg.mva)
  cfg.mva = dml.analysis({ dml.standardizer('verbose',true) ...
                           dml.svm('verbose',true)});
elseif ~isa(cfg.mva,'dml.analysis')
  cfg.mva = dml.analysis(cfg.mva);
end

cv_options = {'mva', cfg.mva, 'type', cfg.cv.type, 'resample', cfg.resample, 'compact', true, 'verbose', true};
if strcmp(cfg.cv.type, 'nfold')
  cv_options = cat(2, cv_options, {'folds', cfg.nfolds});
end
cv = dml.crossvalidator(cv_options{:});

if any(isinf(dat(:)))
  ft_warning('Inf encountered; replacing by zeros');
  dat(isinf(dat(:))) = 0;
end

if any(isnan(dat(:)))
  ft_warning('Nan encountered; replacing by zeros');
  dat(isnan(dat(:))) = 0;
end

% perform everything
cv = cv.train(dat',design');

% extract the statistic of interest
s = cv.statistic(cfg.statistic);
for i=1:length(cfg.statistic)
 stat.statistic.(cfg.statistic{i}) = s{i};
end

% get the model averaged over folds
stat.model = cv.model;

fn = fieldnames(stat.model{1});
if any(ismember(fn,  {'weights', 'primal'})),
  selfn = find(ismember(fn, {'weights', 'primal'}));

  % the mean subtraction is needed only once, but speeds up the covariance
  % computation
  dat = bsxfun(@minus, dat, nanmean(dat,2));
  dat_transp = dat.';
  for j=1:numel(selfn)
    % create the 'encoding' matrix from the weights, as per Haufe 2014.
    %covdat = cov(dat');
    for i=1:length(stat.model)
      i
      W = stat.model{i}.(fn{selfn});

      sW   = size(W);
      sdat = size(dat);
      if sW(2)==sdat(1) && sW(1)~=sdat(1)
        W = transpose(W);
      end

      M    = dat'*W;
      covM = cov(M);
      WcovM = (W/covM)./(size(dat,2)-1); % with the correction term for the covariance computation

      %stat.model{i}.(sprintf('%sinv',fn{selfn})) = covdat*W/covM;
      stat.model{i}.(sprintf('%sinv',fn{selfn})) = dat*(dat_transp*WcovM);

    end
  end
end
fn = fieldnames(stat.model{1}); % update the fieldnames, because some might have been added

fn = fieldnames(stat.model{1}); % may now also contain weightsinv
for i=1:length(stat.model)
  for k=1:length(fn)
    if numel(stat.model{i}.(fn{k}))==prod(cfg.dim)
      stat.model{i}.(fn{k}) = reshape(stat.model{i}.(fn{k}),cfg.dim);
    end
  end

end

% required
stat.trial = [];

% add some stuff to the cfg
cfg.cv = cv;

author: spm


text:
function [MASKSTATS] = canlab_glm_maskstats(DIRS,MASK,varargin)
% Returns a MASKSTATS structure containing data from robfitdir's input subject level
%    SPM analyses.
%
% :Usage:
% ::
%
%     MASKSTATS = canlab_glm_maskstats(robfitdir, mask, [options])
%
% output structure:
% MASKSTATS
%
%   **COV:**
%        subject x covariate matrix (EXPT.cov from robfit design)
%
%   **COVNAME:**
%        names of covariates in COV
%
%   **MASK:**
%        array of structs (1 per mask)
%
%   **MASKFILE:**
%        the filename of the mask used
%
%   **SUB:**
%        struct containing data from subject level images
%
%   **CON:**
%        struct contains data from contrast images
%
%   **NAME:**
%        cell array (1 cell per contrast) of contrast names
%
%   **IMGFILES:**
%        cell array (1 cell per contrast) of character arrays of
%        contrast image filenames
%
%   **(MEASURE):**
%        subject X contrast matrix of measures (see canlab_maskstats)
%
%   **COVxCON.(MEASURE):**
%        arrays of covariate matrix X contrast means correlation
%        results (RHO and P, see help corr())
%
%   **GRP:**
%        struct containing data from group level images
%
%   **BETA:**
%        struct array (1 struct per group level regressor) of data
%        from beta images
%
%   **NAME:**
%        name of group level regressor
%
%   **IMGFILES:**
%        cell array (1 cell per robust directory) of beta image files
%
%   **(MEASURE):**
%        row vector (1 value per robust directory) of measures
%        (see canlab_maskstats)
%
% The following plots are saved in each mask's directory in the plots directory:
%   - contrast means by contrast (means are lines across subjects on x axis)
%   - contrast means by subject (means are dots, lined up along the x axis by contrast)
%   - group level betas (bar plot with group level regressors grouped by
%     subject level contrasts)
%
% If there's more than one regressor in the group level model, for each regressor:
%   - scatter plot of subject level contrast means against group level regressor
%
%
% :Arguments:
%
%   **robfitdir:**
%        a directory in which robfit was run
%        contains robfit directories (e.g., robust0001)
%        preferably contains EXPT.mat or EXPTm.mat
%
%   **mask:**
%        a filename or cell array of filenames of masks to apply to data
%
% :Options:
%
%   **MEASURE OPTIONS:**
%        (see canlab_maskstats) (DEFAULT: mean (within mask's non-zero voxels))
%
%   **'cons', connums:**
%        (vector of contrast numbers)
%        only include data from contrasts as specified by numbers
%
%   **'cons', conname(s):**
%        (string or cell array of strings)
%        only include data from contrasts as specified by name
%
%   **'plots':**
%        make plots
%
%   **'od', dir:**
%        will save plots in dir (DEFAULT: robfitdir/stats_scatterplots)
%
% ..
%    Programmers' notes:
%    to write:
%    save data option
%    save as csv option
%    option to include subject-level beta data
%    option to break input mask into regions?
%    within subjects error bars
%    grouped bar plots?
%    save volInfo (for mask? betas? cons? from fmri_data or spm?)
% ..

DO_PLOTS = false; % SETUP
% set defaults
OP = {};
ALLOPS = false;

applymaskopts = '';

% modify arguments as needed
if ~iscell(MASK), MASK = {MASK}; end

% error checking
if ~exist('DIRS','var') || isempty(DIRS)
    error('No robftidir specified.')
end
if ~exist('MASK','var') || isempty(MASK)
    error('No MASK(s) specified.')
end

% this gets checked/errored in canlab_maskstats
% for i=1:numel(MASK)
%     if ~exist(MASK{i},'file') && ~any(~cellfun('isempty',regexp({'nps','nps_thresh','nps_thresh_smooth'},['^' MASK{i} '$'])))
%         error('No such file: %s',MASK{i})
%     end
% end

% parse varargin
i=1;
while i<=numel(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case {'c' 'cons'}
                i=i+1;
                if ischar(varargin{i})
                    included_connames = varargin(i);
                elseif iscellstr(varargin{i})
                    included_connames = varargin{i};
                elseif isnumeric(varargin{i})
                    included_connums = varargin{i};
                else
                    error('argument to ''c'' must be number, string, or cell string')
                end
            case 'od'
                i=i+1;
                plotdir = varargin{i};
            case 'plots'
                DO_PLOTS = true;
            case {'binmean' 'mean' 'std' 'nonbinmean' 'nonbinstd' 'norm' ...
                  'dot_product' 'centered_dot_product' 'cosine_similarity' 'correlation'}
                OP{end+1} = varargin{i}; %#ok
            case 'all'
                ALLOPS = true;
            case 'l1norm'
                applymaskopts = [applymaskopts ', ''l1norm''']; %#ok
            case 'l2norm'
                applymaskopts = [applymaskopts ', ''l2norm''']; %#ok
            otherwise
                error(['UNRECOGNIZED OPTION: ' varargin{i}])
        end
    elseif iscellstr(varargin{i})
        for j=1:numel(varargin{i})
            switch varargin{i}{j}
                case {'binmean' 'mean' 'std' 'nonbinmean' 'nonbinstd' 'norm' ...
                        'dot_product' 'centered_dot_product' 'cosine_similarity' 'correlation'}
                    OP{end+1} = varargin{i}{j}; %#ok
                case 'all'
                    ALLOPS = true;
                otherwise
                    error('Unrecognized argument %s',varargin{i}{j})
            end
        end
    else
        disp(varargin{i})
        error('Above option UNRECOGNIZED')
    end
    i=i+1;
end

if ALLOPS, OP = {'cosine_similarity' 'dot_product' 'correlation' 'mean' 'std' 'centered_dot_product'}; end

if isempty(OP), OP = {'mean'}; end

%% PREP
if ischar(DIRS)
    DO_GRPLEV = true;
    BASEDOUT = DIRS;

    robdirs = filenames(fullfile(DIRS,'robust[0-9][0-9][0-9][0-9]'),'absolute');
    if isempty(robdirs)
        error('No such directories: %s/robust[0-9][0-9][0-9][0-9]',DIRS)
    end
    fprintf('... LOADING SETUP files from %d robust* directories in %s\n',numel(robdirs),DIRS)
    SETUPS = cell(numel(robdirs),1); allconnames = SETUPS; allconfiles = SETUPS;
    for r = 1:numel(robdirs)
        load(fullfile(robdirs{r},'SETUP'));
        SETUPS{r} = SETUP;
        allconnames{r} = SETUP.name;
        allconfiles{r} = SETUP.files;
        clear SETUP
    end

    % load EXPT
    if exist(fullfile(DIRS,'EXPTm.mat'),'file')
        load(fullfile(DIRS,'EXPTm'));
        EXPT = EXPTm;
        clear EXPTm;
    elseif exist(fullfile(DIRS,'EXPT.mat'),'file')
        load(fullfile(DIRS,'EXPT'));
    end

     % get covariates, if any
    if size(SETUPS{1}.X,2)>1
        fprintf('... GETTING between-subjects covariate(s)\n')
        cov = SETUPS{1}.X(:,2:end);

        covname = {};
        if isfield(EXPT,'covnames')
            covname = EXPT.covnames(2:end);
        end
    else
        cov = {};
        covname = '';
    end
    MASKSTATS.cov = cov;
    MASKSTATS.covname = covname;
elseif iscellstr(DIRS)
    DO_GRPLEV = false;
    BASEDOUT = pwd;

    for i = 1:numel(DIRS)
        clear SPM spmmat
        spmmat = fullfile(DIRS{i},'SPM.mat');
        if isempty(spmmat)
            error('No such file: %s',spmmat)
        end
        load(spmmat);
        tmpconnames{i} = cellstr(char(SPM.xCon.name)); %#ok
        ntmpconnames(i) = size(tmpconnames{i},1); %#ok
        allconfiles{i} = filenames(fullfile(DIRS{i}),'con_*.img'); %#ok
    end

    if size(unique(ntmpconnames))~=1
        error('Not all subject levels have same contrast file names.')
    end
    allconnames = tmpconnames{1};
    clear SPM spmmat tmpconnames
end

if DO_PLOTS && ~exist('plotdir','var')
    plotdir = fullfile(BASEDOUT,'maskstats_plots');
end

%% MAIN
% for r=1:numel(MASK)
%     clear mask
%
%     mask = MASK{r};
%     [ignore maskname] = fileparts(mask);
%     fprintf('... GETTING data for mask: %s\n',mask)

%% SUBJECT LEVEL STAT MAPS (cons)
fprintf('SUBJECT LEVELS\n')

clear lconname lconfiles lcondata

% get filenames
fprintf('Subject level contrasts:\n')
n=0;
lconfiles = {};
if exist('included_connums','var')
    connums = included_connums;
else
    connums = 1:numel(allconfiles);
end
for i = connums
    thisconname = deblank(allconnames{i});
    if exist('included_connames','var') && ~any(~cellfun('isempty',(regexp(thisconname,included_connames))))
        fprintf('EXCLUDED: %s\n',thisconname)
        continue
    else
        n=n+1;
        included_robustdirs(n) = i; %#ok
    end
    % name
    lconname{n} = thisconname; %#ok
    fprintf('%s\n',lconname{n})
    % files
    lconfiles = [lconfiles cellstr(SETUPS{i}.files)]; %#ok
end

% get data
lmaskstats = canlab_maskstats(MASK,lconfiles,OP);

% pack output struct
for m=1:numel(MASK)
    MASKSTATS.mask(m).mask = lmaskstats(m).maskfile;
    MASKSTATS.mask(m).sub.con = lmaskstats(m).stats;
    MASKSTATS.mask(m).sub.con.imgfiles = lmaskstats(m).imgfiles;
    MASKSTATS.mask(m).sub.con.name = lconname;
    if ~isempty(cov)
        for i=1:numel(OP)
            [rho pval] = corr(cov,MASKSTATS.mask(m).sub.con.(OP{i}));
            MASKSTATS.mask(m).sub.covXcon.(OP{i}).rho = rho;
            MASKSTATS.mask(m).sub.covXcon.(OP{i}).pval = pval;
        end
    end
end

% make plots
if DO_PLOTS
    fprintf('... MAKING plots\n')
    for m = 1:numel(MASK)
        [ignore maskname] = fileparts(MASK{m}); %#ok

        dout = fullfile(plotdir,maskname);
        if ~exist(dout,'dir'), mkdir(dout); end

        subticklabels = {};
        for i = 1:size(char(MASKSTATS.mask(m).sub.con.imgfiles{:,1}),1)
            subticklabels{i} = regexprep(deblank(MASKSTATS.mask(m).sub.con.imgfiles{i,1}),'.*/([^/]*)/[^/]*','$1'); %#ok
            subticklabels{i} = sprintf('%s (%d)',subticklabels{i},i); %#ok
        end

        for o=1:numel(OP)
            switch OP{o}
                case {'mean' 'dot_product' 'centered_dot_product' 'cosine_similarity' 'correlation'}
                    % line plot (one line for each contrast, across subjects)

                    figure; plot(MASKSTATS.mask(m).sub.con.(OP{o}));
                    if strcmp(OP{o},'mean')
                        title(sprintf('subject level contrasts mean within %s',maskname),'interpreter','none');
                    else
                        title(sprintf('subject level contrasts %ss with %s',OP{o},maskname),'interpreter','none');
                    end
                    ylabel(sprintf('beta contrast %s',OP{o}),'interpreter','none');
                    xaxis_labeling('subjects',subticklabels)
                    fout = fullfile(dout,sprintf('contrasts_by_subject__%ss.png',OP{o}));
                    fprintf('saving: %s\n',fout);
                    saveas(gcf,fout);
                    close gcf

                    % bar plot (mean across subjects per contrast)
                    evalc(['barplot_columns2(MASKSTATS.mask(m).sub.con.(OP{o}),maskname,'...
                        '''labels'',MASKSTATS.mask(m).sub.con.name,' ...
                        '''ylabel'',OP{o},''plabels'',''robust'',''ind'');']);
                    fout = fullfile(dout,sprintf('group_mean_%s.png',OP{o}));
                    fprintf('saving: %s\n',fout);
                    scn_export_papersetup(500);
                    saveas(gcf,fout);
                    close gcf

                    %%% covariate plots
                    if ~isempty(cov)
                        subdout = fullfile(dout,'scatterplots');
                        if ~exist(subdout,'dir'), mkdir(subdout); end
                        for i = 1:size(cov,2)
                            for j = 1:size(MASKSTATS.mask(m).sub.con.(OP{o}),2)
                                figure;
                                plot(cov(:,i),MASKSTATS.mask(m).sub.con.(OP{o})(:,j),'o');
                                lsline;
                                xlabel(covname{i},'interpreter','none');
                                ylabel(sprintf('beta (meant within %s)',maskname),'interpreter','none');
                                title(sprintf('%s correlation with %s (r=%.4f, p<%.4f)',...
                                    MASKSTATS.mask(m).sub.con.name{j},covname{i},...
                                    MASKSTATS.mask(m).sub.covXcon.(OP{o}).rho(i,j),...
                                    MASKSTATS.mask(m).sub.covXcon.(OP{o}).pval(i,j)),...
                                    'interpreter','none');

                                scon = regexprep(lconname{j},'[()]','');
                                scon = regexprep(scon,'[^A-Za-z0-9_+.-]','');

                                gcov = regexprep(covname{i},'[()]','');
                                gcov = regexprep(gcov,'[^A-Za-z0-9_+.-]','');

                                fout = fullfile(subdout,sprintf('%s_%ss__by__%s.png',scon,OP{o},gcov));
                                fprintf('saving: %s\n',fout)
                                saveas(gcf,fout);
                                close gcf
                            end
                        end
                    end
            end
        end
    end
end

%% GROUP LEVEL STAT MAPS (beta)
if DO_GRPLEV
    fprintf('GROUP LEVELS\n')
    nhcons = size(SETUPS{1}.X,2);
    hbetafiles = {};
    hbetaname = cell(nhcons,1);
    for b = 1:nhcons
        % get name
        if exist('EXPT','var') && isfield(EXPT,'covnames')
            hbetaname{b} = EXPT.covnames{b};
        else
            hbetaname{b} = sprintf('hbeta%04d',b);
        end

        % get list of files
        temphbetafiles = {};
        for i = included_robustdirs
            temphbetafiles{end+1,1} = fullfile(robdirs{i},sprintf('rob_beta_%04d.img',b)); %#ok
        end
        hbetafiles{b} = temphbetafiles; %#ok
    end
    hbetafiles = horzcat(hbetafiles{:});

    % extract data
    hmaskstats = canlab_maskstats(MASK,hbetafiles,OP);

    % pack output struct
    for b = 1:nhcons
        for m=1:numel(MASK)
            for o=1:numel(OP)
                MASKSTATS.mask(m).grp.beta(b).(OP{o}) = hmaskstats(m).stats.(OP{o})(:,b);
            end
            MASKSTATS.mask(m).grp.beta(b).imgfiles = hmaskstats(m).imgfiles;
            MASKSTATS.mask(m).grp.beta(b).name = hbetaname{b};
        end
    end

    if DO_PLOTS
        dout = fullfile(plotdir,maskname);
        if ~exist(dout,'dir'), mkdir(dout); end

        fprintf('... MAKING plots\n')
        for m = 1:numel(MASK)
            conticklabels = {};
            for i = 1:numel(MASKSTATS.mask(m).sub.con.name)
                conticklabels{i} = MASKSTATS.mask(m).sub.con.name{i}; %#ok
            end

            for o = 1:numel(OP)
                X = [];
                for b = 1:nhcons
                    X = [X MASKSTATS.mask(m).grp.beta(b).(OP{o})]; %#ok
                end

                figure; bar(X,'group');
                title(sprintf('group level betas %s with %s',OP{o},maskname),'Interpreter','none');
                ylabel(sprintf('beta %s',OP{o}),'Interpreter','none');
                legend(strvcat(MASKSTATS.mask(m).grp.beta.name)); %#ok
                xaxis_labeling('subject level contrasts',conticklabels);
                fout = fullfile(dout,sprintf('group_level_betas_%s.png',OP{o}));
                fprintf('saving: %s\n',fout);
                saveas(gcf,fout);
                close gcf
            end
        end
    end
end

%% clean up
close all

if DO_PLOTS
    % save MASKSTATS
    for m = 1:numel(MASK)
        [ignore maskname] = fileparts(MASK{m}); %#ok
        dout = fullfile(plotdir,maskname);
        maskstats = MASKSTATS.mask(m); %#ok
        save(fullfile(dout,'raw_data'),'maskstats');
    end
end

end

%% slightly edited code from mathworks.com (search "matlab angled tick labels")
function xaxis_labeling(xaxislabel,ticklabels)

% reduce size of axis to fit labels
pos = get(gca,'Position');
set(gca,'Position',[pos(1), .2, pos(3) .65])

% set tick locations
Xt = 1:numel(ticklabels);
Xl = [0 numel(ticklabels)+1];
set(gca,'XTick',Xt,'Xlim',Xl);

ax = axis; % Current axis limits
axis(axis); % Set the axis limit modes to manual
Yl = ax(3:4); % Y-axis limits

% Place the text labels
t = text(Xt,Yl(1)*ones(1,numel(Xt)),ticklabels,'Interpreter','none');
set(t,'HorizontalAlignment','right','VerticalAlignment','top', ...
    'Rotation',45);

% Remove the default labels
set(gca,'XTickLabel','')

% get the Extent of each text object
for i = 1:length(ticklabels)
    ext(i,:) = get(t(i),'Extent'); %#ok
end

% Determine lower point for alignment of text
LowYPoint = min(ext(:,2));

% Place the axis label
XMidPoint = Xl(1) + abs(diff(Xl))/2;
text(XMidPoint,LowYPoint,xaxislabel,'VerticalAlignment','top','HorizontalAlignment','center');

% adjust width
OP = get(gca,'OuterPosition');
set(gcf,'OuterPosition',[OP(1) OP(2) numel(ticklabels)*50 OP(3)])

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [V, converged, i] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt)
%GAUSSPF  Solves the power flow using a Gauss-Seidel method.
%   [V, CONVERGED, I] = GAUSSPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
%   solves for bus voltages given the full system admittance matrix (for
%   all buses), the complex bus power injection vector (for all buses),
%   the initial vector of complex bus voltages, and column vectors with
%   the lists of bus indices for the swing bus, PV buses, and PQ buses,
%   respectively. The bus voltage vector contains the set point for
%   generator (including ref bus) buses, and the reference angle of the
%   swing bus, as well as an initial guess for remaining magnitudes and
%   angles. MPOPT is a MATPOWER options struct which can be used to
%   set the termination tolerance, maximum number of iterations, and
%   output options (see MPOPTION for details). Uses default options
%   if this parameter is not given. Returns the final complex voltages,
%   a flag which indicates whether it converged or not, and the number
%   of iterations performed.
%
%   See also RUNPF.

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Alberto Borghetti, University of Bologna, Italy
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% default arguments
if nargin < 7
    mpopt = mpoption;
end

%% options
tol     = mpopt.pf.tol;
max_it  = mpopt.pf.gs.max_it;

%% initialize
converged = 0;
i = 0;
V = V0;
Vm = abs(V);

%% set up indexing for updating V
npv = length(pv);
npq = length(pq);

%% evaluate F(x0)
mis = V .* conj(Ybus * V) - Sbus;
F = [   real(mis([pv; pq]));
        imag(mis(pq))   ];

%% check tolerance
normF = norm(F, inf);
if mpopt.verbose > 1
    fprintf('\n it    max P & Q mismatch (p.u.)');
    fprintf('\n----  ---------------------------');
    fprintf('\n%3d        %10.3e', i, normF);
end
if normF < tol
    converged = 1;
    if mpopt.verbose > 1
        fprintf('\nConverged!\n');
    end
end

%% do Gauss-Seidel iterations
while (~converged && i < max_it)
    %% update iteration counter
    i = i + 1;

    %% update voltage
    %% at PQ buses
    for k = pq(1:npq)'
        V(k) =  V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k);
    end

    %% at PV buses
    if npv
        for k = pv(1:npv)'
            Sbus(k) = real(Sbus(k)) + 1j * imag( V(k) .* conj(Ybus(k,:) * V));
            V(k) =  V(k) + (conj(Sbus(k) / V(k)) - Ybus(k,:) * V ) / Ybus(k,k);
%           V(k) = Vm(k) * V(k) / abs(V(k));
        end
        V(pv) = Vm(pv) .* V(pv) ./ abs(V(pv));
    end

    %% evalute F(x)
    mis = V .* conj(Ybus * V) - Sbus;
    F = [   real(mis(pv));
            real(mis(pq));
            imag(mis(pq))   ];

    %% check for convergence
    normF = norm(F, inf);
    if mpopt.verbose > 1
        fprintf('\n%3d        %10.3e', i, normF);
    end
    if normF < tol
        converged = 1;
        if mpopt.verbose
            fprintf('\nGauss-Seidel power flow converged in %d iterations.\n', i);
        end
    end
end

if mpopt.verbose
    if ~converged
        fprintf('\nGauss-Seidel power flow did not converge in %d iterations.\n', i);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function pathS = MC_VGCIR(S0,r,d,T,C,G,M,kappa,eta,...
    lambda,NTime,NSim,NBatches)
intNt = 20;                 % steps in between orginial grid points
allsteps = intNt * NTime;   % grid containing all necessary steps
dT = T / NTime;             % time step
daT = T/allsteps;           % time step
time = 0 : T/NTime : T;     % time for martingale correction

pathS = zeros(NSim,NTime+1,NBatches);% output
lnS = ones(NSim,NTime+1);            % used for each batch
lnS(:,1) = log(S0*exp(-d*T));        % set S(0) dividend adjusted

% precompute constants
psiVG = (-1i)*C*log(G*M/(G*M+(M-G)-1));   % char exp
gamma = sqrt(kappa^2-2*lambda^2*1i*psiVG);% CIR par
denom = ( cosh(0.5*gamma*time) ...        % denom
    + kappa*sinh(0.5*gamma*time)./gamma ).^(2*kappa*eta*lambda^(-2));
% coth is inf at 0
phiCIR(time>0) = kappa^2*eta*time(time>0)*lambda^(-2) ...
         + 2*1i*psiVG./(kappa+gamma.*coth(gamma*time(time>0)/2)) ...
         - log(denom(time>0));            % char func
phiCIR(1) = log(denom(1));                % char func start
omegaT = -phiCIR;                         % martingale correction
omegaT(1) = 0;                            % maringale correction in 0

deg = 4*eta*kappa/lambda^2;
fac1 = 4*kappa*exp(-kappa*daT)/lambda^2/(1-exp(-kappa*daT));
fac2 = lambda^2*(1-exp(-kappa*daT))/4/kappa;

Y = zeros(NSim,NTime);                   % Integrated clock

for l = 1 : NBatches                     % batch loop
   % Generating time change
   yy = ones(NSim,allsteps+1);              % stochastic clock
   for n = 1 : allsteps
       Nvec = poissrnd(0.5 * fac1 * yy(:,n));  % Poissonians
       yy(:,n+1) = fac2 * chi2rnd(deg+2*Nvec); % stochastic time
   end

   for m=1:NTime
        Y(:,m+1) = Y(:,m) + daT * sum(yy(:,1+(m-1)*intNt:m*intNt),2);
   end

   Intensity = C * (Y(:,2:end)-Y(:,1:end-1));           % intensity
   DGam = gamrnd(Intensity,1/M) - gamrnd(Intensity,1/G);% diff gamma proc
   diffomegaT = omegaT(2:end) - omegaT(1:end-1);        % Martingale correction

   for m=2:NTime+1                      % time loop
       lnS(:,m) = lnS(:,m-1) + (r-d)*dT + diffomegaT(m-1) + DGam(:,m-1);
   end
   pathS(:,:,1) = exp(lnS);
end

author: Sable

repo: mcbench-benchmarks


text:
function [ top, btri, bedg, triangle_node, triangle_neighbor ] = swapec ( ...
  i, top, btri, bedg, node_num, node_xy, triangle_num, triangle_node, ...
  triangle_neighbor, work )

%*****************************************************************************80
%
%% SWAPEC swaps diagonal edges until all triangles are Delaunay.
%
%  Discussion:
%
%    The routine swaps diagonal edges in a 2D triangulation, based on
%    the empty circumcircle criterion, until all triangles are Delaunay,
%    given that I is the index of the new vertex added to the triangulation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    Original FORTRAN77 version by Barry Joe,
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, integer I, the index of the new vertex.
%
%    Input/output, integer TOP, the index of the top of the stack.
%    On output, TOP is zero.
%
%    Input/output, integer BTRI, BEDG; on input, if positive, are the
%    triangle and edge indices of a boundary edge whose updated indices
%    must be recorded.  On output, these may be updated because of swaps.
%
%    Input, integer NODE_NUM, the number of node.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of
%    the points.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input/output, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the triangle incidence list.
%    May be updated on output because of swaps.
%
%    Input/output, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle neighbor list;
%    negative values are used for links of the counter-clockwise linked
%    list of boundary edges;  May be updated on output because of swaps.
%
%      LINK = -(3*I + J-1) where I, J = triangle, edge index.
%
%    Workspace, integer WORK(MAXST); on input, entries 1 through TOP
%    contain the indices of initial triangles (involving vertex I)
%    put in stack; the edges opposite I should be in interior;  entries
%    TOP+1 through MAXST are used as a stack.
%

%
%  Determine whether triangles in stack are Delaunay, and swap
%  diagonal edge of convex quadrilateral if not.
%
  x = node_xy(1,i);
  y = node_xy(2,i);

  while ( 1 )

    if ( top <= 0 )
      break
    end

    t = work(top);
    top = top - 1;

    if ( triangle_node(1,t) == i )
      e = 2;
      b = triangle_node(3,t);
    elseif ( triangle_node(2,t) == i )
      e = 3;
      b = triangle_node(1,t);
    else
      e = 1;
      b = triangle_node(2,t);
    end

    a = triangle_node(e,t);
    u = triangle_neighbor(e,t);

    if ( triangle_neighbor(1,u) == t )
      f = 1;
      c = triangle_node(3,u);
    elseif ( triangle_neighbor(2,u) == t )
      f = 2;
      c = triangle_node(1,u);
    else
      f = 3;
      c = triangle_node(2,u);
    end

    swap = diaedg ( x, y, node_xy(1,a), node_xy(2,a), node_xy(1,c), node_xy(2,c), ...
      node_xy(1,b), node_xy(2,b) );

    if ( swap == 1 )

      em1 = i4_wrap ( e - 1, 1, 3 );
      ep1 = i4_wrap ( e + 1, 1, 3 );
      fm1 = i4_wrap ( f - 1, 1, 3 );
      fp1 = i4_wrap ( f + 1, 1, 3 );

      triangle_node(ep1,t) = c;
      triangle_node(fp1,u) = i;
      r = triangle_neighbor(ep1,t);
      s = triangle_neighbor(fp1,u);
      triangle_neighbor(ep1,t) = u;
      triangle_neighbor(fp1,u) = t;
      triangle_neighbor(e,t) = s;
      triangle_neighbor(f,u) = r;

      if ( 0 < triangle_neighbor(fm1,u) )
        top = top + 1;
        work(top) = u;
      end

      if ( 0 < s )

        if ( triangle_neighbor(1,s) == u )
          triangle_neighbor(1,s) = t;
        elseif ( triangle_neighbor(2,s) == u )
          triangle_neighbor(2,s) = t;
        else
          triangle_neighbor(3,s) = t;
        end

        top = top + 1;

        if ( node_num < top )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'SWAPEC - Fatal error!\n' );
          fprintf ( 1, '  Exceeded stacksize.\n' );
          error ( 'SWAPEC - Fatal error!' );
        end

        work(top) = t;

      else

        if ( u == btri & fp1 == bedg )
          btri = t;
          bedg = e;
        end

        l = - ( 3 * t + e - 1 );
        tt = t;
        ee = em1;

        while ( 0 < triangle_neighbor(ee,tt) )

          tt = triangle_neighbor(ee,tt);

          if ( triangle_node(1,tt) == a )
            ee = 3;
          elseif ( triangle_node(2,tt) == a )
            ee = 1;
          else
            ee = 2;
          end

        end

        triangle_neighbor(ee,tt) = l;

      end

      if ( 0 < r )

        if ( triangle_neighbor(1,r) == t )
          triangle_neighbor(1,r) = u;
        elseif ( triangle_neighbor(2,r) == t )
          triangle_neighbor(2,r) = u;
        else
          triangle_neighbor(3,r) = u;
        end

      else

        if ( t == btri & ep1 == bedg )
          btri = u;
          bedg = f;
        end

        l = - ( 3 * u + f - 1 );
        tt = u;
        ee = fm1;

        while ( 0 < triangle_neighbor(ee,tt) )

          tt = triangle_neighbor(ee,tt);

          if ( triangle_node(1,tt) == b )
            ee = 3;
          elseif ( triangle_node(2,tt) == b )
            ee = 1;
          else
            ee = 2;
          end

        end

        triangle_neighbor(ee,tt) = l;

      end

    end

  end

  return
end


text:
function c8_normal_01_test ( )

%*****************************************************************************80
%
%% C8_NORMAL_01_TEST tests C8_NORMAL_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_NORMAL_01_TEST\n' );
  fprintf ( 1, '  C8_NORMAL_01 computes pseudorandom complex values\n' );
  fprintf ( 1, '  normally distributed in the unit circle.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ r, seed ] = c8_normal_01 ( seed );
    fprintf ( 1, '  %6d  (%14f,  %14f)\n', i, real ( r ), imag ( r ) );
  end

  return
end


text:
function R=quat2RotMat(q,handed)
%%QUAT2ROTMAT Turn a unit quaternion into an equivalent rotation matrix.
%             The multiplication rules for the quaternion algebra can be
%             chosen to support standard right-handed quaternion rotation,
%             or non-standard left-handed rotations that some authors use.
%
%INPUTS: q A 4X1 unit quaternion corresponding to the rotation matrix. The
%          quaternion is ordered [cos(theta/2);sin(theta/2)u'] where u is a
%          unit vector for the axis of rotation and theta is the rotation
%          angle about that unit vector according to the specified
%          handedness. The ordering of the elements corresponds to the
%          hypercomplex decomposition q(1)+i*q(2)+j*q(3)+k*q(4), where i,
%          j, and k are all roots of -1.
%   handed The handedness of the quaternion. If omitted, it is assumed that
%          the quaternion is right-handed (the standard). Possible values
%          are:
%          'right' The default if omitted. The quaternion multiplication is
%                  assumed right-handed (standard).
%          'left'  The quaternion multiplication is assumed left-handed.
%                  This is used in someplaces, including the reference from
%                  Shuster, below.
%
%OUTPUTS: R A 3X3 orthonormal rotation matrix.
%
%If q does not have unit magnitude, then R will not be a rotation matrix.
%
%The formula for turning a quaternion into a left-handed rotation matrix is
%given in [1]. However, right-handed quaternion algebra, as used in [2] is
%far more common.
%
%A quaternion form q(1)+i*q(2)+j*q(3)+k*q(4) that obeys right-handed
%multiplication rules supports the following rules for multiplication of i,
%j, and k, where an item in a row is multiplied by an item in the column to
%get the result:
%  i,  j, k
%i -1, k,-j
%j -k,-1, i
%k  j,-i,-1
%On the other hand, left-handed multiplication rules flip the signs of the
%off-diagonal terms:
%  i,  j, k
%i -1,-k, j
%j  k,-1,-i
%k -j, i,-1
%
%REFERENCES:
%[1] M. D. Shuster, "A survey of attitude representations," The Journal of
%    Astronautical Sciences, vol. 41, no. 4, pp. 439-517, Oct. -Dec. 1993.
%[2] Weisstein, Eric W. "Quaternion." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/Quaternion.html
%
%August 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(handed))
    handed='right';
end

switch(handed)
    case 'left'
        %The rotation matrix for a left-handed quaternion algebra, as in
        %[1].
        R=[q(1)^2+q(2)^2-q(3)^2-q(4)^2, 2*(q(2)*q(3)+q(1)*q(4)),    2*(q(2)*q(4)-q(1)*q(3));
           2*(q(2)*q(3)-q(1)*q(4)),     q(1)^2-q(2)^2+q(3)^2-q(4)^2,2*(q(3)*q(4)+q(1)*q(2));
           2*(q(2)*q(4)+q(1)*q(3)),     2*(q(3)*q(4)-q(1)*q(2)),    q(1)^2-q(2)^2-q(3)^2+q(4)^2];
    case 'right'
        %The rotation matrix for a right-handed quaternion algebra is the
        %transpose of the matrix from [1].
        R=[q(1)^2+q(2)^2-q(3)^2-q(4)^2, 2*(q(2)*q(3)-q(1)*q(4)),    2*(q(2)*q(4)+q(1)*q(3));
           2*(q(2)*q(3)+q(1)*q(4)),     q(1)^2-q(2)^2+q(3)^2-q(4)^2,2*(q(3)*q(4)-q(1)*q(2));
           2*(q(2)*q(4)-q(1)*q(3)),     2*(q(3)*q(4)+q(1)*q(2)),    q(1)^2-q(2)^2-q(3)^2+q(4)^2];
    otherwise
        error('Invalid handedness provided.')
end

end

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


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

subsetU = U(:, 1:K);
Z = X * subsetU;

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

end

repo: machine-learning


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function [x_bin_new, x_new_idx] = binary_resampling(binary_models, weights)
% BINARY_RESAMPLING: Function re-samples binary particles based on the
% current weights using systematic resampling.

% Determine the size of x_bin
[n_samples,n_var] = size(binary_models);

% Declare a vector to store new samples and indices
x_bin_new  = zeros(n_samples, n_var);
x_new_idx  = zeros(n_samples, 1);

% Initialize system (u = U[0,1] & v = weights*n_samples & c = v(1))
u = rand;
v = weights*n_samples;
c = v(1);

% Declare iteration counter
counter = 1;

% Resample each particle
for i=1:n_samples

    % Update counter and c value
    while c < u
        counter = counter+1;
        c       = c + v(counter);
    end

    % Assign new binary vector and probability
    x_bin_new(i,:) = binary_models(counter,:);
    x_new_idx(i)   = counter;

    % Update random value
    u = u+1;

end

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function flat = vol2flatAllROIsLevels(volume,flat)
%
% flat = vol2flatAllROIsLevels(volume,flat)
%
% Calls vol2flatROI with the currently selected ROI.  Called from
% callback in xformFlatMenu.
%
% rmk 1/15/99
%
% Modifications:
% djh, 2/2001, replaced globals with local variables
% ras, 10/2004, version for flat level view

if ~volume.selectedROI
  myErrorDlg('Must have a selected ROI in the Volume window before it can be transformed to the Inplane window.');
end

for r=1:length(volume.ROIs)
  volROI = volume.ROIs(r);
  flatROI = vol2flatROILevels(volROI,volume,flat);
  flat = addROI(flat,flatROI,1);
end

author: vistalab

repo: vistasoft


text:
function Output = Bayesian_RPCAmcmc_MarkovDep(X0,Theta0,Row,Col,hyperpara,MCMCpara)
%Bayesian_RPCAmcmc_withMarkovDep: Bayesian robust principle component analysis
%implemented by MCMC, considering Markov dependency on the sparse term in time and space.
%
% Model: X0 = D*(Delta.*Z)*S + S2.*Z2 + E
% ----------------------------------------------------
%
%USAGE: Output = Bayesian_RPCAmcmc_withMarkovDep(X0,Theta0,Row,Col,hyperpara,MCMCpara)
%
%INPUT :
%   X0: P x N, input data matrix. For the video application, every column is a frame of the video.
%
%   Theta0: struct data including the initial parameters of the model. The parameters of the model are initialized as follows,
%         D = Theta0.D: P x K matrix as a dictionary for lowrank learning
%         S = Theta0.S: K x N coefficient matrix
%         Z = Theta0.Z: K x K diagnal binary matrix for rank learning
%         Delta = Theta0.Delta: K x K diagnal matrix for weighting Z.
%         Tao = Theta0.Tao: precision of Delta. Here we set Tao=1.
%         Pi = Theta0.Pi: K x 1 vector. Pi is the probability of Z=1.
%
%         gamma_epsi = Theta0.gamma_epsi:  scalar, precision of the noise. Here we assumed the noise is stationary.
%
%         S2 = Theta0.S2: P x N matrix for Sparse Component learning
%         Z2 = Theta0.Z2: P x N binary matrix for Sparse Component learning
%         gamma_s = Theta0.gamma_s: precision of the S2
%         Pi2 = Theta0.Pi2: P x N . Pi2 is the probability of Z2=1.
%         Note that Pi2 here is a matrix, different from the model structure in function Bayesian_RPCAmcmc.
%
%   Row and Col: Row and column of the video frame.
%
%   hyperpara: MCMC hyperparameters
%     hyperpara.a0: scalar, hyperparameter 1 for Pi [1/K or 1/150 if K<150]
%       hyperpara.b0: scalar, hyperparameter 2 for Pi [1-hyperpara.a0]
%       hyperpara.c0: scalar, hyperparameter 1 for Tao precision [1e-6]
%       hyperpara.d0: scalar, hyperparameter 2 for Tao precision [1e-6]
%       hyperpara.e0: scalar, hyperparameter 1 for nosie precision [1e-6]
%       hyperpara.f0: scalar, hyperparameter 2 for nosie precision [1e-6]
%       hyperpara.g0: scalar, hyperparameter 1 for S2 precision [1e-6]
%       hyperpara.h0: scalar, hyperparameter 2 for S2 precision [1e-6]
%       hyperpara.alpha0: scalar,hyperparameter 1 for Pi2 [0.01*N]
%       hyperpara.beta0: scalar,hyperparameter 2 for Pi2 [0.99*N]
%       hyperpara.alpha1: scalar,hyperparameter 3 for Pi2 [0.99*N]
%       hyperpara.beta1: scalar,hyperparameter 4 for Pi2 [0.01*N]
%
%
%   MCMCpara: MCMC parameters
%       MCMCpara.nBurnin: scalar, number of burnin iterations [2000]
%       MCMCpara.nCollect: scalar, number of collected samples [1000]
%
%OUTPUT:
%   Output: struct data.
%         Output.Lowrank_mean: P x N matrix, mean of the Lowrank Component.
%         Output.Lowrank_std: P x N matrix, std of the Lowrank Component.
%         Output.Sparse_mean: P x N matrix, mean of the Sparse Component.
%         Output.Sparse_std: P x N matrix, std of the Sparse Component.
%         Output.Gamma_epsi_mean: scalar, mean of the noise precision.
%         Output.Gamma_epsi_std: scalar, std of the noise precision.
%         Output.rankL_mean: scalar, mean of the estimated rank of the lowrank Component.
%         Output.rankL_std: scalar, std of the estimated rank of the lowrank Component.
%         Output.NumSparse_mean: scalar, mean of the estimated number of the Sparse Component.
%         Output.NumSparse_std: scalar, std of the estimated number of the Sparse Component.

%--------------------------------------------------------------------------
% References:
% X. Ding, L. He and L. Carin, Bayesian Robust Principal Component Analysis, submitted to IEEE Trans. Image Processing (2010)
%
% Xinghao Ding, Lihan He, ECE, Duke University
% Created: Apr. 12, 2010,
% Last change: Aug. 2, 2010.
%--------------------------------------------------------------------------

% ---------------------
% check input arguments
% ---------------------

% P -- Dimension of the data vector
% N -- Number of samples
[P,N] = size(X0);
% K -- the largest possible rank
K = size(Theta0.D,2);

if nargin<6
    MCMCpara.nBurnin=100;
    MCMCpara.nCollect=100;
end
if nargin<5
    % Hyperparameters
    if K<150
        hyperpara.a0 = 1/150;
    else
        hyperpara.a0 = 1/K;
    end
    hyperpara.b0 = 1-hyperpara.a0;
    hyperpara.c0 = 1e-6;
    hyperpara.d0 = 1e-6;
    hyperpara.e0 = 1e-6;
    hyperpara.f0 = 1e-6;
    hyperpara.g0 = 1e-6;
    hyperpara.h0 = 1e-6;
    hyperpara.alpha0 = 0.01*N;
    hyperpara.beta0 = 0.99*N;
    hyperpara.alpha1 = 0.99*N;
    hyperpara.beta1 = 0.01*N;

end

if isempty(hyperpara)
    % P -- Dimension of the data vector
    % N -- Number of samples
    [P,N] = size(X0);
    % K -- the largest possible rank
    K = size(Theta0.D,2);

    % Hyperparameters
    if K<150
        hyperpara.a0 = 1/150;
    else
        hyperpara.a0 = 1/K;
    end
    hyperpara.b0 = 1-hyperpara.a0;
    hyperpara.c0 = 1e-6;
    hyperpara.d0 = 1e-6;
    hyperpara.e0 = 1e-6;
    hyperpara.f0 = 1e-6;
    hyperpara.g0 = 1e-6;
    hyperpara.h0 = 1e-6;
    hyperpara.alpha0 = 0.01*N;
    hyperpara.beta0 = 0.99*N;
    hyperpara.alpha1 = 0.99*N;
    hyperpara.beta1 = 0.01*N;
end

if isempty(MCMCpara)
   MCMCpara.nBurnin=100;
    MCMCpara.nCollect=100;
end

a0 = hyperpara.a0;
b0 = hyperpara.b0;
c0 = hyperpara.c0;
d0 = hyperpara.d0;
e0 = hyperpara.e0;
f0 = hyperpara.f0;
g0 = hyperpara.g0;
h0 = hyperpara.h0;
alpha0 = hyperpara.alpha0;
beta0 =  hyperpara.beta0;
alpha1 = hyperpara.alpha1;
beta1 =  hyperpara.beta1;

D = Theta0.D;
S = Theta0.S;
Z = Theta0.Z;
Delta = Theta0.Delta;
Tao = Theta0.Tao;
Pi = Theta0.Pi;
gamma_epsi = Theta0.gamma_epsi;
S2 = Theta0.S2;
Z2 = Theta0.Z2;
gamma_s = Theta0.gamma_s;
Pi2 = Theta0.Pi2;

Lowrank_Comp = zeros(P,N);
Sparse_Comp = zeros(P,N);

for iter=1:(MCMCpara.nBurnin+MCMCpara.nCollect)

    X = X0-D*diag(Delta.*Z)*S-S2.*Z2;

    %---------------------------------------------------------------
    % Low-rank component sampling
    %---------------------------------------------------------------

    for k=1:K
        Tao(k) =1;
        X = X+ Delta(k)*Z(k)*D(:,k)*S(k,:);

        %Sample D----------------------------------------------------------
        sigma_Dk = 1/(gamma_epsi*(Delta(k)).^2*Z(k)*(S(k,:)*S(k,:)')+P);
        mu_Dk = gamma_epsi*sigma_Dk* (Delta(k))*Z(k)*X*S(k,:)';
        D(:,k) = mu_Dk + randn(P,1)*sqrt(sigma_Dk);
        %clear sigma_Dk mu_D
        %------------------------------------------------------------------

        %Sample S----------------------------------------------------------

        Dk = D(:,k);
        sigS1 = 1/(1 + gamma_epsi*Z(k)*(Dk'*Dk)*(Delta(k)).^2);
        muS1 = gamma_epsi*sigS1*Z(k)*Delta(k)*D(:,k)'*X;
        S(k,:) = randn(1,N)*sqrt(sigS1) + muS1;
        %------------------------------------------------------------------

        %Sample Delta(k)
        sig_Delta = 1/(Tao(k)+gamma_epsi*(Dk'*Dk)*(S(k,:)*S(k,:)'));
        mu_Delta = gamma_epsi*sig_Delta*D(:,k)'*X*S(k,:)';
        if Z(k) ==1
            %Delta(k) = randn(1)*sqrt(sig_Delta)+mu_Delta;
            Delta(k) = normrnd(mu_Delta,sqrt(sig_Delta));
        else
            %Delta(k) = randn(1);
            Delta(k) =  normrnd(0,sqrt(1/Tao(k)));
        end
        %------------------------------------------------------------------

%        %Sample Z---------------------------------------------------------
        Sk = S(k,:);
        Dk = D(:,k);
        temp =  - 0.5*gamma_epsi*(Delta(k)).^2*(Dk'*Dk)*(Sk*Sk') + gamma_epsi*Delta(k)*Dk'*(X*S(k,:)');
        p1 = exp(temp)*Pi(k);
        p0 = 1-Pi(k);
        Z(k) = rand > p0/(p0+p1);
        %------------------------------------------------------------------

%       %Clappsed Gibbs Sample Z
%         tmpz = log(Pi(k)+eps) - log(1-Pi(k)+eps) + 0.5*log(sig_Delta) + (mu_Delta^2)/(2*sig_Delta) + 0.5*log(Tao(k));
%         if rand < 1/(1+exp(-tmpz))
%             Z(k) = 1; Delta(k) = normrnd(mu_Delta,sqrt(sig_Delta));
%         else
%             Z(k) = 0; Delta(k) =  normrnd(0,sqrt(1/Tao(k)));
%         end

        %sample Pi--------------------------------------------------
        ai = a0 + Z(k);
        bi = b0 + 1 - Z(k);
        Pi(k) =betarnd(ai,bi);
        %------------------------------------------------------------------

        %sample Tao
%         if Z(k) ==1
%             c1 =c0+1/2;
%             d1 =d0 + 0.5*(Delta(k)).^2;
%             Tao(k) = gamrnd(c1,1./d1);
%         else
%             Tao(k) = 1;
%         end
        %------------------------------------------------------------------

        X = X - Delta(k)*Z(k)*D(:,k)*S(k,:);
    end

    %---------------------------------------------------------------
    % Sparse component sampling
    %---------------------------------------------------------------

    X = X + S2.*Z2;

    % Sample S2
    sig_S2 = 1./(gamma_s + gamma_epsi*Z2);
    mu_S2 = gamma_epsi*sig_S2.*Z2.*X;
    S2 = randn(P,N).*sqrt(sig_S2)+mu_S2;

    % Sample Z2
    temp =  exp(-0.5*gamma_epsi*(S2.^2-2*S2.*X));
%     p1 = repmat(Pi2,[1,N]).*temp;
%     p0 = repmat(1- Pi2,[1,N]);
    p1 = Pi2.*temp;
    p0 = 1- Pi2;
    Z2 = rand(P,N)>p0./(p1+p0);

    X = X - S2.*Z2;

    %sample gamma_s
    g1 = g0 + 0.5*P*N;
    tempS2 = S2.*S2;
    h1 = h0 + 0.5*sum(tempS2(:));
    gamma_s = gamrnd(g1,1/h1);

    %Sample Pi2
%     sumZ2 = sum(Z2,2);
%     a2 = a1 + sumZ2;
%     b2 = b1 + N -sumZ2;
%     Pi2 = betarnd(a2,b2);

    %h = fspecial('average');
    h = [0.1,0.1,0.1;0.1,0.2,0.1;0.1,0.1,0.1];
    for n = 1:N
        temp1 = reshape(Z2(:,n),[Row,Col]);
        temp2 = imfilter(double(temp1),h)>0.6;
        tempZ2(:,n) = reshape(temp2,[Row*Col,1]);
    end

    a2 = alpha0 +Z2(:,1);
    b2 = beta0 + 1 - Z2(:,1);
    Pi2(:,1) = betarnd(a2,b2);

    for t = 2:N
%         a2 = alpha0*(~Z2(:,t-1))+alpha1*(Z2(:,t-1)) +Z2(:,t);
%         b2 = beta0*(~Z2(:,t-1))+beta1*(Z2(:,t-1)) +1 - Z2(:,t);
%         a2 = alpha0*(~(Z2(:,t-1).*tempZ2(:,t-1)))+alpha1*(Z2(:,t-1).*tempZ2(:,t-1)) +Z2(:,t);
%         b2 = beta0*(~(Z2(:,t-1).*tempZ2(:,t-1)))+beta1*(Z2(:,t-1).*tempZ2(:,t-1)) +1 - Z2(:,t);
%         a2 = alpha0*(~(tempZ2(:,t-1)))+alpha1*(tempZ2(:,t-1)) +Z2(:,t);
%         b2 = beta0*(~(tempZ2(:,t-1)))+beta1*(tempZ2(:,t-1)) +1 - Z2(:,t);
          a2 = alpha0*(~(tempZ2(:,t)))+alpha1*(tempZ2(:,t)) +Z2(:,t);
          b2 = beta0*(~(tempZ2(:,t)))+beta1*(tempZ2(:,t)) +1 - Z2(:,t);

%         a2 = alpha0*(~(tempZ2(:,t-1).*tempZ2(:,t)))+alpha1*(tempZ2(:,t-1).*tempZ2(:,t)) +Z2(:,t);
%         b2 = beta0*(~(tempZ2(:,t-1).*tempZ2(:,t)))+beta1*(tempZ2(:,t-1).*tempZ2(:,t)) +1 - Z2(:,t);

%         a2 = alpha0*(~(Z2(:,t-1).*tempZ2(:,t)))+alpha1*(Z2(:,t-1).*tempZ2(:,t)) +Z2(:,t);
%         b2 = beta0*(~(Z2(:,t-1).*tempZ2(:,t)))+beta1*(Z2(:,t-1).*tempZ2(:,t)) +1 - Z2(:,t);

        Pi2(:,t) = betarnd(a2,b2);
    end

    %---------------------------------------------------------------
    % noise component sampling
    %---------------------------------------------------------------

    %Sample gamma_epsi
    e1 = e0 + 0.5*P*N;
    f1 = f0 + 0.5*sum(sum(X.^2));
    gamma_epsi = gamrnd(e1,1./f1);

    %------------------------

    % Collect samples
    if iter>MCMCpara.nBurnin
        ii = ceil(iter-MCMCpara.nBurnin);

%         Lowrank_Comp(:,:,ii) = D*diag(Delta.*Z)*S;
%         Sparse_Comp(:,:,ii) = S2.*Z2;
        Lowrank_Comp = Lowrank_Comp + D*diag(Delta.*Z)*S;
        Sparse_Comp = Sparse_Comp + S2.*Z2;

        Gamma_epsi(ii) = gamma_epsi;
        tmpRank(ii) = length(find(Z~=0));
        %tmpNumSparse(ii) = length(find(Sparse_Comp(:,:,ii)~=0));
        %mse_rec(iter) = sum(sum((X0-Lowrank_Comp{ii}-Sparse_Comp{ii}).^2))/(P*N);
    end

end

% Output.Lowrank_mean = mean(Lowrank_Comp,3);
% Output.Lowrank_std = std(Lowrank_Comp,0,3);
% Output.Sparse_mean = mean(Sparse_Comp,3);
% Output.Sparse_std = std(Sparse_Comp,0,3);
Output.Lowrank_mean = Lowrank_Comp/MCMCpara.nCollect;
Output.Sparse_mean = Sparse_Comp/MCMCpara.nCollect;
Output.Gamma_epsi_mean = mean(Gamma_epsi);
Output.Gamma_epsi_std = std(Gamma_epsi);
Output.rankL_mean = mean(tmpRank);
Output.rankL_std = std(tmpRank);
% Output.NumSparse_mean = mean(tmpNumSparse);
% Output.NumSparse_std = std(tmpNumSparse);
end