text:
"""Cail's mcmc.py run wrapper.

Note that this is all set up for Py3k (in argparse) so it might take a bunch of
rewriting to get this down to Python2.
"""

import numpy as np
import argparse
import subprocess as sp
import os
from astropy.io import fits
from collections import OrderedDict
import matplotlib.pyplot as plt
import copy

import fitting
import plotting
import mcmc
from disk_model import debris_disk, raytrace
import aumic_fitting

plt.switch_backend('agg')
run_name = 'run25'

# A dictionary of default param vals (changed in make_best_fits())
param_dict = OrderedDict([
    ('temp_index',        -0.5),
    ('m_disk',            -7.54),
    ('sb_law',             2),
    ('r_in',               20),
    ('d_r',                22),
    ('r_crit',            150.0),
    ('inc',               88.6),
    ('m_star',            0.31),
    ('co_frac',           1e-4),
    ('v_turb',            0.081),
    ('Zq',                70.0),
    ('column_densities', [0.79, 1000]),
    ('abundance_bounds', [50, 500]),
    ('hand',              -1),
    ('rgrid_size',        1000),
    ('zgrid_size',        500),
    ('l_star',            0.09),
    ('scale_factor',      0.031),
    ('pa',                128.49),
    ('mar_starflux',      3.90e-4),
    ('aug_starflux',      1.50e-4),
    ('jun_starflux',      2.20e-4)])

def main():
    """Establish and evaluate some custom argument options.

    When is this actually called?
    """
    parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter, description='''Python commands associated with emcee run25, which has 50 walkers and varies the following parameters:
    1)  disk mass
    2)  surface brightness power law exponent
    3)  scale factor, multiplied by radius to get scale height
    4)  inner radius
    5)  outer radius (really inner radius + dr)
    6)  inclination
    7)  position angle
    8)  march starflux
    9)  august starflux
    10) june starflux
This run is a redo of the previous fiducial model run19. This run uses the new Gaia distance to AU Mic and a new gridding resolution.''')

    parser.add_argument('-r', '--run', action='store_true',
                        help='begin or resume eemcee run.')

    parser.add_argument('-a', '--analyze', action='store_true',
                        help='analyze sampler chain, producing an evolution plot, corner plot, and image domain figure.')

    parser.add_argument('-b', '--burn_in', default=0, type=int,
                        help='number of steps \'burn in\' steps to exclude')

    parser.add_argument('-bf', '--best_fit', action='store_true',
                        help='generate best fit model images and residuals')

    parser.add_argument('-con', '--concise', action='store_true',
                        help='concise best fit')

    parser.add_argument('-c', '--corner', action='store_true',
                        help='generate corner plot')

    parser.add_argument('-cvars', '--corner_vars', default=None, nargs='+',
                        help='concise best fit')

    parser.add_argument('-e', '--evolution', action='store_true',
                        help='generate walker evolution plot.')

    parser.add_argument('-kde', '--kernel_density', action='store_true',
                        help='generate kernel density estimate (kde) of posterior distribution')

    args = parser.parse_args()

    if args.run:
        """Note that to_vary is of form:
                [param name,
                 initial_position_center,
                 initial_position_sigma,
                 (prior low bound, prior high bound)
                 ]
        """
        mcmc.run_emcee(run_name=run_name, nsteps=10000, nwalkers=50,
            lnprob = lnprob,
            to_vary = [
            ('m_disk',            -7.54,        0.06,      (-np.inf, np.inf)),
            ('sb_law',            1,            0.5,       (-np.inf, np.inf)),
            ('scale_factor',      0.031,        0.005,     (0,       np.inf)),
            ('r_in',              18,           5,         (0,       np.inf)),
            ('d_r',               24.3,         5,         (0,       np.inf)),
            ('inc',               88.6,         0.4,       (0,       90.)),
            ('pa',                128.49,       0.07,      (0,       360)),
            ('mar_starflux',      3.9e-4,       0.2e-4,    (0,       np.inf)),
            ('aug_starflux',      1.5e-4,       0.2e-4,    (0,       np.inf)),
            ('jun_starflux',      2.2e-4,       0.2e-4,    (0,       np.inf))])
    else:
        run = mcmc.MCMCrun(run_name, nwalkers=50, burn_in=args.burn_in)
        # old_nsamples = run.groomed.shape[0]
        # run.groomed = run.groomed[run.groomed['r_in'] + run.groomed['d_r'] > 20]
        # print('{} samples removed.'.format(old_nsamples - run.groomed.shape[0]))

        # This was down in the arg analysis but seems separate.
        # Also, aumic_fitting doesn't seem to exist anymore?
        aumic_fitting.label_fix(run)

        # Read the arguments passed and execute them.
        if args.corner_vars:
            cols = list(run.groomed.columns)
            col_indices = [cols.index(col) for col in args.corner_vars]

        if args.analyze or args.best_fit:
            make_best_fits(run, concise=args.concise)

        if args.corner_vars:
            args.corner_vars = run.groomed.columns[col_indices]

        if args.analyze or args.evolution:
            run.evolution()

        if args.analyze or args.kernel_density:
            run.kde()

        if args.analyze or args.corner:
            run.corner(variables=args.corner_vars)

def make_fits(model, disk_params):
    """Run Kevin's code to make a model disk.

    Note that this is debris-disk specific!
    """
    structure_params = disk_params[:-1]
    PA = disk_params[-1]
    freq0 = model.observations[0][0].uvf[0].header['CRVAL4']*1e-9

    # Make a sky-projected debris disk.
    model_disk = debris_disk.Disk(structure_params, obs=[300, 351, 300, 5])
    raytrace.total_model(model_disk,
                         distance=9.725,                # pc
                         imres=0.015,                   # arcsec/pix
                         xnpix=1024,                    # image size in pixels
                         freq0=freq0,                   # obs frequency
                         PA=PA,
                         offs=[0.0, 0.0],               # from image center
                         nchans=1,                      # continum
                         isgas=False,                   # continuum!
                         includeDust=True,              # continuuum!!
                         extra=0,                       # ?
                         modfile=model.root + model.name)

def fix_fits(model, obs, starflux):
    """Open fits and add starflux.

    Seems like this is maybe specific to flair-fixing.
    """
    model_fits = fits.open(model.path + '.fits')
    crpix = int(model_fits[0].header['CRPIX1'])
    model_im = model_fits[0].data[0]
    try:
        model_im[crpix, crpix] = model.crpix_diskflux + starflux

    except AttributeError:
        model.crpix_diskflux = model_im[crpix, crpix]
        model_im[crpix, crpix] += starflux

    model_fits[0].header['CRVAL1'] = obs.ra
    model_fits[0].header['CRVAL2'] = obs.dec

    model_fits.writeto(model.path + '.fits', overwrite=True)

# Define likelehood functions
def lnprob(theta, run_name, to_vary):
    """Calculate the log likelihood of the model based on its chi2.

    For each parameter to be varied, return -inf if the proposed value lies
    outside the bounds for the parameter; if all parameters are OK, return 0.
    """
    # Get the priors and make sure the model value falls within them.
    for i, free_param in enumerate(to_vary):
        lower_bound, upper_bound = free_param[-1]

        if lower_bound < theta[i] < upper_bound:
            param_dict[free_param[0]] = theta[i]
        else:
            return -np.inf

    param_dict_copy = copy.copy(param_dict)
    # Note that original disk mass was an exponent val!
    param_dict_copy['m_disk'] = 10**param_dict_copy['m_disk']
    param_dict_copy['d_r'] += param_dict_copy['r_in']

    # What is this?
    disk_params = param_dict_copy.values()[:-3]
    starfluxes = param_dict_copy.values()[-3:]

    # Make model and the resulting fits image
    model = fitting.Model(observations=aumic_fitting.band6_observations,
                          root=run_name + '/model_files/',
                          name='model' + str(np.random.randint(1e10)))
    make_fits(model, disk_params)

    # I don't understand why this is complicated.
    # Probably something to do with the flares?
    for pointing, starflux in zip(model.observations, starfluxes):
        for obs in pointing:
            fix_fits(model, obs, starflux)
            model.obs_sample(obs)
            model.get_chi(obs)

    model.delete()
    return -0.5 * sum(model.chis)

def make_best_fits(run, concise=False):
    """Docstring.

    Args:
        run (what?): what kind of object is this?
        concise (bool): where does this come up?
    """
    subset_df = run.main  # [run.main['r_in'] < 15]

    # Locate the best fit model from max'ed lnprob.
    max_lnprob = subset_df['lnprob'].max()
    model_params = subset_df[subset_df['lnprob'] == max_lnprob].drop_duplicates()
    print 'Model parameters:\n', model_params.to_string(), '\n\n'

    for param in model_params.columns[:-1]:
        param_dict[param] = model_params[param].values

    param_dict['m_disk'] = 10**param_dict['m_disk']
    param_dict['d_r'] += param_dict['r_in']

    # Why are we picking up the last three?
    disk_params = param_dict.values()[:-3]
    starfluxes = param_dict.values()[-3:]

    # intialize model and make fits image
    print('Making model...')
    model = fitting.Model(observations=aumic_fitting.band6_observations,
                          root=run.name + '/model_files/',
                          name=run.name + '_bestfit')
    make_fits(model, disk_params)

    print('Sampling and cleaning...')
    paths = []
    for pointing, rms, starflux in zip(model.observations, aumic_fitting.band6_rms_values[:-1], starfluxes):
        ids = []
        for obs in pointing:
            fix_fits(model, obs, starflux)

            ids.append('_' + obs.name[12:20])
            model.obs_sample(obs, ids[-1])
            model.make_residuals(obs, ids[-1])

        cat_string1 = ','.join([model.path+ident+'.vis' for ident in ids])
        cat_string2 = ','.join([model.path+ident+'.residuals.vis' for ident in ids])
        paths.append('{}_{}'.format(model.path, obs.name[12:15]))

        sp.call(['uvcat', 'vis={}'.format(cat_string2), 'out={}.residuals.vis'.format(paths[-1])], stdout=open(os.devnull, 'wb'))

        sp.call(['uvcat', 'vis={}'.format(cat_string1), 'out={}.vis'.format(paths[-1])], stdout=open(os.devnull, 'wb'))

        model.clean(paths[-1] + '.residuals', rms, show=False)
        model.clean(paths[-1], rms, show=False)

    cat_string1 = ','.join([path + '.vis' for path in paths])
    cat_string2 = ','.join([path + '.residuals.vis' for path in paths])

    sp.call(['uvcat', 'vis={}'.format(cat_string1), 'out={}_all.vis'.format(model.path)], stdout=open(os.devnull, 'wb'))

    sp.call(['uvcat', 'vis={}'.format(cat_string2), 'out={}_all.residuals.vis'.format(model.path)], stdout=open(os.devnull, 'wb'))

    model.clean(model.path+'_all', aumic_fitting.band6_rms_values[-1], show=False)
    model.clean(model.path+'_all.residuals', aumic_fitting.band6_rms_values[-1], show=False)

    paths.append('{}_all'.format(model.path))

    print('Making figure...')
    if concise:
        # Concise means to just look at the last model (?)
        fig = plotting.Figure(layout=(1, 3),
                              paths=[aumic_fitting.band6_fits_images[-1],
                                     paths[-1] + '.fits',
                                     paths[-1] + '.residuals.fits'],
                              rmses=3*[aumic_fitting.band6_rms_values[-1]],
                              texts=[[[4.6, 4.0, 'Data']],
                                     [[4.6, 4.0, 'Model']],
                                     [[4.6, 4.0, 'Residuals']]
                                     ],
                              title=r'Run 6 Global Best Fit Model & Residuals',
                              savefile=run.name + '/' + run.name + '_bestfit_concise.pdf',
                              show=True)
    else:
        # Note that paths and texts are just long, ugly list comps
        fig = plotting.Figure(layout=(4, 3),
                              paths=[[obs,
                                      path + '.fits',
                                      path + '.residuals.fits']
                                     for obs, path in zip(aumic_fitting.band6_fits_images, paths)],
                              rmses=[3*[rms] for rms in aumic_fitting.band6_rms_values],
                              texts=[[[[4.6, 4.0, date]],
                                      [[4.6, 4.0, 'rms={}'.format(np.round(rms*1e6))]],
                                      None]
                                     for date, rms in zip(['March', 'August', 'June', 'All'],
                                     aumic_fitting.band6_rms_values)
                                     ],
                              title=run.name + r'Global Best Fit Model & Residuals',
                              savefile=run.name+'/' + run.name + '_bestfit_global.pdf')

