blockindexing.r

# Helper function to calulate index price for the representative curve.
# gorl =1 implies all the Rep_KVA values above the given kva value shall be summed up and sent.
# gorl = -1 implies all the Rep_KVA below the passed kva value shall be deducted with the kva and their summation shall be sent.
AreaIntegrator <- function(A,kva,gorl=1)
{
  if(gorl==1)
  {
    B = fn$sqldf("select * from A where Rep_KVA>$kva")
  }
  else if(gorl==-1)
  {
    B = fn$sqldf("select * from A where Rep_KVA<$kva")
  }
  if(length(B$Day)!=0)
  {
    B$Row_no = seq(1,length(B$Day))
    B$ModifiedRep = ifelse((gorl*(B$Rep_KVA-kva)>0),(gorl*(B$Rep_KVA-kva)),0)
    return(sum(B$ModifiedRep))
  }
  else
  {
    return(0.0)
  }
}


# This function rearranges the representative KVA value table in order i.e. from Sunday 00:00:00 to Saturday 12:45:00 
Rearranger <- function(final)
{
  # Rearranging the distinct Days in our sequence
  weekdaylist= c( "Sunday","Monday","Tuesday", "Wednesday","Thursday","Friday","Saturday" )
  yyy =sqldf("select distinct Day from final")
  tutu=NULL
  for (iii in 1:length(yyy$Day))
  {
    day = weekdaylist[iii]
    if(grep(day, yyy$Day)!=0)
    {
      tut=data.frame(Day=day)
    }
    else
      tut=NULL
    tutu= rbind(tutu,tut)
  }
  yyy=tutu

  # Created Dataset based on Daywise grouping
  h3=NULL
  for (h1 in 1:length(yyy$Day))
  {
    h0 = fn$sqldf("select * from final where Day=(select Day from yyy limit $h1-1,1)")
    h2 = fn$sqldf("select * from h0 order by Time")
    plot(h2$Time,h2$Rep_KVA,xlab="Time on Day",ylab="Rep_KVA",las=2)
    lines(h2$Time,h2$Rep_KVA,xlab="Time on Day",ylab="Rep_KVA",las=2)

    # Analysing this weekday for apt pricing plan
    h3 = rbind(h3,h2)
  }

  # Labelling of the data
  h3=transform(Label=paste0(h3$Day," - ",h3$Time),h3)
  h3$Label=as.character(h3$Label)

  # One week representation of Data
  Row_no =seq(1,length(h3$Day))
  h3 = cbind(Row_no,h3)
  plot(h3$Row_no,h3$Rep_KVA)
  lines(h3$Row_no,h3$Rep_KVA)
  return(h3)
}


# Custom function to create an optimised representative 1 week (all time) table for the available data
CentralAnalyser <- function(t)
{
  tt = sqldf("select distinct Time from t")
  tn = as.numeric(sqldf("select count(distinct Day) from t"))
  final=data.frame(  Day=as.Date(character()),
                     Time=character(),
                     Rep_KVA=numeric(),
                     stringsAsFactors=FALSE)

  for (j in 1:length(tt$Time))
  {
    e = fn$sqldf("select * from t where Time=(select Time from tt limit $j-1,1)")
    for (i in 1:tn)
    {
      f = fn$sqldf("select * from e where Date%$tn=$i-1 and Global_active_power!='NA'")
      plot(f$Date,f$KVA,xlab='Day No.',ylab='KVA', main = "Graphical representation of KVA value in 6 months for particular time.")
      lines(f$Date,f$KVA)
      # Finding the best optimized central tendency for the dataset (for particular weekday and particular time)
      nn = as.numeric(sqldf("select count(*) from f"))
      hippo =d = NULL
      for (i1 in seq(1,nn))
      {
        d=NULL
        d1 = as.numeric(fn$sqldf("select KVA from f limit $i1-1,1"))
        for (j1 in 1:nn)
        {
          d2 = as.numeric(fn$sqldf("select KVA from f limit $j1-1,1"))
          d[j1]=ifelse(d2-d1>0,(d2-d1)*.16*.25,0)
        }
        f$hippo[i1]=sum(d)+d1*(nn)*0.1*.25
        #k = k+1
        #f$hippo[k]=hippo[i1]
        #print(paste0("Hippo: for ", d1, " and ", d2,"  = ", hippo[i1]))
      }
      if(!is.na(f$hippo[1]))
      {
        rep_data=sqldf("select * from f where hippo=(select min(hippo)from f) limit 1")
        repkva = as.numeric((sqldf("select KVA from rep_data")))
        reptime = as.character(sqldf("select Time from rep_data"))
        repweekday = weekdays(rep_data$Date)
        final=rbind(final,data.frame(Day=repweekday,Time=reptime,Rep_KVA=repkva))
      } else
      {
        final=rbind(final,data.frame(Day=f$Day[1],Time=f$Time[1],Rep_KVA=0))
      }
    }
  }
  return(final)
}

