Commit 2f6e224c authored by Miroslav Shaltev's avatar Miroslav Shaltev
Browse files

in sa_mc_nomad.py use minimization to draw offsignal parameters for dim >= 5

parent dfc5f039
......@@ -19,6 +19,7 @@ from time import time as clock
from numpy import *
from scipy import *
from scipy.misc import *
from scipy.optimize import *
import random
import string
import fileinput
......@@ -243,6 +244,163 @@ def point_inside_sum(g,d,s,c):
def n3equ(Alpha,Delta):
return {'x':cos(Alpha) * cos(Delta),'y':sin(Alpha) * cos(Delta),'z':sin(Delta)}
def ad_to_n3equ(p):
r = []
tn3 = n3equ(p[0],p[1])
r.append(tn3['x'])
r.append(tn3['y'])
r.append(tn3['z'])
for k in range(2,len(p)):
r.append(p[k])
return r
def min_point_inside_sum(tc,g,d,ts):
c = ad_to_n3equ(tc)
s = ad_to_n3equ(ts)
sumg = 0
for i in range(d):
for j in range(d):
sumg = sumg + g[i][j]*(c[i]-s[i])*(c[j]-s[j])
return sumg
def get_offset_signal(signal,semig,directedsearch,c,bands):
offsignal = {}
point_not_inside = True
s = []
slen = c['spindownorder'] + 1
skybndoff = 0
if not directedsearch:
n3 = n3equ(signal['Alpha'],signal['Delta'])
s.append(n3['x'])
s.append(n3['y'])
s.append(n3['z'])
slen += 3
skybndoff = 2
s.append(signal['Freq'])
s.append(signal['f1dot'])
s.append(signal['f2dot'])
s.append(signal['f3dot'])
while point_not_inside:
tc = []
tmp_alpha = signal['Alpha']
tmp_delta = signal['Delta']
offcnt = 0
if not directedsearch:
tmp_alpha = random.uniform(max(0,signal['Alpha']-bands[0]),min(2.*pi,signal['Alpha']+bands[0]))
tmp_delta = random.uniform(max(-0.5*pi,signal['Delta']-bands[1]),min(0.5*pi,signal['Delta']+bands[1]))
tn3 = n3equ(tmp_alpha,tmp_delta)
tc.append(tn3['x'])
tc.append(tn3['y'])
tc.append(tn3['z'])
offcnt = 3
#print 'slen: ',slen
for i in range(slen-offcnt):
tc.append(random.uniform(s[i+offcnt]-bands[skybndoff+i],s[i+offcnt]+bands[skybndoff+i]))
m = point_inside_sum(semig,len(semig[0]),s,tc)
#print 'this m: ',m,tc
if ( c['mmax'] - m ) > 0 and ( m - c['mmin'] ) > 0:
for i in range(slen-offcnt,8):
tc.append(0)
print 'mmin: %f, m: %f, mmax: %f'%(c['mmin'],m,c['mmax'])
offsignal['Alpha'] = tmp_alpha
offsignal['Delta'] = tmp_delta
offsignal['Freq'] = tc[offcnt + 0]
offsignal['f1dot'] = tc[offcnt + 1]
offsignal['f2dot'] = tc[offcnt + 2]
offsignal['f3dot'] = tc[offcnt + 3]
point_not_inside = False
return offsignal
def get_hdim_offset_signal(signal,semig,directedsearch,c,bands):
offsignal = {}
point_not_inside = True
s = []
slen = c['spindownorder'] + 1
skybndoff = 0
if not directedsearch:
s.append(signal['Alpha'])
s.append(signal['Delta'])
slen += 2
skybndoff = 2
s.append(signal['Freq'])
s.append(signal['f1dot'])
s.append(signal['f2dot'])
s.append(signal['f3dot'])
min_alpha = max(0,signal['Alpha']-bands[0])
max_alpha = min(2.*pi,signal['Alpha']+bands[0])
min_delta = max(-0.5*pi,signal['Delta']-bands[1])
max_delta = min(0.5*pi,signal['Delta']+bands[1])
tc = []
tmp_alpha = signal['Alpha']
tmp_delta = signal['Delta']
offcnt = 0
if not directedsearch:
tmp_alpha = random.uniform(min_alpha,max_alpha)
tmp_delta = random.uniform(min_delta,max_delta)
tc.append(tmp_alpha)
tc.append(tmp_delta)
offcnt = 2
#print 'slen: ',slen
for i in range(slen-offcnt):
tc.append(random.uniform(s[i+offcnt]-bands[skybndoff+i],s[i+offcnt]+bands[skybndoff+i]))
try_to_minimize = True
res_x = []
tryiter = 0
while try_to_minimize:
res = fmin(min_point_inside_sum, tc, args = (semig,len(semig[0]),s))
check_m = point_inside_sum(semig,len(semig[0]),ad_to_n3equ(s),ad_to_n3equ(res))
#print res
if ( c['mmax'] - check_m ) > 0 and ( check_m - c['mmin'] ) > 0:
res_x = res
try_to_minimize = False
else:
tryiter += 1
#print 'try iter: ',tryiter
if tryiter > 10:
tc = []
tmp_alpha = signal['Alpha']
tmp_delta = signal['Delta']
offcnt = 0
if not directedsearch:
tmp_alpha = random.uniform(min_alpha,max_alpha)
tmp_delta = random.uniform(min_delta,max_delta)
tc.append(tmp_alpha)
tc.append(tmp_delta)
offcnt = 2
#print 'slen: ',slen
for i in range(slen-offcnt):
tc.append(random.uniform(s[i+offcnt]-bands[skybndoff+i],s[i+offcnt]+bands[skybndoff+i]))
#print res
check_m = point_inside_sum(semig,len(semig[0]),ad_to_n3equ(s),ad_to_n3equ(res_x))
print res_x
print 'check m: ',check_m
offsignal['Alpha'] = res_x[0]
offsignal['Delta'] = res_x[1]
offsignal['Freq'] = res_x[2]
offsignal['f1dot'] = res_x[3]
offsignal['f2dot'] = res_x[4]
if len(res_x) > 5:
offsignal['f3dot'] = res_x[5]
else:
offsignal['f3dot'] = 0
return offsignal
def create_parameter_file(p):
f = open(p['filename'],'w')
f.write(p['parameters'])
......@@ -1339,54 +1497,13 @@ def main():
print 'bands: ',bands
# draw offsignal parameters
offsignal = {}
point_not_inside = True
s = []
slen = c['spindownorder'] + 1
skybndoff = 0
if not directedsearch:
n3 = n3equ(signal['Alpha'],signal['Delta'])
s.append(n3['x'])
s.append(n3['y'])
s.append(n3['z'])
slen += 3
skybndoff = 2
s.append(signal['Freq'])
s.append(signal['f1dot'])
s.append(signal['f2dot'])
s.append(signal['f3dot'])
while point_not_inside:
tc = []
tmp_alpha = signal['Alpha']
tmp_delta = signal['Delta']
offcnt = 0
if not directedsearch:
tmp_alpha = random.uniform(max(0,signal['Alpha']-bands[0]),min(2.*pi,signal['Alpha']+bands[0]))
tmp_delta = random.uniform(max(-0.5*pi,signal['Delta']-bands[1]),min(0.5*pi,signal['Delta']+bands[1]))
tn3 = n3equ(tmp_alpha,tmp_delta)
tc.append(tn3['x'])
tc.append(tn3['y'])
tc.append(tn3['z'])
offcnt = 3
#print 'slen: ',slen
for i in range(slen-offcnt):
tc.append(random.uniform(s[i+offcnt]-bands[skybndoff+i],s[i+offcnt]+bands[skybndoff+i]))
m = point_inside_sum(semig,len(semig[0]),s,tc)
#print 'this m: ',m,tc
if ( c['mmax'] - m ) > 0 and ( m - c['mmin'] ) > 0:
for i in range(slen-offcnt,8):
tc.append(0)
print 'mmin: %f, m: %f, mmax: %f'%(c['mmin'],m,c['mmax'])
offsignal['Alpha'] = tmp_alpha
offsignal['Delta'] = tmp_delta
offsignal['Freq'] = tc[offcnt + 0]
offsignal['f1dot'] = tc[offcnt + 1]
offsignal['f2dot'] = tc[offcnt + 2]
offsignal['f3dot'] = tc[offcnt + 3]
point_not_inside = False
if not directedsearch and c['spindownorder'] > 1:
offsignal = get_hdim_offset_signal(signal,semig,directedsearch,c,bands)
else:
offsignal = get_offset_signal(signal,semig,directedsearch,c,bands)
print 'offsignal: ',offsignal
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment