text:
import seaborn as sns
import pandas as pd
import cdlib
import matplotlib.pyplot as plt
from typing import Callable

__all__ = [
    "plot_com_properties_relation",
    "plot_com_stat",
    "plot_sim_matrix",
    "plot_scoring",
]

def plot_sim_matrix(
    clusterings: list, scoring: Callable[[object, object], object]
) -> object:
    """
    Plot a similarity matrix between a list of clusterings, using the provided scoring function.

    :param clusterings: list of clusterings to compare
    :param scoring: the scoring function to use
    :return: a ClusterGrid instance

    Example:

    >>> from cdlib import algorithms, viz, evaluation
    >>> import networkx as nx
    >>> g = nx.karate_club_graph()
    >>> coms = algorithms.louvain(g)
    >>> coms2 = algorithms.walktrap(g)
    >>> clustermap = viz.plot_sim_matrix([coms,coms2],evaluation.adjusted_mutual_information)

    """
    forDF = []
    for c in clusterings:
        cID = c.get_description()
        for c2 in clusterings:
            c2ID = c2.get_description()
            forDF.append([cID, c2ID, scoring(c, c2).score])
    df = pd.DataFrame(columns=["com1", "com2", "score"], data=forDF)
    df = df.pivot("com1", "com2", "score")
    return sns.clustermap(df)

def plot_com_stat(
    com_clusters: list, com_fitness: Callable[[object, object, bool], object]
) -> object:
    """
    Plot the distribution of a property among all communities for a clustering, or a list of clusterings (violin-plots)

    :param com_clusters: list of clusterings to compare, or a single clustering
    :param com_fitness: the fitness/community property to use
    :return: the violin-plots

    Example:

    >>> from cdlib import algorithms, viz, evaluation
    >>> import networkx as nx
    >>> g = nx.karate_club_graph()
    >>> coms = algorithms.louvain(g)
    >>> coms2 = algorithms.walktrap(g)
    >>> violinplot = viz.plot_com_stat([coms,coms2],evaluation.size)

    """
    if isinstance(com_clusters, cdlib.classes.clustering.Clustering):
        com_clusters = [com_clusters]

    allVals = []
    allNames = []
    for c in com_clusters:
        prop = com_fitness(c.graph, c, summary=False)
        allVals += prop
        allNames += [c.get_description()] * len(prop)

    ax = sns.violinplot(allNames, allVals, cut=0, saturation=0.5, palette="Set3")
    for tick in ax.get_xticklabels():
        tick.set_rotation(90)

    plt.ylabel("%s" % com_fitness.__name__)
    plt.xlabel("Algorithm")
    plt.tight_layout()

    return ax

def plot_com_properties_relation(
    com_clusters: object,
    com_fitness_x: Callable[[object, object, bool], object],
    com_fitness_y: Callable[[object, object, bool], object],
    **kwargs: dict
) -> object:
    """
    Plot the relation between two properties/fitness function of a clustering

    :param com_clusters:  clustering(s) to analyze (cluster or cluster list)
    :param com_fitness_x: first fitness/community property
    :param com_fitness_y: first fitness/community property
    :param kwargs: parameters for the seaborn lmplot
    :return: a seaborn lmplot

    Example:

    >>> from cdlib import algorithms, viz, evaluation
    >>> import networkx as nx
    >>> g = nx.karate_club_graph()
    >>> coms = algorithms.louvain(g)
    >>> coms2 = algorithms.walktrap(g)
    >>> lmplot = viz.plot_com_properties_relation([coms,coms2],evaluation.size,evaluation.internal_edge_density)
    """
    if isinstance(com_clusters, cdlib.classes.clustering.Clustering):
        com_clusters = [com_clusters]

    for_df = []

    for c in com_clusters:
        x = com_fitness_x(c.graph, c, summary=False)
        y = com_fitness_y(c.graph, c, summary=False)
        for i, vx in enumerate(x):
            for_df.append([c.get_description(), vx, y[i]])

    df = pd.DataFrame(
        columns=["Method", com_fitness_x.__name__, com_fitness_y.__name__], data=for_df
    )
    ax = sns.lmplot(
        x=com_fitness_x.__name__,
        y=com_fitness_y.__name__,
        data=df,
        hue="Method",
        fit_reg=False,
        legend=False,
        x_bins=100,
        **kwargs
    )
    plt.legend(loc="best")
    plt.tight_layout()

    return ax

def plot_scoring(
    graphs: list,
    ref_partitions: object,
    graph_names: list,
    methods: list,
    scoring: Callable[
        [object, object], object
    ] = cdlib.evaluation.adjusted_mutual_information,
    nbRuns: int = 5,
) -> object:
    """
    Plot the scores obtained by a list of methods on a list of graphs.

    :param graphs: list of graphs on which to make computations
    :param ref_partitions: list of reference clusterings corresponding to graphs
    :param graph_names: list of the names of the graphs to display
    :param methods: list of functions that take a graph as input and return a Clustering as output
    :param scoring: the scoring function to use, default anmi
    :param nbRuns: number of runs to do for each method on each graph
    :return: a seaborn lineplot

    Example:

    >>> from cdlib import algorithms, viz, evaluation
    >>> import networkx as nx
    >>> g1 = nx.algorithms.community.LFR_benchmark_graph(1000, 3, 1.5, 0.5, min_community=20, average_degree=5)
    >>> g2 = nx.algorithms.community.LFR_benchmark_graph(1000, 3, 1.5, 0.7, min_community=20, average_degree=5)
    >>> names = ["g1", "g2"]
    >>> graphs = [g1, g2]
    >>> for g in graphs:
    >>>     references.append(NodeClustering(communities={frozenset(g.nodes[v]['community']) for v in g}, graph=g, method_name="reference"))
    >>> algos = [algorithms.crisp_partition.louvain, algorithms.crisp_partition.label_propagation]
    >>> viz.plot_scoring(graphs, references, names, algos, nbRuns=2)

    """
    forDF = []
    for i, g in enumerate(graphs):
        for m in methods:
            for r in range(nbRuns):
                partition = m(g)

                score = scoring(partition, ref_partitions[i]).score
                forDF.append([graph_names[i], score, partition.get_description()])
    df = pd.DataFrame(columns=["graph", "score", "method"], data=forDF)
    ax = sns.lineplot(x="graph", y="score", hue="method", data=df, legend="brief")
    ax.legend(loc="best")
    for tick in ax.get_xticklabels():
        tick.set_rotation(90)
    plt.tight_layout()

    return ax

max_stars_repo_name: xing-lab-pitt/cdlib

max_issues_repo_name: xing-lab-pitt/cdlib

max_forks_repo_name: xing-lab-pitt/cdlib

reason: import networkx


text:
"""Validate the functionality of the `gspits.one_dim.system` module."""

import numpy as np
import pytest
from attr import dataclass
from hypothesis import given
from hypothesis import strategies as stg
from numba import njit
from numpy import pi

from gspits import Partition as Mesh
from gspits.one_dim import BlochState, ExternalPotential, Hamiltonian

@njit
def _ho_potential(domain_mesh: np.ndarray) -> np.ndarray:
    """Evaluate the harmonic potential in the mesh."""
    return 0.5 * domain_mesh ** 2

@dataclass(frozen=True, slots=True)
class _TestHamiltonian:
    """Define a testing Hamiltonian."""

    __slots__ = ()

    @property
    def interaction_factor(self):
        """Bose gas interaction factor."""
        return 4 * pi * 1e-2

    @property
    def external_potential(self):
        """External potential function."""
        return _ho_potential

def test_hamiltonian():
    """Verify subtyping relationships for Hamiltonian instances."""
    hamiltonian = _TestHamiltonian()
    assert isinstance(hamiltonian.external_potential, ExternalPotential)
    assert isinstance(hamiltonian, Hamiltonian)

@given(
    lower_bound=stg.floats(min_value=-2, max_value=-1, allow_nan=False),
    upper_bound=stg.floats(min_value=1, max_value=2, allow_nan=False),
    num_segments=stg.integers(min_value=128, max_value=512),
    wave_vector=stg.floats(
        min_value=-3 * pi, max_value=3 * pi, allow_nan=False
    ),
)
def test_bloch_state(
    lower_bound: float,
    upper_bound: float,
    num_segments: int,
    wave_vector: float,
):
    """Check the basic functionality of `BlochState` instances."""
    mesh = Mesh(
        lower_bound=lower_bound,
        upper_bound=upper_bound,
        num_segments=num_segments,
    )
    bloch_state = BlochState.plane_wave(mesh, wave_vector)
    # A plane wave is normalized, and its norm must be approximately one.
    assert bloch_state.norm == pytest.approx(1, abs=1e-8)
    assert bloch_state.norm == pytest.approx(
        bloch_state.periodic_component.norm, abs=1e-8
    )

ext: py

max_stars_repo_name: oarodriguez/gspits

max_issues_repo_name: oarodriguez/gspits

max_forks_repo_name: oarodriguez/gspits

reason: import numpy,from numpy,from numba


text:
import os
import sys
project_dir = os.path.join(os.path.dirname(__file__), os.pardir, os.pardir)
sys.path.append(project_dir)

import json
import numpy as np
import pandas as pd
from pandas.io.json import json_normalize
import hashlib
import time
import boto3

from sqlalchemy import create_engine, exc, MetaData, select
from sqlalchemy.engine.url import URL
from sqlalchemy.sql import or_, and_
from sqlalchemy.types import JSON as typeJSON

def connect_database(database_dict):

    DATABASE = database_dict

    db_url = URL(**DATABASE)
    engine = create_engine(db_url)
    meta = MetaData()
    meta.bind = engine

    return meta

def tab_raw_data(s3key, s3object):

    def build_raw_df(rec_list):
        df_list = []
        for rec in rec_list:
            rec["rec_string"] = json.dumps(rec, sort_keys=True)
            rec = json_normalize(rec)
            df_list.append(rec)
        raw_data = pd.concat(df_list)
        return raw_data

    #read data file and convert it to a list of dicts.
    file = s3object['Body'].read()
    rec_list = json.loads(file)
    if type(rec_list) is dict:
        rec_list = [rec_list]
    if type(rec_list) is not list:
        raise Exception("It should be a list.")

    raw_data = build_raw_df(rec_list)

    #Get rid of oid (pymongo) if it exists. From now on, object is the same regardless of source.
    raw_data.drop("_id.$oid", axis=1, inplace=True, errors='ignore')

    #clean date fields
    raw_data['startTime'] = pd.to_datetime(raw_data['startTime'].str[:-4])
    raw_data['endTime'] = pd.to_datetime(raw_data['endTime'].str[:-4])
    raw_data['startTime'] = raw_data['startTime'].astype(pd.Timestamp)
    raw_data['endTime'] = raw_data['endTime'].astype(pd.Timestamp)

    #get creation time
    raw_data['date_created'] = pd.to_datetime('today')

    #during creation time, date_updated equals date_created
    raw_data['date_updated'] = raw_data['date_created']

    #get_file_name
    raw_data['file_name'] = s3key

    #get DataFile hash
    raw_data['json_hash'] = raw_data["rec_string"].apply(lambda x: hashlib.sha1(x.encode()).hexdigest())

    return raw_data

