text:
import time
from contextlib import suppress, contextmanager
from astropy import units as u

from panoptes.utils import error
from panoptes.utils.utils import get_quantity_value
from panoptes.utils.time import current_time, wait_for_events, CountdownTimer

from panoptes.pocs.observatory import Observatory
from panoptes.pocs.scheduler.observation.bias import BiasObservation

from huntsman.pocs.utils.logger import get_logger
from huntsman.pocs.guide.bisque import Guide
from huntsman.pocs.archive.utils import remove_empty_directories
from huntsman.pocs.scheduler.observation.dark import DarkObservation
from huntsman.pocs.utils.flats import make_flat_field_sequences, make_flat_field_observation
from huntsman.pocs.utils.flats import get_cameras_with_filter
from huntsman.pocs.utils.safety import get_solar_altaz
from huntsman.pocs.camera.group import CameraGroup, dispatch_parallel
from huntsman.pocs.error import NotTwilightError

class HuntsmanObservatory(Observatory):

    def __init__(self, with_autoguider=True, hdr_mode=False, take_flats=True, logger=None,
                 *args, **kwargs):
        """Huntsman POCS Observatory

        Args:
            with_autoguider (bool, optional): If autoguider is attached, defaults to True.
            hdr_mode (bool, optional): If pics should be taken in HDR mode, defaults to False.
            take_flats (bool, optional): If flat field images should be taken, defaults to True.
            logger (logger, optional): The logger instance. If not provided, use default Huntsman
                logger.
            *args: Parsed to Observatory init function.
            **kwargs: Parsed to Observatory init function.
        """
        if not logger:
            logger = get_logger()
        super().__init__(logger=logger, *args, **kwargs)

        # Make a camera group
        self.camera_group = CameraGroup(self.cameras)

        self._has_hdr_mode = hdr_mode
        self._has_autoguider = with_autoguider

        self.flat_fields_required = take_flats

        # Focusing
        self.last_coarse_focus_time = None
        self.last_coarse_focus_temp = None
        self._coarse_focus_interval = self.get_config('focusing.coarse.interval_hours', 1) * u.hour
        self._coarse_focus_filter = self.get_config('focusing.coarse.filter_name')
        self._coarse_focus_temptol = self.get_config('focusing.coarse.temp_tol_deg', 5) * u.Celsius

        self.last_fine_focus_time = None
        self.last_fine_focus_temp = None
        self._fine_focus_interval = self.get_config('focusing.fine.interval_hours', 1) * u.hour
        self._fine_focus_temptol = self.get_config('focusing.fine.temp_tol_deg', 5) * u.Celsius

        if self.has_autoguider:
            self.logger.info("Setting up autoguider")
            try:
                self._create_autoguider()
            except Exception as e:
                self._has_autoguider = False
                self.logger.warning(f"Problem setting autoguider, continuing without: {e!r}")

        # Hack solution to the observatory not knowing whether it is safe or not
        # This can be overridden when creating the HuntsmanPOCS instance
        self._is_safe = None

    # Properties

    @property
    def has_hdr_mode(self):
        """ Does camera support HDR mode
        Returns:
            bool: HDR enabled, default False
        """
        return self._has_hdr_mode

    @property
    def has_autoguider(self):
        """ Does camera have attached autoguider
        Returns:
            bool: True if has autoguider
        """
        return self._has_autoguider

    @property
    def coarse_focus_required(self):
        """ Return True if we should do a coarse focus. """
        return self._focus_required(coarse=True)

    @property
    def fine_focus_required(self):
        """ Return True if we should do a fine focus. """
        return self._focus_required()

    @property
    def is_past_midnight(self):
        """Check if it's morning, useful for going into either morning or evening flats."""
        # Get the time of the nearest midnight to now
        # If the nearest midnight is in the past, it's the morning
        midnight = self.observer.midnight(current_time(), which='nearest')
        return midnight < current_time()

    @property
    def is_twilight(self):
        """ Return True if it is twilight, else False. """
        return self.is_dark(horizon="twilight_max") and not self.is_dark(horizon="twilight_min")

    @property
    def temperature(self):
        """ Return the ambient temperature. """
        temp = None
        try:
            reading = self.db.get_current("weather")["data"]["ambient_temp_C"]
            temp = get_quantity_value(reading, u.Celsius) * u.Celsius
        except (KeyError, TypeError) as err:
            self.logger.warning(f"Unable to determine temperature: {err!r}")
        return temp

    @property
    def solar_altaz(self):
        """ Return the current solar alt az. """
        return get_solar_altaz(location=self.earth_location, time=current_time())

    # Context managers

    @contextmanager
    def safety_checking(self, *args, **kwargs):
        """ Check safety before and after the code block.
        To be used with a "with" statement, e.g.:
            with self.safety_checking():
                print(x)
        Args:
            *args, **kwargs: Parsed to self._assert_safe
        Raises:
            RuntimeError: If not safe.
        """
        self._assert_safe(*args, **kwargs)
        try:
            yield None
        finally:
            self._assert_safe(*args, **kwargs)

    # Methods

    def initialize(self):
        """Initialize the observatory and connected hardware """
        super().initialize()

        if self.has_autoguider:
            self.logger.debug("Connecting to autoguider")
            self.autoguider.connect()

    def is_safe(self, park_if_not_safe=False, *args, **kwargs):
        """ Return True if it is safe, else False.
        Args:
            *args, **kwargs: Parsed to self._is_safe. See panoptes.pocs.core.POCS.is_safe.
            park_if_not_safe (bool): If True, park if safety fails. Default: False.
        Returns:
            bool: True if safe, else False.
        """
        if self._is_safe is not None:
            return self._is_safe(park_if_not_safe=park_if_not_safe, *args, **kwargs)
        self.logger.warning("Safety function not set. Returning False")
        return False

    def remove_camera(self, cam_name):
        """ Remove a camera from the observatory.
        Args:
            cam_name (str): The name of the camera to remove.
        """
        super().remove_camera(cam_name)
        with suppress(KeyError):
            del self.camera_group.cameras[cam_name]

    def autofocus_cameras(self, coarse=False, filter_name=None, default_timeout=900,
                          blocking=True, **kwargs):
        """ Override autofocus_cameras to update the last focus time and move filterwheels.
        Args:
            coarse (bool, optional): Perform coarse focus? Default False.
            filter_name (str, optional): The filter name to focus with. If None (default), will
                attempt to get from config, by default using the coarse focus filter.
            *args, **kwargs: Parsed to `pocs.observatory.Observatory.autofocus_cameras`.
        Returns:
            threading.Event: The autofocus event.
        """
        focus_type = "coarse" if coarse else "fine"

        # Choose the filter to focus with
        # TODO: Move this logic to the camera level
        if filter_name is None:
            if coarse:
                filter_name = self._coarse_focus_filter
            else:
                try:
                    filter_name = self.current_observation.filter_name
                except AttributeError:
                    filter_name = self._coarse_focus_filter
                    self.logger.warning("Unable to retrieve filter name from current observation."
                                        f" Defaulting to coarse focus filter ({filter_name}).")

        # Asyncronously dispatch autofocus calls
        with self.safety_checking(horizon="focus"):
            events = self.camera_group.autofocus(coarse=coarse, filter_name=filter_name, **kwargs)

        # Wait for sequences to finish
        if blocking:
            timeout = self.get_config(f"focusing.{focus_type}.timeout", default_timeout)
            if not wait_for_events(list(events.values()), timeout=timeout):
                raise error.Timeout(f"Timeout of {timeout} reached while waiting for fine focus.")

        # Update last focus time
        setattr(self, f"last_{focus_type}_focus_time", current_time())

        # Update last focus temperature
        setattr(self, f"last_{focus_type}_focus_temp", self.temperature)

        return events

    def cleanup_observations(self, *args, **kwargs):
        """ Override method to remove empty directories. Called in housekeeping state."""
        super().cleanup_observations(*args, **kwargs)

        self.logger.info("Removing empty directories in images directory.")
        images_dir = self.get_config("directories.images")
        remove_empty_directories(images_dir)

        self.logger.info("Removing empty directories in archive directory.")
        archive_dir = self.get_config("directories.archive")
        remove_empty_directories(archive_dir)

    def take_flat_fields(self, cameras=None, **kwargs):
        """ Take flat fields for each camera in each filter, respecting filter order.
        Args:
            cameras (dict): Dict of cam_name: camera pairs. If None (default), use all cameras.
            **kwargs: Overrides config entries under `calibs.flat`.
        """
        if cameras is None:
            cameras = self.cameras

        # Load the flat field config, allowing overrides from kwargs
        flat_config = self.get_config('calibs.flat', default=dict())
        flat_config.update(kwargs)

        # Specify filter order
        filter_order = flat_config['filter_order'].copy()
        if self.is_past_midnight:  # If it's the morning, order is reversed
            filter_order.reverse()

        # Take flat fields in each filter
        for filter_name in filter_order:

            if not (self.is_safe(horizon="twilight_max") and self.is_twilight):
                raise RuntimeError("Not safe for twilight flats. Aborting.")

            # Get a dict of cameras that have this filter
            cameras_with_filter = get_cameras_with_filter(cameras, filter_name)

            # Go to next filter if there are no cameras with this one
            if not cameras_with_filter:
                self.logger.warning(f'No cameras found with {filter_name} filter.')
                continue

            # Get the flat field observation
            observation = make_flat_field_observation(self.earth_location, filter_name=filter_name)
            observation.seq_time = current_time(flatten=True)

            # Take the flats for each camera in this filter
            self.logger.info(f'Taking flat fields in {filter_name} filter.')
            autoflat_config = flat_config.get("autoflats", {})

            try:
                self._take_autoflats(cameras_with_filter, observation, **autoflat_config)
            # Break out of loop if no longer twilight
            # Catch the error so the state machine keeps running
            except NotTwilightError as err:
                self.logger.warning(f"{err!r}")
                break

        self.logger.info('Finished flat-fielding.')

    def prepare_cameras(self, drop=True, *args, **kwargs):
        """ Make sure cameras are all cooled and ready.
        Args:
            drop (bool): If True, drop cameras that do not become ready in time. Default: True.
            *args, **kwargs: Parsed to self.camera_group.wait_until_ready.
        """
        self.logger.info(f"Preparing {len(self.cameras)} cameras.")

        failed_cameras = self.camera_group.wait_until_ready(*args, **kwargs)

        # Remove cameras that didn't become ready in time
        if drop:
            for cam_name in failed_cameras:
                self.logger.debug(f'Removing {cam_name} from {self} for not being ready.')
                self.remove_camera(cam_name)

    def take_observation_block(self, observation, cameras=None, timeout=60 * u.second,
                               remove_on_error=False, do_focus=True, safety_kwargs=None,
                               do_slew=True):
        """ Macro function to take an observation block.
        This function will perform:
            - slewing (when necessary)
            - fine focusing (when necessary)
            - observation exposures
            - safety checking
        Args:
            observation (Observation): The observation object.
            cameras (dict, optional): Dict of cam_name: camera pairs. If None (default), use all
                cameras.
            timeout (float, optional): The timeout in addition to the exposure time. Default 60s.
            remove_on_error (bool, default False): If True, remove cameras that timeout. If False,
                raise a TimeoutError instead.
            do_slew (bool, optional): If True, do not attempt to slew the telescope. Default
                False.
            **safety_kwargs (dict, optional): Extra kwargs to be parsed to safety function.
        Raises:
            RuntimeError: If safety check fails.
        """
        if cameras is None:
            cameras = self.cameras

        safety_kwargs = {} if safety_kwargs is None else safety_kwargs
        self._assert_safe(**safety_kwargs)

        # Set the sequence time of the observation
        if observation.seq_time is None:
            observation.seq_time = current_time(flatten=True)

        headers = self.get_standard_headers(observation=observation)

        # Take the observation block
        self.logger.info(f"Starting observation block for {observation}")

        # The start new set flag is True before we enter the loop and is set to False
        # immediately inside the loop. This allows the loop to start a new set in case
        # the set_is_finished condition is already satisfied.
        start_new_set = True

        current_field = None
        while (start_new_set or not observation.set_is_finished):

            start_new_set = False  # We don't want to start another set after this one

            # Perform the slew if necessary
            slew_required = (current_field != observation.field) and do_slew
            if slew_required:
                with self.safety_checking(**safety_kwargs):
                    self.slew_to_observation(observation)
                current_field = observation.field

            # Fine focus the cameras if necessary
            focus_required = self.fine_focus_required or observation.current_exp_num == 0
            if do_focus and focus_required:
                with self.safety_checking(**safety_kwargs):
                    self.autofocus_cameras(blocking=True, filter_name=observation.filter_name)

            # Set a common start time for this batch of exposures
            headers['start_time'] = current_time(flatten=True)

            # Start the exposures and get events
            with self.safety_checking(**safety_kwargs):
                events = self.camera_group.take_observation(observation, headers=headers)

            # Wait for the exposures (blocking)
            # TODO: Use same timeout as camera client
            try:
                self._wait_for_camera_events(events, duration=observation.exptime + timeout,
                                             remove_on_error=remove_on_error, **safety_kwargs)
            except error.Timeout as err:
                self.logger.error(f"{err!r}")
                self.logger.warning("Continuing with observation block after error.")

            # Explicitly mark the observation as complete
            with suppress(AttributeError):
                observation.mark_exposure_complete()

            self.logger.info(f"Observation status: {observation.status}")

    def take_dark_observation(self, bias=False, **kwargs):
        """ Take a bias observation block on each camera (blocking).
        Args:
            bias (bool, optional): If True, take Bias observation instead of dark observation.
                Default: False.
            **kwargs: Parsed to `self.take_observation_block`.
        """
        # Move telescope to park position
        if not self.mount.is_parked:
            self.logger.info("Moving telescope to park position for dark observation.")
            self.mount.park()

        # Create the observation
        # Keep the mount where it is since we are just taking darks
        position = self.mount.get_current_coordinates()
        ObsClass = BiasObservation if bias else DarkObservation
        observation = ObsClass(position=position)

        # Dark observations don't care if it's dark or not
        safety_kwargs = {"ignore": ["is_dark"]}

        # Can ignore weather safety if dome is closed
        with suppress(AttributeError):
            if self.dome.is_closed:
                self.logger.warning(f"Ignoring weather safety for {observation}.")
                safety_kwargs["ignore"].append("good_weather")

        # Take the observation (blocking)
        self.take_observation_block(observation, do_focus=False, do_slew=False,
                                    safety_kwargs=safety_kwargs, **kwargs)

    def slew_to_observation(self, observation, min_solar_alt=10 * u.deg):
        """ Slew to the observation field coordinates.
        Args:
            observation (Observation): The observation object.
            min_solar_alt (astropy.Quantity, optional): The minimum solar altitude above which the
                FWs will be moved to their dark positions before slewing.
        """
        self.logger.info(f"Slewing to target coordinates for {observation}.")

        if not self.mount.set_target_coordinates(observation.field.coord):
            raise RuntimeError(f"Unable to set target coordinates for {observation.field}.")

        # Move FWs to dark pos if Sun too high to minimise damage potential
        move_fws = self.solar_altaz.alt > get_quantity_value(min_solar_alt, u.deg) * u.deg

        if move_fws:
            self.logger.warning("Solar altitude above minimum for safe slew. Moving FWs to dark"
                                " positions.")

            # Record curent positions so we can put them back after slew
            # NOTE: These positions could include the dark position so can't use last_light_position
            current_fw_positions = {}
            for cam_name, cam in self.cameras.items():
                if cam.has_filterwheel:
                    current_fw_positions[cam_name] = cam.filterwheel.current_filter

            self.camera_group.filterwheel_move_to(current_fw_positions)

        self.mount.slew_to_target()

        if move_fws:
            self.logger.info("Moving FWs back to last positions.")
            self.camera_group.filterwheel_move_to(current_fw_positions)

    # Private methods

    def _create_autoguider(self):
        guider_config = self.get_config('guider')
        guider = Guide(**guider_config)
        self.autoguider = guider

    def _take_autoflats(
            self, cameras, observation, target_scaling=0.17, scaling_tolerance=0.05, timeout=60,
            bias=32, remove_on_error=False, sleep_time=300, evening_initial_flat_exptime=0.01,
            morning_initial_flat_exptime=1, **kwargs):
        """ Take flat fields using automatic updates for exposure times.
        Args:
            cameras (dict): Dict of camera name: Camera pairs.
            observation: The flat field observation. TODO: Integrate with FlatFieldSequence.
            target_scaling (float, optional): Required to be between [0, 1] so
                target_adu is proportionally between 0 and digital saturation level.
                Default: 0.17.
            scaling_tolerance (float, optional): The minimum precision on the average counts
                required to keep the exposure, expressed as a fraction of the dynamic range.
                Default: 0.05.
            timeout (float): The timeout on top of the exposure time, default 60s.
            bias (int): The bias to subtract from the frames. TODO: Use a real bias image!
            remove_on_error (bool, default False): If True, remove cameras that timeout. If False,
                raise a TimeoutError instead.
            **kwargs: Parsed to FlatFieldSequence.
        """
        # set the initial exposure time
        if self.is_past_midnight:
            initial_exptime = morning_initial_flat_exptime
        else:
            initial_exptime = evening_initial_flat_exptime

        # Create a flat field sequence for each camera
        sequences = make_flat_field_sequences(cameras, target_scaling, scaling_tolerance,
                                              bias, initial_exposure_time=initial_exptime, **kwargs)

        # Loop until sequence has finished
        self.logger.info(f"Starting flat field sequence for {len(self.cameras)} cameras.")
        while True:

            if not self.is_twilight:
                raise NotTwilightError("No longer twilight. Aborting flat fields.")

            # Slew to field
            with self.safety_checking(horizon="twilight_max"):
                self.slew_to_observation(observation)

            # Get standard fits headers
            headers = self.get_standard_headers(observation=observation)

            events = {}
            exptimes = {}
            filenames = {}
            start_times = {}

            # Define function to start the exposures
            def func(cam_name):
                seq = sequences[cam_name]
                camera = cameras[cam_name]

                # Get exposure time, filename and current time
                exptimes[cam_name] = seq.get_next_exptime(past_midnight=self.is_past_midnight)
                filenames[cam_name] = observation.get_exposure_filename(camera)
                start_times[cam_name] = current_time()

                try:
                    events[cam_name] = camera.take_observation(
                        observation, headers=headers, filename=filenames[cam_name],
                        exptime=exptimes[cam_name])
                except error.PanError as err:
                    self.logger.error(f"{err!r}")
                    self.logger.warning("Continuing with flat observation after error.")

            # Start the exposures in parallel
            dispatch_parallel(func, list(cameras.keys()))

            # Wait for the exposures
            self.logger.info('Waiting for flat field exposures to complete.')
            duration = get_quantity_value(max(exptimes.values()), u.second) + timeout
            try:
                self._wait_for_camera_events(events, duration, remove_on_error=remove_on_error,
                                             horizon="twilight_max")
            except error.Timeout as err:
                self.logger.error(f"{err!r}")
                self.logger.warning("Continuing with flat observation after timeout error.")

            # Mark the current exposure as complete
            observation.mark_exposure_complete()

            # Update the flat field sequences with new data
            for cam_name in list(sequences.keys()):

                # Remove sequence for any removed cameras
                if cam_name not in self.cameras:
                    del sequences[cam_name]
                    continue

                # Attempt to update the exposure sequence for this camera.
                # If the exposure failed, use info from the last successful exposure.
                try:
                    sequences[cam_name].update(filename=filenames[cam_name],
                                               exptime=exptimes[cam_name],
                                               time_start=start_times[cam_name])
                except (KeyError, FileNotFoundError) as err:
                    self.logger.warning(f"Unable to update flat field sequence for {cam_name}:"
                                        f" {err!r}")
                # Log sequence status
                status = sequences[cam_name].status
                status["filter_name"] = observation.filter_name
                self.logger.info(f"Flat field status for {cam_name}: {status}")

            # Check if sequences are complete
            if all([s.is_finished for s in sequences.values()]):
                self.logger.info("All flat field sequences finished.")
                break

            # Check if counts are ok
            if self.is_past_midnight:

                # Terminate if Sun is coming up and all exposures are too bright
                if all([s.min_exptime_reached for s in sequences.values()]):
                    self.logger.info(f"Terminating flat sequence for {observation.filter_name}"
                                     f" filter because min exposure time reached.")
                    break

                # Wait if Sun is coming up and all exposures are too faint
                elif all([s.max_exptime_reached for s in sequences.values()]):
                    self.logger.info(f"All exposures are too faint. Waiting for {sleep_time}s")
                    self._safe_sleep(sleep_time, horizon="twilight_max")

            else:
                # Terminate if Sun is going down and all exposures are too faint
                if all([s.max_exptime_reached for s in sequences.values()]):
                    self.logger.info(f"Terminating flat sequence for {observation.filter_name}"
                                     f" filter because max exposure time reached.")
                    break

                # Wait if Sun is going down and all exposures are too bright
                elif all([s.max_exptime_reached for s in sequences.values()]):
                    self.logger.info(f"All exposures are too bright. Waiting for {sleep_time}s")
                    self._safe_sleep(sleep_time, horizon="twilight_max")

    def _wait_for_camera_events(self, events, duration, remove_on_error=False, sleep=1, **kwargs):
        """ Wait for camera events to be set.
        Args:
            events (dict of camera_name: threading.Event): The events to wait for.
            duration (float): The total amount of time to wait for (should include exptime).
            remove_on_error (bool, default False): If True, remove cameras that timeout. If False,
                raise a TimeoutError instead.
            sleep (float): Sleep this long between event checks. Default 1s.
            **kwargs: Parsed to self._assert_safe.
        """
        self.logger.debug(f'Waiting for {len(events)} events with timeout of {duration}.')

        timer = CountdownTimer(duration)
        while not timer.expired():

            # Check safety here
            self._assert_safe(**kwargs)

            # Check if all cameras have finished
            if all([e.is_set() for e in events.values()]):
                break
            time.sleep(sleep)

        # Make sure events are set
        for cam_name, event in events.items():
            if not event.is_set():
                if remove_on_error:
                    self.logger.warning(f"Timeout while waiting for camera event on {cam_name}. "
                                        "Removing from observatory.")
                    self.remove_camera(cam_name)
                else:
                    raise error.Timeout(f"Timeout while waiting for camera event on {cam_name}.")

    def _focus_required(self, coarse=False):
        """ Check if a focus is required based on current conditions.
        Args:
            coarse (bool): If True, check if we need to do a coarse focus. Default: False.
        Returns:
            bool: True if focus required, else False.
        """
        focus_type = "coarse" if coarse else "fine"

        # If a long time period has passed then focus again
        last_focus_time = getattr(self, f"last_{focus_type}_focus_time")
        interval = getattr(self, f"_{focus_type}_focus_interval")

        if last_focus_time is None:  # If we haven't focused yet
            self.logger.info(f"{focus_type} focus required because we haven't focused yet.")
            return True

        if current_time() - last_focus_time > interval:
            self.logger.info(f"{focus_type} focus required because of time difference.")
            return True

        # If there has been a large change in temperature then we need to focus again
        last_focus_temp = getattr(self, f"last_{focus_type}_focus_temp")
        temptol = getattr(self, f"_{focus_type}_focus_temptol")

        if (last_focus_temp is not None) and (self.temperature is not None):
            if abs(last_focus_temp - self.temperature) > temptol:
                self.logger.info(f"{focus_type} focus required because of temperature change.")
                return True

        return False

    def _assert_safe(self, *args, **kwargs):
        """ Raise a RuntimeError if not safe to continue.
        TODO: Raise a custom error type indicating lack of safety.
        Args:
            *args, **kwargs: Parsed to self.is_safe.
        """
        if not self.is_safe(*args, **kwargs):
            raise RuntimeError("Safety check failed!")

    def _safe_sleep(self, duration, interval=1, *args, **kwargs):
        """ Sleep for a specified amount of time while ensuring safety.
        A RuntimeError is raised if safety fails while waiting.
        Args:
            duration (float or Quantity): The time to wait.
            interval (float): The time in between safety checks.
            *args, **kwargs: Parsed to is_safe.
        Raises:
            RuntimeError: If safety fails while waiting.
        """
        self.logger.debug(f"Safe sleeping for {duration}")
        timer = CountdownTimer(duration)
        while not timer.expired():
            self._assert_safe(*args, **kwargs)
            time.sleep(interval)

