text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import os
import argparse
import difflib
import pickle
from pathlib import Path

import numpy as np
import tables

import pyyaks.logger
import Ska.DBI
import Ska.File
from Chandra.Time import DateTime
from ska_helpers.run_info import log_run_info
from ska_helpers.retry import retry_call

from kadi.paths import IDX_CMDS_PATH, PARS_DICT_PATH
from kadi import __version__

MPLOGS_DIR = Path(os.environ['SKA'], 'data', 'mpcrit1', 'mplogs')

MIN_MATCHING_BLOCK_SIZE = 500
BACKSTOP_CACHE = {}
CMDS_DTYPE = [('idx', np.int32),
              ('date', '|S21'),
              ('type', '|S12'),
              ('tlmsid', '|S10'),
              ('scs', np.uint8),
              ('step', np.uint16),
              ('timeline_id', np.uint32),
              ('vcdu', np.int32)]

logger = None  # This is set as a global in main.  Define here for pyflakes.

class UpdatedDict(dict):
    """
    Dict with an ``n_updated`` attribute that gets incremented when any key value is set.
    """
    n_updated = 0

    def __setitem__(self, *args, **kwargs):
        self.n_updated += 1
        super(UpdatedDict, self).__setitem__(*args, **kwargs)

def get_opt(args=None):
    """
    Get options for command line interface to update_cmd_states.
    """
    parser = argparse.ArgumentParser(description='Update HDF5 cmds table')
    parser.add_argument("--mp-dir",
                        help=("MP load directory (default=/data/mpcrit1/mplogs) "
                              "or $SKA/data/mpcrit1/mplogs)"))
    parser.add_argument("--start",
                        help="Start date for update (default=stop-42 days)")
    parser.add_argument("--stop",
                        help="Stop date for update (default=Now+21 days)")
    parser.add_argument("--log-level",
                        type=int,
                        default=10,
                        help='Log level (10=debug, 20=info, 30=warnings)')
    parser.add_argument("--data-root",
                        default='.',
                        help="Data root (default='.')")
    parser.add_argument('--version', action='version',
                        version='%(prog)s {version}'.format(version=__version__))

    args = parser.parse_args(args)
    return args

def fix_nonload_cmds(nl_cmds):
    """
    Convert non-load commands commands dict format from Chandra.cmd_states
    to the values/structure needed here.  A typical value is shown below:
    {'cmd': u'SIMTRANS',               # Needs to be 'type'
    'date': u'2017:066:00:24:22.025',
    'id': 371228,                      # Store as params['nonload_id'] for provenence
    'msid': None,                      # Goes into params
    'params': {u'POS': -99616},
    'scs': None,                       # Set to 0
    'step': None,                      # Set to 0
    'time': 605233531.20899999,        # Ignored
    'timeline_id': None,               # Set to 0
    'tlmsid': None,                    # 'None' if None
    'vcdu': None},                     # Set to -1
    """
    new_cmds = []
    for cmd in nl_cmds:
        new_cmd = {}
        new_cmd['date'] = str(cmd['date'])
        new_cmd['type'] = str(cmd['cmd'])
        new_cmd['tlmsid'] = str(cmd['tlmsid'])
        for key in ('scs', 'step', 'timeline_id'):
            new_cmd[key] = 0
        new_cmd['vcdu'] = -1

        new_cmd['params'] = {}
        new_cmd['params']['nonload_id'] = int(cmd['id'])
        if cmd['msid'] is not None:
            new_cmd['params']['msid'] = str(cmd['msid'])

        # De-numpy (otherwise unpickling on PY3 has problems).
        if 'params' in cmd:
            params = new_cmd['params']
            for key, val in cmd['params'].items():
                key = str(key)
                try:
                    val = val.item()
                except AttributeError:
                    pass
                params[key] = val

        new_cmds.append(new_cmd)

    return new_cmds

def _tl_to_bs_cmds(tl_cmds, tl_id, db):
    """
    Convert the commands ``tl_cmds`` (numpy recarray) that occur in the
    timeline ``tl_id'' to a format mimicking backstop commands from
    Ska.ParseCM.read_backstop().  This includes reading parameter values
    from the ``db``.

    :param tl_cmds: numpy recarray of commands from timeline load segment
    :param tl_id: timeline id
    :param db: Ska.DBI db object

    :returns: list of command dicts
    """
    bs_cmds = [dict((col, row[col]) for col in tl_cmds.dtype.names)
               for row in tl_cmds]
    cmd_index = dict((x['id'], x) for x in bs_cmds)

    # Add 'params' dict of command parameter key=val pairs to each tl_cmd
    for par_table in ('cmd_intpars', 'cmd_fltpars'):
        tl_params = db.fetchall("SELECT * FROM %s WHERE timeline_id %s" %
                                (par_table,
                                 '= %d' % tl_id if tl_id else 'IS NULL'))

        # Build up the params dict for each command in timeline load segment
        for par in tl_params:
            # I.e. cmd_index[par.cmd_id]['params'][par.name] = par.value
            # but create the ['params'] dict as needed.
            if par.cmd_id in cmd_index:
                cmd_index[par.cmd_id].setdefault('params', {})[par.name] = par.value

    return bs_cmds

def get_cmds(start, stop, mp_dir=MPLOGS_DIR):
    """
    Get backstop commands corresponding to the supplied timeline load segments.
    The timeline load segments must be ordered by 'id'.

    Return cmds in the format defined by Ska.ParseCM.read_backstop().
    """
    mp_dir = str(mp_dir)

    # Get timeline_loads within date range.  Also get non-load commands
    # within the date range covered by the timelines.
    server = os.path.join(os.environ['SKA'], 'data', 'cmd_states', 'cmd_states.db3')
    with Ska.DBI.DBI(dbi='sqlite', server=server) as db:
        timeline_loads = db.fetchall("""SELECT * from timeline_loads
                                        WHERE datestop > '{}' AND datestart < '{}'
                                        ORDER BY id"""
                                     .format(start.date, stop.date))

        # Get non-load commands (from autonomous or ground SCS107, NSM, etc) in the
        # time range that the timelines span.
        tl_datestart = min(timeline_loads['datestart'])
        nl_cmds = db.fetchall('SELECT * from cmds where timeline_id IS NULL and '
                              'date >= "{}" and date <= "{}"'
                              .format(tl_datestart, stop.date))

        # Private method from cmd_states.py fetches the actual int/float param values
        # and returns list of dict.
        nl_cmds = _tl_to_bs_cmds(nl_cmds, None, db)
        nl_cmds = fix_nonload_cmds(nl_cmds)
        logger.info(f'Found {len(nl_cmds)} non-load commands between {tl_datestart} : {stop.date}')

    logger.info('Found {} timelines included within {} to {}'
                .format(len(timeline_loads), start.date, stop.date))

    if np.min(np.diff(timeline_loads['id'])) < 1:
        raise ValueError('Timeline loads id not monotonically increasing')

    cmds = []
    orbit_cmds = []
    orbit_cmd_files = set()

    for tl in timeline_loads:
        bs_file = Ska.File.get_globfiles(os.path.join(mp_dir + tl.mp_dir,
                                                      '*.backstop'))[0]
        if bs_file not in BACKSTOP_CACHE:
            bs_cmds = read_backstop(bs_file)
            logger.info('Read {} commands from {}'.format(len(bs_cmds), bs_file))
            BACKSTOP_CACHE[bs_file] = bs_cmds
        else:
            bs_cmds = BACKSTOP_CACHE[bs_file]

        # Process ORBPOINT (orbit event) pseudo-commands in backstop.  These
        # have scs=0 and need to be treated separately since during a replan
        # or shutdown we still want these ORBPOINT to be in the cmds archive
        # and not be excluded by timeline intervals.
        if bs_file not in orbit_cmd_files:
            bs_orbit_cmds = [x for x in bs_cmds if x['type'] == 'ORBPOINT']
            for orbit_cmd in bs_orbit_cmds:
                orbit_cmd['timeline_id'] = tl['id']
                if 'EVENT_TYPE' not in orbit_cmd['params']:
                    orbit_cmd['params']['EVENT_TYPE'] = orbit_cmd['params']['TYPE']
                    del orbit_cmd['params']['TYPE']
            orbit_cmds.extend(bs_orbit_cmds)
            orbit_cmd_files.add(bs_file)

        # Only store commands for this timeline (match SCS and date)
        bs_cmds = [x for x in bs_cmds
                   if tl['datestart'] <= x['date'] <= tl['datestop']
                   and x['scs'] == tl['scs']]

        for bs_cmd in bs_cmds:
            bs_cmd['timeline_id'] = tl['id']

        logger.info('  Got {} backstop commands for timeline_id={} and SCS={}'
                    .format(len(bs_cmds), tl['id'], tl['scs']))
        cmds.extend(bs_cmds)

    orbit_cmds = get_unique_orbit_cmds(orbit_cmds)
    logger.debug('Read total of {} orbit commands'
                 .format(len(orbit_cmds)))

    cmds.extend(nl_cmds)
    cmds.extend(orbit_cmds)

    # Sort by date and SCS step number.
    cmds = sorted(cmds, key=lambda y: (y['date'], y['step']))
    logger.debug('Read total of {} commands ({} non-load commands)'
                 .format(len(cmds), len(nl_cmds)))

    return cmds

def get_unique_orbit_cmds(orbit_cmds):
    """
    Given list of ``orbit_cmds`` find the quasi-unique set.  In the event of a
    replan/reopen or other schedule oddity, it can happen that there are multiple cmds
    that describe the same orbit event.  Since the detailed timing might change between
    schedule runs, cmds are considered the same if the date is within 3 minutes.
    """
    if len(orbit_cmds) == 0:
        return []

    # Sort by (event_type, date)
    orbit_cmds.sort(key=lambda y: (y['params']['EVENT_TYPE'], y['date']))

    uniq_cmds = [orbit_cmds[0]]
    # Step through one at a time and add to uniq_cmds only if the candidate is
    # "different" from uniq_cmds[-1].
    for cmd in orbit_cmds:
        last_cmd = uniq_cmds[-1]
        if (cmd['params']['EVENT_TYPE'] == last_cmd['params']['EVENT_TYPE']
                and abs(DateTime(cmd['date']).secs - DateTime(last_cmd['date']).secs) < 180):
            # Same event as last (even if date is a bit different).  Now if this one
            # has a larger timeline_id that means it is from a more recent schedule, so
            # use that one.
            if cmd['timeline_id'] > last_cmd['timeline_id']:
                uniq_cmds[-1] = cmd
        else:
            uniq_cmds.append(cmd)

    uniq_cmds.sort(key=lambda y: y['date'])

    return uniq_cmds

