text:
#!/usr/bin/env python3
import argparse
import sys
import parallel_mandle as parman
import numpy as np

#    origin=(-1.4011552,0) #easy Feigenbaum(?) point?
#def zoom_loop(start_win_size=2,fin_win_size=0.005,num_steps=2,origin=[-1.4011552,0],output_file="newzoom.gif",numprocs=2):
def zoom_loop(start_win_size,fin_win_size,num_steps,giftime,origin,output_file,numprocs):
#TODOs:
# Add MPI support
# -- 1 frame per cores?
# -- stripe frames across cores?
  gifstep=[]

  for i in np.geomspace(start_win_size,fin_win_size,num=num_steps):

    up_left_corn=[origin[0]-(i/2.0),origin[1]-(i/2.0)]
    gifstep.append(parman.mandlebrot_frame(upper_left_corner=up_left_corn,window_size=i,nprocs=numprocs))

  gifstep[0].save(output_file,format="GIF",append_images=gifstep[1:],save_all=True,duration=giftime,loop=0)

def main():
    #This will eventually hold the argparsing & maybe file-writing.
#Example file open
#  try:
#    with open(file_name,'x') as script_file: #"w" for write, "a" for append, "x" checks for file exist, "r" for read
#      filename.write(script_head+script_imports+script_main)
#  except Exception as err:
#    print("IOError:",err)

  parser=argparse.ArgumentParser()
  parser.add_argument("output_file",type=str,help="Output filename - in GIF format",default="newzoom.gif")
  parser.add_argument("-n","--nprocs",type=int,help="Number of multiprocessing threads to spawn",action="store",default=2)
  parser.add_argument("-d","--dur",type=int,help="Duration of GIF in microseconds",action="store",default=100)
  parser.add_argument("-s","--nsteps",type=int,help="Number of steps in the zoom",action="store",default=2)
  parser.add_argument("-w","--swin",type=float,help="Size of the starting frame",action="store",default=2)
  parser.add_argument("-f","--fwin",type=float,help="Size of the ending frame",action="store",default=0.005)
  parser.add_argument("-o","--origin",type=float,help="Center of the frame in c-space, space delimited",nargs=2,action="store",default=[-1.401152,0])
  args=parser.parse_args()
#  if args.flag:
#    print("Flag turned on")
#  print("This is an argument:",args.string_arg)
# Add args for starting boundary, origin, final boundary, steps
  zoom_loop(start_win_size=args.swin,fin_win_size=args.fwin,origin=args.origin,num_steps=args.nsteps,numprocs=args.nprocs,output_file=args.output_file,giftime=args.dur)

if __name__ == '__main__':
  main()

max_stars_repo_name: XSEDE/nix-container-python-mandle

max_issues_repo_name: XSEDE/nix-container-python-mandle

max_forks_repo_name: XSEDE/nix-container-python-mandle


text:
from keckdrpframework.primitives.base_primitive import BasePrimitive
from kcwidrp.core.bokeh_plotting import bokeh_plot
from kcwidrp.core.kcwi_plotting import get_plot_lims, oplot_slices, \
    set_plot_lims

from bokeh.plotting import figure
import numpy as np
import math
from scipy.interpolate import interpolate
from multiprocessing import get_context
from scipy import signal
import time

def pascal_shift(coefficients=None, x0=None):
    """Shift coefficients to a new reference value (X0)

    This should probably go somewhere else, but will be needed here.
    """
    if not coefficients:
        print("Error, no coefficients for pascal_shift.")
        return None
    if not x0:
        print("Warning, no reference value (x0) supplied")
        return coefficients
    if len(coefficients) == 7:
        usecoeff = list(reversed(coefficients))
        fincoeff = [0.] * 7
    else:
        if len(coefficients) > 7:
            print("Warning - this routine only handles up to 7 coefficients.")
            usecoeff = list(reversed(coefficients[0:7]))
            fincoeff = [0.] * len(coefficients)
        else:
            usecoeff = [0.] * 7
            fincoeff = usecoeff
            for ic, c in enumerate(coefficients):
                usecoeff[len(coefficients) - (ic + 1)] = coefficients[ic]
    # get reference values
    x01 = x0
    x02 = x0 ** 2
    x03 = x0 ** 3
    x04 = x0 ** 4
    x05 = x0 ** 5
    x06 = x0 ** 6
    # use Pascal's Triangle to shift coefficients
    fincoeff[0] = usecoeff[0] - usecoeff[1] * x01 + usecoeff[2] * x02 \
        - usecoeff[3] * x03 + usecoeff[4] * x04 - usecoeff[5] * x05 \
        + usecoeff[6] * x06

    fincoeff[1] = usecoeff[1] - 2.0 * usecoeff[2] * x01 \
        + 3.0 * usecoeff[3] * x02 - 4.0 * usecoeff[4] * x03 \
        + 5.0 * usecoeff[5] * x04 - 6.0 * usecoeff[6] * x05

    fincoeff[2] = usecoeff[2] - 3.0 * usecoeff[3] * x01 \
        + 6.0 * usecoeff[4] * x02 - 10.0 * usecoeff[5] * x03 \
        + 15.0 * usecoeff[6] * x04

    fincoeff[3] = usecoeff[3] - 4.0 * usecoeff[4] * x01 \
        + 10.0 * usecoeff[5] * x02 - 20.0 * usecoeff[6] * x03

    fincoeff[4] = usecoeff[4] - 5.0 * usecoeff[5] * x01 \
        + 15.0 * usecoeff[6] * x02

    fincoeff[5] = usecoeff[5] - 6.0 * usecoeff[6] * x01

    fincoeff[6] = usecoeff[6]
    # Trim if needed
    if len(coefficients) < 7:
        fincoeff = fincoeff[0:len(coefficients)]
    # Reverse for python
    return list(reversed(fincoeff))
    # END: def pascal_shift()

def bar_fit_helper(argument):

    b = argument['b']
    bs = argument['bs']
    # wavelength coefficients
    coefficients = [0., 0., 0., 0., 0.]
    # container for maxima, shifts
    maxima = []
    shifts = []
    # get sub spectrum for this bar
    sub_spectrum = bs[argument['minrow']:argument['maxrow']]
    # now loop over dispersions
    for di, dispersion in enumerate(argument['disps']):
        # populate the coefficients
        coefficients[4] = argument['p0'][b]
        coefficients[3] = dispersion
        cosbeta = dispersion / (argument['PIX'] * argument['ybin']) * \
            argument['rho'] * argument['FCAM'] * 1.e-4
        if cosbeta > 1.:
            cosbeta = 1.
        beta = math.acos(cosbeta)
        coefficients[2] = -(argument['PIX'] * argument['ybin'] /
                            argument['FCAM']) ** 2 * math.sin(beta) / 2. / \
            argument['rho'] * 1.e4
        coefficients[1] = -(argument['PIX'] * argument['ybin'] /
                            argument['FCAM']) ** 3 * math.cos(beta) / 6. / \
            argument['rho'] * 1.e4
        coefficients[0] = (argument['PIX'] * argument['ybin'] /
                           argument['FCAM']) ** 4 * math.sin(beta) / 24. / \
            argument['rho'] * 1.e4
        # what are the min and max wavelengths to consider?
        wl0 = np.polyval(coefficients, argument['xvals'][argument['minrow']])
        wl1 = np.polyval(coefficients, argument['xvals'][argument['maxrow']])
        minimum_wavelength = np.nanmin([wl0, wl1])
        maximum_wavelength = np.nanmax([wl0, wl1])
        # where will we need to interpolate to cross-correlate?
        minrw = [i for i, v in enumerate(argument['refwave'])
                 if v >= minimum_wavelength][0]
        maxrw = [i for i, v in enumerate(argument['refwave'])
                 if v <= maximum_wavelength][-1]
        ref_wave_of_sub_spectrum = argument['refwave'][minrw:maxrw]
        ref_flux_of_sub_spectrum = argument['reflux'][minrw:maxrw]
        # get bell cosine taper to avoid nasty edge effects
        tkwgt = signal.windows.tukey(len(ref_flux_of_sub_spectrum),
                                     alpha=argument['taperfrac'])
        # apply taper to atlas spectrum
        ref_flux_of_sub_spectrum *= tkwgt
        # adjust wavelengths
        waves = np.polyval(coefficients, argument['subxvals'])
        # interpolate the bar spectrum
        obsint = interpolate.interp1d(waves, sub_spectrum, kind='cubic',
                                      bounds_error=False,
                                      fill_value='extrapolate')
        intspec = obsint(ref_wave_of_sub_spectrum)
        # apply taper to bar spectrum
        intspec *= tkwgt
        # get a label
        # cross correlate the interpolated spectrum with the atlas spec
        samples_number = len(ref_wave_of_sub_spectrum)
        offsets_array = np.arange(1 - samples_number, samples_number)
        # Cross-correlate
        crosscorrelation = np.correlate(intspec, ref_flux_of_sub_spectrum,
                                        mode='full')
        # Get central region
        x0c = int(len(crosscorrelation) / 3)
        x1c = int(2 * (len(crosscorrelation) / 3))
        central_crosscorrelation = crosscorrelation[x0c:x1c]
        central_offsets_array = offsets_array[x0c:x1c]
        # Calculate offset
        maxima.append(central_crosscorrelation[
                          central_crosscorrelation.argmax()])
        shifts.append(central_offsets_array[central_crosscorrelation.argmax()])
    # Get interpolations
    int_max = interpolate.interp1d(argument['disps'], maxima, kind='cubic',
                                   bounds_error=False,
                                   fill_value='extrapolate')
    int_shift = interpolate.interp1d(argument['disps'], shifts, kind='cubic',
                                     bounds_error=False,
                                     fill_value='extrapolate')
    xdisps = np.linspace(min(argument['disps']), max(argument['disps']),
                         num=argument['nn'] * 100)
    # get central region
    x0c = int(len(xdisps) / 3)
    x1c = int(2 * (len(xdisps) / 3))
    # get peak values
    central_xdisps = xdisps[x0c:x1c]
    maxima_res = int_max(central_xdisps)
    shifts_res = int_shift(central_xdisps) * argument['refdisp']
    bardisp = central_xdisps[maxima_res.argmax()]
    barshift = shifts_res[maxima_res.argmax()]
    # update coeffs
    coefficients[4] = argument['p0'][b] - barshift
    coefficients[3] = bardisp
    cosbeta = coefficients[3] / (argument['PIX'] * argument['ybin']) * \
        argument['rho'] * argument['FCAM'] * 1.e-4
    if cosbeta > 1.:
        cosbeta = 1.
    beta = math.acos(cosbeta)
    coefficients[2] = -(argument['PIX'] * argument['ybin'] /
                        argument['FCAM']) ** 2 * \
        math.sin(beta) / 2. / argument['rho'] * 1.e4
    coefficients[1] = -(argument['PIX'] * argument['ybin'] /
                        argument['FCAM']) ** 3 * \
        math.cos(beta) / 6. / argument['rho'] * 1.e4
    coefficients[0] = (argument['PIX'] * argument['ybin'] /
                       argument['FCAM']) ** 4 * \
        math.sin(beta) / 24. / argument['rho'] * 1.e4
    shifted_coefficients = pascal_shift(coefficients, argument['x0'])
    print("Bar#: %3d, Cdisp: %.4f" % (b, bardisp))

    # Return results
    return b, shifted_coefficients, coefficients[4], coefficients[3], \
           maxima, bardisp
    # END: def bar_fit_helper()

