text:
import weakref
import numpy as np
from collections import namedtuple

from itertools import count

from ...framework.observer import Observable, Listener

class Planner(object):
    """

    """

    class ValueState(Observable):
        """

        """

        @property
        def value(self):
            return self._value

        def __init__(self, value=None):
            super(Planner.ValueState, self).__init__()

            self._value = value if value is not None else 0.0

        def _set_value(self, value):
            """Changes the value of this ValueState nd sends a notification
            to all of the observers of this ValueState.

            Parameters
            ----------
            value : float
                The new value of this ValueState

            """
            change = value - self._value
            self._value = value
            self.dispatch('value_change', change)

    class DecisionState(ValueState, Listener):
        """

        """
        _instance = object()

        MaxAction = namedtuple('MaxAction', ['action', 'dsa'])

        class Action(Listener):
            """

            """
            _ids = count(0)
            _instance = object()

            ValueStateProbability = namedtuple('ValueStateProbability', ['vs', 'proba'])

            @property
            def q(self):
                """float : The q-value for each state-action."""
                return self._q

            @property
            def mdp(self):
                return self._model

            def __init__(self, token, parent, model):
                if token is not self._instance:
                    raise ValueError("Use 'create' to construct {0}".format(self.__class__.__name__))

                super(Planner.DecisionState.Action, self).__init__()

                self._parent = parent
                """:type: DecisionState"""
                self._model = model
                """StateAction"""

                self._q = 0.0
                self._errorbound = 0.0

                self._successors = {}
                """dict[MDPState, ValueStateProbability]"""

                self._mid = "%s.%s:%i" % (self.__class__.__module__, self.__class__.__name__, next(self._ids))

            def __repr__(self):
                return self._mid

            @classmethod
            def create(cls, parent, model):
                result = cls(cls._instance, parent, model)
                parent._action_values[model._action] = result
                result.compute_successors()
                result.compute_value()
                return result

            def debug(self):
                for s, vsp in self._successors.iteritems():
                    vs = vsp.vs
                    if vsp.vs is None:
                        vs = self._parent
                    print("\t\t%f: %s(%f)" % (vsp.proba, s, vs.value))
                print("\t\tr = %f, Q = %f" % (self._model.reward, self.q))

            def notify(self, event):
                if event.name != 'value_change':
                    return

                change = event.change
                if change < 0:
                    change = -change
                change *= self._parent._planner._gamma
                self._errorbound += change

            # TODO: need to move to Prioritized Sweeping?
            def compute_successors(self):
                model_succs = self._model.successor_probabilities

                # remove erstwhile successors
                for s, vsp in self._successors.items():
                    if s not in model_succs.keys():
                        if vsp.vs is not None:
                            vsp.vs.unsubscribe(self, 'value_change')
                        else:
                            self._parent.unsubscribe(self, 'value_change')
                        del self._successors[s]

                for s, proba in model_succs.iteritems():
                    if s not in self._successors:
                        succ = self._parent._planner._successor_value(s)
                        assert succ
                        succ.subscribe(self, 'value_change', {
                            'func': {
                                'attrib': 'change',
                                'callable': lambda x: proba * x
                            }
                        })

                        # use None reference as a proxy for parent to avoid creating
                        # a cycle of strong references
                        if succ == self._parent:
                            succ = None
                        self._successors[s] = Planner.DecisionState.Action.ValueStateProbability(
                            succ if succ._mid != self._parent._mid else weakref.proxy(succ), proba)
                    elif self._successors[s].proba != proba:
                        # update existing successor
                        vs = self._successors[s].vs
                        if vs is not None:
                            vs.subscribe(self, 'value_change', {
                                'func': {
                                    'attrib': 'change',
                                    'callable': lambda x: proba * x
                                }
                            })
                            self._successors[s] = self._successors[s]._replace(proba=proba)
                        else:
                            self._parent.subscribe(self, 'value_change', {
                                'func': {
                                    'attrib': 'change',
                                    'callable': lambda x: proba * x
                                }
                            })
                            self._successors[s] = self._successors[s]._replace(vs=self._parent, proba=proba)
                self._errorbound = np.infty

            def compute_value(self):
                self._q = 0.0
                for vsp in self._successors.itervalues():
                    succ = vsp.vs if vsp.vs is not None else self._parent
                    self._q += vsp.proba * succ.value

                self._q *= self._parent._planner._gamma
                self._q += self._model.reward
                self._errorbound = 0.0

            def update_value(self):
                if self._errorbound > 0:
                    self.compute_value()

        # -----------------------------
        # DecisionState
        # -----------------------------
        @property
        def state(self):
            return self._model.state

        @property
        def policy_action(self):
            return self._max.action

        @property
        def policy_model(self):
            return self._model

        def __init__(self, token, planner, model):
            if token is not self._instance:
                raise ValueError("Use 'create' to construct {0}".format(self.__class__.__name__))

            super(Planner.DecisionState, self).__init__()

            self._planner = planner
            """:type: Planner"""
            self._model = model
            """:type: StateData"""

            self._action_values = {}
            """dict[MDPAction, Action]"""
            self._max = Planner.DecisionState.MaxAction(None, None)
            """MaxAction"""

            self._errorbound = 0.0

            self._inbox = weakref.WeakSet()
            """:type: set[]"""

            # print("Planner.DecisionState: %s" % self._model.state)

        def __del__(self):
            pass
            # print("Planner.DecisionState.__del__: %s" % self._model.state)
            # for a in self._action_values.itervalues():
            #     for vsp in a._successors.itervalues():
            #         if vsp.vs is not None:
            #             del vsp

        def __str__(self, level=0):
            ret = "\t" * level + repr(self._model.state) + "\n"
            for av in self._action_values.itervalues():
                for vsp in av._successors.itervalues():
                    if vsp.vs._mid != self._mid:
                        ret += vsp.vs.__str__(level + 1)
            return ret

        def __repr__(self):
            return self._mid

        @classmethod
        def create(cls, planner, model, is_completion):
            result = cls(cls._instance, planner, model)
            if is_completion:
                planner._completions[model.state] = result
            planner._nonterminals[model.state] = result
            result._model.subscribe(result, 'mdp_change')
            result.initialize()
            return result

        def debug(self):
            print("State %s" % self.state)
            for a, dsa in self._action_values.iteritems():
                print("\tAction %s" % a)
                dsa.debug()
            print

        def notify(self, event):
            self._inbox.add(event.action)
            self._planner._inbox.add(self)

        def initialize(self):
            """Compute the initial value of this state.

            Note that this method is separate from the constructor, since otherwise
            a cycle in the MDP structure might cause the planner to try to create a
            new DecisionState object for a state in which another DecisionState object
            is still being constructed. This way, we can update the hash table mapping
            states of DecisionState objects after construction but before initialization,
            since the initialization is what causes the cycle.

            """
            for act, model in self._model.state_actions.iteritems():
                if len(self._action_values) <= 0 or act not in self._action_values:
                    # this can lead to the construction and initialization of other
                    # DecisionState objects
                    dsa = type(self).Action.create(weakref.proxy(self), model)
                    assert act in self._action_values
                    assert dsa == self._action_values[act]

                if self._max.dsa is None or self._max.dsa.q < dsa.q:
                    self._max = Planner.DecisionState.MaxAction(act, dsa)

            self._set_value(self._max.dsa.q)

        def propagate_mdp_change(self):
            for i in self._inbox:
                self._action_values[i].compute_successors()

            self._inbox.clear()
            self.propagate_value_change()

        def propagate_value_change(self):
            original_action = self._max.action

            assert self._action_values
            it = self._action_values.iteritems()

            act, dsa = it.next()
            dsa.update_value()
            self._max = Planner.DecisionState.MaxAction(act, dsa)

            while True:
                try:
                    act, dsa = it.next()
                    dsa.update_value()
                    if self._max.dsa.q < dsa.q:
                        self._max = Planner.DecisionState.MaxAction(act, dsa)
                except StopIteration:
                    break

            if self._max.action != original_action:
                # this is a new policy action
                self._planner._outbox[self] = original_action

            if self._max.dsa.q != self.value:
                self._set_value(self._max.dsa.q)

        def propagate_policy_change(self, original):
            if original != self._max.action:
                self.dispatch('policy_change')

    # -----------------------------
    # Planner
    # -----------------------------
    @property
    def mdp(self):
        return self._mdp

    def __init__(self, mdp, terminal, goal, gamma=None):
        self._mdp = mdp
        """:type: MDP"""

        self._terminal = terminal
        """:type: callable"""
        self._goal = goal
        """:type: callable"""

        self._gamma = gamma if gamma is not None else 1.0

        self._nonterminals = weakref.WeakValueDictionary()
        """:type: dict[MDPState, DecisionState] : Each key is the return value of the state()
        method for an MDPStateData. Each datum is the DecisionState object constructed
        with that MDPStateData"""

        self._completions = weakref.WeakValueDictionary()
        """:type: dict[MDPState, ValueState] : Each key is a basis state that appears
        in the value of successors() for sine MDPStateAction object. The data
        pointers are either DecisionState objects in nonterminals or ValueState
        objects representing terminal states"""

        self._inbox = weakref.WeakSet()
        """:type: set[DecisionState] : Includes all DecisionState objects that have not
        executed propagate_value_change since their observed MDPStateData object
        changed"""
        self._outbox = weakref.WeakKeyDictionary()
        """:type: dict[DecisionState, MDPAction] : Each key is a DecisionState object that
        has changed their policy action. The data is the original policy action
        before the first change. (The data can be used to avoid sending policy
        change notifications unnecessarily, when a DecisionState object switches
        back to its original policy action.)"""

    def initialize(self):
        self._mdp.initialize()

    def policy(self, state):
        """Outputs (after possibly computing) the optimal policy for
        the MDP.

        Parameters
        ----------
        state : MDPState
            The state at which to evaluate the optimal policy

        Returns
        -------
        DecisionState:
            An object specifying the optimal child action, as well as
            giving access to the StateData describing that state-action's
            behavior.

        """
        ds = self._policy(self.mdp.state_data(state))
        self.plan()
        return ds

    def update(self, state, act, succ):
        self._mdp.update(state, act, succ)

    def plan(self):
        for ds in self._inbox:
            ds.propagate_mdp_change()
        self._inbox.clear()

        self._propagate_changes()

        for ds, act in self._outbox.iteritems():
            if ds is not None:
                ds.propagate_policy_change(act)
        self._outbox = {}

    def debug(self):
        print("model:")
        for ds in self._nonterminals.itervalues():
            ds.debug()
        print("end model")

        print("value function:")
        self.write_value_function()
        print("end value function")

        print("policy:")
        self.write_policy()
        print("end policy")

    def write_value_function(self):
        for s, ds in self._nonterminals.iteritems():
            print("{0} {1}".format(s, ds.value))

    def write_policy(self):
        policy = {}
        for s, ds in self._nonterminals.iteritems():
            states = policy.setdefault(ds.policy_action, [])
            states.append(ds.state)

        for a, states in policy.iteritems():
            print("# {0}:".format(a))
            for s in states:
                print(s)
        print('\n')

    def _propagate_changes(self):
        pass

    def _policy(self, state_data, is_completion=False):
        """

        Parameters
        ----------
        state_data : MDPState.StateData
            The state at which the policy should be evaluated.

        Returns
        -------
        DecisionState

        """
        try:
            ds = self._nonterminals[state_data.state]
        except KeyError:
            ds = type(self).DecisionState.create(self, state_data, is_completion)
            assert state_data.state in self._nonterminals
            assert ds == self._nonterminals[state_data.state]
        return ds

    def _successor_value(self, successor):
        """

        Parameters
        ----------
        successor : MDPState
            A basis state that appears in teh value of successor_probabilities()
            for some MDPStateAction object

        Returns
        -------
        ValueState

        """
        try:
            vs = self._completions[successor]
        except KeyError:
            if self._terminal(successor):
                vs = Planner.ValueState(self._goal(successor))
            else:
                vs = self._policy(self.mdp.state_data(successor), is_completion=True)
                assert successor in self._completions
                assert vs == self._completions[successor]
        return vs