if __name__ == '__main__':
    main()

max_stars_repo_name: Jonasori/astrocail

max_issues_repo_name: Jonasori/astrocail

max_forks_repo_name: Jonasori/astrocail

reason: import numpy,from astropy


text:
# -*- coding: utf-8 -*-

from pyramid.view import view_config
import networkx as nx
from itertools import chain
from pyramid.httpexceptions import HTTPInternalServerError, HTTPNotFound

from chsdi.models.bod import Catalog
from chsdi.lib.validation import MapNameValidation
from chsdi.lib.filters import filter_by_geodata_staging

# Custom networkx like function
# Inspired by https://github.com/networkx/networkx/blob/master/networkx/readwrite/json_graph/tree.py#L16
def tree_data(G, root, attrs, meta):
    if G.number_of_nodes() != G.number_of_edges() + 1:
        raise TypeError("G is not a tree.")
    if not G.is_directed():
        raise TypeError("G is not directed.")

    id_ = attrs['id']
    children = attrs['children']
    if id_ == children:
        raise nx.NetworkXError('Attribute names are not unique.')

    def add_children(n, G):
        distinct_order_keys = set()
        order_key = lambda x: x['orderKey']
        label_key = lambda x: x['label'].replace(u'Ü', u'U').replace(u'Ä', u'A')
        nbrs = G[n]
        if len(nbrs) == 0:
            return []
        children_ = []
        nbrs = sorted(nbrs)
        for child in nbrs:
            d = dict(chain(G.node[child].items(), [(id_, child)]))
            d['category'] = meta[d['id']]['category']
            d['staging'] = meta[d['id']]['staging']
            d['label'] = meta[d['id']]['label']
            # only for sorting according to 'orderKey'
            if 'orderKey' in meta[d['id']]:
                d['orderKey'] = meta[d['id']]['orderKey']
            if (d['category'] == 'layer'):
                d['layerBodId'] = meta[d['id']]['layerBodId']
            else:
                d['selectedOpen'] = meta[d['id']]['selectedOpen']
            c = add_children(child, G)
            if c:
                d[children] = c
            elif d['category'] != 'layer':
                d[children] = []
            children_.append(d)
            if 'orderKey' in d:
                distinct_order_keys.add(d['orderKey'])

        if len(distinct_order_keys) > 1:
            children_ = sorted(children_, key=order_key)
        else:
            children_ = sorted(children_, key=label_key)

        for d in children_:
            d.pop('orderKey', None)
        return children_

    data = dict(chain(G.node[root].items(), [(id_, root)]))
    data[children] = add_children(root, G)
    return data

def initialize_graph(G, rows, lang):
    root_id = None
    meta = {}
    for r in rows:
        node_parent_id = getattr(r, 'parentId')
        node_id = getattr(r, 'id')
        if node_id not in meta:
            meta[node_id] = r.to_dict(lang)
        G.add_edge(node_parent_id, node_id)
        if getattr(r, 'category') == 'root':
            if (root_id is not None):
                raise HTTPInternalServerError('%s has more than 1 root_id, there can be only one!' % r.topic)
            root_id = node_id
    if root_id is None:
        raise HTTPInternalServerError('%s is missing root_id!' % r.topic)
    return G, meta, root_id

def create_digraph(rows, lang):
    G = nx.DiGraph()
    graph, meta, root_id = initialize_graph(G, rows, lang)
    edge_list = list(nx.bfs_edges(graph, root_id))
    G = nx.DiGraph(edge_list)
    return G, meta, root_id

class CatalogService(MapNameValidation):

    def __init__(self, request):
        self.lang = request.lang
        self.mapName = request.matchdict.get('map')  # The topic
        self.hasMap(request.db, self.mapName)
        self.request = request
        self.staging = request.registry.settings['geodata_staging']

    @view_config(route_name='catalog', renderer='jsonp')
    def catalog(self):
        model = Catalog
        query = self.request.db.query(model)\
            .filter(model.topic.like('%%%s%%' % self.mapName.lower()))
        rows = filter_by_geodata_staging(query, model.staging, self.staging).all()
        if len(rows) == 0:
            raise HTTPNotFound('No catalog with id %s is available' % self.mapName)

        G, meta, root_id = create_digraph(rows, self.lang)

        if len(rows) != len(G.nodes()):
            raise HTTPInternalServerError('Catalog tree for topic %s has unconnected leaves' % self.mapName)
        data = tree_data(G, root_id, {'children': 'children', 'id': 'id'}, meta)
        data['category'] = meta[root_id]['category']
        data['staging'] = meta[root_id]['staging']
        data = {'results': {'root': data}}
        return data

ext: py

reason: import networkx


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

import numpy as np

from examples.seismic import demo_model
from examples.seismic import AcquisitionGeometry
from examples.seismic.acoustic import AcousticWaveSolver
from devito import Function
from examples.seismic import Receiver
from distributed import Client, wait
from ctypes import c_int, c_float, c_bool
from distributed import LocalCluster

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable

from sotb_wrapper import interface

# Serial modeling function
def forward_modeling_single_shot(model, geometry, save=False, dt=4.0):
    solver = AcousticWaveSolver(model, geometry, space_order=4)
    d_obs, u0 = solver.forward(vp=model.vp, save=save)[0:2]
    return d_obs.resample(dt), u0

# Parallel modeling function
def forward_modeling_multi_shots(client, model, geometry, save=False, dt=4.0):

    futures = []
    for i in range(geometry.nsrc):

        # Geometry for current shot
        geometry_i = AcquisitionGeometry(model, geometry.rec_positions,
                                         geometry.src_positions[i, :],
                                         geometry.t0, geometry.tn, f0=geometry.f0,
                                         src_type=geometry.src_type)

        # Call serial modeling function for each index
        futures.append(client.submit(forward_modeling_single_shot, model,
                       geometry_i, save=save, dt=dt))

    # Wait for all workers to finish and collect shots
    wait(futures)
    shots = []
    for i in range(geometry.nsrc):
        shots.append(futures[i].result()[0])

    return shots

# Serial FWI objective function
def fwi_objective_single_shot(model, geometry, d_obs):

    # Devito objects for gradient and data residual
    grad = Function(name="grad", grid=model.grid)
    residual = Receiver(name='rec', grid=model.grid,
                        time_range=geometry.time_axis,
                        coordinates=geometry.rec_positions)
    solver = AcousticWaveSolver(model, geometry, space_order=4)

    # Predicted data and residual
    d_pred, u0 = solver.forward(vp=model.vp, save=True)[0:2]
    dobs_resampled = d_obs.resample(geometry.dt).data[:][0:d_pred.data.shape[0], :]
    residual.data[:] = d_pred.data[:] - dobs_resampled

    # Function value and gradient
    fval = .5*np.linalg.norm(residual.data.flatten())**2
    solver.gradient(rec=residual, u=u0, vp=model.vp, grad=grad)

    # Convert to numpy array and remove absorbing boundaries
    grad_crop = np.array(grad.data[:])[model.nbl:-model.nbl, model.nbl:-model.nbl]

    return fval, grad_crop

# Parallel FWI objective function
def fwi_objective_multi_shots(client, model, geometry, d_obs):

    futures = []
    for i in range(geometry.nsrc):

        # Geometry for current shot
        geometry_i = AcquisitionGeometry(model, geometry.rec_positions,
                                         geometry.src_positions[i, :], geometry.t0,
                                         geometry.tn, f0=geometry.f0,
                                         src_type=geometry.src_type)

        # Call serial FWI objective function for each shot location
        futures.append(client.submit(fwi_objective_single_shot, model, geometry_i,
                       d_obs[i]))

    # Wait for all workers to finish and collect function values and gradients
    wait(futures)
    fval = 0.0
    grad = np.zeros(model.shape)
    for i in range(geometry.nsrc):
        fval += futures[i].result()[0]
        grad += futures[i].result()[1]

    return fval, grad

# Wrapper for scipy optimizer: x is current model in squared slowness [s^2/km^2]
def loss(c, x, model, geometry, d_obs):

    # Convert x to velocity
    v_curr = 1.0/np.sqrt(x.reshape(model.shape))

    # Overwrite current velocity in geometry (don't update boundary region)
    model.update('vp', v_curr.reshape(model.shape))

    # Evaluate objective function
    fval, grad = fwi_objective_multi_shots(c, model, geometry, d_obs)
    return c_float(fval), grad.flatten().astype(np.float32)

def main(c):
    # Set up velocity model
    shape = (101, 101)    # Number of grid points (nx, nz).
    spacing = (10., 10.)  # Grid spacing in m. The domain size is now 1km by 1km.
    origin = (0, 0)  # Need origin to define relative source and receiver locations.
    nbl = 40

    # True model
    model1 = demo_model('circle-isotropic', vp_circle=3.0, vp_background=2.5,
                        origin=origin, shape=shape, spacing=spacing, nbl=nbl)

    # Initial model
    model0 = demo_model('circle-isotropic', vp_circle=2.5, vp_background=2.5,
                        origin=origin, shape=shape, spacing=spacing, nbl=nbl,
                        grid=model1.grid)

    # Set up acquisiton geometry
    t0 = 0.
    tn = 1000.
    f0 = 0.010

    # Set up source geometry, but define 5 sources instead of just one.
    nsources = 5
    src_coordinates = np.empty((nsources, 2))
    src_coordinates[:, 1] = np.linspace(0, model1.domain_size[0], num=nsources)
    src_coordinates[:, 0] = 20.  # Source depth is 20m

    # Initialize receivers for synthetic and imaging data
    nreceivers = 101
    rec_coordinates = np.empty((nreceivers, 2))
    rec_coordinates[:, 1] = np.linspace(spacing[0], model1.domain_size[0] -
                                        spacing[0], num=nreceivers)
    rec_coordinates[:, 0] = 980.    # Receiver depth
    # Set up geometry objects for observed and predicted data
    geometry1 = AcquisitionGeometry(model1, rec_coordinates, src_coordinates, t0,
                                    tn, f0=f0, src_type='Ricker')
    geometry0 = AcquisitionGeometry(model0, rec_coordinates, src_coordinates, t0,
                                    tn, f0=f0, src_type='Ricker')

    # Compute observed data in parallel (inverse crime). In real life we would
    # read the SEG-Y data here.
    d_obs = forward_modeling_multi_shots(c, model1, geometry1, save=False)

    # Box contraints
    vmin = 1.4    # do not allow velocities slower than water
    vmax = 4.0
    lb = [1.0/vmax**2 for _ in range(np.prod(model0.shape))]  # in [s^2/km^2]
    ub = [1.0/vmin**2 for _ in range(np.prod(model0.shape))]  # in [s^2/km^2]
    lb = np.array(lb).astype(np.float32)
    ub = np.array(ub).astype(np.float32)

    # Create an instance of the SEISCOPE optimization toolbox (sotb) Class.
    sotb = interface.sotb_wrapper()

    # Set some fields of the UserDefined derived type in Fortran (ctype structure).
    # parameter initialization
    n = c_int(shape[0]*shape[1])    # dimension
    flag = c_int(0)    			    # first flag
    sotb.udf.conv = c_float(1e-8)   # tolerance for the stopping criterion
    sotb.udf.print_flag = c_int(1)  # print info in output files
    sotb.udf.debug = c_bool(False)  # level of details for output files
    sotb.udf.niter_max = c_int(5)   # maximum iteration number
    sotb.udf.l = c_int(5)

    # Print the derived type.
    print('Hello from Python!')
    print(sotb.udf)

    # Initial guess
    v0 = model0.vp.data[model0.nbl:-model0.nbl, model0.nbl:-model0.nbl]
    X = 1.0 / (v0.reshape(-1).astype(np.float32))**2

    # computation of the cost and gradient associated
    # with the initial guess
    fcost, grad = loss(c, X, model0, geometry0, d_obs)

    while (flag.value != 2 and flag.value != 4):
        sotb.LBFGS(n, X, fcost, grad, flag, lb, ub)

        if (flag.value == 1):
            # compute cost and gradient at point x
            fcost, grad = loss(c, X, model0, geometry0, d_obs)

    # Helpful console writings
    print('END OF TEST')
    print('FINAL iterate is : ', X)
    print('See the convergence history in iterate_LB.dat')

    # Plot FWI result
    vp = 1.0/np.sqrt(X.reshape(model1.shape))
    vmin = 2.4
    vmax = 2.8

    mpl.rcParams['font.size'] = 8.5
    fig, (ax1, ax2) = plt.subplots(1, 2)
    fig.subplots_adjust(wspace=0.25)
    label_format = '{:,.1f}'
    #
    for count, ax in enumerate([ax1, ax2]):
        if count == 0:
            im = ax.imshow(model1.vp.data[nbl:-nbl, nbl:-nbl].T,
                           cmap=plt.cm.cividis, vmin=vmin, vmax=vmax)
        else:
            im = ax.imshow(vp.T, cmap=plt.cm.cividis, vmin=vmin, vmax=vmax)
        ax_divider = make_axes_locatable(ax)
        cax = ax_divider.append_axes("right", size="7%", pad="2%")
        cb = plt.colorbar(im, cax=cax)
        cb.ax.tick_params(labelsize=8)
        #
        ticks_ylabels = (ax.get_yticks()*0.01).tolist()
        ticks_yloc = ax.get_yticks().tolist()
        ticks_xlabels = (ax.get_xticks()*0.01).tolist()
        ticks_xloc = ax.get_xticks().tolist()
        ax.yaxis.set_major_locator(mticker.FixedLocator(ticks_yloc))
        ax.xaxis.set_major_locator(mticker.FixedLocator(ticks_xloc))
        ax.set_yticklabels([label_format.format(y) for y in ticks_ylabels])
        ax.set_xticklabels([label_format.format(x) for x in ticks_xlabels])
        ax.set(xlabel='x (km)', ylabel='Depth (km)')
        ax.label_outer()

    plt.savefig('circle_isotropic_inversion.pdf')

