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

import sys
import time
from scipy import append
#from math import log
from scipy.stats import multivariate_normal,multinomial
import matplotlib.pyplot as plt
#import collections
#import spconavi_path_calculate
import spconavi_read_data
import spconavi_save_data
from __init__ import *
from submodules import *

##Command:
#python3 spcotmhp_dijkstra_execute.py trialname iteration sample init_position_num speech_num initial_position_x initial_position_y waypoint_word
#python3 spcotmhp_dijkstra_execute.py 3LDK_01 1 0 -1 7 100 100 -1

tools     = spconavi_read_data.Tools()
read_data = spconavi_read_data.ReadingData()
save_data = spconavi_save_data.SavingData()
#path_calculate = spconavi_path_calculate.PathPlanner()

#Definition of action (functions in spconavi_read_data)
action_functions = [tools.right, tools.left, tools.up, tools.down, tools.stay] #, migiue, hidariue, migisita, hidarisita]
cost_of_actions  = np.log( np.ones(len(action_functions)) / float(len(action_functions)) ) #[    1/5,    1/5,  1/5,    1/5,    1/5]) #, ,    1,        1,        1,          1]

#GaussMap make (same to the function in spcotmhp_viterbi_path_calculate) #Ito
def PostProbMap_Gauss(CostMapProb,Mu,Sig,map_length,map_width,it): #,IndexMap):
        x,y = np.meshgrid(np.linspace(-10.0,9.92,map_width),np.linspace(-10.0,9.92,map_length))
        pos = np.dstack((x,y))
        #PostProbMap = np.array([ [ PostProb_ij([width, length],Mu,Sig,map_length,map_width, CostMapProb,it) for width in xrange(map_width) ] for length in xrange(map_length) ])
        PostProb=multivariate_normal(Mu[it],Sig[it]).pdf(pos)

        return CostMapProb * PostProb

# The weight of the topoligical-level likelihood in word Se
def TopoWeight(i,j,Connect,w,THETA,psi):
    #THETAを展開
    W, W_index, Myu, S, pi, phi, K, L = THETA

    if (Connect == 1):
        # i_e: i, i_e-1: j, S_e:w
        sum_c  = np.sum([phi[c][i]*(W[c][w]**N_best)*pi[c] for c in range(L)])
        bunbo2 = np.sum([phi[c][i] for c in range(L)])
        weight = psi[i][j] * sum_c / bunbo2
    else:
      weight = 1.0

    return weight

###↓### Sampling of goal candidates ############################################
def EstimateGoal(Otb_B, THETA):
    #THETAを展開
    W, W_index, Mu, S, pi, phi_l, K, L = THETA

    #Prob math func of p(it | Otb_B, THETA) = Σc p(it | phi_c)p(st=Otb_B | Wc)p(c | pi)
    pmf_it = np.ones(K)
    for k in range(K):
        pmf_it[k] = np.sum([ phi_l[c][k] * multinomial.pmf(Otb_B, sum(Otb_B), W[c]) * pi[c] for c in range(L) ])
        bunbo2 = np.sum([phi_l[c][k] for c in range(L)])
        pmf_it[k] = pmf_it[k] / bunbo2

    #argmax_it = np.argmax(pmf_it)

    #Normalization
    pmf_it_n = np.array([pmf_it[i] / float(np.sum(pmf_it)) for i in range(K)])

    #bests_five  = []
    upper_zero_three = np.array([1.0 for i in range(K)])
    for i in range(K):
        if (pmf_it_n[i] < 0.3):
            upper_zero_three[i] = 0.0

    pmf_it = pmf_it * upper_zero_three
    #Normalization2
    pmf_it_n = np.array([pmf_it[i] / float(np.sum(pmf_it)) for i in range(K)])

    #Sampling it from multinomial distribution
    sample_it = np.random.choice(K, size=Sampling_J, replace=True, p=pmf_it_n)
    #random.choices([i for i in range(K)], k=Sampling_J, weights=pmf_it_n)
    #multinomial.rvs(Sampling_J, pmf_it_n, size=1, random_state=None)
    #print(sample_it)

    return sample_it

    """
    goal_candidate = []
    for it in range(K):
        count_it = 0
        while (count_it < sample_it[0][it]):
            goal_candidate += [tools.Map_coordinates_To_Array_index(multivariate_normal.rvs(mean=Myu[it], cov=S[it], size=1, random_state=None))]
            count_it += 1
    #Xt_max = Map_coordinates_To_Array_index( [ Xp[pox.index(max(pox))][0], Xp[pox.index(max(pox))][1] ] ) #[0.0,0.0] ##確率最大の座標候補
    goal_candidate_tuple = [(goal_candidate[j][1], goal_candidate[j][0]) for j in range(Sampling_J)]
    print("Goal candidates:", goal_candidate_tuple)
    return goal_candidate_tuple
    """

###↑### Sampling of goal candidates ############################################

### A star algorithm (by Ryo Ozaki) ############################################
def a_star(start, goal, maze, action_functions, cost_of_actions, PathWeightMap):
    if (maze[goal[0]][goal[1]] != 0):
        print("[ERROR] goal",maze[goal[0]][goal[1]],"is not 0.")

    ###START A*
    open_list = []
    open_list_cost = []
    open_list_key = []
    closed_list = []
    closed_list_cost = []
    closed_list_key = []
    open_list.append(start)
    open_list_cost.append(0)
    open_list_key.append(0 + tools.Manhattan_distance(start, goal))
    OYA = {}
    ko = (0), (0)
    Path = []

    while open_list:
        sorted_idx = np.argsort(open_list_key, kind="stable")
        pop_idx = sorted_idx[0]
        p = open_list.pop(pop_idx)
        p_cost = open_list_cost.pop(pop_idx)
        p_key = open_list_key.pop(pop_idx)
        closed_list.append(p)
        closed_list_cost.append(p_cost)
        closed_list_key.append(p_key)
        if p == goal:
            break
        for act_func, act_cost in zip(action_functions, cost_of_actions):
            q = act_func(p)
            if (int(maze[q]) != 0):
                continue
            q_cost = p_cost - act_cost - np.log(PathWeightMap[q[0]][q[1]])  #current sum cost and action cost
            q_pev = tools.Manhattan_distance(q, goal) * np.log(float(len(action_functions))) #heuristic function
            q_key = q_cost - q_pev

            if q in open_list:
                idx = open_list.index(q)
                key = open_list_key[idx]
                if key > q_key:
                    open_list_key[idx] = q_key
                    open_list_cost[idx] = q_cost
            elif q in closed_list:
                idx = closed_list.index(q)
                key = closed_list_key[idx]
                if key > q_key:
                    closed_list.pop(idx)
                    closed_list_cost.pop(idx)
                    closed_list_key.pop(idx)
                    open_list.append(q)
                    open_list_cost.append(q_cost)
                    open_list_key.append(q_key)
                    #plt.quiver(p[1], p[0], (q[1]-p[1]), (q[0]-p[0]), angles='xy', scale_units='xy', scale=1, color="tab:red")
                    OYA[(q[1], q[0])] = (p[1], p[0])
                    ko = (q[1]), (q[0])
                    #print(ko)
            else:
                open_list.append(q)
                open_list_cost.append(q_cost)
                open_list_key.append(q_key)
                #plt.quiver(p[1], p[0], (q[1]-p[1]), (q[0]-p[0]), angles='xy', scale_units='xy', scale=1, color="tab:red")
                OYA[(q[1], q[0])] = (p[1], p[0])
                ko = (q[1]), (q[0])
                #print(ko)

    #最適経路の決定: ゴールから親ノード（どこから来たか）を順次たどっていく
    #i = len(OYA)
    #for oyako in reversed(OYA):
    ko_origin = ko
    ko = (goal[1], goal[0])
    print(ko,goal)
    #for i in range(p_cost):
    while(ko != (start[1],start[0])):
        #print(OYA[ko])
        try:
            Path = Path + [OYA[ko]]
        except KeyError:
            ko = ko_origin
            Path = Path + [OYA[ko]]
            print("NOT END GOAL.")

        ko = OYA[ko]
        #i = len(Path)
        #print(i, ko)
        #i -= 1

    print(goal,": Total cost using A* algorithm is "+ str(p_cost))
    return Path, p_cost
### A star algorithm (by Ryo Ozaki) ############################################

def ReadDistanceMtrx(outputfile):
    dis  = [ [0.0 for atem in range(K)] for aky in range(K) ]
    c=0
    for line in open(outputfile + 'Astar_SpCoTMHP_distance.csv', 'r'):
        itemList = line[:-1].split(',')
        for i in range(len(itemList)):
            if itemList[i] != "":
              dis[c][i] = float(itemList[i])
        c = c + 1
    return dis

def ReadCostMtrx(outputfile):
    cost  = [ [0.0 for atem in range(K)] for aky in range(K) ]
    c=0
    for line in open(outputfile + 'Astar_SpCoTMHP_cost.csv', 'r'):
        itemList = line[:-1].split(',')
        for i in range(len(itemList)):
            if itemList[i] != "":
              cost[c][i] = float(itemList[i])
        c = c + 1
    return cost

def ReadPathDistance(outputsubfolder,st,gl):
    for line in open(outputsubfolder + "Astar_SpCoTMHP_S" + str(st) + "_G" + str(gl) + "_Distance.csv", 'r'):
        distance = float(line[:-1])

    return int(distance)

def ReadLogLikelihood(outputsubfolder,st,gl):
    LogLikelihood  = [ 0.0 ]
    c=0
    for line in open(outputsubfolder + "Astar_SpCoTMHP_S" + str(st) + "_G" + str(gl) + "_Log_likelihood_sum.csv", 'r'):
        #itemList = line[:-1].split(',')
        #for i in range(len(itemList)):
        #    if itemList[i] != "":
        LogLikelihood.append(float(line[:-1]))
        c = c + 1
    dist = int(ReadPathDistance(outputsubfolder,st,gl))

    return LogLikelihood[dist]

def ReadLogLkhdMtrx(outputsubfolder,CoonectMatricx):
    LogLkhdMtrx  = np.array([ [ 0.0 for gl in range(K)] for st in range(K) ])
    for st in range(K):
        for gl in range(K):
            if (CoonectMatricx[st][gl] == 1) and (st != gl):
                LogLkhdMtrx[st][gl] = ReadLogLikelihood(outputsubfolder,st,gl)
    return LogLkhdMtrx

# おそらく最適探索（Dijkstra）が実装されている
def SearchTopoDijkstra(glf, st, cost): # goal index, start index, cost
    ev   = [ 0.0 for aky in range(K) ]
    node = [ [0,-2] for aky in range(K) ]
    node[st]=[1,-1]

    fin=0
    while fin != 1:
      for i in range(K):
        if node[i][0]==2:
            for k in range(K):
                if cost[i][k]>0:
                    if node[k][0]==0:
                        node[k][0]=1
                        node[k][1]=i
                        ev[k]=ev[i]+cost[i][k]
                    elif node[k][0]==1:
                        if ev[k]>ev[i]+cost[i][k]:
                            ev[k]=ev[i]+cost[i][k]
                            node[k][1]=i
            node[i][0]=3
            break

      minm=100000000000000000
      minm_i=-1
      for i in range(K):
        if node[i][0]==1:
          if minm > ev[i]:
            minm=ev[i]
            minm_i=i
      if minm_i == glf:
         fin=1
         #print("fin")
      else:
       node[minm_i][0]=2

    i=glf
    root=[]
    while i!= -1 :
       for k in range(K):
         if i==k:
           root.insert(0,k)
           i=node[k][1]
    return root

