text:
# Algorithms for Regret Minimization in a Piecewise Stationary MAB with Bernoulli rewards
# julia 1.0.5

# All algorithms take as an input Table, where Table[i,t] is the rewards of arm i at time t
# All algorithms output ChosenArms,ReceivedRewards,episode (=nb of episodes),ChangePoints (=position of the breakpoints)
# (the last one is always empty for algorithm which do not explicitely perform restarts)

using Distributions

### USEFUL FUNCTIONS

function randmax(vector,rank=1)
   # returns an argmax uniformly at random among all maximizers
   # (integer, not CartesianIndex)
   vector=vec(vector)
   Sorted=sort(vector,rev=true)
   m=Sorted[rank]
   Ind=findall(x->x==m,vector)
   index=Ind[floor(Int,length(Ind)*rand())+1]
   return (index)
end

function sampleDiscrete(probas)
   # return a sample from a distribution on {1,...,K} such that P(X=k)=probas[k]
   cumprobas = cumsum(probas)
   if (cumprobas[end]<0.99)
      print("something is wrong!\n")
   end
   cumprobas[end]=1.
   u = rand()
   k = 1
   while (u>cumprobas[k])
      k = k+1
   end
   return k
end

function Normalize(weights)
   # cumpute weights/sum(weights), stable for "big" weights
   K = length(weights)
   probas = zeros(Float32,K)
   for k in 1:K
      exptot = 0
      logweight = log(weights[k])
      for j in 1:K
         if (j!=k)
            exptot = exptot + exp(log(weights[j])-logweight)
         end
      end
      probas[k]=1/(1+exptot)
   end
   return probas
end

function kl(p,q)
   # binary relative entropy
   res=0
   if (p!=q)
      if (p<=0) p = eps() end
      if (p>=1) p = 1-eps() end
      res=(p*log(p/q) + (1-p)*log((1-p)/(1-q)))
   end
   return res
end

function klUp(p,level)
   # KL upper confidence bound:
   # return uM>p such that d(p,uM)=level
   lM = p
   uM = min(min(1,p+sqrt(level/2)),1)
   for j = 1:16
      qM = (uM+lM)/2
      if kl(p,qM) > level
         uM= qM
      else
         lM=qM
      end
   end
   return uM
end

function klSG(p,q)
   # kl-diverge corresponding to a 1/4-subGaussian approximation
   return 2*(p-q)^2
end

function DrawTable(Problem)
   # generate a Table of rewards from a matrix of mean rewards mu[a,t]
   (K,T)=size(Problem)
   Table = Int.(rand(K,T).<Problem)
end

function ComputeCumRegret(Instance,Choices)
   # compute the vector of cumulative regret on a bandit instance based on the successive selections of an algorithm stored in Choices
   K,T = size(Instance)
   regret = 0
   Regret = zeros(Float32,1,T)
   for t in 1:T
      regret+=maximum(Instance[:,t])-Instance[Choices[t],t]
      Regret[t]=regret
   end
   return Regret
end

### ALGORITHMS

## AdSwitch [Auer et al. 19]

function AdSwitch(Table,C=1,subSample=50,subSample2 = 20,verbose="off")
   # table contains a table of rewards for the K arms up to horizon T
   (K,T)=size(Table)
   # index of the episode
   episode = 1
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   ChangePoints = []
   t = 0 # current number of sample
   while (t<T)
      SUMS = zeros(K,T-t,T-t) # array containing the S[a,s,t]
      DRAWS = zeros(K,T-t,T-t) # array containing the N[a,s,t]
      Bad_Gaps = zeros(K) # gaps of the bad arms
      Bad_Means = zeros(K) # means of the bad arms
      tloc = 0 # number of samples in current episode
      TotalNumber=zeros(Int,K) # total number of selections in episode
      # initialize set of good and bad arms
      Good = collect(1:K)
      Bad = []
      # initialize sampling obligations
      Obligations = zeros(Int,K)
      newEpisode = false # should we start a new episode?
      # start the episode
      while (!newEpisode)&&(t<T)
         # form the set of candidate arms to choose from (good arm + sampling obligations)
         Candidate = copy(Good)
         append!(Candidate,[i for i in 1:K if Obligations[i]>0])
         # draw an arm (least sampled among candidates)
         I = randmax(-TotalNumber[Candidate])
         rew = Table[I,t+1]
         TotalNumber[I]=+1
         ChosenArms[t+1]=I
         ReceivedRewards[t+1]=rew
         # update everything
         if (Obligations[I]>0)
            # if the arm was sampled due to obligations, recrease the remaining time to sample
            Obligations[I]=Obligations[I]-1
         end
         t=t+1
         tloc=tloc+1
         for i in 1:K
            if (i!=I)&&(tloc>1)
               DRAWS[i,1:(tloc-1),tloc]=DRAWS[i,1:(tloc-1),tloc-1]
               SUMS[i,1:(tloc-1),tloc]=SUMS[i,1:(tloc-1),tloc-1]
            else
               DRAWS[i,tloc,tloc]=1
               SUMS[i,tloc,tloc]=rew
               if (tloc>1)
                  DRAWS[i,1:(tloc-1),tloc]=DRAWS[i,1:(tloc-1),tloc-1].+1
                  SUMS[i,1:(tloc-1),tloc]=SUMS[i,1:(tloc-1),tloc-1].+rew
               end
            end
         end
         # updating the set of good arms
         for s in 1:tloc
            cand = [a for a in Good if DRAWS[a,s,tloc]>1]
            candmeans = [SUMS[a,s,tloc]/DRAWS[a,s,tloc] for a in cand]
            if (length(candmeans)>0)
               mumax = maximum(candmeans)
               for i in 1:length(cand)
                  arm = cand[i]
                  if mumax - candmeans[i] > sqrt(C*log(T)/(DRAWS[arm,s,tloc]-1))
                     if verbose=="on"
                        print("removed arm $(arm) from Good at t=$(t)!\n")
                     end
                     # remove arm from Good
                     deleteat!(Good, findfirst(isequal(arm), Good))
                     # add arm to bad and store its gap and mean
                     append!(Bad,arm)
                     Bad_Gaps[arm]= mumax - candmeans[i]
                     Bad_Means[arm] = candmeans[i]
                  end
               end
            end
         end
         # perform tests
         if (tloc % subSample ==1)
            # check whether a bad arm has changed
            check=0
            idchange=0
            s=1
            while (s<tloc)&&(check==0)
               for badarm in Bad
                  draws=DRAWS[badarm,s,tloc]
                  if (draws>1)&&(abs(SUMS[badarm,s,tloc]/draws-Bad_Means[badarm]) > Bad_Gaps[badarm]/4 + sqrt(2*log(T)/draws))
                     newEpisode = true
                     check +=1
                     idchange=badarm
                  end
               end
               s = s+1
            end
            if (check==0)
               # check whether a good arm has changed
               s = 1
               while (s < tloc)&&(check==0)
                  for s1 in [j for j in 1:tloc if (j % subSample2 ==1)]
                     for s2 in [j for j in s1:tloc if (j % subSample2 ==1)]
                        for goodarm in Good
                           draws1=DRAWS[goodarm,s1,s2]
                           draws2=DRAWS[goodarm,s,tloc]
                           if (draws1>1)&&(draws2>1)&&(abs(SUMS[goodarm,s1,s2]/draws1 - SUMS[goodarm,s,tloc]/draws2)>sqrt(2*log(T)/draws1)+sqrt(2*log(T)/draws2))
                              newEpisode = true
                              check +=1
                              idchange=goodarm
                           end
                        end
                     end
                  end
                  s = s+1
               end
            end
            if (check>0)
               episode+=1
               ChangePoints=append!(ChangePoints,t)
               if verbose=="on"
                  print("detected a change on $(idchange) at t=$(t)\n\n")
               end
            end
         end
         # possibly add some new sampling obligation
         for badarm in Bad
            i=1
            while (1/(2^i)>= Bad_Gaps[badarm]/16)
               if (rand() < (1/(2^i))*sqrt(episode/(K*T*log(T))))
                  n =floor(Int,2^(2*i+1)*log(T))
                  if (verbose=="on")
                     print("add a sampling obligation for arm $(badarm) at t=$(t) of length $(n)\n")
                  end
                  # update the sampling obligation
                  Obligations[badarm]=max(Obligations[badarm],n)
               end
               i+=1
            end
         end
      end
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

## Variants of GLR-klUCB

# GLRklUCB with Global Restart [Besson et al. 20]
# recommended tuning : alpha0 = sqrt(K*log(T)/T) and delta = 1/sqrt(T)

function GLRklUCBGlobal(Table,beta,klChange=kl,alpha0=0.1,delta=0.05,alphatype="auto",subSample=10,subSample2 = 5,verbose="off")
   # beta, klChange : threshold and kl-divergence function used in the GLRT test
   # parameter of the exploration sequence alpha0 (such that the exploration probability in episode k is alpha0*\sqrt(k)), error probability for the test (delta)
   # allow for a constant exploration probability alpha0 for alphatype="constant"
   (K,T)=size(Table)
   # index of the episode
   episode = 1 # nb of breakpoints detected
   ChangePoints = []
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   t = 0 # current number of sample
   while (t<T)
      SUMS = [[] for a=1:K] # array containing the cumulated sums of rewards
      tloc = 0 # number of samples in current episode
      TotalNumber=zeros(Int,K) # total number of selections in current episode
      TotalSum=zeros(Float32,K) # sum of rewards in current episode
      newEpisode = false # should we start a new episode?
      alpha = alpha0*sqrt(episode) # exploration probability for the episode
      if (alphatype == "constant")
         alpha = alpha0
      end
      ExploRange = floor(Int,K/alpha)
      # start the episode
      while (!newEpisode)&&(t<T)
         # default exploratory action
         I = (t % ExploRange)
         if (I ==0)||(I>K)
            # compute KL-UCB indices if no exploration
            indices = ones(K)
            for i in 1:K
               if (TotalNumber[i]>0)
                  indices[i]=klUp(TotalSum[i]/TotalNumber[i],log(tloc)/TotalNumber[i])
               end
            end
            I = randmax(indices)
         end
         # get the reward
         rew = Table[I,t+1]
         ChosenArms[t+1] = I
         ReceivedRewards[t+1]=rew
         # update everything
         t=t+1
         tloc=tloc+1
         TotalNumber[I]+=1
         TotalSum[I]+=rew
         append!(SUMS[I],TotalSum[I])
         # has the selected arm changed?
         check = 0
         if (tloc % subSample == 1)
            s = 1
            nb = TotalNumber[I]
            sums = SUMS[I]
            while (s < nb)&&(check==0)
               if (s % subSample2 ==1)
                  draw1 = s
                  draw2 = nb-s
                  mu1 = sums[s]/draw1
                  mu2 = (sums[nb]-sums[s])/draw2
                  mu = sums[nb]/nb
                  if (draw1*klChange(mu1,mu)+draw2*klChange(mu2,mu))>beta(nb,delta)
                     newEpisode = true
                     check +=1
                  end
               end
               s=s+1
            end
            if (check>0)
               episode+=1
               ChangePoints=append!(ChangePoints,t)
               if verbose=="on"
                  print("detected a change on $(I) at t=$(t)\n\n")
               end
            end
         end
      end
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

# GLRTklUCB with Local Restarts [Besson et al. 20]
# recommended tuning : alpha0 = sqrt(K*log(T)/T) and delta = 1/sqrt(T)

function GLRklUCBLocal(Table,beta,klChange=kl,alpha0=0.,delta=0.,alphatype="auto",subSample=10,subSample2 = 10,verbose="off")
   # beta, klChange : threshold and kl-divergence function used in the GLRT test
   # parameter of the exploration sequence alpha0 (such that the exploration probability in episode k is alpha0*\sqrt(k)), error probability for the test (delta)
   # allow for a constant exploration probability alpha0 for alphatype="constant"
   (K,T)=size(Table)
   # index of the episode
   episode = 1 # nb of breakpoints detected
   ChangePoints = []
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   t = 0 # current number of sample
   SUMS = [[] for a=1:K] # array containing the cumulated sums of rewards
   tloc = zeros(Int,K) # elapsed time since last restart for each arm
   TotalNumber=zeros(Int,K) # total number of selections since last restart
   TotalSum=zeros(Float32,K) # sum of rewards since last restart
   while (t<T)
      newEpisode = false # should we do a restart on the last arm selects?
      alpha = alpha0*sqrt(episode) # exploration probability for the episode
      if (alphatype == "constant")
         alpha = alpha0
      end
      ExploRange = floor(Int,K/alpha)
      # start the episode
      while (!newEpisode)&&(t<T)
         # default exploratory action
         I = (t % ExploRange)
         if (I ==0)||(I>K)
            # compute KL-UCB indices and choose the argmax if no exploration
            indices = ones(K)
            for i in 1:K
               if (TotalNumber[i]>0)
                  indices[i]=klUp(TotalSum[i]/TotalNumber[i],log(tloc[i])/TotalNumber[i])
               end
            end
            I = randmax(indices)
         end
         # get the reward
         rew = Table[I,t+1]
         ChosenArms[t+1] = I
         ReceivedRewards[t+1]=rew
         # update everything
         t=t+1
         tloc=tloc.+1 # update ALL the local times
         TotalNumber[I]+=1
         TotalSum[I]+=rew
         append!(SUMS[I],TotalSum[I])
         # has the selected arm changed?
         check = 0
         if (tloc[I] % subSample == 1)
            s = 1
            nb = TotalNumber[I]
            sums = SUMS[I]
            while (s < nb)&&(check==0)
               if (s % subSample2 ==1)
                  draw1 = s
                  draw2 = nb-s
                  mu1 = sums[s]/draw1
                  mu2 = (sums[nb]-sums[s])/draw2
                  mu = sums[nb]/nb
                  if (draw1*klChange(mu1,mu)+draw2*klChange(mu2,mu))>beta(nb,delta)
                     newEpisode = true
                     check +=1
                  end
               end
               s=s+1
            end
            if (check>0)
               episode+=1
               ChangePoints=append!(ChangePoints,t)
               if verbose=="on"
                  print("detected a change on $(I) at t=$(t)\n\n")
               end
               # reset the history of arm I only
               tloc[I]=0
               SUMS[I]=[]
               TotalNumber[I]=0
               TotalSum[I]=0
            end
         end
      end
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

