Moved refitting to constrain_mle local function

mnarayan · Jul 25, 2017 · 0b35a0a · 0b35a0a
1 parent a82c28b
commit 0b35a0a
Show file tree

Hide file tree

Showing 4 changed files with 58 additions and 42 deletions.
diff --git a/+estimator/penalized_mle_inverse_covariance.m b/+estimator/penalized_mle_inverse_covariance.m
@@ -64,26 +64,18 @@
             warning('Only lasso_mle supported') 
         end
 
-        if(~isempty(options.W))
-            warning('options.refit is turned on. Overwriting options.W');
-        end 
-        if(ndims(ThetaHat0)==2)       
-            options.W = options.lambda_max*100.*(ThetaHat0==0) + ...
-                                1.0.*(ThetaHat0~=0);
-        else
-           warning(['Need to perform model selection first. ...' ...
-                   ' Cannot refit to a path of estimates']);
-        end
-    end
+        [SigHat ThetaHat options] = constrain_mle(S,ThetaHat0,options);
 
-    switch options.estimator
+    else
 
-    case 'lasso_mle'
-       [SigHat ThetaHat] = lasso_mle(S,options);
-    otherwise
-       warning('Only lasso_mle supported') 
+        switch options.estimator
+
+        case 'lasso_mle'
+           [SigHat ThetaHat] = lasso_mle(S,options);
+        otherwise
+           warning('Only lasso_mle supported') 
+        end
     end
-
 
 end
 
@@ -108,9 +100,28 @@
 end
 
 
-function [SigHat ThetaHat] = constrain_mle(S,options)
+function [SigHat ThetaHat options] = constrain_mle(S,ThetaHat0,options)
+
 
 
+    if(~isempty(options.W))
+        warning('options.refit is turned on. Overwriting options.W');
+    end 
+    if(ndims(ThetaHat0)==2)       
+        options.W = options.lambda_max*100.*(ThetaHat0==0) + ...
+                            1.0.*(ThetaHat0~=0);
+    else
+       warning(['Need to perform model selection first. ...' ...
+               ' Cannot refit to a path of estimates']);
+    end
+
+    switch options.estimator
+
+    case 'lasso_mle'
+       [SigHat ThetaHat] = lasso_mle(S,options);
+    otherwise
+       warning('Only lasso_mle supported') 
+    end
 
 
 end

diff --git a/@GGM/adaptiveWeights.m b/@GGM/adaptiveWeights.m
@@ -3,40 +3,44 @@
 	narginchk(1,2); 
 	p = length(Theta);
 
-	if(nargin==3)
+	if(nargin>=2)
 		Lambda = varargin{1}; 
 	else
 		Lambda = .01*ones(p,p);
 	end
 	nonzeroW = 1.*(abs(Theta)>1e-5); 
 
-	A = diag(diag(Theta)) - Theta; 	
+    laplacian = @(A)(eye(length(A))- ...
+                        diag(1./sqrt(diag(A)))*A*diag(1./sqrt(diag(A))));  
+	A = laplacian(Theta);
 	t_start = cputime;
-	if(exist('mexeig'))
-		[V,D] = mexeig(A);		
-	else 
-		[V,D] = eig(A);
-		D = diag(D);  
-	end	
-	expA = V*diag(exp(D))*inv(V);			 
+    % if(exist('mexeig'))
+    %     [V,D] = mexeig(A);
+    % else
+    %     [V,D] = eig(A);
+    %     D = diag(D);
+    % end
+    % expA = V*diag(exp(D))*inv(V);
+    expA = expm(A);			 
 	timing = cputime-t_start;
 	disp(sprintf('expm took %.4f sec',timing)); 
 
-	% % Convert expA into a distance metric
-	% diagA = diag(expA);
-	% % A_i,i - A_i,j
-	% Gpos = bsxfun(@minus,diagA,expA);
-	% % A_j,j - A_j,i
-	% Gneg = bsxfun(@minus,diagA',expA');
-	% distA = (Gpos + Gneg)/2;
-	% sqrtD = 1./(diagA);
-	% distA = scale_cols(distA,sqrtD);
-	% distA = scale_rows(distA,sqrtD);
+    % % Convert expA into a distance metric
+    diagA = diag(expA);
+    % A_i,i - A_i,j
+    Gpos = bsxfun(@minus,diagA,expA);
+    % A_j,j - A_j,i
+    Gneg = bsxfun(@minus,diagA',expA');
+    distA = (Gpos + Gneg)/2;
+    sqrtD = 1./(diagA);
+    distA = scale_cols(distA,sqrtD);
+    distA = scale_rows(distA,sqrtD);
+    W = distA;
+
+    %distA = 1./abs(expA);    
+    %W = ((nonzeroW==1) + (nonzeroW==0).*(distA)).*Lambda; 
 
-	% distA = 1./abs(expA);		
-	% W = ((nonzeroW==1) + (nonzeroW==0).*(distA)).*Lambda;  
-
-	W = (nonzeroW==0).*Lambda;
+    % W = (nonzeroW==0).*Lambda;
 	W = setdiagLS(W,0); 
 
 end
diff --git a/@GGM/stableEstimate.m b/@GGM/stableEstimate.m
@@ -62,7 +62,7 @@
             end
 			tmpGGM = GGM(self.Data(:,:,1),1,0); 
             [tmpGGM mleresults] = tmpGGM.estimate();
-            mleresults.sparsemle.sparsity'
+            mleresults.sparsemle.sparsity';
 			self.ThetaPath = tmpGGM.ThetaPath;
 			self.Theta = tmpGGM.ThetaPath(:,:,best_lambda);			
 			warning('Change to take in alternative grid of lambdas')		

diff --git a/test/test_mle_inverse_covariance.m → test/test_penalized_mle_inverse_covariance.m b/test/test_mle_inverse_covariance.m → test/test_penalized_mle_inverse_covariance.m
@@ -27,6 +27,7 @@
 
     options.nlambdas = 1;
     options.Lambda = .1;
+    options.refit = true;
     results = estimator.fit(Data,options);