text:
import collections
import numpy as np
import networkx as nx
import itertools
import logging
from networkx.algorithms.flow import shortest_augmenting_path, edmonds_karp, preflow_push
from networkx.algorithms.connectivity import minimum_st_edge_cut
import time
import graph_tool.flow

from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple, Union

from pychunkedgraph.backend import flatgraph_utils

float_max = np.finfo(np.float32).max

def merge_cross_chunk_edges(edges: Iterable[Sequence[np.uint64]],
                            affs: Sequence[np.uint64],
                            logger: Optional[logging.Logger] = None):
    """ Merges cross chunk edges
    :param edges: n x 2 array of uint64s
    :param affs: float array of length n
    :return:
    """
    # mask for edges that have to be merged
    cross_chunk_edge_mask = np.isinf(affs)

    # graph with edges that have to be merged
    cross_chunk_graph = nx.Graph()
    cross_chunk_graph.add_edges_from(edges[cross_chunk_edge_mask])

    # connected components in this graph will be combined in one component
    ccs = nx.connected_components(cross_chunk_graph)

    # Build mapping
    # For each connected component the smallest node id is chosen to be the
    # representative.
    remapping = {}
    mapping_ks = []
    mapping_vs = []

    for cc in ccs:
        nodes = np.array(list(cc))
        rep_node = np.min(nodes)

        remapping[rep_node] = nodes
        mapping_ks.extend(nodes)
        mapping_vs.extend([rep_node] * len(nodes))

    # Initialize mapping with a each node mapping to itself, then update
    # those edges merged to one across chunk boundaries.
    # u_nodes = np.unique(edges)
    # mapping = dict(zip(u_nodes, u_nodes))
    # mapping.update(dict(zip(mapping_ks, mapping_vs)))
    mapping = dict(zip(mapping_ks, mapping_vs))

    # Vectorize remapping
    mapping_vec = np.vectorize(lambda a : mapping[a] if a in mapping else a)
    remapped_edges = mapping_vec(edges)

    # Remove cross chunk edges
    remapped_edges = remapped_edges[~cross_chunk_edge_mask]
    remapped_affs = affs[~cross_chunk_edge_mask]

    return remapped_edges, remapped_affs, mapping, remapping

def merge_cross_chunk_edges_graph_tool(edges: Iterable[Sequence[np.uint64]],
                                       affs: Sequence[np.uint64],
                                       logger: Optional[logging.Logger] = None):
    """ Merges cross chunk edges
    :param edges: n x 2 array of uint64s
    :param affs: float array of length n
    :return:
    """

    # mask for edges that have to be merged
    cross_chunk_edge_mask = np.isinf(affs)

    # graph with edges that have to be merged
    graph, _, _, unique_ids = flatgraph_utils.build_gt_graph(
        edges[cross_chunk_edge_mask], make_directed=True)

    # connected components in this graph will be combined in one component
    ccs = flatgraph_utils.connected_components(graph)

    remapping = {}
    mapping = np.array([], dtype=np.uint64).reshape(-1, 2)

    for cc in ccs:
        print(cc)
        nodes = unique_ids[cc]
        rep_node = np.min(nodes)

        remapping[rep_node] = nodes

        rep_nodes = np.ones(len(nodes), dtype=np.uint64).reshape(-1, 1) * rep_node
        m = np.concatenate([nodes.reshape(-1, 1), rep_nodes], axis=1)

        mapping = np.concatenate([mapping, m], axis=0)

    u_nodes = np.unique(edges)
    u_unmapped_nodes = u_nodes[~np.in1d(u_nodes, mapping)]

    unmapped_mapping = np.concatenate([u_unmapped_nodes.reshape(-1, 1),
                                       u_unmapped_nodes.reshape(-1, 1)], axis=1)
    mapping = np.concatenate([mapping, unmapped_mapping], axis=0)

    sort_idx = np.argsort(mapping[:, 0])
    idx = np.searchsorted(mapping[:, 0], edges, sorter=sort_idx)
    remapped_edges = np.asarray(mapping[:, 1])[sort_idx][idx]

    remapped_edges = remapped_edges[~cross_chunk_edge_mask]
    remapped_affs = affs[~cross_chunk_edge_mask]

    return remapped_edges, remapped_affs, mapping, remapping

def mincut_nx(edges: Iterable[Sequence[np.uint64]], affs: Sequence[np.uint64],
              sources: Sequence[np.uint64], sinks: Sequence[np.uint64],
              logger: Optional[logging.Logger] = None) -> np.ndarray:
    """ Computes the min cut on a local graph
    :param edges: n x 2 array of uint64s
    :param affs: float array of length n
    :param sources: uint64
    :param sinks: uint64
    :return: m x 2 array of uint64s
        edges that should be removed
    """

    time_start = time.time()

    original_edges = edges.copy()

    edges, affs, mapping, remapping = merge_cross_chunk_edges(edges.copy(),
                                                              affs.copy())
    mapping_vec = np.vectorize(lambda a: mapping[a] if a in mapping else a)

    if len(edges) == 0:
        return []

    if len(mapping) > 0:
        assert np.unique(list(mapping.keys()), return_counts=True)[1].max() == 1

    remapped_sinks = mapping_vec(sinks)
    remapped_sources = mapping_vec(sources)

    sinks = remapped_sinks
    sources = remapped_sources

    sink_connections = np.array(list(itertools.product(sinks, sinks)))
    source_connections = np.array(list(itertools.product(sources, sources)))

    weighted_graph = nx.Graph()
    weighted_graph.add_edges_from(edges)
    weighted_graph.add_edges_from(sink_connections)
    weighted_graph.add_edges_from(source_connections)

    for i_edge, edge in enumerate(edges):
        weighted_graph[edge[0]][edge[1]]['capacity'] = affs[i_edge]
        weighted_graph[edge[1]][edge[0]]['capacity'] = affs[i_edge]

    # Add infinity edges for multicut
    for sink_i in sinks:
        for sink_j in sinks:
            weighted_graph[sink_i][sink_j]['capacity'] = float_max

    for source_i in sources:
        for source_j in sources:
            weighted_graph[source_i][source_j]['capacity'] = float_max

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Graph creation: %.2fms" % (dt * 1000))
    time_start = time.time()

    ccs = list(nx.connected_components(weighted_graph))
    for cc in ccs:
        cc_list = list(cc)

        # If connected component contains no sources and/or no sinks,
        # remove its nodes from the mincut computation
        if not np.any(np.in1d(sources, cc_list)) or \
                not np.any(np.in1d(sinks, cc_list)):
            weighted_graph.remove_nodes_from(cc)

    r_flow = edmonds_karp(weighted_graph, sinks[0], sources[0])
    cutset = minimum_st_edge_cut(weighted_graph, sources[0], sinks[0],
                                 residual=r_flow)

    # cutset = nx.minimum_edge_cut(weighted_graph, sources[0], sinks[0], flow_func=edmonds_karp)

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Mincut comp: %.2fms" % (dt * 1000))

    if cutset is None:
        return []

    time_start = time.time()

    edge_cut = list(list(cutset))

    weighted_graph.remove_edges_from(edge_cut)
    ccs = list(nx.connected_components(weighted_graph))

    # assert len(ccs) == 2

    for cc in ccs:
        cc_list = list(cc)
        if logger is not None:
            logger.debug("CC size = %d" % len(cc_list))

        if np.any(np.in1d(sources, cc_list)):
            assert np.all(np.in1d(sources, cc_list))
            assert ~np.any(np.in1d(sinks, cc_list))

        if np.any(np.in1d(sinks, cc_list)):
            assert np.all(np.in1d(sinks, cc_list))
            assert ~np.any(np.in1d(sources, cc_list))

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Splitting local graph: %.2fms" % (dt * 1000))

    remapped_cutset = []
    for cut in cutset:
        if cut[0] in remapping:
            pre_cut = remapping[cut[0]]
        else:
            pre_cut = [cut[0]]

        if cut[1] in remapping:
            post_cut = remapping[cut[1]]
        else:
            post_cut = [cut[1]]

        remapped_cutset.extend(list(itertools.product(pre_cut, post_cut)))
        remapped_cutset.extend(list(itertools.product(post_cut, pre_cut)))

    remapped_cutset = np.array(remapped_cutset, dtype=np.uint64)

    remapped_cutset_flattened_view = remapped_cutset.view(dtype='u8,u8')
    edges_flattened_view = original_edges.view(dtype='u8,u8')

    cutset_mask = np.in1d(remapped_cutset_flattened_view, edges_flattened_view)

    return remapped_cutset[cutset_mask]