class FitCenter(BasePrimitive):
    """ Fit central region"""

    def __init__(self, action, context):
        BasePrimitive.__init__(self, action, context)
        self.logger = context.pipeline_logger
        self.action.args.twkcoeff = []

    def _perform(self):
        """At this point we have the offsets between bars and the approximate
        offset from the reference bar to the atlas spectrum and the approximate
        dispersion.
        """
        self.logger.info("Finding wavelength solution for central region")
        # Are we interactive?
        do_inter = (self.config.instrument.plot_level >= 2)

        # y binning
        y_binning = self.action.args.ybinsize
        # let's populate the 0 points vector
        p0 = self.action.args.cwave + np.array(self.context.bar_offsets) * \
            self.context.prelim_disp - self.action.args.offset_wave
        # next we are going to brute-force scan around the preliminary
        # dispersion for a better solution. We will wander 5% away from it.
        maximum_dispersion_deviation = 0.05  # fraction
        # we will try nn values
        self.logger.info("prelim disp = %.3f, refdisp = %.3f,"
                         " min,max rows = %d, %d" % (self.context.prelim_disp,
                                                     self.action.args.refdisp,
                                                     self.action.args.minrow,
                                                     self.action.args.maxrow))
        number_of_values_to_try = (int(maximum_dispersion_deviation *
                                       abs(self.context.prelim_disp) /
                                   self.action.args.refdisp *
                                       (self.action.args.maxrow -
                                        self.action.args.minrow) / 2.0))
        if number_of_values_to_try < 10:
            number_of_values_to_try = 10
        if number_of_values_to_try > 50:
            number_of_values_to_try = 50
        self.logger.info("N disp. samples: %d" % number_of_values_to_try)
        # dispersions to try
        disps = self.context.prelim_disp * (
                1.0 + maximum_dispersion_deviation *
                (np.arange(0, number_of_values_to_try + 1) -
                 number_of_values_to_try / 2.) *
                2.0 / number_of_values_to_try)
        # values for central fit
        subxvals = self.action.args.xvals[
                   self.action.args.minrow:self.action.args.maxrow]

        # log taperfrac: important!
        self.logger.info("Using TAPERFRAC = %.3f" %
                         self.config.instrument.TAPERFRAC)
        self.action.args.ccddata.header['TAPFRAC'] = (
            self.config.instrument.TAPERFRAC, "taper fraction for central fit")
        # loop over bars and assemble input arguments
        my_arguments = []
        for b, bs in enumerate(self.context.arcs):
            arguments = {
                'b': b,
                'bs': bs,
                'minrow': self.action.args.minrow,
                'maxrow': self.action.args.maxrow,
                'disps': disps,
                'p0': p0,
                'PIX': self.config.instrument.PIX,
                'ybin': y_binning,
                'rho': self.action.args.rho,
                'FCAM': self.config.instrument.FCAM,
                'xvals': self.action.args.xvals,
                'refwave': self.action.args.refwave,
                'reflux': self.action.args.reflux,
                'taperfrac': self.config.instrument.TAPERFRAC,
                'refdisp': self.action.args.refdisp,
                'subxvals': subxvals,
                'nn': number_of_values_to_try, 'x0': self.action.args.x0
            }
            my_arguments.append(arguments)

        twkcoeff = {}
        centwave = []
        centdisp = []

        p = get_context("spawn").Pool()
        results = p.map(bar_fit_helper, list(my_arguments))
        p.close()

        next_bar_to_plot = 0
        for ir, result in enumerate(results):
            b = result[0]
            shifted_coefficients = result[1]
            _centwave = result[2]
            _centdisp = result[3]
            twkcoeff[b] = shifted_coefficients
            centwave.append(_centwave)
            centdisp.append(_centdisp)
            maxima = result[4]
            bardisp = result[5]
            self.logger.info("Central Fit: Bar# %3d, Cdisp %.4f, "
                             "Coefs: %.2f %.4f %13.5e %13.5e" %
                             (b, bardisp, shifted_coefficients[4],
                              shifted_coefficients[3], shifted_coefficients[2],
                              shifted_coefficients[1]))
            if do_inter and ir == next_bar_to_plot:
                # plot maxima
                p = figure(title=self.action.args.plotlabel +
                           "CENTRAL DISPERSION FIT for Bar: %d Slice: %d" %
                           (b, int(b / 5)),
                           plot_width=self.config.instrument.plot_width,
                           plot_height=self.config.instrument.plot_height,
                           x_axis_label="Central dispersion (Ang/px)",
                           y_axis_label="X-Corr Peak Value")
                p.scatter(disps, maxima, color='red', legend_label="Data")
                p.line(disps, maxima, color='blue', legend_label="Data")
                ylim = [min(maxima), max(maxima)]
                p.line([_centdisp, _centdisp], ylim, color='green',
                       legend_label="Fit Disp")
                p.line([self.context.prelim_disp, self.context.prelim_disp],
                       ylim, color='red', legend_label="Calc Disp")
                bokeh_plot(p, self.context.bokeh_session)
                q = input("Next? <int> or <cr>, q to quit: ")
                if 'Q' in q.upper():
                    do_inter = False
                else:
                    try:
                        next_bar_to_plot = int(q)
                    except ValueError:
                        next_bar_to_plot = ir + 1

        self.action.args.twkcoeff = twkcoeff
        # Plot results
        if self.config.instrument.plot_level >= 1:
            # Plot central wavelength
            p = figure(title=self.action.args.plotlabel + "CENTRAL VALUES",
                       x_axis_label="Bar #",
                       y_axis_label="Central Wavelength (A)",
                       plot_width=self.config.instrument.plot_width,
                       plot_height=self.config.instrument.plot_height)
            x = range(len(centwave))
            p.scatter(x, centwave, marker='x', legend_label='bar wave')
            p.line([0, 120], [self.action.args.cwave, self.action.args.cwave],
                   color='red', legend_label='CWAVE')
            xlim = [-1, 120]
            ylim = get_plot_lims(centwave)
            p.xgrid.grid_line_color = None
            oplot_slices(p, ylim)
            p.legend.location = "top_center"
            set_plot_lims(p, xlim=xlim, ylim=ylim)
            bokeh_plot(p, self.context.bokeh_session)
            if self.config.instrument.plot_level >= 2:
                input("Next? <cr>: ")
            else:
                time.sleep(self.config.instrument.plot_pause)
            # Plot central dispersion
            p = figure(title=self.action.args.plotlabel + "CENTRAL VALUES",
                       x_axis_label="Bar #",
                       y_axis_label="Central Dispersion (A)",
                       plot_width=self.config.instrument.plot_width,
                       plot_height=self.config.instrument.plot_height)
            x = range(len(centdisp))
            p.scatter(x, centdisp, marker='x', legend_label='bar disp')
            p.line([0, 120], [self.context.prelim_disp,
                              self.context.prelim_disp], color='red',
                   legend_label='Calc Disp')
            xlim = [-2, 121]
            ylim = get_plot_lims(centdisp)
            p.xgrid.grid_line_color = None
            oplot_slices(p, ylim)
            p.legend.location = "bottom_center"
            set_plot_lims(p, xlim=xlim, ylim=ylim)
            bokeh_plot(p, self.context.bokeh_session)
            if self.config.instrument.plot_level >= 2:
                input("Next? <cr>: ")
            else:
                time.sleep(self.config.instrument.plot_pause)

        log_string = FitCenter.__module__
        self.action.args.ccddata.header['HISTORY'] = log_string
        self.logger.info(log_string)

        return self.action.args
    # END: class FitCenter()

ext: py

reason: import numpy,from scipy


text:
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn import preprocessing
import pylab
import numpy.matlib
import itertools
from line_coding import polar_nrz
import math
import mpmath
from sympy.combinatorics.graycode import gray_to_bin

# m: order of modualtion
# binary bits: a matrix of binary bits. the # of rows is the length of each bin
# gray bits : a matrix of binary bits. the # of rows is the length of each bin

__all__=['binary_bits_to_gray_bits','gray_bits_to_binary_bits','bits_to_str','str_to_bits','gray_str_to_binary_str',
			'binary_str_to_int','bits_to_gray_to_int','bits_to_binary_to_int','int_to_gray_bits','int_to_binary_bits']

def binary_bits_to_gray_bits(binary_bits):
	binary_bits=np.transpose(binary_bits)
	rows,cols=binary_bits.shape
	gray_bits=np.zeros((rows,cols))
	for i in np.arange(rows):
		gray_bits[i,:]=np.hstack((binary_bits[i,0],
			(binary_bits[i,1:cols]^binary_bits[i,0:cols-1])))
	return np.transpose(gray_bits).astype(int)

def gray_bits_to_binary_bits(gray_bits):
	gray_bits=np.transpose(gray_bits)
	rows,cols=gray_bits.shape
	binary_bits=np.zeros((rows,cols)).astype(int)
	for i in np.arange(rows):
		binary_bits[i,0]=gray_bits[i,0]
		for j in np.arange(cols)[1:]:
			binary_bits[i,j]=binary_bits[i,j-1]^gray_bits[i,j]
	return np.transpose(binary_bits).astype(int)

def bits_to_str(bits,m):
	k=int(np.log2(m))
	string=[]
	for i in np.arange(bits.shape[1]):
		b=''
		for j in np.arange(k):
			b+=str(bits[j,i])
		string.append(b)
	return string

def str_to_bits(string,m):
	k=int(np.log2(m))
	bits=np.zeros((k,len(string)))
	for i in np.arange(len(string)):
		for j in np.arange(k):
			bits[j,i]=string[i][j]
	return bits.astype(int)

def gray_str_to_binary_str(gray_str):
	binary_str=[]
	for i in np.arange(len(gray_str)):
		binary_str.append(gray_to_bin(gray_str[i]))
	return binary_str

def binary_str_to_int(binary_str):
	ints=np.array([int(x,2) for x in binary_str])
	return ints

def bits_to_gray_to_int(input_bits,m):

	k=np.log2(m)
	reshape_bits=np.reshape(input_bits,(len(input_bits)/k,k))
	gray_bits=gray_bits_to_binary_bits(np.transpose(reshape_bits))
	gray_str=bits_to_str(gray_bits,m)
	integer=np.array([int(x,2) for x in gray_str])
	return integer

def bits_to_binary_to_int(input_bits,m):
	k=np.log2(m)
	reshape_bits=np.reshape(input_bits,(len(input_bits)/k,k))
	reshape_bits=np.transpose(reshape_bits)
	binary_str=bits_to_str(reshape_bits,m)
	integer=np.array([int(x,2) for x in binary_str])
	return integer

def int_to_gray_bits(input_int,m):
	k=int(np.log2(m))
	binary_str=np.array([format(x,'0%sb'%(k)) for x in input_int])
	binary_bits=str_to_bits(binary_str,m)
	gray_bits=binary_bits_to_gray_bits(binary_bits)
	return gray_bits

def int_to_binary_bits(input_int,m):
	k=int(np.log2(m))
	binary_str=np.array([format(x,'0%sb'%(k)) for x in input_int])
	binary_bits=str_to_bits(binary_str,m)
	return binary_bits

# ...test
# m=8
# k=np.log2(m)
# nb=k*10**3
# input_bits=np.random.randint(2,size=nb)
# integer=bits_to_gray_to_int(input_bits,m)
# integer2=bits_to_binary_to_int(input_bits,m)
# bits=str_to_bits(gray_str,m)
# input_int=np.array([1,2,3])
# b=int_to_gray_bits(input_int,8)
# c=int_to_binary_bits(input_int,8)

reason: import numpy,from sympy,import mpmath


text:
import os
import sys

from gensim.models.doc2vec import Doc2Vec
import numpy
from sklearn.svm import SVC

from preprocessing import normalize, micro_tokenize
import conf

class D2V_feature_extractor(object):
    def __init__(self, d2vfile):
        self.d2v = Doc2Vec.load(d2vfile)

    def fit(self, *args, **kwargs):
        pass

    def transform(self, txts):
        res = []
        for txt in txts:
            # see https://github.com/RaRe-Technologies/gensim/issues/447
            self.d2v.random.seed(conf.SEED)
            v = self.d2v.infer_vector(micro_tokenize(normalize(txt)))
            res.append(v)
        return numpy.vstack(res)

def evaluate(cat, fold, txt_train, txt_test, y_train, y_test):
    d2vmodelfile = os.path.join(conf.D2V_DIR, 'model')
    if not os.path.exists(d2vmodelfile):
        print('Doc2vec model file "%s" not found' % d2vmodelfile)
        print('Did you run train_doc2vec.py?')
        sys.exit(1)

    fe = D2V_feature_extractor(d2vmodelfile)
    predictor = SVC(
        kernel=conf.SVM_KERNEL,
        class_weight=conf.SVM_CLWEIGHT,
        C=conf.SVM_C,
        random_state=conf.SEED,
    )
    X = fe.transform(txt_train)
    predictor.fit(X, y_train)
    X_test = fe.transform(txt_test)
    y_pred = predictor.predict(X_test)

    return y_pred

