text:
using DifferentialEquations, Plots, LaTeXStrings, LinearAlgebra, DiffEqDevTools, SparseArrays
gr()
default(titlefont = 8, legendfontsize = 4, guidefont = 5, tickfont = 5, markersize =2.5)#,  markerstrokewidth = 0) #thickness_scaling=0.8)

setups = [Dict(:alg=>Rodas4()),Dict(:alg=>Rodas4P()),Dict(:alg=>Rodas4P2()),Dict(:alg=>Rodas5()),Dict(:alg=>Rodas5P())]
setups_1 = [Dict(:alg=>Rodas4(;autodiff=false)),Dict(:alg=>Rodas4P(;autodiff=false)),Dict(:alg=>Rodas4P2(;autodiff=false)),Dict(:alg=>Rodas5(;autodiff=false)),Dict(:alg=>Rodas5P(;autodiff=false))]

#---- Problem definitions --------------------------------------------------------------------

title_1 = "Parabol"
nx = 250; dx = 1.0/nx; n = nx-1; x = zeros(n);
M_1 = sparse(1.0I,n,n); J = spzeros(n,n)
#M_1 = Matrix(1.0I, n, n); J = zeros(n,n)
numrun_1 = 10;
for i=1:n
  x[i] = i*dx
end
b1  = zeros(n); b1[1]=3/(4*dx^2); b1[n]=3/(4*dx^2);
for i=1:n
  J[i,i] = -2/(dx^2);
  if i>1 J[i,i-1] = 1/(dx^2); end
  if i<n J[i,i+1] = 1/(dx^2); end
end
p_1 = x, b1, J
function ode_1(du,u,p,t)
  x, b1, J = p
  du[:] = J*u .+ b1/(1+t) .- (x .+ 1/2).*(3/2 .- x)./(t+1).^2 .- 2/(1+t);
  nothing
end
function jac_1(J,u,p,t)
#  x, b1, J = p
  J[:] = p[3]
  nothing
end
function sol_1(t,p)
  x, b1, J = p
  u = (x .+ 1/2).*(3/2 .- x)/(1+t);
end
y0_1 = sol_1(0.0,p_1); tspan_1 = [0.0, 0.1];

title_2 = "Hyperbol"
nx = 250; dx = 1.0/nx; n = nx; x = zeros(n);
M_2 = sparse(1.0I,n,n); J = spzeros(n,n)
#M_2 = Matrix(1.0I, n, n); J = zeros(n,n)
numrun_2 = 10;
for i=1:n
  x[i] = i*dx
end
b1  = zeros(n); b1[1]=-1/dx;
for i=1:n
  J[i,i] = 1/dx;
  if i>1 J[i,i-1] = -1/dx; end
end
p_2 = x, b1, J
function ode_2(du,u,p,t)
  x, b1, J = p
  du[:] = -J*u .- b1/(1+t) .+ (t .- x)./(t .+ 1).^2;
  nothing
end
function jac_2(J,u,p,t)
  J[:] = -p[3]
  nothing
end
function sol_2(t,p)
  x, b1, J = p
  u = (1 .+ x)./(1 .+ t);
end
y0_2 = sol_2(0.0,p_2); tspan_2 = [0.0, 1];

title_3 = "Pendulum Index-1"; title_4 = "Pendulum Index-2"
M_3 = Matrix(1.0I, 5, 5); M_3[5,5] = 0.0; numrun_3 = 20; numrun_4 = 5; p_3 = 1; p_4 = 2;
y0_3 = [2.0, 0.0, 0.0, 0.0, 0.0]; tspan_3 = [0.0, 10.0];
function ode_3(du,u,p,t)
  du[1] = u[2]
  du[2] = u[1]*u[5]
  du[3] = u[4]
  du[4] = -9.81 + u[3]*u[5]
  if p==1
    du[5] = u[2]^2 + u[4]^2 + u[5]*(u[1]^2 + u[3]^2) -9.81*u[3] #-- Index-1
  else
    du[5] = u[1]*u[2] + u[3]*u[4]; #-- Index-2
  end
end
function jac_3(J,y,p,t)
  if p==1
    J[:] = [0 1 0 0 0; y[5] 0 0 0 y[1]; 0 0 0 1 0; 0 0 y[5] 0 y[3]; 2*y[5]*y[1] 2*y[2] 2*y[5]*y[3]-9.81 2*y[4] y[1]^2+y[3]^2] #-- Index-1
  else
    J[:] = [0 1 0 0 0; y[5] 0 0 0 y[1]; 0 0 0 1 0; 0 0 y[5] 0 y[3]; y[2] y[1] y[4] y[3] 0] #-- Index-2
  end
  nothing
end

#-- transistor amplifier from IVP testset
title_5 = "Transistor amplifier"; numrun_5 = 10;
ub=6.0; uf=0.026; alpha=0.99; beta=1.0e-6; r0=1000; r1=9000; r2=9000; r3=9000; r4=9000; r5=9000; r6=9000; r7=9000; r8=9000; r9=9000;
p_5 = ub,uf,alpha,beta,r0,r1,r2,r3,r4,r5,r6,r7,r8,r9
tspan_5 = [0, 0.2];
y0_5 = [0.0, ub/(r2/r1+1.0),  ub/(r2/r1+1.0), ub, ub/(r6/r5+1.0), ub/(r6/r5+1.0), ub, 0.0]
c1=1.0e-6; c2=2.0e-6; c3=3.0e-6; c4=4.0e-6; c5=5.0e-6;
M_5 = zeros(8,8);
M_5[1,1]=-c1; M_5[1,2]=c1; M_5[2,1]=c1; M_5[2,2]=-c1; M_5[3,3]=-c2;
M_5[4,4]=-c3; M_5[4,5]=c3; M_5[5,4]=c3; M_5[5,5]=-c3; M_5[6,6]=-c4;
M_5[7,7]=-c5; M_5[7,8]=c5; M_5[8,7]=c5; M_5[8,8]=-c5;

function ode_5(du,u,p,t)
  ub,uf,alpha,beta,r0,r1,r2,r3,r4,r5,r6,r7,r8,r9 = p
  uet = 0.1*sin(200*pi*t)
  arg = min((u[2]-u[3])/uf,100.0); fac1  = beta*(exp(arg)-1.0)
  arg = min((u[5]-u[6])/uf,100.0); fac2  = beta*(exp(arg)-1.0)
#  fac1  = beta*(exp((u[2]-u[3])/uf)-1.0)
#  fac2  = beta*(exp((u[5]-u[6])/uf)-1.0)
  du[1] = (u[1]-uet)/r0
  du[2] = u[2]/r1 + (u[2]-ub)/r2 + (1.0-alpha)*fac1
  du[3] = u[3]/r3 - fac1
  du[4] = (u[4]-ub)/r4 + alpha*fac1
  du[5] = u[5]/r5 + (u[5]-ub)/r6 + (1.0-alpha)*fac2
  du[6] = u[6]/r7 - fac2
  du[7] = (u[7]-ub)/r8 + alpha*fac2
  du[8] = u[8]/r9
  nothing
end
function jac_5(J,u,p,t)
  ub,uf,alpha,beta,r0,r1,r2,r3,r4,r5,r6,r7,r8,r9 = p
  arg = (u[2]-u[3])/uf; fac1  = beta*exp(arg)/uf
  arg = (u[5]-u[6])/uf; fac2  = beta*exp(arg)/uf
  J .= 0.0
  J[1,1] = 1.0/r0
  J[2,2] = 1.0/r1 + 1.0/r2 + (1.0-alpha)*fac1
  J[2,3] = -(1.0-alpha)*fac1
  J[3,2] = -fac1
  J[3,3] = 1.0/r3 + fac1
  J[4,2] = alpha*fac1
  J[4,3] = -alpha*fac1
  J[4,4] = 1.0/r4
  J[5,5] = 1.0/r5 + 1.0/r6 + (1.0-alpha)*fac2
  J[5,6] = -(1.0-alpha)*fac2
  J[6,5] = -fac2
  J[6,6] = 1.0/r7 + fac2
  J[7,5] = alpha*fac2
  J[7,6] = -alpha*fac2
  J[7,7] = 1.0/r8
  J[8,8] = 1.0/r9
  nothing
end

title_6 = "PV-Network"; numrun_6 = 20
c1 = -3.1037; c2 = 1.0015; c3 = 0.0032; c4 = 1.3984e-09; c5 = 0.4303; c6 = 1.5*0.9562;
R0 = 0.2; R1 = 0.5; C = 4000.0; q_max = 36000.0;
p_6 = c1,c2,c3,c4,c5,c6,R0,R1,C,q_max
y0_6 = [0.0, 11.856598910310167, 0.0, 2.9409008015416687, -2.940900801550821, 0.0, 9000.0]
M_6 = zeros(7,7); M_6[6,6] = 1.0; M_6[7,7] = 1.0;
tspan_6 = [0.0, 36000.0];

function ode_6(dy, y, p, t)
  c1,c2,c3,c4,c5,c6,R0,R1,C,q_max = p
  U0 = y[1]; U1 = y[2]; iV = y[3]; iPV = y[4]; iB = y[5]; uB = y[6]; qB = y[7];
  dy[1] = U0;
  dy[2] = iB + iPV - iV;
  dy[3] = verbraucher(t) - iV*(U1-U0);
  dy[4] = c1 + c2*iPV + c3*(U1-U0) + c4*(exp(c5*iPV+c6*(U1-U0))-1);
  dy[5] = U1-U0 - (ruhespannung(qB/q_max) - uB - R0*iB);
  dy[6] = iB/C - uB/(R1*C);
  dy[7] = -iB;
  nothing
end
function verbraucher(t)
  ts = range(3600.0,36000.0,step = 3600.0)
  return 50.0*einaus(t,ts,60.0)
#  return 50
end
function ruhespannung(soc)
  pp = [6.8072,-10.5555,6.2199,10.2668]; #-- Polynomkoeffizienten
  return ((pp[1]*soc+pp[2])*soc+pp[3])*soc+pp[4]; #-- Polynomauswertung
end
function fstep(t,t0,dauer)
#-- atanh(0.999) = 3.8002, atanh(0.99) = 2.647
   s = 2.647/dauer
   return (tanh(s*(t-t0))+1.0)/2.0
end
function einaus(t,ts,schaltzeit)
   y = 0.0;
   for i=1:2:length(ts)
       y = y + fstep(t,ts[i],schaltzeit);
   end
   for i=2:2:length(ts)
       y = y - fstep(t,ts[i],schaltzeit);
   end
   return y
end

title_7 ="Water tube problem"; numrun_7 = 10;
nu = 1.31e-6; g  = 9.8; rho = 1.0e3; rcrit = 2.3e3; len = 1.0e3; k = 2.0e-4; d = 1.0; b = 2.0e2;
a = pi*d^2/4.0; mu = nu*rho; c = b/(rho*g); v = rho*len/a;
nnodes = 13; nedges = 18;
ANN=zeros(nnodes,nnodes)
ANN[ 1, 2] = 1; ANN[ 2, 3] = 1; ANN[ 2, 6] = 1; ANN[ 3, 4] = 1; ANN[ 3, 5] = 1; ANN[ 4, 5] = 1;
ANN[ 5,10] = 1; ANN[ 6, 5] = 1; ANN[ 7, 4] = 1; ANN[ 7, 8] = 1; ANN[ 8, 5] = 1; ANN[ 8,10] = 1;
ANN[ 9, 8] = 1; ANN[11, 9] = 1; ANN[11,12] = 1; ANN[12, 7] = 1; ANN[12, 8] = 1; ANN[13,11] = 1;
ANE = zeros(nnodes,nedges);
iann = [1,2,2,3,3,4,5,6,7,7,8,8,9,11,11,12,12,13]
jann = [2,3,6,4,5,5,10,5,4,8,5,10,8,9,12,7,8,11]
for i=1:nnodes
  for j = 1:nedges
    if iann[j]==i ANE[i,j] = -1; end
    if jann[j]==i ANE[i,j] = 1; end
  end
end
dein = zeros(nnodes); rghres=zeros(nedges); fdba=zeros(nedges);
p_7 = nu, rho, rcrit, len, k, d, a, mu, c, v, nnodes, nedges, iann, jann, ANE, dein, rghres, fdba
tspan_7 = [0.0,17*3600.0]; neq = 49; y0_7 = zeros(neq);
y0_7[19:36] .= 0.47519404529185289807e-1; y0_7[37:49] .= 109800.0;
M_7 = zeros(neq,neq);
for i=1:18
  M_7[i,i] = 1;
end
i5 = 2*nedges+5; i8 = 2*nedges+8; M_7[i5,i5] = 1; M_7[i8,i8] = 1;

function ode_7(du,u,p,t)
  nu, rho, rcrit, len, k, d, a, mu, c, v, nnodes, nedges, iann, jann, ANE, dein, rghres, fdba = p
  phi = u[1:nedges]; lambda = u[nedges+1:2*nedges]; pp = u[2*nedges+1:end]
  that = t/3600.0;
  dein[1] = -sin(exp(-that)-1.0)*exp(-that)/(3600.0*200);
  dein[13] = -sin(exp(-that)-1.0)*exp(-that)/(3600.0*80);
  dein[10] = -(12*that^3-3*92*that^2+2*720*that)/(1.0e6*3600);
  for n=1:nedges
    pL = pp[iann[n]]; pR = pp[jann[n]];
    rtla = sqrt(lambda[n]); r = abs(phi[n]*d/(nu*a));
    rghres[n] = 1.0/rtla - 1.74 + 2.0*log10(2.0*k/d + 18.7/(max(r,rcrit)*rtla));
    fdba1 = pL - pR - lambda[n]*rho*len*phi[n]^2/(a^2*d);
    fdba0 = pL - pR - 32.0*mu*len*phi[n]/(a*d^2);
    fac0 = 0.9; fac1 = 1.1; alpha = (r/rcrit-fac0)/(fac1-fac0); #-- continous transition
    alpha = max(alpha,0.0); alpha = min(alpha,1.0);
    fdba[n] = ( (1.0 - alpha)*fdba0 + alpha*fdba1 )/v;
#    if r > rcrit
#      fdba[n] = fdba1/v
#    else
#      fdba[n] = fdba0/v
#    end
  end
#--
  netflow = dein .+ ANE*fdba;
  AP = (ANE*phi)/c;  #-- Storage nodes
  netflow[5] = AP[5]; netflow[8] = AP[8];
#--
  du[:]  = [fdba; rghres; netflow];
  nothing
end