def get_idx_cmds(cmds, pars_dict):
    """
    For the input `cmds` (list of dicts), convert to the indexed command format where
    parameters are specified as an index into `pars_dict`, a dict of unique parameter
    values.

    Returns `idx_cmds` as a list of tuples:
       (par_idx, date, time, cmd, tlmsid, scs, step, timeline_id, vcdu)
    """
    idx_cmds = []

    for i, cmd in enumerate(cmds):
        if i % 10000 == 9999:
            logger.info('   Iteration {}'.format(i))

        # Define a consistently ordered tuple that has all command parameter information
        pars = cmd['params']
        keys = set(pars.keys()) - set(('SCS', 'STEP', 'TLMSID'))
        if cmd['tlmsid'] == 'AOSTRCAT':
            # Skip star catalog command because that has many (uninteresting) parameters
            # and increases the file size and load speed by an order of magnitude.
            pars_tup = ()
        else:
            pars_tup = tuple((key.lower(), pars[key]) for key in sorted(keys))

        try:
            par_idx = pars_dict[pars_tup]
        except KeyError:
            # Along with transition to 32-bit idx in #190, ensure that idx=65535
            # never gets used. Prior to #190 this value was being used by
            # get_cmds_from_backstop() assuming that it will never occur as a
            # key in the pars_dict. Adding 65536 allows older versions to work
            # with the new cmds.pkl pars_dict.
            par_idx = len(pars_dict) + 65536
            pars_dict[pars_tup] = par_idx

        idx_cmds.append((par_idx, cmd['date'], cmd['type'], cmd.get('tlmsid'),
                         cmd['scs'], cmd['step'], cmd['timeline_id'], cmd['vcdu']))

    return idx_cmds

def add_h5_cmds(h5file, idx_cmds):
    """
    Add `idx_cmds` to HDF5 file `h5file` of indexed spacecraft commands.
    If file does not exist then create it.
    """
    # Note: reading this file uncompressed is about 5 times faster, so sacrifice file size
    # for read speed and do not use compression.
    h5 = retry_call(tables.open_file, [h5file], {'mode': 'a'},
                    tries=4, delay=1, backoff=4)

    # Convert cmds (list of tuples) to numpy structured array.  This also works for an
    # existing structured array.
    cmds = np.array(idx_cmds, dtype=CMDS_DTYPE)

    # TODO : make sure that changes in non-load commands triggers an update

    try:
        h5d = h5.root.data
        logger.info('Opened h5 cmds table {}'.format(h5file))
    except tables.NoSuchNodeError:
        h5.create_table(h5.root, 'data', cmds, "cmds", expectedrows=2e6)
        logger.info('Created h5 cmds table {}'.format(h5file))
    else:
        date0 = min(idx_cmd[1] for idx_cmd in idx_cmds)
        h5_date = h5d.cols.date[:]
        idx_recent = np.searchsorted(h5_date, date0)
        logger.info('Selecting commands from h5d[{}:]'.format(idx_recent))
        logger.info('  {}'.format(str(h5d[idx_recent])))
        h5d_recent = h5d[idx_recent:]  # recent h5d entries

        # Define the column names that specify a complete and unique row
        key_names = ('date', 'type', 'tlmsid', 'scs', 'step', 'timeline_id', 'vcdu')

        h5d_recent_vals = [tuple(
            row[x].decode('ascii') if isinstance(row[x], bytes) else str(row[x])
            for x in key_names)
            for row in h5d_recent]
        idx_cmds_vals = [tuple(str(x) for x in row[1:]) for row in idx_cmds]

        diff = difflib.SequenceMatcher(a=h5d_recent_vals, b=idx_cmds_vals, autojunk=False)
        blocks = diff.get_matching_blocks()
        logger.info('Matching blocks for existing HDF5 and timeline commands')
        for block in blocks:
            logger.info('  {}'.format(block))
        opcodes = diff.get_opcodes()
        logger.info('Diffs between existing HDF5 and timeline commands')
        for opcode in opcodes:
            logger.info('  {}'.format(opcode))
        # Find the first matching block that is sufficiently long
        for block in blocks:
            if block.size > MIN_MATCHING_BLOCK_SIZE:
                break
        else:
            raise ValueError('No matching blocks at least {} long'
                             .format(MIN_MATCHING_BLOCK_SIZE))

        # Index into idx_cmds at the end of the large matching block.  block.b is the
        # beginning of the match.
        idx_cmds_idx = block.b + block.size

        if idx_cmds_idx < len(cmds):
            # Index into h5d at the point of the first diff after the large matching block
            h5d_idx = block.a + block.size + idx_recent

            if h5d_idx < len(h5d):
                logger.debug('Deleted relative cmds indexes {} .. {}'.format(h5d_idx - idx_recent,
                                                                             len(h5d) - idx_recent))
                logger.debug('Deleted cmds indexes {} .. {}'.format(h5d_idx, len(h5d)))
                h5d.truncate(h5d_idx)

            h5d.append(cmds[idx_cmds_idx:])
            logger.info('Added {} commands to HDF5 cmds table'.format(len(cmds[idx_cmds_idx:])))
        else:
            logger.info('No new timeline commands, HDF5 cmds table not updated')

    h5.flush()
    logger.info('Upated HDF5 cmds table {}'.format(h5file))
    h5.close()

def main(args=None):
    global logger

    opt = get_opt(args)

    logger = pyyaks.logger.get_logger(name='kadi_update_cmds', level=opt.log_level,
                                      format="%(asctime)s %(message)s")

    log_run_info(logger.info, opt)

    # Set the global root data directory.  This gets used in ..paths to
    # construct file names.  The use of an env var is needed to allow
    # configurability of the root data directory within django.
    os.environ['KADI'] = os.path.abspath(opt.data_root)
    idx_cmds_path = IDX_CMDS_PATH()
    pars_dict_path = PARS_DICT_PATH()

    try:
        with open(pars_dict_path, 'rb') as fh:
            pars_dict = pickle.load(fh)
        logger.info('Read {} pars_dict values from {}'.format(len(pars_dict), pars_dict_path))
    except IOError:
        logger.info('No pars_dict file {} found, starting from empty dict'
                    .format(pars_dict_path))
        pars_dict = {}

    if not opt.mp_dir:
        for prefix in ('/', os.environ['SKA']):
            pth = Path(prefix, 'data', 'mpcrit1', 'mplogs')
            if pth.exists():
                opt.mp_dir = str(pth)
                break
        else:
            raise FileNotFoundError('no mission planning directories found (need --mp-dir)')
    logger.info(f'Using mission planning files at {opt.mp_dir}')

    # Recast as dict subclass that remembers if any element was updated
    pars_dict = UpdatedDict(pars_dict)

    stop = DateTime(opt.stop) if opt.stop else DateTime() + 21
    start = DateTime(opt.start) if opt.start else stop - 42

    cmds = get_cmds(start, stop, opt.mp_dir)
    idx_cmds = get_idx_cmds(cmds, pars_dict)
    add_h5_cmds(idx_cmds_path, idx_cmds)

    if pars_dict.n_updated > 0:
        with open(pars_dict_path, 'wb') as fh:
            # Dump as pickle, first converting pars_dict from UpdatedDict to
            # plain dict. The n_updated attribute is not used outside of
            # update_cmds.py. See historical note in NOTES.build for context.
            pickle.dump(dict(pars_dict), fh, protocol=2)
            logger.info('Wrote {} pars_dict values ({} new) to {}'
                        .format(len(pars_dict), pars_dict.n_updated, pars_dict_path))
    else:
        logger.info('pars_dict was unmodified, not writing')

def _coerce_type(val):
    """Coerce the supplied ``val`` (typically a string) into an int or float if
    possible, otherwise as a string.
    """
    try:
        val = int(val)
    except ValueError:
        try:
            val = float(val)
        except ValueError:
            val = str(val)
    return val

def parse_params(paramstr):
    """
    Parse parameters key1=val1,key2=val2,... from ``paramstr``

    Parameter values are cast to the first type (int, float, or str) that
    succeeds.

    :param paramstr: Comma separated string of key=val pairs
    :rtype: dict of key=val pairs
    """
    params = {}
    for opt in paramstr.split(','):
        try:
            key, val = opt.split('=')
            params[key] = val if key == 'HEX' else _coerce_type(val)
        except Exception:
            pass  # backstop has some quirks like blank or '??????' fields

    return params

def read_backstop(filename):
    """
    Read commands from backstop file.

    Create dict with keys as follows for each command.  ``paramstr`` is the
    actual string with comma-separated parameters and ``params`` is the
    corresponding dict of key=val pairs.

    :param filename: Backstop file name
    :returns: list of dict for each command
    """
    bs = []
    for bs_line in open(filename):
        bs_line = bs_line.replace(' ', '')
        date, vcdu, cmd_type, paramstr = [x for x in bs_line.split('|')]
        vcdu = int(vcdu[:-1])  # Get rid of final '0' from '8023268 0' (where space was stripped)
        params = parse_params(paramstr)
        bs.append({'date': date,
                   'type': cmd_type,
                   'params': params,
                   'tlmsid': params.get('TLMSID'),
                   'scs': params.get('SCS'),
                   'step': params.get('STEP'),
                   'vcdu': vcdu
                   })
    return bs

if __name__ == '__main__':
    main()

