x <- 10
x = 10
x <<- 10 # zobrazeni vnitrni promenne funkce


# Urcete hmotnost kapaliny (o hustote rho = 0.887 kg/L) v plechovem sudu o rozmerech d = 0.5 m a v = 1.0 m.
ww = function(d,v,rho){
  V <<- pi*(d/2)^2 * v # [m3]
  m <- V * rho * 1000  # [kg]
  return(m)
}
m = ww(0.5,1.0,0.887); m # [m3]
V  # [m3]


### matice ###################################################################

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)

x12 = rbind(x1,x2)
x12 = cbind(x1,x2)
is.matrix(x12)
nrow(x12)
ncol(x12)

x12[order(x1),]
x12[order(x2),]
x12[,1]
x12[,"x"]

colnames(x12)
colnames(x12)=c("X","Y")
rownames(x12)
rownames(x12) = LETTERS[1:length(x1)] 


### dataframes ##################################################################################

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x12 = rbind(x1,x2,z)
x12 = cbind(x1,x2,z)
nrow(x12)
ncol(x12)
is.data.frame(x12)
is.matrix(x12)

as.data.frame(x12)
is.data.frame(x12)
colnames()
rownames()


### lists ###########################################################################################

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)

x123 = list(x1,x2,x3)
is.list(x123)
names(x123) = c("x1","x2","x3") 

unlist(x123)
do.call(rbind,x123)
do.call(cbind,x123)

stack(x123)


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


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


### vkladani hodnot z konzole

widmark <- function(){ 
  print("Výpocet hladiny alkoholu v krvi podle Widmarka")
  print("Data zadavejte pres konzoli !!!")
  sex<-readline(prompt="Jste muz (m) nebo zena(z)? " )
  if (sex == "m") {BW = 0.58; MR = 0.015} 
  if (sex == "z") {BW = 0.49; MR = 0.017}
  Wt <-readline(prompt="Jaka je Vase hmotnost (v kg)? " )
  Wt = as.numeric(Wt)
  # print("Jeden standardni alkoholicky napoj odpovida 10 g resp. 12.7 ml cisteho ethanolu")
  # SD <-readline(prompt="Kolik jste vypil standardnich alkoholickych napoju? " )
  # prepocitat na piva, vina a panaky
  p0 <-readline(prompt="Kolik sklenic nealkoholickeho piva (0.5 l) jste vypil ? " )
  p10 <-readline(prompt="Kolik sklenic 10 stupnoveho piva (0.5 l) jste vypil ? " )
  p12 <-readline(prompt="Kolik sklenic 12 stupnoveho piva (0.5 l) jste vypil ? " )
  vb <-readline(prompt="Kolik sklenicek bileho vina (0.2 l) jste vypil ? " )
  vc <-readline(prompt="Kolik sklenicek cerveneho vina (0.2 l) jste vypil ? " )
  p30 <-readline(prompt="Kolik panaku 30% alkoholu (0.05 l) jste vypil ? " )
  p45 <-readline(prompt="Kolik panaku 45% alkoholu (0.05 l) jste vypil ? " )
  # SD = as.numeric(SD)
  p0 = as.numeric(p0) * 2 # [g 100% ethanolu]
  p10 = as.numeric(p10) * 16 # [g 100% ethanolu]
  p12 = as.numeric(p12) * 20 # [g 100% ethanolu]
  vb = as.numeric(vb) * 16 # [g 100% ethanolu]
  vc = as.numeric(vc) * 19 # [g 100% ethanolu]
  p30 = as.numeric(p30) * 12 # [g 100% ethanolu]
  p45 = as.numeric(p45) * 18 # [g 100% ethanolu]
  SD = sum(c(p0,p10,p12,vb,vc,p30,p45))/10
  DP <-readline(prompt="Pred kolika hodinami? " )
  DP = as.numeric(DP)
  EBAC = ((0.806*SD*1.2/(BW*Wt)) - MR*DP)*10   # [promile, resp. g/l]
  # if (EBAC <= 0) {EBAC = 0}
  cat("Mate v krvi", round(EBAC,2) , "promile (g/l) alkoholu.")
  # print(paste("Mate v krvi", round(EBAC,2) , "promile (g/l) alkoholu."))
}
widmark()

# https://en.wikipedia.org/wiki/Blood_alcohol_content
# https://en.wikipedia.org/wiki/Standard_drink
# https://www.bezpecnecesty.cz/cz/alkohol-kalkulacka/algoritmus-vypoctu
# https://autobild.pluska.sk/poradca/chcete-vypocitat-hladinu-alkoholu-poradime-vam
#

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