if __name__ == "__main__":
    r'''
    This script demonstrates how we can set up a basic FWI framework
    with gradient-based optimization algorithms from the SEISCOPE
    optimization toolbox (sotb) wrapper. The script is basically a copy
    of the devito FWI tutorial (https://github.com/devitocodes/devito/
    blob/master/examples/seismic/tutorials/04_dask.ipynb) with the addition
    of the sotb wrapper. It uses a simple toy example for validation of the
    code.
    '''
    print('start')
    # Start Dask cluster
    cluster = LocalCluster(n_workers=5, death_timeout=600, asynchronous=False)
    c = Client(cluster)
    main(c)

reason: import numpy


text:
import numpy as np
def rejection_sampling(tri, alpha_boolean, test_data):
    # Might be beneficial to have a seperate section for rejection sampling types
    acceptance_boolean = np.zeros(len(test_data), dtype=bool)

    found_simplices = tri.find_simplex(test_data)
    pass_one_accept = np.argwhere(found_simplices != -1)
    pass_two_accept = pass_one_accept[np.argwhere(alpha_boolean[found_simplices[found_simplices != -1]] == True)]
    acceptance_boolean[pass_two_accept] = 1

    return test_data[acceptance_boolean], test_data[np.invert(acceptance_boolean)], acceptance_boolean


text:
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score as AUC
from sklearn.metrics import mean_absolute_error
from sklearn.decomposition import PCA
from sklearn.preprocessing import LabelEncoder, LabelBinarizer
from sklearn.cross_validation import cross_val_score
from copy import deepcopy

train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')

print train.shape
print test.shape
print train.describe()
print pd.isnull(train).values.any()

# print train.info()
cat_features = list(train.select_dtypes(include=['object']).columns)
print("Number of categorical variables are:", len(cat_features))
print cat_features

cont_features = [cont for cont in list(train.select_dtypes(
                 include=['float64', 'int64']).columns) if cont not in ['loss', 'id']]
print "Continuous: {} features".format(len(cont_features))

id_col = list(train.select_dtypes(include=['int64']).columns)
print "A column of int64: {}".format(id_col)

cat_uniques =[]
for cat in cat_features:
    cat_uniques.append(len(train[cat].unique()))

print cat_uniques

unique_val_categories = pd.DataFrame.from_items([('cat_names',cat_features),('unique_values',cat_uniques)])
print unique_val_categories.head()

fig,(ax1,ax2)= plt.subplots(1,2)
fig.set_size_inches(16,5)
ax1.hist(unique_val_categories.unique_values,bins=50)
ax1.set_title('Amount of categorical features with X distinct values')
ax1.set_xlabel('Distinct values in a feature')
ax1.set_ylabel('Features')

plt.figure(figsize=(16,8))
plt.plot(train['id'], train['loss'])
plt.title('Loss values per id')
plt.xlabel('id')
plt.ylabel('loss')
plt.legend()

print stats.mstats.skew(train['loss']).data
print stats.mstats.skew(np.log(train['loss'])).data

fig, (ax1, ax2) = plt.subplots(1,2)
fig.set_size_inches(16,5)
ax1.hist(train['loss'], bins=50)
ax1.set_title('Train Loss target histogram')
ax1.grid(True)
ax2.hist(np.log(train['loss']), bins=50, color='g')
ax2.set_title('Train Log Loss target histogram')
ax2.grid(True)
# plt.show()

train[cont_features].hist(bins=50,figsize=(16,12))
# plt.show()

for con in cont_features:
    print stats.mstats.skew(train[con]).data

plt.subplots(figsize=(16,9))
correlation_mat = train[cont_features].corr()
sns.heatmap(correlation_mat, annot=True)
plt.show()

max_stars_repo_name: Sharmarajnish/Allstate-Claims-Severity

max_issues_repo_name: Sharmarajnish/Allstate-Claims-Severity

max_forks_repo_name: Sharmarajnish/Allstate-Claims-Severity

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 17 18:00:27 2017