max_stars_repo_name: jzuhone/kadi

max_issues_repo_name: jzuhone/kadi

max_forks_repo_name: jzuhone/kadi


text:
from typing import Union, Tuple

import multiaug
import numpy as np
import scipy.ndimage
from multiaug.augmenters import meta

def _generate_bool_sequence(num: int, random_state: int) -> list:
    return [random_state.choice([True, False], 1)[0] for _ in range(num)]

def _merge_bool_sequences(proposed: list, mask: list) -> list:
    return [x and y for x, y in zip(proposed, mask)]

def _validate_angle(angle: Union[int, float]):
    assert (angle >= 0 and angle <= 360), "Angle not within valid range [0, 360], received {}".format(angle)

class Rotate3d(meta.Augmenter):
    '''
    Apply random rotation to 3D images.

    Parameters
    ----------
    angle : int or float or list
        The maximum angle by which to rotate the image.

            * If 'int' or 'float' then use the angle specified for all axes.
            * If 'list' then must be of format [angle_x, angle_y, angle_z].

    axes : list
        Flags for the x, y, z axes that determine which axes the image can be rotated about.

            * If 'True' then axis can be rotated about
            * If 'False' then axis cannot be rotated about

    interpolation : str
        Interpolation method to use.

            * If 'nearest' then use nearest interpolation.
    '''
    def __init__(self, angle: Union[int, float, list], axes: Tuple[bool, bool, bool] = [True, True, True], interpolation: str = 'nearest'):
        super(Rotate3d, self).__init__()

        if isinstance(angle, list):
            assert len(angle) == 3, "Please specify rotation angle for x, y and z dimensions, only {} angles provided".format(len(angle))
            self.x_max, self.y_max, self.z_max = angle
        elif isinstance(angle, float) or isinstance(angle, int):
            self.x_max = self.y_max = self.z_max = angle
        else:
            raise NotImplementedError("Angle must either be a list of length 3 or a single number, not {}".format(angle))

        _validate_angle(self.x_max)
        _validate_angle(self.y_max)
        _validate_angle(self.z_max)

        self.interpolation = interpolation
        self.active_axes = axes

        if sum(self.active_axes) == 0:
            raise RuntimeError("All axes have been deactivated, please activate atleast one axes for augmentation to take place")

    def apply(self, images: np.ndarray, row_ids: list) -> np.ndarray:
        '''
        Apply transformations to an entire batch.

        Parameters
        ----------
        images : np.ndarray
            Image batch of shape N x H x W x D.

        row_ids : list
            Indices of rows to rotate.

        Returns
        -------
        np.ndarray
            Batch of rotated images.
        '''
        rotated_images = []
        for image in images[row_ids]:
            rotated_images.append(self.apply_to_sample(image))

        return np.array(rotated_images)

    def apply_to_sample(self, image: np.ndarray) -> np.ndarray:
        '''
        Apply transformation to a single image. Randomly samples one or more active axes and applies
        a random rotation about that axes.

        Parameters
        ----------
        image : np.ndarray
            Image to transform.

        Returns
        -------
        np.ndarray
            Rotated image.
        '''
        which_axes = _generate_bool_sequence(3, self.random_state)
        which_axes = _merge_bool_sequences(which_axes, self.active_axes)

        while sum(which_axes) == 0: # at least rotate around one axes
            which_axes = _generate_bool_sequence(3, self.random_state)
            which_axes = _merge_bool_sequences(which_axes, self.active_axes)

        img = image.copy()
        # z-axis
        if which_axes[2]:
            angle = self.random_state.uniform(-self.z_max, self.z_max)
            img = scipy.ndimage.interpolation.rotate(img, angle, mode=self.interpolation, axes=(0, 1), reshape=False)
        # y-axis
        if which_axes[1]:
            angle = self.random_state.uniform(-self.y_max, self.y_max)
            img = scipy.ndimage.interpolation.rotate(img, angle, mode=self.interpolation, axes=(0, 2), reshape=False)
        # x-axis
        if which_axes[0]:
            angle = self.random_state.uniform(-self.x_max, self.x_max)
            img = scipy.ndimage.interpolation.rotate(img, angle, mode=self.interpolation, axes=(1, 2), reshape=False)

        return img

ext: py

max_stars_repo_name: Devin-Taylor/MultiAug

max_issues_repo_name: Devin-Taylor/MultiAug

max_forks_repo_name: Devin-Taylor/MultiAug

reason: import numpy,import scipy


text:
import numpy as np

import tofu as tf

def make_v0_from_v1():

    conf = tf.load_config('AUG-V1')

    lorder = [
        ('SBi', np.r_[range(0, 23)]),
        ('TPLT5', np.r_[range(12, 22), range(0, 6)]),
        ('TPLT4', np.r_[range(14, 23), range(0, 9)]),
        ('TPLT3', np.r_[range(11, 15), range(0, 5)]),
        ('TPLT2', np.r_[range(8, 12), range(0, 4)]),
        ('TPLT1', np.r_[range(14, 17), range(0, 4)]),
        ('TPRT5', np.r_[range(14, 16), range(0, 4)]),
        ('TPRT4', np.r_[range(8, 12), range(0, 4)]),
        ('TPRT3', np.r_[range(9, 15), 0]),
        ('TPRT2', np.r_[range(23, 29), range(0, 4)]),
        ('D2dBu1', np.r_[range(11, 22), 0]),
        ('D2dBu2', np.r_[10, 0]),
        ('D2dBu3', np.r_[10, 0, 1]),
        ('ICRHa', np.r_[range(3, 61), 0]),
        ('D2dBl2', np.r_[range(10, 14), 0]),
        ('D2dBl3', np.r_[range(9, 13), 0, 1]),
        ('D2dBG2', np.r_[range(9, 28), 0, 1]),
        ('D3BG1', np.r_[range(18, 20), 0, 1]),
        ('D3BG10', np.r_[range(11, 18)]),
        ('D2cdomR', np.r_[range(7, 14), 0, 1, 2, 3]),
        ('D2cdome', np.r_[range(5, 14), 0, 1, 2]),
        ('D2cdomL', np.r_[range(5, 14), 0]),
        ('D2cTPic', np.r_[range(7, 12), 0]),
        ('D2cTPi', np.r_[range(12, 18), 0]),
        ('D2cTPib', np.r_[range(7, 23), 0, 1, 2]),
        ('D2ci2', np.r_[range(4, 10), 0, 1]),
        ('D2ci1', np.r_[range(4, 12), 0, 1]),
    ]

    poly = np.concatenate(
        tuple([
            conf.dStruct['dObj']['PFC'][ss].Poly_closed[:, ind]
            for (ss, ind) in lorder
        ]),
        axis=1,
    )

    return poly, lorder

max_stars_repo_name: WinstonLHS/tofu

max_issues_repo_name: WinstonLHS/tofu

reason: import numpy


text:
#!/usr/bin/python

import sys
import os
import numpy as np
import datetime
import pickle

#
# for matplotlib (some problem preventing use of $PYTHONPATH)
#
sys.path.append("/afs/ipp-garching.mpg.de/common/soft/anaconda/amd64_generic/3/5.3.0/lib/python3.7/site-packages")

from matplotlib import pyplot as plotter

sys.path.append("/ptmp/mpp/prieck/atlasDM/dileptonReso/madgraph5atlasval_complete/madgraph5atlasval/excursion/")

import excursion
import excursion.optimize
import excursion.plotting.twodim as exPlot
import excursion.utils

import scan2d as scandetails
#import scan_gl_mzp as scandetails
import getMuLimitBatch

import dmUpperLimitsForTesting
import myTest

def plot_and_save(workdir, index, gps, X, y_list, scandetails, batchsize, doSaveInput = True):

    if doSaveInput:
        outFileXY = open(os.path.join(workdir, "xy_{}.pickle".format(index)), "wb")
        dirXY = { "X" : X,
                  "Y" : y_list,
                  "gps" : gps }
        pickle.dump(dirXY, outFileXY)
        outFileXY.close()
#        print("==================================================")
#        print("GP:")
#        print(gps[0].get_params())
#        print("==================================================")

    figure, axarr = plotter.subplots(1, 2, sharey=True)
    figure.set_size_inches(9.5, 3.5)
    plotter.title('Iteration {}'.format(index))
#    plotter.yscale("log")
    exPlot.plot(axarr, gps, X, y_list, scandetails, figure, batchsize)
    plotter.tight_layout()
    plotter.savefig(os.path.join(workdir,'update_{}.png'.format(str(index).zfill(3))), dpi = 300)
#    plotter.show()

#--------------------------------------------------

def main():

    np.warnings.filterwarnings('ignore')

    scandetails.truth_functions = [ getMuLimitBatch.getMuLimitBatch ]
    scandetails.thresholds      = [ 0.0 ]

    _n_init = 6
    X, y_list, gps = excursion.optimize.init(scandetails, n_init = _n_init, seed = 123)

    print("init X is {}".format(X))
    print("init y_list is {}".format(y_list))

    workdir = datetime.datetime.now().strftime('excursionTest_%Y-%m-%d-%H-%M-%S-2dmulti')
    os.mkdir(workdir)
    plot_and_save(workdir, 0, gps, X, y_list, scandetails, _n_init)

    N_UPDATES = 15

    for iScan in range(N_UPDATES):
        print("X is {}".format(X))
        print("y_list is {}".format(y_list))
        _batchsize = 6
        X_new, acq_vals = excursion.optimize.gridsearch(gps, X, scandetails, batchsize = _batchsize, resampling_frac = 1.0)
        print('Iteration {}. new x: {}'.format(iScan, X_new))
        y_new = [ scandetails.truth_functions[0](X_new) ]
        print('Iteration {}. Evaluted truth function to value: {}'.format(iScan, y_new))
        y_list[0] = np.concatenate( [y_list[0], y_new[0]] )
        X         = np.concatenate( [X, X_new] )
        gps       = [ excursion.get_gp(X, y_list[0]) ]
        plot_and_save(workdir, iScan + 1, gps, X, y_list, scandetails, _batchsize)

    outFile = open(os.path.join(workdir, "upperLimits.txt"), "w")
    for iPoint in range(len(y_list[0])):
        for iCoordinate in range(X.shape[-1]):
            outFile.write("{} ".format(X[iPoint][iCoordinate]))
        outFile.write(" -> {}\n".format(y_list[0][iPoint]))