#### Fyzikalni konstanty  ##########################################################################################################

### Ludolfovo cislo ####
pi
print(pi, digits=20)

### Eulerovo cislo ####
exp(1)
print(exp(1), digits=20) 

library(marelac)
data(Constants)
Constants

Constants$gasCt1
Constants$gasCt2
Constants$gasCt3
R = as.numeric(Constants$gasCt2[1]); R # hodnota
Constants$gasCt2[2] # jednotka


### Vylepseni
mm = "gas"
consts = do.call(rbind,Constants) # list to dataframe
MM = sapply(consts[,3], function(xxx) grepl(mm, xxx, fixed = TRUE))
SMM = consts[which(MM==TRUE),]
SMM
# zadat poradi polozky ve vystupu
SMM[2,] # pro 1 polozku SMM[1]
R = as.numeric(SMM[2,1]); R # hodnota
SMM[2,2] # jednotka


library(constants)
codata
lookup("planck", ignore.case=TRUE)
lookup("avogadro", ignore.case=TRUE)
lookup("faraday", ignore.case=TRUE)
lookup("molar", ignore.case=TRUE)
lookup("molar mass", ignore.case=TRUE)
lookup("light", ignore.case=TRUE)
lookup("boltzmann", ignore.case=TRUE)
lookup("mass", ignore.case=TRUE)
lookup("mass constant", ignore.case=TRUE)
lookup("proton mass", ignore.case=TRUE)
lookup("electron mass", ignore.case=TRUE)

# Univerzalni plynova konstanta
lookup("gas", ignore.case=TRUE)
codata[214,]
R = as.numeric(codata[214,3]) # hodnota
Ru = as.numeric(codata[214,5]) # nejistota
codata[214,4] # jednotka


#### Atomove a molekulove hmotnosti ###################################################################################################################################

library(periodicTable)
data(periodicTable)
periodicTable

library(CHNOSZ)
# atomova hmotnost
AW = mass("Ca"); AW
# molarni hmotnost
MW = mass("CaCO3"); MW
# pocet atomu prvku ve vzorci
ce = count.elements("CaCO3"); ce
as.chemical.formula(ce, drop.zero = TRUE)

library(NORRRM)# archiv
# atomova hmotnost
data(AtomWeight)
AtomWeight
AW = AtomWeight["Si","AWeight"]; AW
el = c("Ca","Si","Al")
AW = AtomWeight[el,'AWeight']; AW
names(AW) = el
AW

# prepocet na oxidy
data(OxiWeight)
OxiWeight
OxiWeight ['Si','OWeight']


library(IsoSpecR)
data(isotopicData)
isotopicData


#### Prevody jednotek ##########################################################################################################

library(scales) # %
percent(0.05)

library(pracma) # uhly vs radiany
deg2rad(120) 
rad2deg(3.14)


library(measurements)
conv_unit_options # pouzivane jednotky

# length
conv_unit(1, from="mm", to="m")
conv_unit(1, from="nm", to="angstrom")

# pressure
conv_unit(101, from="Pa", to="mmHg")
conv_unit(1, from="atm", to="Pa")

# temperature
conv_unit(100, from="C", to="K")

# energy
conv_unit(500, from="cal", to="J")
conv_unit(1, from="erg", to="J")

#volume
conv_unit(100, from="l", to="m3")

# mass
conv_unit(100, from="g", to="mg")
conv_unit(100, from="ug", to="g")

# speed
conv_unit(1, from="km_per_hr", to="m_per_sec")
conv_unit(1, from="light", to="km_per_hr")
conv_unit(1, from="light", to="m_per_sec")


### multiunits
conv_multiunit(x = 1, from="ug / l", to="g / m3")
conv_multiunit(x = 1, from="mg / l", to="kg / m3")
conv_multiunit(x = 1, from="cal / kg", to="J / g")


### Plynný systém mení svuj objem o 1200 ml za konstantního vnejsího tlaku 30 atm.Jakou práci vykoná plyn pri expanzi? Vyjádrete v ruzných energetických jednotkách.
# W = p . dV
dV = 1200 # [ml]
dV = conv_unit(dV, from="ml", to="m3"); dV
p = 30 # [atm] to [Pa]
p = conv_unit(p, from="atm", to="Pa"); p
W = p*dV; W # [J]
W1 = conv_unit(W, from="J", to="cal"); W1 # [cal]
W2 = conv_unit(W, from="J", to="Wsec"); W2 # [W.s]
W3 = conv_unit(W, from="J", to="erg"); W3 # [erg]


