text:
import sortedcontainers
from collections.abc import Iterable

import networkx as nx
from copy import deepcopy
import bidict
import numpy as np

import tnetwork as tn
import numbers
from tnetwork.dyn_graph.dyn_graph import DynGraph

from datetime import datetime, timezone
from tnetwork.dyn_graph.encodings import code_length_SN_M,code_length_SN_E

class DynGraphSN(DynGraph):
    """
    A class to represent dynamic graphs as snapshot sequence.

    Each snapshot is represented as a networkx graph, and is associated to a time step identifier. The time step can be
    an position in the sequence (1,2,3,...) or an arbitrary temporal indicator (year, timestamp...).

    Snpashots are ordered according to their time step identifier using a sorted dictionary (SortedDict).

    """

    def __init__(self, data=None,frequency=1):
        """
        Instanciate a new graph, with or without initial data

        :param data: can be a dictionary {time step:graph} or a list of graph, in which sase time steps are integers starting at 0
        :param frequency: minimal time difference between two observations. Default: 1
        """

        self._snapshots = sortedcontainers.SortedDict()
        self.frequency(frequency)
        if data!=None:
            if isinstance(data,dict):
                self._snapshots = sortedcontainers.SortedDict(data)
            elif isinstance(data,list):
                self._snapshots = sortedcontainers.SortedDict({i:g for i,g in enumerate(data)})
            else:
                raise Exception("data should be a list or a dictionary")

    def start(self):
        """
        Time of the first snapshot

        :return:
        """
        return self.snapshots_timesteps()[0]

    def end(self):
        """
        Time of the last snapshot

        :return:
        """
        return self.snapshots_timesteps()[-1]

    def add_node_presence(self, n, time):
        """
        Add presence for a node at a time

        :param n: node
        :param time: a snapshot time
        """

        if not time in self.snapshots_timesteps():
            self.add_snapshot(time)
        self._snapshots[time].add_node(n)

    def add_nodes_presence_from(self, nodes,times):
        """
        Add nodes for times

        For each node in nodes, add it for each time in times.

        :param nodes: list of nodes, or a single node
        :param times: list of times of same length as node, or a single time
        """
        if not isinstance(nodes,Iterable):
            nodes = [nodes]

        if not isinstance(times, Iterable): #means it is a single time, not list of times
            times = [times]

        for i,node in enumerate(nodes):
            for t in times:
                self.add_node_presence(node, t)

    def remove_node_presence(self, n, time):
        """
        Remove presence for a node at a time

        :param n: node
        :param time: a snapshot time
        """

        self._snapshots[time].remove_node(n)

    def add_interaction(self,u,v,time):
        """
        Add a single interaction at a single time step.

        :param u: first node
        :param v: second node
        :param time: time step identifier
        """

        self.add_interactions_from([(u,v)],[time])

    def add_interactions_from(self, nodePairs, times):
        """
        Add interactions between the provided node pairs for the provided times.

        Add each provided nodePair at each provided time

        :param nodePairs: list of pairs of nodes, or a single pair of nodes as a tuple or set
        :param times: list of times as integer or a single integer
        """
        #note: could be optimized

        #if len(nodePairs)==2 and (isinstance(nodePairs[0],str) or not isinstance(nodePairs[0],Iterable)):
        list_element_example = list(nodePairs)[0]
        if not isinstance(list_element_example, Iterable) or isinstance(list_element_example,str):
            nodePairs = [nodePairs]

        if not isinstance(times,Iterable):
            times = [times]

        for t in times:
            if not t in self._snapshots:
                self.add_snapshot(t)
            self.apply_nx_function(nx.Graph.add_edges_from,start=t,stop=t,ebunch_to_add=nodePairs)

    def remove_interaction(self,u,v,time):
        """
        Remove a single interaction at a single time step.

        Note: it does not remove the node

        :param u: first node
        :param v: second node
        :param time: time step identifier
        """

        self.remove_interactions_from([(u,v)],[time])

    def remove_interactions_from(self, nodePairs, times):
        """
        Remove interactions between the provided node pairs for the provided times.

        If one of the two parameters is a single element,
        will remove the node pair at all provided time steps, or all the node pairs at the provided time step.

        :param nodePairs: list of pairs of nodes, or a single pair of nodes
        :param times: list of times for this node, or a single time
        :return:
        """
        #note: could be optimized

        if len(nodePairs)==2 and not isinstance(nodePairs[0],Iterable):
            nodePairs = [nodePairs]*len(times)

        if not isinstance(times,Iterable):
            times = [times]*len(nodePairs)

        for t in times:
            print("-------",t)
            self.apply_nx_function(nx.Graph.remove_edges_from,start=t,stop=t,ebunch=nodePairs)

    def graph_at_time(self,t):
        """
        Return the graph as it is at time t

        :param t: a time step identifier
        :return: the graph as a networkx graph
        """
        return self.snapshots(t)

    def add_snapshot(self, t=None, graphSN=None):
        """
        Add a snapshot for a time step t

        :param t: the time step identifier. If none, the last one + 1
        :param graphSN: the graph to add (networkx object), if None, add an empty snapshot
        """
        if t==None:
            if len(self.snapshots_timesteps())==0:
                t=0
            else:
                t=self.snapshots_timesteps()[-1]+1
        if graphSN==None:
            graphSN=nx.Graph()
        self._snapshots[t]=graphSN

    def remove_snapshot(self, t):
        """
        Remove a snapshot

        :param t: the time at which to remove a snapshot
        :return:
        """
        del self._snapshots[t]

    def snapshots_timesteps(self):
        """
        Return the list of time steps

        :return: list of time steps
        """
        return list(self._snapshots.keys())

    def apply_nx_function(self,function,start=None,stop=None,**kwargs):
        """
        Apply a networkx function to each snapshot and return the list of result. Parameters of the function to apply can be passed as parameter to this function.
        example
        >>> dg = DynGraphSN.graph_socioPatterns2012()
        >>> dg.apply_nx_function(nx.nodes)
        >>> dg.apply_nx_function(nx.Graph.add_node,node_for_adding="nodeTest")

        :param function: the networkx function
        :return: the list of results for each snapshot
        """

        if start==None:
            start = self.snapshots_timesteps()[0]
        if stop==None:
            stop = self.snapshots_timesteps()[-1]
        to_return = []
        keys = self.snapshots().irange(start, stop, inclusive=(True, True))
        for t in keys:
            theG = self._snapshots[t]
            try:
                answer = function(theG, **kwargs)
            except Exception as e:
                print(str(e))
                answer = Exception
            to_return.append(answer)

        return to_return

    def to_DynGraphIG(self):
        """
        Convert the graph into a DynGraph_IG.

            ##Can be optimized !

        :return:
        """
        sn_duration=self.frequency()

        by_edges = dict()
        by_nodes = dict()
        for t,g in self.snapshots().items():
            for e in g.edges():
                by_edges.setdefault(frozenset(e),[]).append(t)
            for n in g.nodes():
                by_nodes.setdefault(n,[]).append(t)

        by_edges = {tuple(e):tn.Intervals.from_time_list(v,sn_duration) for e,v in by_edges.items()}
        by_nodes = {n:tn.Intervals.from_time_list(v,sn_duration) for n,v in by_nodes.items()}
        return tn.DynGraphIG(start=min(self.snapshots_timesteps()),end=max(self.snapshots_timesteps())+sn_duration,edges=by_edges,nodes=by_nodes,frequency=sn_duration)

    def to_DynGraphLS(self):
        """
        Convert to a linkstream

        Currently, conserve only edges
        :return:
        """
        sn_duration = self.frequency()

        by_edges = dict()
        by_nodes = dict()

        for t, g in self.snapshots().items():
            for e in g.edges():
                by_edges.setdefault(frozenset(e), []).append(t)
            for n in g.nodes():
                by_nodes.setdefault(n, []).append(t)

        by_edges = {tuple(e): v for e, v in by_edges.items()}
        by_nodes = {n: tn.Intervals.from_time_list(v, sn_duration) for n, v in by_nodes.items()}
        return tn.DynGraphLS(start=min(self.snapshots_timesteps()),end=max(self.snapshots_timesteps())+sn_duration,edges=by_edges,nodes=by_nodes,frequency=sn_duration)

    # def to_DynGraphIG(self, sn_duration, convert_time_to_integer=False):
    #     """
    #     Convert the graph into a DynGraph_IG.
    #
    #     ##Can be optimized !
    #
    #     By default, snapshots last from their time ID to the time ID of the next snapshot.
    #     Be careful, for the last snaphsot, we cannot know his duration, therefore, if sn_duration is not provided, it has a default duration equal to the min
    #     of all durations
    #
    #     :param sn_duration: duration of sns, None for automatic behavior
    #     :param convert_time_to_integer: if True, use the snapshot order in the list of SN rather than its time step
    #     :return:
    #     """
    #     toReturn = tn.DynGraphIG()
    #
    #
    #     for i in range(len(self._snapshots)):
    #         if convert_time_to_integer:
    #             current_t=i
    #             tNext=i+1
    #         else:
    #             current_t = self._snapshots.peekitem(i)[0]
    #
    #             if sn_duration!=None:
    #                 tNext = current_t+sn_duration
    #
    #             else:
    #                 if i<len(self._snapshots)-1:
    #                     tNext=self._snapshots.peekitem(i + 1)[0]
    #                 else:
    #                     #computing the min duration to choose as duration of the last period
    #                     dates = self.snapshots_timesteps()
    #
    #                     minDuration = min([dates[i + 1] - dates[i] for i in range(len(dates) - 1)])
    #                     tNext = current_t+minDuration
    #
    #         if (len(self._snapshots.peekitem(i)[1].nodes()))>0:
    #             toReturn.add_nodes_presence_from(self._snapshots.peekitem(i)[1].nodes(), (current_t, tNext))
    #
    #             if len(list(self._snapshots.peekitem(i)[1].edges()))>0:
    #                 toReturn.add_interactions_from(list(self._snapshots.peekitem(i)[1].edges()), (current_t, tNext) )
    #
    #
    #     return toReturn

    def _combine_weighted_graphs(self,graphList, weight=1.0):
        """
        Function to aggregate several graphs into a weighted graph

        :param graphList: enumerable of graphs
        :param weight: default weight
        :return:
        """
        newG = nx.Graph()
        for g in graphList:
            e_weights = nx.get_edge_attributes(g, "weight")
            n_weights = nx.get_node_attributes(g,"weight")

            # add weight of one to unweighted graphs
            if len(e_weights) == 0:
                e_weights = {(u, v): 1 for (u, v) in g.edges()}
            if len(n_weights) == 0:
                n_weights = {n: 1 for n in g.nodes()}

            for n, w in n_weights.items():
                if newG.has_node(n):
                    newG.nodes[n]["weight"] += w
                else:
                    newG.add_node(n, weight=w)

            for (u, v), w in e_weights.items():
                if newG.has_edge(u, v):
                    newG[u][v]["weight"] += w
                else:
                    newG.add_edge(u, v, weight=w)
        return newG

    def cumulated_graph(self,times=None):
        """
        Compute the cumulated graph.

        Return a networkx graph corresponding to the cumulated graph of the given period (whole graph by default)

        :param times: list/set of time steps ID of snapshots to cumulate. Default (None) means all snapshots
        :return: a networkx (weighted) graph
        """

        if times==None:
            snapshots=list(self.snapshots().values())
        else:
            snapshots= [sn for t,sn in self.snapshots().items() if t in times]

        return self._combine_weighted_graphs(snapshots)

    def slice(self,start,end):
        """
        Keep only the selected period

        :param start: time of the beginning of the slice
        :param end: time of the end of the slice
        """

        to_return = tn.DynGraphSN()
        interv = tn.Intervals((start,end))
        for t in list(self._snapshots.keys()):
            if interv.contains_t(t):
                to_return.add_snapshot(t,self._snapshots[t])
        return to_return

    def change_times(self):
        """
        Times of non-empty snapshots

        :return: list of times
        """
        return [t for t,g in self.snapshots().items() if len(g.edges()) > 0]

    def aggregate_sliding_window(self, bin_size=None, shift=None, t_start=None, t_end=None,weighted=True):
        """
        Return a new dynamic graph without modifying the original one, aggregated using sliding windows of the desired size. If Shift is not provided or equal to bin_size, windows are non overlapping.
        If no parameter is provided, creates a single graph aggregating the whole period.
        Yielded graphs are weighted (weight: number of apparition of edges during the period)

        :param bin_size: desired size of bins, in the internal time unit (not necessarily equals to the number of snapshot_affiliations)
        :param shift: time distance (shift) between the start of two successive bins, in the internal time unit (not necessarily number of sn)
        :param t_start: time step to start the binning (default: first)
        :param t_end: time step (not included) to stop the binning (default: last)
        :return: a DynGraph_SN object
        """

        if t_start==None:
            t_start = self.snapshots_timesteps()[0]
        if t_end==None:
            t_end=self.snapshots_timesteps()[-1]

        if bin_size == None:
            bin_size=t_end-t_start

        if shift==None:
            shift=bin_size

        bins = []
        for t in range(t_start, t_end, shift):
            bins.append((t, t + bin_size))

        toReturn = DynGraphSN(frequency=bin_size)
        for (binStart,binEnd) in bins:
            keys = self.snapshots().irange(binStart, binEnd, inclusive=(True, False))
            keys = list(keys)
            if len(keys)>0:
                if weighted:
                    toReturn.add_snapshot(binStart, self._combine_weighted_graphs([self._snapshots[k] for k in keys]))
                else:
                    all_edges = [list(self._snapshots[k].edges()) for k in keys]
                    all_edges = list(set([e for edge_list in all_edges for e in edge_list ]))
                    toReturn.add_snapshot(binStart,nx.Graph(all_edges))
                    #toReturn.add_snapshot(binStart,nx.compose_all([self._snapshots[k] for k in keys]))

            else:
                toReturn.add_snapshot(binStart)

        toReturn.discard_empty_snapshots()
        return toReturn

    def _get_monday_from_calendar_week(self, year, calendar_week):
        monday = datetime.datetime.strptime(f'{year}-{calendar_week}-1', "%Y-%W-%w").date()
        return monday

    def _date_threasholded(self,date,period):
        date = date.replace(second=0)
        if period == "minute":
            return date
        date = date.replace(minute=0)
        if period == "hour":
            return date
        date = date.replace(hour=0)
        if period == "day":
            return date

        if period == "week":
            week = date.isocalendar()[1]
            temp = datetime.combine(self._get_monday_from_calendar_week(date.year, week), datetime.min.time())
            return temp

        date = date.replace(day=0)

        if period == "month":
            return date
        date = date.replace(month=0)
        if period == "year":
            return date

    def aggregate_time_period(self, period, step_to_datetime=datetime.utcfromtimestamp):
        """
        Aggregate graph by time period (day, year, ...)

        Return a new dynamic graph without modifying the original one, aggregated such as all
        Yielded graphs are weighted (weight: number of apparition of edges during the period)

        :param period: either a string (minute,hour,day,week,month,year) or a function returning the timestamp truncated to the start of the desired period
        :param step_to_datetime: function to convert time step to a datetime object, default is utfromtimestamp
        :return: a DynGraph_SN object
        """
        to_return = tn.DynGraphSN()

        if isinstance(period,str):
            period_func = lambda x : self._date_threasholded(x,period)
        else:
            period_func = period

        for t,g in self.snapshots().items():
            t_date = step_to_datetime(t)

            new_t = period_func(t_date)

            #new_t = new_t.timestamp()
            new_t = int(new_t.replace(tzinfo=timezone.utc).timestamp())
            if not new_t in to_return.snapshots():
                to_return.add_snapshot(new_t, g)
            to_return._snapshots[new_t]=self._combine_weighted_graphs([to_return.snapshots(new_t),self.snapshots(t)])
        to_return.frequency(to_return.snapshots_timesteps()[1]-to_return.snapshots_timesteps()[0])
        to_return.discard_empty_snapshots()
        return to_return

    def snapshots(self, t=None):
        """
        Return all snapshots or a particular one

        Default: return a  sorted dictionary, key: the time information, value: a networkx graph.
        If t is provided, return graph at that particular time

        :param t: the time of the snapshot to return
        :return:
        """
        if t==None:
            return self._snapshots
        return self._snapshots[t]

    def node_presence(self, nodes=None):
        """
        Presence time of nodes

         Several usages:

        * If nodes==None (default), return a dict for each note, its existing times
        * If nodes is a single node, return the interval of presence of this node
        * If nodes is a set of nodes, return interval of presence of those nodes as a dictionary

        :param nodes: list of ndoes
        :return: a dictionary, key:node, value: list of time steps
        """

        if isinstance(nodes,str):
            nodes= {nodes}

        toReturn = {}
        for (SNt,g) in self.snapshots().items():
            if nodes==None:
                nodes_this_step = g.nodes()
            else:
                nodes_this_step = g.nodes() & nodes
            for n in nodes_this_step:
                toReturn.setdefault(n,[])
                toReturn[n].append(SNt)
        if nodes!=None and len(nodes)==1:
            return toReturn[list(nodes)[0]]
        return toReturn

    def edge_presence(self, edges=None):
        """
        Presence time of edges

         Several usages:

        * If edges==None (default), return a dict for each edge, its existing times
        * If edges is a set of edges, return interval of presence of those edges as a dictionary

        :param edges: list of edges
        :return: a dictionary, key:edge(pair), value: list of time steps
        """

        if edges!=None:
            list_element_example = list(edges)[0]
            if not isinstance(list_element_example, Iterable) or isinstance(list_element_example, str):
                edges = [edges]

        toReturn = {}
        for (SNt, g) in self.snapshots().items():
            if edges == None:
                edges_this_step = g.edges()
            else:
                edges_this_step = g.edges() & edges
            for e in edges_this_step:
                ef = frozenset([e[0],e[1]])
                toReturn.setdefault(ef, [])
                toReturn[ef].append(SNt)

        if edges != None:
            to_return2 = {}
            for e in edges:
                e2=frozenset(e)
                to_return2[e2]=toReturn[e2]
            if len(to_return2)==1:
                return list(to_return2.values())[0]
            return to_return2

        return toReturn

    def to_tensor(self,always_all_nodes=True):
        """
        Return a tensor representation

        Compute the list of matrices corresponding to each graph, with nodes ordered in a same order
        And the dic of nodes corresponding
        and the list for each sn of nodes
        :param always_all_nodes: if True, even if a node is not active during a snapshot, it is included in the matrix
        :return: 3 elements:(A,B,C) A: list of numpy matrices, B: a bidictionary {node name:node order in the matrix}, C: active node at each step, as a list of list of nodes
        """
        allNodes = list(self.aggregate().nodes().keys())
        nodeIDdict = bidict()
        for i in range(len(allNodes)):
            nodeIDdict[allNodes[i]] = i + 1

        # Create a dynamic network as sequence of nx graphs
        Gs = list(self.snapshots().values())
        nodesPresent = []
        GsMat = []

        nodeIdOrderedList = list(nodeIDdict.keys())
        for g0 in Gs:
            g2 = g0.copy()
            # get nodes of the current graph ordered according to the global order
            filteredOrderedNodes = [x for x in nodeIDdict.keys() if x in g2.nodes]

            # transform to numpy matrix
            if(always_all_nodes):
                g2.add_nodes_presence_from(nodeIdOrderedList)
            GsMat.append(nx.to_numpy_matrix(g2, nodelist=nodeIdOrderedList).tolist())
            nodesPresent.append([nodeIDdict[name] for name in filteredOrderedNodes])
        return(GsMat,nodeIDdict,nodesPresent)

    def last_snapshot(self):
        """
        Return the last snapshot

        :return: the last snapshot as a networkx graph
        """
        return self.snapshots()[self.snapshots_timesteps()[-1]]

    def full_copy(self):
        return deepcopy(self)

    def normalize_to_integers(self,nodes_start_at=1,time_start_at=1):
        """
        Transform time IDs and nodes to Integer

        :return: a new dynamic graph object, a dictionary of nodes {originalID:newID} and a dictionary of times {originalID:newID}
        """
        all_nodes = set()
        to_return = tn.DynGraphSN()
        for g in self.snapshots().values():
            all_nodes.update(set(g.nodes))
        nodes_dict = {v: (i+nodes_start_at) for i, v in enumerate(all_nodes)}

        for i,g in enumerate(self.snapshots().values()):
            to_return.add_snapshot(i+time_start_at,nx.relabel_nodes(g, nodes_dict))

        times = list(self.snapshots().keys())
        time_dict = {i+time_start_at:times[i] for i in range(len(times))}
        nodes_dict_inv = {v:k for k,v in nodes_dict.items()}

        return to_return,nodes_dict_inv,time_dict

    def code_length(self,as_matrix=True,as_edgelist=True):
        to_return =[]
        to_return.append(code_length_SN_M(self))
        to_return.append(code_length_SN_E(self))

        return tuple(to_return)
    # def code_length(self,as_ls=False):
    #     if as_ls:
    #         return self.code_length_ls()
    #     else:
    #         return self.code_length_sn()

    # def code_length_sn(self):
    #     """
    #     Time info should not be repeated
    #     2 versions:
    #     - edges are little repeated => each interaction encoded explicitly
    #     - edges are strongly repeated => matrix form
    #     :return:
    #     """
    #     #we ignore the code of edges identical for all
    #     g_cumulated = self.cumulated_graph()
    #     node_encoding = np.log2(len(g_cumulated.nodes()))
    #     edge_encoding = node_encoding*2
    #     nb_time = len([x for x in self.snapshots().values() if len(x.edges())>0])
    #     time_encoding =  np.log2(nb_time)# +1 for stop
    #
    #     nb_unique_edges = len(g_cumulated.edges())
    #     nb_interactions = sum([len(x.edges()) for x in self.snapshots().values()])
    #
    #     #if "matrix_form":
    #     #matrix form
    #     total_code_matrix = nb_unique_edges*nb_time+edge_encoding*nb_unique_edges+time_encoding*nb_time
    #     print("sn_m: ",edge_encoding,time_encoding,nb_unique_edges,nb_time)
    #     #else:
    #     print("sn_e: ",edge_encoding,time_encoding,nb_interactions,nb_time)
    #     #T1_(N1,N2)_(N3,N4)_STOP_T2_...
    #     total_code_edges = nb_interactions*edge_encoding + nb_time*time_encoding + nb_time*edge_encoding
    #     return total_code_matrix,total_code_edges

    # def _compute_encoding_properties(self):
    #
    #     nb_nodes = len(g_cumulated.nodes())
    #     nb_time = len([x for x in self.snapshots().values() if len(x.edges()) > 0])
    #     #time_encoding = np.log2(len(self.snapshots()))  # +1 for stop
    #     nb_unique_edges = len(g_cumulated.edges())
    #     nb_interactions = sum([len(x.edges()) for x in self.snapshots().values()])
    #     return(nb_nodes,nb_unique_edges,nb_interactions,nb_time)

    def code_length_update(self):
        g_cumulated = self.cumulated_graph()
        node_encoding = np.log2(len(g_cumulated.nodes()))
        edge_encoding = node_encoding * 2
        nb_time = len([x for x in self.snapshots().values() if len(x.edges()) > 0])
        nb_unique_edges = len(g_cumulated.edges())
        #nb_interactions = sum([len(x.edges()) for x in self.snapshots().values()])
        nb_changes = 0
        ordered_graphs = list(self.snapshots().values())
        for i,g in enumerate(ordered_graphs):
            if i==0:
                nb_changes+=len(g.edges())
            else:
                #modified_edges = nx.symmetric_difference(ordered_graphs[i-1],ordered_graphs[i])
                some_edges_1 = {tuple(sorted(edge)) for edge in ordered_graphs[i-1].edges()}
                some_edges_2 = {tuple(sorted(edge)) for edge in ordered_graphs[i].edges()}
                changes= some_edges_1.symmetric_difference(some_edges_2)

                nb_changes+=len(changes)

        time_encoding = np.log2(nb_time)

        #T1_(N1,N2)_(N3,N4)_STOP_T2_CHANGE...
        total_code_edges = nb_changes*edge_encoding + nb_time*time_encoding + nb_time*edge_encoding
        print("updates: ",edge_encoding,time_encoding,nb_changes,nb_unique_edges)

        return total_code_edges

    def stability(self):
        """
        Fraction of successive edge appearances that are in adjacent snasphots
        :return:
        """
        times = self.snapshots_timesteps()
        edge_presence = self.edge_presence()
        successions = 0
        adjacent=0
        for n,presences in edge_presence.items():
            successions += len(presences)-1
            for i in range(len(presences)-1):
                if times.index(presences[i])==times.index(presences[i+1])-1:
                    adjacent+=1
        print(adjacent,successions)
        if successions==0:
            return 1
        return adjacent/successions

    def synchronicity(self):
        """
        Fraction of edges life that appear

        Careful, naive version currently !!!

        for each edge, over its period of existence
        for each other edge, how many appearance, thus how many simultaneous possible
        and how many effective

        -Must take into account that nodes exist over period: if their lifetime do not overlap they are not supposed to be
        synchronous
        -Must take into account that

        :return:
        """

        # times = self.snapshots_timesteps()
        # edge_presence = self.edge_presence()
        # active = 0
        # sn_during_life = 0
        # for n, presences in edge_presence.items():
        #     active += len(presences)
        #     sn_during_life += times.index(presences[-1])- times.index(presences[0])+1
        # print(active, sn_during_life)
        # return active / sn_during_life

        times = self.snapshots_timesteps()
        non_zero = [t for t in times if len(self.snapshots(t).edges())>0]

        edge_presence = self.edge_presence()
        nb_different_edges = len(edge_presence)

        total_interactions = sum([len(x) for x in edge_presence.values()])
        avg_inter_by_sn = total_interactions/len(non_zero)
        #max_possible_inter_by_sn = nb_different_edges
        max_possible_inter_by_sn=total_interactions/max([len(x) for x in edge_presence.values()])
        print("synchro",avg_inter_by_sn,max_possible_inter_by_sn)
        return avg_inter_by_sn/max_possible_inter_by_sn

    def write_interactions(self,filename):
        """
        Write interactions in a file

        Write in corresponding json format

        :param filename:
        :return:
        """
        tn.write_as_SN_E(self,filename)

    def discard_empty_snapshots(self):
        """
        Discard snapshots with no edges

        """
        to_remove=[]
        for t,g in self.snapshots().items():
            if len(g.edges())==0 and len(g.nodes())==0:
                to_remove.append(t)
        for t in to_remove:
            self.remove_snapshot(t)