# おそらく動的計画法（Viterbi）が実装されている
def SearchTopoViterbi(s_i, lk, Step): # start index, likelihoods, T_topo
    lkt=[ [0.0 for i in range(K)] for k in range(Step)]#sum likewood
    rt=[ [-1 for i in range(K)] for k in range(Step)]#oya node save
    node =[0 for i in range(K) ] #
    node[s_i]=1
    ndc=[0 for i in range(K) ]#
    root=[0 for i in range(Step) ]
    mxs=-1#
    mxg=-1
    mx=0
    #print(lk)
    #print(node)

    for t in range(Step):
      for i in range(K):
        if t==0:
          if node[i]==1:
            for j in range(K):
              if lk[i][j]>0:
                 lkt[t][j]=lk[i][j]
                 rt[t][j]=i
                 ndc[j]=1
        else:
          if node[i]==1:
            for j in range(K):
              if lk[i][j]>0:
                if lkt[t][j]<lkt[t-1][i]+lk[i][j]:
                  lkt[t][j]=lkt[t-1][i]+lk[i][j]
                  rt[t][j]=i
                  ndc[j]=1

      for f in range(K):
         if ndc[f]==1 and node[f]==0:
           node[f]=1

    #print(rt)
    for i in range(K):
        if mx<lkt[Step-1][i]:
          mx=lkt[Step-1][i]
          mxg=i
    #print(mxg)

    # Backword path
    for i in range(Step):
        root[Step-i-1]=rt[Step-i-1][mxg]
        mxg=rt[Step-i-1][mxg]

    return root

## パスを読み込み、つなげる (For A*)
def PathConnect(outputsubfolder, dis, Path_ROS, root):
    for v in range(len(root)-1):
        path=[ (0,0) for k in range(int(dis[root[v]][root[v+1]]))]
        #print(dis[root[v]][root[v+1]])
        i=0
        if ( root[v] != root[v+1] ):
            for line in open(outputsubfolder + "Astar_SpCoTMHP_S" + str(root[v]) + '_G' + str(root[v+1]) +  '_Path_ROS.csv', 'r'):
                itemList = line[:-1].split(',')
                path[i] = ( float(itemList[0]) , float(itemList[1]) )
                #path[i][0]=int(po[0])
                #path[i][1]=int(po[1])
                i = i + 1
                #print(i)
            Path_ROS=Path_ROS+path
    return Path_ROS

#################################################
if __name__ == '__main__':
    print("[START] SpCoTMHP (Dijkstra).")

    #Request a folder name for learned parameters.
    trialname = sys.argv[1]

    #Request iteration value
    iteration = sys.argv[2] #1

    #Request sample value
    sample = sys.argv[3] #0

    #Request the index number of the robot initial position
    init_position_num = int(sys.argv[4]) #位置分布インデックス(-1で座標直接指定)

    #Request the file number of the speech instruction
    speech_num = sys.argv[5] #0

    #初期値を指定する場合はTHETA読み込み以降に記載(int(sys.argv[6]) , int(sys.argv[7]))

    #中間地点の単語番号を指定 (未実装：複数指定の場合、コンマ区切りする)
    waypoint_word = sys.argv[8] # -1:中間なし

    WP_list    = waypoint_word[:].split(',')
    print("WP:", WP_list)

    if (int(WP_list[0]) == -1):
        tyukan = 0
    else:
        tyukan = 1

    ##FullPath of folder
    filename = outputfolder_SIG + trialname #+ "/"
    print(filename, iteration, sample)
    outputfile      = filename + navigation_folder #outputfolder + trialname + navigation_folder
    outputsubfolder = outputfile + "astar_node_gauss/"
    outputname_d    = outputfile + "dijkstra_result_wd/"

    Makedir( outputname_d )

    #Read the files of learned parameters  #THETA = [W,W_index,Mu,Sig,Pi,Phi_l,K,L]
    THETA   = read_data.ReadParameters(iteration, sample, filename, trialname)
    W, W_index, Mu, Sig, pi, phi_l, K, L = THETA
    psi   = read_data.ReadPsi(iteration, sample, filename, trialname)

    gridmap = read_data.ReadMap(outputfile)
    map_length, map_width = gridmap.shape

    CostMapProb = read_data.ReadCostMapProb(outputfile)
    #CostMap = read_data.ReadCostMap(outputfile)
    #CostMapProb = (100.0 - CostMap)/100

    #描画
    plt.imshow(gridmap + (40+1)*(gridmap == -1), origin='lower', cmap='binary', vmin = 0, vmax = 100, interpolation='none') #, vmin = 0.0, vmax = 1.0)
    plt.xticks(rotation=90)
    plt.tick_params(axis='x', which='major', labelsize=8)
    plt.tick_params(axis='y', which='major', labelsize=8)
    #plt.xlim([380,800])             #x軸の範囲
    #plt.ylim([180,510])             #y軸の範囲
    plt.xlabel('X', fontsize=10)
    plt.ylabel('Y', fontsize=10)
    #plt.xticks(np.arange(width), np.arange(width))
    #plt.yticks(np.arange(height), np.arange(height))
    plt.gca().set_aspect('equal')

    # 初期位置の設定
    if (init_position_num == -1): #初期値を指定する場合
        start = [int(sys.argv[6]) , int(sys.argv[7])]
        start_inv = [start[1], start[0]]
        print("Start:", start)
    else: #初期値は__init__.pyのリストから選択する場合
        start_inv = Start_Position[init_position_num]
        start = [start_inv[1], start_inv[0]]

    outputname      = outputname_d + "Dijkstra_SpCoTMHP_"+"T"+str(T_topo)+"S"+str(start)+"H"+str(waypoint_word)+"G"+str(speech_num)

    dis  = ReadDistanceMtrx(outputfile)
    cost = ReadCostMtrx(outputfile)
    #lk   = -1.0 * np.array(cost)  # cost is -log likelihood

    #Read CoonectMatricx (psi_setting.csv) file
    CoonectMatricx=[ [0.0 for i in range(K)] for c in range(K) ]
    c = 0
    for line in open(filename + "/" + trialname + '_psi_setting.csv', 'r'):
        itemList = line[:-1].split(',')
        for i in range(len(itemList)):
            if itemList[i] != "":
              CoonectMatricx[c][i] = float(itemList[i])
        c = c + 1

    lk = ReadLogLkhdMtrx(outputsubfolder,CoonectMatricx) # log likelihood

    ### 事前計算として、W_indexごとにトポロジカル重みを計算
    TopoWeight_W_index = [ np.log(np.array([ [TopoWeight(i,j,CoonectMatricx[i][j],w,THETA,psi) for i in range(K)] for j in range(K) ])) for w in range(len(W_index)) ]

    same_lk = np.array([multivariate_normal(Mu[it],Sig[it]).pdf(Mu[it])*0.5 for it in range(K)])
    #same_lk = np.array([CostMapProb[tools.Map_coordinates_To_Array_index(Mu[it])[0]][tools.Map_coordinates_To_Array_index(Mu[it])[1]] for it in range(K)])
    same_lk = np.diag(np.log(same_lk))
    #print("same_lk",same_lk)

    duration_term_lk = np.array([ [-1.0 * dis[i][j] * np.log(T_topo) for j in range(K)] for i in range(K) ]) # log duration prob

    whole_log_lk_inv_W_index =  -1.0 * np.array([TopoWeight_W_index[w] + lk + same_lk + duration_term_lk for w in range(len(W_index))])

    #weight_inv = -1.0 * (TopoWeight_W_index[S_list[t]] + lk + np.diag(same_lk))
    weight_W_index = whole_log_lk_inv_W_index #[ whole_log_lk_inv_W_index[w] - np.min(whole_log_lk_inv_W_index[w]) for w in range(len(W_index))]
    #[np.max(whole_log_lk_inv_W_index[w]) - whole_log_lk_inv_W_index[w] for w in range(len(W_index))]

    if (SAVE_time == 1):
      #Substitution of start time
      start_time = time.time()

    ### 初期位置から一番近くの位置分布を推定（ガウスの重みが一番高い分布のインデックス）
    near_node = np.argmax([multivariate_normal(Mu[k],Sig[k]).pdf(tools.Array_index_To_Map_coordinates(start)) for k in range(K)])

    #ss=Map_coordinates_To_Array_index(Mu[S_i])
    #start=(150,130)

    ###v###近くの位置分布のindexのガウスの平均をゴールとしたA*を実行###v###
    PathWeightMap = PostProbMap_Gauss(CostMapProb,Mu,Sig,map_length,map_width,near_node)

    gg=tools.Map_coordinates_To_Array_index(Mu[near_node])
    #print(gg[0])
    goal_candidate=[[gg[1],gg[0]]]
    #goal_candidate=[[55,135],[89,103]]
    #J = len(goal_candidate)
    J=1
    #if(J != THETA[6]):
    # print("[WARNING] J is not K",J,K)
    p_cost_candidate = [0.0 for j in range(J)]
    Path_candidate   = [[0.0] for j in range(J)]
    #print(goal_candidate)

    ###goal候補ごとにA*を実行
    for gc_index in range(J):
        goal = goal_candidate[gc_index]
        Path, p_cost = a_star(start_inv, goal, gridmap, action_functions, cost_of_actions, PathWeightMap)

        #first_cost = p_cost / float(len(Path))
        p_cost_candidate[gc_index] = p_cost / float(len(Path)+1)
        Path_candidate[gc_index] = Path

    ### select the goal of expected cost
    expect_gc_index = np.argmin(p_cost_candidate)
    Path = Path_candidate[expect_gc_index]
    goal = goal_candidate[expect_gc_index]
    print("Goal:", goal)

    #The moving distance of the path
    Distance = tools.PathDistance(Path)
    print("Path distance using A* algorithm is "+ str(Distance))

    #計算上パスのx,yが逆になっているので直す
    Path_inv = [ (Path[t][1], Path[t][0]) for t in range(len(Path)) ]
    Path_inv.reverse()
    Path_ROS = Path_inv #使わないので暫定的な措置
    #パスを保存
    #SavePath(start, [goal[1], goal[0]], Path_inv, Path_ROS, outputname)
    ###^###近くの位置分布のindexのガウスの平均をゴールとしたA*を実行###^###

    #####Estimate the goal point by spatial concept
    Otb_B      = [int(W_index[i] == Goal_Word[int(speech_num)]) * N_best for i in range(len(W_index))]
    print("BoW:", Otb_B)

    Otb_B_list = [ [int(W_index[i] == Goal_Word[int(WP_list[w])]) * N_best for i in range(len(W_index))] for w in range(len(WP_list))]
    print("WP:BoW_list:", Otb_B_list)
    len_Otb_B_list = len(Otb_B_list)

    #st = near_node
    gl_J_best = [0.0 for j in range(Sampling_J)]
    root_list = [ [] for j in range(Sampling_J)]
    Path_list = [Path_ROS for j in range(Sampling_J)]
    Path_log_likelihood_list = [0.0 for j in range(Sampling_J)]
    S_list = [] # W_indexの番号の配列
    #i = 0

    ## 最終ゴール（位置分布インデックス）の候補リストを推定
    gl_J_best = EstimateGoal(Otb_B, THETA)
    print("gl:",gl_J_best)

    for j in range(Sampling_J):
        i=0
        #tyukan = 0
        st = near_node
        if (tyukan == 1):
            for i in range(len(Otb_B_list)):
                ### 中間地点が指定されている場合、中間地点（がなくなる）までパスプランニング、中間地点がなければゴールへのプランニング ###
                #for j_tyukan in range(Sampling_J):
                #print("tyukan",i)
                ## 次のゴール（位置分布インデックス）を推定
                tyukan_J_best = EstimateGoal(Otb_B_list[i], THETA)
                print("st,tyukan:",st,tyukan_J_best[j])

                ## 次のゴールまでのトポロジカルプランニング
                root = SearchTopoDijkstra(tyukan_J_best[j], st, weight_W_index[np.argmax(Otb_B_list[i])])
                root_list[j] = root_list[j] + root
                print("root_t:",j,root_list[j])

                ## パスの尤度を計算し保存
                for r in range(len(root)):
                    if r == 0:
                        Path_log_likelihood_list[j] += weight_W_index[np.argmax(Otb_B_list[i])][st][root[r]]
                    else:
                        Path_log_likelihood_list[j] += weight_W_index[np.argmax(Otb_B_list[i])][root[r-1]][root[r]]

                ## パスを読み込み、つなげる
                Path_list[j] = PathConnect(outputsubfolder, dis, Path_list[j], root)

                root = root_list[j]
                ## T_topo/(len_Otb_B_list+1))まで最終地点の値を埋める
                for _ in range(int(T_topo/(len_Otb_B_list+1) - len(root))):
                    Path_log_likelihood_list[j] += weight_W_index[np.argmax(Otb_B)][root[-1]][root[-1]]

                st = tyukan_J_best[j]
                #i = i + 1

            ## 中間地点がなければゴールへのプランニングへ
            #if (i >= len(Otb_B_list)):
            #    tyukan = 0

        ## 次のゴールまでのトポロジカルプランニング
        root = SearchTopoDijkstra(gl_J_best[j], st, weight_W_index[np.argmax(Otb_B)])
        root_list[j] = root_list[j] + root
        print("root_j:",j,root_list[j])

        ## パスを読み込み、つなげる
        Path_list[j] = PathConnect(outputsubfolder, dis, Path_list[j], root)

        root = root_list[j]
        ## パスの尤度を計算し保存
        for r in range(len(root)):
            if r == 0:
                Path_log_likelihood_list[j] += weight_W_index[np.argmax(Otb_B)][st][root[r]]
            else:
                Path_log_likelihood_list[j] += weight_W_index[np.argmax(Otb_B)][root[r-1]][root[r]]

        ## T_topoまで最終地点の値を埋める
        for _ in range(T_topo - len(root)):
            Path_log_likelihood_list[j] += weight_W_index[np.argmax(Otb_B)][root[-1]][root[-1]]

    ## 候補のパスから尤度が最大(-log likelihoodなので最小)のものを選ぶ
    saiyou = np.argmin(Path_log_likelihood_list)
    print("saiyou",saiyou,root_list[saiyou])
    Path_ROS = Path_list[saiyou]

    if (SAVE_time == 1):
        #PP終了時刻を保持
        end_pp_time = time.time()
        time_pp = end_pp_time - start_time #end_recog_time
        fp = open( outputname + "_time_pp.txt", 'w')
        fp.write(str(time_pp)+"\n")
        fp.close()

    print("[END] SpCoTMHP (Dijkstra).")
    # print(Path_ROS)

    disp = [tools.PathDistance(Path_ROS)] # tupple in list  #len(Path_ROS)

    for i in range(len(Path_ROS)):
      plt.plot(Path_ROS[i][1], Path_ROS[i][0], "s", color="tab:red", markersize=1)
    plt.savefig(outputname + '_Path.pdf', dpi=300)#, transparent=True
    plt.savefig(outputname + '_Path.png', dpi=300)#, transparent=True

    np.savetxt(outputname + "_fin_Path_ROS.csv", Path_ROS, delimiter=",")
    np.savetxt(outputname + "_fin_Distance.csv", disp, delimiter=",")
    plt.clf()
    print("DONE.",outputname)

