text:
# Copyright (C) 2020 Zurich Instruments
#
# This software may be modified and distributed under the terms
# of the MIT license. See the LICENSE file for details.

import numpy as np
from typing import List

from zhinst.toolkit.control.drivers.base import BaseInstrument, AWGCore, ToolkitError
from zhinst.toolkit.control.node_tree import Parameter
from zhinst.toolkit.control.parsers import Parse
from zhinst.toolkit.interface import DeviceTypes
from zhinst.toolkit.helpers import SequenceType, TriggerMode

MAPPINGS = {
    "result_source": {
        0: "Crosstalk",
        1: "Threshold",
        2: "Rotation",
        4: "Crosstalk Correlation",
        5: "Threshold Correlation",
        7: "Integration",
    },
    "averaging_mode": {0: "Cyclic", 1: "Sequential",},
}

class UHFQA(BaseInstrument):
    """High-level driver for the Zurich Instruments UHFQA Quantum Analyzer.

    Inherits from :class:`BaseInstrument` and adds an :class:`AWGCore` and a
    list of :class:`ReadoutChannels`. They can be accessed as properties of the
    UHFQA.

        >>> import zhinst.toolkit as tk
        >>> ...
        >>> uhfqa = tk.UHFQA("uhfqa", "dev1111")
        >>> uhfqa.setup()
        >>> uhfqa.connect_device()
        >>> ...
        >>> uhfqa.nodetree
        <zhinst.toolkit.tools.node_tree.NodeTree object at 0x0000021091FFC898>
        nodes:
        - stats
        - osc
        - triggers
        - status
        - dio
        - auxin
        - scope
        - system
        - sigins
        - sigouts
        - features
        - auxouts
        - awg
        - qa
        parameters:
        - clockbase

    A measurement using the *AWG Core* could look like this:

        >>> uhfqa.awg.set_sequence_params(...)
        >>> uhfqa.channels[0].enable()
        >>> uhfqa.channels[0].readout_frequency(123e6)
        >>> ...
        >>> uhfqa.arm(length=100, averages=1)   # arm the QA Readout, set length and averages
        >>> uhfqa.awg.run()                     # start the AWG
        >>> uhfqa.awg.wait_done()               # wait until AWG is finished
        >>> result = uhfqa.channels[0].result() # get the result vector

    Arguments:
        name (str): Identifier for the UHFQA.
        serial (str): Serial number of the device, e.g. *'dev1234'*. The serial
            number can be found on the back panel of the instrument.
        discovery: an instance of ziDiscovery

    Attributes:
        awg (:class:`zhinst.toolkit.control.drivers.uhfqa.AWG`):
            A device-specific :class:`AWGCore` for the UHFQA.
        channels (list): A list of ten
            :class:`zhinst.toolkit.control.drivers.uhfqa.ReadoutChannel` s that
            each represent the digital signal processing path on the instrument.
        integration_time (:class:`zhinst.toolkit.control.node_tree.Parameter`):
            The time in seconds used for signal integration. The value must be
            positive. The maximum value when using weighted integration is 4096
            samples or ca. 2.275 us (default: 2.0 us).
        result_source (:class:`zhinst.toolkit.control.node_tree.Parameter`):
            This parameter selects the stage in the signal processing path that
            is used as the source for the QA results. It can be one of
            {`"Crosstalk"`, `"Threshold"`, `"Rotation"`, `"Crosstalk
            Correlation"`, `"Threshold Correlation"`, `"Integration"`}.
        averaging_mode (:class:`zhinst.toolkit.control.node_tree.Parameter`):
            Mode of averaing in the QA Result Acquisition. Either `"Sequential"`
            or `"Cyclic"`.
            `"Sequential"`: The first point of the Result vector is the average
            of the first *N* results where *N* is the value of the *result
            averages* setting. The second point is the average of the next *N*
            results, and so forth.
            `"Cyclic"`: The first point in the Result vector is the average of
            the results *1, M+1, 2M+1, ...*, where *M* is the value of the
            *result length* setting. The second point is the average of the
            results number *2, M+2, 2M+2, ...*, and so forth.

    """

    def __init__(self, name: str, serial: str, discovery=None, **kwargs) -> None:
        super().__init__(name, DeviceTypes.UHFQA, serial, discovery, **kwargs)
        self._awg = AWG(self, 0)
        self._channels = [ReadoutChannel(self, i) for i in range(10)]
        self.integration_time = Parameter(
            self,
            dict(
                Node="qas/0/integration/length",
                Description="The integration time of the QA Integration unit.",
                Type="Double",
                Properties="Read, Write",
                Unit="s",
            ),
            device=self,
            set_parser=Parse.uhfqa_time2samples,
            get_parser=Parse.uhfqa_samples2time,
        )
        self.result_source = Parameter(
            self,
            dict(
                Node="qas/0/result/source",
                Description="The signal source for the QA Results unit. One of {'Crosstalk', 'Threshold', 'Rotation', 'Crosstalk Correlation', 'Threshold Correlation', 'Integration'}.",
                Type="Integer",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self,
            mapping=MAPPINGS["result_source"],
        )
        self.averaging_mode = Parameter(
            self,
            dict(
                Node="qas/0/result/mode",
                Description="Selects the order of the result logger. One of {'Cyclic', 'Sequential'}.",
                Type="Integer",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self,
            mapping=MAPPINGS["averaging_mode"],
        )

    def connect_device(self, nodetree: bool = True) -> None:
        """Connects the device to the data server and initializes the AWG.

        Keyword Arguments:
            nodetree (bool): A flag that specifies if all the parameters from
                the device's nodetree should be added to the object's attributes
                as `zhinst-toolkit` Parameters. (default: True)

        """
        super().connect_device(nodetree=nodetree)
        self.awg._setup()

    def crosstalk_matrix(self, matrix=None):
        """Sets or gets the crosstalk matrix of the UHFQA as a 2D array.

        Keyword Arguments:
            matrix (2D array): The 2D matrix used in the digital signal
                processing path to compensate for crosstalk between the
                different channels. The given matrix can also be a part of the
                entire 10 x 10 matrix. Its maximum dimensions are 10 x 10.
                (default: None)

        Returns:
            If no argument is given the method returns the current crosstalk
            matrix as a 2D numpy array.

        """
        if matrix is None:
            m = np.zeros((10, 10))
            for r in range(10):
                for c in range(10):
                    m[r, c] = self._get(f"qas/0/crosstalk/rows/{r}/cols/{c}")
            return m
        else:
            rows, cols = matrix.shape
            if rows > 10 or cols > 10:
                raise ValueError(
                    f"The shape of the given matrix is {rows} x {cols}. The maximum size is 10 x 10."
                )
            for r in range(rows):
                for c in range(cols):
                    self._set(f"qas/0/crosstalk/rows/{r}/cols/{c}", matrix[r, c])

    def enable_readout_channels(self, channels: List = range(10)) -> None:
        """Enables weighted integration on the specified readout channels.

        Keyword Arguments:
            channels (list): A list of indices of channels to enable.
                (default: range(10))

        """
        for i in channels:
            if i not in range(10):
                raise ValueError(f"The channel index {i} is out of range!")
            self.channels[i].enable()

    def disable_readout_channels(self, channels: List = range(10)) -> None:
        """Disables weighted integration on the specified readout channels.

        Keyword Arguments:
            channels (list): A list of indices of channels to disable.
                (default: range(10))

        """
        for i in channels:
            if i not in range(10):
                raise ValueError(f"The channel index {i} is out of range!")
            self.channels[i].disable()

    def arm(self, length=None, averages=None) -> None:
        """Prepare UHFQA for result acquisition.

        This method enables the QA Results Acquisition and resets the acquired
        points. Optionally, the *result length* and *result averages* can be set
        when specified as keyword arguments. If they are not specified, they are
        not changed.

        Keyword Arguments:
            length (int): If specified, the length of the result vector will be
                set before arming the UHFQA readout. (default: None)
            averages (int): If specified, the result averages will be set before
                arming the UHFQA readout. (default: None)

        """
        if length is not None:
            self._set("qas/0/result/length", int(length))
        if averages is not None:
            self._set("qas/0/result/averages", int(averages))
        self._set("qas/0/result/enable", 1)
        # toggle node value from 0 to 1 for reset
        self._set("qas/0/result/reset", 0)
        self._set("qas/0/result/reset", 1)

    def _init_settings(self):
        settings = [
            ("awgs/0/single", 1),
        ]
        self._set(settings)

    @property
    def awg(self):
        return self._awg

    @property
    def channels(self):
        return self._channels

class AWG(AWGCore):
    """Device-specific AWG Core for UHFQA.

    Inherits from `AWGCore` and adds :mod:`zhinst-toolkit` :class:`Parameters`
    like ouput or gains. This class also specifies sequence specific settings
    for the UHFQA.

    Attributes:
        output1 (:class:`zhinst.toolkit.control.node_tree.Parameter`): The state
            of the output of channel 1. Can be one of {'on', 'off'}.
        output2 (:class:`zhinst.toolkit.control.node_tree.Parameter`): The state
            of the output of channel 2. Can be one of {'on', 'off'}.
        gain1 (:class:`zhinst.toolkit.control.node_tree.Parameter`): Gain of the
            output channel 1. The value must be between -1 and +1 (default: +1).
        gain2 (:class:`zhinst.toolkit.control.node_tree.Parameter`): Gain of the
            output channel 2. The value must be between -1 and +1 (default: +1).

    """

    def __init__(self, parent: BaseInstrument, index: int) -> None:
        super().__init__(parent, index)
        self.output1 = Parameter(
            self,
            dict(
                Node="sigouts/0/on",
                Description="Enables or disables both ouputs of the AWG. Either can be {'1', '0'} or {'on', 'off'}.",
                Type="Integer",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self._parent,
            set_parser=Parse.set_on_off,
            get_parser=Parse.get_on_off,
        )
        self.output2 = Parameter(
            self,
            dict(
                Node="sigouts/1/on",
                Description="Enables or disables both ouputs of the AWG. Either can be {'1', '0'} or {'on', 'off'}.",
                Type="Integer",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self._parent,
            set_parser=Parse.set_on_off,
            get_parser=Parse.get_on_off,
        )
        self.gain1 = Parameter(
            self,
            dict(
                Node="awgs/0/outputs/0/amplitude",
                Description="Sets the gain of the first output channel.",
                Type="Double",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self._parent,
            set_parser=Parse.amp1,
        )
        self.gain2 = Parameter(
            self,
            dict(
                Node="awgs/0/outputs/1/amplitude",
                Description="Sets the gain of the second output channel.",
                Type="Double",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self._parent,
            set_parser=Parse.amp1,
        )

    def _apply_sequence_settings(self, **kwargs):
        if "sequence_type" in kwargs.keys():
            t = SequenceType(kwargs["sequence_type"])
            allowed_sequences = [
                SequenceType.NONE,
                SequenceType.SIMPLE,
                SequenceType.READOUT,
                SequenceType.CW_SPEC,
                SequenceType.PULSED_SPEC,
                SequenceType.CUSTOM,
            ]
            if t not in allowed_sequences:
                raise ToolkitError(
                    f"Sequence type {t} must be one of {[s.value for s in allowed_sequences]}!"
                )
            # apply settings depending on sequence type
            elif t == SequenceType.CW_SPEC:
                self._apply_cw_settings()
            elif t == SequenceType.PULSED_SPEC:
                self._apply_pulsed_settings()
            elif t == SequenceType.READOUT:
                self._apply_readout_settings()
            else:
                self._apply_base_settings()
        # apply settings dependent on trigger type
        if "trigger_mode" in kwargs.keys():
            if TriggerMode(kwargs["trigger_mode"]) == TriggerMode.EXTERNAL_TRIGGER:
                self._apply_trigger_settings()

    def _apply_base_settings(self):
        settings = [
            ("sigouts/*/enables/*", 0),
            ("sigouts/*/amplitudes/*", 0),
            ("awgs/0/outputs/*/mode", 0),
            ("qas/0/integration/mode", 0),
        ]
        self._parent._set(settings)

    def _apply_cw_settings(self):
        settings = [
            ("sigouts/0/enables/0", 1),
            ("sigouts/1/enables/1", 1),
            ("sigouts/0/amplitudes/0", 1),
            ("sigouts/1/amplitudes/1", 1),
            ("awgs/0/outputs/*/mode", 0),
            ("qas/0/integration/mode", 1),
        ]
        self._parent._set(settings)
        self._parent.disable_readout_channels(range(10))

    def _apply_pulsed_settings(self):
        settings = [
            ("sigouts/*/enables/*", 0),
            ("sigouts/*/amplitudes/*", 0),
            ("awgs/0/outputs/*/mode", 1),
            ("qas/0/integration/mode", 1),
        ]
        self._parent._set(settings)
        self._parent.disable_readout_channels(range(10))

    def _apply_readout_settings(self):
        settings = [
            ("sigouts/*/enables/*", 0),
            ("awgs/0/outputs/*/mode", 0),
            ("qas/0/integration/mode", 0),
        ]
        self._parent._set(settings)

    def _apply_trigger_settings(self):
        settings = [
            ("/awgs/0/auxtriggers/*/channel", 0),
            ("/awgs/0/auxtriggers/*/slope", 1),
        ]
        self._parent._set(settings)

    def outputs(self, value=None):
        """Sets both signal outputs simultaneously.

        Keyword Arguments:
            value (tuple): Tuple of values {'on', 'off'} for channel 1 and 2
                (default: {None})

        Returns:
            The state {'on', 'off'} for both outputs if the keyword argument is
            not given.

        """
        if value is None:
            return self.output1(), self.output2()
        else:
            if isinstance(value, tuple) or isinstance(value, list):
                if len(value) != 2:
                    raise ToolkitError(
                        "The values should be specified as a tuple, e.g. ('on', 'off')."
                    )
                return self.output1(value[0]), self.output2(value[1])
            else:
                raise ToolkitError("The value must be a tuple or list of length 2!")

    def update_readout_params(self) -> None:
        """Updates the sequence parameters for 'Simple' sequence with values from the readout channels."""
        if self.sequence_params["sequence_type"] == SequenceType.READOUT:
            freqs = []
            amps = []
            phases = []
            for ch in self._parent.channels:
                if ch.enabled():
                    freqs.append(ch.readout_frequency())
                    amps.append(ch.readout_amplitude())
                    phases.append(ch.phase_shift())
            self.set_sequence_params(
                readout_frequencies=freqs, readout_amplitudes=amps, phase_shifts=phases,
            )
        else:
            raise ToolkitError("AWG Sequence type needs to be 'Readout'")

    def compile(self) -> None:
        """Wraps the 'compile(...)' method of the parent class `AWGCore`."""
        if self.sequence_params["sequence_type"] == SequenceType.READOUT:
            self.update_readout_params()
        super().compile()

class ReadoutChannel:
    """Implements a Readout Channel for UHFQA.

    This class represents the signal processing chain for one of the ten
    :class:`ReadoutChannels` of a UHFQA. One channel is typically used for
    dispersive resonator readout of superconducting qubits.

        >>> ch = uhfqa.channels[0]
        >>> uhfqa.result_source("Threshold")
        >>> ...
        >>> ch.enable()
        >>> ch.readout_frequency(85.6e6)
        >>> ch
        Readout Channel 0:  <zhinst.toolkit.uhfqa.ReadoutChannel object at 0x0000021091FD4F28>
            rotation          : 0.0
            threshold         : 0.0
            readout_frequency : 85600000.0
            readout_amplitude : 1
            phase_shift       : 0
        >>> ch.rotation(123.4)
        >>> ch.threshold(-56.78)
        >>> ...
        >>> ch.result()
        array([0.0, 1.0, 1.0, 1.0, 0.0, ...])

    The readout channel can be enabled with `enable()` which means that the
    weighted integration mode is activated and integration weights are set to
    demodulate the signal at the given readout frequency. If the channel is
    enabled, the readout parameters are also used for signal generation in the
    :class:`AWGCore` if the sequence type is set to *'Readout'*.

    Attributes:
        index (int): The index of the Readout Channel from 1 - 10.
        rotation (:class:`zhinst.toolkit.control.nodetree.Parameter`): The
            rotation applied to the signal in IQ plane. The angle is specified
            in degrees.
        threshold (:class:`zhinst.toolkit.control.nodetree.Parameter`): The
            signal threshold used for state discrimination in the thresholding
            unit.
        result (:class:`zhinst.toolkit.control.nodetree.Parameter`): This
            read-only Parameter holds the result vector for the given readout
            channel as a 1D numpy array.

    """

    def __init__(self, parent: BaseInstrument, index: int) -> None:
        if index not in range(10):
            raise ValueError(f"The channel index {index} is out of range!")
        self._index = index
        self._parent = parent
        self._enabled = False
        self._readout_frequency = 100e6
        self._readout_amplitude = 1
        self._phase_shift = 0
        self.rotation = Parameter(
            self,
            dict(
                Node=f"qas/0/rotations/{self._index}",
                Description="Sets the rotation of the readout channel in degrees.",
                Type="Double",
                Properties="Read, Write",
                Unit="Degrees",
            ),
            device=self._parent,
            set_parser=Parse.deg2complex,
            get_parser=Parse.complex2deg,
        )
        self.threshold = Parameter(
            self,
            dict(
                Node=f"qas/0/thresholds/{self._index}/level",
                Description="Sets the threshold of the readout channel.",
                Type="Double",
                Properties="Read, Write",
                Unit="None",
            ),
            device=self._parent,
        )
        self.result = Parameter(
            self,
            dict(
                Node=f"qas/0/result/data/{self._index}/wave",
                Description="Returns the result vector of the readout channel.",
                Type="Numpy array",
                Properties="Read",
                Unit="None",
            ),
            device=self._parent,
            get_parser=self._average_result,
        )

    @property
    def index(self):
        return self._index

    def enabled(self) -> bool:
        """Returns if weighted integration is enabled."""
        return self._enabled

    def enable(self) -> None:
        """Enables weighted integration for this channel."""
        self._enabled = True
        self._parent._set("qas/0/integration/mode", 0)
        self._parent.integration_time(2e-6)
        self._set_int_weights()

    def disable(self) -> None:
        """Disables weighted integration for this channel."""
        self._enabled = False
        self._reset_int_weights()

    def readout_frequency(self, freq=None):
        """Sets or gets the readout frequency for this channel.

        Readout frequency in Hz of the readout channel. If the AWG
        :class:`SequenceProgram` is of type "Readout", this Parameter is used to
        generate a readout tone at the given readout frequency for all readout
        channels that are enabled. This frequency is also used in the signal
        acquisition for digital demodulation if the readout channel is
        enabled. The frequency must be positive.

        """
        if freq is None:
            return self._readout_frequency
        else:
            Parse.greater0(freq)
            self._readout_frequency = freq
            self._set_int_weights()
            return self._readout_frequency

    def readout_amplitude(self, amp=None):
        """Sets or gets the readout amplitude for this channel.

        The amplitude of the readout pulse is used for signal generation of the
        readout tone if the channel is enabled and if the AWG
        :class:`SequenceProgram` is of type *'Readout'*. (default: 1.0)

        """
        if amp is None:
            return self._readout_amplitude
        else:
            Parse.amp1(amp)
            self._readout_amplitude = amp
            return self._readout_amplitude

    def phase_shift(self, ph=None):
        """Sets or gets the readout phase shift for this channel.

        Additional phase shift in the signal generation between I and Q
        quadtratures. (default: 0)

        """
        if ph is None:
            return self._phase_shift
        else:
            self._phase_shift = Parse.phase(ph)
            return self._phase_shift

    def _reset_int_weights(self):
        node = f"/qas/0/integration/weights/"
        self._parent._set(node + f"{self._index}/real", np.zeros(4096))
        self._parent._set(node + f"{self._index}/imag", np.zeros(4096))

    def _set_int_weights(self):
        self._reset_int_weights()
        l = self._parent._get("qas/0/integration/length")
        freq = self.readout_frequency()
        node = f"/qas/0/integration/weights/{self._index}/"
        self._parent._set(node + "real", self._demod_weights(l, freq, 0))
        self._parent._set(node + "imag", self._demod_weights(l, freq, 90))

    def _demod_weights(self, length, freq, phase):
        if length > 4096:
            raise ValueError(
                "The maximum length of the integration weights is 4096 samples."
            )
        if freq <= 0:
            raise ValueError("This frequency must be positive.")
        clk_rate = 1.8e9
        x = np.arange(0, length)
        y = np.sin(2 * np.pi * freq * x / clk_rate + np.deg2rad(phase))
        return y

    def _average_result(self, result):
        n_averages = self._parent._get("qas/0/result/averages")
        return result / n_averages

    def __repr__(self):
        s = f"Readout Channel {self._index}:  {super().__repr__()}\n"
        s += f"     rotation          : {self.rotation()}\n"
        s += f"     threshold         : {self.threshold()}\n"
        if self._enabled:
            s += f"     readout_frequency : {self.readout_frequency()}\n"
            s += f"     readout_amplitude : {self.readout_amplitude()}\n"
            s += f"     phase_shift       : {self.phase_shift()}\n"
        else:
            s += "      Weighted Integration Disabled\n"
        return s

max_stars_repo_name: ovaar/zhinst-toolkit

max_issues_repo_name: ovaar/zhinst-toolkit

max_forks_repo_name: ovaar/zhinst-toolkit


text:
#!/usr/bin/env python

import os
import cPickle as pickle
from scipy import integrate

def reduce_to_day(df):

    # lambda function for integrated sum
    dayint = lambda x: integrate.trapz(x, dx=1800)*1e-6

    # columns 4-5 are conductance which isn't summed
    up_samps = {clab: 'sum' \
                if i in [2, 9] else dayint \
                if i in [3, 10] else 'mean' \
                for (i, clab) in enumerate(df.columns)}

    # resample based on the up-scaling dict above
    daily_df = df.bfill().resample('D', how=up_samps)

    return daily_df

def main():

    metdata = pickle.load(open(PKLPATH + "hourly/tower_inputs.pkl", 'rb'))

    # daily aggregation here
    daily_met = reduce_to_day(metdata)

    pickle.dump(daily_met, open(PKLPATH + "daily/daily_inputs.pkl", 'wb'))

    return 1

if __name__ == '__main__':

    PKLPATH = os.path.expanduser("~/Savanna/Data/HowardSprings_IAV/pickled/")

    main()

ext: py

reason: from scipy


text:
import requests
import pickle
import gzip
import math
import numpy as np
import torch
import torch.nn.functional as F
from pathlib import Path
from torch import nn, optim
from torch.cuda.amp import autocast
from torch.utils.tensorboard import SummaryWriter
from torch.utils.data import TensorDataset, DataLoader

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

# For reproducibility
torch.manual_seed(0)
torch.backends.cudnn.deterministic = True   # may be slower
torch.backends.cudnn.benchmark = False
np.random.seed(0)

def basic():
    x = torch.tensor([5.5, 3])
    x = x.new_ones((5, 3), dtype=torch.double)
    y = torch.rand(5, 3)
    y += x
    print(y.view(-1, 5).size())     # Resizing

    # The Torch Tensor and NumPy array will share their underlying memory locations
    # (if the Torch Tensor is on CPU), and changing one will change the other
    a = y[: 1]
    b = a.numpy()
    print(a, b)
    val = torch.randn(1).item()     # If you have a one element tensor, use .item() to get the value as a Python number
    a.add_(val)     # in-place
    print(a, b)
    b = torch.from_numpy(b)
    print(b)
    b = b.to(device)
    print(b, b.to("cpu"))

def autograd():
    # Tensor and Function are interconnected and build up an acyclic graph, that encodes a complete
    # history of computation. Each tensor has a .grad_fn attribute that references a Function that has
    # created the Tensor (except for Tensors created by the user - their grad_fn is None).
    x = torch.ones(2, 2, requires_grad=True)    # if requires_grad = True, it starts to track all operations on it
    print(x)
    y = x + 2
    z = y * y * 3
    out = z.mean()
    print(out)
    print(out.grad_fn.next_functions[0][0])
    out.backward()  # all the gradients computed automatically
    print(x.grad)   # d(out)/dx

    # If Tensor isn't a scalar, you need to specify a gradient argument that is a tensor of matching shape
    x = torch.randn(3, requires_grad=True)
    y = x * 2
    while y.data.norm() < 1000:
        y = y * 2
    print(y)
    v = torch.tensor([0.1, 1.0, 0.0001], dtype=torch.float)
    y.backward(v)
    print(x.grad)

    # Stop autograd from tracking history on Tensors
    # either by wrapping the code block in "with torch.no_grad():" or by using .detach()
    with torch.no_grad():
        print((x ** 2).requires_grad)
    x = x.detach()
    print(x.requires_grad)

def neural_networks():
    class MyReLU(torch.autograd.Function):
        @staticmethod
        def forward(ctx, input):
            ctx.save_for_backward(input)  # cache arbitrary objects for use in the backward pass
            return input.clamp(min=0)

        @staticmethod
        def backward(ctx, grad_output):
            input, = ctx.saved_tensors
            grad_input = grad_output.clone()
            grad_input[input < 0] = 0
            return grad_input

    class Net(nn.Module):
        def __init__(self):
            super(Net, self).__init__()
            self.conv1 = nn.Conv2d(1, 6, 3)
            self.pool = nn.MaxPool2d(2, 2)
            self.conv2 = nn.Conv2d(6, 16, 5)
            self.fc1 = nn.Linear(16 * 5 * 5, 120)
            self.fc2 = nn.Linear(120, 84)
            self.fc3 = nn.Linear(84, 10)

        def forward(self, x):
            # You just have to define the forward function, and the backward function
            # is automatically defined for you using autograd
            x = self.pool(MyReLU.apply(self.conv1(x)))
            x = self.pool(MyReLU.apply(self.conv2(x)))
            x = x.view(-1, 16 * 5 * 5)
            x = F.relu(self.fc1(x))
            x = F.relu(self.fc2(x))
            x = self.fc3(x)
            return x

    net = Net()
    net.to(device)
    criterion = nn.CrossEntropyLoss()
    optimizer = optim.SGD(net.parameters(), lr=0.001, momentum=0.9)

    n, bs = 12000, 2000
    data = torch.randn(n, 1, 32, 32).to(device)
    target = torch.randint(10, (n,)).to(device)  # a dummy target
    classes = list('abcdefghij')

    writer = SummaryWriter()
    data0, target0 = data[:10], target[:10]
    class_target0 = [classes[lab] for lab in target0]
    features = data0.view(-1, 32 * 32)
    writer.add_graph(net, data0)
    writer.add_embedding(features, metadata=class_target0)

    for epoch in range(4):
        running_loss = 0.0
        for i in range(0, n, bs):
            inputs, labels = data[i:i+bs], target[i:i+bs]
            optimizer.zero_grad()   # clear the existing gradients though, else gradients will be accumulated
            with autocast():
                outputs = net(inputs)
                loss = criterion(outputs, labels)
            loss.backward()
            optimizer.step()    # does the update

            running_loss += loss.item()
            if i % bs == 0:
                global_step = epoch*n+i+bs
                for name, param in net.named_parameters():
                    writer.add_histogram(name, param.data, global_step)
                    writer.add_histogram(name+'_grad', param.grad, global_step)
                writer.add_scalar('training loss', running_loss/bs, global_step)
                running_loss = 0.0

    print('conv1.bias.grad', net.conv1.bias.grad)
    print('Finished Training')
    writer.close()

    PATH = 'tmp/cnn.pth'
    torch.save(net.state_dict(), PATH)
    net = Net()
    net.load_state_dict(torch.load(PATH))

def mnist():
    class MnistLogistic(nn.Module):
        def __init__(self):
            super().__init__()
            self.weights = nn.Parameter(torch.randn(784, 10) / math.sqrt(784))
            self.bias = nn.Parameter(torch.zeros(10))

        def forward(self, xb):
            return xb @ self.weights + self.bias  # the @ stands for the dot product operation

    class Lambda(nn.Module):
        def __init__(self, func):
            super().__init__()
            self.func = func

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

    class MnistCNN(nn.Module):
        def __init__(self):
            super().__init__()
            self.seq = nn.Sequential(
                nn.Conv2d(1, 16, kernel_size=3, stride=2, padding=1),
                nn.ReLU(),
                nn.Conv2d(16, 16, kernel_size=3, stride=2, padding=1),
                nn.ReLU(),
                nn.Conv2d(16, 10, kernel_size=3, stride=2, padding=1),
                nn.ReLU(),
                nn.AdaptiveAvgPool2d(1),
                Lambda(lambda x: x.view(x.size(0), -1)),
            )

        def forward(self, xb):
            return self.seq(xb)

    class WrappedDataLoader:
        def __init__(self, dl, func):
            self.dl = dl
            self.func = func

        def __len__(self):
            return len(self.dl)

        def __iter__(self):
            batches = iter(self.dl)
            for b in batches:
                yield self.func(*b)

    DATA_PATH = Path("tmp")
    PATH = DATA_PATH / "mnist"
    PATH.mkdir(parents=True, exist_ok=True)
    URL = "http://deeplearning.net/data/mnist/"
    FILENAME = "mnist.pkl.gz"
    if not (PATH / FILENAME).exists():
        content = requests.get(URL + FILENAME).content
        (PATH / FILENAME).open("wb").write(content)
    with gzip.open((PATH / FILENAME).as_posix(), "rb") as f:
        ((x_train, y_train), (x_valid, y_valid), _) = pickle.load(f, encoding="latin-1")
    x_train, y_train, x_valid, y_valid = map(torch.tensor, (x_train, y_train, x_valid, y_valid))
    train_ds, valid_ds = TensorDataset(x_train, y_train), TensorDataset(x_valid, y_valid)

    def get_data(train_ds, valid_ds, bs):
        return DataLoader(train_ds, batch_size=bs, shuffle=True), DataLoader(valid_ds, batch_size=bs * 2)

    def get_model():
        # model = MnistLogistic()
        model = MnistCNN()
        return model, optim.SGD(model.parameters(), lr=lr, momentum=0.9)

    def loss_batch(model, loss_func, xb, yb, opt=None):
        loss = loss_func(model(xb), yb)
        if opt is not None:
            loss.backward()
            opt.step()
            opt.zero_grad()
        return loss.item(), len(xb)

    def fit(epochs, model, loss_func, opt, train_dl, valid_dl):
        for epoch in range(epochs):
            model.train()
            for xb, yb in train_dl:
                loss_batch(model, loss_func, xb, yb, opt)
            model.eval()
            with torch.no_grad():
                losses, nums = zip(*[loss_batch(model, loss_func, xb, yb) for xb, yb in valid_dl])
            val_loss = np.sum(np.multiply(losses, nums)) / np.sum(nums)
            print(epoch, val_loss / len(valid_dl))

    def preprocess(x, y):
        return x.view(-1, 1, 28, 28).to(device), y.to(device)

    bs, lr, epochs = 64, 0.1, 2
    loss_func = F.cross_entropy

    train_dl, valid_dl = get_data(train_ds, valid_ds, bs)
    train_dl, valid_dl = WrappedDataLoader(train_dl, preprocess), WrappedDataLoader(valid_dl, preprocess)
    model, opt = get_model()
    model.to(device)
    fit(epochs, model, loss_func, opt, train_dl, valid_dl)

if __name__ == "__main__":
    # basic()
    # autograd()
    # neural_networks()
    mnist()

reason: import numpy


text:
import numpy as np

JND = 2.3

M_rgb_to_xyz = np.array([[0.4124564, 0.3575761, 0.1804375],
                         [0.2126729, 0.7151522, 0.0721750],
                         [0.0193339, 0.1191920, 0.9503041]])

M_xyz_to_rgb = np.linalg.inv(M_rgb_to_xyz)

# https://en.wikipedia.org/wiki/LMS_color_space#Hunt.2C_RLAB
M_xyz_d65_to_lms = np.array([[ 0.4002, 0.7076, -0.0808],
                             [-0.2263, 1.1653,  0.0457],
                             [ 0     , 0     ,  0.9182]])

M_lms_to_xyz_d65 = np.linalg.inv(M_xyz_d65_to_lms)

# https://en.wikipedia.org/wiki/LMS_color_space#Hunt.2C_RLAB
M_xyz_e_to_lms = np.array([[ 0.38971, 0.68898, -0.07868],
                           [-0.22981, 1.18340,  0.04641],
                           [ 0      , 0      ,  1.     ]])

M_lms_to_xyz_e = np.linalg.inv(M_xyz_e_to_lms)

# http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html - Bradford
M_xyz_e_to_xyz_d55 = np.array([[ 0.9797934, -0.0133624, -0.0096110],
                               [-0.0216865,  1.0235002, -0.0018137],
                               [-0.0037969,  0.0075580,  0.9177289]])

M_xyz_e_to_xyz_d65 = np.array([[ 0.9531874, -0.0265906,  0.0238731],
                               [-0.0382467,  1.0288406,  0.0094060],
                               [ 0.0026068, -0.0030332,  1.0892565]])

M_xyz_e_to_xyz_d75 = np.array([[ 0.9344831, -0.0355691,  0.0508059],
                               [-0.0490066,  1.0307122,  0.0182943],
                               [ 0.0079472, -0.0119897,  1.2304226]])

red_RGB = np.array([1,0,0])
blue_RGB = np.array([0,0,1])
white_RGB = np.array([1,1,1])

red_XYZ = M_rgb_to_xyz @ red_RGB
blue_XYZ = M_rgb_to_xyz @ blue_RGB
white_XYZ = M_rgb_to_xyz @ white_RGB

red_LMS = M_xyz_d65_to_lms @ red_XYZ
blue_LMS = M_xyz_d65_to_lms @ blue_XYZ
white_LMS = M_xyz_d65_to_lms @ white_XYZ

k = np.linspace(0, 1, 101)
zeros = np.zeros_like(k)
ones = np.ones_like(k)

# L-weak
q2p = (1 - k) * (blue_LMS[1] - blue_LMS[0]) / (blue_LMS[1] - blue_LMS[2])
q1p = 1 - k - q2p
SL = np.stack([k, q1p, q2p, zeros, ones, zeros, zeros, zeros, ones], axis=1).reshape(-1, 3, 3)

# M-weak
q2d = (1 - k) * (blue_LMS[0] - blue_LMS[1]) / (blue_LMS[0] - blue_LMS[2])
q1d = 1 - k - q2d
SM = np.stack([ones, zeros, zeros, q1d, k, q2d, zeros, zeros, ones], axis=1).reshape(-1, 3, 3)

# S-weak
q2t = (1 - k) * (red_LMS[0] - red_LMS[2]) / (red_LMS[0] - red_LMS[1])
q1t = 1 - k - q2t
SS = np.stack([ones, zeros, zeros, zeros, ones, zeros, q1t, q2t, k], axis=1).reshape(-1, 3, 3)

# LMS-weak combined
# From regular LMS to deficient LMS
SLMS = np.stack([SL, SM, SS]) # Axes: deficiency (lms), a, i, j

# Cone monochromats
S_mono_cone = np.array([[0, 0, 1],
                        [0, 0, 1],
                        [0, 0, 1]])

# Convert CIE XYZ to CIE L*a*b*
# https://en.wikipedia.org/wiki/CIELAB_color_space#Forward_transformation
def XYZ_to_Lab(XYZ, Xn=1., Yn=1., Zn=1.):
    @np.vectorize
    def f(t):
        delta = 6/29
        f = t**(1/3) if t > delta**3 else (t/(3*delta**2) + 4/29)
        return f

    X, Y, Z = XYZ[...,0], XYZ[...,1], XYZ[...,2]

    Lstar = 116 * f(Y/Yn) - 16
    astar = 500 * (f(X/Xn) - f(Y/Yn))
    bstar = 200 * (f(Y/Yn) - f(Z/Zn))

    Lab = np.stack([Lstar, astar, bstar], axis=-1)
    return Lab

# Calculate CIE Delta E00 distance
# doi 10.1002/col.20070
@np.vectorize
def hue(ap, bs):
    return 0. if ap == bs == 0. else np.rad2deg(np.arctan2(bs, ap)) % 360

@np.vectorize
def delta_hue(h1, h2, Cprime1, Cprime2):
    if Cprime1 * Cprime2 == 0:
        dh = 0
    elif np.abs(h2 - h1) <= 180:
        dh = h2 - h1
    elif h2 - h1 > 180:
        dh = h2 - h1 - 360
    else:
        dh = h2 - h1 + 360
    return dh

@np.vectorize
def hprimebar(hprime1, hprime2, Cprime1, Cprime2):
    if Cprime1 * Cprime2 == 0:
        hpb = hprime1 + hprime2
    elif np.abs(hprime1 - hprime2) <= 180:
        hpb = (hprime1 + hprime2) / 2
    elif np.abs(hprime1 - hprime2) > 180 and hprime1 + hprime2 < 360:
        hpb = (hprime1 + hprime2 + 360) / 2
    else:
        hpb = (hprime1 + hprime2 - 360) / 2
    return hpb

@np.vectorize
def dE00(L1, a1, b1, L2, a2, b2, kL=1., kC=1., kH=1.):
    L, a, b = np.array([L1, L2]), np.array([a1, a2]), np.array([b1, b2])

    C = np.sqrt(a**2 + b**2)
    Cbar = C.mean()
    G = 0.5 * (1 - np.sqrt(Cbar**7 / (Cbar**7 + 25**7)))
    aprime = (1 + G) * a
    Cprime = np.sqrt(aprime**2 + b**2)

    h = hue(aprime, b)

    dL = L2 - L1
    dCprime = Cprime[1] - Cprime[0]

    dh = delta_hue(*h, *Cprime)
    dH = 2 * np.sqrt(Cprime[0] * Cprime[1]) * np.sin(np.deg2rad(dh)/2)

    Lprimebar = L.mean()
    Cprimebar = Cprime.mean()

    hbar = hprimebar(*h, *Cprime)

    T = 1 - 0.17*np.cos(np.deg2rad(hbar-30)) + 0.24*np.cos(np.deg2rad(2*hbar)) + 0.32*np.cos(np.deg2rad(3*hbar+6)) - 0.20*np.cos(np.deg2rad(4*hbar-63))

    dtheta = 30 * np.exp(-((hbar-275)/25)**2)

    Rc = 2*np.sqrt(Cprimebar**7 / (Cprimebar**7 + 25**7))

    SL = 1 + (0.015 * (Lprimebar - 50)**2) / np.sqrt(20 + (Lprimebar - 50)**2)
    SC = 1 + 0.045 * Cprimebar
    SH = 1 + 0.015 * Cprimebar * T
    RT = -np.sin(np.deg2rad(2*dtheta)) * Rc

    dE00 = np.sqrt((dL/(kL * SL))**2 + (dCprime/(kC * SC))**2 + (dH/(kH * SH))**2 + RT * (dCprime/(kC*SC)) * (dH/(kH*SH)))
    return dE00

max_stars_repo_name: burggraaff/cbfu

max_issues_repo_name: burggraaff/cbfu

max_forks_repo_name: burggraaff/cbfu


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

plt.rcParams['figure.figsize'] = [12, 8]
c = [
    "#007fff",
    "#ff3616",
    "#138d75",
    "#7d3c98",
    "#fbea6a"
]  # Blue, Red, Green, Purple, Yellow

lines = sns.color_palette("Accent", 12)
lines = [lines[2], lines[0], lines[4]]
linestyles = ["solid"]*3

sns.set(
    style="ticks",
    rc={
        "font.family": "Arial",
        "font.size": 40,
        "axes.linewidth": 2,
        "lines.linewidth": 3,
    },
    font_scale=2.5,
    palette=sns.color_palette("Set2")
)
set2 = sns.color_palette("Set2", 5)

np.random.seed(2021)

samples = 15
I_test = 1000
noise = 3

total_errors = []
a, b = -10, 10
x = (b - a) * np.random.rand(samples) + a
y = x**3 + ((b - a) * np.random.randn(samples) + a)*noise
x_axis = np.linspace(np.min(x), np.max(x), 200)

a, b = -10, 10
x_test = (b - a) * np.random.rand(I_test) + a
y_test = x_test**3 + ((b - a) * np.random.randn(I_test) + a)*noise
x_test_axis = np.linspace(np.min(x_test), np.max(x_test), 200)

RMSEC, RMSEP = [], []
max_order = 8
px = 1/plt.rcParams['figure.dpi']
fig, axes = plt.subplots(1, 2, figsize=(1848*px, 965*px))
color_idx = 0
color_count = 0
for idx, order in enumerate(np.arange(1, max_order + 1)):
    p = np.polyfit(x, y, order)
    y_func = np.poly1d(p)

    RMSEP.append(np.sqrt(np.mean((y_test - y_func(x_test))**2)))
    RMSEC.append(np.sqrt(np.mean((y - y_func(x))**2)))

    if order in (1, 3, max_order):
        axes[0].plot(
            x_test_axis,
            y_func(x_test_axis),
            label=f"Polynomial order: {order}",
            c=set2[color_count],
            # c=lines[color_idx],
            linestyle=linestyles[color_idx],
            linewidth=30,
            alpha=0.5,
            zorder=max_order - order,
        )
        color_idx += 1
        color_count += 1

    if idx == 0:
        axes[0].scatter(
            x,
            y,
            label="Training set",
            marker="X",
            linewidth=.7,
            c=set2[color_count],
            edgecolor="w",
            s=240,
            zorder=1001,
        )
        # axes[0].scatter(
        #     x,
        #     y,
        #     label="Training set",
        #     marker="x",
        #     c="white",
        #     s=256,
        #     zorder=99,
        # )
        color_count += 1
        axes[0].scatter(
            x_test,
            y_test,
            label="Test set",
            # c=c[0],
            c=set2[color_count],
            s=20,
            # alpha=0.8,
            edgecolor="w",
            linewidth=.4,
            zorder=1000,
        )
        color_count += 1

axes[0].set_xlabel("X")
axes[0].set_ylabel("Y")

axbig = axes[1]

x = np.arange(1, len(RMSEC) + 1)
axbig.bar(
    x - .13, RMSEP, width=0.25, label="Test", facecolor=set2[2], edgecolor="none"
)
axbig.bar(
    x + .13,
    RMSEC,
    width=0.25,
    label="Training",
    facecolor=set2[1],
    edgecolor="none"
)
axbig.axhline(np.min(RMSEP), linestyle="--", c="grey", label="Lowest test")
axbig.set_xlabel("Polynomial order")
axbig.set_ylabel("Error")
axbig.legend(frameon=False, loc="best")
axbig.set_xlim(.7, max_order + .3)
axbig.set_ylim(-.05, np.max(RMSEP) + np.max(RMSEP)*0.05)
labels = np.arange(1, max_order + 1)
axbig.set_xticks(labels)
axbig.set_xticklabels(labels)

handles, labels = axes[0].get_legend_handles_labels()
labels = ["1st order", "3rd order", "8th order", "Training set", "Test set"]
axes[0].legend(handles, labels, loc="lower right")

plt.subplots_adjust(wspace=0.2, hspace=0.0, top=0.95, left=0.06, right=0.99)

fig.savefig("overfitting.pdf")
plt.show()


text:
# Copyright (c) Microsoft. All rights reserved.
# Licensed under the MIT license. See LICENSE file in the project root for full license information.

from PIL import Image, ImageDraw
import numpy as np

train_data_name = "../../data/ch18.train_shape.npz"
test_data_name = "../../data/ch18.test_shape.npz"

def circle(drawObj):
    x0 = np.random.randint(0,14)
    y0 = np.random.randint(0,14)
    x1 = np.random.randint(14,28)
    y1 = np.random.randint(14,28)
    drawObj.ellipse([x0,y0,x1,y1], fill='white', outline='white')

def rectangle(drawObj):
    x0 = np.random.randint(0,14)
    y0 = np.random.randint(0,14)
    x1 = np.random.randint(14,28)
    y1 = np.random.randint(14,28)
    drawObj.rectangle([x0,y0,x1,y1], fill='white', outline='white')

def triangle(drawObj):
    x0 = np.random.randint(0,14)
    y0 = np.random.randint(0,14)
    x1 = np.random.randint(14,28)
    y1 = np.random.randint(0,14)
    x2 = np.random.randint(0,14)
    y2 = np.random.randint(14,28)
    x3 = np.random.randint(14,28)
    y3 = np.random.randint(14,28)
    r = np.random.randint(0,3)
    if r == 0:
        drawObj.polygon([x0,y0,x1,y1,x2,y2], fill='white', outline='white')
    elif r == 1:
        drawObj.polygon([x0,y0,x1,y1,x3,y3], fill='white', outline='white')
    elif r == 2:
        drawObj.polygon([x0,y0,x2,y2,x3,y3], fill='white', outline='white')
    else:
        drawObj.polygon([x1,y1,x2,y2,x3,y3], fill='white', outline='white')

def line(drawObj):
    x0 = np.random.randint(0,28)
    y0 = np.random.randint(0,28)
    x1 = np.random.randint(0,28)
    y1 = np.random.randint(0,28)
    drawObj.line([x0,y0,x1,y1], fill='white', width=2)

def diamond(drawObj):
    x0 = np.random.randint(0,14)
    y0 = np.random.randint(0,14)
    x1 = np.random.randint(14,28)
    y1 = np.random.randint(14,28)
    x2 = (x0+x1)/2
    y2 = (y0+y1)/2
    drawObj.polygon([x2,y0,x1,y2,x2,y1,x0,y2], fill='white', outline='white')

if __name__ == '__main__':
    images_triangle = np.empty((1200,1,28,28))
    for i in range(1200):
        img = Image.new("L", [28,28], "black")
        drawObj = ImageDraw.Draw(img)
        triangle(drawObj)
        images_triangle[i,0] = np.array(img)

    images_circle = np.empty((1200,1,28,28))
    for i in range(1200):
        img = Image.new("L", [28,28], "black")
        drawObj = ImageDraw.Draw(img)
        circle(drawObj)
        images_circle[i,0] = np.array(img)

    images_rectangle = np.empty((1200,1,28,28))
    for i in range(1200):
        img = Image.new("L", [28,28], "black")
        drawObj = ImageDraw.Draw(img)
        rectangle(drawObj)
        images_rectangle[i,0] = np.array(img)

    images_line = np.empty((1200,1,28,28))
    for i in range(1200):
        img = Image.new("L", [28,28], "black")
        drawObj = ImageDraw.Draw(img)
        line(drawObj)
        images_line[i,0] = np.array(img)

    images_diamond = np.empty((1200,1,28,28))
    for i in range(1200):
        img = Image.new("L", [28,28], "black")
        drawObj = ImageDraw.Draw(img)
        diamond(drawObj)
        images_diamond[i,0] = np.array(img)

    images_train = np.empty((5000,1,28,28))
    label_train = np.empty((5000,1))
    images_train[0:1000] = images_circle[0:1000]
    label_train[0:1000] = 0
    images_train[1000:2000] = images_rectangle[0:1000]
    label_train[1000:2000] = 1
    images_train[2000:3000] = images_triangle[0:1000]
    label_train[2000:3000] = 2
    images_train[3000:4000] = images_diamond[0:1000]
    label_train[3000:4000] = 3
    images_train[4000:5000] = images_line[0:1000]
    label_train[4000:5000] = 4

    images_test = np.empty((1000,1,28,28))
    label_test = np.empty((1000,1))
    images_test[0:200] = images_circle[1000:1200]
    label_test[0:200] = 0
    images_test[200:400] = images_rectangle[1000:1200]
    label_test[200:400] = 1
    images_test[400:600] = images_triangle[1000:1200]
    label_test[400:600] = 2
    images_test[600:800] = images_diamond[1000:1200]
    label_test[600:800] = 3
    images_test[800:1000] = images_line[1000:1200]
    label_test[800:1000] = 4

    seed = np.random.randint(0,100)
    np.random.seed(seed)
    np.random.shuffle(images_train)
    np.random.seed(seed)
    np.random.shuffle(label_train)
    np.random.seed(seed)
    np.random.shuffle(images_test)
    np.random.seed(seed)
    np.random.shuffle(label_test)

    np.savez(train_data_name, data=images_train, label=label_train)
    np.savez(test_data_name, data=images_test, label=label_test)

max_stars_repo_name: microsoft/ai-edu

max_issues_repo_name: microsoft/ai-edu

max_forks_repo_name: microsoft/ai-edu


text:
"""
Created on Fri Apr 23 11:34:02 2021
@author: mchini
"""

#%% load basic packages

import numpy as np
from numba import jit

#%% define elephant STTC

@jit(nopython=True)
def run_P(spiketrain_i, spiketrain_j):
    """
    Check every spike in train 1 to see if there's a spike in train 2
    within 0.005
    """
    N2 = len(spiketrain_j)

    # Search spikes of spiketrain_i in spiketrain_j
    # ind will contain index of
    ind = np.searchsorted(spiketrain_j, spiketrain_i)

    # To prevent IndexErrors
    # If a spike of spiketrain_i is after the last spike of spiketrain_j,
    # the index is N2, however spiketrain_j[N2] raises an IndexError.
    # By shifting this index, the spike of spiketrain_i will be compared
    # to the last 2 spikes of spiketrain_j (negligible overhead).
    # Note: Not necessary for index 0 that will be shifted to -1,
    # because spiketrain_j[-1] is valid (additional negligible comparison)
    ind[ind == N2] = N2 - 1

    # Compare to nearest spike in spiketrain_j BEFORE spike in spiketrain_i
    close_left = np.abs(
        spiketrain_j[ind - 1] - spiketrain_i) <= 0.005
    # Compare to nearest spike in spiketrain_j AFTER (or simultaneous)
    # spike in spiketrain_j
    close_right = np.abs(
        spiketrain_j[ind] - spiketrain_i) <= 0.005

    # spiketrain_j spikes that are in [-0.005, 0.005] range of spiketrain_i
    # spikes are counted only ONCE (as per original implementation)
    close = close_left + close_right

    # Count how many spikes in spiketrain_i have a "partner" in
    # spiketrain_j
    return np.count_nonzero(close)

@jit(nopython=True)
def run_T(spiketrain, len_rec):
    """
    Calculate the proportion of the total recording time 'tiled' by spikes.
    """
    N = len(spiketrain)
    time_A = 2 * N * 0.005  # maximum possible time

    if N == 1:  # for just one spike in train
        if spiketrain[0] < 0.005:
            time_A += -0.005 + spiketrain[0]
        if spiketrain[0] + 0.005 > len_rec:
            time_A += -0.005 - spiketrain[0] + len_rec
    else:  # if more than one spike in train
        # Vectorized loop of spike time differences
        diff = np.diff(spiketrain)
        diff_overlap = diff[diff < 2 * 0.005]
        # Subtract overlap
        time_A += -2 * 0.005 * len(diff_overlap) + np.sum(diff_overlap)

        # check if spikes are within 0.005 of the start and/or end
        # if so subtract overlap of first and/or last spike
        if spiketrain[0] < 0.005:
            time_A += spiketrain[0] - 0.005

        if (len_rec - spiketrain[N - 1]) < 0.005:
            time_A += -spiketrain[-1] - 0.005 + len_rec

    T = time_A / len_rec
    return T

@jit(forceobj=True, parallel=True)
def matrix_spike_time_tiling_coefficient(spike_matrix):
    """
    Calculates the Spike Time Tiling Coefficient (STTC) for an entire spike
    matrix. Uses numba to speed up computations.
    It assumes the spike matrix to have 1ms resolution
    """

    # length of recording (assumes spike matrix has 1ms resolution)
    len_rec = np.shape(spike_matrix)[1] / 1000
    num_spike_trains = np.shape(spike_matrix)[0]
    # initialize STTC
    STTC = np.zeros(int(num_spike_trains * (num_spike_trains - 1) / 2))
    idx = 0
    # loop over spike pairs
    for i in range(num_spike_trains):
        # extract spike times
        spiketrain_i = np.array(np.nonzero(spike_matrix[i, :])[0]) / 1000
        for j in range(i + 1, num_spike_trains):
            # extract spike times
            spiketrain_j = np.array(np.nonzero(spike_matrix[i, :])[0]) / 1000
            N1 = len(spiketrain_i)
            N2 = len(spiketrain_j)

            if N1 == 0 or N2 == 0:
                index = np.nan
            else:
                TA = run_T(spiketrain_i, len_rec)
                TB = run_T(spiketrain_j, len_rec)
                PA = run_P(spiketrain_i, spiketrain_j)
                PA = PA / N1
                PB = run_P(spiketrain_j, spiketrain_i)
                PB = PB / N2
                # check if the P and T values are 1 to avoid division by zero
                # This only happens for TA = PB = 1 and/or TB = PA = 1,
                # which leads to 0/0 in the calculation of the index.
                # In those cases, every spike in the train with P = 1
                # is within 0.005 of a spike in the other train,
                # so we set the respective (partial) index to 1.
                if PA * TB == 1:
                    if PB * TA == 1:
                        index = 1.
                    else:
                        index = 0.5 + 0.5 * (PB - TA) / (1 - PB * TA)
                elif PB * TA == 1:
                    index = 0.5 + 0.5 * (PA - TB) / (1 - PA * TB)
                else:
                    index = 0.5 * (PA - TB) / (1 - PA * TB) + \
                    0.5 * (PB - TA) / (1 - PB * TA)
            STTC[idx] = (index)
            idx += 1
    return STTC

max_stars_repo_name: thmspfffr/neonates

max_issues_repo_name: thmspfffr/neonates

max_forks_repo_name: thmspfffr/neonates

reason: import numpy,from numba


text:
"""
Computes tree figure, with size of node proportional to value added

John Stachurski, August 2012.

NOTES: Try delta=1.01 and c(x)=np.exp(x) - 1 vs delta=1.1 and c(x)=np.exp(2 * x) -
1.  Same number of firms but different pattern for value added.  Because of
more convexity in the price function?

"""

from __init__ import *
from rp import RPtree
import networkx as nx
import numpy as np

model1 = { 'delta' : 1.02, 'c' : lambda x: np.exp(4 * x) - 1, 'k' : 2 }
model2 = { 'delta' : 1.20, 'c' : lambda x: np.exp(4 * x) - 1, 'k' : 2 }
model3 = { 'delta' : 1.01, 'c' : lambda x: np.exp(2 * x) - 1, 'k' : 2 }
model4 = { 'delta' : 1.05, 'c' : lambda x: np.exp(2 * x) - 1, 'k' : 2 }
model5 = { 'delta' : 1.01, 'c' : lambda x: np.exp(x) - 1, 'k' : 3 }
model6 = { 'delta' : 1.05, 'c' : lambda x: np.exp(x) - 1, 'k' : 3 }
model7 = { 'delta' : 1.01, 'c' : lambda x: np.exp(2 * x) - 1, 'k' : 3 }
model8 = { 'delta' : 1.15, 'c' : lambda x: np.exp(x**1.01) - 1, 'k' : 6 }
model9 = { 'delta' : 1.15, 'c' : lambda x: np.exp(x**1.00) - 1, 'k' : 6 }

# Set the model
model = model1
fig, ax = plt.subplots()

ps = RPtree(delta=model['delta'], c=model['c'], k=model['k'])

levels = ps.compute_stages()
ts, ls, vs = [], [], []
i = 0
for level in levels:
    t, l, v = level
    ts.append(t)
    ls.append(l)
    vs.append(v)
    print "level:", i, "; sum va:", sum(vs)
    i += 1
# Count the number of levels and firms
num_levels = len(levels)
num_firms = 0
for n in range(num_levels):
    num_firms += ps.k**n
print "Number of firms:", num_firms

## Work out value added for each firm, where firms are
## enumerated top to bottom and left to right in the vertical tree
value_added = []
for i in range(num_levels):
    for j in range(ps.k**i):
        value_added.append(vs[i])
# Now rescale value added relative to first one
value_added = [250 * (v / value_added[0]) for v in value_added]
value_added = np.array(value_added)
# Make the smallest value at least 2% of the largest
temp = value_added.max() * 0.02
value_added = np.where(value_added < temp, temp, value_added)

G = nx.balanced_tree(ps.k, num_levels - 1)
pos = nx.graphviz_layout(G, prog='dot') # twopi or dot
ax.axis('off', frame_on=False)
nx.draw(G,
        pos,
        ax=ax,
        node_size=value_added,
        node_shape='o',
        alpha=0.8,
        node_color="blue",
        with_labels=False)
ax.axes.get_xaxis().set_visible(False)
ax.axes.get_yaxis().set_visible(False)

fig.show()

#plt.savefig('circular_tree.png')
#for i in range(num_firms):
#    x, y = pos[i]
#    x = x * np.random.uniform(low=0.98, high=1.02)
#    y = y * np.random.uniform(low=0.98, high=1.02)
#    pos[i] = (x, y)

reason: import numpy,import networkx


text:
"""
Author: Konstantinos Angelopoulos
Date: 04/02/2020
All rights reserved.
Feel free to use and modify and if you like it give it a star.

Capture background and store it for fabric localization
"""
from pykinect2.PyKinectV2 import *
from pykinect2 import PyKinectV2
from pykinect2 import PyKinectRuntime
import cv2
import numpy as np
import os

# Create Kinect
kinect = PyKinectRuntime.PyKinectRuntime(PyKinectV2.FrameSourceTypes_Depth | PyKinectV2.FrameSourceTypes_Color)
shot = True  # Flag for shot

# Loop until photo taken
while shot is True:
    # wait until frame is returned
    if kinect.has_new_color_frame():
        color_frame = kinect.get_last_color_frame()
        colorImage = color_frame.reshape((kinect.color_frame_desc.Height, kinect.color_frame_desc.Width, 4)).astype(np.uint8)
        cv2.imwrite(os.path.join(os.path.dirname(os.path.realpath(__file__)), 'images/background.png'), colorImage)
        shot = False

max_stars_repo_name: KonstantinosAng/KinectPython

max_issues_repo_name: KonstantinosAng/KinectPython

max_forks_repo_name: KonstantinosAng/KinectPython


text:
import argparse
import numpy as np
import os
import meld_classifier.experiment as exp
import csv

if __name__ == "__main__":

    parser = argparse.ArgumentParser(
        description="Script to rerun failed experiments (only for slurm jobs). Checks for results folders and reruns in their absence"
    )
    parser.add_argument("--params-csv", help="path to params csv file where things failed", default="params.csv")
    parser.add_argument("--run-on-cpu", dest="run_on_cpu", action="store_true")
    args = parser.parse_args()

    # check which lines have no results already
    lines_to_rerun = []
    with open(args.params_csv, "r") as f:
        reader = csv.DictReader(f, delimiter=",")
        for i, row in enumerate(reader):

            experiment_path = row["experiment_path"]
            experiment_name = row["experiment_name"]
            if not os.path.isfile(
                os.path.join(experiment_path, "results", "per_subject_{}_optimal.json".format(experiment_name))
            ):
                lines_to_rerun.append(i)
    # reruns only these lines passing a comma separated list of them
    exp.resubmit_experiments_array(lines_to_rerun, args.params_csv, run_on_cpu=args.run_on_cpu)


text:
# -*- coding: utf-8 -*-
"""Functions for generating test data.

This program code is part of the MultiPy (Multiple Hypothesis Testing in
Python) package.

Author: Tuomas Puoliväli
Email: tuomas.puolivali@helsinki.fi
License: Revised 3-clause BSD
Source: https://github.com/puolival/multipy/blob/master/data.py
Last modified: 25th July 2018.

References:

[1] Neuhaus KL, von Essen R, Tebbe U, Vogt A, Roth M, Riess M, Niederer W,
    Forycki F, Wirtzfeld A, Maeurer W, Limbourg P, Merx W, Haerten K (1992):
    Improved thrombolysis in acute myocardial infarction with front-loaded
    administration of alteplase: Results of the rt-PA-APSAC patency study
    (TAPS). Journal of the American College of Cardiology 19(5):385-391.

[2] Benjamini Y, Hochberg Y (1995): Controlling the false discovery rate:
    A practical and powerful approach to multiple testing. Journal of Royal
    Statistical Society. Series B (Methodological): 57(1):289-300.

[3] Reiss PT, Schwartzman A, Lu F, Huang L, Proal E (2012): Paradoxical
    results of adaptive false discovery rate procedures in neuroimaging
    studies. NeuroImage 63(4):1833-1840.

[4] Efron B (2008): Simultaneous inference: When should hypothesis testing
    problems be combined? The Annals of Applied Statistics 2(1):197-223.

[5] Bennett CM, Wolford GL, Miller MB (2009): The principled control of
    false positives in neuroimaging. Social Cognitive and Affective
    Neuroscience 4(4):417-422.
"""

import numpy as np
from numpy.random import normal, permutation

from scipy.stats import ttest_ind

def neuhaus(permute=False):
    """Function that returns the Neuhaus et al. data that was re-analyzed in
    the classic Benjamini & Hochberg (1995) FDR paper.

    Input arguments:
    permute - If true, the p-values are returned in random order. If false,
              the p-values are returned in ascending order.
    """
    pvals = np.array([0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298,
                      0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590,
                      1.000], dtype='float')
    if (permute):
        m = len(pvals)
        pvals = pvals[permutation(m)]
    return pvals

def two_group_model(N=25, m=1000, pi0=0.1, delta=0.7):
    """A two-group model for generating test data (described in [3] and
    elsewhere). The default input arguments can be used to reproduce the
    result reported by Reiss and colleagues in Figure 2.A.

    Input arguments:
    N     - Number of samples in each group or condition.
    m     - Number of variables
    pi0   - Proportion of null effects among the m variables.
    delta - Location parameter of the non-null part of the distribution of Y,
            which controls the effect size.

    Output arguments:
    tstat - Test statistics (Student's t's)
    pvals - P-values corresponding to the test statistics
    """
    X = normal(loc=0, scale=1, size=(N, m))
    Y = np.hstack([normal(loc=0, scale=1, size=(N, int(pi0*m))),
                   normal(loc=delta, scale=1, size=(N, int(round((1-pi0)*m, 1))))])
    # Two-sample t-test
    tstat, pvals = ttest_ind(X, Y, axis=0)
    return tstat, pvals

def separate_class_model(a_N=25, b_N=25, a_m=500, b_m=500, a_pi0=0.25,
                         b_pi0=0.75, a_delta=0.5, b_delta=0.6):
    """The separate classes model described by Efron [4]. The prefixes
    a_ and b_ refer to the two separate classes. See two_group_model() for
    the definition of the variables."""
    A_tstat, A_pvals = two_group_model(a_N, a_m, a_pi0, a_delta)
    B_tstat, B_pvals = two_group_model(b_N, b_m, b_pi0, b_delta)
    # Combine
    tstats, pvals = (np.hstack([A_tstat, B_tstat]),
                     np.hstack([A_pvals, B_pvals]))
    return tstats, pvals

def square_grid_model(nl=100, sl=60, N=25, delta=0.7, equal_var=True):
    """Bennett et al [5] model.

    Input arguments:
    nl : int
        The side length of the noise region.
    sl : int
        The side length of the signal region.
    N : int
        Sample size in each of the two groups.
    delta : float
        Effect size. Difference in means between the two normal
        distributions.
    equal_var : bool
        Whether to assume equal population variances while performing
        the Student's t-tests. See also scipy.stats.ttest_ind.

    Output arguments:
    P : ndarray
        T-test p-values.
    tstats : ndarray
        The corresponding test statistics.
    X, Y : ndarray
        Samples in the two groups.
    """

    """Generate the noise statistics."""
    X_noise = normal(loc=0, scale=1, size=(nl, nl, N))
    Y_noise = normal(loc=0, scale=1, size=(nl, nl, N))
    noise_tstats, noise_pvals = ttest_ind(
        X_noise, Y_noise, axis=2, equal_var=equal_var)

    """Generate the signal statistics."""
    X_signal = normal(loc=delta, scale=1, size=(sl, sl, N))
    Y_signal = normal(loc=0, scale=1, size=(sl, sl, N))
    signal_tstats, signal_pvals = ttest_ind(
        X_signal, Y_signal, axis=2, equal_var=equal_var)

    """Combine the data so that the desired spatial structure is
    obtained."""
    P = noise_pvals
    d = (nl-sl) // 2
    P[d:(nl-d), d:(nl-d)] = signal_pvals
    tstats = noise_tstats
    tstats[d:(nl-d), d:(nl-d)] = signal_tstats
    X, Y = X_noise, Y_noise
    X[d:(nl-d), d:(nl-d)] = X_signal
    Y[d:(nl-d), d:(nl-d)] = Y_signal
    return P, tstats, X, Y

def spatial_separate_classes_model(delta1, delta2, N=25):
    """Function for simulating data under the spatial separate-classes
    model.

    Input arguments:
    ================
    delta1, delta : float
        The effect sizes at the first and second signal regions
        respectively.
    N : int
        The sample size in each group.

    Output arguments:
    =================
    pvals, tstats : ndarray
        P-values and t-statistics of each test.
    X, Y : ndarray
        The simulated data.
    """
    # TODO: make nl, sl, & N parameters of the function.
    P_a, tstats_a, X_a, Y_a = square_grid_model(nl=45, sl=15, N=N,
                                                delta=delta1, equal_var=True)
    P_b, tstats_b, X_b, Y_b = square_grid_model(nl=45, sl=15, N=N,
                                                delta=delta2, equal_var=True)

    return (np.hstack([P_a, P_b]), np.hstack([tstats_a, tstats_b]),
            np.hstack([X_a, X_b]), np.hstack([Y_a, Y_b]))

def two_class_grid_model():
    """The Bennett et al. [5] model extended for the two separate classes
    case.

    Input arguments:
    nl : int
        The side length of the noise region.
    sl<i> : int
        The side length of the ith signal region.
    N : int
        Sample size in each group.
    delta<i> : float
        Effect size in ith class. Difference in means between the two
        normal distributions.
    """

    P_a, T_a, X_a, Y_a = square_grid_model(nl=45, sl=25, N=25, delta=0.25)
    P_b, T_b, X_b, Y_b = square_grid_model(nl=45, sl=25, N=25, delta=0.50)
    P_c, T_c, X_c, Y_c = square_grid_model(nl=45, sl=25, N=25, delta=0.75)
    P_d, T_d, X_d, Y_d = square_grid_model(nl=45, sl=25, N=25, delta=1.00)

    P = np.vstack([np.hstack([P_a, P_b]), np.hstack([P_c, P_d])])
    T = np.vstack([np.hstack([T_a, T_b]), np.hstack([T_c, T_d])])
    X = np.vstack([np.hstack([X_a, X_b]), np.hstack([X_c, X_d])])
    Y = np.vstack([np.hstack([Y_a, Y_b]), np.hstack([Y_c, Y_d])])

    return P, T, X, Y

max_stars_repo_name: OrionisBio/multipy

max_issues_repo_name: OrionisBio/multipy

max_forks_repo_name: OrionisBio/multipy

reason: import numpy,from numpy,from scipy


text:
import time
import os, sys
import argparse
from argparse import Namespace
from multiprocessing import Process
from typing import Dict, Any, Optional

import torch
from torch import nn
import numpy as np

import dopt
from dopt import NEIOptimizer, Trainer, Server
from dopt.utils import get_output_shape
# from camus_objective import run_kfolds # The objective function

import warnings
warnings.filterwarnings("ignore")

# The configurations
CONFIG = {}
CONFIG["computer_list"] = {
    "acet_update": ['tst008@acet116-lnx-10.bucknell.edu'], # To git pull
    "acet": [
        'tst008@acet116-lnx-11.bucknell.edu',
        'tst008@acet116-lnx-12.bucknell.edu',
        'tst008@acet116-lnx-13.bucknell.edu',
        'tst008@acet116-lnx-14.bucknell.edu',
        'tst008@acet116-lnx-15.bucknell.edu',
#         'tst008@acet116-lnx-16.bucknell.edu',
#         'tst008@acet116-lnx-17.bucknell.edu',
#         'tst008@acet116-lnx-18.bucknell.edu',
#         'tst008@acet116-lnx-1.bucknell.edu',
#         'tst008@acet116-lnx-20.bucknell.edu',
#         'tst008@acet116-lnx-21.bucknell.edu',
    ],
#     "tung-torch": ['tung@jvs008-r1.bucknell.edu']
}
# Commands to run on target machines here
# CONFIG["commands"] = {
#     "acet_update": "cd PycharmProjects/distributed-optimizer/ && git pull" + \
#                    " && module switch python/3.7-2020-05-28" + \
#                    " && export LD_LIBRARY_PATH=/usr/remote/lib:/usr/remote/anaconda-3.7-2020-05-28/lib" + \
#                    " && python3 ~/PycharmProjects/distributed-optimizer/distributed_optimizer.py --run_as trainer",
#     "acet": "sleep 10 && module switch python/3.7-2020-05-28" + \
#                    " && export LD_LIBRARY_PATH=/usr/remote/lib:/usr/remote/anaconda-3.7-2020-05-28/lib" + \
#                    " && python3 ~/PycharmProjects/distributed-optimizer/distributed_optimizer.py --run_as trainer",
#     "tung-torch": "/opt/anaconda/envs/jupyter37/bin/python ~/pj/dopt_v2/distributed_optimizer.py --run_as trainer --data_folder ~/pj/dopt_v2/data/CroppedYale/"
# }

CONFIG["commands"] = {
    "acet_update": "sleep 20 && pkill -9 python && echo 'Hey'",
    "acet": "sleep 20 && pkill -9 python && echo 'Hey'",
    "tung-torch": "sleep 20 && pkill -9 python && echo 'Hey'"
}
CONFIG["server"] = {
    "host": "jvs008-r1.bucknell.edu",
    "port": 15555
}
CONFIG["trainer"] = {
    "username": "tst008",
    "num_constraints": 1
}
CONFIG["optimizer"] = {
    "bounds": {
        "conv1": [3, 16],
        "conv1_kernel": [2, 10],
        "conv2": [20, 32],
        "conv2_kernel": [2, 10],
        "dropout1": [0.1, 1],
        "maxpool1": [2, 10],
        "maxpool2": [2, 10],
        'batch_size': [2, 7],
        'lr': [0.001, 10.0]
    },
    "initial_candidates": None,
    "device": "cpu",
    "seed": 0,
    "filename": "yaleface.json"
}

def get_feasibility(candidate) -> float:
    expected_input_shape = (1, 1, 192, 168)
    conv1_in = round(candidate["conv1"])
    conv1_kernel = round(candidate["conv1_kernel"])
    conv2_in = round(candidate["conv2"])
    conv2_kernel = round(candidate["conv2_kernel"])
    maxpool1_in = round(candidate["maxpool1"])
    maxpool2_in = round(candidate["maxpool2"])
    # Reconstruct the layers for calculation
    conv1 = nn.Conv2d(1, conv1_in, conv1_kernel, 1)
    conv2 = nn.Conv2d(conv1_in, conv2_in, conv2_kernel, 1)
    maxpool1 = nn.MaxPool2d(maxpool1_in)
    maxpool2 = nn.MaxPool2d(maxpool2_in)

    try:
        conv1_out = get_output_shape(conv1, expected_input_shape)
        maxpool1_out = get_output_shape(maxpool1, conv1_out)
        conv2_out = get_output_shape(conv2, maxpool1_out)
        maxpool2_out = get_output_shape(maxpool2, conv2_out)
    except RuntimeError as e:
        if "Output size is too small" in str(e):
            return 2
        else:
            raise e
    fc1_in = np.prod(list(maxpool2_out)) # Flatten

    feasibility = -0.1 # Default is 0, meaning feasible, >0 means not feasible
    if conv1_in > conv2_in:
        feasibility += 1
    if fc1_in > 10**5:
        feasibility += 1
    print(f"Conv1: {get_output_shape(conv1, expected_input_shape)} % {maxpool1_in}, Conv2: {get_output_shape(conv2, maxpool1_out)} % {maxpool2_in}, Linear: {fc1_in}, Feasibility: {feasibility}")
    return feasibility

# Plug in the objective function here
def objective_function(data_folder, candidate):
    feasibility = get_feasibility(candidate)
    if feasibility > 0:
        print("Infeasible!")
        observation = {
            "objective": [0.001, 0.001],
            "constraints": [feasibility]
        }
        return observation
    print("Returning the observation")

    # Simulate input
    input_args = Namespace(
        data_folder=data_folder,
        conv1=round(candidate["conv1"]),
        conv1_kernel=round(candidate["conv1_kernel"]),
        conv2=round(candidate["conv2"]),
        conv2_kernel=round(candidate["conv2_kernel"]),
        dropout1=candidate["dropout1"],
        maxpool1=round(candidate["maxpool1"]),
        maxpool2=round(candidate["maxpool2"]),
        no_cuda=False,
        seed=2,
        batch_size=round(candidate["batch_size"]),
        test_batch_size=1000,
        epochs=1,
        lr=candidate["lr"],
        gamma=0.7,
        log_interval=250, # was 250
        save_model=False,
        num_folds = 5
    )
    time.sleep(60)
    mean = 0
    variance = 0.1
    feasibility = 0
    observation = {
        "objective": [mean, variance],
        "constraints": [feasibility]
    }
    return observation

# Calls start_optimizer and start_trainers simultaneously
def start_server():
    optimizer = NEIOptimizer(
        CONFIG["optimizer"]["filename"],
        CONFIG["optimizer"]["bounds"],
        initial_candidates=CONFIG["optimizer"]["initial_candidates"],
        device=CONFIG["optimizer"]["device"],
        seed=CONFIG["optimizer"]["seed"]
    )
    server = Server(optimizer, CONFIG)
    server.run()

def start_trainers():
    from functools import partial
    of = partial(objective_function, args.data_folder)
    trainer = Trainer(
        of,
        CONFIG["trainer"]["username"],
        CONFIG["server"]["host"],
        CONFIG["server"]["port"],
        num_constraints=CONFIG["trainer"]["num_constraints"]
    )
    trainer.run()

if __name__ == "__main__":
    parser = argparse.ArgumentParser(prog='distributed_optimizer.py',
                                     description='''Optimize objective function of specified by a `Trainer`''')
    parser.add_argument('--run_as', action='store', dest='run_as',
                           help='Specify the role of the machine (server or trainer). Defaults to server',
                           type=str, required=False,
                           default="server")
    parser.add_argument('--data_folder', action='store', dest='data_folder',
                           help='Specify the directory to the data folder (for clients only)',
                           type=str, required=False,
                           default="~/PycharmProjects/summer/data/CroppedYale/")
    args = parser.parse_args()

    # Can modify these code to accomodate more options
    # E.g. Run different Trainers on same task
    if args.run_as == "server":
        start_server()
    elif args.run_as == "trainer":
        start_trainers()

max_stars_repo_name: sontungtran/dopt

max_issues_repo_name: sontungtran/dopt

max_forks_repo_name: sontungtran/dopt


text:
import sqlite3
import pickle
import numpy as np

"""
class for dynamically loading a reaction network
"""

sql_get_reaction = """
    SELECT * FROM reactions WHERE reaction_id = ?;
"""

sql_get_reaction_range = """
    SELECT * FROM reactions WHERE ? <= reaction_id AND reaction_id < ?;
"""

sql_get_redox = """
    SELECT * FROM reactions WHERE is_redox = 1;
"""

sql_get_trajectory = """
    SELECT * FROM trajectories;
"""

sql_get_initial_state = """
    SELECT * FROM initial_state;
"""

class NetworkLoader:

    def __init__(
            self,
            network_database,
            mol_entries_pickle,
            initial_state_database=None
    ):

        self.rn_con = sqlite3.connect(network_database)

        with open(mol_entries_pickle, 'rb') as f:
            self.mol_entries = pickle.load(f)

        cur = self.rn_con.cursor()
        metadata = list(cur.execute("SELECT * FROM metadata"))[0]
        self.number_of_species = metadata[0]
        self.number_of_reactions = metadata[1]

        if initial_state_database:
            self.initial_state_con = sqlite3.connect(initial_state_database)

        self.reactions = {}

    def get_all_redox_reactions(self):
        redox_reactions = []
        cur = self.rn_con.cursor()
        for res in cur.execute(sql_get_redox):
            reaction = {}
            reaction['number_of_reactants'] = res[1]
            reaction['number_of_products'] = res[2]
            reaction['reactants'] = res[3:5]
            reaction['products'] = res[5:7]
            reaction['rate'] = res[7]
            reaction['dG'] = res[8]
            reaction['dG_barrier'] = res[9]
            redox_reactions.append(reaction)

        return redox_reactions

    def get_reactions_in_range(self, lower_bound, upper_bound):
        """
        get range of reactions from database but don't cache them
        """
        cur = self.rn_con.cursor()
        for res in cur.execute(sql_get_reaction_range,
                               (lower_bound, upper_bound)):
            reaction = {}
            reaction['reaction_id'] = res[0]
            reaction['number_of_reactants'] = res[1]
            reaction['number_of_products'] = res[2]
            reaction['reactants'] = res[3:5]
            reaction['products'] = res[5:7]
            reaction['rate'] = res[7]
            reaction['dG'] = res[8]
            reaction['dG_barrier'] = res[9]
            yield reaction

    def index_to_reaction(self, reaction_index):

        """
        this method gets called a lot, so we cache the reactions to
        minimize database interaction
        """

        if reaction_index in self.reactions:
            return self.reactions[reaction_index]

        else:
            print("fetching data for reaction", reaction_index)
            cur = self.rn_con.cursor()
            res = list(
                cur.execute(sql_get_reaction, (reaction_index,))
            )[0]
            reaction = {}
            reaction['number_of_reactants'] = res[1]
            reaction['number_of_products'] = res[2]
            reaction['reactants'] = res[3:5]
            reaction['products'] = res[5:7]
            reaction['rate'] = res[7]
            reaction['dG'] = res[8]
            reaction['dG_barrier'] = res[9]
            self.reactions[reaction_index] = reaction
            return reaction

    def load_trajectories(self):

        cur = self.initial_state_con.cursor()

        # trajectories[seed][step] = (reaction_id, time)
        trajectories = {}
        for row in cur.execute(sql_get_trajectory):
            seed = row[0]
            step = row[1]
            reaction_id = row[2]
            time = row[3]

            if seed not in trajectories:
                trajectories[seed] = {}

            trajectories[seed][step] = (reaction_id, time)

        self.trajectories = trajectories

    def load_initial_state(self):

        cur = self.initial_state_con.cursor()
        initial_state_dict = {}

        for row in cur.execute(sql_get_initial_state):
            initial_state_dict[row[0]] = row[1]

        initial_state_array = np.zeros(
            self.number_of_species,
            dtype=int
        )

        for i in range(self.number_of_species):
            initial_state_array[i] = initial_state_dict[i]

        self.initial_state_dict = initial_state_dict
        self.initial_state_array = initial_state_array

max_stars_repo_name: espottesmith/HiPRGen

max_issues_repo_name: espottesmith/HiPRGen

max_forks_repo_name: espottesmith/HiPRGen