title_8 = "Pollution"; numrun_8 = 20;
k1=.35e0; k2=.266e2; k3=.123e5; k4=.86e-3; k5=.82e-3; k6=.15e5; k7=.13e-3; k8=.24e5; k9=.165e5; k10=.9e4; k11=.22e-1; k12=.12e5; k13=.188e1
k14=.163e5; k15=.48e7; k16=.35e-3; k17=.175e-1; k18=.1e9; k19=.444e12; k20=.124e4; k21=.21e1; k22=.578e1; k23=.474e-1; k24=.178e4; k25=.312e1

p_8 = k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25

y0_8 = zeros(20)
y0_8[2]  = 0.2; y0_8[4]  = 0.04; y0_8[7]  = 0.1; y0_8[8]  = 0.3; y0_8[9]  = 0.01; y0_8[17] = 0.007
tspan_8 = [0.0,60.0]; M_8 = Matrix(1.0I, 20, 20);

function ode_8(dy,y,p,t)
 k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25 = p
 r1  = k1 *y[1]
 r2  = k2 *y[2]*y[4]
 r3  = k3 *y[5]*y[2]
 r4  = k4 *y[7]
 r5  = k5 *y[7]
 r6  = k6 *y[7]*y[6]
 r7  = k7 *y[9]
 r8  = k8 *y[9]*y[6]
 r9  = k9 *y[11]*y[2]
 r10 = k10*y[11]*y[1]
 r11 = k11*y[13]
 r12 = k12*y[10]*y[2]
 r13 = k13*y[14]
 r14 = k14*y[1]*y[6]
 r15 = k15*y[3]
 r16 = k16*y[4]
 r17 = k17*y[4]
 r18 = k18*y[16]
 r19 = k19*y[16]
 r20 = k20*y[17]*y[6]
 r21 = k21*y[19]
 r22 = k22*y[19]
 r23 = k23*y[1]*y[4]
 r24 = k24*y[19]*y[1]
 r25 = k25*y[20]

 dy[1]  = -r1-r10-r14-r23-r24+r2+r3+r9+r11+r12+r22+r25
 dy[2]  = -r2-r3-r9-r12+r1+r21
 dy[3]  = -r15+r1+r17+r19+r22
 dy[4]  = -r2-r16-r17-r23+r15
 dy[5]  = -r3+r4+r4+r6+r7+r13+r20
 dy[6]  = -r6-r8-r14-r20+r3+r18+r18
 dy[7]  = -r4-r5-r6+r13
 dy[8]  = r4+r5+r6+r7
 dy[9]  = -r7-r8
 dy[10] = -r12+r7+r9
 dy[11] = -r9-r10+r8+r11
 dy[12] = r9
 dy[13] = -r11+r10
 dy[14] = -r13+r12
 dy[15] = r14
 dy[16] = -r18-r19+r16
 dy[17] = -r20
 dy[18] = r20
 dy[19] = -r21-r22-r24+r23+r25
 dy[20] = -r25+r24
end

function jac_8(J,y,p,t)
  k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25 = p
      J .= 0.0
      J[1,1]   = -k1-k10*y[11]-k14*y[6]-k23*y[4]-k24*y[19]
      J[1,11]  = -k10*y[1]+k9*y[2]
      J[1,6]   = -k14*y[1]
      J[1,4]   = -k23*y[1]+k2*y[2]
      J[1,19]  = -k24*y[1]+k22
      J[1,2]   = k2*y[4]+k9*y[11]+k3*y[5]+k12*y[10]
      J[1,13]  = k11
      J[1,20]  = k25
      J[1,5]   = k3*y[2]
      J[1,10]  = k12*y[2]

      J[2,4]   = -k2*y[2]
      J[2,5]   = -k3*y[2]
      J[2,11]  = -k9*y[2]
      J[2,10]  = -k12*y[2]
      J[2,19]  = k21
      J[2,1]   = k1
      J[2,2]   = -k2*y[4]-k3*y[5]-k9*y[11]-k12*y[10]

      J[3,1]   = k1
      J[3,4]   = k17
      J[3,16]  = k19
      J[3,19]  = k22
      J[3,3]   = -k15

      J[4,4]   = -k2*y[2]-k16-k17-k23*y[1]
      J[4,2]   = -k2*y[4]
      J[4,1]   = -k23*y[4]
      J[4,3]   = k15

      J[5,5]   = -k3*y[2]
      J[5,2]   = -k3*y[5]
      J[5,7]   = 2k4+k6*y[6]
      J[5,6]   = k6*y[7]+k20*y[17]
      J[5,9]   = k7
      J[5,14]  = k13
      J[5,17]  = k20*y[6]

      J[6,6]   = -k6*y[7]-k8*y[9]-k14*y[1]-k20*y[17]
      J[6,7]   = -k6*y[6]
      J[6,9]   = -k8*y[6]
      J[6,1]   = -k14*y[6]
      J[6,17]  = -k20*y[6]
      J[6,2]   = k3*y[5]
      J[6,5]   = k3*y[2]
      J[6,16]  = 2k18

      J[7,7]   = -k4-k5-k6*y[6]
      J[7,6]   = -k6*y[7]
      J[7,14]  = k13

      J[8,7]   = k4+k5+k6*y[6]
      J[8,6]   = k6*y[7]
      J[8,9]   = k7

      J[9,9]   = -k7-k8*y[6]
      J[9,6]   = -k8*y[9]

      J[10,10] = -k12*y[2]
      J[10,2]  = -k12*y[10]+k9*y[11]
      J[10,9]  = k7
      J[10,11] = k9*y[2]

      J[11,11] = -k9*y[2]-k10*y[1]
      J[11,2]  = -k9*y[11]
      J[11,1]  = -k10*y[11]
      J[11,9]  = k8*y[6]
      J[11,6]  = k8*y[9]
      J[11,13] = k11

      J[12,11] = k9*y[2]
      J[12,2]  = k9*y[11]

      J[13,13] = -k11
      J[13,11] = k10*y[1]
      J[13,1]  = k10*y[11]

      J[14,14] = -k13
      J[14,10] = k12*y[2]
      J[14,2]  = k12*y[10]

      J[15,1]  = k14*y[6]
      J[15,6]  = k14*y[1]

      J[16,16] = -k18-k19
      J[16,4]  = k16

      J[17,17] = -k20*y[6]
      J[17,6]  = -k20*y[17]

      J[18,17] = k20*y[6]
      J[18,6]  = k20*y[17]

      J[19,19] = -k21-k22-k24*y[1]
      J[19,1]  = -k24*y[19]+k23*y[4]
      J[19,4]  = k23*y[1]
      J[19,20] = k25

      J[20,20] = -k25
      J[20,1]  = k24*y[19]
      J[20,19] = k24*y[1]

      return
end

#--- collect problems -----
title_s = [title_1,title_2,title_3,title_4,title_5,title_6,title_7,title_8]
ode_s = [ode_1,ode_2,ode_3,ode_3,ode_5,ode_6,ode_7,ode_8]
y0_s = [y0_1,y0_2,y0_3,y0_3,y0_5,y0_6,y0_7,y0_8]
tspan_s = [tspan_1,tspan_2,tspan_3,tspan_3,tspan_5,tspan_6,tspan_7,tspan_8]
M_s = [M_1,M_2,M_3,M_3,M_5,M_6,M_7,M_8]
Jac_s = [jac_1,jac_2,jac_3,jac_3,jac_5,[],[],jac_8]
param_s = [p_1,p_2,p_3,p_4,p_5,p_6,p_7,p_8]
numrun_s = [numrun_1,numrun_2,numrun_3,numrun_4,numrun_5,numrun_6,numrun_7,numrun_8]

n_probs = length(ode_s); Pl = []

#--- Benchmarks for Rodas-type methods - --
for ip = 1:n_probs
  println(title_s[ip])
  P_1 = plot(); P_2 = plot()
  if Jac_s[ip]==[]
    f = ODEFunction(ode_s[ip], mass_matrix = M_s[ip])
  else
    f = ODEFunction(ode_s[ip], mass_matrix = M_s[ip], jac = Jac_s[ip]);
  end
  tspan = tspan_s[ip]; prob = ODEProblem(f, y0_s[ip], tspan, param_s[ip])
  if ip==4 #-- take test_sol from ip=3
    abstols = 1.0 ./ 10.0 .^ (3:6);  reltols = 1.0 ./ 10.0 .^ (3:6)
  else
    if ip==7
      sol = solve(prob,Rodas5P(;autodiff=false),abstol=1.0e-14,reltol=1.0e-14,maxiters=Int(1e8))
    else
      sol = solve(prob,Rodas5P(),abstol=4.0e-14,reltol=4.0e-14,maxiters=Int(1e8))
    end
    display(plot(sol))
    global test_sol = TestSolution(sol)
    abstols = 1.0 ./ 10.0 .^ (3:9);  reltols = 1.0 ./ 10.0 .^ (3:9)
    if (ip==6)
      abstols = 1.0 ./ 10.0 .^ (7:12);  reltols = 1.0 ./ 10.0 .^ (7:12)
    end
  end
  println("Test_sol ready")
  if ip==7
    wp = WorkPrecisionSet(prob,abstols,reltols,setups_1; save_everystep=true, error_estimate=:l2,appxsol=test_sol,maxiters=Int(1e8),numruns=numrun_s[ip])
  else
    wp = WorkPrecisionSet(prob,abstols,reltols,setups; save_everystep=true, error_estimate=:l2,appxsol=test_sol,maxiters=Int(1e8),numruns=numrun_s[ip])
  end
  display(plot(wp,title=title_s[ip]))
  P_1 = plot(wp)
  if ip==7
    wp = WorkPrecisionSet(prob,abstols,reltols,setups_1; save_everystep=true, error_estimate=:L2, dense_errors=true,appxsol=test_sol,maxiters=Int(1e8),numruns=numrun_s[ip])
  else
    wp = WorkPrecisionSet(prob,abstols,reltols,setups; save_everystep=true, error_estimate=:L2, dense_errors=true,appxsol=test_sol,maxiters=Int(1e8),numruns=numrun_s[ip])
  end
  P_2 = plot(wp)
  display(plot(wp,title=title_s[ip]))
  push!(Pl,plot(P_1,P_2,layout=(1,2),title=title_s[ip]))
end
#plot(Pl[1])
p = plot(Pl[1],Pl[2],layout=(2,1))
savefig(p, "rodas5p_bench_4.pdf")
p = plot(Pl[3],Pl[4],layout=(2,1))
savefig(p, "rodas5p_bench_5.pdf")
p = plot(Pl[5],Pl[6],layout=(2,1))
savefig(p, "rodas5p_bench_6.pdf")
p = plot(Pl[7],Pl[8],layout=(2,1))
savefig(p, "rodas5p_bench_7.pdf")


text:
# setup clean and empty directories

user_dir = abspath(dirname(@__FILE__), "..", "myMatrixDepot")
data_dir = abspath(dirname(@__FILE__),"..","data")

function save_target(path::AbstractString)
    base = basename(path)
    dir = dirname(path)
    target = abspath(dir, string(base, ".saved"))
    while isdir(target)
        target = abspath(target, string(base, ".saved"))
    end
    if isdir(path)
        mv(path, target)
    end
    if basename(path) == "data"
        mkpath(path)
        cp_if("uf_matrices.html", target, path)
        cp_if("mm_matrices.html", target, path)
        cp_if("ss_index.mat", target, path)
    end
    target
end

function cp_if(name::AbstractString, from::AbstractString, to::AbstractString)
    fromname = abspath(from, name)
    toname = abspath(to, name)
    if isfile(fromname) && !MatrixDepot.URL_REDIRECT[]
        cp(fromname, toname)
    end
end

function revert_target(saved::AbstractString, orig::AbstractString)
    savetest = abspath(dirname(orig), string(basename(orig), ".test"))
    if isdir(orig)
        if isdir(savetest)
            mv(savetest, joinpath(orig, basename(savetest)))
        end
        mv(orig, savetest)
    end
    if isdir(saved)
        mv(saved, orig)
    end
end

ext: jl


text:
domain = "example.com"

records = getalldomainrecords!(client, domain)
@testset "List all Domain Records" begin
    @test length(records) == 4
    @test records[1].id == 28448429
    @test records[1].recordtype == "NS"
    @test records[1].name == "@"
    @test records[1].data == "ns1.digitalocean.com"
    @test records[1].priority == nothing
    @test records[1].port == nothing
    @test records[1].ttl == 1800
    @test records[1].weight == nothing
    @test records[1].flags == nothing
    @test records[1].tag == nothing
end;

record = createdomainrecord!(client, domain, "A"; name="www", data="162.10.66.0",
                             ttl=1800)
@testset "Create a new Domain Record" begin
    @test record.id == 28448433
    @test record.recordtype == "A"
    @test record.name == "www"
    @test record.data == "162.10.66.0"
    @test record.priority == nothing
    @test record.port == nothing
    @test record.ttl == 1800
    @test record.weight == nothing
    @test record.flags == nothing
    @test record.tag == nothing
end;

record = getdomainrecord!(client, "example.com", 28448433)
@testset "Retrieve an existing Domain Record" begin
    @test record.id == 28448433
    @test record.recordtype == "A"
    @test record.name == "www"
    @test record.data == "162.10.66.0"
    @test record.priority == nothing
    @test record.port == nothing
    @test record.ttl == 1800
    @test record.weight == nothing
    @test record.flags == nothing
    @test record.tag == nothing
end;

record = updatedomainrecord!(client, "example.com", record, name="blog")
@testset "Update a Domain Record" begin
    @test record.id == 28448433
    @test record.recordtype == "A"
    @test record.name == "blog"
    @test record.data == "162.10.66.0"
    @test record.priority == nothing
    @test record.port == nothing
    @test record.ttl == 1800
    @test record.weight == nothing
    @test record.flags == nothing
    @test record.tag == nothing
end;

@testset "Delete a Domain Record" begin
    @test deletedomainrecord!(client, "example.com", record)
end;


text:
# TODO: we can manage themes like TestImages.jl

const theme_minimal = """
.card-section {
  display: flex;
  flex-direction: row;
  flex-wrap: wrap;
  align-content: space-between;
}

.card:hover{
  box-shadow: 0 4px 8px 0 rgba(0, 0, 0, 0.4), 0 6px 20px 0 rgba(0, 0, 0, 0.1);
}

.card {
    width: 210px;
    max-height: 400px;
    margin: 10px 15px;
    box-shadow: 0 4px 8px 0 rgba(0,0,0,0.2);
    transition: 0.3s;
    border-radius: 5px;
}

.card-text {
    padding: 0 15px;
}

.card-img {
    width: 200px;
    height: 220px;
    padding: 5px;
    box-shadow: 0 2px 4px 0 rgba(0, 0, 0, 0.2);
    transition: 0.3s;
    border-radius: 5px;
}

"""