max_stars_repo_name: aseifert/million-post-corpus

max_issues_repo_name: aseifert/million-post-corpus

max_forks_repo_name: aseifert/million-post-corpus

reason: import numpy


text:
import argparse
import os
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
import numpy as np
import pandas as pd
from tabulate import tabulate

from cgcnn.model import CGCNN
from cgcnn.model import CGCNN_pred
from cgcnn.util import MaskAtom
from cgcnn.dataloader import DataLoaderMasking
from cgcnn.CIFDATA_pred import CIF_Cry_Dataset
from cgcnn.my_logging import *

criterion = nn.CrossEntropyLoss()

def save_checkpoint(state, epoch):
    filename = 'checkingpoint' + '_' + str(epoch) + '_' + '.pth'
    torch.save(state, filename)

def output_training(loss_avg, acc_avg, epoch, log):
    header_1, header_2 = 'Loss | e-stop', 'Acc | TIME'

    train_loss1 = loss_avg['train_loss1'][epoch]
    train_loss2 = loss_avg['train_loss2'][epoch]
    train_loss = loss_avg['train_loss'][epoch]
    train_loss_bond = loss_avg['train_loss_bond'][epoch]

    train_acc1_atom = acc_avg['train_acc1_atom'][epoch]
    train_acc2_atom = acc_avg['train_acc2_atom'][epoch]
    train_acc_atom = acc_avg['train_acc_atom'][epoch]
    train_acc_bond = acc_avg['train_acc_bond'][epoch]

    tab_val = [['Node_prop1', f'{train_loss1:.4f}', f'{train_acc1_atom:.4f}'],
               ['Node_prop2', f'{train_loss2:.4f}', f'{train_acc2_atom:.4f}'],
               ['Node_total', f'{train_loss:.4f}', f'{train_acc_atom:.4f}'],
               ['Bond', f'{train_loss_bond:.4f}', f'{train_acc_bond:.4f}']]

    output = tabulate(tab_val, headers=[f'EPOCH # {epoch}', header_1, header_2], tablefmt='grid_tables')
    print(output)
    log.info("Epoch: {epoch}\t"
             "Node_prop1: {train_acc1_atom}\t"
             "Node_prop2: {train_acc2_atom}\t"
             "Node_total: {train_acc_atom}\t"
             "Bond: {train_acc_bond}\t".format(epoch=epoch,
                                               train_acc1_atom=train_acc1_atom,
                                               train_acc2_atom=train_acc2_atom,
                                               train_acc_atom=train_acc_atom,
                                               train_acc_bond=train_acc_bond))

    return output

def compute_accuracy(pred, target):
    return float(torch.sum(torch.max(pred.detach(), dim=1)[1] == target).cpu().item()) / len(pred)

def train(args, model_list, loader, optimizer_list, device, loss_avg, acc_avg):
    model, linear_pred_atoms, linear_pred_bonds = model_list
    optimizer_model, optimizer_linear_pred_atoms, optimizer_linear_pred_bonds = optimizer_list
    # training mode
    model.train()
    linear_pred_atoms.train()
    linear_pred_bonds.train()
    # loss and accuracy for node and edge respectively
    loss1_accum = 0
    acc1_node_accum = 0
    loss2_accum = 0
    acc2_node_accum = 0
    loss_accum = 0
    acc_node_accum = 0
    loss_edge_accum = 0
    acc_edge_accum = 0
    # train over batches
    for step, batch in enumerate(loader):

        batch = batch.to(device)
        # set the atomic vector of masked atoms to zero vector
        batch.x[batch.masked_atom_indices] = torch.zeros([len(batch.masked_atom_indices), batch.x.size()[1]])
        node_rep = model(batch)
        # predicted vector of masked atoms
        pred_node = linear_pred_atoms(node_rep[batch.masked_atom_indices])
        # predict group vector of masked atoms
        pred_node_group_number = pred_node[:, :19]
        # predict period vector of masked atoms
        pred_node_periodic_number = pred_node[:, 19:28]
        # true group number of masked atoms
        true_node_group_number = batch.mask_node_label[:, :19]
        # true period number of masked atoms
        true_node_periodic_number = batch.mask_node_label[:, 19:28]
        true_node_group_number_idx = [list(true_node_group_number[i]).index(1) for i in range(true_node_group_number.size()[0])]
        true_node_group_number_idx = torch.LongTensor(true_node_group_number_idx)
        true_node_periodic_number_idx = [list(true_node_periodic_number[i]).index(1) for i in range(true_node_periodic_number.size()[0])]
        true_node_periodic_number_idx = torch.LongTensor(true_node_periodic_number_idx)
        # calculate the loss of node task
        loss1 = criterion(pred_node_group_number.double(), true_node_group_number_idx)
        loss2 = criterion(pred_node_periodic_number.double(), true_node_periodic_number_idx)
        loss = loss1 + loss2
        # calculate the accuracy of node task
        acc1_node = compute_accuracy(pred_node_group_number, true_node_group_number_idx)
        acc2_node = compute_accuracy(pred_node_periodic_number, true_node_periodic_number_idx)
        acc_node = (acc1_node + acc2_node) / 2
        acc1_node_accum += acc1_node
        acc2_node_accum += acc2_node
        acc_node_accum += acc_node
        # edge task
        if args.mask_edge:
            masked_edge_index = batch.edge_index[:, batch.connected_edge_indices]
            # take the sum over connected nodes as the representation of masked edges
            edge_rep = node_rep[masked_edge_index[0]] + node_rep[masked_edge_index[1]]
            pred_edge = linear_pred_bonds(edge_rep)
            true_edge_idx = []
            for i in range(batch.mask_edge_label.size()[0]):
                max_val = batch.mask_edge_label[i].max()
                max_index = list(batch.mask_edge_label[i]).index(max_val)
                true_edge_idx.append(max_index)
            true_edge_idx = torch.LongTensor(true_edge_idx)
            # calculate the loss of edge task
            loss_edge = criterion(pred_edge.double(), true_edge_idx)
            loss += loss_edge
            loss_edge_accum += loss_edge
            # calculate the accuracy of edge task
            acc_edge = compute_accuracy(pred_edge, true_edge_idx)
            acc_edge_accum += acc_edge

        optimizer_model.zero_grad()
        optimizer_linear_pred_atoms.zero_grad()
        optimizer_linear_pred_bonds.zero_grad()

        loss.backward()

        optimizer_model.step()
        optimizer_linear_pred_atoms.step()
        optimizer_linear_pred_bonds.step()

        loss1_accum += float(loss1.cpu().item())
        loss2_accum += float(loss2.cpu().item())
        loss_accum += float(loss.cpu().item())
    # average over batches
    loss_avg['train_loss1'].append(loss1_accum / (step + 1))
    loss_avg['train_loss2'].append(loss2_accum / (step + 1))
    loss_avg['train_loss'].append(loss_accum / (step + 1))
    loss_avg['train_loss_bond'].append(loss_edge_accum / (step + 1))
    acc_avg['train_acc1_atom'].append(acc1_node_accum / (step + 1))
    acc_avg['train_acc2_atom'].append(acc2_node_accum / (step + 1))
    acc_avg['train_acc_atom'].append(acc_node_accum / (step + 1))
    acc_avg['train_acc_bond'].append(acc_edge_accum / (step + 1))

    return loss_avg, acc_avg

def main():
    # Training settings
    parser = argparse.ArgumentParser(description='PyTorch implementation of self-supervised training of '
                                                 'graph neural networks')
    parser.add_argument('root', help='The path where you define raw and processed dir')
    parser.add_argument('--device', type=int, default=0,
                        help='which gpu to use if any (default: 0)')
    parser.add_argument('--batch_size', type=int, default=256,
                        help='input batch size for training (default: 256)')
    parser.add_argument('--epochs', type=int, default=100,
                        help='number of epochs to train (default: 100)')
    parser.add_argument('--start-epoch', default=0, type=int, metavar='N',
                        help='manual epoch number (useful on restarts)')
    parser.add_argument('--lr', type=float, default=0.001,
                        help='learning rate (default: 0.001)')
    parser.add_argument('--decay', type=float, default=0,
                        help='weight decay (default: 0)')
    parser.add_argument('--num_layer', type=int, default=5,
                        help='number of GNN message passing layers (default: 5).')
    parser.add_argument('--emb_dim', type=int, default=64,
                        help='embedding dimensions (default: 64)')
    parser.add_argument('--dropout_ratio', type=float, default=0,
                        help='dropout ratio (default: 0)')
    parser.add_argument('--mask_rate', type=float, default=0.15,
                        help='dropout ratio (default: 0.15)')
    parser.add_argument('--mask_edge', type=int, default=1,
                        help='whether to mask edges or not together with atoms')
    parser.add_argument('--output_model_file', type=str, default='', help='filename to output the model')
    parser.add_argument('--model_type', type=str, default="cgcnn", choices=["cgcnn", "gatgnn", "None"])
    parser.add_argument('--seed', type=int, default=0, help="Seed for splitting dataset.")
    parser.add_argument('--num_workers', type=int, default=1, help='number of workers for dataset loading')
    parser.add_argument('--resume', default='', type=str, metavar='PATH',
                        help='path to latest checkpoint (default: none)')
    args = parser.parse_args()

    torch.manual_seed(0)
    np.random.seed(0)
    # device = torch.device("cuda:" + str(args.device)) if torch.cuda.is_available() else torch.device("cpu")
    device = 'cpu'
    if torch.cuda.is_available():
        torch.cuda.manual_seed_all(0)

    print("num layer: %d mask rate: %f mask edge: %d" % (args.num_layer, args.mask_rate, args.mask_edge))

    cwd = os.getcwd()
    log = log_creater(cwd)

    # set up dataset and transform function.
    dataset = CIF_Cry_Dataset(root=args.root,
                              transform=MaskAtom(num_atom_type=92,
                                                 num_edge_type=41,
                                                 mask_rate=args.mask_rate,
                                                 mask_edge=args.mask_edge))

    print("This is length of dataset: ", len(dataset))

    loader = DataLoaderMasking(dataset, batch_size=args.batch_size, shuffle=True, num_workers=args.num_workers)

    # set up models
    if args.model_type == 'cgcnn':
        model = CGCNN(n_conv=args.num_layer,
                      atom_fea_len=args.emb_dim,
                      drop_ratio=args.dropout_ratio).to(device)
        print(model)

    elif args.model_type == 'gatgnn':
        raise NotImplementedError("GATGNN model is not Implemented at present!!!")

    elif args.model_type == "None":
        raise NotImplementedError("Other GNN model is not Implemented at present!!!")

    else:
        raise NotImplementedError("This model is not Implemented !!!")

    linear_pred_atoms = torch.nn.Linear(args.emb_dim, 92).to(device)
    linear_pred_bonds = torch.nn.Linear(args.emb_dim, 41).to(device)

    model_list = [model, linear_pred_atoms, linear_pred_bonds]

    # set up optimizers
    optimizer_model = optim.Adam(model.parameters(), lr=args.lr, weight_decay=args.decay)
    optimizer_linear_pred_atoms = optim.Adam(linear_pred_atoms.parameters(), lr=args.lr, weight_decay=args.decay)
    optimizer_linear_pred_bonds = optim.Adam(linear_pred_bonds.parameters(), lr=args.lr, weight_decay=args.decay)
    optimizer_list = [optimizer_model, optimizer_linear_pred_atoms, optimizer_linear_pred_bonds]

    if args.resume:
        print("=> loading checkpoint '{}'".format(args.resume))
        checkpoint = torch.load(args.resume)
        model.load_state_dict(checkpoint['model_state_dict'])
        linear_pred_atoms.load_state_dict(checkpoint['linear_pred_atoms_state_dict'])
        linear_pred_bonds.load_state_dict(checkpoint['linear_pred_bonds_state_dict'])
        optimizer_model.load_state_dict(checkpoint['optimizer'])
        optimizer_linear_pred_atoms.load_state_dict(checkpoint['optimizer_linear_pred_atoms'])
        optimizer_linear_pred_bonds.load_state_dict(checkpoint['optimizer_linear_pred_bonds'])

    # define dictionaries to record loss and error
    global loss_avg, acc_avg
    loss_avg = {}
    acc_avg = {}
    loss_avg['train_loss1'] = []
    loss_avg['train_loss2'] = []
    loss_avg['train_loss'] = []
    loss_avg['train_loss_bond'] = []
    acc_avg['train_acc1_atom'] = []
    acc_avg['train_acc2_atom'] = []
    acc_avg['train_acc_atom'] = []
    acc_avg['train_acc_bond'] = []

    # start training loop
    for epoch in range(1, args.epochs + 1):
        # loss_avg and acc_avg save average of loss and accuracy of each epoch
        loss_avg, acc_avg = train(args, model_list, loader, optimizer_list, device, loss_avg, acc_avg)

        output_training(loss_avg, acc_avg, epoch - 1, log)

        if epoch % 5 == 0:
            # torch.save(model.state_dict(), args.output_model_file + '_' + str(epoch) + ".pth")
            save_checkpoint({
                'model_state_dict': model.state_dict(),
                'linear_pred_atoms_state_dict': linear_pred_atoms.state_dict(),
                'linear_pred_bonds_state_dict': linear_pred_bonds.state_dict(),
                'optimizer': optimizer_model.state_dict(),
                'optimizer_linear_pred_atoms': optimizer_linear_pred_atoms.state_dict(),
                'optimizer_linear_pred_bonds': optimizer_linear_pred_bonds.state_dict(),
            }, epoch)

    np.save('loss_avg.npy', loss_avg)
    np.save('acc_avg.npy', acc_avg)

