# Introduction to Geochemical modeling
# p.85
m=5;vec=sample(1:10,m,replace = T)
A=t(t(vec))%*%t(vec)+diag(vec)
x0=t(t(c(rep(1,m))))
t=seq(0,10,0.01)
U=eigen(A)$vectors
lambda=diag(eigen(A)$values)
inv_U=solve(U)
U%*%lambda%*%inv_U
t_vec=seq(0,5,0.01)
X=array(0,dim=c(ncol(A),length(t_vec)))
for(t in 1:length(t_vec)){
x=array(0,dim=c(ncol(A),1))
for(j in 1:ncol(A)){
x=x+as.numeric(exp(eigen(A)$values[j]*t_vec[t]))*as.numeric(t(t(U[,j]))%*%t(inv_U[j,])%*%x0)
}
X[,t]=x
}
# p.95
m=5;vec=sample(1:10,m,replace = T)
A=t(t(vec))%*%t(vec)+diag(vec)
x0=t(t(c(rep(1,m))));b=rep(1,m)
U=eigen(A)$vectors
lambda=diag(eigen(A)$values)
inv_U=solve(U)
U%*%lambda%*%inv_U
t_vec=seq(0,5,0.01)
for(t in 1:length(t_vec)){
x=array(0,dim=c(dim(A)))
for(j in 1:ncol(A)){
x=x+as.numeric(exp(eigen(A)$values[j]*t_vec[t]))*(t(t(U[,j]))%*%t(inv_U[j,]))
}
exp_At=x
X[,t]=-solve(A)%*%b+exp_At%*%(x0+solve(A)%*%b)
}
# Introduction to Geochemical Modeling p.150
library(dplyr)
k=0.1;m=5
E=100;N=50
alpha=0.1;k=1-alpha
cost=function(ep_i,a,b){
ni=exp(a*ep_i/k);N=sum(ni)
S=k*sum((ni/N)*log(ni/N))
S_star=S+a*(sum(ep_i*ni)-E)^2+b*(sum(ni)-N)^2
print(S_star)
}
h=0.01
epsilon_i=sample(1:10,m,replace=T)/10
ite=1000;eta=0.0001
for(l in 1:ite){
vec=epsilon_i
for(j in 1:m){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
epsilon_i[j]=epsilon_i[j]-eta*(cost(vec_sub,alpha,beta)-cost(vec,alpha,beta))/h
}
ni=exp(alpha*epsilon_i/k);N=sum(ni)
S=k*sum((ni/N)*log(ni/N));
print(cost(epsilon_i,alpha,beta))
}
# p.308
m=5
x0=t(t(sample(1:100,m,replace=T)))
x=t(t(sample(1:100,m,replace=T)))
vec1=sample(1:10,m,replace=T);vec2=sample(1:10,m,replace=T)
g1=function(y){
return(sum(((y-vec1)^2)))
}
g2=function(y){
return(sum(((y-vec2)^2)))
}
ite=100
h=0.01
for(l in 1:ite){
f1=c();f2=c();x_vec=x
for(i in 1:m){
x_vec_sub=x_vec;x_vec_sub[i]=x_vec_sub[i]+h
dg1=(g1(x_vec_sub)-g1(x_vec))/h
dg2=(g2(x_vec_sub)-g2(x_vec))/h
f1=c(f1,dg1)
f2=c(f2,dg2)
}
f_mat=cbind(f1,f2)
x=x0-f_mat%*%solve(t(f_mat)%*%f_mat)%*%(t(f_mat)%*%(x-x0)-t(t(c(g1(x),g2(x)))))
lambda=solve(t(f_mat)%*%f_mat)%*%t(f_mat)%*%(x-x0)
cost_c2=as.numeric(t(x-x0)%*%(x-x0))-sum(lambda*t(t(c(g1(x),g2(x)))))
print(cost_c2)
}
# p.337
B=t(matrix(c(1,2,0,0,0,0,2,0,0,1,0,2,1,1,0,0,0,0,0,2,0,1,0,0),nrow=4))
P=diag(1,nrow(B))-B%*%solve(t(B)%*%B)%*%t(B)
n=t(t(rep(10,nrow(B))))
C=sample(1:10,nrow(B),replace=T)/10
ite=10000;alpha=0.001
for(l in 1:ite){
n_pre=n
dG_pre=c(C+log(n+1)-rep(1,nrow(B))*as.numeric(log(t(rep(1,nrow(B)))%*%n)))
#n=n+alpha*(P%*%dG)
mu=C+log(n/sum(n));
mu_pre=mu
if(sum(n+alpha*mu>0)==nrow(B)){
n=n+alpha*mu
dG=c(C+log(n+1)-rep(1,nrow(B))*as.numeric(log(t(rep(1,nrow(B)))%*%n)))
}
mu=C+log(n/sum(n));
f=as.numeric(t(P%*%mu_pre)%*%(P%*%mu))
G=sum(c(n)*(C+log(c(n+1))-log(sum(n))))
#print(G)
print(f)
}