## Other approaches based on CPD

# M-(kl)UCB adaptation from [Cao et al. 20]
# recommended tuning: w of order log(2KT^2)/(size of the smallest change)^2 [or other domain-specific tuning]
# b = sqrt(w*log(2*K*T^2)/2)
# alpha = sqrt((nbreaks*K*(2*b+3*\sqrt(w)))/(2*T))

function MklUCB(Table,alpha=0.1,w=150,b=50.,verbose="off")
   # exploration parameter alpha (=gamma in the paper), windows w, threshold b
   (K,T)=size(Table)
   # index of the episode
   episode = 1 # nb of breakpoints detected
   ChangePoints = []
   ExploRange = floor(Int,K/alpha)
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   t = 0 # current number of sample
   while (t<T)
      SUMS = [[] for a=1:K] # array containing the cumulated sums of rewards
      tloc = 0 # number of samples in current episode
      TotalNumber=zeros(Int,K) # total number of selections since last restart
      TotalSum=zeros(Float32,K) # sum of rewards since last restart
      newEpisode = false # should we start a new episode?
      # start the episode
      while (!newEpisode)&&(t<T)
         # default exploratory action
         I = (t % ExploRange)
         if (I ==0)||(I>K)
            # compute KL-UCB indices if no exploration
            indices = ones(K)
            for i in 1:K
               if (TotalNumber[i]>0)
                  indices[i]=klUp(TotalSum[i]/TotalNumber[i],log(tloc)/TotalNumber[i])
               end
            end
            I = randmax(indices)
         end
         # get the reward
         rew = Table[I,t+1]
         ChosenArms[t+1] = I
         ReceivedRewards[t+1]=rew
         # update everything
         t=t+1
         tloc=tloc+1
         TotalNumber[I]+=1
         TotalSum[I]+=rew
         append!(SUMS[I],TotalSum[I])
         # has the selected arm changed?
         nb = TotalNumber[I]
         if (nb >= w)
            sums = SUMS[I]
            mid = sums[nb-floor(Int,w/2)]
            FirstHalf = mid - sums[nb-w+1]
            SecondHalf = sums[nb]-mid
            if (abs(FirstHalf-SecondHalf)> b)
               newEpisode = true
               episode+=1
               ChangePoints=append!(ChangePoints,t)
               if verbose=="on"
                  print("detected a change on $(I) at t=$(t)\n\n")
               end
            end
         end
      end
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

# CUSUM-(kl)UCB adaptation from [Liu et al 2018]
# recommended tuning h = log(T/nbreaks) and alpha=sqrt(nbreaks*log(T/nbreaks)/T)
# M and epsilon require some prior knowledge of the problem

function CUSUMklUCB(Table,alpha=0.05,M=150,epsilon=0.1,h=50,verbose="off")
   # exploration parameter alpha, intial nb of sample M, minimum change epsilon, threshold h
   (K,T)=size(Table)
   # index of the episode
   episode = 1 # nb of breakpoints detected
   ChangePoints = []
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   t = 0 # current number of sample
   TotalNumber=zeros(Int,K) # total number of selections since last restart
   TotalSum=zeros(Float32,K) # sum of rewards since last restart
   # random walks associated to each arm used for stopping
   GPlus = zeros(Float32,K) # first random walk associated to each arm (used in CUSUM test)
   GMinus = zeros(Float32,K) # second random walk associated to each arm
   InitialAverage = zeros(Float32,K) # average of the first M observations from each arm
   newEpisode = false # should we start a new episode?
   while (t<T)
      newEpisode = false # should we start a new episode?
      # start the episode
      while (!newEpisode)&&(t<T)
         I = 0
         if (rand()<alpha)
            # explore uniformly at random
            I = rand(1:K)
         else
            # compute KL-UCB indices and choose the argmax if no exploration
            indices = ones(K)
            nvalid = sum(TotalNumber)
            for i in 1:K
               if (TotalNumber[i]>0)
                  indices[i]=klUp(TotalSum[i]/TotalNumber[i],log(nvalid)/TotalNumber[i])
               end
            end
            I = randmax(indices)
         end
         # get the reward
         rew = Table[I,t+1]
         ChosenArms[t+1] = I
         ReceivedRewards[t+1]=rew
         # update everything
         t=t+1
         TotalNumber[I]+=1
         TotalSum[I]+=rew
         if (TotalNumber[I] == M)
            InitialAverage[I] = TotalSum[I]/TotalNumber[I]
         elseif (TotalNumber[I]>M)
            GPlus[I]=max(0,GPlus[I]+rew-InitialAverage[I]-epsilon)
            GMinus[I]=max(0,GMinus[I]+InitialAverage[I]-rew-epsilon)
         end
         # has the selected arm changed?
         if (GMinus[I]>h)||(GPlus[I]>h)
            newEpisode = true
            episode+=1
            ChangePoints=append!(ChangePoints,t)
            if verbose=="on"
               print("detected a change on $(I) at t=$(t)\n\n")
            end
            # reset the history of arm I only
            TotalNumber[I]=0
            TotalSum[I]=0
            InitialAverage[I]=0
            GPlus[I]=0
            GMinus[I]=0
         end
      end
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

##  EXP3.S [Auer et al. 02]
# recommended tuning alpha=1/T and gamma = min(1,sqrt(K*(nbreak*log(K*T)+exp(1))/((exp(1)-1)*T)))
# (automatic tuning if nbreak is given)

function EXP3S(Table,nbreak=0,gamma=0.1,alpha=0.01)
   # table contains a table of rewards for the K arms up to horizon T
   (K,T)=size(Table)
   episode=1
   ChangePoints = []
   if (nbreak>0)
      # optimized gamma as a function of the nb of breakpoints
      gamma = min(1,sqrt(K*(nbreak*log(K*T)+exp(1))/((exp(1)-1)*T)))
      alpha = 1/T
   end
   Weights=(1/K)*ones(Float32,K) # vector of weights / probability to sample
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   for t in 1:T
      Probas = (1-gamma)*Weights .+ gamma/K
      I = sampleDiscrete(Probas)
      # get the reward
      rew = Table[I,t]
      ChosenArms[t]=I
      ReceivedRewards[t]=rew
      # update the weights
      normrew = rew/Probas[I]
      bonus = (exp(1)*alpha/K)
      for k in 1:K
         if k!=I
            Weights[k]=Weights[k]+bonus
         else
            Weights[k]=(exp(gamma*normrew/K))*Weights[k]+bonus
         end
      end
      # normalization step: storing the normalized weights and working with them leads to the same algorithm
      Weights = Normalize(Weights)
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

## Passively adaptive algorithms

# Discounted-(kl)UCB adaptation of [Kocsis et al. 06,Garivier and Moulines 08]
# recommended tuning: 1- 0.25*sqrt(nbBreaks/T)
# (automatic tuning if nbreak is given)

function DklUCB(Table,nbreak=0,gamma=0.95)
   # table contains a table of rewards for the K arms up to horizon T
   (K,T)=size(Table)
   episode=1
   ChangePoints = []
   if (nbreak>0)
      # optimized gamma as a function of the nb of breakpoints
      gamma = 1- 0.25*sqrt(nbBreaks/T)
   end
   DiscNumber=zeros(Float32,K) # discounted number of selections
   DiscSum=zeros(Float32,K) # discounted sum of rewards
   Time =0 # discounted total number of pulls
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   for t in 1:T
      # compute KL-UCB indices
      indices = ones(K)
      for i in 1:K
         if (DiscNumber[i]>0)
            indices[i]=klUp(DiscSum[i]/DiscNumber[i],log(Time)/DiscNumber[i])
         end
      end
      I = randmax(indices)
      # get the reward
      rew = Table[I,t]
      ChosenArms[t]=I
      ReceivedRewards[t]=rew
      # update everything
      DiscNumber = gamma*DiscNumber
      DiscNumber[I]+=1
      DiscSum = gamma*DiscSum
      DiscSum[I]+=rew
      Time = 1+gamma*Time
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

# SW-(kl)UCB adaptation of [Garivier and Moulines 08]
# recommended tuning:  tau = 2*sqrt(T*log(T)/nbreak)
# (automatic tuning if nbreak is given)

function SWklUCB(Table,nbreak=0,tau=450)
   # table contains a table of rewards for the K arms up to horizon T
   (K,T)=size(Table)
   episode=1
   ChangePoints = []
   if (nbreak>0)
      # optimized tau as a function of the nb of breakpoints
      tau = floor(Int,2*sqrt(T*log(T)/nbreak))
   end
   TabNumber=zeros(Int,T,K) # N[t,k] : nb of selections of arm k in first t rounds
   TabSum=zeros(Float32,T,K) # S[t,k] : sum of rewards from arm k in first t rounds
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   for t in 1:T
      # compute KL-UCB indices
      indices = ones(K)
      NumberLoc = TabNumber[max(t-1,1),:]
      SumLoc = TabSum[max(t-1,1),:]
      if (t>tau)
         NumberLoc = TabNumber[t-1,:]-TabNumber[t-tau,:]
         SumLoc = TabSum[t-1,:] - TabSum[t-tau,:]
      end
      for i in 1:K
         if (NumberLoc[i]>0)
            indices[i]=klUp(SumLoc[i]/NumberLoc[i],log(min(t,tau))/NumberLoc[i])
         end
      end
      I = randmax(indices)
      # get the reward
      rew = Table[I,t]
      ChosenArms[t]=I
      ReceivedRewards[t]=rew
      # update everything
      TabNumber[t,:] = TabNumber[max(t-1,1),:]
      TabNumber[t,I]=TabNumber[t,I]+1
      TabSum[t,:] = TabSum[max(t-1,1),:]
      TabSum[t,I]=TabSum[t,I]+rew
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

# D-TS, investigated by [Raj and Kalyani 2017]
# recommended tuning is gamma=0.75, which is not so robust

function DTS(Table,gamma=0.75,a=1,b=1)
   # discount gamma, Beta(a,b) prior
   (K,T)=size(Table)
   episode=1
   ChangePoints = []
   DiscNumber=zeros(Float32,K) # discounted number of selections
   DiscSum=zeros(Float32,K) # discounted sum of rewards
   Time =0 # discounted total number of pulls
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   for t in 1:T
      # compute KL-UCB indices
      indices = ones(K)
      for i in 1:K
         if (DiscNumber[i]>0)
            indices[i]=rand(Beta(a+DiscSum[i], b+DiscNumber[i]-DiscSum[i]), 1)[1]
         end
      end
      I = randmax(indices)
      # get the reward
      rew = Table[I,t]
      ChosenArms[t]=I
      ReceivedRewards[t]=rew
      # update everything
      DiscNumber = gamma*DiscNumber
      DiscNumber[I]+=1
      DiscSum = gamma*DiscSum
      DiscSum[I]+=rew
      Time = 1+gamma*Time
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

## Variants of kl-UCB [Cappé et al 2013]

function klUCB(Table)
   # table contains a table of rewards for the K arms up to horizon T
   (K,T)=size(Table)
   episode=1
   ChangePoints = []
   # index of the episode
   TotalNumber=zeros(Int,K) # total number of selections
   TotalSum=zeros(Float32,K) # sum of rewards
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   for t in 1:T
      # compute KL-UCB indices
      indices = ones(K)
      for i in 1:K
         if (TotalNumber[i]>0)
            indices[i]=klUp(TotalSum[i]/TotalNumber[i],log(t)/TotalNumber[i])
         end
      end
      I = randmax(indices)
      # get the reward
      rew = Table[I,t]
      ChosenArms[t]=I
      ReceivedRewards[t]=rew
      # update everything
      TotalNumber[I]+=1
      TotalSum[I]+=rew
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end