max_stars_repo_name: evenmarbles/rlpy

max_issues_repo_name: evenmarbles/rlpy

max_forks_repo_name: evenmarbles/rlpy


text:
"""ops.py"""
from __future__ import print_function

import torch
import torch.nn.functional as F

def recon_loss(x, x_recon):
    n = x.size(0)
    loss = F.binary_cross_entropy_with_logits(x_recon, x, size_average=False).div(n)
    return loss

def kl_divergence(mu, logvar):
    kld = -0.5*(1+logvar-mu**2-logvar.exp()).sum(1).mean()
    return kld

def permute_dims(z):
    assert z.dim() == 2

    B, _ = z.size()
    perm_z = []
    for z_j in z.split(1, 1):
        perm = torch.randperm(B).to(z.device)
        perm_z_j = z_j[perm]
        perm_z.append(perm_z_j)

    return torch.cat(perm_z, 1)

def attention_disentanglement(att1, att2):
    """
        Implementation of the equation 5 from author's paper
        ARGS:
            att1: attention map A1
            att2: attention map A2
        Returns:
            Float value representing the loss
    """
    numer = torch.min(att1,att2).sum()
    denomi = torch.add(att1,att2).sum()

    return 2 * (numer / denomi)

def get_cam(image, gcam):
    """ Combines input img and attention map for final colormap """
    gcam = gcam - np.min(gcam)
    gcam = gcam / np.max(gcam)
    h, w, d = image.shape
    gcam = cv2.resize(gcam, (w, h))
    gcam = cv2.applyColorMap(np.uint8(255 * gcam), cv2.COLORMAP_JET)
    gcam = np.asarray(gcam, dtype=np.float) + \
        np.asarray(image, dtype=np.float)
    gcam = 255 * gcam / np.max(gcam)
    gcam = np.uint8(gcam)
    return gcam

########### About GradCam used to compute the attention loss component of the AD-FactorVAE

from collections import OrderedDict

import cv2
import numpy as np
import torch.nn as nn
from torch.autograd import Variable
import os

class PropBase(object):

    def __init__(self, model, target_layer, cuda=True):
        self.model = model
        self.cuda = cuda
        if self.cuda:
            self.model.cuda()
        self.model.eval()
        self.target_layer = target_layer
        self.outputs_backward = OrderedDict()
        self.outputs_forward = OrderedDict()
        self.set_hook_func()

    def set_hook_func(self):
        raise NotImplementedError

    # set the target class as one others as zero. use this vector for back prop added by Lezi
    def encode_one_hot_batch(self, z, mu, logvar, mu_avg, logvar_avg):
        one_hot_batch = torch.FloatTensor(z.size()).zero_()
        return mu

    def forward(self, x):
        self.preds = self.model(x)
        self.image_size = x.size(-1)
        recon_batch, self.mu, self.logvar = self.model(x)
        return recon_batch, self.mu, self.logvar

    # back prop the one_hot signal
    def backward(self, mu, logvar, mu_avg, logvar_avg):
        self.model.zero_grad()
        z = self.model.reparameterize_eval(mu, logvar).cuda()
        one_hot = self.encode_one_hot_batch(z, mu, logvar, mu_avg, logvar_avg)

        if self.cuda:
            one_hot = one_hot.cuda()
        flag = 2
        if flag == 1:
            self.score_fc = torch.sum(F.relu(one_hot.cuda() * mu))
        else:
            self.score_fc = torch.sum(one_hot.cuda())
        self.score_fc.backward(gradient=one_hot, retain_graph=True)

    def get_conv_outputs(self, outputs, target_layer):
        for key, value in outputs.items():
            for module in self.model.named_modules():
                if id(module[1]) == key:
                    if module[0] == target_layer:
                        return value
        raise ValueError('invalid layer name: {}'.format(target_layer))

class GradCamDissen(PropBase):
    """
        Extend the baseline class and implements the FactorVAE pipeline
            with equation (5) from the paper
    """
    def __init__(self, VAE, discrim, target_layer, cuda=True):
        super(GradCamDissen, self).__init__(VAE,target_layer)
        self.D = discrim

    def set_hook_func(self):
        def func_b(module, grad_in, grad_out):
            self.outputs_backward[id(module)] = grad_out[0].cpu()

        def func_f(module, input, f_output):
            self.outputs_forward[id(module)] = f_output

        for module in self.model.named_modules():
            module[1].register_backward_hook(func_b)
            module[1].register_forward_hook(func_f)

    def forward(self, x):
        self.image_size = x.size(-1)
        recon_batch, self.mu, self.logvar, self.z = self.model(x)
        return recon_batch, self.mu, self.logvar, self.z

    def backward(self, mu, logvar, mu_avg, logvar_avg, flag=2):
        self.model.zero_grad()
        one_hot = mu  # self.encode_one_hot_batch(self.z, mu, logvar, mu_avg, logvar_avg)

        if flag == 1:
            self.score_fc = torch.sum(F.relu(one_hot.cuda() * mu))
        else:
            self.score_fc = torch.sum(one_hot.cuda())

        self.score_fc.backward(retain_graph=True)

    def normalize(self, grads):
        l2_norm = torch.sqrt(torch.mean(torch.pow(grads, 2))) + 1e-5
        return grads / l2_norm.item()

    def compute_gradient_weights(self):
        self.grads = self.normalize(self.grads)
        self.map_size = self.grads.size()[2:]
        self.weights = nn.AvgPool2d(self.map_size)(self.grads)

    def generate(self):
        # get gradient
        self.grads = self.get_conv_outputs(
            self.outputs_backward, self.target_layer)

        # compute weithts based on the gradient
        self.compute_gradient_weights()

        # get activation
        self.activiation = self.get_conv_outputs(
            self.outputs_forward, self.target_layer)

        self.weights.volatile = False
        self.activiation = self.activiation[None, :, :, :, :]
        self.weights = self.weights[:, None, :, :, :]

        gcam = F.conv3d(self.activiation, (self.weights.cuda()), padding=0, groups=len(self.weights))
        gcam = gcam.squeeze(dim=0)
        gcam = F.upsample(gcam, (self.image_size, self.image_size), mode="bilinear")
        gcam = torch.abs(gcam)

        return gcam

ext: py

reason: import numpy


text:
"""
Helper classes and functions do illustrate contribution matrices
as segments of circles.
"""
# coding: utf-8

import matplotlib.pyplot as pl

from shapely.geometry import Point, LineString, Polygon, MultiPolygon
from shapely.affinity import translate

import numpy as np
from scipy.optimize import root

import bfmplot as bp
import matplotlib.pyplot as pl

from shapely.geometry import CAP_STYLE, JOIN_STYLE

_BASE_CIRCLE_RESOLUTION = 64

class CircleCaps():

    def __init__(self,r,h,w=1/10,circle_resolution=_BASE_CIRCLE_RESOLUTION):
        #self.y = sorted([amount0, amount1],key=lambda x:-x)
        self.r = r
        self.w = w
        self.h = h
        self.circle_resolution = circle_resolution
        self.compute()

    def compute(self):
        point = Point(0,0)
        self.circle = point.buffer(self.r,resolution=self.circle_resolution)
        r = self.r
        h = self.h
        w = self.w
        box0 = Polygon([(-2*r,h+w), (2*r,h+w), (2*r,h+w+2*r),(-2*r,h+w+2*r)])
        box1 = Polygon([(-2*r,h), (2*r,h), (2*r,h-2*r),(-2*r,h-2*r)])
        self.cap0 = self.circle.intersection(box0)
        self.cap1 = self.circle.intersection(box1)

        filtered_polygons = list(filter(lambda x: x.area > 0,  [self.cap0,self.cap1]))
        self.all = MultiPolygon(filtered_polygons)

    def get_areas(self):
        return (self.cap0.area, self.cap1.area)

    def area(self):
        return self.all.area

