Skip to content
Snippets Groups Projects
Select Git revision
  • 1a7347c926d4c7c29823ffd07510bb7bd431b54e
  • master default
  • mingw_gcc44
  • release_ABP1_012
  • release_ABP1_008
  • release_ABP1_006
  • release_ABP1_007
  • release_ABP1_005
  • release_ABP1_004
  • release_ABP1_003
  • pre_release_0.15
  • release_ABP1_001
  • release_ABP1_002
  • pre_release_0.13
  • pre_release_0.14
  • pre_release_0.11
  • pre_release_0.12
  • pre_release_0.10
  • pre_release_0.09
  • pre_release_0.08
20 results

doxygen_index.h

Blame
  • Forked from einsteinathome / graphicsframework
    Source project has a limited visibility.
    measureComb.py 11.50 KiB
    #!/usr/bin/python
    import math, os, re
    import matplotlib as matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    import matplotlib.cm as cm
    import numpy as np
    
    def combSpectrum(targetDirectory, flag):
    
        # Grant David Meadors
        # g m e a d o r s @  u m i c h . e d u
        # 02012-11-28 (JD 2456260)
        # combSpectrum
        #
        # Scan output log files for parts of the spectrum where feedforward is 
        # adversely affecting Hoft.
        #
        # Input argument: targetDirectory
        # Choose a directory, e.g. "/home/pulsar/public_html/feedforward/diagnostics/LHO/H-H1_AMPS_C02_L2-9531",
        # into which comb files have been written.
        # Input argument: flag
        # Flag whether the directory itself contains comb files ('one') o
        # is a directory of directories containing comb files ('all')
    
        # The results of the frequency comb ratio test are stored in files of the form
        # EleutheriaComb-startGPS-stopGPS.txt
        # and thus can be found be searching for (substitute LLO if appropriate)
        # /home/pulsar/public_html/feedforward/diagnostics/LHO/*/*Comb*
        # In each file, the comb frequencies are listed as the values of the columns in
        # row 5 (using Python indexing, starting from 0), and the ratios as the values
        # in row 6; the differences are in row 7.
    
        if flag == 'all':
        # The comb files are in the directories listed under the target directory.
            highDirectoryList = os.listdir(targetDirectory)
            highDirectoryListDirOnly = sorted([x for x in highDirectoryList if x.find('.') == -1])
            # Because the GPS time 953100000 is used so often for testing,
            # it is temporarily excluded:
            highDirectoryListDirOnly = sorted([x for x in highDirectoryListDirOnly if x.find('9531') == -1])
            combFiles = []
            for x in highDirectoryListDirOnly:
                targetDirectoryFiles = os.listdir(targetDirectory + x)
                for file in targetDirectoryFiles:
                    if file.find('EleutheriaComb') > -1:
                        combFiles.append(x + '/' + file)
            combFiles = sorted(combFiles)
        if flag == 'one':
            # The comb files are in the target directory.
            # Pull a list of all the comb files.
            targetDirectoryFiles = sorted(os.listdir(targetDirectory))
            combFiles = []
            for file in targetDirectoryFiles:
                if file.find('EleutheriaComb') > -1:
                    combFiles.append(file)
            combFiles = sorted(combFiles)
    
        # Now write a function that can read an individual file into memory.
        def combReader(targetDirectory, eachCombFile):
            combFileName = targetDirectory + eachCombFile
            try:
                combObject = open(combFileName)
                combLines = combObject.readlines()
                combObject.close()
            except IOError:
                print 'File not found or accessible; skipping ' +\
                combFileName
            return combLines[5:8] + combLines[11:13]
        # Apply the function to all the comb files in the directory.
        for k, eachCombFile in enumerate(combFiles):
            # Each file will have four properties:
            # (using zero-indexed terminology)
            # start time, read in the GPS time of the file name
            # frequency bins, read from the fifth row of the file (0th of combRaw),
            # ratios, read from the sixth (1st of combRaw),
            # and difference, read from the seventh (2nd of combRaw).
            # Then, after some comment lines,
            # before feedforward in the eleventh line (3rd of combRaw)
            # and after feedforward in the twelfth line (4th of combRaw) 
            if k == 0:
                times = re.search('-(?P<GPS>\d+)-(\d+)\.', eachCombFile)
                timeArray = np.asarray(times.group(1), dtype=np.int32)
                combRaw = combReader(targetDirectory, eachCombFile)
                frequencyArray = np.asarray(str(combRaw[0]).split(), dtype=np.float32)
                frequencyArray = frequencyArray.astype(int)
                ratioArray = np.asarray(str(combRaw[1]).split(), dtype=np.float32)
                differenceArray = np.asarray(str(combRaw[2]).split(), dtype=np.float32)
                beforeArray = np.asarray(str(combRaw[3]).split(), dtype=np.float32)
                afterArray = np.asarray(str(combRaw[4]).split(), dtype=np.float32)
            if k > 0:
                times = re.search('-(?P<GPS>\d+)-(\d+)\.', eachCombFile)
                timeArray = np.vstack([timeArray, np.asarray(times.group(1), dtype=np.int32)])
                combRaw = combReader(targetDirectory, eachCombFile)
                frequencyArray = np.vstack([frequencyArray, np.asarray(str(combRaw[0]).split(), dtype=np.float32)])
                ratioArray = np.vstack([ratioArray, np.asarray(str(combRaw[1]).split(), dtype=np.float32)])
                differenceArray = np.vstack([differenceArray, np.asarray(str(combRaw[2]).split(), dtype=np.float32)])
                beforeArray = np.vstack([beforeArray, np.asarray(str(combRaw[3]).split(), dtype=np.float32)])
                afterArray = np.vstack([afterArray, np.asarray(str(combRaw[4]).split(), dtype=np.float32)])
    
        # Now we are going to plot these arrays.
        # First, set a directory for the output. At least for the time being,
        # that can be the same directory as the target directory.
        # Note the inelegant way from extracting the start time as a label.
        graphTitleRatio = targetDirectory + "EleutheriaPostPlotCombRatio" +\
        '-' + str(timeArray[0]).strip("[").strip("]").strip("'")
        graphTitleDiff = targetDirectory + "EleutheriaPostPlotCombDiff" +\
        '-' + str(timeArray[0]).strip("[").strip("]").strip("'")
        graphTitleBeforeAfter = targetDirectory + "EleutheriaPostPlotCombBeforeAfter" +\
        '-' + str(timeArray[0]).strip("[").strip("]").strip("'")
        graphWholeRun = "EleutheriaWholeRun"
        x, y = np.meshgrid(timeArray, frequencyArray[0])
        extensions = [x[0, 0], x[-1, -1], y[0, 0], y[-1, -1]]
        # Define a helpfully broken color-scheme as found on the 
        # SciPy Matplotlib Cookbook,
        # http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps
        #cdict = {'red': ((0.0, 0.0, 0.0),\
        #(0.5, 1.0, 0.7),\
        #(1.0, 1.0, 1.0)),\
        #'green': ((0.0, 0.0, 0.0),\
        #(0.5, 1.0, 0.0),\
        #(1.0, 1.0, 1.0)),\
        #'blue': ((0.0, 0.0, 0.0),\
        #(0.5, 1.0, 0.0),\
        #(1.0, 0.5, 1.0))}
        #my_cmap = matplotlib.colors.LinearSegmentedColormap('my_colormap', cdict, 256)
        #plt.figure()
        #CSratio = plt.contourf(x.T, y.T, ratioArray, \
        #np.asarray([0.8, 0.85, 0.9, 0.95, \
        #1, 1.05, 1.1, 1.15, 1.2]), cmap=my_cmap)
        #plt.colorbar(CSratio, shrink=0.8, extend='both')
        #plt.xlabel('GPS time (s)')
        #plt.ylabel('Frequency (Hz)')
        #plt.title('Post/pre-filtering Hoft ratio (lower is better)')
        #plt.savefig(graphTitleRatio + 'Contour.png')
        #plt.savefig(graphTitleRatio + 'Contour.pdf')
        #plt.close()
        #plt.clf()
        #fig = plt.figure()
        #ax = fig.add_subplot(111)
        #CSratioIm = ax.imshow(ratioArray.T, origin='lower', \
        #interpolation = 'nearest', extent=extensions, vmin=0.8, vmax=1.2,\
        #cmap=my_cmap) 
        #CSratioIm = fig.colorbar(CSratioIm, shrink = 0.8, extend = 'both')
        #ax.set_aspect('auto')
        #ax.set_xlabel('GPS time (s)')
        #ax.set_ylabel('Frequency (Hz)')
        #ax.set_title('Post/pre-filtering Hoft ratio (lower is better)')
        #fig.savefig(graphTitleRatio + '.png')
        #fig.savefig(graphTitleRatio + '.pdf')
        #plt.clf()
        #fig = plt.figure()
        #ax = fig.add_subplot(111)
        #CSdiffIm = ax.imshow(-differenceArray.T, origin='lower', \
        #interpolation = 'nearest' , extent=extensions, vmin=-2e-24, vmax=2e-24, \
        #cmap=my_cmap)
        #plt.show()
        #CSdiffIm = fig.colorbar(CSdiffIm, shrink = 0.8, extend = 'both')
        #ax.set_aspect('auto')
        #ax.set_xlabel('GPS time (s)')
        #ax.set_ylabel('Frequency (Hz)')
        #ax.set_title('Post - pre Hoft difference (lower is better)')
        #fig.savefig(graphTitleDiff + '.png')
        #fig.savefig(graphTitleDiff + '.pdf')
        #plt.clf()
        #plt.figure()
        #CSdiff = plt.contourf(x.T, y.T, -differenceArray, \
        #np.asarray([-2e-24, -1.5e-24, -1e-24, -5e-25,\
        #0, 5e-25, 1e-24, 1.5e-24, 2e-24]), cmap=my_cmap)
        #plt.colorbar(CSdiff, shrink=0.8, extend='both')
        #plt.xlabel('GPS time (s)')
        #plt.ylabel('Frequency (Hz)')
        #plt.title('Post - pre Hoft difference (lower is better)')
        #plt.savefig(graphTitleDiff + 'Contour.png')
        #plt.savefig(graphTitleDiff + 'Contour.pdf')
        #plt.close()
        #plt.figure()
        #freqList = [25, 27, 29]
        #p0 = plt.scatter(timeArray, -differenceArray[:, 25], color='b')
        #mean0 = np.mean(-differenceArray[:, 25])
        #p1 = plt.scatter(timeArray, -differenceArray[:, 27], color='r')
        #mean1 = np.mean(-differenceArray[:, 27])
        #p2 = plt.scatter(timeArray, -differenceArray[:, 29], color='k')
        #mean2 = np.mean(-differenceArray[:, 29])
        #plt.grid(True)
        #plt.xlabel('GPS time (s)')
        #plt.ylabel('Post - pre Hoft difference (lower is better)')
        #plt.title('Difference vs time by frequency')
        #plt.legend([p0, p1, p2], \
        #[str(frequencyArray[0, freqList[0]]) + ' Hz, mean: ' + str(mean0),\
        #str(frequencyArray[0, freqList[1]]) + ' Hz, mean: ' + str(mean1),\
        #str(frequencyArray[0, freqList[2]]) + ' Hz, mean: ' + str(mean2)])
        #plt.savefig(graphTitleDiff + 'Freq.png')
        #plt.savefig(graphTitleDiff + 'Freq.pdf')
        #plt.close()
        #plt.figure()
        #p0 = plt.scatter(timeArray, beforeArray[:, 27], color='b')
        #mean0 = np.mean(beforeArray[:, 27])
        #numberofthings = np.shape(beforeArray[:,27])[0]
        #mean0array = mean0 * np.ones((np.shape(beforeArray[:, 27])[0], 1))
        #p1 = plt.scatter(timeArray, afterArray[:, 27], color='r')
        #mean1 = np.mean(afterArray[:, 27])
        #mean1array = mean1 * np.ones((np.shape(afterArray[:, 27])[0], 1))
        #p2 = plt.plot(timeArray, mean0array, 'b^')
        #p3 = plt.plot(timeArray, mean1array, 'r+')
        #plt.grid(True)
        #plt.xlabel('GPS time (s)')
        #plt.ylabel('Hoft')
        #plt.title('Before and after feedforward')
        #plt.legend([p0, p1, p2, p3], \
        #[str(frequencyArray[0, freqList[1]]) + ' Hz before, scatterplot', \
        #str(frequencyArray[0, freqList[1]]) + ' Hz before, mean: ' + str(mean0), \
        #str(frequencyArray[0, freqList[1]]) + ' Hz after, scatterplot', \
        #str(frequencyArray[0, freqList[1]]) + ' Hz after, mean: ' + str(mean1), \
        #])
        #plt.savefig(graphTitleBeforeAfter + '.png')
        #plt.savefig(graphTitleBeforeAfter + '.pdf')
        #plt.legend
        #plt.close()
        histBins = 1e-24 * np.arange(50., 101.) 
        plt.figure()
        #p0 = plt.scatter(timeArray, beforeArray[:, 27], color='b')
        #p0 = plt.hist(beforeArray[:, 27], histBins)
        p1 = plt.hist(afterArray[:, 27], histBins)
        plt.grid(True)
        plt.xlabel('Hoft')
        plt.ylabel('Histogram count (1e-24 bins from 5e-23 to 1e-24)')
        plt.title('850 Hz comb histogram after feedforward')
        #plt.legend([p0, p1], ['Before', 'After'])
        plt.savefig(graphWholeRun + 'Histogram.png')
        plt.savefig(graphWholeRun + 'Histogram.pdf')
        plt.close()
            
    # Uncomment below to test on one directory only:
    #combSpectrum('/home/pulsar/public_html/feedforward/diagnostics/LHO/H-H1_AMPS_C02_L2-9326/', 'one')
    #combSpectrum('/home/pulsar/feedforward/2012/10/25/AMPS/files_for_comb/', 'one')
    # Uncomment below to test all directories and produce a whole-run overview.
    grandTarget = '/home/pulsar/public_html/feedforward/diagnostics/LHO/'
    combSpectrum(grandTarget, 'all')
    # Uncomment below to test each directory, making plots one-by-one.
    #grandTarget = '/home/pulsar/public_html/feedforward/diagnostics/LHO'
    #highDirectoryList = os.listdir(grandTarget)
    #highDirectoryListDirOnly = [x for x in highDirectoryList if x.find('.') == -1]
    #for x in highDirectoryListDirOnly:
    #    combSpectrum(grandTarget + '/' + x + '/', 'one')