function OracleklUCB(Table,Breaks=[])
   # requires a vector Breaks giving the positions of the breakpoints
   (K,T)=size(Table)
   ChosenArms = zeros(Int,T)
   ReceivedRewards = zeros(T)
   # episodes and change-points are predetermined
   episode = length(Breaks)+1
   ChangePoints = copy(Breaks)
   # important times
   Times = copy(Breaks)
   append!(Times,T)
   pushfirst!(Times,1)
   for k in 1:episode
      TotalNumber=zeros(Int,K) # total number of selections during episode
      TotalSum=zeros(Float32,K) # sum of rewards
      for t in Times[k]:Times[k+1]
         # compute KL-UCB indices
         indices = ones(K)
         tloc = t -  Times[k] + 1
         for i in 1:K
            if (TotalNumber[i]>0)
               indices[i]=klUp(TotalSum[i]/TotalNumber[i],log(tloc)/TotalNumber[i])
            end
         end
         I = randmax(indices)
         # get the reward
         rew = Table[I,t]
         ChosenArms[t]=I
         ReceivedRewards[t]=rew
         # update everything
         TotalNumber[I]+=1
         TotalSum[I]+=rew
      end
   end
   return ChosenArms,ReceivedRewards,episode,ChangePoints
end


text:
using Clang

# get the headers
VK_DIR = joinpath(@__DIR__, "Vulkan-Headers")
if !isdir(VK_DIR)
    run(`git clone https://github.com/KhronosGroup/Vulkan-Headers.git $VK_DIR`)
else
    run(`git -C $VK_DIR fetch`)
end

run(`git -C $VK_DIR checkout v1.1.114`)

# generate Vulkan bindings
const VK_INCLUDE = joinpath(@__DIR__, "Vulkan-Headers", "include")
const VK_HEADERS = map(x->joinpath(VK_INCLUDE, "vulkan", x), ["vk_platform.h", "vulkan.h", "vulkan_core.h"])
VK_EXTENSIONS = [
    "VK_USE_PLATFORM_ANDROID_KHR",
    "VK_USE_PLATFORM_FUCHSIA",
    "VK_USE_PLATFORM_IOS_MVK",
    "VK_USE_PLATFORM_MACOS_MVK",
    "VK_USE_PLATFORM_METAL_EXT",
    "VK_USE_PLATFORM_VI_NN",
    "VK_USE_PLATFORM_WAYLAND_KHR",
    "VK_USE_PLATFORM_WIN32_KHR",
    "VK_USE_PLATFORM_XCB_KHR",
    "VK_USE_PLATFORM_XLIB_KHR",
    "VK_USE_PLATFORM_XLIB_XRANDR_EXT",
    "VK_USE_PLATFORM_GGP",
    ]

clang_extraargs = String[]
for extension in VK_EXTENSIONS
    push!(clang_extraargs, "-D")
    push!(clang_extraargs, extension)
end

wc = init(; headers = VK_HEADERS,
            output_file = joinpath(@__DIR__, "vk_api.jl"),
            common_file = joinpath(@__DIR__, "vk_common.jl"),
            clang_includes = vcat(VK_INCLUDE, CLANG_INCLUDE),
            clang_args = clang_extraargs,
            header_wrapped = (root, current)->root == current,
            header_library = x->"libvulkan",
            clang_diagnostics = true,
            )
run(wc)

ext: jl


text:
# function for generating data from the DWP-SDE model with the V_{extended} potential function

doc"""
    generate_data(theta, theta_known)

Generates data from the model Z_t = X_t + U_t where
dX_t = -V_extended'(X_t)dt + dWt and U_t follows an UO-process.

# Inputs
* `theta`:true values for the unknown parameters.
* `theta_known`: known parameters.

# Outputs
* `Z_sim`: simualted Z process.
* `dt`: dt for X process.
* `diff_dt`: differeance between dt for X and U process.
"""
function generate_data(theta, theta_known, scale_grid)

  #=
  if length(theta) == 2
    (c,d) = exp(theta)
    (Κ, Γ, A, B, f, g, power1, power2, sigma) = theta_known
  elseif length(theta) == 4
    (A, c, d, g) = exp(theta)
    (Κ, Γ, B, f, power1, power2, sigma) = theta_known
  elseif length(theta) == 3
    (A, c, d) = exp(theta)
    (Κ, Γ, B, f, g, power1, power2, sigma) = theta_known
  elseif length(theta) == 6
    (A,c,d,power1,power2,sigma) = exp(theta)
    (Κ, Γ, B, f, g) = theta_known
  else
    # estiamte all i.e 7 parameters kappa, gamma, A,c,d,g and sigma
    (B,f) = theta_known

    (Κ, Γ, A,B,f,g,power1,power2,sigma) = exp(theta)
  end
  =#

  # set parameters, use this to move between difference parameter combinations easily
  (Κ, Γ, A,A_sign,B,c,d,g,f,power1,power2,sigma) = set_parameters(theta, theta_known, length(theta))

  # grid for the X_t process
  length_of_data = 25000
  #scale_grid = 4 #1  # choose lenght of data
  T = 1*scale_grid
  N = convert(Int64, scale_grid*length_of_data)
  dt_hat = T/N
  dt = dt_hat*(1/4)*10^3*3.5
  diff_dt = 1 #convert(Int64, round(10^-2. / dt))
  x_grid = 0:dt_hat:T

  # grid for the U_t process
  N_U = convert(Int64, scale_grid*length_of_data)
  dt_U = 1
  #dt_U = dt*diff_dt
  u_grid = 0:T/N_U:T

  Z_0 = 24.5 # start in the lower state

  X_sim = zeros(N+1)
  U_sim = zeros(N_U+1)
  dB = zeros(N+1)

  X_sim[1] = Z_0
  U_sim[1] = 0

  for i = 1:length(dB)
    dB[i] = sqrt(dt)*randn()
  end

  j = 1 # index for U_sim
  b_const = sqrt(2*sigma^2 / 2)

  #  should also consider the cache memory!
  for i = 1:N # numerical integration of the X process
    X_sim[i+1] = X_sim[i] - (A*X_sim[i] - f + (power2*abs(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])^(power2 - 1)*sign(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])*(g - (B*power1*abs(c - B*X_sim[i])^(power1 - 1)*sign(c - B*X_sim[i]))/2))/2)*dt + b_const*dB[i]
    if  mod(i, diff_dt) == 0
      U_sim[j+1] = rand(Normal( U_sim[j]*exp(-Κ*dt_U), sqrt( Γ^2*( 1 - exp( -2*Κ*dt_U ) ) ) ),1)[1]
      j = j + 1
    end
  end

  X_thinned = X_sim[1:diff_dt:end]

  Z_sim = X_thinned + U_sim

  return Z_sim, dt, diff_dt, X_thinned

end

# function for generating data from the DWP-SDE model with the V_{extended} potential function

doc"""
    generate_data(theta, theta_known, scale_grid)
Generates data from the model Z_t = X_t + U_t where
dX_t = -V_extended'(X_t)dt + dWt and U_t follows an UO-process.
# Inputs
* `theta`:true values for the unknown parameters.
* `theta_known`: known parameters.
* `scale_grid`: scaling of data set.
# Outputs
* `Z_sim`: simualted Z process.
* `dt`: dt for X process.
* `diff_dt`: differeance between dt for X and U process.
"""
function generate_data(theta,theta_known, scale_grid, dt, dt_U, nbr_sim_steps, start_val)

  # set parameters, use this to move between difference parameter combinations easily
  (Κ, Γ, A,A_sign,B,c,d,g,f,power1,power2,sigma) = set_parameters(theta, theta_known, length(theta))

  # grid for the X_t process
  #length_of_data = 25000
  length_of_data = nbr_sim_steps
  #scale_grid = 4 #1  # choose lenght of data
  T = 1*scale_grid
  N = convert(Int64, scale_grid*length_of_data)

  # grid for X process
  N_dt = 1 #Int(round(dt_U/dt))

  #Z_0 = 24.5 # start in the lower state

  Z_0 =  start_val
  X_sim = zeros(N+1)
  U_sim = zeros(N+1)
  dB = zeros(N_dt)

  X_sim[1] = Z_0
  U_sim[1] = 0

  b_const = sqrt(2*sigma^2 / 2)

  #  should also consider the cache memory!
  for i = 1:N # numerical integration of the X process

    # generare normal dist. noise
    for j = 1:length(N_dt)
      dB[j] = sqrt(dt)*randn()
    end

    # propagate the X process
    for j = 1:N_dt
      X_sim[i+1] = X_sim[i] - (A*X_sim[i] - f + (power2*abs(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])^(power2 - 1)*sign(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])*(g - (B*power1*abs(c - B*X_sim[i])^(power1 - 1)*sign(c - B*X_sim[i]))/2))/2)*dt + b_const*dB[j]
    end

    U_sim[i+1] = rand(Normal( U_sim[i]*exp(-Κ*dt_U), sqrt( Γ^2*( 1 - exp( -2*Κ*dt_U ) ) ) ),1)[1]

  end

  Z_sim = X_sim + U_sim

  return Z_sim, N_dt

end

doc"""
    generate_data(theta, theta_known)

Generates data from the model Z_t = X_t + U_t where
dX_t = -V_extended'(X_t)dt + dWt and U_t follows an UO-process.

# Inputs
* `theta`:true values for the unknown parameters.
* `theta_known`: known parameters.

# Outputs
* `Z_sim`: simualted Z process.
* `dt`: dt for X process.
* `diff_dt`: differeance between dt for X and U process.
"""
function generate_data(theta, theta_known)

  #=
  if length(theta) == 2
    (c,d) = exp(theta)
    (Κ, Γ, A, B, f, g, power1, power2, sigma) = theta_known
  elseif length(theta) == 4
    (A, c, d, g) = exp(theta)
    (Κ, Γ, B, f, power1, power2, sigma) = theta_known
  elseif length(theta) == 3
    (A, c, d) = exp(theta)
    (Κ, Γ, B, f, g, power1, power2, sigma) = theta_known
  elseif length(theta) == 6
    (A,c,d,power1,power2,sigma) = exp(theta)
    (Κ, Γ, B, f, g) = theta_known
  else
    # estiamte all i.e 7 parameters kappa, gamma, A,c,d,g and sigma
    (B,f) = theta_known

    (Κ, Γ, A,B,f,g,power1,power2,sigma) = exp(theta)
  end
  =#

  # set parameters, use this to move between difference parameter combinations easily
  (Κ, Γ, A,B,c,d,g,f,power1,power2,sigma) = set_parameters(theta, theta_known, length(theta))

  thinning_sim = 100

  N = Int64(3.5e6/thinning_sim)
  dt = 0.005*thinning_sim#20*1e-9
  dt_U = 1#20*1e-9
  diff_dt = 1

  X_sim = zeros(N+1)
  U_sim = zeros(N+1)

  X_sim[1] = 35
  U_sim[1] = 0

  dB = zeros(N+1)

  for i = 1:length(dB)
    dB[i] = sqrt(dt)*randn()
  end

  j = 1 # index for U_sim
  b_const = sqrt(2*sigma^2 / 2)

  #  should also consider the cache memory!
  for i = 1:N # numerical integration of the X process
    X_sim[i+1] = X_sim[i] - (A*X_sim[i] - f + (power2*abs(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])^(power2 - 1)*sign(abs(c - B*X_sim[i])^power1/2 - d + g*X_sim[i])*(g - (B*power1*abs(c - B*X_sim[i])^(power1 - 1)*sign(c - B*X_sim[i]))/2))/2)*dt + b_const*dB[i]
    if  mod(i, diff_dt) == 0
      U_sim[j+1] = rand(Normal( U_sim[j]*exp(-Κ*dt_U), sqrt( Γ^2*( 1 - exp( -2*Κ*dt_U ) ) ) ),1)[1]
      j = j + 1
    end
  end

  X_thinned = X_sim[1:diff_dt:end]
  Z_sim = X_thinned + U_sim

  return Z_sim, X_thinned

end

max_stars_repo_path: DWPSDE model/generate_data.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/generate_data.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/generate_data.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
# Definition of the composite type "Resolution"

using Printf

"""
    struct Resolution

`Resolution` objects are needed to perform a number of pixel-related
functions, e.g., convert a direction into a pixel number and vice
versa.

The fields of a `Resolution` object are the following:

- `nside`: the NSIDE parameter
- `nsideTimesTwo`: 2 * NSIDE
- `nsideTimesFour`: 4 * NSIDE
- `numOfPixels`: number of pixels in the map
- `order`: order of the map
- `pixelsPerFace`: number of pixels in each Healpix face
- `ncap`
- `fact2`
- `fact1`
"""
struct Resolution
    nside::Int
    nsideTimesTwo::Int
    nsideTimesFour::Int
    numOfPixels::Int

    order::Int
    pixelsPerFace::Int
    ncap::Int
    fact2::Float64
    fact1::Float64
end