class CircleCapSegments():

    def __init__(self,r,h,x_hi,x_lo,w=1/10,circle_resolution=_BASE_CIRCLE_RESOLUTION):
        #self.y = sorted([amount0, amount1],key=lambda x:-x)
        self.r = r
        self.w = w
        self.h = h
        self.x_lo = x_lo
        self.x_hi = x_hi
        self.circle_resolution = circle_resolution
        self.compute()

    def compute(self):
        point = Point(0,0)
        self.circle = point.buffer(self.r,resolution=self.circle_resolution)
        r = self.r
        h = self.h
        w = self.w
        x_lo = self.x_lo
        x_hi = self.x_hi
        box0 = Polygon([(-2*r,h+w), (2*r,h+w), (2*r,h+w+2*r),(-2*r,h+w+2*r)])
        box1 = Polygon([(-2*r,h), (2*r,h), (2*r,h-2*r),(-2*r,h-2*r)])
        self.cap0 = self.circle.intersection(box0)
        self.cap1 = self.circle.intersection(box1)

        box_lo_left  = Polygon([(-2*r,h+w/2), (x_lo,h+w/2), (x_lo,h-2*r),(-2*r,h-2*r)])
        box_lo_right = Polygon([(x_lo + w,h+w/2), (2*r,h+w/2), (2*r,h-2*r),(x_lo+w,h-2*r)])

        box_hi_left  = Polygon([(-2*r,h+w/2), (x_hi,h+w/2), (x_hi,h+2*r),(-2*r,h+2*r)])
        box_hi_right = Polygon([(x_hi+w,h+w/2), (2*r,h+w/2), (2*r,h+2*r),(x_hi+w,h+2*r)])

        self.seg00 = self.cap0.intersection(box_hi_left)
        self.seg01 = self.cap0.intersection(box_hi_right)
        self.seg10 = self.cap1.intersection(box_lo_left)
        self.seg11 = self.cap1.intersection(box_lo_right)

        filtered_polygons = list(filter(lambda x: x.area > 0, [self.seg00, self.seg01, self.seg10, self.seg11]))
        self.all = MultiPolygon(filtered_polygons)

    def get_areas(self):
        return [ [self.seg00.area, self.seg01.area], [self.seg10.area, self.seg11.area] ]

    def area(self):
        return self.all.area

class CircleCapPresentation():

    def __init__(self,y,r=1,w=1/10,area=None,circle_resolution=_BASE_CIRCLE_RESOLUTION):

        if area is not None:
            self.initial_r = r = np.sqrt(area/np.pi)
            self.target_area = area
        else:
            self.initial_r = r
            self.target_area = np.pi * r**2

        self.y = np.array(y)
        assert(self.y.shape == (2,))
        self.relative_y = self.y/self.y.sum()
        self.target_areas = self.target_area * self.relative_y
        self.w = w
        self.circle_resolution = circle_resolution

    def get_caps(self,r,h,w):
        caps = CircleCaps(r,h,w,circle_resolution=self.circle_resolution)
        return caps

    def get_areas(self,r,h,w):
        caps = self.get_caps(r,h,w)
        return np.array(caps.get_areas())

    def compute(self,tol=1e-3):
        r0 = self.initial_r
        h0 = (self.relative_y[0] - self.relative_y[1])
        #print(r0,h0)
        func = lambda param: self.get_areas(param[0], param[1], self.w) - self.target_areas
        solution = root(func,[r0,h0],tol=tol)
        self.caps = self.get_caps(solution.x[0],solution.x[1],self.w)
        self.r = solution.x[0]

        return self
        #print(solution)

class CircleCapPresentationConstR():

    def __init__(self,y,r=1,w=1/10,circle_resolution=_BASE_CIRCLE_RESOLUTION):

        self.r = r
        self.target_area = np.pi * r**2

        self.y = np.array(y)
        assert(self.y.shape == (2,))
        self.relative_y = self.y/self.y.sum()
        self.rel_target_areas = self.relative_y
        self.w = w
        self.circle_resolution = circle_resolution

    def get_caps(self,r,h,w):
        caps = CircleCaps(r,h,w,circle_resolution=self.circle_resolution)
        return caps

    def get_relative_areas(self,r,h,w):
        caps = self.get_caps(r,h,w)
        areas = np.array(caps.get_areas())
        return areas / areas.sum()

    def compute(self,tol=1e-3):
        h0 = 0

        def func(param):
            rel = self.get_relative_areas(self.r, param[0], self.w)
            trg = self.rel_target_areas
            return [ rel[0] - trg[0] ]

        solution = root(func,[h0],tol=tol)
        self.caps = self.get_caps(self.r,solution.x[0],self.w)
        return self

class CircleCapSegmentPresentation():

    def __init__(self,C,r=1,w=1/10,area=None,circle_resolution=_BASE_CIRCLE_RESOLUTION):

        if area is not None:
            self.initial_r = r = np.sqrt(area/np.pi)
            self.target_area = area
        else:
            self.initial_r = r
            self.target_area = np.pi * r**2

        self.C = np.array(C)
        assert(self.C.shape == (2,2))
        self.relative_C = self.C/self.C.sum()
        self.target_areas = (self.target_area * self.relative_C)
        self.w = w
        self.circle_resolution = circle_resolution

    def get_segs(self,r,h,xhi,xlo,w):
        segs = CircleCapSegments(r,h,xhi,xlo,w,circle_resolution=self.circle_resolution)
        return segs

    def get_areas(self,r,h,xhi,xlo,w):
        segs = self.get_segs(r,h,xhi,xlo,w)
        return np.array(segs.get_areas())

    def compute(self,tol=1e-3):
        r0 = self.initial_r

        h0 = 0
        xhi0 = 0
        xlo0 = 0

        func = lambda p: self.get_areas(p[0], p[1], p[2], p[3], self.w).ravel() - self.target_areas.ravel()
        solution = root(func,[r0,h0,xhi0,xlo0],tol=tol)
        p = solution.x.tolist() + [self.w]
        self.r = p[0]
        self.segs = self.get_segs(*p)

        return self

    def plot(self,
             ax=None,
             upper_color=None,
             lower_color=None,
             brighter_base=2,
             ec='w',
             linewidth=0.5,
             figsize=(6,6),
             ):

        if upper_color is None:
            self.upper_color = bp.epipack[1]
        else:
            self.upper_color = upper_color
        if lower_color is None:
            self.lower_color = bp.epipack[2]
        else:
            self.upper_color = upper_color
        self.upper_brighter = bp.brighter(self.upper_color, brighter_base)
        self.lower_brighter = bp.brighter(self.lower_color, brighter_base)

        if ax is None:
            fig, ax = pl.subplots(1,1,figsize=figsize)

        def _plot(geom, color, ec=None, linewidth=None):
            xs, ys = geom.exterior.xy
            ax.fill(xs, ys, fc=color, ec=ec, linewidth=linewidth)

        for geom, color in [
                (self.segs.seg00, self.upper_color),
                (self.segs.seg01, self.upper_brighter),
                (self.segs.seg10, self.lower_color),
                (self.segs.seg11, self.lower_brighter),
            ]:
            if geom.area > 0:
                _plot(geom, color, ec=ec, linewidth=linewidth)

        return ax

class CircleCapSegmentPresentationConstR(CircleCapSegmentPresentation):

    def __init__(self,C,r=1,w=1/10,circle_resolution=_BASE_CIRCLE_RESOLUTION):

        self.r = r
        self.target_area = np.pi * r**2

        self.C = np.array(C)
        assert(self.C.shape == (2,2))
        self.relative_C = self.C/self.C.sum()
        self.rel_target_areas = self.relative_C
        self.w = w
        self.circle_resolution = circle_resolution

    def get_segs(self,r,h,xhi,xlo,w):
        segs = CircleCapSegments(r,h,xhi,xlo,w,circle_resolution=self.circle_resolution)
        return segs

    def get_relative_areas(self,r,h,xhi,xlo,w):
        segs = self.get_segs(r,h,xhi,xlo,w)
        areas = np.array(segs.get_areas())
        return areas / areas.sum()

    def compute(self,tol=1e-3):
        r0 = self.r

        h0 = 0
        xhi0 = 0
        xlo0 = 0

        def func(p):
            areas = self.get_relative_areas(self.r, p[0], p[1], p[2], self.w).ravel()
            targets = self.rel_target_areas.ravel()
            return areas[1:] - targets[1:]

        solution = root(func,[h0,xhi0,xlo0],tol=tol)

        p = [self.r] + solution.x.tolist() + [self.w]
        self.segs = self.get_segs(*p)

        return self