max_stars_repo_name: a-taniguchi/SpCoTMHP

max_issues_repo_name: a-taniguchi/SpCoTMHP

max_forks_repo_name: a-taniguchi/SpCoTMHP

reason: from scipy


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

"""
from __future__ import division, print_function, unicode_literals, absolute_import

#External Modules------------------------------------------------------------------------------------
import abc
import os
import numpy as np
#External Modules End--------------------------------------------------------------------------------

#Internal Modules------------------------------------------------------------------------------------
from utils.cached_ndarray import c1darray
from utils import utils, InputData
import MessageHandler
from utils import InputData, InputTypes
#Internal Modules End--------------------------------------------------------------------------------

class CheckInterfacePP(InputData.CheckClass):
  """
  Checks that this is an Interface Post Processor of a given type
  """
  def __init__(self, name):
    """
      Creates a CheckInterfacePP class
      @ In, name, string, the method name
      @ Out, None
    """
    self.name = name
    self.__reason = ""

  def check(self, node):
    """
      Checks the node to see if it matches the checkDict
      @ In, node, xml node to check
      @ Out, bool, true if matches
    """
    self.__reason = ""
    passed = "subType" in node.attrib and node.attrib["subType"] == "InterfacedPostProcessor"
    if not passed:
      self.__reason = "subType=InterfacedPostProcessor not in attribs"
    methods = node.findall("method")
    if len(methods) == 1:
      match = methods[0].text == self.name
      if not match:
        self.__reason += ""+repr(methods[0].text)+ "!=" + self.name + " "
      passed = passed and match
    else:
      self.__reason += "wrong number of method blocks "+str(len(methods))+" "
      passed = False
    return passed

  def failCheckReason(self, node):
    """
      returns a string about why the check failed
      @ In, node, xml node to check
      @ Out, string, message for user about why check failed.
    """
    return self.__reason

class PostProcessorInterfaceBase(utils.metaclass_insert(abc.ABCMeta,object),MessageHandler.MessageUser):
  """
    This class is the base interfaced post-processor class
    It contains the three methods that need to be implemented:
      - initialize
      - run
      - readMoreXML
  """

  @classmethod
  def getInputSpecification(cls):
    """
      Method to get a reference to a class that specifies the input data for
      class cls.
      @ In, cls, the class for which we are retrieving the specification
      @ Out, inputSpecification, InputData.ParameterInput, class to use for
        specifying input of cls.
    """
    inputSpecification = InputData.parameterInputFactory(cls.__name__, ordered=False)
    inputSpecification.setCheckClass(CheckInterfacePP("PostProcessorInterfaceBaseClass"))
    inputSpecification.addParam("subType", InputTypes.StringType)
    inputSpecification.addParam("name", InputTypes.StringType)

    return inputSpecification

  def __init__(self, messageHandler):
    """
      Constructor
      @ In, messageHandler, MessageHandler, message handler object
      @ Out, None
    """
    self.type = self.__class__.__name__
    self.name = self.__class__.__name__
    self.messageHandler = messageHandler

    self.transformationSettings = {}   # this dictionary is used to store all the setting required to back transform the data into its original format
                                       # it gets filled in the run method and used in the inverse method

  def initialize(self):
    """
      Method to initialize the Interfaced Post-processor. Note that the user needs to specify two mandatory variables:
       - self.inputFormat:  dataObject that the PP is supposed to receive in input
       - self.outputFormat: dataObject that the PP is supposed to generate in output
      These two variables check that the input and output dictionaries match what PP is supposed to receive and generate
      Refer to the manual on the format of these two dictionaries
      @ In, None
      @ Out, None
    """
    self.inputFormat  = None
    self.outputFormat = None

  def readMoreXML(self,xmlNode):
    """
      Function that reads elements this post-processor will use
      @ In, xmlNode, ElementTree, Xml element node
      @ Out, None
    """
    pass

  def run(self,inputDic):
    """
     Method to post-process the dataObjects
     @ In, inputDic, dict, dictionary which contains the data inside the input DataObject
     @ Out, None
    """
    pass

  def _inverse(self,inputDic):
    """
      Method to perform the inverse of the post-process action
      @ In, inputDic, dict, dictionary which contains the data to be back pre-processed
      @ Out, None
    """
    pass

  def checkGeneratedDicts(self,outputDic):
    """
      Method to check that dictionary generated in def run(self, inputDic) is consistent
      @ In, outputDic, dict, dictionary generated by the run method
      @ Out, True/False, bool, outcome of the outputDic check
    """
    checkInp = self.checkInputFormat(outputDic['data']['input'])
    checkOut = self.checkOutputFormat(outputDic['data']['output'])
    if checkInp and checkOut:
      return True
    else:
      if not checkInp:
        self.raiseAWarning('PP Generation check on Inputs failed!')
      if not checkOut:
        self.raiseAWarning('PP Generation check on Outputs failed!')
      return False

  def checkOutputFormat(self,outputDic):
    """
      This method checks that the generated output part of the generated dictionary is built accordingly to outputFormat
      @ In, outputDic, dict, dictionary generated by the run method
      @ Out, outcome, bool, outcome of the outputDic check (True/False)
    """
    outcome = True
    if isinstance(outputDic,dict):
      if self.outputFormat == 'HistorySet':
        for key in outputDic:
          if isinstance(outputDic[key],dict):
            outcome = outcome and True
          else:
            self.raiseAWarning('Bad PP output type for key:',key,':',type(outputDic[key]),'; should be dict!')
            outcome = False
          for keys in outputDic[key]:
            if isinstance(outputDic[key][keys],(np.ndarray,c1darray)):
              outcome = outcome and True
            else:
              self.raiseAWarning('Bad PP output type for key:',key,keys,':',type(outputDic[key][keys]),'; should be np.ndarray or c1darray!')
              outcome = False
      else:  # self.outputFormat == 'PointSet':
        for key in outputDic:
          if isinstance(outputDic[key],(np.ndarray,c1darray)):
            outcome = outcome and True
          else:
            self.raiseAWarning('Bad PP output type for key:',key,':',type(outputDic[key]),'; should be np.ndarray or c1darray!')
            outcome = False
    else:
      self.raiseAWarning('Bad PP output dict:',type(outputDic),'is not a dict!')
      outcome = False
    return outcome

  def checkInputFormat(self,outputDic):
    """
      This method checks that the generated input part of the generated dictionary is built accordingly to outputFormat
      @ In, outputDic, dict, dictionary generated by the run method
      @ Out, outcome, bool, outcome of the outputDic check (True/False)
    """
    outcome = True
    if isinstance(outputDic,dict):
        for key in outputDic:
          if isinstance(outputDic[key],(np.ndarray,c1darray)):
            outcome = outcome and True
          else:
            self.raiseAWarning('Bad PP output type for key:',key,':',type(outputDic[key]),'; should be np.ndarray or c1darray!')
            outcome = False
    else:
      self.raiseAWarning('Bad PP output dict:',type(outputDic),'is not a dict!')
      outcome = False
    return outcome

  def checkArrayMonotonicity(time):
    """
      This method checks that an array is increasing monotonically
      @ In, time, numpy array, array to be checked
      @ Out, outcome, bool, outcome of the monotonicity check
    """
    outcome = True
    for t in time:
      if t != 0:
        if time[t] > time[t-1]:
          outcome = outcome and True
        else:
          outcome = outcome and False
    return outcome

ext: py

max_stars_repo_name: JiaZhou-PU/raven

max_issues_repo_name: JiaZhou-PU/raven

max_forks_repo_name: JiaZhou-PU/raven

reason: import numpy


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

import collections
import numpy as np
import sonnet as snt
import tensorflow as tf

DatasetTensors = collections.namedtuple('DatasetTensors', ('input',
                                                           'target', 'mask'))

"""
Model for Figure:
python3 -m src.tasks.train --controller=ff --num_read_heads=1 --memory_size=10
    --word_size=3 --batch_size=10 --max_repeats=1 --gpu_usage=0.8 --num_bits=6
    --min_length=5 --max_length=5
    --checkpoint_dir=src/tasks/repeat_copy/checkpoints
    --checkpoint_basename=model-word-size-3.ckpt --checkpoint_interval=10000
    --num_training_iterations=100000
