
%
% BUGS example
%
close all
clear all

curdir=pwd;
n=256; level=4;
wf=[  -0.075765714789341,  -0.029635527645954, ...
  0.497618667632458,   0.803738751805216, ...
  0.297857795605542,  -0.099219543576935, ...
 -0.012603967262261,   0.032223100604071];
ip=10;
is=1;
%precisions tau1, tau2. Tau1 has small precision, large variance
alphatau1 = 0.5; betatau1=0.05; %mean= 0.05/0.5 var=0.05/0.25=0.2
alphatau2=0.05; betatau2 = 0.5; %mean=0.5/0.05  var=0.5/0.0025=200
s=MakeSignal('HeaviSine', n);
subplot(3,1,1);plot(s) 
axis tight
sn=s + 0.7 * randn(size(s));
subplot(3,1,2);plot(sn) 
axis tight


%DD=FWT_PO(sn, level, wf);
DD = dwtr(sn, log2(n)-level, wf);

    finest = log2(n)-1;
    is = std( DD(dyad(finest)));
 
    dataStructu = struct( ...
    'DD', DD , ...
    'inputp', ip,...
    'inputs', is,...
    'at1', alphatau1,...
    'bt1', betatau1,...
     'L',15,...
    'N', n);

init0 = struct( 'tau', 1, 'tau1', 0.1,'theta0',0, 'tau0',10 );
% ===============bugs.txt========================
% model{
%   for (j in 1:N){
%     DD[j] ~ dnorm(theta[j], tau);
%     theta[j] <- p[j] * xi1[j] + (1-p[j]) * xi2[j];
%     xi2[j] ~ ddexp(0, tau1);
% 
% 
%    latent[j]~dcat(pi[1:L])  # L is the truncation point. We need to put a value for L
%    xi1[j] <-theta[latent[j]]
%                 
% 
%      
% 
%     p[j] ~ dbern(r[j]);
%     r[j] ~ dbeta(1, inputp)
%    }
% 
% # Sethuraman's stick-breaking  construction
%   for ( k in 1:L) {
%     r[k]~dbeta(1,alpha)
%     theta[k]~dnorm(theta0,tau0)    #baseline distribution.
%                                              # One can put other distribution also, 
%                                              # like double exponential
%
%    }
% 
% 
%  pi[1]<-r[1]
%    for ( j in 2:(L-1)) {
%         log(pi[j])<-log(r[j])+sum(D1[j,1:j-1])
%          for ( l in 1:j-1) { 
%                D1[j,l]<-log(1-r[l])   
%          }
%    }
%  pi[L]<-1-sum(pi[1:(L-1)])
% 
% 
%   theta0~dnorm(0,0.1)
%   tau0~dgamma(1,1)
%   alpha~dunif(0.4, 5) # based on Ohlssen and 
%   tau ~ dgamma(1, inputs);
%   tau1 ~ dgamma(at1, bt1);
% }

%==============================================
cd('C:\Brani\MyBugs\matbugs\')
[samples, stats, structArray] = matbugs(dataStructu, ...
		fullfile(pwd, 'wshrink1.txt'), ...
		'init', init0, ...
        'openBugs', 1, ...
		'nChains', 1, ...
		'view', 0, 'nburnin', 1000, 'nsamples', 5000, ...
        'monitorParams', {'theta'}, ...
		'Bugdir', 'C:/Program Files/BUGS');
 cd(curdir)
%==============================================

%ss=IWT_PO(stats.mean.theta, level, wf);
ss = idwtr(stats.mean.theta, log2(n)-level, wf);



subplot(3,1,3);plot(ss) 
% hold on
% plot(s,'r-')
axis tight

%stats.DIC.total