class JoinedVectorAndMatrixPresentation():

    def __init__(self,
                 vector_presentation,
                 matrix_presentation,
                 xoff=0.0,
                 yoff=0.0,
                ):

        caps = vector_presentation.caps
        segs = matrix_presentation.segs

        cap_width = 2 * vector_presentation.r
        seg_r = matrix_presentation.r

        #if xoff is None:
        #    xoff = cap_width + seg_r
        if yoff == 0.0 and xoff == 0.0:
            yoff = -1.5*cap_width - seg_r

        self.cap0 = caps.cap0
        self.cap1 = caps.cap1
        self.w = vector_presentation.w

        self.r_caps = vector_presentation.r
        self.r_segs = matrix_presentation.r

        self.seg00 = translate(segs.seg00, xoff=xoff, yoff=yoff)
        self.seg01 = translate(segs.seg01, xoff=xoff, yoff=yoff)
        self.seg10 = translate(segs.seg10, xoff=xoff, yoff=yoff)
        self.seg11 = translate(segs.seg11, xoff=xoff, yoff=yoff)

        self.caps = [self.cap0, self.cap1]
        self.segs = [self.seg00, self.seg01, self.seg10, self.seg11]
        self.seg_matrix = [ [self.seg00, self.seg01], [self.seg10, self.seg11] ]

        filtered_polygons = list(filter(lambda x: x.area > 0, self.caps+self.segs))
        self.all = MultiPolygon(filtered_polygons)

    def plot(self,ax=None,upper_color=None,lower_color=None,brighter_base=2,ec='w',linewidth=0.5,figsize=(6,6)):
        if upper_color is None:
            self.upper_color = bp.epipack[1]
        else:
            self.upper_color = upper_color
        if lower_color is None:
            self.lower_color = bp.epipack[2]
        else:
            self.upper_color = upper_color
        self.upper_brighter = bp.brighter(self.upper_color, brighter_base)
        self.lower_brighter = bp.brighter(self.lower_color, brighter_base)

        if ax is None:
            fig, ax = pl.subplots(1,1,figsize=figsize)

        def _plot(geom, color, ec=None, linewidth=None):
            xs, ys = geom.exterior.xy
            ax.fill(xs, ys, fc=color, ec=ec, linewidth=linewidth)

        for geom, color in [
                (self.cap0, self.upper_color),
                (self.cap1, self.lower_color),
                (self.seg00, self.upper_color),
                (self.seg01, self.upper_brighter),
                (self.seg10, self.lower_color),
                (self.seg11, self.lower_brighter),
            ]:
            if geom.area > 0:
                _plot(geom, color, ec=ec, linewidth=linewidth)

        ax.axis('equal')

        self.ax = ax

        return ax

    def add_arrows_to_plot(self, lw=1.5, symmetrical_arrows=True):

        d = self.w
        arrowwidth = self.w*1.5
        if not symmetrical_arrows:
            left_max = right_max - 2*self.w
            right_max = max(self.r_caps, self.r_segs) + 4*self.w
        else:
            left_max = -max(self.r_caps, self.r_segs) - 2*self.w
            right_max = max(self.r_caps, self.r_segs) + 2*self.w

        first_set = { 'cap': self.cap0, 'seg0': self.seg00, 'seg1': self.seg10, 'max': left_max, 'color': self.upper_brighter }
        second_set = { 'cap': self.cap1, 'seg0': self.seg01, 'seg1': self.seg11, 'max': right_max, 'color': self.lower_brighter }
        capfactors = [+1,-1]
        segfactors = [4,1]
        zorders = [0,0]

        distance_between_circles = 2*self.r_caps

        for cfac, sfac, _set, zorder in zip(capfactors, segfactors, [first_set, second_set], zorders):

            # deal with left arrow#
            _x0, _x1, _y = _find_focused_edge_x(_set['cap'])
            xdist = _x1
            xcap, ycap = _find_mean_x_and_y_of_focused_edge_x(_set['cap'])
            ycap += cfac*d
            xseg0, yseg0 = _find_mean_x_and_y_of_focused_edge_x(_set['seg0'])
            yseg0 += sfac*d
            xseg1, yseg1 = _find_mean_x_and_y_of_focused_edge_x(_set['seg1'])
            yseg1 -= sfac*d

            if symmetrical_arrows:
                xseg0, yseg0 = _set['seg0'].centroid.coords[0]
                xseg1, yseg1 = _set['seg1'].centroid.coords[0]

            xstart0 = xcap-cfac*d

            base_coords = [
                            #(xcap-cfac*xdist*0.85, ycap),
                           # (_set['max'], ycap+0.5*(yseg0-ycap)),
                            (xstart0, ycap),
                            (xstart0, -0.5*distance_between_circles-self.r_caps),
                            (_set['max'], -1.25*distance_between_circles-self.r_caps),
                          ]
            up_coords = [
                            (_set['max'], yseg0),
                            (xseg0, yseg0),
                        ]
            low_coords = [
                            (_set['max'], yseg1),
                            (xseg1, yseg1),
                        ]
            x, y = zip(*(base_coords+up_coords))
            self.ax.plot(x, y, zorder=zorder, color=_set['color'], lw=lw)
            x, y = zip(*(base_coords+low_coords))
            self.ax.plot(x, y, zorder=zorder, color=_set['color'], lw=lw)

            if not symmetrical_arrows and cfac == 1:
                arr_fac = -cfac
            else:
                arr_fac = cfac

            arrow_one = [
                            (_set['max'], yseg1),
                            (_set['max']+arr_fac*arrowwidth*1.5, yseg1),
                            (_set['max'], yseg1+0.5*arrowwidth),
                        ]
            arrow_two = [
                            #(_set['max'], yseg0), #half arrow
                            (_set['max'], yseg0 - 0.5*arrowwidth), #full arrow
                            (_set['max']+arr_fac*arrowwidth*1.5, yseg0),
                            (_set['max'], yseg0+0.5*arrowwidth),
                        ]
            offset = -cfac*arrowwidth*.5 #set this to 0 for full arrow, not half arrow
            arrow_three = [
                            (xstart0+cfac*arrowwidth*.5, ycap,),
                            (xstart0+cfac*arrowwidth*.5, -1.1*self.r_caps),
                            (xstart0, -1.1*self.r_caps-2*arrowwidth),
                            (xstart0+offset, -1.1*self.r_caps),
                            (xstart0+offset, ycap),
                    ]
            for arrow in [arrow_one, arrow_two, arrow_three]:
                xs, ys = zip(*arrow)
                self.ax.fill(xs, ys, fc=_set['color'], ec='None',zorder=-1)

    def add_text_to_plot(self):

        d = self.w
        arrowwidth = self.w*1.5
        right_max = max(self.r_caps, self.r_segs) + 4*self.w

        first_set = { 'cap': self.cap0, 'seg0': self.seg00, 'seg1': self.seg10, 'color': self.upper_brighter }
        second_set = { 'cap': self.cap1, 'seg0': self.seg01, 'seg1': self.seg11, 'color': self.lower_brighter }
        capfactors = [+1,-1]
        segfactors = [4,1]
        zorders = [0,0]

        distance_between_circles = 2*self.r_caps
        labels = {'cap': ['(u)nvacc','(v)acc']}
        labels = {'cap': ['u','v']}

        vas = ['bottom', 'top']

        for i, (cfac, sfac, _set, zorder, va) in enumerate(zip(capfactors, segfactors, [first_set, second_set], zorders, vas)):

            # deal with left arrow#
            _x0, _x1, _y = _find_focused_edge_x(_set['cap'])
            xdist = _x1
            xcap, ycap = _find_mean_x_and_y_of_focused_edge_x(_set['cap'])
            ycap += cfac*d
            xseg0, yseg0 = _find_mean_x_and_y_of_focused_edge_x(_set['seg0'])
            yseg0 += sfac*d
            xseg1, yseg1 = _find_mean_x_and_y_of_focused_edge_x(_set['seg1'])
            yseg1 -= sfac*d

            self.ax.text(xcap, ycap, labels['cap'][i], transform=self.ax.transData, color='w', ha='center',va=va, fontstyle='italic')

def get_circular_vector_and_matrix_presentation(y, C, r=1, w=0.1):

    ypres = CircleCapPresentationConstR(y,r=r,w=w).compute()
    Cpres = CircleCapSegmentPresentation(C, area=ypres.caps.area()*C.sum(), w=w).compute()
    joined = JoinedVectorAndMatrixPresentation(ypres, Cpres)

    return joined

def _get_envelope_xy(geom,which):
    x, y = geom.envelope.exterior.xy
    if not which.lower() in [
                'lower left',
                'upper right',
                'upper left',
                'lower right',
            ]:
        if which == 'lower left':
            return min(x), min(y)
        if which == 'upper left':
            return min(x), max(y)
        if which == 'upper right':
            return max(x), max(y)
        if which == 'lower right':
            return max(x), min(y)
    else:
        raise ValueError(f"Don't know how to interpret {which=}")

def _find_focused_edge_x(geom):
    x, y = geom.exterior.xy
    ndx = np.where(np.diff(y) == 0.)[0][0]
    _x = sorted([ x[ndx], x[ndx+1] ])
    _y = [y[ndx]] * 2
    return _x + [_y[0]]

def _find_focused_edge_y(geom):
    x, y = geom.exterior.xy
    ndx = np.where(np.diff(x) == 0.)[0][0]
    _y = sorted([ y[ndx], y[ndx+1] ])
    _x = [x[ndx]] * 2
    return [_x[0] ] + _y

def _find_mean_x_and_y_of_focused_edge_x(geom):
    x0, x1, y = _find_focused_edge_x(geom)
    return np.mean([x0,x1]), y

def _find_mean_x_and_y_of_focused_edge_y(geom):
    x, y0, y1, = _find_focused_edge_y(geom)
    return x, np.mean([y0,y1]), y

def add_text_to_plot(joined_pres, c):
    pass

if __name__ == "__main__":

    import vaccontrib as vc

    C = vc.covid.get_reduced_vaccinated_susceptible_contribution_matrix_covid([1.,4,4,4,4],variant='alpha')
    K = vc.covid.get_next_generation_matrix_covid([1.,4,4,4,4],variant='alpha')

    C = vc.covid.get_reduced_vaccinated_susceptible_contribution_matrix_covid([1.,1,1,1,1],variant='delta')
    K = vc.covid.get_next_generation_matrix_covid([1.,1,1,1,1],variant='delta')

    C = vc.covid.get_reduced_vaccinated_susceptible_contribution_matrix_covid([6.,6,6,6,6],variant='delta')
    K = vc.covid.get_next_generation_matrix_covid([6.,6,6,6,6],variant='delta')

    C = np.random.rand(2,2)

    y = vc.get_eigenvector(K)
    y = y.sum(axis=0)
    y = np.array([y[0],y[1:].sum()])

    y = np.random.rand(2)
    pres = get_circular_vector_and_matrix_presentation(y, C)

    ax = pres.plot(figsize=(6,6))
    ax.axis('off')
    pres.add_arrows_to_plot()
    pres.add_text_to_plot()
    ax.get_figure().savefig('a.pdf')

    segsonly = CircleCapSegmentPresentation(C,area=C.sum())
    segsonly.compute()
    ax = segsonly.plot()

    pl.show()

max_stars_repo_name: benmaier/vaccontrib

max_issues_repo_name: benmaier/vaccontrib

max_forks_repo_name: benmaier/vaccontrib

reason: import numpy,from scipy


text:
#%% [markdown]
#  # Lung Vasculature Analysis
#  This notebook (.ipynb) is a working project for analyzing lung vasculature. It inculdes three parts:
#  1. converts skeleton analytical output (.xml) into .csv file.
#  2. calulates the length and average thickness of each segment.
#  3. makes two types of plots:
#      1. histogram of each dataset on length and thickness
#      2. average histogram on length and thickness (line plot with error bars)
#

#%%
get_ipython().run_line_magic('load_ext', 'autoreload')
get_ipython().run_line_magic('autoreload', '2')
import os, sys, re, io
import numpy as np
import pandas as pd
from tqdm import tqdm
import time
from core.fileop import DirCheck, ListFiles
import core.mkplot as mkplot

#%% [markdown]
#   ## Part 1:
#   Converting skeleton analytical output (.xml) into .csv file.
#   * Inputs: *.xml
#   * Outputs: *.csv
#   * Dependencies: xml, time, pandas, tqdm </br>
#
#   * *.xml file includes three sheets: nodes, points, and segments.
#   * Warning: the progress bar controled by `tqdm` is not functioning well. It can not overwrite itself and creates multiple lines.
#%% [markdown]
#  ### Functions

#%%
# import dependencies
import xml.etree.ElementTree as etree
from core.msxml import MSXmlReader

# function
def convert_xml_csv(ippath, oppath):
    filelist, fileabslist = ListFiles(ippath, extension='.xml')

    for idx, f in enumerate(filelist):
        filename = f.replace('.xml', '')
        ip = os.path.join(ippath, f)
        op = os.path.join(oppath, filename)

        print(ip)
        print(op)

        # create path
        if filename not in os.listdir(oppath):
            DirCheck(op)

            # convert *.xml to *.csv
            csv_all = MSXmlReader(ip)

            # save each spreadsheet into individual *.csv file
            for key, value in csv_all.items():
                oppath_tmp = os.path.join(op, key + '.csv')
                value.to_csv(oppath_tmp, index = False)

#%% [markdown]
#  ### Execution
#  To run the code, please change `path` to the directory hosts the raw data.