def sep_aji_records(raw_data, aji_type):
    recs = []
    for aji in raw_data[aji_type].iloc[0]:
        aji["rootStartTime"] = raw_data.startTimeMillis.iloc[0]
        df_aji = json_normalize(aji)
        df_aji["rec"] = json.dumps(aji, sort_keys=True)
        recs.append(df_aji)
    df = pd.concat(recs)
    return df

def hash_raw_aji(raw_aji):
    aji_hash = hashlib.sha1(raw_aji.encode()).hexdigest()
    return aji_hash

def align_all_columns(df, col_list):
    all_columns = pd.DataFrame(np.nan, index=[0], columns=col_list)
    df,_ = df.align(all_columns, axis=1)
    return df

def tab_jams(raw_data):
    if ('jams' not in raw_data) or (raw_data.jams.iloc[0] is np.nan):
        print("No jams in this data file.")
        return

    col_dict = {
                'blockingAlertUuid': "blocking_alert_id",
                'startNode': "start_node",
                 'endNode': "end_node",
                 'pubMillis': "pub_millis",
                 'roadType': "road_type",
                 'speedKMH': "speed_kmh",
                 'turnType': "turn_type",
                 }

    other_cols = ['id','city', 'country','delay', 'length',
                  'uuid', 'street', 'level', 'line', 'pub_utc_date']
    col_list = list(col_dict.values())
    col_list = col_list + other_cols

    df_jams = sep_aji_records(raw_data, "jams")
    df_jams = (df_jams
               .rename(columns=col_dict)
               .pipe(align_all_columns, col_list=col_list)
               .assign(pub_utc_date=lambda x: pd.to_datetime(x["pub_millis"], unit='ms'),
                       id=df_jams["rec"].apply(hash_raw_aji)
                      )
              )
    df_jams = df_jams[col_list]

    return df_jams

def tab_irregularities(raw_data):
    if ('irregularities' not in raw_data) or (raw_data.irregularities.iloc[0] is np.nan):
        print("No irregularities in this data file.")
        return

    col_dict = {
                'detectionDateMillis': "detection_date_millis",
                'detectionDate': "detection_date",
                'updateDateMillis': "update_date_millis",
                'updateDate': "update_date",
                'regularSpeed': "regular_speed",
                'delaySeconds': "delay_seconds",
                'jamLevel': "jam_level",
                'driversCount': "drivers_count",
                'alertsCount': "alerts_count",
                'nThumbsUp': "n_thumbs_up",
                'nComments': "n_comments",
                'nImages': "n_images",
                'id': "uuid",
                 }

    other_cols = ['id', 'street', 'city', 'country', 'speed',
                  'seconds', 'length', 'trend', 'type', 'severity', 'line']
    col_list = list(col_dict.values())
    col_list = col_list + other_cols

    df_irregs = sep_aji_records(raw_data, "irregularities")
    df_irregs = (df_irregs
                 .rename(columns=col_dict)
                 .pipe(align_all_columns, col_list=col_list)
                 .assign(detection_utc_date=lambda x: pd.to_datetime(x["detection_date_millis"], unit='ms'),
                         update_utc_date=lambda x: pd.to_datetime(x["update_date_millis"], unit='ms'),
                         id=df_irregs["rec"].apply(hash_raw_aji)
                        )
              )
    df_irregs = df_irregs[col_list]

    return df_irregs

def tab_alerts(raw_data):
    if ('alerts' not in raw_data) or (raw_data.alerts.iloc[0] is np.nan):
        print("No alerts in this data file.")
        return

    col_dict = {
                'pubMillis': "pub_millis",
                'roadType': "road_type",
                'reportDescription': "report_description",
                'reportRating': "report_rating",
                'nThumbsUp': "thumbs_up",
                'jamUuid': "jam_uuid",
                'reportByMunicipalityUser': 'report_by_municipality_user',
                 }

    other_cols = ['id', 'uuid', 'street', 'city', 'country', 'location', 'magvar',
                  'reliability', 'type', 'subtype' ]
    col_list = list(col_dict.values())
    col_list = col_list + other_cols

    df_alerts = sep_aji_records(raw_data, "alerts")
    df_alerts = (df_alerts
                 .rename(columns=col_dict)
                 .pipe(align_all_columns, col_list=col_list)
                 .assign(location=df_alerts.apply(lambda row: {'x': row["location.x"], 'y': row["location.y"]},
                                                  axis=1),
                         pub_utc_date=lambda x: pd.to_datetime(x["pub_millis"], unit='ms'),
                         id=df_alerts["rec"].apply(hash_raw_aji)
                        )

                )
    df_alerts = df_alerts[col_list]

    return df_alerts

#Connection and initial setup
DATABASE = {
'drivername': "postgresql",
'host': "localhost",
'port': 5432,
'username': "tester",
'password': "testmobility123",
'database': "test_mobility",
}

meta = connect_database(DATABASE)

s3 = boto3.client('s3')
bucket="scripted-waze-data-929310922828-test"
#Iterate over all s3 raw data objects
all_data_files = []
paginator = s3.get_paginator('list_objects')
page_iterator = paginator.paginate(Bucket=bucket)
for page in page_iterator:
    all_data_files += [c["Key"] for c in page["Contents"]]

for file in all_data_files:
    #Read raw file
    obj = s3.get_object(Bucket=bucket, Key=file)
    raw_data = tab_raw_data(file, obj)

    i=1
    n = len(raw_data)
    for _, row in raw_data.iterrows():
        start = time.time()
        row = row.to_frame().transpose()
        #Store data_file in database
        col_dict = {"startTimeMillis": "start_time_millis",
                    "endTimeMillis": "end_time_millis",
                    "startTime": "start_time",
                    "endTime": "end_time",
                    "date_created": "date_created" ,
                    "date_updated": "date_updated",
                    "file_name": "file_name",
                    "json_hash": "json_hash",
                    }

        raw_data_tosql = row.rename(columns=col_dict)
        raw_data_tosql = raw_data_tosql[list(col_dict.values())]
        try:
            raw_data_tosql.to_sql(name="data_files", schema="waze", con=meta.bind, if_exists="append", index=False)
        except exc.IntegrityError:
            print("Data file is already stored in the relational database. Stopping process.")
            break

        #Introspect data_file table
        meta.reflect(schema="waze")
        data_files = meta.tables["waze.data_files"]
        datafile_result = (select([data_files.c.id]).where(data_files.c.json_hash == raw_data["json_hash"].iloc[0])
                                                   .execute()
                                                   .fetchall()
                          )
        if len(datafile_result) > 1:
            raise Exception("There should be only one result")

        datafile_id = datafile_result[0][0]

        #Store jams in database
        jams_tosql = tab_jams(row)
        if jams_tosql is not None:
            jams_tosql["datafile_id"] = datafile_id
            jams_tosql.to_sql(name="jams", schema="waze", con=meta.bind, if_exists="append", index=False,
                                  dtype={"line": typeJSON}
                                 )

        #Store alerts in database
        alerts_tosql = tab_alerts(row)
        if alerts_tosql is not None:
            alerts_tosql["datafile_id"] = datafile_id
            alerts_tosql.to_sql(name="alerts", schema="waze", con=meta.bind, if_exists="append", index=False,
                                  dtype={"location": typeJSON}
                                 )

        #Store irregularities in databse
        irregs_tosql = tab_irregularities(row)
        if irregs_tosql is not None:
            irregs_tosql["datafile_id"] = datafile_id
            irregs_tosql.to_sql(name="irregularities", schema="waze", con=meta.bind, if_exists="append", index=False,
                                  dtype={"line": typeJSON}
                                 )

        end = time.time()
        elapsed = str(int(end-start))

        print("Stored DataFile", str(i), "of", str(n),"from", file, "in", elapsed, "seconds.")
        i += 1

max_stars_repo_name: joinvalle/Joinville-Smart-Mobility

max_issues_repo_name: pmj-sepud/Joinville-Smart-Mobility

max_forks_repo_name: pmj-sepud/Joinville-Smart-Mobility

reason: import numpy


text:
from collections import defaultdict
import itertools
import logging
from typing import List, Dict

import numpy as np
import torch
import torch_geometric

from common.torch_helpers import to_numpy

def print_loss_metrics(logger, step, losses=None, metrics=None, input_stats=None):

    logger.info('Step {} validation:'.format(step))
    if losses is not None:
        loss_str = ', '.join(['{:.4g} ({})'.format(losses[k], k) for k in losses if k != 'total'])
        loss_str = 'Losses: ' + loss_str + ' | {:.4g} (total)'.format(losses['total'])
        logger.info(loss_str)

    if input_stats is not None:
        meas_metrics_str = 'Input Measurements have errors: [rot] {:.1f} (mean), {:.1f} (med)'.format(
            input_stats['rot_error_mean'], input_stats['rot_error_med']
        )
        if 'trans_error_mean' in input_stats:
            meas_metrics_str += ' | [trans] errors: {:.2f} | {:.2f}'.format(
                input_stats['trans_error_mean'], input_stats['trans_error_med'])

        logger.info(meas_metrics_str)

    if metrics is not None:  # Training on graphs
        rot_metrics_str = 'Rot: ' + \
                          ', '.join(['{:.1f} ({})'.format(metrics[k] * 100, k[9:])
                                    for k in metrics if k.startswith('rot_ecdf')]) + \
                          ' | {:.3f} (mean), {:.3f} (median)'.format(
                              metrics['rot_error_mean'], metrics['rot_error_med'])
        logger.info(rot_metrics_str)

        if 'trans_error_mean' in metrics:
            trans_metrics_str = 'Trans: ' + \
                                ', '.join(['{:.1f} ({})'.format(metrics[k] * 100, k[11:])
                                          for k in metrics if k.startswith('trans_ecdf')]) + \
                                ' | {:.2f} (mean), {:.2f} (median)'.format(
                                    metrics['trans_error_mean'], metrics['trans_error_med'])
            logger.info(trans_metrics_str)

def compute_metrics(data_batch: torch_geometric.data.batch.Batch, pred_batch, is_rel, only_good):
    if is_rel:
        return compute_rel_metrics(data_batch, pred_batch, check_only_good=only_good)
    else:
        return compute_abs_metrics(data_batch, pred_batch)

def compute_abs_metrics(data_batch: torch_geometric.data.batch.Batch, pred_batch, align=True):
    metrics = {}
    assert data_batch.num_graphs == 1, 'This can only handle batch sizes of 1'

    # Compute statistics for measured transforms
    pairs = data_batch.edge_index
    i, j = pairs
    se3_gt_i = data_batch['t_star'][i]
    se3_gt_j = data_batch['t_star'][j]
    rel_gt = se3_gt_j * se3_gt_i.inv()
    if 'valid_gt' in data_batch:
        valid_pairs = data_batch.valid_gt[i] & data_batch.valid_gt[j]
    else:
        valid_pairs = torch.ones(pairs.shape[1], dtype=torch.bool, device=pairs.device)
    # Compute measurement errors
    measured_errs = rel_gt.compare(data_batch['tij_hat'])
    metrics['measured_rot_errors'] = measured_errs['rot_deg'][valid_pairs]
    if 'trans' in measured_errs:
        metrics['measured_trans_errors'] = measured_errs['trans'][valid_pairs]

    # Evaluate predictions
    T_cw_pred = pred_batch[-1]  # only evaluate final iteration
    T_cw_gt = data_batch.t_star
    if hasattr(data_batch, 'valid_gt'):
        # Remove nodes without valid groundtruth
        T_cw_pred = T_cw_pred[data_batch['valid_gt']]
        T_cw_gt = T_cw_gt[data_batch['valid_gt']]
    if align:
        T_cw_pred = align_predictions(T_cw_gt, T_cw_pred)

    pred_err = T_cw_gt.compare(T_cw_pred)
    for k in pred_err:
        metrics[f'{k}_err_mean'] = torch.mean(pred_err[k])
        metrics[f'{k}_err_med'] = torch.median(pred_err[k])

    return metrics