def mincut_graph_tool(edges: Iterable[Sequence[np.uint64]],
                      affs: Sequence[np.uint64],
                      sources: Sequence[np.uint64],
                      sinks: Sequence[np.uint64],
                      logger: Optional[logging.Logger] = None) -> np.ndarray:
    """ Computes the min cut on a local graph
    :param edges: n x 2 array of uint64s
    :param affs: float array of length n
    :param sources: uint64
    :param sinks: uint64
    :return: m x 2 array of uint64s
        edges that should be removed
    """
    time_start = time.time()

    original_edges = edges

    # Stitch supervoxels across chunk boundaries and represent those that are
    # connected with a cross chunk edge with a single id. This may cause id
    # changes among sinks and sources that need to be taken care of.
    edges, affs, mapping, remapping = merge_cross_chunk_edges(edges.copy(),
                                                              affs.copy())

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Cross edge merging: %.2fms" % (dt * 1000))
    time_start = time.time()

    mapping_vec = np.vectorize(lambda a: mapping[a] if a in mapping else a)

    if len(edges) == 0:
        return []

    if len(mapping) > 0:
        assert np.unique(list(mapping.keys()), return_counts=True)[1].max() == 1

    remapped_sinks = mapping_vec(sinks)
    remapped_sources = mapping_vec(sources)

    sinks = remapped_sinks
    sources = remapped_sources

    # Assemble edges: Edges after remapping combined with edges between sinks
    # and sources
    sink_edges = list(itertools.product(sinks, sinks))
    source_edges = list(itertools.product(sources, sources))

    comb_edges = np.concatenate([edges, sink_edges, source_edges])

    comb_affs = np.concatenate([affs, [float_max, ] *
                                (len(sink_edges) + len(source_edges))])

    # To make things easier for everyone involved, we map the ids to
    # [0, ..., len(unique_ids) - 1]
    # Generate weighted graph with graph_tool
    weighted_graph, cap, gt_edges, unique_ids = \
        flatgraph_utils.build_gt_graph(comb_edges, comb_affs,
                                       make_directed=True)

    sink_graph_ids = np.where(np.in1d(unique_ids, sinks))[0]
    source_graph_ids = np.where(np.in1d(unique_ids, sources))[0]

    if logger is not None:
        logger.debug(f"{sinks}, {sink_graph_ids}")
        logger.debug(f"{sources}, {source_graph_ids}")

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Graph creation: %.2fms" % (dt * 1000))
    time_start = time.time()

    # # Get rid of connected components that are not involved in the local
    # # mincut
    # cc_prop, ns = graph_tool.topology.label_components(weighted_graph)
    #
    # if len(ns) > 1:
    #     cc_labels = cc_prop.get_array()
    #
    #     for i_cc in range(len(ns)):
    #         cc_list = np.where(cc_labels == i_cc)[0]
    #
    #         # If connected component contains no sources and/or no sinks,
    #         # remove its nodes from the mincut computation
    #         if not np.any(np.in1d(source_graph_ids, cc_list)) or \
    #                 not np.any(np.in1d(sink_graph_ids, cc_list)):
    #             weighted_graph.delete_vertices(cc) # wrong

    # Compute mincut
    src, tgt = weighted_graph.vertex(source_graph_ids[0]), \
               weighted_graph.vertex(sink_graph_ids[0])

    res = graph_tool.flow.boykov_kolmogorov_max_flow(weighted_graph,
                                                     src, tgt, cap)

    part = graph_tool.flow.min_st_cut(weighted_graph, src, cap, res)

    labeled_edges = part.a[gt_edges]
    cut_edge_set = gt_edges[labeled_edges[:, 0] != labeled_edges[:, 1]]

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Mincut comp: %.2fms" % (dt * 1000))
    time_start = time.time()

    if len(cut_edge_set) == 0:
        return []

    time_start = time.time()

    # Make sure we did not do something wrong: Check if sinks and sources are
    # among each other and not in different sets
    for i_cc in np.unique(part.a):
        # Make sure to read real ids and not graph ids
        cc_list = unique_ids[np.array(np.where(part.a == i_cc)[0],
                                      dtype=np.int)]

        # if logger is not None:
        #     logger.debug("CC size = %d" % len(cc_list))

        if np.any(np.in1d(sources, cc_list)):
            assert np.all(np.in1d(sources, cc_list))
            assert ~np.any(np.in1d(sinks, cc_list))

        if np.any(np.in1d(sinks, cc_list)):
            assert np.all(np.in1d(sinks, cc_list))
            assert ~np.any(np.in1d(sources, cc_list))

    dt = time.time() - time_start
    if logger is not None:
        logger.debug("Verifying local graph: %.2fms" % (dt * 1000))

    # Extract original ids
    # This has potential to be optimized
    remapped_cutset = []
    for s, t in flatgraph_utils.remap_ids_from_graph(cut_edge_set, unique_ids):

        if s in remapping:
            s = remapping[s]
        else:
            s = [s]

        if t in remapping:
            t = remapping[t]
        else:
            t = [t]

        remapped_cutset.extend(list(itertools.product(s, t)))
        remapped_cutset.extend(list(itertools.product(t, s)))

    remapped_cutset = np.array(remapped_cutset, dtype=np.uint64)

    remapped_cutset_flattened_view = remapped_cutset.view(dtype='u8,u8')
    edges_flattened_view = original_edges.view(dtype='u8,u8')

    cutset_mask = np.in1d(remapped_cutset_flattened_view, edges_flattened_view)

    return remapped_cutset[cutset_mask]