max_stars_repo_name: Yquetzal/tnetwork

max_issues_repo_name: Yquetzal/tnetwork

max_forks_repo_name: Yquetzal/tnetwork

reason: import numpy,import networkx


text:
#!/usr/bin/env python3
# Copyright 2019 Christian Henning, Maria Cervera
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
r"""
A spiking layer module that maintains its own parameters (:mod:`lib.spiking_layer`)
-----------------------------------------------------------------------------------

The module :mod:`lib.spiking_layer` contains the implementation of a single
spiking layer. The goal is to utilize the custom
``Functions`` implemented in module :mod:`lib.spiking_functions`
and to provide a wrapper that takes care of managing the parameters
(:math:`W`) of such a layer. The layers defined here will then be used in
:mod:`lib.snn` to define a multi-layer spiking network.

In biological networks, the electrical activity from a pre-synaptic spike leads
to changes in the membrane potential of a post-synaptic neuron.
In nature, this is a process involving many ions and channels.
Here we will use a simplified model: the leaky integrate-and-fire
model for spiking neurons. Such models are composed of 1) a description of the
dynamics of the membrane potential, and 2) a mechanism for triggering spikes.

In our implementation, the dynamics of the membrane potential, together with
the dynamics of the current and spiking variables, are updated at each timestep
based on a discrete implementation of a set of differential equations described
below.

The ODEs that you will implement sequentially update the membrane potentials
:math:`U`, the auxiliary variable for the alpha function :math:`H`, and the
current :math:`I`. The equation for the membrane potential of neuron :math:`i`
is:

.. math::
    \frac{dU_i}{dt} =& - \frac{1}{\tau_{mem}} \left[ (U_i - U_{rest}) - R I_i
    \right] + S_i(t) \left( U_{rest} - U_{threshold} \right) \\
    :label: eq-U-ODE

Most commonly, post-synaptic currents resulting from spiking inputs from
pre-synaptic neurons are modeled as exponential decay functions, where a spike
causes an instantaneous increase in the post-synaptic membrane potential, which
then decays exponentially with time i.e. :math:`u(t) = e^{-t}`, assuming the
pre-synaptic spike occurs at :math:`t=0`, and a resting potential of 0.
A more biologically plausible model is an alpha-shaped post-synaptic current,
where the post-synaptic current following a pre-synaptic spike has a finite
rise time. In this case we would have :math:`u(t) = te^{-t}`. In this tutorial,
we ask you to implement alpha-shaped post-synaptic currents by filling the
methods :meth:`lib.spiking_layer.update_H` and
:meth:`lib.spiking_layer.update_I`. These are based on the following equations:

.. math::
    \frac{dH_i}{dt} =& - \frac{1}{\tau_{rise}} H_i (t) + \sum_j W_{ij} S_j (t)
    :label: eq-H-ODE

.. math::
    \frac{dI_i}{dt} =& - \frac{1}{\tau_{syn}}  I_i (t) + H_i (t)
    :label: eq-I-ODE

All equations were derived during the tutorial session, and can be found in
the tutorial slides. Note, however, that only a discrete version for the case
of exponential-shaped post-synaptic currents was derived, and not for the
alpha-shaped case, which is the one you need here. Therefore you will need to
figure out how to turn the ODEs into code appropriately.

.. autosummary::

    lib.spiking_layer.SpikingLayer
    lib.spiking_layer.SpikingLayer.update_U
    lib.spiking_layer.SpikingLayer.update_I
    lib.spiking_layer.SpikingLayer.update_H
    lib.spiking_layer.SpikingLayer.forward

.. _`Neftci et al. (2019)`:
    https://arxiv.org/pdf/1901.09948.pdf

"""