function read_cardtheme(theme::AbstractString)
    if theme == "minimal"
        return theme_minimal
    else
        throw(ArgumentError("no theme $(theme) found."))
    end
end


text:
# Stores variables needed to construct the posterior
@kwdef struct DataSet{Td,TCn,TCf,TCf̃,TCϕ,TCn̂,TM,TM̂,TB,TB̂,TD,TG,TP,TL}
    d  :: Td                # data
    Cϕ :: TCϕ               # ϕ covariance
    Cf :: TCf               # unlensed field covariance
    Cf̃ :: TCf̃ = nothing     # lensed field covariance (not always needed)
    Cn :: TCn               # noise covariance
    Cn̂ :: TCn̂ = Cn          # approximate noise covariance, diagonal in same basis as Cf
    M  :: TM  = 1           # user mask
    M̂  :: TM̂  = M           # approximate user mask, diagonal in same basis as Cf
    B  :: TB  = 1           # beam and instrumental transfer functions
    B̂  :: TB̂  = B           # approximate beam and instrumental transfer functions, diagonal in same basis as Cf
    D  :: TD  = IdentityOp  # mixing matrix for mixed parametrization
    G  :: TG  = IdentityOp  # reparametrization for ϕ
    P  :: TP  = 1           # pixelization operator (if estimating field on higher res than data)
    L  :: TL  = alloc_cache(LenseFlow(similar(diag(Cϕ))),d) # a CachedLenseFlow which will be reused for meomry
end

function subblock(ds::DataSet, block)
    DataSet(map(collect(pairs(fields(ds)))) do (k,v)
        @match (k,v) begin
            ((:Cϕ || :G || :L), v)              => v
            (_, L::Union{Nothing,FuncOp,Real})  => L
            (_, L)                              => getindex(L,block)
        end
    end...)
end

function (ds::DataSet)(;θ...)
    DataSet(map(fieldvalues(ds)) do v
        (v isa ParamDependentOp) ? v(;θ...) : v
    end...)
end

function check_hat_operators(ds::DataSet)
    @unpack B̂, M̂, Cn̂, Cf = ds()
    @assert(all([(L isa Scalar) || (L isa typeof(Cf)) || (Cf isa FlatIEBCov && L isa DiagOp{<:FlatIEBFourier}) for L in [B̂,M̂,Cn̂]]),
            "B̂, M̂, Cn̂ should be scalars or the same type as Cf")
end

adapt_structure(to, ds::DataSet) = DataSet(adapt(to, fieldvalues(ds))...)

@doc doc"""
    resimulate(ds::DataSet; f=..., ϕ=...)

Resimulate the data in a given dataset, potentially at a fixed f and/or ϕ (both
are resimulated if not provided)
"""
function resimulate(ds::DataSet; f=simulate(ds.Cf), ϕ=simulate(ds.Cϕ), n=simulate(ds.Cn), f̃=ds.L(ϕ)*f, return_truths=false)
    @unpack M,P,B = ds
    @set! ds.d = M*P*B*f̃ + n
    return_truths ? @namedtuple(ds,f,ϕ,n,f̃) : ds
end

@doc doc"""
    load_sim_dataset

Create a `DataSet` object with some simulated data. E.g.

```julia
@unpack f,ϕ,ds = load_sim_dataset(;
    θpix  = 2,
    Nside = 128,
    pol   = :I,
    T     = Float32
);
```

"""
function load_sim_dataset(;

    # basic configuration
    θpix,
    θpix_data = θpix,
    Nside,
    pol,
    T = Float32,
    storage = Array,

    # noise parameters, or set Cℓn or even Cn directly
    μKarcminT = 3,
    ℓknee = 100,
    αknee = 3,
    Cℓn = nothing,
    Cn = nothing,

    # beam parameters, or set B directly
    beamFWHM = 0,
    B = nothing,

    # mask parameters, or set M directly
    pixel_mask_kwargs = nothing,
    bandpass_mask = LowPass(3000),
    M = nothing,

    # theory
    rfid = 0.05,
    Cℓ = nothing,

    seed = nothing,
    D = nothing,
    G = nothing,
    Nϕ_fac = 2,
    ϕ=nothing, f=nothing, f̃=nothing, Bf̃=nothing, n=nothing, d=nothing, # can override any of these simulated fields
    L = LenseFlow,
    ∂mode = fourier∂
    )

    # the biggest ℓ on the 2D fourier grid
    ℓmax = round(Int,ceil(√2*fieldinfo(Flat(θpix=θpix,Nside=Nside)).nyq)+1)

    # CMB Cℓs
    if Cℓ == nothing
        Cℓ = camb(r=rfid, ℓmax=ℓmax)
    end

    # noise Cℓs (these are non-debeamed, hence beamFWHM=0 below; the beam comes in via the B operator)
    if (Cℓn == nothing)
        Cℓn = noiseCℓs(μKarcminT=μKarcminT, beamFWHM=0, ℓknee=ℓknee, αknee=αknee, ℓmax=ℓmax)
    end

    # some things which depend on whether we chose :I, :P, or :IP
    pol = Symbol(pol)
    S,ks,F,F̂,nF = @match pol begin
        :I  => (S0,  (:TT,),            FlatMap,    FlatFourier,    1)
        :P  => (S2,  (:EE,:BB),         FlatQUMap,  FlatEBFourier,  2)
        :IP => (S02, (:TT,:EE,:BB,:TE), FlatIQUMap, FlatTEBFourier, 3)
        _   => throw(ArgumentError("`pol` should be one of :I, :P, or :IP"))
    end

    # pixelization
    Pix = Flat(Nside=Nside, θpix=θpix, ∂mode=∂mode)
    if (θpix_data == θpix)
        Pix_data = Pix
        P = IdentityOp
    else
        Pix_data = Flat(Nside=Nside÷(θpix_data÷θpix), θpix=θpix_data, ∂mode=∂mode)
        P = FuncOp(
            op  = f -> ud_grade(f, θpix_data, deconv_pixwin=false, anti_aliasing=false),
            opᴴ = f -> ud_grade(f, θpix,      deconv_pixwin=false, anti_aliasing=false)
        )
    end

    # covariances
    Cϕ₀ = Cℓ_to_Cov(Pix,      T, S0, (Cℓ.total.ϕϕ))
    Cfs = Cℓ_to_Cov(Pix,      T, S,  (Cℓ.unlensed_scalar[k] for k in ks)...)
    Cft = Cℓ_to_Cov(Pix,      T, S,  (Cℓ.tensor[k]          for k in ks)...)
    Cf̃  = Cℓ_to_Cov(Pix,      T, S,  (Cℓ.total[k]           for k in ks)...)
    Cn̂  = Cℓ_to_Cov(Pix_data, T, S,  (Cℓn[k]                for k in ks)...)
    if (Cn == nothing); Cn = Cn̂; end
    Cf = ParamDependentOp((mem; r=rfid)->(mem .= Cfs + T(r/rfid)*Cft), similar(Cfs))
    Cϕ = ParamDependentOp((mem; Aϕ=1  )->(mem .= T(Aϕ) .* Cϕ₀), similar(Cϕ₀))

    # data mask
    if (M == nothing)
        M̂ = M = Cℓ_to_Cov(Pix_data, T, S, ((k==:TE ? 0 : 1) * bandpass_mask.diag.Wℓ for k in ks)...)
        if (pixel_mask_kwargs != nothing)
            M = M * Diagonal(F{Pix_data}(repeated(T.(make_mask(Nside÷(θpix_data÷θpix),θpix_data; pixel_mask_kwargs...).Ix),nF)...))
        end
    end

    # beam
    if (B == nothing)
        B̂ = B = Cℓ_to_Cov(Pix, T, S, ((k==:TE ? 0 : 1) * sqrt(beamCℓs(beamFWHM=beamFWHM)) for k=ks)...)
    end

    # D mixing matrix
    if (D == nothing)
        σ²len = T(deg2rad(5/60)^2)
        D = ParamDependentOp(
            function (mem;r=rfid)
                Cfr = Cf(mem,r=r)
                mem .= sqrt(Diagonal(diag(Cfr) .+ σ²len .+ 2*diag(Cn̂)) * pinv(Cfr))
            end,
            similar(Cf())
        )
    end

    # simulate data
    if (seed != nothing); seed!(seed); end
    if (ϕ  == nothing); ϕ  = simulate(Cϕ); end
    if (f  == nothing); f  = simulate(Cf); end
    if (n  == nothing); n  = simulate(Cn); end
    Lϕ = cache(L(ϕ),f)
    if (f̃  == nothing); f̃  = Lϕ*f;         end
    if (Bf̃ == nothing); Bf̃ = B*f̃;          end
    if (d  == nothing); d  = M*P*Bf̃ + n;   end

    # put everything in DataSet
    ds = DataSet(;@namedtuple(d, Cn, Cn̂, Cf, Cf̃, Cϕ, M, M̂, B, B̂, D, P, L=Lϕ)...)

    # with the DataSet created, we can now more conveniently call the quadratic
    # estimate to compute Nϕ if needed for the G mixing matrix
    if (G == nothing)
        Nϕ = quadratic_estimate(ds,(pol in (:P,:IP) ? :EB : :TT)).Nϕ / Nϕ_fac
        G₀ = @. nan2zero(sqrt(1 + 2/($Cϕ()/Nϕ)))
        G = ParamDependentOp((;Aϕ=1)->(@. nan2zero(sqrt(1 + 2/(($(Cϕ(Aϕ=Aϕ))/Nϕ)))/G₀)))
    end
    @set! ds.G = G

    return adapt(storage, @namedtuple(f, f̃, ϕ, n, ds, ds₀=ds(), T, P=Pix, Cℓ, L))

end

###

function load_healpix_sim_dataset(;
    Nside,
    pol,
    gradient_cache,
    T = Float32,
    μKarcminT = 3,
    ℓknee = 100,
    αknee = 3,
    ℓmax_ops = 2Nside,
    ℓmax_data = 3000,
    beamFWHM = 0,
    rfid = 0.05,
    Cℓ = camb(r=rfid, ℓmax=ℓmax_ops),
    Cℓn = nothing,
    Cn = nothing,
    seed = nothing,
    M = nothing,
    B = nothing,
    D = nothing,
    G = nothing,
    ϕ = nothing,
    f = nothing,
    mask_kwargs = nothing,
    L = LenseFlow)

    @assert pol==:T

    # Cℓs
    if (Cℓn == nothing)
        Cℓn = noiseCℓs(μKarcminT, beamFWHM=0, ℓknee=ℓknee, αknee=αknee, ℓmax=ℓmax_ops)
    end
    Cℓf, Cℓf̃ = Cℓ[:f], Cℓ[:f̃]

    Cf = IsotropicHarmonicCov(T.(nan2zero.(Cℓf[:TT][0:ℓmax_ops])), gradient_cache)
    Cf̃ = IsotropicHarmonicCov(T.(nan2zero.(Cℓf̃[:TT][0:ℓmax_ops])), gradient_cache)
    Cn = IsotropicHarmonicCov(T.(nan2zero.(Cℓn[:TT][0:ℓmax_ops])), gradient_cache)
    Cϕ = IsotropicHarmonicCov(T.(nan2zero.(Cℓf[:ϕϕ][0:ℓmax_ops])), gradient_cache)

    P=B=1 #for now

    if (seed != nothing); seed!(seed); end
    if (ϕ==nothing); ϕ = simulate(Cϕ); end
    if (f==nothing); f = simulate(Cf); end
    f̃ = L(ϕ)*f
    n = simulate(Cn)
    d = M*P*B*f̃ + n

    # put everything in DataSet
    ds = DataSet(;@namedtuple(d, Cn, Cf, Cf̃, Cϕ, M)...)

    return @namedtuple(f, f̃, ϕ, n, ds, ds₀=ds())

end


text:
# ------------------------------------------------------------------------------
# No preconditioner
# ------------------------------------------------------------------------------

"""
```julia
IdentityPC()
```

Identity preconditioner to investigate multigrid without smoother

# Returns:
- Identity preconditioner object

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryOperator("mass", 4, 4, mesh);

# preconditioner
identity = IdentityPC(mass);

# verify
println(identity)

# output
identity preconditioner
```
"""
struct IdentityPC <: AbstractPreconditioner
    # data never changes
    operator::Operator

    # inner constructor
    IdentityPC(operator) = new(operator)
end

# printing
# COV_EXCL_START
Base.show(io::IO, preconditioner::IdentityPC) = print(io, "identity preconditioner")
# COV_EXCL_STOP

# COV_EXCL_START
function Base.setproperty!(pc::IdentityPC, f::Symbol, value)
    if f == :operator
        throw(ReadOnlyMemoryError())
    else
        return setfield!(operator, f, value)
    end
end
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# compute symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbols(preconditioner, ω, θ)
```

Compute or retrieve the symbol matrix for a identity preconditioned operator

# Arguments:
- `preconditioner`: Identity preconditioner to compute symbol matrix for
- `ω`:              Smoothing weighting factor array
- `θ`:              Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the identity preconditioner (I)

# Example:
```jldoctest
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryOperator("mass", 4, 4, mesh);

# preconditioner
identity = IdentityPC(mass);

# compute symbols
A = computesymbols(identity, [], []);

# verify
@assert A ≈ I

# output

