text:
module Matcha

import Base: tail

abstract type MatchSteering end
struct Greed{F, T} <: MatchSteering
    x::F
    range::T
end

Greed(x::F, limit::T) where {F, T<:Integer} = Greed{F, Range{T}}(x, 1:limit)
Greed(x::F) where {F} = Greed(x, 1:typemax(Int))

greediness(x) = 1:1
greediness(x::Greed) = x.range

abstract type MatchFunc end

trymatch(ms::MatchSteering, val, history) = trymatch(ms.x, val, history)
function trymatch(f::Union{Function, MatchFunc}, val, history)
    # an atomic pattern can use the match history by having two arguments
    if applicable(f, val, history)
        f(val, history)
    else
        f(val)
    end
end
# a pattern can also be a value
trymatch(val1, val2, history) = val1 == val2

struct History{T, VT, ST}
    buffer::T # optional record of elements for iterators that are volatile. If not volatile, this will be the actual iterator
    matches::Vector{VT} # flattened list of views for each sub pattern match
    last_begin::Ref{ST} # state of last pattern match begin
    env::Dict{Symbol, Any} # environment for more complex saving of matched variables
end
@inline Base.getindex(h::History, i::Integer) = h.matches[i]
@inline Base.getindex(h::History, key::Symbol) = h.env[key]
@inline Base.haskey(h::History, key::Symbol) = haskey(h.env, key)

# trait system
needs_recording(x) = false
view_type(x) = SubArray[]
view_type(x::T) where {T <: AbstractString} = SubString{T}[]
function buffer_type(x)
    needs_recording(x) ? eltype(x)[] : x
end

function History(list, state)
    History(
        buffer_type(list),
        view_type(list),
        Ref(state),
        Dict{Symbol, Any}()
    )
end

"""
Function that walks through `list` and saves `elem` in some way
"""
function next(history::History, list, state)
    elem_state = iterate(list, state)
    elem_state === nothing && return nothing
    elem, state = elem_state
    if needs_recording(history)
        push!(history.buffer, elem)
    end
    elem, state
end

function safe_substring(s::AbstractString, a, b)
    while !isvalid(s, a) a += 1 end
    while !isvalid(s, b) b -= 1 end
    SubString(s, a, b)
end

# seems like copy itself is not generic enough to just use it on any type
_copy(x) = copy(x)
_copy(x::String) = x
_copy(x::Ref) = Ref(x[])

function Base.copy(h::History)
    History(_copy(h.buffer), _copy(h.matches), _copy(h.last_begin), deepcopy(h.env))
end
function view_constructor(pattern::Greed, h::History{X, T, Y}, a, b) where {X, Y, T <: SubArray}
    view_constructor(pattern.x, h, a, b)
end
function view_constructor(pattern, h::History{X, T, Y}, a, b) where {X, Y, T <: SubArray}
    view(h.buffer, a:b)
end
function view_constructor(pattern, h::History{X, T, Y}, a, b) where {X, Y, T <: SubString}
    safe_substring(h.buffer, a, b)
end
function finish_match(pattern, matched, h::History{T, VT, ST}, state) where {T, VT, ST}
    if matched
        push!(h.matches, view_constructor(pattern, h, h.last_begin[], state))
    else
        empty!(h.env)
    end
    matched, h, state
end
function start_match(history, state)
    history.last_begin[] = state
end

function inner_matchat(
        list, last_state,
        patterns::NTuple{N, Any},
        history = History(list, last_state)
    ) where N
    elem_state = next(history, list, last_state)
    elem_state === nothing && return (false, history, last_state, 0)
    elem, state = elem_state
    matches = 0; lastmatchstate = last_state
    start_match(history, last_state)
    pattern = patterns[1]
    while true
        # greed can make one fail, but it depends on the circumstances
        greedrange = greediness(pattern)
        enough = matches in greedrange # we have enough matches when in the range of greed

        # okay lets get matchin'
        matched = trymatch(pattern, elem, history)
        if matched
            matches += 1
            lastmatchstate = last_state
        else
            # we don't have enough matches yet to fail matching, or we don't have any more patterns to match
            if !enough || N == 1
                return finish_match(pattern, enough, history, enough ? lastmatchstate : state)..., matches # we fail or not, depending whether we have enough
            end
            if N > 1
                # okay, we failed but already have enough.
                # The only chance to continue is that next pattern matches
                # this is final, so no copy of history for backtracking needed
                if matches > 0
                    finish_match(pattern, true, history, lastmatchstate)
                end
                return inner_matchat(list, last_state, tail(patterns), history)
            end
        end
        # after match, needs to update enough
        enough = matches in greedrange

        if N > 1 && enough
            # we're in a state were the current pattern can/should stop matching
            # this is where a match of the next pattern could end things!
            if matches == last(greedrange) # we actually are at the last allowed
                finish_match(pattern, true, history, lastmatchstate)
                return inner_matchat(list, state, tail(patterns), history)
            else
                # a copy of history is needed, since we can backtrack
                newbranch = copy(history)
                finish_match(pattern, true, newbranch, lastmatchstate - 1)
                ismatch, history2, state2, n = inner_matchat(list, last_state, tail(patterns), newbranch)
                (ismatch && n > 0) && return true, history2, state2, matches
            end
        elseif N == 1 && matches == last(greedrange) # rest is empty and we have enough -> stop!
            res = finish_match(pattern, true, history, lastmatchstate)
            return res..., matches
        end

        # we matched!
        # But if we have enough and the next pattern starts matching, we must stop here
        last_state = state
        elem_state = next(history, list, state)
        if elem_state === nothing
            # this madness is over!
            # if succesfull or not depends on whether we have enough and no rest!
            success = enough && (N == 1 || all(x-> 0 in greediness(x), tail(patterns)))
            return finish_match(pattern, success, history, success ? lastmatchstate : state)..., matches
        end
        # okay, we're here, meaning we matched in some way and can continue making history!
        elem, state = elem_state
    end
    return false, history, last_state, matches # should be dead code
end

function matchat(
        list, patterns::Tuple,
    )
    matchat(list, 1, patterns)
end
function matchat(
        list, state, patterns::Tuple,
    )
    history = History(list, state)
    matched, hist, state, n = inner_matchat(
        list, state, patterns, history
    )
    matched, hist
end

function matchone(
        list, patterns::Tuple,
    )
    matchone(list, 1, patterns)
end
function matchone(
        list, state, patterns::Tuple,
    )
    history = History(list, state)
    while true
        history = History(list, state)
        match, history, _, n = inner_matchat(list, state, patterns, history)
        match && return true, history.matches
        elem_state = itereate(list, state)
        elem_state !== nothing && break
        elem, state = elem_state
    end
    return false, history.matches
end

function matchitall(
        list, patterns::Tuple,
    )
    matchitall(list, 1, patterns)
end

# TODO find better non clashing name with Base
function matchitall(
        list, state, patterns::Tuple,
    )
    history = History(list, state)
    matches = typeof(history)[]
    while true
        history = History(list, state)
        match, history, _, n = inner_matchat(list, state, patterns, history)
        if match
            push!(matches, history)
        end
        elem_state = iterate(list, state)
        elem_state === nothing && break
        elem, state = elem_state
    end
    return matches
end

function forward(x, elem, state, n)
    for i in 1:n
        elem_state = iterate(x, state)
        elem_state === nothing && break
        elem, state = elem_state
    end
    elem, state
end
slength(x::Union{Tuple, AbstractArray}) = length(x)
slength(x) = 1

@inline firstindex(v::History) = firstindex(v.matches)
@inline lastindex(v::History) = lastindex(v.matches)

@inline firstindex(v::SubArray) = parentindices(v)[1][1]
@inline firstindex(v::Union{Vector, Tuple}) = firstindex(first(v))

@inline lastindex(v::SubArray) = parentindices(v)[1][end]
@inline lastindex(v::Union{Vector, Tuple}) = lastindex(last(v))