@author: lenovo
"""
import pandas
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
import numpy as np

class stock_machine(object):
    def __init__(self):
        self.data = pandas.read_csv('data_spek.csv')
        self.target = pandas.read_csv('target_spek.csv')

    def fitting_frt(self):
        self.forest = RandomForestClassifier(n_estimators=5, random_state=2)
        self.data=self.data.values
        self.target = self.target.values.reshape(len(self.target),)
        self.forest.fit(self.data,self.target)
        RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gint', max_depth=None, max_features='auto', max_leaf_nodes=None,min_samples_leaf=1,min_samples_split=2,min_weight_fraction_leaf=0.0, n_estimators=5,n_jobs=1,oob_score=False, random_state=2, verbose= 2, warm_start=False)
        y_pred_frt = self.forest.predict(self.data)
        return np.mean(y_pred_frt == self.target)

    def fitting_knn(self):
        self.knn =  KNeighborsClassifier(n_neighbors=1)
        self.knn.fit(self.data,self.target)
        KNeighborsClassifier(algorithm='auto',leaf_size=30,metric='minkowski',metric_params=None,n_jobs=1,n_neighbors=1,p=2,weights='uniform')
        y_pred_knn = self.knn.predict(self.data)
        return (np.mean(y_pred_knn == self.target))

    def predict_frt(self,pps,book_prc,divy,pe):
        x_new = np.array([[pps,book_prc,divy,pe]])
        prediction_frt = self.forest.predict(x_new)
        return prediction_frt

    def predict_knn(self,pps,book_prc,divy,pe):
        x_new = np.array([[pps,book_prc,divy,pe]])
        prediction_knn = self.knn.predict(x_new)
        return prediction_knn

if __name__ == '__main__':
    s=stock_machine()
    print("%.2f %%"  %(s.fitting_frt()*100))
    print(s.predict_frt(1.2,1.1,1.1,12.0))


text:
from stuff import *
import sys
from scipy.optimize import minimize
from random import randrange,seed

# Estimate growth advantage of B.1.617.2 over B.1.1.7 by correlating, over LTLAs, the change in R
# (over two pairs of weeks) with the relative prevalence of B.1.617.2 as estimated by Sanger sequencing data.

# Get ltla.csv from https://coronavirus.data.gov.uk/api/v2/data?areaType=ltla&metric=newCasesBySpecimenDate&format=csv
# Get Sanger data from https://covid-surveillance-data.cog.sanger.ac.uk/download/lineages_by_ltla_and_week.tsv

# Get variant numbers from Sanger

apicases=loadcsv("ltla.csv")
sanger=loadcsv("lineages_by_ltla_and_week.2021-05-08.tsv",sep='\t')

# Get case counts c0, c1, c2, c3 for weeks commencing on these dates
# Then estimate old R by (c1/c0)^(mgt/gap) and new R by (c3/c2)^(mgt/gap)
#weeks=['2021-04-11','2021-04-18', '2021-04-25','2021-05-02']
weeks=['2021-04-11','2021-04-18', '2021-04-28','2021-05-05']
#weeks=['2021-04-18','2021-04-25', '2021-04-28','2021-05-05']
sangerdate='2021-05-01'
mgt=5# Mean generation time in days
exclude=set()# set(["E08000001"])

days=[datetoday(w) for w in weeks]
cases={}
for (ltla,date,n) in zip(apicases['areaCode'],apicases['date'],apicases['newCasesBySpecimenDate']):
  if ltla not in cases: cases[ltla]=[0,0,0,0]
  day=datetoday(date)
  for i in range(4):
    if day>=days[i] and day<days[i]+7: cases[ltla][i]+=n

ltlas=set()
sangnum={}
for (date,ltla,var,n) in zip(sanger['WeekEndDate'],sanger['LTLA'],sanger['Lineage'],sanger['Count']):
  if ltla in cases and ltla not in exclude and date==sangerdate:
    if ltla not in sangnum: sangnum[ltla]=[0,0]
    if var=="B.1.617.2": sangnum[ltla][0]+=n
    sangnum[ltla][1]+=n
ltlas=sorted(list(ltla for ltla in sangnum if sangnum[ltla][0]>0))

l=[]
with open('graph','w') as fp:
  for ltla in ltlas:
    (r,n)=sangnum[ltla]
    cc=cases[ltla]
    if cc[0]==0 or cc[1]==0 or cc[2]==0 or cc[3]==0: continue
    gapold=days[1]-days[0]
    gapnew=days[3]-days[2]
    Rold=(cc[1]/cc[0])**(mgt/gapold)
    Rnew=(cc[3]/cc[2])**(mgt/gapnew)
    x=r/n
    y=Rnew/Rold
    # Calculate approx variances of x and y
    p=(r+1)/(n+2);vx=p*(1-p)/n
    vy=y*y*sum(1/c for c in cc)*mgt*((1/gapold+1/gapnew)/2)# This mgt factor on the variance is a fudge
    l.append([x,y,vx,vy,n])
    print("%s %5.3f %5.3f %5.3f %5.3f"%(ltla,x,y,vx,vy),file=fp)

def err(xx,l,cutoff):
  a,b=xx
  tot=0
  for (x,y,vx,vy,n) in l:
    if n>=cutoff:
      tot+=(y-(a+b*x))**2/(vy+b*b*vx)
  return tot

for cutoff in [1,5,10,20,30,50]:
  xx=[1,1]
  bounds=[(0.5,1.5),(-1,2)]
  res=minimize(err,xx,args=(l,cutoff),method="SLSQP",bounds=bounds,options={"maxiter":1000})
  if res.success:
    (a0,b0)=res.x
  else:
    raise RuntimeError(res.message)
  l0=[x for x in l if x[4]>=cutoff]
  n=len(l0)
  # Bootstrap confidence interval
  sl=[]
  nit=1000
  for it in range(nit):
    m=[l0[randrange(n)] for i in range(n)]
    res=minimize(err,xx,args=(m,cutoff),method="SLSQP",bounds=bounds,options={"maxiter":1000})
    if res.success:
      (a,b)=res.x
      sl.append(b)
    else:
      raise RuntimeError(res.message)
  sl.sort()
  conf=0.98
  (bmin,bmax)=sl[int((1-conf)/2*nit)],sl[int((1+conf)/2*nit)]
  print("Cutoff %2d : %3d LTLAs : Extra transmissibility %3.0f%% (%.0f%% CI: %.0f - %.0f%%)"%(cutoff,n,b0*100,conf*100,bmin*100,bmax*100))

reason: from scipy


text:
#   Copyright 2018 Samuel Payne sam_payne@byu.edu
#   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 pandas as pd
import numpy as np
import os
import warnings
import datetime

from cptac.dataset import Dataset
from cptac.dataframe_tools import *
from cptac.exceptions import FailedReindexWarning, PublicationEmbargoWarning, ReindexMapError

class PdcOv(Dataset):

    def __init__(self, version="latest", no_internet=False):
        """Load all of the dataframes as values in the self._data dict variable, with names as keys, and format them properly.

        Parameters:
        version (str, optional): The version number to load, or the string "latest" to just load the latest building. Default is "latest".
        no_internet (bool, optional): Whether to skip the index update step because it requires an internet connection. This will be skipped automatically if there is no internet at all, but you may want to manually skip it if you have a spotty internet connection. Default is False.
        """

        # Set some needed variables, and pass them to the parent Dataset class __init__ function

        # This keeps a record of all versions that the code is equipped to handle. That way, if there's a new data release but they didn't update their package, it won't try to parse the new data version it isn't equipped to handle.
        valid_versions = ["1.0"]

        data_files = {
            "1.0": [
                "clinical.tsv.gz",
                "phosphoproteome.tsv.gz",
                "proteome.tsv.gz",
                "OV_sample_TMT_annotation_UMich_GENCODE34_0315.csv"
            ]
        }

        # Call the parent class __init__ function
        super().__init__(cancer_type="pdcov", version=version, valid_versions=valid_versions, data_files=data_files, no_internet=no_internet, attempt_update_index=False)

        # Load the data into dataframes in the self._data dict
        loading_msg = f"Loading {self.get_cancer_type()} v{self.version()}"
        for file_path in self._data_files_paths: # Loops through files variable

            # Print a loading message. We add a dot every time, so the user knows it's not frozen.
            loading_msg = loading_msg + "."
            print(loading_msg, end='\r')

            path_elements = file_path.split(os.sep) # Get a list of the levels of the path
            file_name = path_elements[-1] # The last element will be the name of the file. We'll use this to identify files for parsing in the if/elif statements below

            if file_name == "clinical.tsv.gz":
                df = pd.read_csv(file_path, sep="\t", index_col=0)
                clin_drop_rows = ['JHU QC', 'PNNL-JHU Ref']
                df = df.drop(clin_drop_rows, axis = 'index') # drop quality control and reference
                self._data["clinical"] = df

            if file_name == "phosphoproteome.tsv.gz":
                df = pd.read_csv(file_path, sep="\t")
                self._data["phosphoproteomics"] = df

            if file_name == "proteome.tsv.gz":
                df = pd.read_csv(file_path, sep="\t")
                self._data["proteomics"] = df

            # This file maps Ov aliquots to patient IDs (case ID with tissue type) and
            # can be found on Box under CPTAC/cptac/pancan/helper_files
            elif file_name == "OV_sample_TMT_annotation_UMich_GENCODE34_0315.csv":
                ov_map = pd.read_csv(file_path, sep = ",", usecols = ['specimen', 'sample'])
                ov_map = ov_map.loc[~ ov_map['sample'].str.contains('JHU', regex = True)] # drop quality control rows
                ov_map = ov_map.set_index('specimen')
                map_dict = ov_map.to_dict()['sample'] # create dictionary with aliquots as keys and patient IDs as values
                self._helper_tables["map_ids"] = map_dict

        print(' ' * len(loading_msg), end='\r') # Erase the loading message
        formatting_msg = f"Formatting {self.get_cancer_type()} dataframes..."
        print(formatting_msg, end='\r')

        # Get dictionary with aliquots as keys and patient IDs as values
        mapping_dict = self._helper_tables["map_ids"]

        # Proteomics
        prot = self._data["proteomics"]
        prot['Patient_ID'] = prot['aliquot_submitter_id'].replace(mapping_dict) # replace aliquot with patient IDs
        prot = prot.drop(['aliquot_submitter_id', 'case_submitter_id'], axis = 'columns')
        prot = prot.set_index('Patient_ID')
        prot.index = prot.index.str.replace('-T$','', regex = True)
        prot.index = prot.index.str.replace('-N$','.N', regex = True)
        self._data["proteomics"] = prot

        # Phosphoproteomics
        phos_df = self._data["phosphoproteomics"]
        phos_df['Patient_ID'] = phos_df['aliquot_submitter_id'].replace(mapping_dict) # replace aliquot with patient IDs
        phos_df = phos_df.set_index('Patient_ID')
        phos_df = phos_df.drop(['aliquot_submitter_id', 'case_submitter_id'], axis = 'columns')
        phos_df.index = phos_df.index.str.replace('-T$','', regex = True)
        phos_df.index = phos_df.index.str.replace('-N$','.N', regex = True)
        phos_df = map_database_to_gene_pdc(phos_df, 'refseq') # map refseq IDs to gene names
        self._data["phosphoproteomics"] = phos_df

        # Sort rows (tumor first then normal) and columns by first level (protein/gene name)
        self._data = sort_all_rows_pancan(self._data)

        print(" " * len(formatting_msg), end='\r') # Erase the formatting message

max_stars_repo_name: old-rob/cptac

max_issues_repo_name: old-rob/cptac

max_forks_repo_name: old-rob/cptac


text:
import glob
import numpy as np
import tensorflow as tf
from PIL import Image

def weight(shape):
    return tf.Variable(tf.truncated_normal(shape, stddev=0.1))

def bias(length):
    return tf.Variable(tf.constant(0.1, shape=[length]))

def layer(input, num_input_channels, filter_size, num_filters, training,
          use_bn=True, use_relu=True, use_pool=True):
    shape = [filter_size, filter_size, num_input_channels, num_filters]
    weights = weight(shape)
    biases = bias(num_filters)

    layer = tf.nn.conv2d(input=input, filter=weights, strides=[1, 1, 1, 1],
                         padding="SAME")
    layer += biases

    if use_bn:
        layer = tf.layers.batch_normalization(layer, training=training)

    if use_relu:
        layer = tf.nn.relu(layer)

    if use_pool:
        layer = tf.nn.max_pool(value=layer, ksize=[1, 2, 2, 1],
                               strides=[1, 2, 2, 1], padding="SAME")

    return layer

def flat_layer(layer):
    layer_shape = layer.get_shape()
    num_features = layer_shape[1:4].num_elements()
    layer_flat = tf.reshape(layer, [-1, num_features])
    return layer_flat, num_features

def fc_layer(input, num_inputs, num_outputs, use_relu=True):
    weights = weight([num_inputs, num_outputs])
    biases = bias(num_outputs)
    layer = tf.matmul(input, weights) + biases
    if use_relu:
        layer = tf.nn.relu(layer)
    return layer

filter_size1 = 3
num_filters1 = 32
filter_size2 = 3
num_filters2 = 64
filter_size3 = 3
num_filters3 = 128
filter_size4 = 3
num_filters4 = 256
fc_size = 128
num_channels = 3
img_size = 128
img_size_flat = img_size * img_size * num_channels
img_shape = (img_size, img_size)
classes = ['maze', 'pig']
num_classes = len(classes)

test_x = tf.placeholder(tf.float32, shape=[None, img_size_flat], name='test_x')
test_x_image = tf.reshape(test_x, [-1, img_size, img_size, num_channels])
test_y_true = tf.placeholder(tf.float32, shape=[None, num_classes],
                             name='test_y_true')
test_y_true_cls = tf.argmax(test_y_true, dimension=1)

test_layer1 = layer(input=test_x_image, num_input_channels=num_channels,
                    filter_size=filter_size1, num_filters=num_filters1,
                    training=False)
test_layer2 = layer(input=test_layer1, num_input_channels=num_filters1,
                    filter_size=filter_size2, num_filters=num_filters2,
                    training=False)
test_layer3 = layer(input=test_layer2, num_input_channels=num_filters2,
                    filter_size=filter_size3, num_filters=num_filters3,
                    training=False)
# test_layer4 = layer(input=test_layer3, num_input_channels=num_filters3,
#                     filter_size=filter_size4, num_filters=num_filters4,
#                     training=False)
test_layer_flat, test_num_features = flat_layer(test_layer3)
test_layer_fc1 = fc_layer(input=test_layer_flat, num_inputs=test_num_features,
                          num_outputs=fc_size, use_relu=True)
test_layer_fc2 = fc_layer(input=test_layer_fc1, num_inputs=fc_size,
                          num_outputs=num_classes, use_relu=False)

test_cross_entropy = tf.nn.softmax_cross_entropy_with_logits(
    logits=test_layer_fc2, labels=test_y_true)
test_y_pred_cls = tf.argmax(test_layer_fc2, dimension=1)

test_correct_prediction = tf.equal(test_y_pred_cls, test_y_true_cls)
test_loss = tf.reduce_mean(test_cross_entropy)
test_accuracy = tf.reduce_mean(tf.cast(test_correct_prediction, tf.float32))

test_path = "data/original/*.png"

test_set = []
test_files = glob.glob(test_path)
for file in test_files:
    img = Image.open(file)
    test_set.append((np.array(img)))
    img.close()

test_labels = []
file_test = "data/original/test_label.txt"
ft = open(file_test, "r")
try:
    for line in ft:
        test_labels.append(line.strip("\n"))
finally:
    ft.close()

test_label = np.zeros((len(test_labels), num_classes))
for i in range(len(test_labels)):
    test_label[i, int(test_labels[i])] = 1

sess = tf.Session()
saver = tf.train.Saver()
if tf.train.latest_checkpoint('ckpts') is not None:
    saver.restore(sess, tf.train.latest_checkpoint('ckpts'))
else:
    assert 'can not find checkpoint folder path!'

test = np.array(test_set).reshape(len(test_set), img_size_flat)
test_los, test_acc = sess.run(
    [test_loss, test_accuracy],
    feed_dict={test_x: test, test_y_true: test_label})
msg = "Test Accuracy: {0:>6.1%}, Loss: {1:.3f}"
print(msg.format(test_acc, test_los))

sess.close()

max_stars_repo_path: CNN_by_Chong Feng/test.py

max_stars_repo_name: shaw-wong/Malmo

max_issues_repo_path: CNN_by_Chong Feng/test.py

max_issues_repo_name: shaw-wong/Malmo

max_forks_repo_path: CNN_by_Chong Feng/test.py

max_forks_repo_name: shaw-wong/Malmo


text:
import matplotlib.pyplot as plt
import numpy as np

from keras.models import Sequential
from keras.layers import LSTM, Dense, Activation
from keras.optimizers import Adam

from sklearn.preprocessing import MinMaxScaler

# Set random seed for reproducibility
np.random.seed(1000)

# Download the dataset from: https://datamarket.com/data/set/22ti/zuerich-monthly-sunspot-numbers-1749-1983#!ds=22ti&display=line
dataset_filename = '<YOUR_PATH>\dataset.csv'

n_samples = 2820
data = np.zeros(shape=(n_samples, ), dtype=np.float32)

sequence_length = 15

if __name__ == '__main__':
    # Load the dataset
    with open(dataset_filename, 'r') as f:
        lines = f.readlines()

    for i, line in enumerate(lines):
        if i == 0:
            continue

        if i == n_samples + 1:
            break

        _, value = line.split(',')
        data[i - 1] = float(value)

    # Scale the dataset between -1 and 1
    mmscaler = MinMaxScaler((-1.0, 1.0))
    data = mmscaler.fit_transform(data.reshape(-1, 1))

    # Show the dataset
    fig, ax = plt.subplots(figsize=(20, 10))

    ax.plot(data)
    ax.grid()
    ax.set_xlabel('Time steps')
    ax.set_ylabel('Monthly sunspots numbers')
    plt.show()

    # Create the train and test sets (rounding to 2800 samples)
    X_ts = np.zeros(shape=(n_samples - sequence_length, sequence_length, 1), dtype=np.float32)
    Y_ts = np.zeros(shape=(n_samples - sequence_length, 1), dtype=np.float32)

    for i in range(0, data.shape[0] - sequence_length):
        X_ts[i] = data[i:i + sequence_length]
        Y_ts[i] = data[i + sequence_length]

    X_ts_train = X_ts[0:2300, :]
    Y_ts_train = Y_ts[0:2300]

    X_ts_test = X_ts[2300:2800, :]
    Y_ts_test = Y_ts[2300:2800]

    # Create the model
    model = Sequential()

    model.add(LSTM(4, stateful=True, batch_input_shape=(20, sequence_length, 1)))

    model.add(Dense(1))
    model.add(Activation('tanh'))

    # Compile the model
    model.compile(optimizer=Adam(lr=0.001, decay=0.0001),
                  loss='mse',
                  metrics=['mse'])

    # Train the model
    model.fit(X_ts_train, Y_ts_train,
              batch_size=20,
              epochs=100,
              shuffle=False,
              validation_data=(X_ts_test, Y_ts_test))

    # Show the result
    fig, ax = plt.subplots(figsize=(20, 10))

    ax.plot(Y_ts_test, label='True values')
    ax.plot(model.predict(X_ts_test, batch_size=20), label='Predicted values')
    ax.grid()
    ax.set_xlabel('Time steps')
    ax.set_ylabel('Monthly sunspots numbers')
    ax.legend()
    plt.show()


text:
import math
import numpy as np
import pygame
import config
from painter import Painter

class PainterExpand(Painter):
    def __init__(self, holder):
        Painter.__init__(self, holder)
        self.clear((0, 0, 0))
        self._led_size = 0
        self._small_surface = None
        self._small_surface_w = None
        self._small_surface_h = None

    def paint(self, wave_data, x, y, output_scroll, output_ennergy, output_spectrum, fft_result):
        Painter.paint(self, wave_data, x, y, output_scroll, output_ennergy, output_spectrum, fft_result)

        w, h = pygame.display.get_surface().get_size()

        #led_data = output_scroll
        led_data = output_ennergy
        #led_data = output_spectrum

        if led_data is not None:
            if self._small_surface is None:
                self._small_surface_h = len(led_data[0])
                self._small_surface_w = int(self._small_surface_h / h * w)
                self._small_surface = pygame.Surface((self._small_surface_w, self._small_surface_h))
                self._small_surface.fill((0, 0, 0))

            alpha_surface = pygame.Surface((self._small_surface_w, self._small_surface_h))
            alpha = 2
            alpha_surface.set_alpha(alpha)
            alpha_surface.fill((0, 0, 0))
            self._small_surface.blit(alpha_surface, (0, 0))

            self.draw_led_points(led_data)

            x_posi = self._small_surface.get_width() // 2
            self._small_surface.blit(self._small_surface, (x_posi + 1, 0), (x_posi, 0, self._small_surface_w - x_posi, self._small_surface_h))
            self._small_surface.blit(self._small_surface, (0, 0), (1, 0, x_posi, self._small_surface_h))

            #large_image = pygame.transform.smoothscale(self._small_surface, (w, h))
            large_image = pygame.transform.scale(self._small_surface, (w, h))

            self._holder.surface.blit(large_image, (0, 0))

            """
            # 这是原来的方法，速度较慢
            self._led_size = h // len(led_data[0])
            x_posi = (w - self._led_size) // 2
            y_posi = (h - self._led_size * len(led_data[0])) // 2

            self._holder.surface.blit(self._holder.surface, (x_posi + self._led_size, 0), (x_posi, 0, w - x_posi, h))
            x = x_posi % self._led_size
            self._holder.surface.blit(self._holder.surface, (x - self._led_size, 0), (x, 0, x_posi - x + self._led_size, h))

            self.clear_with_alpha((0, 0, 0), 2)
            #self.blur(2)

            self.draw_led(led_data, x_posi, y_posi)
            """

    def draw_led(self, led_colors, x_posi, y_posi):
        if led_colors is not None:
            for i in range(len(led_colors[0])):
                pixels = np.clip(led_colors, 0, 255).astype(int)
                pygame.draw.rect(self._holder.surface, (pixels[0][i], pixels[1][i], pixels[2][i]),
                                 (x_posi, y_posi + i*self._led_size, self._led_size, self._led_size))

    def draw_led_points(self, led_colors):
        if led_colors is not None:
            x_posi = self._small_surface.get_width() // 2
            for i in range(len(led_colors[0])):
                pixels = np.clip(led_colors, 0, 255).astype(int)
                self._small_surface.set_at((x_posi, i), (pixels[0][i], pixels[1][i], pixels[2][i]))


text:
import numpy as np
from scipy.sparse import coo_matrix
from shape_function import *
from constitutive import *
from quadrature import *
# import time
from joblib import Parallel, delayed

def get_element_stiffness(ele, coord, connect, E, nu, el_type, problem_type, ngp):
    """Calculates the elemental stiffness matrix
        Inputs:
        1) ele = element number (starts from 0)
        2) coord = coordinate array
        3) connect = connectivity array
        4) E = Young's Modulus
        5) nu = Poisson's Ratio
        6) el_type = type of element used for meshing
            el_type = 0 (4 noded quadrilateral)
            el_type = 1 (8 noded quadrilateral)
        7) problem_type = type of 2D elasticity problem
            problem_type = 0 (plane stress)
            problem_type = 1 (plane strain)
        8) ngp = number of Gauss points in each direction used for integration

        Output:
        K_ele = element stiffness matrix"""

    # Constitutive matrix
    C = constitutive(E, nu, problem_type)
    # Nodes
    node = connect[ele, :]
    # Coordinates of the nodes
    vertex_coord = coord[node, :]
    # Element stiffness matrix
    K_ele = np.zeros((2 * vertex_coord.shape[0], 2 * vertex_coord.shape[0]))
    # Gauss points and weights
    gp, weights = quadrature(ngp)

    for i in range(ngp):
        for j in range(ngp):
            xi = gp[i]
            eta = gp[j]
            shape_func = ShapeFunction(xi, eta, el_type)
            shape_func.compute_Jacobian(vertex_coord)
            J = shape_func.J
            B = shape_func.get_B_matrix()
            K_ele += (B.T @ C @ B) * weights[i] * weights[j] * np.linalg.det(J)

    return K_ele

def global_stiffness(coord, connect, E, nu, el_type, problem_type, ngp):
    """Function to calculate the global stiffness matrix

    Inputs:
    1) coord = coordinate array
    2) connect = connectivity array
    3) E = Young's Modulus
    4) nu = Poisson's ratio
    5) el_type = type of element used for meshing
        el_type = 0 (4 noded quadrilateral)
        el_type = 1 (8 noded quadrilateral)
    6) problem_type = type of 2D elasticity problem
        problem_type = 0 (plane stress)
        problem_type = 1 (plane strain)
    7) ngp = number of Gauss points in each direction used for integration

    Output:
    K_global = global stiffness matrix"""

    n_ele = connect.shape[0]    # No. of elements
    node_perele = connect.shape[1]  # No. of nodes per element
    ndof_perele = 2*node_perele # No. of DOF per element

    # tic = time.perf_counter()
    K_par = Parallel(n_jobs=-1)(delayed(get_element_stiffness)(ele, coord, connect, E, nu, el_type, problem_type, ngp)
                               for ele in range(n_ele))

    # K_temp = np.zeros((n_ele, ndof_perele, ndof_perele))
    # for i in range(K_temp.shape[0]):
    #     K_temp[i, :, :] = K_par[i]

    K_temp = np.array(K_par)

    I = np.zeros((n_ele, K_temp.shape[1] ** 2))
    J = np.zeros((n_ele, K_temp.shape[1] ** 2))

    dof = np.zeros((n_ele, 2 * node_perele))
    for j in range(connect.shape[1]):
        dof[:, 2 * j] = 2 * connect[:, j]
        dof[:, 2 * j + 1] = 2 * connect[:, j] + 1

    for i in range(dof.shape[1]):
        I[:, i*ndof_perele:(i+1)*ndof_perele] = dof[:, i].reshape(dof.shape[0], 1)*np.ones((dof.shape[0], ndof_perele))
        J[:, i*ndof_perele:(i+1)*ndof_perele] = dof

    I = I.flatten()
    J = J.flatten()
    V = K_temp.flatten()
    # Forming the global stiffness matrix in COO format
    K_global = coo_matrix((V, (I, J)), shape=(2 * coord.shape[0], 2 * coord.shape[0]))
    # Converting the global stiffness matrix to LIL format to perform indexing and slicing operations
    K_global = K_global.tocsr()

    return K_global


text:
#%%
import os
import glob
import numpy as np
import scipy as sp
import pandas as pd
import git

# Import library to perform maximum entropy fits
from maxentropy.skmaxent import FeatureTransformer, MinDivergenceModel

# Import matplotlib stuff for plotting
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib as mpl

# Seaborn, useful for graphics
import seaborn as sns

# Import the project utils
import ccutils

# Set PBoC plotting format
ccutils.viz.set_plotting_style()

#%%

# Find home directory for repo
repo = git.Repo("./", search_parent_directories=True)
homedir = repo.working_dir

# Define directories for data and figure
figdir = f'{homedir}/fig/si/'

# %%
# Fit a model p(x) for dice probabilities (x=1,...,6) with the
# single constraint E(X) = 4.5
def first_moment_die(x):
    return np.array(x)

# Put the constraint functions into an array
features = [first_moment_die]
# Write down the constraints (in this case mean of 4.5)
k = np.array([4.5])

# Define the sample space of the die (from 1 to 6)
samplespace = list(range(1, 7))

# Define the minimum entropy
model = MinDivergenceModel(features, samplespace)

# Change the dimensionality of the array
X = np.atleast_2d(k)

# Fit the model
model.fit(X)

# initialize figure
fig, ax = plt.subplots(1, 2, figsize=(6, 2.5), sharex=True, sharey=True)

# Define probability distribution of the "wrong inference"
prob = [0, 0, 0, 0.5, 0.5, 0]
# Plot the "wrong" distribution
ax[0].bar(samplespace, prob)

# Plot the max ent distribution
ax[1].bar(samplespace, model.probdist())

# Label axis
ax[0].set_xlabel("die face")
ax[1].set_xlabel("die face")

ax[0].set_ylabel("probability")

# Set title for plots
ax[0].set_title(r"$\left\langle x \right\rangle = 4.5$")
ax[1].set_title(r"MaxEnt $\left\langle x \right\rangle = 4.5$")

# Add letter label to subplots
plt.figtext(0.1, 0.93, "(A)", fontsize=8)
plt.figtext(0.50, 0.93, "(B)", fontsize=8)

plt.subplots_adjust(wspace=0.05)

plt.savefig(figdir + "figS15.pdf", bbox_inches="tight")

reason: import numpy,import scipy


text:
import requests
import numpy as np

from flask import Blueprint, request, make_response, jsonify
from flask import current_app as app

popular_skills_bp = Blueprint('popular_skills_bp', __name__)

@popular_skills_bp.route('/popular-skills/', methods = ['GET'])
def popular_skills():
    if request.method != 'GET' :
        return make_response('Wrong query performed.', 400)

    paramsReq = {'offset':'0', 'size':'0', 'aggregate':'true'}
    resp = requests.post('https://search.torre.co/opportunities/_search', params = paramsReq)

    if resp.status_code == 400:
        return make_response('Information API of job search is not responding.', 400)

    n = int(request.args.get('size'))
    skills = resp.json()['aggregators']['skill']
    amount = np.array([])
    name = np.array([])

    for skill in skills:

        if len(amount) < n:
            amount = np.append(amount, skill['total'])
            name = np.append(name, skill['value'])

        if (min(amount) < skill['total']) and (len(amount) == n):
            minIndexRes = np.where(amount == min(amount))
            minIndex = minIndexRes[0]
            amount[minIndex[0]] = skill['total']
            name[minIndex[0]] = skill['value']

    sortIndexes = np.argsort(-amount, kind = 'mergesort')

    amountSent = np.array([])
    nameSent = np.array([])
    percentageSent = np.array([])
    totalJobOpportunities = resp.json()['total']

    for index in sortIndexes:
        amountSent = np.append(amountSent, amount[index])
        nameSent = np.append(nameSent, name[index])
        percentageIndex = (amount[index]/totalJobOpportunities)*100
        percentageSent = np.append(percentageSent, round(percentageIndex, 2))

    pop_skills = {'total': amountSent.tolist(), 'value': nameSent.tolist(), 'percentage_jobs': percentageSent.tolist()}
    headers = {'Content-Type': 'applicaton/json', 'Access-Control-Allow-Origin':'https://torre-application-front.herokuapp.com'}

    return make_response(jsonify(pop_skills), 200, headers)


text:
from __future__ import annotations
from os import path
from typing import Sequence, Tuple, Any
from munch import Munch, unmunchify
import inspect
import yaml
import matplotlib.pyplot as plt
import numpy as np
import pathlib
from PIL import Image, ImageDraw
import cv2
from advpipe.log import logger
import eagerpy as ep
from torchvision import transforms

from typing import TYPE_CHECKING
if TYPE_CHECKING:
    import torch
    from typing import Callable, Optional, Iterator, Any
    from advpipe.types import TensorTypeVar

def batches(seq: Iterator[Any], bs: int = 16) -> Iterator[Sequence]:
    batch = []
    while True:
        try:
            batch.append(next(seq))
        except StopIteration:
            break

        if len(batch) == bs:
            yield batch
            batch = []

    if len(batch) > 0:
        yield batch

def tensor_batches(tensor_iter: Iterator[TensorTypeVar], bs: int = 16) -> Iterator[TensorTypeVar]:
    while True:
        batch = []
        try:
            batch.append(ep.astensor(next(tensor_iter)))
        except StopIteration:
            break

        if len(batch) == bs:
            yield ep.concatenate(batch).raw
            batch = []

    if len(batch) > 0:
        yield ep.concatenate(batch).raw

def zfill_align_float(param: float, length: int = 6) -> str:
    return f"{param:.2f}".zfill(length)

def scale_img(img: Image, target_size: int) -> Image:
    """Re-scale image while preserving its aspect ratio"""
    x, y = img.size

    # the shorter side will be scaled to the target_size
    if x < y:
        scale_ratio = target_size / x
    else:
        scale_ratio = target_size / y

    return img.resize((int(x * scale_ratio), int(y * scale_ratio)))

def quantize_img(img: TensorTypeVar, round: bool = True) -> TensorTypeVar:
    """Quantizes tensor image to uint8"""

    ep_img, restore_func = ep.astensor_(img)
    assert ep_img.min() >= 0 and ep_img.max() <= 255
    if ep_img.max() <= 1:
        ep_img *= 255

    # we have to convert it first to numpy, because eagerpy doesn't have round function, and also eagerpy's .astype(dtype) has
    # tensortype-dependent dtype argument
    if round:
        np_img = np.asarray(np.round(ep_img.numpy()), dtype=np.uint8)
    else:
        np_img = np.asarray(ep_img.numpy(), dtype=np.uint8)

    return restore_func(ep.astensor(np_img))    # type: ignore

def mkdir_p(dir_path: str) -> None:
    pathlib.Path(dir_path).mkdir(parents=True, exist_ok=True)

def show_img(img: TensorTypeVar, method: str = "pyplot") -> None:
    np_img = ep.astensor(img).numpy()
    if method == "PIL":
        convert_to_pillow(np_img).show()
    elif method == "opencv":
        cv2.imshow("", np_img[..., ::-1])
        cv2.waitKey(0)
    elif method == "pyplot":
        plt.imshow(np_img, interpolation="bicubic")
        plt.show()
    else:
        raise ValueError("utils.show_img: Unsupported method")

def is_img_filename(img_fn: str) -> bool:
    extensions = [".png", ".jpg", ".jpeg"]
    img_fn = img_fn.lower()
    return any([img_fn.endswith(ext) for ext in extensions])

def load_image_to_numpy(img_path: str) -> np.ndarray:
    image = Image.open(img_path).convert("RGB")
    # convert image to numpy array
    np_img = np.asarray(image, dtype=np.float32) / 255
    return np_img

_to_tensor = transforms.ToTensor()

def load_image_to_torch(img_path: str) -> torch.Tensor:
    image = Image.open(img_path)
    return _to_tensor(image)    # type: ignore

def write_text_to_img(img: Image, text: str, max_lines: int = 20) -> Image:
    """Writes text to the top of the image"""

    font_size = 10
    text = "\n".join(text.split("\n")[:max_lines])
    n_lines = len(text.split("\n"))
    margin = int(n_lines * font_size * 1.5)

    width, height = img.size
    new_img = Image.new("RGB", size=(width + 200, max(height, margin)))
    new_img.paste(img, (0, 0))

    draw = ImageDraw.Draw(new_img)
    draw.text((width + 10, 0), text, (255, 255, 255))

    return new_img

def convert_to_pillow(np_img: np.ndarray) -> Image:
    if np_img.dtype != np.uint8:
        if np_img.min() >= 0 and np_img.max() <= 1:
            np_img = np.asarray(np_img * 255, dtype=np.uint8)
        elif np_img.min() >= 0 and np_img.max() <= 255:
            np_img = np.asarray(np_img, dtype=np.uint8)
        else:
            raise ValueError("Cannot convert numpy array to PIL image: unsupported image format")

    return Image.fromarray(np_img)

# deprecated
def clip_linf(orig_img: np.ndarray, pertubed_img: np.ndarray, epsilon: float = 0.05) -> np.ndarray:
    min_boundary = np.clip(orig_img - epsilon * np.ones_like(orig_img), 0, 1)
    max_boundary = np.clip(orig_img + epsilon * np.ones_like(orig_img), 0, 1)
    return np.clip(pertubed_img, min_boundary, max_boundary)

def load_yaml(yaml_filename: str) -> Munch:
    with open(yaml_filename, 'r') as stream:
        return Munch.fromDict(yaml.safe_load(stream))

def load_yaml_from_str(yaml_str: str) -> Munch:
    return Munch.fromDict(yaml.safe_load(yaml_str))

# TODO: this is quite a bad fuction name
def rel_to_abs_path(relative_path: str) -> str:
    """convert caller's relative path to absolute path"""
    callers_path = inspect.stack()[1].filename
    return path.normpath(path.join(path.dirname(path.abspath(callers_path)), relative_path))

