mirror of
https://asciireactor.com/otho/psdlag-agn.git
synced 2025-10-31 22:38:03 +00:00
236 lines
8.2 KiB
C++
236 lines
8.2 KiB
C++
|
/*
|
||
|
|
* mod.cpp
|
||
|
|
*
|
||
|
|
* Created on: May 31, 2013
|
||
|
|
* Author: azoghbi
|
||
|
|
*/
|
||
|
|
|
||
|
|
#include "inc/mod.hpp"
|
||
|
|
|
||
|
|
mod::mod() {
|
||
|
|
n = 0 ; nfq = 0 ; npar = 0; dt = 0; info = 0;f1=0;f2=0;
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
mod::~mod() {}
|
||
|
|
|
||
|
|
void mod::init( vec fqL , int fac ){
|
||
|
|
|
||
|
|
// --------------- Freq stuff----------------- //
|
||
|
|
int i,j;
|
||
|
|
nfq = fqL.length() - 1;
|
||
|
|
npar = nfq * fac ;
|
||
|
|
FqL.setlength(nfq+1); for( i=0 ; i<=nfq ; i++ ){ FqL[i] = fqL[i]; }
|
||
|
|
// ------------------------------------------ //
|
||
|
|
|
||
|
|
|
||
|
|
// ---------------- Integrals ---------------- //
|
||
|
|
Cfq.resize( nfq );Sfq.resize( nfq );
|
||
|
|
for( i=0 ; i<nfq ; i++ ){ Cfq[i].setlength( n , n ); Sfq[i].setlength( n , n );}
|
||
|
|
_CSfq( fqL , Cfq , Sfq );
|
||
|
|
|
||
|
|
Cfq2.resize( 3 );Sfq2.resize( 3 );
|
||
|
|
for( i=0 ; i<3 ; i++ ){ Cfq2[i].setlength( n , n ); Sfq2[i].setlength( n , n );}
|
||
|
|
vec fqL2;fqL2.setlength(4); fqL2[0] = 1e-1*fqL[0]; fqL2[1] = fqL[0]; fqL2[2] = fqL[nfq]; fqL2[3] = 10*fqL[nfq];
|
||
|
|
f1 = 2; f2 = 0.5;
|
||
|
|
_CSfq( fqL2 , Cfq2 , Sfq2 );
|
||
|
|
// ------------------------------------------ //
|
||
|
|
|
||
|
|
|
||
|
|
// -------------- Containers ---------------- //
|
||
|
|
C.setlength( n , n );
|
||
|
|
Ci.setlength( n , n );
|
||
|
|
I.setlength( n , n );
|
||
|
|
for( i=0 ; i<n ; i++ ){ I(i,i) = 1.0; for( j=0 ; j<i ; j++ ){ I(i,j)=0; I(j,i)=0;}}
|
||
|
|
Cilc.setlength( n );
|
||
|
|
yyT.setlength( n , n );
|
||
|
|
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ yyT(i,j) = lc[i]*lc[j];}}
|
||
|
|
yyTmC.setlength( n , n );
|
||
|
|
dC.setlength( n , n );
|
||
|
|
Cv.resize( npar ); for( i=0 ; i<npar ; i++ ){ Cv[i].setlength( n , n ); }
|
||
|
|
C2.setlength( n , n );
|
||
|
|
C3.setlength( n , n );
|
||
|
|
|
||
|
|
// ------------------------------------------ //
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
void mod::_CSfq( vec fqL , vector<vec2>&cfq , vector<vec2>& sfq ){
|
||
|
|
//fqL[0] *= 1e-6; fqL[nfq]*=2;
|
||
|
|
int i,j,k, nfq=fqL.length()-1;
|
||
|
|
double tt,norm,t1,t2,pi=M_PI,pidt=pi*dt,sdt2,sdt1,cdt2,cdt1,s2dt2,s2dt1,c2dt2,c2dt1,dtmt,dtpt;
|
||
|
|
double st2,st1,ct2,ct1, sm2,sm1,cm2,cm1, sp2,sp1,cp2,cp1;
|
||
|
|
vec w; w.setlength(nfq+1); for( k=0 ; k<=nfq ; k++ ){ w[k] = 2*pi*fqL[k]; }
|
||
|
|
norm = 1./( pidt*pidt );
|
||
|
|
|
||
|
|
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<n ; j++ ){
|
||
|
|
tt = t[j] - t[i];
|
||
|
|
dtmt = dt - tt; dtpt = dt+tt;
|
||
|
|
for( k=0 ; k<nfq ; k++ ){
|
||
|
|
if( tt==0 ){
|
||
|
|
sinecosineintegrals( w[k+1]*dt , sdt2 , cdt2 );
|
||
|
|
sinecosineintegrals( w[k] *dt , sdt1 , cdt1 );
|
||
|
|
t1 = ( cos( w[k+1]*dt ) - 1 ) / ( 2 * fqL[k+1] ) + pidt * sdt2;
|
||
|
|
t1 -= ( cos( w[k] *dt ) - 1 ) / ( 2 * fqL[k] ) + pidt * sdt1;
|
||
|
|
cfq[k](i,j) = norm * t1;
|
||
|
|
sfq[k](i,j) = 0;
|
||
|
|
}else if( fabs(tt)==dt ){
|
||
|
|
sinecosineintegrals( w[k+1]*dt , sdt2 , cdt2 );
|
||
|
|
sinecosineintegrals( w[k] *dt , sdt1 , cdt1 );
|
||
|
|
sinecosineintegrals( 2*w[k+1]*dt , s2dt2 , c2dt2 );
|
||
|
|
sinecosineintegrals( 2*w[k] *dt , s2dt1 , c2dt1 );
|
||
|
|
t1 = -( cos( w[k+1]*dt )*pow( sin( 0.5*w[k+1]*dt ) , 2 ) ) / fqL[k+1];
|
||
|
|
t1 -= -( cos( w[k] *dt )*pow( sin( 0.5*w[k] *dt ) , 2 ) ) / fqL[k] ;
|
||
|
|
t2 = pidt * ( s2dt2 - sdt2 );
|
||
|
|
t2 -= pidt * ( s2dt1 - sdt1 );
|
||
|
|
cfq[k](i,j) = norm * ( t1 + t2 );
|
||
|
|
|
||
|
|
t1 = -( 2*cos( 0.5*w[k+1]*dt )*pow( sin( 0.5*w[k+1]*dt ) , 3 ) ) / fqL[k+1];
|
||
|
|
t1 -= -( 2*cos( 0.5*w[k] *dt )*pow( sin( 0.5*w[k] *dt ) , 3 ) ) / fqL[k] ;
|
||
|
|
t2 = pidt * ( cdt2 - c2dt2 );
|
||
|
|
t2 -= pidt * ( cdt1 - c2dt1 );
|
||
|
|
sfq[k](i,j) = norm * fabs(tt)/tt * ( t1 + t2 );
|
||
|
|
}else{
|
||
|
|
sinecosineintegrals( w[k+1]*tt , st2 , ct2 );
|
||
|
|
sinecosineintegrals( w[k] *tt , st1 , ct1 );
|
||
|
|
sinecosineintegrals( w[k+1]*dtmt , sm2 , cm2 );
|
||
|
|
sinecosineintegrals( w[k] *dtmt , sm1 , cm1 );
|
||
|
|
sinecosineintegrals( w[k+1]*dtpt , sp2 , cp2 );
|
||
|
|
sinecosineintegrals( w[k] *dtpt , sp1 , cp1 );
|
||
|
|
t1 = ( cos( w[k+1]*dtmt ) - 2*cos( w[k+1]*tt ) + cos( w[k+1]*dtpt ) )/(4*fqL[k+1]);
|
||
|
|
t1 -= ( cos( w[k] *dtmt ) - 2*cos( w[k] *tt ) + cos( w[k] *dtpt ) )/(4*fqL[k] );
|
||
|
|
t2 = ( dtmt*sm2 - 2*tt*st2 + dtpt*sp2 ) * (pi/2);
|
||
|
|
t2 -= ( dtmt*sm1 - 2*tt*st1 + dtpt*sp1 ) * (pi/2);
|
||
|
|
cfq[k](i,j) = norm * ( t1 + t2 );
|
||
|
|
|
||
|
|
t1 = ( -sin( w[k+1]*dtmt ) - 2*sin( w[k+1]*tt ) + sin( w[k+1]*dtpt ) )/(4*fqL[k+1]);
|
||
|
|
t1 -= ( -sin( w[k] *dtmt ) - 2*sin( w[k] *tt ) + sin( w[k] *dtpt ) )/(4*fqL[k] );
|
||
|
|
t2 = ( dtmt*cm2 + 2*tt*ct2 - dtpt*cp2 ) * (pi/2);
|
||
|
|
t2 -= ( dtmt*cm1 + 2*tt*ct1 - dtpt*cp1 ) * (pi/2);
|
||
|
|
sfq[k](i,j) = norm * ( t1 + t2 );
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
}}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void mod::_C( vec ){
|
||
|
|
cerr << endl << "mod::_C is an empty function. It should not be called directly" << endl << endl;;
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
void mod::_dC( vec , int ){
|
||
|
|
cerr << endl << "mod::_C is an empty function. It should not be called directly" << endl << endl;;
|
||
|
|
exit(1);
|
||
|
|
}
|
||
|
|
|
||
|
|
double mod::loglikelihood( vec p ){
|
||
|
|
double logdet , chi2 , loglike ; int i;
|
||
|
|
_C( p );
|
||
|
|
if (not spdmatrixcholesky( C , n , false )) {return -1e20;};
|
||
|
|
logdet = 0 ; for( i=0 ; i<n ; i++ ){ logdet += log( C(i,i) ); } logdet *= 2;
|
||
|
|
spdmatrixcholeskysolvem( C , n , false , I , n , info , rep , Ci );
|
||
|
|
rmatrixmv( n , n , Ci , 0 , 0 , 0 , lc , 0 , Cilc , 0 );
|
||
|
|
chi2 = 0 ; for( i=0 ; i<n ; i++ ){ chi2 += Cilc[i]*lc[i] ; }
|
||
|
|
loglike = -0.5 * ( chi2 + logdet + n*log(2*M_PI) );
|
||
|
|
return loglike;
|
||
|
|
}
|
||
|
|
|
||
|
|
void mod::dlikelihood( vec p , double& loglike , vec& grad , vec2& hess ){
|
||
|
|
double logdet , chi2 , res; int i,j,k;
|
||
|
|
_C( p );
|
||
|
|
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ yyTmC(i,j) = yyT(i,j) - C(i,j); yyTmC(j,i) = yyTmC(i,j);}}
|
||
|
|
spdmatrixcholesky( C , n , false );
|
||
|
|
logdet = 0 ; for( i=0 ; i<n ; i++ ){ logdet += log( C(i,i) ); } logdet *= 2;
|
||
|
|
spdmatrixcholeskysolvem( C , n , false , I , n , info , rep , Ci );
|
||
|
|
rmatrixmv( n , n , Ci , 0 , 0 , 0 , lc , 0 , Cilc , 0 );
|
||
|
|
chi2 = 0 ; for( i=0 ; i<n ; i++ ){ chi2 += Cilc[i]*lc[i] ; }
|
||
|
|
loglike = -0.5 * ( chi2 + logdet + n*log(2*M_PI) );
|
||
|
|
|
||
|
|
|
||
|
|
for( i=0 ; i<npar ; i++ ){
|
||
|
|
_dC( p , i );
|
||
|
|
rmatrixgemm( n , n , n , 1.0 , Ci , 0 , 0 , 0 , dC , 0 , 0 , 0 , 0.0 , Cv[i] , 0 , 0 );
|
||
|
|
}
|
||
|
|
|
||
|
|
for( i=0 ; i<npar ; i++ ){
|
||
|
|
rmatrixgemm( n , n , n , 1.0 , Cv[i] , 0 , 0 , 0 , Ci , 0 , 0 , 0 , 0.0 , C2 , 0 , 0 );
|
||
|
|
rmatrixgemm( n , n , n , 1.0 , yyTmC , 0 , 0 , 0 , C2 , 0 , 0 , 0 , 0.0 , C3 , 0 , 0 );
|
||
|
|
res = 0; for( k=0 ; k<n ; k++ ){ res += C3(k,k); }
|
||
|
|
grad[i] = 0.5 * res;
|
||
|
|
for( j=0 ; j<=i ; j++ ){
|
||
|
|
rmatrixgemm( n , n , n , 1.0 , Cv[i] , 0 , 0 , 0 , Cv[j] , 0 , 0 , 0 , 0.0 , C3 , 0 , 0 );
|
||
|
|
res = 0; for( k=0 ; k<n ; k++ ){ res += C3(k,k); }
|
||
|
|
hess(i,j) = 0.5 * res; hess(j,i) = 0.5*res;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void mod::optimize( vec& pars , vec& errs ){
|
||
|
|
int nmax = 200;
|
||
|
|
double tol = 1e-3;
|
||
|
|
|
||
|
|
int i,j,n;
|
||
|
|
double loglike,dpmax;
|
||
|
|
vec tmppars,dpar,grad; vec2 hess,hessi,ii;
|
||
|
|
tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar);
|
||
|
|
hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
|
||
|
|
for( i=0 ; i<npar ; i++ ){
|
||
|
|
tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
|
||
|
|
}
|
||
|
|
|
||
|
|
for( n=1 ; n<= nmax ; n++ ){
|
||
|
|
|
||
|
|
dlikelihood( tmppars , loglike , grad , hess );
|
||
|
|
spdmatrixcholesky( hess , npar , false );
|
||
|
|
spdmatrixcholeskysolvem( hess , npar , false , ii , npar , info , rep , hessi );
|
||
|
|
|
||
|
|
dpmax = -1e20;
|
||
|
|
for( i=0 ; i<npar ; i++ ){
|
||
|
|
dpar[i] = 0; for( j=0 ; j<npar ; j++ ){ dpar[i] += grad[j] * hessi(i,j);}
|
||
|
|
dpmax = max( dpmax , fabs(dpar[i]) );
|
||
|
|
}
|
||
|
|
step_pars( n , dpar , tmppars );
|
||
|
|
|
||
|
|
cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
|
||
|
|
if ( dpmax<tol ) break;
|
||
|
|
}
|
||
|
|
for( i=0 ; i<npar ; i++ ){
|
||
|
|
pars[i] = tmppars[i];
|
||
|
|
errs[i] = sqrt(hessi(i,i));
|
||
|
|
}
|
||
|
|
cout << setfill('-') << setw(40) << "\n";
|
||
|
|
for(i=0;i<npar;i++){cout << pars[i] << " ";} cout << endl;
|
||
|
|
for(i=0;i<npar;i++){cout << errs[i] << " ";} cout << endl;
|
||
|
|
for(i=0;i<npar;i++){cout << grad[i] << " ";} cout << endl;
|
||
|
|
cout << setfill('-') << setw(40) << "\n" << setfill(' ');
|
||
|
|
print_pars( pars, errs );
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
void mod::step_pars( int n , vec& dpar , vec& pars ){
|
||
|
|
for( int i=0 ; i<npar ; i++ ){
|
||
|
|
pars[i] += dpar[i];
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
void mod::print_pars( vec& pars, vec& errs ){
|
||
|
|
cout << setfill('*') << setw(40) << "\n" << setfill(' ');
|
||
|
|
cout << endl << "# "; for(int i=0;i<=nfq;i++){cout << FqL[i] << " ";} cout << endl;
|
||
|
|
for(int i=0;i<npar;i++){ printf( "%3.3e %3.3e\n" , pars[i] , errs[i] );}
|
||
|
|
cout << setfill('*') << setw(40) << "\n" << setfill(' ');
|
||
|
|
}
|
||
|
|
|
||
|
|
void mod::print_pars( vec& pars, vec& errs1 , vec& errs2 ){
|
||
|
|
cout << setfill('*') << setw(40) << "\n" << setfill(' ');
|
||
|
|
cout << endl << "# "; for(int i=0;i<=nfq;i++){cout << FqL[i] << " ";} cout << endl;
|
||
|
|
for(int i=0;i<npar;i++){ printf( "%3.3e %3.3e %3.3e\n" , pars[i] , errs1[i] , errs2[i] );}
|
||
|
|
cout << setfill('*') << setw(40) << "\n" << setfill(' ');
|
||
|
|
}
|
||
|
|
|
||
|
|
void mod::what_pars( int& ip1 , int& ip2 ){
|
||
|
|
ip1 = 0; ip2 = npar;
|
||
|
|
}
|