1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366 |
-
-
-
-
-
-
- #ifndef LOGBARRIERCG_INC
- #define LOGBARRIERCG_INC
-
- #include <Eigen/IterativeLinearSolvers>
- #include <iostream>
- #include <iomanip>
- #include <fstream>
- #include <Eigen/Eigen>
- #include "cg.h"
- #include "bicgstab.h"
- #include "lemma.h"
-
-
- namespace Lemma {
-
- template < typename cgScalar >
- cgScalar PhiB (const cgScalar& mux, const cgScalar& muy, const cgScalar& minVal,
- const cgScalar& maxVal, const VectorXr x) {
- cgScalar phib = std::abs((x.array() - minVal).log().sum()*mux);
-
- return phib;
- }
-
-
-
- template < typename cgScalar >
- cgScalar PhiB2 (const cgScalar& mux, const cgScalar& muy, const cgScalar& minVal,
- const cgScalar& maxVal, const VectorXr x, const int& block,
- const int &nblocks) {
- cgScalar phib = std::abs((x.array() - minVal).log().sum()*mux);
-
- for (int ib=0; ib<nblocks; ++ib) {
-
- phib += cgScalar(block)*std::log(maxVal - x.segment(ib*block, block).sum())*muy;
- }
- return phib;
- }
-
- template < typename cgScalar >
- cgScalar PhiB2 (const cgScalar& minVal, const cgScalar& maxVal, const VectorXr x,
- const int& block, const int &nblocks) {
- cgScalar phib = std::abs((x.array() - minVal).log().sum());
-
- for (int ib=0; ib<nblocks; ++ib) {
-
- phib += cgScalar(block)*std::log(maxVal - x.segment(ib*block, block).sum());
- }
- return phib;
- }
-
- template < typename cgScalar >
- cgScalar PhiB2_NN (const cgScalar& mux, const cgScalar& minVal, const VectorXr x) {
- cgScalar phib = std::abs((x.array() - minVal).log().sum()*mux);
- return phib;
- }
-
-
-
-
- template < typename cgScalar >
- VectorXr LogBarrierCG_4(const MatrixXr &A,
- const VectorXr &b,
- const VectorXr &xr,
- const VectorXr &x0,
- const Eigen::SparseMatrix<cgScalar>& WdTWd,
- const Eigen::SparseMatrix<cgScalar>& WmTWm,
- const cgScalar &minVal,
- const cgScalar &maxVal,
- const Real& lambda) {
-
-
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
-
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
- ATWdTWdA += lambda*WmTWm;
-
- int N = A.cols();
- int M = A.rows();
- cgScalar SIGMA = .25;
-
- cgScalar epsilon = 1e-55;
- Real delta = minVal + 1e-8;
-
- VectorXr x = VectorXr::Zero(N).array() + epsilon;
- VectorXr b_pre = VectorXr::Zero(M);
-
- cgScalar phid = (b_pre - b).norm();
- cgScalar phim = x.norm();
- cgScalar mux = (phid + lambda*phim) * 1e4;
- cgScalar muy = 1e-30;
- cgScalar phib = PhiB(mux, muy, minVal, maxVal, x);
-
-
- int iloop(0);
- do {
-
-
-
- VectorXr X1 = VectorXr::Ones(N).array() / ((x0+x).array()-minVal) ;
- VectorXr Y1 = VectorXr::Ones(N).array() / (maxVal-(x0+x).array()) ;
-
- VectorXr X2 = VectorXr::Ones(N).array() / (((x0+x).array()-minVal)*((x0+x).array()-minVal));
- VectorXr Y2 = VectorXr::Ones(N).array() / ((maxVal-(x0+x).array())*(maxVal-(x0+x).array()));
-
-
-
-
-
-
- VectorXr b2 = (A.transpose()*WdTWd*(b)).array()
-
- + 2.*mux*X1.array() + 2.*muy*Y1.array();
-
-
- MatrixXr A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
- A2.diagonal().array() += mux*X2.array() + muy*Y2.array();
-
-
-
- Eigen::ConjugateGradient<MatrixXr> CG;
-
- CG.compute(A2);
-
- VectorXr ztilde = CG.solve(b2);
-
-
- VectorXr h = ztilde - x;
- cgScalar d(.25);
- for (int ix=0; ix<x.size(); ++ix) {
- d = std::min(d, (cgScalar).925*((x0(ix)+x[ix]-minVal)/std::abs(h[ix])));
- d = std::min(d, (cgScalar).925*(std::abs(x0(ix)+x[ix]-maxVal)/std::abs(h[ix])));
- }
-
- x += d*h;
- std::cout << "\titerations " << CG.iterations() << std::endl;
- std::cout << "\ttolerance " << CG.tolerance() << std::endl;
- std::cout << "\td " << d << std::endl;
-
- for (int i=0; i<x.size(); ++i) {
- if (x(i)+x0(i) < minVal) {
- std::cerr << "overstepp\t" << x(i) << "\t" << x0(i) << "\t" << x(i)+x0(i) << "\t" << minVal << std::endl;
-
- x(i) = minVal + delta;
- exit(1);
- } else if (x(i)+x0(i) > maxVal) {
- std::cout << "overstepp BIG\t" << x(i) << "\t" << x0(i) << "\t" << maxVal<< std::endl;
- x(i) = maxVal - x0(i) - delta;
- exit(1);
- }
- }
-
-
-
- VectorXr s1 = mux * (X2.asDiagonal()*h - 2.*X1);
- VectorXr s2 = muy * (Y2.asDiagonal()*h - 2.*Y1);
-
-
-
-
-
- mux = SIGMA/((cgScalar)(N)) * std::abs( s1.dot(x) ) ;
- muy = SIGMA/((cgScalar)(N)) * std::abs( s2.dot(x) ) ;
-
- b_pre = A*x;
- phid = (b_pre - b).norm();
-
- phim = (WmTWm*(x-xr)).norm();
- phib = PhiB(mux, muy, minVal, maxVal, x0+x);
- ++iloop;
-
- std::cout.precision(14);
- std::cout << iloop << "\t" << phid << "\t" << phim << "\t" << mux << "\t" << muy << "\t" << phib<< std::endl;
-
-
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
- std::cout << "final x\t" << x.array().exp().transpose() << std::endl;
- return x;
-
- }
-
-
-
-
- template < typename cgScalar >
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &b, const cgScalar &minVal,
- const cgScalar &maxVal) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
- MatrixXr ATA = A.transpose()*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 100;
- cgScalar SIGMA = 1e-2;
- cgScalar delta = 1e-4;
-
-
-
-
-
-
-
- cgScalar limsup = 1e10;
- for (int i=0; i<N; ++i) {
- VectorXr Spike = VectorXr::Zero(N);
- Spike(i) = (minVal + maxVal) / 2.;
- limsup = std::min(limsup, (ATA*Spike).array().abs().maxCoeff());
- }
- cgScalar lambda = 1e3*limsup;
- cgScalar mux0 = 1e-1*lambda;
- cgScalar muy0 = 1e-1*lambda;
- cgScalar epsilon = 1e-20;
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- cgScalar phid = (b - (b_pre)).norm();
- cgScalar phim = x.norm();
- cgScalar phib = PhiB(mux0, muy0, minVal, maxVal, x);
- cgScalar mux = mux0;
- cgScalar muy = muy0;
- cgScalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
- std::cout << "Starting CG iterations 7" << std::endl;
-
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
-
-
-
-
- VectorXr WmT_Wm = VectorXr::Ones(N).array()*lambda;
-
- VectorXr Xm1 = x;
- int iter = N*M;
- mux = mux0;
- muy = muy0;
- int iloop(0);
-
- do {
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
-
- MatrixXr A2;
-
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- Y1 = VectorXr::Ones(N).array() / (maxVal-x.array()) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- Y2 = VectorXr::Ones(N).array() / ((maxVal-x.array())*(maxVal-x.array()));
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- b2 = -(A.transpose()*(b_pre-b)).array() -
- (WmT_Wm.array()*x.array()).array() +
- (2.*mux)*X1.array() + (2.*muy)*Y1.array();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- A2 = ATA;
- A2.diagonal().array() += WmT_Wm.array() + mux*X2.array() +
- muy*Y2.array();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
-
-
-
-
-
-
-
-
-
-
-
- tol = 1e-5*phid+mux+muy;
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr ztilde = CG(A2, VectorXr::Zero(N), b2, iter, tol);
-
-
-
-
-
- VectorXr s1, s2;
-
-
-
-
-
- cgScalar d = std::min(1.,0.995*(x.array()/ztilde.array().abs()).minCoeff());
- x += d*ztilde;
-
-
-
-
-
-
-
-
-
-
-
-
- b_pre = A*x;
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phib = PhiB(mux, muy, minVal, maxVal, x);
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phid = (b - (b_pre)).norm();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phim = x.norm();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s1 = mux * (X2.asDiagonal()*(x+ztilde) - 2.*X1);
- mux = SIGMA/((cgScalar)(N)) * std::abs( s1.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s2 = muy * (Y2.asDiagonal()*(x+ztilde) - 2.*Y1);
- muy = SIGMA/((cgScalar)(N)) * std::abs( s2.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
- ++iloop;
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
-
- phireport.precision(12);
- std::cout << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
- phireport << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
-
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .92;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
- template <typename cgScalar>
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &b, const cgScalar &minVal,
- const cgScalar &maxVal, const int& block) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
-
-
-
-
-
-
- MatrixXr ATA = A.transpose()*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 100;
- cgScalar SIGMA = 1e-2;
- cgScalar delta = 1e-4;
-
-
-
-
-
-
-
- cgScalar limsup = 1e10;
- for (int i=0; i<N; ++i) {
- VectorXr Spike = VectorXr::Zero(N);
- Spike(i) = (minVal + maxVal) / 2.;
- limsup = std::min(limsup, (ATA*Spike).array().abs().maxCoeff());
- }
- cgScalar lambda = 1e-6;
- cgScalar mux0 = 1e-1*lambda;
- cgScalar muy0 = 1e-1*lambda;
- cgScalar epsilon = 1e-25;
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- int nblocks = x.size()/block;
- cgScalar phid = (b - (b_pre)).norm();
- cgScalar phim = x.norm();
- cgScalar phib = PhiB2(mux0, muy0, minVal, maxVal, x, block, nblocks);
- cgScalar mux = mux0;
- cgScalar muy = muy0;
- cgScalar tol = 1e-5*phid;
- std::fstream phireport("phimphid.dat", std::ios::out);
- std::cout << "Starting CG iterations 6 " << std::endl;
-
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
-
-
-
-
- VectorXr WmT_Wm = VectorXr::Ones(N).array()*lambda;
-
- VectorXr Xm1 = x;
- int iter = N*M;
- mux = mux0;
- muy = muy0;
- int iloop(0);
-
-
-
-
-
-
-
-
- do {
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
-
- MatrixXr A2;
-
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- for (int ib=0; ib<nblocks; ++ib) {
-
- Y1.segment(ib*block, block) = VectorXr::Ones(block).array() /
- (maxVal - x.segment(ib*block, block).sum());
- }
-
-
-
-
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- for (int ib=0; ib<nblocks; ++ib) {
-
- Y2.segment(ib*block, block) = VectorXr::Ones(block).array() /
- ( (maxVal-x.segment(ib*block, block).sum()) *
- (maxVal-x.segment(ib*block, block).sum()) );
- }
-
-
-
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- b2 = -(A.transpose()*(b_pre-b)).array() -
- (WmT_Wm.array()*x.array()).array() +
- (2.*mux)*X1.array() + (2.*muy)*Y1.array();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- A2 = ATA;
- A2.diagonal().array() += WmT_Wm.array() + mux*X2.array() +
- muy*Y2.array();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
-
-
-
-
-
-
-
-
-
-
-
- tol = 1e-5*phid+mux+muy;
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
- VectorXr ztilde = CG(A2, VectorXr::Zero(N), b2, iter, tol);
-
-
-
-
-
- VectorXr s1, s2;
-
-
-
-
-
- cgScalar d = std::min(1.,0.995*(x.array()/ztilde.array().abs()).minCoeff());
- x += d*ztilde;
-
-
-
-
-
-
-
-
-
-
-
-
- b_pre = A*x;
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phib = PhiB2(mux, muy, minVal, maxVal, x, block, nblocks);
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phid = (b - (b_pre)).norm();
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- phim = x.norm();
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
-
-
- #ifdef LEMMAUSEOMP
- #pragma omp parallel sections
- {
- #endif
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s1 = mux * (X2.asDiagonal()*(x+ztilde) - 2.*X1);
- mux = SIGMA/((cgScalar)(N)) * std::abs( s1.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- #pragma omp section
- #endif
- {
- s2 = muy * (Y2.asDiagonal()*(x+ztilde) - 2.*Y1);
- muy = SIGMA/((cgScalar)(N)) * std::abs( s2.dot(x) ) ;
- }
- #ifdef LEMMAUSEOMP
- }
- #endif
- ++iloop;
- } while (std::abs(phib / (phid+lambda*phim)) > epsilon);
-
-
- phireport.precision(12);
- std::cout << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
- phireport << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << iter << "\t" << iloop << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << x.norm() << "\t" << (b-(A*x)).norm()
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .85;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
- template <typename cgScalar>
- VectorXr LogBarrierCG(const MatrixXr &A, const VectorXr &xr,
- const VectorXr &b,
- const cgScalar &minVal,
- const cgScalar &maxVal, const int& block,
- const Eigen::SparseMatrix<cgScalar>& WdTWd,
- const Eigen::SparseMatrix<cgScalar>& WmTWm, Real lambda0=1e1) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- std::cerr << "A.rows() " << A.rows() << "\n";
- std::cerr << "A.cols() " << A.cols() << "\n";
- std::cerr << "b.size() " << b.size() << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
- int N = A.cols();
- int M = A.rows();
- int MAXITER = 175;
- cgScalar SIGMA = .25;
- cgScalar delta = 1e-4;
-
-
-
-
-
-
-
-
-
- cgScalar lambda = lambda0;
- cgScalar epsilon = 1e-15;
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
- int nblocks = x.size()/block;
- cgScalar phid = (WdTWd*(b - b_pre)).norm();
- cgScalar phim = (WmTWm*(x - xr)).norm();
- cgScalar phib = PhiB2(minVal, maxVal, x, block, nblocks);
- cgScalar mux = (phid + lambda*phim) / phib;
- cgScalar muy = mux;
-
- std::fstream phireport("phimphid.dat", std::ios::out);
- std::cout << "Starting CG iterations 5" << std::endl;
-
- Eigen::ConjugateGradient< MatrixXr, Eigen::Upper > cg;
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
- int iter = N*M;
- mux = (phid + lambda*phim) / phib;
- muy = mux;
- int iloop(0);
- int itertot(0);
- VectorXr h;
- bool cont(true);
- do {
-
-
-
- VectorXr X1(x.size());
- VectorXr Y1(x.size());
- VectorXr X2(x.size());
- VectorXr Y2(x.size());
- VectorXr b2;
- MatrixXr A2;
-
-
-
-
-
-
-
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
-
- for (int ib=0; ib<nblocks; ++ib) {
- Y1.segment(ib*block, block) = VectorXr::Ones(block).array() /
- (maxVal - x.segment(ib*block, block).sum());
- }
-
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
-
- for (int ib=0; ib<nblocks; ++ib) {
- Y2.segment(ib*block, block) = VectorXr::Ones(block).array() /
- ( (maxVal-x.segment(ib*block, block).sum()) *
- (maxVal-x.segment(ib*block, block).sum()) );
- }
-
-
-
-
-
-
-
- b2 = (A.transpose()*WdTWd*(b)).array()
-
- + (2.*mux)*X1.array() + (2.*muy)*Y1.array();
-
- A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
- A2.diagonal().array() += mux*X2.array() + muy*Y2.array();
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- iter = N*M;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- cg.compute(A2);
- VectorXr ztilde;
- ztilde = cg.solveWithGuess(b2, x);
- iter = cg.iterations();
-
-
- ++iloop;
- itertot += iter;
-
-
-
-
-
-
- h = ztilde - x;
-
-
-
- cgScalar d(1.);
- for (int ix=0; ix<x.size(); ++ix) {
- if (h[ix] < 0.) {
- d = std::min(d, (cgScalar).925*(x[ix]/std::abs(h[ix])));
- }
- }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- x += d*h;
-
-
- for (int ib=0; ib<nblocks; ++ib) {
- while (x.segment(ib*block, block).sum() >= maxVal) {
- x.segment(ib*block, block).array() *= .99;
- }
- }
-
- for (int i=0; i<x.size(); ++i) {
- if (x(i) < minVal) {
- x(i) = minVal + delta;
- }
- }
-
- b_pre = A*x;
-
- phib = PhiB2(mux, muy, minVal, maxVal, x, block, nblocks);
- phid = std::sqrt((WdTWd*(b-b_pre)).norm());
- phim = (WmTWm*(x-xr)).norm();
-
-
- VectorXr s1 = mux * (X2.asDiagonal()*(ztilde) - 2.*X1);
- mux = SIGMA/((cgScalar)(N)) * std::abs( s1.dot(x) ) ;
- VectorXr s2 = muy * (Y2.asDiagonal()*(ztilde) - 2.*Y1);
- muy = SIGMA/((cgScalar)(N)) * std::abs( s2.dot(x) ) ;
-
- if ( (std::abs(phib / (phid+lambda*phim)) < epsilon) && phid < 1000. ) {
-
- cont = false;
- }
-
- } while ( cont );
-
-
-
-
- std::cout << std::setw(6) << ii << std::scientific << std::setw(18) << phim << std::setw(18) << phid
- << std::setw(18) << lambda << std::setw(18) << mux << std::setw(18) << muy
- << std::setw(12) << itertot << std::setw(12) << iloop << std::setw(18) << h.norm() << std::endl;
-
- phireport.precision(12);
- phireport << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t"
- << itertot << "\t" << iloop << "\t" << h.norm() << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t" << muy << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .9;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
-
- template <typename cgScalar>
- VectorXr LogBarrierCG_NN(const MatrixXr &A, const VectorXr &xr,
- const VectorXr &b, const cgScalar &minVal,
- const Eigen::SparseMatrix<cgScalar>& WdTWd,
- const Eigen::SparseMatrix<cgScalar>& WmTWm, Real lambda0=1e1, int MAXITER=175) {
-
-
- if (A.rows() != b.size() ) {
- std::cerr << "Matrix A is not aligned with Vector b" << "\n";
- std::cerr << "A.rows() " << A.rows() << "\n";
- std::cerr << "A.cols() " << A.cols() << "\n";
- std::cerr << "b.size() " << b.size() << "\n";
- exit(1);
- }
-
-
- std::fstream obsdata;
- std::string fname = "obsdata.dat";
- obsdata.open(fname.c_str(), std::ios::out);
- obsdata << b << std::endl;
- obsdata.close();
-
-
- MatrixXr ATWdTWdA = A.transpose()*WdTWd*A;
- int N = A.cols();
- int M = A.rows();
-
- cgScalar SIGMA = .25;
- cgScalar delta = 1e-4;
-
-
-
-
-
-
-
-
-
- cgScalar lambda = lambda0;
- cgScalar epsilon = 1e-16;
-
-
- VectorXr x = VectorXr::Zero(N).array() + minVal+delta;
-
-
- VectorXr b_pre = A*x;
-
-
-
- Eigen::ConjugateGradient< MatrixXr, Eigen::Upper, Eigen::DiagonalPreconditioner<Real> > cg;
-
-
- cgScalar phid = (WdTWd*(b - b_pre)).norm();
- cgScalar phim = (WmTWm*(x - xr)).norm();
- cgScalar phib = PhiB2_NN(1., minVal, x);
- cgScalar mux = (phid + lambda*phim) / phib;
-
- std::fstream phireport("phimphid.dat", std::ios::out);
- std::cout << "Starting CG iterations 4" << std::endl;
-
-
- for (int ii=0; ii<MAXITER; ++ii) {
-
- int iter = N*M;
- mux = (phid + lambda*phim) / phib;
- int iloop(0);
- int itertot(0);
- VectorXr h;
- bool cont(true);
- do {
-
-
-
- VectorXr X1(x.size());
- VectorXr X2(x.size());
- VectorXr b2;
- MatrixXr A2;
-
-
-
-
-
-
-
- X1 = VectorXr::Ones(N).array() / (x.array()-minVal) ;
- X2 = VectorXr::Ones(N).array() / ((x.array()-minVal)*(x.array()-minVal));
-
-
- b2 = (A.transpose()*WdTWd*(b)).array()
-
- + (2.*mux)*X1.array();
-
- A2 = ATWdTWdA;
- A2 += lambda*WmTWm;
- A2.diagonal().array() += mux*X2.array();
-
-
-
-
-
-
- iter = N*M;
-
-
-
- cg.compute(A2);
- VectorXr ztilde;
- ztilde = cg.solveWithGuess(b2, x);
-
- iter = cg.iterations();
-
-
- ++iloop;
- itertot += iter;
-
-
-
- h = ztilde - x;
-
-
- cgScalar d(1.);
- for (int ix=0; ix<x.size(); ++ix) {
- if (h[ix] < 0.) {
- d = std::min(d, (cgScalar).925*(x[ix]/std::abs(h[ix])));
- }
- }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- x += d*h;
-
- for (int i=0; i<x.size(); ++i) {
- if (x(i) < minVal) {
- x(i) = minVal + delta;
- }
- }
-
- b_pre = A*x;
-
- phib = PhiB2_NN(mux, minVal, x);
- phid = (WdTWd*(b-b_pre)).norm();
- phim = (WmTWm*(x-xr)).norm();
-
-
- VectorXr s1 = mux * (X2.asDiagonal()*(ztilde) - 2.*X1);
- mux = SIGMA/((cgScalar)(N)) * std::abs( s1.dot(x) ) ;
-
- if ( (std::abs(phib / (phid+lambda*phim)) < epsilon)) {
-
- cont = false;
- }
-
- } while ( cont );
-
-
-
-
- std::cout << std::setw(6) << ii << std::scientific << std::setw(18) << phim << std::setw(18) << phid
- << std::setw(18) << lambda << std::setw(18) << mux
- << std::setw(12) << itertot << std::setw(12) << iloop << std::setw(18) << h.norm() << std::endl;
-
- phireport.precision(12);
- phireport << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t"
- << itertot << "\t" << iloop << "\t" << h.norm() << std::endl;
-
- std::fstream modfile;
- std::string fname = "iteration" + to_string(ii) + ".dat";
- modfile.open( fname.c_str(), std::ios::out);
- modfile << ii << "\t" << phim << "\t" << phid
- << "\t" << lambda << "\t" << mux << "\t" << iter << "\n";
- modfile << x << "\n";
- modfile.close();
-
-
- std::fstream predata;
- fname = "iteration" + to_string(ii) + "pre.dat";
- predata.open(fname.c_str(), std::ios::out);
- predata << b_pre << std::endl;
- predata.close();
-
-
-
- lambda *= .9;
-
- }
- phireport.close();
-
- return x;
- }
-
-
-
-
-
-
- }
-
- #endif
|