def convert_to_absolute_path(module_relative_path: str) -> str:
    """Convert path relative to advpipe's module root directory to its absolute variant"""
    return path.join(get_abs_module_path(), module_relative_path)

def get_abs_module_path() -> str:
    return rel_to_abs_path(".")

class MaxFunctionCallsExceededException(Exception):
    pass

class LossCallCounter:
    def __init__(self, loss_fn: Callable[[np.ndarray], float], max_calls: int):
        self.loss_fn: Callable[[np.ndarray], float] = loss_fn

        self.last_loss_val: float = np.inf
        self.last_img: Optional[np.ndarray] = None

        self.max_calls: int = max_calls    # test comment
        self.i = 0

    def __call__(self, pertubed_image: np.ndarray) -> float:
        if self.i >= self.max_calls:
            msg = f"Max number of function calls exceeded (max_calls={self.max_calls})"
            logger.info(f"LossCallCounter: {msg}")
            raise MaxFunctionCallsExceededException(msg)

        self.i += 1
        self.last_loss_val = self.loss_fn(pertubed_image)
        self.last_img = pertubed_image
        return self.last_loss_val

def config_attr_exists(conf: Munch, attr_name: str) -> bool:
    try:
        _ = conf.__getattr__(attr_name)
        return True
    except AttributeError:
        return False