```
"""
function computesymbols(preconditioner::IdentityPC, ω::Array, θ::Array)
    # return
    return I
end

# ------------------------------------------------------------------------------


text:
# The total board is made up of 9 single boards
mutable struct MetaTac <: Game
  board :: Vector{Int}
  current_player :: Int
  # The focus is either 0 (no focus) or indicates the allowed board
  focus :: Int
  # Stores the Status for the whole game
  status_cache :: Status
  # Stores the status for each region
  region_status_cache :: Vector{Status}
end

function MetaTac()
  MetaTac(zeros(Int, 81), 1, 0, Status(), [Status() for i=1:9])
end

function Base.copy(s :: MetaTac) :: MetaTac
  MetaTac(
    copy(s.board),
    s.current_player,
    s.focus,
    s.status_cache,
    copy(s.region_status_cache)
  )
end

current_player(game :: MetaTac) :: Int = game.current_player

# Returns a list of the Indices of all legal actions
function legal_actions(game :: MetaTac) :: Vector{ActionIndex}
  ActionIndex[ index for index in 1:81 if is_action_legal(game, index) ]
end

# A action is legal, if the following conditions hold:
#   - The game is still active
#   - The region is available
#   - The board position is still empty
function is_action_legal(game :: MetaTac, index :: ActionIndex) :: Bool
  game.board[index] == 0 &&
  is_region_allowed(game, region(index)) &&
  !is_over(game)
end

# Determines, whether a given region can be played in.
function is_region_allowed(game, r) :: Bool
  (game.focus == 0 || game.focus == r) && !is_over(game.region_status_cache[r])
end

# Finds the region 1d index for a given cell 1d index, value is in 1:9
function region(index) :: Int
  col = mod(div(index - 1, 3), 3)
  row = div(index-1, 27)
  1 + col + 3 * row
end

# Finds the inner index for a given cell index, value is in 1:9
function local_index(index) :: Int
  col = mod1(index, 3)
  row = mod(div(index-1, 9), 3)
  col + 3 * row
end

# Copies a region into a new matrix.
function region_view(game, r)
  row = 3 * mod(r-1, 3) + 1
  col = 3 * div(r-1, 3) + 1
  reshape(game.board, (9, 9))[row:row + 2, col:col + 2]
end

function apply_action!(game :: MetaTac, index :: ActionIndex) :: MetaTac
  @assert is_action_legal(game, index) "Action $index is not allowed."
  # Update the board state
  game.board[index] = game.current_player
  game.current_player = -game.current_player

  # Update the status cache
  r = region(index)
  game.region_status_cache[r] = tic_tac_toc_status(region_view(game, r))
  game.status_cache = tic_tac_toc_status(game.region_status_cache)

  # Update the focus. This must be done AFTER the status cache is refreshed.
  game.focus = local_index(index)
  if is_over(game.region_status_cache[game.focus])
    game.focus = 0
  end
  game
end

status(game :: MetaTac) :: Status = game.status_cache

# Returns ( false, * ) => Game is not over yet
# Returns ( true, {-1, 0, 1} ) => Game is over and we return the winner
function tic_tac_toc_status(game :: MetaTac, outer_index) :: Tuple{ Bool, Int }
  start_index = combined_index(outer_index, 1)
  single_board = game.board[start_index:start_index+8]
  tic_tac_toc_status(single_board)
end

# Implements the win condition for the large board
# We need to additionaly verify, that there are legal actions left
function tic_tac_toc_status(board :: Vector{Status}) :: Status
  matrix = reshape(with_default.(board, 0), (3, 3))
  s = tic_tac_toc_status(matrix)
  if is_over(s)
    s
  elseif all(is_over, board) #all(is_over.(board))
    Status(0)
  else
    Status()
  end
end

# Takes a 3 * 3 Array of {1, 0, -1} and returns the status for this board.
function tic_tac_toc_status(board :: Matrix{Int}) :: Status
  # Iterate rows and columns
  for i = 1:3
    if board[1,i] == board[2,i] == board[3,i] != 0
      return Status(board[1, i])
    elseif board[i,1] == board[i,2] == board[i,3] != 0
      return Status(board[i, 1])
    end
  end
  # Diagonals
  if board[1, 1] == board[2, 2] == board[3, 3] != 0
    Status(board[2, 2])
  elseif board[1, 3] == board[2, 2] == board[3, 1] != 0
    Status(board[2, 2])
  elseif all(x -> x != 0, board)
    Status(0)
  else
    Status()
  end
end

policy_length(:: Type{MetaTac}) :: Int = 81

# Size of the data representation of the game
Base.size(:: Type{MetaTac}) :: Tuple{Int, Int, Int} = (9, 9, 3)

# Data representation of the game as layered 2d image
function representation(game :: MetaTac) :: Array{Float32, 3}
  data = zeros(Float32, 81, 3)
  board = game.current_player .* game.board
  data[board .== 1, 1] .= 1
  data[board .== -1, 2] .= 1
  data[legal_actions(game), 3] .= 1
  reshape(data, (9, 9, 3))
end

function augment(game :: MetaTac)

  boards = apply_dihedral_group(reshape(game.board, (9,9)))
  caches = apply_dihedral_group(reshape(game.region_status_cache, (3,3)))

  if game.focus == 0
    focii = fill(0, length(boards))
  else
    focii = apply_dihedral_group(game.focus, (3,3))
  end

  map(boards, caches, focii) do board, cache, focus
    board = reshape(board, (81,))
    cache = reshape(cache, (9,))
    MetaTac(board, game.current_player, focus, game.status_cache, cache)
  end

end

function augment(game :: MetaTac, label :: Vector{Float32})

  matpol = reshape(label[2:end], (9, 9))
  matpols = apply_dihedral_group(matpol)
  labels = [ vcat(label[1], reshape(mp, (81,))) for mp in matpols ]

  augment(game), labels

end

function draw(game :: MetaTac) :: Nothing

  board = reshape(game.board, (9,9))
  symbols = Dict(1 => "X", -1 => "O", 0 => "⋅")

  for i in 1:9
    for j in 1:9
      print(" $(symbols[board[j,i]])")
      if j == 3 || j == 6
        print(" ║")
      end
    end
    println()
    if i == 3 || i == 6
      println("═══════╬═══════╬═══════")
    end
  end
end

max_stars_repo_name: roSievers/jtac


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 794a5356-6e05-44b4-ad9d-375a1fa3613b
begin
	using PlutoUI
	using DataFrames
	using Test
end

# ╔═╡ b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
md"""
## Introduction DataFrames / 8 - Joining DataFrames

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/04_loadsave.ipynb)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ d6993a5b-8158-4ce3-accb-8e8c514229b3
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ de97ded0-3952-4188-9a43-a8f383e1ae65
md"""
### Preparing DataFrames for a join
"""

# ╔═╡ 815474c9-8d95-44c0-ac6c-4df825a61b09
x = DataFrame(ID=[1, 2, 3, 4, missing], name = ["Alice", "Bob", "Conor", "Dave","Zed"])

# ╔═╡ f4362c1b-c787-4b57-9279-00dbac74e445
y = DataFrame(id=[1, 2, 5, 6, missing], age = [21, 22, 23, 24, 99])

# ╔═╡ 3d73280f-89dc-4939-9fdf-11f9cadc3944
md"""
Rules for the on keyword argument:

  - a single Symbol or string if joining on one column with the same name, e.g. `on=:id`
  - a Pair of Symbols or string if joining on one column with different names, e.g. on=:id => :id2`
  - a vector of Symbols or strings if joining on multiple columns with the same name, e.g. `on=[:id1, :id2]`
  - a vector of Pairs of Symbols or strings if joining on multiple columns with the same name, e.g. `on=[:a1=>:a2, :b1=>:b2]`
  - a vector containing a combination of Symbols or strings or Pair of Symbols or strings, e.g. `on=[:a1=>:a2, :b1]`

"""

# ╔═╡ ca9927e5-3ff2-4c73-8206-4d7495668312
md"""
### Standard joins: inner, left, right, outer, semi, anti

"""

# ╔═╡ a1ab6cc4-61e1-47dd-a573-8c8884345e5c
@test_throws ArgumentError innerjoin(x, y, on=:ID=>:id)
## missing is not allowed to join-on by default

# ╔═╡ e155fefb-95ae-469e-bcc0-479bf321b1f9
innerjoin(x, y, on=:ID=>:id, matchmissing=:equal)

# ╔═╡ 5149af4b-4511-443b-abf5-9b53e38e1d45
leftjoin(x, y, on="ID"=>"id", matchmissing=:equal)

# ╔═╡ 0698c48c-b0d8-47a4-a1fb-eb93d0602815
rightjoin(x, y, on=:ID=>:id, matchmissing=:equal)

# ╔═╡ a7cf86b6-11cf-4e24-9d3f-9abec87503a6
outerjoin(x, y, on=:ID=>:id, matchmissing=:equal)

# ╔═╡ aec2bdaa-3e3f-42fb-8747-d9397148983c
semijoin(x, y, on=:ID=>:id, matchmissing=:equal)

# ╔═╡ e2d707ab-13b1-4f85-a688-d4f358a9c1d5
antijoin(x, y, on=:ID=>:id, matchmissing=:equal)

# ╔═╡ a4615e9b-d3d9-443a-9146-e411f9bbd749
md"""
### Cross join
"""

# ╔═╡ 7bc6197d-2c03-471b-b7ac-22f78fa6d383
crossjoin(DataFrame(x=[1, 2]), DataFrame(y=["a", "b", "c"]))  # no `on` argument

# ╔═╡ 79de555b-1858-4fda-adc7-d29c503a0f96

# ╔═╡ 1688d806-c5ef-40dd-8b09-2f9de434023b
md"""
### Complex cases of joins
"""

# ╔═╡ 4b1ba677-fa53-4acd-ba44-4f098e9e4ed5
cx = DataFrame(id1=[1, 1, 2, 2, missing, missing],
              id2=[1, 11, 2, 21, missing, 99],
              name = ["Alice", "Bob", "Conor", "Dave","Zed", "Zoe"])

# ╔═╡ c6a8e4af-63bb-49a6-b810-877a6fe7ab5d
cy = DataFrame(id1=[1, 1, 3, 3, missing, missing],
              id2=[11, 1, 31, 3, missing, 999],
              age = [21, 22, 23, 24, 99, 100])

# ╔═╡ b0015c63-6ee8-422f-95e5-98510b4c893e
innerjoin(cx, cy, on=[:id1, :id2], matchmissing=:equal)  ## joining on two columns

# ╔═╡ 1b8def29-f289-47cb-86ca-0b511e92e193
outerjoin(cx, cy, on=:id1, makeunique=true, indicator=:source, matchmissing=:equal)
## with duplicates all combinations are produced

# ╔═╡ bfecee96-e044-40a6-a2c6-ab4fea31c975
## you can force validation of uniqueness of key on which you join

@test_throws ArgumentError innerjoin(cx, cy, on=:id1, makeunique=true, validate=(true, true),
	matchmissing=:equal)

# ╔═╡ 2ea96ebd-86d1-4def-8404-911d651b6dbf
md"""
mixed on argument for joining on multiple columns
"""

# ╔═╡ 4eae57b0-c23a-47c8-931c-9dc57508a7cc
x₁ = DataFrame(id1=1:6, id2=[1, 2, 1, 2, 1, 2], x1 = 'a':'f')

# ╔═╡ f58a6cbf-38d9-4874-ada7-39812c81a874
y₁ = DataFrame(id1=1:6, ID2=1:6, x2 = 'a':'f')

# ╔═╡ 68b6d030-0c1a-4a7c-913f-7e164eb835e0
innerjoin(x₁, y₁, on=[:id1, :id2=>:ID2])

# ╔═╡ 7434a9b0-9995-4c90-94ca-20778e677ea2
md"""
joining more than two data frames:
"""

# ╔═╡ 8aefe960-02eb-484b-bdf7-ade42d90295d
xs = [DataFrame("id"=>1:6, "v$i"=>((1:6) .+ 10i)) for i in 1:5]

# ╔═╡ e9909c0b-d57b-413b-bccc-0638ea66d43b
innerjoin(xs..., on=:id)   ## also for outerjoin and crossjoin

# ╔═╡ c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e
html"""