"""

def masked_sigmoid_cross_entropy(logits,
                                 target,
                                 mask,
                                 time_average=False,
                                 log_prob_in_bits=False):
    """Add ops to graph which compute the (scalar) NLL of the target sequence.

    The logits parametrize independent bernoulli distributions per time-step
    and per batch element, and irrelevant time/batch elements are masked out by
    the mask tensor.

    Args:
        logits: `Tensor` of activations for which sigmoid(`logits`) gives the
            bernoulli parameter.
        target: time-major `Tensor` of target.
        mask: time-major `Tensor` to be multiplied elementwise with cost T x B
            cost masking out irrelevant time-steps.
        time_average: optionally average over the time dimension (sum by
            default).
        log_prob_in_bits: iff True express log-probabilities in bits (default
            nats).
    Returns:
        A `Tensor` representing the log-probability of the target.
    """
    xent = tf.nn.sigmoid_cross_entropy_with_logits(labels=target,
                                                   logits=logits)
    loss_time_batch = tf.reduce_sum(xent, axis=2)
    loss_batch = tf.reduce_sum(loss_time_batch * mask, axis=0)

    batch_size = tf.cast(tf.shape(logits)[1], dtype=loss_time_batch.dtype)

    if time_average:
        mask_count = tf.reduce_sum(mask, axis=0)
        loss_batch /= (mask_count + np.finfo(np.float32).eps)

    loss = tf.reduce_sum(loss_batch) / batch_size
    if log_prob_in_bits:
        loss /= tf.log(2.)

    return loss

def bitstring_readable(data, batch_size, model_output=None, whole_batch=False):
    """Produce a human readable representation of the sequences in data.

    Args:
        data: data to be visualised
        batch_size: size of batch
        model_output: optional model output tensor to visualize alongside data.
        whole_batch: whether to visualise the whole batch. Only the first
            sample will be visualized if False

    Returns:
        A string used to visualise the data batch
    """
    def _readable(datum):
        return '+' + ' '.join(['-' if x == 0 else '%d' % x for x in datum]) + \
            '+'

    obs_batch = data.input
    targ_batch = data.target

    iterate_over = range(batch_size) if whole_batch else range(1)

    batch_strings = []
    for batch_index in iterate_over:
        obs = obs_batch[:, batch_index, :]
        targ = targ_batch[:, batch_index, :]

        obs_channels = range(obs.shape[1])
        targ_channels = range(targ.shape[1])
        obs_channel_strings = [_readable(obs[:, i]) for i in obs_channels]
        targ_channel_strings = [_readable(targ[:, i]) for i in targ_channels]

        readable_obs = 'Observations:\n' + '\n'.join(obs_channel_strings)
        readable_targ = 'Targets:\n' + '\n'.join(targ_channel_strings)
        strings = [readable_obs, readable_targ]

        if model_output is not None:
            output = model_output[:, batch_index, :]
            output_strings = [_readable(output[:, i]) for i in targ_channels]
            strings.append('Model Output:\n' + '\n'.join(output_strings))

        batch_strings.append('\n\n'.join(strings))

    return '\n' + '\n\n\n\n'.join(batch_strings)

def readable_model(model_state, batch_size, whole_batch=False):
    """Produce a readable representation of the model state.

    Args:
        model_state: A tuple given by the DNC representing the internal state
            of the model.
        batch_size: size of batch.
        whole_batch: Whether to visualize the whole batch. Only the first
            sample will be visualized if False.

    Returns:
        A string used to visualize the model state.
    """
    format_string = \
        "\nRead Weights  = {}\n" \
        "Write Weights = {}\n" \
        "Free Gate     = {}\n" \
        "Alloc. Gate   = {}\n"
    if whole_batch:
        return format_string.format(
            model_state.tape_head.read_weights,
            model_state.tape_head.write_weights,
            model_state.tape_head.free_gate,
            model_state.tape_head.alloc_gate)
    return format_string.format(
        model_state.tape_head.read_weights[-1],
        model_state.tape_head.write_weights[-1],
        model_state.tape_head.free_gate[-1],
        model_state.tape_head.alloc_gate[-1])

def arr_to_str(a):
    """Create more readable array to string."""
    output = "["
    for i in range(len(a)):
        output += str(a[i])
        if i == len(a) - 1:
            break
        output += ", "
    output += "]"
    return output

def figure_data(data, output, model_state, batch_size):
    """Create string with data from the model."""
    input_str = "\n# shape = {}\nobs_str = np.array([\n".format(
        str(data.input.shape))
    for batch_index in range(batch_size):
        obs = data.input[:, batch_index, :]
        for input_row_index in range(obs.shape[0]):
            input_str += "\t" + arr_to_str(obs[input_row_index, :]) + ",\n"
    input_str += "])"

    output_str = "# shape = {}\noutput_str = np.array([\n".format(
        str(output.shape))
    for batch_index in range(batch_size):
        out = output[:, batch_index, :]
        for out_row_index in range(out.shape[0]):
            out_row = out[out_row_index, :]
            output_str += "\t" + arr_to_str(out_row) + ",\n"
    output_str += "])"

    rw_str = "# shape = {}\nrw_str = np.array([\n".format(
        str(model_state['rw'].shape))
    for rw_t in model_state['rw']:
        rw_str += "\t" + arr_to_str(rw_t[0]) + ",\n"
    rw_str += "])"

    ww_str = "# shape = {}\nww_str = np.array([\n".format(
        str(model_state['ww'].shape))
    for ww_t in model_state['ww']:
        ww_str += "\t" + arr_to_str(ww_t) + ",\n"
    ww_str += "])"

    fg_str = "# shape = {}\nfg_str = np.array([\n\t".format(
        str(model_state['fg'].shape))
    count = 0
    for fg_t in model_state['fg']:
        fg_str += str(fg_t[0]) + ", "
        count += 1
        if count == 6:
            count = 0
            fg_str += "\n\t"
    fg_str += "])"

    ag_str = "# shape = {}\nag_str = np.array([\n\t".format(
        str(model_state['ag'].shape))
    count = 0
    for ag_t in model_state['ag']:
        ag_str += str(ag_t[0]) + ", "
        count += 1
        if count == 6:
            count = 0
            ag_str += "\n\t"
    ag_str += "])"
    return input_str + "\n" + output_str + "\n" + rw_str + "\n" + ww_str \
        + "\n" + fg_str + "\n" + ag_str + "\n"

class RepeatCopy(snt.AbstractModule):
    """Sequence data generator for the task of repeating random binary pattern.

    When called, an instance of this class will return a tuple of tensorflow
    ops (obs, targ, mask), representing an input sequence, target sequence, and
    binary mask. Each of these ops produces tensors whose first two dimensions
    represent sequence position and batch index respectively. The value in
    mask[t, b] is equal to 1 iff a prediction about targ[t, b, :] should be
    penalized and 0 otherwise.
    For each realisation from this generator, the observation sequence is
    comprised of I.I.D. uniform-random binary vectors (and some flags).
    The target sequence is comprised of this binary pattern repeated
    some number of times (and some flags). Before explaining in more detail,
    let's examine the setup pictorially for a single batch element:
    ```none
    Note: blank space represents 0.
    time ------------------------------------------>
                +-------------------------------+
    mask:       |0000000001111111111111111111111|
                +-------------------------------+
                +-------------------------------+
    target:     |                              1| 'end-marker' channel.
                |         101100110110011011001 |
                |         010101001010100101010 |
                +-------------------------------+
                +-------------------------------+
    observation:| 1011001                       |
                | 0101010                       |
                |1                              | 'start-marker' channel
                |        3                      | 'num-repeats' channel.
                +-------------------------------+
    ```
    The length of the random pattern and the number of times it is repeated in
    the target are both discrete random variables distributed according to
    uniform distributions whose parameters are configured at construction time.
    The obs sequence has two extra channels (components in the trailing
    dimension) which are used for flags. One channel is marked with a 1 at the
    first time step and is otherwise equal to 0. The other extra channel is
    zero until the binary pattern to be repeated ends. At this point, it
    contains an encoding of the number of times the observation pattern should
    be repeated. Rather than simply providing this integer number directly, it
    is normalised so that a neural network may have an easier time representing
    the number of repetitions internally. To allow a network to be readily
    evaluated on instances of this task with greater numbers of repetitions,
    the range with respect to which this encoding is normalised is also
    configurable by the user. As in the diagram, the target sequence is offset
    to begin directly after the observation sequence; both sequences are padded
    with zeros to accomplish this, resulting in their lengths being equal.
    Additional padding is done at the end so that all sequences in a minibatch
    represent tensors with the same shape.
    """

    def __init__(
        self,
        num_bits=6,
        batch_size=1,
        min_length=1,
        max_length=1,
        min_repeats=1,
        max_repeats=2,
        norm_max=10,
        log_prob_in_bits=False,
        time_average_cost=False,
        name='repeat_copy',
    ):
        """Create an instance of RepeatCopy task.

        Args:
            name: A name for the generator instance (for name scope purposes).
            num_bits: The dimensionality of each random binary vector.
            batch_size: Minibatch size per realization.
            min_length: Lower limit on number of random binary vectors in the
                observation pattern.
            max_length: Upper limit on number of random binary vectors in the
                observation pattern.
            min_repeats: Lower limit on number of times the obervation pattern
                is repeated in targ.
            max_repeats: Upper limit on number of times the observation pattern
                is repeated in targ.
            norm_max: Upper limit on uniform distribution w.r.t which the
                encoding of the number of repetitions presented in the
                observation sequence is normalised.
            log_prob_in_bits: By default, log probabilities are expressed in
                units of nats. If true, express log probabilities in bits.
            time_average_cost: If true, the cost at each time step will be
                divided by the `true`, sequence length, the number of
                non-masked time steps, in each sequence before any subsequent
                reduction over the time and batch dimensions.
        """
        super(RepeatCopy, self).__init__(name=name)

        self._batch_size = batch_size
        self._num_bits = num_bits
        self._min_length = min_length
        self._max_length = max_length
        self._min_repeats = min_repeats
        self._max_repeats = max_repeats
        self._norm_max = norm_max
        self._log_prob_in_bits = log_prob_in_bits
        self._time_average_cost = time_average_cost

    def _normalise(self, val):
        return val / self._norm_max

    def _unnormalise(self, val):
        return val * self._norm_max

    @property
    def time_average_cost(self):
        """The time average cost of this task."""
        return self._time_average_cost

    @property
    def log_prob_in_bits(self):
        """The log-probability in bits."""
        return self._log_prob_in_bits

    @property
    def num_bits(self):
        """The dimensionality of each random binary vector in a pattern."""
        return self._num_bits

    @property
    def output_size(self):
        """The dimensionality of the target tensor."""
        return self._num_bits + 1

    @property
    def batch_size(self):
        """The batch size of this task."""
        return self._batch_size

    def _build(self):
        """Implement build method which adds ops to graph."""
        # short-hand for private fields.
        min_length, max_length = self._min_length, self._max_length
        min_reps, max_reps = self._min_repeats, self._max_repeats
        num_bits = self.num_bits
        batch_size = self.batch_size

        # We reserve one dimension for the num-repeats and one for the
        # start-marker.
        full_obs_size = num_bits + 2
        # We reserve one target dimension for the end-marker.
        full_targ_size = num_bits + 1
        start_end_flag_idx = full_obs_size - 2
        num_repeats_channel_idx = full_obs_size - 1

        # Samples each batch index's sequence length and the number of
        # repeats.
        sub_seq_length_batch = tf.random_uniform(
            [batch_size], minval=min_length, maxval=max_length + 1,
            dtype=tf.int32)
        num_repeats_batch = tf.random_uniform(
            [batch_size], minval=min_reps, maxval=max_reps + 1,
            dtype=tf.int32)

        # Pads all the batches to have the same total sequence length.
        total_length_batch = sub_seq_length_batch * \
            (num_repeats_batch + 1) + 3
        max_length_batch = tf.reduce_max(total_length_batch)
        residual_length_batch = max_length_batch - total_length_batch

        obs_batch_shape = [max_length_batch, batch_size, full_obs_size]
        targ_batch_shape = [max_length_batch, batch_size, full_targ_size]
        mask_batch_trans_shape = [batch_size, max_length_batch]

        obs_tensors = []
        targ_tensors = []
        mask_tensors = []

        # Generates patterns for each batch element independently.
        for batch_index in range(batch_size):
            sub_seq_len = sub_seq_length_batch[batch_index]
            num_reps = num_repeats_batch[batch_index]

            # The observation pattern is a sequence of random binary
            # vectors.
            obs_pattern_shape = [sub_seq_len, num_bits]
            obs_pattern = tf.cast(
                tf.random_uniform(
                    obs_pattern_shape, minval=0, maxval=2, dtype=tf.int32),
                tf.float32)

            # The target pattern is the observation pattern repeated n
            # times. Some reshaping is required to accomplish the tiling.
            targ_pattern_shape = [sub_seq_len * num_reps, num_bits]
            flat_obs_pattern = tf.reshape(obs_pattern, [-1])
            flat_targ_pattern = tf.tile(flat_obs_pattern,
                                        tf.stack([num_reps]))
            targ_pattern = tf.reshape(flat_targ_pattern,
                                      targ_pattern_shape)

            # Expand the obs_pattern to have two extra channels for flags.
            # Concatenate start flag and num_reps flag to the sequence.
            obs_flag_channel_pad = tf.zeros([sub_seq_len, 2])
            obs_start_flag = tf.one_hot(
                [start_end_flag_idx], full_obs_size, on_value=1.,
                off_value=0.)
            num_reps_flag = tf.one_hot(
                [num_repeats_channel_idx],
                full_obs_size,
                on_value=self._normalise(tf.cast(num_reps, tf.float32)),
                off_value=0.)

            # note the concatenation dimensions.
            obs = tf.concat([obs_pattern, obs_flag_channel_pad], 1)
            obs = tf.concat([obs_start_flag, obs], 0)
            obs = tf.concat([obs, num_reps_flag], 0)

            # Now do the same for the targ_pattern (it only has one extra
            # channel).
            targ_flag_channel_pad = tf.zeros([sub_seq_len * num_reps, 1])
            targ_end_flag = tf.one_hot(
                [start_end_flag_idx], full_targ_size, on_value=1.,
                off_value=0.)
            targ = tf.concat([targ_pattern, targ_flag_channel_pad], 1)
            targ = tf.concat([targ, targ_end_flag], 0)

            # Concatenate zeros at end of obs and begining of targ. This
            # aligns them s.t. the target begins as soon as the obs ends.
            obs_end_pad = tf.zeros([sub_seq_len * num_reps + 1,
                                   full_obs_size])
            targ_start_pad = tf.zeros([sub_seq_len + 2, full_targ_size])

            # The mask is zero during the obs and one during the targ.
            mask_off = tf.zeros([sub_seq_len + 2])
            mask_on = tf.ones([sub_seq_len * num_reps + 1])

            obs = tf.concat([obs, obs_end_pad], 0)
            targ = tf.concat([targ_start_pad, targ], 0)
            mask = tf.concat([mask_off, mask_on], 0)

            obs_tensors.append(obs)
            targ_tensors.append(targ)
            mask_tensors.append(mask)

        # End the loop over batch index.
        # Compute how much zero padding is needed to make tensors sequences
        # the same length for all batch elements.
        residual_obs_pad = [
            tf.zeros([residual_length_batch[i], full_obs_size])
            for i in range(batch_size)
        ]
        residual_targ_pad = [
            tf.zeros([residual_length_batch[i], full_targ_size])
            for i in range(batch_size)
        ]
        residual_mask_pad = [
            tf.zeros([residual_length_batch[i]]) for i in range(batch_size)
        ]

        # Concatenate the pad to each batch element.
        obs_tensors = [
            tf.concat([o, p], 0) for o, p in zip(obs_tensors,
                                                 residual_obs_pad)
        ]
        targ_tensors = [
            tf.concat([t, p], 0) for t, p in zip(targ_tensors,
                                                 residual_targ_pad)
        ]
        mask_tensors = [
            tf.concat([m, p], 0) for m, p in zip(mask_tensors,
                                                 residual_mask_pad)
        ]

        # Concatenate each batch element into a single tensor.
        obs = tf.reshape(tf.concat(obs_tensors, 1), obs_batch_shape)
        targ = tf.reshape(tf.concat(targ_tensors, 1), targ_batch_shape)
        mask = tf.transpose(
            tf.reshape(tf.concat(mask_tensors, 0), mask_batch_trans_shape))
        return DatasetTensors(obs, targ, mask)

    def cost(self, output, task_state):
        """Calculate sigmoid cross-entropy cost."""
        targ = task_state.target
        mask = task_state.mask
        return masked_sigmoid_cross_entropy(
            output,
            targ,
            mask,
            time_average=self.time_average_cost,
            log_prob_in_bits=self.log_prob_in_bits)

    def to_string(self, output, task_state, model_state, verbose=False):
        """Return human readable version of data."""
        obs = task_state.input
        unnormalised_num_reps_flag = self._unnormalise(
            obs[:, :, -1:]).round()
        obs = np.concatenate([obs[:, :, :-1], unnormalised_num_reps_flag],
                             axis=2)
        data = task_state._replace(input=obs)
        output_str = bitstring_readable(
            data, self.batch_size, output, whole_batch=verbose)
        if verbose:
            output_str += figure_data(
                data, output, model_state, self.batch_size)
        return output_str

    def process_output(self, output, task_state, model_state):
        """Change output logits to output probabilities."""
        return tf.round(
            tf.expand_dims(task_state.mask, -1) * tf.sigmoid(output))

max_stars_repo_name: derrowap/DNC-TensorFlow

max_issues_repo_name: derrowap/DNC-TensorFlow

max_forks_repo_name: derrowap/DNC-TensorFlow


text:
import os
from pathlib import Path

import imageio
import numpy as np
import SimpleITK as sitk
import h5py
from pytorch3dunet.augment import transforms
from pytorch3dunet.datasets.utils import ConfigDataset, calculate_stats
from pytorch3dunet.unet3d.utils import get_logger, expand_as_one_hot
from pytorch3dunet.datasets.hdf5 import AbstractHDF5Dataset
from pytorch3dunet.datasets.visualize import visualizer
from scipy import ndimage

import re
import yaml
import ipdb
import torch

logger = get_logger('SXTHDataset')

class SXTHDataset(AbstractHDF5Dataset):
    def __init__(self, file_path, phase, slice_builder_config, transformer_config, mirror_padding=(16, 32, 32),
                 raw_internal_path='raw', label_internal_path='label', weight_internal_path=None):
        super().__init__(file_path=file_path,
                         phase=phase,
                         slice_builder_config=slice_builder_config,
                         transformer_config=transformer_config,
                         mirror_padding=mirror_padding,
                         raw_internal_path=raw_internal_path,
                         label_internal_path=label_internal_path,
                         weight_internal_path=weight_internal_path)

    @classmethod
    def create_datasets(cls, dataset_config, phase):
        # prerocess each patient case
        sxth_case_loader = SXTHCaseLoader( dataset_config['original_data_dir']  )
        sxth_case_loader.set_phase(phase)
        cls.original_dir = Path(dataset_config['original_data_dir'])
        cls.train_dir = Path(dataset_config[phase]['file_paths'][0])
        if cls.train_dir.is_dir() == False:
            cls.train_dir.mkdir( parents = True )

        # exist already preprocessed data?
        files = os.listdir(cls.train_dir)
        need_processing = True

        if need_processing:
            # deleted existed files
            for f in files:
                    os.remove(Path(cls.train_dir)/f)

            # We use the following step to carry out the pre-process operation
            # STEP 1: loop all cases for their three axis' voxel space, determin a suitable voxel space
            # STEP 2: for each case,
            #   a. remove dark borders (air) of each slice, crop each slice to human body
            #   b. resample each slice according to the above determined suitable voxel space
            #   c. only keep the slices from the first  to the last slice with mask
            # STEP 3: loop all cases for a suitable path size

            # STEP 1:
            #case_info_dict = {}
            #rt_info_file = open("runtime_info.yml",'w',encoding='utf-8')
            seg_info = []
            foreground_vol = np.array([])
            one_case_visualizer = visualizer()

            for vol,seg,img_info in sxth_case_loader:

                spacing = img_info['Spacing']
                case_id = img_info['Case_id']
                #if phase != 'test':
                #    crop_border = cls.crop_image_only_outside(cls, vol, -1000)
                #    seg = seg[:,crop_border[0]:crop_border[1], crop_border[2]:crop_border[3]]
                #    vol = vol[:,crop_border[0]:crop_border[1], crop_border[2]:crop_border[3]]

                # resample (or re-slice) for isotropic voxel
                new_spacing = [2, 0.8, 0.8]
                resize_factor = np.array(spacing) / new_spacing
                new_real_shape = vol.shape * resize_factor
                new_shape = np.round(new_real_shape)
                real_resize_factor = new_shape / vol.shape
                new_spacing = spacing / real_resize_factor

                #case_info = { 'case_'+str(case_id): np.array( spacing ) }
                #case_info_dict.update( case_info )
                resize_factor = np.array(spacing) / new_spacing
                print(f'\t shape after crop: {vol.shape}')
                print(f'\t resampling image voxel spacing from {spacing} to {new_spacing} ...')
                vol = ndimage.zoom(vol, resize_factor, order=0 )
                print(f'\t shape after resample: {vol.shape}')

                print('\t resampling segmetation voxel ...')
                nclass = 2
                seg_resampled = np.zeros( ( nclass, vol.shape[0], vol.shape[1], vol.shape[2] ), dtype='int32' )
                if phase != 'test':
                    seg_expanded = torch.from_numpy( np.expand_dims(seg.astype('int64'),0) )
                    seg_onehot = expand_as_one_hot( seg_expanded, nclass).squeeze(0).numpy()
                    for ilabel in range( nclass):
                        seg_resampled[ilabel,] = ndimage.zoom( seg_onehot[ilabel,], resize_factor, order= 0)
                    seg = np.argmax( seg_resampled, axis = 0)
		    # collect all fore-ground voxels for further compute  max/min/mean/std value
                    foreground_vol = np.concatenate((foreground_vol, vol[seg>0]),axis = 0)
                #one_case_visualizer.save(vol,seg,'/mnt/sda2/sxth_processed/img')
                print(f'\t after removing slices without masks, image & segmentation data shape is {vol.shape} ')
                print('\t done.')

                # store as a h5d file
                f =  h5py.File(cls.train_dir/'case_{:05d}.h5'.format(case_id),'w')
                f.create_dataset('raw', data = vol)
                if phase != 'test':
                    f.create_dataset('label', data = seg)
                    f.create_dataset('weight', data = np.ones( seg.shape ) )
                f.close()
                seg_info.append(vol.shape)

        # save data to disk
        #save_img = visualizer(vol_folder = dataset_config.get(phase)['file_paths'][0])
        #save_img.save_train_cases( '/mnt/sda2/sxth_processed/img/')
        #ipdb.set_trace()
        return super().create_datasets( dataset_config, phase )

    def _load_files(dir, expand_dims):
        files_data = []
        paths = []
        for file in os.listdir(dir):
            path = os.path.join(dir, file)
            img = np.asarray(imageio.imread(path))
            if expand_dims:
                img = np.expand_dims(img, axis=0)

            files_data.append(img)
            paths.append(path)

        return files_data, paths

    @staticmethod
    def fetch_datasets(input_file_h5, internal_paths):
        return [input_file_h5[internal_path] for internal_path in internal_paths]

    @staticmethod
    def create_h5_file(file_path, internal_paths):
        return h5py.File(file_path, 'r')

    def ds_stats(self):
        # Do not calculate stats on the whole stacks when using lazy loader,
        # they min, max, mean, std should be provided in the config
        logger.info(
            'Using LazyHDF5Dataset. Make sure that the min/max/mean/std values are provided in the loaders config')
        #return 40, 250, 110, 31
        return 0, 600, None, None

    def load_case(self,cid):
        # Resolve location where data should be living
        # data_path = Path(__file__).parent.parent / "data"
        data_path = Path(self.original_dir)
        if not data_path.exists():
            raise IOError(
                    "Data path, {}, could not be resolved".format(str(data_path))
                    )

        # Get case_id from provided cid
        try:
            cid = int(cid)
            case_id = "case_{:05d}".format(cid)
        except ValueError:
            case_id = cid

        # Make sure that case_id exists under the data_path
        case_path = data_path / case_id
        if not case_path.exists():
            raise ValueError(
                    "Case could not be found \"{}\"".format(case_path.name)
                    )
        if (case_path/"imaging.nii.gz").exists() and (case_path/"segmentation.nii.gz").exists():
            vol = nib.load(str(case_path / "imaging.nii.gz"))
            seg = nib.load(str(case_path / "segmentation.nii.gz"))
            return vol, seg
        if (case_path/'imaging.nii.gz').exists() and not (case_path/'segmentation.nii.gz').exists():
            vol = nib.load(str(case_path / "imaging.nii.gz"))
            return vol

    def crop_image_only_outside(self, img, tol):
	# img is 2D or 3D image data
	# tol  is tolerance
        mask = img>tol
        if img.ndim==3:
            mask = mask.all(0)
        m,n = mask.shape
        mask0,mask1 = mask.any(0),mask.any(1)
        col_start,col_end = mask0.argmax(),n-mask0[::-1].argmax()
        row_start,row_end = mask1.argmax(),m-mask1[::-1].argmax()
        return row_start,row_end,col_start,col_end

    def plot_3d(self, image, threshold=-300):
        # Position the scan upright,
        # so the head of the patient would be at the top facing the camera
        p = image.transpose(2,1,0)  #将扫描件竖直放置
        verts, faces = measure.marching_cubes(p, threshold) #Liner推进立方体算法来查找3D体积数据中的曲面。
        fig = plt.figure(figsize=(10, 10))
        ax = fig.add_subplot(111, projection='3d')
        # Fancy indexing: `verts[faces]` to generate a collection of triangles
        mesh = Poly3DCollection(verts[faces], alpha=0.1)  #创建3Dpoly
        face_color = [0.5, 0.5, 1]
        mesh.set_facecolor(face_color)  #设置颜色
        ax.add_collection3d(mesh)
        ax.set_xlim(0, p.shape[0])
        ax.set_ylim(0, p.shape[1])
        ax.set_zlim(0, p.shape[2])
        plt.show()

class SXTHCaseLoader():
    def __init__(self, data_folder):
        self.img_folder = Path(data_folder)/'T2_case1-200'
        self.seg_folder = Path(data_folder)/'ROI_T2'
        self.reader = sitk.ImageSeriesReader()
        if not self.img_folder.exists() or not self.seg_folder.exists():
            raise IOError( "IMG or SEG path, {}, could not be resolved".format(str(img_folder)) )
        self.cases_img = os.listdir( self.img_folder )
        self.cases_seg = os.listdir( self.seg_folder )
        assert len(self.cases_img)==len(self.cases_seg)
        self.n_train = round( len(self.cases_img) * 0.8 )
        self.n_val = len(self.cases_img) - self.n_train

    def __iter__(self):
        assert self.phase == 'train' or self.phase == 'val'
        self.n_cases = self.n_train if self.phase == 'train' else self.n_val
        self.i_case = 0
        return self

    def __next__(self):
        print( f'current case: {self.i_case}/{self.n_cases}' )
        if self.i_case >= self.n_cases:
            raise StopIteration
        img_names = self.reader.GetGDCMSeriesFileNames( str(self.img_path[ self.i_case ] ) )

        self.reader.SetFileNames( img_names )
        image = self.reader.Execute( )
        image_array = sitk.GetArrayFromImage( image )
        info_spacing = image.GetSpacing()
        image_info = {'Case_id': self.i_case,'Spacing': (info_spacing[2], info_spacing[0], info_spacing[1]) }

        seg = sitk.ReadImage( str( self.seg_path[ self.i_case ] ) )
        seg_array = sitk.GetArrayFromImage( seg )

        self.i_case+=1
        return image_array, seg_array, image_info

    def set_phase( self, phase ):
        assert phase == 'train' or phase == 'val'
        self.phase = phase
        self.img_path = []
        self.seg_path = []
        cases_img_selected = self.cases_img[0:self.n_train] if self.phase == 'train' else self.cases_img[self.n_train:]
        for case_name in cases_img_selected:
            idx = re.findall('\d*.?\d', case_name)[0]
            sub_folder = os.listdir( self.img_folder/case_name )[0]
            self.img_path.append( self.img_folder/case_name/sub_folder )
            self.seg_path.append( self.seg_folder/(idx+'hrT2a_Merge.nii') )

reason: import numpy,from scipy


text:
import warnings
from typing import Callable, Optional, List, Union, Dict, Iterable

import cv2
import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn.functional as F
from catalyst.dl import Callback, CallbackOrder, IRunner, CallbackNode
from catalyst.dl.callbacks import TensorboardLogger
from catalyst.contrib.tools.tensorboard import SummaryWriter
from pytorch_toolbelt.utils import render_figure_to_tensor
from pytorch_toolbelt.utils.distributed import all_gather

from ..torch_utils import rgb_image_from_tensor, to_numpy, image_to_tensor

__all__ = [
    "get_tensorboard_logger",
    "ShowPolarBatchesCallback",
    "ShowEmbeddingsCallback",
    "UMAPCallback",
    "draw_binary_segmentation_predictions",
    "draw_semantic_segmentation_predictions",
]

def get_tensorboard_logger(runner: IRunner, tensorboard_callback_name: str = "_tensorboard") -> SummaryWriter:
    tb_callback: TensorboardLogger = runner.callbacks[tensorboard_callback_name]
    if runner.loader_name not in tb_callback.loggers:
        raise RuntimeError(f"Cannot find Tensorboard logger for loader {runner.loader_name}")
    return tb_callback.loggers[runner.loader_name]

class ShowPolarBatchesCallback(Callback):
    """
    Visualize best and worst batch based in metric in Tensorboard
    """

    def __init__(
        self,
        visualize_batch: Callable,
        metric: str = "loss",
        minimize: bool = True,
        min_delta: float = 1e-6,
        targets="tensorboard",
    ):
        """

        :param visualize_batch: Visualization function that must return list of images.
               It's takes two arguments: (batch input, predicted output).
        :param metric:
        :param minimize:
        :param min_delta:
        :param targets: Str 'tensorboard' or 'matplotlib, or ['tensorboard', 'matplotlib']
        """
        super().__init__(CallbackOrder.Logging)
        assert isinstance(targets, (list, str))

        self.best_score = None
        self.best_input = None
        self.best_output = None

        self.worst_score = None
        self.worst_input = None
        self.worst_output = None

        self.target_metric = metric
        self.num_bad_epochs = 0
        self.is_better = None
        self.visualize_batch = visualize_batch
        self.targets = [targets] if isinstance(targets, str) else targets

        if minimize:
            self.is_better = lambda score, best: score <= (best - min_delta)
            self.is_worse = lambda score, worst: score >= (worst + min_delta)
        else:
            self.is_better = lambda score, best: score >= (best + min_delta)
            self.is_worse = lambda score, worst: score <= (worst - min_delta)

    def to_cpu(self, data):
        if isinstance(data, dict):
            return dict((key, self.to_cpu(value)) for (key, value) in data.items())
        if isinstance(data, torch.Tensor):
            return data.detach().cpu()
        if isinstance(data, (list, tuple)):
            return [self.to_cpu(value) for value in data]
        return data

    def on_loader_start(self, runner):
        self.best_score = None
        self.best_input = None
        self.best_output = None

        self.worst_score = None
        self.worst_input = None
        self.worst_output = None

    def on_batch_end(self, runner: IRunner):
        value = runner.batch_metrics.get(self.target_metric, None)
        if value is None:
            warnings.warn(f"Metric value for {self.target_metric} is not available in runner.metrics.batch_values")
            return

        if self.best_score is None or self.is_better(value, self.best_score):
            self.best_score = value
            self.best_input = self.to_cpu(runner.input)
            self.best_output = self.to_cpu(runner.output)

        if self.worst_score is None or self.is_worse(value, self.worst_score):
            self.worst_score = value
            self.worst_input = self.to_cpu(runner.input)
            self.worst_output = self.to_cpu(runner.output)

    def on_loader_end(self, runner: IRunner):
        logger = get_tensorboard_logger(runner)

        if self.best_score is not None:
            best_samples = self.visualize_batch(self.best_input, self.best_output)
            self._log_samples(best_samples, "best", logger, runner.global_batch_step)

        if self.worst_score is not None:
            worst_samples = self.visualize_batch(self.worst_input, self.worst_output)
            self._log_samples(worst_samples, "worst", logger, runner.global_batch_step)

    def _log_samples(self, samples, name, logger, step):
        if "tensorboard" in self.targets:
            for i, image in enumerate(samples):
                logger.add_image(f"{self.target_metric}/{name}/{i}", image_to_tensor(image), step)

        if "matplotlib" in self.targets:
            for i, image in enumerate(samples):
                plt.figure()
                plt.imshow(image)
                plt.tight_layout()
                plt.axis("off")
                plt.show()

class ShowEmbeddingsCallback(Callback):
    def __init__(
        self,
        embedding_key,
        input_key,
        targets_key,
        prefix="embedding",
        mean=(0.485, 0.456, 0.406),
        std=(0.229, 0.224, 0.225),
    ):
        super().__init__(CallbackOrder.Logging)
        self.prefix = prefix
        self.embedding_key = embedding_key
        self.input_key = input_key
        self.targets_key = targets_key
        self.mean = torch.tensor(mean).view((1, 3, 1, 1))
        self.std = torch.tensor(std).view((1, 3, 1, 1))

        self.embeddings = []
        self.images = []
        self.targets = []

    def on_loader_start(self, runner: IRunner):
        self.embeddings = []
        self.images = []
        self.targets = []

    def on_loader_end(self, runner: IRunner):
        logger = get_tensorboard_logger(runner)
        logger.add_embedding(
            mat=torch.cat(self.embeddings, dim=0),
            metadata=self.targets,
            label_img=torch.cat(self.images, dim=0),
            global_step=runner.epoch,
            tag=self.prefix,
        )

    def on_batch_end(self, runner: IRunner):
        embedding = runner.output[self.embedding_key].detach().cpu()
        image = runner.input[self.input_key].detach().cpu()
        targets = runner.input[self.targets_key].detach().cpu().tolist()

        image = F.interpolate(image, size=(256, 256))
        image = image * self.std + self.mean

        self.images.append(image)
        self.embeddings.append(embedding)
        self.targets.extend(targets)

class UMAPCallback(Callback):
    """Visualize embeddings of the classifier using UMAP

    This callback relies on umap-learn package which must be installed beforehand:
    https://github.com/lmcinnes/umap
    """

    def __init__(
        self,
        input_key: str,
        features_key: str,
        output_key: str,
        output_activation: Callable,
        prefix: str = "umap",
        fit_params: Dict = None,
        plot_params: Dict = None,
        loaders: Iterable[str] = ("valid"),
    ):
        """

        Args:
            input_key:
            features_key:
            output_key:
            output_activation:
            prefix:
            fit_params:
            plot_params:
            loaders:
        """
        try:
            import umap
            import umap.plot
        except ImportError as e:
            print(
                "It seems `umap-learn` package is missing."
                "Please install umap (https://github.com/lmcinnes/umap) and restart the script."
            )
            raise e

        super().__init__(CallbackOrder.Metric, CallbackNode.All)
        self.features_key = features_key
        self.prefix = prefix
        self.output_key = output_key
        self.input_key = input_key
        self.output_activation = output_activation
        self.loaders = loaders

        self._plot_params = plot_params or {}
        self._fit_params = fit_params or {}

        self._reset_stats()

    def _reset_stats(self):
        self.features = []
        self.outputs = []
        self.targets = []

    def _add_to_stats(self, features, outputs, targets):
        features = to_numpy(features)
        outputs = to_numpy(outputs)
        targets = to_numpy(targets)

        self.features.extend(features)
        self.outputs.extend(outputs)
        self.targets.extend(targets)

    def _compute_embedings(self):
        features = np.array(self.features)
        features = np.concatenate(all_gather(features))

        targets = np.array(self.targets)
        targets = np.concatenate(all_gather(targets))

        outputs = np.array(self.outputs)
        outputs = np.concatenate(all_gather(outputs))

        import umap

        return umap.UMAP(**self._fit_params).fit(features), targets, outputs

    def _plot_embedings(self, logger, epoch, embeddings, targets, outputs):
        from umap import plot

        fig_gt = plot.points(embeddings, labels=targets, **self._plot_params).figure
        fig_gt = render_figure_to_tensor(fig_gt)
        logger.add_image(f"{self.prefix}/gt/epoch", fig_gt, global_step=epoch)

        fig_pred = plot.points(embeddings, labels=outputs, **self._plot_params).figure
        fig_pred = render_figure_to_tensor(fig_pred)
        logger.add_image(f"{self.prefix}/pred/epoch", fig_pred, global_step=epoch)

    def on_loader_start(self, runner: IRunner):
        """Loader start hook.
        Args:
            runner (IRunner): current runner
        """
        if runner.is_valid_loader:
            self._reset_stats()

    def on_batch_end(self, runner: IRunner):
        """Batch end hook.
        Args:
            runner (IRunner): current runner
        """
        if runner.is_valid_loader:
            self._add_to_stats(
                runner.output[self.features_key].detach(),
                self.output_activation(runner.output[self.output_key].detach()).flatten(),
                runner.input[self.input_key].detach().flatten(),
            )

    def on_loader_end(self, runner: IRunner):
        """Loader end hook.
        Args:
            runner (IRunner): current runner
        """
        if runner.is_valid_loader:
            embeddings, targets, outputs = self._compute_embedings()

            tb_logger = get_tensorboard_logger(runner)
            self._plot_embedings(
                logger=tb_logger, epoch=runner.global_epoch, embeddings=embeddings, targets=targets, outputs=outputs
            )

def draw_binary_segmentation_predictions(
    input: dict,
    output: dict,
    image_key="features",
    image_id_key: Optional[str] = "image_id",
    targets_key="targets",
    outputs_key="logits",
    mean=(0.485, 0.456, 0.406),
    std=(0.229, 0.224, 0.225),
    max_images=None,
    targets_threshold=0.5,
    logits_threshold=0,
    image_format: Union[str, Callable] = "rgb",
) -> List[np.ndarray]:
    """
    Render visualization of model's prediction for binary segmentation problem.
    This function draws a color-coded overlay on top of the image, with color codes meaning:
        - green: True positives
        - red: False-negatives
        - yellow: False-positives

    :param input: Input batch (model's input batch)
    :param output: Output batch (model predictions)
    :param image_key: Key for getting image
    :param image_id_key: Key for getting image id/fname
    :param targets_key: Key for getting ground-truth mask
    :param outputs_key: Key for getting model logits for predicted mask
    :param mean: Mean vector user during normalization
    :param std: Std vector user during normalization
    :param max_images: Maximum number of images to visualize from batch
        (If you have huge batch, saving hundreds of images may make TensorBoard slow)
    :param targets_threshold: Threshold to convert target values to binary.
        Default value 0.5 is safe for both smoothed and hard labels.
    :param logits_threshold: Threshold to convert model predictions (raw logits) values to binary.
        Default value 0.0 is equivalent to 0.5 after applying sigmoid activation
    :param image_format: Source format of the image tensor to conver to RGB representation.
        Can be string ("gray", "rgb", "brg") or function `convert(np.ndarray)->nd.ndarray`.
    :return: List of images
    """
    images = []
    num_samples = len(input[image_key])
    if max_images is not None:
        num_samples = min(num_samples, max_images)

    assert output[outputs_key].size(1) == 1, "Mask must be single-channel tensor of shape [Nx1xHxW]"

    for i in range(num_samples):
        image = rgb_image_from_tensor(input[image_key][i], mean, std)

        if image_format == "rgb":
            pass
        elif image_format == "bgr":
            image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
        elif image_format == "gray":
            image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB)
        elif isinstance(image_format, callable):
            image = image_format(image)

        target = to_numpy(input[targets_key][i]).squeeze(0)
        logits = to_numpy(output[outputs_key][i]).squeeze(0)

        overlay = image.copy()
        true_mask = target > targets_threshold
        pred_mask = logits > logits_threshold

        overlay[true_mask & pred_mask] = np.array(
            [0, 250, 0], dtype=overlay.dtype
        )  # Correct predictions (Hits) painted with green
        overlay[true_mask & ~pred_mask] = np.array([250, 0, 0], dtype=overlay.dtype)  # Misses painted with red
        overlay[~true_mask & pred_mask] = np.array(
            [250, 250, 0], dtype=overlay.dtype
        )  # False alarm painted with yellow
        overlay = cv2.addWeighted(image, 0.5, overlay, 0.5, 0, dtype=cv2.CV_8U)

        if image_id_key is not None and image_id_key in input:
            image_id = input[image_id_key][i]
            cv2.putText(overlay, str(image_id), (10, 15), cv2.FONT_HERSHEY_PLAIN, 1, (250, 250, 250))

        images.append(overlay)
    return images

def draw_semantic_segmentation_predictions(
    input: dict,
    output: dict,
    class_colors: List,
    mode="overlay",
    image_key="features",
    image_id_key="image_id",
    targets_key="targets",
    outputs_key="logits",
    mean=(0.485, 0.456, 0.406),
    std=(0.229, 0.224, 0.225),
    max_images=None,
    image_format: Union[str, Callable] = "rgb",
) -> List[np.ndarray]:
    """
    Render visualization of model's prediction for binary segmentation problem.
    This function draws a color-coded overlay on top of the image, with color codes meaning:
        - green: True positives
        - red: False-negatives
        - yellow: False-positives

    :param input: Input batch (model's input batch)
    :param output: Output batch (model predictions)
    :param class_colors:
    :param mode:
    :param image_key: Key for getting image
    :param image_id_key: Key for getting image id/fname
    :param targets_key: Key for getting ground-truth mask
    :param outputs_key: Key for getting model logits for predicted mask
    :param mean: Mean vector user during normalization
    :param std: Std vector user during normalization
    :param max_images: Maximum number of images to visualize from batch
        (If you have huge batch, saving hundreds of images may make TensorBoard slow)
    :param image_format: Source format of the image tensor to conver to RGB representation.
        Can be string ("gray", "rgb", "brg") or function `convert(np.ndarray)->nd.ndarray`.
    :return: List of images
    """
    assert mode in {"overlay", "side-by-side"}

    images = []
    num_samples = len(input[image_key])
    if max_images is not None:
        num_samples = min(num_samples, max_images)

    for i in range(num_samples):
        image = rgb_image_from_tensor(input[image_key][i], mean, std)

        if image_format == "rgb":
            pass
        elif image_format == "bgr":
            image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
        elif image_format == "gray":
            image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB)
        elif isinstance(image_format, callable):
            image = image_format(image)

        target = to_numpy(input[targets_key][i])
        logits = to_numpy(output[outputs_key][i]).argmax(axis=0)

        if mode == "overlay":
            overlay = image.copy()
            for class_index, class_color in enumerate(class_colors):
                overlay[logits == class_index, :] = class_color

            overlay = cv2.addWeighted(image, 0.5, overlay, 0.5, 0, dtype=cv2.CV_8U)
        elif mode == "side-by-side":
            true_mask = np.zeros_like(image)
            for class_index, class_color in enumerate(class_colors):
                true_mask[target == class_index, :] = class_color

            pred_mask = np.zeros_like(image)
            for class_index, class_color in enumerate(class_colors):
                pred_mask[logits == class_index, :] = class_color

            overlay = np.hstack((image, true_mask, pred_mask)).copy()
        else:
            raise ValueError(mode)

        if image_id_key is not None and image_id_key in input:
            image_id = input[image_id_key][i]
            cv2.putText(overlay, str(image_id), (10, 15), cv2.FONT_HERSHEY_PLAIN, 1, (250, 250, 250))

        images.append(overlay)

    return images

def draw_multilabel_segmentation_predictions(
    input: dict,
    output: dict,
    class_colors: List,
    mode="side-by-side",
    image_key="features",
    image_id_key="image_id",
    targets_key="targets",
    outputs_key="logits",
    mean=(0.485, 0.456, 0.406),
    std=(0.229, 0.224, 0.225),
    max_images=None,
    targets_threshold=0.5,
    logits_threshold=0,
    image_format: Union[str, Callable] = "rgb",
) -> List[np.ndarray]:
    """
    Render visualization of model's prediction for binary segmentation problem.
    This function draws a color-coded overlay on top of the image, with color codes meaning:
        - green: True positives
        - red: False-negatives
        - yellow: False-positives

    :param input: Input batch (model's input batch)
    :param output: Output batch (model predictions)
    :param class_colors:
    :param mode:
    :param image_key: Key for getting image
    :param image_id_key: Key for getting image id/fname
    :param targets_key: Key for getting ground-truth mask
    :param outputs_key: Key for getting model logits for predicted mask
    :param mean: Mean vector user during normalization
    :param std: Std vector user during normalization
    :param max_images: Maximum number of images to visualize from batch
        (If you have huge batch, saving hundreds of images may make TensorBoard slow)
    :param image_format: Source format of the image tensor to conver to RGB representation.
        Can be string ("gray", "rgb", "brg") or function `convert(np.ndarray)->nd.ndarray`.
    :return: List of images
    """
    assert mode in {"overlay", "side-by-side"}

    images = []
    num_samples = len(input[image_key])
    if max_images is not None:
        num_samples = min(num_samples, max_images)

    for i in range(num_samples):
        image = rgb_image_from_tensor(input[image_key][i], mean, std)

        if image_format == "rgb":
            pass
        elif image_format == "bgr":
            image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
        elif image_format == "gray":
            image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB)
        elif isinstance(image_format, callable):
            image = image_format(image)

        target = to_numpy(input[targets_key][i]) > targets_threshold
        logits = to_numpy(output[outputs_key][i]) > logits_threshold

        if mode == "overlay":
            overlay = image.copy()
            for class_index, class_color in enumerate(class_colors):
                overlay[logits[class_index], :] = class_color

            overlay = cv2.addWeighted(image, 0.5, overlay, 0.5, 0, dtype=cv2.CV_8U)
        elif mode == "side-by-side":
            true_mask = image.copy()
            for class_index, class_color in enumerate(class_colors):
                true_mask[target[class_index], :] = class_color

            pred_mask = image.copy()
            for class_index, class_color in enumerate(class_colors):
                pred_mask[logits[class_index], :] = class_color

            true_mask = cv2.addWeighted(image, 0.5, true_mask, 0.5, 0, dtype=cv2.CV_8U)
            pred_mask = cv2.addWeighted(image, 0.5, pred_mask, 0.5, 0, dtype=cv2.CV_8U)
            overlay = np.hstack((true_mask, pred_mask))
        else:
            raise ValueError(mode)

        if image_id_key is not None and image_id_key in input:
            image_id = input[image_id_key][i]
            cv2.putText(overlay, str(image_id), (10, 15), cv2.FONT_HERSHEY_PLAIN, 1, (250, 250, 250))

        images.append(overlay)

    return images

max_stars_repo_name: mohitktanwr/toolkits

max_issues_repo_name: mohitktanwr/toolkits

max_forks_repo_name: mohitktanwr/toolkits


text:
from astropy.coordinates import SkyCoord, AltAz, ICRS, EarthLocation
from astropy.time import Time
import astropy.units as u
import numpy as np

from approximate_coords.almanac_chp_6 import approximate_precesssion
from approximate_coords.altaz_cirs import no_polar_motion
import warnings
from astropy.utils.exceptions import ErfaWarning, AstropyWarning
import matplotlib.pyplot as plt
from time import perf_counter

warnings.filterwarnings('ignore', category=ErfaWarning)
warnings.filterwarnings('ignore', category=AstropyWarning)

coord = SkyCoord(ra='90d', dec='45d')
time = Time('2000-01-01T00:01') + u.Quantity(np.arange(0, 50 * 365, 0.1), u.day)
location = EarthLocation.of_site('Roque de los Muchachos')
print(f'Transforming coordinates for {len(time)} timestamps')

coord_altaz = coord.transform_to(AltAz(obstime=time, location=location))

t0 = perf_counter()
coord_icrs = coord_altaz.transform_to(ICRS)
astropy_time = perf_counter() - t0
print(f'Normal astropy transform took: {astropy_time:.4} s')

with approximate_precesssion():
    with no_polar_motion():
        t0 = perf_counter()
        coord_approx = coord_altaz.transform_to(ICRS)
        approx_time = perf_counter() - t0
        print(f'Approximate transforms took:   {approx_time:.4} s')
        print(f'which is {astropy_time / approx_time:.2f} times faster')

separation = coord_approx.separation(coord).arcsecond
dra = (coord_approx.ra - coord.ra).wrap_at('180d').arcsecond
ddec = (coord_approx.dec - coord.dec).wrap_at('180d').arcsecond

axs = [
    ['time', 'time', 'time'],
    ['sep', 'dra', 'ddec'],
]

fig = plt.figure(constrained_layout=True)
axd = fig.subplot_mosaic(axs)

axd['time'].plot(time.to_datetime(), dra, label=r'$\Delta$ra')
axd['time'].plot(time.to_datetime(), ddec, label=r'$\Delta$dec')
axd['time'].plot(time.to_datetime(), separation, '-', label='sep', color='k')
axd['time'].legend()
axd['time'].set_ylabel('difference / as')

axd['sep'].hist(separation, bins=100, color='k')
axd['sep'].set_xlabel('separation / as')
axd['dra'].hist(dra, bins=100, color='C0')
axd['dra'].set_xlabel(r'$\Delta$ra / as')
axd['ddec'].hist(ddec, bins=100, color='C1')
axd['ddec'].set_xlabel(r'$\Delta$dec / as')

fig.savefig('build/only_precession.pdf')

reason: import numpy,import astropy,from astropy


text:
# Copyright 2020 The Chromium Authors. All rights reserved.
# Use of this source code is governed by a BSD-style license that can be
# found in the LICENSE file.
"""Library for analyzing benchmark results for Instant start."""

import pandas
from scipy import stats

def print_report(runs, model, control='control', experiment='experiment'):
    """Print stats of A/B testing"""
    all_df = pandas.DataFrame(runs, dtype=float)
    report = pandas.DataFrame(
        columns=['Median', 'Diff with control', 'p-value'])
    for metric in sorted(set(all_df['metric_name'])):
        mdf = all_df[all_df['metric_name'] == metric]
        df = pandas.DataFrame()
        for variant in sorted(set(all_df['variant_name'])):
            df[variant] = mdf[mdf['variant_name'] == variant]\
                .value.reset_index(drop=True)

        diff_df = pandas.DataFrame()
        diff_df = df[experiment] - df[control]
        n = len(diff_df)

        row = {}
        row['Median'] = '%.1fms' % df[experiment].median()
        row['Diff with control'] = '%.1fms (%.2f%%)' % (
            diff_df.median(), diff_df.median() / df[experiment].median() * 100)
        row['p-value'] = '%f' % (stats.ttest_rel(df[experiment],
                                                 df[control])[1])
        report = report.append(pandas.Series(data=row, name=metric))
    print('%d samples on %s' % (n, model))
    print(report.sort_values(by='p-value'))

max_stars_repo_name: zealoussnow/chromium

max_issues_repo_name: zealoussnow/chromium

max_forks_repo_name: zealoussnow/chromium


text:
"""Module for testing helper functions."""
from pytest import approx
import numpy as np
import tensorflow as tf
from ..functions import log_loss

def test_log_loss(sample_predictions):
    tf.reset_default_graph()
    session = tf.Session()
    y_true, y_pred = sample_predictions

    def calc_log_loss(true_pred, pred):
        log_loss = (-true_pred * np.log(pred)) - (
            (1 - true_pred) * np.log(1 - pred))
        return log_loss

    calc_log_losses = np.vectorize(calc_log_loss)
    average_loss = np.mean(calc_log_losses(y_true, y_pred))
    test_average_loss = session.run(log_loss(y_true, y_pred, 1))
    session.close()

    assert test_average_loss == approx(average_loss)

max_stars_repo_name: maxbarkhausen/sparseflow

max_issues_repo_name: maxbarkhausen/sparseflow

max_forks_repo_name: maxbarkhausen/sparseflow


text:
from nose.tools import assert_true, assert_false
from operator import eq
import networkx as nx
from networkx.algorithms import isomorphism as iso

def test_categorical_node_match():
    nm = iso.categorical_node_match(['x', 'y', 'z'], [None]*3)
    assert_true(nm(dict(x=1, y=2, z=3), dict(x=1, y=2, z=3)))
    assert_true(not nm(dict(x=1, y=2, z=2), dict(x=1, y=2, z=1)))

class TestGenericMultiEdgeMatch:

    def setup(self):
        self.G1 = nx.MultiDiGraph()
        self.G2 = nx.MultiDiGraph()
        self.G3 = nx.MultiDiGraph()
        self.G4 = nx.MultiDiGraph()
        attr_dict1 = {'id': 'edge1', 'minFlow': 0, 'maxFlow': 10}
        attr_dict2 = {'id': 'edge2', 'minFlow': -3, 'maxFlow': 7}
        attr_dict3 = {'id': 'edge3', 'minFlow': 13, 'maxFlow': 117}
        attr_dict4 = {'id': 'edge4', 'minFlow': 13, 'maxFlow': 117}
        attr_dict5 = {'id': 'edge5', 'minFlow': 8, 'maxFlow': 12}
        attr_dict6 = {'id': 'edge6', 'minFlow': 8, 'maxFlow': 12}
        for attr_dict in [attr_dict1, attr_dict2, attr_dict3, attr_dict4, attr_dict5, attr_dict6]:
            self.G1.add_edge(1 ,2, **attr_dict)
        for attr_dict in [attr_dict5, attr_dict3, attr_dict6, attr_dict1, attr_dict4, attr_dict2]:
            self.G2.add_edge(2, 3, **attr_dict)
        for attr_dict in [attr_dict3, attr_dict5]:
            self.G3.add_edge(3, 4, **attr_dict)
        for attr_dict in [attr_dict6, attr_dict4]:
            self.G4.add_edge(4, 5, **attr_dict)

    def test_generic_multiedge_match(self):
        full_match = iso.generic_multiedge_match(['id', 'flowMin', 'flowMax'], [None]*3, [eq]*3)
        flow_match = iso.generic_multiedge_match(['flowMin', 'flowMax'], [None]*2, [eq]*2)
        min_flow_match = iso.generic_multiedge_match('flowMin', None, eq)
        id_match = iso.generic_multiedge_match('id', None, eq)
        assert_true(flow_match(self.G1[1][2], self.G2[2][3]))
        assert_true(min_flow_match(self.G1[1][2], self.G2[2][3]))
        assert_true(id_match(self.G1[1][2], self.G2[2][3]))
        assert_true(full_match(self.G1[1][2], self.G2[2][3]))
        assert_true(flow_match(self.G3[3][4], self.G4[4][5]))
        assert_true(min_flow_match(self.G3[3][4], self.G4[4][5]))
        assert_false(id_match(self.G3[3][4], self.G4[4][5]))
        assert_false(full_match(self.G3[3][4], self.G4[4][5]))

max_stars_repo_name: MarletteFunding/aws-kube-codesuite

max_issues_repo_name: MarletteFunding/aws-kube-codesuite

max_forks_repo_name: MarletteFunding/aws-kube-codesuite

reason: import networkx,from networkx