#%%
path = 'data'
ipdir = 'raw'
opdir = 'csv'
ippath = os.path.join(path, ipdir)
oppath = os.path.join(path, opdir)
# make dir
DirCheck(oppath)

# convert files in batch
convert_xml_csv(ippath, oppath)

#%% [markdown]
#   ## Part 2:
#   Calulating the length and average thickness of each segment.
#   * Inputs: nodes.csv, points.csv, segments.csv
#   * Outputs: segments_s.csv
#
#   `SegStats` extracts euclidean coordinates and thickness of each point, then calculate the total length and average thickness.
#%% [markdown]
#  ### Functions

#%%
# load dependencies
from core.filamentanalysis import SegStats, PNSCount

path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipdir = 'csv'
ippath = os.path.join(path, ipdir)
img_group = []

# function
def stats_calculator(ippath, oppath):
    imglist = [x for x in os.listdir(ippath) if not x.startswith('.')]

    var = ['df_nodes', 'df_points', 'df_segments']
    counts_combined = []
    names= []

    for img in imglist:
        filelist, fileabslist = ListFiles(os.path.join(ippath, img), extension='.csv')

        df_points = pd.read_csv(os.path.join(ippath, img, 'points.csv'))
        df_segments = pd.read_csv(os.path.join(ippath, img, 'segments.csv'))
        df_nodes = pd.read_csv(os.path.join(ippath, img,'nodes.csv'))

        opfilename = 'segments_s.csv'
        countfilename = 'count.csv'
        countfilename_combined = 'counts_combined.csv'

        if opfilename not in filelist:
            df_segments_s = SegStats(df_points, df_segments)
            df_segments_s.to_csv(os.path.join(oppath, img, opfilename), index = False)

        counts = (PNSCount(df_points, df_nodes, df_segments))
        counts_combined.append(counts)
        names.append(img)

    fileinfo = pd.read_csv(os.path.join('./par', 'lung_file_idx.csv'))
    print(names)

    img_group = []
    img_treatment = []
    for i in names:
        img_group.append(fileinfo[fileinfo['data_filename'] == i]['genotype'].item())
        img_treatment.append(fileinfo[fileinfo['data_filename'] == i]['treatment'].item())

    if countfilename_combined not in imglist:
        df_counts_combined = pd.DataFrame(counts_combined, columns= ['Points', 'Nodes', 'Segments'])
        df_counts_combined['Names'] = names
        df_counts_combined['Genotype'] = img_group
        df_counts_combined['Treatment'] = img_treatment
        df_counts_combined.to_csv(os.path.join(path, countfilename_combined), index = False)

#%% [markdown]
#  ### Execution
#  To run the code, please change `path` to the directory hosts the raw data.

#%%
path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipdir = 'csv'
opdir = 'csv'
ippath = os.path.join(path, ipdir)
oppath = os.path.join(path, opdir)
# make dir
DirCheck(oppath)

# convert files in batch
stats_calculator(ippath, oppath)

#%% [markdown]
#  ## Part 3:
#  Creating two sets of plots:
#  1. histogram of each dataset on length and thickness
#  2. average histogram on length and thickness (line plot with error bars)
#
#   * Inputs: segments_s.csv
#   * Outputs:
#       1. `histo/length/*.png`: frequency - length (µm)
#       2. `histo/thickness/*.png`: frequency - thickness (µm)
#       3. `histo_summary/length.png`: histogram in line plot style
#       4. `histo_summary/thickness.png`: histogram in line plot style
#
#   `SegStats` extracts euclidean coordinates and thickness of each point, then calculate
#  the total length and average thickness.
#
#
#  In the ouputs, the code renames "thickness" to "radius" to avoid confusion. Quotes from
#  Amira User's Manual
#  > As an estimate of the local thickness, the closest distance to the label
#  boundary (boundary distance map) is stored at every point in the *Spatial Graph*.
#  The attribute is named *thickness* and constitutes the *radius* of the circular cross-section
#  of the filament at a given point of the centerline.

#%%
# import depandencies
import matplotlib.pyplot as plt
import matplotlib.style as style
style.use('default')
import scipy.stats as stats
from core.mkplot import GroupImg, FindRange, IndividualHisto
from core.mkplot import make_individul_plots, make_individul_plots_all, make_merged_plots

#%%
path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipdir = 'csv'
opdir1 = 'plot'
opdir2 = 'histogram'
subfolder = ['histo', 'histo_summary']
ippath = os.path.join(path, ipdir)
oppath = os.path.join(path, opdir1, opdir2)
for i in subfolder:
    oppath_sub = os.path.join(oppath, i)
    DirCheck(oppath_sub)

#%%
# load fileinfo
fileinfo = pd.read_csv(os.path.join('./par', 'lung_file_idx.csv'))

columns = {
    'length': {
        'x_label': 'Length (µm)',
        'file_label': 'length',
    },
    'thickness': {
        'x_label': 'Radius (µm)',
        'file_label': 'radius',
    },
}

#%%
# plot individual histogram
make_individul_plots(ippath, oppath, fileinfo, columns)

#%%
path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipfile = 'counts_combined.csv'
ippath_stat = os.path.join(path, ipfile)
ippath_csv = os.path.join(path, 'csv')
oppath_stat = os.path.join(path, 'plot', 'histogram', 'histo_summary')

filestat = pd.read_csv(ippath_stat)
filestat['Group'] = filestat['Genotype'] + ' & ' + filestat['Treatment']
display(filestat)

make_individul_plots_all(ippath_csv, oppath, fileinfo, filestat, columns)

#%% [markdown]
#  Create plots with x-axis in different scales

#%%
# plot merged histogram in counts
make_merged_plots(ippath_csv, oppath, fileinfo, columns, frequency = False, x_max_factor = 0.07)
# plot merged histogram in frequency
make_merged_plots(ippath, oppath, fileinfo, columns, frequency = True, x_max_factor = 0.07)

#%%
# plot merged histogram in counts
make_merged_plots(ippath, oppath, fileinfo, columns, frequency = False, x_max_factor = 0.2)
# plot merged histogram in frequency
make_merged_plots(ippath, oppath, fileinfo, columns, frequency = True, x_max_factor = 0.2)

#%%
# plot merged histogram in counts
make_merged_plots(ippath, oppath, fileinfo, columns, frequency = False, x_max_factor = 1)
# plot merged histogram in frequency
make_merged_plots(ippath, oppath, fileinfo, columns, frequency = True, x_max_factor = 1)

#%% [markdown]
# ## Part 4
# ### 4-1: Plot paird scatter plot
# Module: seaborn

#%%
# from pandas.plotting import scatter_matrix
import seaborn as sns
sns.set(style="ticks", color_codes=True)

path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipfile = 'counts_combined.csv'
ippath_stat = os.path.join(path, ipfile)
ippath_csv = os.path.join(path, 'csv')
oppath_stat = os.path.join(path, 'plot', 'histogram', 'histo_summary')

counts = pd.read_csv(ippath_stat)
counts['Group'] = counts['Genotype'] + ' & ' + counts['Treatment']
display(counts)

#%%
g = sns.PairGrid(counts, vars = ['Points', 'Nodes', 'Segments'], hue = 'Group')
g.map_diag(plt.hist)
g.map_offdiag(plt.scatter)
g.add_legend()
g.savefig(os.path.join(oppath_stat, "pairgrid.png"))

#%% [markdown]
# ### 4-2: Plot histogram comparison
# Module: seaborn

#%%
labels, uniques = pd.factorize(counts['Group'])
for i in uniques:
    print(i)
    tmp_df = counts[counts['Group'] == i]
    data_merge = []
    filename = tmp_df['Names']

    df_all = []
    for index, value in filename.items():
        print(value)
        tmp_df2 = pd.read_csv(os.path.join(ippath_csv, value, 'segments_s.csv'))
        tmp_df2['Filename'] = value
        df_all.append(tmp_df2)

    df_all_con = pd.concat(df_all)
    # display(df_all_con)
    g = sns.FacetGrid(df_all_con, col="Filename")
    g.map(plt.hist, 'thickness', density = True)

#%% [markdown]
# ## Part 5: histogram standardization
#

#%%
from core.mkplot import histo_standardize, make_merged_plots_std

path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipfile = 'counts_combined.csv'
ippath_stat = os.path.join(path, ipfile)
ippath_csv = os.path.join(path, 'csv')

histo_standardize(ippath_csv)

#%% [markdown]
# ## Part 6: Plots with standardize data
#

#%%
path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipfile = 'counts_combined.csv'
ippath_stat = os.path.join(path, ipfile)
ippath_csv = os.path.join(path, 'csv')
oppath = os.path.join(path, 'plot', 'histogram')
DirCheck(oppath)

columns = {
    'length': {
        'x_label': 'Standard Deviation',
        'file_label': 'length',
    },
    'thickness': {
        'x_label': 'Standard Deviation',
        'file_label': 'radius',
    },
}

make_merged_plots_std(ippath_csv, oppath, fileinfo, columns, opdir = 'histo_summary_std', filename = 'segments_s_std.csv', frequency = False, x_max_factor = 1)
make_merged_plots_std(ippath_csv, oppath, fileinfo, columns, opdir = 'histo_summary_std', filename = 'segments_s_std.csv', frequency = True, x_max_factor = 1)

make_merged_plots_std(ippath_csv, oppath, fileinfo, columns, opdir = 'histo_summary_std', filename = 'segments_s_std.csv', frequency = False, x_max_factor = 0.2)
make_merged_plots_std(ippath_csv, oppath, fileinfo, columns, opdir = 'histo_summary_std', filename = 'segments_s_std.csv', frequency = True, x_max_factor = 0.2)

#%% [markdown]
# # Part 6
#     Plot the Points, Nodes, and Segment Count in Bokeh with Holoview

#%%
import numpy as np
import pandas as pd
import holoviews as hv
from holoviews import opts, Cycle

#%%
hv.extension('bokeh')

path = '/Volumes/LaCie_DataStorage/Woo-lungs/2019'
ipfile = 'counts_combined.csv'
ippath = os.path.join(path, ipfile)

counts = pd.read_csv(ippath)

f1 = hv.Scatter((zip(counts.Points.items(), counts.Nodes.items())), ['Points'], ['Nodes'])
f2 = hv.Scatter((zip(counts.Points.items(), counts.Segments.items())), ['Points'], ['Segments'])
f3 = hv.Scatter((zip(counts.Nodes.items(), counts.Points.items())), ['Nodes'], ['Points'])
f4 = hv.Scatter((zip(counts.Nodes.items(), counts.Segments.items())), ['Nodes'], ['Segments'])
f5 = hv.Scatter((zip(counts.Segments.items(), counts.Points.items())), ['Segments'], ['Points'])
f6 = hv.Scatter((zip(counts.Segments.items(), counts.Nodes.items())), ['Segments'], ['Nodes'])