"""

# ╔═╡ Cell order:
# ╟─b3dcf3f4-9d68-11eb-01a6-317fe0e0cd7e
# ╠═794a5356-6e05-44b4-ad9d-375a1fa3613b
# ╟─d6993a5b-8158-4ce3-accb-8e8c514229b3
# ╟─de97ded0-3952-4188-9a43-a8f383e1ae65
# ╠═815474c9-8d95-44c0-ac6c-4df825a61b09
# ╠═f4362c1b-c787-4b57-9279-00dbac74e445
# ╟─3d73280f-89dc-4939-9fdf-11f9cadc3944
# ╟─ca9927e5-3ff2-4c73-8206-4d7495668312
# ╠═a1ab6cc4-61e1-47dd-a573-8c8884345e5c
# ╠═e155fefb-95ae-469e-bcc0-479bf321b1f9
# ╠═5149af4b-4511-443b-abf5-9b53e38e1d45
# ╠═0698c48c-b0d8-47a4-a1fb-eb93d0602815
# ╠═a7cf86b6-11cf-4e24-9d3f-9abec87503a6
# ╠═aec2bdaa-3e3f-42fb-8747-d9397148983c
# ╠═e2d707ab-13b1-4f85-a688-d4f358a9c1d5
# ╟─a4615e9b-d3d9-443a-9146-e411f9bbd749
# ╠═7bc6197d-2c03-471b-b7ac-22f78fa6d383
# ╠═79de555b-1858-4fda-adc7-d29c503a0f96
# ╟─1688d806-c5ef-40dd-8b09-2f9de434023b
# ╠═4b1ba677-fa53-4acd-ba44-4f098e9e4ed5
# ╠═c6a8e4af-63bb-49a6-b810-877a6fe7ab5d
# ╠═b0015c63-6ee8-422f-95e5-98510b4c893e
# ╠═1b8def29-f289-47cb-86ca-0b511e92e193
# ╠═bfecee96-e044-40a6-a2c6-ab4fea31c975
# ╟─2ea96ebd-86d1-4def-8404-911d651b6dbf
# ╠═4eae57b0-c23a-47c8-931c-9dc57508a7cc
# ╠═f58a6cbf-38d9-4874-ada7-39812c81a874
# ╠═68b6d030-0c1a-4a7c-913f-7e164eb835e0
# ╟─7434a9b0-9995-4c90-94ca-20778e677ea2
# ╠═8aefe960-02eb-484b-bdf7-ade42d90295d
# ╠═e9909c0b-d57b-413b-bccc-0638ea66d43b
# ╟─c54cdbf1-7787-4f0f-b9ff-d4389f4a7c5e

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using RunTests

@testmodule TestModuleWithNoTests begin
end


text:
import CUDA
using Knet
using Images
using Distributed
include("./deeper-cut/deeper-cut.config.jl")
include("./deeper-cut/mpii.annotation.reader.jl")
using DataFrames
using CSV

function compute_loss_for_sets(model, named_x_y_tuples)
    loss_dict = Dict()
    for named_x_y_tuple in named_x_y_tuples
        name = named_x_y_tuple[1]
        x = named_x_y_tuple[2]
        y = named_x_y_tuple[3]
        loss_dict["$name-loss"] = model(x, y)
    end
    return loss_dict
end

function compute_accuracy_for_sets(model, named_x_y_tuples)
    accuracy_dict = Dict()
    for named_x_y_tuple in named_x_y_tuples
        name = named_x_y_tuple[1]
        x = named_x_y_tuple[2]
        y = named_x_y_tuple[3]
        current_res = model(x)
        accuracy_dict["$name-accuracy-mae"] = mae(current_res, y)
    end
    return accuracy_dict
end

function compute_mae_for_data(model, data)
    sum = 0
    count = 0
    for (x, y) in data
        current_res = model(x)
        sum += simple_mae(current_res, y)
        count += 1
    end
    return sum / count
end

function simple_mae(x, y)
    sum(abs, (x - y)) / size(y)[end]
end

function simple_mse(x, y)
    sum(abs2, (x - y)) / size(y)[end]
end

function clear_gpu_memory()
    CUDA.reclaim()
    CUDA.memory_status()
    GC.gc(true)
end

function substract_mean_img_from_data!(data, mean_pixel, atype)
    step_size = 2000
    data_length = length(data)
    for i = 1:step_size:data_length
        upper_bound = min(i + step_size, data_length)
        data.x[i:upper_bound] = map(
            element ->
                substract_mean_img_from_data_for_element(element, mean_pixel, atype),
            enumerate(data.x[i:upper_bound]),
        )
    end
end

function substract_mean_img_from_data_for_element(enumerated_img, mean_pixel, atype)
    x = enumerated_img[2]
    count = enumerated_img[1]
    if count % 1024 == 1
        println(count)
    end
    result = Array(x) .- mean_pixel
    return atype(result)
end

function mirror_variant_data(variant_data)
    new_x = deepcopy(variant_data.x)
    new_y = deepcopy(variant_data.y)
    counter = 1
    original_shuffle = variant_data.shuffle
    variant_data.shuffle = false
    for (x, y) in variant_data
        arr_x = Array(x)
        arr_y = Array(y)
        arr_x = reduce_dim(arr_x)
        arr_y = reduce_dim(arr_y)
        flipped_x = fliplr3D(arr_x)
        flipped_y = exchange_symmetric_joint_ids_for_flipped_labels(fliplr3D(arr_y))
        flipped_x = add_dim(flipped_x) |> variant_data.xtype
        flipped_y = add_dim(flipped_y) |> variant_data.ytype
        push!(new_x, flipped_x)
        push!(new_y, flipped_y)
        counter += 1
        if counter % 1024 == 1
            println(counter)
        end
    end
    variant_minibatch(
        new_x,
        new_y;
        shuffle = original_shuffle,
        xtype = variant_data.xtype,
        ytype = variant_data.ytype,
    )
end

function mirror_raw_batch_item(batch_item, is_label, atype)
    arr = Array(batch_item)
    arr = reduce_dim(arr)
    flipped =
        is_label ? exchange_symmetric_joint_ids_for_flipped_labels(fliplr3D(arr)) :
        fliplr3D(arr)
    flipped = add_dim(flipped) |> atype
    return flipped
end

function mirror_data_items_for_batch(data_items)
    mirrored_data_items = map(mirror_data_item, data_items)
    return [data_items..., mirrored_data_items...]
end

function fliplr3D(img)
    copied_img = deepcopy(img)
    y_dim = size(img)[1]
    x_dim = size(img)[2]
    for y = 1:y_dim
        for x = 1:x_dim
            copied_img[y, x_dim-x+1, :] = img[y, x, :]
        end
    end
    return copied_img
end

function fliplr2D(img)
    copied_img = deepcopy(img)
    y_dim = size(img)[1]
    x_dim = size(img)[2]
    for y = 1:y_dim
        for x = 1:x_dim
            copied_img[y, x_dim-x+1] = img[y, x]
        end
    end
    return copied_img
end

function exchange_symmetric_joint_ids_for_flipped_labels(labels)
    copied_labels = deepcopy(labels)
    h, w, c = size(labels)

    for i = 1:global_num_joints
        copied_labels[:, :, i] = labels[:, :, symmetric_joint_dict[i]]
    end

    for i = global_num_joints+1:2:c
        if i < global_num_joints * 3 + 1
            j = (i - global_num_joints) / 2 |> ceil |> Int
            sym_j = symmetric_joint_dict[j]
            new_x = (sym_j - j) * 2 + i
            new_y = (sym_j - j) * 2 + 1 + i
            copied_labels[:, :, i] = labels[:, :, new_x]
            copied_labels[:, :, i+1] = labels[:, :, new_y]
        else
            j = (i - global_num_joints * 3) / 2 |> ceil |> Int
            sym_j = symmetric_joint_dict[j]
            new_x = (sym_j - j) * 2 + i
            new_y = (sym_j - j) * 2 + 1 + i
            copied_labels[:, :, i] = labels[:, :, new_x]
            copied_labels[:, :, i+1] = labels[:, :, new_y]
        end
    end
    return copied_labels
end

function reverse_all_dims(array; until_dim)
    dim_count = array |> size |> length
    for i = 1:dim_count
        if until_dim == i
            break
        end
        array = reverse(array, dims = i)
    end
    return array
end

function write_acc_results_to_csv(filename, acc_results)
    df = DataFrame(
        Ankle1 = Float32[],
        Knee1 = Float32[],
        Hip1 = Float32[],
        Hip2 = Float32[],
        Knee2 = Float32[],
        Ankle2 = Float32[],
        Wrist1 = Float32[],
        Elbow1 = Float32[],
        Shoulder1 = Float32[],
        Shoulder2 = Float32[],
        Elbow2 = Float32[],
        Wrist2 = Float32[],
        Chin = Float32[],
        TopHead = Float32[],
    )
    for i = 1:size(acc_results, 1)
        push!(df, acc_results[i, :])
    end

    CSV.write("$(filename).csv", df)

end


text:
# complex polya

# makie

using CairoMakie
let
        x = y = -2:0.005:2
        f(z) = 1 / (z * (z^2 - z - 1 - 3im))
        fvals = [f(u + 1im * v) for u in x, v in y]
        fvalues = abs.(fvals)
        fargs = angle.(fvals)
        polya(x, y) = Point2f(real(f(x + 1im * y)), -imag(f(x + 1im * y)))

        fig = Figure(resolution = (900, 400))
        axs = [Axis(fig[1, i], aspect = 1) for i in 1:2]
        cmap = :roma
        streamplot!(axs[1], polya, -2 .. 2, -2 .. 2, colormap = (:black, :red),
                gridsize = (40, 40), arrow_size = 6, linewidth = 1)
        pltobj2 = CairoMakie.heatmap!(axs[2], x, y, fargs, colorrange = (-π, π), colormap = cmap)
        CairoMakie.streamplot!(axs[2], polya, -2 .. 2, -2 .. 2, colormap = (:black, :black),
                gridsize = (40, 40), arrow_size = 6, linewidth = 1)
        Colorbar(fig[1, 3], pltobj2, ticks = ([-π, -π / 2, 0, π / 2, π],
                [L"-\pi", L"-\pi/2", L"0", L"\pi/2", L"\pi"]))
        limits!(axs[1], -2, 2, -2, 2)
        limits!(axs[2], -2, 2, -2, 2)
        colsize!(fig.layout, 1, Aspect(1, 1.0))
        colsize!(fig.layout, 2, Aspect(1, 1.0))
        display(fig)
end;

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using Images, ImageView

# crop an image
source_image = load("sample-images/cats-3061372_640.jpg");
size(source_image)
cropped_image = source_image[100:290, 280:540]; # (height, width)
imshow(cropped_image)

# crop an image (option 2)
source_image = load("sample-images/cats-3061372_640.jpg");
cropped_image_view = view(source_image, 100:290, 280:540); # (height, width)
imshow(cropped_image_view)

# resize an image
source_image = load("sample-images/cats-3061372_640.jpg");
resized_image = imresize(source_image, (100, 250)); # (height, width)
imshow(resized_image);

# resize an image (2)
source_image = load("sample-images/cats-3061372_640.jpg");
resized_image = imresize(source_image, (200, 200)); # (height, width)
imshow(resized_image);

# scale an image by percentage
source_image = load("sample-images/cats-3061372_640.jpg");
scale_percentage = 0.6
new_size = trunc.(Int, size(source_image) .* scale_percentage)
scaled_image = imresize(source_image, new_size)
imshow(resized_image);

# scale an image to a specific dimension
source_image = load("sample-images/cats-3061372_640.jpg");
new_width = 200
scale_percentage = new_width / size(source_image)[2]
new_size = trunc.(Int, size(source_image) .* scale_percentage)
scaled_image = imresize(source_image, new_size)
imshow(scaled_image);

# scale/ resize image by two-fold
source_image = load("sample-images/cats-3061372_640.jpg");
resized_image = restrict(source_image, 1); # height
imshow(resized_image);


text:
module Elements
import GaussQuadrature

"""
    lglpoints(::Type{T}, N::Integer) where T <: AbstractFloat

returns the points `r` and weights `w` associated with the `N+1`-point
Gauss-Legendre-Lobatto quadrature rule of type `T`

"""
function lglpoints(::Type{T}, N::Integer) where {T <: AbstractFloat}
    @assert N ≥ 1
    GaussQuadrature.legendre(T, N + 1, GaussQuadrature.both)
end

"""
    lgpoints(::Type{T}, N::Integer) where T <: AbstractFloat

returns the points `r` and weights `w` associated with the `N+1`-point
Gauss-Legendre quadrature rule of type `T`
"""
function lgpoints(::Type{T}, N::Integer) where {T <: AbstractFloat}
    @assert N ≥ 1
    GaussQuadrature.legendre(T, N + 1, GaussQuadrature.neither)
end

"""
    baryweights(r)

returns the barycentric weights associated with the array of points `r`

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
function baryweights(r::AbstractVector{T}) where {T}
    Np = length(r)
    wb = ones(T, Np)

    for j in 1:Np
        for i in 1:Np
            if i != j
                wb[j] = wb[j] * (r[j] - r[i])
            end
        end
        wb[j] = T(1) / wb[j]
    end
    wb
end

"""
    spectralderivative(r::AbstractVector{T},
                       wb=baryweights(r)::AbstractVector{T}) where T

returns the spectral differentiation matrix for a polynomial defined on the
points `r` with associated barycentric weights `wb`

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
function spectralderivative(
    r::AbstractVector{T},
    wb = baryweights(r)::AbstractVector{T},
) where {T}
    Np = length(r)
    @assert Np == length(wb)
    D = zeros(T, Np, Np)

    for k in 1:Np
        for j in 1:Np
            if k == j
                for l in 1:Np
                    if l != k
                        D[j, k] = D[j, k] + T(1) / (r[k] - r[l])
                    end
                end
            else
                D[j, k] = (wb[k] / wb[j]) / (r[j] - r[k])
            end
        end
    end
    D
end

"""
    interpolationmatrix(rsrc::AbstractVector{T}, rdst::AbstractVector{T},
                        wbsrc=baryweights(rsrc)::AbstractVector{T}) where T

returns the polynomial interpolation matrix for interpolating between the points
`rsrc` (with associated barycentric weights `wbsrc`) and `rdst`

Reference:
  Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange Interpolation",
  SIAM Review 46 (2004), pp. 501-517.
  <https://doi.org/10.1137/S0036144502417715>
"""
function interpolationmatrix(
    rsrc::AbstractVector{T},
    rdst::AbstractVector{T},
    wbsrc = baryweights(rsrc)::AbstractVector{T},
) where {T}
    Npdst = length(rdst)
    Npsrc = length(rsrc)
    @assert Npsrc == length(wbsrc)
    I = zeros(T, Npdst, Npsrc)
    for k in 1:Npdst
        for j in 1:Npsrc
            I[k, j] = wbsrc[j] / (rdst[k] - rsrc[j])
            if !isfinite(I[k, j])
                I[k, :] .= T(0)
                I[k, j] = T(1)
                break
            end
        end
        d = sum(I[k, :])
        I[k, :] = I[k, :] / d
    end
    I
end

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.736879506320445 <= q <= 22.217442927624916)

                     add_NL_constraint(m, :(log(1 + exp(-0.3136061325511519 + 0.2821341474705279*log(1 + exp(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))) + 0.4396178028948987*log(1 + exp(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))) + -0.24849525966922226*log(1 + exp(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))) + 0.5972646324288218*log(1 + exp(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))) + -0.8620507270647968*log(1 + exp(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))) + log(1 + exp(-0.9705953045797062 + -0.6259685253929317*log(1 + exp(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))) + 0.41918630672941015*log(1 + exp(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))) + -0.9388772793667424*log(1 + exp(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))) + -0.18849282727638839*log(1 + exp(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))) + -0.7819718693068158*log(1 + exp(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))) + log(1 + exp(0.8285261526161638 + 0.6863822914423507*log(1 + exp(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))) + -0.7060930690780838*log(1 + exp(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))) + -0.2956343589188948*log(1 + exp(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))) + 0.32353098534638214*log(1 + exp(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))) + 0.6963022628505415*log(1 + exp(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))) + log(1 + exp(-0.2246585184074279 + 0.8591320330102441*log(1 + exp(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))) + -0.5037522661683527*log(1 + exp(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))) + 0.9621050832972355*log(1 + exp(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))) + -0.1995700954953512*log(1 + exp(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))) + 0.442251462608958*log(1 + exp(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))) + log(1 + exp(0.5761249271362421 + -0.042724357144423575*log(1 + exp(0.7833824073251248 + 0.1390648465658577*$(x[1]) + -0.6617731040629309*$(x[2]) + -0.5974574624055107*$(x[3]))) + -0.4537006876708345*log(1 + exp(-0.7125949557702831 + 0.9109805514170546*$(x[1]) + 0.42115237771853087*$(x[2]) + 0.7134573642319126*$(x[3]))) + -0.12519521024691826*log(1 + exp(-0.36270963442709325 + 0.07407412170431016*$(x[1]) + -0.6239570881148437*$(x[2]) + -0.6542866380333012*$(x[3]))) + -0.7426009177864223*log(1 + exp(-0.5451890615285553 + -0.19626106399785304*$(x[1]) + 0.5284606976702149*$(x[2]) + 0.2694706847237054*$(x[3]))) + -0.8758186284143128*log(1 + exp(0.35158020081654806 + 0.8300503955389886*$(x[1]) + -0.8477304064925764*$(x[2]) + 0.1840810006485709*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#aproximates liquid volume at a known pressure and t,
#by using isothermal compressibility

"""
    volume_compress(model,p,T,z=SA[1.0];V0=x0_volume(model,p,T,z,phase=:liquid),max_iters=100)

