inference.m

The package includes 3 Matlab-interfaces to the c-code: 1. inference.m An interface to the full

          else
              loopy.rho = rhoMat(args{i+1},adj);
              loopy.trw = 1;
          end
      case 'initMsg'
          if (algorithm~=LOOPY)
              if (algorithm~=GBP)
                  warning(['initMsg argument irrelevant for algorithm ' algorithm_name]);
              else
                  gbp.initMsg = args{i+1};
              end
          else
              loopy.initMsg = args{i+1};
          end
      case 'strategy'
          if (algorithm~=LOOPY)
              warning(['strategy argument irrelevant for algorithm ' algorithm_name]);
          else
              loopy.strategy = args{i+1};
          end
      case 'log'
          if (algorithm~=LOOPY && algorithm~=MEAN_FIELD && algorithm~=GBP)
              warning(['log argument irrelevant for algorithm ' algorithm_name]);
          else
              logsp = args{i+1};
          end
      case 'log_bels'
          if (algorithm~=LOOPY && algorithm~=MEAN_FIELD && algorithm~=GBP)
              warning(['log_bels argument irrelevant for algorithm ' algorithm_name]);
          else
              logbels = args{i+1};
          end
      case 'save_time'
          if (algorithm~=LOOPY)
              warning(['save_time argument irrelevant for algorithm ' algorithm_name]);
          else
              loopy.saveTime = args{i+1};
          end
      case 'trees_per_node'
          if (algorithm~=GBP)
              warning(['trees_per_node parameter irrelevant for algorithm ' algorithm_name]);
          else
              gbp.counting_node = sum(adj,1) - (args{i+1}(:))';
          end
      case 'full'
          if (algorithm~=GBP)
              warning(['full parameter irrelevant for algorithm ' algorithm_name]);
          else
              gbp.full = args{i+1};
          end
      case 'counting_node'
          if (algorithm~=GBP)
              if (algorithm~=LOOPY)
                  warning(['counting_node parameter irrelevant for algorithm ' algorithm_name]);
              else
                  loopy.counting_node = (args{i+1}(:))';
              end
          else
              gbp.counting_node = (args{i+1}(:))';
          end
      case 'alpha'
          if (algorithm~=GBP)
              warning(['alpha parameter irrelevant for algorithm ' algorithm_name]);
          else
              gbp.alpha = args{i+1};
          end
      case 'regions'
          if (algorithm~=GBP)
              warning(['regions argument irrelevant for algorithm ' algorithm_name]);
          else
              gbp.regions = args{i+1};
          end
      case 'all_regs'
          if (algorithm~=GBP)
              warning(['all_regs argument irrelevant for algorithm ' algorithm_name]);
          else
              gbp.all_regs = args{i+1};
          end
      case 'regions_potentials'
          if (algorithm~=GBP)
              warning(['regions_potentials argument irrelevant for algorithm ' algorithm_name]);
          else
              gbp.reg_pot = args{i+1};
          end
      case 'initial_state'
          if (monte_carlo.chosen == 0)
              warning(['initial_state argument irrelevant for algorithm ' algorithm_name]);
          else
              monte_carlo.initial_state = args{i+1};
              monte_carlo.initial_state_given = 1;
          end
      case 'burning_time'
          if (monte_carlo.chosen == 0)
              warning(['burning_time argument irrelevant for algorithm ' algorithm_name]);
          else
              monte_carlo.burning_time = args{i+1};
          end
      case 'sampling_interval'
          if (monte_carlo.chosen == 0)
              warning(['sampling_interval argument irrelevant for algorithm ' algorithm_name]);
          else
              monte_carlo.sampling_interval = args{i+1};
          end
      case 'num_samples'
          if (monte_carlo.chosen == 0)
              warning(['num_samples argument irrelevant for algorithm ' algorithm_name]);
          else
              monte_carlo.num_samples = args{i+1};
          end
      case 'max_iter'
          if (monte_carlo.chosen == 1)
              warning(['max_iter argument irrelevant for algorithm ' algorithm_name]);
          else
              maxIter = args{i+1};
          end
      case 'threshold'
          if (monte_carlo.chosen == 1)
              warning(['max_iter argument irrelevant for algorithm ' algorithm_name]);
          else
              th = args{i+1};
          end
      otherwise,
          error(['invalid argument name ' args{i}]);
  end
end

if (logsp)
    th = log(th);
end

adjCell = adjMat2Cell(adj);
psiCell = [];
if ~isempty(lambda)
	psiCell = psiSqueezeCell(lambda,adjCell,model);