f1 + f2 + f3 + f4 + f5 + f6

#%%
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
from holoviews.operation import gridmatrix
from bokeh.sampledata.iris import flowers
from bokeh.palettes import brewer
import bokeh.models as bmod

counts = pd.read_csv(ippath)
colors = brewer["Spectral"][len(counts.Genotype.unique()) + 1]
colormap = {counts.Genotype.unique()[i] : colors[i] for i in range(len(counts.Genotype.unique()))}
colors = [colormap[x] for x in counts.Genotype]

print(colormap)

iris_ds = hv.Dataset(counts).groupby('Genotype').overlay()

point_grid = gridmatrix(iris_ds, chart_type=hv.Points)

(point_grid).opts(
    opts.Bivariate(bandwidth=0.5, cmap=hv.Cycle(values = colors)),
    opts.Points(size=5, alpha=0.5),
    opts.NdOverlay(batched= False))

#%%

reason: import numpy,import scipy


text:
import logging
from abc import ABC, abstractmethod
from math import log

import numpy as np
from ccmi import CCMI
from scipy.spatial import cKDTree
from sklearn.neighbors import KDTree
from scipy.special import digamma

class Estimator(ABC):

    @abstractmethod
    def estimator_params(self):
        pass

class CMIEstimator(Estimator):

    @abstractmethod
    def compute_cmi(self, x, y, z):
        pass

class MIEstimator(CMIEstimator):

    @abstractmethod
    def compute_mi(self, x, y):
        pass

    # If estimator can do MI, CMI is just difference. Please override if needed.
    def compute_cmi(self, x, y, z):
        yz = np.concatenate([y, z], axis=1)
        mi_x_yz = self.compute_mi(x, yz)
        mi_x_z = self.compute_mi(x, z)
        logging.debug(f'{mi_x_yz} - {mi_x_z}')
        return mi_x_yz - mi_x_z

def _mixed_KSG(x, y, k=3, leafsize=16, use_python=True):
    '''
        Estimate the mutual information I(X;Y) of X and Y from samples {x_i, y_i}_{i=1}^N
        Using *Mixed-KSG* mutual information estimator

        Input: x: 2D array of size N*d_x (or 1D list of size N if d_x = 1)
        y: 2D array of size N*d_y (or 1D list of size N if d_y = 1)
        k: k-nearest neighbor parameter

        Output: one number of I(X;Y)
    '''
    assert len(x) == len(y), "Lists should have same length"
    assert k <= len(x)-1, "Set k smaller than num. samples - 1"
    N = len(x)
    if x.ndim == 1:
        x = x.reshape((N, 1))
    # dx = len(x[0])
    if y.ndim == 1:
        y = y.reshape((N, 1))
    # dy = len(y[0])
    data = np.concatenate((x, y), axis=1)

    tree_xy = KDTree(data, metric='chebyshev', leaf_size=leafsize) if use_python else cKDTree(data, leafsize=leafsize)
    tree_x = KDTree(x, metric='chebyshev', leaf_size=leafsize) if use_python else cKDTree(x, leafsize=leafsize)
    tree_y = KDTree(y, metric='chebyshev', leaf_size=leafsize) if use_python else cKDTree(y, leafsize=leafsize)

    def query_ball(tree, points, r):
        return (tree.query_radius(points, r, count_only=True) if use_python
                else tree.query_ball_point(points, r, p=np.inf, return_length=True))

    eps = 1e-15
    kp = np.full((N,), k)
    knn_dis = tree_xy.query(data, k+1)[0][:, k] if use_python else tree_xy.query(data, k+1, p=np.inf)[0][:, k]
    idxs = np.where(knn_dis == 0)[0]
    if idxs.size:
        kp[idxs] = query_ball(tree_xy, data[idxs], eps)

    nx = query_ball(tree_x, x, np.maximum(knn_dis - eps, eps))
    ny = query_ball(tree_y, y, np.maximum(knn_dis - eps, eps))

    ans = np.mean(digamma(kp)) - np.mean(np.log(nx)) - np.mean(np.log(ny)) + np.log(N)
    return ans

class GKOVEstimator(MIEstimator):

    def __init__(self, *args, **kwargs):
        super().__init__()
        self.k = kwargs.get('k', 3)
        self.use_python = kwargs.get('use_python', True)

    def estimator_params(self):
        return {'k': self.k}

    def compute_mi(self, x, y):
        k = max(1, int(np.sqrt(len(x))/5)) if self.k == 'adaptive' else self.k
        return _mixed_KSG(x, y, k=k, use_python=self.use_python) / log(2)

class CCMIEstimator(MIEstimator):

    def estimator_params(self):
        return {}  # No configurability for now

    def compute_mi(self, x, y):
        model = CCMI(x, y, np.ones((len(x), 1)), tester='Classifier', metric='donsker_varadhan',
                     num_boot_iter=10, h_dim=64, max_ep=20)
        mi = model.get_mi_est(np.concatenate([x, y], axis=1)) / log(2)
        return mi

    def compute_cmi(self, x, y, z):
        model = CCMI(x, y, z, tester='Classifier', metric='donsker_varadhan',
                     num_boot_iter=10, h_dim=64, max_ep=20)
        cmi = model.get_cmi_est() / log(2)
        return cmi

max_stars_repo_name: hpi-epic/mpcsl

max_issues_repo_name: hpi-epic/mpcsl

max_forks_repo_name: hpi-epic/mpcsl


text:
"""
Original ctr_data module developed by T. Trainor and F. Heberling
There is no longer a CTR Data object
The contents of this module are now only responsible for handling correction factors
All plotting and CTR data handling are now done in the wxIntegrator developed by C. Biwer
J. Stubbs April 2015

The key references for how the geometric
corrections are applied include:

* E. Vlieg, J. Appl. Cryst. (1997). 30, 532-543
* C. Schlepuetz et al, Acta Cryst. (2005). A61, 418-425
"""
##############################################################################

import types, copy
import numpy as num
import time

from mathutil import cosd, sind

from active_area import active_area
import gonio_psic

DEBUG = False

##############################################################################

##############################################################################
def image_point_F(scan,point,I='I',Inorm='io',Ierr='Ierr',Ibgr='Ibgr', transm='transm',
                  corr_params={}, preparsed=False):
    """
    compute F for a single scan point in an image scan
    """
    d = {'I':0.0,'Inorm':0.0,'Ierr':0.0,'Ibgr':0.0, 'transm':0.0, 'F':0.0,'Ferr':0.0,
         'ctot':1.0,'alpha':0.0,'beta':0.0}
    d['I']     = scan[I][point]
    d['Inorm'] = scan[Inorm][point]
    d['Ierr']  = scan[Ierr][point]
    d['Ibgr']  = scan[Ibgr][point]
    d['transm'] = scan[transm][point]

    if corr_params == None:
        d['ctot'] = 1.0
        scale = 1.0
    else:
        # compute correction factors
        scale  = corr_params.get('scale')
        if scale == None: scale = 1.
        scale  = float(scale)
        corr = _get_corr(scan,point,corr_params, preparsed)
        if corr == None:
            d['ctot'] = 1.0
        else:
            d['ctot']  = corr.ctot_stationary()
            d['alpha'] = corr.gonio.pangles['alpha']
            d['beta']  = corr.gonio.pangles['beta']

    # compute F
    if d['I'] <= 0.0 or d['Inorm'] <= 0.:
        d['F']    = 0.0
        d['Ferr'] = 0.0
    else:
        scale = scale/d['transm'] * d['ctot']/d['Inorm']
        #scale = scale * d['ctot']/d['Inorm']
        d['F']    = num.sqrt(scale*d['I'])
        d['Ferr'] = 0.5 * scale**0.5 * d['Ierr']/d['I']**0.5
    return d

##############################################################################
def _get_corr(scan,point,corr_params,preparsed=False):
    """
    get CtrCorrection instance
    """
    geom   = corr_params.get('geom','psic')
    beam   = corr_params.get('beam_slits',{})
    det    = corr_params.get('det_slits')
    sample = corr_params.get('sample')
    # get gonio instance for corrections
    if geom == 'psic':
        gonio = gonio_psic.psic_from_spec(scan['G'],preparsed=preparsed)
        _update_psic_angles(gonio,scan,point)
        corr  = CtrCorrectionPsic(gonio=gonio,beam_slits=beam,
                                  det_slits=det,sample=sample)
    else:
        print "Geometry %s not implemented" % geom
        corr = None
    return corr

##############################################################################
def _update_psic_angles(gonio,scan,point,verbose=True):
    """
    given a psic gonio instance, a scandata object
    and a scan point, update the gonio angles...
    """
    try:
        npts = int(scan.dims[0])
    except:
        npts = scan.get('dims', (1,0))[0]
    try:
        scan_name = scan.name
    except:
        scan_name = ''
    #
    try:
      if type(scan['phi']) == types.FloatType:
          phi=scan['phi']
      elif len(scan['phi']) == npts:
          phi=scan['phi'][point]
    except:
        phi=None
    if phi == None and verbose==True:
        print "Warning no phi angle:", scan_name
    #
    try:
        if type(scan['chi']) == types.FloatType:
            chi=scan['chi']
        elif len(scan['chi']) == npts:
            chi=scan['chi'][point]
    except:
        chi = None
    if chi == None and verbose==True:
        print "Warning no chi angle", scan_name
    #
    try:
        if type(scan['eta']) == types.FloatType:
            eta=scan['eta']
        elif len(scan['eta']) == npts:
            eta=scan['eta'][point]
    except:
        eta = None
    if eta == None and verbose==True:
        print "Warning no eta angle", scan_name
    #
    try:
        if type(scan['mu']) == types.FloatType:
            mu=scan['mu']
        elif len(scan['mu']) == npts:
            mu=scan['mu'][point]
    except:
        mu = None
    if mu == None and verbose==True:
        print "Warning no mu angle", scan_name
    #
    try:
        if type(scan['nu']) == types.FloatType:
            nu=scan['nu']
        elif len(scan['nu']) == npts:
            nu=scan['nu'][point]
    except:
        nu = None
    if nu == None and verbose==True:
        print "Warning no nu angle", scan_name
    #
    try:
        if type(scan['del']) == types.FloatType:
            delta=scan['del']
        elif len(scan['del']) == npts:
            delta=scan['del'][point]
    except:
        delta = None
    if delta == None and verbose==True:
        print "Warning no del angle", scan_name
    #
    gonio.set_angles(phi=phi,chi=chi,eta=eta,
                     mu=mu,nu=nu,delta=delta)