max_stars_repo_name: Physarah/huntsman-pocs

max_issues_repo_name: Physarah/huntsman-pocs

max_forks_repo_name: Physarah/huntsman-pocs

reason: from astropy


text:
#from setuptools import setup
from distutils.core import setup, Extension
from Cython.Build import cythonize
from Cython.Distutils import build_ext
from os.path import join, dirname
import numpy

src = ["io_c.pyx","io.c"]
src = [join('PySVF', x) for x in src]

dep = ["allvars.h","proto.h"]
dep = [join('PySVF', x) for x in dep]

#setup(cmdclass={'build_ext':build_ext}, ext_modules=[Extension('io_c', sources=src,depends=dep,language="c++",include_dirs=[numpy.get_include()])])

setup(ext_modules = cythonize('PySVF/*.pyx',language="c++"))

#setup(ext_modules = cythonize('Examples/*.pyx'))#Examples:"hello_word.c"y"rectangle.c"

ext: py

max_stars_repo_name: geralfaro/PySVF

max_issues_repo_name: geralfaro/PySVF

max_forks_repo_name: geralfaro/PySVF

reason: import numpy


text:
import time
import os
from minio import Minio
from minio.error import S3Error
import numpy as np
import random
from sklearn import preprocessing
from common.common_func import gen_unique_str