function matchreplace(f, list, patterns)
    matches = matchitall(list, patterns)
    isempty(matches) && return copy(list)
    result = similar(list, 0)
    elem_state = iterate(list)
    i = 1
    cmatch = matches[i]
    while elem_state !== nothing
        elem, state = elem_state
        last_state = state
        isreplace = i <= length(matches) && last_state == firstindex(cmatch)
        replacements, n = if isreplace
            i += 1
            n = sum(map(slength, cmatch.matches))
            elem, state = forward(list, elem_state, n - 1)
            tmp = f(cmatch)
            r = isa(tmp, Tuple) ? tmp : (tmp,)
            if i <= length(matches)
                cmatch = matches[i]
            end
            r, n
        else
            (elem,), 1
        end
        for r in replacements
            push!(result, r)
        end
        elem_state = iterate(list, state)
    end
    result
end

function matchreplace(f, list::AbstractVector, patterns)
    matches = matchitall(list, patterns)
    isempty(matches) && return copy(list)
    result = similar(list, 0)
    lastidx = 1
    for cmatch in matches
        a, b = firstindex(cmatch), lastindex(cmatch)
        N =  sum(map(slength, cmatch.matches))
        @assert (b - a) == (N - 1) "$a $b $N"
        a < lastidx && continue # ignore matches that go back
        append!(result, view(list, lastidx:a-1))
        lastidx = b + 1
        tmp = f(cmatch)
        r = isa(tmp, Tuple) ? tmp : (tmp,)
        push!(result, r...)
    end
    append!(result, @view list[lastidx:end])
    result
end

alwaysmatch(x) = true
const anything = Greed(alwaysmatch)

export matchat, matchone, matchitall, matchreplace
export Greed
export anything

end # module


text:
## A tool for tokenizing aetherlang expressions

function aetherlang_tokenize(code::String)::Vector{String}
    """Tokenizes an expression, returns Vector{String}. Also makes everything except strings lowercase."""
    tokens::Vector{String} = []
    p_count::UInt = 0     # opened parentheses count
    quotes::Bool = false  # whether '\"' was opened
    special::Bool = false # whether there was a '\\' found inside a string literal
    token::String = ""
    for c::Char in code
        if isspace(c) && p_count == 0 && !quotes
            if token != ""; push!(tokens, token); end
            token = ""
        else
            if c == '\"' && (!quotes || (quotes && !special))
                quotes = !quotes
            elseif c == '(' && !quotes
                p_count += 1
            elseif c == ')' && !quotes
                p_count -= 1
            end
            token *= quotes ? c : lowercase(c)
        end
        special = c == '\\' && quotes
    end
    if quotes || special || p_count != 0
        throw(AetherlangError("Missing closing character. Use '\\' when splitting a line into two:\n(F X X \\\n\tX X)"))
    end
    if token != ""
        push!(tokens, token)
    end
    tokens
end

ext: jl

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
include("loadtrace.jl")

type template
  this_unit::Int
  classMeans::Dict
  tempSumDiff::Array
  POI::Array
  classMeanCovs::Dict
end

#this function constructs templates in a standard manner as presented in
#the semiinal work on template attacks user may add different functions
#to make template in different manner
#this function takes in the output of loadTA_TRAINtrace which is array of TA_TRAINdata type
function buildTemplates(data,cipher_type)
  unit_size = data.unit_size
  num_unit = data.num_unit
  plaintexts = data.plaintexts
  ciphertexts = data.ciphertexts
  keys = data.keys
  traces = data.traces
  res_dir = data.res_dir
  sbox = S_boxes[cipher_type]

  pt_unit = div(num_unit,unit_size)

  #will return this aaray of template
  templates = template[]

  #building templates for byte_n, templates will be
  #I have to build templates for all the bytes
  for byte_n in 1:16
      #based on the HW values which will vary from 0 to 8 hence 9 templates
      #storing the sbox result of the first bytes in the sample
      tempSbox = zeros(Int,length(plaintexts))
      for i in 1:length(plaintexts)
          p_byte = parse(Int,plaintexts[i][(byte_n-1)*pt_unit+1:(byte_n-1)*pt_unit+2],16)
          k_byte = parse(Int,keys[i][(byte_n-1)*pt_unit+1:(byte_n-1)*pt_unit+2],16)
          tempSbox[i] = sbox[p_byte$k_byte+1]
      end
      #storing the HW values of the first byte of the output of the sbox
      #this will decide which trace belongs to which class
      tempHW = zeros(Int,length(plaintexts))
      for i in 1:length(plaintexts)
          tempHW[i] = HW[tempSbox[i]]
      end
      #a dictionary which consistis of elements of the same class bunched together
      #as a 2D array with key equal to their representing HW values
      tempClass = Dict()
      for i in 0:unit_size
          tempClass[i] = []
      end
      #filling the dictionary with the right traces in right classes
      for i in 1:length(plaintexts)
          hw = tempHW[i]
          push!(tempClass[hw],traces[i,:])
      end
      #now I have to calculate the means of traces present in each of the classes
      tempMeans = Dict()
      for i in 0:unit_size
          tempMeans[i] = mean(tempClass[i])
      end
      POI = get_pois(tempMeans, unit_size)[1]
      tempSumDiff = get_pois(tempMeans, unit_size)[2]
      #now for each class I want to build a multivariate model for which
      #I need to define the mean vector and the covariance matrix
      tempMeanCov = Dict()
      for i in 0:unit_size
          meanvec = mean(tempClass[i])[POI]
          traceMat = zeros(length(tempClass[i]),length(traces[1,:]))
          for j in 1:length(tempClass[i])
              tr = tempClass[i][j]
              for k in 1:length(traces[1,:])
                  traceMat[j,k] = tr[k]
              end
          end

          covmat = zeros(length(POI),length(POI))
          for u  in 1:length(POI)
              for v in 1:length(POI)
                  covmat[u,v] = cov(traceMat[:,POI[u]],traceMat[:,POI[v]])
              end
          end
          tempMeanCov[i] = [meanvec, covmat]
      end
      push!(templates,template(byte_n, tempMeans, tempSumDiff, POI, tempMeanCov))
  end
  #might wanna write it down in JLD file
  return templates
end

function get_pois(tempMeans, unit_size)
    tempSumDiff = zeros(length(tempMeans[0]))
    for i in 0:unit_size
        for j in i+1:unit_size
            tempSumDiff += abs(tempMeans[i] - tempMeans[j])
        end
    end
    #now i need to find some POI's from tempSumDiff
    n_pois = 10
    poi_spacing = 5
    tempSumDiff_copy = tempSumDiff
    POI = []
    for i in 1:n_pois
        nextPoi = indmax(tempSumDiff)
        push!(POI,nextPoi)
        inv_left = max(1,nextPoi - poi_spacing)
        inv_right = min(length(tempSumDiff),nextPoi+poi_spacing)
        for j in inv_left:inv_right
            tempSumDiff[j] = 0
        end
    end

    return [POI, tempSumDiff_copy]
end

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
@testset "AlgAss" begin
  include("AlgAss/AbsAlgAss.jl")
  include("AlgAss/AlgAss.jl")
  include("AlgAss/AlgGrp.jl")
  include("AlgAss/AlgMat.jl")
  include("AlgAss/Elem.jl")
  include("AlgAss/Ideal.jl")
end


text:
#include("time_series_opf.jl")

# function to plot the given data files with given plot choice parameters
#
# function parameters
#################################################
# month_st, day_st, period_st, month_en, day_en, period_en: the start and end time for the period of data to plot
# br_arr_raw: the array containg the indexex of subnets to plot (in format of Int64). If using the single for plot_choice, the function only reads in the first index for default.
# folder: the target folder containg the data files for plot
# output_folder: the folder specified to save the output files
# plot_choice::String = "single": default single,
#       "single": optimal result & frequency of a network in a range
#       "compare": compare average value of multiple networks