### Kolik tepla se uvolní pri pruchodu 26.43 coulombu elektrického proudu vodicem pri potenciálovém spádu 2.432 V.
# 64.33 J, 64.33e7 erg, 15.37 cal
# Q = U . I . t = U . q 
q = 26.43 # [C]
U = 2.432 # [V]
Q = U*q; Q # [J]
Q1 = conv_unit(Q, from="J", to="cal"); Q1 # [cal]
Q3 = conv_unit(Q, from="J", to="erg"); Q3 # [erg]


# Prevod casovych jednotek

library(datetime)
xx = as.minute(7658) # 7658 min
as.second(xx)
as.hour(xx)
as.day(xx)


### Pokus trval 7658 minut. Kolik je to dní, hodin a minut?
tt = 7658 # [min]

library(lubridate)
seconds_to_period(tt*60)
library(datetime)
seconds_to_period(as.second(as.minute(tt)))

# opacny vypocet
library(datetime)
dur <- duration(week = 0, day = 5, hours = 7, minutes = 38)
as.numeric(dur, "minutes")


#### nasobky jednotek SI
library(sitools)
f2si(10000)
f2si(0.023, unit="l")
f2si(3.5e-5, unit="mol")
numbers <- c(1e5, 3.5e-12, 0.004) 
f2si(numbers, unit="g")


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


#### Nacitani externich dat ##########################################################################################################

pisky<- read.table("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.txt", header=T)
pisky<- read.csv("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.csv", header=T)

library(openxlsx)
shn = getSheetNames("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx")
shn
pisky <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "DVPisky", startRow = 1, colNames = TRUE, rowNames = TRUE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
summary(pisky)
pisky = pisky[,-1]


#### NA and data imputation #####################################################################################################

### vektory

x = pisky[,"Sr"]
length(x)
x[4] = NA
any(is.na(x))
is.na(x)

library(jwutil) # rel. podil hodnot ve vektoru
propIsNa(x)

xx = na.omit(x)
length(xx)


### Matice a dataframes

x = pisky
x[5,5]<- NA
any(is.na(x))
is.na(x)

library(ryouready)
count_na(x, along = 1)
count_na(x, along = 2)
# along: 1 = rows, 2=columns

library(jwutil) # rel. podil hodnot ve sloupcich
propNaPerField(x)

library(agricolae)
delete.na(x, alternative="less")
delete.na(x, alternative="greater")
# "less"    = vynechava sloupce
# "greater" = vynechava radky


### Imputace

library(DescTools)
ZeroIfNA(x)
NAIfZero(x)
#
BlankIfNA(x, blank="")
BlankIfNA(x, blank="UDL")
NAIfBlank(x)

library(jwutil) # NA na 0
zero_na(x)
zero_na(x, na_ish = TRUE)


library(DescTools) # nahrazeni medianem sloupce
x5 = Impute(x[,5], FUN = function(x) median(x, na.rm = TRUE))
rbind(x[,5],x5,pisky[,5])
x[,5] = x5 # nahrazeni

library(compositions) # nahrazeni nasobkem DL
x0 = ZeroIfNA(x) # library(DescTools)
dl=rep(0.1,ncol(x)) # vektor DL
xdl = zeroreplace(x0[,5],d = 0.1,a = 1/2)
x[,5] = xdl # nahrazeni

library(na.tools)
all_na(x)
any_na(x)
is_na(x)
which_na(x[,5])
which_na(x[5,])

xdl = na.replace(x[,5], .na=0.05)
x[,5] = xdl # nahrazeni

na.constant(x[,5], 0.05)
na.max(x[,5])
na.min(x[,5])
na.mean(x[,5])
na.median(x[,5])
na.quantile(x[,5], prob=0.4 )
na.mode(x[,5])
na.most_freq(x[,5])
na.inf(x[,5])
na.neginf(x[,5])
na.true(x[,5])
na.false(x[,5])
na.zero(x[,5])

set.seed(123)
xb = na.bootstrap(x[,5])
xr = na.resample(x[,5])
rbind(x[,5],xb,xr,pisky[,5])



### ZAOKROUHLOVANI ##########################################################################################

# Zaokrouhlovani na dany pocet desetinnych mist #####
round(1234.567,2)
round(123.456,digits=2)

library(guf)
round_something(123.456, decimals = 2)
round_something("123.456", decimals = 2)

library(PKNCA)
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "e")
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "x10^")
roundString(3141.59265, digits = 3, sci_range = Inf, sci_sep = "e")

# Pocet desetinnych mist daneho cisla
decimalplaces <- function(x) { # error due to rounding double precision floating point numbers -- replacing the x %% 1
  if (abs(x - round(x)) > .Machine$double.eps^0.5) {
    nchar(strsplit(sub('0+$', '', as.character(x)), ".", fixed = TRUE)[[1]][[2]])
  } else { return(0) } }