def mincut(edges: Iterable[Sequence[np.uint64]],
           affs: Sequence[np.uint64],
           sources: Sequence[np.uint64],
           sinks: Sequence[np.uint64],
           logger: Optional[logging.Logger] = None) -> np.ndarray:
    """ Computes the min cut on a local graph
    :param edges: n x 2 array of uint64s
    :param affs: float array of length n
    :param sources: uint64
    :param sinks: uint64
    :return: m x 2 array of uint64s
        edges that should be removed
    """

    return mincut_graph_tool(edges=edges, affs=affs, sources=sources,
                             sinks=sinks, logger=logger)

max_stars_repo_name: perlman/PyChunkedGraph

max_issues_repo_name: perlman/PyChunkedGraph

max_forks_repo_name: perlman/PyChunkedGraph

reason: import numpy,import networkx,from networkx


text:
### Voxel-wise calculation of flip angle deviation with the double angle method on spin echo data
#
# Author: Francesco Grussu, University College London
#		    CDSQuaMRI Project
#		   <f.grussu@ucl.ac.uk> <francegrussu@gmail.com>
#
# Code released under BSD Two-Clause license
#
# Copyright (c) 2019 University College London.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# The views and conclusions contained in the software and documentation are those
# of the authors and should not be interpreted as representing official policies,
# either expressed or implied, of the FreeBSD Project.

### Load useful modules
import argparse, os, sys
import nibabel as nib
import numpy as np
import warnings

def DASE(*argv):
	''' Calculate B1 map from two spin echo images using the "double angle" method

	    INTERFACES
	    DASE(se1_nifti,se2_nifti,theta,se_type,b1_output)
	    DASE(se1_nifti,se2_nifti,theta,se_type,b1_output,mask_nifti)

	    PARAMETERS
	    - se1_nifti: path of a Nifti file storing the spin echo image acquired with excitation flip angle theta
	    - se2_nifti: path of a Nifti file storing the spin echo image acquired with excitation flip angle 2 x theta
	    - theta: path of a text file storing the flip angle (in degrees)
	    - se_type: a string specifying the type of spin echo sequence used (set it to "theta_2theta" for a spin echo
                       sequence like "theta-TE/2-2theta-TE/2-acquisition", as in Insko et al, JMR 1993;
                       or set it to "theta_180" for a spin echo sequence like "theta-TE/2-180-TE/2-acquisition",
                       as in Stollberger et al, MRM 1996)
	    - b1_output:  path of the Nifti file that will store the 3D B1 map (saved as a
			  double-precision floating point image FLOAT64); such an output map is such that

				AF = B1 x NF

			  above, NF is the nominal flip angle prescribed to the scanner;
			  AF is the actual flip angle (true flip angle played out by the scanner in each voxel);
			  B1 is the voxel-wise B1 map (deviation between actual and nominal flip angle)

	    - mask_nifti: path of a Nifti file storting a mask (B1 map will be calculated only where
			  mask_nifti equals 1; 0 will be set in the B1 output map otherwise)

	    Dependencies (Python packages): nibabel, numpy (other than standard library)

	    References:
		        "Mapping of the radiofrequency field",
		        Insko EK and Bolinger L, Journal of Magnetic Resonance Series A, 103:82-85 (1993)

			"Imaging of the active B1 field in vivo",
                        Stollberger R and Wach P, Magnetic Resonance in Medicine, 35(2): 246-251 (1996)

	    Author: Francesco Grussu, University College London
		    CDSQuaMRI Project
		   <f.grussu@ucl.ac.uk> <francegrussu@gmail.com>'''

	### Get input parametrs
	Nargv = len(argv)
	scan1_nifti = argv[0]
	scan2_nifti = argv[1]
	theta_file = argv[2]
	se_type = argv[3]
	b1_output = argv[4]

	### Load scan with original flip angle and check for consistency
	print('    ... loading input data')
	try:
		scan1_obj = nib.load(scan1_nifti)
	except:
		print('')
		print('ERROR: the 3D spin echo NIFTI {} does not exist or is not in NIFTI format. Exiting with 1.'.format(scan1_nifti))
		print('')
		sys.exit(1)
	scan1_data = scan1_obj.get_fdata()
	imgsize = scan1_data.shape
	imgsize = np.array(imgsize)
	if imgsize.size!=3:
		print('')
		print('ERROR: the 3D spin echo NIFTI {} is not a 3D NIFTI. Exiting with 1.'.format(scan1_nifti))
		print('')
		sys.exit(1)

	### Load scan acquired with doubled flip angle and check for consistency
	try:
		scan2_obj = nib.load(scan2_nifti)
	except:
		print('')
		print('ERROR: the 3D spin echo NIFTI {} does not exist or is not in NIFTI format. Exiting with 1.'.format(scan2_nifti))
		print('')
		sys.exit(1)
	scan2_data = scan2_obj.get_fdata()
	scan1_header = scan1_obj.header
	scan1_affine = scan1_header.get_best_affine()
	scan1_dims = scan1_obj.shape
	scan2_header = scan2_obj.header
	scan2_affine = scan2_header.get_best_affine()
	scan2_dims = scan2_obj.shape
	scan2_size = scan2_data.shape
	scan2_size = np.array(scan2_size)
	if scan2_size.size!=3:
		print('')
		print('ERROR: the 3D spin echo NIFTI {} is not a 3D NIFTI. Exiting with 1.'.format(scan2_nifti))
		print('')
		sys.exit(1)
	elif ( (np.sum(scan1_affine==scan2_affine)!=16) or (scan1_dims[0]!=scan2_dims[0]) or (scan1_dims[1]!=scan2_dims[1]) or (scan1_dims[2]!=scan2_dims[2]) ):
		print('')
		print('ERROR: the geometries of the the two spin echo files {} and {} do not match. Exiting with 1.'.format(scan1_nifti,scan2_nifti))
		print('')
		sys.exit(1)

	### Load the flip angle file

	# Make sure data exists and makes sense
	try:
		thetavalue = np.loadtxt(theta_file)
		thetavalue = np.array(thetavalue,'float64')
		thetavalue_size = thetavalue.size
	except:
		print('')
		print('ERROR: the flip angle file {} does not exist or is not a numeric text file. Exiting with 1.'.format(theta_file))
		print('')
		sys.exit(1)

	# Make sure we get only one scalar value
	if thetavalue_size!=1:
		print('')
		print('ERROR: the flip angle file {} contains more than one entry or is empty. Exiting with 1.'.format(theta_file))
		print('')
		sys.exit(1)
	else:
		theta_nominal = thetavalue

	## Select the sequence type
	if se_type=='theta_2theta':
		isseq_theta180 = 0

	elif se_type=='theta_180':
		isseq_theta180 = 1

	else:
		print('')
		print('ERROR: the sequence type "{}" is not recognised. Choose either "theta_2theta" or "theta_180". Exiting with 1.'.format(se_type))
		print('')
		sys.exit(1)

	### Deal with optional arguments: mask
	if Nargv==6:
		got_mask = True
		mask_nifti = argv[4]
		try:
			mask_obj = nib.load(mask_nifti)
		except:
			print('')
			print('ERROR: the mask file {} does not exist or is not in NIFTI format. Exiting with 1.'.format(mask_nifti))
			print('')
			sys.exit(1)
		# Make sure that the mask geometry matches that of the other files
		mask_data = mask_obj.get_fdata()
		mask_size = mask_data.shape
		mask_size = np.array(mask_size)
		mask_header = mask_obj.header
		mask_affine = mask_header.get_best_affine()
		mask_dims = mask_obj.shape
		if mask_size.size!=3:
			print('')
			print('WARNING: the mask file {} is not a 3D NIFTI file. Ignoring mask...'.format(mask_nifti))
			print('')
			mask_data = np.ones(imgsize[0:3],'float64')
		elif ( (np.sum(scan1_affine==mask_affine)!=16) or (scan1_dims[0]!=mask_dims[0]) or (scan1_dims[1]!=mask_dims[1]) or (scan1_dims[2]!=mask_dims[2]) ):
			print('')
			print('WARNING: the geometry of the the mask file {} does not match that of the spin echo data. Ignoring mask...'.format(mask_nifti))
			print('')
			mask_data = np.ones(imgsize[0:3],'float64')
		else:
			mask_data = np.array(mask_data,'float64')
			# Make sure mask is a binary file
			mask_data[mask_data>1] = 1
			mask_data[mask_data<0] = 0

	else:
		got_mask = False

	### Calculate B1 map
	print('    ... calculating B1 map')
	scan1_data = np.array(scan1_data,'float64')
	scan2_data = np.array(scan2_data,'float64')
	warnings.filterwarnings('ignore')    # Ignore warnings - these are going to happen in the background for sure
	Sratio = scan1_data/scan2_data
	if isseq_theta180==0:
		theta_measured = np.arccos( 0.5*( (1.0/Sratio)**(1/3) ) )  # Sequence is theta - TE/2 - 2 x theta - TE/2, as Insko EK and Bolinger L
	elif isseq_theta180==1:
		theta_measured = np.arccos( 0.5*(1.0/Sratio) ) # Sequence is theta - TE/2 - 180 - TE/2, as Insko EK and Bolinger L
	else:
		print('')
		print('ERROR: the sequence type "{}" is not recognised. Choose either "theta_2theta" or "theta_180". Exiting with 1.'.format(se_type))
		print('')
		sys.exit(1)

	b1_map = theta_measured/(np.pi*theta_nominal/180.0)
	b1_map[np.isnan(b1_map)] = 0.0   # Remove NaNs
	b1_map[np.isinf(b1_map)] = 0.0   # Remove Inf
	print('    ... substituting any nan and inf values with 0')
	# Use mask if required
	if got_mask==True:
		b1_map[mask_data==0]=0

	### Save output
	print('    ... saving output file')
	buffer_header = scan1_obj.header
	buffer_header.set_data_dtype('float64')   # Make sure we save quantitative map as a float64
	b1_obj = nib.Nifti1Image(b1_map,scan1_obj.affine,buffer_header)
	nib.save(b1_obj, b1_output)

	### Done
	print('')