end
clear lambda;
if isa(local,'sparse_cell')
    local = full(local);
end

belpairs = [];
switch algorithm,
    case LOOPY
        if (logsp == 0 && loopy.trw ~= 0)
            epsilon = 10^(-64);
            for i=1:length(psiCell)
                if iscell(psiCell{i})
                    for n=1:length(psiCell{i})
                        psiCell{i}{n} = psiCell{i}{n} / max(max(psiCell{i}{n}));
                        psiCell{i}{n} = psiCell{i}{n} .^ (1/loopy.rho{i}(n));
                        psiCell{i}{n}(find(psiCell{i}{n}<epsilon)) = epsilon;
                    end
                else
                    psiCell{i} = psiCell{i} ./ loopy.rho{i};
                end
            end
        end
        clear adj;
        if (nargout > 3)
            [bel, converged, belpairs_cform, msgs] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,...
                maxIter,th,logsp,logbels,loopy.sumOrMax,loopy.strategy,loopy.trw,loopy.rho,loopy.initMsg,loopy.saveTime,loopy.counting_node);
            belpairs = cform2SparseCell(belpairs_cform, adjCell);
        else
            if (nargout > 2)
                [bel, converged, belpairs_cform] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,...
                    maxIter,th,logsp,logbels,loopy.sumOrMax,loopy.strategy,loopy.trw,loopy.rho,loopy.initMsg,loopy.saveTime,loopy.counting_node);
                belpairs = cform2SparseCell(belpairs_cform, adjCell);
            else
                [bel, converged] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,maxIter,th,...
                    logsp,logbels,loopy.sumOrMax,loopy.strategy,loopy.trw,loopy.rho,loopy.initMsg,loopy.saveTime,loopy.counting_node);
            end
        end
    case GBP
        if (isempty(gbp.counting_node)==1)
            if gbp.trw>0
                gbp.counting_node = zeros(1,size(adj,1));
            end
        else
            if (isempty(gbp.regions)==0)
                warning('since regions are defined, the trees_per_node parameter will be ignored');
                gbp.counting_node = [];
            end
        end
        % if regions are not defined - pairs of neighbours are taken
        % (+isolated nodes as seperate regions)
        if (isempty(gbp.regions)==1)
            [I,J]=find(triu(adj));
            for i=1:length(I)
                gbp.regions{i} = [I(i), J(i)];
            end
            isolated = find(sum(adj)==0);
            for i=1:length(isolated)
                gbp.regions{end+1} = isolated(i);
            end
        end

        gbp.full = gbp.full || ~gbp.trw;

        clear adj;
        if (nargout > 3)
            [bel, converged, belpairs_cform,msgs] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,maxIter,th, ...
                logsp,logbels,gbp.regions,gbp.all_regs,gbp.sumOrMax,gbp.alpha,gbp.trw,gbp.counting_node,gbp.full,...
		gbp.initMsg,gbp.reg_pot,0);

            belpairs = cform2SparseCell(belpairs_cform, adjCell);
        else
            if (nargout > 2)
                [bel, converged, belpairs_cform] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,maxIter,th, ...
                	logsp,logbels,gbp.regions,gbp.all_regs,gbp.sumOrMax,gbp.alpha,gbp.trw,gbp.counting_node,gbp.full,...
			gbp.initMsg,gbp.reg_pot,0);

                belpairs = cform2SparseCell(belpairs_cform, adjCell);
            else
                [bel, converged] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,maxIter,th, ...
                	logsp,logbels,gbp.regions,gbp.all_regs,gbp.sumOrMax,gbp.alpha,gbp.trw,gbp.counting_node,gbp.full,...
			gbp.initMsg,gbp.reg_pot,0);
            end
        end
    case MEAN_FIELD
        clear adj;
        [bel, converged] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model,maxIter,th,logsp,logbels);
    case {GIBBS, WOLFF, SWENDSEN_WANG, METROPOLIS}
        if (monte_carlo.initial_state_given == 0)
            num_nodes = length(local);
            monte_carlo.initial_state = zeros(1,num_nodes);
            for i=1:num_nodes
                vi = length(local{i});
                q = ones(1,vi) / vi;
                monte_carlo.initial_state(i) = chooseInteger(q);
            end
        end
        clear adj;
        [bel, converged] = c_inference(adjCell,psiCell,local,algorithm-1,temperature,model, ...
            monte_carlo.initial_state,monte_carlo.burning_time, ...
            monte_carlo.sampling_interval,monte_carlo.num_samples);
    otherwise,
        error(['invalid algorithm ' algorithm_name]);
end