function?[weightsusedsigma2errbarscountalphaindexbasis]?=?BCS_fast_rvm（PHItsigma2etaadaptiveoptimalscale）
%------------------------------------------------------------------
%?The?BCS?algorithm?for?the?following?paper:
%?“Bayesian?Compressive?Sesning“?（Preprint?2007）.?The?algorithm?
%?adopts?from?the?fast?RVM?algorithm?[Tipping?&?Faul?2003].
%?Coded?by:?Shihao?Ji?ECE?Duke?University
%?last?change:?Jan.?2?2007
%?You?are?suggested?to?use?mt_CS.m?for?improved?robustness
%------------------------------------------------------------------
%?Input?for?BCS:
%???PHI:?projection?matrix
%???t:???CS?measurements
%???sigma2:?initial?noise?variance
%??????If?measurement?noise?exists?and/or?w?is?not?truely?sparse?
%?????????????then?sigma2?=?std（t）^2/1e2?（suggested）
%??????If?no?measurement?noise?and?w?is?truely?sparse
%?????????????then?sigma2?=?std（t）^2/1e6?（suggested）
%??????This?term?is?in?fact?not?updated?in?the?implementation?to?allow?
%??????the?fast?algorithm.?For?this?reason?you?are?recommended?to?use
%??????mt_CS.m?in?which?the?noise?variance?is?marginalized.
%???eta:?threshold?for?stopping?the?algorithm?（suggested?value:?1e-8）
%?Input?for?Adaptive?CS:
%???adaptive:?generate?basis?for?adpative?CS??（default:?0）
%???optimal:?use?the?rigorous?implementation?of?adaptive?CS??（default:?1）
%???scale:?diagonal?loading?parameter?（default:?0.1）
%?Output:
%???weights:??sparse?weights
%???used:?????the?positions?of?sparse?weights
%???sigma2:???re-estimated?noise?variance
%???errbars:??one?standard?deviation?around?the?sparse?weights
%???basis:????if?adaptive==1?then?basis?=?the?next?projection?vector
%
if?nargin?????adaptive?=?0;
end
if?nargin?????optimal?=?1;
end
if?nargin?????scale?=?0.1;
end

%?find?initial?alpha
[NM]?=?size（PHI）;
PHIt?=?PHI‘*t;
PHI2?=?sum（PHI.^2）‘;
ratio?=?（PHIt.^2）./PHI2;
[maxrindex]?=?max（ratio）;
alpha?=?PHI2（index）/（maxr-sigma2）;
%?compute?initial?mu?Sig?S?Q
phi?=?PHI（:index）;
Hessian?=?alpha?+?phi‘*phi/sigma2;
Sig?=?1/Hessian;
mu?=?Sig*PHIt（index）/sigma2;
left?=?PHI‘*phi/sigma2;
S?=?PHI2/sigma2-Sig*left.^2;
Q?=?PHIt/sigma2-Sig*PHIt（index）/sigma2*left;
%
for?count?=?1:10000
????s?=?S;?q?=?Q;
????s（index）?=?alpha.*S（index）./（alpha-S（index））;
????q（index）?=?alpha.*Q（index）./（alpha-S（index））;
????theta?=?q.^2-s;

????%?choice?the?next?alpha?that?maximizes?marginal?likelihood
????ml?=?-inf*ones（1M）;
????ig0?=?find（theta>0）;
????%?index?for?re-estimate
????[irefoowhich]?=?intersect（ig0index）;
????if?~isempty（ire）
????????Alpha?=?s（ire）.^2./theta（ire）;
????????delta?=?（alpha（which）-Alpha）./（Alpha.*alpha（which））;
????????ml（ire）?=?Q（ire）.^2.*delta./（S（ire）.*delta+1）-log（1+S（