# Run the module as a script when required
if __name__ == "__main__":

	### Print help and parse arguments
	parser = argparse.ArgumentParser(description='Voxel-wise calculation of B1 map (flip angle deviation) using the "double angle" method based on two spin echo acquisitions. Dependencies (Python packages): nibabel, numpy (other than standard library). References: "Mapping of the radiofrequency field", Insko EK and Bolinger L, Journal of Magnetic Resonance Series A, 103:82-85 (1993); "Imaging of the active B1 field in vivo", Stollberger R and Wach P, Magnetic Resonance in Medicine, 35(2): 246-251 (1996). Author: Francesco Grussu, University College London, CDSQuaMRI Project. Email: <francegrussu@gmail.com> <f.grussu@ucl.ac.uk>.')
	parser.add_argument('nifti1_file', help='3D Nifti file storing the spin echo image acquired with flip angle theta')
	parser.add_argument('nifti2_file', help='3D Nifti file storing the spin echo image acquired with flip angle 2 x theta')
	parser.add_argument('flipangle_file', help='path of a text file storting the flip angle theta (in degrees)')
	parser.add_argument('out_file', help='3D Nifti file that will store the B1 map as double-precision floating point (FLOAT64), such that the actual flip angle (AF) is related to a nominal flip angle (NF) as AF = B1 x NF')
	parser.add_argument('--refocus', metavar='<type>', default="theta_180", help='type of refocussing pulse of the spin echo sequence (set it to "theta_2theta" for a spin echo sequence like "theta-TE/2-2theta-TE/2-acquisition", as in Insko et al, JMR 1993; or set it to "theta_180" for a spin echo sequence like "theta-TE/2-180-TE/2-acquisition", as in Stollberger et al, MRM 1996; default "theta_180")')
	parser.add_argument('--mask', metavar='<file>', help='mask in Nifti format where 1 flags voxels where the flip angle map is required, 0 where is not')
	args = parser.parse_args()

	### Get input arguments
	spinecho1file = args.nifti1_file
	spinecho2file = args.nifti2_file
	fafile = args.flipangle_file
	outfile = args.out_file
	refocustype = args.refocus
	maskfile = args.mask

	### Deal with optional arguments
	if isinstance(maskfile, str)==1:
	    # A mask for flip angle calculation has been provided
	    maskrequest = True
	else:
	    # A mask for flip angle calculation has not been provided
	    maskrequest = False

	print('')
	print('********************************************************************')
	print('     B1 mapping with the double angle method on spin echo data      ')
	print('********************************************************************')
	print('')
	print('Called on 3D Nifti files: {} (spin echo images for flip angle theta) and {} (spin echo images for flip angle 2 x theta)'.format(spinecho1file,spinecho2file))
	print('Flip angle file: {}'.format(fafile))
	print('Refocusing pulse thype: {}'.format(refocustype))
	print('Output B1 map file: {}'.format(outfile))
	print('')

	if (maskrequest==False):
		DASE(spinecho1file, spinecho2file, fafile, refocustype, outfile)
	else:
		DASE(spinecho1file, spinecho2file, fafile, refocustype, outfile, maskfile)

	### Done
	print('Processing completed.')
	print('')
	sys.exit(0)

ext: py

max_stars_repo_name: fragrussu/MyRelax

max_issues_repo_name: fragrussu/MyRelax

max_forks_repo_name: fragrussu/MyRelax

reason: import numpy


text:
import json
import os
import argparse
from functools import partial
from multiprocessing.pool import Pool
from os import cpu_count

import cv2
from cv2 import fillPoly
from shapely import wkt
import numpy as np
from shapely.geometry import mapping
from tqdm import tqdm