#--------------------------------------------------

if __name__ == "__main__":
    main()

max_stars_repo_name: irinaespejo/excursion

max_issues_repo_name: irinaespejo/excursion

max_forks_repo_name: irinaespejo/excursion


text:
#! /usr/bin/env python3
#
#  Copyright 2019 California Institute of Technology
#
#  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.
#
# ISOFIT: Imaging Spectrometer Optimal FITting
# Author: David R Thompson, david.r.thompson@jpl.nasa.gov
#

from os.path import split, abspath
import numpy as np
import scipy
from scipy.ndimage.filters import gaussian_filter1d
from spectral.io import envi

from isofit.core.common import expand_path, json_load_ascii

def instrument_model(config):
    """."""

    hdr_template = """ENVI
    samples = {samples}
    lines   = {lines}
    bands   = 1
    header offset = 0
    file type = ENVI Standard
    data type = 4
    interleave = bsq
    byte order = 0
    """

    config = json_load_ascii(config, shell_replace=True)
    configdir, configfile = split(abspath(config))

    infile = expand_path(configdir, config['input_radiance_file'])
    outfile = expand_path(configdir, config['output_model_file'])
    flatfile = expand_path(configdir, config['output_flatfield_file'])
    uniformity_thresh = float(config['uniformity_threshold'])

    infile_hdr = infile + '.hdr'
    img = envi.open(infile_hdr, infile)
    inmm = img.open_memmap(interleave='bil', writable=False)
    X = np.array(inmm[:, :, :], dtype=np.float32)
    nr, nb, nc = X.shape

    FF, Xhoriz, Xhorizp, use_ff = _flat_field(X, uniformity_thresh)
    np.array(FF, dtype=np.float32).tofile(flatfile)
    with open(flatfile+'.hdr', 'w') as fout:
        fout.write(hdr_template.format(lines=nb, samples=nc))

    C, Xvert, Xvertp, use_C = _column_covariances(X, uniformity_thresh)
    cshape = (C.shape[0], C.shape[1]**2)
    out = np.array(C, dtype=np.float32).reshape(cshape)
    mdict = {'columns': out.shape[0], 'bands': out.shape[1],
             'covariances': out, 'Xvert': Xvert, 'Xhoriz': Xhoriz,
             'Xvertp': Xvertp, 'Xhorizp': Xhorizp, 'use_ff': use_ff,
             'use_C': use_C}
    scipy.io.savemat(outfile, mdict)

def _high_frequency_vert(X, sigma=4.0):
    """."""

    nl, nb, nr = X.shape
    Xvert = X.copy()
    for r in range(nr):
        for b in range(nb):
            filt = gaussian_filter1d(Xvert[:, b, r], sigma, mode='nearest')
            Xvert[:, b, r] = X[:, b, r] - filt
    return Xvert

def _low_frequency_horiz(X, sigma=4.0):
    """."""

    nl, nb, nr = X.shape
    Xhoriz = X.copy()
    for l in range(nl):
        for b in range(nb):
            Xhoriz[l, b, :] = gaussian_filter1d(
                Xhoriz[l, b, :], sigma, mode='nearest')
    return Xhoriz

def _flat_field(X, uniformity_thresh):
    """."""

    Xhoriz = _low_frequency_horiz(X, sigma=4.0)
    Xhorizp = _low_frequency_horiz(X, sigma=3.0)
    nl, nb, nc = X.shape
    FF = np.zeros((nb, nc))
    use_ff = np.ones((X.shape[0], X.shape[2])) > 0
    for b in range(nb):
        xsub = Xhoriz[:, b, :]
        mu = xsub.mean(axis=0)
        dists = abs(xsub - mu)
        thresh = np.percentile(dists.flatten(), 90.0)
        use = dists < thresh
        FF[b, :] = ((xsub*use).sum(axis=0)/use.sum(axis=0)) / \
            ((X[:, b, :]*use).sum(axis=0)/use.sum(axis=0))
        use_ff = np.logical_and(use_ff, use)
    return FF, Xhoriz, Xhorizp, np.array(use_ff)

def _column_covariances(X, uniformity_thresh):
    """."""

    Xvert = _high_frequency_vert(X, sigma=4.0)
    Xvertp = _high_frequency_vert(X, sigma=3.0)
    models = []
    use_C = []
    for i in range(X.shape[2]):
        xsub = Xvert[:, :, i]
        mu = xsub.mean(axis=0)
        dists = np.sqrt(pow((xsub - mu), 2).sum(axis=1))
        thresh = np.percentile(dists, 95.0)
        use = dists < thresh
        C = np.cov(xsub[use, :], rowvar=False)
        [U, V, D] = scipy.linalg.svd(C)
        V[V < 1e-8] = 1e-8
        C = U.dot(np.diagflat(V)).dot(D)
        models.append(C)
        use_C.append(use)
    return np.array(models), Xvert, Xvertp, np.array(use_C).T

max_stars_repo_name: jfahlen/isofit

max_issues_repo_name: jfahlen/isofit

max_forks_repo_name: serbinsh/isofit

reason: import numpy,import scipy,from scipy


text:
# Copyright (c) Microsoft Corporation. All rights reserved.
# Licensed under the MIT License.

import numpy as np
import torch

from .fm_stage_wrapper import FM_Stage_Wrapper
from .modules.spm_tracker_cm import CM_Stage

from ..utils import (generate_anchors_on_response_map, rpn_logit_to_prob,
                     load_weights, crop_with_boxes)
from ..utils import boxes as ubox