def align_predictions(rot1, rot2, method='median', mask=None):
    """Align two sets of rotations, similar to what is done in the rotation
    averaging papers

    Notes:
        1. Only works for SO3
        2. (Update 2 June 2021) The rotation conventions are the same as that
           used in NeuRoRA/Chatterjee code. Specifically they assume that the
           rotations maps from world to camera, i.e. p_c = rot1 * p_w.

    Args:
        rot1: Groundtruth rotations, mapping points from the world frame to
          the camera frame
        rot2: Predicted rotations
        method: Either 'mean' or 'median'
        mask: If provided, will ignore entries corresponding to mask values of 0

    Returns:
        Aligned rot2
    """

    Transform = type(rot1)

    # Chooses a random frame for initial alignment
    if mask is None:
        j = np.random.randint(rot1.shape[0])
        rot1a = rot1 * rot1[j:j+1].inv()
        rot2a = rot2 * rot2[j:j+1].inv()
    else:
        j = np.random.choice(to_numpy(torch.nonzero(mask)[:, 0]))
        rot1a = rot1 * rot1[j:j+1].inv()
        rot2a = rot2 * rot2[j:j+1].inv()

    d, count = float('inf'), 0
    while d > 1e-5 and count < 20:
        if mask is None:
            align_w = (rot2a.inv() * rot1a).log()
        else:
            align_w = (rot2a[mask].inv() * rot1a[mask]).log()
        if method == 'mean':
            align_w = torch.mean(align_w, dim=0, keepdim=True)
        elif method == 'median':
            align_w = torch.median(align_w, dim=0, keepdim=True).values
        else:
            raise AssertionError('Alignment method not recognized')
        d = torch.norm(align_w)
        correction = Transform.exp(align_w)
        rot2a = rot2a * correction
        count += 1

    aligned = rot2a * rot1[j:j+1]
    return aligned

def compute_rel_metrics(data_batch: torch_geometric.data.batch.Batch, pred_batch, check_only_good):
    """Accumulate statistics

    Args:
        data_batch: Validation/Test data batch
        pred_batch: network outputs
    """
    metrics = {}
    batch_size = data_batch.num_graphs
    predicted_transforms = pred_batch[-1]  # Only evaluate last step

    # Check if predicted rotations are normalize, if not normalize them
    assert predicted_transforms.is_valid(), 'Predicted transforms are not valid SO(3)/SE(3)'

    # Compute statistics for measured transforms
    pairs = data_batch.edge_index
    i, j = pairs
    se3_gt_i = data_batch['t_star'][i]
    se3_gt_j = data_batch['t_star'][j]
    rel_gt = se3_gt_j * se3_gt_i.inv()
    if 'valid_gt' in data_batch:
        valid_pairs = data_batch.valid_gt[i] & data_batch.valid_gt[j]
    else:
        valid_pairs = torch.ones(pairs.shape[1], dtype=torch.bool, device=pairs.device)
    # Compute measurement errors
    measured_errs = rel_gt.compare(data_batch['tij_hat'])
    metrics['measured_rot_errors'] = measured_errs['rot_deg'][valid_pairs]
    if 'trans' in measured_errs:
        metrics['measured_trans_errors'] = measured_errs['trans'][valid_pairs]

    # Compute prediction errors
    # First compute which pairs to evaluate
    if not check_only_good:
        # Check all pairise edges
        pairs_to_eval = []
        for iBatch in range(batch_size):
            # Compute all pairwise combinations
            in_batch = np.nonzero(to_numpy(data_batch.batch == iBatch))[0]
            pairs_to_eval.append(torch.tensor(list(itertools.combinations(in_batch, 2))))
        pairs_to_eval = torch.cat(pairs_to_eval, dim=0)
    else:
        # Hack to compute metrics just on good edges
        pairs_to_eval = data_batch.edge_index[:, data_batch['edge_sigma'] < 0.05].transpose(0, 1)

    # Then evaluate the pairs
    i, j = pairs_to_eval[:, 0], pairs_to_eval[:, 1]
    se3_gt_i = data_batch['t_star'][i]
    se3_gt_j = data_batch['t_star'][j]
    se3_pred_i = predicted_transforms[i]
    se3_pred_j = predicted_transforms[j]
    rel_gt = se3_gt_j * se3_gt_i.inv()
    rel_pred = se3_pred_j * se3_pred_i.inv()
    pred_errs = rel_gt.compare(rel_pred)

    metrics['rot_errors'] = pred_errs['rot_deg']
    if 'trans' in pred_errs:
        metrics['trans_errors'] = pred_errs['trans']

    # Handles pairs without groundtruth
    if 'valid_gt' in data_batch:
        valid_pairs = data_batch.valid_gt[i] & data_batch.valid_gt[j]
        metrics['rot_errors'] = pred_errs['rot_deg'][valid_pairs]
        if 'trans' in pred_errs:
            metrics['trans_errors'] = pred_errs['trans'][valid_pairs]

    return metrics

def get_aggregated_metrics(metrics, is_rel):
    if is_rel:
        return get_aggregated_metrics_rel(metrics)
    else:
        return get_aggregated_metrics_abs(metrics)

def get_aggregated_metrics_rel(metrics):
    rot_thresh_deg = [3, 5, 10, 30, 45]
    trans_thresh = [0.05, 0.1, 0.25, 0.5, 0.75]

    # rot errors
    input_stats = {
        'rot_error_mean': torch.mean(metrics['measured_rot_errors']),
        'rot_error_med': torch.median(metrics['measured_rot_errors']),
    }
    aggr_metrics = {
        'rot_error_mean': torch.mean(metrics['rot_errors']),
        'rot_error_med': torch.median(metrics['rot_errors'])
    }
    for t in rot_thresh_deg:
        aggr_metrics['rot_ecdf_{}'.format(t)] = torch.mean((metrics['rot_errors'] < t).float())

    # trans errors
    if 'trans_errors' in metrics:
        input_stats.update({
            'trans_error_mean': torch.mean(metrics['measured_trans_errors']),
            'trans_error_med': torch.median(metrics['measured_trans_errors']),
        })
        aggr_metrics.update({
            'trans_error_mean': torch.mean(metrics['trans_errors']),
            'trans_error_med': torch.median(metrics['trans_errors'])
        })
        for t in trans_thresh:
            aggr_metrics['trans_ecdf_{}'.format(t)] = torch.mean((metrics['trans_errors'] < t).float())

    return aggr_metrics, input_stats

def get_aggregated_metrics_abs(metrics):
    # rot errors
    input_stats = {
        'rot_error_mean': torch.mean(metrics['measured_rot_errors']),
        'rot_error_med': torch.median(metrics['measured_rot_errors']),
    }
    aggr_metrics = {
        'rot_error_mean': torch.mean(metrics['rot_deg_err_mean']),
        'rot_error_med': torch.mean(metrics['rot_deg_err_med'])
    }

    # trans errors
    if 'trans_errors' in metrics:
        input_stats.update({
            'trans_error_mean': torch.mean(metrics['measured_trans_errors']),
            'trans_error_med': torch.median(metrics['measured_trans_errors']),
        })
        aggr_metrics.update({
            'trans_error_mean': torch.mean(metrics['trans_err_mean']),
            'trans_error_med': torch.mean(metrics['trans_err_medi='])
        })

    return aggr_metrics, input_stats

class Accumulator(object):
    """Generic Accumulator"""
    def __init__(self):
        self._data = defaultdict(float)  # Weighted sum
        self.reset()

    def reset(self):
        """Clears accumulated statistics
        """
        self._data = defaultdict(list)  # Stores metrics separately then concat

    def update(self, data):
        for k in data:
            self._data[k].append(data[k])

    def get_aggregated(self):
        """Return concatenated data
        """

        if len(self._data) == 0 or len(self._data[next(iter(self._data.keys()))]) == 0:
            raise AssertionError('You need to call update() at least once before calling this')
        self._concat_all()

        data = {k: self._data[k][0] for k in self._data}
        return data

    def _concat_all(self):
        for k in self._data:
            if self._data[k][0].dim() == 0:
                self._data[k] = [torch.stack(self._data[k], dim=0)]
            else:
                self._data[k] = [torch.cat(self._data[k], dim=0)]

max_stars_repo_name: yewzijian/MultiReg

max_issues_repo_name: yewzijian/MultiReg

max_forks_repo_name: yewzijian/MultiReg


text:
import os
import time
from datetime import datetime
import numpy as np
import numpy.random as npr
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

import tensorflow as tf
print(tf.__version__)
from tensorflow.python.client import device_lib
device_lib.list_local_devices()

from tensorflow.contrib.eager.python import tfe
from scipy.stats import mode
from scipy.io import loadmat
from scipy.special import expit

from tensorflow.python.keras.constraints import nonneg

import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors

import argparse

# always reload dependency code
import importlib
import utils
import idgmmvae
import gmmvae
importlib.reload(utils)
importlib.reload(idgmmvae)
importlib.reload(gmmvae)
from utils import *
from idgmmvae import *
from gmmvae import *

config = tf.ConfigProto()
config.gpu_options.allow_growth = True  # dynamically grow the memory used on the GPU
config.log_device_placement = True  # to log device placement (on which device the operation ran)
# (nothing gets printed in Jupyter, only if you run it standalone)
config.gpu_options.per_process_gpu_memory_fraction = 0.90
config.allow_soft_placement = True

# enable eager mode
tf.enable_eager_execution(config=config)

# set random seed
timenow = int(datetime.utcnow().strftime('%Y%m%d%H%M%S%f'))
randseed = int(timenow % 10000)
print("Random seed", randseed)
tf.set_random_seed(randseed)
np.random.seed(randseed)

parser = argparse.ArgumentParser()

parser.add_argument('-dd', '--discretez_dim', \
        type=int, default=10)
parser.add_argument('-hd', '--h_dim', \
        type=int, default=512)
parser.add_argument('-cd', '--contiz_dim', \
        type=int, default=64)
parser.add_argument('-od', '--out_dim', \
        type=int, default=784)
parser.add_argument('-opt', '--optimizer', \
        type=str, default="adam")
parser.add_argument('-lr', '--lr', \
        type=float, default=0.001)
parser.add_argument('-mmt', '--momentum', \
        type=float, default=0.5)
parser.add_argument('-ep', '--num_epochs', \
        type=int, default=1000)
parser.add_argument('-l2', '--l2reg', \
        type=float, default=0.001)
parser.add_argument('-nfd', '--nf_dim', \
        type=int, default=128)