# Determining the limiting values for the block size beyond which block size (height) is not econoomical.
MinKvaCalculation <- function(tob)
{
  library(sqldf)
  nom=as.numeric(sqldf("select count(*) from tob"))
  hr=nom/4
  divtocover=hr*4
  mink=as.numeric(sqldf("select Min(Rep_KVA) from tob"))
  maxk=as.numeric(sqldf("select Max(Rep_KVA) from tob"))
  mkva=(c(floor(mink),ceiling(maxk)))
  p=mkva[[1]]
  ond=0
  # Plus 1 excluded
  while(p>=mkva[[1]] && p<=(mkva[[2]]))
  {
    if((p*divtocover*.10+AreaIntegrator(tob,kva=p)*.16)<=AreaIntegrator(tob,kva=0)*.16)
    {
      ond[1]=AreaIntegrator(tob,kva=p,gorl=-1)
      p=p+1
    }
    else{
      ond[1]=AreaIntegrator(tob,kva=p,gorl=-1)
    if(ond==0)
        p= p +1
    else
      break}
  }

  return(p)
}

# Finding the possible feasible solutions for eq. k1+k2+k3+.+kx.=n.
# where, n = MinKvaValCalculation
# x = number of plans that the user added.
# SOlution : (n+x-1)!/((x-1)!(n!))
CombinationFinder <- function(tub,uc)
{
  maxsizeofblock = MinKvaCalculation(tub)
  ul = length(uc$Hour)-1
  pc=factorial(ul+maxsizeofblock)/(factorial(ul)*factorial(maxsizeofblock))
  pc=(as.numeric(pc))
  return(pc)
}

# Gives all possible block plan possible irrespective of the weekday and time been declared as an array.
ArrayCreator <- function(hour,set)
{
  
  # No of hour variation
  hv=24-hour+1
  
  # No of day variation
  dv=8-set
  
  # No of divisions in each block 
  m = as.numeric(hour*4)
  arr=array(0,dim=c(hv*dv,672))
  
  for( i in 1:dv)
  {  
    for(j in 1:hv)
    { 
      ll=(j-1)*4+(i-1)*96
      k=1
      while (1<=k && k<=set && ll<=672)
      {
        arr[j+(i-1)*hv,(ll+1):(ll+m)]=1
        ll=ll+96
        k=k+1
      }
    }
  }
  return(arr)
}