def get_config_attr(conf: Munch, attr_name: str, default_val: Any) -> Any:
    """Returns config attribute value if it exists, otherwise returns default value"""
    val = default_val
    try:
        val = conf.__getattr__(attr_name)
    except AttributeError:
        pass
    return val

def serialize_config(conf: Munch) -> str:
    unmunched = unmunchify(conf)
    return yaml.dump(unmunched)    # type: ignore


text:
import os
import time
import random
import numpy as np
import torch
import metaworld

from collections import OrderedDict

import maml_trpo.utils.pytorch_utils as ptu

from maml_trpo.envs import AsyncVectorEnv
from maml_trpo.policies import GaussianPolicy, CategoricalPolicy
from maml_trpo.critics import ValueFun
from maml_trpo.samplers import Sampler, MetaWorldTaskSampler
from maml_trpo.utils import Logger
from maml_trpo.meta_learner import MetaLearner

def init_metaworld_env(benchmark, seed, num_workers):
    """
    initializes metaworld env across num_workers

    Args:
    benchmark : Metworld Benchmark(one of ML1, ML10, ML45 from metaworld)
    seed : random seed
    num_workers : number of parallel workers

    Returns
    MetaWorld Env distributed across workers
    """
    _, env_cls = list(benchmark.train_classes.items())[0]
    seeds = [ seed * (i+1) for i in range(num_workers)]
    env = AsyncVectorEnv([env_cls for _ in range(num_workers)], seeds)
    return env

def set_random_seed(seed):
    """
    Sets random seed for reproducibility of experiments
    """
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)

def set_benchmark(benchmark_type):
    """
    initializes metaworld benchmark
    Args:
    benchmark_type : str input one of ('ml1', 'ml10', 'ml45')

    Returns
    MetaWorld Benchmark
    """
    if benchmark_type == 'ml1':
        # pick-place-v2 is chosen arbitarily, can choose any other env
        return metaworld.ML1('pick-place-v2')
    elif benchmark_type == 'ml10':
        return metaworld.ML10()
    elif benchmark_type == 'ml45':
        return metaworld.ML45()