parser.add_argument('-nb', '--num_bijectors', \
        type=int, default=2)
parser.add_argument('-gc', '--grad_clip', \
        type=float, default=5.)
parser.add_argument('-nl', '--gen_model_hidden_layers', \
        type=int, default=1)
parser.add_argument('-dt', '--data', \
        type=str, default="mnist", choices=['mnist', 'omniglot', 'fashionmnist', 'pinwheel'])
parser.add_argument('-md', '--model', \
        type=str, default='idvae', choices=['idvae', 'vae', 'laggingvae', 'betavae', 'li_et_al'])
parser.add_argument('--betavae', type=float, default=0.2, help='the beta value of betavae')
parser.add_argument('-ss', '--subsample', \
        type=int, default=20000)
# parameter for Li et al. 2019
parser.add_argument("--fb", type=int, default=0,
                 help="0: no fb; 1: fb; 2: max(target_kl, kl) for each dimension")
parser.add_argument("--target_kl", type=float, default=7.,
                 help="target kl of the free bits trick, value 7 or 3 in Li et al.")

args, unknown = parser.parse_known_args()

discretez_dim = args.discretez_dim
h_dim = args.h_dim

momentum = args.momentum
lr = args.lr

num_epochs = args.num_epochs

batch_size = 100

nf_dim = args.nf_dim

l2reg = args.l2reg
regularizer = tf.keras.regularizers.l2(l2reg)

num_bijectors = args.num_bijectors
MODEL = 'RealNVP'
USE_BATCHNORM = False

grad_clip = args.grad_clip

qy_num_hidden_layers = 2
qz_num_hidden_layers = 2
gen_model_hidden_layers = args.gen_model_hidden_layers

dataname = args.data

train_data, train_labels, test_data, test_labels = load_data(dataname)

train_idx = npr.choice(train_data.shape[0], args.subsample, replace=False)
train_data = train_data[train_idx,:]
train_labels = train_labels[train_idx,:]

optimizer = set_optimizer(args.optimizer, lr, momentum)

contiz_dim = args.contiz_dim # the same as data dimension
out_dim = train_data.shape[1]

outfilename = 'out_'+dataname+'_'+args.model+str(timenow)

checkpoint_path = '../../res/'+dataname+'_gmm_' + args.model + 'lr' + str(args.lr) + 'l2reg' + str(args.l2reg) + 'seed' + str(timenow)+'/cp-{itr:04d}.ckpt'
checkpoint_dir = os.path.dirname(checkpoint_path)
if not os.path.exists(checkpoint_dir):
    os.makedirs(checkpoint_dir)

device = '/cpu:0' if tfe.num_gpus() == 0 else '/gpu:0'

# with tf.device(device):
# build model and optimizer
if args.model == 'idvae':
        kl_weight = 1.0
        model = idgmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers, gen_model_hidden_layers, USE_BATCHNORM)
elif args.model == 'vae':
        kl_weight = 1.0
        model = gmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers,  gen_model_hidden_layers, kl_weight=kl_weight, USE_BATCHNORM=USE_BATCHNORM)
elif args.model == 'betavae':
        kl_weight = args.betavae
        model = gmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers,  gen_model_hidden_layers, kl_weight=kl_weight, USE_BATCHNORM=USE_BATCHNORM)
elif args.model == 'li_et_al':
        kl_weight = 1.

        ae_model = gmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers,  gen_model_hidden_layers, kl_weight=0., USE_BATCHNORM=USE_BATCHNORM)

        vae_model = gmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers,  gen_model_hidden_layers, kl_weight=1., fb=2., target_kl=args.target_kl, USE_BATCHNORM=USE_BATCHNORM)
elif args.model == 'laggingvae':
        kl_weight = 1.
        vae_model1 = gmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers, gen_model_hidden_layers, kl_weight=0., USE_BATCHNORM=USE_BATCHNORM)
        vae_model2 = gmmvae(discretez_dim, contiz_dim, h_dim, out_dim, regularizer, num_bijectors, nf_dim, qy_num_hidden_layers, qz_num_hidden_layers,  gen_model_hidden_layers, kl_weight=1., USE_BATCHNORM=USE_BATCHNORM)
for outfile in [outfilename+'.log', outfilename+'_au_discrete.log', outfilename+'_au_conti.log']:

    with open(outfile, 'a') as f:
        f.write("starting")
        f.write("\n"+dataname+"_"+args.model+"\n")
        f.write("\noptimizer")
        f.write(str(args))
        f.write("\nl2reg" + str(l2reg))
        f.write("\nkl_weight" + str(kl_weight))
        f.write("\nnf_dim" + str(nf_dim))
        f.write("\nnf_model" + str(MODEL))
        f.write("\nnum_bijectors" + str(num_bijectors))
        f.write("\nbatch_norm" + str(USE_BATCHNORM))
        f.write("\ngrad_clip" + str(grad_clip))
        f.write("\nrand_seed" + str(randseed))
        f.write("\ngen_model_hidden_layers" + str(gen_model_hidden_layers))
        f.write("\n\n\n")

if args.model == 'idvae':
        trained_model = \
                train_idgmmvae(model, train_data, train_labels, test_data, test_labels, optimizer, batch_size, num_epochs, outfilename, device, grad_clip, checkpoint_path)
elif args.model == 'vae':
        trained_model = \
                train_gmmvae(model, train_data, train_labels, test_data, test_labels, optimizer, batch_size, num_epochs, outfilename, device, grad_clip, checkpoint_path)
elif args.model == 'betavae':
        trained_model = \
                train_gmmvae(model, train_data, train_labels, test_data, test_labels, optimizer, batch_size, num_epochs, outfilename, device, grad_clip, checkpoint_path)

elif args.model == 'li_et_al':
        print("initializing with autoencoder")
        trained_ae_model = \
                train_gmmvae(ae_model, train_data, train_labels, test_data, test_labels, optimizer, batch_size, 100, outfilename, device, grad_clip, checkpoint_path)
        latest = tf.train.latest_checkpoint(checkpoint_dir)
        vae_model.load_weights(latest)
        print("fitting variational autoencoder with KL truncation")
        trained_vae_model = \
                train_gmmvae(vae_model, train_data, train_labels, test_data, test_labels, optimizer, batch_size, num_epochs, outfilename, device, grad_clip, checkpoint_path)
elif args.model == 'laggingvae':
        print("perform lagging inference (aggressivesly optimize inference network")
        trained_vae_model1 = \
                train_gmmvae(vae_model1, train_data, train_labels, test_data, test_labels, optimizer, batch_size, 100, outfilename, device, grad_clip, checkpoint_path, lagging=1)

        latest = tf.train.latest_checkpoint(checkpoint_dir)
        print("perform standard inference")
        vae_model2.load_weights(latest)
        trained_vae_model2 = \
                train_gmmvae(vae_model2, train_data, train_labels, test_data, test_labels, optimizer, batch_size, num_epochs, outfilename, device, grad_clip, checkpoint_path, lagging=0)

max_stars_repo_name: yixinwang/lidvae-public

max_issues_repo_name: yixinwang/lidvae-public

max_forks_repo_name: yixinwang/lidvae-public

reason: import numpy,from scipy


text:
from __future__ import print_function, division, absolute_import

from numba import utils
from numba.bytecode import ByteCodeInst, CustomByteCode
from collections import defaultdict

def lift_loop(bytecode, dispatcher_factory):
    """Lift the top-level loops.

    Returns (outer, loops)
    ------------------------
    * outer: ByteCode of a copy of the loop-less function.
    * loops: a list of ByteCode of the loops.
    """
    outer = []
    loops = []
    # Discover variables references
    outer_rds, outer_wrs = find_varnames_uses(bytecode, iter(bytecode))
    # Separate loops and outer
    separate_loops(bytecode, outer, loops)

    # Prepend arguments as negative bytecode offset
    for a in bytecode.pysig.parameters:
        outer_wrs[a] = [-1] + outer_wrs[a]

    dispatchers = []
    outerlabels = set(bytecode.labels)
    outernames = list(bytecode.co_names)

    for loop in loops:
        args, rets = discover_args_and_returns(bytecode, loop, outer_rds,
                                               outer_wrs)

        disp = insert_loop_call(bytecode, loop, args,
                                outer, outerlabels, rets,
                                dispatcher_factory)
        dispatchers.append(disp)

    # Build outer bytecode
    codetable = utils.SortedMap((i.offset, i) for i in outer)
    outerbc = CustomByteCode(func=bytecode.func,
                             func_qualname=bytecode.func_qualname,
                             is_generator=bytecode.is_generator,
                             pysig=bytecode.pysig,
                             filename=bytecode.filename,
                             co_names=outernames,
                             co_varnames=bytecode.co_varnames,
                             co_consts=bytecode.co_consts,
                             co_freevars=bytecode.co_freevars,
                             table=codetable,
                             labels=outerlabels & set(codetable.keys()))

    return outerbc, dispatchers

@utils.total_ordering
class SubOffset(object):
    """The loop-jitting may insert bytecode between two bytecode but we
    cannot guarantee that there is enough integral space between two offsets.
    This class workaround the problem by introducing a fractional part to the
    offset.
    """
    def __init__(self, val, sub=1):
        assert sub > 0, "fractional part cannot be <= 0"
        self.val = val
        self.sub = sub

    def next(self):
        """Helper method to get the next suboffset by incrementing the
        fractional part only
        """
        return SubOffset(self.val, self.sub + 1)

    def __add__(self, other):
        """Adding to a suboffset will only increment the fractional part.
        The integral part is immutable.
        """
        return SubOffset(self.val, self.sub + other)

    def __hash__(self):
        return hash((self.val, self.sub))

    def __lt__(self, other):
        """Can only compare to SubOffset or int
        """
        if isinstance(other, SubOffset):
            if self.val < other.val:
                return self
            elif self.val == other.val:
                return self.sub < other.sub
            else:
                return False
        elif isinstance(other, int):
            return self.val < other
        else:
            return NotImplemented

    def __eq__(self, other):
        if isinstance(other, SubOffset):
            return self.val == other.val and self.sub == other.sub
        elif isinstance(other, int):
            # Can never be equal to a integer by definition
            return False
        else:
            return NotImplemented

    def __repr__(self):
        """Print like a floating-point by it is not one at all.
        """
        return "{0}.{1}".format(self.val, self.sub)