if __name__ == "__main__":
    main()


text:
import re
import sys
import os
import time
import tempfile
from struct import pack
import copy
import math
import random
from subprocess import Popen, PIPE, call
import numpy as num
#import matplotlib
#matplotlib.use('PDF')
#import matplotlib.pylab as lab
import pickle
import scipy.optimize

pjoin = os.path.join

def kiwi_aux_dir():
    ieq_home = os.getenv('KIWI_HOME')
    if ieq_home is None:
         sys.exit('KIWI_HOME environment variable not set')
    d = pjoin(ieq_home, 'aux')
    if not os.path.isdir(d):
        sys.exit('directory not found: "%s"' % d)
    return d

def dump(data, filename):
    f = open(filename, 'w')
    pickle.dump( data, f )
    f.close()

class Config:
    seismosizer_prog  = "minimizer"
    source_info_prog  = "source_info"
    gfdb_info_prog = "gfdb_info"
    plotstyle = [ "1p/0/0/200", "1p/200/0/0", "1p/0/150/50", "1p/226/113/0", "1p/226/0/113" ]
    verbose = 0
    component_names = { 'a':'R@-+@-',
                        'c':'R@--@-',
                        'r':'T@-+@-',
                        'l':'T@--@-',
                        'd':'Z@-+@-',
                        'u':'Z@--@-',
                        'n':'N',
                        'e':'E',
                        's':'S',
                        'w':'W' }
    earthradius = 6371.*1000.