# Customize printing
Base.show(io::IO, r::Resolution) = @printf(io, "Healpix resolution(NSIDE = %d)", r.nside)

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

"""
    Resolution(nside) -> Resolution

Create a `Resolution` object, given a value for `NSIDE`.
"""
function Resolution(nside)
    (1 ≤ nside ≤ NSIDE_MAX) || throw(DomainError(nside))
    # The expression (nside & (nside - 1)) == 0 is a quick check for
    # detecting if nside is a power of two or not
    (nside & (nside - 1) == 0) || throw(DomainError())

    order          = ilog2(nside)
    pixelsPerFace  = nside^2
    numOfPixels    = 12pixelsPerFace
    ncap           = 2 * nside * (nside - 1)
    fact2          = 4 / numOfPixels
    fact1          = 2 * nside * fact2

    result = Resolution(nside,
                        2nside,
                        4nside,
                        numOfPixels,
                        order,
                        pixelsPerFace,
                        ncap,
                        fact2,
                        fact1)

end

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

"""
    nsideok(nside::Integer) -> Bool

Check whether `nside` is a valid `NSIDE` parameter.
"""
nsideok(nside::Integer) = (nside > 0) && ((nside) & (nside - 1) == 0)

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

"""
    nside2npix(nside::Integer) -> Integer

Return the number of pixels for a Healpix map with the specified
`NSIDE` value. If `NSIDE` is not an integer power of two, the function
throws a `DomainError` exception.
"""
function nside2npix(nside::Integer)
    (nside > 0) || throw(DomainError(nside, "`NSIDE` is not a positive number"))

    nsidelog2 = round(Int, log2(nside))
    (2^nsidelog2 == nside) || throw(DomainError(nside, "`NSIDE` is not an integer power of two"))

    12(nside^2)
end

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

"""
    npix2nside(npix::Integer) -> Integer

Given the number of pixels in a Healpix map, return the `NSIDE`
resolution parameter. If the number is invalid, throw a `DomainError`
exception.
"""
function npix2nside(npix::Integer)
    (npix % 12 == 0) || throw(DomainError(npix, "Invalid number of pixels"))

    square_root = sqrt(npix / 12)
    (square_root^2 == npix / 12) || throw(DomainError(npix, "Invalid number of pixels"))

    convert(Int, round(square_root))
end

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

"""
    nside2pixarea(nside::Integer) -> Real

Return the solid angle of a pixel in a map with the specified `NSIDE` parameter.
The result is expressed in steradians.
"""
nside2pixarea(nside::Integer) = 4π / nside2npix(nside)

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

"""
    nside2resol(nside::Integer) -> Real

Return the approximate resolution of a map with the specified `NSIDE`. The
resolution is expressed in radians, and it is the square root of the pixel
size.
"""
nside2resol(nside::Integer) = sqrt(nside2pixarea(nside))

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


text:
module TestHistograms

include("../../../examples/histograms.jl")
using Test

function test_one_to_ten()
    indices = CuVector(1:10)
    h = countints_svector_functional(indices, Val(10))
    @test collect(h) == fill(1, 10)
end

function test_rand()
    indices = (floor(Int, x * 10) + 1 for x in CUDA.rand(2^30))
    h = countints_svector_functional(indices, Val(10))
    h = collect(h)
    @test h ./ h[1] ≈ fill(1, 10) atol = 0.01
end

end  # module


text:
#########################
# Main types for meshes #
#########################
"""
    Node{dim, T}

A `Node` is a point in space.

# Fields
- `x::Vec{dim,T}`: stores the coordinates
"""
struct Node{dim,T}
    x::Vec{dim,T}
end
Node(x::NTuple{dim,T}) where {dim,T} = Node(Vec{dim,T}(x))
getcoordinates(n::Node) = n.x

abstract type AbstractCell{dim,N,M} end
"""
    Cell{dim,N,M} <: AbstractCell{dim,N,M}
A `Cell` is a sub-domain defined by a collection of `Node`s as it's vertices.
However, a `cell` is not defined by the nodes but rather by the global node ids.
The parameter `dim` refers here to the geometrical/ambient dimension, i.e. the dimension of the `nodes` in the grid and **not** the topological dimension of the cell.
A `Cell` has `N` nodes and `M` faces.

# Fields
- `nodes::Ntuple{N,Int}`: N-tuple that stores the node ids. The ordering defines a cell's and its subentities' orientations.
"""
struct Cell{dim,N,M} <: AbstractCell{dim,N,M}
    nodes::NTuple{N,Int}
end
nfaces(c::C) where {C<:AbstractCell} = nfaces(typeof(c))
nfaces(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = M
nnodes(c::C) where {C<:AbstractCell} = nnodes(typeof(c))
nnodes(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = N

# Typealias for commonly used cells
const Line  = Cell{1,2,2}
const Line2D = Cell{2,2,1}
const Line3D = Cell{3,2,0}
const QuadraticLine = Cell{1,3,2}

const Triangle = Cell{2,3,3}
const QuadraticTriangle = Cell{2,6,3}

const Quadrilateral = Cell{2,4,4}
const Quadrilateral3D = Cell{3,4,1}
const QuadraticQuadrilateral = Cell{2,9,4}

const Tetrahedron = Cell{3,4,4}
const QuadraticTetrahedron = Cell{3,10,4}

const Hexahedron = Cell{3,8,6}

"""
A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh
"""
struct CellIndex
    idx::Int
end

"""
A `FaceIndex` wraps an (Int, Int) and defines a local face by pointing to a (cell, face).
"""
struct FaceIndex <: BoundaryIndex
    idx::Tuple{Int,Int} # cell and side
end

"""
A `EdgeIndex` wraps an (Int, Int) and defines a local edge by pointing to a (cell, edge).
"""
struct EdgeIndex <: BoundaryIndex
    idx::Tuple{Int,Int} # cell and side
end

"""
A `VertexIndex` wraps an (Int, Int) and defines a local vertex by pointing to a (cell, vert).
"""
struct VertexIndex <: BoundaryIndex
    idx::Tuple{Int,Int} # cell and side
end

abstract type AbstractGrid{dim} end

"""
    Grid{dim, C<:AbstractCell, T<:Real} <: AbstractGrid}

A `Grid` is a collection of `Cells` and `Node`s which covers the computational domain, together with Sets of cells, nodes and faces.
There are multiple helper structures to apply boundary conditions or define subdomains. They are gathered in the `cellsets`, `nodesets`,
`facesets`, `edgesets` and `vertexsets`.

# Fields
- `cells::Vector{C}`: stores all cells of the grid
- `nodes::Vector{Node{dim,T}}`: stores the `dim` dimensional nodes of the grid
- `cellsets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of cell ids
- `nodesets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of global node ids
- `facesets::Dict{String,Set{FaceIndex}}`: maps a `String` to a `Set` of `Set{FaceIndex} (global_cell_id, local_face_id)`
- `edgesets::Dict{String,Set{EdgeIndex}}`: maps a `String` to a `Set` of `Set{EdgeIndex} (global_cell_id, local_edge_id`
- `vertexsets::Dict{String,Set{VertexIndex}}`: maps a `String` key to a `Set` of local vertex ids
- `boundary_matrix::SparseMatrixCSC{Bool,Int}`: optional, only needed by `onboundary` to check if a cell is on the boundary, see, e.g. Helmholtz example
"""
mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
    cells::Vector{C}
    nodes::Vector{Node{dim,T}}
    # Sets
    cellsets::Dict{String,Set{Int}}
    nodesets::Dict{String,Set{Int}}
    facesets::Dict{String,Set{FaceIndex}}
    edgesets::Dict{String,Set{EdgeIndex}}
    vertexsets::Dict{String,Set{VertexIndex}}
    # Boundary matrix (faces per cell × cell)
    boundary_matrix::SparseMatrixCSC{Bool,Int}
end

function Grid(cells::Vector{C},
              nodes::Vector{Node{dim,T}};
              cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              facesets::Dict{String,Set{FaceIndex}}=Dict{String,Set{FaceIndex}}(),
              edgesets::Dict{String,Set{EdgeIndex}}=Dict{String,Set{EdgeIndex}}(),
              vertexsets::Dict{String,Set{VertexIndex}}=Dict{String,Set{VertexIndex}}(),
              boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
    return Grid(cells, nodes, cellsets, nodesets, facesets, edgesets, vertexsets, boundary_matrix)
end

##########################
# Grid utility functions #
##########################
@inline getdim(grid::AbstractGrid{dim}) where {dim} = dim
"""
    getcells(grid::AbstractGrid)
    getcells(grid::AbstractGrid, v::Union{Int,Vector{Int}}
    getcells(grid::AbstractGrid, setname::String)

Returns either all `cells::Collection{C<:AbstractCell}` of a `<:AbstractGrid` or a subset based on an `Int`, `Vector{Int}` or `String`.
Whereas the last option tries to call a `cellset` of the `grid`. `Collection` can be any indexable type, for `Grid` it is `Vector{C<:AbstractCell}`.
"""
@inline getcells(grid::AbstractGrid) = grid.cells
@inline getcells(grid::AbstractGrid, v::Union{Int, Vector{Int}}) = grid.cells[v]
@inline getcells(grid::AbstractGrid, setname::String) = grid.cells[collect(getcellset(grid,setname))]
"Returns the number of cells in the `<:AbstractGrid`."
@inline getncells(grid::AbstractGrid) = length(grid.cells)
"Returns the celltype of the `<:AbstractGrid`."
@inline getcelltype(grid::AbstractGrid) = eltype(grid.cells)
@inline getcelltype(grid::AbstractGrid, i::Int) = typeof(grid.cells[i])

"""
    getnodes(grid::AbstractGrid)
    getnodes(grid::AbstractGrid, v::Union{Int,Vector{Int}}
    getnodes(grid::AbstractGrid, setname::String)

Returns either all `nodes::Collection{N}` of a `<:AbstractGrid` or a subset based on an `Int`, `Vector{Int}` or `String`.
The last option tries to call a `nodeset` of the `<:AbstractGrid`. `Collection{N}` refers to some indexable collection where each element corresponds
to a Node.
"""
@inline getnodes(grid::AbstractGrid) = grid.nodes
@inline getnodes(grid::AbstractGrid, v::Union{Int, Vector{Int}}) = grid.nodes[v]
@inline getnodes(grid::AbstractGrid, setname::String) = grid.nodes[collect(getnodeset(grid,setname))]
"Returns the number of nodes in the grid."
@inline getnnodes(grid::AbstractGrid) = length(grid.nodes)
"Returns the number of nodes of the `i`-th cell."
@inline nnodes_per_cell(grid::AbstractGrid, i::Int=1) = nnodes(grid.cells[i])

"""
    getcellset(grid::AbstractGrid, setname::String)

Returns all cells as cellid in a `Set` of a given `setname`.
"""
@inline getcellset(grid::AbstractGrid, setname::String) = grid.cellsets[setname]
"""
    getcellsets(grid::AbstractGrid)

Returns all cellsets of the `grid`.
"""
@inline getcellsets(grid::AbstractGrid) = grid.cellsets

"""
    getnodeset(grid::AbstractGrid, setname::String)

Returns all nodes as nodeid in a `Set` of a given `setname`.
"""
@inline getnodeset(grid::AbstractGrid, setname::String) = grid.nodesets[setname]
"""
    getnodesets(grid::AbstractGrid)

Returns all nodesets of the `grid`.
"""
@inline getnodesets(grid::AbstractGrid) = grid.nodesets

"""
    getfaceset(grid::AbstractGrid, setname::String)

Returns all faces as `FaceIndex` in a `Set` of a given `setname`.
"""
@inline getfaceset(grid::AbstractGrid, setname::String) = grid.facesets[setname]
"""
    getfacesets(grid::AbstractGrid)

Returns all facesets of the `grid`.
"""
@inline getfacesets(grid::AbstractGrid) = grid.facesets

"""
    getedgeset(grid::AbstractGrid, setname::String)

Returns all edges as `EdgeIndex` in a `Set` of a given `setname`.
"""
@inline getedgeset(grid::AbstractGrid, setname::String) = grid.edgesets[setname]
"""
    getedgesets(grid::AbstractGrid)

Returns all edge sets of the grid.
"""
@inline getedgesets(grid::AbstractGrid) = grid.edgesets

"""
    getedgeset(grid::AbstractGrid, setname::String)

Returns all vertices as `VertexIndex` in a `Set` of a given `setname`.
"""
@inline getvertexset(grid::AbstractGrid, setname::String) = grid.vertexsets[setname]
"""
    getvertexsets(grid::AbstractGrid)

Returns all vertex sets of the grid.
"""
@inline getvertexsets(grid::AbstractGrid) = grid.vertexsets

n_faces_per_cell(grid::Grid) = nfaces(eltype(grid.cells))

"""
    function compute_vertex_values(grid::AbstractGrid, f::Function)
    function compute_vertex_values(grid::AbstractGrid, v::Vector{Int}, f::Function)
    function compute_vertex_values(grid::AbstractGrid, set::String, f::Function)

Given a `grid` and some function `f`, `compute_vertex_values` computes all nodal values,
 i.e. values at the nodes,  of the function `f`.
The function implements two dispatches, where only a subset of the grid's node is used.

```julia
    compute_vertex_values(grid, x -> sin(x[1]) + cos([2]))
    compute_vertex_values(grid, [9, 6, 3], x -> sin(x[1]) + cos([2])) #compute function values at nodes with id 9,6,3
    compute_vertex_values(grid, "right", x -> sin(x[1]) + cos([2])) #compute function values at nodes belonging to nodeset right
```

"""
@inline function compute_vertex_values(nodes::Vector{Node{dim,T}}, f::Function) where{dim,T}
    map(n -> f(getcoordinates(n)), nodes)
end

@inline function compute_vertex_values(grid::AbstractGrid, f::Function)
    compute_vertex_values(getnodes(grid), f::Function)
end

@inline function compute_vertex_values(grid::AbstractGrid, v::Vector{Int}, f::Function)
    compute_vertex_values(getnodes(grid, v), f::Function)
end

@inline function compute_vertex_values(grid::AbstractGrid, set::String, f::Function)
    compute_vertex_values(getnodes(grid, set), f::Function)
end

# Transformations
"""
    transform!(grid::Abstractgrid, f::Function)

Transform all nodes of the `grid` based on some transformation function `f`.
"""
function transform!(g::AbstractGrid, f::Function)
    c = similar(g.nodes)
    for i in 1:length(c)
        c[i] = Node(f(g.nodes[i].x))
    end
    copyto!(g.nodes, c)
    g
end

# Sets

_check_setname(dict, name) = haskey(dict, name) && throw(ArgumentError("there already exists a set with the name: $name"))
_warn_emptyset(set, name) = length(set) == 0 && @warn("no entities added to the set with name: $name")

"""
    addcellset!(grid::AbstractGrid, name::String, cellid::Union{Set{Int}, Vector{Int}})
    addcellset!(grid::AbstractGrid, name::String, f::function; all::Bool=true)

Adds a cellset to the grid with key `name`.
Cellsets are typically used to define subdomains of the problem, e.g. two materials in the computational domain.
The `MixedDofHandler` can construct different fields which live not on the whole domain, but rather on a cellset.

```julia
addcellset!(grid, "left", Set((1,3))) #add cells with id 1 and 3 to cellset left
addcellset!(grid, "right", x -> norm(x[1]) < 2.0 ) #add cell to cellset right, if x[1] of each cell's node is smaller than 2.0
```
"""
function addcellset!(grid::AbstractGrid, name::String, cellid::Union{Set{Int},Vector{Int}})
    _check_setname(grid.cellsets,  name)
    cells = Set(cellid)
    _warn_emptyset(cells, name)
    grid.cellsets[name] = cells
    grid
end

function addcellset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true)
    _check_setname(grid.cellsets, name)
    cells = Set{Int}()
    for (i, cell) in enumerate(getcells(grid))
        pass = all
        for node_idx in cell.nodes
            node = grid.nodes[node_idx]
            v = f(node.x)
            all ? (!v && (pass = false; break)) : (v && (pass = true; break))
        end
        pass && push!(cells, i)
    end
    _warn_emptyset(cells, name)
    grid.cellsets[name] = cells
    grid
end

"""
    addfaceset!(grid::AbstractGrid, name::String, faceid::Union{Set{FaceIndex},Vector{FaceIndex}})
    addfaceset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true)

Adds a faceset to the grid with key `name`.
A faceset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id, local_face_id)`.
Facesets are used to initialize `Dirichlet` structs, that are needed to specify the boundary for the `ConstraintHandler`.

```julia
addfaceset!(gird, "right", Set(((2,2),(4,2))) #see grid manual example for reference
addfaceset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0) #see incompressible elasticity example for reference
```
"""
addfaceset!(grid::Grid, name::String, set::Union{Set{FaceIndex},Vector{FaceIndex}}) =
    _addset!(grid, name, set, grid.facesets)
addedgeset!(grid::Grid, name::String, set::Union{Set{EdgeIndex},Vector{EdgeIndex}}) =
    _addset!(grid, name, set, grid.edgesets)
addvertexset!(grid::Grid, name::String, set::Union{Set{VertexIndex},Vector{VertexIndex}}) =
    _addset!(grid, name, set, grid.vertexsets)
function _addset!(grid::AbstractGrid, name::String, _set, dict::Dict)
    _check_setname(dict, name)
    set = Set(_set)
    _warn_emptyset(set, name)
    dict[name] = set
    grid
end

addfaceset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) =
    _addset!(grid, name, f, Ferrite.faces, grid.facesets, FaceIndex; all=all)
addedgeset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) =
    _addset!(grid, name, f, Ferrite.edges, grid.edgesets, EdgeIndex; all=all)
