function OUT=GibbsSample(objfun,IN)

%GIBBSSAMPLE		Gibbs Sampling inversion routine
%                   employing the Heat Bath algorithm.
%
%USAGE              OUT=GibbsSample(objfun,IN)
%
%INPUT              objfun: a string containting the name of the objective
%                       function
%                   IN: the Input structure (see below for description of the
%                       required fields for 'IN')
%OUTPUT             OUT: the Output structure (see below for description)
%
%EJEMPLO           1) short-run mode: SR_tst=GibbsSample('okada_insar_obj',QT_CORR_4763_9773_DSC)
%                   2) heat-bath mode: GS_tst=GibbsSample('okada_insar_obj',QT_CORR_4763_9773_DSC)
%
%DESCRIPTION
%   Gibbs Sampling w/Heat Bath is based on work described in 3 principal papers with many referenced therein:
%
%   1) 'Rapid Determination of the Critical Temperature in Simulated Annealing Inversion', Basu, A., Frazer, N.; Science, 1990
%   2) 'Freeze Bath Inversion for Estimation of Geoacoustic Parameters'; Chapman,N.R, Jaschke, L.; Problems in Underwater Acoustics; 2001
%   3) 'Importance reweighting reduces dependence on temperature in Gibbs samplers: an application to the coseismic geodetic inverse
%       problem'; Brooks, B.A., Frazer, L.N.; Geophysical Journal International; 2005
%
%   The 'IN' structure is divided into 2 fields: '.GS', which contains fields specific to the Gibbs Sampling algorithm and
%   '.OBJ', which contains fields specific to the objective function.
%
%   The '.GS' field  must contain the following fields:
%       .num_subrange:        number of subranges to divide range of each parameter (default = 25)
%       .ts_alg:              a string designating the sampling temperature algorithm. Highly
%                               recommended to use 'ShortRun' as I have not used 'FreezeProb' for
%                               a year or so and so it will probably not work with these later
%                               versions of GibbsSample. 'ShortRun' is what is discussed in the
%                               Brooks & Frazer (2005) paper and it is a robust method to determine
%                               a range of sampling temperatures as
%                               discussed in Brooks&Frazer (2005).
%       .ts_num:              the number of temperature steps during the
%                               'burn-in' towards sampling temperature
%       .i_save:              the number of iterations between save
%                               commands
%       .num_bins:            the number of bins used in defining marginal
%                               distribution plots (default = 30)
%       .ts (optional):       the sampling temperature. If .ts is defined
%                               then the routine will run in Heat-Bath mode. If not, then the
%                               routine will go into sampling temperature mode according to
%                               .ts_alg.
%
%       ShortRun.init_T:      initial temp.
%       ShortRun.final_T:     final temp.
%       ShortRun.num_T:       number of temp. steps
%       ShortRun.num_R:       number of random starting models (default = 5)
%       ShortRun.num_sweep:   number of short-run sweeps  (default =25)


%
%   The '.OBJ' field comprises all of the fields which must be passed to the objective function.
%   The user is responsible for the correct creation of these fields, though, '.OBJ' must contain two fields:
%       .MB:  an nx2 array containing the model parameter bounds in the format [min max].
%       .xlabels: an nx1 cell string array of labels for plotting such as -
%                   {'dx','dy','l','w','d','strike','dip','u1','u2'} for
%                   the 9 parameter okada dislocation problem
%
% see also FuzzySample, Sweep, HeatBath, ShortRun, FreezeProb, RobustTc, RapidTc
%
% BAB   4/6/05

%-------------
%CHECK & ASSIGN FIELDS
%-------------
check_GS_fields(IN)
GS=IN.GS;
OBJ=IN.OBJ;

%-------------------------------------------------
%SET RANDOM NUMBER GENERATOR
%-------------------------------------------------
rand('state',13); %set to the 13th state so that the results are exactly reproducible%rand('state',sum(100*clock))

