text:
from __future__ import print_function
import argparse
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torchvision import datasets, transforms

from pytorch_acdc.layers import StackedLinearACDC

from obj import PyTorchObjective
from scipy import optimize
from tqdm import tqdm

from layers import ShuffleNetLinear, StackedAFDF, GroupedStackedACDC, GroupedStackedAFDF

class Net(nn.Module):
    def __init__(self, original):
        super(Net, self).__init__()
        #self.model = StackedLinearACDC(784, 784, 32, original=original)
        #self.model = nn.Linear(784,10)
        self.model = ShuffleNetLinear(784, 784, 4, 784//4)

    def forward(self, x):
        n, c, h, w = x.size()
        x = x.view(n, c*h*w)
        x = self.model(x)
        return F.log_softmax(x[:,:10], dim=1)

def train(args, device, train_loader):

    # make an objective as a class
    class Objective(nn.Module):
        def __init__(self):
            super(Objective, self).__init__()
            self.model = Net(original=False)
            self.model.train()

        def forward(self):
            data, target = train_loader[0][0], train_loader[0][1]
            data, target = data.to(device), target.to(device)
            self.model.to(device)
            output = self.model(data)
            return F.nll_loss(output, target)

    objective = Objective()
    obj = PyTorchObjective(objective)
    maxiter = 1000
    with tqdm(total=maxiter) as pbar:
        def verbose(*args):
            pbar.update(1)
        xL = optimize.minimize(obj.fun, obj.x0, method='CG', jac=obj.jac,
                callback=verbose, options={'maxiter':maxiter})
        #xL = optimize.minimize(obj.fun, obj.x0, method='BFGS', jac=obj.jac,
        #        callback=verbose, options={'maxiter':maxiter})
    # make sure the final value is loaded to the model
    _ = obj.fun(xL.x)
    return obj.f.model

def test(args, model, device, test_loader):
    model.eval()
    test_loss = 0
    correct = 0
    with torch.no_grad():
        for data, target in test_loader:
            data, target = data.to(device), target.to(device)
            output = model(data)
            test_loss += F.nll_loss(output, target, reduction='sum').item() # sum up batch loss
            pred = output.max(1, keepdim=True)[1] # get the index of the max log-probability
            correct += pred.eq(target.view_as(pred)).sum().item()

    test_loss /= len(test_loader.dataset)
    print('\nTest set: Average loss: {:.4f}, Accuracy: {}/{} ({:.0f}%)\n'.format(
        test_loss, correct, len(test_loader.dataset),
        100. * correct / len(test_loader.dataset)))

def main():
    # Training settings
    parser = argparse.ArgumentParser(description='PyTorch MNIST Example')
    parser.add_argument('--batch-size', type=int, default=1000, metavar='N',
                        help='input batch size for training (default: 64)')
    parser.add_argument('--test-batch-size', type=int, default=1000, metavar='N',
                        help='input batch size for testing (default: 1000)')
    parser.add_argument('--epochs', type=int, default=10, metavar='N',
                        help='number of epochs to train (default: 10)')
    parser.add_argument('--lr', type=float, default=0.001, metavar='LR',
                        help='learning rate (default: 0.001)')
    parser.add_argument('--no-cuda', action='store_true', default=False,
                        help='disables CUDA training')
    parser.add_argument('--seed', type=int, default=1, metavar='S',
                        help='random seed (default: 1)')
    parser.add_argument('--log-interval', type=int, default=10, metavar='N',
                        help='how many batches to wait before logging training status')
    parser.add_argument('--original', action='store_true', default=False,
                        help='uses original ACDC parameterisation')
    args = parser.parse_args()
    use_cuda = not args.no_cuda and torch.cuda.is_available()

    torch.manual_seed(args.seed)

    device = torch.device("cuda" if use_cuda else "cpu")

    #kwargs = {'num_workers': 1, 'pin_memory': True} if use_cuda else {}
    torch.manual_seed(1)
    train_loader = torch.utils.data.DataLoader(
        datasets.MNIST('data', train=True, download=True,
                       transform=transforms.Compose([
                           transforms.ToTensor(),
                           transforms.Normalize((0.1307,), (0.3081,))
                       ])),
        batch_size=args.batch_size, shuffle=True)
    train_loader = [next(iter(train_loader))] # train on only one minibatch
    test_loader = torch.utils.data.DataLoader(
        datasets.MNIST('data', train=False, transform=transforms.Compose([
                           transforms.ToTensor(),
                           transforms.Normalize((0.1307,), (0.3081,))
                       ])),
        batch_size=args.test_batch_size, shuffle=True)

    model = train(args, device, train_loader)
    test(args, model, device, test_loader)

if __name__ == '__main__':
    main()

max_stars_repo_name: gngdb/pytorch-acdc

max_issues_repo_name: gngdb/pytorch-acdc

max_forks_repo_name: gngdb/pytorch-acdc

reason: from scipy


text:
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from sklearn import svm
import matplotlib.gridspec as gridspec

def train_linear_svm(X, y, c, sigma):
    svmAux = svm.SVC(kernel= "linear", C=c)
    svmAux.fit(X,y.ravel())
    return svmAux
def train_rbf_svm(X, y, c, sigma):
    svmAux = svm.SVC(kernel= "rbf", C=c, gamma= 1/(2 * sigma**2))
    svmAux.fit(X,y.ravel())
    return svmAux
def create_dictionary():
    return {"linear":train_linear_svm, "rbf":train_rbf_svm}

def visualize_boundary_mod(X, y, svm, fig):
    x1 = np.linspace(X[:, 0].min(), X[:, 0].max(), 100)
    x2 = np.linspace(X[:, 1].min(), X[:, 1].max(), 100)
    x1, x2 = np.meshgrid(x1, x2)
    yp = svm.predict(np.array([x1.ravel(), x2.ravel()]).T).reshape(x1.shape)
    pos = (y == 1).ravel()
    neg = (y == 0).ravel()
    fig.scatter(X[pos, 0], X[pos, 1], color='black', marker='+')
    fig.scatter(
    X[neg, 0], X[neg, 1], color='yellow', edgecolors='black', marker='o')
    fig.contour(x1, x2, yp)

def print_versions(X, y, file_name, C, kern, sigma=None):
    dic = create_dictionary()
    fig = plt.figure(figsize=(10,20))
    gs = gridspec.GridSpec(nrows=5, ncols=2)
    for i in range(C.shape[0]):
        aux = fig.add_subplot(gs[i//2, i%2])
        aux.set_title("C = " + str(C[i]))
        svmAux = dic[kern](X, y, C[i], sigma)
        visualize_boundary_mod(X, y, svmAux, aux)
    fig.savefig(file_name)

def visualize_boundary(X, y, svm, file_name):
    x1 = np.linspace(X[:, 0].min(), X[:, 0].max(), 100)
    x2 = np.linspace(X[:, 1].min(), X[:, 1].max(), 100)
    x1, x2 = np.meshgrid(x1, x2)
    yp = svm.predict(np.array([x1.ravel(), x2.ravel()]).T).reshape(x1.shape)
    pos = (y == 1).ravel()
    neg = (y == 0).ravel()
    plt.figure()
    plt.scatter(X[pos, 0], X[pos, 1], color='black', marker='+')
    plt.scatter(
    X[neg, 0], X[neg, 1], color='yellow', edgecolors='black', marker='o')
    plt.contour(x1, x2, yp)
    plt.savefig(file_name)
    plt.close()

C = np.linspace(0,100,11,dtype=int)
C = np.arange(0, 50, 5)
C[-1] = 100
C[0] = 1

data1 = loadmat("ex6data1.mat")
X1 = data1['X']
Y1 = data1['y']
print_versions(X1, Y1, "data1", C, "linear")

data2 = loadmat("ex6data2.mat")
X2 = data2['X']
Y2 = data2['y']
print_versions(X2, Y2, "data2", C, "rbf", 0.1)

values = np.array([0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30])
data3 = loadmat("ex6data3.mat")
X3 = data3['X']
Y3 = data3['y']
X3val = data3['Xval']
Y3val = data3['yval']
m = X3val.shape[0]

pairs = np.empty((64,2))
errors = np.empty(64)
for i in range(8):
    for j in range(8):
        new = np.array([values[i],values[j]])
        pairs[8*i + j] = new
        aux = train_rbf_svm(X3, Y3, new[0], new[1])
        H = aux.predict(X3val)
        error = np.sum((H - Y3val.ravel())**2)*(1/(2*m))
        errors[8*i + j] = error
print(errors)
opt = pairs[np.argmin(errors)]
print(opt)
aux = train_rbf_svm(X3, Y3, opt[0], opt[1])
visualize_boundary(X3, Y3, aux, "data3")

ext: py

reason: import numpy,from scipy


text:
#!/usr/bin/env python

"""
    Utility to compute optical properties along a sampled file.

    Ed Nowottnick, January, 2015.

    Modified to return lidar intensive properties
    P. Castellanos June 2017

"""

import os
import MieObs_
from types import *
from netCDF4 import Dataset
from mieobs import VNAMES, getAOPext, getAOPscalar, getAOPint
from numpy import zeros, arange, array, ones, zeros, interp, isnan, ma, NaN, squeeze, transpose, shape, asarray
from math import pi, sin, cos, asin, acos

from types           import *
from optparse        import OptionParser
from datetime        import datetime, timedelta
from dateutil.parser import parse         as isoparser
from csv             import DictReader

from MAPL           import Config, eta
from MAPL.constants import *
from pyobs          import ICARTT
from pyobs.sgp4     import getTrack

# Generic Lists of Varnames and Units
VNAMES_DU = ['DU001','DU002','DU003','DU004','DU005']
VNAMES_SS = ['SS001','SS002','SS003','SS004','SS005']
VNAMES_BC = ['BCPHOBIC','BCPHILIC']
VNAMES_OC = ['OCPHOBIC','OCPHILIC']
VNAMES_SU = ['SO4']
MieVarsNames = ['ext','scatext','backscat','aback_sfc','aback_toa','depol','ext2back','tau','ssa','g']
MieVarsUnits = ['km-1','km-1','km-1 sr-1','sr-1','sr-1','unitless','sr','unitless','unitless','unitless']
MieVarsLongNames = ['total aerosol extinction','scattering extinction','total aerosol backscatter',
                    'attenuated aerosol backscatter at the surface','attenuated aerosol backscatter at TOA',
                    'depolarization ratio','extinction to backscatter ratio','aerosol optical depth',
                    'single scattering albedo','assymetry parameter']

IntVarsUnits = {'area':'m2 m-3',
                'vol': 'm3 m-3',
                'refi':'unitless',
                'refr':'unitless',
                'reff':'m'}

IntVarsLongNames = {'area':'aerosol cross sectional area',
                    'vol' :'aerosol volume',
                    'refi':'imaginary refractive index',
                    'refr':'real refractive index',
                    'reff':'aerosol effective raidus'}

class MieCalc(object):
    pass
    """
    Generic container for Variables
    """

#---
def getVars(inFile):
    """
#    Parse input file, create variable dictionary
#    """

    Vars       = MieCalc()
    file       = Dataset(inFile)
    names      = file.variables.keys()
    MIENAMES   = names
    for n, name in enumerate(MIENAMES):
        var = file.variables[name]
        if (name == 'trjLon') or (name == 'stnLon'):
            name = 'LONGITUDE'
        if (name == 'trjLat') or (name == 'stnLat'):
            name = 'LATITUDE'
        name = name.upper()
        size = len(var.shape)
        if size == 3:
            setattr(Vars,name,var[:,:,:])
        if size == 2:
            setattr(Vars,name,var[:,:])
        if size == 1:
            setattr(Vars,name,var[:])
    return Vars

#---
def computeMie(Vars, channel, varnames, rcFile, options):
    """
#    Computes optical quantities and combines into a dictionary
#   """

    #STN Sampled?
    if options.station:
        NAMES = VNAMES + ['PS','DELP','RH','AIRDENS']
        nstn = len(Vars.STATION)
        nobs = len(Vars.TIME)

        for v in range(nstn):
            VarsIn = MieCalc()
            for n, name in enumerate(NAMES):
                Var = getattr(Vars,name)
                size = len(Var.shape)
                if (size == 2):
                    #1D Variables, ex. PS
                    setattr(VarsIn,name,Var[v,:])
                if size == 3:
                    #2D Variables, ex. DU001
                    setattr(VarsIn,name,Var[v,:,:])

            if (v==0):
                tau,ssa,g = getAOPscalar(VarsIn,channel,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
                ext,sca,backscat,aback_sfc,aback_toa,depol = getAOPext(VarsIn,channel,I=None,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
                ext2back = ones(backscat.shape)*MAPL_UNDEF
                I = backscat > 0
                ext2back[I] = ext[I]/backscat[I]
                MieVars = {"ext":[ext],"scatext":[sca],"backscat":[backscat],"aback_sfc":[aback_sfc],"aback_toa":[aback_toa],"depol":[depol],"ext2back":[ext2back],"tau":[tau],"ssa":[ssa],"g":[g]}

                if options.intensive:
                    vol, area, refr, refi, reff = getAOPint(VarsIn,channel,I=None,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
                    MieVars['vol']  = [vol]
                    MieVars['area'] = [area]
                    MieVars['refr'] = [refr]
                    MieVars['refi'] = [refi]
                    MieVars['reff'] = [reff]
            else:
                tau,ssa,g = getAOPscalar(VarsIn,channel,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
                ext,sca,backscat,aback_sfc,aback_toa,depol = getAOPext(VarsIn,channel,I=None,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
                ext2back = ones(backscat.shape)*MAPL_UNDEF
                I = backscat > 0
                ext2back[I] = ext[I]/backscat[I]
                MieVars['ext'].append(ext)
                MieVars['scatext'].append(sca)
                MieVars['backscat'].append(backscat)
                MieVars['aback_sfc'].append(aback_sfc)
                MieVars['aback_toa'].append(aback_toa)
                MieVars['depol'].append(depol)
                MieVars['ext2back'].append(ext2back)
                MieVars['tau'].append(tau)
                MieVars['ssa'].append(ssa)
                MieVars['g'].append(g)

                if options.intensive:
                    vol, area, refr, refi, reff = getAOPint(VarsIn,channel,I=None,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
                    MieVars['vol'].append(vol)
                    MieVars['area'].append(area)
                    MieVars['refr'].append(refr)
                    MieVars['refi'].append(refi)
                    MieVars['reff'].append(reff)

    #TRJ Sampled?
    else:
        tau,ssa,g = getAOPscalar(Vars,channel,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
        ext,sca,backscat,aback_sfc,aback_toa,depol = getAOPext(Vars,channel,I=None,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
        ext2back = ones(backscat.shape)*MAPL_UNDEF
        I = backscat > 0
        ext2back[I] = ext[I]/backscat[I]
        MieVars = {"ext":[ext],"scatext":[sca],"backscat":[backscat],"aback_sfc":[aback_sfc],"aback_toa":[aback_toa],"depol":[depol],"ext2back":[ext2back],"tau":[tau],"ssa":[ssa],"g":[g]}
        if options.intensive:
            vol, area, refr, refi, reff  = getAOPint(Vars,channel,I=None,vnames=varnames,vtypes=varnames,Verbose=True,rcfile=rcFile)
            MieVars['vol']  = [vol]
            MieVars['area'] = [area]
            MieVars['refr'] = [refr]
            MieVars['refi'] = [refi]
            MieVars['reff'] = [reff]

    return MieVars

#---
def writeNC ( stations, lons, lats, tyme, isotimeIn, MieVars, MieVarsNames, MieVarsLongNames,
              MieVarsUnits, inFile, outFile, options, zlib=False):
    """
    Write a NetCDF file with sampled GEOS-5 variables along the satellite track
    described by (lon,lat,tyme).
    """
    km = 72

    # Open NC file
    # ------------
    nc = Dataset(outFile,'w',format=options.format)

    # Set global attributes
    # ---------------------
    nc.title = 'GEOS-5 Sampled Aerosol Optical Properties File'
    nc.institution = 'NASA/Goddard Space Flight Center'
    nc.source = 'Global Model and Assimilation Office'
    nc.history = 'Created from sampled GEOS-5 collections'
    nc.references = 'n/a'
    nc.comment = 'This file contains sampled GEOS-5 aerosol optical properties.'
    nc.contact = 'Ed Nowottnick <edward.p.nowottnick@nasa.gov>'
    nc.Conventions = 'CF'
    nc.inFile = inFile

    # Create dimensions
    # -----------------
    nt = nc.createDimension('time',len(tyme))
    if options.station:
        ns = nc.createDimension('station',len(stations))
    ls = nc.createDimension('ls',19)
    if km>0:
        nz = nc.createDimension('lev',km)
    x = nc.createDimension('x',1)
    y = nc.createDimension('y',1)

    if options.station:
        # Station names
        # -------------
        stnName_ = nc.createVariable('stnName','S1',('station','ls'),zlib=zlib)
        stnName_.long_name = 'Station Names'
        stnName_.axis = 'e'
        stnName_[:] = stations[:]

    # Coordinate variables
    # --------------------
    time = nc.createVariable('time','i4',('time',),zlib=zlib)
    time.long_name = 'Time'
    t0 = tyme[0]
    isot0 = isotimeIn[0]
    date0 = ''.join(isot0[:10])
    time0 = ''.join(isot0[-8:])
    time.units = 'seconds since '+date0+' '+time0
    time[:] = tyme
    if km > 0: # pressure level not supported yet
        lev = nc.createVariable('lev','f4',('lev',),zlib=zlib)
        lev.long_name = 'Vertical Level'
        lev.units = 'km'
        lev.positive = 'down'
        lev.axis = 'z'
        lev[:] = range(1,km+1)

    # Add fake dimensions for GrADS compatibility
    # -------------------------------------------
    x = nc.createVariable('x','f4',('x',),zlib=zlib)
    x.long_name = 'Fake Longitude for GrADS Compatibility'
    x.units = 'degrees_east'
    x[:] = zeros(1)
    y = nc.createVariable('y','f4',('y',),zlib=zlib)
    y.long_name = 'Fake Latitude for GrADS Compatibility'
    y.units = 'degrees_north'
    y[:] = zeros(1)
    if options.station:
        e = nc.createVariable('station','i4',('station',),zlib=zlib)
        e.long_name = 'Station Ensemble Dimension'
        e.axis = 'e'
        e.grads_dim = 'e'
        e[:] = range(len(stations))

    # Lat/Lon Coordinates
    # ----------------------
    if options.station:
        lon = nc.createVariable('longitude','f4',('station',),zlib=zlib)
        lon.long_name = 'Longitude'
        lon.units = 'degrees_east'
        lon[:] = lons[:]
        lat = nc.createVariable('latitude','f4',('station',),zlib=zlib)
        lat.long_name = 'Latitude'
        lat.units = 'degrees_north'
        lat[:] = lats[:]
    else:
        lon = nc.createVariable('longitude','f4',('time',),zlib=zlib)
        lon.long_name = 'Longitude'
        lon.units = 'degrees_east'
        lon[:] = lons[:]
        lat = nc.createVariable('latitude','f4',('time',),zlib=zlib)
        lat.long_name = 'Latitude'
        lat.units = 'degrees_north'
        lat[:] = lats[:]

    # Time in ISO format if so desired
    # ---------------------------------
    isotime = nc.createVariable('isotime','S1',('time','ls'),zlib=zlib)
    isotime.long_name = 'Time (ISO Format)'
    isotime[:] = isotimeIn[:]

    # Write each variable
    # --------------------------------------------------
    for n, name in enumerate(MieVarsNames):

        var = squeeze(MieVars[name])
        size = len(var.shape)
        if options.station:
            if size == 2:
                var = asarray([var])
                size = len(var.shape)
        if size == 3:
            dim = ('station','time','lev')
        if size == 2:
            dim = ('time','lev')
        if size == 1:
            dim = ('time')
        this = nc.createVariable(name,'f4',dim,zlib=zlib)
        this.standard_name = name
        this.long_name = MieVarsLongNames[n]
        this.units = MieVarsUnits[n]
        this.missing_value = MAPL_UNDEF
        if options.station:
            this[:] = transpose(var,(0,2,1))
        else:
            this[:] = transpose(var)

    if options.intensive:
        for name in IntVarsUnits:
            var = squeeze(MieVars[name])
            size = len(var.shape)
            if options.station:
                if size == 2:
                    var = asarray([var])
                    size = len(var.shape)
            if size == 3:
                dim = ('station','time','lev')
            if size == 2:
                dim = ('time','lev')
            if size == 1:
                dim = ('time')
            this = nc.createVariable(name,'f4',dim,zlib=zlib)
            this.standard_name = name
            this.long_name     = IntVarsLongNames[name]
            this.units = IntVarsUnits[name]
            this.missing_value = MAPL_UNDEF
            if options.station:
                this[:] = transpose(var,(0,2,1))
            else:
                this[:] = transpose(var)

    # Close the file
    # --------------
    nc.close()

    if options.verbose:
        print " <> wrote %s file %s"%(options.format,options.outFile)

#------------------------------------ M A I N ------------------------------------

if __name__ == "__main__":

    format = 'NETCDF3_CLASSIC'
    inFile  = 'trj_sampler.nc'
    outFile = 'ext_sampler.nc'
    channel = (532)
    rcFile = 'Aod3d_532nm.rc'
    intensive = False

#   Parse command line options
#   --------------------------
    parser = OptionParser()

    parser.add_option("-i", "--input", dest="inFile", default=inFile,
              help="Sampled input file")

    parser.add_option("-o", "--output", dest="outFile", default=outFile,
              help="Output file containing optical properties")

    parser.add_option("-r", "--rc", dest="rcFile", default=rcFile,
              help="Resource file pointing to optical tables")

    parser.add_option("-f", "--format", dest="format", default=format,
              help="Output file format: one of NETCDF4, NETCDF4_CLASSIC, NETCDF3_CLASSIC, NETCDF3_64BIT or EXCEL (default=%s)"%format )

    parser.add_option("-c", "--channel", dest="channel", default=channel,
              help="Channel for Mie calculation")

    parser.add_option("-I", "--intensive",default=intensive,
                      action="store_true", dest="intensive",
                      help="return intensive variables")

    parser.add_option("-v", "--verbose",
                      action="store_true", dest="verbose",
                      help="Verbose mode")

    parser.add_option("--du",
                      action="store_true", dest="dust",
                      help="Dust Only")

    parser.add_option("--ss",
                      action="store_true", dest="seasalt",
                      help="Seasalt Only")

    parser.add_option("--su",
                      action="store_true", dest="sulfate",
                      help="Sulfate Only")

    parser.add_option("--bc",
                      action="store_true", dest="bcarbon",
                      help="Black Carbon Only")

    parser.add_option("--oc",
                      action="store_true", dest="ocarbon",
                      help="Organic Carbon Only")

    parser.add_option("--stn",
                      action="store_true", dest="station",
                      help="Input File is from stn_sampler.py")

    (options, args) = parser.parse_args()

    # Create consistent file name extension
    # -------------------------------------
    name, ext = os.path.splitext(options.outFile)
    if ext.upper() == '.XLS':
        options.format = 'EXCEL'
    if 'NETCDF4' in options.format:
        options.outFile = name + '.nc4'
    elif 'NETCDF3' in options.format:
        options.outFile = name + '.nc'
    elif 'EXCEL' in options.format:
        options.outFile = name + '.xls'
    else:
        raise ValueError, 'invalid extension <%s>'%ext

    # Get Variables
    # --------------------------
    Vars = getVars(options.inFile)

    # Run Mie Calculator and Write Output Files
    # --------------------------
    if options.station:
        StnNames = Vars.STNNAME
    else:
        StnNames = ''

    channelIn = float(options.channel)
    MieVars = computeMie(Vars,channelIn,VNAMES,options.rcFile,options)
    writeNC(StnNames,Vars.LONGITUDE,Vars.LATITUDE,Vars.TIME,Vars.ISOTIME,
            MieVars,MieVarsNames,MieVarsLongNames,MieVarsUnits,options.inFile,options.outFile,options)

    if options.dust:
        outFile = options.outFile+'.dust'
        MieVars = computeMie(Vars,channelIn,VNAMES_DU,options.rcFile,options)
        writeNC(StnNames,Vars.LONGITUDE,Vars.LATITUDE,Vars.TIME,Vars.ISOTIME,
                MieVars,MieVarsNames,MieVarsLongNames,MieVarsUnits,options.inFile,outFile,options)

    if options.seasalt:
        outFile = options.outFile+'.ss'
        MieVars = computeMie(Vars,channelIn,VNAMES_SS,options.rcFile,options)
        writeNC(StnNames,Vars.LONGITUDE,Vars.LATITUDE,Vars.TIME,Vars.ISOTIME,
                MieVars,MieVarsNames,MieVarsLongNames,MieVarsUnits,options.inFile,outFile,options)

    if options.sulfate:
        outFile = options.outFile+'.su'
        MieVars = computeMie(Vars,channelIn,VNAMES_SU,options.rcFile,options)
        writeNC(StnNames,Vars.LONGITUDE,Vars.LATITUDE,Vars.TIME,Vars.ISOTIME,
                MieVars,MieVarsNames,MieVarsLongNames,MieVarsUnits,options.inFile,outFile,options)

    if options.bcarbon:
        outFile = options.outFile+'.bc'
        MieVars = computeMie(Vars,channelIn,VNAMES_BC,options.rcFile,options)
        writeNC(StnNames,Vars.LONGITUDE,Vars.LATITUDE,Vars.TIME,Vars.ISOTIME,
                MieVars,MieVarsNames,MieVarsLongNames,MieVarsUnits,options.inFile,outFile,options)

    if options.ocarbon:
        outFile = options.outFile+'.oc'
        MieVars = computeMie(Vars,channelIn,VNAMES_OC,options.rcFile,options)
        writeNC(StnNames,Vars.LONGITUDE,Vars.LATITUDE,Vars.TIME,Vars.ISOTIME,
                MieVars,MieVarsNames,MieVarsLongNames,MieVarsUnits,options.inFile,outFile,options)

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps

reason: from numpy


text:
import numpy as np
from scipy.signal import find_peaks
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy import ndimage

class peakAnalysisTools(object):

    def __init__(self, camera):
        self.camera = camera
        self.grid = plt.GridSpec(2, 2, wspace=0.4, hspace=0.3)
        self.fig = plt.figure()
        self.prominence = 200

    def measureIntensities(self):
        data = {"xAmplitudes": [], "yAmplitudes": []}
        #data['xAmplitudes'] = list(self.summedFunctionx[self.peaksx].astype(int))
        #data['yAmplitudes'] = list(self.summedFunctiony[self.peaksy].astype(int))
        data['xAmplitudes'] = self.summedFunctionx[self.peaksx].tolist()
        data['yAmplitudes'] = self.summedFunctiony[self.peaksy].tolist()
        return data

    def getPlot(self, top, bottom, left, right, rotation, cutSizeX, cutSizeY):

        grayimg = self.camera.get_image()
        # There was an error retrieving image
        if grayimg is None:
            return

        plt.clf()
        grayimg = ndimage.rotate(grayimg, int(rotation))
        grayimg = grayimg[left:right, top:bottom]
        cutx = self.fig.add_subplot(self.grid[0, 0])
        cutx.title.set_text('Channel 1')
        cuty = self.fig.add_subplot(self.grid[1, 0])
        cuty.title.set_text('Channel 0')
        image = self.fig.add_subplot(self.grid[:, 1])
        image.imshow(grayimg, cmap='gray', vmin=0, vmax=255)
        dim = grayimg.shape
        x0 = int(dim[0]/2)
        x1 = int(dim[0]/2)+cutSizeX
        y0 = int(dim[1]/2)
        y1 = int(dim[1]/2)+cutSizeY
        print ("left to right cut: (transpose)" + str(x0) + " " + str(x1))
        print ("top to bottom cut: " + str(y0) + " " + str(y1))
        image.axvline(x0, color='r')
        image.axvline(x1, color='r')
        image.axhline(y0, color='r')
        image.axhline(y1, color='r')
        self.summedFunctionx = np.sum(grayimg[:][y0:y1], axis=0)
        grayimg = np.transpose(grayimg)
        self.summedFunctiony = np.sum(grayimg[x0:x1][:], axis=0)

        self.peaksx, properties = find_peaks(self.summedFunctionx, prominence=(self.prominence, None))
        self.peaksy, properties = find_peaks(self.summedFunctiony, prominence=(self.prominence, None))

        cutx.plot(self.summedFunctionx)
        cutx.plot(self.peaksx, self.summedFunctionx[self.peaksx], "x")
        cuty.plot(self.summedFunctiony)
        cuty.plot(self.peaksy, self.summedFunctiony[self.peaksy], "x")
        cutx.set_ylim([0, 2000])
        cuty.set_ylim([0, 2000])
        return self.fig


text:
import unittest

import numpy as np
import pandas as pd

from dsbox.ml.feature_engineering.timeseries import Shifter, RollingWindower, DistributionTransformer, \
    np_rolling_agg_window, np_rolling_agg_abs_deviation_window
from dsbox.ml.outliers import median_absolute_deviation
from pandas.util.testing import assert_frame_equal, assert_numpy_array_equal, assert_series_equal

class TestComputeShift(unittest.TestCase):
    def test_shift_ts(self):
        # given
        df_ts = pd.DataFrame({'data': [0.0, 1.0, 2.0, 3.0, 4.0]},
                             index=[pd.Timestamp('20130101 09:00:00'),
                                    pd.Timestamp('20130101 09:00:02'),
                                    pd.Timestamp('20130101 09:00:03'),
                                    pd.Timestamp('20130101 09:00:05'),
                                    pd.Timestamp('20130101 09:00:06')])

        # when
        shifter = Shifter(shifts=[1, 2], prefix='t-')
        df_shifted = shifter.transform(df_ts)

        # then
        df_expected = pd.DataFrame(
            {'t-1_data': [np.nan, 0, 1.0, 2.0, 3.0],
             't-2_data': [np.nan, np.nan, 0.0, 1.0, 2.0]},
            index=[pd.Timestamp('20130101 09:00:00'),
                   pd.Timestamp('20130101 09:00:02'),
                   pd.Timestamp('20130101 09:00:03'),
                   pd.Timestamp('20130101 09:00:05'),
                   pd.Timestamp('20130101 09:00:06')]
        )

        self.assertTrue(df_expected.equals(df_shifted))

    def test_shift_several_features(self):
        # given
        df = pd.DataFrame({'x': [0.0, 1.0, 2.0, 3.0, 4.0],
                           'y': [0.0, 0.0, 2.0, 2.0, 1.0]})

        # when
        shifter = Shifter(shifts=[1, 2])
        df_shifted = shifter.transform(df)

        # then
        df_expected = pd.DataFrame(
            {'x_1': [np.nan, 0.0, 1.0, 2.0, 3.0],
             'x_2': [np.nan, np.nan, 0.0, 1.0, 2.0],
             'y_1': [np.nan, 0.0, 0.0, 2.0, 2.0],
             'y_2': [np.nan, np.nan, 0.0, 0.0, 2.0],
             }
        )

        self.assertTrue(df_expected.equals(df_shifted[df_expected.columns]))

    def test_shift_groupby(self):
        # given
        df = pd.DataFrame({'data': [0.0, 1.0, 2.0, 3.0, 4.0],
                           'category': ['A', 'A', 'A', 'B', 'B']})

        # when
        shifter = Shifter(shifts=[1, 2])
        df_shifted = shifter.transform(df.groupby('category'))

        # then
        df_expected = pd.DataFrame(
            {'data_1': [np.nan, 0.0, 1.0, np.nan, 3.0],
             'data_2': [np.nan, np.nan, 0.0, np.nan, np.nan]
             }
        )

        self.assertTrue(df_expected.equals(df_shifted[df_expected.columns]))

    def test_shift_empty_df(self):
        # given
        df = pd.DataFrame()

        # when
        shifter = Shifter()
        df_shifted = shifter.transform(df)

        # then
        df_expected = pd.DataFrame()

        self.assertTrue(df_expected.equals(df_shifted[df_expected.columns]))

class TestRollingWindower(unittest.TestCase):
    def test_compute_moving_average(self):
        df = pd.DataFrame({'data': [0, 1, 2, 3, 4], 'data_bis': [3, 1, 0, 4, 1]})

        roller = RollingWindower(windows=[3], min_periods=1)

        # compute potential differences -> test will fail if a difference is present
        df_mean = pd.DataFrame({'mean_3_data': [0, 0.5, 1, 2, 3], 'mean_3_data_bis': [3, 2, 4 / 3.0, 5 / 3.0, 5 / 3.0]})
        df_test = df_mean.sub(roller.transform(df))
        self.assertEqual(df_test.sum().sum(), 0)

    def test_compute_custom_function(self):
        df = pd.DataFrame({'data': [0, 1, 2, 3, 4], 'data_bis': [3, 1, 0, 4, 1]})

        roller = RollingWindower(operation=median_absolute_deviation, windows=[3], min_periods=1)

        # compute potential differences -> test will fail if a difference is present
        df_mad_true = pd.DataFrame({'mad_3_data': [0, 0.5, 1, 1, 1], 'mad_3_data_bis': [0, 1, 1, 1, 1]})
        df_mad = roller.transform(df)
        df_test = df_mad_true.sub(df_mad)
        self.assertEqual(df_test.sum().sum(), 0)

    def test_compute_moving_average_with_datetime(self):
        df_ts = pd.DataFrame({'data': [0, 1, 2, 3, 4]},
                             index=[pd.Timestamp('20130101 09:00:00'),
                                    pd.Timestamp('20130101 09:00:02'),
                                    pd.Timestamp('20130101 09:00:03'),
                                    pd.Timestamp('20130101 09:00:05'),
                                    pd.Timestamp('20130101 09:00:06')])

        df_mean_ts = pd.DataFrame(
            {'mean_5s_data': [0, 0.5, 1, 2, 2.5],
             'mean_2s_data': [0, 1, 1.5, 3, 3.5]},
            index=[pd.Timestamp('20130101 09:00:00'),
                   pd.Timestamp('20130101 09:00:02'),
                   pd.Timestamp('20130101 09:00:03'),
                   pd.Timestamp('20130101 09:00:05'),
                   pd.Timestamp('20130101 09:00:06')]
        )

        roller = RollingWindower(windows=['5s', '2s'])

        df_test = df_mean_ts.sub(roller.transform(df_ts))
        self.assertEqual(df_test.sum().sum(), 0)

class TestDistribution(unittest.TestCase):
    def test_distribution_transformer_should_produce_df_distribution_per_column(self):
        # given
        df = pd.DataFrame({'sales': [3, 10, 12, 23, 48, 19, 21]})

        # when
        distrib_transformer = DistributionTransformer(3)
        df_distrib = distrib_transformer.fit_transform(df)

        # then
        df_expected = pd.DataFrame({'sales_bin_1': [0, 1, 1, 2, 1, 2, 2],
                                    'sales_bin_2': [0, 0, 0, 0, 1, 0, 0],
                                    'sales_bin_3': [1, 0, 0, 0, 0, 0, 0],
                                    'sales_bin_4': [0, 1, 2, 1, 1, 1, 1]})

        assert_frame_equal(df_expected, df_distrib)

    def test_distribution_transformer_should_produce_df_quantile_per_column(self):
        # given
        df = pd.DataFrame({'sales': [3, 10, 12, 23, 48, 19, 21]})

        # when
        distrib_transformer = DistributionTransformer(3, quantiles=[0., 0.25, 0.5, 0.75, 1.])
        df_distrib = distrib_transformer.fit_transform(df)

        # then
        df_expected = pd.DataFrame({'sales_quantile_1': [3., 3, 3, 10, 12, 19, 19],
                                    'sales_quantile_2': [3, 4.75, 6.5, 11, 17.5, 21, 20],
                                    'sales_quantile_3': [3, 6.5, 10, 12, 23, 23, 21],
                                    'sales_quantile_4': [3, 8.25, 11, 17.5, 35.5, 35.5, 34.5],
                                    'sales_quantile_5': [3., 10, 12, 23, 48, 48, 48],
                                    })

        assert_frame_equal(df_expected, df_distrib)

class TestNumpyRolling(unittest.TestCase):
    def test_np_rolling_mean_window_on_np_array(self):
        # given
        array = np.array([0, 1, 2, 3, 4])

        # when
        array_rol_mean = np_rolling_agg_window(array)

        # then
        array_expected = np.array([np.nan, np.nan, 1, 2, 3])

        assert_numpy_array_equal(array_expected, array_rol_mean)

    def test_np_rolling_median_window_on_np_array(self):
        # given
        array = np.array([0, 1, 2, 3, 0])

        # when
        array_rol_median = np_rolling_agg_window(array, agg_func=np.nanmedian)

        # then
        array_expected = np.array([np.nan, np.nan, 1, 2, 2])

        assert_numpy_array_equal(array_expected, array_rol_median)

    def test_np_rolling_mean_window_on_pd_dataframe(self):
        # given
        df = pd.DataFrame({'values': [0, 1, 2, 3, 4]})

        # when
        df_rol_mean = df.apply(np_rolling_agg_window)

        # then
        df_expected = pd.DataFrame({'values': [np.nan, np.nan, 1, 2, 3]})

        assert_frame_equal(df_expected, df_rol_mean)

    def test_np_rolling_mad_window_on_np_array(self):
        # given
        array = np.array([0, 1, 2, 3, 0])

        # when
        array_rol_mad = np.round(np_rolling_agg_abs_deviation_window(array), 2)

        # then
        array_expected = np.array([np.nan, np.nan, 0.67, 0.67, 1.11])

        assert_numpy_array_equal(array_expected, array_rol_mad)

    def test_np_rolling_mad_window_on_pd_dataframe(self):
        # given
        df = pd.DataFrame({'values': [0, 1, 2, 3, 0]})

        # when
        df_rol_mad = df.apply(np_rolling_agg_abs_deviation_window).round(2)

        # then
        df_expected = pd.DataFrame({'values': [np.nan, np.nan, 0.67, 0.67, 1.11]})

        assert_frame_equal(df_expected, df_rol_mad)

    def test_np_rolling_mad_window_on_pd_dataframe_using_groupby(self):
        # given
        df = pd.DataFrame({
            'group': ['A'] * 5 + ['B'] * 5,
            'values': [0, 1, 2, 3, 4, 0, 1, 2, 3, 0]})

        # when
        serie_rol_group_mad = df.groupby('group')['values'].apply(
            np_rolling_agg_abs_deviation_window).explode().reset_index(drop=True).apply(np.round, decimals=2)

        # then
        serie_expected = pd.Series(name='values',
                                   data=[np.nan, np.nan, 0.67, 0.67, 0.67, np.nan, np.nan, 0.67, 0.67, 1.11])

        assert_series_equal(serie_expected, serie_rol_group_mad)

max_stars_repo_name: Pandinosaurus/dsbox

max_issues_repo_name: Pandinosaurus/dsbox

max_forks_repo_name: Pandinosaurus/dsbox

reason: import numpy


text:
import os
import sys

import numpy as np

import qcdb
sys.path.insert(1, os.path.join(sys.path[0], '..'))
from .utils import *
from .addons import *

#! conventional and density-fitting mp2 test of mp2 itself and setting scs-mp2

def check_mp2(return_value, is_df, is_5050=False):
    if is_df:
        ref       =    -76.0167614256151865
        mp2ss     =     -0.0527406422061238
        mp2os     =     -0.1562926850310142
        mp2corl   =     -0.2090333272371381
        mp2tot    =    -76.2257947528523232
        scscorl   =     -0.2051314361059251
        scstot    =    -76.2218928617211162
        a5050corl =      0.5 * (mp2ss + mp2os)
        a5050tot  =      a5050corl + ref
    else:
        ref       =   -76.01678947133706
        mp2ss     =    -0.05268120425816
        mp2os     =    -0.15637564436589
        mp2corl   =    -0.20905684862405
        mp2tot    =   -76.22584631996111
        scscorl   =    -0.20521117465845
        scstot    =   -76.22200064599551
        a5050corl =     0.5 * (mp2ss + mp2os)
        a5050tot  =     a5050corl + ref

    assert compare_values(ref, qcdb.get_variable('SCF TOTAL ENERGY'), 5, 'mp2 ref')
    assert compare_values(mp2ss, qcdb.get_variable('MP2 SAME-SPIN CORRELATION ENERGY'), 5, 'mp2 ss')
    assert compare_values(mp2os, qcdb.get_variable('MP2 OPPOSITE-SPIN CORRELATION ENERGY'), 5, 'mp2 os')
    assert compare_values(mp2corl, qcdb.get_variable('MP2 CORRELATION ENERGY'), 5, 'mp2 corl')
    assert compare_values(mp2tot, qcdb.get_variable('MP2 TOTAL ENERGY'), 5, 'mp2 tot')
    assert compare_values(scscorl, qcdb.get_variable('SCS-MP2 CORRELATION ENERGY'), 5, 'mp2 scscorl')
    assert compare_values(scstot, qcdb.get_variable('SCS-MP2 TOTAL ENERGY'), 5, 'mp2 scstot')
    if is_5050:
        assert compare_values(a5050corl, qcdb.get_variable('custom SCS-MP2 CORRELATION ENERGY'), 5, 'mp2 scscorl')
        assert compare_values(a5050tot, qcdb.get_variable('custom SCS-MP2 TOTAL ENERGY'), 5, 'mp2 scstot')

    assert compare_values(ref, qcdb.get_variable('CURRENT REFERENCE ENERGY'), 5, 'mp2 ref')
    assert compare_values(mp2corl, qcdb.get_variable('CURRENT CORRELATION ENERGY'), 5, 'mp2 corl')
    assert compare_values(mp2tot, qcdb.get_variable('CURRENT ENERGY'), 5, 'mp2 tot')
    assert compare_values(mp2tot, return_value, 5, 'mp2 return')

def h2o():
    return """
O
H 1 1.0
H 1 1.0 2 90.0
"""

@pytest.fixture
def h2o_y():
    return """
H 0 0 1
O 0 0 0
H 1 0 0
no_com
no_reorient
"""

@pytest.fixture
def h2o_z():
    return """
O 0 0 0
H 0 1 0
H 1 0 0
no_com
no_reorient
"""

dip_z = np.array([   1.52054865,  1.52054865,  0.])
grad_z = np.array([[-0.04159183, -0.04159183, -0.        ],
                   [ 0.01821547,  0.02337636,  0.        ],
                   [ 0.02337636,  0.01821547,  0.        ]])
dip_y = np.array([   1.52054865, 0.,  1.52054865])
grad_y = np.array([[ 0.01821547, 0.,  0.02337636],
                   [-0.04159183, 0., -0.04159183],
                   [ 0.02337636, 0.,  0.01821547]])

@using_cfour
def test_2_conv_mp2(h2o_z):
    qcdb.set_molecule(h2o_z)
    qcdb.set_options({
        'basis': 'cc-pvdz',
        'psi4_mp2_type': 'conv'
    })

    print('   Testing mp2 (conv) ...')
    val = qcdb.energy('c4-mp2')
    check_mp2(val, is_df=False)

@using_cfour
def test_4_conv_scs_mp2(h2o_z):
    qcdb.set_molecule(h2o_z)
    qcdb.set_options({
        'basis': 'cc-pvdz',
        'mp2_os_scale': 1.2,
        'mp2_ss_scale': 0.33333333333333333,
        'psi4_mp2_type': 'conv',
    })

    print('   Testing explicit scs mp2 (conv) ...')
    val = qcdb.energy('c4-mp2')
    check_mp2(val, is_df=False)

@using_cfour
def test_scale(h2o_z):
    h2o_z = qcdb.Molecule(h2o_z)
    qcdb.set_options({
        'basis': 'cc-pvdz',
        'cfour_spin_scal': 'on',
        'cfour_reference': 'uhf',
        'cfour_calc_level': 'mp2',
        'cfour_deriv_level': 'first',
        'CFOUR_DIFF_TYPE': 'relaxed',
    })

    val, jrec = qcdb.energy('c4-cfour', return_wfn=True, molecule=h2o_z)
    dip = np.array([float(jrec['qcvars']['CURRENT DIPOLE {}'.format(i)].data) for i in 'XYZ'])

    tnm = sys._getframe().f_code.co_name
    assert compare_arrays(grad_z, jrec['qcvars']['CURRENT GRADIENT'].data, 5, tnm + ' grad')
    assert compare_arrays(dip_z, dip, 5, tnm + ' dipole')

@using_cfour
def test_scale_2(h2o_y):
    h2o_y = qcdb.Molecule(h2o_y)
    qcdb.set_options({
        'basis': 'cc-pvdz',
        'cfour_spin_scal': 'on',
        'cfour_reference': 'uhf',
        'cfour_calc_level': 'mp2',
        'cfour_deriv_level': 'first',
        'CFOUR_DIFF_TYPE': 'relaxed',
    })

    val, jrec = qcdb.energy('c4-cfour', return_wfn=True, molecule=h2o_y)
    dip = np.array([float(jrec['qcvars']['CURRENT DIPOLE {}'.format(i)].data) for i in 'XYZ'])

    tnm = sys._getframe().f_code.co_name
    assert compare_arrays(grad_y, jrec['qcvars']['CURRENT GRADIENT'].data, 5, tnm + ' grad')
    assert compare_arrays(dip_y, dip, 5, tnm + ' dipole')

@using_cfour
def test_6_conv_custom_scs_mp2(h2o_z):
    qcdb.set_molecule(h2o_z)
    qcdb.set_options({
        'basis': 'cc-pvdz',
        'mp2_os_scale':  0.5,
        'mp2_ss_scale':  0.5,
        'psi4_mp2_type': 'conv',
    })

    print('   Testing user-def scs mp2 (conv) ...')
    val = qcdb.energy('c4-mp2')
    check_mp2(val, is_df=False, is_5050=True)
#    assert False

max_stars_repo_name: nuwandesilva/qcdb

max_issues_repo_name: nuwandesilva/qcdb

max_forks_repo_name: nuwandesilva/qcdb


text:
#!/usr/bin/env python

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

plt.style.use('bmh')
gfs = (6, 6)
fs = 15
c = 2.

def main():
    fig, ax = plt.subplots(figsize=gfs)

    xs = np.linspace(0., 10., 500)
    rvs = stats.lomax(c).pdf(xs)
    ax.plot(xs, rvs, label='$\lambda(x) = c/x$')

    rvs = stats.expon(scale=1/c).pdf(xs)
    ax.plot(xs, rvs, label='$\lambda(x) = c$')

    rvs = stats.rayleigh(scale=np.sqrt(1/c)).pdf(xs)
    ax.plot(xs, rvs, label='$\lambda(x) = cx$')

    ax.set_xlim(0, 10)
    ax.set_ylim(0, 2)
    ax.set_xlabel('$X$ (failure rv)', fontsize=fs)
    ax.set_ylabel('$p(x)$', fontsize=fs)
    ax.tick_params(labelsize=fs)
    ax.legend(fontsize=fs)

    plt.tight_layout()
    plt.savefig('../hf2pdf-examples.png', bbox_inches='tight')
    plt.close()

    # logspace now
    fig, ax = plt.subplots(figsize=gfs)

    xs = np.linspace(0.1, 100., 1000)
    rvs = stats.lomax(c).pdf(xs)
    ax.plot(xs, rvs, label='$\lambda(x) = c/x$')

    rvs = stats.expon(scale=1/c).pdf(xs)
    ax.plot(xs, rvs, label='$\lambda(x) = c$')

    rvs = stats.rayleigh(scale=np.sqrt(1/c)).pdf(xs)
    ax.plot(xs, rvs, label='$\lambda(x) = cx$')

    ax.set_xlabel('$X$ (failure rv)', fontsize=fs)
    ax.set_ylabel('$p(x)$', fontsize=fs)
    ax.tick_params(labelsize=fs)
    ax.legend(fontsize=fs)
    ax.set_ylim(10**-9,)
    ax.set_yscale('log')
    ax.set_xscale('log')

    plt.tight_layout()
    plt.savefig('../hf2pdf-examples-loglog.png', bbox_inches='tight')
    plt.close()

if __name__ == "__main__":
    main()


text:
"""
    Copyright (C) 2015 Rocco Meli

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

from basis import *
from integrals import *

import numpy.linalg as la

def S_overlap(basis):
    """
    Compute overlap matrix S.

    INPUT:
        BASIS: basis set
    OUTPUT:
        S: Overlap matrix
    """

    # Size of the basis set
    K = basis.K

    # List of basis functions
    B = basis.basis()

    S = np.zeros((K,K))

    for i,b1 in enumerate(B):
        for j,b2 in enumerate(B):
            for a1,d1 in zip(b1["a"],b1["d"]):
                for a2,d2 in zip(b2["a"],b2["d"]):
                    R1 = b1["R"]
                    R2 = b2["R"]

                    tmp = d1.conjugate()*d2
                    tmp *= overlap(b1["lx"],b1["ly"],b1["lz"],b2["lx"],b2["ly"],b2["lz"],a1,a2,R1,R2)

                    S[i,j] +=  tmp

    return S

def X_transform(S):
    """
    Compute the transformation matrix X using canonical orthogonalization.

    INPUT:
        S: Overlap matrix
    OUTPUT:
        X: Transformation matrix

    Source:
        Modern Quantum Chemistry
        Szabo and Ostlund
        Dover
        1989
    """

    s, U = la.eig(S)

    s = np.diag(s**(-1./2.))

    X = np.dot(U,s)

    return X

def T_kinetic(basis):
    """
    Compute kinetic matrix T.

    INPUT:
        BASIS: basis set
    OUTPUT:
        T: Kinetic matrix
    """
    # Size of the basis set
    K = basis.K

    # List of basis functions
    B = basis.basis()

    T = np.zeros((K,K))

    for i,b1 in enumerate(B):
        for j,b2 in enumerate(B):
            for a1,d1 in zip(b1["a"],b1["d"]):
                for a2,d2 in zip(b2["a"],b2["d"]):
                    R1 = b1["R"]
                    R2 = b2["R"]

                    tmp = d1.conjugate()*d2
                    tmp *= kinetic(b1["lx"],b1["ly"],b1["lz"],b2["lx"],b2["ly"],b2["lz"],a1,a2,R1,R2)

                    T[i,j] +=  tmp

    return T

def V_nuclear(basis,atom):
    """
    Compute nuclear-electron potential energy matrix Vn.

    INPUT:
        BASIS: basis set
        ATOM: atom specifications (position and charge)
    OUTPUT:
        VN: Nuclear-attraction matrix for atom ATOM
    """
    # Size of the basis set
    K = basis.K

    # List of basis functions
    B = basis.basis()

    # Nuclear coordinates
    Rn = atom.R

    # Nuclear charge
    Zn = atom.Z

    Vn = np.zeros((K,K))

    for i,b1 in enumerate(B):
        for j,b2 in enumerate(B):
            for a1,d1 in zip(b1["a"],b1["d"]):
                for a2,d2 in zip(b2["a"],b2["d"]):
                    R1 = b1["R"]
                    R2 = b2["R"]

                    tmp = d1.conjugate()*d2
                    tmp *= nuclear(b1["lx"],b1["ly"],b1["lz"],b2["lx"],b2["ly"],b2["lz"],a1,a2,R1,R2,Rn,Zn)

                    Vn[i,j] +=  tmp

    return Vn

def H_core(basis,molecule):
    """
    Compute core Hamiltonian (sum of T and all the VN)

    INPUT:
        BASIS: basis set
        MOLECULE: molecule, collection of atom objects
    OUTPUT:
        (T + VN): Core Hamitlonian
    """
    T = T_kinetic(basis)

    print("Kinetic energy")
    print(T)

    # Size of the basis set
    K = basis.K

    Vn = np.zeros((K,K))

    Vnn  = np.zeros((K,K))

    for atom in molecule:
        Vnn = V_nuclear(basis,atom)

        print("Nuclear attraction Vn")
        print(Vnn)

        Vn += Vnn

    print("Total nuclear attraction matrix")
    print(Vn)

    return T + Vn

def P_density(C,N):
    """
    Compute dansity matrix.

    INPUT:
        C: Matrix of coefficients
        N: Number of electrons
    OUTPUT:
        P: density matrix

    Source:
        Modern Quantum Chemistry
        Szabo and Ostlund
        Dover
        1989
    """

    # Size of the basis set
    K = C.shape[0]

    P = np.zeros((K,K))

    for i in range(K):
        for j in range(K):
            for k in range(int(N/2)): #TODO Only for RHF
                P[i,j] += 2 * C[i,k] * C[j,k].conjugate()

    return P

def G_ee(basis,molecule,P,ee):
    """
    Compute core Hamiltonian matrix.

    INPUT:
        BASIS: Basis set.
        MOLECULE: Collection of atoms
        P: Density matrix
        EE: Two-electron integrals
    OUTPUT:
        G: Electron-electron interaction matrix
    """

    # Size of the basis set
    K = basis.K

    G = np.zeros((K,K))

    for i in range(K):
        for j in range(K):
            for k in range(K):
                for l in range(K):
                    G[i,j] += P[k,l] * (ee[i,j,k,l]  - 0.5 * ee[i,l,k,j])

    return G

if __name__ == "__main__":

    """
    Results compared with

        Modern Quantum Chemistry
        Szabo and Ostlund
        Dover
        1989

    and

        The Mathematica Journal
        Evaluation of Gaussian Molecular Integrals
        I. Overlap Integrals
        Minhhuy Hô and Julio Manuel Hernández-Pérez
        2012

    and

        The Mathematica Journal
        Evaluation of Gaussian Molecular Integrals
        II. Kinetic-Energy Integrals
        Minhhuy Hô and Julio Manuel Hernández-Pérez
        2013

    and

        The Mathematica Journal
        Evaluation of Gaussian Molecular Integrals
        III. Nuclear-Electron attraction Integrals
        Minhhuy Hô and Julio Manuel Hernández-Pérez
        2014
    """

    # H2
    H2 = [Atom("H",(0,0,0),1,["1s"]),Atom("H",(0,0,1.4),1,["1s"])]

    # Create the basis set
    sto3g_H2 = STO3G(H2)

    # Compute matrices
    S_H2 = S_overlap(sto3g_H2)
    T_H2 = T_kinetic(sto3g_H2)
    Vn1_H2 = V_nuclear(sto3g_H2,H2[0])
    Vn2_H2 = V_nuclear(sto3g_H2,H2[1])
    H_core_H2 = H_core(sto3g_H2,H2)

    print("###########")
    print("H2 molecule")
    print("###########")

    print("\nOverlap matrix S:")
    print(S_H2)

    print("\nKinetic matrix T:")
    print(T_H2)

    print("\nElectron-nucleus interaction " + H2[0].name + " :")
    print(Vn1_H2)

    print("\nElectron-nucleus interaction " + H2[1].name + " :")
    print(Vn2_H2)

    print("\nCore Hamiltonian:")
    print(H_core_H2)

    # HeH+
    HeH = [Atom("H",(0,0,0),1,["1s"]),Atom("He",(0,0,1.4632),2,["1s"])]

    # Create the basis set
    sto3g_HeH = STO3G(HeH)

    # Compute matrices
    S_HeH = S_overlap(sto3g_HeH)
    T_HeH = T_kinetic(sto3g_HeH)
    Vn1_HeH = V_nuclear(sto3g_HeH,HeH[0])
    Vn2_HeH = V_nuclear(sto3g_HeH,HeH[1])
    H_core_HeH = H_core(sto3g_HeH,HeH)

    print("\n\n\n")
    print("############")
    print("HeH molecule")
    print("############")

    print("\nOverlap matrix S:")
    print(S_HeH)

    print("\nKinetic matrix T:")
    print(T_HeH)

    print("\nElectron-nucleus interaction " + HeH[0].name + " :")
    print(Vn1_HeH)

    print("\nElectron-nucleus interaction " + HeH[1].name + " :")
    print(Vn2_HeH)

    print("\nCore Hamiltonian:")
    print(H_core_HeH)

    # H2O
    H2O = [   Atom("H",(0,+1.43233673,-0.96104039),1,["1s"]),
                Atom("H",(0,-1.43233673,-0.96104039),1,["1s"]),
                Atom("O",(0,0,0.24026010),8,["1s","2s","2p"])]

    sto3g_H2O = STO3G(H2O)

    # Overlap matrix
    S_H2O = S_overlap(sto3g_H2O)
    T_H2O = T_kinetic(sto3g_H2O)
    Vn1_H2O = V_nuclear(sto3g_H2O,H2O[0])
    Vn2_H2O = V_nuclear(sto3g_H2O,H2O[1])
    Vn3_H2O = V_nuclear(sto3g_H2O,H2O[2])

    Vn_H2O = Vn1_H2O + Vn2_H2O + Vn3_H2O

    H_core_H2O = H_core(sto3g_H2O,H2O)

    print("\n\n\n")
    print("############")
    print("H2O molecule")
    print("############")

    print("\nOverlap matrix S:")
    print(S_H2O)

    print("\nKinetic matrix T:")
    print(T_H2O)

    print("\nTotal electron-nucleus interaction:")
    print(Vn_H2O)

    print("\nCore hamiltonian:")
    print(H_core_H2O)

max_stars_repo_name: RMeli/Hartree-Fock

max_issues_repo_name: RMeli/Hartree-Fock

max_forks_repo_name: RMeli/Hartree-Fock


text:
# Test Dingle
import os
import unittest
from unittest.mock import patch
import pandas as pd
import numpy as np
from NN_Predictor import NN_Predictor

class Test_NN_Predictor(unittest.TestCase):

    def setUp(self):
        phony_csv = np.array([["field1,field2,field3,field4,class"],
                             ["0,1,0,1,A"],
                             ["1,0,1,0,B"],
                             ["1,0,1,0,A"],
                             ["0,1,0,1,B"],
                             ["1,0,1,0,C"]])
        self.ersatz_data = pd.DataFrame(phony_csv)

    def test_init(self):
        with self.assertRaises(Exception, msg="A data set must be supplied to"
                               + "create a Predictor."):
            _ = NN_Predictor()
        with self.assertRaises(FileNotFoundError):
            _ = NN_Predictor("Not A File")
        with patch('pandas.read_csv') as mock_reader:
            _ = NN_Predictor("file.csv")
            mock_reader.assert_called_with(f"{os.getcwd()}/file.csv")

    def test_split_data(self):
        print(self.ersatz_data)


text:
#!/usr/bin/env python3

import tensorflow.keras as kr
import numpy as np
import matplotlib.pyplot as plt

def preprocess(images):
	images = images.reshape(images.shape[0], 28, 28, 1)
	return images / 255.0

dataset = kr.datasets.mnist
(train_images, train_labels), (test_images, test_labels) = dataset.load_data()
train_images = preprocess(train_images)
test_images = preprocess(test_images)

model = kr.Sequential([
	kr.layers.Conv2D(16, 5, activation='relu', input_shape=(28, 28, 1)),
	kr.layers.Conv2D(32, 5, activation='relu'),
	kr.layers.MaxPooling2D(pool_size=(2, 2)),
	kr.layers.Flatten(),
	kr.layers.Dense(128, activation='relu'),
	kr.layers.Dense(128, activation='relu'),
	kr.layers.Dropout(0.2),
	kr.layers.Dense(10, activation='softmax')
])
model.compile(optimizer='adam', loss='sparse_categorical_crossentropy',
	metrics=['accuracy'])
model.fit(train_images, train_labels, epochs=5)
model.save('mnist-model.h5')

max_stars_repo_name: callumobrien-bluefruit/tensorflow-stuff

max_issues_repo_name: callumobrien-bluefruit/tensorflow-stuff

max_forks_repo_name: callumobrien-bluefruit/tensorflow-stuff


text:
import argparse
import os

import torch
from pytorch_toolbelt.utils import fs

from xview.dataset import get_datasets
from xview.inference import model_from_checkpoint, run_inference_on_dataset_oof
import numpy as np

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("model", type=str, nargs="+")
    parser.add_argument("-o", "--output-dir", type=str, default=None)
    parser.add_argument("--fast", action="store_true")
    parser.add_argument("--tta", type=str, default=None)
    parser.add_argument("-b", "--batch-size", type=int, default=1, help="Batch Size during training, e.g. -b 64")
    parser.add_argument("-w", "--workers", type=int, default=0, help="")
    parser.add_argument("-dd", "--data-dir", type=str, default="data", help="Data directory")
    parser.add_argument("--size", default=1024, type=int)
    parser.add_argument("--fold", default=None, type=int)
    parser.add_argument("--no-save", action="store_true")
    parser.add_argument("--fp16", action="store_true")
    parser.add_argument("--activation", default="model", type=str)
    parser.add_argument("--align", action="store_true")

    args = parser.parse_args()

    fp16 = args.fp16
    activation = args.activation

    average_score = []
    average_dmg = []
    average_loc = []

    for model_checkpoint in args.model:
        model_checkpoint = fs.auto_file(model_checkpoint)
        checkpoint = torch.load(model_checkpoint)

        print("Model        :", model_checkpoint)
        print(
            "Metrics      :",
            checkpoint["epoch_metrics"]["valid"]["weighted_f1"],
            checkpoint["epoch_metrics"]["valid"]["weighted_f1/localization_f1"],
            checkpoint["epoch_metrics"]["valid"]["weighted_f1/damage_f1"],
        )

        workers = args.workers
        data_dir = args.data_dir
        fast = args.fast
        tta = args.tta
        no_save = args.no_save
        image_size = args.size or checkpoint["checkpoint_data"]["cmd_args"]["size"]
        batch_size = args.batch_size or checkpoint["checkpoint_data"]["cmd_args"]["batch_size"]
        fold = args.fold or checkpoint["checkpoint_data"]["cmd_args"]["fold"]
        align = args.align

        print("Image size :", image_size)
        print("Fold       :", fold)
        print("Align      :", align)
        print("Workers    :", workers)
        print("Save       :", not no_save)

        output_dir = None
        if not no_save:
            output_dir = args.output_dir or os.path.join(
                os.path.dirname(model_checkpoint), fs.id_from_fname(model_checkpoint) + "_oof_predictions"
            )
            print("Output dir :", output_dir)

        # Load models
        model, info = model_from_checkpoint(model_checkpoint, tta=tta, activation_after=None, report=False)
        print(info)
        _, valid_ds, _ = get_datasets(data_dir=data_dir, image_size=(image_size, image_size), fast=fast, fold=fold, align_post=align)

        score, localization_f1, damage_f1, damage_f1s = run_inference_on_dataset_oof(
            model=model,
            dataset=valid_ds,
            output_dir=output_dir,
            batch_size=batch_size,
            workers=workers,
            save=not no_save,
            fp16=fp16
        )

        average_score.append(score)
        average_dmg.append(damage_f1)
        average_loc.append(localization_f1)

        print("Score        :", score)
        print("Localization :", localization_f1)
        print("Damage       :", damage_f1)
        print("Per class    :", damage_f1s)
        print()

    print("Average")
    if len(average_score) > 1:
        print("Score        :", np.mean(average_score), np.std(average_score))
        print("Localization :", np.mean(average_loc), np.std(average_loc))
        print("Damage       :", np.mean(average_dmg), np.std(average_dmg))

if __name__ == "__main__":
    torch.backends.cudnn.benchmark = True
    main()


text:
# Copyright (c) Facebook, Inc. and its affiliates. All Rights Reserved
"""
TridentNet Training Script.

This script is a simplified version of the training script in detectron2/tools.
"""
import argparse
import os
import sys
import torch
import tqdm
import cv2
import numpy as np
sys.path.append('detectron2')

import detectron2.utils.comm as comm
from detectron2.checkpoint import DetectionCheckpointer
from detectron2.data import build_detection_test_loader, build_detection_train_loader
from detectron2.config import get_cfg
from detectron2.engine import DefaultTrainer, default_setup, launch
from detectron2.evaluation import COCOEvaluator, verify_results
from detectron2.structures import Instances

from utils.utils import mkdir, save_features
from utils.extract_utils import get_image_blob, save_bbox, save_roi_features_by_gt_bbox, save_roi_features
from models import add_config
from models.bua.box_regression import BUABoxes

def setup(args):
    """
    Create configs and perform basic setups.
    """
    cfg = get_cfg()
    add_config(args, cfg)
    cfg.merge_from_file(args.config_file)
    cfg.merge_from_list(args.opts)
    cfg.freeze()
    default_setup(cfg, args)
    return cfg

def main():
    parser = argparse.ArgumentParser(description="PyTorch Object Detection2 Inference")
    parser.add_argument(
        "--config-file",
        default="configs/bua-caffe/extract-bua-caffe-r101.yaml",
        metavar="FILE",
        help="path to config file",
    )

    parser.add_argument("--mode", default="caffe", type=str, help="bua_caffe, ...")

    parser.add_argument('--out-dir', dest='output_dir',
                        help='output directory for features',
                        default="features")
    parser.add_argument('--image-dir', dest='image_dir',
                        help='directory with images',
                        default="image")
    parser.add_argument('--gt-bbox-dir', dest='gt_bbox_dir',
                        help='directory with gt-bbox',
                        default="bbox")
    parser.add_argument(
        "--resume",
        action="store_true",
        help="whether to attempt to resume from the checkpoint directory",
    )
    parser.add_argument(
        "opts",
        help="Modify config options using the command-line",
        default=None,
        nargs=argparse.REMAINDER,
    )

    args = parser.parse_args()

    cfg = setup(args)

    MIN_BOXES = cfg.MODEL.BUA.EXTRACTOR.MIN_BOXES
    MAX_BOXES = cfg.MODEL.BUA.EXTRACTOR.MAX_BOXES
    CONF_THRESH = cfg.MODEL.BUA.EXTRACTOR.CONF_THRESH

    model = DefaultTrainer.build_model(cfg)
    DetectionCheckpointer(model, save_dir=cfg.OUTPUT_DIR).resume_or_load(
        cfg.MODEL.WEIGHTS, resume=args.resume
    )
    # Extract features.
    imglist = os.listdir(args.image_dir)
    num_images = len(imglist)
    print('Number of images: {}.'.format(num_images))
    model.eval()

    for im_file in tqdm.tqdm(imglist):
        if os.path.exists(os.path.join(args.output_dir, im_file.split('.')[0]+'.npz')):
            continue
        im = cv2.imread(os.path.join(args.image_dir, im_file))
        if im is None:
            print(os.path.join(args.image_dir, im_file), "is illegal!")
            continue
        dataset_dict = get_image_blob(im, cfg.MODEL.PIXEL_MEAN)
        # extract roi features
        if cfg.MODEL.BUA.EXTRACTOR.MODE == 1:
            attr_scores = None
            with torch.set_grad_enabled(False):
                if cfg.MODEL.BUA.ATTRIBUTE_ON:
                    boxes, scores, features_pooled, attr_scores = model([dataset_dict])
                else:
                    boxes, scores, features_pooled = model([dataset_dict])
            save_roi_features(args, cfg, im_file, im, dataset_dict, boxes, scores, features_pooled, attr_scores)
            # extract bbox only
        elif cfg.MODEL.BUA.EXTRACTOR.MODE == 2:
            with torch.set_grad_enabled(False):
                boxes, scores = model([dataset_dict])
            save_bbox(args, cfg, im_file, im, dataset_dict, boxes, scores)
            # extract roi features by gt bbox
        elif cfg.MODEL.BUA.EXTRACTOR.MODE == 3:
            if not os.path.exists(os.path.join(args.gt_bbox_dir, im_file.split('.')[0]+'.npz')):
                continue
            bbox = torch.from_numpy(np.load(os.path.join(args.gt_bbox_dir, im_file.split('.')[0]+'.npz'))['bbox']) * dataset_dict['im_scale']
            proposals = Instances(dataset_dict['image'].shape[-2:])
            proposals.proposal_boxes = BUABoxes(bbox)
            dataset_dict['proposals'] = proposals

            attr_scores = None
            with torch.set_grad_enabled(False):
                if cfg.MODEL.BUA.ATTRIBUTE_ON:
                    boxes, scores, features_pooled, attr_scores = model([dataset_dict])
                else:
                    boxes, scores, features_pooled = model([dataset_dict])

            save_roi_features_by_gt_bbox(args, cfg, im_file, im, dataset_dict, boxes, scores, features_pooled, attr_scores)

if __name__ == "__main__":
    main()


text:
#!/usr/bin/env python3
# Copyright 2004-present Facebook. All Rights Reserved.

import argparse
import logging
import json
import numpy as np
import os
import trimesh
import csv

import asdf
import asdf.workspace as ws
import glob
import re

def compute_chamfer_distance(chamfer_dist_file):
    chamfer_distance = []
    with open(chamfer_dist_file, newline='') as csvfile:
        spamreader = csv.reader(csvfile, delimiter=' ', quotechar='|')
        for idx, row in enumerate(spamreader):
            if idx>0:
                chamfer_distance.append(float(row[-1]))
    print("avg chamfer distance: ", np.mean(np.array(chamfer_distance)))

def evaluate(experiment_directory, checkpoint, data_dir, mode, specs):

    with open(specs["TestSplit"], "r") as f:
        split = json.load(f)

    chamfer_results = []

    for dataset in split:
        for class_name in split[dataset]:
            if mode=='recon_testset':
                all_names = sorted(glob.glob(os.path.join( experiment_directory, 'Results_recon_testset', checkpoint, 'Meshes', dataset, '*.ply' )))
            elif mode=='recon_testset_ttt':
                all_names = sorted(glob.glob(os.path.join( experiment_directory, 'Results_recon_testset_ttt', checkpoint, 'Meshes', dataset, '*.ply' )))
            elif mode=='generation':
                all_names = sorted(glob.glob(os.path.join( experiment_directory, 'Results_generation', checkpoint, 'Meshes', dataset, '*.ply' )))
            elif mode=='generation_ttt':
                all_names = sorted(glob.glob(os.path.join( experiment_directory, 'Results_generation_ttt', checkpoint, 'Meshes', dataset, '*.ply' )))

            logging.debug(
                'Num of files to be evaluated: {}'.format(len(all_names))
            )
            for instance_name in all_names:
                reconstructed_name = re.split('/', instance_name)[-1]
                if class_name=='laptop':
                    #instance_name = reconstructed_name[:7] + reconstructed_name[-8:-4]
                    if 'generation' in mode:
                        instance_name = glob.glob(os.path.join('data/SurfaceSamples/rbo/laptop/' , reconstructed_name[:4]+'*'+'{:04d}'.format(int(reconstructed_name[5:9])-20)+'.npz'))[0]
                        instance_name = re.split('/', instance_name)[-1][:-4]
                    else:
                        instance_name = reconstructed_name[:-4]
                else:
                    raise Exception("no such class")

                logging.debug(
                    "evaluating " + os.path.join(dataset, class_name, reconstructed_name)
                )

                if mode=='recon_testset':
                    reconstructed_mesh_filename = ws.get_recon_testset_mesh_filename(
                        experiment_directory, checkpoint, dataset, class_name, instance_name
                    )
                elif mode=='recon_testset_ttt':
                    reconstructed_mesh_filename = ws.get_recon_testset_ttt_mesh_filename(
                        experiment_directory, checkpoint, dataset, class_name, instance_name
                    )
                elif mode=='generation':
                    reconstructed_mesh_filename = ws.get_generation_mesh_filename(
                        experiment_directory, checkpoint, dataset, class_name, reconstructed_name[:-4]
                    )
                elif mode=='generation_ttt':
                    reconstructed_mesh_filename = ws.get_generation_ttt_mesh_filename(
                        experiment_directory, checkpoint, dataset, class_name, reconstructed_name[:-4]
                    )

                logging.debug(
                    'reconstructed mesh is "' + reconstructed_mesh_filename + '"'
                )

                ground_truth_samples_filename = os.path.join(
                    data_dir,
                    "SurfaceSamples",
                    dataset,
                    class_name,
                    instance_name + ".npz",
                )

                logging.debug(
                    "ground truth samples are " + ground_truth_samples_filename
                )

                #ground_truth_points = trimesh.load(ground_truth_samples_filename)
                ground_truth_points = np.load(ground_truth_samples_filename)
                ground_truth_points = ground_truth_points["pos"]
                reconstruction = trimesh.load(reconstructed_mesh_filename)

                chamfer_dist = asdf.metrics.chamfer.compute_depth_chamfer(
                    ground_truth_points,
                    reconstruction,
                )

                logging.debug("chamfer distance: " + str(chamfer_dist))

                chamfer_results.append(
                    (os.path.join(dataset, class_name, instance_name), chamfer_dist)
                )

    chamfer_dist_file = os.path.join(ws.get_evaluation_dir(experiment_directory+'/Evaluation_with_Latent/', checkpoint, True), "chamfer.csv")

    with open(os.path.join(chamfer_dist_file),"w") as f:
        f.write("shape, chamfer_dist\n")
        for result in chamfer_results:
            f.write("{}, {}\n".format(result[0], result[1]))

    return chamfer_dist_file

if __name__ == "__main__":

    arg_parser = argparse.ArgumentParser(description="Evaluate a DeepSDF autodecoder")
    arg_parser.add_argument(
        "--experiment",
        "-e",
        dest="experiment_directory",
        required=True,
        help="The experiment directory. This directory should include experiment specifications in "
        + '"specs.json", and logging will be done in this directory as well.',
    )
    arg_parser.add_argument(
        "--checkpoint",
        "-c",
        dest="checkpoint",
        default="latest",
        help="The checkpoint to test.",
    )
    arg_parser.add_argument(
        "--data",
        "-d",
        dest="data_source",
        default="data",
        help="The data source directory.",
    )
    arg_parser.add_argument(
        "--mode",
        "-m",
        required=True,
        help="choose from recon_testset | inter_testset | genration",
    )

    asdf.add_common_args(arg_parser)
    args = arg_parser.parse_args()
    asdf.configure_logging(args)

    specs_filename = os.path.join(args.experiment_directory, "specs.json")
    if not os.path.isfile(specs_filename):
        raise Exception(
            'The experiment directory does not include specifications file "specs.json"'
        )
    specs = json.load(open(specs_filename))

    chamfer_dist_file = evaluate(
        args.experiment_directory,
        args.checkpoint,
        args.data_source,
        args.mode,
        specs,
    )

    compute_chamfer_distance(chamfer_dist_file)

max_stars_repo_name: JitengMu/A-SDF

max_issues_repo_name: JitengMu/A-SDF

max_forks_repo_name: JitengMu/A-SDF


text:
import numpy as np
from scipy.interpolate import interp1d
from IC_functions import particles_to_gadget_massrange

def run(argv):

  # Let's just fix these.
  H0 = 0.1 # km/s/kpc
  OmegaM = 0.3089
  G = 4.3022682e-6 # (km/s)^2 kpc/Msun
  h = 0.6774
  rho = 3*H0**2/(8*np.pi*G) * OmegaM # comoving matter density

  if len(argv) < 3:
    print('python script.py <particle count> <lognormal mass center (in Msun/h),spread (e-folds)>')
    return 1

  NP = int(argv[1])

  mass_info = argv[2].split(',')
  mass_center = float(mass_info[0])
  mass_spread = float(mass_info[1])

  # lognormal masses
  mass = np.exp(np.random.normal(np.log(mass_center),mass_spread,NP))

  # set box size to fit the mass content
  mass_mean = np.mean(mass)
  BoxSize = (np.sum(mass)/rho)**(1./3)

  # uniformly distribute positions
  pos = np.random.uniform(0,BoxSize,(3,NP))

  # zero velocities
  # if they were nonzero, we'd have to divide by sqrt(a)
  vel = np.zeros((3,NP))

  # write to snapshot
  particles_to_gadget_massrange('init.hdf5',pos,vel,mass)

  # these go in the parameter file
  print('Omega0                   %lg'%OmegaM)
  print('OmegaLambda              %lg'%(1-OmegaM))
  print('HubbleParam              %lg'%h) # h
  print('BoxSize                  %lg'%BoxSize)
  print('UnitLength_in_cm         3.085678e21') # 1 kpc/h
  print('UnitMass_in_g            1.989e33') # 1 Msun/h
  print('UnitVelocity_in_cm_per_s 1e5') # 1 km/s
  print('SofteningComovingClass1  %lg     ; or so'%(BoxSize/NP**(1./3) * 0.03))

if __name__ == '__main__':
  from sys import argv
  run(argv)


text:
import cv2
import numpy as np
import pandas as pd

im_num = 50
im_path = r'C:\Users\Lab\GoogleDrive\HPNet\pytorch\logs\ds_015\MODEL_007_BATCH_008_EPOCH_100_NORM\results\{:05d}.png'.format(im_num)
im = cv2.imread(im_path)
h = im.shape[0]  # Height
w = im.shape[1]  # Width
a_ratio = 9 / 16  # Aspect ratio
theta = 35  # Angle of view [deg]

"""
Data shape:
+-----+-----+-----+
|  x  |  y  |  z  |
+-----+-----+-----+
x: Right direction of an image
y: Upward of an image
z:
"""

# Import data as pandas.DataFrame
lbl_path = r'C:\Users\Lab\GoogleDrive\Database\ds_015\val_015.csv'
pred_path = r'C:\Users\Lab\GoogleDrive\HPNet\pytorch\logs\ds_015\MODEL_007_BATCH_008_EPOCH_100_NORM\pred_015.csv'
df_lbl = pd.read_csv(lbl_path)
df_pred = pd.read_csv(pred_path)

x_t, y_t, z_t = df_lbl.iloc[im_num-1, :3]
z_t = - z_t
x_pred, y_pred, z_pred = df_pred.iloc[im_num-1, :3]

def img_w(x_):  # Image width
    global w
    return w * (1 + x_) / 2

def img_h(y_):  # Image height
    global h
    return h * (1 - y_) / 2

r_t = z_t * np.tan(np.deg2rad(theta))
r_pred = z_pred * np.tan(np.deg2rad(theta))

x_t = x_t / r_t
y_t = y_t / (r_t * a_ratio)
x_pred = x_pred / r_pred
y_pred = y_pred / (r_pred * a_ratio)

x_t = img_w(x_t).round().astype(int)
y_t = img_h(y_t).round().astype(int)
x_pred = img_w(x_pred).round().astype(int)
y_pred = img_h(y_pred).round().astype(int)

print(x_t, y_t, x_pred, y_pred)

cv2.drawMarker(im, (x_t, y_t), color=(0, 0, 255), markerType=cv2.MARKER_TILTED_CROSS, thickness=2)
cv2.drawMarker(im, (x_pred, y_pred), color=(0, 255, 0), markerType=cv2.MARKER_TILTED_CROSS, thickness=2)

cv2.imshow('image', im)
cv2.waitKey(0)
cv2.destroyAllWindows()


text:
# -*- coding: utf-8 -*-
#process--->1.load data(X:list of lint,y:int). 2.create session. 3.feed data. 4.training (5.validation) ,(6.prediction)
"""

train the model(transformer) with data enhanced by pre-training of two tasks.
default hyperparameter is d_model=512,h=8,d_k=d_v=64(big). if you have a small data set or want to train a
small model, use d_model=128,h=8,d_k=d_v=16(small), or d_model=64,h=8,d_k=d_v=8(tiny).
"""
#import sys
#reload(sys)
#sys.setdefaultencoding('utf8')
import tensorflow as tf
import numpy as np
from model.transfomer_model import TransformerModel
from data_util_hdf5 import create_or_load_vocabulary,load_data_multilabel,assign_pretrained_word_embedding,set_config
import os
from evaluation_matrix import *
from model.config_transformer import Config
#configuration
FLAGS=tf.app.flags.FLAGS

# you can change as you like
tf.app.flags.DEFINE_boolean("test_mode",False,"whether it is test mode. if it is test mode, only small percentage of data will be used")
tf.app.flags.DEFINE_string("data_path","./data/","path of traning data.")
tf.app.flags.DEFINE_string("training_data_file","./data/bert_train.txt","path of traning data.") #./data/cail2018_bi.json
tf.app.flags.DEFINE_string("valid_data_file","./data/bert_test.txt","path of validation data.")
tf.app.flags.DEFINE_string("test_data_file","./data/bert_test.txt","path of validation data.")
tf.app.flags.DEFINE_integer("d_model", 64, "dimension of model") # 512-->128-->64
tf.app.flags.DEFINE_integer("num_layer", 6, "number of layer")
tf.app.flags.DEFINE_integer("num_header", 8, "number of header")
tf.app.flags.DEFINE_integer("d_k", 8, "dimension of k") # 64-->16-->8
tf.app.flags.DEFINE_integer("d_v", 8, "dimension of v") # 64-->16-->8

# below hyperparameter you can use default one, seldom change
tf.app.flags.DEFINE_string("ckpt_dir","./checkpoint_transformer/","checkpoint location for the model") #save to here, so make it easy to upload for test
tf.app.flags.DEFINE_string("tokenize_style","word","checkpoint location for the model")
tf.app.flags.DEFINE_integer("vocab_size",50002,"maximum vocab size.")
tf.app.flags.DEFINE_float("learning_rate",0.0001,"learning rate") #0.001
tf.app.flags.DEFINE_integer("batch_size", 64, "Batch size for training/evaluating.") # 32-->128
tf.app.flags.DEFINE_integer("decay_steps", 1000, "how many steps before decay learning rate.") # 32-->128
tf.app.flags.DEFINE_float("decay_rate", 1.0, "Rate of decay for learning rate.") #0.65
tf.app.flags.DEFINE_float("dropout_keep_prob", 0.9, "percentage to keep when using dropout.") #0.65
tf.app.flags.DEFINE_integer("sequence_length",200,"max sentence length")#400
tf.app.flags.DEFINE_boolean("is_training",True,"is training.true:tranining,false:testing/inference")
tf.app.flags.DEFINE_integer("num_epochs",30,"number of epochs to run.")
tf.app.flags.DEFINE_integer("process_num",3,"number of cpu used")
tf.app.flags.DEFINE_integer("validate_every", 1, "Validate every validate_every epochs.") #
tf.app.flags.DEFINE_boolean("use_pretrained_embedding",False,"whether to use embedding or not.")#
tf.app.flags.DEFINE_string("word2vec_model_path","./data/Tencent_AILab_ChineseEmbedding_100w.txt","word2vec's vocabulary and vectors") # data/sgns.target.word-word.dynwin5.thr10.neg5.dim300.iter5--->data/news_12g_baidubaike_20g_novel_90g_embedding_64.bin--->sgns.merge.char
tf.app.flags.DEFINE_integer("sequence_length_lm",10,"max sentence length of language model")
tf.app.flags.DEFINE_boolean("is_fine_tuning",False,"is_finetuning.ture:this is fine-tuning stage")

def main(_):
    vocab_word2index, label2index= create_or_load_vocabulary(FLAGS.data_path,FLAGS.training_data_file,FLAGS.vocab_size,test_mode=FLAGS.test_mode,tokenize_style=FLAGS.tokenize_style,model_name='transfomer')
    vocab_size = len(vocab_word2index);print("cnn_model.vocab_size:",vocab_size);num_classes=len(label2index);print("num_classes:",num_classes)
    train,valid, test= load_data_multilabel(FLAGS.data_path,FLAGS.training_data_file,FLAGS.valid_data_file,FLAGS.test_data_file,vocab_word2index,label2index,FLAGS.sequence_length,
                                            process_num=FLAGS.process_num,test_mode=FLAGS.test_mode,tokenize_style=FLAGS.tokenize_style,model_name='transfomer')
    train_X, train_Y= train
    valid_X, valid_Y= valid
    test_X,test_Y = test
    print("Test_mode:",FLAGS.test_mode,";length of training data:",train_X.shape,";valid data:",valid_X.shape,";test data:",test_X.shape,";train_Y:",train_Y.shape)
    # 1.create session.
    gpu_config=tf.ConfigProto()
    gpu_config.gpu_options.allow_growth=True
    with tf.Session(config=gpu_config) as sess:
        #Instantiate Model
        config=set_config(FLAGS,num_classes,vocab_size)
        model=TransformerModel(config)
        #Initialize Save
        saver=tf.train.Saver()
        if os.path.exists(FLAGS.ckpt_dir+"checkpoint"):
            print("Restoring Variables from Checkpoint.")
            saver.restore(sess,tf.train.latest_checkpoint(FLAGS.ckpt_dir))
            #for i in range(2): #decay learning rate if necessary.
            #    print(i,"Going to decay learning rate by half.")
            #    sess.run(model.learning_rate_decay_half_op)
        else:
            print('Initializing Variables')
            sess.run(tf.global_variables_initializer())
            if FLAGS.use_pretrained_embedding:
                vocabulary_index2word={index:word for word,index in vocab_word2index.items()}
                assign_pretrained_word_embedding(sess, vocabulary_index2word, vocab_size,FLAGS.word2vec_model_path,model.embedding,config.d_model) # assign pretrained word embeddings
        curr_epoch=sess.run(model.epoch_step)
        # 2.feed data & training
        number_of_training_data=len(train_X)
        batch_size=FLAGS.batch_size
        iteration=0
        score_best=-100
        f1_score=0
        for epoch in range(curr_epoch,FLAGS.num_epochs):
            loss_total, counter =  0.0, 0
            for start, end in zip(range(0, number_of_training_data, batch_size),range(batch_size, number_of_training_data, batch_size)):
                iteration=iteration+1
                if epoch==0 and counter==0:
                    print("trainX[start:end]:",train_X[start:end],"train_X.shape:",train_X.shape)
                feed_dict = {model.input_x: train_X[start:end],model.input_y:train_Y[start:end],model.dropout_keep_prob: FLAGS.dropout_keep_prob}
                current_loss,lr,l2_loss,_=sess.run([model.loss_val,model.learning_rate,model.l2_loss,model.train_op],feed_dict)
                loss_total,counter=loss_total+current_loss,counter+1
                if counter %30==0:
                    print("Learning rate:%.5f\tLoss:%.3f\tCurrent_loss:%.3f\tL2_loss%.3f\t"%(lr,float(loss_total)/float(counter),current_loss,l2_loss))
                if start!=0 and start%(3000*FLAGS.batch_size)==0:
                    loss_valid, f1_macro_valid, f1_micro_valid= do_eval(sess, model, valid,num_classes,label2index)
                    f1_score_valid=((f1_macro_valid+f1_micro_valid)/2.0)*100.0
                    print("Valid.Epoch %d ValidLoss:%.3f\tF1_score_valid:%.3f\tMacro_f1:%.3f\tMicro_f1:%.3f\t" % (epoch, loss_valid, f1_score_valid, f1_macro_valid, f1_micro_valid))

                    # save model to checkpoint
                    if f1_score_valid>score_best:
                        save_path = FLAGS.ckpt_dir + "model.ckpt"
                        print("going to save check point.")
                        saver.save(sess, save_path, global_step=epoch)
                        score_best=f1_score_valid
            #epoch increment
            print("going to increment epoch counter....")
            sess.run(model.epoch_increment)

            # 4.validation
            print(epoch,FLAGS.validate_every,(epoch % FLAGS.validate_every==0))
            if epoch % FLAGS.validate_every==0:
                loss_valid,f1_macro_valid2,f1_micro_valid2=do_eval(sess,model,valid,num_classes,label2index)
                f1_score_valid2 = ((f1_macro_valid2 + f1_micro_valid2) / 2.0) #* 100.0
                print("Valid.Epoch %d ValidLoss:%.3f\tF1 score:%.3f\tMacro_f1:%.3f\tMicro_f1:%.3f\t"% (epoch,loss_valid,f1_score_valid2,f1_macro_valid2,f1_micro_valid2))
                #save model to checkpoint
                if f1_score_valid2 > score_best:
                    save_path=FLAGS.ckpt_dir+"model.ckpt"
                    print("going to save check point.")
                    saver.save(sess,save_path,global_step=epoch)
                    score_best = f1_score_valid2
            if (epoch == 2 or epoch == 4 or epoch == 6 or epoch == 9 or epoch == 13):
                for i in range(1):
                    print(i, "Going to decay learning rate by half.")
                    sess.run(model.learning_rate_decay_half_op)

        # 5.最后在测试集上做测试，并报告测试准确率 Testto 0.0
        loss_test, f1_macro_test, f1_micro_test=do_eval(sess, model, test,num_classes, label2index)
        f1_score_test=((f1_macro_test + f1_micro_test) / 2.0) #* 100.0
        print("Test.Epoch %d TestLoss:%.3f\tF1_score:%.3f\tMacro_f1:%.3f\tMicro_f1:%.3f\t" % (epoch, loss_test, f1_score_test,f1_macro_test, f1_micro_test))
        print("training completed...")

          #sess,model,valid,iteration,num_classes,label2index
def do_eval(sess,model,valid,num_classes,label2index):
    """
    do evaluation using validation set, and report loss, and f1 score.
    :param sess:
    :param model:
    :param valid:
    :param num_classes:
    :param label2index:
    :return:
    """
    number_examples=valid[0].shape[0]
    valid_x,valid_y=valid
    print("number_examples:",number_examples)
    eval_loss,eval_counter=0.0,0
    batch_size=FLAGS.batch_size
    label_dict=init_label_dict(num_classes)
    eval_macro_f1, eval_micro_f1 = 0.0,0.0
    for start,end in zip(range(0,number_examples,batch_size),range(batch_size,number_examples,batch_size)):
        feed_dict = {model.input_x: valid_x[start:end],model.input_y:valid_y[start:end],model.dropout_keep_prob: 1.0}
        curr_eval_loss, logits= sess.run([model.loss_val,model.logits],feed_dict) # logits：[batch_size,label_size]
        #compute confuse matrix
        label_dict=compute_confuse_matrix_batch(valid_y[start:end],logits,label_dict,name='bright')
        eval_loss=eval_loss+curr_eval_loss
        eval_counter=eval_counter+1
    #compute f1_micro & f1_macro
    f1_micro,f1_macro=compute_micro_macro(label_dict) #label_dict is a dict, key is: an label,value is: (TP,FP,FN). where TP is number of True Positive
    compute_f1_score_write_for_debug(label_dict,label2index)
    return eval_loss/float(eval_counter+small_value),f1_macro,f1_micro

if __name__ == "__main__":
    tf.app.run()

max_stars_repo_name: yyht/bert

max_issues_repo_name: yyht/bert

max_forks_repo_name: yyht/bert


text:
'''
Flockers
=============================================================
A Mesa implementation of Craig Reynolds's Boids flocker model.
Uses numpy arrays to represent vectors.
'''

import numpy as np
import random

from mesa import Model, Agent
from mesa.space import ContinuousSpace
from mesa.time import RandomActivation

class Boid(Agent):
    '''
    A Boid-style flocker agent.

    The agent follows three behaviors to flock:
        - Cohesion: steering towards neighboring agents.
        - Separation: avoiding getting too close to any other agent.
        - Alignment: try to fly in the same direction as the neighbors.

    Boids have a vision that defines the radius in which they look for their
    neighbors to flock with. Their speed (a scalar) and heading (a unit vector)
    define their movement. Separation is their desired minimum distance from
    any other Boid.
    '''
    def __init__(self, unique_id, pos, speed=5, heading=None,
            vision=5, separation=1):
        '''
        Create a new Boid flocker agent.

        Args:
            unique_id: Unique agent identifyer.
            pos: Starting position
            speed: Distance to move per step.
            heading: numpy vector for the Boid's direction of movement.
            vision: Radius to look around for nearby Boids.
            separation: Minimum distance to maintain from other Boids.
        '''
        self.unique_id = unique_id
        self.pos = pos
        self.speed = speed
        if heading is not None:
            self.heading = heading
        else:
            self.heading = np.random.random(2)
            self.heading /= np.linalg.norm(self.heading)
        self.vision = vision
        self.separation = separation

    def cohere(self, neighbors):
        '''
        Return the vector toward the center of mass of the local neighbors.
        '''
        center = np.array([0.0, 0.0])
        for neighbor in neighbors:
            center += np.array(neighbor.pos)
        return center / len(neighbors)

    def separate(self, neighbors):
        '''
        Return a vector away from any neighbors closer than separation dist.
        '''
        my_pos = np.array(self.pos)
        sep_vector = np.array([0, 0])
        for neighbor in neighbors:
            their_pos = np.array(neighbor.pos)
            dist = np.linalg.norm(my_pos - their_pos)
            if dist < self.separation:
                sep_vector -= np.int64(their_pos - my_pos)
        return sep_vector

    def match_heading(self, neighbors):
        '''
        Return a vector of the neighbors' average heading.
        '''
        mean_heading = np.array([0, 0])
        for neighbor in neighbors:
            mean_heading += np.int64(neighbor.heading)
        return mean_heading / len(neighbors)

    def step(self, model):
        '''
        Get the Boid's neighbors, compute the new vector, and move accordingly.
        '''

        neighbors = model.space.get_neighbors(self.pos, self.vision, False)
        if len(neighbors) > 0:
            cohere_vector = self.cohere(neighbors)
            separate_vector = self.separate(neighbors)
            match_heading_vector = self.match_heading(neighbors)
            self.heading += (cohere_vector + separate_vector
                + match_heading_vector)
            self.heading /= np.linalg.norm(self.heading)
        new_pos = np.array(self.pos) + self.heading * self.speed
        new_x, new_y = new_pos
        model.space.move_agent(self, (new_x, new_y))

class BoidModel(Model):
    '''
    Flocker model class. Handles agent creation, placement and scheduling.
    '''
    N = 100
    width = 100
    height = 100

    def __init__(self, N, width, height, speed, vision, separation):
        '''
        Create a new Flockers model.

        Args:
            N: Number of Boids
            width, height: Size of the space.
            speed: How fast should the Boids move.
            vision: How far around should each Boid look for its neighbors
            separtion: What's the minimum distance each Boid will attempt to
                       keep from any other
        '''
        self.N = N
        self.vision = vision
        self.speed = speed
        self.separation = separation
        self.schedule = RandomActivation(self)
        self.space = ContinuousSpace(width, height, True,
            grid_width=10, grid_height=10)
        self.make_agents()
        self.running = True

    def make_agents(self):
        '''
        Create N agents, with random positions and starting headings.
        '''
        for i in range(self.N):
            x = random.random() * self.space.x_max
            y = random.random() * self.space.y_max
            pos = (x, y)
            heading = np.random.random(2) * 2 - np.array((1, 1))
            heading /= np.linalg.norm(heading)
            boid = Boid(i, pos, self.speed, heading, self.vision, self.separation)
            self.space.place_agent(boid, pos)
            self.schedule.add(boid)

    def step(self):
        self.schedule.step()

max_stars_repo_name: joshainglis/mesa

max_issues_repo_name: joshainglis/mesa

max_forks_repo_name: joshainglis/mesa


text:
import math

import mpmath
import scipy.integrate as integrate

def f(x):
 return mpmath.sec(x)

def f_a(x):
 return mpmath.sec(x) * mpmath.tan(x)

def f_b(x):
 return mpmath.cos(x)

def f_c(x):
 return mpmath.sec(x)

def f_d(x):
 return mpmath.log(mpmath.cos(x))

if __name__ == "__main__":
 i = integrate.quad(f, 0, 1)
 print(i)
 print(integrate.quad(f_a, 0, 1))
 print(integrate.quad(f_b, 0, 1))
 print(integrate.quad(f_c, 0, 1))
 print(integrate.quad(f_d, 0, 1))

reason: import scipy,import mpmath