GNUPLOT：繪製文本文件中由參數給出的分佈總和

Question

假設我有一個帶有任意數量行的文本文件，其中每行給出一組定義函數的參數（比如說（x，y）位置和sigma （可能性不等）的2D高斯）。例如，在這種情況下，文本文件可能包含：

100 112 3 4 
97 38 8 9 
88 79 3 9 
    ... 
    ... 
102 152 9 5 

我想（用pm3d）的SUM由文本文件中定義的所有發行的繪製。如何做到這一點？

Answer 1

我想繪製（使用pm3d）由文本文件定義的所有分佈的總和。如何做到這一點？

我認爲不能在原生gnuplot中完成，至少不能以我知道的任何理智方式完成。這種數字處理並不真正是gnuplot設計的目的。當您的示例數據運行

的Python，然而，使得它非常可行的......

#!/usr/bin/env python2 

import math 
import numpy 
import os 

class Gaussian(object): 
    '''A 2D gaussian function (normalized), defined by 
    mx: x mean position 
    my: y mean position 
    sx: sigma in x 
    sy: sigma in y''' 

    def __init__(self, mx, my, sx, sy): 
     self.mx = float(mx) 
     self.my = float(my) 
     self.sx = float(sx) 
     self.sy = float(sy) 

    def value(self,x,y): 
     '''Evaluates the value of a Gaussian at a certain point (x,y)''' 
    prefactor = (1.0/(self.sx*self.sy*2*math.pi)) 
    ypart = math.exp(-(x - self.mx)**2/(2*self.sx**2)) 
    xpart = math.exp(-(y - self.my)**2/(2*self.sy**2)) 
    return prefactor*ypart*xpart 

    def getLimits(self, devs): 
     '''Finds the upper and lower x and y limits for the distribution, 
     defined as points a certain number of deviations from the mean.''' 
     xmin = self.mx - devs*self.sx 
     xmax = self.mx + devs*self.sx 
     ymin = self.my - devs*self.sy 
     ymax = self.my + devs*self.sy 

     return (xmin, xmax, ymin, ymax) 

def readGaussians(filename): 
    '''reads in gaussian parameters from an input file in the format 
    mx my sx sy 

    This makes some assumptions about how perfect the input file will be''' 
    gaussians = [] 
    with open(filename, 'r') as f: 
     for line in f.readlines(): 
      (mx, my, sx, sy) = line.split() 
      gaussians.append(Gaussian(mx, my, sx, sy)) 

    return gaussians 

def getPlotLimits(gaussians, devs): 
    '''finds the x and y limits of the field of gaussian distributions. 
    Sets the boundary a set number of deviations from the mean''' 
    # get the limits from the first gaussian and work from there 
    (xminlim, xmaxlim, yminlim, ymaxlim) = gaussians[0].getLimits(devs) 
    for i in range(1, len(gaussians)): 
     (xmin, xmax, ymin, ymax) = gaussians[i].getLimits(devs) 
     if (xmin < xminlim): 
      xminlim = xmin 
     if (xmax > xmaxlim): 
      xmaxlim = xmax 
     if (ymin < yminlim): 
      yminlim = ymin 
     if (ymax > ymaxlim): 
      ymaxlim = ymax 

    return (xminlim, xmaxlim, yminlim, ymaxlim) 

def makeDataFile(gaussians, limits, res, outputFile): 
    '''makes a data file of x,y coordinates to be plotted''' 
    xres = res[0] 
    yres = res[1] 

    # make arrays of x and y values 
    x = numpy.linspace(limits[0], limits[1], xres) 
    y = numpy.linspace(limits[2], limits[3], yres) 
    # initialize grid of z values 
    z = numpy.zeros((xres, yres)) 

    # compute z value at each x, y point 
    for i in range(len(x)): 
     for j in range(len(y)): 
      for gaussian in gaussians: 
       z[i][j] += gaussian.value(x[i], y[j]) 

    # write out the matrix 
    with open(outputFile, 'w') as f: 
     for i in range(len(x)): 
      for j in range(len(y)): 
       f.write('%f %f %f\n' % (x[i], y[j], z[i][j])) 
      f.write('\n') 

def makePlotFile(limits, outputFile): 
    '''makes a plot file for gnuplot''' 
    with open('plot.plt', 'w') as f: 
     f.write("#!/usr/bin/env gnuplot\n") 
     f.write("set terminal png font 'Courier'\n") 
     f.write("set output 'gaussians.png'\n") 
     f.write("set xr [%f:%f]\n" % (limits[0], limits[1])) 
     f.write("set yr [%f:%f]\n" % (limits[2], limits[3])) 
     f.write("set pm3d map\n") 
     f.write("plot '%s' with image\n" % outputFile) 

    # make plot file executable 
    os.system('chmod 755 plot.plt') 

    # plot 
    os.system('./plot.plt') 

# name of input file 
inputFile = 'data.dat' 
# name of output (processed data) 
outputFile = 'gaussians.dat' 
# number of x and y points in plot 
res = (100, 100) 
# deviations from the mean by which to define Gaussian limits 
devs = 3 

# read in the gaussians from the data file 
print 'reading in data...' 
gaussians = readGaussians(inputFile) 

# find the plot limits 
limits = getPlotLimits(gaussians, devs) 

# make the gaussian data file 
print 'computing data for plotting...' 
makeDataFile(gaussians, limits, res, outputFile) 

# make the plot file 
print 'plotting...' 
makePlotFile(limits, outputFile) 

該腳本生成以下輸出。