Page MenuHomeHEPForge

No OneTemporary

diff --git a/UserSpace/Cpp/Axion/Axion.cpp b/UserSpace/Cpp/Axion/Axion.cpp
index ccfc5be..29bd4cd 100644
--- a/UserSpace/Cpp/Axion/Axion.cpp
+++ b/UserSpace/Cpp/Axion/Axion.cpp
@@ -1,146 +1,146 @@
#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include"src/Axion/AxionSolve.hpp"
// #define printPoints // print the entire evolution of the axion angle, its derivative, and other quantities
// #define printPeaks // print these quantities at the peaks of the oscillation
// #define smallPrint //print the results along with string on names
#define resultPrint //print only the results
#ifndef LONG
#define LONG
#endif
#ifndef LD
#define LD LONG double
#endif
int main(int argc, char **argv){
if(argc!=18){
- std::cerr<<"usage: AxionExample.run theta_i f_a umax TSTOP ratio_ini N_convergence_max convergence_lim inputFile\n \
+ std::cerr<<"usage: Axion.run theta_i f_a umax TSTOP ratio_ini N_convergence_max convergence_lim inputFile\n \
initial_step_size minimum_step_size maximum_step_size absolute_tolerance relative_tolerance beta \n \
fac_max fac_min maximum_No_steps\n";
std::cerr<<"With:\n\n \
theta_i: initial angle\n\n \
fa: PQ scale in GeV (the temperature dependent mass is defined as m_a^2(T) = \\chi(T)/f^2)\n\n \
umax: if u>umax the integration stops (rempember that u=log(a/a_i))\n\n \
TSTOP: if the temperature drops below this, integration stops\n\n \
ratio_ini: integration starts when 3H/m_a<~ratio_ini (this is passed to AxionEOM,\n \
in order to make the interpolations start at this point)\n\n \
N_convergence_max and convergence_lim: integration stops when the relative difference\n \
between two consecutive peaks is less than convergence_lim for N_convergence_max\n \
consecutive peaks\n\n \
inputFile: file that describes the cosmology. the columns should be: u T[GeV] logH\n\n \
-----------Optional arguments------------------------\n\n \
initial_stepsize: initial step the solver takes.\n\n \
maximum_stepsize: This limits the sepsize to an upper limit.\n\n \
minimum_stepsize: This limits the sepsize to a lower limit.\n\n \
absolute_tolerance: absolute tolerance of the RK solver.\n\n \
relative_tolerance: relative tolerance of the RK solver.\n\n \
Note:\n \
Generally, both absolute and relative tolerances should be 1e-8.\n\n \
In some cases, however, one may need more accurate result (eg if f_a is extremely high,\n \
the oscillations happen violently, and the ODE destabilizes). Whatever the case, if the\n \
tolerances are below 1e-8, long doubles *must* be used.\n\n \
beta: controls how agreesive the adaptation is. Generally, it should be around but less than 1.\n\n \
fac_max, fac_min: the stepsize does not increase more than fac_max, and less than fac_min.\n\n \
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must\n \
tweak them in order to avoid instabilities.\n\n \
maximum_No_steps: maximum steps the solver can take Quits if this number is reached even if integration\n \
is not finished.\n";
return 1;
}
int ar=0;
//model parameters
LD theta_i = atof(argv[++ar]) ;
LD fa = atof(argv[++ar]);
// solver parameters
LD umax = atof(argv[++ar]); //u at which the integration stops
LD TSTOP = atof(argv[++ar]); // temperature at which integration stops
LD ratio_ini=atof(argv[++ar]); // 3H/m_a at which integration begins (should be larger than 500 or so)
// stopping conditions.
// integration stops after the adiabatic invariant hasn't changed
// more than convergence_lim% for N_convergence_max consecutive peaks
unsigned int N_convergence_max=atoi(argv[++ar]);
LD convergence_lim=atof(argv[++ar]);
//file in which the cosmology is defined. the columns should be : u T[GeV] logH
std::string inputFile=argv[++ar];
/*options for the solver (These variables are optional. Yoou can use the Axion class without them).*/
LD initial_step_size=atof(argv[++ar]); //initial step the solver takes.
LD minimum_step_size=atof(argv[++ar]); //This limits the sepsize to an upper limit.
LD maximum_step_size=atof(argv[++ar]); //This limits the sepsize to a lower limit.
LD absolute_tolerance=atof(argv[++ar]); //absolute tolerance of the RK solver
LD relative_tolerance=atof(argv[++ar]); //relative tolerance of the RK solver
LD beta=atof(argv[++ar]); //controls how agreesive the adaptation is. Generally, it should be around but less than 1.
/*
the stepsize does not increase more than fac_max, and less than fac_min.
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
tweak them in order to avoid instabilities.
*/
LD fac_max=atof(argv[++ar]);
LD fac_min=atof(argv[++ar]);
unsigned int maximum_No_steps=atoi(argv[++ar]); //maximum steps the solver can take Quits if this number is reached even if integration is not finished.
// instance of Axion
mimes::Axion<LD> Ax(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max,convergence_lim,inputFile,
initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance, relative_tolerance, beta,
fac_max, fac_min, maximum_No_steps);
// Solve the Axion EOM
Ax.solveAxion();
//get the most important results
#ifdef resultPrint
std::cout<<std::setprecision(16)
<<theta_i<<" "<< fa<<" "<<Ax.theta_osc<<" "<<Ax.T_osc<<" "<<Ax.relic<<"\n";
#endif
#ifdef smallPrint
std::cout<<std::setprecision(25)
<<"theta_i="<<theta_i<<" "<<"f_a="<< fa<<" GeV\n"<<"theta_osc~="<<Ax.theta_osc<<" "
<<"T_osc~="<<Ax.T_osc<<"GeV \n"
<<"Omega h^2="<<Ax.relic<<"\n";
#endif
// print all the points
#ifdef printPoints
std::cout<<"---------------------points:---------------------\n";
std::cout<<"a/a_i\tT [GeV]\ttheta\tzeta\trho_a [GeV^4]"<<std::endl;
for(int i=0; i<Ax.pointSize; ++i ){
for(int j=0; j<5; ++j){
std::cout<<Ax.points[i][j];
if(j==4){std::cout<<"\n";}else{std::cout<<"\t";}
}
}
#endif
// print all the peaks
#ifdef printPeaks
std::cout<<"---------------------peaks:---------------------\n";
std::cout<<"a/a_i\tT [GeV]\ttheta\tzeta\trho_a [GeV^4]\tadiabatic_inv [GeV]"<<std::endl;
for(int i=0; i<Ax.peakSize; ++i ){
for(int j=0; j<6; ++j){
std::cout<<Ax.peaks[i][j];
if(j==5){std::cout<<"\n";}else{std::cout<<"\t";}
}
}
#endif
return 0;
}
diff --git a/UserSpace/Cpp/Axion/Axion.sh b/UserSpace/Cpp/Axion/Axion.sh
deleted file mode 100755
index c577573..0000000
--- a/UserSpace/Cpp/Axion/Axion.sh
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/bin/bash
-
-#This script reads a file (which is passed as argument) with the inputs for Axion.run,
-# and prints the result.
-# Also, in stderr, print the time it took (is seconds) to run it.
-
-time (cat $1 | xargs ./Axion.run) 2> .time
-
-sed 's/m/ /g; s/\,/./g; s/s//g' .time | awk '$1 ~ /real/ {print 60*$2+$3}' >&2
-
-
-rm -f .time
-
diff --git a/UserSpace/Cpp/AxionMass/AxionMass.cpp b/UserSpace/Cpp/AxionMass/AxionMass.cpp
index 37c148b..27526e9 100644
--- a/UserSpace/Cpp/AxionMass/AxionMass.cpp
+++ b/UserSpace/Cpp/AxionMass/AxionMass.cpp
@@ -1,74 +1,74 @@
#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include"src/AxionMass/AxionMass.hpp"
#include"src/misc_dir/path.hpp"
//-- Get a number (length) of log_10-spaced points from 10^min to 10^max. --//
template<class LD>
void logspace(LD min, LD max, int length, std::vector<LD> &X ){
for(int i = 0; i<length ; ++i){
X.push_back(pow( 10, min + i*(max-min)/( length-1 )));
}
}
#ifndef LONG
#define LONG
#endif
#ifndef LD
#define LD LONG double
#endif
int main(int argc, char **argv){
if(argc!=4){
- std::cout<<"usage: AxionMassExample.run fa minT maxT\n";
+ std::cout<<"usage: AxionMass.run fa minT maxT\n";
std::cout<<"With:\n \
fa: PQ scale.\n \
minT: minimum interpolation tempareture.\n \
maxT: maximum interpolation tempareture.\n\n \
Beyond these limits, everything is taken constant, so one can use minT=1e-5 GeV and maxT=1e3 with good accuracy.\n";
return 0;
}
LD fa=atof(argv[1]);
LD minT=atof(argv[2]);
LD maxT=atof(argv[3]);
mimes::AxionMass<LD> axM(chi_PATH,minT,maxT);
if(minT==maxT){
std::cout<<axM.ma2(minT,fa)<<"\t"<<axM.ma2(maxT,fa)<<"\n";
return 0;
}
// print N points
unsigned int N=50;
// take logarithmically spaced points
std::vector<LD> T;
logspace<LD>(std::log10(minT), std::log10(maxT),N,T);
std::cout<<"-----Interpolation-----\n";
std::cout<<"T[GeV]\tm_a^2[GeV^2]\tdm_a^2dT[GeV]\n";
for(unsigned int i = 0; i < N; ++i){
std::cout<<T[i]<<"\t"<<axM.ma2(T[i],fa)<<"\t"<<axM.dma2dT(T[i],fa)<<"\n";
}
std::cout<<"-----Approximations-----\n";
std::cout<<"T[GeV]\tm_{a,approx}^2[GeV^2]\tdm_a^2dT_{approx}[GeV]\n";
for(unsigned int i = 0; i < N; ++i){
std::cout<<T[i]<<"\t"<<axM.ma2_approx(T[i],fa)<<"\t"<<axM.dma2dT_approx(T[i],fa)<<"\n";
}
return 0;
}
diff --git a/UserSpace/Cpp/Cosmo/Cosmo.cpp b/UserSpace/Cpp/Cosmo/Cosmo.cpp
index 90c5479..c396ba8 100644
--- a/UserSpace/Cpp/Cosmo/Cosmo.cpp
+++ b/UserSpace/Cpp/Cosmo/Cosmo.cpp
@@ -1,61 +1,61 @@
#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include"src/Cosmo/Cosmo.hpp"
#include"src/misc_dir/path.hpp"
//-- Get a number (length) of log_10-spaced points from 10^min to 10^max. --//
template<class LD>
void logspace(LD min, LD max, int length, std::vector<LD> &X ){
for(int i = 0; i<length ; ++i){
X.push_back(pow( 10, min + i*(max-min)/( length-1 )));
}
}
#ifndef LONG
#define LONG
#endif
#ifndef LD
#define LD LONG double
#endif
int main(int argc, char **argv){
if(argc!=3){
- std::cout<<"usage: CosmoExample.run minT maxT\n";
+ std::cout<<"usage: Cosmo.run minT maxT\n";
std::cout<<"With:\n \
minT: minimum interpolation tempareture.\n \
maxT: maximum interpolation tempareture.\n\n \
Beyond these limits, everything is taken constant, so one can use minT=1e-5 GeV and maxT=3e3 with good accuracy.\n";
return 0;
}
LD minT=atof(argv[1]);
LD maxT=atof(argv[2]);
mimes::Cosmo<LD> cosmo(cosmo_PATH,minT,maxT);
// print N points
unsigned int N=50;
// take logarithmically spaced points
std::vector<LD> T;
logspace<LD>(std::log10(minT), std::log10(maxT),N,T);
std::cout<<"T[GeV]\th_eff\tg_eff\tdh_effdT[GeV^-1]\tdg_effdT[GeV^-1]\tdh\t";
std::cout<<"H[GeV]\ts[GeV^3]\n";
for(unsigned int i = 0; i < N; ++i){
std::cout<<T[i]<<"\t"<<cosmo.heff(T[i])<<"\t"<<cosmo.geff(T[i])
<<"\t"<<cosmo.dheffdT(T[i])<<"\t"<<cosmo.dgeffdT(T[i])<<"\t"
<<"\t"<<cosmo.dh(T[i])<<"\t"<<cosmo.Hubble(T[i])<<"\t"<<cosmo.s(T[i])<<"\n";
}
return 0;
}
diff --git a/UserSpace/scan/AxionScan/AxionScan.py b/UserSpace/scan/AxionScan/AxionScan.py
index 904dc78..7175a6e 100644
--- a/UserSpace/scan/AxionScan/AxionScan.py
+++ b/UserSpace/scan/AxionScan/AxionScan.py
@@ -1,47 +1,47 @@
###-----------------------------------------------------------------------------###
###-------------------Example of interfacePy.ScanScript.Scan-------------------###
###-----------------------------------------------------------------------------###
#------------------------------------Note:------------------------------------#
###########-Scan: scans for the entire theta_i and f_a one provides-###########
#-----------------------------------------------------------------------------#
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../../../src'))
from interfacePy.ScanScript import Scan
from numpy import logspace,linspace,log10,pi
scan=Scan(
cpus=8,
- table_fa= logspace(10,19,100),
- table_theta_i=[1e-2],
+ table_fa= logspace(10,19,50),
+ table_theta_i=linspace(0.1,2,50),
umax=500,
TSTOP=1e-4,
ratio_ini=1e3,
N_convergence_max=5,
convergence_lim=1e-2,
inputFile="../../InputExamples/RDinput.dat",
# inputFile="../../InputExamples/MatterInput.dat",
# inputFile="../../InputExamples/KinationInput.dat",
PathToCppExecutable=r"../../Cpp/Axion/Axion.run",
break_after=60*60*3,
break_time=60,
break_command='',
initial_step_size=1e-2,
minimum_step_size=1e-8,
maximum_step_size=1e-2,
absolute_tolerance=1e-8,
relative_tolerance=1e-8,
beta=0.9,
fac_max=1.2,
fac_min=0.8,
maximum_No_steps=int(1e7)
)
scan.run()
\ No newline at end of file
diff --git a/configure.sh b/configure.sh
index aab2ecb..9f56682 100755
--- a/configure.sh
+++ b/configure.sh
@@ -1,77 +1,80 @@
#!/bin/bash
source Definitions.mk
-for file in $(find . -regex ".*\.sh");do
+# change permissions to all .sh files
+for file in $(find . -type f -regex ".*\.sh");do
echo "allow $file to be executed"
chmod +x $file
done
+# since FormatFile.sh can change directories in future versions, let's find it first
+formatF=$(find . -type f -regex ".*FormatFile.sh")
for dataFile in $cosmoDat $axMDat $anFDat;do
if [ -f "$dataFile" ]; then
echo "format $dataFile"
- bash src/FormatFile.sh $dataFile
+ eval $formatF $dataFile
else
echo "$dataFile does not exist! Make sure that you provide valid paths in Definitions.mk"
exit
fi
done
# ---------these are needed for python and c++---------------- #
echo "make lib directory"
mkdir "lib" 2> /dev/null
echo "make exec directory"
mkdir "exec" 2> /dev/null
echo "make src/misc_dir directory"
mkdir "src/misc_dir" 2> /dev/null
PathHead=src/misc_dir/path.hpp
PathHeadPy=src/misc_dir/path.py
PathTypePy=src/misc_dir/type.py
echo "write type in $PathTypePy"
echo "from ctypes import c_"$LONGpy"double as cdouble" > $PathTypePy
echo "define _PATH_=$PWD in $PathHeadPy"
echo "_PATH_=\"$PWD\" "> $PathHeadPy
echo "write paths in $PathHead"
echo "#ifndef PATHS_HEAD
#define PATHS_HEAD
#define cosmo_PATH \"$PWD/$cosmoDat\"
#define chi_PATH \"$PWD/$axMDat\"
#define anharmonic_PATH \"$PWD/$anFDat\"
#define PWD \"$PWD\"
#endif
">$PathHead
echo "Run \"bash configure.sh License\" in order to read the License, or read the file named LICENSE."
if test "$1" = "License"
then
echo -e "License:"
cat LICENSE
echo ""
echo ""
echo ""
fi
echo -e "\033[1;5;35m
__ __ _ __ __ _____
| \/ | (_) | \/ | / ____|
| \ / | _ | \ / | ___ | (___
| |\/| | | | | |\/| | / _ \ \___ \
| | | | | | | | | | | __/ ____) |
|_| |_| |_| |_| |_| \___| |_____/
"
echo -e "\033[0;97m
You can run \"make\" to compile everything. After that, you will find several examples in UserSpace,
and you can run the executables in exec/ in order to see if the code actually works.
"
diff --git a/src/interfacePy/ScanScript/Scan.py b/src/interfacePy/ScanScript/Scan.py
index e49ef3c..f048a9c 100644
--- a/src/interfacePy/ScanScript/Scan.py
+++ b/src/interfacePy/ScanScript/Scan.py
@@ -1,205 +1,205 @@
# Scan for axion in radiation dominated Uinverse.
#
#------------------Note:------------------#
# Scan scans for the entire theta_i and f_a one provides.
#-----------------------------------------#
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../'))
sysPath.append(osPath.join(osPath.dirname(__file__), '../../src'))
from misc_dir.path import _PATH_
from numpy import array, float64, savetxt
from subprocess import check_output as subprocess_check_output
from sys import stdout as sys_stdout
from time import time
from datetime import datetime
from datetime import timedelta
from os import remove as os_remove
from os import system as os_system
from pathlib import Path
-parallelScan=_PATH_+r"/src/interfacePy/ScanScript/parallel_scan.sh"
+parallelScan=_PATH_+r"/src/util/parallel_scan.sh"
class Scan:
def __init__(self,cpus,table_fa,table_theta_i,umax,TSTOP,ratio_ini,N_convergence_max,convergence_lim,inputFile,
PathToCppExecutable, break_after=0,break_time=60,break_command='',
initial_step_size=1e-2, minimum_step_size=1e-8, maximum_step_size=1e-2,
absolute_tolerance=1e-8, relative_tolerance=1e-8,
beta=0.9, fac_max=1.2, fac_min=0.8, maximum_No_steps=int(1e7)):
'''
scan for different values of fa (in table_fa) table_theta_i. The result file is timecoded
(so it would be difficult to write over it), and the columns correspond to:
theta_i fa [GeV] theta_osc T_osc [GeV] relic (Omega h^2)
Comment: If the scan exits before it finishes, it will continue from the point it stopped.
In order to start from the beginning, delete the file "count.dat".
cpus: number of points to run simultaneously (No of cpus available).
table_fa: table of fa to scan
table_theta_i: table of theta_i to scan
umax: if u>umax the integration stops (rempember that u=log(a/a_i))
TSTOP: if the temperature drops below this, integration stops
ratio_ini: integration starts when 3H/m_a<~ratio_ini (this is passed to AxionEOM,
in order to make the interpolations start at this point)
N_convergence_max and convergence_lim: integration stops when the relative difference
between two consecutive peaks is less than convergence_lim for N_convergence_max
consecutive peaks
inputFile: file that describes the cosmology. the columns should be: u T[GeV] logH
PathToCppExecutable: path to an executable that takes "theta_i fa umax TSTOP ratio_ini N_convergence_max convergence_lim inputFile"
and prints "theta_i fa theta_osc T_osc relic"
break_after,break_time: take a break after break_after seconds for break_time seconds
(optional) break_command: before it takes a break, run this system command (this may be a script to send the results
to an e-mail, or back them up)
-----------Optional arguments------------------------
initial_stepsize: initial step the solver takes.
maximum_stepsize: This limits the sepsize to an upper limit.
minimum_stepsize: This limits the sepsize to a lower limit.
absolute_tolerance: absolute tolerance of the RK solver
relative_tolerance: relative tolerance of the RK solver
Note:
Generally, both absolute and relative tolerances should be 1e-8.
In some cases, however, one may need more accurate result (eg if f_a is extremely high,
the oscillations happen violently, and the ODE destabilizes). Whatever the case, if the
tolerances are below 1e-8, long doubles *must* be used.
beta: controls how agreesive the adaptation is. Generally, it should be around but less than 1.
fac_max, fac_min: the stepsize does not increase more than fac_max, and less than fac_min.
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
tweak them in order to avoid instabilities.
maximum_No_steps: maximum steps the solver can take Quits if this number is reached even if integration
is not finished.
'''
self.cpus=cpus
self.Table_fa=table_fa
self.table_theta_i=table_theta_i
self.umax=umax
self.TSTOP=TSTOP
self.ratio_ini=ratio_ini
self.Npeaks=N_convergence_max
self.conv=convergence_lim
self.inputFile=inputFile
self.PathToCppExecutable=PathToCppExecutable
self.break_after=break_after
self.break_time=break_time
self.break_command=break_command
self.initial_step_size=initial_step_size
self.minimum_step_size=minimum_step_size
self.maximum_step_size= maximum_step_size
self.absolute_tolerance=absolute_tolerance
self.relative_tolerance=relative_tolerance
self.beta=beta
self.fac_max=fac_max
self.fac_min=fac_min
self.maximum_No_steps=maximum_No_steps
self.FileDate=datetime.now()
self.FileName = "{}".format(self.FileDate.strftime('%d-%m-%Y_%H-%M-%S'))
self._p = Path(self.FileName+'.dat')
self.in_file="in.xtx"
self.count_file="count.xtx"
def run_batch(self):
'''
run a batch.
'''
# get the result
points=subprocess_check_output( [parallelScan, self.PathToCppExecutable, str(self.cpus),self.in_file]).decode(sys_stdout.encoding)
points=points.split('\n')
points=array([array(i.split(),float64) for i in points[:-1] ])
points=array(sorted(points, key=lambda arr: arr[0]))
File= self._p.open('ab')
savetxt(File,points)
File.close()
def run(self):
try:
with open(self.count_file,'r') as _:
last_batch=int(_.read())
except:
last_batch=0
Total_batches=len(self.Table_fa)-last_batch
totalT=0
meanT=0
ETA=0
batch=0
sleepTimer=0
for fa in self.Table_fa[last_batch:]:
time0=time()
file=open(self.in_file , 'w')
for theta_i in self.table_theta_i :
file.write( '{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12} {13} {14} {15} {16} \n'.
format(theta_i, fa, self.umax, self.TSTOP, self.ratio_ini, self.Npeaks, self.conv, self.inputFile,
self.initial_step_size, self.minimum_step_size, self.maximum_step_size, self.absolute_tolerance, self.relative_tolerance, self.beta,
self.fac_max, self.fac_min, self.maximum_No_steps) )
file.close()
self.run_batch()
batch+=1
totalT+=time()-time0
meanT=totalT/batch
if self.break_after>0:
sleepTimer+=time()-time0
if sleepTimer>self.break_after:
os_system(self.break_command)
print("taking a break")
time.sleep(self.break_time)
sleepTimer=0
ETA=meanT*(Total_batches-batch)
print('======================================== \n',
'Completed batches: ',batch,' out of ', Total_batches ,'\n',
'Running for: {0:.50s}'.format( str(timedelta(seconds=totalT)) ), '\n',
'ETA: {0:.50s}'.format(str(timedelta(seconds=ETA) ) ),'\n',
'Time per batch: {0:.3f} sec'.format(meanT),
'\n========================================' )
with open(self.count_file,'w') as _:
_.write(str(batch+last_batch))
os_remove(self.in_file)
os_remove(self.count_file)
print('Done!')
\ No newline at end of file
diff --git a/src/interfacePy/ScanScript/ScanObs.py b/src/interfacePy/ScanScript/ScanObs.py
index fcc8258..72e1c2e 100644
--- a/src/interfacePy/ScanScript/ScanObs.py
+++ b/src/interfacePy/ScanScript/ScanObs.py
@@ -1,261 +1,261 @@
# Scan for axion in radiation dominated Uinverse.
#
#------------------Note:------------------#
# The scan proceeds as follows:
# The for each value of $f_a$, we find $\theta_i^{\rm approx}$ such that
# $\Omega h^2 = 0.12$, assuming that $sin(\theta_i^{\rm approx}) \approx \theta_i^{\rm approx}$.
# Then, we scan values of $\theta_i$ close to $\theta_i^{\rm approx}$.
# This results in a range of values of $\theta_i$ for each $f_a$,
# with $\Omega h^2$ close to the observed value.
#-----------------------------------------#
from sys import path as sysPath
from os import path as osPath
sysPath.append(osPath.join(osPath.dirname(__file__), '../'))
sysPath.append(osPath.join(osPath.dirname(__file__), '../../src'))
from misc_dir.path import _PATH_
from interfacePy.Axion import Axion
from interfacePy.Cosmo import h_hub,T0,rho_crit,s
from interfacePy.AxionMass import ma2
from numpy import pi, sqrt, abs, array, float64, argmin, savetxt, linspace
from subprocess import check_output as subprocess_check_output
from sys import stdout as sys_stdout
from time import time
from datetime import datetime
from datetime import timedelta
from os import remove as os_remove
from os import system as os_system
from pathlib import Path
-parallelScan=_PATH_+r"/src/interfacePy/ScanScript/parallel_scan.sh"
+parallelScan=_PATH_+r"/src/util/parallel_scan.sh"
class ScanObs:
def __init__(self,cpus,table_fa,len_theta,umax,TSTOP,ratio_ini,N_convergence_max,convergence_lim,inputFile,
PathToCppExecutable, relic_obs,relic_err_up,relic_err_low,break_after,break_time,break_command='',
initial_step_size=1e-2, minimum_step_size=1e-8, maximum_step_size=1e-2,
absolute_tolerance=1e-8, relative_tolerance=1e-8,
beta=0.9, fac_max=1.2, fac_min=0.8, maximum_No_steps=int(1e7)):
'''
scan for different values of fa (in table_fa) and find the theta_i closer to relic_obs.
The result file is timecoded (so it would be difficult to write over it), and the columns correspond to
theta_i fa [GeV] theta_osc T_osc [GeV] relic (Omega h^2)
Comment 1: The way it works is the following:
1) we solve for theta_i=1e-5.
2) based on this, rescale it in order to find an appropriate theta_i such that Omegah^2 = relic_obs
3) use this rescaled theta_i as initial point, and scan for len_theta between
np.linspace(min([theta_i*0.85,1]),min([theta_i*1.2,pi]),len_theta)
4) finally, the point with Omegah^2 closer to relic_obs is stored.
Comment 2: If the scan exits before it finishes, it will continue from the point it stopped.
In order to start from the beginning, delete the file "count.dat".
cpus: number of points to run simultaneously (No of cpus available).
table_fa: table of fa to scan
len_theta: number of points to search for theta closest to relic_obs. This search happens in
np.linspace(min([theta_i*0.85,1]),min([theta_i*1.2,pi]),self.len_theta), with theta_i is the angle
that results in Omega h^2=relic_obs assuming theta_i<<1.
umax: if u>umax the integration stops (rempember that u=log(a/a_i))
TSTOP: if the temperature drops below this, integration stops
ratio_ini: integration starts when 3H/m_a<~ratio_ini (this is passed to AxionEOM,
in order to make the interpolations start at this point)
N_convergence_max and convergence_lim: integration stops when the relative difference
between two consecutive peaks is less than convergence_lim for N_convergence_max
consecutive peaks
inputFile: file that describes the cosmology. the columns should be: u T[GeV] logH
PathToCppExecutable: path to an executable that takes "theta_i fa umax TSTOP ratio_ini N_convergence_max convergence_lim inputFile"
and prints "theta_i fa theta_osc T_osc relic"
relic_obs: central value of the relic (example, relic_obs=0.12)
relic_err_up: upper error bar around relic_obs (example, relic_err_up=0.001)
relic_err_low: lower error bar around relic_obs (example, relic_err_low=0.001)
break_after,break_time: take a break after break_after seconds for break_time seconds
(optional) break_command: before it takes a break, run this system command (this may be a script to send the results
to an e-mail, or back them up)
-----------Optional arguments------------------------
initial_stepsize: initial step the solver takes.
maximum_stepsize: This limits the sepsize to an upper limit.
minimum_stepsize: This limits the sepsize to a lower limit.
absolute_tolerance: absolute tolerance of the RK solver
relative_tolerance: relative tolerance of the RK solver
Note:
Generally, both absolute and relative tolerances should be 1e-8.
In some cases, however, one may need more accurate result (eg if f_a is extremely high,
the oscillations happen violently, and the ODE destabilizes). Whatever the case, if the
tolerances are below 1e-8, long doubles *must* be used.
beta: controls how agreesive the adaptation is. Generally, it should be around but less than 1.
fac_max, fac_min: the stepsize does not increase more than fac_max, and less than fac_min.
This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must
tweak them in order to avoid instabilities.
maximum_No_steps: maximum steps the solver can take Quits if this number is reached even if integration
is not finished.
'''
self.cpus=cpus
self.Table_fa=table_fa
self.len_theta=len_theta
self.umax=umax
self.TSTOP=TSTOP
self.ratio_ini=ratio_ini
self.Npeaks=N_convergence_max
self.conv=convergence_lim
self.inputFile=inputFile
self.PathToCppExecutable=PathToCppExecutable
self.relic_obs=relic_obs
self.relic_err_up,self.relic_err_low = relic_err_up,relic_err_low
self.break_after=break_after
self.break_time=break_time
self.break_command=break_command
self.initial_step_size=initial_step_size
self.minimum_step_size=minimum_step_size
self.maximum_step_size= maximum_step_size
self.absolute_tolerance=absolute_tolerance
self.relative_tolerance=relative_tolerance
self.beta=beta
self.fac_max=fac_max
self.fac_min=fac_min
self.maximum_No_steps=maximum_No_steps
self.FileDate=datetime.now()
self.FileName = "{}".format(self.FileDate.strftime('%d-%m-%Y_%H-%M-%S'))
self._p = Path(self.FileName+'.dat')
self.in_file="in.xtx"
self.count_file="count.xtx"
def run_batch(self):
'''
run a batch.
'''
# get the result
points=subprocess_check_output( [parallelScan, self.PathToCppExecutable, str(self.cpus),self.in_file]).decode(sys_stdout.encoding)
points=points.split('\n')
points=array([array(i.split(),float64) for i in points[:-1] ])
points=array(sorted(points, key=lambda arr: arr[0]))
absDiff=argmin(abs(points[:,4]-self.relic_obs))
_=points[absDiff]
if _[4]<=self.relic_obs+self.relic_err_up and _[4]>=self.relic_obs-self.relic_err_low:
File= self._p.open('ab')
savetxt(File,array([_]))
File.close()
def run(self):
try:
with open(self.count_file,'r') as _:
last_batch=int(_.read())
except:
last_batch=0
Total_batches=len(self.Table_fa)-last_batch
totalT=0
meanT=0
ETA=0
batch=0
sleepTimer=0
for fa in self.Table_fa[last_batch:]:
time0=time()
########################---find theta_i (assuming theta_i<<1) such that Omega h^=0.12---########################
theta_small=1e-3
ax=Axion(theta_small,fa,self.umax,self.TSTOP,self.ratio_ini,self.Npeaks,self.conv,self.inputFile)
ax.solve()
ax.getPeaks()
T=ax.T_peak[-1]
theta=ax.theta_peak[-1]
relic=ax.relic
theta_obs=sqrt(theta**2/relic*self.relic_obs)
relic=s(T0)/s(T)*0.5*sqrt(ma2(0,1)*ma2(T,1))*theta_obs**2*h_hub**2/rho_crit
theta_small_i=theta_small*theta_obs/theta
if theta_small_i < theta_small:
Table_theta_i = array([theta_small_i])
else:
if theta_small_i>pi:
theta_small_i=pi
Table_theta_i=linspace(min([theta_small_i*0.85,1]),min([theta_small_i*1.2,pi]),self.len_theta)
del ax
########################---END---########################
file=open(self.in_file , 'w')
for theta_i in Table_theta_i :
file.write( '{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11} {12} {13} {14} {15} {16} \n'.
format(theta_i, fa, self.umax, self.TSTOP, self.ratio_ini, self.Npeaks, self.conv, self.inputFile,
self.initial_step_size, self.minimum_step_size, self.maximum_step_size, self.absolute_tolerance, self.relative_tolerance, self.beta,
self.fac_max, self.fac_min, self.maximum_No_steps) )
file.close()
self.run_batch()
batch+=1
totalT+=time()-time0
meanT=totalT/batch
if self.break_after>0:
sleepTimer+=time()-time0
if sleepTimer>self.break_after:
os_system(self.break_command)
print("taking a break")
time.sleep(self.break_time)
sleepTimer=0
ETA=meanT*(Total_batches-batch)
print('======================================== \n',
'Completed batches: ',batch,' out of ', Total_batches ,'\n',
'Running for: {0:.50s}'.format( str(timedelta(seconds=totalT)) ), '\n',
'ETA: {0:.50s}'.format(str(timedelta(seconds=ETA) ) ),'\n',
'Time per batch: {0:.3f} sec'.format(meanT),
'\n========================================' )
with open(self.count_file,'w') as _:
_.write(str(batch+last_batch))
os_remove(self.in_file)
os_remove(self.count_file)
print('Done!')
\ No newline at end of file
diff --git a/src/FormatFile.sh b/src/util/FormatFile.sh
similarity index 100%
rename from src/FormatFile.sh
rename to src/util/FormatFile.sh
diff --git a/src/interfacePy/ScanScript/parallel_scan.sh b/src/util/parallel_scan.sh
similarity index 100%
rename from src/interfacePy/ScanScript/parallel_scan.sh
rename to src/util/parallel_scan.sh
diff --git a/src/util/timeit.sh b/src/util/timeit.sh
new file mode 100755
index 0000000..6930f14
--- /dev/null
+++ b/src/util/timeit.sh
@@ -0,0 +1,42 @@
+#!/bin/bash
+
+# This script executes the first argument and prints the time it took (is seconds) to run it.
+#
+
+# check if the number of arguments is correct
+expectedArgs=2
+if [ $# -ne $expectedArgs ]; then
+ echo "Please provide $expectedArgs arguments."
+ echo "The first should be a path to an executable, the second should be a file that contains its arguments (each argument should be in different line)."
+ exit
+fi
+
+# check if the first argument can be executed
+if ! command -v $1 &>/dev/null ; then
+ echo "$1 is not a valid executable."
+ echo "Please provide a valid path to an executable."
+ exit
+fi
+
+
+# check if the file exists
+if [ ! -f $2 ]; then
+ echo "File $2 does not exist."
+ echo "Please provide a valid path to a file."
+ exit
+fi
+
+# check if the file is readable
+if [ ! -r $2 ]; then
+ echo "File $2 is not readable."
+ echo "Please provide a valid path to a readable file."
+ exit
+fi
+
+# run the executable ($1) passing as arguments the contents of the file ($2)
+time (xargs -a $2 $1) 2> .time.mimes
+
+# print the time it took in seconds
+perl -pe 's/m/ /g; s/\,/./g; s/s//g' .time.mimes | awk '$1 ~ /real/ {print 60*$2+$3}' >&2
+rm -f .time.mimes
+

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 11:07 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242807
Default Alt Text
(36 KB)

Event Timeline