text:
import numpy as np
import os
import scipy.io
import scipy.interpolate
import falco.utils
from falco.utils import _spec_arg
import collections

_influence_function_file = os.path.join(
    os.path.dirname(os.path.abspath(__file__)), "influence_dm5v2.mat")

def falco_gen_dm_poke_cube(dm, mp, dx_dm, flagGenCube=True):
    #Compute sampling of the pupil. Assume that it is square.
    dm.dx_dm = dx_dm
    dm.dx = dx_dm

    # Compute coordinates of original influence function
    Ninf0 = len(dm.inf0) #Number of points across the influence function at its native resolution
    #x_inf0 = (-(Ninf0-1)/2:(Ninf0-1)/2)*dm.dx_inf0 # True for even- or odd-sized influence function maps as long as they are centered on the array.
    x_inf0 = np.arange(-(Ninf0-1)//2, (Ninf0-1)//2)*dm.dx_inf0

    Ndm0 = falco.utils.ceil_even( Ninf0 + (dm.Nact - 1)*(dm.dm_spacing/dm.dx_inf0) ) #Number of points across the DM surface at native influence function resolution
    dm.NdmMin = falco.utils.ceil_even( Ndm0*(dm.dx_inf0/dm.dx))+2 #Number of points across the (un-rotated) DM surface at new, desired resolution.
    dm.Ndm = falco.utils.ceil_even( max(np.abs([np.sqrt(2)*np.cos(np.radians(45-dm.zrot)), np.sqrt(2)*np.sin(np.radians(45-dm.zrot))]))*Ndm0*(dm.dx_inf0/dm.dx))+2 #Number of points across the array to fully contain the DM surface at new, desired resolution and z-rotation angle.

    Xinf0,Yinf0 = np.meshgrid(x_inf0, x_inf0)

    #Compute list of initial actuator center coordinates (in actutor widths).
    dm.Xact, dm.Yact = np.meshgrid(np.arange(dm.Nact)-dm.xc, np.arange(dm.Nact)-dm.yc) # in actuator widths
    x_vec = dm.Xact.reshape(-1)
    y_vec = dm.Yact.reshape(-1)

    dm.NactTotal = len(x_vec) #Total number of actuators in the 2-D array

    sa = np.sin(np.radians(dm.xtilt))
    ca = np.cos(np.radians(dm.xtilt))
    sb = np.sin(np.radians(dm.ytilt))
    cb = np.cos(np.radians(dm.ytilt))
    sg = np.sin(np.radians(-dm.zrot))
    cg = np.cos(np.radians(-dm.zrot))

    Mrot = [[cb * cg, sa * sb * cg - ca * sg, ca * sb * cg + sa * sg, 0.0], \
            [cb * sg, sa * sb * sg + ca * cg, ca * sb * sg - sa * cg, 0.0], \
            [    -sb,                sa * cb,                ca * cb, 0.0], \
            [    0.0,                    0.0,                    0.0, 1.0]]

    dm.xy_cent_act = []
    for iact in range(dm.NactTotal):
        xyzVals = [x_vec[iact], y_vec[iact], 0, 1]
        xyzValsRot = np.dot(Mrot, xyzVals)
        dm.xy_cent_act.append(xyzValsRot[:2])

    dm.xy_cent_act = np.array(dm.xy_cent_act).T[::-1]

    N0 = max(dm.inf0.shape)
    Npad = falco.utils.ceil_odd( np.sqrt(2)*N0 )
    inf0pad = np.zeros((Npad,Npad))

    fromx = int(np.ceil(Npad/2)-np.floor(N0/2))
    tox = int(np.ceil(Npad/2)+np.floor(N0/2)) + 1
    inf0pad[fromx:tox, fromx:tox] = dm.inf0

    ydim,xdim = inf0pad.shape

    xd2  = np.floor(xdim / 2) + 1
    yd2  = np.floor(ydim / 2) + 1
    cx   = np.arange(xdim) - xd2
    cy   = np.arange(ydim) - yd2
    Xs0, Ys0 = np.meshgrid(cx, cy)

    xsNew = 0*Xs0
    ysNew = 0*Ys0

    for ii in range(xdim):
        for jj in range(ydim):
            xyzVals = [Xs0[ii,jj], Ys0[ii,jj], 0, 1]
            xyzValsRot = np.dot(Mrot, xyzVals)
            xsNew[ii,jj] = xyzValsRot[0]
            ysNew[ii,jj] = xyzValsRot[1]

    # Calculate the interpolated DM grid (set extrapolated values to 0.0)
    dm.infMaster = scipy.interpolate.griddata((xsNew.reshape(-1),ysNew.reshape(-1)), inf0pad.reshape(-1), (Xs0.reshape(-1),Ys0.reshape(-1)), "cubic")
    dm.infMaster[np.isnan(dm.infMaster)] = 0
    dm.infMaster = dm.infMaster.reshape(inf0pad.shape)

    #Crop down the influence function until it has no zero padding left
    infSum = np.sum(dm.infMaster)
    infDiff = 0
    counter = 2
    while abs(infDiff) <= 1e-7:
        counter += 2
        Ninf0pad = len(dm.infMaster)-counter #Number of points across the rotated, cropped-down influence function at the original resolution
        infDiff = infSum - np.sum( dm.infMaster[1+counter//2:1-counter//2,1+counter//2:1-counter//2] ) #Subtract an extra 2 to negate the extra step that overshoots.

    counter -= 2
    Ninf0pad = len(dm.infMaster)-counter
    infMaster2 = dm.infMaster[1+counter//2:1-counter//2,1+counter//2:1-counter//2]

    dm.infMaster = infMaster2
    Npad = Ninf0pad

    x_inf0 = np.arange(-(Npad-1)//2, (Npad+1)//2)*dm.dx_inf0 # True for even- or odd-sized influence function maps as long as they are centered on the array.
    Xinf0,Yinf0 = np.meshgrid(x_inf0,x_inf0)

    #Compute the size of the postage stamps.
    dm.Nbox = falco.utils.ceil_even(Ninf0pad*dm.dx_inf0//dx_dm) # Number of points across the influence function in the pupil file's spacing. Want as even
    #Also compute their padded sizes for the angular spectrum (AS) propagation between P2 and DM1 or between DM1 and DM2
    Nmin = falco.utils.ceil_even( max(mp.sbp_center_vec)*max(np.abs([mp.d_P2_dm1, mp.d_dm1_dm2, mp.d_P2_dm1+mp.d_dm1_dm2]))//dx_dm**2 ) #Minimum number of points across for accurate angular spectrum propagation
    dm.NboxAS = max([dm.Nbox, Nmin]) #Uses a larger array if the max sampling criterion for angular spectrum propagation is violated

    ## Pad the pupil to at least the size of the DM(s) surface(s) to allow all actuators to be located outside the pupil. (Same for both DMs)

    #Find actuator farthest from center:
    dm.r_cent_act = np.sqrt(dm.xy_cent_act[0]**2 + dm.xy_cent_act[1]**2)
    dm.rmax = max(np.abs(dm.r_cent_act))
    NpixPerAct = dm.dm_spacing/dx_dm
    dm.NdmPad = falco.utils.ceil_even( ( dm.NboxAS + 2*(1+ (np.max(np.abs(dm.xy_cent_act))+0.5)*NpixPerAct)) ) # DM surface array padded by the width of the padded influence function to prevent indexing outside the array. The 1/2 term is because the farthest actuator center is still half an actuator away from the nominal array edge.

    #Compute coordinates (in meters) of the full DM array
    if dm.centering ==  "pixel":
        dm.x_pupPad = np.arange(-dm.NdmPad/2, dm.NdmPad/2)*dx_dm # meters, coords for the full DM arrays. Origin is centered on a pixel
    else:
        dm.x_pupPad = np.arange(-(dm.NdmPad-1)/2, (dm.NdmPad+1)/2)*dx_dm # meters, coords for the full DM arrays. Origin is centered between pixels for an even-sized array

    dm.y_pupPad = dm.x_pupPad

    #DM: (use NboxPad-sized postage stamps)

    if flagGenCube:
        #-Find the locations of the postage stamps arrays in the larger pupilPad array
        dm.xy_cent_act_inPix = dm.xy_cent_act*(dm.dm_spacing/dx_dm) # Convert units to pupil-file pixels
        dm.xy_cent_act_inPix = dm.xy_cent_act_inPix + 0.5 #For the half-pixel offset if pixel centered.

        dm.xy_cent_act_box = np.round(dm.xy_cent_act_inPix) # Center locations of the postage stamps (in between pixels), in actuator widths
        dm.xy_cent_act_box_inM = dm.xy_cent_act_box*dx_dm # now in meters
        dm.xy_box_lowerLeft = dm.xy_cent_act_box + (dm.NdmPad-dm.Nbox)/2 + 1 # indices of pixel in lower left of the postage stamp within the whole pupilPad array

        #Starting coordinates (in actuator widths) for updated influence function. This is interpixel centered, so do not translate!
        dm.x_box0 = np.arange(-(dm.Nbox-1)/2, (dm.Nbox+1)/2)*dx_dm
        dm.Xbox0,dm.Ybox0 = np.meshgrid(dm.x_box0, dm.x_box0) #meters, interpixel-centered coordinates for the master influence function

        #Limit the actuators used to those within 1 actuator width of the pupil
        r_cent_act_box_inM = np.sqrt(dm.xy_cent_act_box_inM[0]**2 + dm.xy_cent_act_box_inM[1]**2)
        #Compute and store all the influence functions:
        dm.inf_datacube = np.zeros((dm.Nbox,dm.Nbox,dm.NactTotal)) #initialize array of influence function "postage stamps"

        for iact in range(dm.NactTotal): #dm.Nact^2
            Xbox = dm.Xbox0 - (dm.xy_cent_act_inPix[0,iact]-dm.xy_cent_act_box[0,iact])*dx_dm
            Ybox = dm.Ybox0 - (dm.xy_cent_act_inPix[1,iact]-dm.xy_cent_act_box[1,iact])*dx_dm

            #For consistency, the MATLAB data was generated also with 'cubic' interpolation method in falco_gen_dm_poke_cube (dm.inf_datacube(:,:,iact) = interp2(Xinf0,Yinf0,dm.infMaster,dm.Xbox,dm.Ybox,'cubic',0);)
            #The result is still very slightly off because of non-identical interpolation algorithm
            dm.inf_datacube[:,:,iact] = scipy.interpolate.griddata((Xinf0.reshape(-1),Yinf0.reshape(-1)), dm.infMaster.reshape(-1), (Xbox.reshape(-1),Ybox.reshape(-1)), "cubic", fill_value=0).reshape(Xbox.shape)

    dm.act_ele = np.arange(dm.NactTotal)

class DeformableMirrorParameters:
    class _base_dm:
        def __init__(self, **kwargs):
            self.dm_spacing = _spec_arg("dm_spacing", kwargs, 0.001)
            self.dx_inf0 = _spec_arg("dx_inf0", kwargs, 0.0001)
            self.xc = _spec_arg("xc", kwargs, 23.5)
            self.dummy = _spec_arg("dummy", kwargs, 1)
            self.ytilt = _spec_arg("ytilt", kwargs, 0)
            self.xtilt = _spec_arg("xtilt", kwargs, 0)
            self.Nact = _spec_arg("Nact", kwargs, 48)
            self.edgeBuffer = _spec_arg("edgeBuffer", kwargs, 1)
            self.maxAbsV = _spec_arg("maxAbsV", kwargs, 125)
            self.yc = _spec_arg("yc", kwargs, 23.5)
            self.zrot = _spec_arg("zrot", kwargs, 0)
            self.VtoH = _spec_arg("VtoH", kwargs, 1e-9*np.ones((1, self.Nact)))
            self.inf0 = _spec_arg("inf0", kwargs, scipy.io.loadmat(
                _influence_function_file)["inf0"])

    def __init__(self, **kwargs):
        self.dm_weights = _spec_arg("dm_weights", kwargs, [1, 1, 1, 1, 1, 1, 1, 1, 1])
        self.dm_ind = _spec_arg("dm_ind", kwargs, [1, 2])
        self.maxAbsdV = _spec_arg("maxAbsdV", kwargs, 30)
        self.dm1 = _spec_arg("dm1", kwargs, self._base_dm())
        self.dm2 = _spec_arg("dm2", kwargs, self._base_dm())

    def init_ws(self, mp):
        self.flagZYX = False

        self.dm1.centering = mp.centering;
        self.dm2.centering = mp.centering;

        assert(self.dm_ind == [1,2])#Anything else is not implemented yet
        self.dm1.compact = self.dm1
        falco_gen_dm_poke_cube(self.dm1, mp, mp.P2.full.dx, False)
        falco_gen_dm_poke_cube(self.dm1.compact, mp, mp.P2.compact.dx);

        self.dm2.compact = self.dm2
        self.dm2.dx = mp.P2.full.dx
        self.dm2.compact.dx = mp.P2.compact.dx

        falco_gen_dm_poke_cube(self.dm2, mp, mp.P2.full.dx, False)
        falco_gen_dm_poke_cube(self.dm2.compact, mp, mp.P2.compact.dx)

        #Initial DM voltages
        self.dm1.V = np.zeros(self.dm1.Nact)
        self.dm2.V = np.zeros(self.dm2.Nact)

        #Peak-to-Valley DM voltages
        self.dm1.Vpv = np.zeros((mp.Nitr,1))
        self.dm2.Vpv = np.zeros((mp.Nitr,1))

        #First delta DM settings are zero (for covariance calculation in Kalman filters or robust controllers)
        self.dm1.dV = np.zeros((self.dm1.Nact,self.dm1.Nact)) # delta voltage on DM1
        self.dm2.dV = np.zeros((self.dm2.Nact,self.dm2.Nact)) # delta voltage on DM2

        #Store the DM commands
        self.dm1.Vall = np.zeros((self.dm1.Nact,self.dm1.Nact,mp.Nitr+1))
        self.dm2.Vall = np.zeros((self.dm2.Nact,self.dm2.Nact,mp.Nitr+1))

        #Initialize the number of actuators used.
        self.dm1.Nele = len(self.dm1.act_ele)
        self.dm2.Nele = len(self.dm2.act_ele)

        ## Array Sizes for Angular Spectrum Propagation with FFTs

        #Compact Model: Set nominal DM plane array sizes as a power of 2 for angular spectrum propagation with FFTs
        NdmPad = 2**np.ceil(1 + np.log2(max(self.dm1.compact.NdmPad, self.dm2.compact.NdmPad)))

        while (NdmPad < min(mp.sbp_center_vec)*abs(mp.d_dm1_dm2)/mp.P2.full.dx**2) or (NdmPad < min(mp.sbp_center_vec)*abs(mp.d_P2_dm1)/mp.P2.compact.dx**2): #Double the zero-padding until the angular spectrum sampling requirement is not violated
            NdmPad = 2*NdmPad

        self.compact = collections.namedtuple("_compact", "NdmPad")
        self.compact.NdmPad = NdmPad

        #Full Model: Set nominal DM plane array sizes as a power of 2 for angular spectrum propagation with FFTs
        NdmPad = 2**np.ceil(1 + np.log2(max(self.dm1.NdmPad, self.dm2.NdmPad)))
        while (NdmPad < min(mp.lam_array)*abs(mp.d_dm1_dm2)/mp.P2.full.dx**2) or (NdmPad < min(mp.lam_array)*abs(mp.d_P2_dm1)/mp.P2.full.dx**2): #Double the zero-padding until the angular spectrum sampling requirement is not violated
            NdmPad = 2*NdmPad

        self.full = collections.namedtuple("_full", "NdmPad")
        self.full.NdmPad = NdmPad

max_stars_repo_name: leonidprinceton/falco-python

max_issues_repo_name: leonidprinceton/falco-python

max_forks_repo_name: leonidprinceton/falco-python

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
#
# Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
# SPDX-License-Identifier: Apache-2.0

import numpy as np
import torch
import torch.nn as nn

from torch.utils.data import DataLoader

from dgllife.data import PubChemBioAssayAromaticity
from dgllife.model import load_pretrained
from dgllife.utils import EarlyStopping, Meter, RandomSplitter

from utils import set_random_seed, collate_molgraphs, load_model

def regress(args, model, bg):
    bg = bg.to(args['device'])
    atom_feats, bond_feats = bg.ndata.pop('hv'), bg.edata.pop('he')
    atom_feats, bond_feats = atom_feats.to(args['device']), bond_feats.to(args['device'])
    return model(bg, atom_feats, bond_feats)

def run_a_train_epoch(args, epoch, model, data_loader,
                      loss_criterion, optimizer):
    model.train()
    train_meter = Meter()
    for batch_id, batch_data in enumerate(data_loader):
        smiles, bg, labels, masks = batch_data
        labels, masks = labels.to(args['device']), masks.to(args['device'])
        prediction = regress(args, model, bg)
        loss = (loss_criterion(prediction, labels) * (masks != 0).float()).mean()
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        train_meter.update(prediction, labels, masks)
    total_score = np.mean(train_meter.compute_metric(args['metric_name']))
    print('epoch {:d}/{:d}, training {} {:.4f}'.format(
        epoch + 1, args['num_epochs'], args['metric_name'], total_score))

def run_an_eval_epoch(args, model, data_loader):
    model.eval()
    eval_meter = Meter()
    with torch.no_grad():
        for batch_id, batch_data in enumerate(data_loader):
            smiles, bg, labels, masks = batch_data
            labels = labels.to(args['device'])
            prediction = regress(args, model, bg)
            eval_meter.update(prediction, labels, masks)
        total_score = np.mean(eval_meter.compute_metric(args['metric_name']))
    return total_score

def main(args):
    args['device'] = torch.device("cuda: 0") if torch.cuda.is_available() else torch.device("cpu")
    set_random_seed(args['random_seed'])

    dataset = PubChemBioAssayAromaticity(smiles_to_graph=args['smiles_to_graph'],
                                         node_featurizer=args.get('node_featurizer', None),
                                         edge_featurizer=args.get('edge_featurizer', None))
    train_set, val_set, test_set = RandomSplitter.train_val_test_split(
            dataset, frac_train=args['frac_train'], frac_val=args['frac_val'],
            frac_test=args['frac_test'], random_state=args['random_seed'])
    train_loader = DataLoader(dataset=train_set,
                              batch_size=args['batch_size'],
                              shuffle=True,
                              collate_fn=collate_molgraphs)
    val_loader = DataLoader(dataset=val_set,
                            batch_size=args['batch_size'],
                            collate_fn=collate_molgraphs)
    test_loader = DataLoader(dataset=test_set,
                             batch_size=args['batch_size'],
                             collate_fn=collate_molgraphs)

    if args['pre_trained']:
        args['num_epochs'] = 0
        model = load_pretrained(args['exp'])
    else:
        model = load_model(args)
        loss_fn = nn.MSELoss(reduction='none')
        optimizer = torch.optim.Adam(model.parameters(), lr=args['lr'],
                                     weight_decay=args['weight_decay'])
        stopper = EarlyStopping(mode=args['mode'], patience=args['patience'])
    model.to(args['device'])

    for epoch in range(args['num_epochs']):
        # Train
        run_a_train_epoch(args, epoch, model, train_loader, loss_fn, optimizer)

        # Validation and early stop
        val_score = run_an_eval_epoch(args, model, val_loader)
        early_stop = stopper.step(val_score, model)
        print('epoch {:d}/{:d}, validation {} {:.4f}, best validation {} {:.4f}'.format(
            epoch + 1, args['num_epochs'], args['metric_name'], val_score,
            args['metric_name'], stopper.best_score))

        if early_stop:
            break

    if not args['pre_trained']:
        stopper.load_checkpoint(model)
    test_score = run_an_eval_epoch(args, model, test_loader)
    print('test {} {:.4f}'.format(args['metric_name'], test_score))

if __name__ == "__main__":
    import argparse

    from configure import get_exp_configure

    parser = argparse.ArgumentParser(description='Aromaticity Prediction')
    parser.add_argument('-m', '--model', type=str, choices=['AttentiveFP'],
                        default='AttentiveFP', help='Model to use')
    parser.add_argument('-p', '--pre-trained', action='store_true',
                        help='Whether to skip training and use a pre-trained model')
    args = parser.parse_args().__dict__
    args['dataset'] = 'Aromaticity'
    args['exp'] = '_'.join([args['model'], args['dataset']])
    args.update(get_exp_configure(args['exp']))

    main(args)

ext: py

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Tests for `et_ppmd` package."""

import pytest
import numpy as np
import et_ppmd.grid

def test_build():
    grid = et_ppmd.grid.Grid(cell_size=1.0, wx=2., )
    assert grid.cl.shape[0] == 2 # number of cells in x-direction
    assert grid.cl.shape[1] == 2 # number of cells in y-direction

    rx = np.array([0.5, 1.5, 0.5, 1.5,])
    ry = np.array([0.5, 0.5, 1.5, 1.5,])
    for build_method in [grid.build_simple, grid.build]:
        build_method(rx, ry, linearise=False)
        print(f'{build_method}\n{grid}')

        # all cells contain 1 atom
        assert np.all(grid.cl[:,:,0] == 1)
        # cell (0,0) contains atom 0
        assert np.all(grid.cl[0,0,1] == 0)
        # cell (1,0) contains atom 1
        assert np.all(grid.cl[1,0,1] == 1)
        # cell (0,1) contains atom 2
        assert np.all(grid.cl[0,1,1] == 2)
        # cell (1,1) contains atom 3
        assert np.all(grid.cl[1,1,1] == 3)

        grid.clear()

def test_build_random():
    max_atoms_per_cell = 40
    grid1 = et_ppmd.grid.Grid(cell_size=0.5, wx=1.0, max_atoms_per_cell=max_atoms_per_cell)
    grid2 = et_ppmd.grid.Grid(cell_size=0.5, wx=1.0, max_atoms_per_cell=max_atoms_per_cell)
    grid1.clear(full=True) # full clear, otherwise the test will fail because the entries in
    grid2.clear(full=True) # the cell lists that are not use.

    rx = np.random.random((100,))
    ry = np.random.random((100,))

    grid1.build_simple(rx, ry, linearise=False)
    print(grid1)

    grid2.build(rx, ry, linearise=False)
    print(grid2)

    # all entries should be equal since the order of atoms is the same.
    assert np.all(grid1.cl == grid2.cl)

def test_linearised():
    n_atoms = 100
    if n_atoms == 4: # a small test that is easy to debug
        rx = np.array([0.25, 0.75, 0.25, 0.75,])
        ry = np.array([0.25, 0.25, 0.75, 0.75,])
        max_atoms_per_cell = 2
    else:
        rx = np.random.random((n_atoms,))
        ry = np.random.random((n_atoms,))
        # estimate the maximum number of atoms per cell, this may accidentally be too small
        # because we are using random positions. The average number of atoms per cell is
        # the n_atoms/4 (there are 4 cells)
        max_atoms_per_cell = int((n_atoms/4)*1.5)

    grid1 = et_ppmd.grid.Grid(cell_size=0.5, wx=1.0, max_atoms_per_cell=max_atoms_per_cell)
    grid2 = et_ppmd.grid.Grid(cell_size=0.5, wx=1.0, max_atoms_per_cell=max_atoms_per_cell)

    grid1.build_simple(rx, ry, linearise=False)
    print(grid1)

    grid2.build(rx, ry, linearise=False)
    print(grid2)
    grid2.linearise()
    print(grid2)

    for l in range(grid1.n):
        for k in range(grid1.m):
            assert np.all(grid1.cell_list(k,l) == grid2.cell_list(k,l))

# ==============================================================================
# The code below is for debugging a particular test in eclipse/pydev.
# (normally all tests are run with pytest)
# ==============================================================================
if __name__ == "__main__":
    the_test_you_want_to_debug = test_linearised

    the_test_you_want_to_debug()
    print("-*# finished #*-")
# ==============================================================================


text:
# -*- coding: utf-8 -*-
from __future__ import print_function

import json
import sys
import numpy
import numpy as np
from time import time
from itertools import islice

def numpy_array_from_list_or_numpy_array(vectors):
    """
    Returns numpy array representation of argument.

    Argument maybe numpy array (input is returned)
    or a list of numpy vectors.
    """
    # If vectors is not a numpy matrix, create one
    if not isinstance(vectors, numpy.ndarray):
        V = numpy.zeros((vectors[0].shape[0], len(vectors)))
        for index in range(len(vectors)):
            vector = vectors[index]
            V[:, index] = vector
        return V

    return vectors

def perform_pca(A):
    """
    Computes eigenvalues and eigenvectors of covariance matrix of A.
    The rows of a correspond to observations, the columns to variables.
    """
    # First subtract the mean
    M = (A-numpy.mean(A, axis=0)).T
    # Get eigenvectors and values of covariance matrix
    return numpy.linalg.eig(numpy.cov(M))

def perform_online_pca(data, dimension=None):
    if dimension is None:
        dimension = data[0].shape

    total = 0
    mean = np.zeros(dimension, dtype=np.float64)
    comoment = np.zeros((dimension, dimension), dtype=np.float64)

    for element in data:
        total += len(element)
        last_mean = mean
        mean = mean + np.sum(element - mean, axis=0, dtype=np.float64)/total
        comoment += np.dot((element-mean).T, (element-last_mean))

    cov = comoment / (total-1)

    eigenvalues, eigenvectors = np.linalg.eig(cov)
    ordered_idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[ordered_idx]
    eigenvectors = eigenvectors[:, ordered_idx]

    return mean, (eigenvalues, eigenvectors)

PY2 = sys.version_info[0] == 2
if PY2:
    bytes_type = str
else:
    bytes_type = bytes

def want_string(arg, encoding='utf-8'):
    if isinstance(arg, bytes_type):
        rv = arg.decode(encoding)
    else:
        rv = arg
    return rv

def chunk(sequence, n):
    """ Yield successive n-sized chunks from sequence. """
    for i in xrange(0, len(sequence), n):
        yield sequence[i:i + n]

def ichunk(sequence, n):
    """ Yield successive n-sized chunks from sequence. """
    sequence = iter(sequence)
    chunk = list(islice(sequence, n))
    while len(chunk) > 0:
        yield chunk
        chunk = list(islice(sequence, n))

def load_dict_from_json(path):
    try:
        with open(path, "r") as json_file:
            return json.loads(json_file.read())
    except:
        print("infos.json is corrupted")

    return {}

def save_dict_to_json(path, dictionary):
    data = json.dumps(dictionary, indent=4, separators=(',', ': '))
    with open(path, "w") as json_file:
        json_file.write(data)

class Timer():
    def __init__(self, txt):
        self.txt = txt

    def __enter__(self):
        self.start = time()
        print(self.txt + "... ", end="")
        sys.stdout.flush()

    def __exit__(self, type, value, tb):
        print("{:.2f} sec.".format(time()-self.start))

max_stars_repo_name: MarcCote/NearPy

max_issues_repo_name: MarcCote/NearPy

max_forks_repo_name: MarcCote/NearPy


text:
# Copyright 2020 Huawei Technologies Co., Ltd
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import numpy as np

import mindspore as ms
import mindspore.nn as nn
from mindspore import Tensor
from mindspore import context
from mindspore.common.api import _executor
from mindspore.common.parameter import Parameter
from mindspore.ops import composite as C
from mindspore.ops import operations as P
from mindspore.common.initializer import initializer
from mindspore.nn import TrainOneStepCell, Momentum
from tests.ut.python.ops.test_math_ops import VirtualLoss

grad_all = C.GradOperation(get_all=True)

class NetWithLoss(nn.Cell):
    def __init__(self, network):
        super(NetWithLoss, self).__init__()
        self.loss = VirtualLoss()
        self.network = network

    def construct(self, x):
        predict = self.network(x)
        return self.loss(predict)

class GradWrap(nn.Cell):
    def __init__(self, network):
        super(GradWrap, self).__init__()
        self.network = network

    def construct(self, x):
        return grad_all(self.network)(x)

def test_unique_column_split():
    class Net(nn.Cell):
        def __init__(self):
            super().__init__()
            self.unique = P.Unique().shard(((1,),))
            self.relu = P.ReLU()
            self.mul = P.Mul()
            self.embedding_lookp = P.GatherV2().shard(((1, 8), (1,)))
            self.embedding_table = Parameter(initializer('normal', [2000, 128]),
                                             name='embedding_table')
            self.gatherv2 = P.GatherV2().shard(((1, 8), (1,)))
            self.reshape = P.Reshape()
            self.matmul = P.MatMul()
            self.mul_weight = Parameter(Tensor(np.full([32, 64, 1], 0.5, dtype=np.float32)), name="mul_weight")

        def construct(self, indices):
            indices_flatten = self.reshape(indices, (-1,))
            unique_id, unique_idx = self.unique(indices_flatten)
            unique_id_weight = self.embedding_lookp(self.embedding_table, unique_id, 0)
            weight_flatten = self.gatherv2(unique_id_weight, unique_idx, 0)
            weight = self.reshape(weight_flatten, (32, 64, 128))
            vx = self.mul(weight, self.mul_weight)
            return vx

    size = 8
    context.set_auto_parallel_context(device_num=size, global_rank=0, parallel_mode="auto_parallel")
    x = Tensor(np.ones([32, 64]), dtype=ms.int32)
    net = Net()
    optimizer = Momentum(net.trainable_params(), learning_rate=0.1, momentum=0.9)
    train_net = TrainOneStepCell(net, optimizer)
    train_net.set_auto_parallel()
    train_net.set_train()
    _executor.compile(train_net, x)

def test_unique_row_split():
    class Net(nn.Cell):
        def __init__(self):
            super().__init__()
            self.unique = P.Unique().shard(((1,),))
            self.relu = P.ReLU()
            self.mul = P.Mul()
            self.embedding_lookp = P.GatherV2().shard(((8, 1), (1,)))
            self.embedding_table = Parameter(initializer('normal', [2000, 128]),
                                             name='embedding_table')
            self.gatherv2 = P.GatherV2().shard(((1, 1), (8,)))
            self.reshape = P.Reshape()
            self.matmul = P.MatMul()
            self.mul_weight = Parameter(Tensor(np.full([32, 64, 1], 0.5, dtype=np.float32)), name="mul_weight")

        def construct(self, indices):
            indices_flatten = self.reshape(indices, (-1,))
            unique_id, unique_idx = self.unique(indices_flatten)
            unique_id_weight = self.embedding_lookp(self.embedding_table, unique_id, 0)
            weight_flatten = self.gatherv2(unique_id_weight, unique_idx, 0)
            weight = self.reshape(weight_flatten, (32, 64, 128))
            vx = self.mul(weight, self.mul_weight)
            return vx

    size = 8
    context.set_auto_parallel_context(device_num=size, global_rank=0, parallel_mode="stand_alone")
    x = Tensor(np.ones([32, 64]), dtype=ms.int32)
    net = Net()
    optimizer = Momentum(net.trainable_params(), learning_rate=0.1, momentum=0.9)
    train_net = TrainOneStepCell(net, optimizer)
    train_net.set_auto_parallel()
    train_net.set_train()
    _executor.compile(train_net, x)

max_stars_repo_name: dongkcs/mindspore

max_issues_repo_name: dongkcs/mindspore

max_forks_repo_name: dongkcs/mindspore


text:
from itertools import izip

from sklearn.base import RegressorMixin
from sklearn.cross_validation import train_test_split
from sklearn.datasets import make_regression
from sklearn.grid_search import GridSearchCV
from sklearn.linear_model import PassiveAggressiveRegressor
from sklearn.metrics import mean_absolute_error
from sklearn.metrics.pairwise import rbf_kernel, \
    polynomial_kernel
from sklearn.preprocessing import StandardScaler
from sklearn.utils import check_array
import numpy as np

from pa import PAEstimator, KernelPAEstimator, linear_kernel_local

__author__ = 'desouza'

class PARegressor(PAEstimator, RegressorMixin):
    def __init__(self, C=1, epsilon=0.01, loss="pai", n_iter=1,
                 fit_intercept=False):
        super(PARegressor, self).__init__(loss, C, n_iter, fit_intercept)

        self.epsilon = epsilon

    def partial_fit(self, X, y):
        # check_arrays(X, y)

        if self.loss == "pa":
            reg_func = self._pa
        elif self.loss == "pai":
            reg_func = self._pai
        elif self.loss == "paii":
            reg_func = self._paii

        # inits coefficients
        if self.coef_ is None:
            self.coef_ = np.zeros(X.shape[0])

        for _ in xrange(self.n_iter):
            y_hat_t = np.dot(self.coef_, X)
            discrepancy_t = np.abs(y - y_hat_t)
            loss_t = 0
            if discrepancy_t > self.epsilon:
                loss_t = discrepancy_t - self.epsilon

            tau_t = reg_func(loss_t, X)

            self.coef_ = self.coef_ + (np.sign(y - y_hat_t) * tau_t * X)

        return self

    def fit(self, X, y):
        check_array(X, y)

        for x_i, y_i in izip(X, y):
            self.partial_fit(x_i, y_i)

        return self

    def predict(self, X):
        X = check_array(X)

        # print self.coef_.shape
        # print X.shape
        y_preds = np.dot(self.coef_, X.T)
        return y_preds

class KernelPARegressor(KernelPAEstimator, RegressorMixin):
    def __init__(self, kernel="linear", gamma=0.01, C=1, epsilon=0.01,
                 loss="pa", n_iter=1,
                 fit_intercept=False):
        super(KernelPARegressor, self).__init__(kernel, gamma, loss, C, n_iter,
                                                fit_intercept)

        self.epsilon = epsilon

    def predict(self, X):
        X = check_array(X)

        kern_func = linear_kernel_local
        if self.kernel == "rbf":
            kern_func = rbf_kernel
        elif self.kernel == "poly":
            kern_func = polynomial_kernel

        support = np.array(self.support_)
        alphas = np.array(self.alphas_)

        phi = kern_func(support, X, self.gamma)
        margins = np.dot(phi.T, alphas)

        y_preds = margins

        return y_preds

    def fit(self, X, y):
        check_array(X, y)

        for x_i, y_i in izip(X, y):
            self.partial_fit(x_i, y_i)

        return self

    def partial_fit(self, X, y):

        reg_func = self._pai
        if self.loss == "pa":
            reg_func = self._pa
        elif self.loss == "paii":
            reg_func = self._paii

        kern_func = linear_kernel_local
        if self.kernel == "rbf":
            kern_func = rbf_kernel
        elif self.kernel == "poly":
            kern_func = polynomial_kernel

        # computes kernel and prediction
        if len(self.support_) > 0:
            # kernel function (phi) value
            # computes phi of support vectors against current example t
            support_t = np.array(self.support_)
            phi_t = kern_func(support_t, X, self.gamma)
            # takes the prediction
            alphas = np.array(self.alphas_)
            y_hat_t = np.dot(alphas.T, phi_t.ravel())

        else:  # support has 0 elements
            # then prediction is zero
            y_hat_t = 0

        # computes loss
        discrepancy_t = np.abs(y - y_hat_t)
        loss_t = 0
        if discrepancy_t > self.epsilon:
            loss_t = discrepancy_t - self.epsilon
        # loss_t = max(0, np.abs(y - y_hat_t) - self.epsilon)
        # loss_t = max(0, np.abs(y_hat_t - y) - self.epsilon)

            # computes weight for new support vector and updates the model
            kern_t = np.ravel(kern_func(X, X, self.gamma))
            # print kern_t

            tau_t = reg_func(loss_t, kern_t)
            alpha_t = np.sign(y - y_hat_t) * tau_t

            if alpha_t != 0:
                self.support_.append(X)
                self.alphas_.append(alpha_t)

        return self

def main_qeda_1000():
    # boston = load_boston()
    # X = boston.data
    # y = boston.target
    # X_train, X_test, y_train, y_test = train_test_split(X, y,
    # test_size=0.3, random_state=40)

    X_train = np.loadtxt(
        '/home/desouza/Projects/qe_da/domains_large_bb17/sample-nomarkup-1000'
        '/stl_40/features/1.0-it.src.csv',
        delimiter=',')
    y_train = np.loadtxt(
        '/home/desouza/Projects/qe_da/domains_large_bb17/sample-nomarkup-1000'
        '/stl_40/labels/1.0-it.src.hter')

    X_test = np.loadtxt(
        '/home/desouza/Projects/qe_da/domains_large_bb17/sample-nomarkup-1000'
        '/stl_40/features/it.tgt.csv',
        delimiter=',')
    y_test = np.loadtxt(
        '/home/desouza/Projects/qe_da/domains_large_bb17/sample-nomarkup-1000'
        '/stl_40/labels/it.tgt.hter')

    # X_train_bias = np.column_stack(tuple([X_train, np.ones(X_train.shape[
    # 0])]))
    # X_test_bias = np.column_stack(tuple([X_test, np.ones(X_test.shape[0])]))

    param_grid = {
        "C": [0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0],
        "epsilon": [0.0001, 0.001, 0.01, 0.1],
        "loss": ["pa", "pai", "paii"]}

    search = GridSearchCV(PARegressor(), param_grid,
        scoring='mean_absolute_error', n_jobs=8, iid=True,
        refit=True, cv=5, verbose=1)

    search.fit(X_train, y_train)
    # est = PARegressor(C=0.4, epsilon=0.01, loss="paii")

    # for i in range(10):
    # for x_i, y_i in izip(X_train_bias, y_train):
    # for x_i, y_i in izip(X_train, y_train):
    # est.partial_fit(x_i, y_i)
    # est.fit(X_train, y_train)

    # y_preds = est.predict(X_test_bias)
    # y_preds = est.predict(X_test)
    y_preds = search.predict(X_test)

    # print est.coef_
    print search.best_estimator_.coef_
    print search.best_params_

    print y_preds.shape

    mae_my_pa = mean_absolute_error(y_preds, y_test)
    print "My PA MAE = %2.4f" % mae_my_pa

    param_grid = {
        "C": [0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0],
        "epsilon": [0.0001, 0.001, 0.01, 0.1],
        "loss": ["epsilon_insensitive", "squared_epsilon_insensitive"]}

    search = GridSearchCV(PassiveAggressiveRegressor(), param_grid,
        scoring='mean_absolute_error', n_jobs=8, iid=True,
        refit=True, cv=5, verbose=1)
    # est_pa = PassiveAggressiveRegressor()
    # est_pa.partial_fit(X_train, y_train)
    search.fit(X_train, y_train)
    # y_preds_pa = est_pa.predict(X_test)
    y_preds_pa = search.predict(X_test)
    mae_sk_pa = mean_absolute_error(y_preds_pa, y_test)
    print "Sklearn PA MAE = %2.4f" % mae_sk_pa

    # print est_pa.coef_
    print search.best_estimator_.coef_
    print search.best_params_

def main():
    X, y, coef = make_regression(1000, 200, 10, 1, noise=0.05, coef=True,
                                 random_state=42)

    # X = np.column_stack((X, np.ones(X.shape[0])))

    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3,
                                                        random_state=42)

    # sca = StandardScaler()
    # sca.fit(X_train)
    # X_train = sca.transform(X_train)
    # X_test = sca.transform(X_test)

    # print X.shape
    # print y.shape
    # print coef.shape

    param_grid = {
        "C": [0.0000001, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 10,
              100, 1000],
        "epsilon": [0.0001, 0.001, 0.01, 0.1]}

    param_grid_kern = {
        "C": [0.0000001, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 10,
              100, 1000],
        "epsilon": [0.0001, 0.001, 0.01, 0.1],
        "gamma": [0.00001, 0.0001, 0.001, 0.01, 0.1, 1, 10, 100]}
    # "loss": ["pa", "pai", "paii"]}}

    my_pa = PARegressor(loss="paii", C=1, epsilon=0.001, n_iter=1,
                        fit_intercept=False)
    #
    # search = GridSearchCV(my_pa, param_grid,
    #                       scoring='mean_absolute_error', n_jobs=8, iid=True, refit=True, cv=5,
    #                       verbose=1)
    # search.fit(X_train, y_train)
    # print search.best_params_

    my_pa.fit(X_train, y_train)
    print my_pa.coef_

    # y_preds = search.predict(X_test)
    y_preds = my_pa.predict(X_test)

    mae_my_pa = mean_absolute_error(y_test, y_preds)
    print "My PA MAE = %2.4f" % mae_my_pa

    my_kpa_linear = KernelPARegressor(kernel="linear", loss="paii", C=1, epsilon=0.001, n_iter=1, fit_intercept=False)
    my_kpa_linear.fit(X_train, y_train)
    print "alphas", len(my_kpa_linear.alphas_), my_kpa_linear.alphas_
    y_preds = my_kpa_linear.predict(X_test)
    mae_kpa_linear = mean_absolute_error(y_test, y_preds)
    print "My KPA linear MAE = %2.4f" % mae_kpa_linear

    my_kpa_rbf = KernelPARegressor(kernel="rbf", loss="paii", gamma=0.001, C=1, epsilon=0.001, n_iter=1, fit_intercept=False)
    # search = GridSearchCV(my_kpa_rbf, param_grid_kern,
    #                       scoring='mean_absolute_error', n_jobs=8, iid=True, refit=True, cv=5,
    #                       verbose=1)
    # search.fit(X_train, y_train)

    my_kpa_rbf.fit(X_train, y_train)
    print "alphas", len(my_kpa_rbf.alphas_), my_kpa_rbf.alphas_
    print "support", len(my_kpa_rbf.support_)
    # print "alphas", len(search.best_estimator_.alphas_)  # , my_kpa_rbf.alphas_
    # print "support", len(search.best_estimator_.support_)
    # print search.best_params_
    y_preds = my_kpa_rbf.predict(X_test)
    # y_preds = search.predict(X_test)
    mae_my_kpa = mean_absolute_error(y_test, y_preds)
    print "My Kernel PA MAE = %2.4f" % mae_my_kpa

    # print search.best_estimator_
    # print np.corrcoef(search.best_estimator_.coef_, coef)

    # param_grid = {
    # "C": [0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 10,
    #           100, 1000, 10000],
    #     "epsilon": [0.0001, 0.001, 0.01, 0.1],
    #     # "loss": ["epsilon_insensitive", "squared_epsilon_insensitive"]}
    #     "loss": ["squared_epsilon_insensitive"]}

    # search = GridSearchCV(PassiveAggressiveRegressor(fit_intercept=True),
    # param_grid, scoring='mean_absolute_error', n_jobs=8, iid=True,
    # refit=True, cv=5, verbose=1)
    # search.fit(X_train, y_train)

    sk_pa = PassiveAggressiveRegressor(loss="squared_epsilon_insensitive", C=1,
                                       epsilon=0.001, n_iter=1,
                                       fit_intercept=False,
                                       warm_start=True)
    for i in xrange(X_train.shape[0]):
        # for x_i, y_i in zip(X_train, y_train):
        x = np.array(X_train[i], ndmin=2)
        y = np.array(y_train[i], ndmin=1)
        # print x.shape
        # print y
        sk_pa.partial_fit(x, y)

    # sk_pa.fit(X_train, y_train)

    # y_preds = search.predict(X_test)
    y_preds = sk_pa.predict(X_test)
    mae_sk_pa = mean_absolute_error(y_preds, y_test)
    print "Sklearn PA MAE = %2.4f" % mae_sk_pa
    # print search.best_estimator_
    # print np.corrcoef(search.best_estimator_.coef_, coef)

    # en = ElasticNet(fit_intercept=True)
    # en.fit(X_train, y_train)
    # y_preds = en.predict(X_test)
    # mae_en = mean_absolute_error(y_test, y_preds)
    # print "Sklearn EN MAE = %2.4f" % mae_en

    # svr = SVR(kernel="linear")
    # svr.fit(X_train, y_train)
    # y_preds = svr.predict(X_test)
    # mae_svr = mean_absolute_error(y_test, y_preds)
    # print "Sklearn SVR MAE = %2.4f" % mae_svr

    # sgd = SGDRegressor(loss="epsilon_insensitive", penalty="l2",
    # epsilon=0.001)
    # # sgd.fit(X_train, y_train)
    # for i in xrange(X_train.shape[0]):
    # # for x_i, y_i in zip(X_train, y_train):
    # x = np.array(X_train[i], ndmin=2)
    #     y = np.array(y_train[i], ndmin=1)
    #     # print x.shape
    #     # print y
    #     sgd.partial_fit(x, y)
    #
    # y_preds = sgd.predict(X_test)
    # mae_sgd = mean_absolute_error(y_test, y_preds)
    # print "Sklearn SGD MAE = %2.4f" % mae_sgd

if __name__ == "__main__":
    main()
    # main_qeda_1000()

max_stars_repo_name: jsouza/pamtl

max_issues_repo_name: jsouza/pamtl

max_forks_repo_name: jsouza/pamtl


text:
import sys
import numpy
import gzip
import math
import itertools
import time
import codecs
import argparse
import pickle
from index_data import DataProcessor
from onto_attention import OntoAttentionLSTM
from keras.models import Model, model_from_yaml
from keras.layers import Input, Dropout, LSTM, Dense
from keras.layers.wrappers import TimeDistributed
from keras.callbacks import EarlyStopping
from keras_extensions import HigherOrderEmbedding

class SentenceModel(object):
  def __init__(self, word_dim=50, num_senses=2, num_hyps=5):
    self.dp = DataProcessor(word_syn_cutoff=num_senses, syn_path_cutoff=num_hyps)
    self.num_hyps = num_hyps
    self.num_senses = num_senses
    self.numpy_rng = numpy.random.RandomState(12345)
    self.word_dim = word_dim
    self.model = None

  def read_sentences(self, tagged_sentences, sentlenlimit=None, test=False, remove_singletons=False):
    num_sentences = len(tagged_sentences)
    all_words = []
    all_pos_tags = []
    maxsentlen = 0
    for tagged_sentence in tagged_sentences:
      words = []
      pos_tags = []
      # Expects each token to be a "_" separated combination of word and POS tag.
      tagged_words = tagged_sentence.split(" ")
      if sentlenlimit is not None:
        tagged_words = tagged_words[:sentlenlimit]
      for word_tag in tagged_words:
        parts = word_tag.split("_")
        tag = parts[-1]
        word = "_".join(parts[:-1]).lower()
        words.append(word)
        pos_tags.append(tag)
      if len(words) > maxsentlen:
        maxsentlen = len(words)
      all_words.append(words)
      all_pos_tags.append(pos_tags)
    if not sentlenlimit:
      sentlenlimit = maxsentlen
    C_ind = numpy.zeros((num_sentences, sentlenlimit, self.num_senses, self.num_hyps), dtype='int32')
    S_ind = numpy.zeros((num_sentences, sentlenlimit), dtype='int32')
    for i, (words, pos_tags) in enumerate(zip(all_words, all_pos_tags)):
      sentlen = len(words)
      # test=True locks the word and syn index dicts. No new keys will be added
      word_inds, syn_inds = self.dp.index_sentence(words, pos_tags, test=test, remove_singletons=remove_singletons)
      S_ind[i][-sentlen:] = word_inds
      for j in range(sentlen):
        sense_syn_ind = syn_inds[j]
        sense_syn_ind_len = len(sense_syn_ind)
        for k, syn_ind in enumerate(sense_syn_ind):
          C_ind[i][-sentlen+j][-sense_syn_ind_len+k][-len(syn_ind):] = syn_ind
    return S_ind, C_ind

  def _make_one_hot(self, word_inds, vec_size):
    onehot = numpy.zeros((word_inds.shape + (vec_size,)))
    for inds in itertools.product(*[numpy.arange(s) for s in word_inds.shape]):
      onehot[inds+(word_inds[inds],)] = 1
    return onehot

  def _factor_target_indices(self, Y_inds, vocab_size=None, base=2):
    if vocab_size is None:
      vocab_size = len(self.dp.word_index)
    print >>sys.stderr, "Factoring targets of vocabulary size: %d"%(vocab_size)
    num_vecs = int(math.ceil(math.log(vocab_size)/math.log(base))) + 1
    base_inds = []
    div_Y_inds = Y_inds
    print >>sys.stderr, "Number of factors: %d"%num_vecs
    for i in range(num_vecs):
      new_inds = div_Y_inds % base
      if i == num_vecs - 1:
        if new_inds.sum() == 0:
          # Most significant "digit" is a zero. Omit it.
          break
      base_inds.append(new_inds)
      div_Y_inds = numpy.copy(div_Y_inds/base)
    base_vecs = [self._make_one_hot(base_inds_i, base) for base_inds_i in base_inds]
    return base_vecs

  def train(self, S_ind, C_ind, use_onto_lstm=True, use_attention=True, num_epochs=20,  hierarchical=False, base=2):
    # Predict next word from current synsets
    X = C_ind[:,:-1] if use_onto_lstm else S_ind[:,:-1] # remove the last words' hyps in all sentences
    Y_inds = S_ind[:,1:] # remove the first words in all sentences
    if hierarchical:
      train_targets = self._factor_target_indices(Y_inds, base=base)
    else:
      train_targets = [self._make_one_hot(Y_inds, Y_inds.max() + 1)]
    length = Y_inds.shape[1]
    lstm_outdim = self.word_dim

    num_words = len(self.dp.word_index)
    num_syns = len(self.dp.synset_index)
    input = Input(shape=X.shape[1:], dtype='int32')
    embed_input_dim = num_syns if use_onto_lstm else num_words
    embed_layer = HigherOrderEmbedding(name='embedding', input_dim=embed_input_dim, output_dim=self.word_dim, input_shape=X.shape[1:], mask_zero=True)
    sent_rep = embed_layer(input)
    reg_sent_rep = Dropout(0.5)(sent_rep)
    if use_onto_lstm:
      lstm_out = OntoAttentionLSTM(name='sent_lstm', input_dim=self.word_dim, output_dim=lstm_outdim, input_length=length, num_senses=self.num_senses, num_hyps=self.num_hyps, return_sequences=True, use_attention=use_attention)(reg_sent_rep)
    else:
      lstm_out = LSTM(name='sent_lstm', input_dim=self.word_dim, output_dim=lstm_outdim, input_length=length, return_sequences=True)(reg_sent_rep)
    output_nodes = []
    # Make one node for each factored target
    for target in train_targets:
      node = TimeDistributed(Dense(input_dim=lstm_outdim, output_dim=target.shape[-1], activation='softmax'))(lstm_out)
      output_nodes.append(node)

    model = Model(input=input, output=output_nodes)
    print >>sys.stderr, model.summary()
    early_stopping = EarlyStopping()
    precompile_time = time.time()
    model.compile(loss='categorical_crossentropy', optimizer='adam')
    postcompile_time = time.time()
    print >>sys.stderr, "Model compilation took %d s"%(postcompile_time - precompile_time)
    model.fit(X, train_targets, nb_epoch=num_epochs, validation_split=0.1, callbacks=[early_stopping])
    posttrain_time = time.time()
    print >>sys.stderr, "Training took %d s"%(posttrain_time - postcompile_time)
    concept_reps = model.layers[1].get_weights()
    self.model = model
    return concept_reps

  def test(self, vocab_size, use_onto_lstm, S_ind_test=None, C_ind_test=None, hierarchical=False, base=2, oov_list=None):
    X_test = C_ind_test[:,:-1] if use_onto_lstm else S_ind_test[:,:-1] # remove the last words' hyps in all sentences
    Y_inds_test = S_ind_test[:,1:]
    if hierarchical:
      test_targets = self._factor_target_indices(Y_inds_test, vocab_size, base=base)
    else:
      test_targets = [self._make_one_hot(Y_inds_test, vocab_size)]
    print >>sys.stderr, "Evaluating model on test data"
    test_loss = self.model.evaluate(X_test, test_targets)
    print >>sys.stderr, "Test loss: %.4f"%test_loss
    if oov_list is not None:
      oov_inds = [self.dp.word_index[w] for w in oov_list]
      non_oov_Y_inds = numpy.copy(Y_inds_test)
      for ind in oov_inds:
	non_oov_Y_inds[non_oov_Y_inds == ind] = 0
      non_oov_test_targets = self._factor_target_indices(non_oov_Y_inds, vocab_size, base=base)
      non_oov_test_loss = self.model.evaluate(X_test, non_oov_test_targets)
      print >>sys.stderr, "Non-oov test loss: %.4f"%non_oov_test_loss
    factored_test_preds = [-((numpy.log(pred) * target).sum(axis=-1)) for pred, target in zip(self.model.predict(X_test), test_targets)]
    test_preds = sum(factored_test_preds)
    #non_null_probs = []
    #for test_pred, inds in zip(test_preds, Y_inds_test):
    #  wanted_probs = []
    #  for tp, ind in zip(test_pred, inds):
    #    if ind != 0:
    #      wanted_probs.append(tp)
    #  non_null_probs.append(wanted_probs)
    #return non_null_probs
    return test_preds

  def get_attention(self, C_ind):
    if not self.model:
      raise RuntimeError, "Model not trained!"
    model_embedding = None
    model_weights = None
    for layer in self.model.layers:
      if layer.name.lower() == "embedding":
        model_embedding = layer
      if layer.name.lower() == "sent_lstm":
        model_lstm = layer
    if model_embedding is None or model_lstm is None:
      raise RuntimeError, "Did not find expected layers"
    lstm_weights = model_lstm.get_weights()
    embedding_weights = model_embedding.get_weights()
    embed_in_dim, embed_out_dim = embedding_weights[0].shape
    att_embedding = HigherOrderEmbedding(input_dim=embed_in_dim, output_dim=embed_out_dim, weights=embedding_weights)
    onto_lstm = OntoAttentionLSTM(input_dim=embed_out_dim, output_dim=embed_out_dim, input_length=model_lstm.input_length, num_senses=self.num_senses, num_hyps=self.num_hyps, use_attention=True, return_attention=True, weights=lstm_weights)
    att_input = Input(shape=C_ind.shape[1:], dtype='int32')
    att_sent_rep = att_embedding(att_input)
    att_output = onto_lstm(att_sent_rep)
    att_model = Model(input=att_input, output=att_output)
    att_model.compile(optimizer='adam', loss='mse') # optimizer and loss are not needed since we are not going to train this model.
    C_att = att_model.predict(C_ind)
    print >>sys.stderr, "Got attention values. Input, output shapes:", C_ind.shape, C_att.shape
    return C_att

if __name__ == '__main__':
  argparser = argparse.ArgumentParser(description="Model sentences using ontoLSTM")
  argparser.add_argument('--train_file', type=str, help="One sentence per line, POS tagged")
  argparser.add_argument('--test_file', type=str, help="Test file for which attention values will be printed. One sentence per line, POS tagged")
  argparser.add_argument('--dim', type=int, help="Word/synset dimensionality", default=50)
  argparser.add_argument('--num_senses', type=int, help="Number of senses per word if using OntoLSTM (default 2)", default=2)
  argparser.add_argument('--num_hyps', type=int, help="Number of hypernyms per sense if using OntoLSTM (default 5)", default=5)
  argparser.add_argument('--use_onto_lstm', help="If this flag is not set, will use traditional LSTM", action='store_true')
  argparser.add_argument('--use_attention', help="Use attention in ontoLSTM. If this flag is not set, will use average concept representations", action='store_true')
  argparser.add_argument('--show_attention', help="Print attention values for sentences in the test file (printed in <testfilename>.att_out)", action='store_true')
  argparser.add_argument('--hierarchical', type=int, help="Use hierarchical softmax with the provided number of classes per factor")
  argparser.add_argument('--synset_embedding_output', type=str, help="Print learned synset representations in the given file")
  argparser.add_argument('--num_epochs', type=int, help="Number of epochs (default 20)", default=20)
  args = argparser.parse_args()
  if args.show_attention and (not args.use_attention or not args.use_onto_lstm):
    raise RuntimeError, "Use OntoLSTM with attention to print attention values of the test file"
  sm = SentenceModel(word_dim=args.dim, num_senses=args.num_senses, num_hyps=args.num_hyps)
  hierarchical = False
  base = 2
  if args.hierarchical is not None:
    hierarchical = True
    base = args.hierarchical
  model_name_prefix = "sent_model_ontolstm=%s_att=%s_senses=%d_hyps=%d"%(str(args.use_onto_lstm), str(args.use_attention), args.num_senses, args.num_hyps)
  do_test = False
  do_train = False

  if args.train_file is not None:
    ts = [x.strip() for x in codecs.open(args.train_file, "r", "utf-8").readlines()]
    print >>sys.stderr, "Reading training data"
    S_ind, C_ind = sm.read_sentences(ts, remove_singletons=True)
    print >>sys.stderr, "Target vocab size after reading training data: %d"%(len(sm.dp.word_index))
    singletons = []
    for w, n in sm.dp.word_index.items():
      if n == 1:
        singletons.append(w)
    print >>sys.stderr, "Number of singletons: %d. Here are some:"%(len(singletons)), singletons[:5]
    _, train_sent_len, _, _ = C_ind.shape
    do_train = True
  else:
    print >>sys.stderr, "Loading stored model"
    sm.model = model_from_yaml(open("%s.yaml"%model_name_prefix).read(), custom_objects={"HigherOrderEmbedding": HigherOrderEmbedding, "OntoAttentionLSTM": OntoAttentionLSTM})
    print >>sys.stderr, sm.model.summary()
    sm.model.load_weights("%s.h5"%model_name_prefix)
    sm.model.compile(optimizer='adam', loss='categorical_crossentropy')
    dataproc_pkl_file = open("%s_dataproc.pkl"%model_name_prefix)
    sm.dp = pickle.load(dataproc_pkl_file)
    print >>sys.stderr, "Target vocab size of loaded index: %d"%(len(sm.dp.word_index))
    print sm.model.get_input_shape_at(0)
    train_sent_len = sm.model.get_input_shape_at(0)[1] + 1 #because the input shape will be one less than the sentence length

  if args.test_file is not None:
    print >>sys.stderr, "Reading test data"
    ts_test = [x.strip() for x in codecs.open(args.test_file, "r", "utf-8").readlines()]
    #lock_index = False if do_train else True
    lock_index = True
    S_ind_test, C_ind_test = sm.read_sentences(ts_test, sentlenlimit=train_sent_len, test=lock_index)
    print >>sys.stderr, "Target vocab size after reading test data: %d"%(len(sm.dp.word_index))
    do_test = True

  if do_train:
    concept_reps = sm.train(S_ind, C_ind, use_onto_lstm=args.use_onto_lstm, use_attention=args.use_attention, num_epochs=args.num_epochs, hierarchical=hierarchical, base=base)
    model_yaml_string = sm.model.to_yaml()
    open("%s.yaml"%model_name_prefix, "w").write(model_yaml_string)
    sm.model.save_weights("%s.h5"%model_name_prefix, overwrite=True)
    dataproc_pkl_file = open("%s_dataproc.pkl"%model_name_prefix, "w")
    pickle.dump(sm.dp, dataproc_pkl_file)

  if do_test:
    vocab_size = len(sm.dp.word_index)
    test_probs = sm.test(vocab_size, args.use_onto_lstm, S_ind_test, C_ind_test, hierarchical, base)
    outfile = open("%s.out"%model_name_prefix, "w")
    for nnp in test_probs:
      print >>outfile, " ".join(["%.4f"%p for p in nnp])
    outfile.close()
  if args.synset_embedding_output is not None:
    concrepfile = open(args.synset_embedding_output, "w")
    for syn in sm.dp.synset_index:
      print >>concrepfile, syn, " ".join(["%s"%x for x in concept_reps[0][sm.dp.synset_index[syn]]])
  rev_synset_ind = {ind: syn for (syn, ind) in sm.dp.synset_index.items()}
  if args.show_attention:
    if args.test_file is None:
      raise RuntimeError, "Test file not provided!"
    outfile_name = args.test_file.split("/")[-1] + ".att_out"
    outfile = open(outfile_name, "w")
    C_att = sm.get_attention(C_ind_test)
    for i, (sent, sent_c_inds, sent_c_atts) in enumerate(zip(ts_test, C_ind_test, C_att)):
      print >>outfile, "SENT %d: %s"%(i, sent)
      words = sent.split()[:train_sent_len]
      word_id = 0
      for word_c_inds, word_c_atts in zip(sent_c_inds, sent_c_atts):
        if word_c_inds.sum() == 0:
          continue
        sense_id = 0
        #print >>outfile, "Attention for %s"%(words[word_id])
        best_sense = ""
        max_sense_prob = 0.0
        for s_h_ind, s_h_att in zip(word_c_inds, word_c_atts):
          if sum(s_h_ind) == 0:
            continue
          #print >>outfile, "\nSense %d"%(sense_id)
          sense_id += 1
          if s_h_att[-1] > max_sense_prob:
            max_sense_prob = s_h_att[-1]
            best_sense = rev_synset_ind[s_h_ind[-1]]
          #for h_ind, h_att in zip(s_h_ind, s_h_att):
          #  print >>outfile, rev_synset_ind[h_ind], h_att
        print >>outfile, words[word_id], best_sense
        word_id += 1
        #print >>outfile
      print >>outfile

max_stars_repo_name: pdasigi/onto-lstm

max_issues_repo_name: StevenLOL/onto-lstm

max_forks_repo_name: StevenLOL/onto-lstm


text:
import argparse
import sys
import torch
import os
import numpy as np
import time

from torch.utils import data

from datasets.diagnostic_dataset import DiagnosticDataset
from enums.dataset_type import DatasetType
from helpers.metadata_helper import get_metadata_properties
from helpers.train_helper import TrainHelper
from models.diagnostic_ensemble import DiagnosticEnsemble
from metrics.average_ensemble_meter import AverageEnsembleMeter

import matplotlib.pyplot as plt

header = '  Time Epoch Iteration    Progress (%Epoch) | Loss-Avg  Acc-Avg | Loss-Color Loss-Shape Loss-Size Loss-PosH Loss-PosW | Acc-Color Acc-Shape Acc-Size Acc-PosH Acc-PosW |    Dataset | Best'
log_template = ' '.join(
    '{:>6.0f},{:>5.0f},{:>9.0f},{:>5.0f}/{:<5.0f} {:>7.0f}%,| {:>8.6f} {:>7.6f} | {:>10.6f} {:>10.6f} {:>9.6f} {:>9.6f} {:>9.6f} | {:>9.6f} {:>9.6f} {:>8.6f} {:>8.6f} {:>8.6f} | {:>10s} | {:>4s}'.split(','))

start_time = time.time()

def parse_arguments(args):
    # Training settings
    parser = argparse.ArgumentParser(
        description="Training Sender/Receiver Agent on a task"
    )
    parser.add_argument(
        "--max-epochs",
        type=int,
        default=1000,
        metavar="N",
        help="number of batch epochs to train (default: 1k)",
    )
    parser.add_argument(
        "--messages-seed",
        type=int,
        required=True,
        metavar="S",
        help="The seed which was used for sampling messages"
    )
    parser.add_argument(
        "--training-seed",
        type=int,
        default=13,
        help="The seed that will be used to initialize the current model (default: 13)"
    )
    parser.add_argument(
        "--embedding-size",
        type=int,
        default=32,
        metavar="N",
        help="embedding size for embedding layer (default: 64)",
    )
    parser.add_argument(
        "--hidden-size",
        type=int,
        default=64,
        metavar="N",
        help="hidden size for hidden layer (default: 64)",
    )
    parser.add_argument(
        "--batch-size",
        type=int,
        default=1024,
        metavar="N",
        help="input batch size for training (default: 1024)",
    )
    parser.add_argument(
        "--max-length",
        type=int,
        default=10,
        metavar="N",
        help="max sentence length allowed for communication (default: 10)",
    )
    parser.add_argument(
        "--vocab-size",
        type=int,
        default=25,
        metavar="N",
        help="Size of vocabulary (default: 25)",
    )
    parser.add_argument(
        "--lr",
        type=float,
        default=0.01,
        metavar="N",
        help="Adam learning rate (default: 0.01)",
    )
    parser.add_argument(
        "--device",
        type=str,
        help="Device to be used. Pick from none/cpu/cuda. If default none is used automatic check will be done")

    parser.add_argument(
        "--training-sample-count",
        type=int,
        default=None,
        help="How many models should be sampled while training. Default is None, meaning that all models will be used")

    args = parser.parse_args(args)

    return args

def print_results(accuracies_meter: AverageEnsembleMeter, losses_meter: AverageEnsembleMeter, epoch, max_epochs, dataset: str, new_best:bool):
    print(log_template.format(
        time.time()-start_time,
        epoch,
        epoch,
        1 + epoch,
        max_epochs,
        100. * (1+epoch) / max_epochs,
        losses_meter.avg,
        accuracies_meter.avg,
        losses_meter.averages[0],
        losses_meter.averages[1],
        losses_meter.averages[2],
        losses_meter.averages[3],
        losses_meter.averages[4],
        accuracies_meter.averages[0],
        accuracies_meter.averages[1],
        accuracies_meter.averages[2],
        accuracies_meter.averages[3],
        accuracies_meter.averages[4],
        dataset,
        "BEST" if new_best else ""
    ))

def perform_iteration(model: DiagnosticEnsemble, dataloader, batch_size: int, device, sample_count=5):
    accuracies_meter = AverageEnsembleMeter(number_of_values=5)
    losses_meter = AverageEnsembleMeter(number_of_values=5)

    for messages, properties in dataloader:
        messages = messages.long().to(device)
        properties = properties.long().to(device)

        current_accuracies, current_losses = model.forward(messages, properties, sample_count)

        accuracies_meter.update(current_accuracies)
        losses_meter.update(current_losses)

    return accuracies_meter, losses_meter

def generate_unique_name(length, vocabulary_size, seed):
    result = f'max_len_{length}_vocab_{vocabulary_size}_seed_{seed}'
    return result

def generate_model_name(length, vocabulary_size, messages_seed, training_seed):
    result = f'max_len_{length}_vocab_{vocabulary_size}_msg_seed_{messages_seed}_train_seed_{training_seed}'
    result += '.p'
    return result

def save_model(path, model):
    torch_state = {'model_state_dict': model.state_dict(),
                   'optimizer0_state_dict': model.optimizers[0].state_dict(),
                   'optimizer1_state_dict': model.optimizers[1].state_dict(),
                   'optimizer2_state_dict': model.optimizers[2].state_dict(),
                   'optimizer3_state_dict': model.optimizers[3].state_dict(),
                   'optimizer4_state_dict': model.optimizers[4].state_dict()}

    torch.save(torch_state, path)

def load_model(path, model):
    checkpoint = torch.load(path)
    model.load_state_dict(checkpoint['model_state_dict'])

    for i, optimizer in enumerate(model.optimizers):
        if checkpoint[f'optimizer{i}_state_dict']:
            optimizer.load_state_dict(checkpoint[f'optimizer{i}_state_dict'])

def initialize_model(args, device, checkpoint_path: str = None):
    diagnostic_model = DiagnosticEnsemble(
        num_classes_by_model=[3, 3, 2, 3, 3],
        batch_size=args.batch_size,
        device=device,
        embedding_size=args.embedding_size,
        num_hidden=args.hidden_size,
        learning_rate=args.lr)

    if checkpoint_path:
        load_model(checkpoint_path, diagnostic_model)

    return diagnostic_model

def baseline(args):
    args = parse_arguments(args)

    if args.device:
        device = torch.device(args.device)
    else:
        device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

    inference_path = os.path.join('data', 'inference')
    if not os.path.exists(inference_path):
        os.mkdir(inference_path)

    unique_name = generate_unique_name(
        length=args.max_length,
        vocabulary_size=args.vocab_size,
        seed=args.messages_seed)

    model_name = generate_model_name(
        length=args.max_length,
        vocabulary_size=args.vocab_size,
        messages_seed=args.messages_seed,
        training_seed=args.training_seed)

    model_path = os.path.join(inference_path, model_name)

    train_helper = TrainHelper(device)
    train_helper.seed_torch(args.training_seed)

    train_dataset = DiagnosticDataset(unique_name, DatasetType.Train)
    train_dataloader = data.DataLoader(train_dataset, num_workers=1, pin_memory=True, shuffle=True, batch_size=args.batch_size)

    validation_dataset = DiagnosticDataset(unique_name, DatasetType.Valid)
    validation_dataloader = data.DataLoader(validation_dataset, num_workers=1, pin_memory=True, shuffle=False, batch_size=args.batch_size)

    test_dataset = DiagnosticDataset(unique_name, DatasetType.Test)
    test_dataloader = data.DataLoader(test_dataset, num_workers=1, pin_memory=True, shuffle=False, batch_size=args.batch_size)

    diagnostic_model = initialize_model(args, device)

    diagnostic_model.to(device)

    # Setup the loss and optimizer

    print(header)

    best_accuracy = -1.

    for epoch in range(args.max_epochs):

        # TRAIN
        diagnostic_model.train()
        perform_iteration(diagnostic_model, train_dataloader, args.batch_size, device, sample_count=args.training_sample_count)

        # VALIDATION

        diagnostic_model.eval()
        validation_accuracies_meter, validation_losses_meter = perform_iteration(diagnostic_model, validation_dataloader, args.batch_size, device)

        new_best = False
        if validation_accuracies_meter.avg > best_accuracy:
            new_best = True
            best_accuracy = validation_accuracies_meter.avg
            save_model(model_path, diagnostic_model)

        print_results(validation_accuracies_meter, validation_losses_meter, epoch, args.max_epochs, "validation", new_best)

    best_model = initialize_model(args, device, model_path)
    best_model.eval()

    test_accuracies_meter, test_losses_meter = perform_iteration(best_model, test_dataloader, args.batch_size, device)
    print_results(test_accuracies_meter, test_losses_meter, epoch, args.max_epochs, "test", False)

if __name__ == "__main__":
    baseline(sys.argv[1:])

max_stars_repo_name: gautierdag/sygnal

max_issues_repo_name: gautierdag/sygnal

max_forks_repo_name: gautierdag/sygnal


text:
import math
from typing import TypeVar, Optional, Iterator
import os
import torch
from torch.utils.data.distributed import DistributedSampler
from torch.utils.data.dataset import Dataset
import torch.distributed as dist
import numpy as np
import random

T_co = TypeVar('T_co', covariant=True)

# This sampler is implemented for multi-GPU training, modified from the DistributedSampler of pytorch.
# The strategy follows VIGOR
class DistributedMiningSampler(DistributedSampler[T_co]):
    def __init__(self, dataset: Dataset, num_replicas: Optional[int] = None,
                 rank: Optional[int] = None, shuffle: bool = False,
                 seed: int = 0, drop_last: bool = False, batch_size: int = 128, mode = 'similarity', dim=1000, save_path = None) -> None:
        super(DistributedMiningSampler, self).__init__(dataset, num_replicas, rank, shuffle, seed, drop_last)
        self.dim = dim
        self.batch_size = batch_size * self.num_replicas
        self.queue_length = len(self.dataset)
        self.current_size = len(self.dataset) // self.batch_size * self.batch_size
        self.current_indices = np.arange(self.current_size)
        self.queue_size = 1 # for computing moving average, not used in this implementation
        self.queue = np.zeros([self.queue_length, self.queue_size, self.dim, 2])
        self.queue_ptr = 0
        self.queue_counter = np.zeros(self.queue_length,dtype=np.int)
        self.save_path = save_path
        self.mining_start = 1
        self.mining_pool_size = min(40000, self.queue_length)
        self.mining_save_size = 100
        self.choice_pool = range(self.mining_save_size)
        self.mining_save = np.zeros([self.queue_length, self.mining_save_size],dtype=int)

        self.mode = mode

    def update(self, data_sat, data_grd, indexes):
        data_sat_norm = data_sat / np.linalg.norm(data_sat, axis=1, keepdims=True)
        data_grd_norm = data_grd / np.linalg.norm(data_grd, axis=1, keepdims=True)
        batch_size = data_sat.shape[0]
        # writing in distributed training style, complicated. Update the queue according to the previous index.
        for j in range(self.num_replicas):
            index_j = self.indices_out[j:self.current_size:self.num_replicas]

            for i in range(batch_size // self.num_replicas):
                index = index_j[self.queue_ptr + i]
                assert index == indexes[i + j * (batch_size // self.num_replicas)]
                self.queue[index, self.queue_counter[index] % self.queue_size, :, 0] = \
                data_sat_norm[i + j * (batch_size // self.num_replicas)]
                self.queue[index, self.queue_counter[index] % self.queue_size, :, 1] = \
                    data_grd_norm[i + j * (batch_size // self.num_replicas)]
                self.queue_counter[index] += 1
        self.queue_ptr = (self.queue_ptr + batch_size // self.num_replicas)

    def generate_indices_sim(self):
        self.queue_ptr = 0

        random.seed(7 + self.epoch)
        self.current_indices = np.arange(self.current_size)
        random.shuffle(self.current_indices)

        if self.epoch >= self.mining_start:
            assert self.mining_pool_size <= self.queue_length
            mining_pool = np.array(random.sample(range(self.queue_length), self.mining_pool_size),dtype=int)
            product_train = np.matmul(self.queue[:,:,:,1].mean(axis=1), np.transpose(self.queue[mining_pool,:,:,0].mean(axis=1)))
            product_index = np.argsort(product_train, axis=1)
            ranking = np.zeros(product_train.shape[0])
            # update mining pool
            for i in range(product_train.shape[0]):
                ranking[i] = product_train.shape[1] - 1 - np.where(mining_pool[product_index[i]] == i)[0]
                self.mining_save[i, :] = mining_pool[product_index[i, -self.mining_save_size:]]
            # randomly sample first half
            ori_list = self.current_indices[:self.current_size//2]
            self.current_indices = []
            # global hard mining for the other half
            for i in range(self.current_size//self.batch_size):
                index_s = i * (self.batch_size//2)
                index_e = index_s + min(self.batch_size//2, self.current_size//2 - index_s)
                self.current_indices.extend(ori_list[index_s:index_e])
                hard_list = []
                for j in range(index_s, index_e):
                    idx = random.choice(self.mining_save[ori_list[j]])
                    # random sampling until no overlap in the batch
                    while idx in ori_list[index_s:index_e] or idx in hard_list:
                        idx = random.choice(self.mining_save[ori_list[j]])
                    hard_list.append(idx)
                self.current_indices.extend(hard_list)
        self.current_indices = np.array(self.current_indices, dtype=int)
        assert len(self.current_indices) == self.current_size
        print('sampler updated!')

    def update_epoch(self):
        self.generate_indices_sim()
        if self.rank == 0:
            np.save(os.path.join(self.save_path,'queue.npy'), self.queue)
            np.save(os.path.join(self.save_path,'queue_counter.npy'), self.queue_counter)

    def load(self, path):
        self.queue_counter = np.load(os.path.join(path,'queue_counter.npy'))
        self.queue = np.load(os.path.join(path,'queue.npy'))

    def __iter__(self) -> Iterator[T_co]:
        if self.shuffle:
            # deterministically shuffle based on epoch and seed
            g = torch.Generator()
            g.manual_seed(self.seed + self.epoch)
            indices = torch.randperm(len(self.current_indices), generator=g).tolist()  # type: ignore
        else:
            indices = list(range(len(self.current_indices)))  # type: ignore

        if not self.drop_last:
            # add extra samples to make it evenly divisible
            padding_size = self.current_size - len(indices)
            if padding_size <= len(indices):
                indices += indices[:padding_size]
            else:
                indices += (indices * math.ceil(padding_size / len(indices)))[:padding_size]
        else:
            # remove tail of data to make it evenly divisible.
            indices = indices[:self.current_size]
        assert len(indices) == self.current_size

        # subsample
        self.indices_out = self.current_indices[indices].tolist()
        indices = indices[self.rank:self.current_size:self.num_replicas]
        # assert len(indices) == self.num_samples
        # print(indices)
        indices_out = self.current_indices[indices].tolist()
        # print(self.rank, len(indices), len(indices_out))

        return iter(indices_out)

class DistributedMiningSamplerVigor(DistributedSampler[T_co]):
    def __init__(self, dataset: Dataset, num_replicas: Optional[int] = None,
                 rank: Optional[int] = None, shuffle: bool = False,
                 seed: int = 0, drop_last: bool = False, batch_size: int = 128, mode = 'similarity', dim=1000, save_path=None) -> None:
        super(DistributedMiningSamplerVigor, self).__init__(dataset, num_replicas, rank, shuffle, seed, drop_last)
        self.dim = dim
        self.batch_size = batch_size * self.num_replicas
        self.queue_length = max(dataset.train_data_size, len(dataset.train_sat_cover_list))
        self.current_size = len(self.dataset) // self.batch_size * self.batch_size
        self.current_indices = np.arange(self.current_size)
        self.queue_size = 1 # for computing moving average, not used in this implementation
        self.queue = np.zeros([self.queue_length, self.queue_size, self.dim, 2])
        self.queue_ptr = 0
        self.queue_counter = np.zeros(self.queue_length,dtype=np.int)
        self.save_path = save_path
        self.mining_start = 1
        self.mining_pool_size = min(40000, len(dataset.train_sat_cover_list))
        self.mining_save_size = 100
        self.choice_pool = range(self.mining_save_size)
        self.mining_save = np.zeros([self.queue_length, self.mining_save_size],dtype=int)
        self.mode = mode
        # raise Exception

    def update(self, data_sat, data_grd, indexes):
        data_sat_norm = data_sat / np.linalg.norm(data_sat, axis=1, keepdims=True)
        data_grd_norm = data_grd / np.linalg.norm(data_grd, axis=1, keepdims=True)
        batch_size = data_sat.shape[0]
        # writing in distributed training style, complicated. Update the queue according to the previous index.
        for j in range(self.num_replicas):
            index_j = self.indices_out[j:self.current_size:self.num_replicas]

            for i in range(batch_size // self.num_replicas):
                index = index_j[self.queue_ptr + i] %len(self.dataset.train_sat_cover_list)
                assert indexes[i + j * (batch_size // self.num_replicas)] in self.dataset.train_sat_cover_dict[self.dataset.train_sat_cover_list[index]]
                self.queue[index, self.queue_counter[index] % self.queue_size, :, 0] = \
                data_sat_norm[i + j * (batch_size // self.num_replicas)]
                self.queue[indexes[i + j * (batch_size // self.num_replicas)], self.queue_counter[index] % self.queue_size, :, 1] = \
                    data_grd_norm[i + j * (batch_size // self.num_replicas)]
                self.queue_counter[index] += 1
        self.queue_ptr = (self.queue_ptr + batch_size // self.num_replicas)

    def generate_indices_sim(self):
        self.queue_ptr = 0

        random.seed(7 + self.epoch)
        self.current_indices = np.arange(self.current_size) %len(self.dataset.train_sat_cover_list)
        random.shuffle(self.current_indices)

        if self.epoch >= self.mining_start:
            assert self.mining_pool_size <= self.queue_length
            mining_pool = np.array(random.sample(range(len(self.dataset.train_sat_cover_list)), self.mining_pool_size),dtype=int)
            product_train = np.matmul(self.queue[:,:,:,1].mean(axis=1), np.transpose(self.queue[mining_pool,:,:,0].mean(axis=1)))
            product_index = np.argsort(product_train, axis=1)
            # update mining pool
            for i in range(product_train.shape[0]):
                self.mining_save[i, :] = mining_pool[product_index[i, -self.mining_save_size:]]
            # randomly sample the first half
            ori_list = self.current_indices[:self.current_size//2]
            self.current_indices = []
            # global hard mining for the other half
            for i in range(self.current_size//self.batch_size):
                index_s = i * (self.batch_size//2)
                index_e = index_s + min(self.batch_size//2, self.current_size//2 - index_s)
                self.current_indices.extend(ori_list[index_s:index_e])
                hard_list = []
                for j in range(index_s, index_e):
                    grd_id = random.choice(self.dataset.train_sat_cover_dict[self.dataset.train_sat_cover_list[ori_list[j]]])
                    idx = int(random.choice(self.mining_save[grd_id]))
                    # keep random sampling until there is no overlap in the batch, hard coded as VIGOR is complicated
                    while True:
                        flag = False
                        for grd_idx in self.dataset.train_sat_cover_dict[self.dataset.train_sat_cover_list[idx]]:
                            if not self.dataset.check_overlap(ori_list[index_s:index_e],
                                                          grd_idx) or not self.dataset.check_overlap(hard_list, grd_idx):
                                flag = True
                        if flag:
                            idx = random.choice(self.mining_save[grd_id])
                        else:
                            break
                    hard_list.append(idx)
                self.current_indices.extend(hard_list)
        self.current_indices = np.array(self.current_indices, dtype=int)
        assert len(self.current_indices) == self.current_size
        print('sampler updated!')

    def update_epoch(self):
        # if self.epoch >= self.mining_start:
        self.generate_indices_sim()
        if self.rank == 0:
            np.save(os.path.join(self.save_path, 'queue.npy'), self.queue)
            np.save(os.path.join(self.save_path, 'queue_counter.npy'), self.queue_counter)

    def load(self, path):
        self.mining_start = 0
        self.queue_counter = np.load(os.path.join(path, 'queue_counter.npy'))
        self.queue = np.load(os.path.join(path, 'queue.npy'))

    def __iter__(self) -> Iterator[T_co]:
        if self.shuffle:
            # deterministically shuffle based on epoch and seed
            g = torch.Generator()
            g.manual_seed(self.seed + self.epoch)
            indices = torch.randperm(len(self.current_indices), generator=g).tolist()  # type: ignore
        else:
            indices = list(range(len(self.current_indices)))  # type: ignore

        if not self.drop_last:
            # add extra samples to make it evenly divisible
            padding_size = self.current_size - len(indices)
            if padding_size <= len(indices):
                indices += indices[:padding_size]
            else:
                indices += (indices * math.ceil(padding_size / len(indices)))[:padding_size]
        else:
            # remove tail of data to make it evenly divisible.
            indices = indices[:self.current_size]
        assert len(indices) == self.current_size

        # subsample
        self.indices_out = self.current_indices[indices].tolist()
        indices = indices[self.rank:self.current_size:self.num_replicas]
        # assert len(indices) == self.num_samples
        # print(indices)
        indices_out = self.current_indices[indices].tolist()
        # print(self.rank, len(indices), len(indices_out))

        return iter(indices_out)


text:
import numpy

from openff.recharge.grids import GridGenerator, GridSettings
from openff.recharge.utilities.openeye import smiles_to_molecule

def test_generate_fcc_grid():

    # Build a simple case monatomic test case.
    oe_molecule = smiles_to_molecule("[Ar]")
    conformer = numpy.array([[0.0, 0.0, 0.0]])

    # Select grid the grid settings so the corners of the FCC
    # lattice should be cut off.
    grid_settings = GridSettings(
        type="fcc",
        spacing=numpy.sqrt(2) * 1.8 / 2.0,
        inner_vdw_scale=0.9,
        outer_vdw_scale=1.1,
    )

    grid = GridGenerator.generate(oe_molecule, conformer, grid_settings)

    # Ensure that the grid encompasses the correct number of grid points.
    assert len(grid) == 24

max_stars_repo_name: openforcefield/openff-recharge

max_issues_repo_name: openforcefield/openff-recharge

max_forks_repo_name: openforcefield/openff-recharge


text:
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
from pprint import pprint
from datetime import datetime
from datetime import timedelta
import pickle
import copy
from pyhdf.SD import SD, SDC, SDS
from operator import add

##sampling times
sampling_times_file = 'C:/Users/Sarah Hanna/Documents/Data/WHI long term record/GOES-Chem/Junwei_runs/WHI_SP2_6h_rBC_mass_concs.txt'
sampling_times = []
with open(sampling_times_file,'r') as f:
	f.readline()
	for line in f:
		newline = line.split()
		sampling_date = newline[0]
		sampling_time = newline[1]
		sampling_datetime = datetime(int(sampling_date[0:4]),int(sampling_date[5:7]),int(sampling_date[8:10]),int(sampling_time[0:2]))
		sampling_datetime = datetime(int(sampling_date[0:4]),int(sampling_date[5:7]),int(sampling_date[8:10])) #without hour - date only
		sampling_times.append(sampling_datetime)

#GC data

GC_data = {}

data_dir = 'C:/Users/Sarah Hanna/Documents/Data/WHI long term record/GOES-Chem/Junwei_runs/no_biomass/'  # ['default','Vancouver_emission','wet_scavenging','no_biomass','All_together']

lat = 20 #20 corresponds to 50deg
lon = 7 #7 corresponds to -122.5deg
level = 9 #1-47 #9 is closest to WHI avg P (WHI 95% CI = 770-793)

molar_mass_BC = 12.0107 #in g/mol
ng_per_g = 10**9
R = 8.3144621 # in m3*Pa/(K*mol)
GEOS_Chem_factor = 10**-9
i=0
pressure = []

start_hour = 3
end_hour = 15

os.chdir(data_dir)
for file in os.listdir(data_dir):
	if file.endswith('.hdf'):

		file_year = int(file[2:6])
		file_month = int(file[6:8])
		file_day = int(file[8:10])
		file_hour = int(file[11:13])
		file_datetime = datetime(file_year,file_month,file_day,file_hour)

		#if (file_year == 2009 and file_month in [7,8]) or  (file_year == 2010 and file_month in [6,7]) or  (file_year == 2012 and file_month in [4,5]):

		if start_hour <= file_hour < end_hour:  #ignore any times not in the 2000-0800 PST window (0400-1600 UTC)
			hdf_file = SD(file, SDC.READ)

			#pprint(hdf_file.datasets())

			#pressures = hdf_file.select('PEDGE-$::PSURF')
			#pressure.append(pressures[level,lat,lon])
			#lats = hdf_file.select('LAT')
			#lons = hdf_file.select('LON')
			#print lats[lat], lons[lon]

			if start_hour <= file_hour < (start_hour+6):
				period_midtime = datetime(file_year,file_month,file_day,23) - timedelta(days=1) #this is the early night period of 2000-0200 PST (mid time is 2300 of the previous day when converting from UTC to PST)

			if (start_hour+6) <= file_hour < end_hour:
				period_midtime = datetime(file_year,file_month,file_day,05) #this is the late night period of 0200-0800 PST

			##switch for 24-hr averaging
			period_midtime = datetime(file_year,file_month,file_day) #this is the late night period of 0200-0800 PST

			hydrophilic_BC = hdf_file.select('IJ-AVG-$::BCPI') #3d conc data in ppbv (molBC/molAIR)
			hydrophobic_BC = hdf_file.select('IJ-AVG-$::BCPO')

			total_BC_ppbv = (hydrophilic_BC[level,lat,lon]) + (hydrophobic_BC[level,lat,lon])

			#hydrophilic_BC.__del__()
			#hydrophobic_BC.__del__()

			Factor =  (1000 * 1e2 * 1e6 * 1e-9) / (8.31 * 273)

			BC_conc_ngm3 = total_BC_ppbv*(molar_mass_BC*ng_per_g*GEOS_Chem_factor*(101325/(R*273)))  #101325/(R*273) corrects to STP
			temp_BC = total_BC_ppbv * Factor * 12 *1000

			if period_midtime in sampling_times:  #this excludes BB times already
				if period_midtime in GC_data:
					GC_data[period_midtime].append(temp_BC)
				else:
					GC_data[period_midtime] = [temp_BC]

			hdf_file.end()
			i+=1

##save stats to file
#os.chdir('C:/Users/Sarah Hanna/Documents/Data/WHI long term record/GOES-Chem/')
#file = open('WHI_long_term_GC_JW_stats_Vancouver_emission.txt', 'w')
#file.write('datetime' + '\t' +  'mean_bc_conc_ug/m3' +'\n')
#temp = []
#for period_midtime in GC_data:
#	mean_BC_conc = np.nanmean(GC_data[period_midtime])
#	row = [datetime(period_midtime.year,period_midtime.month,period_midtime.day), period_midtime.hour, mean_BC_conc]
#	temp.append(row)
#
#temp.sort()
#
#for row in temp:
#	line = '\t'.join(str(x) for x in row)
#	file.write(line + '\n')
#
#file.close()

#print np.mean(pressure)

GC_6h_all_non_BB = []
errs = []
list_of_all_hourly_concs = []

for period_midtime in GC_data:

	GC_concs = GC_data[period_midtime]
	mean_BC_conc = np.nanmean(GC_concs)
	GC_6h_all_non_BB.append(mean_BC_conc)

	for value in GC_concs:
		list_of_all_hourly_concs.append(value)

	level_std_err  = np.std(GC_data[period_midtime])   #standard error of the mean
	rel_err = level_std_err/mean_BC_conc

	errs.append(rel_err)

print np.percentile(GC_6h_all_non_BB, 10),np.percentile(GC_6h_all_non_BB, 50), np.percentile(GC_6h_all_non_BB, 90), np.mean(GC_6h_all_non_BB)
print np.percentile(list_of_all_hourly_concs, 10),np.percentile(list_of_all_hourly_concs, 50), np.percentile(list_of_all_hourly_concs, 90), np.mean(list_of_all_hourly_concs)
#print sum(errs)/len(errs)

####plots

fig = plt.figure()
bin_number_all_FT = 40
UL_all_FT = 300
bin_range_all_FT = (0,UL_all_FT)

ax1 = plt.subplot2grid((2,1), (0,0), colspan=1)
ax2 = plt.subplot2grid((2,1), (1,0), colspan=1)

#SP2
ax1.hist(GC_6h_all_non_BB,bins = bin_number_all_FT, range = bin_range_all_FT)
ax1.xaxis.set_visible(True)
ax1.yaxis.set_visible(True)
ax1.set_ylabel('frequency')
ax1.xaxis.tick_top()
ax1.xaxis.set_label_position('top')
ax1.xaxis.set_ticks(np.arange(0, UL_all_FT, 50))
ax1.axvline(np.percentile(GC_6h_all_non_BB,10), color= 'black', linestyle = '--')
ax2.set_xlabel('6h rBC mass concentration (ng/m3 - STP)')

#GC
ax2.hist(list_of_all_hourly_concs,bins = bin_number_all_FT, range = bin_range_all_FT, color = 'green')
ax2.xaxis.set_visible(True)
ax2.yaxis.set_visible(True)
ax2.set_ylabel('frequency')
ax2.xaxis.set_ticks(np.arange(0, UL_all_FT, 50))
ax2.axvline(np.percentile(list_of_all_hourly_concs,10), color= 'black', linestyle = '--')
ax2.set_xlabel('1h rBC mass concentration (ng/m3 - STP)')

#plt.savefig('histograms - GEOS-Chem and measurements - all non-BB FT - 6h - JW_data.png',bbox_inches='tight')

plt.show()


text:
# -*- coding: utf-8 -*-
"""
Created on Jan 24 2019

@author: Angelo Rossi, Ralph Florent
"""
# Import relevant libraries
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline

import numpy as np
import pandas as pd # we'll use pandas to load data
from scipy.interpolate import interp1d

# We'll be looking at phylosiclicates (clays), as they have a fairly complex spectra with many features:
# In the data_samples folder we have two spectra for Montmorillonite
PATH = '../assets/data/spectral_data/'
specSlPhyMontm = pd.read_csv(PATH + 'FRTC596_537-220-5x5_Al-phyllo_montm.txt', skipinitialspace=True, sep=' ',
                    skiprows=3, header=None)

relabMontm = pd.read_csv(PATH + 'montmorillonite-1292F35-RELAB.txt', skipinitialspace=True, sep=' ',
                      skiprows=3, header=None)
# also we have two spectra for Nontronite
specFePhyNontr = pd.read_csv(PATH + 'FRTC596_537-220-5x5_Fe-phyllo_nontr.txt', skipinitialspace=True, sep=' ',
                         skiprows=3, header=None)

relabNontr = pd.read_csv(PATH + 'nontronite-CBJB26-RELAB.txt', skipinitialspace=True, sep=' ',
                         skiprows=3, header=None)

specFePhyNontr[:5] # take a look at what our data looks like

# This worked, lets fix all dataframes now:
specSlPhyMontm = specSlPhyMontm.rename(columns={0:'wavelength',1:'specSlPhyMontm'})
relabMontm     = relabMontm.rename(columns={0:'wavelength',1:'relabMontm'})
specFePhyNontr = specFePhyNontr.rename(columns={0:'wavelength',1:'specFePhyNontr'})
relabNontr     = relabNontr.rename(columns={0:'wavelength',1:'relabNontr'})

#lets check:
relabNontr[:5]

#lets plot:
ax1 = specSlPhyMontm.plot.line(0, figsize=(15,5), title="Montmorillonite")
relabMontm.plot.line(0, ax=ax1) # ax=ax plots all on the same axis

ax2 = specFePhyNontr.plot.line(0, figsize=(15,5), title="Nontronite")
relabNontr.plot.line(0, ax=ax2)

ax3 = relabMontm.plot.line(0, figsize=(15,5), title="Montmorillonite vs Nontronite")
relabNontr.plot.line(0, ax=ax3)

plt.show()

#lets put this into a function so we can reuse it later:
def removeCont(pSample):
    pSampleLineX=[pSample[0][0],pSample[0][-1]]
    pSampleLineY=[pSample[1][0],pSample[1][-1]]
    pSampleLine=[pSampleLineX,pSampleLineY]
    finterp = interp1d(pSampleLine[0],pSampleLine[1])#create interploation function
    return pSample[1]-finterp(pSample[0])

# Lets try this on a bigger fake dataset:
pivot = lambda sample: [[a[0] for a in sample],[a[1] for a in sample]]

sample = np.asarray([list(x) for x in zip(range(32),[x[0] for x in np.random.random((32,1)).tolist()])])
pSample=pivot(sample)
plt.plot(pSample[0],pSample[1], '-',pSample[0],removeCont(pSample),'-')
plt.legend(['original', 'continuum removed'], loc='best')
plt.show()

#we can now easily find a maximum value above the continuum for points other than first and last:
max(removeCont(pSample)[1:-1])

#furthermore, we can ask numpy for an index of max:
maxIndex = np.argmax(removeCont(pSample)[1:-1]) +1
maxIndex

#Armed with this information we can repeat the previous step for a subset:
plt.plot(pSample[0][:maxIndex],pSample[1][:maxIndex], '-',pSample[0][:maxIndex],
         removeCont([pSample[0][:maxIndex],pSample[1][:maxIndex]]),'-')
plt.legend(['original', 'continuum removed'], loc='best')
plt.show()

# From henceforth we could do this recursively - that is to use the output of a function as an input to itself:
def getMaxima(pSample):
    def getMaximaInner(innerSample):
        contRem=removeCont(innerSample)
        #print(contRem)
        maxIndex=np.argmax(contRem)
        #print(maxIndex)
        maxVal=contRem[maxIndex]
        maxLoc=innerSample[0][maxIndex]
        if len(contRem)>2 and maxVal>contRem[0] and maxVal>contRem[-1]: # check that the maximum is more than edges
            maxLocArray.append(maxLoc)
            #print(maxLoc)
            subsetLeft=[innerSample[0][:maxIndex+1],innerSample[1][:maxIndex+1]]
            #print(subsetLeft[0])
            subsetRight=[innerSample[0][maxIndex:],innerSample[1][maxIndex:]]
            #print(subsetRight[0])
            getMaximaInner(subsetLeft)
            getMaximaInner(subsetRight)
    maxLocArray=[] #initialize array to store a list of points on a convex hull
    getMaximaInner(pSample)
    maxLocArray.sort()
    return [pSample[0][0]]+maxLocArray+[pSample[0][-1]]

#maxList=getMaxima([pSample[0][:5],pSample[1][:5]])
maxList=getMaxima(pSample)
print(maxList)

hull=[maxList,[x[1] for x in sample if x[0] in maxList]]
plt.plot(pSample[0],pSample[1], '-',hull[0],hull[1],'-')
plt.show()

#Now lets try it with a real dataset:
sample = np.asarray(relabMontm)
pSample = pivot(sample)
maxList = getMaxima(pSample)
print(maxList)
hull = [maxList,[x[1] for x in sample if x[0] in maxList]]
plt.plot(pSample[0],pSample[1], '-',hull[0],hull[1],'-')

# You may notice that the formation of a convex hull is distorted
# by a long row of zeros at the end of the data sample
# Lets remove all zeros:
cleanSample=[value for value in sample if value[1]>0]
pSample=pivot(cleanSample)
maxList=getMaxima(pSample)
#print(maxList)
pHull=[maxList,[x[1] for x in sample if x[0] in maxList]]
plt.plot(pSample[0],pSample[1], '-',pHull[0],pHull[1],'-')
plt.legend(['spectrum', 'convex hull'])
plt.show()

# Next we can subtact the convex hull from our data, in a manner similar to how we subtracted continuum ealier
def removeHull(pSample,pHull):
    finterp = interp1d(pHull[0],pHull[1])#create interploation function
    return pSample[1]-finterp(pSample[0])

hullRemoved = removeHull(pSample,pHull)
plt.plot(pSample[0],pSample[1],'-',pSample[0],hullRemoved, '-')
plt.legend(['spectrum', 'convex hull removed'])
plt.show()

# we can easily find indices for these values:
splitInd=[pSample[0].index(x) for x in pHull[0]]
print(splitInd)

#then we can split the array along the indices using list comprehension:
splitSample=[[pSample[0][splitInd[i]:splitInd[i+1]],hullRemoved[splitInd[i]:splitInd[i+1]]]
             for i in range(len(splitInd)-1) if splitInd[i+1]-splitInd[i]>2]
for s in splitSample:
    plt.plot(s[0],s[1],'-')
plt.show()

# Finding local minima is then straightforward:
listMinimaX=[x[0][np.argmin(x[1])] for x in np.asarray(splitSample)]
print(listMinimaX)

# we can use list comprehension again to get corresponding Y-values
listMinimaY=[pSample[1][pSample[0].index(x)] for x in listMinimaX]
print(listMinimaY)

#And we can combine the two and plot the minima on a graph:
splitSample=[[pSample[0][splitInd[i]:splitInd[i+1]],pSample[1][splitInd[i]:splitInd[i+1]]]
             for i in range(len(splitInd)-1) if splitInd[i+1]-splitInd[i]>2]
for s in splitSample:
    plt.plot(s[0],s[1],'-')
plt.plot(listMinimaX,listMinimaY,'x',color='black')
plt.show()

# We've now identified some deep absorption bands, a some shallow that are probably noise.
# We can filter out the shallow ones by appling a threshold:
# First get the band depths with hull removed:
listMinimaYhullRemoved=[hullRemoved[pSample[0].index(x)] for x in listMinimaX]
print(listMinimaYhullRemoved)

# Now apply a threshold:
threshold=0.05
listMinimaSigX=[q[0] for q in list(zip(listMinimaX,listMinimaYhullRemoved)) if q[1]<-threshold]
listMinimaSigYhullRemoved=[q[1] for q in list(zip(listMinimaX,listMinimaYhullRemoved)) if q[1]<-threshold]
listMinimaSigY=[pSample[1][pSample[0].index(x)] for x in listMinimaSigX]

#then we can split the array along the indices using list comprehension:
splitSample=[[pSample[0][splitInd[i]:splitInd[i+1]],hullRemoved[splitInd[i]:splitInd[i+1]]]
             for i in range(len(splitInd)-1) if splitInd[i+1]-splitInd[i]>2]
plt.figure(figsize=(15,5)) #make larger figure
for s in splitSample:
    plt.plot(s[0],s[1],'-')
for xc in listMinimaSigX:
    plt.axvline(x=xc,color='black')
plt.show()

# Finally, lets see where do these band depths plot on original spectra:
plt.figure(figsize=(15,5))
splitSample=[[pSample[0][splitInd[i]:splitInd[i+1]],pSample[1][splitInd[i]:splitInd[i+1]]]
             for i in range(len(splitInd)-1) if splitInd[i+1]-splitInd[i]>2]
for s in splitSample:
    plt.plot(s[0],s[1],'-')
plt.plot(listMinimaSigX,listMinimaSigY,'o',color='black')
plt.show()

print("adsorption band center wavelenghts are:")
for item in listMinimaSigX:
    print(item, "micrometers")

max_stars_repo_name: ralflorent/owm

max_issues_repo_name: ralflorent/owm

max_forks_repo_name: ralflorent/owm

reason: import numpy,from scipy


text:
import h5py
import numpy as np
from PostProcessing_HDF_Data.Process_LBEPR_Table_Class import Process_LBEPR_Table

def CreateHDF5Dataset_ChChGen(hdf5_file):
    """ Create an array which you will append to the existing 'LBEPR' table.
        shape = len(lbepr, 3); [child_1, child_2, generation]
    """

    # Call the class & its relevant functions:
    call = Process_LBEPR_Table(hdf5_file=hdf5_file, channel="GFP")
    children = call.FindBothChildrenOfParent()
    generation = call.FindGenerationalDepth()

    # Create a list of lists with data, return in np.array form:
    maxi_list = []
    for child_list, gen in zip(children, generation):
        child_list.append(gen)
        maxi_list.append(child_list)

    chchtime = np.array(maxi_list, dtype=np.int32)
    return chchtime

def AppendChChTime_ToLBEPR(hdf5_file):
    """ "LBEPR" shape (X, 5)    ->    "LBEPRChChGen" shape (X, 8) """

    chchtime_array = CreateHDF5Dataset_ChChGen(hdf5_file=hdf5_file)

    # Open the HDF5 file:
    f = h5py.File(hdf5_file, 'a')
    lbepr_array = list(f["tracks"]["obj_type_1"]["LBEPR"])

    lst = []
    for item, add in zip(lbepr_array, chchtime_array):
        lst.append(np.append(item, add))
    full_array = np.array(lst, dtype=np.int32)

    if len(lbepr_array) != full_array.shape[0]:
        raise ValueError("The shapes of LBEPR and LBEPRChChGen arrays are different.")

    # Append to HDF5:
    if "LBEPRChChGen" in list(f["tracks"]["obj_type_1"]):
        del f["tracks"]["obj_type_1"]["LBEPRChChGen"]

    grp = f["tracks"]["obj_type_1"]
    grp.create_dataset(name="LBEPRChChGen", data=full_array)

    if f.__bool__():
        f.close()

def DeleteLBEPR_KeepLBEPRChChGen(hdf5_file):
    """ Not necessary to call. But just in case you don't want data duplicates. """

    # Open the HDF5 file:
    f = h5py.File(hdf5_file, 'a')
    lbepr = f["tracks"]["obj_type_1"]["LBEPR"]
    lbeprchchgen = f["tracks"]["obj_type_1"]["LBEPRChChGen"]

    # Delete only if replicas of each other:
    if lbepr.shape[0] == lbeprchchgen.shape[0]:
        if lbepr.shape[1] == 5 and lbeprchchgen.shape[1] == 8:
            del f["tracks"]["obj_type_1"]["LBEPR"]

max_stars_repo_name: The-Kristina/CellComp

max_stars_repo_licenses: CNRI-Python RSA-MD Xnet Net-SNMP X11

max_issues_repo_name: The-Kristina/CellComp

max_issues_repo_licenses: CNRI-Python RSA-MD Xnet Net-SNMP X11

max_forks_repo_name: The-Kristina/CellComp

max_forks_repo_licenses: CNRI-Python RSA-MD Xnet Net-SNMP X11