def generate_localization_polygon(json_path, out_dir):
    with open(json_path, "r") as f:
        annotations = json.load(f)
    h = annotations["metadata"]["height"]
    w = annotations["metadata"]["width"]
    mask_img = np.zeros((h, w), np.uint8)
    out_filename = os.path.splitext(os.path.basename(json_path))[0] + ".png"
    for feat in annotations['features']['xy']:
        feat_shape = wkt.loads(feat['wkt'])
        coords = list(mapping(feat_shape)['coordinates'][0])
        fillPoly(mask_img, [np.array(coords, np.int32)], (255))
    cv2.imwrite(os.path.join(out_dir, out_filename), mask_img)

def generate_damage_polygon(json_path, out_dir):
    with open(json_path, "r") as f:
        annotations = json.load(f)

    h = annotations["metadata"]["height"]
    w = annotations["metadata"]["width"]
    mask_img = np.zeros((h, w), np.uint8)

    damage_dict = {
        "no-damage": 1,
        "minor-damage": 2,
        "major-damage": 3,
        "destroyed": 4,
        "un-classified": 255
    }
    out_filename = os.path.splitext(os.path.basename(json_path))[0] + ".png"
    for feat in annotations['features']['xy']:
        feat_shape = wkt.loads(feat['wkt'])
        coords = list(mapping(feat_shape)['coordinates'][0])
        fillPoly(mask_img, [np.array(coords, np.int32)], damage_dict[feat['properties']['subtype']])
    cv2.imwrite(os.path.join(out_dir, out_filename), mask_img)

if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument('--input',
                        default="/home/selim/datasets/xview/train/",
                        help='Path to parent dataset directory "xBD"')
    args = parser.parse_args()
    out_dir = os.path.join(args.input, "masks")
    in_dir = os.path.join(args.input, "labels")
    pre_images = [os.path.join(in_dir, f) for f in os.listdir(in_dir) if '_pre_' in f]
    post_images = [os.path.join(in_dir, f) for f in os.listdir(in_dir) if '_post_' in f]

    pool = Pool(cpu_count())
    with tqdm(total=len(pre_images)) as pbar:
        for i, v in enumerate(pool.imap_unordered(partial(generate_localization_polygon, out_dir=out_dir), pre_images)):
            pbar.update()
    with tqdm(total=len(post_images)) as pbar:
        for i, v in enumerate(pool.imap_unordered(partial(generate_damage_polygon, out_dir=out_dir), post_images)):
            pbar.update()


text:
import pickle
import collections
import numpy
import redis

from threepio import logger
from django.db.models import (
        Avg, ExpressionWrapper,
        F, Q, fields)
from django.utils import timezone
from dateutil import rrule
from core.models import (
    Instance, InstanceStatusHistory
)

METRICS_CACHE_DURATION = 4*24*60*60  # 4 days (persist over the weekend)

def get_application_metrics(application, now_time=None, read_only=False):
    """
    Skip image metrics on end-dated applications
    Otherwise look through the cache to find application metrics
    """
    metrics = {}
    if application.end_date:
        return metrics
    metrics = _get_application_metrics(application, interval=rrule.DAILY, now_time=now_time, read_only=read_only)
    metrics = _get_application_metrics(application, interval=rrule.WEEKLY, now_time=now_time, read_only=read_only)
    metrics = _get_application_metrics(application, interval=rrule.MONTHLY, now_time=now_time, read_only=read_only)
    return metrics

def _get_application_metrics(application, interval=rrule.MONTHLY, day_limit=120,
                             now_time=None, force=False, read_only=False):
    metrics = collections.OrderedDict()
    if not interval:
        interval = rrule.MONTHLY
    redis_cache = redis.StrictRedis()
    key = "metrics-application-%s-interval-%s-limited-to-%s" % (application.id, rrule.FREQNAMES[interval], day_limit)
    try:
        if redis_cache.exists(key) and not force:
            pickled_object = redis_cache.get(key)
            metrics = pickle.loads(pickled_object)
        elif not read_only:
            metrics = calculate_application_metrics(application, interval, day_limit, now_time=now_time)
            pickled_object = pickle.dumps(metrics)
            redis_cache.set(key, pickled_object)
            redis_cache.expire(key, METRICS_CACHE_DURATION)
    except:
        logger.exception("Unexpected errror in application metrics")
    return metrics

def calculate_application_metrics(application, interval=rrule.MONTHLY,
                                  day_limit=None, sum_datapoints=True, now_time=None):
    """
    From start_date of Application to now/End-date of application
      - Create a timeseries by splitting by 'interval'
      - Query for metrics datapoints
      - Return the timeseries + datapoints
    """
    if not now_time:
        now_time = timezone.now()
    end_date = application.end_date or now_time
    start_date = end_date - timezone.timedelta(days=day_limit)
    timeseries = _generate_time_series(start_date, end_date, interval)
    all_instance_ids = application.versions.values_list('machines__instance_source__instances', flat=True)
    application_metrics = collections.OrderedDict()
    for idx, ts in enumerate(timeseries):
        interval_start = ts
        interval_key = interval_start.strftime("%x %X")
        if sum_datapoints:
            interval_start = start_date
        if idx == len(timeseries)-1:
            interval_end = end_date
        else:
            interval_end = timeseries[idx+1]
        all_instances = Instance.objects\
            .filter(id__in=all_instance_ids)\
            .filter(start_date__gt=interval_start, start_date__lt=interval_end)
        all_histories = InstanceStatusHistory.objects.filter(
            instance__in=all_instances)
        per_interval_metrics = calculate_instance_metrics_for_interval(
            all_instances, all_histories, interval_start, interval_end)
        application_metrics[interval_key] = per_interval_metrics
    return application_metrics

# Alternative calculation method, drilling down per-version rather than per-application.
def calculate_detailed_application_metrics(application, interval=rrule.MONTHLY):
    """
    From start_date of Application to now/End-date of application
      - Create a timeseries by splitting by 'interval'
      - Query for metrics datapoints
      - Return the timeseries + datapoints
    """
    now_time = timezone.now()
    the_beginning = application.start_date
    the_end = application.end_date or now_time
    application_metrics = {}
    for version in application.versions.all():
        per_version_metrics = calculate_application_version_metrics(
            version, the_beginning, the_end, interval=interval)
        application_metrics[version.name] = per_version_metrics
    return application_metrics

def calculate_application_version_metrics(version, start_date, end_date, interval=rrule.MONTHLY):
    timeseries = _generate_time_series(start_date, end_date, interval)
    per_version_metrics = collections.OrderedDict()
    for idx, ts in enumerate(timeseries):
        interval_start = ts
        if idx == len(timeseries)-1:
            interval_end = end_date
        else:
            interval_end = timeseries[idx+1]
        interval_key = interval_start.strftime("%x %X")
        interval_metrics = calculate_application_version_metrics_for_interval(version, interval_start, interval_end)
        per_version_metrics[interval_key] = interval_metrics
    return per_version_metrics

