diff --git a/inference/dimmwitted/src/sinco/AUTHORS b/inference/dimmwitted/src/sinco/AUTHORS deleted file mode 100644 index 7f5201903..000000000 --- a/inference/dimmwitted/src/sinco/AUTHORS +++ /dev/null @@ -1,6 +0,0 @@ -The author: Katya Scheinberg, -address: IBM T.J. Watson Research Center, - 1101 Kitchawan Rd., Yorktown Heights, - 10598, NY -email: katyascheinberg@gmail.com - \ No newline at end of file diff --git a/inference/dimmwitted/src/sinco/COPYING b/inference/dimmwitted/src/sinco/COPYING deleted file mode 100644 index c9990a7ea..000000000 --- a/inference/dimmwitted/src/sinco/COPYING +++ /dev/null @@ -1,213 +0,0 @@ -Common Public License Version 1.0 - -THE ACCOMPANYING PROGRAM IS PROVIDED UNDER THE TERMS OF THIS COMMON PUBLIC -LICENSE ("AGREEMENT"). ANY USE, REPRODUCTION OR DISTRIBUTION OF THE PROGRAM -CONSTITUTES RECIPIENT'S ACCEPTANCE OF THIS AGREEMENT. - -1. DEFINITIONS - -"Contribution" means: - - a) in the case of the initial Contributor, the initial code and -documentation distributed under this Agreement, and - - b) in the case of each subsequent Contributor: - - i) changes to the Program, and - - ii) additions to the Program; - - where such changes and/or additions to the Program originate from and are -distributed by that particular Contributor. A Contribution 'originates' from a -Contributor if it was added to the Program by such Contributor itself or anyone -acting on such Contributor's behalf. Contributions do not include additions to -the Program which: (i) are separate modules of software distributed in -conjunction with the Program under their own license agreement, and (ii) are not -derivative works of the Program. - -"Contributor" means any person or entity that distributes the Program. - -"Licensed Patents " mean patent claims licensable by a Contributor which are -necessarily infringed by the use or sale of its Contribution alone or when -combined with the Program. - -"Program" means the Contributions distributed in accordance with this Agreement. - -"Recipient" means anyone who receives the Program under this Agreement, -including all Contributors. - -2. GRANT OF RIGHTS - - a) Subject to the terms of this Agreement, each Contributor hereby grants -Recipient a non-exclusive, worldwide, royalty-free copyright license to -reproduce, prepare derivative works of, publicly display, publicly perform, -distribute and sublicense the Contribution of such Contributor, if any, and such -derivative works, in source code and object code form. - - b) Subject to the terms of this Agreement, each Contributor hereby grants -Recipient a non-exclusive, worldwide, royalty-free patent license under Licensed -Patents to make, use, sell, offer to sell, import and otherwise transfer the -Contribution of such Contributor, if any, in source code and object code form. -This patent license shall apply to the combination of the Contribution and the -Program if, at the time the Contribution is added by the Contributor, such -addition of the Contribution causes such combination to be covered by the -Licensed Patents. The patent license shall not apply to any other combinations -which include the Contribution. No hardware per se is licensed hereunder. - - c) Recipient understands that although each Contributor grants the licenses -to its Contributions set forth herein, no assurances are provided by any -Contributor that the Program does not infringe the patent or other intellectual -property rights of any other entity. Each Contributor disclaims any liability to -Recipient for claims brought by any other entity based on infringement of -intellectual property rights or otherwise. As a condition to exercising the -rights and licenses granted hereunder, each Recipient hereby assumes sole -responsibility to secure any other intellectual property rights needed, if any. -For example, if a third party patent license is required to allow Recipient to -distribute the Program, it is Recipient's responsibility to acquire that license -before distributing the Program. - - d) Each Contributor represents that to its knowledge it has sufficient -copyright rights in its Contribution, if any, to grant the copyright license set -forth in this Agreement. - -3. REQUIREMENTS - -A Contributor may choose to distribute the Program in object code form under its -own license agreement, provided that: - - a) it complies with the terms and conditions of this Agreement; and - - b) its license agreement: - - i) effectively disclaims on behalf of all Contributors all warranties and -conditions, express and implied, including warranties or conditions of title and -non-infringement, and implied warranties or conditions of merchantability and -fitness for a particular purpose; - - ii) effectively excludes on behalf of all Contributors all liability for -damages, including direct, indirect, special, incidental and consequential -damages, such as lost profits; - - iii) states that any provisions which differ from this Agreement are offered -by that Contributor alone and not by any other party; and - - iv) states that source code for the Program is available from such -Contributor, and informs licensees how to obtain it in a reasonable manner on or -through a medium customarily used for software exchange. - -When the Program is made available in source code form: - - a) it must be made available under this Agreement; and - - b) a copy of this Agreement must be included with each copy of the Program. - -Contributors may not remove or alter any copyright notices contained within the -Program. - -Each Contributor must identify itself as the originator of its Contribution, if -any, in a manner that reasonably allows subsequent Recipients to identify the -originator of the Contribution. - -4. COMMERCIAL DISTRIBUTION - -Commercial distributors of software may accept certain responsibilities with -respect to end users, business partners and the like. While this license is -intended to facilitate the commercial use of the Program, the Contributor who -includes the Program in a commercial product offering should do so in a manner -which does not create potential liability for other Contributors. Therefore, if -a Contributor includes the Program in a commercial product offering, such -Contributor ("Commercial Contributor") hereby agrees to defend and indemnify -every other Contributor ("Indemnified Contributor") against any losses, damages -and costs (collectively "Losses") arising from claims, lawsuits and other legal -actions brought by a third party against the Indemnified Contributor to the -extent caused by the acts or omissions of such Commercial Contributor in -connection with its distribution of the Program in a commercial product -offering. The obligations in this section do not apply to any claims or Losses -relating to any actual or alleged intellectual property infringement. In order -to qualify, an Indemnified Contributor must: a) promptly notify the Commercial -Contributor in writing of such claim, and b) allow the Commercial Contributor to -control, and cooperate with the Commercial Contributor in, the defense and any -related settlement negotiations. The Indemnified Contributor may participate in -any such claim at its own expense. - -For example, a Contributor might include the Program in a commercial product -offering, Product X. That Contributor is then a Commercial Contributor. If that -Commercial Contributor then makes performance claims, or offers warranties -related to Product X, those performance claims and warranties are such -Commercial Contributor's responsibility alone. Under this section, the -Commercial Contributor would have to defend claims against the other -Contributors related to those performance claims and warranties, and if a court -requires any other Contributor to pay any damages as a result, the Commercial -Contributor must pay those damages. - -5. NO WARRANTY - -EXCEPT AS EXPRESSLY SET FORTH IN THIS AGREEMENT, THE PROGRAM IS PROVIDED ON AN -"AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, EITHER EXPRESS OR -IMPLIED INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OR CONDITIONS OF TITLE, -NON-INFRINGEMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Each -Recipient is solely responsible for determining the appropriateness of using and -distributing the Program and assumes all risks associated with its exercise of -rights under this Agreement, including but not limited to the risks and costs of -program errors, compliance with applicable laws, damage to or loss of data, -programs or equipment, and unavailability or interruption of operations. - -6. DISCLAIMER OF LIABILITY - -EXCEPT AS EXPRESSLY SET FORTH IN THIS AGREEMENT, NEITHER RECIPIENT NOR ANY -CONTRIBUTORS SHALL HAVE ANY LIABILITY FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING WITHOUT LIMITATION LOST -PROFITS), HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, -STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY -OUT OF THE USE OR DISTRIBUTION OF THE PROGRAM OR THE EXERCISE OF ANY RIGHTS -GRANTED HEREUNDER, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. - -7. GENERAL - -If any provision of this Agreement is invalid or unenforceable under applicable -law, it shall not affect the validity or enforceability of the remainder of the -terms of this Agreement, and without further action by the parties hereto, such -provision shall be reformed to the minimum extent necessary to make such -provision valid and enforceable. - -If Recipient institutes patent litigation against a Contributor with respect to -a patent applicable to software (including a cross-claim or counterclaim in a -lawsuit), then any patent licenses granted by that Contributor to such Recipient -under this Agreement shall terminate as of the date such litigation is filed. In -addition, if Recipient institutes patent litigation against any entity -(including a cross-claim or counterclaim in a lawsuit) alleging that the Program -itself (excluding combinations of the Program with other software or hardware) -infringes such Recipient's patent(s), then such Recipient's rights granted under -Section 2(b) shall terminate as of the date such litigation is filed. - -All Recipient's rights under this Agreement shall terminate if it fails to -comply with any of the material terms or conditions of this Agreement and does -not cure such failure in a reasonable period of time after becoming aware of -such noncompliance. If all Recipient's rights under this Agreement terminate, -Recipient agrees to cease use and distribution of the Program as soon as -reasonably practicable. However, Recipient's obligations under this Agreement -and any licenses granted by Recipient relating to the Program shall continue and -survive. - -Everyone is permitted to copy and distribute copies of this Agreement, but in -order to avoid inconsistency the Agreement is copyrighted and may only be -modified in the following manner. The Agreement Steward reserves the right to -publish new versions (including revisions) of this Agreement from time to time. -No one other than the Agreement Steward has the right to modify this Agreement. -IBM is the initial Agreement Steward. IBM may assign the responsibility to serve -as the Agreement Steward to a suitable separate entity. Each new version of the -Agreement will be given a distinguishing version number. The Program (including -Contributions) may always be distributed subject to the version of the Agreement -under which it was received. In addition, after a new version of the Agreement -is published, Contributor may elect to distribute the Program (including its -Contributions) under the new version. Except as expressly stated in Sections -2(a) and 2(b) above, Recipient receives no rights or licenses to the -intellectual property of any Contributor under this Agreement, whether -expressly, by implication, estoppel or otherwise. All rights in the Program not -expressly granted under this Agreement are reserved. - -This Agreement is governed by the laws of the State of New York and the -intellectual property laws of the United States of America. No party to this -Agreement will bring a legal action under this Agreement more than one year -after the cause of action arose. Each party waives its rights to a jury trial in -any resulting litigation. diff --git a/inference/dimmwitted/src/sinco/INSTALL b/inference/dimmwitted/src/sinco/INSTALL deleted file mode 100644 index 2596cc13c..000000000 --- a/inference/dimmwitted/src/sinco/INSTALL +++ /dev/null @@ -1,3 +0,0 @@ -SINCO software comex with Matlab mex interface. To compile -run compile_sinco comand (file provided) in Matlab and then -use the appropriate calling line in your Matlab code. \ No newline at end of file diff --git a/inference/dimmwitted/src/sinco/a.out b/inference/dimmwitted/src/sinco/a.out deleted file mode 100755 index cb84057d9..000000000 Binary files a/inference/dimmwitted/src/sinco/a.out and /dev/null differ diff --git a/inference/dimmwitted/src/sinco/compile_sinco.m b/inference/dimmwitted/src/sinco/compile_sinco.m deleted file mode 100644 index 963ad38e7..000000000 --- a/inference/dimmwitted/src/sinco/compile_sinco.m +++ /dev/null @@ -1,27 +0,0 @@ -% Script to compile mex files for SINCO using MATLAB interface -% -% -% Caution!!!! Only tested for Linux right now. - -% Test for mex extension to determine matlab version and platform -switch mexext - case {'mexw32'} % Win32 MATLAB after 7.1 - files='sinco.cpp sinco_mex.cpp'; -% libs=[matlabroot,'\extern\lib\win32\microsoft\libmwlapack.lib']; - switches='-v -O -DWIN32'; - eval(['mex ',switches,' ',files,' ','''']); - case {'dll'} % Win32 MATLAB before 7.1 - files='sinco.cpp sinco_mex.cpp'; -% libs=[matlabroot,'\extern\lib\win32\microsoft\msvc60\libmwlapack.lib']; - switches='-v -O -DWIN32'; - eval(['mex ',switches,' ',files,' ','''']); - case {'mexmac'}% Macintosh using the VecLib framework - files='cinco.cpp sinco_mex.cpp'; - switches='-v -O -Dmac'; - eval(['mex ',switches,' ',files,' ']); - otherwise % All other platforms - files='sinco.cpp sinco_mex.cpp'; - switches='-v -O -Dlinuxp'; - eval(['mex ',switches,' ',files,' ']); -end -disp(' ......................... Done Compiling Source .........................') diff --git a/inference/dimmwitted/src/sinco/sinco.cpp b/inference/dimmwitted/src/sinco/sinco.cpp deleted file mode 100644 index 4d1d3395d..000000000 --- a/inference/dimmwitted/src/sinco/sinco.cpp +++ /dev/null @@ -1,437 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -#include "sinco.hpp" - -/* This is the C++ implementation to solve the sparse inverse covariance - problem of the form - -max _C[ K*log(det(C)-tr(AC)-\lambda ||S.*C||_1] - -K and lambda are positive scalars, A is a given symmetric matrix, S is a given nonnegative matrix -(not necessarily symmetric), S.*C is a elemwhise product and ||*||_1 is the sum -of all absolute values. - - */ - -double funvalue_update(double alpha, double K, double Wij, double Wii, double Wjj, - double Aij, double Sij,double Sji, double update) { - double detratio, fchange; - /* This routine computes the change in the objective function value when a step of length\ alpha is made in the direction e_ie_j^T+e_je_i^T */ -if (update ==1) { - detratio=(1+alpha*Wij)*(1+alpha*Wij-alpha*alpha*Wii*Wjj/(1+alpha*Wij)); - fchange=K*(log(detratio))-2*alpha*Aij-alpha*(Sij+Sji); - } - else { - detratio=(1-alpha*Wij)*(1-alpha*Wij-alpha*alpha*Wii*Wjj/(1-alpha*Wij)); - fchange=K*(log(detratio))+2*alpha*Aij-alpha*(Sij+Sji); - } - return fchange; -} - - -void invupdate(Matrix& Gradp, Matrix& Gradpp, vector& UpInd, Matrix& W, const double theta, const double K, const int p, int i, int j, int update) -{ - /* This routine computes the change in the dual matrix W (inverse of C) when a step of length\ alpha is made in the direction e_ie_j^T+e_je_i^T */ -double a,b,c,d, wij, wjj, wii; - int ii, jj, upindind; -Matrix UpW(p,p); -wii=W(i,i); -wjj=W(j,j); -wij=W(i,j); -if (update==1){ - d=theta*theta*(wii*wjj-wij*wij)-1-2*theta*wij; - a=-(1+theta*wij)/d; - b=theta*wjj/d; - c=theta*wii/d; - for (int ii=p-1; ii>=0; ii--) { - for (int jj=ii; jj>=0; jj--){ - upindind=ii*p+jj; - UpW(ii,jj)=-theta*(a*W(i,ii)*W(j,jj)+ c*W(j,ii)*W(j,jj)+b*W(i,ii)*W(i,jj)+a*W(j,ii)*W(i,jj)); - if (fabs(UpW(ii,jj))>zerotol){UpInd[upindind]=1;} - else {UpInd[ii*p+jj]=0;} - UpW(jj,ii)=UpW(ii,jj); - UpInd[jj*p+ii]=UpInd[upindind]; - } - } -} - else{ - d=theta*theta*(-wii*wjj+wij*wij)+1-2*theta*wij; - a=(1-theta*wij)/d; - b=theta*wjj/d; - c=theta*wii/d; - - for (int ii=p-1; ii>=0; ii--) { - for (int jj=ii; jj>=0; jj--) { - upindind=ii*p+jj; - UpW(ii,jj)=theta*(a*W(i,ii)*W(j,jj)+ c*W(j,ii)*W(j,jj)+b*W(i,ii)*W(i,jj)+a*W(j,ii)*W(i,jj)); - if (fabs(UpW(ii,jj))>zerotol){UpInd[upindind]=1;} - else {UpInd[ii*p+jj]=0;} - UpW(jj,ii)=UpW(ii,jj); - UpInd[jj*p+ii]=UpInd[upindind]; - } - } -} -W.sum(UpW); -Gradp.sumwithscal(K,UpW); -Gradpp.sumwithscal(-K,UpW); - -} - - - -double findposstep(double K, double Wij, double Wii, double Wjj, - double Aij, double Sij, double Sji, int update) { - /* This routine computes the optimal length of a step of length in the direction e_ie_j^T+e_je_i^T for the positive component of C */ - double aux1, aux2, aux3, aux4, aux5; - double a, b, c, D, alpha, alpha1, alpha2; - - if ( update ==1 ) { - aux1=2*(K*Wij-Aij)-Sji-Sij; - aux2=(Wii*Wjj-Wij*Wij); - aux3=2*Wij; - aux4=-2*K*(Wii*Wjj+Wij*Wij); - aux5=2*K*Wij*aux2; - } - else { - aux1=2*(-K*Wij+Aij)-Sij-Sji; - aux2=(-Wii*Wjj+Wij*Wij); - aux3=2*Wij; - aux4=2*K*(Wii*Wjj+Wij*Wij); - aux5=-2*K*Wij*aux2; - } - a=aux2*aux1-aux5; - b=-aux3*aux1-aux4; - c=-update*aux1; - if (fabs(a)>zerotol) { - D=b*b-4*a*c; - if (D<0) { - printf ( "negative discriminant, \n"); - abort; - } - alpha1=fmin((-b-sqrt(D))/(2*a),(-b+sqrt(D))/(2*a)); - alpha2=fmax((-b-sqrt(D))/(2*a),(-b+sqrt(D))/(2*a)); - if (alpha1>=0){ alpha=alpha1;} - else { - if (alpha2>=0){ alpha=alpha2;} - else { - printf ( "unbounded direction!!!, \n"); - abort; - } - } - } - else if (-c/b>0) - { alpha=-c/b;} - else { - printf ( "unbounded direction!!!,\n"); - abort; - } - - /* Check correctness -double theta=alpha; - double fprime; - if ( update ==1 ) { - fprime=(K*Wij+K*Wij-Aij-Aij-Sij-Sji)- - K*theta/(theta*theta*(Wii*Wjj-Wij*Wij)-1-2*theta*Wij)*(-2*(Wii*Wjj+Wij*Wij) - +2*theta*Wij*(Wii*Wjj-Wij*Wij)); - } - else { - fprime=(-K*Wij-K*Wij+Aij+Aij-Sij-Sji)-K*theta/(theta*theta*(-Wii*Wjj+Wij*Wij) - +1-2*theta*Wij)*(2*(Wii*Wjj+Wij*Wij) - +2*theta*Wij*(Wii*Wjj-Wij*Wij)); - } - if (fabs(fprime) > 10e-6){ - printf("\n bad gradient"); - fflush(stdout); - abort; - } */ - return alpha; -} - - -double findnegstep(double K, bool diag, double Wij, double Wii, double Wjj, - double Aij, double Sij, double Sji, double Cpij, double Cppij, int update) { -/* This routine computes the optimal length of a step of length in the direction e_ie_j^T+e_je_i^T for the negative component of C */ - double aux1, aux2, aux3, aux4, aux5; - double a, b, c, D, maxstep, alpha, alpha1, alpha2; - - if (update ==1) { - aux1=2*(K*Wij-Aij)-Sji-Sij; - aux2=(Wii*Wjj-Wij*Wij); - aux3=2*Wij; - aux4=-2*K*(Wii*Wjj+Wij*Wij); - aux5=2*K*Wij*aux2; - maxstep=Cpij;} - else { - aux1=2*(-K*Wij+Aij)-Sij-Sji; - aux2=(-Wii*Wjj+Wij*Wij); - aux3=2*Wij; - aux4=2*K*(Wii*Wjj+Wij*Wij); - aux5=-2*K*Wij*aux2; - maxstep=Cppij; -} -if (diag) { - maxstep=maxstep/2; -} -a=aux2*aux1-aux5; -b=-aux3*aux1-aux4; -c=-update*aux1; - if (fabs(a)>zerotol) { - D=b*b-4*a*c; - if (D<0) { - printf ( "negative discriminant, \n"); - abort; - } - double sqrootD=sqrt(D); - alpha1=fmin(((-b-sqrootD)/(2*a)),((-b+sqrootD)/(2*a))); - alpha2=fmax(((-b-sqrootD)/(2*a)),((-b+sqrootD)/(2*a))); - if ( alpha2<0 ){ - alpha=fmax(alpha2, -maxstep);} - else if ( alpha1<0 ){ - alpha=fmax(alpha1, -maxstep);} - else { - alpha=-maxstep; - } - } - else if (-c/b<0){ - alpha=fmax(-c/b, -maxstep);} - else { - alpha=-maxstep; - } - /* Check correctness - -double theta=alpha; - double fprime; - if ( update ==1 ) { - fprime=(K*Wij+K*Wij-Aij-Aij-Sij-Sji)- - K*theta/(theta*theta*(Wii*Wjj-Wij*Wij)-1-2*theta*Wij)*(-2*(Wii*Wjj+Wij*Wij) - +2*theta*Wij*(Wii*Wjj-Wij*Wij)); - } - else { - fprime=(-K*Wij-K*Wij+Aij+Aij-Sij-Sji)-K*theta/(theta*theta*(-Wii*Wjj+Wij*Wij) - +1-2*theta*Wij)*(2*(Wii*Wjj+Wij*Wij) - +2*theta*Wij*(Wii*Wjj-Wij*Wij)); - } - if ((abs(theta) 10e-6)){ - printf("\n bad gradient"); - fflush(stdout); - abort; - } */ - return alpha; -} - - -void ICS::sinco(SOL& ICS_sol, const SOL& ICS_start, double tol ){ - /*This is the main routine for computing the sparse inverse covariance. - ICS class contains the problem data, that is A, K, p, lambda and S. - The input for this routine is a starting point ICS_start which contains - intial C, initial W=C^{-1} and initial objective function value. The output -of the routine is final C, W and function value. Tolerance is set in tol and -if set too high may significantly slow down the solver, since convergence - is slow in the end. */ - vector UpInd; - double alphamax, funmax, K, fnew; - int imax, jmax, updatemax, iter=0; - UpInd.resize(p*p, 1); - bool stop=false; - double f=ICS_start.f; - Steps AllSteps(p,p); - K=N/2; - ICS_sol.W=ICS_start.W; - ICS_sol.C=ICS_start.C; - Matrix Cp(p,p,0); - Matrix Cpp(p,p,0); - Matrix Gradp(p,p,0); - Matrix Gradpp(p,p,0); - for (int i=p-1; i>=0; i--) { - for (int j=i; j>=0; j--) { - if (ICS_start.C(i,j)>0) { - Cp(i,j)=ICS_start.C(i,j); - Cp(j,i)=ICS_start.C(j,i); } - else { - Cpp(i,j)=-ICS_start.C(i,j); - Cpp(j,i)=-ICS_start.C(j,i); - } - } - } - Gradp.sumwithscal(-1*lambda, S); - Gradp.sumwithscal(-1, A); - Gradp.sumwithscal(K,ICS_sol.W); - Gradpp.sumwithscal(-1*lambda, S); - Gradpp.sum(A); - Gradpp.sumwithscal(-1*K,ICS_sol.W); - - - while (!stop && (iter<100000)) { - // %find the largest positive derivative - alphamax=0; - funmax=f; - iter+=1; - // printf("\n iteration [%4d] [%6f]", iter, f); - // fflush(stdout); - - for (int i=p-1; i>=0; i--) { - for (int j=i; j>=0; j--) { - if (UpInd[i*p+j] || UpInd[i*p+i] || UpInd[j*p+j]){ - AllSteps(i,j).alpha=0; - if ((Gradp(i,j)+Gradp(j,i) > tol) && (Cpp(i,j)<=tol)){ - AllSteps(i,j).update=1; - AllSteps(i,j).alpha=findposstep(K, ICS_sol.W(i,j),ICS_sol.W(i,i), - ICS_sol.W(j,j), A(i,j), lambda*S(i,j), - lambda*S(j,i), AllSteps(i,j).update); - } - else if ((Gradpp(i,j)+Gradpp(j,i) > tol) && (Cp(i,j)<=tol)) { - AllSteps(i,j).update=-1; - AllSteps(i,j).alpha=findposstep(K, ICS_sol.W(i,j),ICS_sol.W(i,i), - ICS_sol.W(j,j), A(i,j), lambda*S(i,j), - lambda*S(j,i), AllSteps(i,j).update); - } - else if ((Gradp(i,j)+Gradp(j,i)<-tol) && (Cp(i,j)> tol )) { - AllSteps(i,j).update=1; - AllSteps(i,j).alpha=findnegstep(K, (i==j), ICS_sol.W(i,j),ICS_sol.W(i,i), - ICS_sol.W(j,j), A(i,j), lambda*S(i,j), - lambda*S(j,i), Cp(i,j), Cpp(i,j), - AllSteps(i,j).update); - } - else if ((Gradpp(i,j)+Gradpp(j,i)<-tol) && (Cpp(i,j)> tol)) { - AllSteps(i,j).update=-1; - AllSteps(i,j).alpha=findnegstep(K, (i==j), ICS_sol.W(i,j),ICS_sol.W(i,i), - ICS_sol.W(j,j), A(i,j), lambda*S(i,j), - lambda*S(j,i), Cp(i,j), Cpp(i,j), AllSteps(i,j).update); - } - if ( fabs(AllSteps(i,j).alpha)> tol) { - AllSteps(i,j).fchange=funvalue_update(AllSteps(i,j).alpha, - K,ICS_sol.W(i,j),ICS_sol.W(i,i), - ICS_sol.W(j,j), A(i,j), lambda*S(i,j), - lambda*S(j,i), AllSteps(i,j).update); - } - else { - AllSteps(i,j).fchange=0; - } - } - fnew=f+AllSteps(i,j).fchange; - if (fnew > funmax ) { - funmax=fnew; - imax=i; - jmax=j; - alphamax=AllSteps(i,j).alpha; - updatemax=AllSteps(i,j).update; - } - } - } - /* for (int i=0; i>test.A(i,j); - } - } - /* myfile.close(); */ - /* for (int i=0; i -#include -#include -#include -#include -#include -#include -#include - -#define zerotol 1e-10 -using namespace std; - - -class Matrix{ -public: - Matrix(){} - Matrix(int d1, int d2, double elm=0){ - M.resize(d1*d2, elm); - dim1=d1; - dim2=d2; - } - Matrix(const Matrix& A){ - M=A.M; - dim1=A.dim1; - dim2=A.dim2; - } - Matrix& operator=(const Matrix& A){ - if (this !=&A){ - M=A.M; - dim1=A.dim1; - dim2=A.dim2; - } - return *this; - } - - inline double operator()(int i, int j) const { - return M[i*dim2+j]; - } - inline double& operator()(int i, int j){ - return M[i*dim2+j]; - } - - double dotmultiply (Matrix A){ - double S=0; - for (int i=dim1*dim2-1; i>=0; i--) { - S+=A.M[i]*M[i]; - } - return S; - } - void importMat(double* A, int d1, int d2) { - M.resize(d1*d2); - for (int i=d1*d2-1; i>=0; i--) { - M[i]=A[i]; - } - dim1=d1; - dim2=d2; - } - - void exportMat(double* A) { - for (int i=dim1*dim2-1; i>=0; i--) { - A[i]=M[i]; - } - } - - - - /* double dotmultiply (Matrix& A){ - double S=0; - for (int i=0; i=0; i--) { - M[i]+=A.M[i]; - } - } - void sumwithscal(double N, Matrix& A){ - for (int i=dim1*dim2-1; i>=0; i--) { - M[i]+=N*A.M[i]; - } - } - void addidentity(double N){ - for (int i=0; i=0; i--) { - S+=fabs(M[i]*A.M[i]); - } - return S; - } - /* double normS (Matrix& A){ - double S=0; - for (int i=0; i M; - int dim1; - int dim2; -}; - -class Step{ -public: - double alpha; - double fchange; - int update; -}; - - -class Steps{ -public: - Steps(int d1, int d2){ - Stepmat.resize(d1*d2); - dim1=d1; - dim2=d2; - } - inline Step& operator()(int i, int j){ - return Stepmat[i*dim2+j]; - } - vector< Step > Stepmat; - int dim1; - int dim2; -}; - -class SOL{ -public: - Matrix C; - Matrix W; - double f; -}; - -/* Class ICS is a class of "Inverse Covariance Selection" problems, it has the optimization function, the data A, lambda, S, p and N. */ - -class ICS{ -public: - void sinco(SOL& ICS_sol, const SOL& ICS_start, double tol ); - Matrix A; - Matrix S; - double lambda; - int p; - int N; - -}; - - diff --git a/inference/dimmwitted/src/sinco/sinco_mex.cpp b/inference/dimmwitted/src/sinco/sinco_mex.cpp deleted file mode 100644 index 210ff35d1..000000000 --- a/inference/dimmwitted/src/sinco/sinco_mex.cpp +++ /dev/null @@ -1,56 +0,0 @@ -// MEX wrapper for the function computing sparse inv. covariance - -// function [Csol, Wsol, fsol]=sinco(Cstart, fstart, Wstart, A, S, lambda, p, K); - -// Calling sequence: void ICS::sinco(SOL& ICS_sol, const SOL& ICS_start, double tol ){ - -// Last Modified: Katya Scheinber Feb 2009 -// - - -#include "sinco.hpp" -#include "mex.h" - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - double *Cstart, *Wstart, *Csol, *Wsol, *A, *S, *fsol; - double fstart, K, lambda, tol; - int p; - ICS prob; - SOL Init_Sol; - SOL Opt_Sol; - - Cstart = mxGetPr(prhs[0]); - p = mxGetN(prhs[0]); - Wstart = mxGetPr(prhs[2]); - A = mxGetPr(prhs[3]); - S = mxGetPr(prhs[4]); - fstart = mxGetScalar(prhs[1]); - lambda = mxGetScalar(prhs[5]); - K = mxGetScalar(prhs[6]); - tol = mxGetScalar(prhs[7]); - - plhs[0] = mxCreateDoubleMatrix(p, p, mxREAL); - plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL); - plhs[1] = mxCreateDoubleMatrix(p, p, mxREAL); - Csol = mxGetPr(plhs[0]); - fsol = mxGetPr(plhs[2]); - Wsol = mxGetPr(plhs[1]); - - prob.p=p; - prob.S.importMat(S, p, p); - prob.A.importMat(A, p, p); - prob.lambda=lambda; - prob.N=K; - Init_Sol.C.importMat(Cstart, p, p); - Init_Sol.W.importMat(Wstart, p, p); - Init_Sol.f=fstart; - - prob.sinco(Opt_Sol, Init_Sol, tol ); - - Opt_Sol.C.exportMat(Csol); - Opt_Sol.W.exportMat(Wsol); - *fsol=Opt_Sol.f; - -} -