decimalplaces(3.14) 
decimalplaces(3.1415) 
decimalplaces(3.14159265)


#### Zaokrouhlovani na dany pocet platnych cislic #########
signif(1234.567,4)
signif(123.456,digits=4)

library(PKNCA)
signifString(3141.59265, digits = 1, sci_range = 0, sci_sep = "e")
signifString(3141.59265, digits = 1, sci_range = 0,sci_sep = "x10^")
signifString(3141.59265, digits = 1, sci_range = Inf, sci_sep = "e")

# Urceni poctu platnych cislic daneho cisla
sigdig <- function(x) length(gregexpr("[[:digit:]]", as.character(x))[[1]]) 

sigdig(3.14) 
sigdig(3.1415) 
sigdig(3.14159265)


### Relativni chyba zaokrouhleneho cisla
X = 0.053450
n = 2 # zvoleny pocet platnych cislic
library(PKNCA)
# nacist funkci sigdig()
A = as.numeric(substring(signifString(X, digits = sigdig(X), sci_range = 0, sci_sep = "e"), 1, 1))
s = 1/(2*A*10^(n-1)); s
library(scales)
percent(s)
library(jwutil)
percentize(s)


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

### if - else ###

### Vypocítejte vazebnou energii atomového jádra pomoci Bethe - Weiszackerovy rovnice
#     B = 14.0 - 13.1*A^(2/3) + 0.585*Z*(Z-1)/A^(1/3) - (18.1*(A-2*Z)^2)/A + C/A
# A - nukleonove cislo, Z - atomove cislo, C - konstanta, pro jadra se sudym poctem protonu a neutronu je rovna 132, 
# pro jadra s lichym poctem protonu a neutronu je rovna -132, pro lichy pocet protonu a sudy pocet neutronu a naopak je rovna 0.

# alpha particle
A = 4 
Z = 2

# U-235
A = 235 
Z = 92

N = A-Z
if (N%%2==0 & Z%%2==0){C = 132} else if (N%%2!=0 & Z%%2!=0){C = -132} else {C = 0}
B = 14.0 - 13.1*A^(2/3) + 0.585*Z*(Z-1)/A^(1/3) - (18.1*(A-2*Z)^2)/A + C/A
B # [MeV]
B/A # energie na 1 nukleon [MeV]

# polomer jadra: R = R0 * A^(1/3)
R0 = 1.4e-13 # [cm]
R = R0 * A^(1/3)
R = R*10e-15 
R   # [m]


### Anorganická kvalitativní analýza - 1. skupina

KAA1 <- function(){ 
  print("Kvalitativni analyza - 1. skupina iontu (Ag, Pb, Hg)")
  print("Odpovidejte a (ano) nebo n (ne)")
  print("Pozor. Nez zacnete, prepnete se do konzole.")
  an <- readline(prompt="Chcete pokracovat? " )
  if (an == "n") {print("Tak snad nekdy jindy ...")
  } else {print("Takze zaciname ...")  
    an <- readline(prompt="Vznikla po pridani HCl bila srazenina? " )
    if (an == "n") {print("Neni pritomen zadny prvek 1. skupiny.")
    } else {print("Je pritomno olovo, stribro, rtut, nebo jejich kombinace.")
      an <- readline(prompt="Rozpousti se srazenina zcela v horke vode? " )
      if (an == "a") {print("Je pritomno pouze olovo.")
      } else {
        an <- readline(prompt="Rozpousti se srazenina alespon castecne? " )  
        if (an == "a") {print("Je pritomno olovo, ale take stribro a/nebo rtut.")
        } else {print("Olovo neni pritomno, jsou pritomny stribro a/nebo rtut.")}
        an <- readline(prompt="Rozpousti se srazenina zcela v amoniaku? " )
        if (an == "a") {print("Je pritomno stribro, neni pritomna rtut.")
        } else {
          an <- readline(prompt="Rozpousti se srazenina alespon castecne? " )  
          if (an == "a") {print("Je pritomno stribro, a take rtut.")
          } else {print("Neni pritomno stribro, je pritomna rtut.")}
          an <- readline(prompt="Doslo pridavkem amoniaku ke zcernani srazeniny? " )
          if (an == "a") {print("Rtut je potvrzena.")
          } else {print("To je divne. Nekde se stala chyba ...")}
        } 
      }   
    }
    print("Hotovo. Nashledanou priste.")
  }
}
KAA1()



### while cyklus (while loop) ###

