#showing maps
------------------------------------------------------------------
Plot=function(H,main=NULL){
image(-H,zlim=c(-0.006,0),col=heat.colors(1000),main=main)}


PlotMAP=function(MAP,HD=0,main=NULL){
lev=c(0.1,0.25,0.5,1,1.5,2,3,4,5)
if(!HD){Plot(MAP,main); contour(MAP*1600,levels=lev,add=TRUE)}
if(HD){kernel_smoothMAP(MAP,11,lev,main)}}


kernel_smoothMAP=function(M,R,lev,main){
M=M/sum(M)
pencil=crea_tondi_SMAll(R)
n=ncol(M)
M2=matrix(0,5*n,5*n)
for(i in 1:(5*n)){for(j in 1:(5*n)){M2[i,j]=M[(i-1)%/%5+1,(j-1)%/%5+1]}}
M3=M2*0
for(i in R:(5*n-R+1)){for(j in R:(5*n-R+1)){
M3[i,j]=sum(M2[(i-R+1):(i+R-1),(j-R+1):(j+R-1)]*pencil)}}
Plot(M3,main=main); mM3=max(M3)
contour(M3*1600,add=T,levels=lev)}


crea_tondi_SMAll=function(R){ 
d=2*R-1
M=matrix(0,d,d)
for(i in 1:d){
for(j in 1:d){ 
x=i-R
y=j-R
M[i,j]=sqrt(x^2+y^2)}}
M=(sign(-M+R-0.499999)+1)/2
M=M/sum(M)
return(M)}


.................................................................
EXAMPLE:
PlotMAP(vkernel12); PlotMAP(vkernel12, HD=1)
PlotMAP(vkernel3)

PlotMAP(vents12) 
PlotMAP(vents3)

PlotMAP(vfaults,HD=1); PlotMAP(vfracts,HD=1); PlotMAP(maskTOT)







#creating syntetic experts and seed responses
------------------------------------------------------------------
EXAMPLE: Q=4 experts, k=10 seed, h=4 target
Q=4; k=10, h=4

syntSEED=round(runif(k,10,250)); syntSEED
..................................................................

syntEXP=function(Err,Inc,flag=1,i=1){
Bias=runif(k,-Err,Err)/100
Mexp=matrix(0,k,3)
Mexp[,2]=syntSEED*(1+Bias)
Mexp[,1]=Mexp[,2]*(1-Inc/100)
Mexp[,3]=Mexp[,2]*(1+Inc/100)
if(flag){plot(syntSEED,seq(1,k),pch=4,xlim=c(0,max(Mexp)),main=c('Expert',i),xlab='Seed Responses', ylab='Seed Questions' )
points(Mexp[,2],seq(1,10),col='red',pch=16)
points(Mexp[,1],seq(1,10),pch=1)
points(Mexp[,3],seq(1,10),pch=1)
for(i in 1:k){
lines(seq(1,1000)/1000*Mexp[,3],rep(i,1000))
lines(seq(1,1000)/1000*(Mexp[i,3]-Mexp[i,1])+Mexp[i,1],rep(i,1000),lwd=2)}}
return(round(Mexp))}


createEXP=function(vErr,vInc,Q){
Mexp=NULL
for(i in 1:Q){
print(c('Expert',i)); flush.console()
print(c(vErr[i], vInc[i])); flush.console()}
for(i in Q:1){
dev.new(); Mexp[[i]]=syntEXP(vErr[i],vInc[i],1,i)}
return(Mexp)}

.....................................................................
EXAMPLE: 
vErr=c(30,30,50,50); vInc=c(50,30,30,10)

MexpS=createEXP(vErr,vInc,Q)







#calculating ERF weight from the seed responses
--------------------------------------------------------------------------------
ERFweights=function(Mexp,syntSEED){
pERF=numeric(Q)
p_single=numeric(k)
for(i in 1:Q){
for(j in 1:k){
p_single[j]=ERFsingle(syntSEED[j],Mexp[[i]][j,1],Mexp[[i]][j,2],Mexp[[i]][j,3])}
pERF[i]=mean(p_single)}
return(100*pERF/sum(pERF))}


ERFsingle=function(x,a,b,c){
A=min(max(0.95*x,a),c)
B=max(min(1.05*x,c),a)
p=((A-c)^2-(B-c)^2)/((c-a)*(c-b))
if(A<b){
p=1-(B-c)^2/((c-a)*(c-b))-(A-a)^2/((b-a)*(c-a))
if(B<b){
p=((B-a)^2-(A-a)^2)/((b-a)*(c-a))}}
p}


........................................................................
EXAMPLE:
ERFweights(MexpS,syntSEED)
#do the same with EXCALIBUR and Classical Method








#defining target responses and plotting them
--------------------------------------------------------------------------
EXAMPLE:
#VAR-UNIF, VENTS-STRUCT, EPOI/II-EPOIII, FAULTS-FRACT
# INC 10% e 25%

