mirror of
https://asciireactor.com/otho/psdlag-agn.git
synced 2024-11-21 20:05:06 +00:00
removed references to psdlag src
This commit is contained in:
caes
parent
bf88001dc8
commit
1ba53f4334
Binary file not shown.
Binary file not shown.
@ -1,36 +0,0 @@
|
||||
/*
|
||||
* def.hpp
|
||||
*
|
||||
* Created on: May 24, 2013
|
||||
* Author: abduz
|
||||
*/
|
||||
|
||||
#ifndef DEF_HPP_
|
||||
#define DEF_HPP_
|
||||
|
||||
#include <alglib/ap.h>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
|
||||
using namespace alglib;
|
||||
using namespace std;
|
||||
|
||||
typedef real_1d_array vec;
|
||||
typedef real_2d_array vec2;
|
||||
typedef integer_1d_array ivec;
|
||||
|
||||
|
||||
class lcurve {
|
||||
public:
|
||||
int len; double dt; vec lc,lce,t;
|
||||
lcurve(int n=1){t.setlength(n); lc.setlength(n); lce.setlength(n) ;len=n;dt=1.0;}
|
||||
lcurve(vec& t_,vec& lc_,vec& lce_,double dt_){t=t_;lc=lc_;lce=lce_;len=t.length();dt=dt_;}
|
||||
double demean(){double m=0;int i;for(i=0;i<len;i++){m+=lc[i];} m/=len;for(i=0;i<len;i++){lc[i]-=m;}; return m;}
|
||||
friend ostream& operator<<(ostream&o,lcurve&l){ o << setprecision(12);
|
||||
for(int i=0;i<l.len;i++){o << l.t[i] << " " << l.lc[i] << " " << l.lce[i] << endl;}
|
||||
return(o);}
|
||||
};
|
||||
|
||||
|
||||
|
||||
#endif /* DEF_HPP_ */
|
||||
@ -1,56 +0,0 @@
|
||||
/*
|
||||
* lag.hpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#ifndef LAG_HPP_
|
||||
#define LAG_HPP_
|
||||
|
||||
#include "mod.hpp"
|
||||
|
||||
class lag : public mod {
|
||||
|
||||
int n1;
|
||||
vec p1,p2;
|
||||
ivec icx,iphi;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
|
||||
public:
|
||||
lag( lcurve , lcurve , vec , vec );
|
||||
virtual ~lag();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
};
|
||||
|
||||
// ----------------------------------- //
|
||||
|
||||
class lag10 : public mod {
|
||||
|
||||
int n1;
|
||||
vec p1,p2;
|
||||
ivec icx,iphi;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
|
||||
public:
|
||||
lag10( lcurve , lcurve , vec , vec );
|
||||
virtual ~lag10();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
void what_pars( int& , int& );
|
||||
};
|
||||
|
||||
#endif /* LAG_HPP_ */
|
||||
@ -1,58 +0,0 @@
|
||||
/*
|
||||
* lag.hpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#ifndef LAG_RMS_HPP_
|
||||
#define LAG_RMS_HPP_
|
||||
|
||||
#include "mod.hpp"
|
||||
|
||||
class lag_rms : public mod {
|
||||
|
||||
int n1;
|
||||
vec p1,p2;
|
||||
ivec icx,iphi;
|
||||
double mean1,mean2,mean1sq,mean2sq,mean12;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
|
||||
public:
|
||||
lag_rms( lcurve , lcurve , vec , vec );
|
||||
virtual ~lag_rms();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
};
|
||||
|
||||
// ----------------------------------- //
|
||||
|
||||
class lag10_rms : public mod {
|
||||
|
||||
int n1;
|
||||
vec p1,p2;
|
||||
ivec icx,iphi;
|
||||
double mean1,mean2,mean1sq,mean2sq,mean12;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
|
||||
public:
|
||||
lag10_rms( lcurve , lcurve , vec , vec );
|
||||
virtual ~lag10_rms();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
void what_pars( int& , int& );
|
||||
};
|
||||
|
||||
#endif /* LAG_HPP_ */
|
||||
@ -1,58 +0,0 @@
|
||||
/*
|
||||
* main.hpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#ifndef MAIN_HPP_
|
||||
#define MAIN_HPP_
|
||||
|
||||
|
||||
#include <fstream>
|
||||
#include <cstring>
|
||||
#include <sstream>
|
||||
#include <vector>
|
||||
|
||||
#include "def.hpp"
|
||||
#include "psd.hpp"
|
||||
#include "psd_rms.hpp"
|
||||
#include "lag.hpp"
|
||||
#include "lag_rms.hpp"
|
||||
#include "psdlag.hpp"
|
||||
#include "psdlag_rms.hpp"
|
||||
#include "mcmc.hpp"
|
||||
|
||||
|
||||
void usage();
|
||||
void do_work( char* );
|
||||
void readLC( vector<vector<lcurve> >&, string , int , int , int , bool );
|
||||
|
||||
|
||||
double mcmc_lag10( vec& x ,void*ptr ){
|
||||
int i,hnp,np = x.length();hnp=np/2; double phi,twopi=2*M_PI;
|
||||
for(i=0;i<hnp;i++){
|
||||
if(x[i]>7){x[i]=7;}if(x[i]<-7){x[i]=-7;}
|
||||
phi = x[i+hnp]+M_PI;
|
||||
while(phi>twopi){phi -= twopi;}
|
||||
while(phi<0){phi += twopi;}
|
||||
x[i+hnp] = phi-M_PI;
|
||||
}
|
||||
Mod<lag10_rms> *mod = (Mod<lag10_rms>*) ptr; double logl=mod->loglikelihood(x);
|
||||
return logl;
|
||||
}
|
||||
|
||||
double mcmc_psdlag10( vec& x ,void*ptr ){
|
||||
int i,hnp,np = x.length();hnp=np/2; double phi,twopi=2*M_PI;
|
||||
for(i=0;i<hnp;i++){
|
||||
if(x[i]>7){x[i]=7;}if(x[i]<-7){x[i]=-7;}
|
||||
phi = x[i+hnp]+M_PI;
|
||||
while(phi>twopi){phi -= twopi;}
|
||||
while(phi<0){phi += twopi;}
|
||||
x[i+hnp] = phi-M_PI;
|
||||
}
|
||||
Mod<psdlag10_rms> *mod = (Mod<psdlag10_rms>*) ptr; double logl=mod->loglikelihood(x);
|
||||
return logl;
|
||||
}
|
||||
|
||||
#endif /* MAIN_HPP_ */
|
||||
@ -1,35 +0,0 @@
|
||||
/*
|
||||
* mcmc.hpp
|
||||
*
|
||||
* Created on: May 14, 2013
|
||||
* Author: abduz
|
||||
*/
|
||||
|
||||
#ifndef MCMC_HPP_
|
||||
#define MCMC_HPP_
|
||||
|
||||
#include <alglib/ap.h>
|
||||
#include <alglib/alglibmisc.h>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
using namespace alglib;
|
||||
using namespace std;
|
||||
|
||||
typedef real_1d_array vec;
|
||||
typedef real_2d_array vec2;
|
||||
typedef integer_1d_array ivec;
|
||||
typedef integer_2d_array ivec2;
|
||||
|
||||
class mcmc {
|
||||
int npar; double(*loglikelihood)(vec&,void*); void*ptr;
|
||||
hqrndstate rnd;
|
||||
public:
|
||||
int nrun,nburn,nwk,ncheck; double avalue;
|
||||
mcmc(int np,double (*f)(vec&,void*),void *p);
|
||||
void setseed(int s1,int s2 ){hqrndseed(s1,s2,rnd);}
|
||||
void run( vec& , vec& , const char*fname="mcmc.dat");
|
||||
virtual ~mcmc();
|
||||
};
|
||||
|
||||
#endif /* MCMC_HPP_ */
|
||||
289
src/inc/mod.hpp
289
src/inc/mod.hpp
@ -1,289 +0,0 @@
|
||||
/*
|
||||
* mod.hpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#ifndef MOD_HPP_
|
||||
#define MOD_HPP_
|
||||
|
||||
#include "def.hpp"
|
||||
#include <alglib/specialfunctions.h>
|
||||
#include <alglib/solvers.h>
|
||||
|
||||
class mod {
|
||||
|
||||
void _CSfq( vec , vector<vec2>& , vector<vec2>& );
|
||||
|
||||
public:
|
||||
|
||||
int n,nfq,npar;
|
||||
double dt,f1,f2;
|
||||
vec lc,lce,t,FqL;
|
||||
vector<vec2> Cfq,Sfq,Cv,Cfq2,Sfq2;
|
||||
|
||||
vec2 C,Ci,I,yyT,yyTmC,dC,C2,C3;
|
||||
vec Cilc;
|
||||
|
||||
ae_int_t info; densesolverreport rep;
|
||||
|
||||
mod();
|
||||
virtual ~mod();
|
||||
|
||||
virtual void _C( vec );
|
||||
virtual void _dC( vec , int );
|
||||
virtual void step_pars( int , vec& , vec& );
|
||||
virtual void print_pars( vec& , vec& );
|
||||
virtual void print_pars( vec& , vec& , vec& );
|
||||
virtual void what_pars( int&, int& );
|
||||
|
||||
void setlc(){ lc.setlength(n); lce.setlength(n); t.setlength(n);}
|
||||
void init( vec , int );
|
||||
double loglikelihood( vec );
|
||||
void dlikelihood( vec , double& , vec& , vec2& );
|
||||
void optimize( vec& , vec& );
|
||||
};
|
||||
|
||||
// ----------------------------------------- //
|
||||
// ---------- mod container ---------------- //
|
||||
// ----------------------------------------- //
|
||||
|
||||
template <class M>
|
||||
class Mod {
|
||||
|
||||
int nmod,nfq;
|
||||
vector<M> mods;
|
||||
|
||||
public:
|
||||
int npar;
|
||||
Mod( vector<lcurve> inlc , vec fqL ){
|
||||
int i;
|
||||
nmod = inlc.size();
|
||||
for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( inlc[i] , fqL ) ); }
|
||||
npar = mods[0].npar;
|
||||
nfq = mods[0].nfq;
|
||||
}
|
||||
Mod( vector<lcurve> lc1, vector<lcurve> lc2 , vec fqL , vec pars ){
|
||||
int i;
|
||||
nmod = lc1.size();
|
||||
for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( lc1[i] , lc2[i] , fqL , pars ) ); }
|
||||
npar = mods[0].npar;
|
||||
nfq = mods[0].nfq;
|
||||
}
|
||||
Mod( vector<lcurve> lc1, vector<lcurve> lc2 , vec fqL ){
|
||||
int i;
|
||||
nmod = lc1.size();
|
||||
for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( lc1[i] , lc2[i] , fqL ) ); }
|
||||
npar = mods[0].npar;
|
||||
nfq = mods[0].nfq;
|
||||
}
|
||||
virtual ~Mod(){}
|
||||
|
||||
double loglikelihood( vec p ){
|
||||
double l=0; for( int i=0 ; i<nmod ; i++ ){ l += mods[i].loglikelihood(p); }
|
||||
return l;
|
||||
}
|
||||
|
||||
void dlikelihood( vec p , double& loglike , vec& grad , vec2& hess ){
|
||||
int i,j,nm;vec g; vec2 h; double l;g.setlength(npar);h.setlength(npar,npar);
|
||||
loglike = 0; for( i=0 ; i<npar ; i++ ){ grad[i]=0; for(j=0;j<npar;j++){ hess(i,j) = 0;}}
|
||||
|
||||
for( nm=0 ; nm<nmod ; nm++ ){
|
||||
mods[nm].dlikelihood( p , l , g , h );
|
||||
loglike += l;
|
||||
for( i=0 ; i<npar ; i++ ){ grad[i] += g[i]; for(j=0;j<npar;j++){ hess(i,j) += h(i,j);}}
|
||||
}
|
||||
}
|
||||
|
||||
void step_pars( int n, vec& dpar , vec& pars ){ mods[0].step_pars( n , dpar , pars );}
|
||||
void print_pars( vec& pars , vec& errs ){ mods[0].print_pars( pars , errs );}
|
||||
void print_pars( vec& pars , vec& errs1 , vec& errs2 ){ mods[0].print_pars( pars , errs1 , errs2 );}
|
||||
void what_pars( int& ip1 , int& ip2 ){ mods[0].what_pars( ip1 , ip2 );}
|
||||
|
||||
double optimize( vec& pars , vec& errs ){
|
||||
int nmax = 400;
|
||||
double tol = 1e-3;
|
||||
|
||||
ae_int_t info; densesolverreport rep;
|
||||
int i,j,n,sing=0,np;
|
||||
double loglike,dpmax,dl,prevl=-1e20;
|
||||
vec tmppars,dpar,grad,prevp; vec2 hess,hessi,ii;
|
||||
tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar); prevp.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;}
|
||||
prevp[i] = pars[i];
|
||||
}
|
||||
|
||||
for( n=1 ; n<= nmax ; n++ ){
|
||||
|
||||
dlikelihood( tmppars , loglike , grad , hess );
|
||||
// ---------- //
|
||||
if( loglike != loglike ){ for( i=0 ; i<npar ; i++ ){ dpar[i]/=10;tmppars[i] = prevp[i] + dpar[i]; } continue;}
|
||||
dl = loglike - prevl;
|
||||
// ---------- //
|
||||
vector<int> ic,iv; sing = 0;
|
||||
for( i=0; i <npar ; i++ ){ if( fabs(grad[i]) < 1e-5 ){ sing=1;ic.push_back(i); }else { iv.push_back(i);}}
|
||||
//if( sing == 1 and iv.size()!=0 ){
|
||||
if( sing == 1 ){
|
||||
np = iv.size();
|
||||
vec g; vec2 h,hi,iii; g.setlength(np); h.setlength(np,np); hi.setlength(np,np);iii.setlength(np,np);
|
||||
for( i=0 ; i<np ; i++ ){
|
||||
iii(i,i) = 1; h(i,i) = hess(iv[i],iv[i]);
|
||||
g[i] = grad[iv[i]]; for(j=0;j<i;j++){ h(i,j) = hess(iv[i],iv[j]);h(j,i)=h(i,j);iii(i,j)=0;iii(j,i)=0;}
|
||||
}
|
||||
spdmatrixcholesky( h , np , false );
|
||||
spdmatrixcholeskysolvem( h , np , false , iii , np , info , rep , hi );
|
||||
|
||||
dpmax = -1e20;
|
||||
for( i=0 ; i<np ; i++ ){
|
||||
dpar[iv[i]] = 0; for( j=0 ; j<np ; j++ ){ dpar[iv[i]] += g[j] * hi(i,j);}
|
||||
dpmax = max( dpmax , fabs(dpar[iv[i]]) );
|
||||
}
|
||||
for( i=0 ; i<int(ic.size()); i++ ){ dpar[ic[i]] = 0;}
|
||||
}else {
|
||||
// --------- //
|
||||
|
||||
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]) );
|
||||
}
|
||||
}
|
||||
for(i=0;i<npar;i++){prevp[i] = tmppars[i];}
|
||||
prevl = loglike;
|
||||
if ( dl < 2 and n>10 ){for(i=0;i<npar;i++){dpar[i]/=2;}}
|
||||
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 );
|
||||
return loglike;
|
||||
}
|
||||
|
||||
double opti( vec pars , vec errs ,int k ){
|
||||
int nmax = 400;
|
||||
double tol = 1e-3;
|
||||
|
||||
ae_int_t info; densesolverreport rep;
|
||||
int i,j,n,np;
|
||||
double loglike,dpmax,dl,prevl=-1e20;
|
||||
vec tmppars,dpar,grad,prevp; vec2 hess,hessi,ii;
|
||||
tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar); prevp.setlength(npar);
|
||||
hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
|
||||
vector<int> ic,iv;
|
||||
vec g; vec2 h,hi,iii;
|
||||
|
||||
dlikelihood( pars , loglike , grad , hess );
|
||||
if( fabs(grad[k])<1e-5 ) { return loglike;}
|
||||
|
||||
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;}
|
||||
prevp[i] = pars[i];
|
||||
if(i==k or fabs(grad[i])<1e-5 ){ic.push_back(i);}else{ iv.push_back(i); }
|
||||
}
|
||||
np = npar - ic.size();
|
||||
g.setlength(np); h.setlength(np,np); hi.setlength(np,np);iii.setlength(np,np);
|
||||
|
||||
|
||||
for( n=1 ; n<= nmax ; n++ ){
|
||||
|
||||
dlikelihood( tmppars , loglike , grad , hess );
|
||||
// ---------- //
|
||||
if( loglike != loglike ){ for( i=0 ; i<npar ; i++ ){ tmppars[i] = prevp[i] + dpar[i]/10; } continue;}
|
||||
dl = loglike - prevl;
|
||||
// ---------- //
|
||||
|
||||
|
||||
for( i=0 ; i<np ; i++ ){
|
||||
iii(i,i) = 1; h(i,i) = hess(iv[i],iv[i]);
|
||||
g[i] = grad[iv[i]]; for(j=0;j<i;j++){ h(i,j) = hess(iv[i],iv[j]);h(j,i)=h(i,j);iii(i,j)=0;iii(j,i)=0;}
|
||||
}
|
||||
spdmatrixcholesky( h , np , false );
|
||||
spdmatrixcholeskysolvem( h , np , false , iii , np , info , rep , hi );
|
||||
|
||||
dpmax = -1e20;
|
||||
for( i=0 ; i<np ; i++ ){
|
||||
dpar[iv[i]] = 0; for( j=0 ; j<np ; j++ ){ dpar[iv[i]] += g[j] * hi(i,j);}
|
||||
dpmax = max( dpmax , fabs(dpar[iv[i]]) );
|
||||
}
|
||||
for( i=0 ; i<int(ic.size()); i++ ){ dpar[ic[i]] = 0;}
|
||||
|
||||
for(i=0;i<npar;i++){prevp[i] = tmppars[i];}
|
||||
prevl = loglike;
|
||||
if ( dl < 2 and n>10 ){for(i=0;i<npar;i++){dpar[i]/=2;}}
|
||||
step_pars( n , dpar , tmppars );
|
||||
|
||||
cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
|
||||
if ( dpmax<tol ) break;
|
||||
}
|
||||
return loglike;
|
||||
}
|
||||
|
||||
void errors_avg( vec& pars , vec& errs ){
|
||||
vec errs1,errs2;errs1.setlength(npar);errs2.setlength(npar);
|
||||
errors( pars , errs1 , 1 );
|
||||
errors( pars , errs2 , -1 );
|
||||
for(int i=0;i<npar;i++){errs[i] = (errs1[i]+errs2[i])/2;}
|
||||
}
|
||||
|
||||
void errors( vec& pars , vec& errs1 , vec& errs2 ){
|
||||
errors( pars , errs1 , 1 );
|
||||
errors( pars , errs2 , -1 );
|
||||
}
|
||||
|
||||
void errors( vec& pars , vec& errs , int sign=1){
|
||||
double tol = 1e-2 , delchi = 1;
|
||||
int i,ip,ip1,ip2;
|
||||
double pu,pd,phalf=0,bestlike,tmplike,dl;
|
||||
vec pars0,errs0,tmpp; pars0.setlength(npar);errs0.setlength(npar);tmpp.setlength(npar);
|
||||
bestlike = optimize( pars , errs );
|
||||
for(i=0;i<npar;i++){pars0[i] = pars[i];errs0[i]=errs[i];}
|
||||
|
||||
what_pars( ip1 , ip2 );
|
||||
|
||||
for( ip=ip1 ; ip<ip2 ; ip++ ){
|
||||
cout << endl << "******* par " << ip << " *******" << endl << endl;
|
||||
if( errs0[ip]>10 ) { continue ;}
|
||||
pd = pars[ip];
|
||||
for( i=0 ;i<npar ; i++ ){tmpp[i] = pars0[i]; }
|
||||
|
||||
i=2; dl = 0;
|
||||
while( dl < delchi ){
|
||||
pu = pars[ip] + sign*i*errs0[ip];
|
||||
tmpp[ip] = pu;
|
||||
tmplike = opti( tmpp , errs , ip );
|
||||
dl = 2*(bestlike - tmplike);
|
||||
i++;
|
||||
}
|
||||
while( fabs(dl-delchi)>tol ){
|
||||
phalf = (pu+pd)/2.0;
|
||||
tmpp[ip] = phalf;
|
||||
tmplike = opti( tmpp , errs , ip );
|
||||
dl = 2*(bestlike - tmplike);
|
||||
cout << "+++ " << bestlike << " " << tmplike << " " << pars[ip] << " "<< phalf << " " << dl << endl;
|
||||
if( dl<delchi ){ pd = phalf;}else{ pu = phalf;}
|
||||
}
|
||||
errs0[ip] = sign* ( phalf - pars0[ip] );
|
||||
}
|
||||
for( i=0 ; i<npar ; i++ ){ pars[i] = pars0[i]; errs[i] = errs0[i];}
|
||||
print_pars( pars0 , errs0 );
|
||||
//optimize( pars0 , errs );
|
||||
}
|
||||
};
|
||||
|
||||
#endif /* MOD_HPP_ */
|
||||
@ -1,37 +0,0 @@
|
||||
/*
|
||||
* psd.hpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#ifndef PSD_HPP_
|
||||
#define PSD_HPP_
|
||||
|
||||
#include "mod.hpp"
|
||||
|
||||
class psd: public mod {
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
|
||||
public:
|
||||
psd( lcurve , vec );
|
||||
virtual ~psd();
|
||||
};
|
||||
|
||||
// ---------------------- //
|
||||
|
||||
class psd10: public mod {
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
|
||||
public:
|
||||
psd10( lcurve , vec );
|
||||
virtual ~psd10();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
};
|
||||
|
||||
#endif /* PSD_HPP_ */
|
||||
@ -1,41 +0,0 @@
|
||||
/*
|
||||
* psd.hpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#ifndef PSD_RMS_HPP_
|
||||
#define PSD_RMS_HPP_
|
||||
|
||||
#include "mod.hpp"
|
||||
|
||||
class psd_rms: public mod {
|
||||
|
||||
double mean,mean2;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
|
||||
public:
|
||||
psd_rms( lcurve , vec );
|
||||
virtual ~psd_rms();
|
||||
};
|
||||
|
||||
// ---------------------- //
|
||||
|
||||
class psd10_rms: public mod {
|
||||
|
||||
double mean,mean2;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
|
||||
public:
|
||||
psd10_rms( lcurve , vec );
|
||||
virtual ~psd10_rms();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
};
|
||||
|
||||
#endif /* PSD_RMS_HPP_ */
|
||||
@ -1,54 +0,0 @@
|
||||
/*
|
||||
* psdlag.hpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "mod.hpp"
|
||||
|
||||
#ifndef PSDLAG_HPP_
|
||||
#define PSDLAG_HPP_
|
||||
|
||||
class psdlag : public mod {
|
||||
|
||||
int n1;
|
||||
ivec ip1,ip2,icx,iphi;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dP1( vec , int );
|
||||
void _dP2( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
public:
|
||||
psdlag( lcurve , lcurve , vec );
|
||||
virtual ~psdlag();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
};
|
||||
|
||||
class psdlag10 : public mod {
|
||||
|
||||
int n1;
|
||||
ivec ip1,ip2,icx,iphi;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dP1( vec , int );
|
||||
void _dP2( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
public:
|
||||
psdlag10( lcurve , lcurve , vec );
|
||||
virtual ~psdlag10();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
void what_pars( int&, int& );
|
||||
};
|
||||
|
||||
#endif /* PSDLAG_HPP_ */
|
||||
@ -1,56 +0,0 @@
|
||||
/*
|
||||
* psdlag.hpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "mod.hpp"
|
||||
|
||||
#ifndef PSDLAG_RMS_HPP_
|
||||
#define PSDLAG_RMS_HPP_
|
||||
|
||||
class psdlag_rms : public mod {
|
||||
|
||||
int n1;
|
||||
ivec ip1,ip2,icx,iphi;
|
||||
double mean1,mean2,mean1sq,mean2sq,mean12;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dP1( vec , int );
|
||||
void _dP2( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
public:
|
||||
psdlag_rms( lcurve , lcurve , vec );
|
||||
virtual ~psdlag_rms();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
};
|
||||
|
||||
class psdlag10_rms : public mod {
|
||||
|
||||
int n1;
|
||||
ivec ip1,ip2,icx,iphi;
|
||||
double mean1,mean2,mean1sq,mean2sq,mean12;
|
||||
|
||||
void _C( vec );
|
||||
void _dC( vec , int );
|
||||
void _dP1( vec , int );
|
||||
void _dP2( vec , int );
|
||||
void _dCx( vec , int );
|
||||
void _dPhi( vec , int );
|
||||
|
||||
public:
|
||||
psdlag10_rms( lcurve , lcurve , vec );
|
||||
virtual ~psdlag10_rms();
|
||||
|
||||
void step_pars( int , vec& , vec& );
|
||||
void print_pars( vec& , vec& );
|
||||
void what_pars( int&, int& );
|
||||
};
|
||||
|
||||
#endif /* PSDLAG_RMS_HPP_ */
|
||||
227
src/lag.cpp
227
src/lag.cpp
@ -1,227 +0,0 @@
|
||||
/*
|
||||
* lag.cpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/lag.hpp"
|
||||
|
||||
lag::lag( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
lc1.demean(); lc2.demean();
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 2 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
p1[i] = pars[i];
|
||||
p2[i] = pars[i+nfq];
|
||||
icx[i] = i;
|
||||
iphi[i] = i+nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
lag::~lag() {}
|
||||
|
||||
void lag::_C( vec p ){
|
||||
int i,j,k; double d;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void lag::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dCx( p , k );
|
||||
}else{ _dPhi( p , k-nfq );}
|
||||
}
|
||||
|
||||
void lag::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void lag::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=p[k],phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
|
||||
pars[i] += dpar[i]/((n<10)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
if( pars[i]<0 ){ pars[i]*=-1; pars[i+nfq] += M_PI;}
|
||||
while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
|
||||
while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
|
||||
}
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
|
||||
|
||||
// ****************************************** //
|
||||
|
||||
|
||||
lag10::lag10( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
lc1.demean(); lc2.demean();
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 2 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
p1[i] = pow(10,pars[i]);
|
||||
p2[i] = pow(10,pars[i+nfq]);
|
||||
icx[i] = i;
|
||||
iphi[i] = i+nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
lag10::~lag10() {}
|
||||
|
||||
void lag10::_C( vec p ){
|
||||
int i,j,k; double d; for(i=0;i<nfq;i++){p[i] = pow(10,p[i]);}
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); }
|
||||
d += f1*p1[0]*Cfq2[0](i,i) + f2*p1[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); }
|
||||
d += f1*p1[0]*Cfq2[0](i,j) + f2*p1[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); }
|
||||
d += f1*p2[0]*Cfq2[0](i,i) + f2*p2[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); }
|
||||
d += f1*p2[0]*Cfq2[0](i,j) + f2*p2[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); }
|
||||
d += f1*p[icx[0]] * ( Cfq2[0](i,j)*cos(p[iphi[0]]) - Sfq2[0](i,j)*sin(p[iphi[0]]) );
|
||||
d += f2*p[icx[nfq-1]] * ( Cfq2[2](i,j)*cos(p[iphi[nfq-1]]) - Sfq2[2](i,j)*sin(p[iphi[nfq-1]]) );
|
||||
C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void lag10::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dCx( p , k );
|
||||
}else{ _dPhi( p , k-nfq );}
|
||||
}
|
||||
|
||||
void lag10::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[k]);
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void lag10::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=pow(10,p[k]),phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag10::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
||||
//pars[i] += dpar[i];
|
||||
pars[i] += dpar[i]/((n<5)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag10::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
|
||||
while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
|
||||
}
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
|
||||
void lag10::what_pars( int& ip1 , int& ip2 ){
|
||||
ip1 = nfq; ip2 = npar;
|
||||
}
|
||||
230
src/lag_rms.cpp
230
src/lag_rms.cpp
@ -1,230 +0,0 @@
|
||||
/*
|
||||
* lag.cpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/lag_rms.hpp"
|
||||
|
||||
lag_rms::lag_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
mean1 = lc1.demean(); mean2 = lc2.demean();
|
||||
mean1sq = mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 2 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
p1[i] = pars[i] * mean1sq;
|
||||
p2[i] = pars[i+nfq] * mean2sq;
|
||||
icx[i] = i;
|
||||
iphi[i] = i+nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
lag_rms::~lag_rms() {}
|
||||
|
||||
void lag_rms::_C( vec p ){
|
||||
int i,j,k; double d;
|
||||
for( i=0 ; i<nfq ; i++ ){ p[icx[i]] *= mean12;}
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void lag_rms::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dCx( p , k );
|
||||
}else{ _dPhi( p , k-nfq );}
|
||||
}
|
||||
|
||||
void lag_rms::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void lag_rms::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=p[k]*mean12,phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag_rms::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
|
||||
pars[i] += dpar[i]/((n<10)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag_rms::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
if( pars[i]<0 ){ pars[i]*=-1; pars[i+nfq] += M_PI;}
|
||||
while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
|
||||
while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
|
||||
}
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
|
||||
|
||||
// ****************************************** //
|
||||
|
||||
|
||||
lag10_rms::lag10_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
mean1 = lc1.demean(); mean2 = lc2.demean();
|
||||
mean1sq = mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 2 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
p1[i] = pow(10,pars[i]) * mean1sq;
|
||||
p2[i] = pow(10,pars[i+nfq]) * mean2sq;
|
||||
icx[i] = i;
|
||||
iphi[i] = i+nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
lag10_rms::~lag10_rms() {}
|
||||
|
||||
void lag10_rms::_C( vec p ){
|
||||
int i,j,k; double d; for(i=0;i<nfq;i++){p[icx[i]] = pow(10,p[i]) * mean12;}
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); }
|
||||
d += f1*p1[0]*Cfq2[0](i,i) + f2*p1[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); }
|
||||
d += f1*p1[0]*Cfq2[0](i,j) + f2*p1[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); }
|
||||
d += f1*p2[0]*Cfq2[0](i,i) + f2*p2[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); }
|
||||
d += f1*p2[0]*Cfq2[0](i,j) + f2*p2[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); }
|
||||
d += f1*p[icx[0]] * ( Cfq2[0](i,j)*cos(p[iphi[0]]) - Sfq2[0](i,j)*sin(p[iphi[0]]) );
|
||||
d += f2*p[icx[nfq-1]] * ( Cfq2[2](i,j)*cos(p[iphi[nfq-1]]) - Sfq2[2](i,j)*sin(p[iphi[nfq-1]]) );
|
||||
C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void lag10_rms::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dCx( p , k );
|
||||
}else{ _dPhi( p , k-nfq );}
|
||||
}
|
||||
|
||||
void lag10_rms::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[k])*mean12;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void lag10_rms::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=pow(10,p[k])*mean12,phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag10_rms::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
||||
//pars[i] += dpar[i];
|
||||
pars[i] += dpar[i]/((n<5)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void lag10_rms::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
|
||||
while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
|
||||
}
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
|
||||
void lag10_rms::what_pars( int& ip1 , int& ip2 ){
|
||||
ip1 = nfq; ip2 = npar;
|
||||
}
|
||||
221
src/main.cpp
221
src/main.cpp
@ -1,221 +0,0 @@
|
||||
/*
|
||||
* main.cpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/main.hpp"
|
||||
|
||||
/**************************************
|
||||
* Main function
|
||||
* Checks if there is an input and
|
||||
* raises an error.
|
||||
* otherwise, call do_work(fname)
|
||||
**************************************/
|
||||
int main( int argc , char* argv[] ){
|
||||
if ( argc < 2 ){
|
||||
cerr << "** Error ** No input file given." << endl;
|
||||
usage();
|
||||
}else {
|
||||
do_work( argv[1] );
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************
|
||||
* Main function that does the work
|
||||
* INPUT: char* fname: input file name
|
||||
* RETURN: void
|
||||
**************************************/
|
||||
|
||||
void do_work( char* fname ){
|
||||
|
||||
/************ Reading the input file ***********/
|
||||
ifstream fp(fname);
|
||||
string line,mcmcfile,colon; stringstream ss;
|
||||
int i,nfiles,nfq,mode,npar,bin1,bin2,fit_type,nrun,nburn,nwk,ref; bool strict;
|
||||
getline(fp,line);ss.str(line); ss >> nfiles; ss.clear();
|
||||
vector<string> files(nfiles);ivec secL;secL.setlength(nfiles);for(i=0;i<nfiles;i++){
|
||||
getline(fp,line); ss.str(line); ss >> files[i] >> secL[i]; ss.clear();}
|
||||
getline(fp,line); ss.str(line); ss >> nfq;vec fqL;fqL.setlength(nfq+1);
|
||||
for(i=0;i<=nfq;i++){ss >> fqL[i];} ss.clear();
|
||||
getline(fp,line); ss.str(line); ss >> mode; ss.clear();
|
||||
npar = (mode==0)?4*nfq:nfq; vec pars; pars.setlength(npar);
|
||||
getline(fp,line); ss.str(line); for(i=0;i<npar;i++){ss >> pars[i];} ss.clear();
|
||||
getline(fp,line); ss.str(line); ss >> ref >> colon >> bin2; ss.clear();
|
||||
if(colon.compare(0,1,":") == 0 ){
|
||||
bin1 = atoi(colon.substr(1).c_str());
|
||||
}else{
|
||||
bin1 = ref; ref = -10; bin2 = atoi(colon.c_str());
|
||||
}
|
||||
if(mode!=0){bin1 = (mode==-1)?-1:(mode-1);bin2=-1;}
|
||||
getline(fp,line); ss.str(line); ss >> fit_type; ss.clear();
|
||||
getline(fp,line); ss.str(line); ss >> i; strict = i!=0; ss.clear();
|
||||
getline(fp,line); ss.str(line); ss >> nrun >> nburn >> nwk >> mcmcfile; ss.clear();
|
||||
fp.close();
|
||||
/********* END Reading the input file **********/
|
||||
|
||||
/********** Read the light curves **************/
|
||||
vector<vector<lcurve> > LC,LC_ref;
|
||||
if ( ref == -10 ){
|
||||
for(i=0;i<nfiles;i++) readLC( LC , files[i] , secL[i] , bin1 , bin2 , strict );
|
||||
}else {
|
||||
for(i=0;i<nfiles;i++) {
|
||||
if( i!=ref ) { readLC( LC , files[i] , secL[i] , bin2 , -1 , strict ); }
|
||||
}
|
||||
readLC( LC_ref , files[ref] , secL[ref] , bin1 , -1 , strict );
|
||||
}
|
||||
/******** END Read the light curves ************/
|
||||
|
||||
if( mode > 0 or mode==-1 ){
|
||||
vec errs; errs.setlength(nfq);
|
||||
vector<lcurve> lc1; for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC[i][0]); }
|
||||
Mod<psd10_rms> p1( lc1 , fqL );
|
||||
|
||||
p1.optimize( pars , errs );
|
||||
if( fit_type==1 ) p1.errors( pars , errs );
|
||||
}else if( mode == 0 ){
|
||||
vec ec,e1,e2,errs; ec.setlength(2*nfq);e1.setlength(nfq);e2.setlength(nfq);errs.setlength(4*nfq);
|
||||
vector<lcurve> lc1,lc2;
|
||||
if( ref==-10 ){ // standard way of getting reference
|
||||
for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC[i][0]); lc2.push_back( LC[i][1]); }
|
||||
}else{
|
||||
for( i=0 ; i<int(LC_ref.size()) ; i++ ){ lc2.push_back( LC[i][0]);}
|
||||
for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC_ref[i][0]);}
|
||||
}
|
||||
Mod<psd10_rms> p1( lc1 , fqL ), p2( lc2 , fqL );
|
||||
vec pars1,pars2;pars1.setlength(nfq);pars2.setlength(nfq);for(i=0;i<nfq;i++){pars1[i]=pars[i];pars2[i]=pars[i+nfq];}
|
||||
|
||||
p1.optimize( pars1 , e1 );
|
||||
p2.optimize( pars2 , e2 );
|
||||
if( fit_type==1 ) {
|
||||
p1.errors( pars1 , e1 );
|
||||
p2.errors( pars2 , e2 );
|
||||
}
|
||||
|
||||
vec pc; pc.setlength(2*nfq);
|
||||
for(i=0;i<nfq;i++){pc[i]=pars[i+2*nfq];pc[i+nfq]=pars[i+3*nfq];pars[i]=pars1[i];pars[i+nfq]=pars2[i];pc[i]=(pars1[i]+pars2[i])/2;}
|
||||
Mod<lag10_rms> l( lc1 , lc2 , fqL , pars );
|
||||
l.optimize( pc , ec );
|
||||
if( fit_type==1 ) l.errors( pc , ec );
|
||||
if( fit_type==2 ) {
|
||||
mcmc mc( 2*nfq , mcmc_lag10 , (void*)&l );
|
||||
mc.nrun=nrun; mc.nburn=nburn; mc.nwk=nwk;
|
||||
mc.run( pc , ec , mcmcfile.c_str() );
|
||||
}
|
||||
|
||||
for(i=0;i<nfq;i++){ pars[i+2*nfq]=pc[i]; pars[i+3*nfq]=pc[i+nfq];}
|
||||
for(i=0;i<nfq;i++){ errs[i] = e1[i]; errs[i+nfq] = e2[i]; errs[i+2*nfq]=ec[i]; errs[i+3*nfq]=ec[i+nfq];}
|
||||
|
||||
Mod<psdlag10_rms> pl( lc1 , lc2 , fqL );
|
||||
if( fit_type==3 ) pl.errors( pars , errs );
|
||||
if( fit_type==4 ) {
|
||||
mcmc mc( 4*nfq , mcmc_psdlag10 , (void*)&pl );
|
||||
mc.nrun=nrun; mc.nburn=nburn; mc.nwk=nwk;
|
||||
mc.run( pars , errs , mcmcfile.c_str() );
|
||||
}
|
||||
|
||||
pl.print_pars( pars , errs );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**************************************
|
||||
* Usage function
|
||||
* prints a sample input text
|
||||
* INPUT: none
|
||||
* RETURN: void
|
||||
**************************************/
|
||||
void usage(){
|
||||
cerr << setw(40) << setfill('*') << "\n";
|
||||
cout << setw(30) << left << "1" << "# of lc" << endl;
|
||||
cout << setw(30) << left << "lc.dat 0" << "name1 secL1" << endl;
|
||||
cout << setw(30) << left << "3 1e-5 1e-4 5e-4 1e-3" << "nfq fqL1 fqL2.." << endl;
|
||||
cout << setw(30) << left << "1" << "mode? 0:lag, n:psd" << endl;
|
||||
cout << setw(30) << left << "1 1 1 1 1 1 1 1 1 1 1 1" << "init. pars" << endl;
|
||||
cout << setw(30) << left << "0 1" << "bin1, bin2,-1: all" << endl;
|
||||
cout << setw(30) << left << "0" << "fit_type" << endl;
|
||||
cout << setw(30) << left << "0" << "strict" << endl;
|
||||
cout << setw(30) << left << "100 50 50 mcmc.dat" << "nrun,nburn,nwk,file"<< endl;
|
||||
cerr << setw(40) << setfill('*') << "\n";
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/***************************************
|
||||
* READ LIGHTCURVES
|
||||
* Read a standard light curve file that has a description at the top
|
||||
* the first line is something like: # dT nlc b0 b1 b2 ... b_lc (b0 .. are bins boundaries)
|
||||
* Then the content with: time lc1 err1 lc2 err2 ... etc
|
||||
*
|
||||
* INPUT:
|
||||
* LC : a vector of vectors of light curves, where the first dimension is the number of
|
||||
* light curves (i.e. segments), this is incremented in this function when reading a new file.
|
||||
* The other dimension is either 1 for psd and 2 for lags, which contains either one light curve
|
||||
* or two for the lag. Passed by reference
|
||||
* fname : name of the light curve file
|
||||
* secL : the section length to split the light curves to.
|
||||
* b1 : the bin1 to read for the first light curve. 0-based, -1 for all
|
||||
* b2 : the bin2 to read, 0-based, -1 for none (.i.e just psd)
|
||||
* strict : a boolean on whether to be strict or not in segmenting the light curve
|
||||
* e.g. for secL=220: 0 gives 100,100 and 1 gives 100,120
|
||||
*
|
||||
**************************************/
|
||||
void readLC( vector<vector<lcurve> >&LC , string fname , int secL , int b1 , int b2 , bool strict ){
|
||||
|
||||
/* Initialize file stream, define some variables */
|
||||
ifstream fp(fname.c_str()); string line,sdum; stringstream ss;
|
||||
int i,j,nlc,n,nsec,ns,sl; double dt;
|
||||
|
||||
/* read dt from the description line */
|
||||
getline(fp,line); ss.str(line); ss >> sdum >> dt >> nlc; ss.clear();
|
||||
|
||||
/* get the number of lines */
|
||||
n=0;while(getline(fp,line)){if(line!="")n++;} fp.clear(); fp.seekg(0,ios::beg);getline(fp,line);
|
||||
|
||||
/* calculate section lengths */
|
||||
if(secL==0){secL=n;}
|
||||
nsec = n/secL; vector<int> secl(nsec,secL); if(not strict) secl[nsec-1] += n-nsec*secL;
|
||||
|
||||
/* Start looping through the number of sections */
|
||||
for(ns=0;ns<nsec;ns++){
|
||||
sl = secl[ns]; vector<lcurve> Lc;
|
||||
vec t; vector<vec> lc(nlc),lce(nlc); t.setlength( sl ); for(i=0;i<nlc;i++){lc[i].setlength(sl);lce[i].setlength(sl);}
|
||||
|
||||
// Read the first chunk
|
||||
for( i=0 ; i<sl ; i++ ){
|
||||
getline(fp,line);ss.str(line);ss>>t[i];for(j=0;j<nlc;j++){ss>>lc[j][i]>>lce[j][i];} ss.clear();
|
||||
}
|
||||
|
||||
// if we want the whole light curves in bin1
|
||||
if( b1==-1 ){
|
||||
vec lc_,lce_; lc_.setlength( sl );lce_.setlength( sl );for(i=0;i<sl;i++){lc_[i]=0;lce_[i]=0;}
|
||||
for( i=0 ; i<sl ; i++ ){
|
||||
for( j=0 ; j<nlc ; j++ ){ if(j==b2){continue;} lc_[i] += lc[j][i]; lce_[i] += lce[j][i]*lce[j][i]; }
|
||||
lce_[i] = sqrt(lce_[i]);
|
||||
}
|
||||
Lc.push_back(lcurve(t,lc_,lce_,dt));
|
||||
|
||||
// else: i.e. light curve 1 is bin1
|
||||
}else {
|
||||
Lc.push_back(lcurve(t,lc[b1],lce[b1],dt));
|
||||
}
|
||||
// Do we need a second light curve?
|
||||
if(b2!=-1){Lc.push_back(lcurve(t,lc[b2],lce[b2],dt));}
|
||||
|
||||
// push the resulting vector<lcurve> to LC
|
||||
LC.push_back(Lc);
|
||||
}
|
||||
fp.close();
|
||||
}
|
||||
13
src/makefile
13
src/makefile
@ -1,13 +0,0 @@
|
||||
|
||||
# original
|
||||
#libdir='/eos/azoghbi/soft/usr/lib'
|
||||
#incdir='/eos/azoghbi/soft/usr/include'
|
||||
# most systems
|
||||
#libdir='/home/caes/science/psdlag-agn/src/inc'
|
||||
#incdir='/home/caes/science/psdlag-agn/src/inc'
|
||||
# thor
|
||||
libdir='/home/research/xaw5719/psdlag-agn/src/inc'
|
||||
alglibdir='/home/research/xaw5719/psdlag-agn/src/inc/alglib'
|
||||
incdir='/home/research/xaw5719/psdlag-agn/src/inc'
|
||||
psdlag:
|
||||
g++ *cpp inc/alglib/*cpp -o psdlag -O3 -Wall -I${incdir} -L${libdir} -L${alglibdir} -I${alglibdir}
|
||||
73
src/mcmc.cpp
73
src/mcmc.cpp
@ -1,73 +0,0 @@
|
||||
/*
|
||||
* mcmc.cpp
|
||||
*
|
||||
* Created on: May 14, 2013
|
||||
* Author: abduz
|
||||
*/
|
||||
|
||||
#include "inc/mcmc.hpp"
|
||||
|
||||
mcmc::mcmc(int np,double (*f)(vec&,void*),void *p) {
|
||||
npar = np;
|
||||
loglikelihood = f;
|
||||
ptr = p;
|
||||
hqrndrandomize(rnd);
|
||||
nrun = 100; nburn = 100; nwk = 100; ncheck = 10;avalue = 2.0;
|
||||
}
|
||||
|
||||
mcmc::~mcmc() {}
|
||||
|
||||
void mcmc::run( vec& p, vec& pe, const char*fname){
|
||||
int i,j,k,nacc,nr,n,nsum;
|
||||
vec2 pars; pars.setlength(nwk,npar);
|
||||
vec prob,currp,bestp,sum,sum2;
|
||||
ivec2 compj;
|
||||
double bestlike,tmplike=-1e20,z,q,a1=1./sqrt(avalue),a2=sqrt(avalue);;
|
||||
prob.setlength(nwk);currp.setlength(npar);bestp.setlength(npar);sum.setlength(npar);sum2.setlength(npar);
|
||||
compj.setlength(nwk,nwk-1);
|
||||
nsum=0;for(i=0;i<npar;i++){sum[i]=0;sum2[i]=0;}
|
||||
ofstream fout(fname);
|
||||
|
||||
for(k=0;k<nwk;k++){j=0;for(i=0;i<k;i++){compj[k][j]=i;j++;}for(i=k+1;i<nwk;i++){compj[k][j]=i;j++;}}
|
||||
|
||||
bestlike = -1e20;
|
||||
for(k=0;k<nwk;k++){
|
||||
//pars[k][0] = p[0];currp[0] = p[0];
|
||||
for(i=0;i<npar;i++){pars[k][i] = p[i]*(1.0+hqrndnormal(rnd)*0.05); currp[i]=pars[k][i];}
|
||||
prob[k] = loglikelihood( currp , ptr );
|
||||
if(bestlike<prob[k]){bestlike=prob[k];for(i=0;i<npar;i++){bestp[i]=pars[k][i];}}
|
||||
}
|
||||
|
||||
nr = 0;nacc = 0;
|
||||
for(n=1;n<=(nrun+nburn);n++){
|
||||
for(k=0;k<nwk;k++){
|
||||
nr++;j = compj[k][hqrnduniformi(rnd,nwk-1)];
|
||||
z = hqrnduniformr(rnd)*(a2-a1) + a1; z = z*z;
|
||||
for(i=0;i<npar;i++){currp[i]=pars[j][i] + z*(pars[k][i]-pars[j][i]);}
|
||||
tmplike = loglikelihood(currp,ptr);
|
||||
q = pow(z,npar-1)*exp(tmplike-prob[k]);
|
||||
if( hqrnduniformr(rnd)<=q and tmplike>-1e10 ){
|
||||
nacc++;for(i=0;i<npar;i++){pars[k][i]=currp[i];} prob[k] = tmplike;
|
||||
if(bestlike<tmplike){bestlike=tmplike;for(i=0;i<npar;i++){bestp[i]=currp[i];}}
|
||||
if(n>nburn) {
|
||||
for(i=0;i<npar;i++){
|
||||
fout << currp[i] << " ";
|
||||
sum[i] += currp[i];sum2[i] += currp[i]*currp[i]; nsum++;
|
||||
} fout << tmplike << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
if(n%ncheck==0){
|
||||
q = 1.0*nacc/nr; nr=0;nacc=0;
|
||||
cout << setw(8) << n << setw(12) << q << setw(12) << tmplike << setw(12) << bestlike << endl;
|
||||
}
|
||||
}
|
||||
for(i=0;i<npar;i++){
|
||||
p[i]=bestp[i];
|
||||
pe[i] = sqrt((sum2[i] - sum[i]*sum[i]/nsum)/(nsum+1));
|
||||
}cout << endl;
|
||||
|
||||
|
||||
fout.close();
|
||||
|
||||
}
|
||||
236
src/mod.cpp
236
src/mod.cpp
@ -1,236 +0,0 @@
|
||||
/*
|
||||
* 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 = 1.5; f2 = 0.6;
|
||||
f1 = 0.5; f2 = .0;
|
||||
_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;
|
||||
}
|
||||
105
src/psd.cpp
105
src/psd.cpp
@ -1,105 +0,0 @@
|
||||
/*
|
||||
* psd.cpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/psd.hpp"
|
||||
|
||||
psd::psd( lcurve inlc , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n = inlc.len;
|
||||
dt = inlc.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
inlc.demean();
|
||||
int i;
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
lc[i] = inlc.lc[i];
|
||||
lce[i] = inlc.lce[i]*inlc.lce[i];
|
||||
t[i] = inlc.t[i];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 1 );
|
||||
// ------------------------------------------ //
|
||||
}
|
||||
|
||||
psd::~psd() {}
|
||||
|
||||
void psd::_C( vec p ){
|
||||
int i,j,k; double d;
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
|
||||
d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
|
||||
d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psd::_dC( vec p , int k ){
|
||||
int i,j;
|
||||
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
}
|
||||
|
||||
|
||||
// ********************************************** //
|
||||
// ********************************************** //
|
||||
|
||||
psd10::psd10( lcurve inlc , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n = inlc.len;
|
||||
dt = inlc.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
inlc.demean();
|
||||
int i;
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
lc[i] = inlc.lc[i];
|
||||
lce[i] = inlc.lce[i]*inlc.lce[i];
|
||||
t[i] = inlc.t[i];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 1 );
|
||||
// ------------------------------------------ //
|
||||
}
|
||||
|
||||
psd10::~psd10() {}
|
||||
|
||||
void psd10::_C( vec p ){
|
||||
int i,j,k; double d; for(i=0;i<npar;i++){p[i]=pow(10,p[i]);}
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
|
||||
d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
|
||||
d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psd10::_dC( vec p , int k ){
|
||||
int i,j; double d = log(10)*pow(10,p[k]);
|
||||
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = d*Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
}
|
||||
|
||||
void psd10::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
||||
pars[i] += dpar[i];
|
||||
}
|
||||
}
|
||||
108
src/psd_rms.cpp
108
src/psd_rms.cpp
@ -1,108 +0,0 @@
|
||||
/*
|
||||
* psd.cpp
|
||||
*
|
||||
* Created on: May 31, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/psd_rms.hpp"
|
||||
|
||||
psd_rms::psd_rms( lcurve inlc , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n = inlc.len;
|
||||
dt = inlc.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
mean = inlc.demean();
|
||||
mean2 = mean*mean;
|
||||
int i;
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
lc[i] = inlc.lc[i];
|
||||
lce[i] = inlc.lce[i]*inlc.lce[i];
|
||||
t[i] = inlc.t[i];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 1 );
|
||||
// ------------------------------------------ //
|
||||
}
|
||||
|
||||
psd_rms::~psd_rms() {}
|
||||
|
||||
void psd_rms::_C( vec p ){
|
||||
int i,j,k; double d;
|
||||
for( i=0 ; i<npar ; i++ ){ p[i] *= mean2; }
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
|
||||
d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
|
||||
d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psd_rms::_dC( vec p , int k ){
|
||||
int i,j;
|
||||
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
}
|
||||
|
||||
|
||||
// ********************************************** //
|
||||
// ********************************************** //
|
||||
|
||||
psd10_rms::psd10_rms( lcurve inlc , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n = inlc.len;
|
||||
dt = inlc.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
mean = inlc.demean();
|
||||
mean2 = mean*mean;
|
||||
int i;
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
lc[i] = inlc.lc[i];
|
||||
lce[i] = inlc.lce[i]*inlc.lce[i];
|
||||
t[i] = inlc.t[i];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 1 );
|
||||
// ------------------------------------------ //
|
||||
}
|
||||
|
||||
psd10_rms::~psd10_rms() {}
|
||||
|
||||
void psd10_rms::_C( vec p ){
|
||||
int i,j,k; double d; for(i=0;i<npar;i++){p[i]=pow(10,p[i])*mean2;}
|
||||
for( i=0 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
|
||||
d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
|
||||
d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
|
||||
d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
|
||||
C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psd10_rms::_dC( vec p , int k ){
|
||||
int i,j; double d = log(10)*pow(10,p[k])*mean2;
|
||||
for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = d*Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
}
|
||||
|
||||
void psd10_rms::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
||||
pars[i] += dpar[i];
|
||||
}
|
||||
}
|
||||
252
src/psdlag.cpp
252
src/psdlag.cpp
@ -1,252 +0,0 @@
|
||||
/*
|
||||
* psdlag.cpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/psdlag.hpp"
|
||||
|
||||
psdlag::psdlag( lcurve lc1, lcurve lc2 , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
lc1.demean(); lc2.demean();
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 4 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
ip1[i] = i;
|
||||
ip2[i] = i+nfq;
|
||||
icx[i] = i+2*nfq;
|
||||
iphi[i] = i+3*nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
psdlag::~psdlag() {}
|
||||
|
||||
void psdlag::_C( vec p ){
|
||||
int i,j,k; double d;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dP1( p , k );
|
||||
}else if( k<2*nfq ){ _dP2( p , k%nfq );
|
||||
}else if( k<3*nfq ){ _dCx( p , k%nfq );
|
||||
}else{ _dPhi( p , k%nfq );}
|
||||
}
|
||||
|
||||
void psdlag::_dP1( vec p , int k ){
|
||||
int i,j;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag::_dP2( vec p , int k ){
|
||||
int i,j;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=p[icx[k]],phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
|
||||
pars[i] += dpar[i]/((n<10)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void psdlag::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
if( pars[i+2*nfq]<0 ){ pars[i+2*nfq]*=-1; pars[i+3*nfq] += M_PI;}
|
||||
while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
|
||||
while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
}
|
||||
|
||||
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ //
|
||||
|
||||
psdlag10::psdlag10( lcurve lc1, lcurve lc2 , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
lc1.demean(); lc2.demean();
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 4 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
ip1[i] = i;
|
||||
ip2[i] = i+nfq;
|
||||
icx[i] = i+2*nfq;
|
||||
iphi[i] = i+3*nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
psdlag10::~psdlag10() {}
|
||||
|
||||
void psdlag10::_C( vec p ){
|
||||
int i,j,k; double d; for(i=0;i<3*nfq;i++){ p[i] = pow(10,p[i]);}
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dP1( p , k );
|
||||
}else if( k<2*nfq ){ _dP2( p , k%nfq );
|
||||
}else if( k<3*nfq ){ _dCx( p , k%nfq );
|
||||
}else{ _dPhi( p , k%nfq );}
|
||||
}
|
||||
|
||||
void psdlag10::_dP1( vec p , int k ){
|
||||
int i,j; double p1 = log(10)*pow(10,p[k]);
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = p1 * Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10::_dP2( vec p , int k ){
|
||||
int i,j; double p2 = log(10) * pow(10,p[ip2[k]]);
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = p2 * Cfq[k](i,j); dC(j,i) = dC(i,j);}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]], cx=log(10)*pow(10,p[icx[k]]);
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=pow(10,p[icx[k]]),phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
||||
//pars[i] += dpar[i];
|
||||
pars[i] += dpar[i]/((n<5)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void psdlag10::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
|
||||
while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
|
||||
}
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
|
||||
void psdlag10::what_pars( int& ip1 , int& ip2 ){
|
||||
ip1 = 3*nfq; ip2 = 4*nfq;
|
||||
}
|
||||
@ -1,255 +0,0 @@
|
||||
/*
|
||||
* psdlag.cpp
|
||||
*
|
||||
* Created on: Jun 1, 2013
|
||||
* Author: azoghbi
|
||||
*/
|
||||
|
||||
#include "inc/psdlag_rms.hpp"
|
||||
|
||||
psdlag_rms::psdlag_rms( lcurve lc1, lcurve lc2 , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
mean1 = lc1.demean(); mean2 = lc2.demean();
|
||||
mean1sq = mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 4 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
ip1[i] = i;
|
||||
ip2[i] = i+nfq;
|
||||
icx[i] = i+2*nfq;
|
||||
iphi[i] = i+3*nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
psdlag_rms::~psdlag_rms() {}
|
||||
|
||||
void psdlag_rms::_C( vec p ){
|
||||
int i,j,k; double d;
|
||||
for( i=0 ; i<nfq ; i++ ){ p[ip1[i]] *= mean1sq; p[ip2[i]] *= mean2sq; p[icx[i]] *= mean12;}
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag_rms::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dP1( p , k );
|
||||
}else if( k<2*nfq ){ _dP2( p , k%nfq );
|
||||
}else if( k<3*nfq ){ _dCx( p , k%nfq );
|
||||
}else{ _dPhi( p , k%nfq );}
|
||||
}
|
||||
|
||||
void psdlag_rms::_dP1( vec p , int k ){
|
||||
int i,j;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag_rms::_dP2( vec p , int k ){
|
||||
int i,j;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag_rms::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag_rms::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=p[icx[k]]*mean12,phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag_rms::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
|
||||
pars[i] += dpar[i]/((n<10)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void psdlag_rms::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
if( pars[i+2*nfq]<0 ){ pars[i+2*nfq]*=-1; pars[i+3*nfq] += M_PI;}
|
||||
while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
|
||||
while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
}
|
||||
|
||||
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ //
|
||||
|
||||
psdlag10_rms::psdlag10_rms( lcurve lc1, lcurve lc2 , vec fqL ) {
|
||||
|
||||
// ----------- initial parameters ------------ //
|
||||
n1 = lc1.len;
|
||||
n = n1 + lc2.len;
|
||||
dt = lc1.dt;
|
||||
// ------------------------------------------ //
|
||||
|
||||
|
||||
// ----------- light curve setup ------------ //
|
||||
setlc();
|
||||
mean1 = lc1.demean(); mean2 = lc2.demean();
|
||||
mean1sq = mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
|
||||
int i;
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
lc[i] = lc1.lc[i];
|
||||
lce[i] = lc1.lce[i]*lc1.lce[i];
|
||||
t[i] = lc1.t[i];
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
lc[i] = lc2.lc[i-n1];
|
||||
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
||||
t[i] = lc2.t[i-n1];
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
// ----------- parent initilizer ------------ //
|
||||
init( fqL , 4 );
|
||||
// ------------------------------------------ //
|
||||
|
||||
// --------- constants and indices ---------- //
|
||||
ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
||||
for( i=0 ; i<nfq ; i++ ){
|
||||
ip1[i] = i;
|
||||
ip2[i] = i+nfq;
|
||||
icx[i] = i+2*nfq;
|
||||
iphi[i] = i+3*nfq;
|
||||
}
|
||||
// ------------------------------------------ //
|
||||
|
||||
}
|
||||
|
||||
psdlag10_rms::~psdlag10_rms() {}
|
||||
|
||||
void psdlag10_rms::_C( vec p ){
|
||||
int i,j,k; double d; for(i=0;i<nfq;i++){ p[ip1[i]] = pow(10,p[ip1[i]])*mean1sq; p[ip2[i]] = pow(10,p[ip2[i]])*mean2sq; p[icx[i]] = pow(10,p[icx[i]])*mean12;}
|
||||
for( i=0 ; i<n1 ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
}
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
||||
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
|
||||
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
||||
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10_rms::_dC( vec p , int k ){
|
||||
if( k<nfq ){ _dP1( p , k );
|
||||
}else if( k<2*nfq ){ _dP2( p , k%nfq );
|
||||
}else if( k<3*nfq ){ _dCx( p , k%nfq );
|
||||
}else{ _dPhi( p , k%nfq );}
|
||||
}
|
||||
|
||||
void psdlag10_rms::_dP1( vec p , int k ){
|
||||
int i,j; double p1 = log(10)*pow(10,p[k])*mean1sq;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = p1 * Cfq[k](i,j); dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10_rms::_dP2( vec p , int k ){
|
||||
int i,j; double p2 = log(10) * pow(10,p[ip2[k]])*mean2sq;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = p2 * Cfq[k](i,j); dC(j,i) = dC(i,j);}
|
||||
for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10_rms::_dCx( vec p , int k ){
|
||||
int i,j; double phi=p[iphi[k]], cx=log(10)*pow(10,p[icx[k]])*mean12;
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10_rms::_dPhi( vec p , int k ){
|
||||
int i,j; double cx=pow(10,p[icx[k]])*mean12,phi=p[iphi[k]];
|
||||
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
||||
for( i=n1 ; i<n ; i++ ){
|
||||
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
||||
for( j=0 ; j<n1 ; j++ ){
|
||||
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
||||
}
|
||||
}
|
||||
|
||||
void psdlag10_rms::step_pars( int n , vec& dpar , vec& pars ){
|
||||
for( int i=0 ; i<npar ; i++ ){
|
||||
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
||||
//pars[i] += dpar[i];
|
||||
pars[i] += dpar[i]/((n<5)?10:1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void psdlag10_rms::print_pars( vec& pars , vec& errs ){
|
||||
for( int i=0 ; i<nfq ; i++ ){
|
||||
while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
|
||||
while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
|
||||
}
|
||||
mod::print_pars( pars , errs );
|
||||
}
|
||||
|
||||
void psdlag10_rms::what_pars( int& ip1 , int& ip2 ){
|
||||
ip1 = 3*nfq; ip2 = 4*nfq;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user