### Vypocitejte molekulove hmotnosti alkanu C1 - C12
n = 1
Mhc = NULL
while(n <= 12){
  hc = n*12.011 + (2*n + 2)*1.008
  Mhc = c(Mhc,hc)
  n = n+1
}
names(Mhc) = c(1:12)
Mhc


### switch
switch(2,"red","green","blue")
switch(1,"red","green","blue")

# Jednoducha kalkulacka pro scitani, odcitani, nasobeni a deleni
simcalc = function(){
  add <- function(x, y) { return(x + y)}
  subtract <- function(x, y) { return(x - y) }
  multiply <- function(x, y) { return(x * y) }
  divide <- function(x, y) { return(x / y)}
  # take input from the user
  num1 = as.integer(readline(prompt="Enter first number: "))
  choice = as.character(readline(prompt="Enter choice[+, -, *, /]: "))
  num2 = as.integer(readline(prompt="Enter second number: "))
  result <- switch(choice, "+"=add(num1, num2), "-"=subtract(num1, num2), "*"=multiply(num1, num2), "/"=divide(num1, num2))
  print(paste(num1, choice, num2, "=", result))
}
simcalc()
# s do.call()
# grafy s moznosti vyberu kombinace promennych

# https://www.datamentor.io/r-programming/switch-function/
# https://www.datamentor.io/r-programming/examples/simple-calculator/
# https://stackoverflow.com/questions/10393508/how-to-use-the-switch-statement-in-r-functions
# https://stackoverflow.com/questions/31156957/use-of-switch-in-r-to-replace-vector-values
# 


### iterace, rovnice, integraly - presnost



### for cyklus (for loop) ###


### Vypocitejte molekulove hmotnosti alkanu C1 - C12

Mhc = NULL
for(n in c(1:12)){
  hc = n*12.011 + (2*n + 2)*1.008
  Mhc = c(Mhc,hc) 
}
names(Mhc) = c(1:12)
Mhc


### Vypocitejte stredni kvadratickou rychlost molekul danych plynu v rozmezi teplot od 150 K do 500 K s krokem 50 K.

el = c("H2","He", "O2", "Kr")
M = c(2.0159, 4.0026, 31.9988, 83.8) # [g / mol]
Ts = seq(from=150, to=500, by=50) # [K]
R = 8314.34 # [J / kmol K]
MKV = NULL
for(tt in Ts){
  mkv= rep(0, length(M))
  for(ii in c(1:length(M))){
    mkv[ii] = (3*R*tt/M[ii])^(1/2)
  }
  MKV = rbind(MKV, mkv)   
}
colnames(MKV) = el
rownames(MKV) = Ts
MKV
plot(0,0,xlim=c(min(Ts),max(Ts)),ylim = c(min(MKV),max(MKV)), type ="n",xlab="T [K]",ylab="str. kv. rychlost [m/s]")
for(ii in c(1:length(M))){points(Ts,MKV[,ii], type="p",col=ii,pch=16)}
legend("topleft",legend=el,col=c(1:length(M)),pch=16,cex=.8)



### apply funkce ###

# https://www.r-bloggers.com/r-tutorial-on-the-apply-family-of-functions/
# http://www.datasciencemadesimple.com/apply-function-r/


### lapply a sapply (vektory, seznamy) ###

### Vypocitejte molekulove hmotnosti alkanu C1 - C12.
# lapply
Mhc = lapply(c(1:12), function(n){n*12.011 + (2*n + 2)*1.008})
Mhc = unlist(Mhc)
names(Mhc) = c(1:12)
Mhc
# sapply
Mhc = sapply(c(1:12), function(n){n*12.011 + (2*n + 2)*1.008})
names(Mhc) = c(1:12)
Mhc


### tapply (matice, dataframes) ###

### Effectiveness of Insect Sprays: The counts of insects in agricultural experimental units treated with different insecticides. 
data(InsectSprays)
InsectSprays
tapply(InsectSprays$count,InsectSprays$spray,mean)
tapply(InsectSprays$count,InsectSprays$spray,sd)
tapply(InsectSprays$count,InsectSprays$spray,length)
tapply(InsectSprays$count,InsectSprays$spray,boxplot.stats)

# by (x, INDICES, FUN, ...)
by(as.data.frame(mtcars), list(cyl,vs), FUN=mean)

# aggregate(x, by, FUN, …, simplify = TRUE, drop = TRUE)
# aggregate(formula, data, FUN, …, subset, na.action = na.omit)
attach(mtcars)
aggdata <-aggregate(mtcars, by=list(cyl,vs), FUN=mean, na.rm=TRUE)
print(aggdata)