Main routine to calculate a volume, given a pressure, temperature, composition and intitial volume guess. each step is taken by locally aproximating the EoS as an isothermal compressibility process.
The new volume is calculated by the following recurrence formula:

```julia
v[i+1] = v[i]*exp(β[i]*(p-p(v[i])))
```

In the liquid root region, the iterations follow `v0 < v[i] < v[i+1] < v(p)`, allowing the calculation of the liquid root without entering the metastable region.
"""
function volume_compress(model,p,T,z=SA[1.0];V0=x0_volume(model,p,T,z,phase=:liquid),max_iters=100)
    _0 = zero(p+T+first(z))
    _nan = _0/_0
    logV0 = log(V0)
    lb_v = lb_volume(model,z)
    function logstep(_V)
        _V < log(lb_v) && return zero(_V)/zero(_V)
        _V = exp(_V)
        _p,dpdV = p∂p∂V(model,_V,T,z)
        _Δ = (p-_p)/(_V*dpdV)
        return _Δ
    end

    function f_fixpoint(_V)
        Δ = logstep(_V)
        _V + sign(Δ)*abs(Δ)^(one(Δ)+Δ)
    end

    Δ0 = logstep(logV0)
    Δ00 = Δ0
    Δ1 = logstep(logV0+Δ0)
    logV1 = logV0
    iter = 0

    #convergence hill. we do Successive substitution until we land on a
    #good initial point for aitken acceleration
    for i in 1:20
        if abs(Δ00) > abs(Δ1)
            break
        end
        if !isfinite(Δ1)
            return _nan
        end
        Δ0 = Δ1
        Δ1 = logstep(logV1+Δ0)
        logV1 = logV1 + Δ1
        iter +=1
    end
    if iter == 20
        return _nan
    end
    res = @nan(Solvers.fixpoint(f_fixpoint,logV1,Solvers.SSFixPoint(),rtol = 1e-12,max_iters=max_iters),_nan)

    return exp(res)
end

"""
    volume_virial(model::EoSModel,p,T,z=SA[1.0])
    volume_virial(B::Real,p,T,z=SA[1.0])

Calculates an aproximation to the gas volume at specified pressure, volume and composition, by aproximating:

```julia

Z(v) ≈ 1 + B(T)/v
```
where `Z` is the compressibility factor and `B` is the second virial coefficient.
If `B>0`, (over the inversion temperature) returns `NaN`. If the solution to the problem is complex (`Z = 1 + B/v` implies solving a quadratic polynomial), returns `-2*B`.

If you pass an `EoSModel` as the first argument, `B` will be calculated from the EoS at the input `T`. You can provide your own second virial coefficient instead of a model.

"""
function volume_virial end

function volume_virial(model::EoSModel,p,T,z=SA[1.0])
    B = second_virial_coefficient(model,T,z)
    return volume_virial(B,p,T,z)
end

function volume_virial(B::Real,p,T,z=SA[1.0])
    _0 = zero(B)
    B > _0 && return _0/_0
    a = p/(R̄*T)
    b = -1
    c = -B
    Δ = b*b-4*a*c
    n = sum(z)
    if Δ <= 0
        #virial approximation could not be calculated
        #return value at spinodal
        return -2*B
    end
    #only the left root has physical meaning
    return (-b + sqrt(b*b-4*a*c))/(2*a)
end

#(z = pV/RT)
#(RT/p = V/z)
"""
    volume(model::EoSModel,p,T,z=SA[1.0];phase=:unknown,threaded=true)

calculates the volume (m³) of the compound modelled by `model` at a certain pressure,temperature and moles.

`phase` is a Symbol that determines the initial volume root to look for:

- If `phase =:liquid` it will return the volume of the phase using a liquid initial point.

- If `phase =:vapor` it will return the volume of the phase using a gas initial point.

The default is `phase =:unknown`. with this, both liquid and volume roots will be calculated, and
the phase with the least amount of energy is returned.

The calculation of both volume roots can be calculated in serial (`threaded=false`) or in parallel (`threaded=true`)

"""
function volume(model::EoSModel,p,T,z=SA[1.0];phase=:unknown,threaded=true)

    fp(_V) = log(pressure(model,_V,T,z)/p)

#Threaded version
    phase = Symbol(phase)
    if phase != :unknown
        V0 = x0_volume(model,p,T,z,phase=phase)
        #return Solvers.ad_newton(fp,Vg0)
        return volume_compress(model,p,T,z,V0=V0)
    end

    if threaded
            Vg0 = x0_volume(model,p,T,z,phase=:v)
            Vl0 = x0_volume(model,p,T,z,phase=:l)
            _Vg = Threads.@spawn volume_compress(model,$p,$T,$z;V0=Vg0)
            _Vl = Threads.@spawn volume_compress(model,$p,$T,$z;V0=$Vl0)
            Vg = fetch(_Vg)
            Vl = fetch(_Vl)
    else
            Vg0 = x0_volume(model,p,T,z,phase=:v)
            Vl0 = x0_volume(model,p,T,z,phase=:l)
            Vg =  volume_compress(model,p,T,z,V0=Vg0)
            Vl =  volume_compress(model,p,T,z,V0=Vl0)
    end

# Serial version
#=
    Vg0 = volume_virial(model,p,T,z)
    Vg =Solvers.ad_newton(fp,Vg0,rtol=1e-08)

    Vl =  volume_compress(model,p,T,z)
=#
    function gibbs(V)
        _df,f =  ∂f(model,V,T,z)
        dV,dt = _df
        if abs((p+dV)/p) > 0.03
            return Inf
        else
            return f  +p*V
        end
    end
    #this catches the supercritical phase as well

    isnan(Vl) && return Vg
    isnan(Vg) && return Vl
    (isnan(Vl) & isnan(Vg)) && error("Failed to converge to a root")

    gg = gibbs(Vg)
    gl = gibbs(Vl)
    #@show Vg,Vl
    #@show gg,gl
    return ifelse(gg<gl,Vg,Vl)

end

function volume(model::IdealModel,p,T,z=SA[1.0];phase=:unknown,threaded=false)
    return sum(z)*R̄*T/p
end

export volume


text:
using ForwardDiff
using PyPlot
using NLopt
using Ipopt
using JuMP
import ForwardDiff: value
using FixedSizeArrays
import SoftRobots
reload("SoftRobots")

function potential_energy{T}(robot, state::SoftRobots.SoftRobotState{T})
    energy = zero(T)
    for edge in robot.edges
        energy += 0.5 * edge.stiffness *
         (norm(state.positions[edge.parent] - state.positions[edge.child]) -
          edge.rest_length)^2
    end
    energy
end

Base.convert{T}(::Type{SoftRobots.SoftRobotState{T}}, x::Vector{T}) = SoftRobots.SoftRobotState(SoftRobots.Point{3, T}[x[i:i+2]
    for i in 1:3:(3*length(robot.nodes)-2)],
    [SoftRobots.Point{3, T}(0) for i in 1:length(robot.nodes)])

potential_energy{T}(robot, x::Vector{T}) = potential_energy(robot, convert(SoftRobots.SoftRobotState{T}, x))

draw{T}(robot, x::Vector{T}) = SoftRobots.draw(robot, convert(SoftRobots.SoftRobotState{T}, x))

value(x::Real) = x

function tower()
    k = 10.
    b = 4.0
    l = 1.0
    m = 0.1
    nodes = SoftRobots.PointMass[]
    positions = SoftRobots.Point{3, Float64}[]
    for z = 0:1:3
        bulge = 1 - 0.1*(z - 1.5)^2
        for x = -0.5:0.5
            for y = -0.5:0.5
                push!(nodes, SoftRobots.PointMass(m))
                push!(positions, SoftRobots.Point{3, Float64}(x*bulge, y*bulge, z))
            end
        end
    end
    num_nodes = length(nodes)
    edges = SoftRobots.DampedSpring[]
    for i = 1:num_nodes
        for j = (i+1):num_nodes
            rest_length = norm(positions[i] - positions[j])
            push!(edges, SoftRobots.DampedSpring(i, j, k, b, rest_length))
        end
    end
    faces = SoftRobots.convex_hull(positions)
    robot = SoftRobots.SoftRobot(nodes, edges, faces)
    velocities = [SoftRobots.Point{3, Float64}(0) for n in robot.nodes]
    state = SoftRobots.SoftRobotState(positions, velocities)
    robot, state
end

robot, state = tower()
SoftRobots.draw(robot, state)
point = [0.25; 0.0]

function potential_energy_fixed_base{T}(robot, state::SoftRobots.SoftRobotState{T})
    energy = zero(T)
    for edge in robot.edges
        energy += 0.5 * edge.stiffness *
         (norm(state.positions[edge.parent] - state.positions[edge.child]) -
          edge.rest_length)^2
    end
    for position in state.positions[1:4]
        energy += 0.5 * 1000 * position[3]^2
    end
    energy
end

function cost{T}(x::Vector{T})
    state = convert(SoftRobots.SoftRobotState{T}, x)
    energy = potential_energy_fixed_base(robot, state)
    SoftRobots.update_barrier!(robot, state)
#     v = SoftRobots.evaluate(state.barrier, point)
#     if v > 0
#         energy += 0.5 * 1000 * v^2
#     end
    energy
end

x0 = collect(destructure(state.positions))
f = cost
g = ForwardDiff.gradient(f)
x = copy(x0)
energies = []
step = 1e-4
@time for i = 1:1
    gi = g(x)
    x -= step * gi
    draw(robot, map(value, x))
    push!(energies, f(x))
end


text:
# Geometry which displays points at given (x, y) positions.
immutable PointGeometry <: Gadfly.GeometryElement
    tag::Symbol

    function PointGeometry(; tag::Symbol=empty_tag)
        new(tag)
    end
end

const point = PointGeometry

function element_aesthetics(::PointGeometry)
    [:x, :y, :size, :color, :shape]
end

# Generate a form for a point geometry.
#
# Args:
#   geom: point geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::PointGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(RGBA{Float32}[theme.default_color])
    default_aes.size = Measure[theme.default_point_size]
    aes = inherit(aes, default_aes)

    lw_hover_scale = 10
    lw_ratio = theme.line_width / aes.size[1]

    aes_x, aes_y = concretize(aes.x, aes.y)

    ctx = context()
    if aes.shape != nothing
        xs = Array(eltype(aes_x), 0)
        ys = Array(eltype(aes_y), 0)
        cs = Array(eltype(aes.color), 0)
        size = Array(eltype(aes.size), 0)
        shape_max = maximum(aes.shape)
        if shape_max > length(theme.shapes)
            error("Too many values for the shape aesthetic. Define more shapes in Theme.shapes")
        end

        for shape in 1:maximum(aes.shape)
            for (x, y, c, sz, sh) in Compose.cyclezip(aes.x, aes.y, aes.color,
                                                      aes.size, aes.shape)
                if sh == shape
                    push!(xs, x)
                    push!(ys, y)
                    push!(cs, c)
                    push!(size, sz)
                end
            end
            compose!(ctx, (context(), theme.shapes[shape](xs, ys, size), fill(cs)))
            empty!(xs)
            empty!(ys)
            empty!(cs)
            empty!(size)
        end
    else
        compose!(ctx,
            circle(aes.x, aes.y, aes.size, geom.tag),
            fill(aes.color))
    end
    compose!(ctx, linewidth(theme.highlight_width))

    if aes.color_key_continuous != nothing && aes.color_key_continuous
        compose!(ctx,
            stroke(map(theme.continuous_highlight_color, aes.color)))
    else
        stroke_colors =
            Gadfly.pooled_map(RGBA{Float32}, theme.discrete_highlight_color, aes.color)
        classes =
            Gadfly.pooled_map(Compat.ASCIIString,
                c -> svg_color_class_from_label(escape_id(aes.color_label([c])[1])),
                aes.color)

        compose!(ctx, stroke(stroke_colors), svgclass(classes))
    end

    return compose!(context(order=4), svgclass("geometry"), ctx)
end


text:
using MCMC

println("    Testing basic HMC constructors...")

HMC()
HMC(20)
HMC(0.75)
HMC(20, 0.75)
HMC(init=20)
HMC(scale=0.75)
HMC(init=20, scale=0.75)
mctuner = EmpMCTuner(0.85)
HMC(mctuner)
HMC(20, mctuner)
HMC(0.75, mctuner)
HMC(20, 0.75, mctuner)
HMC(init=20, tuner=mctuner)
HMC(scale=0.75, tuner=mctuner)
HMC(init=20, scale=0.75, tuner=mctuner)

println("    Testing that HMC sampler can run...")

mcmodel = model(v->-dot(v,v), grad=v->-2v, init=ones(3))
mcchain = run(mcmodel, HMC(0.75), SerialMC(steps=5000, burnin=1000))
# mcchain = run(mcmodel, HMC(0.75, mctuner), SerialMC(steps=5000, burnin=1000))


text:
using WebIO
using Test

@testset "Scope id deprecations" begin
    @test_deprecated Scope("manual-ids-are-deprecated")
    @test_deprecated Scope(; id="still-deprecated")

    # Make sure that non-deprecated arguments are still handled as expected.
    # This surfaced in https://github.com/JuliaGizmos/Interact.jl/issues/315.
    @test_deprecated !isempty(Scope("foo"; imports=["foo.js"]).imports)
end


text:
# R formula to julia Formula

# a special wrapper for rcopy specialized for formulas
function rcopy_formula(s::Ptr{S}) where S<:Sxp
    r = rcopy(s)
    if isa(r, Number) && isinteger(r)
        Int(r)
    else
        r
    end
end

function rcopy(::Type{Expr}, l::Ptr{LangSxp})
    op = rcopy(Symbol, l[1])
    if op == :(:)
        op = :&
    end
    args = cdr(l)
    if op == Symbol("(")
        f = rcopy(Expr, l[2])
    else
        cargs = []
        for (t, s) in pairs(args)
            if isNull(t)
                push!(cargs, rcopy_formula(s))
            else
                push!(cargs, Expr(:(=), rcopy_formula(t), rcopy_formula(s)))
            end
        end
        f = Expr(:call, op, cargs...)
    end
    # unwind these operators
    if op in (:+, :*, :&) && isa(f.args[2], Expr) && f.args[2].args[1] == op
        f = Expr(:call, op, f.args[2].args[2:end]..., f.args[3])
    end
    f
end

function rcopy(::Type{FormulaTerm}, l::Ptr{LangSxp})
    expr = rcopy(Expr, l)
    @eval StatsModels.@formula($expr)
end

# julia Formula to R formula

# formula
function sexp_formula(e::Expr)
    e.head == :call || error("invalid formula object")
    op = e.args[1]
    if op == :&
        op = :(:)
    end
    if length(e.args) > 3 && op in (:+, :*, :(:))
        rlang_p(op, sexp_formula(Expr(e.head, e.args[1:end-1]...)), sexp_formula(e.args[end]))
    else
        rlang_p(op, map(sexp_formula, e.args[2:end])...)
    end
end
sexp_formula(t::StatsModels.TupleTerm) = sexp_formula(Expr(:call, :+, t...))
sexp_formula(e::Symbol) = sexp(SymSxp, e)
sexp_formula(n::Integer) = sexp(RClass{:numeric}, Float64(n))
sexp_formula(n::Number) = sexp(n)
sexp_formula(t::ConstantTerm) = sexp(float(t.n))
sexp_formula(t::FunctionTerm) = sexp_formula(t.exorig)
sexp_formula(i::InteractionTerm) = sexp_formula(Expr(:call, :(:), i.terms...))
sexp_formula(t::Term) = sexp(SymSxp, t.sym)

# R formula objects
function sexp(::Type{RClass{:formula}}, f::FormulaTerm)
    s = protect(rlang_p(:~, sexp_formula(f.lhs), sexp_formula(f.rhs)))
    try
        setattrib!(s, Const.ClassSymbol, "formula")
        setattrib!(s, ".Environment", Const.GlobalEnv)
    finally
        unprotect(1)
    end
    s
end


text:
function add_note!(a_schedule, song, note, time, key)
    add!(
        a_schedule,
        Map(Scale(song.volume), Map(song.wave, Cycles(key * Rational(note.interval)))),
        time,
        song.make_envelope(note.beats * song.beat_duration)...,
    )
end

function add_chord!(a_schedule, song, chord, time, key)
    modulation = chord.modulation
    key = key * Rational(modulation.interval)
    for note in chord.notes
        add_note!(a_schedule, song, note, time, key)
    end
    time + modulation.beats * song.beat_duration, key
end

function add_chord!(a_schedule, song, chord_vectors::Vector, time, key)
    for chords in chord_vectors
        key, time = add_chord!(a_schedule, song, chords, time, key)
    end
    time, key
end

function AudioSchedule(song::Song; chords = song.chords, initial_key = song.initial_key)
    time = 0s
    a_schedule = AudioSchedule(; sample_rate = song.sample_rate)
    key = initial_key
    for chord in chords
        time, key = add_chord!(a_schedule, song, chord, time, key)
    end
    a_schedule
end

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
module SNN

using Reexport
@reexport using Utils

const Int = Int32
const Float = Float32
# srand(1000)

include("units.jl")
include("main.jl")
include("utils.jl")

include("neuron/if.jl")
include("neuron/if2.jl")
include("neuron/noisy_if.jl")
include("neuron/poisson.jl")
include("neuron/iz.jl")
include("neuron/hh.jl")
include("neuron/rate.jl")

include("synapse/spiking_synapse.jl")
include("synapse/rate_synapse.jl")
include("synapse/pinning_full_synapse.jl")
# include("synapse/pinning_synapse.jl")
include("synapse/fl_full_synapse.jl")
# include("synapse/fl_synapse.jl")

isdefined(Main, :Plots) && include("plot.jl")

end


text:
using Documenter, ASTInterpreter2

makedocs(
    modules = [ASTInterpreter2],
    clean = false,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "ASTInterpreter2.jl",
    authors = "Keno Fischer, Tim Holy, and others",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "ast.md",
        "internals.md",
        "dev_reference.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaDebug/ASTInterpreter2.jl.git",
)


text:
module RunMetrics

export import_strava
export import_garmin

include("run_summary.jl")
include("tcx.jl")
include("stats.jl")
include("stats2.jl")
include("file_imports.jl")

end # module


text:
import REPL
using REPL.TerminalMenus
using Plots

options = ["sin", "cos", "sinh", "cosh", "exp"]
menu = MultiSelectMenu(options)

function main()
    choices = request("Select the function you like to plot:", menu)
    p = plot()
    if length(choices) > 0
        println("plot")
        x = -π:0.01:π
        for i in choices
            fname = options[i]
            f = eval(Meta.parse(fname))
            println("func name = $f")
            plot!(p, f, x)
        end
    else
        println("Menu canceled")
    end
    return p
end

main() |> display

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Copyright (c) 2016, Jernej Kovacic
# Licensed under the MIT "Expat" License.
#
# See LICENSE.md for more details.

"""
Exception thrown when any input argument is not valid.
"""
type SpecfunUndefinedError <: Exception end

"""
Exception thrown when the internal algorithm does not converge.
"""
type NoConvergenceError <: Exception end


text:
###############################################################################
#
# Plant Hydraulic system
#
###############################################################################
"""
    abstract type AbstractPlantOrganism{FT}