class Trainer:
    """
    Trainer class to run train and evaluation experiments
    """

    def __init__(self, params):

        self.params = params
        self.benchmark = set_benchmark(self.params['benchmark_type'])

        set_random_seed(self.params['seed'])

        self.env = init_metaworld_env(self.benchmark, self.params['seed'], self.params['num_workers'])
        self.eval_env = init_metaworld_env(self.benchmark, self.params['seed'], 1)

        self.task_sampler = MetaWorldTaskSampler(self.benchmark, mode='train')
        self.test_task_sampler = MetaWorldTaskSampler(self.benchmark, mode='test')

        self.sampler = Sampler(self.env, num_workers=self.params['num_workers'])
        self.render_sampler = Sampler(self.eval_env, num_workers=1)

        obs_space, ac_space = self.env.get_spaces()

        self.device = ptu.init_gpu(self.params['use_gpu'], self.params['which_gpu'])

        self.policy = GaussianPolicy(
            in_size = int(np.prod(obs_space.shape)),
            out_size = int(np.prod(ac_space.shape)),
            hidden_units=[100, 100],
            activation='tanh',
            device=self.device)

        self.value_fn = ValueFun(
            in_size = int(np.prod(obs_space.shape)),
            out_size=1,
            hidden_units=[32, 32],
            activation='tanh',
            device=self.device)

        self.logger = Logger(self.params['logdir'])

        self.meta_train_itrs = self.params['n_iters']
        self.meta_train_batch_size = self.params['meta_batch_size']

        # meta-train and meta-test envs have same number of task variations in Metaworld benchmark
        # Eval-itrs is to evaluate for each task variation across meta-train and meta-test envs
        self.eval_itrs =  int(len(self.benchmark.test_tasks)/len(self.benchmark.test_classes))

        self.metalearner = MetaLearner(
            policy = self.policy,
            value_fn = self.value_fn,
            inner_lr = self.params['inner_lr'],
            num_adapt_steps = self.params['num_adapt_steps'],
            episodes_per_task = self.params['episodes_per_task'],
            gamma = self.params['gamma'],
            gae_lam = self.params['gae_lam'],
            max_kl = self.params['max_kl_increment'],
            max_backtracks = self.params['max_backtracks'],
            cg_steps = self.params['cg_steps'],
            cg_damping =self.params['cg_damping'],
            device = self.device,
            restore = self.params['restore'],
            ckpt_path = self.params['ckpt_path'])

        self.start_itr = self.metalearner.itr

        self.meta_eval_freq =  params['eval_freq']

    def meta_train(self):

        if self.start_itr!=0:
            self.start_itr+= 1

        for itr in range(self.start_itr, self.meta_train_itrs):

            print(f"At Iteration {itr}")

            logs = OrderedDict()

            tasks = self.task_sampler.sample_tasks(self.meta_train_batch_size)
            train_logs = self.metalearner.train_step(tasks, self.sampler)
            logs.update(train_logs)

            if itr==0 or ((itr+1)% self.meta_eval_freq)==0:
                meta_batch_size = len(self.benchmark.test_classes)
                test_tasks = self.test_task_sampler.sample_tasks(meta_batch_size)
                test_logs = self.metalearner.evaluate_step(test_tasks, self.sampler)
                logs.update(test_logs)

                self.metalearner.save('{}/maml_ckpt_itr_{}.pt'.format(self.params['logdir'], itr))

            for key, value in logs.items():
                print('{} : {}'.format(key, value))
                self.logger.log_scalar(value, key, itr)

            self.logger.flush()

    def meta_evaluate(self, eval_on='train_tasks', render=False):

        if eval_on == 'train_tasks':
            task_sampler = self.task_sampler
            meta_batch_size = len(self.benchmark.train_classes)
            prefix ='mt_tr'

        else:
            task_sampler = self.test_task_sampler
            meta_batch_size = len(self.benchmark.test_classes)
            prefix ='mt_ts'

        if render==True:
            sampler = self.render_sampler
        else:
            sampler = self.sampler

        for itr in range(self.eval_itrs):
            tasks = task_sampler.sample_tasks(meta_batch_size, shuffle=False)
            logs = self.metalearner.evaluate_step(tasks, sampler, log_videofn=self.logger.log_videos,
                                                                render=render, prefix=prefix,video_itr=itr)
            for key, value in logs.items():
                print('{} : {}'.format(key, value))
                self.logger.log_scalar(value, key, itr)

            self.logger.flush()

def main():

    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument('--train', action='store_true')
    parser.add_argument('--evaluate', action='store_true')
    parser.add_argument('--eval_on', type=str, choices=['train_tasks', 'test_tasks'], default='meta_train',
                                    help='test meta-trained policy on either meta_train/meta_test tasks')
    parser.add_argument('--render', action='store_true')

    parser.add_argument('--benchmark_type', type=str, default='ml10', choices=['ml1', 'ml10', 'ml45'])
    parser.add_argument('--exp_name', type=str, default='todo')
    parser.add_argument('--n_iters', '-n', type=int, default=300)
    parser.add_argument('--num_workers', '-n_w', type=int, default=1)

    parser.add_argument('--use_gae', action='store_true')
    parser.add_argument('--gae_lam', type=float, default=0.95)
    parser.add_argument('--gamma', type=float, default=0.99)

    parser.add_argument('--meta_batch_size', '-mb', type=int, default=20) # tasks collected per meta train iteration
    parser.add_argument('--episodes_per_task', type=int, default=10) # rollouts collected per a single task

    parser.add_argument('--num_adapt_steps', type =int, default=1)

    parser.add_argument('--inner_lr', type=float, default=1e-4, help='learning_rate for inner-loop/adaptation steps')

    # TRPO update params
    parser.add_argument('--cg_steps', type=int, default=10, help='num cg steps')
    parser.add_argument('--cg_damping', type=float, default=1e-1, help='damping to use in cg')
    parser.add_argument('--max_backtracks', type=int, default=10, help='no of backtracks for line search')
    parser.add_argument('--max_kl_increment', type=float, default=1e-2,
                                           help='max_kl divergence between old and new policy')

    parser.add_argument('--seed', type=int, default=1)
    parser.add_argument('--use_gpu', '-ngpu', action='store_true')
    parser.add_argument('--which_gpu', '-gpu_id', default=0)
    parser.add_argument('--scalar_log_freq', type=int, default=1)
    parser.add_argument('--eval_freq', type=int, default=10)

    parser.add_argument('--restore', action='store_true')
    parser.add_argument('--ckpt_path', type=str, help='checkpoint path to restore from')

    args = parser.parse_args()

    # convert to dictionary
    params = vars(args)

    ##################################
    ### CREATE DIRECTORY FOR LOGGING
    ##################################

    data_path = os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../data')

    if not (os.path.exists(data_path)):
        os.makedirs(data_path)

    logdir = args.benchmark_type + '_' + args.exp_name + '_' + time.strftime("%d-%m-%Y-%H-%M-%S")
    logdir = os.path.join(data_path, logdir)
    params['logdir'] = logdir
    if not(os.path.exists(logdir)):
        os.makedirs(logdir)

    ###################
    ### RUN TRAINING
    ###################

    trainer = Trainer(params)

    if args.train:
        trainer.meta_train()

    elif args.evaluate:
        trainer.meta_evaluate(params['eval_on'], params['render'])

if __name__ == "__main__":
    main()

max_stars_repo_name: adityabingi/maml-trpo-metaworld

max_issues_repo_name: adityabingi/maml-trpo-metaworld

max_forks_repo_name: adityabingi/maml-trpo-metaworld


text:
#
# Copyright 2018, 2020 Lars Pastewka
#           2019-2020 Antoine Sanner
#           2016 km7219@lsdf-28-131.scc.kit.edu
#           2016 Till Junge
#
# ### MIT license
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
"""
Standard program for a repulsive contact sim. The idea is to use this
to profile and accelerate PyCo
"""

import numpy as np

from ContactMechanics import FreeFFTElasticHalfSpace
from SurfaceTopography.Special import make_sphere
from ContactMechanics.Factory import make_system
import sys

#import matplotlib.pyplot as plt

# geometry
res = int(sys.argv[1])
print(res)
size = 65e-6
radius = 500e-6

# material
e_combo = 4e6

# system

surface = make_sphere(radius, (res, res), (size, size))
system = make_system(substrate, interaction, surface)

# minimize
offset_max = 1e-6
disp = None
for offset in np.linspace(0, offset_max, 10):
    print("offset = {}".format(offset))
    result = system.minimize_proxy(offset, initial_displacements=disp)
    disp = system.disp

print(result.success)

max_stars_repo_name: ComputationalMechanics/ContactMechanics

max_issues_repo_name: ComputationalMechanics/ContactMechanics

max_forks_repo_name: ComputationalMechanics/ContactMechanics


text:
from collections import OrderedDict

import numpy as np
import torch
from torch import nn as nn
import torch.nn.functional as F
import torch.optim as optim
import itertools

import rlkit.torch.utils.pytorch_util as ptu
from rlkit.core.trainer import Trainer
from rlkit.core.eval_util import create_stats_ordered_dict

class SoftActorCritic(Trainer):
    """
    version that:
        - uses reparameterization trick
        - has two Q functions and a V function
    TODO: Recently in rlkit there is a version which only uses two Q functions
    as well as an implementation of entropy tuning but I have not implemented
    those
    """

    def __init__(
        self,
        policy,
        qf1,
        qf2,
        vf,
        reward_scale=1.0,
        discount=0.99,
        alpha=1.0,
        policy_lr=1e-3,
        qf_lr=1e-3,
        vf_lr=1e-3,
        soft_target_tau=1e-2,
        policy_mean_reg_weight=1e-3,
        policy_std_reg_weight=1e-3,
        optimizer_class=optim.Adam,
        beta_1=0.9,
        **kwargs,
    ):
        self.policy = policy
        self.qf1 = qf1
        self.qf2 = qf2
        self.vf = vf
        self.reward_scale = reward_scale
        self.discount = discount
        self.soft_target_tau = soft_target_tau
        self.policy_mean_reg_weight = policy_mean_reg_weight
        self.policy_std_reg_weight = policy_std_reg_weight

        self.target_vf = vf.copy()
        self.eval_statistics = None

        self.policy_optimizer = optimizer_class(
            self.policy.parameters(), lr=policy_lr, betas=(beta_1, 0.999)
        )
        self.qf1_optimizer = optimizer_class(
            self.qf1.parameters(), lr=qf_lr, betas=(beta_1, 0.999)
        )
        self.qf2_optimizer = optimizer_class(
            self.qf2.parameters(), lr=qf_lr, betas=(beta_1, 0.999)
        )
        self.vf_optimizer = optimizer_class(
            self.vf.parameters(), lr=vf_lr, betas=(beta_1, 0.999)
        )

        self.alpha = alpha

    def train_step(self, batch):
        # q_params = itertools.chain(self.qf1.parameters(), self.qf2.parameters())
        # v_params = itertools.chain(self.vf.parameters())
        # policy_params = itertools.chain(self.policy.parameters())

        rewards = self.reward_scale * batch["rewards"]
        terminals = batch["terminals"]
        obs = batch["observations"]
        actions = batch["actions"]
        next_obs = batch["next_observations"]

        """
        QF Loss
        """
        # Only unfreeze parameter of Q
        # for p in itertools.chain(self.qf1.parameters(), self.qf2.parameters()):
        #     p.requires_grad = True
        # for p in self.vf.parameters():
        #     p.requires_grad = False
        # for p in self.policy.parameters():
        #     p.requires_grad = False
        self.qf1_optimizer.zero_grad()
        self.qf2_optimizer.zero_grad()
        q1_pred = self.qf1(obs, actions)
        q2_pred = self.qf2(obs, actions)
        target_v_values = self.target_vf(
            next_obs
        )  # do not need grad || it's the shared part of two calculation
        q_target = (
            rewards + (1.0 - terminals) * self.discount * target_v_values
        )  ## original implementation has detach
        q_target = q_target.detach()
        qf1_loss = 0.5 * torch.mean((q1_pred - q_target) ** 2)
        qf2_loss = 0.5 * torch.mean((q2_pred - q_target) ** 2)

        # freeze parameter of Q
        # for p in itertools.chain(self.qf1.parameters(), self.qf2.parameters()):
        #     p.requires_grad = False

        """
        VF Loss
        """
        # Only unfreeze parameter of V
        # for p in itertools.chain(self.qf1.parameters(), self.qf2.parameters()):
        #     p.requires_grad = False
        # for p in self.vf.parameters():
        #     p.requires_grad = True
        # for p in self.policy.parameters():
        #     p.requires_grad = True  ##
        self.vf_optimizer.zero_grad()
        v_pred = self.vf(obs)
        # Make sure policy accounts for squashing functions like tanh correctly!
        policy_outputs = self.policy(obs, return_log_prob=True)
        # in this part, we only need new_actions and log_pi with no grad
        new_actions, policy_mean, policy_log_std, log_pi = policy_outputs[:4]
        q1_new_acts = self.qf1(obs, new_actions)
        q2_new_acts = self.qf2(obs, new_actions)  ## error
        q_new_actions = torch.min(q1_new_acts, q2_new_acts)
        v_target = q_new_actions - self.alpha * log_pi
        v_target = v_target.detach()
        vf_loss = 0.5 * torch.mean((v_pred - v_target) ** 2)

        qf1_loss.backward()
        qf2_loss.backward()
        vf_loss.backward()

        self.qf1_optimizer.step()
        self.qf2_optimizer.step()

        self.vf_optimizer.step()

        """
        Policy Loss
        """
        # for p in itertools.chain(self.qf1.parameters(), self.qf2.parameters()):
        #     p.requires_grad = False
        # for p in self.vf.parameters():
        #     p.requires_grad = False
        # for p in self.policy.parameters():
        #     p.requires_grad = True
        new_actions, policy_mean, policy_log_std, log_pi = policy_outputs[:4]
        q1_new_acts = self.qf1(obs, new_actions)
        q2_new_acts = self.qf2(obs, new_actions)  ## error
        q_new_actions = torch.min(q1_new_acts, q2_new_acts)

        self.policy_optimizer.zero_grad()
        policy_loss = torch.mean(self.alpha * log_pi - q_new_actions)  ##
        mean_reg_loss = self.policy_mean_reg_weight * (policy_mean ** 2).mean()
        std_reg_loss = self.policy_std_reg_weight * (policy_log_std ** 2).mean()
        policy_reg_loss = mean_reg_loss + std_reg_loss
        policy_loss = policy_loss + policy_reg_loss
        policy_loss.backward()
        self.policy_optimizer.step()

        """
        Update networks
        """
        # unfreeze all -> initial states
        # for p in itertools.chain(self.qf1.parameters(), self.qf2.parameters()):
        #     p.requires_grad = True
        # for p in self.vf.parameters():
        #     p.requires_grad = True
        # for p in self.policy.parameters():
        #     p.requires_grad = True

        # unfreeze parameter of Q
        # for p in itertools.chain(self.qf1.parameters(), self.qf2.parameters()):
        #     p.requires_grad = True

        self._update_target_network()

        """
        Save some statistics for eval
        """
        if self.eval_statistics is None:
            """
            Eval should set this to None.
            This way, these statistics are only computed for one batch.
            """
            self.eval_statistics = OrderedDict()
            self.eval_statistics["Reward Scale"] = self.reward_scale
            self.eval_statistics["QF1 Loss"] = np.mean(ptu.get_numpy(qf1_loss))
            self.eval_statistics["QF2 Loss"] = np.mean(ptu.get_numpy(qf2_loss))
            self.eval_statistics["VF Loss"] = np.mean(ptu.get_numpy(vf_loss))
            self.eval_statistics["Policy Loss"] = np.mean(ptu.get_numpy(policy_loss))
            self.eval_statistics.update(
                create_stats_ordered_dict(
                    "Q1 Predictions",
                    ptu.get_numpy(q1_pred),
                )
            )
            self.eval_statistics.update(
                create_stats_ordered_dict(
                    "Q2 Predictions",
                    ptu.get_numpy(q2_pred),
                )
            )
            self.eval_statistics.update(
                create_stats_ordered_dict(
                    "V Predictions",
                    ptu.get_numpy(v_pred),
                )
            )
            self.eval_statistics.update(
                create_stats_ordered_dict(
                    "Log Pis",
                    ptu.get_numpy(log_pi),
                )
            )
            self.eval_statistics.update(
                create_stats_ordered_dict(
                    "Policy mu",
                    ptu.get_numpy(policy_mean),
                )
            )
            self.eval_statistics.update(
                create_stats_ordered_dict(
                    "Policy log std",
                    ptu.get_numpy(policy_log_std),
                )
            )

    @property
    def networks(self):
        return [
            self.policy,
            self.qf1,
            self.qf2,
            self.vf,
            self.target_vf,
        ]

    def _update_target_network(self):
        ptu.soft_update_from_to(self.vf, self.target_vf, self.soft_target_tau)

    def get_snapshot(self):
        return dict(
            qf1=self.qf1,
            qf2=self.qf2,
            policy=self.policy,
            vf=self.vf,
            target_vf=self.target_vf,
            policy_optimizer=self.policy_optimizer,
            qf1_optimizer=self.qf1_optimizer,
            qf2_optimizer=self.qf2_optimizer,
            vf_optimizer=self.vf_optimizer,
        )

    def load_snapshot(self, snapshot):
        self.qf1 = snapshot["qf1"]
        self.qf2 = snapshot["qf2"]
        self.policy = snapshot["policy"]
        self.vf = snapshot["vf"]
        self.target_vf = snapshot["target_vf"]
        self.policy_optimizer = snapshot["policy_optimizer"]
        self.qf1_optimizer = snapshot["qf1_optimizer"]
        self.qf2_optimizer = snapshot["qf2_optimizer"]
        self.vf_optimizer = snapshot["self.vf_optimizer"]

    def get_eval_statistics(self):
        return self.eval_statistics

    def end_epoch(self):
        self.eval_statistics = None