##############################################################################
class CtrCorrectionPsic:
    """
    Data point operations / corrections for Psic geometry

    Notes:
    ------
    All correction factors are defined such that the
    measured data is corrected by multiplying times
    the correction:
      Ic  = Im*ct
    where
      Im = Idet/Io = uncorrected (measured) intensity

    In other words we use the following formalism:
      Im = (|F|**2)* prod_i(Xi)
    where Xi are various (geometric) factors that
    influence the measured intensity.  To get the
    structure factor:
      |F| = sqrt(Im/prod_i(Xi)) = sqrt(Im* ct)
    and
      ct = prod_i(1/Xi) = prod_i(ci)
      ci = 1/Xi

    If there is an error or problem in the routine for a specific
    correction factor, (e.g. divide by zero), the routine should
    return a zero.  This way the corrected data is zero'd....

    * The correction factors depend on the goniometer geometry
      gonio = gonio_psic.Psic instance

    * The slits settings are needed.  Note if using a large area detector
      you may pass det_slits = None and just spill off will be computed
        beam_slits = {'horz':.6,'vert':.8}
        det_slits = {'horz':20.0,'vert':10.5}
      these are defined wrt psic phi-frame:
        horz = beam/detector horz width (total slit width in lab-z,
               or the horizontal scattering plane)
        vert = detector vert hieght (total slit width in lab-x,
               or the vertical scattering plane)

    * A sample description is needed.
      sample = {}
        sample['dia'] = is taken as the diameter of a round sample
                        mounted on center. if dia<=0 then we either use
                        the polygon description or ignore the sample.
        sample['polygon'] = [[1.,1.], [.5,1.5], [-1.,1.],
                             [-1.,-1.],[0.,.5],[1.,-1.]]
        sample['angles']  = {'phi':108.0007,'chi':0.4831}

        polygon = [[x,y,z],[x,y,z],[x,y,z],....]
                  is a list of vectors that describe the shape of
                  the sample.  They should be given in general lab
                  frame coordinates.

        angles = {'phi':0.,'chi':0.,'eta':0.,'mu':0.}
                 are the instrument angles at which the sample
                 vectors were determined.

    Note: the lab frame coordinate systems is defined such that:
    x is vertical (perpendicular, pointing to the ceiling of the hutch)
    y is directed along the incident beam path
    z make the system right handed and lies in the horizontal scattering plane
    (i.e. z is parallel to the phi axis)

    The center (0,0,0) of the lab frame is the rotation center of the instrument.

    If the sample vectors are given at the flat phi and chi values and with
    the correct sample hieght (sample Z set so the sample surface is on the
    rotation center), then the z values of the sample vectors will be zero.
    If 2D vectors are passed we therefore assume these are [x,y,0].  If this
    is the case then make sure:
    angles = {'phi':flatphi,'chi':flatchi,'eta':0.,'mu':0.}

    The easiest way to determine the sample coordinate vectors is to take a picture
    of the sample with a camera mounted such that is looks directly down the omega
    axis and the gonio angles set at the sample flat phi and chi values and
    eta = mu = 0. Then find the sample rotation center and measure the position
    of each corner (in mm) with up being the +x direction, and downstream
    being the +y direction.

    Note this routine does not correct for attenuation factors.

    """
    def __init__(self,gonio=None,beam_slits={},det_slits=None,sample={}):
        """
        Initialize

        Parameters:
        -----------
        * gonio is a goniometer instance used for computing reciprocal
          lattice indicies and psuedo angles from motor angles
        * beam_slits are dictionary defining the incident beam aperature
        * det_slits are a dictionary defining the detector aperature
        * sample is a dictionary describing the sample geometry
        (see the instance documentation for more details)
        """
        self.gonio      = gonio
        if self.gonio.calc_psuedo == False:
            self.gonio.calc_psuedo = True
            self.gonio._update_psuedo()
        self.beam_slits = beam_slits
        self.det_slits  = det_slits
        self.sample     = sample
        # fraction horz polarization
        self.fh         = 1.0

    ##########################################################################
    def ctot_stationary(self,plot=False,fig=None):
        """
        correction factors for stationary measurements (e.g. images)
        """
        cp = self.polarization()
        cl = self.lorentz_stationary()
        ca = self.active_area(plot=plot,fig=fig)
        ct = (cp)*(cl)*(ca)
        if plot == True:
            print "Correction factors (mult by I)"
            print "   Polarization=%f" % cp
            print "   Lorentz=%f" % cl
            print "   Area=%f" % ca
            print "   Total=%f" % ct
        return ct

    ##########################################################################
    def lorentz_stationary(self):
        """
        Compute the Lorentz factor for a stationary (image)
        measurement.  See Vlieg 1997

        Measured data is corrected for Lorentz factor as:
          Ic  = Im * cl
        """
        beta  = self.gonio.pangles['beta']
        cl = sind(beta)
        return cl

    ##########################################################################
    def polarization(self,):
        """
        Compute polarization correction factor.

        For a horizontally polarized beam (polarization vector
        parrallel to the lab-frame z direction) the polarization
        factor is normally defined as:
           p = 1-(cos(del)*sin(nu))^2
        For a beam with mixed horizontal and vertical polarization:
           p = fh( 1-(cos(del)*sin(nu))^2 ) + (1-fh)(1-sin(del)^2)
        where fh is the fraction of horizontal polarization.

        Measured data is corrected for polarization as:
          Ic  = Im * cp = Im/p
        """
        fh    = self.fh
        delta = self.gonio.angles['delta']
        nu    = self.gonio.angles['nu']
        p = 1. - ( cosd(delta) * sind(nu) )**2.
        if fh != 1.0:
            p = fh * c_p + (1.-fh)*(1.0 - (sind(delta))**2.)
        if p == 0.:
            cp = 0.
        else:
            cp = 1./p

        return cp

    ##########################################################################
    def active_area(self,plot=False,fig=None):
        """
        Compute active area correction (c_a = A_beam/A_int**2)

        Use to correct scattering data for area effects,
        including spilloff, A_int/A_beam and normailization
        to unit surface area (1/A_beam), i.e.
            Ic = Im * ca = Im/A_ratio
            A_ratio = A_int/(A_beam**2)
        where
            A_int = intersection area (area of beam on sample
                    viewed by detector)
            A_beam = total beam area
        """
        if self.beam_slits == {} or self.beam_slits == None:
            print "Warning beam slits not specified"
            return 1.0
        alpha = self.gonio.pangles['alpha']
        beta  = self.gonio.pangles['beta']
        if plot == True:
            print 'Alpha = ', alpha, ', Beta = ', beta
        if alpha < 0.0:
            print 'alpha is less than 0.0'
            return 0.0
        elif beta < 0.0:
            print 'beta is less than 0.0'
            return 0.0

        # get beam vectors
        bh = self.beam_slits['horz']
        bv = self.beam_slits['vert']
        beam = gonio_psic.beam_vectors(h=bh,v=bv)

        # get det vectors
        if self.det_slits == None:
            det = None
        else:
            dh = self.det_slits['horz']
            dv = self.det_slits['vert']
            det  = gonio_psic.det_vectors(h=dh,v=dv,
                                          nu=self.gonio.angles['nu'],
                                          delta=self.gonio.angles['delta'])
        # get sample poly
        if type(self.sample) == types.DictType:
            sample_dia    = self.sample.get('dia',0.)
            sample_vecs   = self.sample.get('polygon',None)
            sample_angles = self.sample.get('angles',{})
            #
            if sample_vecs != None and sample_dia <= 0.:
                sample = gonio_psic.sample_vectors(sample_vecs,
                                                   angles=sample_angles,
                                                   gonio=self.gonio)
            elif sample_dia > 0.:
                sample = sample_dia
            else:
                sample = None
        else:
            sample = self.sample

        # compute active_area
        (A_beam,A_int) = active_area(self.gonio.nm,ki=self.gonio.ki,
                                     kr=self.gonio.kr,beam=beam,det=det,
                                     sample=sample,plot=plot,fig=fig)
        if A_int == 0.:
            ca = 0.
        else:
            ca = A_beam/(A_int**2)

        return ca

##############################################################################


text:
import torch
import torch.nn as nn
import numpy as np
from models import ModelBuilder
from torch.autograd import Variable
import torch.optim as optim
import torch.backends.cudnn as cudnn
from torch.utils.data import DataLoader
import os
import random
import nibabel as nib
from utils import AverageMeter
from distutils.version import LooseVersion
import argparse
from dataset import ValDataset_full as ValDataset
import SimpleITK as sitk

# save prediction results in the format of online submission
def visualize_result(name, pred, args):
    _, name1, _, name2 = name.split("/")
    _, _, _, _, _, name3, _ = name2.split(".")
    pred = sitk.GetImageFromArray(pred)
    sitk.WriteImage(pred, args.result + "/VSD"+"."+ str(name1) + '.'+ str(name3)+ '.mha')

def accuracy(pred, mask, label):
    # columns in score is (# pred, # label, pred and label)
    score = np.zeros([3,3])

    # compute Enhance score (label==4) in the first line
    score[0,0] = np.count_nonzero(pred * mask == 4)
    score[0,1] = np.count_nonzero(label == 4)
    score[0,2] = np.count_nonzero(pred * mask * label == 16)

    # compute Core score (label == 1,3,4) in the second line
    pred[pred > 2] = 1
    label[label > 2] = 1
    score[1,0] = np.count_nonzero(pred * mask == 1)
    score[1,1] = np.count_nonzero(label == 1)
    score[1,2] = np.count_nonzero(pred * mask * label == 1)

    # compute Whole score (all labels) in the third line
    pred[pred > 1] = 1
    label[label > 1] = 1
    score[2,0] = np.count_nonzero(pred * mask == 1)
    score[2,1] = np.count_nonzero(label == 1)
    score[2,2] = np.count_nonzero(pred * mask * label == 1)
    return score