# https://www.r-bloggers.com/aggregate-a-powerful-tool-for-data-frame-in-r/
# https://stackoverflow.com/questions/9809166/count-number-of-rows-within-each-group
# https://www.r-bloggers.com/how-to-aggregate-data-in-r/
# https://stats.stackexchange.com/questions/8225/how-to-summarize-data-by-group-in-r
# https://datascienceplus.com/aggregate-data-frame-r/
#

# mean pomoci sapply a lapply (pres indexy)
xx = as.character(unique(InsectSprays$spray))
insp = sapply(xx, function(x){mean(InsectSprays$count[which(InsectSprays$spray==x)])})
insp
insp = lapply(xx, function(x){mean(InsectSprays$count[which(InsectSprays$spray==x)])})
names(insp) = xx
insp = unlist(insp)
insp
# mean pomoci sapply a lapply (pres seznamy)
xx = split(InsectSprays$count, InsectSprays$spray, drop = FALSE)
names(xx) = as.character(unique(InsectSprays$spray))
insp = lapply(names(xx),function(x){mean(xx[[x]])})
names(insp) = names(xx)
unlist(insp)
insp = sapply(names(xx),function(x){mean(xx[[x]])})
insp


### apply a pbapply (matice, dataframes) ###

n<-1000
mat<-matrix(rnorm(n),ncol=10,nrow=100)

# apply()
x<-apply(mat,1,sum)
x<-apply(mat,2,sum)

# pbapply()
library(pbapply) # pridava progress bar v konzoli
x<-pbapply(mat,1,sum)
x<-pbapply(mat,2,sum)

# https://www.r-bloggers.com/pbapply/
# https://www.r-bloggers.com/how-to-add-pbapply-to-r-packages-2/
# 

library(purrr)

library(doBy)
library(AggregateR)


###  Red wine characteristics and quality: dataset related to wine samples from north Portugal.
library(live)
data(wine)
wine
apply(wine, 2, mean)
apply(wine, 2, sd)
apply(wine, 2, length)
apply(wine, 2, boxplot.stats)

library(pbapply)
pbapply(wine, 2, mean)
pbapply(wine, 2, sd)
pbapply(wine, 2, length)
pbapply(wine, 2, boxplot.stats)

# mean pomoci sapply a lapply
wn = sapply(c(1:length(wine[1,])), function(x){mean(wine[,x])})
names(wn) = colnames(wine)
wn
#
wn = lapply(c(1:length(wine[1,])), function(x){mean(wine[,x])})
names(wn) = colnames(wine)
wn

# https://aosmith.rbind.io/2018/06/05/a-closer-look-at-replicate-and-purrr/
# http://www.datasciencemadesimple.com/repeat-and-replicate-function-in-r/
# 
# 
# 


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




library(scales) # %
scales::percent(0.1)
scales::percent(-0.1)
xx =c(1:100)/100
scales::percent(xx)

library(formattable)
x <- c(0.23, 0.95, 0.3)
formattable::percent(x)

# https://stackoverflow.com/questions/7145826/how-to-format-a-number-as-percentage-in-r



library(pracma) # uhly vs radiany
deg2rad(120) 
rad2deg(3.14)


library(sitools)
f2si(10000)
f2si(0.023, unit="l")
f2si(3.5e-5, unit="mol")
numbers <- c(1e5, 3.5e-12, 0.004) 
f2si(numbers, unit="g")


library(quantities)

library(units)


library(convertr)
get_conversion_table()
convert(1:20, "kg", "g") 
convert(1:20, "galUK/min.ft2", "kft/s")


library(constants)

# https://stackoverflow.com/questions/7214781/converting-units-in-r


## vysledna hodnota v zakladnich jednotkach

centi <- 0.01
# objem krychle o delce strany 10 cm (v m3)
volume <- (10 * centi)^3 
volume
# objem kvadru 10 cm x 20 cm x 30 cm (v m3)
volume <- 10*centi * 20*centi * 30*centi
volume

mili = 0.001
# molarni koncentrace roztoku NaCl (v mol/l)
m = 100 * mili # 100 mg NaCl
M = 58.44      # 58.44 g/mol
V = 250 * mili # 250 ml
c = m/(M*V); c



library(measurements)
conv_unit_options # pouzivane jednotky

# length
conv_unit(1, from="mm", to="m")
conv_unit(1, from="nm", to="angstrom")

# pressure
conv_unit(101, from="Pa", to="mmHg")
conv_unit(1, from="atm", to="Pa")

# temperature
conv_unit(100, from="C", to="K")

# energy
conv_unit(500, from="cal", to="J")
conv_unit(1, from="erg", to="J")

#volume
conv_unit(100, from="l", to="m3")