# Please fill out identification details in
#   lib/spiking_layer.py
#   lib/spiking_functions.py
#   theory_question.txt

# Student name : Anian Ruoss
# Student ID   : 16-934-002
# Email address: anruoss@ethz.ch

import lib.spiking_functions as sf
import numpy as np
import torch
import torch.nn as nn

class SpikingLayer(nn.Module):
    r"""Implements a single spiking layer.

    The :class:`lib.SpikingLayer` class contains all the parameters and
    variables necessary to implement a single spiking layer. It will be a
    submodule of the ``torch.nn.Module`` instance named :class:`lib.snn.SNN`.
    You will use the class :class:`lib.SpikingLayer` to define your hidden
    layer and your output layer.

    Attributes:
        tau_mem (float): Membrane time constant.
        tau_syn (float): Synaptic time constant.
        tau_rise (float): Rise synaptic time constant.
        u_rest (float): Resting membrane potential.
        u_threshold (float): Firing threshold.
        R (float): Membrane resistance.
        delta_t (float): The time step duration.
        gamma (float): Decay rate for the post-synaptic current :math:`I`.
        beta (float): Decay rate for the membrane potential :math:`U`.
        phi (float): Decay rate for the auxiliary current variable :math:`H`.
        weights (torch.nn.Parameter): The weight matrix :math:`W` of the layer.
        compute_spikes (func): Spike non-linearity function.

    Args:
        in_features (int): Size of the pre-synaptic layer.
        out_features (int): Size of the current layer.
        args (argparse.Namespace): The command-line arguments.

    """

    def __init__(self, in_features, out_features, args):
        nn.Module.__init__(self)

        # Store parameters of the layer.
        self.tau_mem = args.tau_mem
        self.tau_syn = args.tau_syn
        self.tau_rise = args.tau_rise
        self.u_rest = args.u_rest
        self.u_threshold = args.u_threshold
        self.R = args.R
        self.delta_t = args.delta_t

        # Calculate the decay scales.
        self.gamma = float(np.exp(-args.delta_t/args.tau_syn))
        self.beta = float(np.exp(-args.delta_t/args.tau_mem))
        self.phi = float(np.exp(-args.delta_t/args.tau_rise))

        # Define which spike nonlinearity function to use.
        self.compute_spikes = sf.spike_function

        # Create and initialize weights.
        self._weights = nn.Parameter(torch.Tensor(out_features, in_features),
                                     requires_grad=True)
        nn.init.normal_(self._weights, mean=0.0,
                        std=args.weight_scale/np.sqrt(in_features))

    @property
    def weights(self):
        r"""Getter for read-only attribute :attr:`weights`."""
        return self._weights

    def update_U(self, U, I, S):
        r"""Updates the membrane potential.

        Updates the membrane potential given the current and past states
        of the network as specified in eq. :eq:`eq-U-ODE`.

        Note that the :class:`lib.SpikingLayer` has an attribute called
        :math:`beta` that will be useful here.

        Args:
            U: The membrane potential :math:`U`.
            I: The post-synaptic current :math:`I`.
            S: The spiking activity :math:`S`.

        Returns:
            The updated membrane potential of the neurons in the layer.

        """
        return self.beta * (U - self.u_rest) + self.delta_t * self.R * I / self.tau_mem - S * (self.u_threshold - self.u_rest) + self.u_rest

    def update_H(self, H, inputs):
        r"""Updates the state of the auxiliary variable for alpha-shaped
        post-synaptic current.

        Implementation of eq. :eq:`eq-H-ODE`. Please note that the pre-synaptic
        inputs have already been multiplied by the weights at the beginning of
        the forward method. Therefore, you should only make sure you give the
        correct inputs argument to this function.

        Note that the :class:`lib.SpikingLayer` has an attribute called
        :math:`phi` that will be useful here.

        Args:
            H: The auxiliary variable for the alpha-shaped post-synaptic
                currents :math:`H`.
            inputs: The inputs (weighted spikes) to the layer in the current
                time step.

        Returns:
            The updated auxiliary current variable for the neurons in the
            layer.

        """
        return self.phi * H + inputs

    def update_I(self, I, H):
        r"""Updates the post-synaptic current.

        Updates the post-synaptic current given the current and past states
        of the network as specified in eq. :eq:`eq-I-ODE`.

        Note that the :class:`lib.SpikingLayer` has an attribute called
        :math:`gamma` that will be useful here.

        Args:
            I: The post-synaptic current :math:`I`.
            H: The auxiliary variable for the alpha-shaped post-synaptic
                currents :math:`H`.

        Returns:
            The updated post-synaptic current of the neurons in the layer.

        """
        return self.gamma * I + self.delta_t * H

    def forward(self, X):
        r"""Computes the output activation of a spiking layer.

        This method computes the membrane potential and spiking activity of the
        current layer across all time steps given the pre-synaptic spiking
        activity. For this, the state of the layer is
        updated time step by time step; i.e. the post-synaptic current,
        membrane potential and spiking activity are computed in each time step.
        The states are updated based on the computational graph provided in
        Figure 2 in `Neftci et al. (2019)`_, and when filling in the missing
        lines you should pay extra attention and make sure that the values
        you provide to the update methods belong to the right time step
        according to this computational graph.

        Note that since we deal with alpha-shaped post-synaptic currents here
        (and not exponential decay post-synaptic currents), the computational
        graph has an extra variable :math:`H` that is updated based on the
        inputs in the previous time step, and its own value in the previous
        time step. :math:`H` in a given time step is then used to compute the
        post-synaptic current :math:`I` in the following time step. Notice that
        for this extra equation, the :class:`lib.SpikingLayer` class has an
        attribute :math:`phi` that governs the decay rate of the variable
        :math:`H`. For further discussion see :mod:`lib.spiking_layer`.

        Args:
            X: The spiking activity of the previous layer.

        Returns:
            (tuple): Tuple containing:

            - **U**: The membrane potential for all neurons of the layer in all
              time steps.
            - **S**: The spiking activity of all neurons of the layer in all
              time steps.

        """

        device = X.device
        dtype = torch.float

        num_h = self.weights.shape[0]  # number of neurons in the layer
        N = X.shape[1]
        batch_size = X.shape[0]

        # Compute the total input to the layer in all time steps
        inputs = torch.einsum("abc,cd->abd", (X, self.weights.t()))

        # Create lists to store the network states in each time step
        U = [[]]*N  # membrane potential
        S = [[]]*N  # spiking activity
        H = [[]]*N  # auxiliary current state
        I = [[]]*N  # post-synaptic current

        # Initialize states
        U[0] = self.u_rest * torch.ones((batch_size, num_h), device=device,
                                        dtype=dtype)
        S[0] = torch.zeros_like(U[0])
        H[0] = torch.zeros_like(U[0])
        I[0] = torch.zeros_like(U[0])

        # Compute hidden layer state at each time step
        for n in range(1, N):

            # Compute the post-synaptic current
            H[n] = self.update_H(H[n - 1], inputs[:, n - 1])
            I[n] = self.update_I(I[n - 1], H[n - 1])

            # Compute the membrane potential
            U[n] = self.update_U(U[n - 1], I[n - 1], S[n - 1])

            # Compute the spiking activity
            S[n] = self.compute_spikes(U[n] - self.u_threshold)

        U = torch.stack(U, dim=1)
        S = torch.stack(S, dim=1)

        return U, S