def test(test_loader, model, args):
    # switch to evaluate mode
    model.eval()

    # initialization
    num_ignore = 0
    num_images = int(len(test_dir)/args.num_input)
    # columns in score is (# pred, # label, pred and label)
    # lines in score is (Enhance, Core, Whole)
    dice_score = np.zeros([num_images, 3]).astype(float)

    for i, sample in enumerate(test_loader):
        image = sample['images'].float().cuda()
        label = sample['labels'].long().cuda()
        mask = sample['mask'].long().cuda()
        name = sample['name'][0]

        with torch.no_grad():
            image = Variable(image)
            label = Variable(label)
            mask = Variable(mask)

            # The dimension of out should be in the dimension of B,C,H,W,D
            out = model(image)
            out = out.permute(0,2,3,4,1)
            out = torch.max(out, 4)[1]

            # start index corresponding to the model
            start = [14, 14, 14]
            center = out.size()[1:]
            prediction = torch.zeros(label.size())
            prediction[:, start[0]:start[0]+center[0], start[1]: start[1]+center[1], start[2]: start[2]+center[2]] = out

            # make the prediction corresponding to the center part of the image
            prediction = prediction.contiguous().view(-1).cuda()
            label = label.contiguous().view(-1)
            mask = mask.contiguous().view(-1)

        # compute the dice score
        score_per_image = accuracy(prediction.data.cpu().numpy(), mask.data.cpu().numpy(), label.data.cpu().numpy())

        if np.sum(score_per_image[0,:]) == 0 or np.sum(score_per_image[1,:]) == 0 or np.sum(score_per_image[2,:]) == 0:
            num_ignore += 1
            continue
        # compute the Enhance, Core and Whole dice score
        dice_score_per = [2 * np.sum(score_per_image[k,2]) / (np.sum(score_per_image[k,0]) + np.sum(score_per_image[k,1])) for k in range(3)]
        print('Image: %d, Enhance score: %.4f, Core score: %.4f, Whole score: %.4f' % (i, dice_score_per[0], dice_score_per[1], dice_score_per[2]))

        dice_score[i, :] = dice_score_per

        if args.visualize:
            vis = out.data.cpu().numpy()[0]
            vis = np.swapaxes(vis, 0, 2).astype(dtype=np.uint8)
            visualize_result(name, vis, args)

    count_image = num_images - num_ignore
    dice_score = dice_score[:count_image,:]
    mean_dice = np.mean(dice_score, axis=0)
    std_dice = np.std(dice_score, axis=0)
    print('Evalution Done!')
    print('Enhance score: %.4f, Core score: %.4f, Whole score: %.4f, Mean Dice score: %.4f' % (mean_dice[0], mean_dice[1], mean_dice[2], np.mean(mean_dice)))
    print('Enhance std: %.4f, Core std: %.4f, Whole std: %.4f, Mean Std: %.4f' % (std_dice[0], std_dice[1], std_dice[2], np.mean(std_dice)))

def main(args):
    # import network architecture
    builder = ModelBuilder()
    model = builder.build_net(
            arch=args.id,
            num_input=args.num_input,
            num_classes=args.num_classes,
            num_branches=args.num_branches,
            padding_list=args.padding_list,
            dilation_list=args.dilation_list)
    model = torch.nn.DataParallel(model, device_ids=list(range(args.num_gpus))).cuda()
    cudnn.benchmark = True

    if args.resume:
        if os.path.isfile(args.resume):
            print("=> Loading checkpoint '{}'".format(args.resume))
            checkpoint = torch.load(args.resume)
            state_dict = checkpoint['state_dict']
            model.load_state_dict(state_dict)
            print("=> Loaded checkpoint (epoch {})".format(checkpoint['epoch']))
        else:
            raise Exception("=> No checkpoint found at '{}'".format(args.resume))

    tf = ValDataset(test_dir, args)
    test_loader = DataLoader(tf, batch_size=args.batch_size, shuffle=args.shuffle, num_workers=args.num_workers, pin_memory=False)
    test(test_loader, model, args)

if __name__ == '__main__':
    assert LooseVersion(torch.__version__) >= LooseVersion('0.3.0'), \
        'PyTorch>=0.3.0 is required'

    parser = argparse.ArgumentParser()
    # Model related arguments
    parser.add_argument('--id', default='AFN1',
                        help='a name for identitying the model. Choose from the following options: AFN1-6, Basic, ASPP_c, ASPP_s.')
    parser.add_argument('--padding_list', default=[0,4,8,12], nargs='+', type=int,
                        help='list of the paddings in the parallel convolutions')
    parser.add_argument('--dilation_list', default=[2,6,10,14], nargs='+', type=int,
                        help='list of the dilation rates in the parallel convolutions')
    parser.add_argument('--num_branches', default=4, type=int,
                        help='the number of parallel convolutions in autofocus layer')

    # Path related arguments
    parser.add_argument('--test_path', default='datalist/test_list.txt',
                        help='txt file of the name of test data')
    parser.add_argument('--root_path', default='./',
                        help='root directory of data')
    parser.add_argument('--ckpt', default='./saved_models',
                        help='folder to output checkpoints')

    # Data related arguments
    parser.add_argument('--crop_size', default=[75,75,75], nargs='+', type=int,
                        help='crop size of the input image (int or list)')
    parser.add_argument('--center_size', default=[47,47,47], nargs='+', type=int,
                        help='the corresponding output size of the input image (int or list)')
    parser.add_argument('--num_classes', default=5, type=int,
                        help='number of classes')
    parser.add_argument('--num_input', default=5, type=int,
                        help='number of input image for each patient plus the mask')
    parser.add_argument('--num_workers', default=16, type=int,
                        help='number of data loading workers')
    parser.add_argument('--normalization', default=True, type=bool,
                        help='normalizae the data before running the test')
    parser.add_argument('--shuffle', default=False, type=bool,
                        help='if shuffle the data in test')
    parser.add_argument('--mask', default=True, type=bool,
                        help='if have the mask')

    # test related arguments
    parser.add_argument('--num_gpus', default=4, type=int, help='number of GPUs to use')
    parser.add_argument('--batch_size', default=1, type=int,
                        help='test batch size')
    parser.add_argument('--test_epoch', default=400, type=int,
                        help='epoch to start test.')
    parser.add_argument('--visualize', action='store_true',
                        help='save the prediction result as 3D images')
    parser.add_argument('--result', default='./result',
                        help='folder to output prediction results')
    parser.add_argument('--num_round', default=None, type=int,
                        help='restore the models from which run')

    args = parser.parse_args()
    print("Input arguments:")
    for key, value in vars(args).items():
        print("{:16} {}".format(key, value))

    test_file = open(args.test_path, 'r')
    test_dir = test_file.readlines()

    if not args.num_round:
        args.ckpt = os.path.join(args.ckpt, args.id)
    else:
        args.ckpt = os.path.join(args.ckpt, args.id, str(args.num_round))

    args.result = os.path.join(args.result, args.id)
    if not os.path.isdir(args.result):
        os.makedirs(args.result)

    assert len(args.padding_list) == args.num_branches, \
        '# parallel convolutions should be the same as the length of padding list'

    assert len(args.dilation_list) == args.num_branches, \
        '# parallel convolutions should be the same as # dilation rates'

    assert isinstance(args.crop_size, (int, list))
    if isinstance(args.crop_size, int):
        args.crop_size = [args.crop_size, args.crop_size, args.crop_size]

    assert isinstance(args.center_size, (int, list))
    if isinstance(args.center_size, int):
        args.center_size = [args.center_size, args.center_size, args.center_size]

    # do the test on a series of models
    args.resume = args.ckpt + '/' + str(args.test_epoch) + '_checkpoint.pth.tar'
    main(args)


text:
import json

import numpy as np
from fastapi import BackgroundTasks, FastAPI, HTTPException
from fastapi.middleware.cors import CORSMiddleware
from typing import List
import pandas as pd
from itertools import groupby
from pyjarowinkler import distance
from sklearn.cluster import AffinityPropagation
import os

from fastapi.encoders import jsonable_encoder
from fastapi.responses import JSONResponse

app = FastAPI()

origins = [
    "http://localhost:8080",
]

app.add_middleware(
    CORSMiddleware,
    allow_origins=["*"],
    allow_credentials=True,
    allow_methods=["*"],
    allow_headers=["*"],
)

#grouped_names = SqliteDict('./grouped_names.sqlite', encode=json.dumps, autocommit=True)
#intelligently_grouped_names = SqliteDict('./intelligently_grouped_names.sqlite', encode=json.dumps, autocommit=True)

grouped_names = dict()
intelligently_grouped_names = dict()

path = os.getcwd() + "/input"
input_file_path = f"{path}/names.csv"
input_names_df = pd.read_csv(input_file_path, index_col=None, header=None, names=['Name'])
names:List[str] = input_names_df['Name'].tolist()

@app.on_event("startup")
async def startup_event():
    names.sort()
    names_grouped_on_first_prefix = {j:tuple(i) for j, i in groupby(names,lambda a: a.split('_')[0])}

    for group_name, values in names_grouped_on_first_prefix.items():
        grouped_names[group_name] = values

def intelligent_clustering():
    names_np_array = np.asarray(names)
    jaro_winkler_similarity = -1*np.array([[distance.get_jaro_distance(w1, w2, winkler=True, scaling=0.1) for w1 in names_np_array] for w2 in names_np_array])
    affinity_propogation_model = AffinityPropagation(random_state=None)
    affinity_propogation_model.fit(jaro_winkler_similarity)

    for cluster_id in np.unique(affinity_propogation_model.labels_):
        cluster_name = names_np_array[affinity_propogation_model.cluster_centers_indices_[cluster_id]]
        cluster_values = np.unique(names_np_array[np.nonzero(affinity_propogation_model.labels_==cluster_id)])
        intelligently_grouped_names[cluster_name] = list(cluster_values)

@app.post("/perform-intelligent-grouping/")
async def intelligent_grouping(background_tasks: BackgroundTasks):
    background_tasks.add_task(intelligent_clustering)
    return {"message": "Affinity propagation based clustering in background"}

@app.get("/")
def intelligent_grouping_status():
    return {"Intellingent_cluster_status": len(intelligently_grouped_names)}

@app.get("/groups/{group_name}")
async def get_group_items(group_name: str):
    return grouped_names.get(group_name)

@app.get("/group")
async def get_all_groups():
    return grouped_names.items()

@app.get("/intelligent_group")
async def get_all_intelligent_groups():
    return intelligently_grouped_names.items()

@app.get("/names/{name}")
async def get_group_name(name: str):
    group_name = [grp_name for grp_name, list_of_names in grouped_names.items() if any(name in individual_name for individual_name in list_of_names)]
    if len(group_name) == 0:
        raise HTTPException(status_code=404, detail="Name not found in any group")
    return group_name