minio = "10.96.1.23:9000"  # TODO update hardcode
bucket_name = "yanliang-bulk-load"  # TODO update hardcode

data_source = "/tmp/bulk_load_data"

BINARY = "binary"
FLOAT = "float"

def gen_file_prefix(row_based=True, auto_id=True, prefix=""):
    if row_based:
        if auto_id:
            return f"{prefix}row_auto"
        else:
            return f"{prefix}row_cust"
    else:
        if auto_id:
            return f"{prefix}col_auto"
        else:
            return f"{prefix}col_cust"

def entity_suffix(rows):
    if rows // 1000000 > 0:
        suffix = f"{rows // 1000000}m"
    elif rows // 1000 > 0:
        suffix = f"{rows // 1000}k"
    else:
        suffix = f"{rows}"
    return suffix

def gen_float_vectors(nb, dim):
    vectors = [[random.random() for _ in range(dim)] for _ in range(nb)]
    vectors = preprocessing.normalize(vectors, axis=1, norm='l2')
    return vectors.tolist()

def gen_binary_vectors(nb, dim):
    # binary: each int presents 8 dimension
    # so if binary vector dimension is 16，use [x, y], which x and y could be any int between 0 to 255
    vectors = [[random.randint(0, 255) for _ in range(dim)] for _ in range(nb)]
    # vectors = preprocessing.normalize(vectors, axis=1, norm='l2')
    return vectors