class Minimizer:
    """This is a wrapper to the seismogram calculation part of a minimizer process (see minimizer.f90),
       allowing on-the-fly conversion of the seismosizer output files to vtk files.
       Communication with the child seismosizer process is done via pipes,
       connected to stdin and stdout of the child process.
       Furthermore, at startup, it queries source_info (source_info.f90) for
       information about the possible source model parameterizations, which it
       stores in self.params"""

    commands = ['set_database',
                'set_local_interpolation',
                'set_spacial_undersampling',
                'set_receivers',
                'switch_receiver',
                'set_ref_seismograms',
                'set_source_location',
                'set_source_crustal_thickness_limit',
                'get_source_crustal_thickness',
                'set_source_params',
                'set_source_params_mask',
                'set_source_subparams',
                'set_effective_dt',
                'set_misfit_method',
                'set_misfit_filter',
                'set_misfit_taper',
                'set_synthetics_factor',
                'minimize_lm',
                'output_seismograms',
                'output_seismogram_spectra',
                'output_source_model',
                'get_source_subparams',
                'get_global_misfit',
                'get_misfits',
                'get_principal_axes',
                'output_distances',
                'output_cross_correlations',
                'shift_ref_seismogram',
                'autoshift_ref_seismogram',
                'set_verbose']

    def __init__(self):
        """start minimizer program"""

        cmd = [Config.source_info_prog]

        self.tempdir = tempfile.mkdtemp("","minimizer-")
        self.tempfilebase = self.tempdir + "/s"

        cmd = [Config.seismosizer_prog]

        try:
            p = Popen( cmd, stdin=PIPE, stdout=PIPE )
            self.p = p
            self.to_p = p.stdin
            self.from_p = p.stdout
        except:
            sys.exit("cannot start minimizer")

    def close(self):
        self.to_p.close()
        self.from_p.close()

    def do(self, command_name, *args):
        """Put commands to minimizer and return the results"""

        command  = command_name+' '+' '.join( [str(arg) for arg in args ] )

        answer = None
        if Config.verbose >= 2:
            sys.stderr.write( command+"\n" )

        self.to_p.write(command+"\n")
        self.to_p.flush()
        retval = self.from_p.readline().rstrip()

        if Config.verbose >= 3:
            sys.stderr.write( retval+"\n" )

        if retval.endswith('nok'):
            sys.exit("minimizer failed doing command: "+ command)

        if retval.endswith('ok >'):
            answer = self.from_p.readline().rstrip()
            if Config.verbose >= 3:
                sys.stderr.write( answer+"\n" )

        return answer

    def set_receivers( self, lat=None, lon=None, components=None, filename=None ):
        """Set minimizers receiver list to given lists or from file"""
        names = []
        if filename != None:
            file = open(filename, "r")
            lat = []
            lon = []
            components = []
            for line in file:
                toks = line.split()
                if len(toks) >= 3:
                    lat.append(float(toks[0]))
                    lon.append(float(toks[1]))
                    components.append(toks[2])
                if len(toks) == 4:
                    names.append(toks[3])
            file.close()

        receiverfn = self.tempdir+"/receivers"

        self.lat = lat
        self.lon = lon
        self.components = components
        self.receiver_names = names
        self.nreceivers = len(lat)

        file = open(receiverfn, "w")
        for i in xrange(len(lat)):
            file.write( ' '.join( (str(lat[i]), str(lon[i]), components[i]) ) + "\n" )

        file.close()

        self.do("set_receivers", receiverfn)
        os.remove( receiverfn )

    def set_source( self, source ):
        self.do("set_source_params", source )

    def set_source_receiver_circle( self, nreceivers, epidist, comps='ard' ):
        """set source to pole and layout n receivers around it at distance epidist"""

        self.do("set_source_location", 90, 0, 0)
        lat = []
        lon = []
        compo = []

        epidist = 5. # degrees
        for i in xrange(nreceivers):
            lat.append( 90. - epidist )
            lon.append( i*360./nreceivers )
            compo.append( comps )

        self.set_receivers( lat, lon, compo )

    def set_source_receiver_zigzag( self, nreceivers, epidist, comps='ard' ):
        """set source to pole and layout n receivers around it at distance epidist"""

        self.do("set_source_location", 90, 0, 0)
        lat = []
        lon = []
        compo = []

        epidist = 5. # degrees
        for i in xrange(nreceivers):
            if i % 2:
                lat.append( 90. - epidist )
            else:
                lat.append( 90. - epidist*1.5 )
            lon.append( i*360./nreceivers )
            compo.append( comps )

        self.set_receivers( lat, lon, compo )

    def set_source_receiver_random( self, nreceivers, mindist, maxdist, comps='ard' ):
        """set source to pole and layout n receivers around it at distance epidist"""

        self.do("set_source_location", 90, 0, 0)
        lat = []
        lon = []
        compo = []

        epidist = 5. # degrees
        for i in xrange(nreceivers):
            epidist = ( (maxdist-mindist)*random.random()+mindist ) * 360. / (2.0 * math.pi * Config.earthradius)
            lat.append( 90. - epidist )
            lon.append( random.random()*360. )
            compo.append( comps )

        self.set_receivers( lat, lon, compo )

    def set_synthetic_reference( self, source ):
        """Calculate seismograms for given source and use these as reference"""

        tempfnbase = self.tempdir + "/reference"
        self.set_source( source )
        self.do("output_seismograms", tempfnbase, "mseed", "synthetics", "plain")
        self.do("set_ref_seismograms", tempfnbase, "mseed")

    def get_misfits_for_source( self, source, source_has_been_set=False ):
        """Calculate misfit for given source"""

        if not source_has_been_set:
            self.set_source( source )
        values = [ float(x) for x in self.do_get_misfits().split() ]
        misfits_by_r = []
        norms_by_r = []
        for irec in xrange(len(self.components)):
            misfits_by_c = []
            norms_by_c = []
            for icomp in xrange(len(self.components[irec])):
                misfits_by_c.append( values.pop(0) )
                norms_by_c.append( values.pop(0) )

            misfits_by_r.append(misfits_by_c)
            norms_by_r.append(norms_by_c)

        return source, misfits_by_r, norms_by_r

    def get_misfits_for_sources( self, sources,  prefunc=None ):
        results = []
        for source in sources:
            if prefunc is not None:
                source = prefunc(self,source)
            results.append( self.get_misfits_for_source( source ) )

        return results

    def global_misfit( self, misfits_by_r, norms_by_r ):
        smisfit = 0.
        for misfit_by_c, norms_by_c in zip(misfits_by_r, norms_by_r):
            sm = 0.
            sn = 0.
            for misfit, norm  in zip(misfit_by_c, norms_by_c):
                if norm > 0.:
                    sm += misfit**2
                    sn += norm**2

            if sn > 0.0:
                smisfit += sm / sn

        return math.sqrt(smisfit)

    def get_best_factor_for_source( self, source ):

        def func(factor):
            self.do_set_synthetics_factor( factor )
            misfits_by_r, norms_by_r = self.get_misfits_for_source(source, source_has_been_set=True )[1:]
            return self.global_misfit(misfits_by_r, norms_by_r)

        self.set_source(source)
        result = scipy.optimize.brent(func, brack=(0.,1.), tol=1.0e-4)

        self.do_set_synthetics_factor(1.)
        return result

    def get_grid_minimum( self, sm_grid ):
        smg_grid = []
        for (source, misfits, norms) in sm_grid:
            glob_misfit = self.global_misfit( misfits, norms )
            smg_grid.append( (source, glob_misfit) )

        misfit, s =  min( [(xm,xs) for (xs,xm) in smg_grid ] )
        return s, misfit

    def get_grid_maximum( self, sm_grid ):
        smg_grid = []
        for (source, misfits, norms) in sm_grid:
            glob_misfit = self.global_misfit( misfits, norms )
            smg_grid.append( (source, glob_misfit) )

        misfit, s =  max( [(xm,xs) for (xs,xm) in smg_grid ] )
        return s, misfit

    def get_grid_minimum_bootstrap( self, sm_grid ):

        # randomization
        nrecs = len(sm_grid[0][1])
        sample = []
        for i in xrange(nrecs):
            sample.append( random.randint(0,nrecs-1) )

        # construct misfit and norm vectors for this randomization
        smb_grid = []
        for (source, misfits, norms) in sm_grid:
            bmisfits = []
            bnorms = []
            for irec in xrange(nrecs):
                bmisfits.append( misfits[sample[irec]] )
                bnorms.append( norms[sample[irec]] )
            smb_grid.append((source, bmisfits, bnorms))

        # get minimum in the same way as usual
        return self.get_grid_minimum( smb_grid )

    def grid_bootstrap( self, sm_grid, niter,  sourceparams=None ):
        sources = []
        for i in xrange(niter):
            source, misfit = self.get_grid_minimum_bootstrap( sm_grid )
            sources.append( source )

        if sourceparams==None:
            sourceparams = sources[0].params.keys()

        meansource = Source(sources[0].sourcetype)
        meansourceparams = {}
        stddevsourceparams = {}
        bresults = {}
        for param in sourceparams:
            param_results = num.zeros(niter, dtype=num.float)
            param_histogram = {}
            for i, source in enumerate(sources):
                value = source.params[param]
                param_results[i] = value
                param_histogram[value] = param_histogram.get(value, 0) + 1

            meansourceparams[param] = param_results.mean()
            stddevsourceparams[param] = param_results.std()
            bresults[param] = param_results

        return meansourceparams, stddevsourceparams, bresults

    def get_misfit_for_source( self, source ):
        """Calculate misfit for given source"""

        self.set_source( source )
        misfit = float(self.do("get_global_misfit"))
        return misfit

    def get_misfit_for_sources( self, sources, misfit_setup_functions=None ):
        """Calculate misfit for many sources"""

        results = []
        for source in sources:
            self.set_source( source )
            if misfit_setup_functions:
                for setup in misfit_setup_functions:
                    setup()
                    misfit = float(self.do_get_global_misfit())
                    results.append( ((source,setup), misfit) )
            else:
                misfit = float(self.do_get_global_misfit())
                results.append( (source, misfit) )

        return results

    def brute_force_minimize(self,
                    base_source,
                    param_ranges,
                    source_constraints=None,
                    prefunc=None,
                    misfit_method='l2norm',
                    bootstrap_iterations=1000,
                    histogram_filename='histogram-%s.pdf',
                    result_filename='results-%s.txt',
                    misfit_filename='misfit-%s.txt'):

        sourcetype = base_source.sourcetype
        params = []
        pmin = {}
        pmax = {}
        pinc = {}
        plen = {}
        for param, vmin, vmax, vinc in param_ranges:
            params.append( param )
            pmin[param] = vmin
            pmax[param] = vmax
            pinc[param] = vinc
            plen[param] = int(round((vmax-vmin)/vinc))+1

        self.do_set_misfit_method( misfit_method )

        if prefunc is not None: base_source=prefunc(self,base_source)
        base_misfit = self.global_misfit(
                          *self.get_misfits_for_source(base_source)[1:] )

        sources = base_source.make_source_grid(
            param_ranges,
            source_constraints=source_constraints)

        results = self.get_misfits_for_sources( sources, prefunc=prefunc )

        opt_source, misfit = self.get_grid_minimum( results )
        worst_source, worst_misfit = self.get_grid_maximum( results )

        mean, stddev, bresults = self.grid_bootstrap( results,
                                                        bootstrap_iterations,
                                                        params )

        for param in params:

            # plot histogram
            count = {}
            for r in bresults[param]:
                count[r] = count.get(r,0) + 1

            counted = ([],[])
            for r in sorted(count.keys()):
                counted[0].append(r)
                counted[1].append(float(count[r]) / bootstrap_iterations)

            unit = source_infos(sourcetype)[param].unit
            lab.clf()
            lab.bar(counted[0], counted[1],
                facecolor=(0.74,0.81,0.74),
                align='center',
                width=pinc[param]*0.7 )
            lab.xlabel('%s [%s]' % (param,unit) )
            lab.ylabel('probability of result')
            lab.xlim(pmin[param]-pinc[param]/2.,pmax[param]+pinc[param]/2.)
            fn = histogram_filename % param
            lab.savefig(fn)
            dump( [counted[0], counted[1], param, unit, pmin[param], pmax[param], pinc[param] ], fn+'.data' )

            fout = open(result_filename % param, 'w')
            fout.write( '--- Grid search and bootstrap results for parameter: %s ---\n' % param )
            fout.write( 'Best value (using all stations):         %g\n' % opt_source.params[param] )
            fout.write( 'Mean value of bootstrap results:         %g\n' % mean[param] )
            fout.write( 'Standard deviation of bootstrap results: %g\n' % stddev[param] )
            fout.write( 'Misfit method used:                      %s\n' % misfit_method )
            fout.write( 'Base misfit value:                       %g\n' % base_misfit )
            fout.write( 'Worst misfit value:                      %g\n' % worst_misfit )
            fout.write( 'Best misfit value:                       %g\n' % misfit )
            fout.write( 'Base source:  %s\n' % str(base_source) )
            fout.write( 'Worst source: %s\n' % str(worst_source) )
            fout.write( 'Best source:  %s\n' % str(opt_source) )
            fout.close()

            # plot misfit cross-section
            aparams = num.zeros(len(results), dtype=num.float)
            amisfits = num.zeros(len(results), dtype=num.float)
            for i, (source, misfits, norms) in enumerate(results):
                aparams[i] = source.params[param]
                amisfits[i] = self.global_misfit( misfits, norms )

            lab.clf()
            xunit = source_infos(sourcetype)[param].unit
            if xunit == 'm':
                xunit = 'km'
                aparams /= 1000.

            lab.scatter( aparams, amisfits, facecolor=(0.74,0.81,0.74) )
            lab.xlabel('%s [%s]' % (param,xunit) )
            lab.ylabel('norm. misfit')
            fn = misfit_filename % param
            lab.savefig(fn)

            dump( [aparams, amisfits, param, xunit], fn+'.data' )

        for ixparam in xrange(len(params)):
            for iyparam in xrange(ixparam+1,len(params)):

                # plot 2D histograms
                xparam = params[ixparam]
                yparam = params[iyparam]

                xunit = source_infos(sourcetype)[xparam].unit
                yunit = source_infos(sourcetype)[yparam].unit
                xunitfactor = 1.
                yunitfactor = 1.
                if xunit == 'm':
                    xunit = 'km'
                    xunitfactor = 1000.

                if yunit == 'm':
                    yunit = 'km'
                    yunitfactor = 1000.

                count = {}

                for r,s in zip(bresults[xparam],bresults[yparam]):
                    count[(r,s)] = count.get((r,s),0) + 1

                counted = ([],[],[])
                for (r,s) in sorted(count.keys()):
                    counted[0].append(r)
                    counted[1].append(s)
                    counted[2].append(count[(r,s)])

                lab.clf()
                lab.scatter(num.array(counted[0])/xunitfactor, num.array(counted[1])/yunitfactor, counted[2], facecolor=(0.74,0.81,0.74))
                lab.xlabel('%s [%s]' % (xparam,xunit) )
                lab.ylabel('%s [%s]' % (yparam,yunit) )
                fn = histogram_filename % '-'.join((xparam,yparam))
                lab.savefig(fn)
                dump( [counted, xunitfactor, yunitfactor, xparam, xunit, yparam, yunit ], fn+'.data' )

                # plot 2D misfit cross-section
                axparams = num.linspace(pmin[xparam], pmax[xparam], plen[xparam])
                ayparams = num.linspace(pmin[yparam], pmax[yparam], plen[yparam])

                amisfits = num.zeros((plen[yparam],plen[xparam]),dtype=num.float) - 1.0

                for i, (source, misfits, norms) in enumerate(results):
                    ix = round((source.params[xparam] - pmin[xparam])/pinc[xparam])
                    iy = round((source.params[yparam] - pmin[yparam])/pinc[yparam])
                    if amisfits[iy,ix] < 0.:
                        amisfits[iy,ix] = self.global_misfit( misfits, norms )
                    else:
                        amisfits[iy,ix] = min(amisfits[iy,ix], self.global_misfit( misfits, norms ))

                amisfits_masked = num.ma.masked_less( amisfits, 0.0 )

                lab.clf()
                lab.contourf( axparams/xunitfactor, ayparams/yunitfactor, amisfits_masked )
                lab.xlabel('%s [%s]' % (xparam,xunit) )
                lab.ylabel('%s [%s]' % (yparam,yunit) )
                fn = misfit_filename % '-'.join((xparam,yparam))
                lab.savefig(fn)

                dump( [axparams, ayparams, amisfits_masked, xunitfactor, yunitfactor, xparam, xunit, yparam, yunit ], fn+'.data' )

        return opt_source, misfit

    def get_distazi( self ):
        fn = self.tempdir + '/distances'
        self.do_output_distances( fn )
        f = open(fn,'r')
        distances_deg = []
        distances_m = []
        azimuths = []
        for line in f:
            toks = [ float(x) for x in line.split() ]
            distances_deg.append( toks[0] )
            distances_m.append( toks[1] )
            azimuths.append( toks[2] )
        f.close()
        return distances_deg, distances_m, azimuths

    def make_seismograms( self, sources, fnbase ):
        """Calculate seismograms for many sources."""

        for i,source in enumerate(sources):
            self.set_source( source )
            self.do("output_seismograms", fnbase+"-"+str(i), "table", "synthetics", "plain")

    def show_seismograms_gmt( self, which_set="synthetics", which_processing="plain" ):
        tempfnbase = self.tempdir + "/seismogram"
        outfile = tempfnbase+".pdf"
        self.do("output_seismograms ",tempfnbase, "table", which_set, which_processing)
        fns = self.file_names( tempfnbase, "table" )
        cmd = ["autoplot", " ".join(fns), outfile, "--ydistrib", "--fit", "--width=10" ]
        os.system( " ".join(cmd) )
        return outfile

    def show_seismograms_gmt_new_dgg( self, comp='d', which_processing='plain',
                                    scaling='per-component', addopts=[], plotstyle=None,
                                    timeranges=None, ncols=2, phases=None, filenamepattern=None):

        outfile     = filenamepattern
        tempfnbase1 = self.tempdir + "/reference"
        tempfnbase2 = self.tempdir + "/synthetic"
        dummyfile   = self.tempdir + "/dummy"

        extension = 'table'

        f = open(dummyfile, 'w')
        f.write( '0 0\n')
        f.close()

        if plotstyle is None:
            plotstyle = Config.plotstyle

        nrecs = len(self.components)

        self.do("output_seismograms", tempfnbase1, "table", "references", which_processing)
        self.do("output_seismograms", tempfnbase2, "table", "synthetics", which_processing)

        # layout
        ncomps = len(comp)
        outfiles = []
        nrows = int(math.ceil(nrecs/float(ncols)))
        frames = []
        for irow in xrange(nrows):
            for icol in xrange(ncols):
                for icomp in xrange(ncomps):
                    frames.append( ( "--xnlayout=" + str(1),
                                     "--xilayout=" + str(1),
                                     "--ynlayout=" + str(ncomps),
                                     "--yilayout=" + str(icomp + 1),
                                     "--yspacing=0.1"  ) )

        # gather ranges acording to scaling method
        if scaling == 'per-component':
            yranges = []
            for icomp in xrange(ncomps):
                infiles = []
                infiles.extend(  self.file_names( tempfnbase1, extension, components=comp[icomp] )  )
                infiles.extend(  self.file_names( tempfnbase2, extension, components=comp[icomp] )  )
                ymin, ymax = self.autoplot_get_yrange( infiles, addopts )
                yranges.append( (ymin,ymax) )

        if scaling == 'all-the-same':
            yranges = []
            infiles = []
            infiles.extend(  self.file_names( tempfnbase1, extension, components=comp )  )
            infiles.extend(  self.file_names( tempfnbase2, extension, components=comp )  )
            ymin, ymax = self.autoplot_get_yrange( infiles, addopts )
            yranges.append( (ymin,ymax) )

        if scaling == 'per-receiver':
            yranges = []
            infiles = []
            for irec in xrange(nrecs):
                infiles = []
                infiles.extend(  self.file_names( tempfnbase1, extension, components=comp, receivers=[ irec ] )  )
                infiles.extend(  self.file_names( tempfnbase2, extension, components=comp, receivers=[ irec ] )  )
                ymin, ymax = self.autoplot_get_yrange( infiles, addopts )
                yranges.append( (ymin,ymax) )

        if phases:
            distances_deg, distances_m, azi = self.get_distazi()

        # plot it
        for irec in xrange(nrecs):
            currentoutfile = outfile % (irec+1)
            for icomp in xrange(ncomps):
                if comp[icomp] in self.components[irec]:
                    f1 = tempfnbase1 + "-" + str(irec+1) + "-" + comp[icomp] + "." + extension
                    f2 = tempfnbase2 + "-" + str(irec+1) + "-" + comp[icomp] + "." + extension
                else:
                    f1 = dummyfile
                    f2 = dummyfile

                cmd = [ "autoplot",
                    f1+"["+plotstyle[0]+"]", f2+"["+plotstyle[1]+"]",
                    currentoutfile, "--fit",
                     "--ylabel="+Config.component_names[comp[icomp]],
                    "--yannotevery=0", "--yfunit=off", "--ylabelofs=0.15" ]
                cmd.extend( frames[irec*ncomps+icomp] )
                cmd.extend( addopts )
                yrange = None
                if scaling == 'per-component':
                    yrange = yranges[icomp]
                if scaling == 'all-the-same':
                    yrange = yranges[0]
                if scaling == 'per-receiver':
                    yrange = yranges[irec]
                if yrange:
                    cmd.append("--yrange=%g/%g" % yrange )
                axes = 'eW'
                if icomp == 0:
                    axes = axes + 'S'
                if icomp == ncomps-1:
                    axes = axes + 'n'
                    title = "Receiver %i" % (irec+1)
                    if len(self.receiver_names) == nrecs:
                        title += ': ' + self.receiver_names[irec]
                    cmd.append( "--title='%s'" % title )
                    if phases:
                        anotations = open(currentoutfile+'.anot','w')
                        for phase in phases:
                            arrivaltime = phase(distances_m[irec])
                            if arrivaltime is not None:
                                anotations.write( "%f %f %s\n" % (arrivaltime, yrange[0]/2, phase.name) )
                        anotations.close()

                cmd.append( "--axes="+axes )
                if icomp > 0:
                    cmd.append( "-O" )
                if icomp < ncomps-1:
                    cmd.append( "-K" )

                if icomp == 0:
                    cmd.extend( ["--xlabel=time", "--xunit=s"] )

                if not timeranges is None:
                    cmd.append( "--xrange=%f/%f" % timeranges[irec] )

                call(cmd)

                if (os.path.isfile(currentoutfile+'.anot')):
                    os.remove(currentoutfile+'.anot')

            outfiles.append(currentoutfile)

        cmd = ['pdfjoin', '--outfile', outfile % 'all' ]
        cmd.extend(outfiles)
        call(cmd)

        return outfile

    def autoplot_get_yrange(self, infiles, addopts=[] ):
        cmd = [ "autoplot" ]
        cmd.extend( infiles )
        cmd.extend( [ "nofile.pdf", "--fit", "--printenv" ] )
        cmd.extend( addopts )
        plotenv = os.popen( " ".join(cmd), "r")
        penv = {}
        for line in  plotenv:
            k,v = line.strip().split("=")
            if (re.match( r'[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?$', v ) ):
                v = float(v)
                penv[k] = v

        return penv['ymin'], penv['ymax']

    def file_names( self, basefn, extension, components=None, receivers=None ):
        fns = []

        if receivers is None:
            receivers = range(len(self.components))

        for irec in receivers:
            if components == None:
                components = self.components[irec]
            for comp in components:
                if comp in self.components[irec]:
                    fns.append( basefn + "-" + str(irec+1) + "-" + comp + "." + extension)

        return fns

    def show_grid_results_1d_gmt( self, sourcetype, paramname, traces, plotstyle=None, setupname='rel.misfit' ):

        if plotstyle is None:
            plotstyle = Config.plotstyle

        xyfns = []
        for i, trace in enumerate(traces):
            tempfnbase = self.tempdir + "/xy" + str(i)
            xyfn = tempfnbase+".table"
            file = open(xyfn, "w")

            for point in trace:
                (x, y) = point
                file.write( str(x)+" "+str(y)+"\n" )

            file.close()
            xyfns.append( xyfn + "[" + plotstyle[i] + "]" )

        outfilename = self.tempdir+"/"+paramname+".pdf"

        unit = source_infos(sourcetype)[paramname].unit

        cmd = [ "autoplot", ' '.join(xyfns), outfilename, "--fit",
                "--width=4", "--xlabel="+paramname, "--ylabel="+setupname,
                "--yapproxticks=3" ]

        if unit == 'm':
            cmd.extend( [ "--xfunit=km", "--xexp=3" ] )
        elif unit == 'm/s':
            cmd.extend( [ "--xfunit=km", "--xexp=3" ] )
        else:
            cmd.extend( [ "--xunit="+unit ] )

        os.system( " ".join(cmd) )
        return outfilename

    def make_misfit_cross_sections_1d( self, gfdb, s, steps, params, misfit_setup_functions=None ):

        results = []
        for param in params:
            if Config.verbose >= 1:
                sys.stderr.write( "-- making misfit cross-section for: " + param + "\n" )

            vmin = source_infos(s.sourcetype)[param].soft_min
            vmax = source_infos(s.sourcetype)[param].soft_max

            # exceptions:
            if param == 'depth':
                depth_ext = math.sqrt( 0.25*(s.params['length-a']+s.params['length-b'])**2 +
                                        0.25*(s.params['width'])**2 )

                if vmin < gfdb.firstz + depth_ext:
                    vmin = gfdb.firstz + depth_ext

                if vmax > gfdb.firstz + (gfdb.nz-1) * gfdb.dz - depth_ext:
                    vmax = gfdb.firstz + (gfdb.nz-1) * gfdb.dz - depth_ext

            if param == 'radius':
                vmax_limit = s.params['depth']
                if vmax > vmax_limit:
                    vmax = vmax_limit

            if param == 'length-a':
                vmax_limit = 2.0*math.sqrt(s.params['depth']**2 - 0.25*s.params['width']**2 ) - s.params['length-b']
                if vmax > vmax_limit:
                    vmax = vmax_limit

            if param == 'length-b':
                vmax_limit = 2.0*math.sqrt(s.params['depth']**2 - 0.25*s.params['width']**2 ) - s.params['length-a']
                if vmax > vmax_limit:
                    vmax = vmax_limit

            if param == 'width':
                vmax_limit = 2.0*math.sqrt(s.params['depth']**2 - 0.25*(s.params['length-a'] + s.params['length-b'])**2 )
                if vmax > vmax_limit:
                    vmax = vmax_limit

            if param == 'moment':
                vmax = s.params['moment']*2.
                vmin = s.params['moment']/2.

            vdelta = (vmax - vmin) / float(steps)
            grid = s.make_source_grid( [ ( param, vmin, vmax, vdelta )] )
            results_this_param = self.get_misfit_for_sources( grid, misfit_setup_functions )

            for result in results_this_param:
                if misfit_setup_functions:
                    ((source,setup), misfit) = result
                    results.append( ((param,source,setup), misfit) )
                else:
                    (source, misfit) = result
                    results.append( ((param,source), misfit) )

        return results

# add commands as methods
def gen_do_method(command):
    def func(self,*args):
        return self.do( func.command, *args )

    func.command = command
    return func

for command in Minimizer.commands:
    method = gen_do_method(command)
    setattr( Minimizer, 'do_'+command, method )

class Source:
    def __init__(self, sourcetype='bilateral', sourceparams=None, sourceparams_str=None):

        self.sourcetype = sourcetype
        self.params = {}

        if not sourceparams_str is None:
            sourceparams_float = [ float(s) for s in sourceparams_str.split() ]
            for i,sparam in enumerate(param_names(self.sourcetype)):
                self.params[sparam] = sourceparams_float[i]

        if sourceparams is None:
            sourceparams = {}

        for sparam in param_names(sourcetype):
            self.params[sparam] = source_infos(sourcetype)[sparam].default

        self.set_params( sourceparams )

    def set_params(self, sourceparams):
        for sparam in sourceparams.keys():
            if sparam in param_names(self.sourcetype):
                self.params[sparam] = sourceparams[sparam]
            else:
                raise Exception('invalid source parameter: "%s" for source type: "%s"'% (sparam, self.sourcetype) )

    def __str__(self):

        values = []
        for sparam in param_names(self.sourcetype):
            values.append( self.params[sparam] )

        return self.sourcetype + ' ' + ' '.join( [ str(f) for f in values ] )

    def get_params_as_list(self):
        values = []
        for sparam in param_names(self.sourcetype):
            values.append( self.params[sparam] )
        return values

    def set_params_from_list(self, values):
        for sparam, value in zip( param_names(self.sourcetype), values):
            self.params[sparam] = value

    def make_source_grid( self, sourceparams, irecurs=0, paramlist=None, sourceslist=None, source_constraints=None ):
        """create a grid of sources by """

        if paramlist is None:
            paramlist = [0.] * len(sourceparams)

        if sourceslist is None:
            sourceslist = []

        key, vmin, vmax, vstep = sourceparams[ irecurs ]
        n = int(round((vmax-vmin)/vstep))+1
        for i in xrange(n):
            v = vmin + i * vstep
            paramlist[irecurs] = v
            if irecurs < len(sourceparams)-1:
                self.make_source_grid( sourceparams, irecurs+1, paramlist, sourceslist, source_constraints=source_constraints )
            else:
                v = {}
                for i,elem in enumerate(sourceparams):
                    v[elem[0]] = paramlist[i]

                s = copy.deepcopy(self)
                s.set_params(v)
                if source_constraints == None or source_constraints(s):
                    sourceslist.append(s)

        return sourceslist

    def make_source_randomize( self, sourceparams, nsources ):
        """create a grid of sources by """

        sourceslist = []
        for isource in xrange(nsources):
            v = {}
            for iparam in xrange(len(sourceparams)):
                key, vmin, vmax = sourceparams[ iparam ]
                v[key] = random.uniform(vmin, vmax)

            s = copy.deepcopy(self)
            s.set_params(v)
            sourceslist.append(s)

        return sourceslist

class Phase:
    def __init__(self,name,filename=None):

        self.name = name
        if filename is None:
            filename = os.path.join(kiwi_aux_dir(), 'phases', name)

        f = open(filename,'r')
        self.ref_points = []
        for line in f:
            distance, time = [float(x) for x in line.split()]
            self.ref_points.append( (distance, time) )

    def __call__(self, distance):
        for (low,high) in zip( self.ref_points[0:-1],self.ref_points[1:len(self.ref_points)]):
            if low[0] <= distance <= high[0]:
                return low[1] + (distance-low[0])/(high[0]-low[0])*(high[1]-low[1])
        return None
        #raise Exception("distance %f out of range [%f,%f] for phase %s" %
        #                 (distance, self.ref_points[0][0], self.ref_points[-1][0], self.name) )

class SourceInfo:
    info = None
    info_flat = None

def source_types():
    """Get available source types."""

    if SourceInfo.info is None:
        SourceInfo.info, SourceInfo.info_flat = get_source_infos()

    return SourceInfo.info.keys()

def source_infos( sourcetype ):
        """get some information about possible sources"""

        if SourceInfo.info is None:
            SourceInfo.info, SourceInfo.info_flat = get_source_infos()

        return SourceInfo.info[sourcetype]

def source_infos_flat( sourcetype ):
        """get some information about possible sources"""

        if SourceInfo.info is None:
            SourceInfo.info, SourceInfo.info_flat = get_source_infos()

        return SourceInfo.info_flat[sourcetype]

def param_names( sourcetype ):
        """returns param names in correct order for given source type"""

        if SourceInfo.info is None:
            SourceInfo.info, SourceInfo.info_flat = get_source_infos()

        return SourceInfo.info_flat[sourcetype]["name"]

def get_source_infos():
        """get some information about possible sources in minimizer by asking source_info"""
        info = {}
        info_flat = {}

        class Entry:
            pass

        # get avail. source types
        cmd = [ Config.source_info_prog ]

        source_info = os.popen( ' '.join(cmd), 'r' )

        for line in source_info:
            if re.match(r'\s*source types: ', line):
                sourcetypes = re.sub(r'\s*source types: ','',line).split()

        # get parameter names for source types
        for sourcetype in sourcetypes:
            cmd = [ Config.source_info_prog, sourcetype ]

            source_info = os.popen( ' '.join(cmd), 'r' )

            key_translate = { 'parameter names: ': 'name',
                              'parameter units: ': 'unit',
                              'parameter hard min: ': 'hard_min',
                              'parameter hard max: ': 'hard_max',
                              'parameter soft min: ': 'soft_min',
                              'parameter soft max: ': 'soft_max',
                              'parameter defaults: ': 'default' }

            params_flat = {}
            try:
              for line in source_info:
                  # string fields
                  for key in ['parameter names: ',
                              'parameter units: ']:

                      if re.match(r'\s*'+key, line):
                          pars = re.sub(r'\s*'+key, '', line).split()
                          pkey = key_translate[key]
                          params_flat[pkey] = pars

                  for key in ['parameter hard min: ',
                              'parameter hard max: ',
                              'parameter soft min: ',
                              'parameter soft max: ',
                              'parameter defaults: ' ]:

                      if re.match(r'\s*'+key, line):
                          pars = [ float(s) for s in re.sub(r'\s*'+key, '', line).split() ]
                          pkey = key_translate[key]
                          params_flat[pkey] = pars
            except:
                print 'uuuups'

            source_info.close()
            info_flat[sourcetype] = params_flat

            params = {}
            for i,pname in enumerate(params_flat['name']):
                entry = Entry()
                for key in params_flat.keys():
                    setattr(entry, key, params_flat[key][i])
                params[pname] = entry

            info[sourcetype] = params

        return info, info_flat