Hierachy of AbstractPlantOrganism
- [`GrassLikeOrganism`](@ref)
- [`PalmLikeOrganism`](@ref)
- [`TreeLikeOrganism`](@ref)
"""
abstract type AbstractPlantOrganism{FT<:AbstractFloat} end

"""
    mutable struct GrassLikeOrganism{FT}

A plant hydraulic system like a grass, which contains multiple root layers, and
multiple canopy layers. No trunk or branch system applies.

# Fields
$(TYPEDFIELDS)
"""
mutable struct GrassLikeOrganism{FT} <: AbstractPlantOrganism{FT}
    # structre information
    "Root Layers"
    n_root  ::Int
    "Canopy Layers"
    n_canopy::Int

    # Arrays of roots and leaves
    "Roots system"
    roots ::Array{RootHydraulics{FT},1}
    "Leaves"
    leaves::Array{LeafHydraulics{FT},1}

    # Root and canopy index in Soil and Atmosphere
    "Corresponding soil layer per root layer"
    root_index_in_soil ::Array{Int,1}
    "Corresponding air layer per canopy layer"
    canopy_index_in_air::Array{Int,1}

    # containers for root flows
    "Conductances for each layer at given flow"
    cache_k::Array{FT,1}
    "Pressure for each layer at given flow"
    cache_p::Array{FT,1}
    "Flow rate"
    cache_q::Array{FT,1}
end

"""
    mutable struct PalmLikeOrganism{FT}

A plant hydraulic system like a palm, which contains multiple root layers, one
trunk, and multiple canopy layers. No branch system applies.

# Fields
$(TYPEDFIELDS)
"""
mutable struct PalmLikeOrganism{FT} <: AbstractPlantOrganism{FT}
    # structre information
    "Root Layers"
    n_root  ::Int
    "Canopy Layers"
    n_canopy::Int

    # Arrays of roots and leaves
    "Roots system"
    roots ::Array{RootHydraulics{FT},1}
    "Trunk"
    trunk ::StemHydraulics{FT}
    "Leaves"
    leaves::Array{LeafHydraulics{FT},1}

    # Root and canopy index in Soil and Atmosphere
    "Corresponding soil layer per root layer"
    root_index_in_soil ::Array{Int,1}
    "Corresponding air layer per canopy layer"
    canopy_index_in_air::Array{Int,1}

    # containers for root flows
    "Conductances for each layer at given flow"
    cache_k::Array{FT,1}
    "Pressure for each layer at given flow"
    cache_p::Array{FT,1}
    "Flow rate"
    cache_q::Array{FT,1}
end

"""
    mutable struct TreeLikeOrganism{FT}

A plant hydraulic system like a tree, which contains multiple root layers, one
trunk, and multiple branch and canopy layers.

# Fields
$(TYPEDFIELDS)
"""
mutable struct TreeLikeOrganism{FT} <: AbstractPlantOrganism{FT}
    # structre information
    "Root Layers"
    n_root  ::Int
    "Canopy Layers"
    n_canopy::Int

    # Arrays of roots and leaves
    "Roots system"
    roots ::Array{RootHydraulics{FT},1}
    "Trunk"
    trunk ::StemHydraulics{FT}
    "Branch system"
    branch::Array{StemHydraulics{FT},1}
    "Leaves"
    leaves::Array{LeafHydraulics{FT},1}

    # Root and canopy index in Soil and Atmosphere
    "Corresponding soil layer per root layer"
    root_index_in_soil ::Array{Int,1}
    "Corresponding air layer per canopy layer"
    canopy_index_in_air::Array{Int,1}

    # containers for root flows
    "Conductances for each layer at given flow"
    cache_k::Array{FT,1}
    "Pressure for each layer at given flow"
    cache_p::Array{FT,1}
    "Flow rate"
    cache_q::Array{FT,1}
end

"""
    mutable struct TreeSimple{FT}

A plant hydraulic system with one root, one stem, and one leaf for testing
    purpose

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct TreeSimple{FT} <: AbstractPlantOrganism{FT}
    # Arrays of roots and leaves
    "Root"
    root::RootHydraulics{FT} = RootHydraulics{FT}()
    "Stem"
    stem::StemHydraulics{FT} = StemHydraulics{FT}()
    "Leaf"
    leaf::LeafHydraulics{FT} = LeafHydraulics{FT}()

    # Local container for tree information
    "Relative hydraulic conductance"
    krs::Array{FT,1} = ones(FT,4)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
################################################################################
#
# NfOrd/Ideal/Arithmetic.jl : Arithmetic for ideals in orders of absolute
#                             number fields
#
# This file is part of Hecke.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann
#  Copyright (C) 2015, 2016, 2017 Claus Fieker
#
################################################################################

################################################################################
#
#  Ideal addition / GCD
#
################################################################################

doc"""
***
    +(x::NfOrdIdl, y::NfOrdIdl)

> Returns $x + y$.
"""
function +(x::NfOrdIdl, y::NfOrdIdl)
  d = degree(order(x))
  H = vcat(basis_mat(x), basis_mat(y))
  g = gcd(minimum(x), minimum(y))
  if isone(g)
    return ideal(order(x), g)
  end
  H = sub(_hnf_modular_eldiv(H, g, :lowerleft), (d + 1):2*d, 1:d)
  return ideal(order(x), H)
end

################################################################################
#
#  Intersection / LCM
#
################################################################################

doc"""
    intersection(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl
    lcm(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl

> Returns $x \cap y$.
"""
function intersection(x::NfOrdIdl, y::NfOrdIdl)
  d = degree(order(x))
  H = vcat(basis_mat(x), basis_mat(y))
  K = _kernel(H)
  g = lcm(minimum(x),minimum(y))
  return ideal(order(x), _hnf_modular_eldiv(sub(K, 1:d, 1:d)*basis_mat(x), g, :lowerleft))
end

doc"""
    intersection(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl
    lcm(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl

> Returns $x \cap y$.
"""
lcm(x::NfOrdIdl, y::NfOrdIdl) = intersection(x, y)

################################################################################
#
#  Ideal multiplication
#
################################################################################

doc"""
    *(x::NfOrdIdl, y::NfOrdIdl)

> Returns $x \cdot y$.
"""
function *(x::NfOrdIdl, y::NfOrdIdl)
  if ismaximal_known(order(x)) && ismaximal(order(x))
    return mul_maximal(x, y)
  else
    return mul_gen(x, y)
  end
end

function mul_gen(x::NfOrdIdl, y::NfOrdIdl)
  order(x) != order(y) && error("Not compatible")
  O = order(x)
  d = degree(O)
  l = minimum(x, Val{false})*minimum(y, Val{false})
  z = fmpz_mat(degree(O)*degree(O), degree(O))
  z.base_ring = FlintZZ
  X = basis(x, Val{false})
  Y = basis(y, Val{false})
  t = O()
  for i in 1:d
    for j in 1:d
      mul!(t, X[i], Y[j])
      assure_has_coord(t)
      for k in 1:d
        z[(i - 1)*d + j, k] = t.elem_in_basis[k]
      end
    end
  end
  # This is a d^2 x d matrix
  return ideal(O, sub(_hnf_modular_eldiv(z, l, :lowerleft),
                      (d*(d - 1) + 1):d^2, 1:d))
end

# using the 2-normal representation
function prod_via_2_elem_normal(a::NfOrdIdl, b::NfOrdIdl)
  @hassert :NfOrd 1 has_2_elem_normal(a)
  @hassert :NfOrd 1 has_2_elem_normal(b)
  O = order(a)
  a1 = a.gen_one
  b1 = b.gen_one
  m = lcm(a1, b1)
  e, f = ppio(m, a1)
  if f == 1
    a2 = a.gen_two
  else
    # we need to become normal at m, we are normal at a
    # higher powers in a are fine
    # CRT: the 2nd gen of a needs to stay the same at a
    # and should be  1 at f

    #a2 = a.gen_two*f + a1^2
    if !false
      mul!(e, a1, a1)
      a2 = f*a.gen_two
      add!(a2, a2, e)
    else
      g, x, y = gcdx(f, e)

      mul!(e, a1, a1)
      mul!(x, x, f)
      a2 = x*a.gen_two
      mul!(y, y, e)
      add!(a2, a2, y)
     end
    # now (a1, a2) should be m-normal for a
  end

  e, f = ppio(m, b1)
  if f == 1
    b2 = b.gen_two
  else

    #b2 = b.gen_two*f + b1^2
    if !false #Carlo's new method
      mul!(e, b1, b1)
      b2 = f*b.gen_two
      add!(b2, b2, e)
     else
      g, x, y = gcdx(f, e)

      #b2 = b.gen_two*f*x + y*b1^2
      mul!(e, b1, b1)
      mul!(x, x, f)
      b2 = x*b.gen_two
      mul!(y, y, e)
      add!(b2, b2, y)
    end
  end
  C = ideal(O, a1*b1, a2*b2)
  C.norm = norm(a) * norm(b)
#
#CF: too expensive, need norm_mod to compute the norm only modulo...
#
#  if C.norm != gcd(C.gen_one^degree(O), FlintZZ(norm(C.gen_two)))
#    println("a:", a)
#    println("b:", b)
#    println("C:", C)
#    @hassert :NfOrd 1 gcd(a1^degree(O), norm(a2)) == norm(a)
#    @hassert :NfOrd 1 gcd(b1^degree(O), norm(b2)) == norm(b)
##    assert(false)
#  end

  if has_minimum(a) && has_minimum(b)
    ma = minimum(a, Val{false})
    mb = minimum(b, Val{false})

    if gcd(ma, mb) == 1
      C.minimum = ma * mb
    end

    # Otherwise, I don't know the
    # correct value.
  end

  C.gens_normal = m
  C.gens_weakly_normal = true #for the 1-ideal: it will not be normal...

  return C
