Skip to content

Commit

Permalink
Merge pull request #64 from siyuan-chen/fdjump
Browse files Browse the repository at this point in the history 
add FDJUMP parameters to the designmmatrix
  • Loading branch information
@mattpitkin
mattpitkin authored Sep 5, 2024
2 parents 60c76f3 + 6a6a880 commit a6098e3
Show file tree
Hide file tree
Showing 2 changed files with 73 additions and 8 deletions.
78 changes: 70 additions & 8 deletions libstempo/libstempo.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,7 @@ cdef extern from "tempo2.h":
enum: param_LAST
enum: param_ZERO
enum: param_JUMP
enum: param_FDJUMP
enum: MAX_T2EFAC
enum: MAX_T2EQUAD
enum: MAX_TNEF
Expand Down Expand Up @@ -211,7 +212,9 @@ cdef extern from "tempo2.h":
int noWarnings
double fitChisq
int nJumps
char fjumpID[16]
double jumpVal[MAX_JUMPS]
# char jumpSAT[MAX_JUMPS]
int fitJump[MAX_JUMPS]
double jumpValErr[MAX_JUMPS]
char *binaryModel
Expand All @@ -236,6 +239,15 @@ cdef extern from "tempo2.h":
double ne_sw_ifuncE[MAX_IFUNC]
int ne_sw_ifuncN

# new parameters for fdjumps
int nfdJumps
char ffdjumpID[16]
double fdjumpVal[MAX_JUMPS]
int fdjumpIdx[MAX_JUMPS]
int fitfdJump[MAX_JUMPS]
double fdjumpValErr[MAX_JUMPS]
char fdjump_log

# noise parameters follow

# T2EFAC
Expand Down Expand Up @@ -325,6 +337,9 @@ cdef extern from "t2fit-stub.h":

double t2FitFunc_zero(pulsar *psr,int ipsr,double x,int ipos,param_label label,int k)
void t2UpdateFunc_zero(pulsar *psr,int ipsr,param_label label,int k,double val,double err)

double t2FitFunc_fdjump(pulsar *psr,int ipsr,double x,int ipos,param_label label,int k)
void t2UpdateFunc_fdjump(pulsar *psr,int ipsr,param_label label,int k,double val,double err)

void t2fit_fillOneParameterFitInfo(pulsar* psr,param_label fit_param,const int k,FitInfo& OUT)

Expand Down Expand Up @@ -354,6 +369,7 @@ cdef class tempopar:
cdef public int subct

cdef int _isjump
cdef int _isfdjump
cdef void *_val
cdef void *_err
cdef int *_fitFlag
Expand All @@ -370,7 +386,7 @@ cdef class tempopar:

property val:
def __get__(self):
if not self._isjump:
if not self._isjump and not self._isfdjump:
return self._unitify(get_longdouble_as_scalar((<long double*>self._val)[0]))
else:
return self._unitify(float((<double*>self._val)[0]))
Expand All @@ -384,7 +400,7 @@ cdef class tempopar:
elif self.unit and isinstance(value,Quantity):
value = value.to(self.unit).value

if not self._isjump:
if not self._isjump and not self._isfdjump:
if not self._paramSet[0]:
self._paramSet[0] = 1

Expand All @@ -394,7 +410,7 @@ cdef class tempopar:

property err:
def __get__(self):
if not self._isjump:
if not self._isjump and not self._isfdjump:
return self._unitify(get_longdouble_as_scalar((<long double*>self._err)[0]))
else:
return self._unitify(float((<double*>self._err)[0]))
Expand All @@ -403,7 +419,7 @@ cdef class tempopar:
if self.unit and isinstance(value,Quantity):
value = value.to(self.unit).value

if not self._isjump:
if not self._isjump and not self._isfdjump:
set_longdouble(<long double*>self._err,value)
else:
(<double*>self._err)[0] = value
Expand All @@ -414,8 +430,9 @@ cdef class tempopar:

def __set__(self,value):
if value:
if not self._isjump and not self._paramSet[0]:
self._paramSet[0] = 1
if not self._isjump and not self._isfdjump:
if not self._paramSet[0]:
self._paramSet[0] = 1

self._fitFlag[0] = 1
else:
Expand All @@ -424,13 +441,13 @@ cdef class tempopar:
# note that paramSet is not always respected in tempo2
property set:
def __get__(self):
if not self._isjump:
if not self._isjump and not self._isfdjump:
return True if self._paramSet[0] else False
else:
return True

def __set__(self,value):
if not self._isjump:
if not self._isjump and not self._isfdjump:
if value:
self._paramSet[0] = 1
else:
Expand All @@ -442,6 +459,10 @@ cdef class tempopar:
def __get__(self):
return True if self._isjump else False

property isfdjump:
def __get__(self):
return True if self._isfdjump else False

def __str__(self):
if self.set:
return 'tempo2 parameter %s (%s): %s +/- %s' % (self.name,'fitted' if self.fit else 'not fitted',repr(self.val),repr(self.err))
Expand Down Expand Up @@ -473,6 +494,7 @@ cdef create_tempopar(parameter par,int ct,int subct,int eclCoord,object units):
newpar.timescale = None

newpar._isjump = 0
newpar._isfdjump = 0

newpar._val = &par.val[subct]
newpar._err = &par.err[subct]
Expand All @@ -499,6 +521,7 @@ cdef create_tempojump(pulsar *psr,int ct,object units):
newpar.name = 'JUMP{0}'.format(ct)

newpar._isjump = 1
newpar._isfdjump = 0

newpar._val = &psr.jumpVal[ct]
newpar._err = &psr.jumpValErr[ct]
Expand All @@ -509,6 +532,33 @@ cdef create_tempojump(pulsar *psr,int ct,object units):

return newpar

cdef create_tempofdjump(pulsar *psr,int ct,int fddmct,object units):
cdef tempopar newpar = tempopar.__new__(tempopar)

# TO DO: proper units
if units:
newpar.unit = u.dimensionless_unscaled
newpar.timescale = None
else:
newpar.unit = None
newpar.timescale = None

if psr.fdjumpIdx[ct] == -2:
newpar.name = 'FDJUMPDM{0}'.format(fddmct)
else:
newpar.name = 'FDJUMP{0}'.format(ct)

newpar._isjump = 0
newpar._isfdjump = 1

newpar._val = &psr.fdjumpVal[ct]
newpar._err = &psr.fdjumpValErr[ct]
newpar._fitFlag = &psr.fitfdJump[ct]

newpar.ct = param_FDJUMP
newpar.subct = ct

return newpar

# TODO: check if consistent with new API
cdef class GWB:
Expand Down Expand Up @@ -850,6 +900,12 @@ cdef class tempopulsar:
for ct in range(1,self.psr[0].nJumps+1): # jump 1 in the array not used...
newpar = create_tempojump(&self.psr[0],ct,self.units)
self.pardict[newpar.name] = newpar

for ct in range(1,self.psr[0].nfdJumps+1): # jump 1 in the array not used...
if self.psr[0].fdjumpIdx[ct] == -2:
fddmct += 1
newpar = create_tempofdjump(&self.psr[0],ct,fddmct,self.units)
self.pardict[newpar.name] = newpar

# the designmatrix plugin also adds extra parameters for sinusoidal whitening
# but they don't seem to be used in the EPTA analysis
Expand Down Expand Up @@ -1713,6 +1769,12 @@ cdef class tempopulsar:
fitinfo.paramDerivs[fitinfo.nParams] = t2FitFunc_jump
fitinfo.updateFunctions[fitinfo.nParams] = t2UpdateFunc_jump
fitinfo.nParams = fitinfo.nParams + 1
elif self[par].isfdjump:
fitinfo.paramIndex[fitinfo.nParams] = param_FDJUMP
fitinfo.paramCounters[fitinfo.nParams] = self[par].subct
fitinfo.paramDerivs[fitinfo.nParams] = t2FitFunc_fdjump
fitinfo.updateFunctions[fitinfo.nParams] = t2UpdateFunc_fdjump
fitinfo.nParams = fitinfo.nParams + 1
else:
t2fit_fillOneParameterFitInfo(&self.psr[0],self[par].ct,self[par].subct,fitinfo)
# the function already increases nParams
Expand Down
3 changes: 3 additions & 0 deletions libstempo/t2fit-stub.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,9 @@ void t2UpdateFunc_zero(struct pulsar *psr,int ipsr,param_label label,int k,doubl
double t2FitFunc_jump(struct pulsar *psr,int ipsr,double x,int ipos,param_label label,int k);
void t2UpdateFunc_jump(struct pulsar *psr,int ipsr,param_label label,int k,double val,double err);

double t2FitFunc_fdjump(struct pulsar *psr,int ipsr,double x,int ipos,param_label label,int k);
void t2UpdateFunc_fdjump(struct pulsar *psr,int ipsr,param_label label,int k,double val,double err);

// defined in t2fit.C, not declared in t2fit.h

void t2fit_fillOneParameterFitInfo(struct pulsar *psr,param_label fit_param,const int k,FitInfo& OUT);

0 comments on commit a6098e3

Please sign in to comment.