class SPMTrackerInference(object):

    def __init__(self, model_device, img_norm, cfg):

        self.model_device = model_device
        self.cfg = cfg

        # create and initialize CM stage & FM stage
        self.cm = CM_Stage(cfg)
        load_weights(self.cm, self.cfg.MODEL.CM_WEIGHTS_FILE, verbose=False)
        self.cm.to(self.model_device)
        self.cm.eval()
        self.fm_wrapper = FM_Stage_Wrapper(cfg, model_device)

        self.z_size = self.cfg.INFERENCE.Z_SIZE
        self.x_size = self.cfg.INFERENCE.X_SIZE

        self.anchors = generate_anchors_on_response_map(
            z_size=self.z_size,
            z_response_size=self.cm.z_res_size,
            x_response_size=self.cm.x_res_size,
            model_stride=self.cm.stride,
            ratios=self.cfg.MODEL.ANCHOR_RATIOS,
            scales=self.cfg.MODEL.ANCHOR_SCALES)

        self.norm = img_norm

    def init_tracker(self, first_frame_img, current_box):
        """ Init the tracker (called in the first frame) """
        with torch.no_grad():
            search_img_tensor, z_box, x_box, scale_x = self.generate_search_image(first_frame_img, current_box)
            z_out = self.cm(z_image=search_img_tensor, x_image=None)

            # init fine matching stage
            z_box_in_search_img = \
                ubox.image_to_search_region(np.array([0.] + current_box, dtype=np.float32), x_box, scale_x)
            z_rois = torch.tensor(z_box_in_search_img).type_as(search_img_tensor).view(1, 5)
            z_feats = self.fm_wrapper.fm.feat_extractor(z_out, z_rois)
            self.fm_wrapper.set_template(z_feats)

    def predict_next_frame(self, img_tensor, current_box, vis=False):
        """ Predict the tracking result in given image
        Args:
            img_tensor: torch.Tensor in shape of [1, 3, H, W]
            current_box: list, in order of [x1, y1, x2, y2]
            vis: boolean, visualize tracking results or not (for debug)
        Returns:
            results: list in order of [x1, y1, x2, y2]
        """
        with torch.no_grad():
            search_img_tensor, z_box, x_box, scale_x = self.generate_search_image(img_tensor, current_box)
            det_boxes, x_out = self.generate_rpn_boxes(search_img_tensor, x_box, scale_x)

        proposed_boxes = self.generate_proposal_boxes(det_boxes, current_box)

        fm_scores, fm_boxes = self.fm_wrapper.get_similarity(x_out, proposed_boxes[:, 5:9])
        win_scores = self.calculate_window_scores(proposed_boxes)

        if self.cfg.INFERENCE.ENABLE_FM_BBOX_REFINE:
            fm_boxes_img = ubox.search_region_to_image(fm_boxes, x_box, scale_x)

            # calculate the weights
            v_influence = self.cfg.INFERENCE.FM_BBOX_INFLUENCE
            sum_scores = np.clip(proposed_boxes[:, 4] + fm_scores * v_influence, a_min=0.00001, a_max=99999.)
            _w = (proposed_boxes[:, 4] / sum_scores)[:, np.newaxis]

            proposed_boxes[:, :4] = proposed_boxes[:, :4] * _w + fm_boxes_img[:, :4] * (1. - _w)
            proposed_boxes[:, 5:9] = proposed_boxes[:, 5:9] * _w + fm_boxes[:, :4] * (1. - _w)

        cm_scores = proposed_boxes[:, 4]
        pred_scores = cm_scores * (1 - self.cfg.INFERENCE.FM_SCORE_INFLUENCE) + \
                      fm_scores * self.cfg.INFERENCE.FM_SCORE_INFLUENCE

        final_scores = pred_scores * (1 - self.cfg.INFERENCE.WINDOWS_INFLUENCE) + \
                       win_scores * self.cfg.INFERENCE.WINDOWS_INFLUENCE

        max_score_idx = np.argmax(final_scores)

        best_box = proposed_boxes[max_score_idx, :4].tolist()
        best_score = pred_scores[max_score_idx]

        if best_score < self.cfg.INFERENCE.MAX_SCORE_THRESHOLD:
            best_box = current_box
        else:
            lr = best_score * self.cfg.INFERENCE.LR
            best_box = self.linear_inter(best_box, current_box, lr, img_tensor.size(2), img_tensor.size(3))

        return best_box

    def generate_rpn_boxes(self, img_tensor, x_crop_box, scale_x):
        """ Forward SiamRPN network and get the predicted boxes.
        Args:
            img_tensor: torch.Tensor in shape of [1, 3, x_size, x_size], search image tensor
        Returns:
            det_boxes: numpy.ndarray in shape of [x_res_size * x_res_size * num_anchors, 9]. The first four values
                       indicate the coordinates in source image domain. The following one value denotes the rpn score.
                       the last four values are coordinates in search region domain.
            x_out: dict. Convolutional features of search image.

        """
        z_out, x_out, rpn_cls_logit, rpn_bbox_delta = self.cm(z_image=None, x_image=img_tensor)
        rpn_cls_prob = \
            rpn_logit_to_prob(rpn_cls_logit)
        rpn_cls_prob = rpn_cls_prob[:, rpn_cls_prob.size(1) // 2:, :, :]

        rpn_bbox_delta = rpn_bbox_delta.data.cpu().numpy()[0].transpose((1, 2, 0)).reshape(-1, 4)
        rpn_score = rpn_cls_prob.data.cpu().numpy()[0].transpose((1, 2, 0)).flatten()
        rpn_bbox_pred = ubox.bbox_transform(self.anchors, rpn_bbox_delta)

        # transform to original image domain
        rpn_bbox_pred_img = ubox.search_region_to_image(rpn_bbox_pred, x_crop_box, scale_x)

        # combine the coordinates in two domains
        det_boxes = np.hstack((rpn_bbox_pred_img, rpn_score[:, np.newaxis], rpn_bbox_pred))

        return det_boxes, x_out

    def generate_proposal_boxes(self, rpn_boxes, current_box):
        """ Propose some candidate boxes.
        Args:
            rpn_boxes: numpy.ndarray in shape of [N, 9], see 'generate_rpn_boxes' function for more details.
            current_box: list in order of [x1, y1, x2, y2]
        Returns:
            proposed_boxes: numpy.ndarray in shape of [K, 9], the proposal results
        """

        num_near_candidates = self.cfg.INFERENCE.NUM_NEAR_CANDIDATES
        num_candidates = self.cfg.INFERENCE.NUM_CANDIDATES

        current_box_xyxy = np.array(current_box, dtype=np.float32)

        # Step 1. Find the candidates which are closest to the current box
        if num_near_candidates > 0:
            overlaps = ubox.bbox_overlaps(np.ascontiguousarray(rpn_boxes[:, :4], dtype=np.float32),
                                          np.array(current_box_xyxy, dtype=np.float32).reshape(1, 4)).flatten()
            is_near_box = (overlaps > 0.80)
            near_inds = np.where(is_near_box)[0]

            # if there exist a few boxes which is close to current box
            # we will select the one with the highest rpn score.
            if len(near_inds) > num_near_candidates:
                ov_sorted_inds = near_inds[np.argsort(rpn_boxes[near_inds, 4])[::-1]]
                near_boxes = rpn_boxes[ov_sorted_inds[:num_near_candidates], :]

                other_inds = np.where(~is_near_box)[0]
                other_inds = np.concatenate((ov_sorted_inds[num_near_candidates:], other_inds))
                proposed_boxes = rpn_boxes[other_inds, :]
            else:
                ov_sorted_inds = np.argsort(overlaps)[::-1]
                near_boxes = rpn_boxes[ov_sorted_inds[:num_near_candidates], :]
                proposed_boxes = rpn_boxes[ov_sorted_inds[num_near_candidates:], :]
        else:
            proposed_boxes = rpn_boxes

        # Step 2. filter some box whose score is too low.
        max_score = np.max(proposed_boxes[:, 4])
        kept_inds = np.where(proposed_boxes[:, 4] >= self.cfg.INFERENCE.SCORE_THRESH_FACTOR * max_score)[0]
        proposed_boxes = proposed_boxes[kept_inds, :]

        # Step 3. apply for NMS
        if proposed_boxes.shape[0] > 1:
            kept_inds = ubox.nms(np.ascontiguousarray(proposed_boxes[:, :5], dtype=np.float32),
                                 self.cfg.INFERENCE.NMS_THRESHOLD)
            proposed_boxes = proposed_boxes[kept_inds, :]
            sort_inds = np.argsort(proposed_boxes[:, 4])[::-1]
            proposed_boxes = proposed_boxes[sort_inds, :]

            _k = min(proposed_boxes.shape[0], num_candidates - num_near_candidates)
            proposed_boxes = proposed_boxes[:_k, :]

        if num_near_candidates > 0:
            proposed_boxes = np.vstack((proposed_boxes, near_boxes))

        return proposed_boxes

    def calculate_window_scores(self, det_boxes):
        """ Add cosine windows penalty.
        Args:
            det_boxes: numpy.ndarray
        """
        det_boxes_on_x_xcycwh = ubox.xyxy_to_xcycwh(det_boxes[:, 5:9])
        center_dist = det_boxes_on_x_xcycwh[:, 0:2] / (self.x_size - 1)
        center_dist = np.maximum(np.minimum(center_dist, 1.0), 0.0)
        win_scores_2d = 0.5 - 0.5 * np.cos(2. * np.pi * center_dist)
        win_scores = win_scores_2d[:, 0] * win_scores_2d[:, 1]
        return win_scores

    def linear_inter(self, det_box, current_box, lr, im_height, im_width):
        """ To make tracking more stable, the final result is a linear combination between
            predicted result and last frame result.

        Args:
            det_box: list in order of [x1, y1, x2, y2], predicted result
            current_box: list in order of [x1, y1, x2, y2], predicted result
            lr: linear combination weight
        Results:
            refine_box: list in order of [x1, y1, x2, y2]

        """
        _, _, current_w, current_h = ubox.xyxy_to_xcycwh(current_box)
        refine_box = ubox.xyxy_to_xcycwh(det_box)
        refine_box[2] = refine_box[2] * lr + current_w * (1 - lr)
        refine_box[3] = refine_box[3] * lr + current_h * (1 - lr)

        refine_box[2] = max(10, refine_box[2])
        refine_box[3] = max(10, refine_box[3])
        refine_box[0] = max(0, min(refine_box[0], im_width - (refine_box[2]) / 2.))
        refine_box[1] = max(0, min(refine_box[1], im_height - (refine_box[3]) / 2.))

        refine_box = ubox.xcycwh_to_xyxy(refine_box)

        return refine_box

    def generate_search_image(self, img_tensor, current_box):
        # generate search regions
        z_box, x_box, scale_x = ubox.generate_search_region_boxes(
            boxes=current_box,
            z_size=self.z_size,
            x_size=self.x_size,
            context_amount=self.cfg.MODEL.Z_CONTEXT_AMOUNT)
        search_img_tensor = \
            crop_with_boxes(img_tensor, x_box, self.x_size, self.x_size)

        search_img_tensor = self.norm(search_img_tensor)

        return search_img_tensor, z_box, x_box, scale_x

max_stars_repo_name: songheony/SPM-Tracker

max_issues_repo_name: songheony/SPM-Tracker

max_forks_repo_name: songheony/SPM-Tracker


text:
import unittest
import filterdesigner.IIRDesign as IIRDesign
import scipy.signal as signal
import numpy as np

class TestButtord(unittest.TestCase):

    def setUp(self):
        self.f1 = 0.2
        self.f2 = 0.3
        self.f3 = [0.2, 0.5]
        self.f4 = [0.1, 0.6]
        self.Rp = 1
        self.Rs = 1

    def test_buttord_1(self):
        # Test case for lowpass filter
        self.assertTrue(np.all(IIRDesign.buttord(self.f1, self.f2, self.Rp, self.Rs) == signal.buttord(self.f1, self.f2, self.Rp, self.Rs, analog=False, fs=2)))

    def test_buttord_2(self):
        # Test case for highpass filter
         self.assertTrue(np.all(IIRDesign.buttord(self.f2, self.f1, self.Rp, self.Rs) == signal.buttord(self.f2, self.f1, self.Rp, self.Rs, analog=False, fs=2)))

    def test_buttord_3(self):
        # Test case for bandpass filter
        ORD = IIRDesign.buttord(self.f3, self.f4, self.Rp, self.Rs)
        ord = signal.buttord(self.f3, self.f4, self.Rp, self.Rs, analog=False, fs=2)
        self.assertTrue((ORD[0] == ord[0]) and (ORD[1] == ord[1]).all())

    def test_buttord_4(self):
        # Test case for bandstop filter
        ORD = IIRDesign.buttord(self.f4, self.f3, self.Rp, self.Rs)
        ord = signal.buttord(self.f4, self.f3, self.Rp, self.Rs, analog=False, fs=2)
        self.assertTrue((ORD[0] == ord[0]) and (ORD[1] == ord[1]).all())

    def test_buttord_5(self):
        # Test case for analog filter
        self.assertTrue(np.all(IIRDesign.buttord(60, 75, self.Rp, self.Rs, zs='s') == signal.buttord(60, 75, self.Rp, self.Rs, analog=True, fs=None)))

    def test_buttord_6(self):
        # Test case for Exception
        with self.assertRaises(ValueError):
            IIRDesign.buttord(self.f1, self.f2, self.Rp, self.Rs, zs='x')

max_stars_repo_name: Yuki-F-HCU/filterdesigner

max_issues_repo_name: Yuki-F-HCU/filterdesigner

max_forks_repo_name: Yuki-F-HCU/filterdesigner


text:
# 파넬백 광학 흐름 변환으로 영상을 만들어주는 파일

import cv2
import numpy as np

cap = cv2.VideoCapture("./data/squat_wrong1.mp4")

fourcc = cv2.VideoWriter_fourcc('M', 'J', 'P', 'G')
out = cv2.VideoWriter('./data/optical_squat_wrong1.avi', fourcc, 20.0, (int(cap.get(3)), int(cap.get(4))))

ret, frame1 = cap.read()
prvs = cv2.cvtColor(frame1, cv2.COLOR_BGR2GRAY)
hsv = np.zeros_like(frame1)
hsv[..., 1] = 255
i = 0
prediction_str = ""
repetitions = 0
up = 0
down = 0
no_move = 0
current_move = 0
initial = -1
while (cap.isOpened()):
    i += 1

    ret, frame2 = cap.read()
    if not (ret): break
    next = cv2.cvtColor(frame2, cv2.COLOR_BGR2GRAY)

    flow = cv2.calcOpticalFlowFarneback(prvs, next, None, 0.5, 3, 15, 3, 5, 1.2, 0)

    mag, ang = cv2.cartToPolar(flow[..., 0], flow[..., 1])
    hsv[..., 0] = ang * 180 / np.pi / 2
    hsv[..., 2] = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX)
    rgb = cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)

    cv2.imshow('Mediapipe Feed', rgb)
    if cv2.waitKey(10) & 0xFF == ord('q'):
        break
    out.write(rgb)
    prvs = next