addvertexset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) =
    _addset!(grid, name, f, Ferrite.vertices, grid.vertexsets, VertexIndex; all=all)
function _addset!(grid::AbstractGrid, name::String, f::Function, _ftype::Function, dict::Dict, _indextype::Type; all::Bool=true)
    _check_setname(dict, name)
    _set = Set{_indextype}()
    for (cell_idx, cell) in enumerate(getcells(grid))
        for (face_idx, face) in enumerate(_ftype(cell))
            pass = all
            for node_idx in face
                v = f(grid.nodes[node_idx].x)
                all ? (!v && (pass = false; break)) : (v && (pass = true; break))
            end
            pass && push!(_set, _indextype(cell_idx, face_idx))
        end
    end
    _warn_emptyset(_set, name)
    dict[name] = _set
    grid
end

"""
    addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int},Set{Int}})
    addnodeset!(grid::AbstractGrid, name::String, f::Function)

Adds a `nodeset::Dict{String, Set{Int}}` to the `grid` with key `name`. Has the same interface as `addcellset`.
However, instead of mapping a cell id to the `String` key, a set of node ids is returned.
"""
function addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int},Set{Int}})
    _check_setname(grid.nodesets, name)
    grid.nodesets[name] = Set(nodeid)
    _warn_emptyset(grid.nodesets[name], name)
    grid
end

function addnodeset!(grid::AbstractGrid, name::String, f::Function)
    _check_setname(grid.nodesets, name)
    nodes = Set{Int}()
    for (i, n) in enumerate(getnodes(grid))
        f(n.x) && push!(nodes, i)
    end
    grid.nodesets[name] = nodes
    _warn_emptyset(grid.nodesets[name], name)
    grid
end

"""
    getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::Int)
Fills the vector `x` with the coordinates of a cell, defined by its cell id.
"""
@inline function getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::Int) where {dim,T}
    #@assert length(x) == N
    @inbounds for i in 1:length(x)
        x[i] = grid.nodes[grid.cells[cell].nodes[i]].x
    end
end
@inline getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::CellIndex) where {dim, T} = getcoordinates!(x, grid, cell.idx)
@inline getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, face::FaceIndex) where {dim, T} = getcoordinates!(x, grid, face.idx[1])

"""
    getcoordinates(grid::AbstractGrid, cell)
Return a vector with the coordinates of the vertices of cell number `cell`.
"""
@inline function getcoordinates(grid::AbstractGrid, cell::Int)
    # TODO pretty ugly, worth it?
    dim = typeof(grid.cells[cell]).parameters[1]
    T = typeof(grid).parameters[3]
    nodeidx = grid.cells[cell].nodes
    return [grid.nodes[i].x for i in nodeidx]::Vector{Vec{dim,T}}
end
@inline getcoordinates(grid::AbstractGrid, cell::CellIndex) = getcoordinates(grid, cell.idx)
@inline getcoordinates(grid::AbstractGrid, face::FaceIndex) = getcoordinates(grid, face.idx[1])

# Iterate over cell vector
function Base.iterate(c::Vector{Cell{dim,N}}, state = 1) where {dim, N}
    if state > length(c)
        return nothing
    else
        return (CellIndex(state), state + 1)
    end
end

function Base.show(io::IO, ::MIME"text/plain", grid::Grid)
    print(io, "$(typeof(grid)) with $(getncells(grid)) ")
    typestrs = sort!(collect(Set(celltypes[typeof(x)] for x in grid.cells)))
    str = join(io, typestrs, '/')
    print(io, " cells and $(getnnodes(grid)) nodes")
end

const celltypes = Dict{DataType, String}(Cell{1,2,2}  => "Line",
                                         Cell{2,2,2}  => "2D-Line",
                                         Cell{3,2,0}  => "3D-Line",
                                         Cell{1,3,2}  => "QuadraticLine",
                                         Cell{2,3,3}  => "Triangle",
                                         Cell{2,6,3}  => "QuadraticTriangle",
                                         Cell{2,4,4}  => "Quadrilateral",
                                         Cell{3,4,1}  => "3D-Quadrilateral",
                                         Cell{2,9,4}  => "QuadraticQuadrilateral",
                                         Cell{3,4,4}  => "Tetrahedron",
                                         Cell{3,10,4} => "QuadraticTetrahedron",
                                         Cell{3,8,6}  => "Hexahedron",
                                         Cell{3,20,6} => "Cell{3,20,6}")

# Functions to uniquely identify vertices, edges and faces, used when distributing
# dofs over a mesh. For this we can ignore the nodes on edged, faces and inside cells,
# we only need to use the nodes that are vertices.
# 1D: vertices
faces(c::Union{Line,QuadraticLine}) = (c.nodes[1], c.nodes[2])
vertices(c::Union{Line,Line2D,Line3D,QuadraticLine}) = (c.nodes[1], c.nodes[2])
# 2D: vertices, faces
faces(c::Line2D) = ((c.nodes[1],c.nodes[2]),)
vertices(c::Union{Triangle,QuadraticTriangle}) = (c.nodes[1], c.nodes[2], c.nodes[3])
faces(c::Union{Triangle,QuadraticTriangle}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]))
vertices(c::Union{Quadrilateral,Quadrilateral3D,QuadraticQuadrilateral}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
faces(c::Union{Quadrilateral,QuadraticQuadrilateral}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))
# 3D: vertices, edges, faces
edges(c::Line3D) = ((c.nodes[1],c.nodes[2]),)
vertices(c::Union{Tetrahedron,QuadraticTetrahedron}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
edges(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]), (c.nodes[1],c.nodes[4]), (c.nodes[2],c.nodes[4]), (c.nodes[3],c.nodes[4]))
faces(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2],c.nodes[3]), (c.nodes[1],c.nodes[2],c.nodes[4]), (c.nodes[2],c.nodes[3],c.nodes[4]), (c.nodes[1],c.nodes[4],c.nodes[3]))
vertices(c::Union{Hexahedron,Cell{3,20,6}}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6], c.nodes[7], c.nodes[8])
edges(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]), (c.nodes[5],c.nodes[6]), (c.nodes[6],c.nodes[7]), (c.nodes[7],c.nodes[8]), (c.nodes[8],c.nodes[5]), (c.nodes[1],c.nodes[5]), (c.nodes[2],c.nodes[6]), (c.nodes[3],c.nodes[7]), (c.nodes[4],c.nodes[8]))
faces(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[4],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[6],c.nodes[5]), (c.nodes[2],c.nodes[3],c.nodes[7],c.nodes[6]), (c.nodes[3],c.nodes[4],c.nodes[8],c.nodes[7]), (c.nodes[1],c.nodes[5],c.nodes[8],c.nodes[4]), (c.nodes[5],c.nodes[6],c.nodes[7],c.nodes[8]))
edges(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))
faces(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2],c.nodes[3],c.nodes[4]),)

# random stuff
default_interpolation(::Union{Type{Line},Type{Line2D},Type{Line3D}}) = Lagrange{1,RefCube,1}()
default_interpolation(::Type{QuadraticLine}) = Lagrange{1,RefCube,2}()
default_interpolation(::Type{Triangle}) = Lagrange{2,RefTetrahedron,1}()
default_interpolation(::Type{QuadraticTriangle}) = Lagrange{2,RefTetrahedron,2}()
default_interpolation(::Union{Type{Quadrilateral},Type{Quadrilateral3D}}) = Lagrange{2,RefCube,1}()
default_interpolation(::Type{QuadraticQuadrilateral}) = Lagrange{2,RefCube,2}()
default_interpolation(::Type{Tetrahedron}) = Lagrange{3,RefTetrahedron,1}()
default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{3,RefTetrahedron,2}()
default_interpolation(::Type{Hexahedron}) = Lagrange{3,RefCube,1}()
default_interpolation(::Type{Cell{3,20,6}}) = Serendipity{3,RefCube,2}()

boundaryfunction(::Type{FaceIndex}) = Ferrite.faces
boundaryfunction(::Type{EdgeIndex}) = Ferrite.edges
boundaryfunction(::Type{VertexIndex}) = Ferrite.vertices