if __name__ == '__main__':
    pass

ext: py

max_stars_repo_name: anianruoss/LDABNN

max_issues_repo_name: anianruoss/LDABNN

max_forks_repo_name: anianruoss/LDABNN

reason: import numpy


text:
from geo3d import Vector, Point, express_vector_in_frame
import numpy as np
import pytest

def test_vector_from_array():
    p = Vector([1, 2, 3])
    assert p[2] == 3

def test_vector_from_vector():
    p = Vector(Vector([1, 2, 3]))
    assert p[2] == 3

def test_vector_from_point():
    p = Vector(Point([1, 2, 3]))
    assert p[2] == 3

def test_vector_from_tuple():
    p = Vector((1, 2, 3))
    assert p[2] == 3

def test_vector_from_np():
    p = Vector(np.array([1, 2, 3]))
    assert p[2] == 3

def test_vector_from_np_bare_nocopy():
    a = np.array([1, 2, 3])
    p = Vector.from_array(a, copy=False)
    a[2] = 7
    assert p[2] == 7

def test_vector_from_np_bare_copy():
    a = np.array([1, 2, 3])
    p = Vector.from_array(a, copy=True)
    a[2] = 7
    assert p[2] == 3

def test_express_vector_in_frame(example_frames):
    fa, fb, fc = example_frames
    v0 = (1, 3, 0)
    v1 = (2.82843, -1.41421, 0)

    res0 = Vector(v0).express_in_frame(
        fa, original_frame=fb
    )
    assert res0.as_array() == pytest.approx(v1, abs=1e-4)
    assert isinstance(res0._a, np.ndarray)

    assert express_vector_in_frame(
        v0, fa, original_frame=fb
    ).as_array() == pytest.approx(v1, abs=1e-4)

    assert Vector(v0).express_in_frame(fa, original_frame=fb) == Vector(
        v0
    ).express_in_frame(fa.express_in_frame(fb))


text:
# v2 -- impact parameter; Gaussian priors on a, p, rp
# v3 -- binning, timing, TDB correction deprecated, save results,
#		resolved issue of sigma > 1, naming scheme

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid.anchored_artists import AnchoredText
from scipy.stats import binned_statistic
import batman
import emcee
import corner
import os, sys, time

'''
# ADJUST THINGS IN THIS SECTION. BE VERY CAREFUL!
##########################################################
targ_name = 'K2-79'
epic = 210402237
save_results = True	# Save plots?
show_plots = False	# Show plots?
annotate_plot = True

bjd_offset = 57710.		# Makes nicer plots

# Define initial parameters
param_names = ["$T_0$", "$R_p/R_*$", "$a/R_*$", "$b$", "$offset$", "$\sigma_f$"]
t0_i = 57718.93		# BJD (look at light curve)
b_i = 0.1			# degrees
offset_i = 0.005	# offset added to model
sigma_i = 0.001		# flux uncertainty

# From Crossfield et al. 2016
rp_i = 0.0261		# Rp/R*
s_rp = 0.0012
per_i = 10.99573	# days
s_per = 0.00070
Rstar = 1.28		# Rsun
s_Rstar = 0.12
Mstar = 1.06		# Msun
s_Mstar = 0.04

# From Claret et al. 2012/13
u1 = 0.0764		# LD coeff 1
u2 = 0.1004		# LD coeff 2

# Number of Bins
nbins = 100

# MCMC parameters
nsteps = 10000
burn_in = 2500
ndim = 6
nwalkers = 100
##########################################################
'''

# ADJUST THINGS IN THIS SECTION. BE VERY CAREFUL!
##########################################################
targ_name = 'K2-36'
epic = 201713348
save_results = True	# Save plots?
show_plots = False	# Show plots?
annotate_plot = True

bjd_offset = 57460.		# Makes nicer plots

# Define initial parameters
param_names = ["$T_0$", "$R_p/R_*$", "$a/R_*$", "$b$", "$offset$", "$\sigma_f$"]
t0_i = 57464.45		# BJD (look at light curve)
b_i = 0.5			# impact parameter
offset_i = 0.00005	# offset added to model
sigma_i = 0.001		# flux uncertainty

# From Crossfield et al. 2016
rp_i = 0.0323		# Rp/R*
s_rp = 0.0028
per_i = 5.34072		# days
s_per = 0.00011
Rstar = 0.75		# Rsun
s_Rstar = 0.02
Mstar = 0.81		# Msun
s_Mstar = 0.02

# From Claret et al. 2012/13
u1 = 0.0965		# LD coeff 1
u2 = 0.1144		# LD coeff 2

# Number of Bins
nbins = 500

# MCMC parameters
nsteps = 10000
burn_in = 2500
ndim = 6
nwalkers = 100
##########################################################

############## THE CODE STARTS HERE ##############
##################################################

start = time.time()

a_i = (((per_i/365.)**2) / Mstar)**(1./3.) * Rstar * 215.035216		# semimajor axis from period and stellar params
s_a = a_i * np.sqrt((s_Rstar / Rstar)**2 + ((2*s_per) / (3*per_i))**2 + (s_Mstar / (3*Mstar))**2)

# Priors.
def lnprior(theta, t0_i=t0_i, mu_rp=rp_i, sigma_rp=s_rp, mu_per=per_i, sigma_per=s_per, mu_a=a_i, sigma_a=s_a):
	t0, rp, a, b, offset, sigma = theta
	if (t0_i-0.05 < t0 < t0_i+0.05) \
	and (0. < rp) \
	and (b < a) \
	and (0. <= b < 1.) \
	and (-0.1 < offset < 0.1) \
	and (0.0001 < sigma < 0.1):
		Pr_rp = (1. / (np.sqrt(2.*np.pi)*sigma_rp)) * np.exp((-1. / (2. * sigma_rp**2.)) * (rp - mu_rp)**2.)		# Gaussian priors on Rp and a
		Pr_a = (1. / (np.sqrt(2.*np.pi)*sigma_a)) * np.exp((-1. / (2. * sigma_a**2.)) * (a - mu_a)**2.)
		return np.log(Pr_rp + Pr_a)
	return -np.inf

# Import data from text file
data = []
with open('PhotData//%s_Spitzer_unbin_TDB.txt' % targ_name) as file:
	for line in file:
		data_str = line.split()
		data_flt = [float(i) for i in data_str]
		data.append(data_flt)
data = np.asarray(data)
data = data[data[:,0].argsort()]

t_all = []
f_all = []

for i in data[:,0]:
	t_all.append(i)
for i in data[:,1]:
	f_all.append(i)

def bin_data(t,f,nbins):
	binnedTimes = []
	binnedData = binned_statistic(t,f,bins=nbins)
	bin_edges = binnedData[1]
	for n in range(len(bin_edges) - 1):
		time = (bin_edges[n] + bin_edges[n+1]) / 2.
		binnedTimes.append(time)
	t_binned = np.asarray(binnedTimes)
	f_binned = binnedData[0]
	return t_binned, f_binned

DataBinned = bin_data(t_all,f_all,nbins)

t = DataBinned[0]
f = DataBinned[1]

# Define log of the likelihood function.
def lnlike(theta, x, y, per=per_i, u1=u1, u2=u2):
	t0, rp, a, b, offset, sigma = theta
	# Set up transit parameters.
	params = batman.TransitParams()
	params.t0 = t0
	params.per = per
	params.rp = rp
	params.a = a
	params.inc = np.arccos(b / a) * (180. / np.pi)
	params.ecc = 0.
	params.w = 0.
	params.u = [u1,u2]
	params.limb_dark = 'quadratic'
	# Initialize the transit model.
	m_init = batman.TransitModel(params, x)
	model = m_init.light_curve(params) + offset
	inv_sigma2 = 1.0 / (sigma**2)
	return -0.5*(np.sum((y-model)**2*inv_sigma2 - np.log(2*np.pi*inv_sigma2)))

# Define log of probability function.
def lnprob(theta, x, y):
	lp = lnprior(theta)
	if not np.isfinite(lp):
		return -np.inf
	return lp + lnlike(theta, x, y)

initial_params = t0_i, rp_i, a_i, b_i, offset_i, sigma_i

# Initialize walkers around maximum likelihood.
pos = [initial_params + 1e-5*np.random.randn(ndim) for i in range(nwalkers)]

# Set up sampler.
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(t, f))

# Run MCMC for n steps and display progress bar.
width = 50
for i, result in enumerate(sampler.sample(pos, iterations=nsteps)):
	n = int((width+1) * float(i) / nsteps)
	sys.stdout.write("\r{}[{}{}]{}".format('sampling... ', '#' * n, ' ' * (width - n), ' (%s%%)' % str(100. * float(i) / nsteps)))
sys.stdout.write("\n")
print 'Sampling complete!'

samples = sampler.chain

# Plot walkers as a function of step number.
walker_fig, ax = plt.subplots(ndim, sharex=True, figsize=(10, 8))
ax[-1].set_xlabel('Iteration')
for n in range(ndim):
    for k in range(nwalkers):
        ax[n].plot(samples[k, :, n], color='k', alpha=0.1, lw=1)
    ax[n].set_ylabel((param_names)[n], fontsize=9)
    ax[n].margins(0, None)
    ax[n].axvline(burn_in, color='b', alpha=0.5, lw=1, ls='--')

# Discard burn-in.
samples = samples[:, burn_in:, :].reshape((-1, ndim))

# Final params and uncertainties based on the 16th, 50th, and 84th percentiles of the samples in the marginalized distributions.
t0_mcmc, rp_mcmc, a_mcmc, b_mcmc, offset_mcmc, sigma_mcmc = map(
	lambda v: (v[1], v[2]-v[1], v[1]-v[0]), zip(*np.percentile(samples, [16, 50, 84], axis=0)))

# create array for flux uncert
ferr = []
for i in f:
	ferr.append(sigma_mcmc[0])
ferr = np.asarray(ferr)

# Plot the final transit model.
params_final = batman.TransitParams()
params_final.t0 = t0_mcmc[0]
params_final.per = per_i
params_final.rp = rp_mcmc[0]
params_final.a = a_mcmc[0]
params_final.inc = np.arccos(b_mcmc[0] / a_mcmc[0]) * (180. / np.pi)
params_final.ecc = 0.
params_final.w = 0.
params_final.u = [u1,u2]
params_final.limb_dark = "quadratic"
tl = np.linspace(min(t),max(t),5000)
m = batman.TransitModel(params_final, tl)
m_res = batman.TransitModel(params_final, t)

f_final = m.light_curve(params_final)
final_fig, ax = plt.subplots(figsize=(10,8))
ax.set_title('%s (EPIC %s)' % (targ_name, str(epic)))
ax.errorbar(t-bjd_offset,f - offset_mcmc[0],yerr=ferr,fmt='k.',capsize=0,alpha=0.4,zorder=1)
ax.plot(tl-bjd_offset, f_final, 'r-',alpha=0.8,lw=3,zorder=2)
if annotate_plot == True:
	ant = AnchoredText('$T_0 = %s^{+%s}_{-%s}$ \n $R_p/R_* = %s^{+%s}_{-%s}$' % (2400000+round(t0_mcmc[0],4),round(t0_mcmc[1],4),
	round(t0_mcmc[2],4),round(rp_mcmc[0],4),round(rp_mcmc[1],4),round(rp_mcmc[2],4)), prop=dict(size=11), frameon=True, loc=3)
	ant.patch.set_boxstyle('round,pad=0.,rounding_size=0.2')
	ax.add_artist(ant)
ax.set_xlabel("$BJD_{TDB} - 24%s$" % str(bjd_offset))
ax.set_ylabel("Relative Flux")
ax.legend(('BATMAN','Spitzer'), loc=2)

f_res = m_res.light_curve(params_final) + offset_mcmc[0]
resid_fig, axr = plt.subplots(figsize=(10,8))
axr.errorbar(t-bjd_offset,f-f_res,yerr=ferr,fmt='k.',capsize=0,alpha=0.4,zorder=1)
axr.plot(t-bjd_offset, np.zeros(len(t)), 'r-',alpha=0.8,lw=3,zorder=2)
axr.set_xlabel("$BJD_{TDB} - 24%s$" % str(bjd_offset))
axr.set_ylabel("Residual")

corn_fig = corner.corner(samples, labels=param_names)

fit_params = [t0_mcmc, rp_mcmc, a_mcmc, b_mcmc, offset_mcmc, sigma_mcmc]

exec_time = time.time() - start

if save_results == True:
	save_to = 'fig//Spitzer_plots//%s//V3//bin%s' % (targ_name,nbins)
	if not os.path.exists(save_to):
		os.makedirs(save_to)

	final_fig.savefig('%s//%s_%sb_%ss_MCMCfit_Spitzer.png' % (save_to,targ_name,nbins,nsteps), bbox_inches='tight')
	corn_fig.savefig('%s//%s_%sb_%ss_corner_Spitzer.png' % (save_to,targ_name,nbins,nsteps), bbox_inches='tight')
	walker_fig.savefig('%s//%s_%sb_%ss_walkers_Spitzer.png' % (save_to,targ_name,nbins,nsteps), bbox_inches='tight')
	resid_fig.savefig('%s//%s_%sb_%ss_residuals_Spitzer.png' % (save_to,targ_name,nbins,nsteps), bbox_inches='tight')

	with open('%s//%s_%sb_%ss_LOG_Spitzer.txt' % (save_to,targ_name,nbins,nsteps), 'w') as newfile:
		newfile.write(
			'K2_ID' + '\t' + targ_name + '\n'
			'EPIC' + '\t' + str(epic) + '\n'
			'data' + '\t' + 'Spitzer' + '\n'
			'\n'
			'nbins' + '\t' + str(nbins) + '\n'
			'nsteps' + '\t' + str(nsteps) + '\n'
			'nwalk' + '\t' + str(nwalkers) + '\n'
			'burnin' + '\t' + str(burn_in) + '\n'
			'\n'
			'per' + '\t' + str(per_i) + '\t' + str(s_per) + '\t' + '[d]' + '\t' + '(fixed)' + '\n'
			'Rstar' + '\t' + str(Rstar) + '\t' + str(s_Rstar) + '\t' + '[Rsun]' + '\t' + '(fixed)' + '\n'
			'Mstar' + '\t' + str(Mstar) + '\t' + str(s_Mstar) + '\t' + '[Msun]' + '\t' + '(fixed)' + '\n'
			'u1' + '\t' + str(u1) + '\t' + '(fixed)' + '\n'
			'u2' + '\t' + str(u2) + '\t' + '(fixed)' + '\n'
			'\n'
			)
		for i,param in enumerate(param_names):
			newfile.write(
				param[1:len(param)-1] + '\t' +
				str(fit_params[i][0]) + '\t' +
				str(fit_params[i][1]) + '\t' +
				str(fit_params[i][2]) + '\n'
				)
		newfile.write('\n\n\n' + 'runtime' + '\t' + str(round(exec_time,2)) + '\t' + '[s]')

print 't0: ',t0_mcmc
print 'rp: ',rp_mcmc
print 'a: ',a_mcmc
print 'b: ',b_mcmc
print 'offset: ',offset_mcmc
print 'sigma: ',sigma_mcmc
print 'runtime: ',exec_time

if show_plots == True:
	plt.show()

reason: import numpy,from scipy


text:
import numpy as np
import robot_sim._kinematics.jlchain as jl
import robot_sim.robots.robot_interface as ri
import rrt
import rrt_star as rrts

class XYBot(ri.RobotInterface):

    def __init__(self, pos=np.zeros(3), rotmat=np.eye(3), name='XYBot'):
        super().__init__(pos=pos, rotmat=rotmat, name=name)
        self.jlc = jl.JLChain(homeconf=np.zeros(2), name='XYBot')
        self.jlc.jnts[1]['type'] = 'prismatic'
        self.jlc.jnts[1]['loc_motionax'] = np.array([1, 0, 0])
        self.jlc.jnts[1]['loc_pos'] = np.zeros(3)
        self.jlc.jnts[1]['motion_rng'] = [-2.0, 15.0]
        self.jlc.jnts[2]['type'] = 'prismatic'
        self.jlc.jnts[2]['loc_motionax'] = np.array([0, 1, 0])
        self.jlc.jnts[2]['loc_pos'] = np.zeros(3)
        self.jlc.jnts[2]['motion_rng'] = [-2.0, 15.0]
        self.jlc.reinitialize()

    def fk(self, component_name='all', jnt_values=np.zeros(2)):
        if component_name != 'all':
            raise ValueError("Only support hnd_name == 'all'!")
        self.jlc.fk(jnt_values)

    def rand_conf(self, component_name='all'):
        if component_name != 'all':
            raise ValueError("Only support hnd_name == 'all'!")
        return self.jlc.rand_conf()

    def get_jntvalues(self, component_name='all'):
        if component_name != 'all':
            raise ValueError("Only support hnd_name == 'all'!")
        return self.jlc.get_jnt_values()

    def is_collided(self, obstacle_list=[], otherrobot_list=[]):
        for (obpos, size) in obstacle_list:
            dist = np.linalg.norm(np.asarray(obpos) - self.get_jntvalues())
            if dist <= size / 2.0:
                return True  # collision
        return False  # safe

# ====Search Path with RRT====
obstacle_list = [
    ((5, 5), 3),
    ((3, 6), 3),
    ((3, 8), 3),
    ((3, 10), 3),
    ((7, 5), 3),
    ((9, 5), 3),
    ((10, 5), 3),
    ((10, 0), 3),
    ((10, -2), 3),
    ((10, -4), 3),
    ((15, 5), 3),
    ((15, 7), 3),
    ((15, 9), 3),
    ((15, 11), 3),
    ((0, 12), 3),
    ((-1, 10), 3),
    ((-2, 8), 3)
]  # [x,y,size]
# Set Initial parameters
robot = XYBot()
rrt_s = rrt.RRT(robot)
rrts_s = rrts.RRT(robot)
import time

start_conf = np.array([15, 0])
goal_conf = np.array([5, -2.5])
start_conf = np.array([0, 0])
goal_conf = np.array([5, 10])

total_t = 0
for i in range(200):
    tic = time.time()
    path = rrt_s.plan(start_conf=start_conf, goal_conf=goal_conf,
                       obstacle_list=obstacle_list,
                       ext_dist=1, rand_rate=70, max_time=1000,
                       max_iter=100000,
                       component_name='all', smoothing_iterations=0,
                       animation=False)
    toc = time.time()
    total_t = total_t + toc - tic
print("rrt costs:", total_t)

total_t = 0
for i in range(200):
    tic = time.time()
    path = rrts_s.plan(start_state=start_conf, goal_conf=goal_conf,
                       obstacle_list=obstacle_list,
                       ext_dist=1, rand_rate=70, max_time=1000,
                       max_iter=100000,
                       component_name='all', smoothing_iterations=0,
                       animation=False)
    toc = time.time()
    total_t = total_t + toc - tic
print("rrts costs:", total_t)

max_stars_repo_name: takuya-ki/wrs

max_issues_repo_name: takuya-ki/wrs

max_forks_repo_name: takuya-ki/wrs


text:
from __future__ import print_function

import argparse
import os
import shutil
import time
import random
import cv2
import imageio

import numpy as np

import pickle
import scipy.misc

import torch
import torch.nn as nn
import torch.backends.cudnn as cudnn
import torch.optim as optim
import torch.utils.data as data

import timecycle as tc
import davis_test as davis

from scipy.ndimage.morphology import binary_dilation,generate_binary_structure
from torch.autograd import Variable

params = {}
params['filelist'] = '/home/ajabri/data/davis/val2017.txt'
# params['batchSize'] = 24
params['imgSize'] = 320
params['cropSize'] = 320
params['videoLen'] = 8
params['offset'] = 0
params['sideEdge'] = 80
params['predFrames'] = 1

def str_to_bool(v):
    if v.lower() in ('yes', 'true', 't', 'y', '1'):
        return True
    elif v.lower() in ('no', 'false', 'f', 'n', '0'):
        return False
    else:
        raise argparse.ArgumentTypeError('Boolean value expected.')

# Parse arguments
parser = argparse.ArgumentParser(description='PyTorch ImageNet Training')

# Datasets
parser.add_argument('-j', '--workers', default=4, type=int, metavar='N',
                    help='number of data loading workers (default: 4)')
parser.add_argument('--resume', default='', type=str, metavar='PATH',
                    help='path to latest checkpoint (default: none)')

parser.add_argument('--manualSeed', type=int, help='manual seed')
parser.add_argument('-e', '--evaluate', dest='evaluate', action='store_true',
                    help='evaluate model on validation set')
parser.add_argument('--pretrained', default='', type=str, metavar='PATH',
                    help='use pre-trained model')
#Device options
parser.add_argument('--gpu-id', default='0,1,2,3', type=str,
                    help='id(s) for CUDA_VISIBLE_DEVICES')
parser.add_argument('--batchSize', default=1, type=int,
                    help='batchSize')
parser.add_argument('--T', default=1.0, type=float,
                    help='temperature')

parser.add_argument('--topk_vis', default=20, type=int,
                    help='topk_vis')

parser.add_argument('--videoLen', default=8, type=int,
                    help='predict how many frames away')

parser.add_argument('--cropSize', default=320, type=int,
                    help='predict how many frames away')

parser.add_argument('--save-path', default='', type=str)

args = parser.parse_args()
state = {k: v for k, v in args._get_kwargs()}

params['batchSize'] = state['batchSize']
print('batchSize: ' + str(params['batchSize']) )

params['videoLen'] = state['videoLen']
print('videoLen: ' + str(params['videoLen']) )

params['cropSize'] = state['cropSize']
print('cropSize: ' + str(params['cropSize']) )
params['imgSize'] = state['cropSize']

# Use CUDA
os.environ['CUDA_VISIBLE_DEVICES'] = args.gpu_id
use_cuda = torch.cuda.is_available()

print(args.gpu_id)

# Random seed
if args.manualSeed is None:
    args.manualSeed = random.randint(1, 10000)
random.seed(args.manualSeed)
torch.manual_seed(args.manualSeed)
if use_cuda:
    torch.cuda.manual_seed_all(args.manualSeed)

def partial_load(pretrained_dict, model):
    model_dict = model.state_dict()

    # 1. filter out unnecessary keys
    pretrained_dict = {k: v for k, v in pretrained_dict.items() if k in model_dict}
    # 2. overwrite entries in the existing state dict
    model_dict.update(pretrained_dict)
    # 3. load the new state dict
    model.load_state_dict(model_dict)

def process_labels(lbls, height_dim, width_dim):
    # processing labels
    t00 = time.time()

    lbls = lbls[0].data.cpu().numpy()
    print(lbls.shape)

    lbls_new = []

    lbl_set = []
    lbl_set.append(np.zeros(3).astype(np.uint8))
    count_lbls = []
    count_lbls.append(0)

    for i in range(lbls.shape[0]):
        nowlbl = lbls[i].copy()
        if i == 0:
            for j in range(nowlbl.shape[0]):
                for k in range(nowlbl.shape[1]):

                    pixellbl = nowlbl[j, k, :].astype(np.uint8)

                    flag = 0
                    for t in range(len(lbl_set)):
                        if lbl_set[t][0] == pixellbl[0] and lbl_set[t][1] == pixellbl[1] and lbl_set[t][2] == pixellbl[2]:
                            flag = 1
                            count_lbls[t] = count_lbls[t] + 1
                            break

                    if flag == 0:
                        lbl_set.append(pixellbl)
                        count_lbls.append(0)

        lbls_new.append(nowlbl)

    print('lbls', time.time() - t00)

    import pdb; pdb.set_trace()

    lbl_set_temp = []
    for i in range(len(lbl_set)):
        if count_lbls[i] > 10:
            lbl_set_temp.append(lbl_set[i])

    lbl_set = lbl_set_temp
    print(lbl_set)
    print(count_lbls)

    import pdb; pdb.set_trace()
    t01 = time.time()

    lbls_resize = np.zeros((lbls.shape[0], lbls.shape[1], lbls.shape[2], len(lbl_set)))
    lbls_resize2 = np.zeros((lbls.shape[0], height_dim, width_dim, len(lbl_set)))

    for i in range(lbls.shape[0]):
        nowlbl = lbls[i].copy()
        for j in range(nowlbl.shape[0]):
            for k in range(nowlbl.shape[1]):

                pixellbl = nowlbl[j, k, :].astype(np.uint8)
                for t in range(len(lbl_set)):
                    if lbl_set[t][0] == pixellbl[0] and lbl_set[t][1] == pixellbl[1] and lbl_set[t][2] == pixellbl[2]:
                        lbls_resize[i, j, k, t] = 1

    for i in range(lbls.shape[0]):
        lbls_resize2[i] = cv2.resize(lbls_resize[i], (height_dim, width_dim))

    t02 = time.time()
    print(t02 - t01, 'relabel', t01-t00, 'label')

    return lbl_set, lbls_resize2, lbls_new

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

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

def main():
    global best_loss
    start_epoch = 0  # start from epoch 0 or last checkpoint epoch

    val_loader = torch.utils.data.DataLoader(
        davis.DavisSet(params, is_train=False),
        batch_size=int(params['batchSize']), shuffle=False,
        num_workers=args.workers, pin_memory=True)

    model = tc.TimeCycle()
    model = Wrap(model, 'forward_affinity')

    model = torch.nn.DataParallel(model).cuda()

    cudnn.benchmark = False
    print('    Total params: %.2fM' % (sum(p.numel() for p in model.parameters())/1000000.0))

    # Load checkpoint.
    if os.path.isfile(args.resume):
        print('==> Resuming from checkpoint..')
        checkpoint = torch.load(args.resume)
        partial_load(checkpoint['state_dict'], model)
        del checkpoint

    model.eval()

    if not os.path.exists(args.save_path):
        os.makedirs(args.save_path)

    print('\Testing')
    test_loss = test(val_loader, model, 1, use_cuda)

def test(val_loader, model, epoch, use_cuda):

    save_path = args.save_path + '/'
    save_file = '%s/list.txt' % save_path
    fileout = open(save_file, 'w')

    end = time.time()

    for batch_idx, (imgs_total, lbls, meta) in enumerate(val_loader):

        n_context = params['videoLen']
        finput_num = n_context

        # measure data loading time
        imgs_total = torch.autograd.Variable(imgs_total.cuda())

        bs = imgs_total.size(0)
        total_frame_num = imgs_total.size(1)
        channel_num = imgs_total.size(2)
        height_len  = imgs_total.size(3)
        width_len   = imgs_total.size(4)

        assert(bs == 1)

        folder_paths = meta['folder_path']
        print('total_frame_num: ' + str(total_frame_num))

        height_dim = int(params['cropSize'] / 8)
        width_dim  = int(params['cropSize'] / 8)

        lbl_set, lbls_resize2, lbls_new = process_labels(lbls, height_dim, width_dim)

        ##################################################################
        # Print the images
        ##################################################################

        imgs_set = imgs_total.data
        imgs_set = imgs_set.cpu().numpy()
        imgs_set = imgs_set[0]
        mean=[0.485, 0.456, 0.406]
        std=[0.229, 0.224, 0.225]

        imgs_toprint = []

        # ref image
        for t in range(imgs_set.shape[0]):
            img_now = imgs_set[t]

            for c in range(3):
                img_now[c] = img_now[c] * std[c]
                img_now[c] = img_now[c] + mean[c]

            img_now = img_now * 255
            img_now = np.transpose(img_now, (1, 2, 0))
            img_now = cv2.resize(img_now, (img_now.shape[0] * 2, img_now.shape[1] * 2) )

            imgs_toprint.append(img_now)

            imname  = save_path + str(batch_idx) + '_' + str(t) + '_frame.jpg'
            scipy.misc.imsave(imname, img_now)

        for t in range(n_context):
            nowlbl = lbls_new[t]
            imname  = save_path + str(batch_idx) + '_' + str(t) + '_label.jpg'
            scipy.misc.imsave(imname, nowlbl)

        ##################################################################
        # Compute image features
        ##################################################################

        ##################################################################
        # Compute correlation features
        ##################################################################

        now_batch_size = 4
        imgs_stack = []

        im_num = total_frame_num - n_context
        corrfeat2_set = []

        imgs_tensor = torch.Tensor(now_batch_size, finput_num, 3, params['cropSize'], params['cropSize']).cuda()
        target_tensor = torch.Tensor(now_batch_size, 1, 3, params['cropSize'], params['cropSize']).cuda()

        t03 = time.time()

        for iter in range(0, im_num, now_batch_size):

            print(iter)

            startid = iter
            endid   = iter + now_batch_size

            if endid > im_num:
                endid = im_num

            now_batch_size2 = endid - startid

            for i in range(now_batch_size2):

                imgs = imgs_total[:, iter + i + 1: iter + i + n_context, :, :, :]
                imgs2 = imgs_total[:, 0, :, :, :].unsqueeze(1)
                imgs = torch.cat((imgs2, imgs), dim=1)

                imgs_tensor[i] = imgs
                target_tensor[i, 0] = imgs_total[0, iter + i + n_context]

            corrfeat2_now = model(imgs_tensor, target_tensor)
            corrfeat2_now = corrfeat2_now.view(now_batch_size, n_context, corrfeat2_now.size(1), corrfeat2_now.size(2), corrfeat2_now.size(3))

            for i in range(now_batch_size2):
                corrfeat2_set.append(corrfeat2_now[i].data.cpu().numpy())

        t04 = time.time()
        print(t04-t03, 'model forward', t03-t02, 'image prep')

        ##################################################################
        # Label propagation
        ##################################################################
        topk_vis = args.topk_vis

        for iter in range(total_frame_num - n_context):

            if iter % 10 == 0:
                print(iter)

#             imgs = imgs_total[:, iter + 1: iter + n_context, :, :, :]
#             imgs2 = imgs_total[:, 0, :, :, :].unsqueeze(1)
#             imgs = torch.cat((imgs2, imgs), dim=1)

            corrfeat2   = corrfeat2_set[iter]
            corrfeat2   = torch.from_numpy(corrfeat2)

            out_frame_num = int(finput_num)
            height_dim = corrfeat2.size(2)
            width_dim = corrfeat2.size(3)

            corrfeat2 = corrfeat2.view(corrfeat2.size(0), height_dim, width_dim, height_dim, width_dim)
            corrfeat2 = corrfeat2.data.cpu().numpy()

            vis_ids_h = np.zeros((corrfeat2.shape[0], height_dim, width_dim, topk_vis)).astype(np.int)
            vis_ids_w = np.zeros((corrfeat2.shape[0], height_dim, width_dim, topk_vis)).astype(np.int)

            t05 = time.time()

            atten1d  = corrfeat2.reshape(corrfeat2.shape[0], height_dim * width_dim, height_dim, width_dim)
            ids = np.argpartition(atten1d, -topk_vis, axis=1)[:, -topk_vis:]
            # ids = np.argsort(atten1d, axis=1)[:, -topk_vis:]

            hid = ids // width_dim
            wid = ids % width_dim

            vis_ids_h = wid.transpose(0, 2, 3, 1)
            vis_ids_w = hid.transpose(0, 2, 3, 1)

            t06 = time.time()

            img_now = imgs_toprint[iter + n_context]

            predlbls = np.zeros((height_dim, width_dim, len(lbl_set)))
            # predlbls2 = np.zeros((height_dim * width_dim, len(lbl_set)))

            for t in range(finput_num):

                tt1 = time.time()

                h, w, k = np.meshgrid(np.arange(height_dim), np.arange(width_dim), np.arange(topk_vis), indexing='ij')
                h, w = h.flatten(), w.flatten()

                hh, ww = vis_ids_h[t].flatten(), vis_ids_w[t].flatten()

                if t == 0:
                    lbl = lbls_resize2[0, hh, ww, :]
                else:
                    lbl = lbls_resize2[t + iter, hh, ww, :]

                np.add.at(predlbls, (h, w), lbl * corrfeat2[t, ww, hh, h, w][:, None])

            t07 = time.time()
            # print(t07-t06, 'lbl proc', t06-t05, 'argsorts')

            predlbls = predlbls / finput_num

            for t in range(len(lbl_set)):
                nowt = t
                predlbls[:, :, nowt] = predlbls[:, :, nowt] - predlbls[:, :, nowt].min()
                predlbls[:, :, nowt] = predlbls[:, :, nowt] / predlbls[:, :, nowt].max()

            lbls_resize2[iter + n_context] = predlbls

            predlbls_cp = predlbls.copy()
            predlbls_cp = cv2.resize(predlbls_cp, (params['imgSize'], params['imgSize']))
            predlbls_val = np.zeros((params['imgSize'], params['imgSize'], 3))

            ids = np.argmax(predlbls_cp[:, :, 1 : len(lbl_set)], 2)

            predlbls_val = np.array(lbl_set)[np.argmax(predlbls_cp, axis=-1)]
            predlbls_val = predlbls_val.astype(np.uint8)
            predlbls_val2 = cv2.resize(predlbls_val, (img_now.shape[0], img_now.shape[1]), interpolation=cv2.INTER_NEAREST)

            # activation_heatmap = cv2.applyColorMap(predlbls, cv2.COLORMAP_JET)
            img_with_heatmap =  np.float32(img_now) * 0.5 + np.float32(predlbls_val2) * 0.5

            imname  = save_path + str(batch_idx) + '_' + str(iter + n_context) + '_label.jpg'
            imname2  = save_path + str(batch_idx) + '_' + str(iter + n_context) + '_mask.png'

            scipy.misc.imsave(imname, np.uint8(img_with_heatmap))
            scipy.misc.imsave(imname2, np.uint8(predlbls_val))

    fileout.close()

if __name__ == '__main__':
    main()

max_stars_repo_name: ajabri/vision

max_issues_repo_name: ajabri/vision

max_forks_repo_name: ajabri/vision

reason: import numpy,import scipy,from scipy


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

"""
Generate figure showing the distribution of memorization values for models
trained with two different learning rates.

Author: G.J.J. van den Burg
License: See LICENSE file.
Copyright: 2021, The Alan Turing Institute

"""

import argparse
import numpy as np

from fitter import Fitter

from analysis_utils import dict2tex

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--lr3-file",
        help="NumPy compressed archive with results for learning rate 1e-3",
        required=True,
    )
    parser.add_argument(
        "--lr4-file",
        help="NumPy compressed archive with results for learning rate 1e-4",
        required=True,
    )
    parser.add_argument(
        "-o", "--output", help="Output file to write to (.tex)", required=True
    )
    return parser.parse_args()

def make_tex(M3, M4):
    m3 = M3[:, -1]
    m4 = M4[:, -1]
    assert len(m3) == len(m4) == 60_000
    del M3, M4

    tex = []

    pgfplotsset = """
    \\pgfplotsset{compat=newest,%
        /pgf/declare function={
    		gauss(\\x) = 1/sqrt(2*pi) * exp(-0.5 * \\x * \\x);%
    		johnson(\\x,\\a,\\b) = \\b/sqrt(\\x * \\x + 1) * gauss(\\a+\\b*ln(\\x+sqrt(\\x*\\x+1)));%
    		johnsonsu(\\x,\\a,\\b,\\loc,\\scale) = johnson((\\x - \\loc)/\\scale,\\a,\\b)/\\scale;%
    	},
    }
    """

    tex.append("\\documentclass[10pt,preview=true]{standalone}")
    # LuaLaTeX
    tex.append(
        "\\pdfvariable suppressoptionalinfo \\numexpr1+2+8+16+32+64+128+512\\relax"
    )
    tex.append("\\usepackage[utf8]{inputenc}")
    tex.append("\\usepackage[T1]{fontenc}")
    tex.append("\\usepackage{pgfplots}")
    tex.append(pgfplotsset)

    tex.append("\\definecolor{MyBlue}{HTML}{004488}")
    tex.append("\\definecolor{MyYellow}{HTML}{DDAA33}")
    tex.append("\\definecolor{MyRed}{HTML}{BB5566}")

    tex.append("\\begin{document}")
    tex.append("\\begin{tikzpicture}")

    fontsize = "\\normalsize"
    bins = 100
    xmin = -15
    xmax = 35
    ymin = 0
    ymax = 0.12

    axis_opts = {
        "xmin": xmin,
        "xmax": xmax,
        "ymin": ymin,
        "ymax": ymax,
        "scale only axis": None,
        "xlabel": "Memorization score",
        "ylabel": "Density",
        "width": "6cm",
        "height": "8cm",
        "ytick": "{0, 0.05, 0.10}",
        "yticklabels": "{0, 0.05, 0.10}",
        "xlabel style": {"font": fontsize},
        "ylabel style": {"font": fontsize},
        "xticklabel style": {"font": fontsize},
        "yticklabel style": {"font": fontsize},
        "legend pos": "north east",
        "legend style": {"font": fontsize},
        "legend cell align": "left",
    }

    tex.append(f"\\begin{{axis}}[{dict2tex(axis_opts)}]")

    thickness = "ultra thick"

    hist_plot_opts = {
        "forget plot": None,
        "draw": "none",
        "fill opacity": 0.3,
        "hist": {
            "bins": bins,
            "density": "true",
            "data min": xmin,
            "data max": xmax,
        },
    }
    line_plot_opts = {
        "domain": f"{xmin}:{xmax}",
        "samples": 201,
        "mark": "none",
        "solid": None,
        thickness: None,
    }
    vert_plot_opts = {
        "forget plot": None,
        "densely dashed": None,
        thickness: None,
    }

    ms = [m3, m4]
    labels = ["$\\eta = 10^{-3}$", "$\\eta = 10^{-4}$"]
    colors = ["MyBlue", "MyYellow"]

    for m, label, color in zip(ms, labels, colors):
        hist_plot_opts["fill"] = color
        line_plot_opts["draw"] = color
        vert_plot_opts["draw"] = color

        tex.append(
            f"\\addplot [{dict2tex(hist_plot_opts)}] table[y index=0] {{%"
        )
        tex.append("data")
        for v in m:
            tex.append(str(v.item()))
        tex.append("};")

        f = Fitter(
            m, distributions=["johnsonsu"], xmin=xmin, xmax=xmax, bins=bins
        )
        f.fit()
        params = f.get_best()
        a, b, loc, scale = params["johnsonsu"]

        tex.append(f"\\addplot [{dict2tex(line_plot_opts)}] {{%")
        tex.append(f"johnsonsu(x, {a}, {b}, {loc}, {scale})")
        tex.append("};")
        tex.append(f"\\addlegendentry{{{label}}}")

        tex.append(f"\\addplot [{dict2tex(vert_plot_opts)}] coordinates {{%")
        tex.append(
            f"({np.quantile(m, 0.95)}, {ymin}) ({np.quantile(m, 0.95)}, {ymax})"
        )
        tex.append("};")

    tex.append("\\end{axis}")
    tex.append("\\end{tikzpicture}")
    tex.append("\\end{document}")
    return tex

def main():
    args = parse_args()

    bmnist3 = np.load(args.lr3_file)
    bmnist4 = np.load(args.lr4_file)

    M3 = bmnist3["M"]
    M4 = bmnist4["M"]

    tex = make_tex(M3, M4)
    with open(args.output, "w") as fp:
        fp.write("\n".join(tex))

if __name__ == "__main__":
    main()

max_stars_repo_name: alan-turing-institute/memorization

max_issues_repo_name: alan-turing-institute/memorization

max_forks_repo_name: alan-turing-institute/memorization


text:
#!/usr/bin/env python
from __future__ import division

import matplotlib as mpl
mpl.use('Agg')

import matplotlib.pyplot as plt
plt.style.use('seaborn-paper')
#plt.style.use('fivethirtyeight')
import numpy as np

mpl.rcParams.update({'text.usetex': True,
                     'lines.linewidth': 2.5,
                     'font.size': 16,
                     'xtick.labelsize': 'x-small',
                     'ytick.labelsize': 'x-small',
                     'axes.labelsize': 'x-small',
                     'axes.titlesize': 'large',
                     'axes.grid': True,
                     'grid.alpha': 0.73,
                     'lines.markersize': 12,
                     'legend.borderpad': 0.2,
                     'legend.fancybox': True,
                     'legend.fontsize': 13,
                     'legend.framealpha': 0.7,
                     'legend.handletextpad': 0.1,
                     'legend.labelspacing': 0.2,
                     'legend.loc': 'best',
                     'savefig.dpi': 240,
                     'pdf.compression': 9})

lambduh = 532e-9
L = np.linspace(-lambduh/2, lambduh/2, 3300)

k = 2*np.pi / lambduh

r1 = np.sqrt(0.99)
t1 = np.sqrt(1 - r1**2)
r2 = np.sqrt(0.9)
t2 = np.sqrt(1 - r2**2)

tcav1 = t1 * t1 / (1 - r1*r1*np.exp(-1j*2*k*L) )
tcav2 = t2 * t2 / (1 - r2*r2*np.exp(-1j*2*k*L) )

Tcav1 = np.abs(tcav1)**2
Tcav2 = np.abs(tcav2)**2

fig = plt.figure(772)

plt.semilogy(L*1e9, Tcav1, c = 'xkcd:Royal Blue',
    alpha = 0.7, rasterized = True, label = r'$R = 0.99$')
plt.semilogy(L*1e9, Tcav2, c = 'xkcd:Tomato',
    alpha = 0.7, rasterized = True, label = r'$R = 0.90$')

#ax[0].text(5, 2e-1, r'$|Z(\omega)|$')
#ax[0].text(1.1e6, 48e3, r'$L = 250~\mathrm{nH}$')
#plt.text(8, 0.5, r'$n_{air}   = 1$')
#plt.text(8, 0.6, r'$n_{glass} = 1.45$')
#plt.axvline(t[m]*180/np.pi, color='xkcd:Green', lw=4, alpha=0.5)
#plt.text(t[m]*180/np.pi, .9, "Brewster's angle", fontsize="x-small")
#plt.text(t[m]*180/np.pi, .8, "  = " + str(round(t[m]*180/np.pi,2)) + " deg", fontsize="x-small")

plt.xlabel('Length [nm]')
plt.ylabel(r'Transmission')
#plt.title('Fresnel Reflection (Air to Glass)')
plt.grid(True)
#plt.xlim([0,90])
#plt.ylim([0,1])
#ax[1].yaxis.set_ticks(np.arange(-180, 181, 45))
plt.legend()
plt.savefig("FabryPerot.pdf", bbox_inches='tight')


text:
# Copyright (c) 2021 Alix Routhier-Lalonde. Licence included in root of package.

import numpy as np
from flask import Flask, abort
from src.hyperparameters import Train
from src.utils import get_base_path
from src.data import Mode, Dataset
import pandas_market_calendars as mcal
from flask_cors import CORS
from hashlib import sha256
import json
import datetime
import os
import torch
import sys
from src import predict_stock

class Backend:
    def __init__(self):
        # creating the models
        self.model_daily = load_model(1)
        self.model_weekly = load_model(5)
        self.model_monthly = load_model(22)

        self.app = Flask(__name__)
        CORS(self.app)

        @self.app.route("/hello")
        def hello():
            """
            tester route
            :return: 'Hello World!'
            """
            return "Hello World!"

        @self.app.route("/predict/<code>/<int:num_days>", methods=['GET'])
        def predict(code, num_days):
            """
            Predicts stock prices for the next `num_days`
            :param code: the stock code
            :param num_days: the number of days to predict
            :return: the prediction data
            """
            data = predict_stock(code, num_days)
            if data is None:
                abort(404)
            data = data.squeeze().tolist()
            returned_data = []
            if not isinstance(data, list):
                data = [data]

            # getting dates
            nyse = mcal.get_calendar('NYSE')
            num_days_dates = len(data) * 2
            if num_days_dates < 6:
                num_days_dates = 6
            current_date = datetime.date.today() + datetime.timedelta(days=1)
            max_date = current_date + datetime.timedelta(days=num_days_dates)
            open_days = nyse.valid_days(start_date=str(
                current_date), end_date=str(max_date))
            if num_days == 1:
                returned_data = [{
                    "date": str(open_days[0]).split()[0],
                    "price": data
                }]
            else:
                for index, d in enumerate(data[-1]):
                    returned_data.append({
                        "date": str(open_days[index]).split()[0],
                        "price": d
                    })
            return json.dumps(returned_data)

        @self.app.route("/portfolios", methods=['GET'])
        def get_portfolios():
            """
            Will return all portfolio data in the backend
            If no portfolios data exists, will return a 404 code
            :return: the portfolio data
            """
            file = os.path.join(
                get_base_path(), "backend/data/portfolios.json")
            if os.path.exists(file):
                with open(file) as json_file:
                    data = json.load(json_file)
                    return json.dumps(data)
            abort(404)

        @self.app.route("/portfolios/<portfolio_id>", methods=['GET'])
        def get_portfolio(portfolio_id):
            """
            Given a portfolio id, it will return the portfolio data associated with that id
            If that portfolio cannot be found, it will return a 404 code
            :param portfolio_id: portfolio's id
            :return: the portfolio data
            """
            file = os.path.join(
                get_base_path(), "backend/data/portfolios.json")
            if os.path.exists(file):
                with open(file) as json_file:
                    data = json.load(json_file)
                    portfolios = data['portfolios']
                    if portfolio_id in portfolios:
                        return json.dumps(portfolios[portfolio_id])
                    else:
                        abort(404)
            return abort(404)

        @self.app.route("/portfolios/<name>", methods=['POST'])
        def add_portfolio(name):
            """
            Will add a portfolio to all the backend's portfolios.
            If no portfolio data exists in the backend, it will initialize it and create a file for it
            If the portfolio already exists, it will return all the current portfolios with a 208 code.
            :param name: the portfolio's name
            :return: all backend portfolios
            """
            file = os.path.join(
                get_base_path(), "backend/data/portfolios.json")
            portfolio_id = str(sha256(name.encode('utf-8')).hexdigest())
            if os.path.exists(file):
                with open(file) as json_file:
                    data = json.load(json_file)
                    portfolios = data['portfolios']
                    if portfolio_id in portfolios:
                        return json.dumps(data), 208
                    else:
                        new_portfolio = {
                            "id": portfolio_id,
                            "name": name,
                            "stocks": []
                        }
                        portfolios[portfolio_id] = new_portfolio
            else:
                data = {
                    "portfolios": {
                        portfolio_id: {
                            "id": portfolio_id,
                            "name": name,
                            "stocks": []
                        }
                    },
                    "stocks": {}
                }

            with open(file, 'w') as outfile:
                json.dump(data, outfile)
            return json.dumps(data)

        @self.app.route("/portfolios/<portfolio_id>", methods=['DELETE'])
        def delete_portfolio(portfolio_id):
            """
            Given a portfolio id, it delete the portfolio attached to it
            If that portfolio cannot be found, it will return a 404 code
            :param portfolio_id: portfolio's id
            :return: the portfolio data
            """
            file = os.path.join(
                get_base_path(), "backend/data/portfolios.json")
            if os.path.exists(file):
                with open(file) as json_file:
                    data = json.load(json_file)
                    portfolios = data['portfolios']
                    stocks = data['stocks']
                    if portfolio_id in portfolios:
                        portfolio = portfolios[portfolio_id]
                        stock_ids = portfolio['stocks']
                        del portfolios[portfolio_id]
                        for stock_id in stock_ids:
                            condition = False
                            for portfolio in portfolios.values():
                                if stock_id in portfolio['stocks']:
                                    condition = True
                            if not condition and stock_id in stocks:
                                del stocks[stock_id]
                        with open(file, 'w') as outfile:
                            json.dump(data, outfile)
                        return json.dumps(data)
                    else:
                        abort(404)
            return abort(404)

        @self.app.route("/portfolios/<portfolio_id>/stocks/<name>/<int:mode>", methods=['POST'])
        def add_stock_to_portfolio(portfolio_id, name, mode):
            """
            Will add a stock to a portfolio.
            If the stock already exists in the portfolio, it will return all the current portfolios with a 208 code.
            If the portfolio cannot be found, a 404 code will be returned
            :param portfolio_id: the portfolio's id
            :param name: the stock's code
            :param mode: the number of predicted days wanted
            :return: the whole portfolios data
            """
            file = os.path.join(
                get_base_path(), "backend/data/portfolios.json")
            if os.path.exists(file):
                with open(file) as json_file:
                    data = json.load(json_file)
                    portfolios = data['portfolios']
                    stocks = data['stocks']
                    if portfolio_id in portfolios:
                        stock_id = str(
                            sha256((name + str(mode)).encode('utf-8')).hexdigest())
                        if stock_id in portfolios[portfolio_id]["stocks"] and stock_id in stocks:
                            return json.dumps(data), 208
                        else:
                            new_stock = {
                                "id": stock_id,
                                "name": name,
                                "mode": mode
                            }
                            if stock_id not in stocks:
                                stocks[stock_id] = new_stock
                            if stock_id not in portfolios[portfolio_id]["stocks"]:
                                portfolios[portfolio_id]["stocks"].append(
                                    stock_id)
                            with open(file, 'w') as outfile:
                                json.dump(data, outfile)
                            data['new_stock'] = new_stock
                            return json.dumps(data)
                    else:
                        abort(404)
            else:
                abort(404)

        @self.app.route("/portfolios/<portfolio_id>/stocks/<stock_id>", methods=['DELETE'])
        def delete_stock_from_portfolio(portfolio_id, stock_id):
            """
            Will delete a stock from a portfolio.
            If that stock doesn't exist in that portfolio, will return a 404 error
            :param portfolio_id: the portfolio's id
            :param stock_id: the stock's id
            :return: the whole portfolios data
            """
            file = os.path.join(
                get_base_path(), "backend/data/portfolios.json")
            if os.path.exists(file):
                with open(file) as json_file:
                    data = json.load(json_file)
                    portfolios = data['portfolios']
                    stocks = data['stocks']
                    if portfolio_id in portfolios:
                        if stock_id in portfolios[portfolio_id]["stocks"]:
                            portfolios[portfolio_id]["stocks"].remove(stock_id)

                            # purging in the stocks
                            condition = False
                            for portfolio in portfolios.values():
                                if stock_id in portfolio['stocks']:
                                    condition = True
                            if not condition and stock_id in stocks:
                                del stocks[stock_id]
                            with open(file, 'w') as outfile:
                                json.dump(data, outfile)
                            return json.dumps(data)
                        else:
                            abort(404)
                    else:
                        abort(404)
            else:
                abort(404)

    def run(self):
        self.app.run(host='0.0.0.0', port=8000, threaded=False)

    def predict_data(self, code: str, num_days):
        """
        Predicts stock prices for the next num_days days
        :param code: the stock's code
        :param num_days: the number of days to predict for
        :return: None if the model cannot be found, otherwise predicted stock data
        """
        # get the data and model
        dataset = Dataset(code, mode=num_days, y_flag=True)
        dataset.transform_to_torch()
        if num_days == 1:
            model = self.model_daily
        elif num_days == 5:
            model = self.model_weekly
        elif num_days == 22:
            model = self.model_monthly
        else:
            return None
        if model is None:
            return None

        # predict
        data = torch.from_numpy(np.array(dataset.prediction_data)).float()
        if torch.cuda.is_available():
            data = data.to(device='cuda')
        with torch.no_grad():
            predicted = model(data)

        # re-transforming to numpy
        predicted = predicted.detach().cpu().numpy()
        return predicted
        # return dataset.normalizer.inverse_transform(predicted)

def load_model(num_days: int):
    """
    Will load a model from the number of days it has been trained for
    :param num_days: the number of days the model has been trained for
    :return: the model
    """
    model = Train.model(6, Train.hidden_dim, Train.num_dim,
                        Train.dropout, num_days)

    # if we have a gpu available, we store the models there
    gpu = torch.cuda.is_available()
    if gpu:
        model.to('cuda')
    file = os.path.abspath(os.path.join(
        get_base_path(), f'src/model/models/model-{num_days}.hdf5'))
    print(file)
    if not os.path.exists(file):
        return None
    model.load_state_dict(torch.load(file))
    model.eval()
    return model

if __name__ == '__main__':
    app = Backend()
    app.run()

max_stars_repo_name: SgtMilk/stock-prediction

max_issues_repo_name: SgtMilk/stock-prediction

max_forks_repo_name: SgtMilk/stock-prediction


text:
import random as rnd
import json
from sympy import pretty, Symbol
import math

def elementosListaEhDistinta(lista):
    for indiceLista in range(len(lista)):
        for indiceListaComparacao in range(len(lista)):
            if indiceLista == 4:
                return True
            elif lista[indiceLista] == lista[indiceListaComparacao] and indiceLista != indiceListaComparacao:
                return False

enunciado = [None]*200
k = 0
while k < 200:
    questoes = open("./src/json/Gerador3/questao{}.json".format(k+1), 'w')

    expoenteConcentracao = rnd.randint(1,1000)
    resposta = round((expoenteConcentracao ** 2)/3, 2)

    e = Symbol('e')

    listLetra = ["A","B","C","D","E"]
    questaoCerta = rnd.choice(listLetra)
    questaoInvertida = rnd.choice(listLetra)

    # Verifica se a letra com a resposta certa e invertida são iguais
    repetida = True
    while repetida == True:
        if questaoCerta == questaoInvertida:
            questaoCerta = rnd.choice(listLetra)
            questaoInvertida = rnd.choice(listLetra)
        else: repetida = False

    listAlternativas = [0,0,0,0,0]
    isCorrect = ['','','','','']
    howGenerated = ['','','','','']

    # Insere a resposta certa na letra escolhida para ser certa, uma letra recebe a questão invertida e o resto recebe números aleatórios
    while(elementosListaEhDistinta(listAlternativas) == False):
        possuiQuestaoInvertida = rnd.randint(0,1)

        for numLetra in range(0, 5):
            if questaoCerta == listLetra[numLetra]:
                listAlternativas[numLetra] = resposta
                isCorrect[numLetra] = "Sim"
                howGenerated[numLetra] = "nenhum"
            elif questaoInvertida == listLetra[numLetra] and possuiQuestaoInvertida == 1:
                ehNegativoPositivo = rnd.randint(0,1)

                if ehNegativoPositivo == 0:

                    listAlternativas[numLetra] = str(round((expoenteConcentracao ** 3)/2, 2)).replace(".",",")
                    isCorrect[numLetra] = "Nao"
                    howGenerated[numLetra] = "invertida e positiva"
                else:

                    listAlternativas[numLetra] = str(round(-(expoenteConcentracao ** 3)/2, 2)).replace(".",",")
                    isCorrect[numLetra] = "Nao"
                    howGenerated[numLetra] = "invertida e negativa"
            else:
                ehNegativoPositivo = rnd.randint(0,1)

                if (ehNegativoPositivo == 0):
                    listAlternativas[numLetra] = str(round((rnd.randint(1,expoenteConcentracao) ** 2)/3, 2)).replace(".",",")
                    isCorrect[numLetra] = "Nao"
                    howGenerated[numLetra] = "gerada aleatoriamente e positiva"
                else:
                    listAlternativas[numLetra] = str(round(-(rnd.randint(1,expoenteConcentracao) ** 2)/3, 2)).replace(".",",")
                    isCorrect[numLetra] = "Nao"
                    howGenerated[numLetra] = "gerada aleatoriamente e negativa"

    # Cria a variável que será convertida em um arquivo json
    dados = {
        'logaritmo' : [
            {
                'base' : 10,
                'logaritmando' : 1 / 10 ** expoenteConcentracao,
                'expoente' : expoenteConcentracao,
                'resultado' : resposta
            }
        ],
        'respostas' : [
            {
                'letra': 'a)',
                'resposta': listAlternativas[0],
                'correta': isCorrect[0],
                'tipoerro': howGenerated[0]
            },
            {
                'letra': 'b)',
                'resposta': listAlternativas[1],
                'correta': isCorrect[1],
                'tipoerro': howGenerated[1]
            },
            {
                'letra': 'c)',
                'resposta': listAlternativas[2],
                'correta': isCorrect[2],
                'tipoerro': howGenerated[2]
            },{
                'letra': 'd)',
                'resposta': listAlternativas[3],
                'correta': isCorrect[3],
                'tipoerro': howGenerated[3]
            },{
                'letra': 'e)',
                'resposta': listAlternativas[4],
                'correta': isCorrect[4],
                'tipoerro': howGenerated[4]
            }
        ],
        'atributosquestao': [
            {
                'enunciado': 'O pH de uma solução é definido por pH = log(1⁄H+), em que H+ é a concentração de hidrogênio em íons-gama por litro de solução. Determine o quadrado sobre três do pH de uma solução tal que H+ = 1,0 × 10{}.'.format(expoenteConcentracao),
                'corretaspossiveis': listAlternativas[isCorrect.index("Sim")],
                'corretas': isCorrect.count("Sim"),
                'aleatoriapositiva': howGenerated.count("gerada aleatoriamente e positiva"),
                'aleatorianegativa': howGenerated.count("gerada aleatoriamente e negativa"),
                'invertidapositiva': howGenerated.count("invertida e positiva"),
                'invertidanegativa': howGenerated.count("invertida e negativa"),
                'respostascorretas': listLetra[isCorrect.index("Sim")]
            }
        ]
    }

    # Verifica os enunciados
    if dados['atributosquestao'][0]['enunciado'] in enunciado:
        continue
    else:

    # Armazena os enunciados
        enunciado[k] = dados['atributosquestao'][0]['enunciado']

    # Cria o arquivo JSON
        print("\nquestao {}\n".format(k+1),json.dumps(dados))
        json.dump(dados, questoes, indent=4)
        k = k + 1

questoes.close()

max_stars_repo_name: oismaelazevedo/PIBIC-EM

max_issues_repo_name: oismaelazevedo/PIBIC-EM

max_forks_repo_name: oismaelazevedo/PIBIC-EM

reason: from sympy