# This returns a default block to be taken based on user's input (if the weekday or the time is not defined default value taken is Monday and default time is 6:00 am in the morning)
ArrayCreatorNew <- function(hour,set,weekday,time)
{
    # No of hour variation
    hv=24-hour+1

    # No of day variation
    dv=8-set

    # No of divisions in each block 
    m = as.numeric(hour*4)
    arr=array(0,dim=c(hv*dv,672))

    for( i in 1:dv)
    {
      for(j in 1:hv)
      {
        ll=(j-1)*4+(i-1)*96
        k=1
        while (1<=k && k<=set && ll<=672)
        {
          arr[j+(i-1)*hv,(ll+1):(ll+m)]=1
          ll=ll+96
          k=k+1
        }
      }
    }

        v = dim(arr)[1]/dv


    v1 = seq(1,v)
    v2 = seq(v+1,2*v)
    v3 = seq(2*v+1,3*v)

    ii =-2
    jj = -2

    if(set==7)
      ii = -1

    if (is.na(time))
      jj = 0
    else
      jj = as.integer(time) + 1

    ii = ifelse (is.na(weekday),-1,(ifelse ((weekday=='Sunday'),1,ifelse ((weekday=='Monday'),2,ifelse ((weekday=='Tuesday'),3,0)))))

    if(ii==-2 && set ==5)
    {
      ii = 2
      weekday = 'Monday'
    }
    if (jj==0 && hour !=24)
    {
      jj = 7
      time=7-1
    }

    if(ii==1)
      ii=v1
    else
    {
        if(ii==2)
            ii=v2
        else
        {
            if(ii==3)
              ii=v3
        }
    }

  u = 0
  u1 = 0

  if(length(ii)>1)
  {
    u = array((arr[ii,]),dim=c(length(ii),672))
    u1 = array((u[jj,]),dim=c(1,672))
    return(u1)
  }
  else
  {
    if(ii==-1)
    {
      if (jj==0)
        u = array((arr[,]),dim=c(1,672))
      else
      {
        u = array((arr[jj,]),dim=c(1,672))
      }
      return (u)
    }
    else
    {
      u = array((arr[ii,]),dim=c(length(ii),672))

      if(jj ==0)
        return(u)
      else
      {
        u1 = array((u[jj,]),dim=c(1,672))
        return(u1)
      }
    }
  }

}

# This result gives us the possible combinations the block may have based on the different plans the user have provided.( A cross relation between certain plans are includede in this. eg. Combination with plan1 and plan2 block to be taken. )
# Basically this results an array of all feasible clock structue combinations with different plans as units.
ProductFinder <- function(co)
{
  l = length(co$Hour)
  arow=c()
  for (i in 1:l)
  {
    hr = co$Hour[i]
    st = co$Set[i]
    arow[i]=dim(ArrayCreator(hr,st))[1]
  }
  product=prod(arow)
  return(product)
}

FinalProcessing <- function(teb,co)
{
  library(abind)
  ab=NULL
  l = length(co$Hour)
  arow=itneed=c()
  for (i in 1:l)
  {
    hr = co$Hour[i]
    st = co$Set[i]
    arow[i]=dim(ArrayCreator(hr,st))[1]
  }
  product = ProductFinder(co)
  for (i in 1:l)
  {
    itneed[i]=product/arow[i]
  }
  for (o in 1:length(itneed))
  {
    p = ArrayCreator(co$Hour[o],co$Set[o])
    xo = do.call("rbind", replicate(itneed[o], p, simplify = FALSE))
    ab=abind(ab,xo,along =1)
  }
  return(ab)
}

# To process the result for user specific input giving result non optimized result.
FinalProcessingNew <- function(teb,co)
{
  library(abind)
  ab = NULL
  l = length(co$Hour)
  arow = c()

  for (o in 1:l)
  {
    p = ArrayCreatorNew(co$Hour[o],co$Set[o],as.character(co$Day[o]),as.character(co$Time24[o]))
    ab=abind(ab,p,along =1)
  }
  return (ab)
}

# This function will prepare the complete block structure by summing block heights of each plan one by one together and finally calculate the BlockPrice, Index Price for each combination hence total price for each possibility.
SummationArray <- function(h3,k,arr,co)
{
  Prod=ProductFinder(co)
  sm=matrix(0,nrow=Prod,ncol=672)
  for(i in 1:Prod)
  {
    for(j in 1: length(co$Hour))
    {
      sm[i,]= (sm[i,]+arr[i+Prod*(j-1),]*k[j])
    }
  }
  BSC=matrix(0,nrow=1,ncol=673)
  loo=matrix(h3$Rep_KVA,nrow=1,ncol=672)
  loo = do.call("rbind", replicate(Prod, loo, simplify = FALSE))
  t = matrix(0,nrow=Prod,ncol=672)
  for( u in 1:Prod)
  {
    for (j in 1:672)
    {
      p = (loo[u,j]-sm[u,j])
      if(p>=0 && (!is.na(p)))
      {
        t[u,j] = loo[u,j]-sm[u,j]
      }
      else t[u,j] =0
    }
  }
  BPrice=rowSums(sm, na.rm=TRUE)*.1*.15
  IPrice=rowSums(t,na.rm=TRUE)*.16*.15
  Price=(BPrice+IPrice)
  minpindex=match(min(Price),Price)
  BSC[1,1:672]=sm[minpindex,]
  BSC[1,673]=min(Price)

  return(BSC)
}