def gen_row_based_json_file(row_file, str_pk, multi_scalars, float_vect, rows, dim, autoid):
    with open(row_file, "w") as f:
        f.write("{")
        f.write("\n")
        f.write('"rows":[')
        f.write("\n")
        for i in range(rows):
            if i > 0:
                f.write(",")
                f.write("\n")
            # pk fields
            if not autoid:
                if str_pk:
                    f.write('{"uid":"' + str(gen_unique_str()) + '",')
                else:
                    f.write('{"uid":' + str(i) + ',')
            else:
                f.write('{')

            # scalar fields
            if multi_scalars:
                f.write('"int_scalar":' + str(random.randint(-999999, 9999999)) + ',')
                f.write('"float_scalar":' + str(random.random()) + ',')
                f.write('"string_scalar":"' + str(gen_unique_str()) + '",')
                f.write('"bool_scalar":' + str(random.choice(["true", "false"])) + ',')

            # vector field
            if float_vect:
                vectors = gen_float_vectors(1, dim)
            else:
                vectors = gen_binary_vectors(1, (dim//8))
            f.write('"vectors":' + ",".join(str(x) for x in vectors) + "}")
        f.write("\n")
        f.write("]")
        f.write("\n")
        f.write("}")
        f.write("\n")

def gen_column_base_json_file(col_file, str_pk, float_vect, multi_scalars, rows, dim, autoid):
    with open(col_file, "w") as f:
        f.write("{")
        f.write("\n")
        # pk columns
        if not autoid:
            if str_pk == "str_pk":
                f.write('"uid":["' + ',"'.join(str(gen_unique_str()) + '"' for i in range(rows)) + '],')
                f.write("\n")
            else:
                f.write('"uid":[' + ",".join(str(i) for i in range(rows)) + "],")
                f.write("\n")

        # scalar columns
        if multi_scalars:
            f.write('"int_scalar":[' + ",".join(str(random.randint(-999999, 9999999)) for i in range(rows)) + "],")
            f.write("\n")
            f.write('"float_scalar":[' + ",".join(str(random.random()) for i in range(rows)) + "],")
            f.write("\n")
            f.write('"string_scalar":["' + ',"'.join(str(gen_unique_str()) + '"' for i in range(rows)) + '],')
            f.write("\n")
            f.write('"bool_scalar":[' + ",".join(str(random.choice(["true", "false"])) for i in range(rows)) + "],")
            f.write("\n")

        # vector columns
        if float_vect:
            vectors = gen_float_vectors(rows, dim)
            f.write('"vectors":[' + ",".join(str(x) for x in vectors) + "]")
            f.write("\n")
        else:
            vectors = gen_binary_vectors(rows, (dim//8))
            f.write('"vectors":[' + ",".join(str(x) for x in vectors) + "]")
            f.write("\n")

        f.write("}")
        f.write("\n")

def gen_vectors_in_numpy_file(dir, vector_type, rows, dim, num):
    # vector columns
    if vector_type == FLOAT:
        vectors = gen_float_vectors(rows, dim)
    else:
        vectors = gen_binary_vectors(rows, (dim // 8))

    suffix = entity_suffix(rows)
    # print(vectors)
    arr = np.array(vectors)
    path = f"{dir}/{dim}d_{suffix}_{num}"
    if not os.path.isdir(path):
        os.mkdir(path)
    file = f'{path}/vectors_{dim}d_{suffix}.npy'
    np.save(file, arr)

def gen_scalars_in_numpy_file(dir, vector_type, rows, dim, num, start):
    # scalar columns
    if vector_type == FLOAT:
        data = [random.random() for i in range(rows)]
    elif vector_type == "int":
        data = [i for i in range(start, start + rows)]

    suffix = entity_suffix(rows)
    path = f"{dir}/{dim}d_{suffix}_{num}"
    arr = np.array(data)
    file = f'{path}/uid.npy'
    np.save(file, arr)

def gen_json_file_name(row_based, rows, dim, auto_id, str_pk, float_vector, multi_scalars, file_num):
    suffix = entity_suffix(rows)
    scalars = "only"
    if multi_scalars:
        scalars = "multi_scalars"
    vt = FLOAT
    if not float_vector:
        vt = BINARY
    pk = ""
    if str_pk:
        pk = "str_pk_"
    prefix = gen_file_prefix(row_based=row_based, auto_id=auto_id)
    return f"{prefix}_{pk}{vt}_vectors_{scalars}_{dim}d_{suffix}_{file_num}.json"

def gen_json_files(row_based, rows, dim, auto_id, str_pk, float_vector, multi_scalars, file_nums, force=False):
    """
    row_based: Boolean
        generate row based json file if True
        generate column base json file if False
    rows: entities of data
    dim: dim of vector data
    auto_id: Boolean
        generate primary key data if False, else not
    str_pk: Boolean
         generate string as primary key if True, else generate INT64 as pk
    float_vector: Boolean
        generate float vectors if True, else binary vectors
    multi_scalars: Boolean
        only generate vector data (and pk data depended on auto_switches) if False
        besides vector data, generate INT, STRING, BOOLEAN, etc scalar data if True
    file_nums: file numbers that to be generated
    """

    # gen json files
    for i in range(file_nums):
        file_name = gen_json_file_name(row_based=row_based, rows=rows,dim=dim,
                                       auto_id=auto_id,str_pk=str_pk,float_vector=float_vector,
                                       multi_scalars=multi_scalars, file_num=i)
        file = f"{data_source}/{file_name}"
        if not os.path.exists(file) or force:
            if row_based:
                gen_row_based_json_file(row_file=file, str_pk=str_pk,
                                        float_vect=float_vector, multi_scalars=multi_scalars,
                                        rows=rows, dim=dim, autoid=auto_id)
            else:
                gen_column_base_json_file(col_file=file, str_pk=str_pk,
                                          float_vect=float_vector, multi_scalars=multi_scalars,
                                          rows=rows, dim=dim, autoid=auto_id)

def gen_npy_files():
    # # gen numpy files
    # uid = 0
    # for i in range(file_nums):
    #     gen_vectors_in_numpy_file(data_dir, FLOAT, rows=rows_list[0], dim=dim_list[0], num=i)
    #     gen_scalars_in_numpy_file(data_dir, "int", rows=rows_list[0], dim=dim_list[0], num=i, start=uid)
    #     uid += rows_list[0]
    pass

def copy_files_to_bucket(client, r_source, bucket_name, force=False):
    # check 'xxx' bucket exist.
    found = client.bucket_exists(bucket_name)
    if not found:
        print(f"Bucket {bucket_name} not found, create it.")
        client.make_bucket(bucket_name)

    # copy files from root source folder
    os.chdir(r_source)
    onlyfiles = [f for f in os.listdir(r_source) if
                 os.path.isfile(os.path.join(r_source, f))]
    for file in onlyfiles:
        if not file.startswith("."):
            found = False
            try:
                result = client.stat_object(bucket_name, file)
                found = True
            except S3Error as exc:
                pass

            if force:
                res = client.fput_object(bucket_name, file, f"{r_source}/{file}")
                print(res.object_name)
            elif not found:
                res = client.fput_object(bucket_name, file, f"{r_source}/{file}")
                print(res.object_name)

    # copy subfolders
    sub_folders = [f.name for f in os.scandir(r_source) if f.is_dir()]
    for sub_folder in sub_folders:
        if sub_folder not in ["backup", "tested"]:
            source = f"{r_source}/{sub_folder}"
            os.chdir(source)
            onlyfiles = [f for f in os.listdir(source) if
                         os.path.isfile(os.path.join(source, f))]
            for file in onlyfiles:
                if not file.startswith("."):
                    found = False
                    try:
                        result = client.stat_object(bucket_name, f"{sub_folder}/{file}")
                        found = True
                    except S3Error as exc:
                        pass

                    if force:
                        res = client.fput_object(bucket_name, f"{sub_folder}/{file}", f"{source}/{file}")
                        print(res.object_name)
                    elif not found:
                        res = client.fput_object(bucket_name, f"{sub_folder}/{file}", f"{source}/{file}")
                        print(res.object_name)

def copy_files_to_minio(host, bucket_name, access_key="minioadmin", secret_key="minioadmin", secure=False):
    client = Minio(
        host,
        access_key=access_key,
        secret_key=secret_key,
        secure=secure,
    )
    try:
        # TODO: not copy all the files, just copy the new generated files
        copy_files_to_bucket(client, r_source=data_source, bucket_name=bucket_name, force=False)
    except S3Error as exc:
        print("error occurred.", exc)

def parpar_bulk_load_data(json_file, row_based, rows, dim, auto_id, str_pk, float_vector, multi_scalars, file_nums, force=False):
    if json_file:
        gen_json_files(row_based=row_based, rows=rows, dim=dim,
                       auto_id=auto_id, str_pk=str_pk, float_vector=float_vector,
                       multi_scalars=multi_scalars, file_nums=file_nums, force=force)

        copy_files_to_minio(host=minio, bucket_name=bucket_name)
    else:
        # TODO: for npy files
        # gen_npy_files()
        # copy()
        pass

# if __name__ == '__main__':
#     gen_json_files(row_based=True, rows=10,
#                    dim=4, auto_id=False, str_pk=False,
#                    float_vector=True, multi_scalars=False, file_nums=2)
#
#     copy_files_to_minio(host=minio, bucket_name=bucket_name)

max_stars_repo_name: letian-jiang/milvus

max_issues_repo_name: letian-jiang/milvus

max_forks_repo_name: letian-jiang/milvus


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 14 14:54:18 2020

@author: djrob
"""
import os
from pathlib import Path
import pandas as pd
from IPython.display import display
import os
from pprint import pprint
from pathlib import Path
from pmutt.io.excel import read_excel
import numpy as np
from scipy.interpolate import Rbf, interp1d, BSpline, make_interp_spline
import matplotlib
import matplotlib.pyplot as plt
from scipy import interpolate

excel_path = r'C:\Users\djrob\Desktop\Work\Fall_2020\openMKM v Chemkin\Ammonia\ammonia_transient_data.xlsx'
chemkin_cov_data = read_excel(io=excel_path, sheet_name='chemkin_cov')
openmkm_cov_data = read_excel(io=excel_path, sheet_name='omkm_cov')
chemkin_mass_data = read_excel(io=excel_path, sheet_name='chemkin_mass')
openmkm_mass_data = read_excel(io=excel_path, sheet_name='omkm_mass')
chemkin_time = []
chemkin_conv = []
openmkm_time = []
openmkm_conv = []
chemkin_ru_s = []
chemkin_ru_t = []
openmkm_ru_s = []
openmkm_ru_t = []
omkm_H_s = []
omkm_H_t = []
omkm_NH_s = []
omkm_NH_t = []
chemkin_H_s = []
chemkin_H_t = []
chemkin_NH_s = []
chemkin_NH_t = []
omkm_H2 = []
chemkin_H2 = []
omkm_N2 = []
chemkin_N2 = []

for i in range(0, 1002):
    chemkin_time.append(chemkin_cov_data[i]['Time [s]'])
    chemkin_conv.append(chemkin_mass_data[i]['Conv %'])
    chemkin_ru_s.append(chemkin_cov_data[i]['RU(S)'])
    chemkin_ru_t.append(chemkin_cov_data[i]['RU(T)'])
    chemkin_H_s.append(chemkin_cov_data[i]['H(S)'])
    chemkin_H_t.append(chemkin_cov_data[i]['H(T)'])
    chemkin_NH_s.append(chemkin_cov_data[i]['NH(S)'])
    chemkin_NH_t.append(chemkin_cov_data[i]['NH(T)'])
    chemkin_H2.append(chemkin_mass_data[i]['H2'])
    chemkin_N2.append(chemkin_mass_data[i]['N2'])

for i in range(0, 329):
    openmkm_time.append(openmkm_mass_data[i]['t(s)'])
    openmkm_conv.append(openmkm_mass_data[i]['Conv %'])
    openmkm_ru_s.append(openmkm_cov_data[i]['RU(S)'])
    openmkm_ru_t.append(openmkm_cov_data[i]['RU(T)'])
    omkm_H_s.append(openmkm_cov_data[i]['H(S)'])
    omkm_H_t.append(openmkm_cov_data[i]['H(T)'])
    omkm_NH_s.append(openmkm_cov_data[i]['NH(S)'])
    omkm_NH_t.append(openmkm_cov_data[i]['NH(T)'])
    omkm_H2.append(openmkm_mass_data[i]['H2'])
    omkm_N2.append(openmkm_mass_data[i]['N2'])

interpol_openmkm = make_interp_spline(np.array(openmkm_time), np.array(openmkm_conv),2)
interpol_step = make_interp_spline(np.array(openmkm_time), np.array(openmkm_ru_s),2)
interpol_terr = make_interp_spline(np.array(openmkm_time), np.array(openmkm_ru_t),2)
h2_interpol = make_interp_spline(openmkm_time, omkm_H2, 2)
n2_interpol = make_interp_spline(openmkm_time, omkm_N2, 2)
step_interpol = make_interp_spline(openmkm_time, openmkm_ru_s,2)
terr_interpol = make_interp_spline(openmkm_time, openmkm_ru_t, 2)
h_s_interpol = make_interp_spline(openmkm_time, omkm_H_s, 2)
h_t_interpol = make_interp_spline(openmkm_time, omkm_H_t, 2)
nh_s_interpol = make_interp_spline(openmkm_time, omkm_NH_s, 2)
nh_t_interpol = make_interp_spline(openmkm_time, omkm_NH_t, 2)

'''ABSOLUTE MEAN AND MAX DIFF MASS'''
conv_diff = np.array(chemkin_conv) - np.array(interpol_openmkm(chemkin_time))
conv_avg_diff = np.average(conv_diff)
H2_diff = np.array(chemkin_H2) - np.array(h2_interpol(chemkin_time))
H2_avg_diff = np.average(H2_diff)
N2_diff = np.array(chemkin_N2) - np.array(n2_interpol(chemkin_time))
N2_avg_diff = np.average(N2_diff)
overall_mass_diff = np.concatenate((conv_diff, H2_diff, N2_diff), axis = 0)
max_mass_mean_diff = max(abs(np.array(overall_mass_diff)))
print("Maximum Difference for Mass: " + str(max_mass_mean_diff))
print('Absolute Mean for Mass: ' + str(np.average(abs(overall_mass_diff))))

'''ABSOLUTE MEAN AND MAX DIFF COV'''
ru_s_diff = np.array(chemkin_ru_s) - np.array(interpol_step(chemkin_time))
ru_s_avg_diff = np.average(ru_s_diff)
ru_t_diff = np.array(chemkin_ru_t) - np.array(interpol_terr(chemkin_time))
ru_t_avg_diff = np.average(ru_t_diff)
h_s_diff = np.array(chemkin_H_s) - np.array(h_s_interpol(chemkin_time))
h_s_avg_diff = np.average(h_s_diff)
h_t_diff = np.array(chemkin_H_t) - np.array(h_t_interpol(chemkin_time))
h_t_avg_diff = np.average(h_t_diff)
nh_s_diff = np.array(chemkin_NH_s) - np.array(nh_s_interpol(chemkin_time))
nh_s_avg_diff = np.average(nh_s_diff)
nh_t_diff = np.array(chemkin_NH_t) - np.array(nh_t_interpol(chemkin_time))
nh_t_avg_diff = np.average(nh_t_diff)
overall_cov_diff = np.concatenate((ru_s_diff, ru_t_diff, h_s_diff, h_t_diff, nh_s_diff, nh_t_diff), axis = 0)
max_cov_mean_diff = max(abs(np.array(overall_cov_diff)))
print("Maximum Difference for Cov: " + str(max_cov_mean_diff))
print('Absolute Mean for Cov: ' + str(np.average(abs(overall_cov_diff))))
plt.rc('font', size=28)

plt.figure('Conversion')
omkm = plt.plot(openmkm_time, openmkm_conv, 'k', linewidth=3, label = 'oMKM')
chemkin = plt.plot(chemkin_time, chemkin_conv, '--r', linewidth=3, label = 'Chemkin')
#plt.plot(chemkin_time, interpol_openmkm(chemkin_time), '--b', linewidth=3)
#plt.xticks(fontsize=32)
plt.xlim([0, 25])
plt.ylim([0, 0.65])
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('Conversion %', fontsize=40)
plt.legend(loc ='lower right')

plt.figure('Step Vacancies')
plt.plot(chemkin_time, chemkin_ru_s, '--r', linewidth=3, label = 'Chemkin')
plt.plot(openmkm_time, openmkm_ru_s, 'k', linewidth=3, label = 'oMKM')
plt.xlim([0, 5])
plt.ylim([0.3, 0.5])
#plt.xticks(fontsize=32)
plt.yticks(np.arange(.30, .51, step = 0.05))
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('RU(S) Vacancies', fontsize=40)
plt.legend(loc ='lower right')

plt.figure('Terrace Vacancies')
plt.plot(chemkin_time, chemkin_ru_t, '--r', linewidth=3)
plt.plot(openmkm_time, openmkm_ru_t, 'k', linewidth=3)
plt.xlim([0, 8])
plt.ylim([0, 0.5])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('RU(T) Vacancies', fontsize=40)

plt.figure('H(S) Coverage')
plt.plot(chemkin_time, chemkin_H_s, '--r', linewidth=3)
plt.plot(openmkm_time, omkm_H_s, 'k', linewidth=3)
plt.xlim([0, 10])
plt.ylim([0.00, 0.35])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('H(S) Coverage', fontsize=40)

plt.figure('H(T) Coverage')
plt.plot(chemkin_time, chemkin_H_t, '--r', linewidth=3)
plt.plot(openmkm_time, omkm_H_t, 'k', linewidth=3)
plt.xlim([0, 15])
plt.ylim([0, 0.3])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('H(T) Coverage', fontsize=40)

plt.figure('NH(S) Coverage')
plt.plot(chemkin_time, chemkin_NH_s, '--r', linewidth=3)
plt.plot(openmkm_time, omkm_NH_s, 'k', linewidth=3)
plt.xlim([0, 8])
plt.ylim([0, 0.50])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('NH(S) Coverage', fontsize=40)

plt.figure('NH(T) Coverage')
plt.plot(chemkin_time, chemkin_NH_t, '--r', linewidth=3, label = 'Chemkin')
plt.plot(openmkm_time, omkm_NH_t, 'k', linewidth=3, label = 'oMKM')
plt.xlim([0, 10])
plt.ylim([0.2, 0.7])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('NH(T) Coverage', fontsize=40)
plt.legend(loc ='upper right')

plt.figure('H2 Mass')
plt.plot(chemkin_time, chemkin_H2, '--r', linewidth=3)
plt.plot(openmkm_time, omkm_H2, 'k', linewidth=3)
plt.xlim([0, 50])
plt.ylim([0, .12])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('H2 Mass Fraction', fontsize=40)

plt.figure('N2 Mass')
plt.plot(chemkin_time, chemkin_N2, '--r', linewidth=3)
plt.plot(openmkm_time, omkm_N2, 'k', linewidth=3)
plt.xlim([0, 25])
plt.ylim([0, 0.55])
#plt.xticks(fontsize=32)
#plt.yticks(fontsize=32)
plt.xlabel('Time[s]', fontsize=40)
plt.ylabel('N2 Mass Fraction', fontsize=40)

max_stars_repo_name: skasiraj/Parameter-Estimation-BO

max_issues_repo_name: skasiraj/Parameter-Estimation-BO

max_forks_repo_name: skasiraj/Parameter-Estimation-BO

reason: import numpy,from scipy


text:
"""
TAKEN FROM lsdnumbatools, temporary measure before the latter is available
this provides a numba jitclass to loop through node effiviently

B.G.
"""
import numpy as np
import numba as nb
import math as m
from numba.experimental import jitclass

@nb.njit(inline='always')
def _get_neighbouring_nodes_with_checks(output,adder,node):
	i = 0
	for add in adder:
		if(add != -2147483648):
			output[i] = node + add
		else:
			output[i] = -2147483648
		i += 1

# specs for the jitclass which needs to be statically typed
spec = [
		('node_type', nb.int16[:]),
		('nx', nb.int32),
		('ny', nb.int32),
		('nxy', nb.int32),
		('dx', nb.float64),
		('dy', nb.float64),
		('dxy', nb.float64),
		('nodata', nb.int32),
		('neighbourer', nb.int32[:,:]),
		('lengthener', nb.float64[:]),
	]
	#min, max int32 = -2,147,483,648 2,147,483,647

@jitclass(spec)
class iterasator(object):
	"""
		The iterasator class is made to return neighboring informations from a node.
		It takes care of boundary conditions for example, where the geometric neighbors are not the right ones.
		Date: November 2020
		Authors: B.G.
	"""
	def __init__(self,node_type,nx,ny, dx, dy):
		"""
		Initialises the iterasator.
		node_type (1D numpy array of int): vectorised array of node type
		nx (int): number of cols
		ny (int): number of rows
		Date: November 2020
		Authors: B.G.
		"""
		# The nodes type are documented in the file node_graph.py
		self.node_type = node_type
		# Number of cols
		self.nx = nx
		# Number of rows
		self.ny = ny

		#number of elements
		self.nxy = nx*ny

		# No data value for the neighbourer
		self.nodata = -2147483648

		self.dx = dx
		self.dy = dy
		self.dxy = m.sqrt(dx**2 + dy**2)
		# The neighbourer is the most important (and irksome to code) structure of this class:
		# It hardcodes the window looping through vectorised neighbours for each cases
		# I adopted this approach rather than padding no data to be as flexible as possible
		self.neighbourer = np.array([
[ - nx - 1, - nx, - nx + 1, -1, 1, nx - 1, nx, nx + 1], # neighbourer 0: classic internal D8 node and all the neighbors are returned
[ - 1, - nx, - nx + 1, nx - 1, 1, 2 * nx - 1, nx, nx + 1], # neighbourer 1: periodic D8 node at the Western Boundary
[ - 1 + ny * nx,  ny * nx, ny * nx + 1, -1, 1, nx - 1, nx, nx + 1], # neighbourer 2: periodic D8 node at the Northern Boundary
[ - nx - 1, - nx, - 2 * nx + 1, -1, - nx + 1, nx - 1, nx, 1], # neighbourer 3: periodic D8 node at the Eastern Boundary
[ - nx - 1, - nx, - nx + 1, -1, 1, - ny * nx - 1, - ny * nx, - ny * nx + 1], # neighbourer 4: periodic D8 node at the Southern Boundary
[ self.nodata, self.nodata, self.nodata, -1, 1, nx - 1, nx, nx + 1], # neighbourer 5: classic D8 node Northern boundary
[ self.nodata, - nx, - nx + 1, self.nodata, 1, self.nodata, nx, nx + 1], # neighbourer 6: classic D8 node Western boundary
[ - nx - 1, - nx, - nx + 1, -1, 1, self.nodata, self.nodata, self.nodata], # neighbourer 7: classic D8 node Southern boundary
[ - nx - 1, - nx, self.nodata, -1, self.nodata, nx - 1, nx, self.nodata], # neighbourer 8: classic D8 node Eastern boundary
[ self.nodata,  self.nodata, self.nodata, self.nodata, 1, self.nodata, nx, nx + 1], # neighbourer 9: classic D8 node at the Northern-western corner Boundary
[ self.nodata, - nx, - nx + 1, self.nodata, 1, self.nodata, self.nodata, self.nodata], # neighbourer 10: classic D8 node at the Southwestern corner Boundary
[ - nx - 1, - nx, self.nodata, -1, self.nodata, self.nodata, self.nodata, self.nodata], # neighbourer 11: classic D8 node at the Southeastern corner Boundary
[ self.nodata, self.nodata, self.nodata, -1, self.nodata, nx - 1, nx, self.nodata], # neighbourer 12: classic D8 node at the Northeastern corner Boundary
# To add for the future, not urgent:
# - Cyclic boundaries cause why not
# - specific corner boundary (e.g. NW with W periodic and NW with N periodic)
# - Internal periodic? (i.e. goes straight to the next periodic in the interested direction)
		],dtype = nb.int32)

		# The lengthener is associated with the neighbourer to return the length to the neighbour associated to each one
		self.lengthener = np.array([self.dxy,self.dy,self.dxy,self.dx,self.dx,self.dxy,self.dy,self.dxy], nb.float64)

	def node2rowcol(self,node):
		"""
		Converts vectorised node ID to row index and col index
		"""
		r = node // self.nx
		c = node % self.nx
		return r,c

	def get_neighbouring_nodes(self,node):
		"""
		return all the neighboring nodes of a diven node index
		"""
		output = np.full(8,self.nodata,dtype = nb.int32)
		# case 1: active, classic internal node, keeping it test-free to optimise code as this will be the overwhelming majority of cases
		# THE ONLY REQUIREMENT IS THAT NO BOUNDARY NODES ARE LABELLED WITH 1
		if(self.node_type[node] == 1):
			i = 0
			for add in self.neighbourer[0]:
				output[i] = node + add
				i += 1
		# North row (in array point of view)
		elif(node < self.nx):
			# classic border bourndary
			if(self.node_type[node] == 0):
				# N case
				if(node > 0 and node < self.nx - 1):
					_get_neighbouring_nodes_with_checks(output,self.neighbourer[5],node)
				# NW case
				elif (node == 0):
					_get_neighbouring_nodes_with_checks(output,self.neighbourer[9],node)
				# NE case
				elif (node == self.nx - 1):
					_get_neighbouring_nodes_with_checks(output,self.neighbourer[12],node)
			# periodic boundaries
			if(self.node_type[node] == 2):
				_get_neighbouring_nodes_with_checks(output,self.neighbourer[2],node)
		# South row
		elif(node >= (self.nx * self.ny) - self.nx):
			# classic border bourndary
			if(self.node_type[node] == 0):
				# Middle of the S row
				if(node > (self.nx * self.ny - self.nx) and node < (self.nx * self.ny - self.nx) - 1):
					_get_neighbouring_nodes_with_checks(output,self.neighbourer[7],node)
				# SW classic boulder
				elif (node == (self.nx * self.ny - self.nx)):
					_get_neighbouring_nodes_with_checks(output,self.neighbourer[10],node)
					# SE classic boundary
				elif (node == (self.nx * self.ny - self.nx) - 1):
					_get_neighbouring_nodes_with_checks(output,self.neighbourer[11],node)
			elif(self.node_type[node] == 2):
				_get_neighbouring_nodes_with_checks(output,self.neighbourer[4],node)
	    # first Column
		elif(node % self.nx == 0):
			if(self.node_type[node] == 0):
				_get_neighbouring_nodes_with_checks(output,self.neighbourer[6],node)
			elif(self.node_type[node] == 2):
				_get_neighbouring_nodes_with_checks(output,self.neighbourer[1],node)
	    # Last Column
		elif(node % self.nx - 1 == 0):
			if(self.node_type[node] == 0):
				_get_neighbouring_nodes_with_checks(output,self.neighbourer[8],node)
			elif(self.node_type[node] == 2):
				_get_neighbouring_nodes_with_checks(output,self.neighbourer[3],node)

		return output

#End of fucking file (I am extremely frustrated at Sublime text which delete the last empty rows of a file automatically and wants you to then work using obly the absolute bottom of your screen. I am not alone right?)

max_stars_repo_name: fastscape-lem/fastscape-litho

max_issues_repo_name: fastscape-lem/fastscape-litho

max_forks_repo_name: fastscape-lem/fastscape-litho

reason: import numpy,import numba,from numba


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np
from src.cosine_annealing import CosineAnnealingWithRestartsLR
from src.adamw import AdamW

class ScheduledOptimizer:
    """
    If otype is adam, then optimizer will take regular step to decay learning rate
    If otype is sgd, then optimizer will take irregular step to decay learning rate, which can be specified
    If otype is sgdr, then optimizer will take cosine annealing with restart to decay learning rate
    If otype is adamw, then optimizer will take regular step to decay learning rate
    """

    def __init__(self, parameters, otype, lr, lr_scheduler_step_size, lr_gamma,
                 lr_clip=1e-5, lr_weight_decay=1e-3, lr_sgd_momentum=0.9, T_mult=2):
        if otype == "adam":
            self.optimizer = torch.optim.Adam(parameters, lr=lr, weight_decay=lr_weight_decay)
            self.scheduler = torch.optim.lr_scheduler.StepLR(self.optimizer, step_size=lr_scheduler_step_size, gamma=lr_gamma)
        elif otype == "sgd":
            self.optimizer = torch.optim.SGD(parameters, lr=lr, momentum=lr_sgd_momentum, weight_decay=lr_weight_decay)
            self.scheduler = torch.optim.lr_scheduler.MultiStepLR(self.optimizer, milestones=lr_scheduler_step_size, gamma=lr_gamma)
        elif otype == "sgdr":
            self.optimizer = torch.optim.SGD(parameters, lr=lr, momentum=lr_sgd_momentum, weight_decay=lr_weight_decay)
            self.scheduler = CosineAnnealingWithRestartsLR(self.optimizer, T_max=lr_scheduler_step_size,
                                                           eta_min=lr_clip, last_epoch=-1, T_mult=T_mult)
        elif otype == "adamw":
            self.optimizer = AdamW(parameters, lr=lr, weight_decay=lr_weight_decay)
            self.scheduler = torch.optim.lr_scheduler.StepLR(self.optimizer, step_size=lr_scheduler_step_size, gamma=lr_gamma)

        self.lr_clip = lr_clip

    def zero_grad(self):
        self.optimizer.zero_grad()

    def step(self):
        self.optimizer.step()

    def scheduler_step(self):
        self.scheduler.step()
        for param_group in self.optimizer.param_groups:
            param_group['lr'] = max(param_group['lr'], self.lr_clip)

class BNMomentumScheduler:
    def __init__(self, net, init_momentum=0.1, step_size=40, gamma=0.5, clip=0.01):
        self.step_size = step_size
        self.gamma = gamma
        self.clip = clip
        self.epoch = 0
        self.net = net

        for m in self.net.modules():
            if isinstance(m, nn.BatchNorm3d) or isinstance(m, nn.BatchNorm1d) or isinstance(m, nn.BatchNorm2d):
                m.momentum = init_momentum

    def scheduler_step(self):
        self.epoch += 1
        if self.epoch % self.step_size == 0:
            for m in self.net.modules():
                if isinstance(m, nn.BatchNorm3d) or isinstance(m, nn.BatchNorm1d) or isinstance(m, nn.BatchNorm2d):
                    m.momentum = max(m.momentum*self.gamma, self.clip)

def label_smooth_loss(log_prob, label, confidence=0.9):
    """
    :param log_prob: log probability
    :param label: one hot encoded
    :param confidence: we replace one (in the one hot) with confidence. 0 <= confidence <= 1.
    :return:
    """
    N = log_prob.size(0)
    C = log_prob.size(1)
    smoothed_label = torch.full(size=(N, C), fill_value=(1-confidence) / (C - 1)).to(log_prob)
    smoothed_label.scatter_(dim=1, index=torch.unsqueeze(label, dim=1), value=confidence)
    loss = - torch.sum(log_prob * smoothed_label) / N
    return loss

def cal_label_weight(label):
    num_label = np.unique(label).shape[0]
    label_weight = np.zeros(num_label).astype(np.float32)
    for i in range(num_label):
        label_weight[i] = (label == i).sum()
    label_weight = label_weight / label_weight.sum()
    label_weight = 1 / np.log(1 + label_weight)
    label_weight = label_weight / label_weight.sum()
    label_weight = torch.from_numpy(label_weight)
    return label_weight

class Conv1dBatchRelu(nn.Module):
    def __init__(self, in_channel, out_channel, kernel_size=1, stride=1, padding=0, batch_norm=True):
        super(Conv1dBatchRelu, self).__init__()
        if batch_norm:
            self.block = nn.Sequential(
                nn.Conv1d(in_channel, out_channel, kernel_size=kernel_size, stride=stride, padding=padding),
                nn.BatchNorm1d(out_channel),
                nn.ReLU(inplace=True))
        else:
            self.block = nn.Sequential(
                nn.Conv1d(in_channel, out_channel, kernel_size=kernel_size, stride=stride, padding=padding),
                nn.ReLU(inplace=True))

    def forward(self, x):
        x = self.block(x)
        return x

class Conv2dBatchRelu(nn.Module):
    def __init__(self, in_channel, out_channel, kernel_size=3, stride=1, padding=1, dilation=1, batch_norm=True):
        super(Conv2dBatchRelu, self).__init__()
        if batch_norm:
            self.block = nn.Sequential(
                nn.Conv2d(in_channel, out_channel, kernel_size=kernel_size, stride=stride, padding=padding, dilation=dilation),
                nn.BatchNorm2d(out_channel),
                nn.ReLU(inplace=True))
        else:
            self.block = nn.Sequential(
                nn.Conv2d(in_channel, out_channel, kernel_size=kernel_size, stride=stride, padding=padding, dilation=dilation),
                nn.ReLU(inplace=True))

    def forward(self, x):
        x = self.block(x)
        return x

class Conv3dBatchRelu(nn.Module):
    def __init__(self, in_channel, out_channel, kernel_size=3, stride=1, padding=1, dilation=1, batch_norm=True):
        super(Conv3dBatchRelu, self).__init__()
        if batch_norm:
            self.block = nn.Sequential(
                nn.Conv3d(in_channel, out_channel, kernel_size=kernel_size, stride=stride, padding=padding, dilation=dilation),
                nn.BatchNorm3d(out_channel),
                nn.ReLU(inplace=True))
        else:
            self.block = nn.Sequential(
                nn.Conv3d(in_channel, out_channel, kernel_size=kernel_size, stride=stride, padding=padding, dilation=dilation),
                nn.ReLU(inplace=True))

    def forward(self, x):
        x = self.block(x)
        return x

class CatDilationConv2dBatchRelu(nn.Module):
    def __init__(self, in_channel, out_channel, dilation=(1, 2), batch_norm=True):
        super(CatDilationConv2dBatchRelu, self).__init__()
        self.dilation_num = len(dilation)
        assert out_channel % self.dilation_num == 0
        self.block = nn.ModuleList([Conv2dBatchRelu(in_channel, int(out_channel/self.dilation_num), batch_norm=batch_norm,
                                                    kernel_size=3, stride=1, padding=dilation[i], dilation=dilation[i])
                                    for i in range(self.dilation_num)])

    def forward(self, x):
        x = torch.cat([self.block[i](x) for i in range(self.dilation_num)], dim=1)
        return x

class CatDilationConv3dBatchRelu(nn.Module):
    def __init__(self, in_channel, out_channel, dilation=(1, 2), batch_norm=True):
        super(CatDilationConv3dBatchRelu, self).__init__()
        self.dilation_num = len(dilation)
        assert out_channel % self.dilation_num == 0
        self.block = nn.ModuleList([Conv3dBatchRelu(in_channel, int(out_channel/self.dilation_num), batch_norm=batch_norm,
                                                    kernel_size=3, stride=1, padding=dilation[i], dilation=dilation[i])
                                    for i in range(self.dilation_num)])

    def forward(self, x):
        x = torch.cat([self.block[i](x) for i in range(self.dilation_num)], dim=1)
        return x

class LinearBatchRelu(nn.Module):
    def __init__(self, in_feature, out_feature, batch_norm=True):
        super(LinearBatchRelu, self).__init__()
        if batch_norm:
            self.block = nn.Sequential(
                nn.Linear(in_features=in_feature, out_features=out_feature),
                nn.BatchNorm1d(out_feature),
                nn.ReLU(inplace=True))
        else:
            self.block = nn.Sequential(
                nn.Linear(in_features=in_feature, out_features=out_feature),
                nn.ReLU(inplace=True))

    def forward(self, x):
        return self.block(x)

def one_hot_encoding(x):
    one_hot = np.zeros((len(x), x.max()+1))
    one_hot[np.arange(len(x)), x] = 1
    return one_hot

def torch_one_hot_encoding(x):
    x = torch.LongTensor(x).view(-1, 1)
    one_hot = torch.zeros((len(x), x.max() + 1))
    one_hot.scatter_(dim=1, index=x, value=1)
    return one_hot

if __name__ == '__main__':
    x = torch.tensor([1,2,1,3])
    print(torch_one_hot_encoding(x))

max_stars_repo_name: DIUDIUDIUDIUDIU/AdvectiveNet

max_issues_repo_name: DIUDIUDIUDIUDIU/AdvectiveNet

max_forks_repo_name: DIUDIUDIUDIUDIU/AdvectiveNet


text:
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from math import *
import time
# User-friendly machine learning library
# Front end for TensorFlow
import tensorflow.keras
# Different methods from Keras needed to create an RNN
# This is not necessary but it shortened function calls
# that need to be used in the code.
from tensorflow.keras import datasets, layers, models
from tensorflow.keras.layers import Input
from tensorflow.keras import regularizers
from tensorflow.keras.models import Model, Sequential
#from tensorflow.keras.layers.core import Dense, Activation
from tensorflow.keras.layers import Dense, SimpleRNN, LSTM, GRU
# For timing the code
from timeit import default_timer as timer
# For plotting
import matplotlib.pyplot as plt

# Define Analytical, Euler-Cromer, and Velocity-Verlet methods of solving
def analytical(k,m,x0,v0,dt,tfinal):
    t = np.arange(0,tfinal+dt,dt)
    v = -x0 * np.sin(t) + v0 * np.cos(t)
    x = x0 * np.cos(t) + v0 * np.sin(t)
    K = 1/2 *m*v**2
    U = 1/2 *k*x**2
    return x, v, K, U, t

def euler_cromer(k,m,x0,v0,dt,tfinal):
    n = np.ceil(tfinal/dt)\
    # Set up arrays
    t = np.zeros(n)
    v = np.zeros(n)
    x = np.zeros(n)
    K = np.zeros(n)
    U = np.zeros(n)
    # Define Initial Conditions
    x[0] = x0
    v[0] = v0
    K[0] = 1/2 *m*v0**2
    U[0] = 1/2 *k*x0**2

    # Integrate using the Euler-Cromer Method
    for i in range(n-1):
        a = -x[i]
        v[i+1] = v[i] + dt*a
        x[i+1] = x[i] + dt*v[i+1]
        K[i+1] = 1/2 *m*v[i+1]**2
        U[i+1] = 1/2 *k*x[i+1]**2
        t[i+1] = t[i] + dt
    return x, v, K, U, t

def velocity_verlet(k,m,x0,v0,dt,tfinal):
    n = np.ceil(tfinal/dt)
    # Set up arrays
    t = np.zeros(n)
    v = np.zeros(n)
    x = np.zeros(n)
    K = np.zeros(n)
    U = np.zeros(n)
    # Define Initial Conditions
    x[0] = x0
    v[0] = v0
    K[0] = 1/2 *m*v0**2
    U[0] = 1/2 *k*x0**2

    # Integrate using the Velocity-Verlet Method
    for i in range(n-1):
        a = -x[i]
        x[i+1] = x[i] + dt*v[i] + dt**2 /2*a
        a1 = -x[i+1]
        v[i+1] = v[i] + dt/2*(a+a1)
        K[i+1] = 1/2 *m*v[i+1]**2
        U[i+1] = 1/2 *k*x[i+1]**2
        t[i+1] = t[i] + dt
    return x, v, K, U, t

# Define Constants
dt = 0.01
tfinal = 50
x0 = 1
v0 = 0
m = 1
k = 1

# Call the Integration Function
ax, av, aK, aU, at = analytical(k,m,x0,v0,dt,tfinal)
ecx, ecv, ecK, ecU, ect = euler_cromer(k,m,x0,v0,dt,tfinal)
vvx, vvv, vvK, vvU,vvt = velocity_verlet(k,m,x0,v0,dt,tfinal)

# Plots
fig, axes = plt.subplots(1,3, figsize = (15,5))
fig.suptitle("System of an Undamped Spring", fontsize=14, y = 1.05)

axes[0].plot(at, ax, label = "Analytical Method")
axes[0].plot(ect, ecx, label = "Euler-Cromer Method")
axes[0].plot(vvt, vvx, label = "Velocity-Verlet Method")
axes[0].set_title("Position as a Function of Dimensionless Time")
axes[0].set_xlabel("Dimensionless Time")
axes[0].set_ylabel("Position")

axes[1].plot(at, av, label = "Analytical Method")
axes[1].plot(ect, ecv, label = "Euler-Cromer Method")
axes[1].plot(vvt, vvv, label = "Velocity-Verlet Method")
axes[1].set_title("Velocity as a Function of Dimensionless Time")
axes[1].set_xlabel("Dimensionless Time")
axes[1].set_ylabel("Velocity")

axes[2].plot(at, aU+aK, label = "Analytical Method")
axes[2].plot(ect, ecU+ecK, label = "Euler-Cromer Method")
axes[2].plot(vvt, vvU+vvK, label = "Velocity-Verlet Method")
axes[2].set_title("Energy as a Function of Dimensionless Time")
axes[2].set_xlabel("Dimensionless Time")
axes[2].set_ylabel("Total Energy")

plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

plt.tight_layout()

def damp(gamma,m,x0,v0,DeltaT,tfinal):
    n = np.ceil(tfinal/DeltaT)
    # Set up arrays
    t = np.zeros(n)
    v = np.zeros(n)
    r = np.zeros(n)
    # Define Initial Conditions
    r[0] = x0
    v[0] = v0
    # Integrate over using the Velocity Verlet Method
    for i in range(n-1):
        a = -r[i] - 2*gamma*v[i]
        r[i+1] = r[i] + DeltaT*v[i] + DeltaT**2 /2 *a
        a1 = -r[i+1] - 2*gamma*v[i]
        v[i+1] = v[i] + DeltaT/2*(a+a1)
        t[i+1] = t[i] + DeltaT

    return t,r,v

# Define Constants
under_gamma = 0.1
crit_gamma = 1
over_gamma = 2
m = 1
x0 = 1
v0 = 0
dt  = 0.1
tfinal = 50

# Call the Integration Function
t, under_r, under_v = damp(under_gamma,m,x0,v0,dt,tfinal)
t, crit_r, crit_v = damp(crit_gamma,m,x0,v0,dt,tfinal)
t, over_r, over_v = damp(over_gamma,m,x0,v0,dt,tfinal)

# Plots
fig, axes = plt.subplots(1,2, figsize = (15,5))
fig.suptitle("System of a Damped Spring", fontsize=14, y = 1.05)

axes[0].plot(t, under_r, label = "Under Damping")
axes[0].plot(t, crit_r, label = "Critical Damping")
axes[0].plot(t, over_r, label = "Over Damping")
axes[0].set_title("Position as a Function of Dimensionless Time")
axes[0].set_xlabel("Dimensionless Time")
axes[0].set_ylabel("Position")

axes[1].plot(t, under_v, label = "Under Damping")
axes[1].plot(t, crit_v, label = "Critical Damping")
axes[1].plot(t, over_v, label = "Over Damping")
axes[1].set_title("Velocity as a Function of Dimensionless Time")
axes[1].set_xlabel("Dimensionless Time")
axes[1].set_ylabel("Velocity")

plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

plt.tight_layout()

def forced(gamma,r0,v0,F0,omega,d,DeltaT,tfinal):
    n = np.ceil(tfinal/DeltaT)
    # Set up arrays
    t = np.zeros(n)
    v = np.zeros(n)
    r = np.zeros(n)
    # Define Initial Conditions
    r[0] = r0
    v[0] = v0

    # Integrate using the 4th-Order RK Method
    for i in range(n-1):
        t[i+1] = t[i] + DeltaT

        Force = (-r[i] - 2*gamma*v[i] - F0*np.cos(omega*t[i]-d))*m
        k1x = DeltaT*v[i]
        k1v = DeltaT*Force

        vv = v[i]+k1v*0.5*DeltaT
        rr = r[i]+k1x*0.5*DeltaT
        Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m
        k2x = DeltaT*vv
        k2v = DeltaT*Force

        vv = v[i]+k2v*0.5*DeltaT
        rr = r[i]+k2x*0.5*DeltaT
        Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m
        k3x = DeltaT*vv
        k3v = DeltaT*Force

        vv = v[i]+k3v*DeltaT
        rr = r[i]+k3x*DeltaT
        Force = (-rr - 2*gamma*vv - F0*np.cos(omega*(t[i]+DeltaT*0.5)-d))*m
        k4x = DeltaT*vv
        k4v = DeltaT*Force

        r[i+1] = r[i]+(k1x+2*k2x+2*k3x+k4x)/6.
        v[i+1] = v[i]+(k1v+2*k2v+2*k3v+k4v)/6.

    return t,r,v

# Define Constants
gamma = 0.1
r0 = 1
v0 = 0
F0 = 3
omega = 3
d = 0
dt = 0.1
tfinal = 100

# Call the Integration Function
t, r, v = forced(gamma,r0,v0,F0,omega,d,dt,tfinal)

# Plots
fig, axes = plt.subplots(1,2, figsize = (15,5))
fig.suptitle("System of a Damped Spring with a Driving Force", fontsize=14, y = 1.05)

axes[0].plot(t, r)
axes[0].set_title("Position as a Function of Dimensionless Time")
axes[0].set_xlabel("Dimensionless Time")
axes[0].set_ylabel("Position")

axes[1].plot(t, v)
axes[1].set_title("Velocity as a Function of Dimensionless Time")
axes[1].set_xlabel("Dimensionless Time")
axes[1].set_ylabel("Velocity")

plt.tight_layout()

max_stars_repo_name: anacost/MachineLearning

max_issues_repo_name: anacost/MachineLearning

max_forks_repo_name: anacost/MachineLearning


text:
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import Pipeline
from sklearn import svm
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import SGDClassifier
from sklearn.model_selection import cross_val_predict
from sklearn.preprocessing import StandardScaler
from skimage.feature import hog
from skimage.io import imread
from skimage.transform import rescale
import skimage
import joblib
from skimage import io
import numpy as np
import re
from sklearn import svm, metrics
from skimage import feature, filters

class ImageLoader:
    def __init__(self):
        self.classifer = None
        ## TODO: Change folder
        self.folder = '../app/data/'

    def imread_convert(self, f):
        return io.imread(f).astype(np.uint8)

    def load_data_from_folder(self, dir):
        # read all images into an image collection
        ic = io.ImageCollection(self.folder + dir + "*.bmp",
                                load_func=self.imread_convert)

        #create one large array of image data
        data = io.concatenate_images(ic)

        #extract labels from image names
        labels = np.array(ic.files)
        for i, f in enumerate(labels):
            m = re.search("_", f)
            labels[i] = (f[len(dir):m.start()]).split('/')[-1]
            #labels[i] = f[len(dir):m.start()]

        return(data,labels)

class GaussianBlurTransformer(BaseEstimator, TransformerMixin):
    """
    Perform Guassian blurring of greyscale images
    """

    def __init__(self):
        pass

    def fit(self, X, y=None):
        """returns itself"""
        return self

    def transform(self, X, y=None):
        """perform the transformation and return an array"""
        return np.array([filters.gaussian(img) for img in X])

class RGB2GrayTransformer(BaseEstimator, TransformerMixin):
    """
    Convert an array of RGB images to grayscale
    """

    def __init__(self):
        pass

    def fit(self, X, y=None):
        """returns itself"""
        return self

    def transform(self, X, y=None):
        """perform the transformation and return an array"""
        return np.array([skimage.color.rgb2gray(img) for img in X])

class HogTransformer(BaseEstimator, TransformerMixin):
    """
    Expects an array of 2d arrays (1 channel images)
    Calculates hog features for each img
    """

    def __init__(self, y=None, orientations=9,
                 pixels_per_cell=(8, 8),
                 cells_per_block=(3, 3), block_norm='L2-Hys'):
        self.y = y
        self.orientations = orientations
        self.pixels_per_cell = pixels_per_cell
        self.cells_per_block = cells_per_block
        self.block_norm = block_norm

    def fit(self, X, y=None):
        return self

    def transform(self, X, y=None):

        def local_hog(X):
            return hog(X,
                       orientations=self.orientations,
                       pixels_per_cell=self.pixels_per_cell,
                       cells_per_block=self.cells_per_block,
                       block_norm=self.block_norm)

        try: # parallel
            return np.array([local_hog(img) for img in X])
        except:
            return np.array([local_hog(img) for img in X])

HOG_pipeline = Pipeline([
    ('grayify', RGB2GrayTransformer()),
    ('blurify', GaussianBlurTransformer()),
    ('hogify', HogTransformer(
        pixels_per_cell=(8, 8),
        cells_per_block=(2,2),
        orientations=9,
        block_norm='L2-Hys')
    ),
    ('scalify', StandardScaler()),
    ('classify', SGDClassifier(random_state=42, max_iter=1000, tol=1e-3))
])

param_grid = [
    {'hogify__orientations': [9],
    'hogify__cells_per_block': [(3, 3)],
    'hogify__pixels_per_cell': [(8, 8), (14, 14)]},
    {'hogify__orientations': [9],
     'hogify__cells_per_block': [(3, 3)],
     'hogify__pixels_per_cell': [(14, 14)],
     'classify': [
         SGDClassifier(random_state=42, max_iter=1000, tol=1e-3),
         svm.SVC(kernel='linear')]}
]

max_stars_repo_name: baconandchips/myImageClassifier

max_issues_repo_name: baconandchips/myImageClassifier

max_forks_repo_name: baconandchips/myImageClassifier


text:
import numpy as np
class MyTuple:
    '''
    The basic object representing the input variable 'w'
    represents the core of our AD calculator.  An instance
    of this class is a tuple containining one function/derivative
    evaluation of the variable 'w'.  Because it is meant to
    represent the simple variable 'w' the derivative 'der' is
    preset to 1.  The value 'val' can be set to 0 by default.
    '''
    def __init__(self,**kwargs):
        # variables for the value (val) and derivative (der) of our input function
        self.val = 0
        self.der = 1

        # re-assign these default values
        if 'val' in kwargs:
            self.val = kwargs['val']
        if 'der' in kwargs:
            self.der = kwargs['der']

# our implementation of the sinusoid rule from Table 1
def sin(a):
    # Create output evaluation and derivative object
    b = MyTuple()

    # Produce new function value
    b.val = np.sin(a.val)

    # Produce new derivative value - we need to use the chain rule here!
    b.der = np.cos(a.val)*a.der

    # Return updated object
    return b

# create instance of our function to differentiate - notice this uses our homemade sine function not numpy's
g = lambda w: sin(w)

# initialize our AutoDiff object at each point
a1 = MyTuple(val = 0); a2 = MyTuple(val = 0.5)

# evaluate
result1 = g(a1); result2 = g(a2)

# print results
#print ('function value at ' + str(0) + ' = ' + str(result1.val))

#print ('derivaive value at ' + str(0) + ' = ' + str(result1.der))
#print ('function value at ' + str(0.5) + ' = ' + str(result2.val))
#print ('derivaive value at ' + str(0.5) + ' = ' + str(result2.der))

print(sin(MyTuple(val = 5)).der)


text:
# Copyright 2019-2020 ETH Zurich and the DaCe authors. All rights reserved.
import dace
from dace.transformation import optimizer
from dace.transformation.dataflow import GPUTransformMap
import numpy as np
import pytest

@dace.program
def prog1(A: dace.float32[32], B: dace.float32[32]):
    @dace.map
    def work1(i: _[0:32]):
        a << A[i]
        b >> B[i]
        b = a * 2.0

@dace.program
def prog2(A: dace.float32[32], B: dace.float32[32]):
    @dace.map
    def work2(i: _[0:32]):
        a << A[i]
        b >> B[i]
        b = a / 2.0

######################################
@pytest.mark.gpu
def test_multiprogram():
    print('Multi-program CUDA test')

    A = np.random.rand(32).astype(np.float32)
    B = np.random.rand(32).astype(np.float32)
    C = np.random.rand(32).astype(np.float32)

    s1 = prog1.to_sdfg()
    s1.apply_transformations(GPUTransformMap)

    s2 = prog2.to_sdfg()
    s2.apply_transformations(GPUTransformMap)

    s1func = s1.compile()
    s2func = s2.compile()

    s1func(A=A, B=B)
    s2func(A=B, B=C)

    diff = np.linalg.norm(A - C)

    print('Difference:', diff)
    assert diff <= 1e-5

if __name__ == '__main__':
    test()

max_stars_repo_name: tobiasholenstein/dace

max_issues_repo_name: tobiasholenstein/dace

max_forks_repo_name: tobiasholenstein/dace


text:
import numpy as np
import scipy
import scipy.optimize

from . import units, utils

def center(coordinates, masses):
    """
    Given coordinates and masses, return coordinates with center of mass at 0.
    Also works to remove COM translational motion.

    Args:
        coordinates ({nparticle, ndim} ndarray): xyz (to compute COM) or velocities (COM velocity)
        masses ({nparticle,} array_like): masses

    Returns:
        centered_coordinates ({nparticle, ndim} ndarray): xyz with COM at origin or velocities with COM velocity removed
    """
    com_coordinates = utils.compute_center_of_mass(coordinates, masses)
    new_coordinates = np.array(coordinates) - com_coordinates
    return new_coordinates

def angular_momentum_components(coordinates, velocities, masses):
    """
    Computes the minimal change in 3D momentum vectors to remove the angular
    momentum by solving

    min. || \Delta p ||^2_2
    s.t. x \cross (p + \Delta p) = 0

    Args:
        coordinates ({nparticle, 3} ndarray): xyz coordinates
        velocities ({nparticle, 3} ndarray): particle velocities
        masses ({nparticle,} array_like): masses

    Returns:
        dv ({nparticle, 3} ndarray): minimal perturbation to velocities that have no angular momentum, such that v + dv has no angular momentum
    """
    mass_r = np.reshape(masses, (-1, 1))
    linear_momentum = velocities * mass_r
    angular_constraint_mat = np.reshape(np.transpose(utils.cross_product_matrix(coordinates), (1, 0, 2)), (3, -1)) # 3 x (N x 3)
    constraint_vec = - np.dot(angular_constraint_mat, np.reshape(linear_momentum, (-1, 1))) # 3 x 1
    dp = np.linalg.lstsq(angular_constraint_mat, constraint_vec, rcond=None)[0]
    dv = np.reshape(dp, (-1, 3)) / mass_r
    return dv

def initialize_centered(coordinates, velocities, masses):
    """
    Initialize coordinates to have COM at the origin and velocities to have COM
    translation and angular momentum removed by minimizing

    min. || \Delta p ||^2_2
    s.t. x \cross (p + \Delta p) = 0
         \sum_i (p_i + \Delta p_i) = 0
         \sum_i x_i m_i = 0

    Args:
        coordinates ({nparticle, 3} ndarray): particle coordinates
        velocities ({nparticle, 3} ndarray): particle velocities
        masses ({nparticle,} array_like): masses

    Returns:
        centered_coordinates ({nparticle, 3} ndarray): particle coordinates
        centered_velocities ({nparticle, 3} ndarray): particle velocities
    """
    center_coordinates = center(coordinates, masses)
    center_velocities = center(velocities, masses)
    mass_r = np.reshape(masses, (-1, 1))
    linear_momentum = center_velocities * mass_r
    linear_constraint_mat = np.hstack([np.eye(3) for i in range(np.shape(center_coordinates)[0])]) # 3 x (N x 3)
    angular_constraint_mat = np.reshape(np.transpose(utils.cross_product_matrix(center_coordinates), (1, 0, 2)), (3, -1)) # 3 x (N x 3)
    total_constraint_mat = np.vstack([linear_constraint_mat, angular_constraint_mat]) # 6 x (N x 3)
    constraint_vec = - np.dot(total_constraint_mat, np.reshape(linear_momentum, (-1, 1))) # 6 x 1
    dp = np.linalg.lstsq(total_constraint_mat, constraint_vec, rcond=None)[0]
    dv = np.reshape(dp, (-1, 3)) / mass_r
    center_velocities = center_velocities + dv
    return center_coordinates, center_velocities

def rescale_velocity(velocities, masses, kT):
    """
    Rescale velocities to correspond to a specific temperature.

    Args:
        velocities ({nparticle, ndim} ndarray): particle velocities
        masses ({nparticle,} ndarray): particle masses
        kT (float): Temperature in energy units

    Returns:
        rescaled_velocities ({nparticle, ndim} ndarray): particle velocities after scaling to target kT
    """
    alpha = np.sqrt(kT / utils.compute_temperature(velocities, masses))
    rescaled_velocities = velocities * alpha
    return rescaled_velocities

def boltzmann(kT, masses, dim=3, recenter=True):
    """
    Take sample in dim dimensions for each masses at each temperature.

    Args:
        kT (float): Temperature in energy units
        masses ({nparticle,} ndarray): particle masses
        dim (int): number of dimensions to sample
        recenter (bool): True to remove center of mass motion

    Returns:
        velocities ({nparticle, ndim} ndarray): particle velocities sampled from boltzmann distribution at energy kT
    """
    mass_r = np.reshape(masses, (-1, 1))
    factor = np.sqrt(kT/mass_r)
    velocities = np.random.normal(scale=factor, size=(mass_r.shape[0], dim))
    if recenter:
        velocities = center(velocities, mass_r)
    return velocities

def sin_velocity(masses, Qs, kT):
    """
    Samples velocities and auxiliary velocities from Boltzmann distribution and then rescales the joint velocity and auxiliary velocity to fit the isokinetic Nose Hoover constraint.

    Args:
        masses ({nparticle,} ndarray): masses of nonauxiliary dofs
        Qs ({2, L, nparticle, ndim} ndarray): auxiliary masses
        kT (float): temperature in energy units

    Returns:
        ({nparticle, ndim} ndarray, {2, L, nparticle, ndim} ndarray): tuple of velocity and auxiliary velocity arrays
    """
    _, Lint, natom, ndim = np.shape(Qs)
    factor = np.sqrt(kT/Qs)
    velocities = boltzmann(kT, masses, dim=ndim)
    aux_v = np.random.normal(scale=factor)
    L = float(Lint)
    for i in range(natom):
        for a in range(ndim):
            v = np.concatenate((velocities[i, a], aux_v[0, :, i, a]), axis=None)
            m = np.concatenate((masses[i] / L, 1.0/(L+1.0)*Qs[0, :, i, a]), axis=None)
            curr_kT = np.sum(v**2 * m)
            alpha = np.sqrt(kT / curr_kT)
            velocities[i, a] *= alpha
            aux_v[0, :, i, a] *= alpha

    Lhalf = Lint//2
    # The first aux velocity should not change sign during dynamics, so we
    # initialize half of them to each sign (if there are an even number of them)
    if Lhalf*2 == Lint:
        aux_v[0, :Lhalf, :, :] =  np.abs(aux_v[0, :Lhalf, :, :])
        aux_v[0, Lhalf:, :, :] = -np.abs(aux_v[0, Lhalf:, :, :])

    return velocities, aux_v

def boltzmann_range(kT, masses, ppf_range):
    """
    Take single sample of the Boltzmann distribution within a probability range.

    Args:
        kT (float): Temperature in energy units
        masses ({nparticle,} ndarray): particle masses
        ppf_range ({2, nparticle} array_like): percent range to sample in probability distribution

    Returns:
        speeds ({nparticle,} ndarray): speeds sampled from boltzmann distribution in given probability range
    """
    ppf = np.random.uniform(ppf_range[0], ppf_range[1])
    factor = np.sqrt(kT/masses)
    sample = scipy.stats.maxwell.ppf(ppf, scale=factor)
    return sample

def bimol(
        coordinates1,
        coordinates2,
        masses1,
        masses2,
        impact_kT=10000.0 * units.K_TO_AU,
        vibration_kT=600.0 * units.K_TO_AU,
        nosc=100,
        freq=3000.0 * units.INV_CM_TO_AU,
        impact_ppf_range=(0.4,0.6),
        tolerance=4.0,
    ):
    """
    Function to go from two molecular coordinates to a total configuration for bimolecular collision

    Args:
        coordinates1 ({nparticle1, 3} ndarray): coordinates of molecule 1
        coordinates2 ({nparticle2, 3} ndarray): coordinates of molecule 2
        masses1 ({nparticle1,} ndarray): masses of atoms in molecule 1
        masses2 ({nparticle2,} ndarray): masses of atoms in molecule 2
        impact_kT (float): temperature at which to sample the impact in energy units
        vibration_kT (float): temperature at which to sample molecular vibrations in energy units
        nosc (float): number of oscillations at given frequency, used to compute distance
        freq (float): vibrational frequency to use in energy units
        impact_ppf_range ({2,} array_like): percents to sample of Boltzmann cumulative distribution function for impact velocity
        tolerance (float): minimum usable distance between two molecules

    Returns:
        coordinates_tot ({nparticle1+nparticle2, 3} ndarray): coordinates of two molecules randomly rotated, separated on x axis
        velocities_tot ({nparticle1+nparticle2, 3} ndarray): velocities drawn from boltzmann distributions pointed at one another on x axis
    """
    # Convert to atomic units
    mass_tot = np.concatenate([masses1, masses2])
    natom1 = len(masses1)
    natom2 = len(masses2)

    # Get randomly rotated molecules
    coordinates1_rot = np.dot(center(coordinates1, masses1), utils.random_rotation_matrix())
    coordinates2_rot = np.dot(center(coordinates2, masses2), utils.random_rotation_matrix())

    # View each molecule as rigid body of their own mass at the given temperature,
    # sample boltzmann distribution of velocities within given percentage bounds
    sample_velocities1 = boltzmann_range(impact_kT, np.sum(masses1), impact_ppf_range)
    sample_velocities2 = boltzmann_range(impact_kT, np.sum(masses2), impact_ppf_range)
    speed_tot = sample_velocities1 + sample_velocities2

    # Compute time to impact based on number of oscillations at the given frequency
    impact_time = nosc / freq

    # Use velocities and time to impact to compute minimal distance
    impact_dist = speed_tot * impact_time

    # Place first molecule with center of mass zero, second molecule on 1d line with minimal distance
    min_dist = impact_dist + np.max(coordinates1_rot[:, 0]) - np.min(coordinates2_rot[:, 0])
    coordinates2_rot[:, 0] = coordinates2_rot[:, 0] + min_dist

    # Ensure molecules are beyond tolerance
    while np.min(utils.pairwise_dist(coordinates1_rot, coordinates2_rot)) < tolerance:
        coordinates2_rot[:, 0] += 0.2

    # Construct total coordinate matrix, center total mass at zero
    coordinates_tot = np.vstack([coordinates1_rot, coordinates2_rot])
    coordinates_tot = center(coordinates_tot, mass_tot)

    # Assign velocities
    velocities_tot = np.zeros_like(coordinates_tot)
    velocities_tot[:natom1, 0] = sample_velocities1
    velocities_tot[natom1:, 0] = -sample_velocities2

    # Add random boltzmann energy to each atom (may be different temperature)
    velocities_tot[:natom1, :] += boltzmann(vibration_kT, masses1)
    velocities_tot[natom1:, :] += boltzmann(vibration_kT, masses2)

    # Remove COM motion and total angular momentum
    coordinates_tot, velocities_tot = initialize_centered(coordinates_tot, velocities_tot, mass_tot)

    return coordinates_tot, velocities_tot

#def gaussian_x(x, a=1.0, x0=0.0, k0=1.0):
#    """
#    Gaussian wave packet of width a, centered at x0, with momentum k0
#    """
#    return ((a * np.sqrt(np.pi)) ** (-0.5)
#            * np.exp(-0.5 * ((x - x0) * 1. / a) ** 2 + 1j * x * k0))
#
#def gaussian_k(k, a=1.0, x0=0.0, k0=1.0):
#    """
#    Fourier transform of gaussian_x(x)
#    """
#    return ((a / np.sqrt(np.pi))**0.5
#            * np.exp(-0.5 * (a * (k - k0)) ** 2 - 1j * (k - k0) * x0))

max_stars_repo_name: jeford/mdprop

max_issues_repo_name: jeford/mdprop

max_forks_repo_name: jeford/mdprop

reason: import numpy,import scipy


text:
import logging
import os
import warnings
from datetime import datetime
import json
import copy
import pickle
import time

import numpy as np

logging.EXP = logging.INFO - 1
logging.TEST = 25

def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1, length = 100, fill = '█', printEnd = "\r"):
    """
    Call in a loop to create terminal progress bar
    @params:
        iteration   - Required  : current iteration (Int)
        total       - Required  : total iterations (Int)
        prefix      - Optional  : prefix string (Str)
        suffix      - Optional  : suffix string (Str)
        decimals    - Optional  : positive number of decimals in percent complete (Int)
        length      - Optional  : character length of bar (Int)
        fill        - Optional  : bar fill character (Str)
        printEnd    - Optional  : end character (e.g. "\r", "\r\n") (Str)
    """
    try:
        stage = printProgressBar.stage
    except:
        printProgressBar.stage = 0.
        stage = printProgressBar.stage
    ratio = iteration / float(total)
    if total > 10 and ratio < stage - 1.0e-5:
        return
    printProgressBar.stage += 0.1
    percent = ("{0:." + str(decimals) + "f}").format(100 * ratio)
    print(prefix, percent + "% Complete" + suffix)
    # Print New Line on Complete
    if iteration == total:
        print()

def log_init(log_prefix, level=logging.WARN):
    """
    If level<=EXP, Initialze the log file with the current time as ID.
    Otherwise only retrun ID
    """
    log = logging.getLogger()  # root logger
    for hdlr in log.handlers[:]:  # remove all old handlers
        log.removeHandler(hdlr)

    date = datetime.now()
    id = log_prefix + "-" + date.strftime('%Y%m%d-%H%M%S')
    if os.path.isfile("logging/" + id + '.log'):
        warnings.warn("ID repeatition.")
        id += "r"
    if level == logging.EXP:
        logging.basicConfig(filename="logging/" + id + '.log', level=level)
    else:
        logging.basicConfig(level=level)
    logging.info("ID: " + id)
    logging.info("Date: " + date.strftime('%Y-%m-%d %H:%M:%S') + "\n")
    return id

def log_params(parameters):
    paramstr = "Parameters:\n"
    for name, value in parameters.items():
        paramstr += name + " " + str(value) + "\n"
    logging.info(paramstr)

def log_finish(id, parameters, remark=""):
    """
    Save the parameters and
    """
    parameters = copy.deepcopy(parameters)
    parameters["ID"] = id
    date = datetime.now()
    parameters["finishing date"] = str(date.date())
    parameters["finishing time"] = str(date.time())
    parameters["remark"] = remark
    if logging.getLogger().getEffectiveLevel() == logging.EXP:
        logging.info("Saving parameters into record...")
        with open("logging/logging_record.json", "a") as write_file:
            try:
                json.dump(parameters, write_file)
                write_file.write("\n")
            except Exception as err:
                logging.exception(
                    "Failed to save heading data into the record.\n" +
                    str(parameters))

    logging.info("Simulation finished.")

def save_data(id, data):
    if logging.getLogger().getEffectiveLevel() != logging.EXP:
        return
    logging.info("Saving data...")
    try:
        data_address = "data/" + id + ".pickle"
        outfile = open(data_address, 'wb')
        pickle.dump(data, outfile)
        outfile.close()
        logging.info("Data saved at " + '"' + data_address + '"')
    except Exception as err:
        logging.exception("Failed to save data")

def load_data(id):
    infile = open("data//" + id + ".pickle", "rb")
    data = pickle.load(infile)
    infile.close()
    return data

def create_iter_kwargs(parameters):

    from itertools import product
    names = parameters.keys()
    values = [value if isinstance(value, list) else [value]
                for value in list(parameters.values())]

    kwarg_list = []
    for param_list in product(*values):
        param_dict = {}
        for i, key in enumerate(names):
            value = param_list[i]
            param_dict[key] = value
        kwarg_list.append(param_dict)
    return kwarg_list

def find_record_id(ID):
    patterns = {"ID":ID}
    result = find_record_patterns(patterns)
    if not result:
        return None
    else:
        return result

def find_record_patterns(patterns):
    with open("logging/logging_record.json") as log:
        record_list = log.readlines()
    matched_records = []
    if patterns is not None:
        for record in record_list:
            record = json.loads(record)
            match = True
            for key in patterns.keys():
                try:
                    if record[key] != patterns[key]:
                        match = False
                        continue
                except KeyError:
                    match = False
                    continue
            if match:
                matched_records.append(record)
    for parameters in matched_records:
        try:
            parameters["protocol"] = tuple(parameters["protocol"])
        except KeyError:
            pass
    if len(matched_records) == 1:
        return matched_records[0]
    else:
        return matched_records

def mytimeit(func):
    """
    Decorator for automatically calculate the elapsed time.
    """
    def inner(*arg, **kwargs):
        start_time = time.time()
        result = func(*arg, **kwargs)
        end_time = time.time()
        print("Elapse time: {}\n".format(end_time - start_time))
        return result
    return inner

def clear_junk_data():
    for filename in os.listdir("./logging"):
        if filename.endswith(".log"):
            id_end = filename.rfind('.')
            ID = filename[:id_end]
            result = find_record_id(ID)
            if result is None:
                try:
                    os.remove("./logging/" + ID + ".log")
                except:
                    pass
    for filename in os.listdir("./data"):
        if filename.endswith(".pickle"):
            id_end = filename.rfind('.')
            ID = filename[:id_end]
            result = find_record_id(ID)
            if result is None:
                try:
                    os.remove("./data/" + ID + ".pickle")
                except:
                    pass

max_stars_repo_name: BoxiLi/repeater-cut-off-optimization

max_issues_repo_name: BoxiLi/repeater-cut-off-optimization

max_forks_repo_name: BoxiLi/repeater-cut-off-optimization