def get_gfdb_infos( gfdbpath ):

    class Gfdb:
        pass

    gfdb_infos_str = {}
    cmd = [ Config.gfdb_info_prog, gfdbpath ]
    gfdb_info = os.popen( ' '.join(cmd), 'r' )
    for line in gfdb_info:
        k,v = line.split('=')
        gfdb_infos_str[k] = v.strip()

    gfdb_infos = Gfdb()
    for k in [ 'dt', 'dx', 'dz', 'firstx', 'firstz' ]:
        setattr(gfdb_infos, k, float( gfdb_infos_str[k] ))

    for k in [ 'nchunks', 'nx', 'nz', 'ng' ]:
        setattr(gfdb_infos, k, int( gfdb_infos_str[k] ))

    return gfdb_infos

def table_to_bin(ifn, ofn):
    """make binary file which vtk can understand from tabular output from seismosizer"""

    i = open(ifn)
    o = open(ofn,"w")
    for line in i:
        vals = line.split()
        nvals = len(vals)
        val = vals[-1]
        for ival in range(nvals,4):
            vals.append(val)

        data = pack("ffff", *([float(x) for x in vals[:4]]))
        o.write(data)

def psm_info_to_vtk(infofilename, outfilenamebase):
    """convert PSM info file to VTK format files"""
    sections = set(["center","outline","rupture","slip","eikonal-grid"])
    i = open(infofilename)
    atsec = ''
    points = []
    for line in i:
        sline = line.strip()
        if sline == '':    # at a section end
            if atsec != '':
                ofn = outfilenamebase+ "-" + atsec +".vtk"
                psm_info_to_vtk_section( atsec, points, ofn )
            atsec = ''
            points = []
            continue
        if sline in sections:
            atsec = sline
            continue
        if atsec != '':
            points.append( sline.split() )
    if atsec != '':
        psm_info_to_vtk_section( atsec, points )

    i.close()

def psm_info_to_vtk_section(atsec, points, filename ):
    """called by PsmInfoToVtk() for every section <atsec> in the PSM infofile
        with point data in <points>, this then writes a vtk file for this section."""
    npoints = len(points)
    vtk_head = """# vtk DataFile Version 3.1
generated by minimizer.py
ASCII
DATASET POLYDATA
"""

    vtk_head_ugr = """# vtk DataFile Version 3.1
generated by minimizer.py
ASCII
DATASET UNSTRUCTURED_GRID
"""

    vtk_head_sgr = """# vtk DataFile Version 3.1
generated by minimizer.py
ASCII
DATASET STRUCTURED_GRID
"""

    o = open(filename,"wb")

    def vecstr(v):
        return ' '.join([str(e) for e in v])

    if atsec == "outline":
        o.write(vtk_head)
        o.write("POINTS %i FLOAT\n" % (npoints*2))
        for p in points:
            o.write(vecstr(p[0:3]) + "\n")
        for p in points:
            o.write(vecstr(p[0:2]) + " 0\n")
        o.write("\nPOLYGONS 2 %i\n" % ((npoints+1)*2))
        o.write(str(npoints)+" ")
        o.write(vecstr(range(0,npoints)) + "\n")
        o.write(str(npoints) + " ")
        o.write(vecstr(range(npoints,npoints*2)) + "\n")
        o.write("\nLINES 2 %i\n" % ((npoints+1)*2))
        o.write(str(npoints)+" ")
        o.write(vecstr(range(0,npoints)) + "\n")
        o.write(str(npoints) + " ")
        o.write(vecstr(range(npoints,npoints*2)) + "\n")

    if atsec == "center":
        o.write(vtk_head)
        o.write("POINTS 2 FLOAT\n")
        o.write(vecstr(points[0][0:2]) + " 0\n" )
        o.write(vecstr(points[0]) + "\n")
        o.write("\nLINES 1 3\n")
        o.write("2 0 1\n")

    if atsec == "rupture" or atsec == "slip":
        o.write(vtk_head_ugr)
        o.write("POINTS %i FLOAT\n" % (npoints/2))
        for i in range(0,npoints,2):
            o.write(vecstr(points[i])+"\n")
        o.write("\nPOINT_DATA %i\n" % (npoints/2))
        o.write("VECTORS "+atsec+"vector FLOAT\n")
        for i in xrange(1,npoints,2):
            o.write(vecstr(points[i])+"\n")

    if atsec == "eikonal-grid":
        o.write(vtk_head_sgr)
        gridsize = (int(points[0][0]), int(points[0][1]), 1)
        o.write("DIMENSIONS %i %i %i\n" % gridsize)
        o.write("POINTS %i FLOAT\n" % (gridsize[0]*gridsize[1]))
        for i in xrange(1,npoints):
            o.write(vecstr(points[i][0:3])+"\n")
        o.write("\nPOINT_DATA %i\n" % (gridsize[0]*gridsize[1]))
        o.write("SCALARS rupturetime FLOAT 1\n")
        o.write("LOOKUP_TABLE default\n")

        itimecol = 3
        if len(points[1]) > 4:
            itimecol = 5
        for i in xrange(1,npoints):
            o.write('%s\n' % points[i][itimecol])

    o.close()

max_stars_repo_name: emolch/kiwi

max_issues_repo_name: alirezaniki/kiwi

max_forks_repo_name: emolch/kiwi

reason: import numpy,import scipy


text:
#!/usr/bin/env python3

"""
Example used in the readme. In this example a Bell state is made and then measured.

## Running this script using the "lightweight" CK infrastructure to import QISKit library...

# 1) on local simulator:
    ck virtual env --tags=lib,qiskit --shell_cmd=quantum_coin_flip.py

# 2) on remote simulator (need the API Token from IBM QuantumExperience) :
    CK_IBM_BACKEND=ibmq_qasm_simulator ck virtual `ck search env:* --tags=qiskit,lib`  `ck search env:* --tags=ibmqx,login` --shell_cmd=quantum_coin_flip.py

# 3) on remote quantum hardware (need the API Token from IBM QuantumExperience) :
    CK_IBM_BACKEND=ibmqx4 ck virtual `ck search env:* --tags=qiskit,lib`  `ck search env:* --tags=ibmqx,login` --shell_cmd=quantum_coin_flip.py

"""
import sys
import os

# We don't know from where the user is running the example,
# so we need a relative position from this file path.
# TODO: Relative imports for intra-package imports are highly discouraged.
# http://stackoverflow.com/a/7506006
sys.path.append(os.path.join(os.path.dirname(__file__), '../..'))

from qiskit import QuantumProgram, QISKitError, available_backends, register

try:
    import Qconfig
    register(Qconfig.APItoken, Qconfig.config["url"], verify=False,
                      hub=Qconfig.config["hub"],
                      group=Qconfig.config["group"],
                      project=Qconfig.config["project"])
except:
    print("""
            WARNING: There's no connection with IBMQuantumExperience servers.
            cannot test I/O intesive tasks, will only test CPU intensive tasks
            running the jobs in the local simulator
            """)

available_backends = available_backends()
print("The backends available for use are: {}\n".format(available_backends))
backend = os.environ.get('CK_IBM_BACKEND', 'local_qasm_simulator')

email   = os.environ.get('CK_IBM_API_EMAIL', 'N/A')
print("User email: {}\n".format(email))

timeout = int( os.environ.get('CK_IBM_TIMEOUT', 120) )
shots   = int( os.environ.get('CK_IBM_REPETITION', 10) )

# Create a QuantumProgram object instance.
Q_program = QuantumProgram()

try:
    # Create a Quantum Register called "qr" with 2 qubits.
    qr = Q_program.create_quantum_register("qr", 2)
    # Create a Classical Register called "cr" with 2 bits.
    cr = Q_program.create_classical_register("cr", 2)
    # Create a Quantum Circuit called "qc" with the Quantum Register "qr"
    # and the Classical Register "cr".
    qc = Q_program.create_circuit("bell", [qr], [cr])

    # Add an H gate to qubit 0, putting this qubit in superposition.
    qc.h(qr[0])
    # Add a CX gate to control qubit 0 and target qubit 1, putting
    # the qubits in a Bell state.
    qc.cx(qr[0], qr[1])

    # Add a Measure gate to observe the state.
    qc.measure(qr, cr)

    # Compile and execute the Quantum Program using the given backend.
    result = Q_program.execute(["bell"], backend=backend, shots=shots, seed=1, timeout=timeout)

    # Show the results.
    print(result)       # 'COMPLETED'

    q_execution_time = result.get_data().get('time')
    if q_execution_time:
        print("Quantum execution time: {} sec".format(q_execution_time) )

    print(result.get_data("bell"))

except QISKitError as ex:
    print('Error in the circuit! {}'.format(ex))

########################### Save output to CK format. ##############################

import json
import numpy as np

# See https://stackoverflow.com/questions/26646362/numpy-array-is-not-json-serializable
#
class NumpyEncoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.ndarray):
            return obj.tolist()
        elif isinstance(obj, np.bool_):
            return bool(obj)
        elif isinstance(obj, np.complex):
            return obj.real         # if you care about the imaginary part, try (obj.real, obj.imag)
        return json.JSONEncoder.default(self, obj)

output_dict = {
    'backends': available_backends,
    'email':    email,
    'result':   result.get_data("bell"),
}

formatted_json  = json.dumps(output_dict, cls=NumpyEncoder, sort_keys = True, indent = 4)
#print(formatted_json)

with open('tmp-ck-timer.json', 'w') as json_file:
    json_file.write( formatted_json )

max_stars_repo_name: ctuning/qiskit

max_issues_repo_name: ctuning/qiskit

max_forks_repo_name: ctuning/qiskit


text:
import numpy as np

from gips.gistmodel.fitting import MC_fitter

from gips.gistmodel._numerical_ext import gist_functional_6p_ext
from gips.gistmodel._numerical_ext import gist_functional_5p_ext
from gips.gistmodel._numerical_ext import gist_functional_4p_ext
from gips.gistmodel._numerical_ext import gist_restraint_ext
from gips.gistmodel._numerical_ext import merge_casedata_ext
from gips.gistmodel._numerical_ext import pair_difference_ext
from gips.utils.misc import parms_error

from gips import FLOAT
from gips import DOUBLE

MODE=3