def insert_loop_call(bytecode, loop, args, outer, outerlabels, returns,
                     dispatcher_factory):
    endloopoffset = loop[-1].next
    # Accepted. Create a bytecode object for the loop
    args = tuple(args)

    lbc = make_loop_bytecode(bytecode, loop, args, returns)

    # Generate dispatcher for this inner loop, and append it to the
    # consts tuple.
    disp = dispatcher_factory(lbc)
    disp_idx = len(bytecode.co_consts)
    bytecode.co_consts += (disp,)

    # Insert jump to the end
    insertpt = SubOffset(loop[0].next)
    jmp = ByteCodeInst.get(loop[0].offset, 'JUMP_ABSOLUTE', insertpt)
    jmp.lineno = loop[0].lineno
    insert_instruction(outer, jmp)

    outerlabels.add(outer[-1].next)

    # Prepare arguments
    loadfn = ByteCodeInst.get(insertpt, "LOAD_CONST", disp_idx)
    loadfn.lineno = loop[0].lineno
    insert_instruction(outer, loadfn)

    insertpt = insertpt.next()
    for arg in args:
        loadarg = ByteCodeInst.get(insertpt, 'LOAD_FAST',
                                   bytecode.co_varnames.index(arg))
        loadarg.lineno = loop[0].lineno
        insert_instruction(outer, loadarg)
        insertpt = insertpt.next()

    # Call function
    assert len(args) < 256
    call = ByteCodeInst.get(insertpt, "CALL_FUNCTION", len(args))
    call.lineno = loop[0].lineno
    insert_instruction(outer, call)

    insertpt = insertpt.next()

    if returns:
        # Unpack arguments
        unpackseq = ByteCodeInst.get(insertpt, "UNPACK_SEQUENCE",
                                  len(returns))
        unpackseq.lineno = loop[0].lineno
        insert_instruction(outer, unpackseq)
        insertpt = insertpt.next()

        for out in returns:
            # Store each variable
            storefast = ByteCodeInst.get(insertpt, "STORE_FAST",
                                      bytecode.co_varnames.index(out))
            storefast.lineno = loop[0].lineno
            insert_instruction(outer, storefast)
            insertpt = insertpt.next()
    else:
        # No return value
        poptop = ByteCodeInst.get(insertpt, "POP_TOP", None)
        poptop.lineno = loop[0].lineno
        insert_instruction(outer, poptop)
        insertpt = insertpt.next()

    jmpback = ByteCodeInst.get(insertpt, 'JUMP_ABSOLUTE',
                               endloopoffset)

    jmpback.lineno = loop[0].lineno
    insert_instruction(outer, jmpback)

    return disp

def insert_instruction(insts, item):
    i = find_previous_inst(insts, item.offset)
    insts.insert(i, item)

def find_previous_inst(insts, offset):
    for i, inst in enumerate(insts):
        if inst.offset > offset:
            return i
    return len(insts)

def make_loop_bytecode(bytecode, loop, args, returns):
    # Add return None
    co_consts = tuple(bytecode.co_consts)
    if None not in co_consts:
        co_consts += (None,)

    if returns:
        for out in returns:
            # Load output
            loadfast = ByteCodeInst.get(loop[-1].next, "LOAD_FAST",
                                         bytecode.co_varnames.index(out))
            loadfast.lineno = loop[-1].lineno
            loop.append(loadfast)

        # Build tuple
        buildtuple = ByteCodeInst.get(loop[-1].next, "BUILD_TUPLE",
                                    len(returns))
        buildtuple.lineno = loop[-1].lineno
        loop.append(buildtuple)

    else:
        # Load None
        load_none = ByteCodeInst.get(loop[-1].next, "LOAD_CONST",
                                     co_consts.index(None))
        load_none.lineno = loop[-1].lineno
        loop.append(load_none)

    # Return TOS
    return_value = ByteCodeInst.get(loop[-1].next, "RETURN_VALUE", 0)
    return_value.lineno = loop[-1].lineno
    loop.append(return_value)

    # Function name
    loop_qualname = bytecode.func_qualname + ".__numba__loop%d__" % loop[0].offset

    # Code table
    codetable = utils.SortedMap((i.offset, i) for i in loop)

    # Custom bytecode object
    lbc = CustomByteCode(func=bytecode.func,
                         func_qualname=loop_qualname,
                         # Enforced in separate_loops()
                         is_generator=False,
                         pysig=bytecode.pysig,
                         arg_count=len(args),
                         arg_names=args,
                         filename=bytecode.filename,
                         co_names=bytecode.co_names,
                         co_varnames=bytecode.co_varnames,
                         co_consts=co_consts,
                         co_freevars=bytecode.co_freevars,
                         table=codetable,
                         labels=bytecode.labels)

    return lbc

def stitch_instructions(outer, loop):
    begin = loop[0].offset
    i = find_previous_inst(outer, begin)
    return outer[:i] + loop + outer[i:]

def remove_from_outer_use(inneruse, outeruse):
    for name in inneruse:
        inuse = inneruse[name]
        outuse = outeruse[name]
        outeruse[name] = sorted(list(set(outuse) - set(inuse)))

def discover_args_and_returns(bytecode, insts, outer_rds, outer_wrs):
    """
    Basic analysis for args and returns
    This completely ignores the ordering or the read-writes.

    outer_rds and outer_wrs are modified

    Note:
    An invalid argument and return set will likely to cause a RuntimeWarning
    in the dataflow analysis due to mismatch in stack offset.
    """
    rdnames, wrnames = find_varnames_uses(bytecode, insts)

    # Remove all local use from the set
    remove_from_outer_use(rdnames, outer_rds)
    remove_from_outer_use(wrnames, outer_wrs)

    # Return every variables that are written inside the loop and read
    # afterwards
    rets = set()
    for name, uselist in wrnames.items():
        if name in outer_rds:
            endofloop = insts[-1].offset

            # Find the next read
            for nextrd in outer_rds[name]:
                if nextrd > endofloop:
                    break
            else:
                nextrd = None

            # Find the next write
            for nextwr in outer_wrs[name]:
                if nextwr > endofloop:
                    break
            else:
                nextwr = None

            # If there is a read but no write OR
            # If the next use is a read, THEN
            # it is a return value
            if nextrd is not None and (nextwr is None or nextwr > nextrd):
                rets.add(name)

    # Make variables arguments if they are read before defined before the loop.
    # Since we can't tell if things are conditionally defined here,
    # We will have to be more conservative.
    args = set()
    firstline = insts[0].offset
    for name in rdnames.keys():
        outer_write = outer_wrs[name]
        # If there exists a definition before the start of the loop
        # for a variable read in side the loop.
        if any(i < firstline for i in outer_write):
            args.add(name)

    # Make variables arguments if it is being returned but defined before the
    # loop.
    for name in rets:
        if any(i < insts[0].offset for i in outer_wrs[name]):
            args.add(name)

    # Re-add the arguments back to outer_rds
    for name in args:
        outer_rds[name] = sorted(set(outer_rds[name]) | set([firstline]))

    # Re-add the arguments back to outer_wrs
    for name in rets:
        outer_wrs[name] = sorted(set(outer_wrs[name]) | set([firstline]))

    return args, rets

def find_varnames_uses(bytecode, insts):
    rdnames = defaultdict(list)
    wrnames = defaultdict(list)
    for inst in insts:
        if inst.opname == 'LOAD_FAST':
            name = bytecode.co_varnames[inst.arg]
            rdnames[name].append(inst.offset)
        elif inst.opname == 'STORE_FAST':
            name = bytecode.co_varnames[inst.arg]
            wrnames[name].append(inst.offset)
    return rdnames, wrnames

def separate_loops(bytecode, outer, loops):
    """
    Separate top-level loops from the function

    Stores loopless instructions from the original function into `outer`.
    Stores list of loop instructions into `loops`.
    Both `outer` and `loops` are list-like (`append(item)` defined).
    """
    # XXX When an outer loop is rejected, there may be an inner loop
    # which would still allow lifting.
    endloop = None
    cur = None
    for inst in bytecode:
        if endloop is None:
            if inst.opname == 'SETUP_LOOP':
                cur = [inst]
                # Python may set the end of loop to the final jump destination
                # when nested in a if-else.  We need to scan the bytecode to
                # find the actual end of loop
                endloop = _scan_real_end_loop(bytecode, inst)
            else:
                outer.append(inst)
        else:
            cur.append(inst)
            if inst.next == endloop:
                for inst in cur:
                    if inst.opname in ['RETURN_VALUE', 'YIELD_VALUE',
                                       'BREAK_LOOP']:
                        # Reject if return, yield or break inside loop
                        outer.extend(cur)
                        break
                else:
                    loops.append(cur)
                endloop = None

def _scan_real_end_loop(bytecode, setuploop_inst):
    """Find the end of loop.
    Return the instruction offset.
    """
    start = setuploop_inst.next
    end = start + setuploop_inst.arg
    offset = start
    depth = 0
    while offset < end:
        inst = bytecode[offset]
        depth += inst.block_effect
        if depth < 0:
            return inst.next
        offset = inst.next

max_stars_repo_name: wangyum/anaconda

max_issues_repo_name: wangyum/anaconda

max_forks_repo_name: wangyum/anaconda

reason: from numba


text:
import sys
sys.path.append('../codes')
from options import get_options
import os
import networkx as nx
from networkx import topological_sort
import torch as th
import pickle
import torch.nn as nn
from random import shuffle

class Dataset(th.utils.data.Dataset):
    def __init__(self, data):
        # print("init dataset, len=", len(data))
        # print(torch.max(data))
        self.data = data

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

    def __getitem__(self, idx):
        #print(self.data[idx][1].shape)
        return (self.data[idx][0], self.data[idx][1])

class EVCNN(nn.Module):
    def __init__(self, in_dim, hidden_dim):
        super().__init__()
        self.bn = nn.BatchNorm1d(in_dim)
        self.cnn = nn.Sequential(
            nn.Conv1d(hidden_dim, 64, 8), nn.ReLU(), nn.MaxPool1d(2), nn.Dropout(0.25),
        )
        self.fc = nn.Sequential(
            nn.Linear(int(64 * (in_dim - 8) / 2), 32),
            nn.ReLU(),
            nn.Dropout(0.5),
            nn.Linear(32, 4),
        )

    def forward(self, x):
        # x = self.bn(x)
        # print(x)
        x = self.cnn(x)
        # print(x)
        x = x.view(x.size(0), -1)
        # print(x.size())
        x = self.fc(x)
        # print(x)
        # input()
        return x

def top_k(evs,k):
    evs2 = th.clone(evs)
    evs2[evs2>0]=1
    count = th.sum(evs2,dim=1)
    _,order = th.sort(count)
    topk = reversed(order)[:k]
    return evs[topk]

def construct_ev(graph):
    evs = None
    max_t= 0
    nxg = graph.to_networkx()
    for i in  range(graph.number_of_nodes()):
        ev = th.clone(graph.ndata["ntype"][i]).detach()
        l = nx.single_source_shortest_path_length(nxg, i, cutoff=2)
        for n, d in l.items():
            ev += graph.ndata["ntype"][n] * (1 if d == 1 else 0.5)
        l = nx.single_source_shortest_path_length(
            nxg.reverse(copy=False), i, cutoff=2
        )
        for n, d in l.items():
            ev += graph.ndata["ntype"][n] * (1 if d == 1 else 0.5)
        max_t = max(max_t, th.max(ev).item())
        ev = th.clamp(ev, 0, 64).unsqueeze(0)
        if evs is None:
            evs = ev
        else:
            evs = th.cat((evs,ev),dim=0)
    return evs
def build_feature(graph,p,k):
    evs = construct_ev(graph)
    feature = None
    order = list(topological_sort(graph.to_networkx()))
    #print(order)
    group_size = int(len(order)/p)
    groups = []
    if len(order)<p*k:
        added = 0
        feature = evs
        while len(feature)<p*k:
            feature = th.cat((feature,evs[:p*k-len(feature)]))
        #feature = th.cat((evs,evs[:p*k-len(order)]))
        #print(len(feature))
    else:
        for i in range(p):
            if i ==p-1:
                group = order[i*group_size:]
            else:
                group = order[i*group_size:(i+1)*group_size]
            group_evs = evs[group]
            ktop_evs = top_k(group_evs,k)
            if feature is None:
                feature = ktop_evs
            else:
                feature = th.cat((feature,ktop_evs))
    if feature.shape[0]!=120:
        print(graph)
        exit()
    return feature

