function [ R ] = experimentLHCRF( splits, params, H )

R = cell(1,numel(H));
best_accuracy = -1;
best_idx = -1;

for i=1:numel(H)
    params.nbHiddenStates = H{i}; 
    R{i}.params = params;
    
    % Create toolbox
    matHCRF('createToolbox', params.modelType, params.graphType, ...
        params.nbViews, params.nbHiddenStates, params.rawFeatureIndex);
    matHCRF('setOptimizer',params.optimizer);
    matHCRF('setInferenceEngine',params.inferenceEngine);
    matHCRF('set','regularizationL2',params.regFactorL2);
    matHCRF('set','maxIterations',params.maxIterOpt);
    matHCRF('set','randomSeed', params.seed); 

    % Train
    matHCRF('setData',splits.seqs_train,[],splits.labels_train);
    id=tic(); matHCRF('train'); time = toc(id);

    [R{i}.model,R{i}.features] = matHCRF('getModel');
    matHCRF('clearToolbox');

    % Load the model 
    matHCRF('createToolbox', R{i}.params.modelType, R{i}.params.graphType, ...
        R{i}.params.nbViews, R{i}.params.nbHiddenStates, ...
        R{i}.params.rawFeatureIndex);
    matHCRF('setOptimizer',R{i}.params.optimizer);
    matHCRF('setInferenceEngine',R{i}.params.inferenceEngine);
    matHCRF('initToolbox');
    matHCRF('setModel',R{i}.model,R{i}.features);

    % Test on validatation split
    pYstar = cell(1,numel(splits.labels_valid));
    for j=1:numel(splits.labels_valid), 
        matHCRF('setData',splits.seqs_valid(j),[],splits.labels_valid(j));
        matHCRF('test');
        ll = matHCRF('getResults');
        pYstar{j} = ll{1};
    end
    matHCRF('clearToolbox');
        
    [~,Ystar] = max(cell2mat(pYstar)); Ystar = Ystar-1;
    accuracy = sum(Ystar==splits.labels_valid)/numel(splits.labels_valid);
    if accuracy > best_accuracy
        best_idx = i;
        best_accuracy = accuracy;
    end    
    %fprintf('H=[%d %d], accuracy = %f, time = %.2f mins\n', ...
    %    R{i}.params.nbHiddenStates(1), R{i}.params.nbHiddenStates(2), accuracy, time/60);
        
    R{i}.accuracy_valid = accuracy;
    R{i}.time = time;
end

% Load the best model
 matHCRF('createToolbox', R{best_idx}.params.modelType, ...
     R{best_idx}.params.graphType, R{best_idx}.params.nbViews, ...
     R{best_idx}.params.nbHiddenStates, R{best_idx}.params.rawFeatureIndex);
matHCRF('setOptimizer',R{best_idx}.params.optimizer);
matHCRF('setInferenceEngine',R{best_idx}.params.inferenceEngine);
matHCRF('initToolbox');
matHCRF('setModel',R{best_idx}.model,R{best_idx}.features);

% Test 
pYstar = cell(1,numel(splits.labels_test));
for j=1:numel(splits.labels_test), 
    matHCRF('setData',splits.seqs_test(j),[],splits.labels_test(j));
    matHCRF('test');
    ll = matHCRF('getResults');
    pYstar{j} = ll{1};
end
matHCRF('clearToolbox');

% Report the result
[~,Ystar] = max(cell2mat(pYstar)); Ystar = Ystar-1;
accuracy = sum(Ystar==splits.labels_test)/numel(splits.labels_test);
fprintf('Best> Linked-HCRF [%d %d] accuracy = %f, time = %.2f mins\n', ...
    R{best_idx}.params.nbHiddenStates(1), R{best_idx}.params.nbHiddenStates(2),...
    accuracy, R{best_idx}.time/60);
end