def _generate_time_series(start_date, end_date, interval, limit=None):
    time_series = list(rrule.rrule(interval, dtstart=start_date, until=end_date))
    if not limit:
        return time_series
    return time_series[-1*limit:]

def calculate_instance_metrics_for_interval(all_instances, all_histories, start_date, end_date):
    total_count = all_instances.count()
    active_count = all_instances.filter(Q(instancestatushistory__status__name='active')).distinct().count()
    # build_stats = _stats_for_instance_history(all_histories, 'build', end_date)
    # networking_stats = _stats_for_instance_history(all_histories, 'networking', end_date)
    # deploying_stats = _stats_for_instance_history(all_histories, 'deploying', end_date)
    time_specific_metrics = {
            "active": active_count,
            "total": total_count,
            # "build": build_stats,
            # "networking": networking_stats,
            # "deploying": deploying_stats,
        }
    return time_specific_metrics

def calculate_application_version_metrics_for_interval(application_version, interval_start, interval_end):
    provider_metrics = {}
    for pm in application_version.machines.all():
        all_instances = pm.instance_source.instances.all()
        all_instances = all_instances.filter(start_date__gt=interval_start, start_date__lt=interval_end)
        all_histories = InstanceStatusHistory.objects.filter(
            instance__in=all_instances)
        interval_metrics = calculate_instance_metrics_for_interval(
            all_instances, all_histories, interval_start, interval_end)
        provider_metrics[pm.instance_source.provider.location] = interval_metrics
    return provider_metrics

def get_historical_average(instance_history_list, status_name, limit=None):
    filtered_history_list = instance_history_list.filter(status__name=status_name).distinct()
    tmp_duration = ExpressionWrapper(F('end_date') - F('start_date'), output_field=fields.DurationField())
    tmp_average_time = filtered_history_list.annotate(
            tmp_duration=tmp_duration).aggregate(Avg('tmp_duration'))['tmp_duration__avg']
    return tmp_average_time

def _stats_for_instance_history(histories, status_name, end_date):
    filtered_histories = histories.filter(status__name=status_name)
    status_in_secs = map(lambda x: (x.force_end_date(end_date) - x.start_date).total_seconds(), filtered_histories)
    median_status_time = numpy.median(status_in_secs) if status_in_secs != 0 else 0
    avg_status_time = numpy.mean(status_in_secs) if status_in_secs != 0 else 0
    std_dev_status_time = numpy.std(status_in_secs) if status_in_secs != 0 else 0
    status_stats = {
        'average': avg_status_time,
        'std_deviation': std_dev_status_time,
        'median': median_status_time,
    }
    return status_stats


text:
from ch02.main import records
from ch02.main import print

from pandas import *
from matplotlib import pyplot as plt

import numpy as np

frame = DataFrame(records)
print(frame['tz'].value_counts()[:10])

clean_tz = frame['tz'].fillna('Missing')
clean_tz[clean_tz == ''] = 'Unknown'

tz_counts = clean_tz.value_counts()

print(tz_counts[:10])
tz_counts[:10].plot(kind = 'barh', rot = 0)

cframe = frame[frame.a.notnull()]
operating_system = np.where(cframe['a'].str.contains('Windows'), 'Windows', 'Not Windows')
print(operating_system[:5])

by_tz_os = cframe.groupby(['tz', operating_system])
print(by_tz_os)
agg_counts = by_tz_os.size().unstack().fillna(0)
print(agg_counts[:10])

indexer = agg_counts.sum(1).argsort()
print(indexer[:10])
count_subset = agg_counts.take(indexer)[-10:]
print(count_subset)
count_subset.plot(kind = 'barh', stacked = True)

normed_subset = count_subset.div(count_subset.sum(1), axis = 0)
normed_subset.plot(kind = 'barh', stacked = True)

plt.show()

max_stars_repo_name: QianYuXiang/UsePythonForDataAnalysis

max_issues_repo_name: QianYuXiang/UsePythonForDataAnalysis

max_forks_repo_name: QianYuXiang/UsePythonForDataAnalysis


text:
import os
import time
import sys
import math
import cmath
import glob
import numpy as np

# donde cada uno de los puntajes Spropiedad son funciones tipo Heavside
# que asignan 1 si el valor de tal propiedad cumple las propiedades de
# Lipinski, Ghose y Veber o 0 en caso contrario:

def Heavside (up, down, value):
    boolean = 0
    #
    if((value > down ) and (value < up)):
        boolean = 1

    return boolean

def minmax (array,title):
    max_=np.amax(array)
    min_=np.amin(array)
    desv=np.std(array)
    avg=np.average(array)
    perc=np.percentile(array, 90)
    print (title)
    percQ1=np.percentile(array, 25)
    percQ3=np.percentile(array, 75)
    print("Max: "+str(max_))
    print("Q3: "+str(percQ3))
    print("Avg: "+str(avg))
    print("Q1: "+str(percQ1))
    print("Min: "+str(min_))
    print("DS : "+str(desv))
    print("90perc: "+str(perc)+"\n")
    return avg

def sumalista(listaNumeros):
    laSuma = 0
    for i in listaNumeros:
        laSuma = laSuma + i
    return laSuma