# Summation of blocks plans based on user specific information.
SummationArrayNew <- function(h3,k,arr,co)
{
  sm=matrix(0,nrow=length(k),ncol=672)
  for(i in 1:length(k))
  {
      sm[i,]= (sm[i,]+arr[i,]*k[i])
  }

  BSC=matrix(0,nrow=1,ncol=673)
  loo=matrix(h3$Rep_KVA,nrow=1,ncol=672)
  loo = do.call("rbind", replicate(dim(sm)[1], loo, simplify = FALSE))
  t = matrix(0,nrow=dim(sm)[1],ncol=672)
  for (i in 1:length(k))
  {
    for (j in 1:672)
    {
      p = (loo[i,j]-sm[i,j])
      if(p>=0 && (!is.na(p)))
      {
        t[i,j] = loo[i,j]-sm[i,j]
      }
      else t[i,j] =0
    }
  }

  BPrice=rowSums(sm, na.rm=TRUE)*.1*.15
  IPrice=rowSums(t,na.rm=TRUE)*.16*.15

  Price=(BPrice+IPrice)
  minpindex=match(min(Price),Price)
  BSC[1,1:672]=sm[minpindex,]
  BSC[1,673]=min(Price)

  return(BSC)
}

# This includes the instructions to be taken in the beginning i.e., to fetch data , prepare it for processing and process it based on our choice.
# Basically an integration of the entire functions together
BlockIndexingInitialInstructions <- function()
{
    library(sqldf)

    # Reading content from file
    x1<-read.table("//home//innovator//Documents//Internship//Data//household_power_consumption.txt",header=TRUE,sep=';',dec = '.',stringsAsFactors=FALSE)

    # Performing necessary conversions
    x1$Date<- as.Date(x1$Date,"%d/%m/%Y")
    x1$Time  =  format(as.POSIXlt(x1$Time, format = "%H:%M:%S"), format="%H:%M:%S")
    x1$Global_active_power = as.numeric(x1$Global_active_power)
    x1$Global_reactive_power = as.numeric(x1$Global_reactive_power)

    # Defining KVA
    x1$KVA = sqrt(x1$Global_active_power*x1$Global_active_power+x1$Global_reactive_power*x1$Global_reactive_power)

    # Creating Day column
    x1$Day = weekdays(x1$Date)


    ## Creating sample pools for six month data
    # The starting date
    first<- as.Date('16-12-2008',format = '%d-%m-%Y')
    # The Ending Date
    second <- as.Date('30-05-2009',format = '%d-%m-%Y')

    y1 = fn$sqldf("select * from x1 where Date>=$first and Date<=$second")

    ## Creating sample pools of data at time 0, 15min, 30min, 45min
    z1  = fn$sqldf("select * from y1 where substr(Time,4,2)='00' or substr(Time,4,2)='15' or substr(Time,4,2)='30' or substr(Time,4,2)='45'")

    final2 = NULL
    final1=CentralAnalyser(z1)
    final2 = Rearranger(final1)

    #UserChoiceTables
    choice = NULL
    print(paste0("Enter number of blocks to be taken of your choice. "))
    block1= as.integer(readline())
    for (g in 1:block1)
    {
    print(paste0("Enter hour plan to be taken for block choice ",g))
    choice$Hour[g] = as.numeric(readline())
    print(paste0("Enter set plan to be taken for block choice ",g))
    choice$Set[g] =  as.numeric(readline())
    if(choice$Hour[g]!=24)
    {
        print(paste0("Enter time plan to be taken for block choice ",g))
        choice$Time[g] =  as.character(readline())
    }
    else
    {
        choice$Time[g] =  NA
    }
    if (choice$Set[g]==5)
    {
        print(paste0("Enter Start Day plan to be taken  for block choice ",g))
        choice$Day[g] =  as.character(readline())
    }
    else
    {
        choice$Day[g] =  NA
    }
    }

    choice=data.frame(choice)

    #Updating price column
    choice$Price = choice$Hour*choice$Set*.10*0.4

    for (i in 1:length(choice$Hour))
    {
      choice$Day[i] = ifelse(as.character  (choice$Day[i])=="","NA",as.character(choice$Day[i]))
      choice$Time[i] = ifelse(as.character  (choice$Time[i])=="","NA",as.character(choice$Time[i]))
    }

    #Updating time in 24 format column
    for (iii in 1:length(choice$Time))
    {
    T = as.character(choice$Time[iii])
    if(is.na(T) || is.null(T))
    {
    choice$Time24[iii]=T
    }
    else
    {
        time1=as.numeric(strsplit(T,':')[[1]][1])
        time2=(substring(strsplit(T,':')[[1]][2],1,2))
        if(length(grep('p', g, ignore.case = TRUE)))
        {time1 = as.numeric(time1 +12)}
        else { if(!(time1==10 || time1==11 || time1==12))
        time1 = paste0("0",time1)}
        choice$Time24[iii]=time1
    }
    }

    co = choice
    # Calculating Value for the Optimized Result
    optim = BlockIndexingOptimized(final2,co)
    print(optim)
    # Calculating Value for the Non-Optimized Result
    nonoptim = BlockIndexingOptimized(final2,co)
    print(nonoptim)

}