# all network indexes:
# [1 2 3 4 5 6 8 9 10 11 12 13 14 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
# 43 44 45 46 47 49 50 51 52 53 54 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
# 83 84 85 87 88 89 90 91 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
# 118 119]

# the function of plotting results of multifrequency-opf on different netwoks and in different periods
function multi_plot_new(month_st, day_st, period_st, month_en, day_en, period_en, br_arr_raw, folder, output_folder,
      plot_choice::String = "single")

      # time range
      start_t = DateTime(2020,month_st,day_st,period_st)
      end_t = DateTime(2020,month_en,day_en,period_en)
      range_t = start_t:Hour(1):end_t

      all_range = DateTime(2020,1,1,1):Hour(1):DateTime(2020,12,31,24)
      part_range = DateTime(2020, month_st, day_st, period_st):Hour(1):DateTime(2020,month_en, day_en, period_en)
      # num_br: the number of branches we want to plot
      num_br = max(size(br_arr_raw, 1), size(br_arr_raw, 2))

      # br_result: the optimization result of the subnets (temporary storage)
      # res_summary_dict: the final data structure of optimization results of the subnets
      # br_arr: array of the indexes of all the subnets into string format (for further indexing format)
      br_result = Dict()
      res_summary_dict = Dict()
      br_arr = Array{String, 1}()
      for i = 1:num_br
            push!(br_arr, string(br_arr_raw[i]))
      end

      base_info = JSON.parsefile("D:/github/VariableFrequencyOPF/time_series_out_modified/res_summary_base_(1.1.1__12.31.24).json")
      # base_info_part = Dict()
      # for i = 1:size(part_range, 1)
      #       base_info_part[i] = deepcopy(base_info)
      # end
      for br_ind = 1:num_br
            # input result summary from specified input folder (output folder of the optimization function)
            #dict_info = JSON.parsefile("$(folder)/res_summary_br$(br_arr[br_ind])_(1.1.1__12.31.24).json")
            #dict_info = JSON.parsefile("$(folder)/res_summary_br$(br_arr[br_ind])_(7.15.1_7.15.24).json")
            dict_info = JSON.parsefile("$(folder)/res_summary_br$(br_arr[br_ind])_(1.1.1__12.31.24).json")
            br_result[br_ind] = dict_info
            # the result dictionary can be of different data structure
            res_summary_dict[br_arr[br_ind]] = deepcopy(br_result[br_ind])

            # if (plot_choice == "comp")
            #       # modify all the types of the fields to be Float64
            #       for i in 1:size(all_range, 1)
            #             # convert the optimal value to Float64
            #             if typeof(res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][20][1]) == String
            #                   res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][20][1] =
            #                        parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][20][1])
            #             end
            #
            #             # # convert the optimal frequency to Float64
            #             # if typeof(res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][2][1]) == String
            #             #       res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][2][1] =
            #             #            parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][2][1])
            #             # end
            #       end

            # normalize the data formats
            # loop through the locations of results from the output file of optimization function
            # ( the data format need to be changed for different )
            for i in 1:size(part_range, 1)
            #elseif (plot_choice == "single")

                  # convert the optimal value to Float64
                  # if typeof(res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1]) == String
                  #       res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1] =
                  #            parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1])
                  # end

                  # convert the optimal frequency to Float64
                  if typeof(res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1]) == String
                        res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1] =
                             parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1])
                            #FIXME parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["f"][1])
                  end

                  # if typeof(res_summary_dict[br_arr[br_ind]]["$(i)"]["f"][1]) == String
                  #       res_summary_dict[br_arr[br_ind]]["$(i)"]["f"][1] =
                  #            parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["f"][1])
                  #           #FIXME parse(Float64, res_summary_dict[br_arr[br_ind]]["$(i)"]["f"][1])
                  # end
            end

      end

      # change the type of the object for base case

      for i in 1:size(all_range, 1)
            if typeof(base_info["$(i)"]["objective"][1]) == String
                  base_info["$(i)"]["objective"][1] =
                       parse(Float64, base_info["$(i)"]["objective"][1])
            end
      end

      # the dictionaries to store output values (objective value/frequency/size), "_n" represents the array of points
      # where the solver did not find objective values
      objval = Dict()
      objval_n = Dict()
      obj_freq = Dict()
      obj_freq_n = Dict()
      size_obj = Array{DateTime, 1}()
      size_obj_n = Array{DateTime, 1}()
      objval_base = Array{Float64, 1}()

      # this is the single array that store the data for single branch for comparison to base case
      objval_all = Dict()
      size_obj_all = Array{DateTime, 1}()
      #objfreq_base = Array{Float64, 1}()

      #if (plot_choice == "comp")
            for br_ind = 1:num_br
                  objval[br_ind] = Array{Float64, 1}()
                  objval_n[br_ind] = Array{Float64, 1}()
                  obj_freq[br_ind] = Array{Float64, 1}()
                  obj_freq_n[br_ind] = Array{Float64, 1}()
                  objval_all[br_ind] = Array{Float64, 1}()
            end
      #end
      #
      # for i in 1:size(range_t, 1)
      #       if ((res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][22][1] in ["LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"])
      #          && (res_summary_dict[br_arr[br_ind]]["$(i)"]["columns"][22][2] in ["LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"]))
      #          push!(size_obj[br_ind], range_t[i])
      #       else
      #          push!(size_obj_n[br_ind], range_t[i])
      #       end
      # end

      # find the starting time for the indexing below based on input time range and the whole year time range
      start_t_ind = findfirst(isequal(start_t), part_range)
      start_all_ind = findfirst(isequal(start_t), all_range)
      println("$(start_t_ind)")
      println("$(start_all_ind)")

      # for plot average value comparison:

      # here divide each time index into feasible array and infeasible array
      # index start from the start time above and loop for the size of the range we want to plot for
      for i in start_t_ind:(start_t_ind + size(range_t, 1) - 1)
            push!(size_obj_all, range_t[i-start_t_ind + 1])
            if ((res_summary_dict[br_arr[1]]["$(i)"]["status"] in ["LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"]))
               push!(size_obj, range_t[i-start_t_ind + 1])
            else
               push!(size_obj_n, range_t[i-start_t_ind + 1])
               println("infeasible$(i), $(res_summary_dict[br_arr[1]]["$(i)"]["status"])")
            end
      end

      # store data for base case

      for i in start_t_ind:(start_t_ind + size(range_t, 1) - 1)
            all_ind = start_all_ind + i - 1
            println(all_ind)
            push!(objval_base, base_info["$(all_ind)"]["objective"][1])
      end

      # store the values of objective values
      for br_ind = 1:num_br
            for i in start_t_ind:(start_t_ind + size(range_t, 1) - 1)
                  #println("absd$(i)")
                  push!(objval_all[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1])

                  # if ((res_summary_dict[br_arr[br_ind]]["$(i)"]["status"] in ["LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"]))
                  #    push!(objval[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1])
                  #    push!(obj_freq[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1])
                  #    #push!(size_obj[br_ind], range_t[i])
                  # else
                  #    push!(objval_n[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1])
                  #    push!(obj_freq_n[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1])
                  #    #push!(size_obj_n[br_ind], range_t[i])
                  # end

                  if ((res_summary_dict[br_arr[br_ind]]["$(i)"]["status"] in ["LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"]))
                     push!(objval[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1])
                     push!(obj_freq[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1])
                     #push!(size_obj[br_ind], range_t[i])
                  else
                     push!(objval_n[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["objective"][1])
                     push!(obj_freq_n[br_ind], res_summary_dict[br_arr[br_ind]]["$(i)"]["frequency (Hz)"][1])
                     #push!(size_obj_n[br_ind], range_t[i])
                  end
            end
      end

      # string = JSON.json(res_summary_dict)
      # open("test.json", "w") do f
      #     write(f, string)
      # end

      #size_obj_n = size(range_t, 1) - size_obj
      # below are the parts of plotting
      #println(range_t)
      #println(objval)

      # plot if the plot choice is "single", which means we plot the average optimal values and frequency
      # versus the branch numbers.
      if (plot_choice == "single")
            gr()

            # plot for the optimal cost

            # p = bar(size_obj, objval[1],  color = :blue,  linecolor = :transparent, title = "optimal result: br$(br_arr[1])",
            #           xlabel = "time/h", ylabel = "optimal cost", titlefontsize=20,
            #           right_margin = 10px)
            #          #xtickfontsize=18,ytickfontsize=18,
            # plot!(size_obj_n, objval_n[1], color = :red, seriestype = :scatter)
            # if ((num_br - 1) > 0)
            #       for br_ind = 2:num_br
            #             plot!(size_obj, objval[br_ind], color = :green, seriestype = :line)
            #             plot!(size_obj_n, objval_n[br_ind], color = :yellow, seriestype = :line)
            #       end
            # end

            # plot for the optimal frequency
            p2 = plot(size_obj, obj_freq[1], color = :blue, seriestype = :line, title = "optimal frequency: br$(br_arr[1])",
                     xlabel = "time/h", ylabel = "optimal frequency", titlefontsize=20,
                     right_margin = 10px, labels = ["feasible" "infeasible"])
            plot!(size_obj_n, obj_freq_n[1], color = :red, seriestype = :scatter,
                  markersize = 2.5)
            if ((num_br - 1) > 0)
                  for br_ind = 2:num_br
                        plot!(size_obj, obj_freq[br_ind], color = :green, seriestype = :line)
                        plot!(size_obj_n, obj_freq_n[br_ind], color = :blue, seriestype = :line)
                  end
            end

            # get the data difference compared to base case and store to array
            base_diff = Array{Float64, 1}()
            base_diff_percent = Array{Float64, 1}()
            objval_all_mod = Array{Float64, 1}()

            number_bins = 30
            total_cost_improve = 0
            total_cost_base = 0

            # FIXME get the max improve rate
            max_improve = 0;

            for i = 1:size(range_t, 1)
                  push!(base_diff, objval_base[i] - objval_all[1][i])
                  push!(objval_all_mod, objval_all[1][i]-25000)
                  push!(base_diff_percent, (objval_base[i] - objval_all[1][i])/objval_base[i])
                  total_cost_improve += (objval_base[i] - objval_all[1][i]);
                  if ((objval_base[i] - objval_all[1][i]) < 0)
                        println("time point $(i) is not right.")
                  end

                  # calculate the total cost of the base case
                  total_cost_base += objval_base[i];
            end

            max_base_diff_percent = maximum(base_diff_percent)
            min_base_diff_percent = minimum(base_diff_percent)
            println("The maximum percent of difference is $(max_base_diff_percent)")

            total_improve_rate = total_cost_improve/total_cost_base
            #bins_use=range(min_base_diff_percent,max_base_diff_percent,length=number_bins)
            bins=range(min_base_diff_percent,max_base_diff_percent,length=number_bins)

            base_comp_data = [base_diff objval_all[1]]
            base1 = JSON.json(base_comp_data)
            open("time_series_out/opt_res_plot_single_10_(1.1.1__1.1.3).json", "w") do f
              write(f, base1)
            end
            # zeros(Float64, size(range_t, 1), 2)
            # for i = 1:size(range_t, 1)
            #       push!(base_diff, objval_base[i] - objval_all[1][i])
            #       base_comp_data[i][1] = objval_all[1][i]
            #       base_comp_data[i][2] = base_diff[i]
            # end

            p1 = groupedbar(base_comp_data, bar_position = :stack, bar_width=0.7, linecolor = :transparent,
                        title = "optimal value VS base: br$(br_arr[1])",
                      xlabel = "time/h", ylabel = "opt val diff", titlefontsize=20,
                      right_margin = 10px)
            # p3 = plot(size_obj_all, base_diff, color = :yellow, seriestype = :bar, title = "optimal value VS base",
            #          xlabel = "time/h", ylabel = "optimal value difference", right_margin = 10px, labels = ["feasible" "infeasible"])

            # combine two plots together on a single figure

            p4 = bar(size_obj_all, base_diff_percent, line_color = :transparent, title = "Rate of improvement: br$(br_arr[1])",
                        linecolor = :transparent, xlabel = "time/h", ylabel = "rate improved", titlefontsize=20,
                        right_margin = 10px)

            p5 = histogram(base_diff_percent, bins = bins, linecolor = :transparent, yaxis=(:log10), title = "Num of hours: br$(br_arr[1])",
                        xlabel = "rate of improvements (cost reduceed/base cost)", ylabel = "amount of hours", titlefontsize=20,
                        right_margin = 10px)

            p_result = plot(p1, p2, layout = (2,1), legend = false)

            savefig(p_result, "$(output_folder)/opt_res_plot_single_$(br_arr[1])_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).pdf")

            p_improve = plot(p4, p5, layout = (2,1), legend = false)

            savefig(p_improve, "$(output_folder)/opt_improve_single_$(br_arr[1])_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).pdf")

            println("the current total improvement rate is $(total_improve_rate)")
      end

      # plot if the plot choice is "compare", which means we plot the average optimal values and frequency
      # versus the time steps (the comparison between different ).
      if (plot_choice == "compare")
            #objval_ave_dic = Dict()
            objval_ave = Array{Float64, 1}()
            opt_f_ave = Array{Float64, 1}()
            # objval_ave_dic[1] = Array{String, 1}()
            # objval_ave_dic[2] = Array{Float64, 1}()

            # the base value to help make the plot easy to analyze ( show the differences between the different cases)
            base_val_plot = 25000
            for br_ind = 1:num_br
                  # for i in 1:size()
                  # push!(objval_ave_dic[1], br_arr[br_ind])
                  # push!(objval_ave_dic[2], sum(objval[br_ind])/size(objval[br_ind], 1))
                  push!(objval_ave, (sum(objval[br_ind])/size(objval[br_ind], 1)) - base_val_plot)
                  push!(opt_f_ave, (sum(obj_freq[br_ind])/size(obj_freq[br_ind], 1)))
                  #end
            end
            #sort(collect(objval_ave_dic), by=x->x[2])
            ind = sortperm(objval_ave)
            println(ind)

            objval_order = Array{Float64, 1}()
            subnet_order = Array{Int64, 1}()
            for i = 1:num_br
                  push!(subnet_order, br_arr_raw[ind[i]])
                  push!(objval_order, objval_ave[ind[i]])
            end
            println(subnet_order)
            println(objval_order)
            println("result is above")

            net_order = JSON.json(subnet_order)
            val_order = JSON.json(objval_order)
            open("$(output_folder)/opt_data_comp_$(num_br)_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).json", "w") do f
                  write(f, net_order)
                  write(f, "\n")
                  write(f, val_order)
            end
            # string = JSON.json(res_summary_dict)
            # open("test.json", "w") do f
            #     write(f, string)
            # end

            gr()
            # p1 = plot(br_arr[ind], objval_ave[ind],  color = :blue, seriestype = :bar, title = "partial optimal result plot ($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en))",
            #          xlabel = "subnet number", ylabel = "opt_val - $(base_val_plot)", xtickfontsize=18,ytickfontsize=18,titlefontsize=20,right_margin = 10px, legend = false)
            # p2 = plot(br_arr[ind], opt_f_ave[ind], color = :blue, seriestype = :bar, title = "partial optimal frequency plot ($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en))",
            #          xlabel = "subnet number", ylabel = "optimal f", xtickfontsize=18,ytickfontsize=18,titlefontsize=20, right_margin = 10px, legend = false)
            #
            # p_all = plot(p1, p2, layout = (2,1), legend = false)
            # savefig(p_all, "$(output_folder)/opt_res_plot_multi_comp_$(num_br)_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).pdf")

            # get the data difference compared to base case and store to array
            base_diff = Dict()
            base_diff_percent = Dict()
            base_diff_percent_posi = Dict()
            base_diff_percent_all = Dict()
            objval_all_mod = Dict()

            for br_ind = 1:num_br
                  base_diff[br_ind] = Array{Float64, 1}()
                  base_diff_percent[br_ind] = Array{Float64, 1}()
                  base_diff_percent_posi[br_ind] = Array{Float64, 1}()
                  base_diff_percent_all[br_ind] = Array{Float64, 1}()
                  objval_all_mod[br_ind] = Array{Float64, 1}()
            end

            number_bins = 30
            total_cost_improve = 0
            total_cost_base = 0

            # FIXME get the max improve rate
            max_improve = 0;
            for br_ind = 1:num_br
                  for i = 1:size(range_t, 1)
                        push!(base_diff[br_ind], objval_base[i] - objval_all[br_ind ][i])
                        push!(objval_all_mod[br_ind], objval_all[br_ind ][i]-25000)
                        push!(base_diff_percent[br_ind], (objval_base[i] - objval_all[br_ind ][i])/objval_base[i])
                        if ((objval_base[i] - objval_all[br_ind ][i])/objval_base[i])>0
                              push!(base_diff_percent_posi[br_ind], (objval_base[i] - objval_all[br_ind ][i])/objval_base[i])
                              push!(base_diff_percent_all[br_ind], (objval_base[i] - objval_all[br_ind ][i])/objval_base[i])
                        else
                              push!(base_diff_percent_all[br_ind], 0.0)
                        end
                        #total_cost_improve += (objval_base[i] - objval_all[br_ind ][i]);
                        if ((objval_base[i] - objval_all[br_ind][i]) < 0)
                              println("br$(br_ind), time point $(i) is not right.")
                        end

                        # calculate the total cost of the base case
                        #total_cost_base += objval_base[i];
                  end
            end

            max_base_diff_percent = 0;
            for br_ind = 1:num_br
                  for i = 1:size(base_diff_percent[br_ind], 1)
                        if max_base_diff_percent < base_diff_percent[br_ind][i]
                              max_base_diff_percent = base_diff_percent[br_ind][i]
                        else
                              #
                        end
                  end
            end
            # min_base_diff_percent = minimum(base_diff_percent)
            println("The maximum percent of difference is $(max_base_diff_percent)")

            #total_improve_rate = total_cost_improve/total_cost_base
            #bins_use=range(min_base_diff_percent,max_base_diff_percent,length=number_bins)
            bins=range(0,max_base_diff_percent,length=number_bins)

            # base_comp_data = Dict()
            # for br_ind = 1:num_br
            #       push!(base_comp_data, [base_diff[br_ind] objval_all[1][br_ind]])
            # end
            # base1 = JSON.json(base_diff_percent)
            # open("opt_res_plot_single_10_(1.1.1__12.31.24).json", "w") do f
            #   write(f, base1)
            # end
            # zeros(Float64, size(range_t, 1), 2)
            # for i = 1:size(range_t, 1)
            #       push!(base_diff, objval_base[i] - objval_all[1][i])
            #       base_comp_data[i][1] = objval_all[1][i]
            #       base_comp_data[i][2] = base_diff[i]
            # end

            p3 = bar(size_obj_all, base_diff_percent_all[1], line_color = :transparent, title = "rate of improvement: br$(br_arr[1])",
                       linecolor = :transparent, xlabel = "time/h", ylabel = "rate improved", titlefontsize=20,
                       right_margin = 10px, )
            p4 = bar(size_obj_all, base_diff_percent_all[2], line_color = :transparent, title = "rate of improvement: br$(br_arr[2])",
                       linecolor = :transparent, xlabel = "time/h", ylabel = "rate improved", titlefontsize=20,
                       right_margin = 10px)
            p5 = bar(size_obj_all, base_diff_percent_all[3], line_color = :transparent, title = "rate of improvement: br$(br_arr[3])",
                       linecolor = :transparent, xlabel = "time/h", ylabel = "rate improved", titlefontsize=20,
                       right_margin = 10px)
            p6 = bar(size_obj_all, base_diff_percent_all[4], line_color = :transparent, title = "rate of improvement: br$(br_arr[4])",
                        linecolor = :transparent, xlabel = "time/h", ylabel = "rate improved", titlefontsize=20,
                        right_margin = 10px)

            p7 = histogram(base_diff_percent_posi[1], bins = bins, linecolor = :transparent, yaxis=(:log10), title = "br$(br_arr[1])",
                       xlabel = "rate of improvements (cost reduceed/base cost)", ylabel = "amount of hours", titlefontsize=20,
                       right_margin = 10px, xlim = max_base_diff_percent)
            p8 = histogram(base_diff_percent_posi[2], bins = bins, linecolor = :transparent, yaxis=(:log10), title = "br$(br_arr[2])",
                       xlabel = "rate of improvements (cost reduceed/base cost)", ylabel = "amount of hours", titlefontsize=20,
                       right_margin = 10px, xlim = max_base_diff_percent)
            p9 = histogram(base_diff_percent_posi[3], bins = bins, linecolor = :transparent, yaxis=(:log10), title = "br$(br_arr[3])",
                       xlabel = "rate of improvements (cost reduceed/base cost)", ylabel = "amount of hours", titlefontsize=20,
                       right_margin = 10px, xlim = max_base_diff_percent)
            p10 = histogram(base_diff_percent_posi[4], bins = bins, linecolor = :transparent, yaxis=(:log10), title = "br$(br_arr[4])",
                       xlabel = "rate of improvements (cost reduceed/base cost)", ylabel = "amount of hours", titlefontsize=20,
                       right_margin = 10px, xlim = max_base_diff_percent)
            #Plots.link_axes!()
            p_first = plot(p3, p4, p5, p6, layout = (4,1), link = :all,legend = false, size = (500, 800))
            savefig(p_first, "$(output_folder)/opt_improve_plot_comp_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).pdf")

            # title = "amount of hours with each rate of improvements"
            p_second = plot(p7, p8, p9, p10, layout = (4,1), link = :x,legend = false, size = (500, 800) )
            savefig(p_second, "$(output_folder)/opt_rate_plot_comp_($(month_st).$(day_st).$(period_st)__$(month_en).$(day_en).$(period_en)).pdf")
            #for br_ind = 1:num_br
            #end
      end

      println("The objective value of the current period")
end


text:
module DarcyTests

using Test
using Gridap
import Gridap: ∇, divergence
using LinearAlgebra

u(x) = VectorValue(2*x[1],x[1]+x[2])

divergence(::typeof(u)) = (x) -> 3

p(x) = x[1]-x[2]

∇p(x) = VectorValue(1,-1)

∇(::typeof(p)) = ∇p

f(x) = u(x) + ∇p(x)

domain = (0,1,0,1)
partition = (4,4)
order = 1
model = CartesianDiscreteModel(domain,partition)

V = FESpace(model,ReferenceFE(raviart_thomas,Float64,order),conformity=:Hdiv,
  dirichlet_tags=[5,6])

Q = FESpace(model,ReferenceFE(lagrangian,Float64,order); conformity=:L2)

U = TrialFESpace(V,u)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

trian = Triangulation(model)
degree = 2
dΩ = Measure(trian,degree)

neumanntags = [7,8]
btrian = BoundaryTriangulation(model,tags=neumanntags)
degree = 2*(order+1)
dΓ = Measure(btrian,degree)
nb = get_normal_vector(btrian)

a((u, p),(v, q)) = ∫( u⋅v - (∇⋅v)*p + q*(∇⋅u) )*dΩ

b(( v, q)) = ∫( v⋅f + q*(∇⋅u))*dΩ - ∫((v⋅nb)*p )*dΓ

op = AffineFEOperator(a,b,X,Y)
xh = solve(op)
uh, ph = xh

eu = u - uh
ep = p - ph

l2(v) = sqrt(sum(∫(v⋅v)*dΩ))
h1(v) = sqrt(sum(∫(v*v + ∇(v)⋅∇(v))*dΩ))

eu_l2 = l2(eu)
ep_l2 = l2(ep)
ep_h1 = h1(ep)

tol = 1.0e-9
@test eu_l2 < tol
@test ep_l2 < tol
@test ep_h1 < tol

end # module


text:
@testset "z" begin
    @test ScientificTest.z(2) == 6
    @test ScientificTest.z(3) == 9
end


text:
"Generic PDDL planning domain."
@kwdef mutable struct GenericDomain <: Domain
    name::Symbol # Name of domain
    requirements::Dict{Symbol,Bool} = Dict() # PDDL requirements used
    typetree::Dict{Symbol,Vector{Symbol}} = Dict() # Types and their subtypes
    datatypes::Dict{Symbol,Type} = Dict() # Non-object data types
    constants::Vector{Const} = [] # List of constants
    constypes::Dict{Const,Symbol} = Dict() # Types of constants
    predicates::Dict{Symbol,Signature} = Dict() # Dictionary of predicates
    functions::Dict{Symbol,Signature} = Dict() # Dictionary of function declarations
    funcdefs::Dict{Symbol,Any} = Dict() # Dictionary of function definitions
    axioms::Dict{Symbol,Clause} = Dict() # Axioms / derived predicates
    actions::Dict{Symbol,Action} = Dict() # Action definitions
    _extras::Dict{Symbol,Any} # Extra fields
end

function GenericDomain(name::Symbol, header::Dict{Symbol,Any}, body::Dict{Symbol,Any})
    h_extras = filter(item -> !(first(item) in fieldnames(GenericDomain)), header)
    b_extras = filter(item -> !(first(item) in fieldnames(GenericDomain)), body)
    extras = merge!(h_extras, b_extras)
    header = filter(item -> first(item) in fieldnames(GenericDomain), header)
    axioms = Clause[get(body, :axioms, []); get(body, :deriveds, [])]
    body = filter(item -> first(item) in fieldnames(GenericDomain), body)
    body[:axioms] = Dict(ax.head.name => ax for ax in axioms)
    body[:actions] = Dict(act.name => act for act in body[:actions])
    return GenericDomain(;name=name, _extras=extras, header..., body...)
end

Base.getproperty(d::GenericDomain, s::Symbol) =
    hasfield(GenericDomain, s) ? getfield(d, s) : d._extras[s]

Base.setproperty!(d::GenericDomain, s::Symbol, val) =
    hasfield(GenericDomain, s) && s != :_extras ?
        setfield!(d, s, val) : setindex!(d._extras, val, s)

function Base.propertynames(d::GenericDomain, private=false)
    if !private
        tuple(fieldnames(GenericDomain)..., keys(d._extras)...)
    else
        tuple(filter(f -> f != :_extras, fieldnames(GenericDomain))...,
              keys(d._extras)...)
    end
end

Base.copy(domain::GenericDomain) = deepcopy(domain)

get_name(domain::GenericDomain) = domain.name

get_requirements(domain::GenericDomain) = domain.requirements

get_typetree(domain::GenericDomain) = domain.typetree

get_datatypes(domain::GenericDomain) = domain.datatypes

get_constants(domain::GenericDomain) = domain.constants

get_constypes(domain::GenericDomain) = domain.constypes

get_predicates(domain::GenericDomain) = domain.predicates

get_functions(domain::GenericDomain) = domain.functions

get_funcdefs(domain::GenericDomain) = domain.funcdefs

get_fluents(domain::GenericDomain) = merge(domain.predicates, domain.functions)

get_axioms(domain::GenericDomain) = domain.axioms

get_actions(domain::GenericDomain) = domain.actions


text:
module TestTikZWiringDiagrams

using Test

using Catlab.Doctrines
using Catlab.WiringDiagrams
using Catlab.Graphics
import Catlab.Graphics: TikZ
import Catlab.Graphics.TikZWiringDiagrams.Defaults

# We can't test that the pictures look right, but we can test that they exist!
is_pic(obj) = isa(obj, TikZ.Picture)

A, B = Ob(FreeCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)
@test is_pic(to_tikz(f))
@test is_pic(to_tikz(compose(f,g)))
@test is_pic(to_tikz(compose(id(A), f, id(B))))

A, B = Ob(FreeDaggerCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)
@test is_pic(to_tikz(dagger(f)))
@test is_pic(to_tikz(dagger(compose(f,g))))

A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)
@test is_pic(to_tikz(otimes(f,g)))
@test is_pic(to_tikz(compose(otimes(f,g),braid(B,A))))