class mode3(MC_fitter):

    def __init__(self, gdatarec_dict,
                    gdata_dict,
                    ref_energy=-11.108,
                    parms=6,
                    pairs=False,
                    radiusadd=[0.,3.],
                    softness=1.0,
                    softcut=2.0,
                    boundsdict=None,
                    pairlist=None,
                    exclude=None,
                    scaling=2.0,
                    select=None,
                    decomp_E=False,
                    decomp_S=False,
                    verbose=False):

        super(mode3, self).__init__(gdatarec_dict=gdatarec_dict,
                                    gdata_dict=gdata_dict,
                                    ref_energy=ref_energy,
                                    mode=MODE,
                                    radiusadd=radiusadd,
                                    softness=softness,
                                    softcut=softcut,
                                    exclude=exclude,
                                    scaling=scaling,
                                    verbose=verbose)

        self.pairs                = pairs
        self.parms                = parms
        self._parms               = parms
        self._gist_functional_ext = None
        self.boundsdict           = boundsdict
        self.pairlist             = pairlist

        if self.pairs:
            self.set_pairs()

        self.set_selection(select)
        self.set_functional()
        self.set_bounds()
        self.set_step()
        self.set_x0()

        self.w      = self.w.astype(DOUBLE)
        self.w_cplx = self.w_cplx.astype(DOUBLE)
        self.w_lig  = self.w_lig.astype(DOUBLE)

    def gist_functional(self, x):

        ### &PyArray_Type, &E,
        ### &PyArray_Type, &S,
        ### &PyArray_Type, &g,
        ### &PyArray_Type, &vol,
        ### &PyArray_Type, &ind,
        ### &PyArray_Type, &x,
        ### &PyArray_Type, &dx,
        ### &PyArray_Type, &fun,
        ### &PyArray_Type, &grad,
        ### &verbose

        ### x[0] = E_aff
        ### x[1] = e_co
        ### x[2] = S_aff
        ### x[3] = s_co
        ### x[4] = g_co
        ### x[5] = C

        _x      = np.zeros(self.parms, dtype=DOUBLE)
        _x[:-1] = x[:-1]

        ### Make sure all types are DOUBLE
        if not self.pairs:
            if not self._gist_functional_ext(self.E, self.S, self.g, self.vol,
                                            self.ind_rec, _x, self._dx, self._calc_data,
                                            self._gradients, 0, int(self.anal_grad)):
                raise ValueError("Something went wrong in gist functional calculation.")

        if not self._gist_functional_ext(self.E_cplx, self.S_cplx, self.g_cplx, self.vol_cplx,
                                        self.ind_rec_cplx, _x, self._dx, self._calc_data_cplx,
                                        self._gradients_cplx, 0, int(self.anal_grad)):
            raise ValueError("Something went wrong in gist functional calculation.")

        if not self._gist_functional_ext(self.E_lig, self.S_lig, self.g_lig, self.vol_lig,
                                        self.ind_rec_lig, _x, self._dx, self._calc_data_lig,
                                        self._gradients_lig, 0, int(self.anal_grad)):
            raise ValueError("Something went wrong in gist functional calculation.")

        ### &PyArray_Type, &x,
        ### &PyArray_Type, &xmin,
        ### &PyArray_Type, &xmax,
        ### &k,
        ### &restraint,
        ### &PyArray_Type, &restraint_grad

        if self.anal_boundary:
            self._restraint = gist_restraint_ext(x,
                                                self.xmin,
                                                self.xmax,
                                                self.kforce_f,
                                                self.kforce,
                                                self._restraint_grad)

    def _f_process(self, x):

        __doc___= """
        returns the squared sum of residuals
        objective function is the free energy
        """

        self.gist_functional(x)

        self._f[:] = 0.
        if self.anal_grad:
            self._g[:] = 0.
            ### Complex and Ligand contributions

            ### &PyArray_Type, &source,
            ### &PyArray_Type, &assign,
            ### &PyArray_Type, &factor,
            ### &PyArray_Type, &assign_factor

            if self.pairs:
                _f          = merge_casedata_ext(self._calc_data_cplx, self.ind_case_cplx, self.w_cplx, self.ind_case_cplx)
                _f         -= merge_casedata_ext(self._calc_data_lig,  self.ind_case_lig,  self.w_lig,  self.ind_case_lig)
                self._f[:] += pair_difference_ext(_f, self.pairidx)

                for i in range(self.parms-1):
                    _g            = merge_casedata_ext(self._gradients_cplx[:,i], self.ind_case_cplx, self.w_cplx, self.ind_case_cplx)
                    _g           -= merge_casedata_ext(self._gradients_lig[:, i], self.ind_case_lig,  self.w_lig,  self.ind_case_lig)
                    self._g[:,i] += pair_difference_ext(_g, self.pairidx)

            else:
                self._f[:]  = merge_casedata_ext(self._calc_data_cplx, self.ind_case_cplx, self.w_cplx, self.ind_case_cplx)
                self._f[:] -= merge_casedata_ext(self._calc_data, self.ind_case, self.w, self.ind_rec)
                self._f[:] -= merge_casedata_ext(self._calc_data_lig,  self.ind_case_lig,  self.w_lig,  self.ind_case_lig)

                for i in range(self.parms):
                    self._g[:,i]  = merge_casedata_ext(self._gradients_cplx[:,i], self.ind_case_cplx, self.w_cplx, self.ind_case_cplx)
                    self._g[:,i] -= merge_casedata_ext(self._gradients[:,i], self.ind_case, self.w, self.ind_rec)
                    self._g[:,i] -= merge_casedata_ext(self._gradients_lig[:, i], self.ind_case_lig,  self.w_lig,  self.ind_case_lig)

            self._f[:]   += x[-1]
            self._g[:,-1] = 1

        else:

            if self.pairs:
                ### Complex and Ligand contributions
                _f          = merge_casedata_ext(self._calc_data_cplx, self.ind_case_cplx, self.w_cplx, self.ind_case_cplx)
                _f         -= merge_casedata_ext(self._calc_data_lig,  self.ind_case_lig,  self.w_lig,  self.ind_case_lig)
                self._f[:] += pair_difference_ext(_f, self.pairidx)

            else:
                ### Receptor contributions
                self._f[:]  = merge_casedata_ext(self._calc_data_cplx, self.ind_case_cplx, self.w_cplx, self.ind_case_cplx)
                self._f[:] -= merge_casedata_ext(self._calc_data, self.ind_case,  self.w, self.ind_rec)
                self._f[:] -= merge_casedata_ext(self._calc_data_lig,  self.ind_case_lig,  self.w_lig,  self.ind_case_lig)

            self._f[:] += x[-1]

    def set_bounds(self):

        __doc__ = """
        Ensures that we don't run out of bounds during MC steps.
        """

        self.xmin = np.zeros(self._parms, dtype=DOUBLE)
        self.xmax = np.zeros(self._parms, dtype=DOUBLE)

        self._restraint_grad = np.zeros(self._parms, dtype=DOUBLE)
        self._restraint      = 0.

        self.kforce_f        = np.zeros(self._parms, dtype=DOUBLE)

        _E = np.min([np.min(self.E_cplx),np.min(self.E),np.min(self.E_lig)]),\
             np.max([np.max(self.E_cplx),np.max(self.E),np.max(self.E_lig)])

        _S = np.min([np.min(self.S_cplx),np.min(self.S),np.min(self.S_lig)]),\
             np.max([np.max(self.S_cplx),np.max(self.S),np.max(self.S_lig)])

        _g = np.min([np.min(self.g_cplx),np.min(self.g),np.min(self.g_lig)]),\
             np.max([np.max(self.g_cplx),np.max(self.g),np.max(self.g_lig)])

        if isinstance(self.boundsdict, dict):
            self.xmin[-1], self.xmax[-1] = self.boundsdict['C'][0], self.boundsdict['C'][1] ### C
        else:
            self.xmin[-1], self.xmax[-1] = -10.          , 10.                              ### C
        self.kforce_f[-1] = 1.

        if self.parms==6:

            if isinstance(self.boundsdict, dict):
                self.xmin[0], self.xmax[0] = self.boundsdict['E'][0],    self.boundsdict['E'][1]    ### E_aff
                self.xmin[1], self.xmax[1] = self.boundsdict['e_co'][0], self.boundsdict['e_co'][1] ### e_co
                self.xmin[2], self.xmax[2] = self.boundsdict['S'][0],    self.boundsdict['S'][1]    ### S_aff
                self.xmin[3], self.xmax[3] = self.boundsdict['s_co'][0], self.boundsdict['s_co'][1] ### s_co
                self.xmin[4], self.xmax[4] = self.boundsdict['g_co'][0], self.boundsdict['g_co'][1] ### g_co
            else:
                self.xmin[0], self.xmax[0] = -10       , 10.        ### E_aff
                self.xmin[1], self.xmax[1] = np.min(_E), np.max(_E) ### e_co
                self.xmin[2], self.xmax[2] = -10.      , 10.        ### S_aff
                self.xmin[3], self.xmax[3] = np.min(_S), np.max(_S) ### s_co
                self.xmin[4], self.xmax[4] = 1.        , np.max(_g) ### g_co
            self.kforce_f[0] = 1.
            self.kforce_f[1] = 10.
            self.kforce_f[2] = 1.
            self.kforce_f[3] = 10.
            self.kforce_f[4] = 10.

        elif self.parms==5:

            if isinstance(self.boundsdict, dict):
                self.xmin[0], self.xmax[0] = self.boundsdict['E'][0],    self.boundsdict['E'][1]    ### Aff
                self.xmin[1], self.xmax[1] = self.boundsdict['e_co'][0], self.boundsdict['e_co'][1] ### e_co
                self.xmin[2], self.xmax[2] = self.boundsdict['s_co'][0], self.boundsdict['s_co'][1] ### s_co
                self.xmin[3], self.xmax[3] = self.boundsdict['g_co'][0], self.boundsdict['g_co'][1] ### g_co
            else:
                self.xmin[0], self.xmax[0] = -10       , 10.        ### Aff
                self.xmin[1], self.xmax[1] = np.min(_E), np.max(_E) ### e_co
                self.xmin[2], self.xmax[2] = np.min(_S), np.max(_S) ### s_co
                self.xmin[3], self.xmax[3] = 1.        , np.max(_g) ### g_co
            self.kforce_f[0] = 1.
            self.kforce_f[1] = 10.
            self.kforce_f[2] = 10.
            self.kforce_f[3] = 10.

        elif self.parms==4:

            if isinstance(self.boundsdict, dict):
                self.xmin[0], self.xmax[0] = self.boundsdict['e_co'][0], self.boundsdict['e_co'][1] ### e_co
                self.xmin[1], self.xmax[1] = self.boundsdict['s_co'][0], self.boundsdict['s_co'][1] ### s_co
                self.xmin[2], self.xmax[2] = self.boundsdict['g_co'][0], self.boundsdict['g_co'][1] ### g_co
            else:
                self.xmin[0], self.xmax[0] = np.min(_E), np.max(_E) ### e_co
                self.xmin[1], self.xmax[1] = np.min(_S), np.max(_S) ### s_co
                self.xmin[2], self.xmax[2] = 1.        , np.max(_g) ### g_co
            self.kforce_f[0] = 10.
            self.kforce_f[1] = 10.
            self.kforce_f[2] = 10.

    def set_step(self):

        self.steps = np.zeros(self._parms, dtype=DOUBLE)

        self.steps[-1] = 1.0

        if self.parms==6:
            self.steps[0] = 1.
            self.steps[1] = 2.0
            self.steps[2] = 1.
            self.steps[3] = 2.0
            self.steps[4] = 2.0

        elif self.parms==5:
            self.steps[0] = 1.
            self.steps[1] = 2.0
            self.steps[2] = 2.0
            self.steps[3] = 2.0

        elif self.parms==4:
            self.steps[0] = 2.0
            self.steps[1] = 2.0
            self.steps[2] = 2.0

        else:
            parms_error(self.parms, self._parms)

    def set_functional(self):

        ### Note, all arrays which are passed to the functionals (such as
        ### gist_functional_6p_ext), must be DOUBLE (i.e. 32bit floating
        ### point type in C). This will not checked within the C routine
        ### (but should be implemented at some point ...).

        if self.pairs:
            self._exp_data = pair_difference_ext(self.dg.astype(DOUBLE), self.pairidx)
            self._f        = np.zeros(self.N_pairs, dtype=DOUBLE)
            self._g        = np.zeros((self.N_pairs, self._parms), dtype=DOUBLE)
        else:
            self._exp_data = np.copy(self.dg.astype(DOUBLE))
            self._f        = np.zeros(self.N_case, dtype=DOUBLE)
            self._g        = np.zeros((self.N_case, self._parms), dtype=DOUBLE)

        self._gradients      = np.zeros((self.N_pos,  self.parms), dtype=DOUBLE)
        self._gradients_cplx = np.zeros((self.N_cplx, self.parms), dtype=DOUBLE)
        self._gradients_lig  = np.zeros((self.N_lig,  self.parms), dtype=DOUBLE)
        self._calc_data      = np.zeros(self.N_pos,  dtype=DOUBLE)
        self._calc_data_cplx = np.zeros(self.N_cplx, dtype=DOUBLE)
        self._calc_data_lig  = np.zeros(self.N_lig,  dtype=DOUBLE)
        self._dx             = 0.00000001

        if self.parms == 6:
            self._gist_functional_ext = gist_functional_6p_ext
        elif self.parms == 5:
            self._gist_functional_ext = gist_functional_5p_ext
        elif self.parms == 4:
            self._gist_functional_ext = gist_functional_4p_ext

        else:
            parms_error(self.parms, self._parms)

max_stars_repo_name: accsc/gips

max_issues_repo_name: accsc/gips

max_forks_repo_name: accsc/gips