print("Video Generated")
out.release()
cap.release()
cv2.destroyAllWindows()

max_stars_repo_path: create data(OpticalFlow_Farneback).py

max_stars_repo_name: Y-greatigr/Muscle-Loss-Prevention

max_issues_repo_path: create data(OpticalFlow_Farneback).py

max_issues_repo_name: Y-greatigr/Muscle-Loss-Prevention

max_forks_repo_path: create data(OpticalFlow_Farneback).py

max_forks_repo_name: Y-greatigr/Muscle-Loss-Prevention


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

import os
from typing import List, Tuple
import pandas as pd
import numpy as np
import datetime as dt
import py_scripts.transform
import pickle

MODELSDIR = r'../models'
DATASET = r'../data/sim_ts_limpo.csv'

class modelo_produtos:
    """ encapsulador para os modelos de cada produto """

    def __init__(self, dataset: str = DATASET, modelsdir: str = MODELSDIR):
        self.dataset = dataset
        self.modelsdir = modelsdir

        # importanto dados limpos
        ts_raw = pd.read_csv(self.dataset)
        tsd, self.tswide = py_scripts.transform.pipeline(ts_raw)

        # produtos
        self.produtos = self.tswide.columns

        # importar modelos
        self.modelo, self.serie_treino = self.importar_modelos()

        self.fat_total = self.tswide.sum(axis = 'columns')

    def importar_modelos(self):

        modelo = {}
        serie_treino = {}

        for produto in self.produtos:
            produto_  = produto.split('_')[0]

            picklefile = fr'produto_{produto_}.model'
            picklefn = os.path.join(os.path.abspath(self.modelsdir), picklefile)

            with open(picklefn, 'rb') as modelo_arq:
                unpickler = pickle.Unpickler(modelo_arq)
                modelo_dict = unpickler.load()
                modelo[produto] = modelo_dict['modelo']
                serie_treino[produto] = modelo_dict['serie_treino']

        return modelo, serie_treino

    def get_models(self):
        return self.modelo, self.serie_treino

    def get_test_begin(self, produtos: List or None = None):

        if produtos is None:
            produtos = self.produtos

        serie_treino_prods = { p: s for p, s in self.serie_treino.items() if p in produtos }

        train_end = pd.Series(
            [ v.index[-1] for v in serie_treino_prods.values() ],
            index = produtos
        )

        test_start = train_end + pd.offsets.MonthBegin(1)

        return test_start

    def get_all_test_begin(self):

        train_end = max([ v.index[-1] for v in self.serie_treino.values() ])

        test_start = train_end + pd.offsets.MonthBegin(1)

        return test_start

    def predict(self, n_periods: int, return_conf_int: bool = False,
        predict_array: bool = True,
        *args, **kwards
    ) -> pd.Series or pd.DataFrame or np.array or Tuple[np.array, np.array]:

        if n_periods <= 0:
            raise ValueError('Can only predict forward!')

        # construimos um dataframe onde ficarão as predições individuais
        preds = pd.DataFrame([], columns = self.produtos)

        if return_conf_int:
            colsmult = pd.MultiIndex.from_product((self.produtos, ['lb', 'ub']))
            preds_ci = pd.DataFrame([], columns = colsmult)

        # obtemos a maior data em todos os conjuntos de treino

        max_train_right_bound = max([ v.index[-1] for v in self.serie_treino.values() ])

        for produto in self.produtos:
            # maior data de cada conjunto de treino
            train_right_bound = self.serie_treino[produto].index[-1]

            # construimos o índice de datas do conjunto de teste de cada produto:
            # range entre mês após o último contido no conjunto de treino e
            idx_test = pd.date_range(
                start = train_right_bound + dt.timedelta(days = 1),
                end = max_train_right_bound + pd.offsets.MonthBegin(n_periods), freq = 'MS')

            # geramos a predição para cada produto. Essa predição vem em um np.array
            # como queremos o intervalo de confiança, o resultado da função é uma tupla com
            #    - o array da predição média
            #    - um array com duas colunas contendo o lower bound e o upper bound
            arr_pred_all = self.modelo[produto].predict(n_periods = idx_test.shape[0], return_conf_int = return_conf_int)

            # primeiro trataremos das médias
            if return_conf_int:
                arr_pred = arr_pred_all[0]
            else:
                arr_pred = arr_pred_all

            # convertemos o array para Series
            pred = pd.Series(arr_pred, index = idx_test)
            pred.name = 'predicted_mean'

            # adicionamos a Series ao DataFrame `preds`
            preds[produto] = pred

            # agora trabalharemos nos bounds
            if return_conf_int:
                arr_pred_ci = arr_pred_all[1]

                pred_ci = pd.DataFrame(
                    arr_pred_ci,
                    columns = pd.MultiIndex.from_product(((produto, ), ('lb', 'ub'))),
                    index = idx_test
                )

                preds_ci[pred_ci.columns] = pred_ci

        preds_series = preds.dropna().sum(axis = 'columns')
        preds_series.name = 'predicted_mean'

        if return_conf_int:
            fat_test = pd.DataFrame([])

            fat_test['predicted_mean'] = preds_series

            fat_test['lb'] = preds_ci.loc[:, (slice(None), 'lb')].dropna().sum(axis = 'columns')
            fat_test['ub'] = preds_ci.loc[:, (slice(None), 'ub')].dropna().sum(axis = 'columns')

            if predict_array:
                return (
                    fat_test['predicted_mean'].values,
                    fat_test[['lb', 'ub']].values
                )
            else:
                return fat_test

        else:
            if predict_array:
                return preds_series.values
            else:
                return preds_series

    def __str__(self):
        totalstr = 'Modelos:'

        tamanho_campo = max([ len(produto) for produto in self.produtos ]) + 4

        for produto, modelo in self.modelo.items():
            totalstr += f"\n{produto:>{tamanho_campo}s}: {modelo}"

        return totalstr

    def __repr__(self):
        return str(self.modelo)

max_stars_repo_name: flimao/case-previsao-faturamento

max_issues_repo_name: flimao/case-previsao-faturamento

max_forks_repo_name: flimao/case-previsao-faturamento


text:
import numpy as np
from utilities.objective_functions import TSPObjectiveFunction
from perm_LQUBO.form_new_LQUBO import NewLQUBO
from perm_LQUBO.next_perm import NextPerm
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from scipy import stats
ax = plt.figure().gca()

class ObjectiveFunction:
    def __init__(self, objective_function=None):
        if objective_function:
            self.objective_function = objective_function
        else:
            raise AttributeError('Objective function missing.')

class LocalQUBOQuality(ObjectiveFunction):
    """
    For a specified objective function and LQUBO type, LocalQUBOQuality class has 2 useful visual functions that can
    be run from the command line.

    """
    def __init__(self,
                 objective_function=None,
                 num_points=None):
        super().__init__(objective_function=objective_function)

        # Initialize switch network:
        # The default behavior here is to choose the smaller of either permutation or
        # sorting networks for the given input size.
        self.n_problem = self.objective_function.n

        if num_points:
            self.num_points = num_points
        else:
            self.num_points = 200

        # Initialize random bitstring
        self.p = np.random.permutation(self.n_problem)
        self.v = self.objective_function(self.p)
        self.n_qubo = len(self.p) - 1

        if type(objective_function) == TSPObjectiveFunction:
            self.OF_type = 'TSP'
        else:
            self.OF_type = 'QAP'

        self.form_qubo = NewLQUBO(objective_function=self.objective_function).form_lqubo(self.p)

        self.qubo = np.zeros((self.n_qubo, self.n_qubo))

        for i in range(self.n_qubo):
            for j in range(self.n_qubo):
                if i <= j:
                    self.qubo[i][j] = self.form_qubo[(i, j)]

    def q_q_plot(self, hamming_dist):

        """
        Given a specified hamming dist (int), q_q_plot will produce a Q-Q scatter plot of the LQUBO change in objective
        function versus the actual change in objective function.
        """

        if hamming_dist:
            hamming_dist = hamming_dist
        else:
            raise AttributeError('Hamming dist missing.')

        def swap(input_list, first_entry, second_entry):
            input_list[first_entry], input_list[second_entry] = input_list[second_entry], input_list[first_entry]
            return input_list

        def random_binary(binary):
            for index in range(hamming_dist):
                swap(binary, index, np.random.randint(self.n_qubo))
            return binary

        # delta q of specified hamming dist
        hd = []
        x = []
        y = []

        for i in range(self.num_points):
            hd.append(np.zeros(self.n_qubo))

        for i in range(len(hd)):
            for j in range(hamming_dist):
                hd[i][j] = 1
            # only swap first n_qubo qubits to preserve total weight qubit
            random_binary(hd[i][:self.n_qubo])

        for vec in hd:
            p_new = NextPerm(current_perm=self.p, lqubo_result=vec).next_perm()
            v_new = self.objective_function(p_new)
            x.append(v_new - self.v)
            y.append(np.matmul(np.matmul(vec, self.qubo), np.transpose(vec)))

        plt.scatter(x, y, label='HD = {}'.format(hamming_dist))

        plt.xlabel("true delta obj")
        plt.ylabel("local qubo delta obj")
        plt.title('n = {} {}'.format(self.n_problem, self.OF_type))

        plt.axhline(y=0, color='k')
        plt.axvline(x=0, color='k')
        plt.legend(loc='upper left')
        plt.show()

    def plot_r_squared(self):

        """
        plot_r_squared function will take objective function from the base class and produce multiple scatter plots to
        retrieve the R squared value for a specified hamming distance.  Then for a hamming dist from 1 to 20 it will
        produce a plot of the R squared val versus the input hamming dist.
        """

        def swap(input_list, first_entry, second_entry):
            input_list[first_entry], input_list[second_entry] = input_list[second_entry], input_list[first_entry]
            return input_list

        def random_binary(binary):
            for index in range(hamming_dist):
                swap(binary, index, np.random.randint(self.n_qubo))
            return binary

        r_val = []
        hamming_weight = []
        for hamming_dist in range(self.n_problem - 1):

            hamming_dist = hamming_dist + 1
            hd = []
            x = []
            y = []

            for i in range(self.num_points):
                hd.append(np.zeros(self.n_qubo))

            for i in range(len(hd)):
                for j in range(hamming_dist):
                    hd[i][j] = 1
                random_binary(hd[i])

            for vec in hd:
                p_new = NextPerm(current_perm=self.p, lqubo_result=vec).next_perm()
                v_new = self.objective_function(p_new)
                x.append(v_new - self.v)
                y.append(np.matmul(np.matmul(vec, self.qubo), np.transpose(vec)))

            slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
            r_val.append(r_value**2)
            hamming_weight.append(hamming_dist)

        plt.plot(hamming_weight, r_val, label='n = {} {}'.format(self.n_problem, self.OF_type))

        plt.scatter(hamming_weight, r_val)
        plt.xticks(hamming_weight)
        plt.xlabel("Hamming dist")
        plt.ylabel("R Squared")
        plt.title('Hamming dist of LQUBO vs R Squared of scatter plot')
        plt.legend(loc='best')
        ax.xaxis.set_major_locator(MaxNLocator(integer=True))
        plt.show()