for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex)
    @eval begin
        #Constructor
        ($INDEX)(a::Int, b::Int) = ($INDEX)((a,b))

        Base.getindex(I::($INDEX), i::Int) = I.idx[i]

        #To be able to do a,b = faceidx
        Base.iterate(I::($INDEX), state::Int=1) = (state==3) ?  nothing : (I[state], state+1)

        #For (cellid, faceidx) in faceset
        Base.in(v::Tuple{Int, Int}, s::Set{$INDEX}) = in($INDEX(v), s)
    end
end


text:
using InteractionNetworkModels, StructuredDistances, Plots, Distributions

# Testing sampler
# ---------------

E = [[1,2,1,2],
    [1,2,1],
    [3,4,3],
    [3,4],
    [1,2],
    [1,2,1],
    [1,2,3],
    [4,5],
    [7,8]]
d = MatchingDistance(FastLCS(101))
K_inner, K_outer = (DimensionRange(1,10), DimensionRange(1,25))
model = SIM(
    E, 3.5,
    d,
    1:10,
    K_inner, K_outer
)

mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=1, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.8, 1, model.K_inner.u),
    burn_in=2000, lag=75, init=InitRandIns(10)
)

@time test = mcmc_sampler(
    model, desired_samples=5000, burn_in=0, lag=1
)
plot(test)
summaryplot(test)

@time mcmc_out = mcmc_sampler(
    model,
    desired_samples=50,
    lag=500,
    burn_in=10000
)
plot(mcmc_out)
summaryplot(mcmc_out)