# Gives result for the optimized case where user's plan details about just Hour (24,16,8) and Set (5,7) is considered
BlockIndexingOptimized<- function(final2,co)
{
    library(rPython)
    python.load('multichoose.py')
    balls = MinKvaCalculation(tob=final2)
    boxes=length(co$Hour)
    op = CombinationFinder(tub=final2,uc=co)
    com= python.get(paste0("multichoose(",boxes,",",balls,")"))

    BSC=matrix(0,nrow=op,ncol=673)
    #pro = ProductFinder(co=choice)
    #n1 = power(n+1,length(choice$Hour))

    ouo=FinalProcessing(teb=c(final2),co=co)

    # Creating loop
    for( i in 1:op)
    {
      l=SummationArray(h3=final2,k=com[[i]],arr=ouo,co=co)
      BSC[i,] = l
    }

    Rows = seq(1:672)
    blockstruct = data.frame(Rows)
    MinP = min(BSC[,673])
    blockstruct$X = BSC[match(MinP,BSC[,673]),1:672]
    par(mfrow=c(3,1))
    plot(final2$Row_no,final2$Rep_KVA)
    lines(final2$Row_no,final2$Rep_KVA)
    barplot(space=10,blockstruct$Row,height=c(blockstruct$X))
    barplot(space=10,final2$Row_no,height=final2$Rep_KVA)
    print(paste0("Total Pay Amount for the optimized graph: $ ",MinP))

    return(MinP)

}

# Gives result for the Non Optimized Plan considering all the user's choices and default values in case user input is unavailable.
BlockIndexingNonOptimized<- function(final2,co)
{
    library(rPython)
    python.load('multichoose.py')
    balls = MinKvaCalculation(tob=final2)
    boxes=length(co$Hour)
    op = CombinationFinder(tub=final2,uc=co)
    com= python.get(paste0("multichoose(",boxes,",",balls,")"))

    BSC=matrix(0,nrow=op,ncol=673)

    # For optimized Result

    ouo=FinalProcessingNew(teb=final2,co=co)

    # Creating loop
    for( i in 1:op)
    {
      l=SummationArrayNew(h3=final2,k=com[[i]],arr=ouo,co=co)
      BSC[i,] = l
    }

    Rows = seq(1:672)
    blockstruct = data.frame(Rows)
    MinP = min(BSC[,673])
    blockstruct$X = BSC[match(MinP,BSC[,673]),1:672]
    par(mfrow=c(3,1))
    plot(final2$Row_no,final2$Rep_KVA)
    lines(final2$Row_no,final2$Rep_KVA)
    barplot(space=10,blockstruct$Row,height=blockstruct$X)
    barplot(space=10,final2$Row_no,height=final2$Rep_KVA)
    print(paste0("Total Pay Amount for the non-optimized graph: $ ",MinP))

    return(MinP)
}

multichoose.py

def multichoose(n,k):
    if k < 0 or n < 0: return "Error"
    if not k: return [[0]*n]
    if not n: return []
    if n == 1: return [[k]]
    return [[0]+val for val in multichoose(n-1,k)] + \
        [[val[0]+1]+val[1:] for val in multichoose(n,k-1)]