EXP=NULL
EXP[[1]]=matrix(0,3,4)
EXP[[1]][,1]=c(70,90,100)
EXP[[1]][,2]=c(80,90,95)
EXP[[1]][,3]=c(25,50,80)
EXP[[1]][,4]=c(20,50,75)

EXP[[2]]=matrix(0,3,4)
EXP[[2]][,1]=c(75,80,95)
EXP[[2]][,2]=c(70,80,90)
EXP[[2]][,3]=c(5,25,40)
EXP[[2]][,4]=c(10,25,45)

EXP[[3]]=matrix(0,3,4)
EXP[[3]][,1]=c(75,85,95)
EXP[[3]][,2]=c(65,75,90)
EXP[[3]][,3]=c(65,80,90)
EXP[[3]][,4]=c(55,75,90)

EXP[[4]]=matrix(0,3,4)
EXP[[4]][,1]=c(75,80,95)
EXP[[4]][,2]=c(60,75,90)
EXP[[4]][,3]=c(10,30,40)
EXP[[4]][,4]=c(60,75,85)
..............................................................................

plotEXPt=function(EXP,i=0){if(i){EXP=EXP[[i]]}; if(!i){i='DM'}
plot(EXP[2,],seq(h,1),pch=16,xlim=c(0,100),ylim=c(0.5,h),col='red',main=c('Expert',i),xlab='Target Responses',ylab='Target Questions')
points(EXP[1,],seq(h,1),pch=1)
points(EXP[3,],seq(h,1),pch=1)
for(i in 1:h){
lines(seq(1,1000)/1000*100,rep(i,1000))
lines(seq(1,1000)/1000*(EXP[3,h+1-i]-EXP[1,h+1-i])+EXP[1,h+1-i],rep(i,1000),lwd=2)}}

...............................................................................
EXAMPLE: 
plotEXPt(EXP,1)










#function generating MAP coefficients and producing maps
-----------------------------------------------------------------------------------
createMAP=function(EXPv,flag=1,Ckernel=0,Csubzones=0,HD=0,MC=0){if(MC==0){EXPv=EXPv[2,]/100}
P=numeric(5)
P[1]=1-EXPv[1]
P[2]=EXPv[1]*EXPv[2]*EXPv[3]
P[3]=EXPv[1]*EXPv[2]*(1-EXPv[3])
P[4]=EXPv[1]*(1-EXPv[2])*EXPv[4]
P[5]=EXPv[1]*(1-EXPv[2])*(1-EXPv[4])
if(!flag){return(c(P[2:5],P[1])*100)}
if(flag){
if(!Ckernel){MAP=maskTOT*P[1]+vents12*P[2]+vents3*P[3]+vfaults*P[4]+vfracts*P[5]}
if(Ckernel){MAP=maskTOT*P[1]+vkernel12*P[2]+vkernel3*P[3]+vfaults*P[4]+vfracts*P[5]}}
if(Csubzones>0){Prob=calcSubZones(MAP,Csubzones);return(Prob*100)}
if(flag==2){PlotMAP(MAP,HD)}
if(flag==1){return(MAP)}}

.....................................................................................
EXAMPLE
createMAP(EXP[[1]],flag=0)
createMAP(EXP[[1]],flag=2)
createMAP(EXP[[1]],flag=2,Csubzones=8)
createMAP(EXP[[1]],flag=2,HD=1)
createMAP(EXP[[1]],flag=2,Ckernel=1)








#producing confidence distributions for coefficients and probability maps of a single expert
--------------------------------------------------------------------------------------
MonteCarloCoeff=function(EXP,N){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); coeff=matrix(0,N,5)
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
coeff[i,]=createMAP(EXPv,0,MC=1)}
coln=c('EpoI_II','EpoIII','Faults','Fracts','Homog'); rown=c('5%ile','mean','95%ile')
a=matrix(0,3,5,dimnames=list(rown,coln))
for(i in 1:5){a[,i]=createQuantiles(coeff[,i])}
return(a)}


MonteCarloMap=function(EXP,N,Ckernel=0,Csubzones=0,HD=0){
EXP=EXP/100
coef=matrix(0,2,h); EXPv=numeric(h); MAPS=NULL; 
if(Csubzones){Prob=numeric(N)}
for(i in 1:h){
coef[,i]=NewRap(EXP[1,i],EXP[2,i],EXP[3,i])
for(j in 1:2){coef[j,i]=min(max(coef[j,i],0),1)}}
for(i in 1:N){for(j in 1:h){
EXPv[j]=rtrianINN(coef[1,j],EXP[2,j],coef[2,j])}
if(i%%100==0){print(i); flush.console()}
MAPS[[i]]=createMAP(EXPv,1,Ckernel,MC=1)
if(Csubzones>0){Prob[i]=calcSubZones(MAPS[[i]],Csubzones)}}
if(Csubzones>0){P=createQuantiles(Prob); return(P*100)}
MAPS=seekMAPS(MAPS,N)
PlotMAP(MAPS[[1]],HD,main='5%ile') 
dev.new(); PlotMAP(MAPS[[3]],HD,main='95%ile')  
dev.new(); PlotMAP(MAPS[[2]],HD,main='mean') 
return(MAPS)}