data = mcmc_out.sample
E_prior = SIM(E, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
γ_prior = Uniform(0.5,7.0)

posterior = SimPosterior(data, E_prior, γ_prior)

# Construct posterior sampler
posterior_sampler = SimIexSplitMerge(
    mcmc_sampler,
    len_dist=TrGeometric(0.9,K_inner.l,K_inner.u),
    ν_ed=1, ν_td=1,
    β=0.7, ε=0.3
)

# Mode Conditional

E_init = [[1,2,1,3,4,3],[7,8,4,5,1,2,1],[10,1]]
@time posterior_out = posterior_sampler(
    posterior,
    # 4.9,
    desired_samples=100, lag=1, burn_in=0,
    S_init=E_init, γ_init=2.8,
    aux_init_at_prev=true,
);
plot(
    posterior_out, E, size=(600,400),
    xlabel=["" "Sample Index"],
    ylabel=["Distance to True Mode" "γ"],
    yguidefontsize=[9 9],
    xguidefontsize=9
)
posterior_out.S_sample

posterior_out
# Testing iMCMC proposal functions
# --------------------------------

aux_mcmc = SimMcmcInsertDelete()
mcmc = SimIexSplitMerge(aux_mcmc)

# Standard proposal generation functions

# Split
S_curr = [rand(1:3, 3) for i in 1:5]
S_prop = deepcopy(S_curr)

imcmc_multi_split_prop_sample!(S_curr, S_prop, mcmc)

S_curr
S_prop

# Merge
S_curr = [rand(1:3, 3) for i in 1:5]
S_prop = deepcopy(S_curr)

imcmc_multi_merge_prop_sample!(S_curr, S_prop, mcmc)

S_curr
S_prop

# When we have length one paths
S_curr = [rand(1:3, 1) for i in 1:5]
push!(S_curr, [1,2])
push!(S_curr,[1,2,3,2])
S_prop = deepcopy(S_curr)

imcmc_special_multi_split_prop_sample!(S_curr, S_prop, mcmc)

S_curr
S_prop

S_curr = deepcopy(S_prop)

imcmc_special_multi_merge_prop_sample!(S_curr, S_prop, mcmc)

S_curr
S_prop


text:
abstract type TDFunctional{T} end
Base.eltype(x::TDFunctional{T}) where {T} = T

function quaddata(exc::TDFunctional, testrefs, timerefs, testels, timeels)

    testqd = quadpoints(testrefs, testels, (2,))
    timeqd = quadpoints(timerefs, timeels, (10,))

    testqd, timeqd

end

function quaddata(excitation::TDFunctional,
        test_refspace, time_refspace::DiracBoundary,
        test_elements, time_elements)

    test_quad_data = quadpoints(test_refspace, test_elements, (2,))

    test_quad_data, nothing
end

function quadrule(exc::TDFunctional, testrefs, timerefs, p, τ, r, ρ, qd)

    MultiQuadStrategy(
        qd[1][1,p],
        SingleQuadStrategy2(
            qd[2][1,r]
        )
    )

end

function quadrule(exc::TDFunctional, testrefs, timerefs::DiracBoundary, p, τ, r, ρ, qd)

    MultiQuadStrategy(
        qd[1][1,p],
        nothing
    )

end

function assemble(exc::TDFunctional, testST)
    testfns = spatialbasis(testST)
    timefns = temporalbasis(testST)
    Z = zeros(eltype(exc), numfunctions(testfns), numfunctions(timefns))
    store(v,m,k) = (Z[m,k] += v)
    assemble!(exc, testST, store)
    return Z
end

function assemble(exc::TDFunctional, testST::StagedTimeStep)
    stageCount = length(testST.c);
    spatialBasis = testST.spatialBasis;
    Nt = testST.Nt;
    Δt = testST.Δt;
    Z = zeros(eltype(exc), numfunctions(spatialBasis) * stageCount, Nt);
    for i = 1:stageCount
        store(v,m,k) = (Z[(m-1)*stageCount+i,k] += v);
        tbsd = TimeBasisDeltaShifted(timebasisdelta(Δt, Nt), testST.c[i]);
        assemble!(exc, spatialBasis ⊗ tbsd, store);
    end
    return Z;
end

function assemble!(exc::TDFunctional, testST, store)

    testfns = spatialbasis(testST)
    timefns = temporalbasis(testST)

    testrefs = refspace(testfns)
    timerefs = refspace(timefns)

    testels, testad = assemblydata(testfns)
    timeels, timead = assemblydata(timefns)

    qd = quaddata(exc, testrefs, timerefs, testels, timeels)

    z = zeros(eltype(exc), numfunctions(testrefs), numfunctions(timerefs))
    for p in eachindex(testels)
        τ = testels[p]
        for r in eachindex(timeels)
            ρ = timeels[r]

            fill!(z, 0)
            qr = quadrule(exc, testrefs, timerefs, p, τ, r, ρ, qd)
            momintegrals!(z, exc, testrefs, timerefs, τ, ρ, qr)

            for i in 1 : numfunctions(testrefs)
                for d in 1 : numfunctions(timerefs)

                    v = z[i,d]
                    for (m,a) in testad[p,i]
                        for (k,c) in timead[r,d]
                            store(a*c*v, m, k)
                        end
                    end

                end
            end
        end
    end
end

abstract type NumQuadStrategy end

mutable struct MultiQuadStrategy{P,R} <: NumQuadStrategy
    quad_points::P
    inner_rule::R
end

# TODO: consolidate with the existing definition of SingleQuadStrategy
mutable struct SingleQuadStrategy2{P} <: NumQuadStrategy
    quad_points::P
end

timequadrule(qr::MultiQuadStrategy, p) = qr.inner_rule

function momintegrals!(z, exc::TDFunctional, testrefs, timerefs, τ, ρ, qr)

    for p in qr.quad_points
        x = p.point
        w = p.weight
        f = p.value
        dx = w

        # try
        #     @assert ρ.vertices[1][1] <= cartesian(x)[1] <= ρ.vertices[2][1]
        # catch
        #     @show ρ.vertices[1][1]
        #     @show cartesian(x)[1]
        #     @show ρ.vertices[2][1]
        #     error("")
        # end

        tqr = timequadrule(qr,p)
        timeintegrals!(z, exc, testrefs, timerefs, x, ρ, dx, tqr, f)

    end

end

function timeintegrals!(z, exc::TDFunctional, testrefs, timerefs, testpoint, timeelement, dx, qr, f)

    for p in qr.quad_points
        t = p.point
        w = p.weight
        U = p.value
        dt = w #* jacobian(t) # * volume(timeelement)

        for i in 1 : numfunctions(testrefs)
            for k in 1 : numfunctions(timerefs)
                z[i,k] += dot(f[i][1]*U[k], exc(testpoint,t)) * dt * dx
            end
        end
    end
end

function timeintegrals!(z, exc::TDFunctional,
        spacerefs, timerefs::DiracBoundary,
        testpoint, timeelement,
        dx, qr, testvals)

        # since timeelement uses barycentric coordinates,
        # the first/left vertex has coords u = 1.0!
        testtime = neighborhood(timeelement, point(0.0))
        @assert cartesian(testtime)[1] ≈ timeelement.vertices[2][1]

        for i in 1 : numfunctions(spacerefs)
            z[i,1] += dot(testvals[i][1], exc(testpoint, testtime)) * dx
        end
end


text:
using Requests
using URIParser
using Gadfly

## geturl(line) = strip(split(line, ",")[4], '\"')

## function getallurls(filename)
##     lines = split(readall(filename), "\n")[2:end-1]
##     urls = Array(String, length(lines))
##     for i=1:length(lines)
##         urls[i] = "http://" * geturl(lines[i])
##     end
##     return urls
## end

## function resolvehost(url)
##     h = URI(url).host
##     ip = string(Base.getaddrinfo(h))
##     replace(url, h, ip)
## end

## function resolvehosts(urls)
##     result = String[]
##     @sync for url in urls
##         @async try
##             resolved = resolvehost(strip(url))
##             push!(result, resolved)
##         catch
##             println("Can't resolve $url")
##         end
##     end
##     return result
## end

## function get_stream(url::String)
##     Requests.open_stream(URI(url), Dict(), "", "GET")
## end

# @time urls = resolvehosts(open(readlines, expanduser("../domains.txt"))[1:5000])

function fetch_urls(urls; timeout=Inf)
    statuses = -ones(length(urls))
    times = -ones(length(urls))
    @sync for i=1:length(urls)
        @async begin
            try
                t = @elapsed resp = get(urls[i]; timeout=timeout)
                statuses[i] = resp.status
                times[i] = t
                print("OK\n")
            catch
                print("Error\n")
            end
        end
    end
    return (statuses, times)
end

function main()
    @time urls = [strip(l) for l in
                  open(readlines, expanduser("~/Dropbox/ips.txt"))][1:100]
    @time (statuses, times) = fetch_urls(urls, timeout=2)
    plot(x=[i for i=1:length(times)], y=times)
end


text:
module UtilsTests

using Test

@testset "Helpers" begin include("HelpersTests.jl") end

@testset "CachedArrays" begin include("CachedArraysTests.jl") end

@testset "CachedSubVectors" begin include("CachedSubVectorsTests.jl") end

@testset "CachedValues" begin include("CachedValuesTests.jl") end

end # module


text:
module SimpleLockFiles

import Random: randstring

include("core.jl")

export SimpleLockFile
export lock_path, rand_lkid
export acquire_lock, force_unlock, unlock
export lock

end


text:
using DataFrames, XLSX
wb = "C://Users//onmyw//Downloads//Craig's Design and Landscaping Services Dec 10, 2021//General_ledger.xlsx"
ws = "General Ledger"
rng = "A5:J452"
x = DataFrame(XLSX.readtable(wb,ws)...)
names(x)

function qb(;wb=missing,ws=missing,rng=missing,teststate=false)
    x = DataFrame(XLSX.readtable(wb,ws)..., header=false)

XLSX.writetable(wb, x, overwrite=true, sheetname=ws,anchor_cell=outcell())
names


text:
include("lib_for_robot.jl")
function task12(r::Robot,n::Int)
    x, y = go_to_bottom_left_corner_return_coords_shift_of_start(r)

    all = show_path_after_moving(r, Nord)+1
    while isborder(r, Sud) == false
        move!(r, Sud)
    end

    side=Nord
    t=0
    while isborder(r, Ost) == false
        if n <= t < 2n
            for i in 1:n
                move!(r,side)
            end
        end

        direction_mark(r,side,n,all)

        t = t + 1
        if t > 2n
            t -= 2n
        end

        move!(r, Ost)
        side = HorizonSide(mod(Int(side) + 2, 4))
        show_path_after_moving(r, side)
        side = HorizonSide(mod(Int(side) + 2, 4))
    end
    if n <= t < 2n
        for i in 1:n
            move!(r, side)
        end
    end
    direction_mark(r, side, n, all)

    go_to_bottom_left_corner_return_coords_shift_of_start(r)
    move_n(r, Ost, x)
    move_n(r, Nord, y)
end


text:
"""
Unit tests for methods defined in utils.jl.

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

# Standard library
using Test

# External packages
using ZMQ: Socket

# Boid.jl package
using Boid

# --- Test fixtures

# --- Constructor tests

@testset "utils.jl: construct_ipc_url()" begin
    # --- ipc_dir does not start with "ipc://"

    # ipc_dir is a relative path
    ipc_dir = "path/to/ipc/dir"
    ipc_name = "ipc-name"
    ipc_url = construct_ipc_url(ipc_dir, ipc_name)
    @test ipc_url == "ipc://path/to/ipc/dir/ipc-name"

    # ipc_dir is an absolute path
    ipc_dir = "/path/to/ipc/dir"
    ipc_name = "ipc-name"
    ipc_url = construct_ipc_url(ipc_dir, ipc_name)
    @test ipc_url == "ipc:///path/to/ipc/dir/ipc-name"

    # --- ipc_dir starts with "ipc://"

    # ipc_dir is a relative path
    ipc_dir = "ipc://path/to/ipc/dir"
    ipc_name = "ipc-name"
    ipc_url = construct_ipc_url(ipc_dir, ipc_name)
    @test ipc_url == "ipc://path/to/ipc/dir/ipc-name"

    # ipc_dir is an absolute path
    ipc_dir = "ipc:///path/to/ipc/dir"
    ipc_name = "ipc-name"
    ipc_url = construct_ipc_url(ipc_dir, ipc_name)
    @test ipc_url == "ipc:///path/to/ipc/dir/ipc-name"

    # --- error cases

    ipc_dir = "invalid://ipc:///path/to/ipc/dir"
    ipc_name = "ipc-name"
    expected_message = """Invalid `ipc_dir` "$ipc_dir"."""
    @test_throws(ArgumentError(expected_message),
                 construct_ipc_url(ipc_dir, ipc_name))
end


text:
import JuMP: DiffableFunction

function frank_wolfe(p::CombiProblem{DiffableFunction})
  # check only constraint is one combinatorial constraint
  @assert length(p.convex_constraints==0) && length(p.combi_constraints==1)
  combiset = p.combi_constraints[1]
  n = combiset.dim
  x = zeros(n)
  for k=1:100 # maxiters
    g = grad(p.objective, x)
    h = fenchel(combiset, g)
    alpha = 1/k # stepsize
    x = (1-alpha)*x + alpha*h
  end
  return x
end


text:
export fetchval

const SHOULD_SEND = Dict(:PATCH => true, :POST => true, :PUT => true)

"""
A wrapper around a response from the REST API. Every function which wraps a Discord REST
API endpoint returns a `Future` which will contain a value of this type. To retrieve the
`Response` from the `Future`, use `fetch`. See also: [`fetchval`](@ref).

## Fields
- `val::Nullable{T}`: The object contained in the HTTP response. For example, for a
  call to [`get_channel_message`](@ref), this value will be a [`Message`](@ref).
- `ok::Bool`: The state of the request. If `true`, then it is safe to access `val`.
- `http_response::Nullable{HTTP.Messages.Response}`: The underlying HTTP response, if
  a request was made.
- `exception::Nullable{Exception}`: The caught exception, if one is thrown.

## Examples
Multiple API calls which immediately return `Future`s and can be awaited:
```julia
futures = map(i -> create_message(c, channel_id; content=string(i)), 1:10);
other_work_here()
resps = fetch.(futures)
```
Skipping error checks and returning the value directly:
```julia
guild = fetchval(create_guild(c; name="foo"))
```
"""
struct Response{T}
    val::Nullable{T}
    ok::Bool
    http_response::Nullable{HTTP.Messages.Response}
    exception::Nullable{Exception}
end

# HTTP response with no body.
function Response{Nothing}(r::HTTP.Messages.Response)
    return Response{Nothing}(nothing, r.status < 300, r, nothing)
end

# HTTP response with body (maybe).
function Response{T}(c::Client, r::HTTP.Messages.Response) where T
    if T == Nothing
        return Response{Nothing}(r)
    end
    r.status == 204 && return Response{T}(nothing, true, r, nothing)
    r.status >= 300 && return Response{T}(nothing, false, r, nothing)

    val = if HTTP.header(r, "Content-Type") == "application/json"
        symbolize(JSON3.read(String(copy(r.body)), T))
    else
        copy(r.body)
    end

    return Response{T}(val, true, r, nothing)
end

# HTTP request with no expected response body.
function Response(
    c::Client,
    method::Symbol,
    endpoint::AbstractString;
    body="",
    kwargs...
)
    return Response{Nothing}(c, method, endpoint; body=body, kwargs...)
end

# HTTP request.
function Response{T}(
    c::Client,
    method::Symbol,
    endpoint::AbstractString,
    cachetest::Function=alwaystrue;
    headers=Dict(),
    body=HTTP.nobody,
    kwargs...
) where T
    f = Future()
    @spawn begin
        kws = (
            channel=cap("channels", endpoint),
            guild=cap("guilds", endpoint),
            message=cap("messages", endpoint),
            user=cap("users", endpoint),
        )

        # Retrieve the value from the cache, maybe.
        @debug "Preparing Request: Looking in cache" Time=now()
        if c.use_cache && method === :GET
            val = get(c.state, T; kws...)
            if val !== nothing
                try
                    # Only use the value if it satisfies the cache test.
                    if cachetest(val) === true
                        put!(f, Response{T}(val, true, nothing, nothing))
                        return
                    end
                catch
                    # TODO: Log this?
                end
            end
        end
        # Prepare the request.
        url = "$DISCORD_API/v$(c.version)$endpoint"
        @debug "Preparing Request: Did not find in cache, preparing request" Time=now()
        isempty(kwargs) || (url *= "?" * HTTP.escapeuri(kwargs))
        headers = merge(Dict(
            "User-Agent" => "Discord.jl $D1ZKORD_JL_VERSION",
            "Content-Type" => "application/json",
            "Authorization" => c.token,
        ), Dict(headers))
        args = [method, url, headers]
        @debug "Preparing Request: Type of request: $T" Time=now()
        if get(SHOULD_SEND, method, false)
            push!(args, headers["Content-Type"] == "application/json" ? JSON3.write(body) : body)
        end

        @debug "Preparing Request: Enqueuing request" Time=now()
        # Queue a job to be run within the rate-limiting constraints.
        enqueue!(c.limiter, method, endpoint) do
            @debug "Enqueued job running"
            http_r = nothing

            try
                # Make an HTTP request, and generate a Response.
                # If we hit an upstream rate limit, return the response immediately.
                http_r = HTTP.request(args...; status_exception=false)
                @debug "Sent http request" Time=now()
                if http_r.status - 200 >= 100
                    @warn "Got an unexpected response" To=args[2] Code=http_r.status Body=http_r Sent=args[4]
                end
                http_r.status == 429 && return http_r
                r = Response{T}(c, http_r)

                if r.ok && r.val !== nothing
                    # Cache the value, and also replace any missing data
                    # in the response with data from the cache.
                    r = @set r.val = put!(c.state, r.val; kws...)
                end

                # Store the successful Response to the Future.
                put!(f, r)
                @debug "Got response"
            catch e
                kws = logkws(c; conn=undef, exception=(e, catch_backtrace()))
                @error "Creating response failed" kws...
                put!(f, Response{T}(nothing, false, http_r, e))
            end

            @debug "Request completed" url=args[2]
            return http_r
        end
    end
    @debug "Returning future" Time=now()
    return f
end

Base.eltype(::Response{T}) where T = T
Base.eltype(::Type{Response{T}}) where T = T

"""
    fetchval(f::Future{Response{T}}) -> Nullable{T}

Shortcut for `fetch(f).val`: Fetch a [`Response`](@ref) and return its value. Note that
there are no guarantees about the response's success and the value being returned, and it
discards context that can be useful for debugging, such as HTTP responses and caught
exceptions.
"""
fetchval(f::Future) = fetch(f).val

# Capture an ID from a string.
function cap(path::AbstractString, s::AbstractString)
    m = match(Regex("/$path/(\\d+)"), s)
    return m === nothing ? nothing : parse(Snowflake, first(m.captures))
end

include(joinpath("endpoints", "endpoints.jl"))
include(joinpath("crud", "crud.jl"))


text:
module TestShow

include("preamble.jl")

if VERSION >= v"1.6-"
    ==ᵣ = ==ₛ
else
    ==ᵣ = ==
end

xforms = [
    Cat(),
    Count(),
    Count(2),
    Count(3, 4),
    GetIndex([0]),
    # GetIndex{true}([0]),
    Iterated(sqrt, 0.1),
    # Iterated(sqrt, 0.1, Number),
    Map(sin),
    Filter(isfinite),
    Scan(*),
    ZipSource(Filter(isfinite)),
    ZipSource(opcompose(Filter(isfinite), Map(sin))),
    opcompose(ZipSource(opcompose(Filter(isfinite), Map(sin))), Map(sum)),
    let xf = ZipSource(opcompose(Filter(isfinite), Map(sin)))
        opcompose(xf, Map(sum), xf)
    end,
    let xf = opcompose(Map(first), Map(last))
        xf = opcompose(ZipSource(opcompose(ZipSource(xf), Map(identity))), xf)
        opcompose(xf, ZipSource(xf), xf)
    end,
    # opcompose(Zip(Map(sin), Map(cos), Map(tan)), Map(prod)),
]

@testset "smoke test summary(xf)" begin
    for xf in xforms
        @test sprint(summary, xf) isa String
    end
end

# Note: Use `ENV["JULIA_DEBUG"] = "Main"` to enable debugging
@testset "eval(show(xf))" begin
    @testset "$(summary(xf))" for xf in xforms
        code = sprint(show, xf)
        @debug """
        show(xf) =
        $code
        """
        VERSION < v"1.5-beta" && occursin("⨟", code) && continue
        xf2 = include_string(@__MODULE__, code)
        @test xf == xf2
    end
end
@testset "eval(show(text/plain, xf))" begin
    @testset "$(summary(xf))" for xf in xforms
        code = sprint(show, "text/plain", xf)
        @debug """
        show("text/plain", xf) =
        $code
        """
        VERSION < v"1.5-beta" && occursin("⨟", code) && continue
        xf2 = include_string(@__MODULE__, code)
        @test xf == xf2
    end
end

@testset "Reduced" begin
    @test sprint(show, Reduced(1)) == "Transducers.Reduced(1)"
    @test sprint(show, Reduced(1); context = :limit => true) == "Reduced(1)"
end

@testset "Completing(+)" begin
    rf = Completing(+)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.Completing(+)"
    @test sprint(print, rf; kw...) == "Transducers.Completing(+)"
    @test sprint(show, rf; kw...) == "Transducers.Completing(+)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.Completing(+)"
end

@testset "Completing{Any}(+)" begin
    rf = Completing{Any}(+)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.Completing{Any}(+)"
    @test sprint(print, rf; kw...) == "Transducers.Completing{Any}(+)"
    @test sprint(show, rf; kw...) == "Transducers.Completing{Any}(+)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.Completing{Any}(+)"
end

@testset "TeeRF(min, max)" begin
    rf = TeeRF(min, max)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.TeeRF(min, max)"
    @test sprint(print, rf; kw...) == "Transducers.TeeRF(min, max)"
    @test sprint(show, rf; kw...) == "Transducers.TeeRF(min, max)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.TeeRF(min, max)"
end

@testset "TeeRF{2,Tuple{Any,Any}}((min, max))" begin
    rf = TeeRF{2,Tuple{Any,Any}}((min, max))
    kw = (; context = :module => Base)
    @test repr(rf; kw...) ==ᵣ "Transducers.TeeRF{2,Tuple{Any,Any}}((min, max))"
    @test sprint(print, rf; kw...) ==ᵣ "Transducers.TeeRF{2,Tuple{Any,Any}}((min, max))"
    @test sprint(show, rf; kw...) ==ᵣ "Transducers.TeeRF{2,Tuple{Any,Any}}((min, max))"
    @test sprint(show, "text/plain", rf; kw...) ==ᵣ
          "Transducers.TeeRF{2,Tuple{Any,Any}}((min, max))"
end

@testset "ProductRF(min, max)" begin
    rf = ProductRF(min, max)
    kw = (; context = :module => Base)
    @test repr(rf; kw...) == "Transducers.ProductRF(min, max)"
    @test sprint(print, rf; kw...) == "Transducers.ProductRF(min, max)"
    @test sprint(show, rf; kw...) == "Transducers.ProductRF(min, max)"
    @test sprint(show, "text/plain", rf; kw...) == "Transducers.ProductRF(min, max)"
end

@testset "ProductRF{2,Tuple{Any,Any}}((min, max))" begin
    rf = ProductRF{2,Tuple{Any,Any}}((min, max))
    kw = (; context = :module => Base)
    @test repr(rf; kw...) ==ᵣ "Transducers.ProductRF{2,Tuple{Any,Any}}((min, max))"
    @test sprint(print, rf; kw...) ==ᵣ "Transducers.ProductRF{2,Tuple{Any,Any}}((min, max))"
    @test sprint(show, rf; kw...) ==ᵣ "Transducers.ProductRF{2,Tuple{Any,Any}}((min, max))"
    @test sprint(show, "text/plain", rf; kw...) ==ᵣ
          "Transducers.ProductRF{2,Tuple{Any,Any}}((min, max))"
end

end  # module


text:
# Python version created in 2012.
# Ported to julia in January, 2018.
#
# Author: Rich Plevin
#
# Copyright (c) 2012-2018. The Regents of the University of California (Regents)
# and Richard Plevin. See the file COPYRIGHT.txt for details.
#
# Implements the Latin Hypercube Sampling technique as described by Iman and Conover, 1982,
# including correlation control both for no correlation or for a specified rank correlation
# matrix for the sampled parameters. Original python version was heavily modified from
# http://nullege.com/codes/show/src@m@o@model-builder-HEAD@Bayes@lhs.py
#
import StatsBase
using Statistics
using LinearAlgebra

"""
    rank_corr_coef(m::Matrix{Float64})

Take a 2-D array of values and produce a array of rank correlation
coefficients representing the rank correlations pairs of columns.
"""
function rank_corr_coef(m::Matrix{Float64})
    cols = size(m, 2)
    corrCoef = Matrix(1.0I, cols, cols)    # identity matrix

    for i in 1:cols
        for j in (i + 1):cols
            corr = corspearman(m[:, i], m[:, j])
            corrCoef[i, j] = corrCoef[j, i] = corr
        end
    end

    return corrCoef
end

"""
    _gen_rank_values(params::Int, trials::Int, corrmatrix::Matrix{Float64})

Generate a data set of 'trials' ranks for 'params'
parameters that obey the given correlation matrix.

params: number of parameters.
trials: number of trials.
corrmatrix: rank correlation matrix for parameters.
corrmatrix[i,j] denotes the rank correlation between parameters
i and j.

Output is a Matrix with 'trials' rows and 'params' columns.
The i'th column represents the ranks for the i'th parameter.
"""
function _gen_rank_values(params::Int, trials::Int, corrmatrix::Matrix{Float64})
    # Create van der Waarden scores
    strata = collect(1.0:trials) / (trials + 1)
    vdwScores = quantile.(Normal(), strata)

    S = zeros(trials, params)
    for i in 1:params
        shuffle!(vdwScores)
        S[:, i] = vdwScores
    end

    P = Matrix(cholesky(corrmatrix, Val(false)).L)
    E = rank_corr_coef(S)
    Q = Matrix(cholesky(E, Val(false)).L)
    final = (S * inv(Q)') * P'

    ranks = zeros(Int, trials, params)
    for i in 1:params
        ranks[:, i] = ordinalrank(final[:, i])
    end

    return ranks
end

"""
    _get_percentiles(trials::Int)

Generate a list of 'trials' values, one from each of 'trials' equal-size
segments from a uniform distribution. These are used with an RV's ppf
(percent point function = inverse cumulative function) to retrieve the
values for that RV at the corresponding percentiles.
"""
function _get_percentiles(trials::Int)
    segmentSize = 1.0 / trials
    points = rand(Uniform(), trials) * segmentSize + collect(0:trials-1) * segmentSize
    return points
end

"""
    lhs(rvlist::Vector{RandomVariable}, trials::Int; corrmatrix::Union{Matrix{Float64},Nothing}=nothing, asDataFrame::Bool=true)

Produce an array or DataFrame of 'trials' rows of values for the given parameter
list, respecting the correlation matrix 'corrmatrix' if one is specified, using Latin
Hypercube (stratified) sampling.

The values in the i'th column are drawn from the ppf function of the i'th parameter
from rvlist, and each columns i and j are rank correlated according to corrmatrix[i,j].

rvlist: (list of rv-like objects representing parameters) Only requirement
       on parameter objects is that they must implement the ppf function.

trials: (int) number of trials to generate for each parameter.

corrmatrix: a numpy matrix representing the correlation between the parameters.
       corrmatrix[i,j] should give the correlation between the i'th and j'th
       entries of rvlist.

columns: (None or list(str)) Column names to use to return a DataFrame.

skip: (list of params)) Parameters to process later because they are
       dependent on other parameter values (e.g., they're "linked"). These
       cannot be correlated.

Returns DataFrame with `trials` rows of values for the `rvlist`.
"""
function lhs(rvlist::Vector{RandomVariable}, trials::Int;
             corrmatrix::Union{Matrix{Float64},Nothing}=nothing,
             asDataFrame::Bool=true)

    num_rvs = length(rvlist)
    ranks = corrmatrix === nothing ? nothing : _gen_rank_values(num_rvs, trials, corrmatrix)

    samples = zeros(trials, num_rvs)

    for (i, rv) in enumerate(rvlist)
        values = quantile.(Ref(rv.dist), _get_percentiles(trials))  # extract values from the RV for these percentiles

        if corrmatrix === nothing
            shuffle!(values)           # randomize the stratified samples
        else
            indices = ranks[:, i]
            values = values[indices]   # reorder to respect correlations
        end

        samples[:, i] = values
    end

    return asDataFrame ? DataFrame(samples, map(rv->rv.name, rvlist)) : samples
end

function lhs!(mcs::MonteCarloSimulation; corrmatrix::Union{Matrix{Float64},Nothing}=nothing)
    # TBD: verify that any correlated values are actual distributions, not stored vectors?

    trials = mcs.trials
    rvdict = mcs.rvdict
    num_rvs = length(rvdict)
    rvlist = mcs.dist_rvs

    samples = lhs(rvlist, mcs.trials, corrmatrix=corrmatrix, asDataFrame=false)

    for (i, rv) in enumerate(rvlist)
        dist = rv.dist
        name = rv.name
        values = samples[:, i]
        rvdict[name] = RandomVariable(name, SampleStore(values))
    end
    return nothing
end

# TBD: handle this outside lhs since it applies to any linked variables (or drop support for this...)
"""
    lhs_amend!(df::DataFrame, rvlist::Vector{RandomVariable}, trials::Int)

Amend the DataFrame with LHS data by adding columns for the given parameters.
This allows "linked" parameters to refer to the values of other parameters.

df: Generated by prior call to LHS or something similar.
rvlist: The random variables to add to the df
trials: the number of trials to generate for each parameter
"""
function lhs_amend!(df::DataFrame, rvlist::Vector{RandomVariable}, trials::Int, sampling::SamplingOptions=LHS)
    for rv in rvlist
        if sampling == LHS
            values = quantile.(rv.dist, _get_percentiles(trials))  # extract values from the RV for these percentiles
            shuffle!(values)                                       # randomize the stratified samples

        else # sampling == Random
            values = rv.dist.rand(trials)
        end

        df[rv.name] = values
    end
    return nothing
end

"""
    correlation_matrix(mcs::MonteCarloSimulation)

Return a Matrix holding the correlations between random variables
as indicated in the MonteCarloSimulation, or nothing if no correlations
have been defined.

TBD: if needed, compute correlation matrix only for correlated
     RVs, leaving all uncorrelated RVs alone.
"""
function correlation_matrix(mcs::MonteCarloSimulation)
    if length(mcs.corrlist) == 0
        return nothing
    end

    rvdict = mcs.rvdict

    # create a mapping of names to RV position in list
    names = Dict([(rv.name, i) for (i, rv) in enumerate(values(rvdict))])

    count = length(rvdict)
    corrmatrix = Matrix(1.0I, count, count)

    for corr in mcs.corrlist
        n1 = corr.name1
        n2 = corr.name2

        # We don't support correlation between stored samples
        # if rvdict[n1].dist isa SampleStore || rvdict[n2].dist isa SampleStore
        #     error("Correlations with SampleStores is not supported ($n1, $n2)")
        # end

        i = names[n1]
        j = names[n2]

        corrmatrix[i, j] = corrmatrix[j, i] = corr.value
    end

    return corrmatrix
end


text:
using Distributed
@show nworkers() myid()

@everywhere function hoge()
    sleep(1)
    return myid()
end

@everywhere function hoge3d(x, y, z)
    sleep(1)
    return x*y*z
end

function main()
    N = 10
    a = zeros(N)
    vs = [@spawn hoge() for i in 1:N]
    vv = [fetch(s) for s in vs]
    println(vv)
    # fetchを呼ぶとwaitもする
    vs = [@spawn hoge3d(x, y, z) for x in 1:N for y in 1:N for z in 1:N]
    vv = [fetch(s) for s in vs]
    println(sum(vv))
end

@time main()

max_stars_repo_name: Soyukke/julia-paralell-examples

max_issues_repo_name: Soyukke/julia-paralell-examples

max_forks_repo_name: Soyukke/julia-paralell-examples


text:
using GPUArrays, Base.Test
using GPUArrays: JLArray
using GPUArrays.TestSuite

@testset "Julia reference implementation:" begin
    TestSuite.run_tests(JLArray)
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module VisCaptions

export vis_runtype_caps, vis_sim_caps, vis_runinfo_caps

using RLESUtils, Obj2Dict
using SISLES.GenerativeModel

using ...DefineSave
using ...SaveHelpers

#runtype captions
function vis_runtype_caps(d::TrajLog)
    typ = get_run_type(d)
    if typ == "ONCE"
        cap = "Encounter. "
    elseif typ == "MCBest"
        cap = "MC. N=$(get_mcbest_n(d)). "
    elseif typ == "MCTS"
        cap = "MCTS. N=$(get_mcts_iterations(d)). "
    else
        warn("vis_captions::vis_runtype_caps: No such run type! ")
        cap = ""
    end
    cap
end

#sim parameter captions.  TODO: make this more robust
function vis_sim_caps(d::TrajLog)
    #ACASX_GM...for now
    return ""
end

#runinfo captions
function vis_runinfo_caps(d::TrajLog)
  r = round(get_reward(d), 2)
  nmac = get_nmac(d)
  vmd = round(get_vmd(d), 2)
  hmd = round(get_hmd(d), 2)
  mdt = get_md_time(d)
  return "R=$r. vmd=$vmd. hmd=$hmd. md-time=$mdt. NMAC=$nmac. "
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using Test
using ITensors.LazyApply: LazyApply, Add, Mul, ∑, ∏, α, materialize

@testset "LazyApply general functionality" begin
  @test materialize(∏([1, 2, Add(3, 4)])) == prod([1, 2, 3 + 4])
  @test ∏([1, 2, Add(3, 4)]) isa ∏
  @test materialize(3 * ∏([1, 2, Add(3, 4)])) == 3 * prod([1, 2, 3 + 4])
  @test materialize(exp(∏([1, 2, ∑([3, 4])]))) == exp(prod([1, 2, sum([3 + 4])]))
  @test materialize(2 * ∑([1, 2, ∏([3, 4])])) == 2 * sum([1, 2, prod([3, 4])])
  @test 2 * ∑([1, 2, ∏([3, 4])]) == ∑([2, 4, 2∏([3, 4])])
  @test 2 * ∑([1, 2, ∏([3, 4])]) isa ∑
  @test 2∑(["X", "Y"]) == ∑([Mul(2, "X"), Mul(2, "Y")])
  @test materialize(∑() + 3 + 4) == sum([3, 4])
  @test ∑() + 3 + 4 isa ∑
  @test materialize(∑([1, 2, 3]) + ∑([4, 5, 6])) == sum([1, 2, 3, 4, 5, 6])
  @test ∑([1, 2, 3]) + ∑([4, 5, 6]) isa ∑
  @test materialize(Add(1, 2) + Add(3, 4)) == 1 + 2 + 3 + 4
  @test Add(1, 2) + Add(3, 4) == Add(1, 2, 3, 4)
  @test Add(1, 2) + Add(3, 4) isa Add
  @test materialize(2 * Add(1, 2)) == 2 * (1 + 2)
  @test 2 * Add(1, 2) isa Add
  @test materialize(3 + Add(1, 2)) == 3 + 1 + 2
  @test 3 + Add(1, 2) isa Add
  @test materialize(2 * ∏([1, 2])) == 2 * prod([1, 2])
  @test 2 * ∏([1, 2]) isa α
  @test 2 * ∏([1, 2]) isa α{<:∏}
  @test 2 * ∏([1, 2]) isa α{∏{Int}}
  @test ∏([1, 2]) + ∏([3, 4]) == ∑([∏([1, 2]), ∏([3, 4])])
  @test ∏([1, 2]) + ∏([3, 4]) isa ∑
  @test materialize(∑(∏([1, 2]) + ∏([3, 4]))) == sum([prod([1, 2]), prod([3, 4])])
  @test ∏([1, 2]) + ∏([3, 4]) == ∑([∏([1, 2]), ∏([3, 4])])
  @test ∏([1, 2]) + ∏([3, 4]) isa ∑
  @test ∏([1, 2]) - ∏([3, 4]) isa ∑
  @test materialize(∏(["X", "Y", "Z"])) == "XYZ"
  @test ∏(["X", "Y", "Z"]) isa ∏
  @test materialize(∏() * "X" * "Y" * "Z") == "XYZ"
  @test ∏() * "X" * "Y" * "Z" == ∏(["X", "Y", "Z"])
  @test ∏() * "X" * "Y" * "Z" isa ∏
  @test 2∏() * "X" * "Y" == 2∏(["X", "Y"])
  @test 2∏() * "X" * "Y" isa α{<:∏}
  @test 2∏() * "X" * "Y" isa α{∏{String}}
  @test 2∏() * "X" * "Y" isa α{∏{String},Int}
  @test 2∏(["X"]) * 3∏(["Y"]) == 6∏(["X", "Y"])
end