A, B = Ob(FreeBiproductCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)
@test is_pic(to_tikz(compose(mcopy(A), otimes(f,f), mmerge(B))))
@test is_pic(to_tikz(compose(create(A), f, g, delete(A))))

A, B = Ob(FreeCompactClosedCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)
h = Hom(:h, otimes(A,B), otimes(A,B))
@test is_pic(to_tikz(dunit(A)))
@test is_pic(to_tikz(dcounit(A)))
trace = compose(
  f,
  otimes(dunit(A), id(B)),
  otimes(id(dual(A)), h),
  otimes(braid(dual(A),A), id(B)),
  otimes(dcounit(A), id(B)),
  g,
)
@test is_pic(to_tikz(trace))

end


text:
# Save Gradcam

"""
    save_gradcam(gradient, original_image_path, grad_file_name, heatmap_file_name, combined_file_name)

Function to save the gradient calculated using `gradcam` as a `gray image`, `gradient heatmap` and the `combined heatmap` with different transparency.

# Arguments
 - `gradient`: Gradient calculated using the `viz_gradcam` function
 - `original_image_path`: Path to the original image file, for which the gradient is calculated.
 - `grad_file_name`: Filename with path, to which the gradient image is to be saved as a `gray` image
 - `heatmap_file_name`: Filename to save the gradient heatmap
 - `combined_file_name`: Filename to save the gradient the combined heatmap with original image.

"""
function save_gradcam(gradient, original_image_path, grad_file_name, heatmap_file_name, combined_file_name)
  gradient = max.(gradient, zero(gradient))
  gradient = permutedims((gradient .- minimum(gradient))/(maximum(gradient) - minimum(gradient)), [3, 2, 1])
  img = Gray.(colorview(RGB{eltype(gradient)},gradient))
  img = imresize(img, 224, 224)
  display(img)
  h1 = plot(heatmap(channelview(float.(augment(img, FlipY()))), color = :rainbow))
  try
    display(h1)
  catch
    @info("The heatmap could not be displayed. The file will be saved at $heatmap_file_name")
  end
  mapped = float.(channelview(imresize(load(original_image_path), (224, 224)))) * 0.9 .+ reshape(float.(img), 1, 224, 224) * 2
  mapped .-= minimum(mapped)
  mapped ./= maximum(mapped)
  h2 = plot(heatmap(channelview(float.(Gray.(augment(colorview(RGB{eltype(mapped)}, mapped), FlipY())))), color=:rainbow))
  try
    display(h2)
  catch
    @info("The heatmap could not be displayed. The file will be saved at $combined_file_name")
  end
  save(grad_file_name, img)
  try
    savefig(h1, heatmap_file_name)
    savefig(h2, combined_file_name)
  catch
    @info("Encountered Errors while trying to save heatmaps. So the heatmaps are being returned")
    @info("They need to be saved manually")
    (h1, h2)
  end