max_stars_repo_name: seangholson/lqubo

max_issues_repo_name: seangholson/lqubo

max_forks_repo_name: seangholson/QAP-Quantum-Computing

reason: import numpy,from scipy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Date    : 2015-12-19 02:09:53
# @Author  : Gefu Tang (tanggefu@gmail.com)
# @Link    : https://github.com/primetang/pylsd
# @Version : 0.0.1

import ctypes
import os
import sys
import random
import numpy as np

def load_lsd_library():

    root_dir = os.path.abspath(os.path.dirname(__file__))

    libnames = ['linux/liblsd.so']
    libdir = 'lib'
    if sys.platform == 'win32':
        if sys.maxsize > 2 ** 32:
            libnames = ['win32/x64/lsd.dll', 'win32/x64/liblsd.dll']
        else:
            libnames = ['win32/x86/lsd.dll', 'win32/x86/liblsd.dll']
    elif sys.platform == 'darwin':
        libnames = ['darwin/liblsd.dylib']

    while root_dir != None:
        for libname in libnames:
            try:
                lsdlib = ctypes.cdll[os.path.join(root_dir, libdir, libname)]
                return lsdlib
            except Exception as e:
                pass
        tmp = os.path.dirname(root_dir)
        if tmp == root_dir:
            root_dir = None
        else:
            root_dir = tmp

    # if we didn't find the library so far, try loading without
    # a full path as a last resort
    for libname in libnames:
        try:
            # print "Trying",libname
            lsdlib = ctypes.cdll[libname]
            return lsdlib
        except:
            pass

    return None

lsdlib = load_lsd_library()
if lsdlib == None:
    raise ImportError('Cannot load dynamic library. Did you compile LSD?')

max_stars_repo_name: TEAM-QUANDO/pylsd

max_issues_repo_name: TEAM-QUANDO/pylsd

max_forks_repo_name: TEAM-QUANDO/pylsd


text:
import numpy as np
import lqr_control as control

A = np.array([[1.01,0.01,0],[0.01,1.01,0.01],[0,0.01,1.01]])
B = np.identity(3)
Q = np.identity(3)
R1 = np.identity(3)*100
R2 = np.identity(3)*1000

x0 = np.array([[1],[0],[-1]]) #(0,0,0) is stable
u0 = np.zeros((3,1)) # default to 0 init for starting action

# number of time steps to simulate
T = 100

K_1, _, _ = control.dlqr(A,B,Q,R1)

K_2, _, _ = control.dlqr(A,B,Q,R2)

x_1, u_1 = control.simulate_discrete(A,B,K_1,x0,u0,T)
x_2, u_2 = control.simulate_discrete(A,B,K_2,x0,u0,T)

control.plot_states(x_1, 'State Temps', R1)
control.plot_states(x_2, 'State Temps', R2)


text:
'''
Encapsulates mesh for grasping operations
Author: Jeff Mahler
'''
import copy
import logging
import numpy as np
import os
import sys

class Mesh3D(object):
    """
    A Mesh is a three-dimensional shape representation

    Params:
       vertices:  (list of 3-lists of float)
       triangles: (list of 3-lists of ints)
       normals:   (list of 3-lists of float)
       metadata:  (dictionary) data like category, etc
       pose:      (tfx pose)
       scale:     (float)
       component: (int)
    """
    def __init__(self, vertices, triangles, normals):
        self.vertices_ = vertices
        self.triangles_ = triangles
        self.normals_ = normals

    @property
    def vertices(self):
        return self.vertices_

    @property
    def triangles(self):
        return self.triangles_

    @property
    def normals(self):
        return self.normals_

    @property
    def colors(self):
        return

    def compute_normals(self):
        """ Get vertex normals from triangles"""
        vertex_array = np.array(self.vertices_)
        tri_array = np.array(self.triangles_)
        self.normals_ = []
        for i in range(len(self.vertices_)):
            inds = np.where(tri_array == i)
            first_tri = tri_array[inds[0][0],:]
            t = vertex_array[first_tri, :]
            v0 = t[1,:] - t[0,:]
            v1 = t[2,:] - t[0,:]
            v0 = v0 / np.linalg.norm(v0)
            v1 = v1 / np.linalg.norm(v1)
            n = np.cross(v0, v1)
            n = n / np.linalg.norm(n)
            self.normals_.append(n.tolist())

    def tri_centers(self):
        """ Return a list of the triangle centers as 3D points """
        centers = []
        for tri in self.triangles_:
            v = np.array([self.vertices_[tri[0]], self.vertices_[tri[1]], self.vertices_[tri[2]]])
            center = np.mean(v, axis=0)
            centers.append(center.tolist())
        return centers

    def tri_normals(self):
        """ Return a list of the triangle normals """
        vertex_array = np.array(self.vertices_)
        self.normals_ = []
        for tri in self.triangles_:
            v0 = vertex_array[tri[0],:]
            v1 = vertex_array[tri[1],:]
            v2 = vertex_array[tri[2],:]
            n = np.cross(v1 - v0, v2 - v0)
            n = n / np.linalg.norm(n)
            self.normals_.append(n.tolist())

    def center_vertices_avg(self):
        """ Re-center vertices at average vertex """
        vertex_array = np.array(self.vertices_)
        centroid = np.mean(vertex_array, axis = 0)
        vertex_array_cent = vertex_array - centroid
        self.vertices_ = vertex_array_cent.tolist()

    def center_vertices_bb(self):
        """ Re-center vertices at center of bounding box """
        vertex_array = np.array(self.vertices_)
        min_vertex = np.min(vertex_array, axis = 0)
        max_vertex = np.max(vertex_array, axis = 0)
        centroid = (max_vertex + min_vertex) / 2
        vertex_array_cent = vertex_array - centroid
        self.vertices_ = vertex_array_cent.tolist()

    def rescale(self, scale_factor):
        """
        Rescales the vertex coordinates by scale factor
        Params:
           scale_factor: (float) the scale factor
        Returns:
           Nothing. Modified the mesh in place (for now)
        """
        vertex_array = np.array(self.vertices_)
        scaled_vertices = scale_factor * vertex_array
        self.vertices_ = scaled_vertices.tolist()


text:
# Copyright 2019 Jetperch LLC
#
# 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 logging
import numpy as np
import pyqtgraph as pg
from PySide2 import QtWidgets
from .histogram_config_widget import Ui_Dialog
from . import plugin_helpers

log = logging.getLogger(__name__)

PLUGIN = {
    'name': 'Histogram',
    'description': 'Histogram for current/voltage data',
}

_NORMALIZATIONS = {
    'Frequency Distribution': ('count', 'Sample Count'),
    'Discrete Probability Distribution': ('unity', 'Probability'),
    'Probability Density Distribution': ('density', 'Probability Density'),
}