if __name__=='__main__':
    # Time
    t1=time.time()

    totalfiles = []

    list_of_files = glob.glob('./QSAR-2D*.txt.csv')
    FO = open("01DataBase_lin_test.csv", 'w')
    FO.write("MolWeight\tMolLogP\tHbondA\tHbondD\tPolarSA\tRbonds\tMolarRefractivity\tnAtoms\tSynthAcces\tAmesMutagenic\tSeq\tSMILES\n")
    for file_name in list_of_files:
        print ("File: " + str(file_name))
        FI = open(file_name, 'r')
        for line in FI:
            line = line.replace('\n', '').replace('\r', '')
            totalfiles.append(line)
            FO.write(line + "\n")
        FI.close()
    FO.close()

    # Lists
    ListMolWeight        =[]
    ListMolLogP          =[]
    ListHbondA           =[]
    ListHbondD           =[]
    ListPolarSA          =[]
    ListRbonds           =[]
    ListMolarRefractivity=[]
    ListnAtoms           =[]
    #
    ListSynthAcces       =[]
    ListAmesMutagenic    =[]
    #
    ListScoreADME        =[]
    ListScoreTotal       =[]
    #
    ListSeq              =[]
    ##print("MolWeight\tMolLogP\tHbondA\tHbondD\tPolarSA\tRbonds\tMolarRefractivity\tnAtoms\tSynthAcces\tAmesMutagenic\n")
    for data in totalfiles:
        parameter=data.split(sep="\t",maxsplit=11)
        MolWeight        =float(parameter[0])
        MolLogP          =float(parameter[1])
        HbondA           =int(parameter[2])
        HbondD           =int(parameter[3])
        PolarSA          =float(parameter[4])
        Rbonds           =int(parameter[5])
        MolarRefractivity=float(parameter[6])
        nAtoms           =int(parameter[7])
        #
        ListMolWeight.append(MolWeight)
        ListMolLogP.append(MolLogP)
        ListHbondA.append(HbondA)
        ListHbondD.append(HbondD)
        ListPolarSA.append(PolarSA)
        ListRbonds.append(Rbonds)
        ListMolarRefractivity.append(MolarRefractivity)
        ListnAtoms.append(nAtoms)
        #
        SynthAcces       =float(parameter[8])
        AmesMutagenic    =float(parameter[9])
        ListSynthAcces.append(SynthAcces)
        ListAmesMutagenic.append(AmesMutagenic)
        #
        seq              =str(parameter[10])
        ListSeq.append(seq)

    # Promedio
    # Min Max
    minmax(ListMolWeight,"Weight")
    minmax(ListMolLogP,"LogP")
    minmax(ListHbondA,"HbondA")
    minmax(ListHbondD,"HbondD")
    minmax(ListPolarSA,"TPSA")
    minmax(ListRbonds,"Rbonds")
    minmax(ListMolarRefractivity,"MR")
    minmax(ListnAtoms,"N°")
    #
    avgSynthAcces=minmax(ListSynthAcces,"SA")
    avgAmesMutagenic=minmax(ListAmesMutagenic,"Ames")
    #
    #Proceso de testeo
    for i in range(len(ListnAtoms)):
        booleanMW  = Heavside(float(875.985) , float(160),ListMolWeight[i])
        booleanLogP= Heavside(float(5), float(-4.88),ListMolLogP[i])
        booleanHBA = Heavside(float(15), float(0),ListHbondA[i])
        booleanHBD = Heavside(float(14), float(0),ListHbondD[i])
        booleanTPSA= Heavside(float(388.029), float(160),ListPolarSA[i])
        booleanRB  = Heavside(float(27), float(0),ListRbonds[i])
        booleanMR  = Heavside(float(225), float(40),ListMolarRefractivity[i])
        booleanX   = Heavside(float(62), float(0),ListnAtoms[i])
        #
        # Propiedades fisicoquímicas ScoreADME
        sumaADME  = sumalista([booleanMW,booleanLogP,booleanHBA,booleanHBD,booleanTPSA,booleanRB,booleanMR,booleanX])
        ScoreADME = (sumaADME/8)
        #
        ListScoreADME.append(ScoreADME)
        #print (str(sumaADME) + " " + str(ScoreADME)  )

    F2 = open("ScoreTotal_lineal_test.csv", 'w')
    for i in range(len(ListScoreADME)):
        _ScoreADME       =float(ListScoreADME[i])
        #
        _SynthAcces_dM   =float(ListSynthAcces[i])
        _AmesMutagenic_dM=float(ListAmesMutagenic[i])
        avgSynthAcces    =4.626848289181031
        avgAmesMutagenic =0.3092776990414759
        #
        ScoreSA = (1 + (math.log(avgSynthAcces/_SynthAcces_dM)))
        ScoreAM = (1 + (math.log(_AmesMutagenic_dM/avgAmesMutagenic)))
        #
        ScoreTotal = (_ScoreADME + ScoreSA + ScoreAM)
        ListScoreTotal.append(ScoreTotal)
        F2.write( str(ListSeq[i])+ "\t" + str(ScoreTotal) + "\n")

    F2.close()

    print ("## Final ScoreTotal ##")
    minmax(ListScoreTotal,"Score Total")

    print ("Fin programa")
    t2=time.time()
    resTime1=(t2-t1)
    print (' Reading took %.2f seconds. \n' % resTime1 )


text:
# -*- coding: utf-8 -*-
import os
import sys
import shutil
import logging

import numpy as np
import torch
import torchvision
from torch.utils.tensorboard import SummaryWriter

base_dir = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
sys.path.append(base_dir)
import skeleton
import resnet

LOGGER = logging.getLogger(__name__)

def parse_args():
    import argparse
    parser = argparse.ArgumentParser(description="fast converge cifar10")

    parser.add_argument('--dataset-base', type=str, default='./data')
    parser.add_argument('--dataset-name', type=str, default='cifar10')
    parser.add_argument('--batch', type=int, default=128)
    parser.add_argument('--epoch', type=int, default=200)
    parser.add_argument('--result-path', type=str, default=None)

    parser.add_argument('--lr', type=float, default=0.01)
    parser.add_argument('--hypergrad', action='store_true')
    parser.add_argument('--hypergrad-multi', action='store_true')

    parser.add_argument('--download', action='store_true')
    parser.add_argument('--seed', type=lambda x: int(x, 0), default=None)

    parser.add_argument('--log-filename', type=str, default='')
    parser.add_argument('--debug', action='store_true')
    return parser.parse_args()