end

# GradCAM

normalize_grad(grad) = grad / (sqrt(mean(abs2.(grad))) + eps())

function save_gradient(x, tracker)
    tracker["grad"] = x |>cpu
    x
end

function tracked_model_gc(x, model, layer, tracker)
    x = model[1:layer](x) |> gpu
    tracker["fwd"] = x |> cpu
    model[layer+1:end](Zygote.hook(x->save_gradient(x, tracker), x)) |>gpu
end

"""
    viz_gradcam(img, model, layer; top_k=1)

Layer attribution method that computes the gradients of the target output with respect to the given `layer`, averages for each output channel (dimension 2 of output), and multiplies the average gradient for each channel by the layer activations. The results are summed over all channels.

It returns an array of tuples of (`gradCAM`, `prediction_probability`, `class_index`) for `top_k` predicted classes.
"""
function viz_gradcam(img, model, layer; top_k=1)
    model_ = model[1:end-1]
    tracker = Dict()
    set_training_false()
    preds, back = Zygote.pullback(x -> tracked_model_gc(x, model_, layer, tracker), img)
    set_training_true()
    probs = softmax(preds)
    prob, inds = get_topk(probs, k=top_k)
    grads = []
    for (i, idx) in enumerate(inds)
        back(one_hot_encode(preds, idx))
        weights = reshape(maximum(normalize_grad(tracker["grad"]), dims=[1, 2, 4]), 1, 1, size(tracker["fwd"])[3], 1)
        cam = ones(size(tracker["fwd"])[1:2])
        cam += dropdims(sum(weights .* tracker["fwd"],dims=[3, 4]), dims=(3, 4))
        push!(grads, (cam, probs[i], idx))
    end
    grads