class Histogram:

    def __init__(self):
        self._cfg = None
        self.hist = np.asarray([])
        self.bin_edges = np.asarray([])
        self._win = None
        self._label = None
        self._bg = None

    def run_pre(self, data):
        rv = HistogramDialog().exec_()
        if rv is None:
            return 'Cancelled'
        self._cfg = rv

    def run(self, data):
        norm = _NORMALIZATIONS[self._cfg['norm']][0]
        signal = self._cfg['signal']
        num_bins = self._cfg['num_bins']

        hist, bin_edges = plugin_helpers.calculate_histogram(data, num_bins, signal)
        self.hist, _bin_edges = plugin_helpers.normalize_hist(hist, bin_edges, norm)
        self.bin_edges = _bin_edges[:-1]

    def run_post(self, data):
        if self.hist.size == 0 or self.bin_edges.size == 0:
            log.exception('Histogram is empty')
            return

        self._win = pg.GraphicsLayoutWidget(show=True, title='Histogram')
        p = self._win.addPlot(row=1, col=0)
        p.getAxis('left').setGrid(128)
        p.getAxis('bottom').setGrid(128)

        self._label = pg.LabelItem(justify='right')
        width = self.bin_edges[1] - self.bin_edges[0]
        brushes = [(128, 128, 128)] * len(self.bin_edges)
        self._bg = pg.BarGraphItem(
            x0=self.bin_edges, height=self.hist, width=width, brushes=brushes)
        self.prev_hover_index = 0

        p.addItem(self._bg)
        self._win.addItem(self._label, row=0, col=0)

        left_axis_label = _NORMALIZATIONS[self._cfg['norm']][1]
        p.setLabels(left=left_axis_label, bottom=self._cfg['signal'])
        p.setXRange(self.bin_edges[0], self.bin_edges[-1], padding=0.05)
        p.setYRange(np.nanmin(self.hist), np.nanmax(self.hist), padding=0.05)

        def mouseMoved(evt):
            pos = evt[0]
            if p.sceneBoundingRect().contains(pos):
                mousePoint = p.vb.mapSceneToView(pos)
                xval = mousePoint.x()
                index = np.searchsorted(self.bin_edges, xval) - 1
                if index >= 0 and index < len(self.bin_edges):
                    self._label.setText(
                        "{}={:.5f},   {}: {:.5f}".format(
                            self._cfg['signal'], mousePoint.x(), left_axis_label, self.hist[index])
                    )
                    self._bg.opts['brushes'][self.prev_hover_index] = (128, 128, 128)
                    self._bg.opts['brushes'][index] = (213, 224, 61)
                    self.prev_hover_index = index
                    self._bg.drawPicture()

        self.proxy = pg.SignalProxy(
            p.scene().sigMouseMoved, rateLimit=60, slot=mouseMoved)

        self._win.closeEvent = lambda evt: data.on_tool_finished()
        return True

class HistogramDialog(QtWidgets.QDialog):
    def __init__(self):
        QtWidgets.QDialog.__init__(self)
        self.ui = Ui_Dialog()
        self.ui.setupUi(self)
        self.ui.signal.addItem("current")
        self.ui.signal.addItem("voltage")
        self.ui.signal.addItem("power")
        for name in _NORMALIZATIONS.keys():
            self.ui.normalization.addItem(name)

    def exec_(self):
        if QtWidgets.QDialog.exec_(self) == 1:
            num_bins = int(self.ui.num_bins.value())
            signal = str(self.ui.signal.currentText())
            norm = str(self.ui.normalization.currentText())
            return {
                'num_bins': num_bins,
                'signal': signal,
                'norm': norm
            }
        else:
            return None

def plugin_register(api):
    """Register the example plugin.

    :param api: The :class:`PluginServiceAPI` instance.
    :return: True on success any other value on failure.
    """
    api.range_tool_register('Analysis/Histogram', Histogram)
    return True


text:
import util
import util.cast
import numpy as np
import ctypes

def get_particle_hsml( x, y, z, DesNgb=32, Hmax=0. , **kwargs ):
    x=util.cast.fcor(x); y=util.cast.fcor(y); z=util.cast.fcor(z); N=util.cast.checklen(x);
    ok=(util.cast.ok_scan(x) & util.cast.ok_scan(y) & util.cast.ok_scan(z)); x=x[ok]; y=y[ok]; z=z[ok];
    if(Hmax==0.):
        dx=np.max(x)-np.min(x); dy=np.max(y)-np.min(y); dz=np.max(z)-np.min(z); ddx=np.max([dx,dy,dz]);
        Hmax=5.*ddx*(np.float(N)**(-1./3.)); ## mean inter-particle spacing

    ## load the routine we need
    exec_call=util.dir.c_routines_dir()+'/StellarHsml/starhsml.so'
    h_routine=ctypes.cdll[exec_call];

    h_out_cast=ctypes.c_float*N; H_OUT=h_out_cast();
    ## main call to the hsml-finding routine
    h_routine.stellarhsml( ctypes.c_int(N),
                          util.cast.vfloat(x),
                          util.cast.vfloat(y),
                          util.cast.vfloat(z),
                          ctypes.c_int(DesNgb),
                          ctypes.c_float(Hmax),
                          ctypes.byref(H_OUT) )

    ## now put the output arrays into a useful format
    h = np.ctypeslib.as_array(H_OUT);
    return h;

def get_gas_density_around_stars( x_gas, y_gas, z_gas, m_gas, x_star, y_star, z_star, DesNgb=32, Hmax=0. , **kwargs ):

    x_gas=util.cast.fcor(x_gas); y_gas=util.cast.fcor(y_gas); z_gas=util.cast.fcor(z_gas);
    m_gas=util.cast.fcor(m_gas);

    x_star=util.cast.fcor(x_star); y_star=util.cast.fcor(y_star); z_star=util.cast.fcor(z_star);
    N_star=util.cast.checklen(x_star);

    ok=(util.cast.ok_scan(x_gas) & util.cast.ok_scan(y_gas) & util.cast.ok_scan(z_gas) & util.cast.ok_scan(m_gas));
    x_gas=x_gas[ok]; y_gas=y_gas[ok]; z_gas=z_gas[ok]; m_gas=m_gas[ok];
    N_gas=util.cast.checklen(x_gas);

    if(Hmax==0.):
        dx=np.max(x_gas)-np.min(x_gas); dy=np.max(y_gas)-np.min(y_gas); dz=np.max(z_gas)-np.min(z_gas); ddx=np.max([dx,dy,dz]);
        Hmax=5.*ddx*(np.float(N_gas)**(-1./3.)); ## mean inter-particle spacing

    ## load the routine we need
    exec_call=util.dir.c_routines_dir()+'/StellarGasDensity/stargasdensity.so'
    h_routine=ctypes.cdll[exec_call];

    h_out_cast=ctypes.c_float*N_star; H_OUT=h_out_cast();
    ## main call to the hsml-finding routine
    h_routine.stellargasdensity( ctypes.c_int(N_gas), ctypes.c_int(N_star),
                          util.cast.vfloat(x_gas),  util.cast.vfloat(y_gas),    util.cast.vfloat(z_gas), util.cast.vfloat(m_gas),
                          util.cast.vfloat(x_star),  util.cast.vfloat(y_star),    util.cast.vfloat(z_star),
                          ctypes.c_int(DesNgb),
                          ctypes.c_float(Hmax),
                          ctypes.byref(H_OUT) )

                          ## now put the output arrays into a useful format
    h = np.ctypeslib.as_array(H_OUT);
    return h;

def get_gas_temperature_around_stars( x_gas, y_gas, z_gas, m_gas, rho_gas, u_gas, x_star, y_star, z_star, DesNgb=32, Hmax=0. , **kwargs ):

    x_gas=util.cast.fcor(x_gas); y_gas=util.cast.fcor(y_gas); z_gas=util.cast.fcor(z_gas);
    m_gas=util.cast.fcor(m_gas);
    rho_gas=util.cast.fcor(rho_gas);
    u_gas=util.cast.fcor(u_gas);

    x_star=util.cast.fcor(x_star); y_star=util.cast.fcor(y_star); z_star=util.cast.fcor(z_star);
    N_star=util.cast.checklen(x_star);

    ok=(util.cast.ok_scan(x_gas) & util.cast.ok_scan(y_gas) & util.cast.ok_scan(z_gas) & util.cast.ok_scan(m_gas));
    x_gas=x_gas[ok]; y_gas=y_gas[ok]; z_gas=z_gas[ok]; m_gas=m_gas[ok]; rho_gas=rho_gas[ok]; u_gas=u_gas[ok]
    N_gas=util.cast.checklen(x_gas);

    if(Hmax==0.):
        dx=np.max(x_gas)-np.min(x_gas); dy=np.max(y_gas)-np.min(y_gas); dz=np.max(z_gas)-np.min(z_gas); ddx=np.max([dx,dy,dz]);
        Hmax=5.*ddx*(np.float(N_gas)**(-1./3.)); ## mean inter-particle spacing

    ## load the routine we need
    exec_call=util.dir.c_routines_dir()+'/StellarGasTemperature/stargasavg.so'
    h_routine=ctypes.cdll[exec_call];

    h_out_cast=ctypes.c_float*N_star; H_OUT=h_out_cast();
    ## main call to the hsml-finding routine
    h_routine.stellargastemperature( ctypes.c_int(N_gas), ctypes.c_int(N_star),
                          util.cast.vfloat(x_gas),  util.cast.vfloat(y_gas),    util.cast.vfloat(z_gas), util.cast.vfloat(m_gas), util.cast.vfloat(rho_gas), util.cast.vfloat(u_gas),
                          util.cast.vfloat(x_star),  util.cast.vfloat(y_star),    util.cast.vfloat(z_star),
                          ctypes.c_int(DesNgb),
                          ctypes.c_float(Hmax),
                          ctypes.byref(H_OUT) )

                          ## now put the output arrays into a useful format
    h = np.ctypeslib.as_array(H_OUT);
    return h;


text:
"""
Created on April 15, 2018

@author: Alejandro Molina
@author: Antonio Vergari
"""
from spn.algorithms.Sampling import add_leaf_sampling
from spn.structure.leaves.histogram.Histograms import (
    Histogram
)

import numpy as np

from spn.structure.leaves.parametric.utils import get_scipy_obj_params
import logging

logger = logging.getLogger(__name__)

def sample_histogram_node(node, n_samples, data, rand_gen):
    assert isinstance(node, Histogram)
    assert n_samples > 0
    # sample the value at each bin according to the densities of each bin
    # print(len(node.bin_repr_points))
    X = rand_gen.choice(np.array(node.bin_repr_points), p=node.densities, size=n_samples)

    return X

def add_histogram_sampling_support():
    add_leaf_sampling(Histogram, sample_histogram_node)

max_stars_repo_name: minimrbanana/SPFlow

max_issues_repo_name: minimrbanana/SPFlow

max_forks_repo_name: minimrbanana/SPFlow