# mass
conv_unit(100, from="g", to="mg")
conv_unit(100, from="ug", to="g")

# speed
conv_unit(1, from="km_per_hr", to="m_per_sec")
conv_unit(1, from="light", to="km_per_hr")
conv_unit(1, from="light", to="m_per_sec")


### Plynný systém mení svuj objem o 1200 ml za konstantního vnejsího tlaku 30 atm.Jakou práci vykoná plyn pri expanzi? Vyjádrete v ruzných energetických jednotkách.
# W = p . dV
dV = 1200 # [ml]
dV = conv_unit(dV, from="ml", to="m3"); dV
p = 30 # [atm] to [Pa]
p = conv_unit(p, from="atm", to="Pa"); p
W = p*dV; W # [J]
W1 = conv_unit(W, from="J", to="cal"); W1 # [cal]
W2 = conv_unit(W, from="J", to="Wsec"); W2 # [W.s]
W3 = conv_unit(W, from="J", to="erg"); W3 # [erg]


### Kolik tepla se uvolní pri pruchodu 26.43 coulombu elektrického proudu vodicem pri potenciálovém spádu 2.432 V.
# 64.33 J, 64.33e7 erg, 15.37 cal
# Q = U . I . t = U . q 
q = 26.43 # [C]
U = 2.432 # [V]
Q = U*q; Q # [J]
Q1 = conv_unit(Q, from="J", to="cal"); Q1 # [cal]
Q3 = conv_unit(Q, from="J", to="erg"); Q3 # [erg]


### 


### multiunits

conv_multiunit(x = 1, from="ug / l", to="g / m3")
conv_multiunit(x = 1, from="mg / l", to="kg / m3")

conv_multiunit(x = 1, from="cal / kg", to="J / g")


conv_dim(x, x_unit, trans, trans_unit, y, y_unit)
conv_dim(x=60, x_unit="ug", trans_unit="mg_per_m3", y=100, y_unit="l")
conv_dim(x_unit="ug", trans=0.6, trans_unit="mg_per_m3", y=100, y_unit="l")

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

library(NISTunits)

ShortenLongNames('kilogram per cubic meter')
ShortenTitles('degree Fahrenheit (something long here)')


NISTsecTOhour(10)
NISTsecTOmin(10)

NISTmillimeterOfMercuryConvtnlTOpascal(10)
NISTmillimeterOfWaterConvtnlTOpascal(10)
NISTatmosphereStdTOpascal(10)
NISTatmosphereTechnicalTOpascal(10)
NISTbarTOpascal(10)

NISTliterTOcubMeter(10)

NISTkTOdegC(10)

NISTkTOdegCinterval(10)

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

library(udunits2)

# This function takes udunits compatible strings and determines whether or not it is possible to convert between them.
ud.are.convertible("miles", "km") # TRUE 
ud.are.convertible("grams", "kilograms") # TRUE 
ud.are.convertible("celsius", "grams") # FALSE 

# Determine whether a unit string is parseable and recognized by the udunits library.
ud.is.parseable("K") # TRUE 
ud.is.parseable(" K ") # TRUE 
ud.is.parseable("miles") # TRUE 

# 
x <- seq(10) 
ud.convert(x, "miles", "km")

x <- c(-40, 0, 100) 
ud.convert(x, "celsius", "degree_fahrenheit")
ud.convert(x, "celsius", "kelvin")


ud.get.name("kg") 
ud.get.symbol("kilogram")

# Retrieve the udunits name or symbol from the database for a given units string.
units.to.display <- c("celsius", # has no name, messed up symbol (maybe a bug in R?)
                      "kg",
                      "hr", # has no symbol 
                      "K", "degrees", "m", "ohm")
lapply(units.to.display, function(xx){cbind(print(ud.get.name(xx)), print(ud.get.symbol(xx))) })


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


library(prettyunits)


## Datum a cas
library(lubridate)
library(datetime)


Sys.Date()
Sys.time()
Sys.timezone(location = TRUE)

# stopky
library(pracma)
tic(gcFirst=FALSE) # start
toc(echo=TRUE)     # stop


library(lubridate)
today(tzone = "")
now(tzone = "")

# Prevod casovych jednotek
library(datetime)
as.second(as.minute(7658))
as.hour(as.minute(7658))
as.day(as.minute(7658))
as.week(as.minute(7658))

# Prestupny rok
library(climtrends)
IsLeapYear(1620)

# Posledni den mesice
library(climtrends)
LastDayOfTheMonth(5,1945)


###
library(datetime)
dur <- duration(second = 3, minute = 1.5, hour = 2, day = 1.5, week = 0) 
as.numeric(dur, "hours") 
as.numeric(dur, "minutes")