end

# Guided GradCAM

"""
    viz_guidedgradcam(img, model, layer; top_k=1)

It computes the element-wise product of guided backpropagation attributions (evaluated using `viz_guidedbackprop`) with resized (layer) GradCAM attributions (computed using `viz_gradcam`). GradCAM attributions are computed with respect to a given `layer`, and attributions are resized to match the input size.

It returns an array of tuples of (`guidedgradcam`, `prediction_probability`, `class_index`) for `top_k` predicted classes.

Also see: [`viz_guidedbackprop`](@ref), [`viz_guidedgradcam`](@ref).
"""
function viz_guidedgradcam(img, model, layer; top_k=1)
    # Ideally we should not be doing the forward pass twice
    gcam_grads = viz_gradcam(img, model, layer, top_k=1)
    backprop_grads = viz_guidedbackprop(img, model; top_k=1)
    grads = []
    for i in 1:top_k
      normalized_camgrads = gcam_grads[i][1] .- minimum(gcam_grads[i][1])
      normalized_camgrads ./= maximum(normalized_camgrads)
      normalized_camgrads = reshape(channelview(float.(imresize(Gray.(normalized_camgrads), 224, 224))), 224, 224, 1, 1)
      normalized_bpropgrads = backprop_grads[i][1][1] .- minimum(backprop_grads[i][1][1])
      normalized_bpropgrads ./= maximum(normalized_bpropgrads)
      push!(grads, (normalized_camgrads .* normalized_bpropgrads, gcam_grads[i][2], gcam_grads[i][3]))
    end
    grads