%----------------------------
%CHECK MODEL PARAMETER BOUNDS
%-----------------------------
if (find(OBJ.MB(:,2)-OBJ.MB(:,1) < 0))
    error('Model parameter bounds must be in format [min max]')
end

%---------------------------------------------------------
%EITHER DO A FULL HB RUN OR DETERMINE SAMPLING TEMPERATURE
%---------------------------------------------------------

if isfield(GS,'ts') %if there is a predefined sampling temperature then go into Heat-Bath Mode
    
    %******************
    %*** HEAT-BATH ***
    %******************
    fprintf('\n\nMODE: HeatBath \n');

    %check to make sure the saving increment is greater than the
    %burn-in time
    if GS.i_save <= GS.ts_num
        error('The saving increment (i_save) must be greater than the burn-in time (tc_num)')
    end

    if GS.ts < 1
        error('T samp < 1, check defintion of Ts')
    end
    
    if isfield(GS,'num_bins')==1
    else %DEFAULT VALUES
        GS.num_bins=30;
    end

    %assign 'burn-in' cooling schedule
    c=GS.ts*log10(1+GS.ts_num);     %'tc_num' is the number of cooling steps for the burn-in time
    t=1:GS.ts_num;
    T=1000*exp(-((t-1)./2000));
    T=T-(min(T)-GS.ts);
    %T=1e5*c./log10(1+t);
    %plot(t,T)
    %hold on
    %plot(t,T,'r')
    %pause
    T(GS.ts_num:GS.ts_num+GS.num_it)=GS.ts;
    t=1:length(T);

    %create directory structure for iterative saving inside of Heat Bath
    [year,month,day,hour,minute,second]=datevec(now);   %get the date
    dir_str=['TS' num2str(GS.ts) '_' num2str(year) num2str(month) num2str(day) num2str(hour) num2str(minute)];
    command_str=['mkdir ' dir_str];
    system(command_str);
    input_name=inputname(2);

    %Heat-Bath
    HeatBath(objfun,IN,T,dir_str,input_name);

else

    %*****************************
    %*** SAMPLING TEMPERATURE ***
    %****************************
    if isfield(GS,'ts_alg')
    else
        GS.ts_alg='ShortRun';
        fprintf('\n\n Sampling Temperature Algorithm Assigned to default: Short Run \n');
    end

    switch GS.ts_alg
        case 'ShortRun'
            fprintf('\n\nMODE: Short Run \n');
            
            %---------------------------------------------------------------------
            %SHORT-RUN PARAMETERS
            %---------------------------------------------------------------------
            if isfield(GS,'ShortRun')==1
            else %DEFAULT VALUES
                IN.GS.ShortRun.init_T=1e5;
                IN.GS.ShortRun.final_T=1e-1;
                IN.GS.ShortRun.num_T=10;
                IN.GS.ShortRun.num_R=5;
                IN.GS.ShortRun.num_sweep=25;
                GS=IN.GS;
            end
            
            %Temperature profile
            T=logspace(log10(GS.ShortRun.final_T),log10(GS.ShortRun.init_T),GS.ShortRun.num_T);
            T=fliplr(T);
            
            %Short-Run method of Basu and Frazer (1990)
            [ShortRun_dat.Ebar,ShortRun_dat.M,ShortRun_dat.E,ShortRun_dat.T,ShortRun_dat.R]=ShortRun(objfun,IN,T);	
            OUT.ShortRun=ShortRun_dat;
            OUT.date=date;
            OUT.IN=IN;
            OUT.input_name=inputname(2);
            subplot(2,3,6)
            sr_plot_env(OUT)
            OUT.ShortRun.BEST=sr_best_mod(OUT);
            return
        case 'FreezeProb'
            GS.T=T_freezeprob(IN);
            fprintf('\n\nMODE: FreezeProb \n');
            [T,p]=FreezeProb(objfun,IN); %Chapman (2001)
            OUT.T=T;
            OUT.p=p;
            OUT.date=date;
            return
        otherwise
            error('must specify .ts_alg field: ShortRun, FreezeProb')
    end
end