def load_data(path,p,k):
    if os.path.exists(path):
        with open(path,'rb') as f:
            train_data,val_data = pickle.load(f)
    else:
        with open(os.path.join('../data/global_new', 'train_data.pkl'), 'rb') as f:
            train_dataset = pickle.load(f)
        with open(os.path.join('../data/global_new', 'val_data.pkl'), 'rb') as f:
            val_dataset = pickle.load(f)

        train_data = []
        val_data = []
        train_labels =  th.tensor([],dtype=th.long)
        val_labels = th.tensor([],dtype=th.long)
        remove_train = 0
        remove_val = 0
        for i,(label, g, _, _) in enumerate(train_dataset):

            # if g.number_of_nodes() < p*k:
            #     remove_train += 1
            #     continue

            feature = build_feature(g, p, k)
            train_data.append((label, feature))
            print(i,feature.shape,feature, '\n')

        for i, (label, g, _, _) in enumerate(val_dataset):

            # if g.number_of_nodes() < p*k:
            #     remove_val += 1
            #     continue

            feature = build_feature(g, p, k)
            val_data.append((label, feature))
            print(i, feature.shape,feature, '\n')

        print(remove_train, remove_val)
        os.makedirs('../data/evcnn/',exist_ok=True)
        with open(path, 'wb') as f:
            pickle.dump((train_data, val_data), f)

    return train_data,val_data

def validate(val_dataloader, device, model,loss, options):
    total_num, total_loss, correct, fn, fp, tn, tp = 0, 0.0, 0, 0, 0, 0, 0

    with th.no_grad():
        for i, (labels,features)  in enumerate(val_dataloader):
            features = features.to(device)
            labels = labels.to(device)
            label_hats = model(features)
            predict_labels = th.argmax(nn.functional.softmax(label_hats, 1), dim=1)
        # print('ground-truth labels:',labels.shape,labels)
        # print('predict labels:',predict_labels.shape,predict_labels)
            val_loss = loss(label_hats, labels)
        #print('val loss:', val_loss.item())
            total_num += len(labels)

            correct += (
                    predict_labels == labels
            ).sum().item()
            #errors = th.tensor(range(0, len(predict_labels)))[predict_labels != labels]
            #print('errors:', errors)
            #print(correct,total_num)
            val_acc = correct / total_num

    return [val_loss, val_acc]

def train():
    options = get_options()
    p,k = 40,3

    train_data,val_data = load_data('../data/evcnn/{}.pkl'.format(options.datapath),p,k)
    device = th.device("cuda"  if th.cuda.is_available() else "cpu")
    #shuffle(train_dataset)

    if os.path.exists(options.model_saving_dir):
        with open(os.path.join(options.model_saving_dir,'model.pkl'),'rb') as f:
            options,model = pickle.load(f)
    else:
        os.makedirs(options.model_saving_dir, exist_ok=True)
        if options.datapath=='data2':
            p,k=40,2
        else:
            p,k=40,3
        model =EVCNN(16,p*k)
        with open(os.path.join(options.model_saving_dir,'model.pkl'),'wb') as f:
            pickle.dump((options,model),f)
        with open(os.path.join(options.model_saving_dir, 'res.txt'),'wb') as f:
            pass
    model.to(device)

    Loss = nn.CrossEntropyLoss()
    optim = th.optim.Adam(
        model.parameters(),
        options.learning_rate, weight_decay=options.weight_decay
    )
    model.train()

    max_acc = 0
    # for lb,ft in train_data:
    #     print(ft.shape)
    # exit()
    shuffle(train_data)
    train_data = train_data[:int(len(val_data)*options.train_percent)]
    print('num train:',len(train_data))
    print('num_val:',len(val_data))
    train_dataset = Dataset(train_data)
    train_dataloader = th.utils.data.DataLoader(train_dataset, batch_size=options.batch_size)
    val_dataset = Dataset(val_data)
    val_dataloader = th.utils.data.DataLoader(val_dataset, batch_size=len(val_dataset))
    print('Start training ')
    for epoch in range(options.num_epoch):
        total_num, total_loss, correct, fn, fp, tn, tp = 0, 0.0, 0, 0, 0, 0, 0
        shuffle(train_data)
        labels = th.tensor([], dtype=th.long).to(device)
        label_hats = None
        for i, (labels,features)  in enumerate(train_dataloader):
            #print(labels.shape,features.shape)
            #labels = th.cat((labels, th.tensor([label], dtype=th.long).to(device)))
            features = features.to(device)
            labels = labels.to(device)
            label_hats = model(features)
            #print(label_hats.shape)
            #print(labels,label_hats)
            # if label_hats is None:
            #     label_hats = label_hat.unsqueeze(0)
            # else:
            #     label_hats = th.cat((label_hats,label_hat),dim=0)

            predict_labels = th.argmax(nn.functional.softmax(label_hats, 1), dim=1)
            # print('ground-truth labels:',labels.shape,labels)
            # print('predict labels:',predict_labels.shape,predict_labels)
            train_loss = Loss(label_hats, labels)
            #print('loss:', train_loss.item())
            total_num += len(labels)
            total_loss += train_loss.item() * len(labels)
            correct += (
                    predict_labels == labels
            ).sum().item()
            optim.zero_grad()
            train_loss.backward()
            # print(model.GCN1.layers[0].attn_n.grad)
            optim.step()

        Train_loss = total_loss / total_num
        Train_acc = correct / len(train_data)
        # print("epoch[{:d}]".format(epoch))
        # print("  train:")
        # # print("\ttp:", tp, " fp:", fp, " fn:", fn, " tn:", tn, " precision:", round(Train_precision,3))
        # print("\tloss:{:.8f}, acc:{:.3f}, ".format(Train_loss, Train_acc))
        val_loss, val_acc = validate(val_dataloader, device, model,Loss, options)
        #if not os.path.exists(os.path.join(options.model_saving_dir, 'res.txt'))
        with open(os.path.join(options.model_saving_dir, 'res.txt'), 'a') as f:
            f.write(str(round(Train_loss, 8)) + " " + str(round(Train_acc, 3))  + "\n")
            f.write(str(round(val_loss.item(), 3)) + " " + str(round(val_acc, 3)) + "\n")
            f.write('\n')

        judgement = val_acc > max_acc
        #judgement = True
        if judgement:
           print("epoch {}".format(epoch))
           print("\tval loss:{:.3f}, acc:{:.3f}".format(val_loss, val_acc))
           max_acc = val_acc
           print("Saving model.... ", os.path.join(options.model_saving_dir))
           if os.path.exists(options.model_saving_dir) is False:
              os.makedirs(options.model_saving_dir)
           with open(os.path.join(options.model_saving_dir, 'model.pkl'), 'wb') as f:
              parameters = options
              pickle.dump((parameters, model), f)
           print("Model successfully saved")

if __name__ == "__main__":
    seed = 1234
    # th.set_deterministic(True)
    th.manual_seed(seed)
    th.cuda.manual_seed(seed)
    train()

max_stars_repo_name: ZeayW/graph-contrastive-learning

max_issues_repo_name: ZeayW/graph-contrastive-learning

max_forks_repo_name: ZeayW/graph-contrastive-learning

reason: import networkx,from networkx


text:
# ---
# jupyter:
#   jupytext:
#     formats: py:percent
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.11.1
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# ![SkeideLab and MPI CBS logos](../misc/header_logos.png)
#
# # Notebook #06: Fail-Safe N Analysis
#
# *Created April 2021 by Alexander Enge* ([enge@cbs.mpg.de](mailto:enge@cbs.mpg.de))
#
# This notebook describes another type of robustness check for our meta-analytic results. Here, we ask what happens with our inferences under the (not unlikely) assumption that there are additional studies with inconsistent findings that we've failed to include because they never got published (also called the "file drawer" effect). As for the jackknife analysis, we cannot directly test or correct for the presence of this bias. We can, however, perform simulations to estimate the impact that such a bias would have *if it existed*. The logic is simple and has been used for classical (i.e., non-neuroimaging) meta-analyses for decades (Rosenthal et al., 1979, *Psychol Bull*): Just create fictional studies with non-significant effect sizes and see how many of these can be added to your meta-analysis until the statistical significance of your effect size is overturned. If this number (called fail-safe N or FSN) exceeds any plausible guess of studies that *could* actually be in the file drawer, it would speak to the robustness of the meta-analytic results. Acar and colleagues (2018; *PLOS ONE*) have extended this logic to ALE-based meta-analysis of neuroimaging experiments: Just create new (fictional) experiments with inconsistent results and see how many of these can added to the original ALE until the significance of your clusters is overturned. Here we'll perform such an analysis for our four ALEs of interest (main analysis, knowledge, relatedness, and objects).
#
# As always, we start by loading the relevant packages.

# %%
import random
from glob import glob
from math import sqrt
from os import makedirs, path
from re import sub
from shutil import copy
from sys import argv

import numpy as np
import pandas as pd
from nibabel import save
from nilearn import image, reporting
from nimare import correct, io, meta, utils
from scipy.stats import norm

# %% [markdown]
# Our FSN simulations will heavily rely on the generation of so called "null experiments," i.e., fictional experiments with their peaks randomly distributed across the brain. We'll start by writing a helper function for creating these. It takes as its input a "real" ALE data set (in the form of a Sleuth text file, see Notebook #1 and [this example](http://www.brainmap.org/ale/foci2.txt)). It then creates a desired number (`k_null`) of null experiments that are similar to the real experiments in terms of sample sizes and number of peak coordinates. However, the location of these peaks is drawn randomly from all voxels within our gray matter template. For these experiments, we know that the null hypothesis (i.e., no spatial convergence) is true, thus providing a testing ground for the file drawer effect.

# %%
# Define function to generate a new data set with k null studies added
def generate_null(
    text_file="peaks.txt",
    space="ale_2mm",
    k_null=100,
    random_seed=None,
    output_dir="./",
):

    # Load NiMARE's gray matter template
    temp = utils.get_template(space=space, mask="gm")

    # Extract possible MNI coordinates for all gray matter voxels
    x, y, z = np.where(temp.get_fdata() == 1.0)
    within_mni = image.coord_transform(x=x, y=y, z=z, affine=temp.affine)
    within_mni = np.array(within_mni).transpose()

    # Read the original Sleuth file into a NiMARE data set
    dset = io.convert_sleuth_to_dataset(text_file, target=space)

    # Set a random seed to make the results reproducible
    if random_seed:
        random.seed(random_seed)

    # Resample numbers of subjects per experiment based on the original data
    nr_subjects_dset = [n[0] for n in dset.metadata["sample_sizes"]]
    nr_subjects_null = random.choices(nr_subjects_dset, k=k_null)

    # Resample numbers of peaks per experiment based on the original data
    nr_peaks_dset = dset.coordinates["study_id"].value_counts().tolist()
    nr_peaks_null = random.choices(nr_peaks_dset, k=k_null)

    # Create random peak coordinates
    idx_list = [
        random.sample(range(len(within_mni)), k=k_peaks) for k_peaks in nr_peaks_null
    ]
    peaks_null = [within_mni[idx] for idx in idx_list]

    # Copy original experiments to the destination Sleuth file
    makedirs(output_dir, exist_ok=True)
    text_file_basename = path.basename(text_file)
    null_file_basename = sub(
        pattern=".txt", repl="_plus_k" + str(k_null) + ".txt", string=text_file_basename
    )
    null_file = output_dir + "/" + null_file_basename
    copy(text_file, null_file)

    # Append all the null studies to the Sleuth file
    f = open(null_file, mode="a")
    for i in range(k_null):
        f.write(
            "\n// nullstudy"
            + str(i + 1)
            + "\n// Subjects="
            + str(nr_subjects_null[i])
            + "\n"
        )
        np.savetxt(f, peaks_null[i], fmt="%.3f", delimiter="\t")
    f.close()

    # Read the new Sleuth file and return it as a NiMARE data set
    dset_null = io.convert_sleuth_to_dataset(null_file, target=space)
    return dset_null