max_stars_repo_name: FineArtz/ILSwiss

max_issues_repo_name: FineArtz/ILSwiss

max_forks_repo_name: FineArtz/ILSwiss


text:
import os
import sys
import time
import numpy as np
from threading import Thread as thread
import subprocess

from matplotlib import pyplot as plt

def run_interaction_process():
    monkey_runner_path = '/home/vishal/Android/Sdk/tools/bin/monkeyrunner'
    event_path = '/home/vishal/Github_umass/efficient_cv/utility/main_runner/activity.py'
    out = subprocess.run([monkey_runner_path, event_path],
                         stdout=subprocess.DEVNULL)

class screen_interaction_runner:
    def __init__(self):
        self._is_running = True

    def terminate(self):
        self._is_running = False

    def run(self):
        while self._is_running:
            run_interaction_process()

def parse_encoded_out(s):
    """
    Format: lines with three characters represents
            (time_stamp, average_load, exp_load, raw_load)
            lines with two charatacters represents
            (time_stamp, runtime)
            The first line and any blank line will be discarded
    """
    load_timestamp = []
    avg_loads = []
    exp_loads = []
    raw_loads = []
    runtime_timestamps = []
    runtimes = []
    s_lines = s.splitlines()
    # print(s_lines)

    for index, line in enumerate(s_lines):
        if index == 0 or index == 1:
            continue
        if len(line) == 0:
            continue
        words = line.split()
        if len(words) == 4:
            load_timestamp.append(int(words[0]))
            avg_loads.append(float(words[1]))
            exp_loads.append(float(words[2]))
            raw_loads.append(float(words[3]))
        else:
            runtime_timestamps.append(int(words[0]))
            runtimes.append(float(words[1]))
    return (load_timestamp, avg_loads, exp_loads, raw_loads,
            runtime_timestamps, runtimes)

def parse_encoded_out_with_exit(s):
    load_timestamp = []
    avg_loads = []
    exp_loads = []
    raw_loads = []
    runtime_timestamps = []
    runtimes = []
    exits = []
    s_lines = s.splitlines()
    # print(s_lines)

    for index, line in enumerate(s_lines):
        if index == 0 or index == 1:
            continue
        if len(line) == 0:
            continue
        words = line.split()
        if len(words) == 4:
            load_timestamp.append(int(words[0]))
            avg_loads.append(float(words[1]))
            exp_loads.append(float(words[2]))
            raw_loads.append(float(words[3]))
        elif len(words) == 1:
            exits.append(int(words[0]))
        else:
            runtime_timestamps.append(int(words[0]))
            runtimes.append(float(words[1]))
    return (load_timestamp, avg_loads, exp_loads, raw_loads,
            runtime_timestamps, runtimes, exits)

def get_plots(stdout):
    start_index = 5
    encoded_out = str(stdout, encoding='utf-8')
    parsed_out = parse_encoded_out(encoded_out)
    (load_timestamp, avg_loads, exp_loads, raw_loads,
            runtime_timestamps, runtimes) = parsed_out

    # First dump all the date, and lets worry about it layer
    np.save("load_timestamp.npy", np.array(load_timestamp))
    np.save("avg_loads.npy", np.array(avg_loads))
    np.save("exp_loads.npy", np.array(exp_loads))
    np.save("raw_loads.npy", np.array(raw_loads))
    np.save("runtime_timestamps.npy", np.array(runtime_timestamps))
    np.save("runtimes.npy", np.array(runtimes))

    plt.plot(load_timestamp[start_index:], avg_loads[start_index:],
             label='avg_loads')
    plt.plot(load_timestamp[start_index:], exp_loads[start_index:],
             label='exp_loads')
    plt.plot(load_timestamp[start_index:], raw_loads[start_index:],
             label='instant_loads')
    #plt.legend(loc="upper right")
    #plt.show()

    #plt.clf()
    #plt.cla()
    #plt.close()
    runtimes = runtimes[start_index:]
    runtimes_arr = np.array(runtimes, dtype = np.float)
    diff = np.max(runtimes_arr) - np.min(runtimes_arr)
    runtimes_scaled = list((runtimes_arr - np.min(runtimes_arr))/diff)
    # print(runtime_timestamps)
    # print(runtimes)
    # print(runtimes_scaled)
    plt.plot(runtime_timestamps[start_index:], runtimes_scaled,
             label='runtimes', marker='o')
    plt.legend(loc="upper right")
    plt.show()

    plt.clf()
    plt.cla()
    plt.close()

def get_plots_dynamic(stdout, exit = 0):
    start_index = 5
    encoded_out = str(stdout, encoding='utf-8')
    if exit == 0:
        parsed_out = parse_encoded_out(encoded_out)
        (load_timestamp, avg_loads, exp_loads, raw_loads,
                runtime_timestamps, runtimes) = parsed_out
    else:
        parsed_out = parse_encoded_out_with_exit(encoded_out)
        (load_timestamp, avg_loads, exp_loads, raw_loads,
            runtime_timestamps, runtimes, exits) = parsed_out

    # First dump all the date, and lets worry about it layer
    np.save("load_timestamp.npy", np.array(load_timestamp))
    np.save("avg_loads.npy", np.array(avg_loads))
    np.save("exp_loads.npy", np.array(exp_loads))
    np.save("raw_loads.npy", np.array(raw_loads))
    np.save("runtime_timestamps.npy", np.array(runtime_timestamps))
    np.save("runtimes.npy", np.array(runtimes))
    if exit == 1:
        np.save("exits.npy", np.array(exits))
        print(exits)

    plt.plot(load_timestamp[start_index:], avg_loads[start_index:],
             label='avg_loads')
    plt.plot(load_timestamp[start_index:], exp_loads[start_index:],
             label='exp_loads')
    plt.plot(load_timestamp[start_index:], raw_loads[start_index:],
             label='instant_loads')
    #plt.legend(loc="upper right")
    #plt.show()

    #plt.clf()
    #plt.cla()
    #plt.close()
    runtimes = runtimes[start_index:]
    runtimes_arr = np.array(runtimes, dtype = np.float)
    diff = np.max(runtimes_arr) - np.min(runtimes_arr)
    runtimes_scaled = list((runtimes_arr - np.min(runtimes_arr))/diff)
    # print(runtime_timestamps)
    # print(runtimes)
    # print(runtimes_scaled)
    plt.plot(runtime_timestamps[start_index:], runtimes_scaled,
             label='runtimes', marker='o')
    plt.legend(loc="upper right")
    plt.show()

    plt.clf()
    plt.cla()
    plt.close()

def main():
    adb_command = '/home/vishal/Android/Sdk/platform-tools/adb'
    # program_file = '/data/local/tmp/milestone3/final_app'
    program_file = '/data/local/tmp/final_app_static/final_app_static'
    model_file = '/data/local/tmp/milestone3/mobilenet2_1.tflite'
    nr_runs = 1000
    alpha = 0.33
    print_selection = 1

    runner = screen_interaction_runner()
    t = thread(target = runner.run)
    t.start()

    time.sleep(20) # Wait

    print("Activity has started")
    # out = subprocess.Popen([adb_command, 'shell', program_file,
    #                         model_file, str(nr_runs)],
    #                         stdout=subprocess.PIPE,
    #                         stderr=subprocess.STDOUT)
    out = subprocess.Popen([adb_command, 'shell', program_file,
                            str(nr_runs), str(alpha), str(print_selection)],
                            stdout=subprocess.PIPE,
                            stderr=subprocess.STDOUT)
    stdout, _ = out.communicate()
    print("Now stopping the interaction")

    try:
        get_plots_dynamic(stdout, print_selection)
    except Exception as err:
        print(err)
        print(stdout)
    finally:
        pass

    runner.terminate()
    t.join()

    # Calculate mean and variance

    # Compute more statistics with 99th percentile and 95th percentile

if __name__ == "__main__":
    main()