seekMAPS=function(MAPS,N){ 
n=length(MAPS[[1]])
QMAP=matrix(0,n,3)
MAPmat=matrix(0,n,N)
for(j in 1:N){
MAPmat[,j]=matrix(MAPS[[j]],n,1)}
for(i in 1:n){QMAP[i,]=createQuantiles(MAPmat[i,])}
QMAPf=NULL
for(i in 1:3){
QMAPf[[i]]=matrix(QMAP[,i],100,100)}
return(QMAPf)}

createQuantiles=function(vect,flag=0){
a=numeric(3)
a[1]=quantile(vect,0.05)
if(!flag){a[2]=mean(vect)} 
if(flag){a[2]=quantile(vect,0.5)}
a[3]=quantile(vect,0.95)
return(a)}



#measuring probability inside sub-sectors
calcSubZones=function(M,z){
Prob=0
if(z>0){Prob=sum(M*Z_partit[[z]])}
return(Prob)}



#testing with triangular distributions
rtrian=function(a,b,c){
if(a==c){R=b}
else {R=rtrianINN(NewRap(a,b,c)[1],b,NewRap(a,b,c)[2])}}

rtrianINN=function(a,b,c){
u=runif(1)
if(u<((b-a)/(c-a))){x=sqrt(u*(b-a)*(c-a))+a}
else{x=c-sqrt((1-u)*(c-a)*(c-b))}
x}

NewRap=function(a,b,c){
x0=a-(c-a)/6
y0=c+(c-a)/6
x=c(x0,y0)
for(i in 1:5){
x=x-crossprod(InverseJac(x[1],x[2],a,b,c),FunRap(x[1],x[2],a,b,c))}
x}

InverseJac=function(x,y,a,b,c){
A=1.9*x+0.05*y-2*a+0.05*b
B=0.05*(x-b)
C=0.05*(y-b)
D=1.9*y+0.05*x-2*c+0.05*b
M=matrix(0,2,2)
M[1,1]=D
M[2,2]=A
M[1,2]=-B
M[2,1]=-C
M=M/(A*D-B*C)
M}

FunRap=function(x,y,a,b,c){
A1=(a-x)^2-0.05*(y-x)*(b-x)
A2=(y-c)^2-0.05*(y-x)*(y-b)
c(A1,A2)}

...................................................................................
EXAMPLE
MonteCarloCoeff(EXP[[1]],10^4)
outp=MonteCarloMap(EXP[[1]],10^4)
outp=MonteCarloMap(EXP[[1]],10^4,1)
MonteCarloMap(EXP[[1]],10^4,Csubzones=8)








#define the ERF decision maker (quantile pooling and triangular dist.)
--------------------------------------------------------------------------------------
DM_ERF=function(EXP,Mexp,syntSEED){
pERF=ERFweights(Mexp,syntSEED)/100
DM=EXP[[1]]*0
for(i in 1:h){
DM=DM+EXP[[i]]*pERF[i]}
return(DM)}

............................................................................................
EXAMPLE:
DMerf=DM_ERF(EXP,Mexp,syntSEED)

plotEXPt(DMerf)

EXAMPLE
MonteCarloCoeff(DMerf,10^4)
outp=MonteCarloMap(DMerf,10^4)
outp=MonteCarloMap(DMerf,10^4,HD=1)
MonteCarloMap(DMerf,10^4,Csubzones=8)






#testing with max entropy samples 
--------------------------------------------------------------------------------------------
max_entropy=function(a,b,c,rA,rB){
x=runif(1)
y=runif(1,b,c)
if(x>0.95){y=runif(1,c,rB)}
if(x<0.5){y=runif(1,a,b)
if(x<0.05){y=runif(1,rA,a)}}
return(y)}







#define EW decision maker, or assuming arbitrary weights (linear pooling and max entropy dist.)
--------------------------------------------------------------------------------------------
#10% overshooting, 10^5 MC samples

DM_linear=function(EXP,P=rep(1/Q,Q),N=10^5){
coln=c('5%ile','mean','95%ile')
quan=matrix(0,h,3,dimnames=list(NULL,coln))
vexp=matrix(0,Q,3)
for(i in 1:h){
for(j in 1:Q){vexp[j,]=EXP[[j]][,i]}
rA=min(vexp[,1]); rB=max(vexp[,3])
R=rB-rA; rA=rA-R/10; rB=rB+R/10
test=numeric(N)
for(j in 1:N){estr=ceiling(runif(1)*Q); test[j]=max_entropy(vexp[estr,1],vexp[estr,2],vexp[estr,3],rA,rB)}
quan[i,]=createQuantiles(test,1)
print(i); flush.console()}
return(quan)}

.........................................................
EXAMPLE
DM_linear(EXP)
DM_linear(EXP,c(0.5,0.5,0,0))
#DM_linear with Classical weights