end


text:
using Logging
using SIAN

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

# to run this, one should comment the NFkB part in the benchmark file
include("benchmarks.jl")
NUM_RUNS = 5

index_to_run = 1

bmark = benchmarks[index_to_run]
@info "Processing $(bmark[:name])"
ode = bmark[:ode]

total_time = 0

for i in 1:(NUM_RUNS + 1)
    @info "Run numer $i"
    rtime = @elapsed identifiability_ode(ode, get_parameters(ode))
    @info "Runtime: $rtime"
    if i > 1
        global total_time += rtime
    end
end

@info total_time
@info "Average: $(total_time / NUM_RUNS)"


text:
using Term.progress
import Term.progress: AbstractColumn, getjob, get_columns
import Term: install_term_logger, uninstall_term_logger

using ProgressLogging

@testset "\e[34mProgress - jobs" begin

    pbar = ProgressBar()

    j1 = addjob!(pbar; description="test", N=10)

    @test j1.id == 1
    @test j1.N == 10
    @test j1.i == 1
    @test j1.description == "test"
    @test j1.started == true
    @test typeof(j1.columns) == Vector{AbstractColumn}

    update!(j1)
    @test j1.i == 2
    update!(j1; i=10)
    @test j1.i == 12

    @test getjob(pbar, j1.id).id == j1.id

    removejob!(pbar, j1)
    @test length(pbar.jobs) == 0
end

@testset "\e[34mProgress basic" begin

    @test_nothrow begin
        pbar = ProgressBar()
        with(pbar) do
            job = addjob!(pbar; N=10)
            for i in 1:10
                update!(job)
                sleep(.001)
            end
        end
    end

    @test_nothrow begin
        pbar2 = ProgressBar(; transient=true)
        with(pbar2) do
            job = addjob!(pbar2; N=10)
            job2 = nothing
            for i in 1:10
                if i == 50
                    job2 = addjob!(pbar2; N=10)
                end

                if i >= 50
                    update!(job2)
                end

                if i == 75
                    removejob!(pbar2, job2)
                end

                update!(job)
                sleep(.001)
            end
        end
    end

    @test_nothrow begin
        @track for i in 1:10
            sleep(0.001)
        end
    end
end

@testset "\e[34mProgress columns" begin

    for colinfo in (:minimal, :default, :spinner, :detailed)
        pbar = ProgressBar(columns=colinfo)
        @test pbar.columns == get_columns(colinfo)
        @test typeof(pbar.columns) == Vector{DataType}

        job = addjob!(pbar; N= colinfo == :spinner ? nothing : 10)
        @test typeof(job.columns) == Vector{AbstractColumn}

        with(pbar) do
            for i in 1:10
                update!(job)
                sleep(0.01)
            end
        end
    end

    colkwargs = Dict(:DescriptionColumn=>Dict(:style=>"red"))
    pbar = ProgressBar(;columns_kwargs=colkwargs)
    job = addjob!(pbar; N=10)
    @test job.columns[1].segments[1].text == "\e[31mRunning...\e[39m"
end

# @testset "\e[34mProgress ProgressLogging" begin

#     install_term_logger()

#     @test_nothrow begin
#         @progress "inner... $i" for j in  1:10
#             sleep(0.01)
#         end
#     end
# end


text:
# # Finding boundaries

# The output of a `wombling` operation can be used to pick boundaries, *i.e.*
# areas where the values on the grid transition sharply. We will illustrate this
# with a simple example of a three-patch landscape.

using SpatialBoundaries
using StatsPlots
default(; dpi=500, size=(600, 600), aspectratio=1, c=:batlow, frame=:box)

# Let's create a landscape with two values:

A = rand(Float64, 200, 150);
A[1:80, 1:85] .+= 5.0;
A[110:end, 130:end] .+= 3.0;

# We can check out what this landscape looks likes:

heatmap(A)

# We can apply wombling to this landscape, assuming that all cells have the same
# size:

W = wombling(A);

# Let's look at the rate of change:

heatmap(W.m; c=:nuuk)

# Picking the boundaries is done by passing the wombling output to the
# `boundaries` function, with a specific threshold giving the proportion of
# points that should be retained as part of the boundaries. Checking what the
# effect of this threshold is would be a good idea:

thresholds = LinRange(0.0, 0.2, 200)
patches = [length(boundaries(W; threshold=t)) for t in thresholds]

plot(thresholds, log1p.(patches), aspectratio=:none)
xaxis!("Threshold", (0., 0.2))
yaxis!("log(boundary patches + 1)", (0., 9.))

# Let's eyeball this as a 0.01, and see how the patches are distributed.

# Another way we can look at the boundaries is to see *when* a patch is
# considered to be a boundary. To do so we will create an empty matrix, and fill
# each position with the lowest threshold at which it is considered to be a
# boundary:

b = similar(W.m)

for t in reverse(LinRange(0.0, 1.0, 200))
    b[boundaries(W; threshold=t)] .= t
end

heatmap(b, c=:tofino, clim=(0,1))

# This also suggests that we will get well delineated patches for low values of
# the threshold.

B = boundaries(W; threshold=0.01);

