Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add FDJUMP parameters to the designmmatrix #64

Merged
merged 4 commits into from
Sep 5, 2024
Merged

add FDJUMP parameters to the designmmatrix #64

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
725933d
add FDJUMP parameters to the designmmatrix
siyuan-chen Jun 19, 2024
fa29b2c
Merge remote-tracking branch 'vallis/master' into fdjump
siyuan-chen Jul 24, 2024
5923908
separate FDJUMP and FDJUMPDM in parameter names
siyuan-chen Sep 3, 2024
6a6a880
merge vallis/master
siyuan-chen Sep 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 70 additions & 8 deletions libstempo/libstempo.pyx
View file
Open in desktop
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]
mattpitkin marked this conversation as resolved.
Show resolved Hide resolved
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]
mattpitkin marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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);
Loading