# %% [markdown]
# With this helper function, we can go on to implement the actual FSN simulation. This function will be rather long-winded and complex, but the overall logic is simple: Take the Sleuth file from the original ALE analysis, generate a large number of null experimetns, and add them iteratively to the analysis. At each step, re-estimate the ALE and record which voxel have remained statistically significant and which have not. Based on this, each (initially singificant) voxel will be assigned an FSN value which is equal to the highest number of null experiments that we were able to add before the voxel failed to reach significant for the first time. Because these simulations are going to take a really (!) long time, we implement two stopping rules: We either stop if there are no more significant voxels left *or* if we've added a sufficiently large number of null studies for our purpose (e.g., five times times the number of original studies in the meta-analysis).

# %%
# Define function to compute the FSN for all voxels from a Sleuth file
def compute_fsn(
    text_file="peaks.txt",
    space="ale_2mm",
    voxel_thresh=0.001,
    cluster_thresh=0.01,
    n_iters=1000,
    k_max_factor=5,
    random_ale_seed=None,
    random_null_seed=None,
    output_dir="./",
):

    # Let's show the user what we are doing
    print("\nCOMPUTING FSN FOR " + text_file + " (seed: " + str(random_null_seed) + ")")

    # Set random seed for original ALE if requested
    if random_ale_seed:
        np.random.seed(random_ale_seed)

    # Recreate the original ALE analysis
    ale = meta.cbma.ALE()
    corr = correct.FWECorrector(
        method="montecarlo", voxel_thresh=voxel_thresh, n_iters=n_iters
    )
    dset_orig = io.convert_sleuth_to_dataset(text_file=text_file, target=space)
    res_orig = ale.fit(dset_orig)
    cres_orig = corr.transform(res_orig)

    # Extract the original study IDs
    ids_orig = dset_orig.ids.tolist()

    # Create a new data set with a large number null studies added
    k_max = len(ids_orig) * k_max_factor
    dset_null = generate_null(
        text_file=text_file,
        space=space,
        k_null=k_max,
        random_seed=random_null_seed,
        output_dir=output_dir,
    )

    # Create thresholded cluster mask
    img_fsn = cres_orig.get_map("z_level-cluster_corr-FWE_method-montecarlo")
    cluster_thresh_z = norm.ppf(1 - cluster_thresh / 2)
    img_fsn = image.threshold_img(img_fsn, threshold=cluster_thresh_z)
    img_fsn = image.math_img("np.where(img > 0, 1, 0)", img=img_fsn)

    # Create cluster-thresholded z map
    img_z = cres_orig.get_map("z")
    img_z = image.math_img("img1 * img2", img1=img_fsn, img2=img_z)

    # Iteratively add null studies up to our pre-defined maximum
    for k in range(1, k_max):

        # Print message
        print("Computing ALE for k = " + str(k) + " null studies added...")

        # Create a new data set with k null studies added
        ids_null = ["nullstudy" + str(x) + "-" for x in range(1, k + 1)]
        ids = ids_orig + ids_null
        dset_k = dset_null.slice(ids)

        # Compute the ALE
        res_k = res = ale.fit(dset_k)
        cres_k = corr.transform(result=res_k)

        # Create a thresholded cluster mask
        img_k = cres_k.get_map("z_level-cluster_corr-FWE_method-montecarlo")
        img_k = image.threshold_img(img_k, threshold=cluster_thresh_z)
        img_k = image.math_img("np.where(img > 0, 1, 0)", img=img_k)

        # Use this to update the per-voxel FSN - this is a bit hack-ish: On a voxel-by-
        # voxel basis, we increase the value by 1 if and only if the voxel has remained
        # significant. As soon as it has failed to reach significance once, we never
        # increase FSN any further. This is handeled by comparing the current FSN to
        # the current value of k.
        count = str(k + 1)
        formula = "np.where(img_fsn + img_k == " + count + ", img_fsn + 1, img_fsn)"
        img_fsn = image.math_img(formula, img_fsn=img_fsn, img_k=img_k)

        # Quit as soon as there are no significant clusters left in the current map
        if not np.any(img_k.get_fdata()):
            print("No more significant voxels - terminating\n")
            break

    # Save the FSN map that we've created in the loop
    filename_img = path.basename(text_file).replace(".txt", "_fsn.nii.gz")
    save(img_fsn, filename=output_dir + "/" + filename_img)

    # Extract the FSN values at the original cluster peaks
    tab_fsn = reporting.get_clusters_table(img_z, stat_threshold=0, min_distance=1000)
    inv_affine = np.linalg.inv(img_z.affine)
    x, y, z = [np.array(tab_fsn[col]) for col in ["X", "Y", "Z"]]
    x, y, z = image.coord_transform(x=x, y=y, z=z, affine=inv_affine)
    x, y, z = [arr.astype("int") for arr in [x, y, z]]
    tab_fsn["FSN"] = img_fsn.get_fdata()[x, y, z]

    # Save this cluster table with the new FSN column
    filename_tab = path.basename(text_file).replace(".txt", "_fsn.tsv")
    tab_fsn.to_csv(output_dir + "/" + filename_tab, sep="\t", index=False)

    return img_fsn, tab_fsn

# %% [markdown]
# Ideally, we want to perform all of this multiple times for different (random) filedrawers. Otherwise, the resulting FSN values would hinge a lot on the random patterns of these specific null experiments. However, doing all of these iterative simulations multiple times is extremly computationally expensive. We therefore wrote the next bit of the notebook in a way so that it can be run in parallel on our high performance computing (HPC) cluster. For this, we would call this notebook as a Python script from the command line and need to provide it with two additional parameters: The name of the original ALE analysis for which we want to compute the FSN and the number of different filedrawers we want to estimate (so we can always compute multiple filedrawers in parallel). If you don't happen to have access to an HPC and want to try the out the simulations directly withing the notebook, simply uncomment the two lines of code to define `prefixes` and `nr_filedrawers` locally.

# %%
# Get which FSN analyses to perform from the command line
prefixes = argv[1].split(",")

# Get number of filedrawers per analysis from the command line
nr_filedrawers = int(argv[2])

# # Or define them here for debugging
# prefixes = ["all", "knowledge", "relatedness", "objects"]
# nr_filedrawers = 10

# List the filenames of the Sleuth text files
text_files = ["../results/ale/" + prefix + ".txt" for prefix in prefixes]

# Create output directory names
output_dirs = ["../results/fsn/" + prefix + "/" for prefix in prefixes]

# Create random seeds for filedrawers
random_null_seeds = random.sample(range(1000), k=nr_filedrawers)
filedrawers = ["filedrawer" + str(seed) for seed in random_null_seeds]

# %% [markdown]
# The actual simulations are happening here: For each of the input text file (one for each original ALE analysis), we compute one FSN map for a number of different filedrawers (in our paper, we used ten). The results for each file drawer will be stored in separate folder which also indicates the random seed value (e.g., `166`) that was used to generate the null experiments (e.g., `results/fsn/all/filedrawer166`).

# %%
# Use our function to compute multiple filedrawers for each text file
_ = [
    [
        compute_fsn(
            text_file=text_file,
            space="ale_2mm",
            voxel_thresh=0.001,
            cluster_thresh=0.01,
            n_iters=1000,
            k_max_factor=5,
            random_ale_seed=1234,
            random_null_seed=random_null_seed,
            output_dir=output_dir + filedrawer,
        )
        for random_null_seed, filedrawer in zip(random_null_seeds, filedrawers)
    ]
    for text_file, output_dir in zip(text_files, output_dirs)
]

# %% [markdown]
# Once the simulation are done, we need to aggregate the results that we've obtained for each analysis across multiple file drawers. Remember that for each file drawer, we've already stored an FSN map as well as a table that contains the FSN value for each of the original cluster peaks. Here we just average these maps and perform some summary statistics on the tables, such as computing the mean FSN and its 95% confidence interval across the multiple file drawers (in our case, ten).

# %%
# Compute mean FSN across filedrawers
for prefix in prefixes:

    # Read FSN maps from all filedrawers
    fnames_maps = glob(
        "../results/fsn/" + prefix + "/filedrawer*/" + prefix + "_fsn.nii.gz"
    )
    imgs_fsn = [image.load_img(fname) for fname in fnames_maps]

    # Average and save
    img_mean = image.mean_img(imgs_fsn)
    fname_img_mean = "../results/fsn/" + prefix + "/" + prefix + "_mean_fsn.nii.gz"
    save(img_mean, fname_img_mean)

    # Read FSN tables from all filedrawers
    fnames_tabs = glob(
        "../results/fsn/" + prefix + "/filedrawer*/" + prefix + "_fsn.tsv"
    )
    tabs_fsn = [pd.read_csv(fname, delimiter="\t") for fname in fnames_tabs]
    tab_fsn = pd.concat(tabs_fsn)

    # Compute summary statistics
    agg = tab_fsn.groupby("Cluster ID")["FSN"].agg(["mean", "count", "std"])

    # Compute confidence intervals
    ci_level = 0.05
    z_crit = abs(norm.ppf(ci_level / 2))
    agg["se"] = [std / sqrt(count) for std, count in zip(agg["std"], agg["count"])]
    agg["ci_lower"] = agg["mean"] - z_crit * agg["se"]
    agg["ci_upper"] = agg["mean"] + z_crit * agg["se"]

    # Save summary statistics
    fname_agg = "../results/fsn/" + prefix + "/" + prefix + "_mean_fsn.csv"
    agg.to_csv(fname_agg, float_format="%.3f")


text:
import json
import numpy as np
from cv2 import VideoWriter, VideoWriter_fourcc
import cv2
import librosa
import os
from PIL import Image
from modules.background import Background
from modules.widgets import *
from modules.tools import GetUniqueID
import moviepy.editor as mp
from math import ceil