# In the following figure, cells identified as candidate boundaries are marked
# in white:

heatmap(A)
scatter!([(reverse(x.I)) for x in B], leg=false, msw=0, c=:white)

# We can see that the boundaries of the patches have been well identified!


text:
"""
    SolverTrace(num_steps, i, print_interval, io, progress, columns, callbacks)

Represents a solver trace for `num_steps` iterations, `i` being the
current one. Each time [`next!`](@ref) is called, `i` is incremented
by one, and if a multiple of `print_interval` is reached, a line is
printed in the trace. Optionally (default), a `ProgressMeter.Progress`
is shown as well, providing progress information between solver trace
printouts.
"""
mutable struct SolverTrace{P,TCs,Callbacks,IOt<:IO}
    num_steps::Int
    i::Int
    print_interval::Int
    io::IOt
    progress::P
    columns::TCs
    callbacks::Callbacks
end

"""
    SolverTrace(num_steps[, column=CurrentStep(num_steps), columns...;
                            io=stdout, num_printouts=10, progress_meter=true,
                            print_interval, kwargs...])

Construct a [`SolverTrace`](@ref) for `num_steps` iterations, with the
`columns` that are to appear in the solver trace. The solver trace is
printed to `io`, `num_printouts` lines will appear; alternatively, the
`print_interval` can be specified directly. Below the solver trace, a
progressbar is displayed, unless `!progress_meter`. `kwargs...` are
passed on to `ProgressMeter.Progress`.
"""
function SolverTrace(num_steps::Int,
                     column::TraceColumn = CurrentStep(num_steps),
                     columns::TraceColumn...;
                     io::IO=stdout,
                     num_printouts::Integer=min(num_steps,10),
                     progress_meter::Bool=true,
                     callbacks=(),
                     print_interval = num_printouts > num_steps ? 1 : num_steps÷num_printouts,
                     kwargs...)
    columns = (column,columns...)

    progress = if progress_meter
        Progress(num_steps; kwargs...)
    else
        nothing
    end
    SolverTrace(num_steps, 0, print_interval, io, progress, columns, callbacks)
end

clear_current_line(io::IO=stdout) = print(io, "\r\u1b[K")

"""
    print_header(s)

Print the associated header of the [`SolverTrace`](@ref) `s`.
"""
print_header(s::SolverTrace) =
    println(s.io,
            crayon"underline",
            join(map(column -> column.header, s.columns), " "),
            crayon"reset")

print_header(::Nothing) = nothing

Formatting.format(s::SolverTrace) =
    join(map(column -> format(column, s.i), s.columns), " ")

next!(::Nothing; kwargs...) = nothing
next!(progress::Progress; kwargs...) = ProgressMeter.next!(progress; kwargs...)

"""
    next!(fun, s::SolverTrace; kwargs...)

Signal to the [`SolverTrace`](@ref) `s` that the next iteration has
been performed, increasing the internal counter `s.i`. If a multiple
of `s.print_interval` is reached, a line in the solver trace is
printed, formatting each and every trace column, the callback function
`fun` is called, as well as all the registered callbacks.

`kwargs` are passed on to `ProgressMeter.next!` to allow
e.g. `showvalues` functionality.
"""
function next!(fun::Function, s::SolverTrace; kwargs...)
    s.i += 1
    if s.i%s.print_interval == 0 || s.i==1
        fun()
        # Carriage return & clearing has to be done on stderr since
        # ProgressMeter prints to that stream; this is important on
        # Windows.
        isnothing(s.progress) || clear_current_line(stderr)
        println(s.io, format(s))
        foreach(f -> f(s.i), s.callbacks)
    end
    next!(s.progress; kwargs...)
end
next!(s::SolverTrace; kwargs...) = next!(() -> nothing, s; kwargs...)

"""
    ColumnSeparator(fmt, header)

Simple separator between columns, that by default simply prints a
pipe.
"""
struct ColumnSeparator <: TraceColumn
    fmt::FormatExpr
    header::String
end

ColumnSeparator() = ColumnSeparator(FormatExpr("│"), " ")
(c::ColumnSeparator)(i::Integer) = ()

export SolverTrace, print_header, ColumnSeparator


text:
#### Woodbury matrices ####
# This type provides support for the Woodbury matrix identity
type Woodbury{T} <: AbstractMatrix{T}
    A
    U::Matrix{T}
    C
    Cp
    V::Matrix{T}
    tmpN1::Vector{T}
    tmpN2::Vector{T}
    tmpk1::Vector{T}
    tmpk2::Vector{T}

    function Woodbury(A, U::Matrix{T}, C, V::Matrix{T})
        N = size(A, 1)
        k = size(U, 2)
        if size(A, 2) != N || size(U, 1) != N || size(V, 1) != k || size(V, 2) != N
            throw(DimensionMismatch(""))
        end
        if k > 1
            if size(C, 1) != k || size(C, 2) != k
                throw(DimensionMismatch(""))
            end
        end
        Cp = inv(inv(C) + V*(A\U))
        # temporary space for allocation-free solver
        tmpN1 = Array(T, N)
        tmpN2 = Array(T, N)
        tmpk1 = Array(T, k)
        tmpk2 = Array(T, k)
        # don't copy A, it could be huge
        new(A, copy(U), copy(C), Cp, copy(V), tmpN1, tmpN2, tmpk1, tmpk2)
    end
end

Woodbury{T}(A, U::Matrix{T}, C, V::Matrix{T}) = Woodbury{T}(A, U, C, V)
Woodbury{T}(A, U::Vector{T}, C, V::Matrix{T}) = Woodbury{T}(A, reshape(U, length(U), 1), C, V)

size(W::Woodbury) = size(W.A)

function show(io::IO, W::Woodbury)
    println(io, summary(W), ":")
    print(io, "A:\n", W.A)
    print(io, "\nU:\n")
    print_matrix(io, W.U)
    if isa(W.C, Matrix)
        print(io, "\nC:\n")
        print_matrix(io, W.C)
    else
        print(io, "\nC: ", W.C)
    end
    print(io, "\nV:\n")
    print_matrix(io, W.V)
end

full{T}(W::Woodbury{T}) = convert(Matrix{T}, W)
convert{T}(::Type{Matrix{T}}, W::Woodbury{T}) = full(W.A) + W.U*W.C*W.V

function similar(W::Woodbury, T, dims::Dims)
    (length(dims) == 2 && dims[1] == dims[2] ) || throw(DimensionMismatch(""))
    n = size(W, 1)
    k = size(W.U, 2)
    Woodbury{T}(similar(W.A), Array(T, n, k), Array(T, k, k), Array(T, k, n))
end

copy(W::Woodbury) = Woodbury(W.A, W.U, W.C, W.V)

## Woodbury matrix routines ##

*(W::Woodbury, B::StridedVecOrMat)=W.A*B + W.U*(W.C*(W.V*B))

function \(W::Woodbury, R::StridedVecOrMat)
    AinvR = W.A\R
    AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
end

det(W::Woodbury)=det(W.A)*det(W.C)/det(W.Cp)

function A_ldiv_B!(W::Woodbury, B::AbstractVecOrMat)
    length(B) == size(W, 1) || throw(DimensionMismatch("Vector length $(length(B)) must match matrix size $(size(W,1))"))
    copy!(W.tmpN1, B)
    Alu = lufact(W.A) # Note. This makes an allocation (unless A::LU). Alternative is to destroy W.A.
    A_ldiv_B!(Alu, W.tmpN1)
    A_mul_B!(W.tmpk1, W.V, W.tmpN1)
    A_mul_B!(W.tmpk2, W.Cp, W.tmpk1)
    A_mul_B!(W.tmpN2, W.U, W.tmpk2)
    A_ldiv_B!(Alu, W.tmpN2)
    for i = 1:length(W.tmpN2)
        @inbounds B[i] = W.tmpN1[i] - W.tmpN2[i]
    end
    B
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia