nInstances = 100; % Number of training examples
p = 15; % Number of variables
lambda = 500; % Strength of regularizer

% Make vector containing group memberships (note: Group 0 is reserved for non-regularized variables
groups = zeros(p,1);
groups(1:3:p) = [1:p/3]';
groups(2:3:p) = [1:p/3]';
groups(3:3:p) = [1:p/3]';
nGroups = length(unique(groups(groups>0)));

% Generate Data
X = randn(nInstances,p); % Generate Random Features
wTrue = randn(p,1).*groups; % Generate Random Weights, where higher groups have higher weight multiplier
y = X*wTrue + randn(nInstances,1); % Generate labels as linear combination plus random noise
funObj_sub = @(w)GaussianLoss(w,X,y); % We will be minimizing the Least Squares lost

% Code to optimize funObj_sub with L1-regularization
lambdaVect = lambda*ones(p,1); % Vector of regularization parameters (one for each variable)
funObj = @(w)nonNegGrad(w,lambdaVect,funObj_sub); % Surrogate objective function that divides variables into non-negative parts
wPosNeg = minConF_BC(funObj,zeros(2*p,1),zeros(2*p,1),inf(2*p,1)); % Solve non-negatively constrained optimization problem
w = wPosNeg(1:p)-wPosNeg(p+1:end) % Put non-negative parts together to get final anwser

% Code to optimize funObj_sub with Group L1-regularization (L1-regularization of group 2-norms)

% Make Auxiliary variable Objective and Projection Function
funObj = @(w)auxGroupLoss(w,groups,lambda,funObj_sub);
funProj = @(w)auxGroupL2Proj(w,groups);
 
% Solve
wAlpha = minConF_SPG(funObj,zeros(p+nGroups,1),funProj);
w = wAlpha(1:p) % Remove auxiliary variables

% Code to optimize funObj_sub with Group L1-regularization (L1-regularization of group infinity-norms)

% Make Auxiliary variable Objective and Projection Function
funObj = @(w)auxGroupLoss(w,groups,lambda,funObj_sub);
funProj = @(w)auxGroupLinfProj(w,groups);
 
% Solve
wAlpha = minConF_SPG(funObj,zeros(p+nGroups,1),funProj);
w = wAlpha(1:p) % Remove auxiliary variables

nTrain = 500; % number of examples to use for training
nTest = 1000; % number of examples to use for testing
nFeatures = 10; % number of features for each node
nNodes = 10; % number of nodes
nStates = 2; % number of states that each node can take

edgeProb = .5; % probability of each edge being included in the graph
edgeType = 1; % set to 0 to make the edge features normally distributed
ising = 0; % set to 1 to use ising potentials

trainType = 'pseudo'; % set to 'pseudo' for pseudo-likelihood, 'loopy' for loopy belief propagation, 'exact' for 'exact inference
testType = 'exact'; % set to 'loopy' to test with loopy belief propagation, 'exact' for exact inference
structureSeed = 0; % change this to generate different structures
trainSeed = 0; % vary seed from 0:9 to get paper results
useMex = 1; % use mex files in UGM to speed things up

% Regularization Parameters
lambdaNode = 10;
lambdaEdge = 10;

%% Generate data
nInstances = nTrain+nTest;
[y,adjTrue,X] = UGM_generate(nInstances,nFeatures,nNodes,edgeProb,nStates,ising,tied,useMex,edgeType);
perm = randperm(nInstances);
trainNdx = perm(1:nTrain);
testNdx = perm(nTrain+1:end);

%% Fixed Structure Methods 
edgePenaltyType = 'L2';

% Empty Structure
type = 'Fixed: Empty';
adjInit = zeros(nNodes);
example_UGMlearnSub;

%% Generative L1 Methods

% Regular L1
type = 'Generative L1: L1-L1';
edgePenaltyType = 'L1'; % Train with L1 regularization on the edges
Xold = X;
X = zeros(nInstances,0,nNodes); % Train while ignoring features
adjInit = fixed_Full(nNodes);
example_UGMlearnSub
 
adjInit = adjFinal;
X = Xold;
edgePenaltyType = 'L2';
subDisplay = display;
example_UGMlearnSub

%% Discriminative L1 Methods
type = 'Discriminative L1: L1-L1';
edgePenaltyType = 'L1'; % Train with L1 regularization on the edges
adjInit = fixed_Full(nNodes);
example_UGMlearnSub