def make_video(audio_file, widgetstring, background=None, fps=24, width=1280, height=720):
    # sets the name of the destination file
    video_id = GetUniqueID("GENERATED_VIDEOS", ".mp4")
    filename = video_id + ".mp4"
    # loads the audio in time series and its sample rate
    time_series, sample_rate = librosa.load(audio_file)
    # gets the absolutae values of the stft operation
    stft = np.abs(librosa.stft(time_series, hop_length=512, n_fft=2048 * 4))
    # convert the amplitude based diagram into a db-based one
    spectrogram = librosa.amplitude_to_db(stft, ref=np.max)
    # gets the Discrete Fourier Transform sample frequencies
    frequencies = librosa.core.fft_frequencies(n_fft=2048 * 4)
    times = librosa.core.frames_to_time(np.arange(
        spectrogram.shape[1]), sr=sample_rate, hop_length=512, n_fft=2048 * 4)  # get each frame time's in seconds
    time_index_ratio = len(times) / times[len(times) - 1]
    frequencies_index_ratio = len(
        frequencies) / frequencies[len(frequencies) - 1]
    # estimates tempo
    bpm, beats = librosa.beat.beat_track(time_series, sample_rate)
    first_time = list(librosa.frames_to_time(beats, sr=sample_rate))[0]

    background = cv2.cvtColor(cv2.resize(np.array(Image.open(background) if background is not None else np.random.randint(
        255, 256, (width, height, 3), dtype=np.uint8), dtype=np.uint8), (width, height)), cv2.COLOR_RGB2BGR)
    backgroundEffect = False

    fourcc = VideoWriter_fourcc(*'avc1')
    video = VideoWriter(os.path.join(
        "static", "GENERATED_VIDEOS", filename), fourcc, float(fps), (width, height))

    base_freqs = frequencies
    analyzed_freqs = [frequencies[i]
                      for i in range(0, len(base_freqs), len(base_freqs)//26)]
    while len(analyzed_freqs) > 26:
        analyzed_freqs.pop()

    widgetsSrc = {
        "circularBeatVisualizer": BeatVisualizer,
        "circularBeatVisualizer_default": {
            "color": (0, 0, 0),
            "first_time": first_time,
            "bpm": bpm,
            "x": 0,
            "y": 0,
            "max_effect_strenght": 50,
            "min_height": 100,
            "thickness": 14,
            "decreasing_step": 1,
            "delay_tolerance": 0.1,
        },
        "freqVolumeBars": FreqVisualizerGroup,
        "freqVolumeBars_default": {
            "s_width": width,
            "s_height": height,
            "direction": "down",
            "freqvolumetype": "rectangle",
            "color": (0, 0, 0),
            "bpm": bpm,
            "freqs": analyzed_freqs,
            "x": 0,
            "y": 0,
            "circle_y_gap": 2,
            "max_effect_strenght": 50,
            "min_height": 10,
            "thickness": 10,
            "delay_tolerance": 0.1,
            "width": width//26,
            "max_height": 300,
            "min_decibel": -80,
            "max_decibel": 0,
        },
        "colorWoofer": ColorWoofer,
        "colorWoofer_default": {
            "color": (0, 0, 0),
            "bpm": bpm,
            "freqs": [freq for freq in frequencies][1:6],
            "x": 0,
            "y": 0,
            "thickness": 10,
            "width": 50,
            "height": 100,
            "min_decibel": 300,
            "max_decibel": 10000
        },
    }

    widgets = []
    query = json.loads(widgetstring)["widgets"]

    for widget in query:
        params = {}
        if widget["name"] == "blurryBackground":
            if backgroundEffect:
                continue
            backgroundEffect = True
            background = Background(**{
                "image": background,
                "height": height,
                "width": width,
                "bpm": bpm,
                "first_time": first_time,
                "max_blur_effect_strenght": 5,
            })
            continue
        elif not widget["name"].endswith("_default") and widget["name"] in widgetsSrc.keys():
            toHave = list(widgetsSrc[widget["name"]+"_default"].keys())
        else:
            break
        for param in toHave:
            try:
                params[param] = widget[param]
            except KeyError:
                params[param] = widgetsSrc[widget["name"]+"_default"][param]
        params["name"] = widget["name"]
        widgets.append(widgetsSrc[widget["name"]](**params))

    duration = librosa.get_duration(filename=audio_file)
    freqN = len(frequencies)-1
    timesN = len(times)-1

    def get_decibel(target_time, freq):
        return spectrogram[min(int(freq * frequencies_index_ratio), freqN)][min(int(target_time * time_index_ratio), timesN)]

    def get_ms(i):
        return (i/fps)*1000

    def get_s(i):
        return i/fps

    i = 0
    getTicksLastFrame = 0

    for i in range(ceil(duration*fps)):
        t = get_ms(i)
        ts = get_s(i)

        if backgroundEffect:
            frame = background.get_frame(time=ts)
        else:
            frame = background.copy()

        deltaTime = (t - getTicksLastFrame) / 1000.0
        getTicksLastFrame = t

        for widget in widgets:
            if widget.name == "circularBeatVisualizer":
                widget.draw(dt=deltaTime, time=ts, frame=frame)
            elif widget.name == "freqVolumeBars":
                widget.draw(dt=deltaTime, db=[get_decibel(
                    t/1000, freq) for freq in widget.freqs], frame=frame)
            elif widget.name == "colorWoofer":
                widget.draw(db=[get_decibel(ts, freq)
                                for freq in widget.freqs], frame=frame)
            elif widget.name == "text":
                widget.draw(frame)

        video.write(frame)

    video.release()
    filename2 = video_id + "_audio.mp4"
    videoClip = mp.VideoFileClip(os.path.join(
        "static", "GENERATED_VIDEOS", filename))
    audioClip = mp.AudioFileClip(audio_file)
    final = videoClip.set_audio(audioClip)
    final.write_videofile(os.path.join(
        "static", "GENERATED_VIDEOS", filename2), fps=fps, codec="libx264", audio_codec="libvorbis", bitrate="6000k")
    return filename2

max_stars_repo_name: ldvlpr-dev/AudioVisualizer-Release

max_issues_repo_name: ldvlpr-dev/AudioVisualizer-Release

max_forks_repo_name: ldvlpr-dev/AudioVisualizer-Release


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

import os, sys
import numpy as np

import tensorflow as tf
import util as util

def read_doc_keys(fname):
    keys = set()
    with open(fname) as f:
        for line in f:
            keys.add(line.strip())
    return keys

if __name__ == "__main__":
  config = util.initialize_from_env()
  model = util.get_model(config)
  saver = tf.train.Saver()
  log_dir = config["log_dir"]
  with tf.Session() as session:
    model.restore(session)
    trigger_ids = np.array([[int(x) for x in sys.argv[2:]]]) #[ 1109 13307 12680]
    # Make sure eval mode is True if you want official conll results
    model.evaluate(session, official_stdout=True, eval_mode=True, trigger_token_ids=trigger_ids)

max_stars_repo_name: JieyuZhao/coref

max_issues_repo_name: JieyuZhao/coref

max_forks_repo_name: JieyuZhao/coref


text:
# Copyright 2018 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
""" Decode stimulus for a given stimulus.

s_r_embedding_multiple_retina_refractor --logtostderr --stim_layers='1, 5, 1, 3, 128, 1, 3, 128, 1, 3, 128, 1, 3, 128, 2, 3, 128, 2, 3, 1, 1' --resp_layers='3, 128, 1, 3, 128, 1, 3, 128, 1, 3, 128, 2, 3, 128, 2, 3, 1, 1' --batch_norm=True --save_folder='/cns/in-d/home/bhaishahster/end_to_end_refrac_2' --save_suffix='_stim-resp_wn_nsem' --is_test=1 --taskid=24
"""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from absl import app
from absl import flags
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

import os.path
import tensorflow as tf
import tensorflow.contrib.slim as slim
from tensorflow.python.platform import gfile
import numpy as np, h5py,numpy
import scipy.io as sio
from numpy.random import RandomState
import pickle

FLAGS = tf.app.flags.FLAGS

def predict_stimulus(resp_batch, sr_graph, resp): # responses[iretina]
  '''Predict responses'''

  resp_iters  = []

  n_cells = resp['responses'].shape[1]
  t_len = resp_batch.shape[0]

  stim_np = np.zeros((t_len, 80, 40, 30)).astype(np.float32)
  step_sz = 0.001
  eps = 1e-1
  dist_prev = np.inf
  d_log = []
  # theta_np_log = []

  grad_stim = tf.gradients(sr_graph.d_s_r_pos, sr_graph.stim_tf)[0]

  for iiter in range(200):
    feed_dict ={sr_graph.stim_tf: stim_np.astype(np.float32),
                sr_graph.anchor_model.responses_tf: resp_batch,
                sr_graph.anchor_model.map_cell_grid_tf: resp['map_cell_grid'],
                sr_graph.anchor_model.cell_types_tf: resp['ctype_1hot'],
                sr_graph.anchor_model.mean_fr_tf: resp['mean_firing_rate']}

    dist_np, grad_stim_np = sr_graph.sess.run([sr_graph.d_s_r_pos, grad_stim], feed_dict=feed_dict)

    # update theta
    stim_np = stim_np - step_sz * grad_stim_np

    # theta_np_log += [theta_np]
    if np.sum(np.abs(dist_prev - dist_np)) < eps:
      break
    print(iiter, np.sum(dist_np), np.sum(np.abs(dist_prev - dist_np)))
    dist_prev = dist_np
    d_log += [np.sum(dist_np)]

  return stim_np

def plot_decoding(stimulus_target, stim_decode, n_targets=10):

  batch_size = stimulus_target.shape[0]

  plt.figure()
  plt_tm = np.random.randint(0, batch_size, n_targets)
  n_cols = 5
  for itarget in range(n_targets):
    plt.subplot(n_targets, n_cols, itarget * n_cols + 1)
    plt.imshow(stimulus_target[plt_tm[itarget], :, :, 4].T, cmap='gray', interpolation='nearest')
    if itarget == 0:
      plt.title('Target, t-4')
    plt.axis('off')

    plt.subplot(n_targets, n_cols, itarget * n_cols + 2)
    s_target = stimulus_target[plt_tm[itarget], :, :, 4].T
    s_target_blur =scipy.ndimage.filters.gaussian_filter(s_target, 2)
    plt.imshow(s_target_blur, cmap='gray', interpolation='nearest')
    if itarget == 0:
      plt.title('Smoothened target, t-4')
    plt.axis('off')

    plt.subplot(n_targets, n_cols, itarget * n_cols + 3)
    plt.imshow(stim_decode[plt_tm[itarget], :, :, 4].T, cmap='gray', interpolation='nearest')
    if itarget == 0:
      plt.title('decoded t-4')
    plt.axis('off')

    plt.subplot(n_targets, n_cols, itarget * n_cols + 4)
    plt.imshow(stim_decode[plt_tm[itarget], :, :, 7].T, cmap='gray', interpolation='nearest')
    if itarget == 0:
      plt.title('decoded t-7')
    plt.axis('off')

    plt.subplot(n_targets, n_cols, itarget * n_cols + 5)
    plt.imshow(stim_decode[plt_tm[itarget], :, :, 10].T, cmap='gray', interpolation='nearest')
    if itarget == 0:
      plt.title('decoded t-10')
    plt.axis('off')
  plt.show()

max_stars_repo_name: googlearchive/rgc-models

max_issues_repo_name: google/rgc-models

max_forks_repo_name: google/rgc-models

reason: import numpy,from numpy,import scipy