def dataloaders(base, name, download, batch_size, device):
    num_class = 10 if name == 'cifar10' else 100
    dataset = torchvision.datasets.CIFAR10 if name == 'cifar10' else torchvision.datasets.CIFAR100

    train_dataset = dataset(
        root=base + '/cifar10',
        train=True,
        download=download
    )
    test_dataset = dataset(
        root=base + '/cifar10',
        train=False,
        download=download
    )

    post_transform = torchvision.transforms.Compose([
        torchvision.transforms.ToTensor(),
        torchvision.transforms.Normalize(
            mean=(0.4914, 0.4822, 0.4465),
            std=(0.2471, 0.2435, 0.2616)
        ),
    ])

    train_dataloader = skeleton.data.FixedSizeDataLoader(
        skeleton.data.TransformDataset(
            skeleton.data.prefetch_dataset(
                skeleton.data.TransformDataset(
                    train_dataset,
                    transform=torchvision.transforms.Compose([
                        skeleton.data.transforms.Pad(4),
                        post_transform
                    ]),
                    index=0
                ),
                num_workers=16
            ),
            transform=torchvision.transforms.Compose([
                skeleton.data.transforms.TensorRandomCrop(32, 32),
                skeleton.data.transforms.TensorRandomHorizontalFlip(),
                skeleton.data.transforms.Cutout(8, 8)
            ]),
            index=0
        ),
        steps=None,  # for prefetch using infinit dataloader
        batch_size=batch_size,
        num_workers=32,
        pin_memory=False,
        drop_last=True,
        shuffle=True,
        # sampler=skeleton.data.StratifiedSampler(train_dataset.targets)
    )

    test_dataloader = torch.utils.data.DataLoader(
        skeleton.data.prefetch_dataset(
            skeleton.data.TransformDataset(
                test_dataset,
                transform=torchvision.transforms.Compose([
                    post_transform,
                ]),
                index=0
            ),
            num_workers=16
        ),
        batch_size=batch_size // 2,
        shuffle=False,
        num_workers=0,
        pin_memory=False,
        drop_last=False
    )

    train_dataloader = skeleton.data.PrefetchDataLoader(train_dataloader, device=device, half=False)
    test_dataloader = skeleton.data.PrefetchDataLoader(test_dataloader, device=device, half=False)
    return num_class, int(len(train_dataset) // batch_size), train_dataloader, test_dataloader

def main():
    timer = skeleton.utils.Timer()

    args = parse_args()
    log_format = '[%(asctime)s] [%(levelname)s] [%(filename)s:%(lineno)03d] %(message)s'
    level = logging.DEBUG if args.debug else logging.INFO
    if not args.log_filename:
        logging.basicConfig(level=level, format=log_format, stream=sys.stderr)
    else:
        logging.basicConfig(level=level, format=log_format, filename=args.log_filename)

    torch.backends.cudnn.benchmark = True
    if args.seed is not None:
        skeleton.utils.set_random_seed_all(args.seed, deterministic=False)

    if args.result_path is not None:
        if os.path.exists(args.result_path):
            LOGGER.warning('remove exists result path at %s', args.result_path)
            shutil.rmtree(args.result_path)
        os.makedirs(args.result_path + '/models', exist_ok=True)
        writers = {
            'train': SummaryWriter(args.result_path + '/train'),
            'test': SummaryWriter(args.result_path + '/test'),
        }
    LOGGER.debug('args: %s', args)

    epochs = args.epoch
    batch_size = args.batch
    device = torch.device('cuda', 0)

    num_class, steps_per_epoch, train_loader, test_loader = dataloaders(
        args.dataset_base, args.dataset_name, args.download, batch_size, device
    )
    train_iter = iter(train_loader)

    model = resnet.resnet18(num_classes=num_class).to(device=device)
    criterion = torch.nn.CrossEntropyLoss(reduction='mean')
    metrics = skeleton.nn.Accuracy(1)

    optimizer = skeleton.optim.ScheduledOptimizer(
        [p for p in model.parameters() if p.requires_grad],
        skeleton.optim.HypergradSGD,
        steps_per_epoch=steps_per_epoch,
        lr=args.lr,
        momentum=0.9,
        weight_decay=5e-4,
        # nesterov=True,
        hypergrad_lr=(1e-2 * batch_size / 128) if args.hypergrad else 0,
        hypergrad_momentum=0.9 if args.hypergrad and not args.hypergrad_multi else 0,
        update_multiple=args.hypergrad_multi
    )

    class ModelLoss(torch.nn.Module):
        def __init__(self, model, criterion):
            super(ModelLoss, self).__init__()
            self.model = model
            self.criterion = criterion

        def forward(self, inputs, targets=None):
            logits = self.model(inputs)
            if targets is not None:
                loss = self.criterion(logits, targets)
                return logits, loss
            return logits
    model = ModelLoss(model, criterion)

    # warmup
    torch.cuda.synchronize()
    model.train()
    for _ in range(2):
        inputs, targets = next(train_iter)
        logits, loss = model(inputs, targets)
        loss.backward()
        model.zero_grad()
    torch.cuda.synchronize()
    timer('init')

    # train
    for epoch in range(epochs):
        model.train()
        train_metrics = {'loss': [], 'accuracy': []}
        timer('init', reset_step=True)
        for step in range(steps_per_epoch):
            inputs, targets = next(train_iter)
            logits, loss = model(inputs, targets)

            loss.sum().backward()
            accuracy = metrics(logits, targets)
            train_metrics['loss'].append(loss.detach())
            train_metrics['accuracy'].append(accuracy.detach())

            optimizer.update()
            optimizer.step()
            optimizer.zero_grad()
            LOGGER.debug(
                '[%02d] [%04d/%04d] [train] lr:%.5f loss:%.5f',
                epoch,
                step,
                steps_per_epoch,
                optimizer.get_learning_rate(),
                float(loss.detach())
            )
        timer('train')
        LOGGER.debug('lr:%s', ', '.join(['%.3f' % float(lr) for lr in optimizer.param_groups[0]['lr_per_layers']]))

        model.eval()
        test_metrics = {'loss': [], 'accuracy': []}
        with torch.no_grad():
            for inputs, targets in test_loader:
                logits, loss = model(inputs, targets)

                accuracy = metrics(logits, targets)

                test_metrics['loss'].append(loss.detach())
                test_metrics['accuracy'].append(accuracy.detach())
        timer('test')
        LOGGER.info(
            '[%02d] train loss:%.3f accuracy:%.3f test loss:%.3f accuracy:%.3f lr:%.5f',
            epoch,
            np.average([t.cpu().numpy() for t in train_metrics['loss']]),
            np.average([t.cpu().numpy() for t in train_metrics['accuracy']]),
            np.average([t.cpu().numpy() for t in test_metrics['loss']]),
            np.average([t.cpu().numpy() for t in test_metrics['accuracy']]),
            optimizer.get_learning_rate(),
        )

        if args.result_path is not None:
            state = {
                'epoch': epoch,
                'model': model.state_dict(),
                'optimizer': optimizer.state_dict()
            }
            torch.save(state, args.result_path + '/models/epoch%04d.pth' % (epoch))
            for key, value in train_metrics.items():
                writers['train'].add_scalar('metrics/{}'.format(key), np.average([t.cpu().numpy() for t in value]), global_step=epoch)
            for key, value in optimizer.params().items():
                if key == 'lr':
                    continue
                writers['train'].add_scalar('optimizer/{}'.format(key), value, global_step=epoch)
            writers['train'].add_scalar('optimizer/lr'.format(key), optimizer.get_learning_rate(), global_step=epoch)

            for key, value in test_metrics.items():
                writers['test'].add_scalar('metrics/{}'.format(key), np.average([t.cpu().numpy() for t in value]), global_step=epoch)

if __name__ == '__main__':
    main()


text:
import numpy, pylab
import moments
import OutOfAfrica

data = moments.Spectrum.from_file('YRI_CEU_CHB.fs')

params = numpy.array([2.10065897, 0.25066579, 0.22247642, 3.05297944,
					  0.09022469, 5.82773903, 3.79104318, 0.25730946,
					  0.12569788, 1.07182332, 0.36429414, 0.1108222,
					  0.07072507])

model = OutOfAfrica.OutOfAfrica(params, data.sample_sizes)

fig = pylab.figure(1)
fig.clear()
moments.Plotting.plot_1d_comp_multinom(model.marginalize([1,2]),
                                    data.marginalize([1,2]))
fig.savefig('1d_comp.pdf')

fig = pylab.figure(2)
fig.clear()
moments.Plotting.plot_single_2d_sfs(data.marginalize([2]), vmin=1)
fig.savefig('2d_single.pdf')

fig = pylab.figure(3)
fig.clear()
moments.Plotting.plot_2d_comp_multinom(model.marginalize([2]),
                                    data.marginalize([2]), vmin=1)
fig.savefig('2d_comp.pdf')

fig = pylab.figure(4, figsize=(8,10))
fig.clear()
moments.Plotting.plot_3d_comp_multinom(model, data, vmin=1)
fig.savefig('3d_comp.pdf')

pylab.show()

max_stars_repo_name: grahamgower/moments

max_issues_repo_name: grahamgower/moments

max_forks_repo_name: grahamgower/moments