span <- interval(ymd("2009-01-01"), ymd("2009-08-01")) 
as.duration(span) 

span <- interval(ymd_hms("2020-03-21 08:16:00"), ymd_hms("2020-03-21 12:03:00")) 
as.numeric(span, "hours") 
as.numeric(span, "minutes")


###

dtm <- strptime(c("1.1.2010 11:35"), format = "%d.%m.%Y %H:%M", tz = "CET")
dtm
str(dtm)

hrs <- function(u) {
  x <- u * 3600
  return(x)
}
dtm + hrs(3)

mns <- function(m) {
  x <- m * 60
  return(x)
}
dtm + hrs(3) + mns(10)

dtm + 3*3600 + 10*60   # 3 hod. 10 min.

dtm - 111  # 111 sec.

seq(as.Date("2000-01-01"), as.Date("2010-12-31"), by = "1 day") 


### formaty data

# default format is yyyy-mm-dd
mydates <- as.Date(c("2007-06-22", "2004-02-13"))
# convert dates to character data
strDates <- as.character(dates) 

dates <- c("05/27/84", "07/07/05")
betterDates <- as.Date(dates,format = "%m/%d/%y"); betterDates

dates <- c("May 27 1984", "July 7 2005")
betterDates <- as.Date(dates, format = "%B %d %Y"); betterDates



### Pokus zacal v 8 hodin 16 minut a skoncil ve 12 hodin 3 minuty. Kolik minut trval pokus?
library(lubridate)
date1 <- ymd_hms("2020-03-21 08:16:00") 
date2 <- ymd_hms("2020-03-21 12:03:00") 
dft1 = difftime(date1, date2, tz="Europe/Prague", units = "mins") # base
dft1
# units = c("auto", "secs", "mins", "hours", "days", "weeks")
abs(as.numeric(dft1)) # [min]
span <- interval(ymd_hms("2020-03-21 08:16:00"), ymd_hms("2020-03-21 12:03:00")) 
(dft2 = as.numeric(span, "minutes"))


### Pokus trval 7658 minut. Kolik je to dní, hodin a minut?
library(lubridate)
tt = 7658 # [min]
seconds_to_period(tt*60)
library(datetime)
seconds_to_period(as.second(as.minute(tt)))

# opacny vypocet
library(datetime)
dur <- duration(week = 0, day = 5, hours = 7, minutes = 38)
as.numeric(dur, "minutes")


### Bitva na Bílé hore probehla 8. 11. 1620. Kolik dní/mesícu/let uplynulo od bitvy?
today(tzone = "")
span <- interval(ymd("1620-11-08"), ymd(today(tzone = ""))) 
dur = as.duration(span) 
as.numeric(dur, "days")
as.numeric(dur, "months")
as.numeric(dur, "years")


### 
# "2007-9-29 23:00:00", "2007-9-30 23:00:00"
temp = c(0.280,0.208,-0.264,-0.480,-0.708,-0.714,-0.498,-0.216,0.126,0.574,1.042,1.086,1.820, 4.570,3.808,7.400,5.572,5.402,6.288,4.966,5.180,2.380,4.710, 5.366,4.766)
t1 = c("2007-9-29 23:00:00", "2007-9-30 00:00:00")
t2 = sapply(c(1:9), function(xx) paste("2007-9-30 0",xx,":00:00", sep=""))
t3 = sapply(c(10:23), function(xx) paste("2007-9-30 ",xx,":00:00", sep=""))
tt = c(t1,t2,t3)
df <- strptime(tt, format = "%Y-%m-%d %H:%M:%S")
plot(df,temp,xlab="time [h:m]",ylab="Temperature [°C]")


### V plasti duplikatoru zaznamenali pri jeho nahrivani zaznamenány teploty: v case 18 h 0 m  78.5 °C, v case 19 h 0 m 81.3 °C a v case 20 h 0 min 81.4 °C. Jaká byla teplota v case 18 h 30 m a 19 h 30 m?
temp = c(78.5, 81.3, 81.4) # [°C]
tt = c("18:00","19:00","20:00")
tt2 = c("18:30","19:30")
df <- strptime(tt, format = "%H:%M")
df2 <- strptime(tt2, format = "%H:%M")
plot(df,temp,xlab="time [h:m]",ylab="Temperature [°C]",type ="b")

tt = c(18.0, 19.0, 20.0)
tt2 = c(18.5, 19.5)
avt = approx(x=tt, y=temp, xout=tt2, method="linear")
res = cbind(avt$`x`,avt$y); colnames(res)=c("time [h]", "temperature [°C]")
res