end

# using the 2-weak-normal representation
function prod_via_2_elem_weakly(a::NfOrdIdl, b::NfOrdIdl)
  @hassert :NfOrd 1 has_2_elem(a)
  @hassert :NfOrd 1 has_2_elem(b)

  O = order(a)
  K = nf(O)
  bas = basis(O)
  n = degree(O)

  norm_c = norm(a) * norm(b)        # we ARE in the maximal order case
  norm_int_c = norm_c
  mod_c = fmpz(1)

  if has_minimum(a)
    mod_c *= minimum(a)
  else
    mod_c *= norm(a)
  end

  if has_minimum(b)
    mod_c *= minimum(b)
  else
    mod_c *= norm(b)
  end

  # Is this a good idea? Bad magic constants

  if mod_c > 250
    r = -500:500  # should still have enough elements...
  else
    r = -Int(div(mod_c, 2)):Int(div(mod_c, 2))
  end

  @vprint :NfOrd 1 "a: $a \nb: $b"
  @vprint :NfOrd 1 "using basis: $bas"

  gen = O()
  gen2 = O()
  t = O()
  s = O()
  u = O()

  cnt = 0
  while true
    cnt += 1
    if cnt % 20 == 0
      assure_2_normal(a)
      assure_2_normal(b)
      return a*b
    end
    rand!(t, O, r) # Nemo.rand_into!(bas, r, t)
    r2 = rand(r)
    rand!(s, O, r) # Nemo.rand_into!(bas, r, s)
    r4 = rand(r)
    mul!(t, t, a.gen_two)      # Nemo.mult_into!(t, a.gen_two, t)
    add!(t, t, r2*a.gen_one)   # Nemo.add_into!(t, r2*a.gen_one, t)
    mul!(s, s, b.gen_two)      # Nemo.mult_into!(s, b.gen_two, s)
    add!(s, s, r4*b.gen_one)   # Nemo.add_into!(s, r4*b.gen_one, s)
    mul!(u, s, t)              # Nemo.mult_into!(s, t, u)
    add!(gen, u, gen)          # Nemo.add_into!(u, gen, gen)
#    gen2 += (r1*K(a.gen_two) + r2*a.gen_one) *
#           (r3*K(b.gen_two) + r4*b.gen_one)
    gen = mod(gen, mod_c^2)    # = element_reduce_mod(gen, O, FlintZZ(mod_c)^2)

    if gcd(norm(gen), norm_int_c^2) == norm_int_c # should be ^n, but for
                                                  # the test ^2 is sufficient
      break
    end
  end

  @vprint :NfOrd 1 "prod_via_2_elem: used $cnt tries\n"

  c = ideal(O, norm_int_c, gen)

  c.norm = norm_c

  if has_minimum(a) && has_minimum(b) && gcd(minimum(a), minimum(b)) == 1
    c.minimum = minimum(a) * minimum(b)
                    # otherwise, I don't know the correct value
  end

  c.gens_weakly_normal = 1

  return c
end

# dispatching
doc"""
***
  *(x::NfMaxIdl, y::NfOrdIdl)

> Returns the ideal x*y.
"""
function mul_maximal(x::NfOrdIdl, y::NfOrdIdl)
  if x.iszero == 1 || y.iszero == 1
    z = ideal(order(x), zero_matrix(FlintZZ, degree(order(x)), degree(order(x))))
    z.iszero = 1
    return z
  end
  if isone(x)
    return y
  elseif isone(y)
    return x
  end
  if has_2_elem_normal(x) && has_2_elem_normal(y)
    return prod_via_2_elem_normal(x, y)
  end
  if has_2_elem(x) && has_2_elem(y)
    return prod_via_2_elem_weakly(x, y)
  end
  # Fall back to the generic algorithm _mul(::NfOrdIdl, ::NfOrdIdl)
  # Could also use invoke
  return mul_gen(x, y)
end

################################################################################
#
#  Addition
#
################################################################################

# Falls back to generic case +(::NfOrd, ::NfOrd)
#for ideals in the maximal order, the gcd is well defined...

doc"""
***
  gcd(A::NfOrdIdl, B::NfOrdIdl) -> NfOrdIdl
> The gcd or sum (A+B).
"""
function gcd(A::NfOrdIdl, B::NfOrdIdl)
  return A+B
end

function gcd_into!(A::NfOrdIdl, B::NfOrdIdl, C::NfOrdIdl)
  return C+B
end

#TODO: write a ppio version that allows for p-powers as well
doc"""
***
  gcd(A::NfOrdIdl, p::fmpz) -> NfOrdIdl
> The gcd or sum (A+ pO).
"""
function gcd(A::NfOrdIdl, p::fmpz)
  if isdefined(A, :minimum)
    if gcd(A.minimum, p) == 1
      return ideal(order(A), fmpz(1))
    end
  end
  if isdefined(A, :norm)
    if gcd(A.norm, p) == 1
      return ideal(order(A), fmpz(1))
    end
  end
  if has_2_elem(A)
    g = gcd(A.gen_one, p)
    @assert !isone(g)
    return ideal(order(A), g, A.gen_two)
  else
    @assert isdefined(A, :basis_mat)
    return A + ideal(order(A), p)
  end
end

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

# multiplication by fmpz, using two normal presentation
function prod_by_int_2_elem_normal(A::NfOrdIdl, a::fmpz)
  @assert has_2_elem(A) && has_2_elem_normal(A)

  # <a,a> is a a-normal presentation
  # we need to have a common presentation for A and a though

  m = lcm(a, A.gen_one)

  e, f = ppio(m, A.gen_one)
  if f == 1
    a2 = A.gen_two
  else
                           # we need to become normal at m, we are normal at a
                           # higher powers in a are fine
                           # CRT: the 2nd gen of a needs to stay the same at a
                           # and should be  1 at f
    a2 = A.gen_two*f + A.gen_one^2 # now (a1, a2) should be m-normal for a
  end

  B = NfOrdIdl(A.gen_one*a, a2*a)
  B.gens_normal = m

  if has_minimum(A)
    B.minimum = A.minimum * a
  end

  if has_norm(A)
    B.norm = A.norm * a^degree(order(A))
  end

  @assert has_2_elem(B) && has_2_elem_normal(B)
  return B
end

function prod_by_int_2_elem(A::NfOrdIdl, a::fmpz)
  @assert has_2_elem(A)

  B = NfOrdIdl(A.gen_one*a, A.gen_two*a)
  B.gens_normal = A.gens_normal

  if has_minimum(A)
    B.minimum = A.minimum * a
  end

  if has_norm(A)
    B.norm = A.norm * a^degree(order(A))
  end

  @assert has_2_elem(B)
  return B
end

doc"""
***
    *(x::NfOrdIdl, y::fmpz) -> NfOrdIdl

> Returns the ideal $x \cdot y$.
"""
function *(x::NfOrdIdl, y::fmpz)
  if ismaximal_known(order(x)) && ismaximal(order(x))
    return mul_maximal(x, y)
  else
    return mul_gen(x, y)
  end
end

function mul_maximal(x::NfOrdIdl, y::fmpz)
  if y == 1 || y == -1
    return x
  end

  if has_2_elem(x)
    if has_2_elem_normal(x)
      return prod_by_int_2_elem_normal(x,y)
# this does not make any sense since prod_by_int_2_elem(x, y) works only if has_2_elem_norm(x) == true
#    else
#      return prod_by_int_2_elem(x,y)
    end
  end

  return x*ideal(order(x), y)
end

*(x::fmpz, y::NfOrdIdl) = y * x

doc"""
***
    *(x::NfOrdIdl, y::Integer) -> NfOrdIdl

> Returns the ideal $x \cdot y$.
"""
*(x::NfOrdIdl, y::Integer) = x * fmpz(y)

*(x::Integer, y::NfOrdIdl) = y * x

function *(x::NfOrdElem, y::NfOrdIdl)
  parent(x) != order(y) && error("Orders of element and ideal must be equal")
  return ideal(parent(x), x) * y
end

*(x::NfOrdIdl, y::NfOrdElem) = y * x

function mul_gen(x::NfOrdIdl, y::fmpz)
  z = ideal(order(x), basis_mat(x)*y)
  if isdefined(x, :princ_gen)
    z.princ_gen = x.princ_gen * y
  end
  if isdefined(x, :princ_gen_special)
    z.princ_gen_special = (2, 0, x.princ_gen_special[x.princ_gen_special[1] + 1] * y)
  end
  return z
end

################################################################################
#
#  Idempotents
#
################################################################################

doc"""
    idempotents(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdElem, NfOrdElem

> Returns a tuple `(e, f)` consisting of elements `e in x`, `f in y` such that
> `1 = e + f`.
>
> If the ideals are not coprime, an error is raised.
"""
function idempotents(x::NfOrdIdl, y::NfOrdIdl)
  !(order(x) === order(y)) && error("Parent mismatch")
  O = order(x)
  d = degree(O)

  mx = minimum(x)
  my = minimum(y)

  g, u, v = gcdx(mx, my)
  if isone(g)
    z = O(u*mx)
    @hassert :NfOrd 2 z in x
    @hassert :NfOrd 2 (1 - z) in y
    return z, 1 - z
  end

  # form the matrix
  #
  # ( 1 |  1  | 0 )
  # ( 0 | M_x | I )
  # ( 0 | M_y | 0 )

  V = O.tidempotents

  u = elem_in_basis(one(O))

  V[1, 1] = 1

  for i in 1:d
    V[1, i + 1] = u[i]
  end

  _copy_matrix_into_matrix(V, 2, 2, basis_mat(x))
  _copy_matrix_into_matrix(V, 2 + d, 2, basis_mat(y))

  for i in 1:d
    V[1 + i, d + 1 + i] = 1
  end

  H = hnf(V) # upper right

  for i in 2:(1 + d)
    if H[1, i] != 0
      error("Ideals are not coprime")
    end
  end

  basisofx = basis(x, Val{false})

  z = basisofx[1]*H[1, d + 2]

  for i in 2:d
    z = z + basisofx[i]*H[1, d + 1 + i]
  end

  @hassert :NfOrd 2 -z in x
  @hassert :NfOrd 2 1 + z in y

  ccall((:fmpz_mat_zero, :libflint), Void, (Ptr{fmpz_mat}, ), &V)

  return -z, 1 + z
end

################################################################################
#
#  crt
#
################################################################################
doc"""
   crt(r1::NfOrdElem, i1::NfOrdIdl, r2::NfOrdElem, i2::NfOrdIdl) -> NfOrdElem
> Find $x$ s.th $x \equiv r1 \bmod i1$ and $x \equiv r2 \bmod i2$
using (((idempotents)))
"""
function crt(r1::NfOrdElem, i1::NfOrdIdl, r2::NfOrdElem, i2::NfOrdIdl)
  u, v = idempotents(i1, i2)
  return r1*v + r2*u
end

function crt(a::Array{NfOrdElem, 1}, I::Array{NfOrdIdl, 1})
  if length(a) == 1
    return a[1]
  end
  if length(a) == 2
    return crt(a[1], I[1], a[2], I[2])
  end
  A = [crt(a[2*i-1], I[2*i-1], a[2*i], I[2*i]) for i=1:div(length(a), 2)]
  B = [I[2*i-1]*I[2*i] for i=1:div(length(a), 2)]
  if isodd(length(a))
    push!(A, a[end])
    push!(B, I[end])
  end
  return crt(A, B)
end
################################################################################
#
#  Division
#
################################################################################

divexact(A::NfOrdIdl, b::Integer) = divexact(A, fmpz(b))
doc"""
***
    divexact(A::NfOrdIdl, y::fmpz) -> NfOrdIdl
    divexact(A::NfOrdIdl, y::Integer) -> NfOrdIdl

> Returns $A/y$ assuming that $A/y$ is again an integral ideal.
"""
#TODO: write a divexact! to change the ideal?
#  difficult due to Julia's inability to unset entries...
function divexact(A::NfOrdIdl, b::fmpz)
  zk = order(A)
  if has_2_elem(A)
    B = ideal(zk, divexact(A.gen_one, b), divexact(A.gen_two, b))
    if isdefined(A, :gens_normal)
      B.gens_normal = A.gens_normal
    end
    B.gens_weakly_normal = A.gens_weakly_normal
    if has_basis_mat(A)
      B.basis_mat = divexact(A.basis_mat, b)
    end
    if false && has_basis_mat_inv(A)
      error("not defined at all")
      B.basis_mat_inv = b*A.basis_mat_inv
    end
  else
    B = ideal(zk, divexact(A.basis_mat, b))
    if false && has_basis_mat_inv(A)
      error("not defined at all")
      B.basis_mat_inv = b*A.basis_mat_inv
    end
  end
  if has_minimum(A)
    B.minimum = divexact(A.minimum, b)
  end
  if has_norm(A)
    B.norm = divexact(A.norm, b^degree(zk))
  end
  if has_princ_gen(A)
    B.princ_gen = divexact(A.princ_gen, b)
  end
  #TODO princ_gen_special missing
  return B
end

doc"""
***
    divexact(A::NfOrdIdl, B::NfOrdIdl) -> NfOrdIdl

> Returns $AB^{-1}$ assuming that $AB^{-1}$ is again an integral ideal.
"""
function divexact(A::NfOrdIdl, B::NfOrdIdl)
  # It is assumed that B divides A, that is, A \subseteq B
  t_prod = 0.0
  t_simpl = 0.0
  t_b_mat = 0.0
  t_2_elem_weak = 0.0
  t_2_elem = 0.0

  if norm(A) == norm(B)
    return ideal(order(A), one(FlintZZ), order(A)(1))
  else
    t_prod += @elapsed I = A*inv(B)
    t_simpl += @elapsed simplify_exact!(I)
    #t_b_mat += @elapsed B = basis_mat(I)
    I.den != 1 && error("Division not exact")
    #I = ideal(order(A), B.num)
    #t_2_elem_weak += @elapsed _assure_weakly_normal_presentation(I)
    #t_2_elem += @elapsed assure_2_normal(I)

    #println("  computation for product: $t_prod")
    #println("  simplification         : $t_simpl")
    #println("  basis matrix           : $t_b_mat")
    #println("  2 weak presentation    : $t_2_elem_weak")
    #println("  2 elem presentation    : $t_2_elem")

    return I.num

  end
end