text:
Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
Centroid_ODE(z,w)=conj( 1 ./(1+z.^2)+1 ./(1+w.^2) )
Ham_numerator(Q,P)=((Q[1]+Q[2])^2+(P[1]+P[2])^2)*((Q[1]-Q[2])^2+(P[1]-P[2])^2)
Ham_denominator(Q,P)=(Q[1]^2+(P[2]-1)^2)*(Q[1]^2+(P[2]+1)^2)*(Q[2]^2+(P[1]-1)^2)*(Q[2]^2+(P[1]+1)^2)
Aref_Hamiltonian(Q,P,t)=1/2*sqrt(Ham_numerator(Q,P)/Ham_denominator(Q,P));
Hamiltonian_Dimer(q,p,t)=4sqrt(abs( (q[1]^2+p[1]^2)*(q[2]^2+p[2]^2) )/
(
        (2-2sqrt(2)*(p[1]+p[2])+(p[1]+p[2])^2+(q[1]-q[2])^2 )*
        (2+2sqrt(2)*(p[1]+p[2])+(p[1]+p[2])^2+(q[1]-q[2])^2 )*
        (2-2sqrt(2)*(p[1]-p[2])+(p[1]-p[2])^2+(q[1]+q[2])^2 )*
        (2+2sqrt(2)*(p[1]-p[2])+(p[1]-p[2])^2+(q[1]+q[2])^2 )

));
Hamiltonian_Dimer_Backwards(q,p,t)=-4sqrt(abs( (q[1]^2+p[1]^2)*(q[2]^2+p[2]^2) )/
(
        (2-2sqrt(2)*(p[1]+p[2])+(p[1]+p[2])^2+(q[1]-q[2])^2 )*
        (2+2sqrt(2)*(p[1]+p[2])+(p[1]+p[2])^2+(q[1]-q[2])^2 )*
        (2-2sqrt(2)*(p[1]-p[2])+(p[1]-p[2])^2+(q[1]+q[2])^2 )*
        (2+2sqrt(2)*(p[1]-p[2])+(p[1]-p[2])^2+(q[1]+q[2])^2 )

));
# function Eq_of_M(du,u,p,t)
#     du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[5]=0
#     return
# end

# function Eq_of_M_Lagrangian_Descriptors(du,u,p,t)
#     du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[5]=sqrt(abs(du[1])^2+ abs(du[2])^2+ abs(du[3])^2+ abs(du[4])^2)
#     return
# end

function Eq_of_M_Lagrangian_Descriptors(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=abs(du[1])+ abs(du[2])+ abs(du[3])+ abs(du[4])
    return
end

# function Eq_of_M_Lagrangian_Descriptors(du,u,p,t)
#     du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[5]=(sqrt(abs(du[1]))+ sqrt(abs(du[2]))+ sqrt(abs(du[3]))+ sqrt(abs(du[4])))^2
# end

# function Eq_of_M_Lagrangian_Descriptors_Backwards(du,u,p,t)
#     du[1]=-real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[2]=-imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[3]=-real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[4]=-imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[5]=(sqrt(abs(du[1]))+ sqrt(abs(du[2]))+ sqrt(abs(du[3]))+ sqrt(abs(du[4])))^2
#     return
# end

# function Eq_of_M_Lagrangian_Descriptors_Backwards(du,u,p,t)
#     du[1]=-real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[2]=-imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[3]=-real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[4]=-imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[5]=sqrt(abs(du[1])^2+ abs(du[2])^2+ abs(du[3])^2+ abs(du[4])^2)
#     return
# end

# function Eq_of_M_Lagrangian_Descriptors_Backwards(du,u,p,t)
#     du[1]=-real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[2]=-imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
#     du[3]=-real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[4]=-imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
#     du[5]=sqrt(abs(du[1])^2+ abs(du[2])^2+ abs(du[3])^2+ abs(du[4])^2)
#     return
# end

function Eq_of_M_Lagrangian_Descriptors_Backwards(du,u,p,t)
    du[1]=-real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=-imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=-real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=-imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=abs(du[1])+ abs(du[2])+ abs(du[3])+ abs(du[4])
    return
end

function Eq_of_M_LAB_FRAME(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=real(Centroid_ODE(u[1]+im*u[2],u[3]+im*u[4]))
    du[6]=imag(Centroid_ODE(u[1]+im*u[2],u[3]+im*u[4]))
    return
end;

function Yfind(Q,P,H)
    try
        Y=sqrt( (1+(-1)*H*((-1)+P^2)^2)^(-1)*(P^2+(-1)*Q^2+H*((-1)+P^2)
        ^2*((-1)+Q^2)+sqrt(H*((-1)+P^2)^4+(-4)*((-1)+H*((-1)+P^2)^2)
        *((-1)*P^2+H*((-1)+P^2)^2)*Q^2 ))   )
     catch
        Y=zeros(0)
    end
end
function Yfind_Aref(Q,P,H_Aref)
    H=(2*H_Aref)^2
    try
        Y=sqrt( (1+(-1)*H*((-1)+P^2)^2)^(-1)*(P^2+(-1)*Q^2+H*((-1)+P^2)
        ^2*((-1)+Q^2)+sqrt(H*((-1)+P^2)^4+(-4)*((-1)+H*((-1)+P^2)^2)
        *((-1)*P^2+H*((-1)+P^2)^2)*Q^2 ))   )
     catch
        Y=zeros(0)
    end
end
using Roots
function Q1_find_dimer(Q2,P2,H_dimer)
    Q1_to_find(Q1)=Hamiltonian_Dimer([Q1,Q2],[0,P2],1)-H_dimer
    try
        Q1=find_zeros(Q1_to_find,0,4,maxeval=100,maxfnevals=300,tol=1e-15)
     catch
        Q1=zeros(0)
    end
end
function P1_find_dimer_second(Q2,P2,H_dimer)
    P1_to_find(P1)=Hamiltonian_Dimer([0,Q2],[P1,P2],1)-H_dimer
    try
        P1=find_zero(P1_to_find,1.4,maxeval=10000,maxfnevals=300,tol=1e-15)
     catch
        P1=zeros(0)
    end
end
function P1_find_dimer(Q,P,H)
    coeff=(-16)*H^2 + 32*H^2*P^2 + (-24)*H^2*P^4 +8*H^2*P^6 + (-1)*H^2*P^8 + (-32)*H^2*Q^2 + 16*H^2*P^2*Q^2 + 8*H^2*P^4*Q^2 + (-4)*H^2*P^6*Q^2 + (-24)*H^2*Q^4 + (-8)*H^2*P^2*Q^4 + (-6)*H^2*P^4*Q^4 + (-8)*H^2*Q^6 + (-4)*H^2*P^2*Q^6 + (-1)*H^2*Q^8
    coeff2=(32*H^2 + 16*P^2 + (-16)*H^2*P^2 + (-8)*H^2*P^4 + 4*H^2*P^6 + 16*Q^2 + 16*H^2*Q^2 + (-80)*H^2*P^2*Q^2 + 4*H^2*P^4*Q^2 + (-8)*H^2*Q^4 + (-4)*H^2*P^2*Q^4 + (-4)*H^2*Q^6)
    coeff4=((-24)*H^2 + (-8)*H^2*P^2 + (-6)*H^2*P^4 + 8*H^2*Q^2 +  4*H^2*P^2*Q^2 + (-6)*H^2*Q^4)
    coeff6=(8*H^2 + 4*H^2*P^2 + (-4)*H^2*Q^2)
    coeff8=(-1)*H^2;
    poynomial4=[coeff,coeff2,coeff4,coeff6,coeff8];
    p1=[]
    sol=roots(poynomial4,polish=true)
    for k in 1:length(sol)
        if abs(imag(sol[k]))<1e-8
            push!(p1,sqrt(abs(sol[k])))
        end
    end
    if length(p1)==2
        return p1[2]
    else
        return zeros(0)
    end
end

function P1_poly(Q,P,H)
    coeff=(-16)*H^2 + 32*H^2*P^2 + (-24)*H^2*P^4 +8*H^2*P^6 + (-1)*H^2*P^8 + (-32)*H^2*Q^2 + 16*H^2*P^2*Q^2 + 8*H^2*P^4*Q^2 + (-4)*H^2*P^6*Q^2 + (-24)*H^2*Q^4 + (-8)*H^2*P^2*Q^4 + (-6)*H^2*P^4*Q^4 + (-8)*H^2*Q^6 + (-4)*H^2*P^2*Q^6 + (-1)*H^2*Q^8
    coeff2=(32*H^2 + 16*P^2 + (-16)*H^2*P^2 + (-8)*H^2*P^4 + 4*H^2*P^6 + 16*Q^2 + 16*H^2*Q^2 + (-80)*H^2*P^2*Q^2 + 4*H^2*P^4*Q^2 + (-8)*H^2*Q^4 + (-4)*H^2*P^2*Q^4 + (-4)*H^2*Q^6)
    coeff4=((-24)*H^2 + (-8)*H^2*P^2 + (-6)*H^2*P^4 + 8*H^2*Q^2 +  4*H^2*P^2*Q^2 + (-6)*H^2*Q^4)
    coeff6=(8*H^2 + 4*H^2*P^2 + (-4)*H^2*Q^2)
    coeff8=(-1)*H^2;
    poynomial4=[coeff,coeff2,coeff4,coeff6,coeff8];
    p1=[]
    sol=roots(poynomial4,polish=true)
    for k in 1:length(sol)
        if abs(imag(sol[k]))<1e-8
            push!(p1,sqrt(abs(sol[k])))
        end
    end
    if length(p1)==2
        return minimum(p1)
    else
        return zeros(0)
    end
end
function P1_poly_min(Q,P,H)
    coeff=(-16)*H^2 + 32*H^2*P^2 + (-24)*H^2*P^4 +8*H^2*P^6 + (-1)*H^2*P^8 + (-32)*H^2*Q^2 + 16*H^2*P^2*Q^2 + 8*H^2*P^4*Q^2 + (-4)*H^2*P^6*Q^2 + (-24)*H^2*Q^4 + (-8)*H^2*P^2*Q^4 + (-6)*H^2*P^4*Q^4 + (-8)*H^2*Q^6 + (-4)*H^2*P^2*Q^6 + (-1)*H^2*Q^8
    coeff2=(32*H^2 + 16*P^2 + (-16)*H^2*P^2 + (-8)*H^2*P^4 + 4*H^2*P^6 + 16*Q^2 + 16*H^2*Q^2 + (-80)*H^2*P^2*Q^2 + 4*H^2*P^4*Q^2 + (-8)*H^2*Q^4 + (-4)*H^2*P^2*Q^4 + (-4)*H^2*Q^6)
    coeff4=((-24)*H^2 + (-8)*H^2*P^2 + (-6)*H^2*P^4 + 8*H^2*Q^2 +  4*H^2*P^2*Q^2 + (-6)*H^2*Q^4)
    coeff6=(8*H^2 + 4*H^2*P^2 + (-4)*H^2*Q^2)
    coeff8=(-1)*H^2;
    poynomial4=[coeff,coeff2,coeff4,coeff6,coeff8];
    p1=[]
    sol=roots(poynomial4,polish=true)
    for k in 1:length(sol)
        if abs(imag(sol[k]))<1e-5
            push!(p1,sqrt(abs(sol[k])))
        end
    end
    if !isempty(p1)
        return minimum(p1)
    else
        return zeros(0)
    end
end
function P1_poly_max(Q,P,H)
    coeff=(-16)*H^2 + 32*H^2*P^2 + (-24)*H^2*P^4 +8*H^2*P^6 + (-1)*H^2*P^8 + (-32)*H^2*Q^2 + 16*H^2*P^2*Q^2 + 8*H^2*P^4*Q^2 + (-4)*H^2*P^6*Q^2 + (-24)*H^2*Q^4 + (-8)*H^2*P^2*Q^4 + (-6)*H^2*P^4*Q^4 + (-8)*H^2*Q^6 + (-4)*H^2*P^2*Q^6 + (-1)*H^2*Q^8
    coeff2=(32*H^2 + 16*P^2 + (-16)*H^2*P^2 + (-8)*H^2*P^4 + 4*H^2*P^6 + 16*Q^2 + 16*H^2*Q^2 + (-80)*H^2*P^2*Q^2 + 4*H^2*P^4*Q^2 + (-8)*H^2*Q^4 + (-4)*H^2*P^2*Q^4 + (-4)*H^2*Q^6)
    coeff4=((-24)*H^2 + (-8)*H^2*P^2 + (-6)*H^2*P^4 + 8*H^2*Q^2 +  4*H^2*P^2*Q^2 + (-6)*H^2*Q^4)
    coeff6=(8*H^2 + 4*H^2*P^2 + (-4)*H^2*Q^2)
    coeff8=(-1)*H^2;
    poynomial4=[coeff,coeff2,coeff4,coeff6,coeff8];
    p1=[]
    sol=roots(poynomial4,polish=true)
    for k in 1:length(sol)
        if abs(imag(sol[k]))<1e-5
            push!(p1,sqrt(abs(sol[k])))
        end
    end
    if !isempty(p1)
        return maximum(p1)
    else
        return zeros(0)
    end
end
function Dimer_Eq_of_M(du, u,p, t)
    p=2
    du[1]=-((-4)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4]) +(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(   2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+  (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)  *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[3]^9+24*u[3]^6*   u[4]*u[1]*u[2]+4*u[3]^7*((-4)+(-2)*u[4]^2+3*u[1]^2+2*u[2]^2)+8*   u[3]^4*u[4]*u[1]*u[2]*(4+(-2)*u[4]^2+7*u[1]^2+2*u[2]^2)+4*u[3]^3*   u[1]^2*(12+3*u[4]^4+3*u[1]^4+(-4)*u[2]^2+3*u[2]^4+2*u[1]^2*((-2)+   u[2]^2)+(-2)*u[4]^2*((-2)+u[1]^2+u[2]^2))+8*u[4]*u[1]^3*u[2]*(u[4]^4+(   u[1]^2+(-1)*u[2]^2)^2+2*u[4]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+u[1]^2+   u[2]^2))+2*u[3]^5*(12+3*u[4]^4+9*u[1]^4+(-4)*u[2]^2+3*u[2]^4+8*   u[1]^2*((-2)+u[2]^2)+(-2)*u[4]^2*((-2)+4*u[1]^2+u[2]^2))+(-8)*   u[3]^2*u[4]*u[1]*u[2]*(u[4]^4+(-5)*u[1]^4+u[2]^4+(-8)*u[1]^2*(2+u[2]^2)   +(-4)*(3+u[2]^2)+2*u[4]^2*(2+4*u[1]^2+u[2]^2))+u[3]*((-1)*u[4]^8+   3*u[1]^8+8*u[1]^2*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*(2+u[2]^2)^4+(   -4)*u[4]^6*((-2)+2*u[1]^2+u[2]^2)+(-2)*u[1]^4*(20+36*u[2]^2+5*   u[2]^4)+(-2)*u[4]^4*(12+5*u[1]^4+(-4)*u[2]^2+3*u[2]^4+4*u[1]^2*((   -2)+u[2]^2))+4*u[4]^2*((-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+u[1]^4*(   18+23*u[2]^2)+2*u[1]^2*(4+20*u[2]^2+u[2]^4)))));
    du[2]=-(4*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*  2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+  u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2])   ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+   (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)   *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*((-3)*u[4]^9+(-24)*   u[3]*u[4]^6*u[1]*u[2]+8*u[3]*u[4]^4*u[1]*u[2]*((-4)+2*u[3]^2+(-2)*   u[1]^2+(-7)*u[2]^2)+4*u[4]^7*(4+2*u[3]^2+(-2)*u[1]^2+(-3)*u[2]^2)+  u[4]*((u[3]^4+2*u[3]^2*((-2)+u[1]^2)+(2+u[1]^2)^2)^2+8*(u[3]^6+  u[3]^4*((-2)+u[1]^2)+(-1)*((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[3]^2*  (4+20*u[1]^2+u[1]^4))*u[2]^2+2*(20+5*u[3]^4+36*u[1]^2+5*u[1]^4+(-2)*u[3]^2*(18+23*u[1]^2))*u[2]^4+(-3)*u[2]^8)+(-4)*u[4]^3*  u[2]^2*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+2*((-2)+u[1]^2)*u[2]^2+   3*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+u[2]^2))+(-8)*u[3]*u[1]*u[2]^3*(   u[3]^4+(u[1]^2+(-1)*u[2]^2)^2+2*u[3]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+   u[1]^2+u[2]^2))+(-2)*u[4]^5*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+8*((  -2)+u[1]^2)*u[2]^2+9*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+4*u[2]^2))+8*  u[3]*u[4]^2*u[1]*u[2]*((-12)+u[3]^4+(-4)*u[1]^2+u[1]^4+(-8)*(2+u[1]^2)   *u[2]^2+(-5)*u[2]^4+2*u[3]^2*(2+u[1]^2+4*u[2]^2))));
    du[3]=-(4*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*( 2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+( -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)*((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[1]^9+24*u[3]*u[4]*u[1]^6*u[2]+4*u[1]^7*(4+3*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+8*u[3]*u[4]* u[1]^4*u[2]*((-4)+7*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+4*u[3]^2*u[1]^3*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[2]^2))+2*u[1]^5*(12+9*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+8*u[3]^2*(2+u[4]^2+(-1)* u[2]^2))+8*u[3]*u[4]*u[1]^2*u[2]*(5*u[3]^4+8*u[3]^2*((-2)+u[4]^2+(-1) *u[2]^2)+4*(3+(-1)*u[4]^2+u[2]^2)+(-1)*(u[4]^2+u[2]^2)^2)+8*u[3]^3*u[4]*u[2]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[2]^2)+2*u[3]^2*(2+( -1)*u[4]^2+u[2]^2)+(u[4]^2+u[2]^2)^2)+u[1]*(3*u[3]^8+(-1)*(u[4]^4+2* u[4]^2*((-2)+u[2]^2)+(2+u[2]^2)^2)^2+(-2)*u[3]^4*(20+5*u[4]^4+ 36*u[2]^2+5*u[2]^4+(-2)*u[4]^2*(18+23*u[2]^2))+8*u[3]^2*(u[4]^6+ u[4]^4*((-2)+u[2]^2)+(-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*u[4]^2*(4+20*u[2]^2+u[2]^4)))));
    du[4]=-((-4)*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*((-8)*u[3]*u[4]^3*u[1]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[1]^2)+2*u[3]^2*(2+(-1)*u[4]^2+u[1]^2)+(u[4]^2+u[1]^2)^2)+(u[3]^8+(-3)*u[4]^8+8*u[4]^2*((-2)+u[1]^2)*(2+u[1]^2)^2+(2+u[1]^2)^4+4*u[3]^6*((-2)+(-2)*u[4]^2+u[1]^2)+2*u[3]^4*(12+8*u[4]^2+5*u[4]^4+(-4)*(1+u[4]^2)*u[1]^2+3*u[1]^4)+2*u[4]^4*(20+36*u[1]^2+5*u[1]^4)+4*u[3]^2*(((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[4]^4*(18+23*u[1]^2)+2*u[4]^2*(4+20*u[1]^2+u[1]^4)))*u[2]+8*u[3]*u[4]*u[1]*(u[3]^4+(-5)*u[4]^4+u[1]^4+8*u[4]^2*(2+u[1]^2)+(-4)*(3+u[1]^2)+2*u[3]^2*(2+(-4)*u[4]^2+u[1]^2))*u[2]^2+( -4)*u[4]^2*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[1]^2+3*u[1]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[1]^2))*u[2]^3+(-8)*u[3]*u[4]* u[1]*(2*u[3]^2+7*u[4]^2+(-2)*(2+u[1]^2))*u[2]^4+(-2)*(12+3*u[3]^4+ 9*u[4]^4+(-4)*u[1]^2+3*u[1]^4+u[3]^2*(4+8*u[4]^2+(-2)*u[1]^2)+(-8) *u[4]^2*((-2)+u[1]^2))*u[2]^5+(-24)*u[3]*u[4]*u[1]*u[2]^6+(-4)*(4+ 2*u[3]^2+3*u[4]^2+(-2)*u[1]^2)*u[2]^7+(-3)*u[2]^9)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+ (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+( u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2) *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2));
    du[5]=(abs(du[1])^(1/p))+ (abs(du[2])^(1/p))+ (abs(du[3])^(1/p))+ (abs(du[4])^(1/p)) ;
    du[6]=0;
end

# function Dimer_Eq_of_M(du, u,p, t)
#     du[1]=-((-4)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4]) +(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(   2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+  (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)  *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[3]^9+24*u[3]^6*   u[4]*u[1]*u[2]+4*u[3]^7*((-4)+(-2)*u[4]^2+3*u[1]^2+2*u[2]^2)+8*   u[3]^4*u[4]*u[1]*u[2]*(4+(-2)*u[4]^2+7*u[1]^2+2*u[2]^2)+4*u[3]^3*   u[1]^2*(12+3*u[4]^4+3*u[1]^4+(-4)*u[2]^2+3*u[2]^4+2*u[1]^2*((-2)+   u[2]^2)+(-2)*u[4]^2*((-2)+u[1]^2+u[2]^2))+8*u[4]*u[1]^3*u[2]*(u[4]^4+(   u[1]^2+(-1)*u[2]^2)^2+2*u[4]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+u[1]^2+   u[2]^2))+2*u[3]^5*(12+3*u[4]^4+9*u[1]^4+(-4)*u[2]^2+3*u[2]^4+8*   u[1]^2*((-2)+u[2]^2)+(-2)*u[4]^2*((-2)+4*u[1]^2+u[2]^2))+(-8)*   u[3]^2*u[4]*u[1]*u[2]*(u[4]^4+(-5)*u[1]^4+u[2]^4+(-8)*u[1]^2*(2+u[2]^2)   +(-4)*(3+u[2]^2)+2*u[4]^2*(2+4*u[1]^2+u[2]^2))+u[3]*((-1)*u[4]^8+   3*u[1]^8+8*u[1]^2*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*(2+u[2]^2)^4+(   -4)*u[4]^6*((-2)+2*u[1]^2+u[2]^2)+(-2)*u[1]^4*(20+36*u[2]^2+5*   u[2]^4)+(-2)*u[4]^4*(12+5*u[1]^4+(-4)*u[2]^2+3*u[2]^4+4*u[1]^2*((   -2)+u[2]^2))+4*u[4]^2*((-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+u[1]^4*(   18+23*u[2]^2)+2*u[1]^2*(4+20*u[2]^2+u[2]^4)))));
#     du[2]=-(4*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*  2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+  u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2])   ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+   (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)   *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*((-3)*u[4]^9+(-24)*   u[3]*u[4]^6*u[1]*u[2]+8*u[3]*u[4]^4*u[1]*u[2]*((-4)+2*u[3]^2+(-2)*   u[1]^2+(-7)*u[2]^2)+4*u[4]^7*(4+2*u[3]^2+(-2)*u[1]^2+(-3)*u[2]^2)+  u[4]*((u[3]^4+2*u[3]^2*((-2)+u[1]^2)+(2+u[1]^2)^2)^2+8*(u[3]^6+  u[3]^4*((-2)+u[1]^2)+(-1)*((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[3]^2*  (4+20*u[1]^2+u[1]^4))*u[2]^2+2*(20+5*u[3]^4+36*u[1]^2+5*u[1]^4+(-2)*u[3]^2*(18+23*u[1]^2))*u[2]^4+(-3)*u[2]^8)+(-4)*u[4]^3*  u[2]^2*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+2*((-2)+u[1]^2)*u[2]^2+   3*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+u[2]^2))+(-8)*u[3]*u[1]*u[2]^3*(   u[3]^4+(u[1]^2+(-1)*u[2]^2)^2+2*u[3]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+   u[1]^2+u[2]^2))+(-2)*u[4]^5*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+8*((  -2)+u[1]^2)*u[2]^2+9*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+4*u[2]^2))+8*  u[3]*u[4]^2*u[1]*u[2]*((-12)+u[3]^4+(-4)*u[1]^2+u[1]^4+(-8)*(2+u[1]^2)   *u[2]^2+(-5)*u[2]^4+2*u[3]^2*(2+u[1]^2+4*u[2]^2))));
#     du[3]=-(4*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*( 2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+( -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)*((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[1]^9+24*u[3]*u[4]*u[1]^6*u[2]+4*u[1]^7*(4+3*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+8*u[3]*u[4]* u[1]^4*u[2]*((-4)+7*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+4*u[3]^2*u[1]^3*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[2]^2))+2*u[1]^5*(12+9*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+8*u[3]^2*(2+u[4]^2+(-1)* u[2]^2))+8*u[3]*u[4]*u[1]^2*u[2]*(5*u[3]^4+8*u[3]^2*((-2)+u[4]^2+(-1) *u[2]^2)+4*(3+(-1)*u[4]^2+u[2]^2)+(-1)*(u[4]^2+u[2]^2)^2)+8*u[3]^3*u[4]*u[2]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[2]^2)+2*u[3]^2*(2+( -1)*u[4]^2+u[2]^2)+(u[4]^2+u[2]^2)^2)+u[1]*(3*u[3]^8+(-1)*(u[4]^4+2* u[4]^2*((-2)+u[2]^2)+(2+u[2]^2)^2)^2+(-2)*u[3]^4*(20+5*u[4]^4+ 36*u[2]^2+5*u[2]^4+(-2)*u[4]^2*(18+23*u[2]^2))+8*u[3]^2*(u[4]^6+ u[4]^4*((-2)+u[2]^2)+(-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*u[4]^2*(4+20*u[2]^2+u[2]^4)))));
#     du[4]=-((-4)*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*((-8)*u[3]*u[4]^3*u[1]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[1]^2)+2*u[3]^2*(2+(-1)*u[4]^2+u[1]^2)+(u[4]^2+u[1]^2)^2)+(u[3]^8+(-3)*u[4]^8+8*u[4]^2*((-2)+u[1]^2)*(2+u[1]^2)^2+(2+u[1]^2)^4+4*u[3]^6*((-2)+(-2)*u[4]^2+u[1]^2)+2*u[3]^4*(12+8*u[4]^2+5*u[4]^4+(-4)*(1+u[4]^2)*u[1]^2+3*u[1]^4)+2*u[4]^4*(20+36*u[1]^2+5*u[1]^4)+4*u[3]^2*(((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[4]^4*(18+23*u[1]^2)+2*u[4]^2*(4+20*u[1]^2+u[1]^4)))*u[2]+8*u[3]*u[4]*u[1]*(u[3]^4+(-5)*u[4]^4+u[1]^4+8*u[4]^2*(2+u[1]^2)+(-4)*(3+u[1]^2)+2*u[3]^2*(2+(-4)*u[4]^2+u[1]^2))*u[2]^2+( -4)*u[4]^2*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[1]^2+3*u[1]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[1]^2))*u[2]^3+(-8)*u[3]*u[4]* u[1]*(2*u[3]^2+7*u[4]^2+(-2)*(2+u[1]^2))*u[2]^4+(-2)*(12+3*u[3]^4+ 9*u[4]^4+(-4)*u[1]^2+3*u[1]^4+u[3]^2*(4+8*u[4]^2+(-2)*u[1]^2)+(-8) *u[4]^2*((-2)+u[1]^2))*u[2]^5+(-24)*u[3]*u[4]*u[1]*u[2]^6+(-4)*(4+ 2*u[3]^2+3*u[4]^2+(-2)*u[1]^2)*u[2]^7+(-3)*u[2]^9)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+ (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+( u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2) *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2));
#     du[5]=sqrt(abs(du[1])^2+ abs(du[2])^2+ abs(du[3])^2+ abs(du[4])^2)
#     du[6]=0;
# end

function Dimer_Eq_of_M_Backwards(du, u,p, t)
    p=2
    du[1]=((-4)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4]) +(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(   2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+  (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)  *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[3]^9+24*u[3]^6*   u[4]*u[1]*u[2]+4*u[3]^7*((-4)+(-2)*u[4]^2+3*u[1]^2+2*u[2]^2)+8*   u[3]^4*u[4]*u[1]*u[2]*(4+(-2)*u[4]^2+7*u[1]^2+2*u[2]^2)+4*u[3]^3*   u[1]^2*(12+3*u[4]^4+3*u[1]^4+(-4)*u[2]^2+3*u[2]^4+2*u[1]^2*((-2)+   u[2]^2)+(-2)*u[4]^2*((-2)+u[1]^2+u[2]^2))+8*u[4]*u[1]^3*u[2]*(u[4]^4+(   u[1]^2+(-1)*u[2]^2)^2+2*u[4]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+u[1]^2+   u[2]^2))+2*u[3]^5*(12+3*u[4]^4+9*u[1]^4+(-4)*u[2]^2+3*u[2]^4+8*   u[1]^2*((-2)+u[2]^2)+(-2)*u[4]^2*((-2)+4*u[1]^2+u[2]^2))+(-8)*   u[3]^2*u[4]*u[1]*u[2]*(u[4]^4+(-5)*u[1]^4+u[2]^4+(-8)*u[1]^2*(2+u[2]^2)   +(-4)*(3+u[2]^2)+2*u[4]^2*(2+4*u[1]^2+u[2]^2))+u[3]*((-1)*u[4]^8+   3*u[1]^8+8*u[1]^2*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*(2+u[2]^2)^4+(   -4)*u[4]^6*((-2)+2*u[1]^2+u[2]^2)+(-2)*u[1]^4*(20+36*u[2]^2+5*   u[2]^4)+(-2)*u[4]^4*(12+5*u[1]^4+(-4)*u[2]^2+3*u[2]^4+4*u[1]^2*((   -2)+u[2]^2))+4*u[4]^2*((-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+u[1]^4*(   18+23*u[2]^2)+2*u[1]^2*(4+20*u[2]^2+u[2]^4)))));
    du[2]=(4*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*  2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+  u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2])   ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+   (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)   *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*((-3)*u[4]^9+(-24)*   u[3]*u[4]^6*u[1]*u[2]+8*u[3]*u[4]^4*u[1]*u[2]*((-4)+2*u[3]^2+(-2)*   u[1]^2+(-7)*u[2]^2)+4*u[4]^7*(4+2*u[3]^2+(-2)*u[1]^2+(-3)*u[2]^2)+  u[4]*((u[3]^4+2*u[3]^2*((-2)+u[1]^2)+(2+u[1]^2)^2)^2+8*(u[3]^6+  u[3]^4*((-2)+u[1]^2)+(-1)*((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[3]^2*  (4+20*u[1]^2+u[1]^4))*u[2]^2+2*(20+5*u[3]^4+36*u[1]^2+5*u[1]^4+(-2)*u[3]^2*(18+23*u[1]^2))*u[2]^4+(-3)*u[2]^8)+(-4)*u[4]^3*  u[2]^2*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+2*((-2)+u[1]^2)*u[2]^2+   3*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+u[2]^2))+(-8)*u[3]*u[1]*u[2]^3*(   u[3]^4+(u[1]^2+(-1)*u[2]^2)^2+2*u[3]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+   u[1]^2+u[2]^2))+(-2)*u[4]^5*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+8*((  -2)+u[1]^2)*u[2]^2+9*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+4*u[2]^2))+8*  u[3]*u[4]^2*u[1]*u[2]*((-12)+u[3]^4+(-4)*u[1]^2+u[1]^4+(-8)*(2+u[1]^2)   *u[2]^2+(-5)*u[2]^4+2*u[3]^2*(2+u[1]^2+4*u[2]^2))));
    du[3]=(4*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*( 2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+( -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)*((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[1]^9+24*u[3]*u[4]*u[1]^6*u[2]+4*u[1]^7*(4+3*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+8*u[3]*u[4]* u[1]^4*u[2]*((-4)+7*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+4*u[3]^2*u[1]^3*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[2]^2))+2*u[1]^5*(12+9*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+8*u[3]^2*(2+u[4]^2+(-1)* u[2]^2))+8*u[3]*u[4]*u[1]^2*u[2]*(5*u[3]^4+8*u[3]^2*((-2)+u[4]^2+(-1) *u[2]^2)+4*(3+(-1)*u[4]^2+u[2]^2)+(-1)*(u[4]^2+u[2]^2)^2)+8*u[3]^3*u[4]*u[2]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[2]^2)+2*u[3]^2*(2+( -1)*u[4]^2+u[2]^2)+(u[4]^2+u[2]^2)^2)+u[1]*(3*u[3]^8+(-1)*(u[4]^4+2* u[4]^2*((-2)+u[2]^2)+(2+u[2]^2)^2)^2+(-2)*u[3]^4*(20+5*u[4]^4+ 36*u[2]^2+5*u[2]^4+(-2)*u[4]^2*(18+23*u[2]^2))+8*u[3]^2*(u[4]^6+ u[4]^4*((-2)+u[2]^2)+(-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*u[4]^2*(4+20*u[2]^2+u[2]^4)))));
    du[4]=((-4)*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*((-8)*u[3]*u[4]^3*u[1]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[1]^2)+2*u[3]^2*(2+(-1)*u[4]^2+u[1]^2)+(u[4]^2+u[1]^2)^2)+(u[3]^8+(-3)*u[4]^8+8*u[4]^2*((-2)+u[1]^2)*(2+u[1]^2)^2+(2+u[1]^2)^4+4*u[3]^6*((-2)+(-2)*u[4]^2+u[1]^2)+2*u[3]^4*(12+8*u[4]^2+5*u[4]^4+(-4)*(1+u[4]^2)*u[1]^2+3*u[1]^4)+2*u[4]^4*(20+36*u[1]^2+5*u[1]^4)+4*u[3]^2*(((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[4]^4*(18+23*u[1]^2)+2*u[4]^2*(4+20*u[1]^2+u[1]^4)))*u[2]+8*u[3]*u[4]*u[1]*(u[3]^4+(-5)*u[4]^4+u[1]^4+8*u[4]^2*(2+u[1]^2)+(-4)*(3+u[1]^2)+2*u[3]^2*(2+(-4)*u[4]^2+u[1]^2))*u[2]^2+( -4)*u[4]^2*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[1]^2+3*u[1]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[1]^2))*u[2]^3+(-8)*u[3]*u[4]* u[1]*(2*u[3]^2+7*u[4]^2+(-2)*(2+u[1]^2))*u[2]^4+(-2)*(12+3*u[3]^4+ 9*u[4]^4+(-4)*u[1]^2+3*u[1]^4+u[3]^2*(4+8*u[4]^2+(-2)*u[1]^2)+(-8) *u[4]^2*((-2)+u[1]^2))*u[2]^5+(-24)*u[3]*u[4]*u[1]*u[2]^6+(-4)*(4+ 2*u[3]^2+3*u[4]^2+(-2)*u[1]^2)*u[2]^7+(-3)*u[2]^9)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+ (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+( u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2) *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2));
    du[5]=(abs(du[1])^(1/p))+ (abs(du[2])^(1/p))+ (abs(du[3])^(1/p))+ (abs(du[4])^(1/p)) ;
    du[6]=0;
end

# function Dimer_Eq_of_M_Backwards(du, u,p, t)
#     du[1]=((-4)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4]) +(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(   2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+  (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)  *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[3]^9+24*u[3]^6*   u[4]*u[1]*u[2]+4*u[3]^7*((-4)+(-2)*u[4]^2+3*u[1]^2+2*u[2]^2)+8*   u[3]^4*u[4]*u[1]*u[2]*(4+(-2)*u[4]^2+7*u[1]^2+2*u[2]^2)+4*u[3]^3*   u[1]^2*(12+3*u[4]^4+3*u[1]^4+(-4)*u[2]^2+3*u[2]^4+2*u[1]^2*((-2)+   u[2]^2)+(-2)*u[4]^2*((-2)+u[1]^2+u[2]^2))+8*u[4]*u[1]^3*u[2]*(u[4]^4+(   u[1]^2+(-1)*u[2]^2)^2+2*u[4]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+u[1]^2+   u[2]^2))+2*u[3]^5*(12+3*u[4]^4+9*u[1]^4+(-4)*u[2]^2+3*u[2]^4+8*   u[1]^2*((-2)+u[2]^2)+(-2)*u[4]^2*((-2)+4*u[1]^2+u[2]^2))+(-8)*   u[3]^2*u[4]*u[1]*u[2]*(u[4]^4+(-5)*u[1]^4+u[2]^4+(-8)*u[1]^2*(2+u[2]^2)   +(-4)*(3+u[2]^2)+2*u[4]^2*(2+4*u[1]^2+u[2]^2))+u[3]*((-1)*u[4]^8+   3*u[1]^8+8*u[1]^2*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*(2+u[2]^2)^4+(   -4)*u[4]^6*((-2)+2*u[1]^2+u[2]^2)+(-2)*u[1]^4*(20+36*u[2]^2+5*   u[2]^4)+(-2)*u[4]^4*(12+5*u[1]^4+(-4)*u[2]^2+3*u[2]^4+4*u[1]^2*((   -2)+u[2]^2))+4*u[4]^2*((-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+u[1]^4*(   18+23*u[2]^2)+2*u[1]^2*(4+20*u[2]^2+u[2]^4)))));
#     du[2]=(4*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*  2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+  u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2])   ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+   (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)   *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*((-3)*u[4]^9+(-24)*   u[3]*u[4]^6*u[1]*u[2]+8*u[3]*u[4]^4*u[1]*u[2]*((-4)+2*u[3]^2+(-2)*   u[1]^2+(-7)*u[2]^2)+4*u[4]^7*(4+2*u[3]^2+(-2)*u[1]^2+(-3)*u[2]^2)+  u[4]*((u[3]^4+2*u[3]^2*((-2)+u[1]^2)+(2+u[1]^2)^2)^2+8*(u[3]^6+  u[3]^4*((-2)+u[1]^2)+(-1)*((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[3]^2*  (4+20*u[1]^2+u[1]^4))*u[2]^2+2*(20+5*u[3]^4+36*u[1]^2+5*u[1]^4+(-2)*u[3]^2*(18+23*u[1]^2))*u[2]^4+(-3)*u[2]^8)+(-4)*u[4]^3*  u[2]^2*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+2*((-2)+u[1]^2)*u[2]^2+   3*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+u[2]^2))+(-8)*u[3]*u[1]*u[2]^3*(   u[3]^4+(u[1]^2+(-1)*u[2]^2)^2+2*u[3]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+   u[1]^2+u[2]^2))+(-2)*u[4]^5*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+8*((  -2)+u[1]^2)*u[2]^2+9*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+4*u[2]^2))+8*  u[3]*u[4]^2*u[1]*u[2]*((-12)+u[3]^4+(-4)*u[1]^2+u[1]^4+(-8)*(2+u[1]^2)   *u[2]^2+(-5)*u[2]^4+2*u[3]^2*(2+u[1]^2+4*u[2]^2))));
#     du[3]=(4*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*( 2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+( -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)*((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[1]^9+24*u[3]*u[4]*u[1]^6*u[2]+4*u[1]^7*(4+3*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+8*u[3]*u[4]* u[1]^4*u[2]*((-4)+7*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+4*u[3]^2*u[1]^3*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[2]^2))+2*u[1]^5*(12+9*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+8*u[3]^2*(2+u[4]^2+(-1)* u[2]^2))+8*u[3]*u[4]*u[1]^2*u[2]*(5*u[3]^4+8*u[3]^2*((-2)+u[4]^2+(-1) *u[2]^2)+4*(3+(-1)*u[4]^2+u[2]^2)+(-1)*(u[4]^2+u[2]^2)^2)+8*u[3]^3*u[4]*u[2]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[2]^2)+2*u[3]^2*(2+( -1)*u[4]^2+u[2]^2)+(u[4]^2+u[2]^2)^2)+u[1]*(3*u[3]^8+(-1)*(u[4]^4+2* u[4]^2*((-2)+u[2]^2)+(2+u[2]^2)^2)^2+(-2)*u[3]^4*(20+5*u[4]^4+ 36*u[2]^2+5*u[2]^4+(-2)*u[4]^2*(18+23*u[2]^2))+8*u[3]^2*(u[4]^6+ u[4]^4*((-2)+u[2]^2)+(-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*u[4]^2*(4+20*u[2]^2+u[2]^4)))));
#     du[4]=((-4)*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*((-8)*u[3]*u[4]^3*u[1]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[1]^2)+2*u[3]^2*(2+(-1)*u[4]^2+u[1]^2)+(u[4]^2+u[1]^2)^2)+(u[3]^8+(-3)*u[4]^8+8*u[4]^2*((-2)+u[1]^2)*(2+u[1]^2)^2+(2+u[1]^2)^4+4*u[3]^6*((-2)+(-2)*u[4]^2+u[1]^2)+2*u[3]^4*(12+8*u[4]^2+5*u[4]^4+(-4)*(1+u[4]^2)*u[1]^2+3*u[1]^4)+2*u[4]^4*(20+36*u[1]^2+5*u[1]^4)+4*u[3]^2*(((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[4]^4*(18+23*u[1]^2)+2*u[4]^2*(4+20*u[1]^2+u[1]^4)))*u[2]+8*u[3]*u[4]*u[1]*(u[3]^4+(-5)*u[4]^4+u[1]^4+8*u[4]^2*(2+u[1]^2)+(-4)*(3+u[1]^2)+2*u[3]^2*(2+(-4)*u[4]^2+u[1]^2))*u[2]^2+( -4)*u[4]^2*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[1]^2+3*u[1]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[1]^2))*u[2]^3+(-8)*u[3]*u[4]* u[1]*(2*u[3]^2+7*u[4]^2+(-2)*(2+u[1]^2))*u[2]^4+(-2)*(12+3*u[3]^4+ 9*u[4]^4+(-4)*u[1]^2+3*u[1]^4+u[3]^2*(4+8*u[4]^2+(-2)*u[1]^2)+(-8) *u[4]^2*((-2)+u[1]^2))*u[2]^5+(-24)*u[3]*u[4]*u[1]*u[2]^6+(-4)*(4+ 2*u[3]^2+3*u[4]^2+(-2)*u[1]^2)*u[2]^7+(-3)*u[2]^9)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+ (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+( u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2) *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2));
#     du[5]=sqrt(abs(du[1])^2+ abs(du[2])^2+ abs(du[3])^2+ abs(du[4])^2)
# end

function Dimer_Eq_of_M_hits(du, u,p, t)
  du[1]=-((-4)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4]) +(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(   2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+  (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)  *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[3]^9+24*u[3]^6*   u[4]*u[1]*u[2]+4*u[3]^7*((-4)+(-2)*u[4]^2+3*u[1]^2+2*u[2]^2)+8*   u[3]^4*u[4]*u[1]*u[2]*(4+(-2)*u[4]^2+7*u[1]^2+2*u[2]^2)+4*u[3]^3*   u[1]^2*(12+3*u[4]^4+3*u[1]^4+(-4)*u[2]^2+3*u[2]^4+2*u[1]^2*((-2)+   u[2]^2)+(-2)*u[4]^2*((-2)+u[1]^2+u[2]^2))+8*u[4]*u[1]^3*u[2]*(u[4]^4+(   u[1]^2+(-1)*u[2]^2)^2+2*u[4]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+u[1]^2+   u[2]^2))+2*u[3]^5*(12+3*u[4]^4+9*u[1]^4+(-4)*u[2]^2+3*u[2]^4+8*   u[1]^2*((-2)+u[2]^2)+(-2)*u[4]^2*((-2)+4*u[1]^2+u[2]^2))+(-8)*   u[3]^2*u[4]*u[1]*u[2]*(u[4]^4+(-5)*u[1]^4+u[2]^4+(-8)*u[1]^2*(2+u[2]^2)   +(-4)*(3+u[2]^2)+2*u[4]^2*(2+4*u[1]^2+u[2]^2))+u[3]*((-1)*u[4]^8+   3*u[1]^8+8*u[1]^2*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*(2+u[2]^2)^4+(   -4)*u[4]^6*((-2)+2*u[1]^2+u[2]^2)+(-2)*u[1]^4*(20+36*u[2]^2+5*   u[2]^4)+(-2)*u[4]^4*(12+5*u[1]^4+(-4)*u[2]^2+3*u[2]^4+4*u[1]^2*((   -2)+u[2]^2))+4*u[4]^2*((-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+u[1]^4*(   18+23*u[2]^2)+2*u[1]^2*(4+20*u[2]^2+u[2]^4)))));
  du[2]=-(4*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*  2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+  u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*   u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2)   *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2])   ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+   (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(   -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(   u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)   *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*((-3)*u[4]^9+(-24)*   u[3]*u[4]^6*u[1]*u[2]+8*u[3]*u[4]^4*u[1]*u[2]*((-4)+2*u[3]^2+(-2)*   u[1]^2+(-7)*u[2]^2)+4*u[4]^7*(4+2*u[3]^2+(-2)*u[1]^2+(-3)*u[2]^2)+  u[4]*((u[3]^4+2*u[3]^2*((-2)+u[1]^2)+(2+u[1]^2)^2)^2+8*(u[3]^6+  u[3]^4*((-2)+u[1]^2)+(-1)*((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[3]^2*  (4+20*u[1]^2+u[1]^4))*u[2]^2+2*(20+5*u[3]^4+36*u[1]^2+5*u[1]^4+(-2)*u[3]^2*(18+23*u[1]^2))*u[2]^4+(-3)*u[2]^8)+(-4)*u[4]^3*  u[2]^2*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+2*((-2)+u[1]^2)*u[2]^2+   3*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+u[2]^2))+(-8)*u[3]*u[1]*u[2]^3*(   u[3]^4+(u[1]^2+(-1)*u[2]^2)^2+2*u[3]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+   u[1]^2+u[2]^2))+(-2)*u[4]^5*(12+3*u[3]^4+(-4)*u[1]^2+3*u[1]^4+8*((  -2)+u[1]^2)*u[2]^2+9*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+4*u[2]^2))+8*  u[3]*u[4]^2*u[1]*u[2]*((-12)+u[3]^4+(-4)*u[1]^2+u[1]^4+(-8)*(2+u[1]^2)   *u[2]^2+(-5)*u[2]^4+2*u[3]^2*(2+u[1]^2+4*u[2]^2))));
  du[3]=-(4*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*( 2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(u[4]^2+u[2]^2)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+( -1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+(u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2)*((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2)*(3*u[1]^9+24*u[3]*u[4]*u[1]^6*u[2]+4*u[1]^7*(4+3*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+8*u[3]*u[4]* u[1]^4*u[2]*((-4)+7*u[3]^2+2*u[4]^2+(-2)*u[2]^2)+4*u[3]^2*u[1]^3*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[2]^2))+2*u[1]^5*(12+9*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[2]^2+3*u[2]^4+8*u[3]^2*(2+u[4]^2+(-1)* u[2]^2))+8*u[3]*u[4]*u[1]^2*u[2]*(5*u[3]^4+8*u[3]^2*((-2)+u[4]^2+(-1) *u[2]^2)+4*(3+(-1)*u[4]^2+u[2]^2)+(-1)*(u[4]^2+u[2]^2)^2)+8*u[3]^3*u[4]*u[2]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[2]^2)+2*u[3]^2*(2+( -1)*u[4]^2+u[2]^2)+(u[4]^2+u[2]^2)^2)+u[1]*(3*u[3]^8+(-1)*(u[4]^4+2* u[4]^2*((-2)+u[2]^2)+(2+u[2]^2)^2)^2+(-2)*u[3]^4*(20+5*u[4]^4+ 36*u[2]^2+5*u[2]^4+(-2)*u[4]^2*(18+23*u[2]^2))+8*u[3]^2*(u[4]^6+ u[4]^4*((-2)+u[2]^2)+(-1)*((-2)+u[2]^2)*(2+u[2]^2)^2+(-1)*u[4]^2*(4+20*u[2]^2+u[2]^4)))));
  du[4]=-((-4)*(u[3]^2+u[1]^2)*(2+u[3]^2+(-2)*2^(1/2)*u[4]+u[4]^2+2*u[3]*((-1)*2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+2*u[3]*(2^(1/2)+u[4])+(u[1]+(-1)*u[2])^2)^(-2)*((-8)*u[3]*u[4]^3*u[1]*(u[3]^4+4*((-3)+u[4]^2+(-1)*u[1]^2)+2*u[3]^2*(2+(-1)*u[4]^2+u[1]^2)+(u[4]^2+u[1]^2)^2)+(u[3]^8+(-3)*u[4]^8+8*u[4]^2*((-2)+u[1]^2)*(2+u[1]^2)^2+(2+u[1]^2)^4+4*u[3]^6*((-2)+(-2)*u[4]^2+u[1]^2)+2*u[3]^4*(12+8*u[4]^2+5*u[4]^4+(-4)*(1+u[4]^2)*u[1]^2+3*u[1]^4)+2*u[4]^4*(20+36*u[1]^2+5*u[1]^4)+4*u[3]^2*(((-2)+u[1]^2)*(2+u[1]^2)^2+(-1)*u[4]^4*(18+23*u[1]^2)+2*u[4]^2*(4+20*u[1]^2+u[1]^4)))*u[2]+8*u[3]*u[4]*u[1]*(u[3]^4+(-5)*u[4]^4+u[1]^4+8*u[4]^2*(2+u[1]^2)+(-4)*(3+u[1]^2)+2*u[3]^2*(2+(-4)*u[4]^2+u[1]^2))*u[2]^2+( -4)*u[4]^2*(12+3*u[3]^4+4*u[4]^2+3*u[4]^4+(-2)*(2+u[4]^2)*u[1]^2+3*u[1]^4+2*u[3]^2*(2+u[4]^2+(-1)*u[1]^2))*u[2]^3+(-8)*u[3]*u[4]* u[1]*(2*u[3]^2+7*u[4]^2+(-2)*(2+u[1]^2))*u[2]^4+(-2)*(12+3*u[3]^4+ 9*u[4]^4+(-4)*u[1]^2+3*u[1]^4+u[3]^2*(4+8*u[4]^2+(-2)*u[1]^2)+(-8) *u[4]^2*((-2)+u[1]^2))*u[2]^5+(-24)*u[3]*u[4]*u[1]*u[2]^6+(-4)*(4+ 2*u[3]^2+3*u[4]^2+(-2)*u[1]^2)*u[2]^7+(-3)*u[2]^9)*(2+u[3]^2+2* u[3]*(2^(1/2)+(-1)*u[4])+(-2)*2^(1/2)*u[4]+u[4]^2+(u[1]+u[2])^2)^(-2) *(2+u[3]^2+2*2^(1/2)*u[4]+u[4]^2+(-2)*u[3]*(2^(1/2)+u[4])+(u[1]+u[2]) ^2)^(-2)*((u[3]^2+u[1]^2)*(2+(-2)*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+ (u[1]+(-1)*u[2])^2)^(-1)*(2+2*2^(1/2)*(u[3]+u[4])+(u[3]+u[4])^2+(u[1]+(-1)*u[2])^2)^(-1)*(u[4]^2+u[2]^2)*(2+2*2^(1/2)*(u[3]+(-1)*u[4])+( u[3]+(-1)*u[4])^2+(u[1]+u[2])^2)^(-1)*(2+(u[3]+(-1)*u[4])^2+2*2^(1/2) *((-1)*u[3]+u[4])+(u[1]+u[2])^2)^(-1))^(-1/2));
  du[5]=0;
  du[6]=0;
end


text:
using JuMP
using Gurobi
using MathProgBase

# solver_Outputflag_check check if user didn't turn Original solver OutputFlag off
function solver_Outputflag_check(solver::ASCIIString)
  if contains(solver , "Any[<:OutputFlag,0>]") == false
    if contains(solver , "GurobiSolver")
      return "GurobiSolver"
    ###elseif other solver can be used in JuMP
    end
  end
end

# getFullsizeModelData: getting the data of full-size model which is entered by users
function getFullsizeModelData(m::JuMP.Model)
  # parameter
  # input:
  # m: the model entered by users using JuMP
  # output:
  # fullsizeObjSense: the objective sense of full-size model (Max or Min)
  # fullsizeObjCoeff: the coefficients of objective function of full-size model
  # fullsizeConstLB: the lower bound of each constraint of full-size model, which is a vector
  # fullsizeConstUB: the upper bound of each constraint of full-size model, which is a vector
  # fullsizeConstMatrix: the constraint matrix of full-size model
  # fullsizeVarLB: the lower bound of each varible of full-size model, which is a vector
  # fullsizeVarUB: the upper bound of each varible of full-size model, which is a vector
  # fullsizeNumVar: number of variables of full-size model
  # fullsizeNumConst: number of constraints of full-size model

  fullsizeObjSense = m.objSense
  fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB = JuMP.prepProblemBounds(m)
  fullsizeConstMatrix = JuMP.prepConstrMatrix(m)
  fullsizeVarLB = m.colLower
  fullsizeVarUB = m.colUpper

  fullsizeNumVar = m.numCols
  fullsizeNumConst = length(fullsizeConstLB)

  return fullsizeObjSense, fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeNumVar, fullsizeNumConst
end

# getBigger: getting the top n largest elements in vector v
function largestn(v::Vector{Float64}, n::Int64)
  # function getBigger(v::Vector{Float64}, n::Int)
  # parameters
  # input
  # v: the vector which is sorted to find top n largest elements
  # n: the number of elements which have to be finded from vector v
  # output
  # TopN_Idx: top n largest elements, which are the indices w.r.t the vector v
  # TopN_Value: values of top n largest elements

  # initialization
  TopN_Value = v[1:n]  # the top n largest value in vector v
  TopN_Idx = collect(1:n) # the top n largest elements in vector v

  minValue, minIdx = findmin(TopN_Value)

  for i = n+1:length(v)
    if v[i] > minValue
      TopN_Value[minIdx] = v[i]
      TopN_Idx[minIdx] = i

      minValue, minIdx = findmin(TopN_Value)
    end
  end

  return TopN_Value, TopN_Idx
end

# smallestn: getting the top n smallest elements in vector v
function smallestn(v::Vector{Float64}, n::Int64)
  # parameters
  # input
  # v: the vector which is sorted to find top n smallest elements
  # n: the number of elements which have to be finded from vector v
  # output
  # BottomN_Idx: top n smallest elements, which are the indices w.r.t the vector v
  # BottomN_Value: values of top n smallest elements

  # initialization
  BottomN_Value = v[1:n] # the top n smallest value in vector v
  BottomN_Idx = collect(1:n) # the top n smallest elements in vector v

  maxValue, maxIdx = findmax(BottomN_Value)

  for i = n+1:length(v)
    if v[i] < maxValue
      BottomN_Value[maxIdx] = v[i]
      BottomN_Idx[maxIdx] = i

      maxValue, maxIdx = findmax(BottomN_Value)
    end
  end

  return BottomN_Value, BottomN_Idx
end

ext: jl


text:
# We use the max norm as default distance. In the Distance.jl package
# the max norm is called the Chebyshev distance.
default_knn_constructor(data) = KDTree(data, Chebyshev())

# Classic Kozachenko-Leonenko k-nearest neighbor differential entropy estimator from
# samples from a continuous random variable. X should be a d*N dimensional
# matrix where each column is a d-dimensional sample from the distribution.
#
# Example: X = [1.3 3.7 5.1 2.4] if the distribution is 1-dimensional and
# we have 4 samples.
#
# Differential entropy (also referred to as continuous entropy) is a concept in
# information theory that extends the idea of (Shannon) entropy, a measure of
# average surprisal of a random variable, to continuous probability distributions.
# For more information see: http://en.wikipedia.org/wiki/Differential_entropy
#
# Note that differential entropy can be negative, unlike discrete entropy.
#
# This is based on the entropy() function in NPEET version 1.1 which, in turn,
# states that it is from Kozachenko and Leonenko as described in the paper:
# A. Kraskov, H. Stogbauer, and P. Grassberger. Estimating mutual information.
# Phys. Rev. E, 69:066138, Jun 2004.
function differential_entropy{T <: Number}(X::Array{T, 2}, k=3, base=2, knnConstructor = default_knn_constructor)

  # Get dimensions
  d, N = size(X)

  # Check that arguments are valid
  if k >= N
    error("k ($k) must be smaller than the number of samples ($N)")
  end

  # Add small noise to break degeneracy. This can happen if samples have limited
  # resolution and just happen to be exactly the same even if the thing measured
  # would not be the same. By adding a very small noise term we avoid any numerical
  # problems because of this.
  xprim = add_noise(X, d, N, 1e-10)

  # Build kNN tree/graph/solution.
  knn_tree = knnConstructor(xprim)

  # Get the k-nearest neighbor distance per point
  knn_dists = zeros(N)
  for i in 1:N
    knn_dists[i] = knn_nearest_distance(xprim, knn_tree, xprim[:,i], k)
  end

  # and estimate the entropy:
  constant = digamma(N) - digamma(k) + d*log(2)
  return( (constant + d * mean(log(knn_dists))) / log(base) )

end


text:
module TestArrayOfComponents

println("\nTestArrayOfComponents: Demonstrating the handling of arrays of components")

using Modia
using Modia.Electric

# Desired:
#   using ModiaMath: plot
#   using Test
#
# In order that these packages need not to be defined in the user environment, they are included via Modia:
using Modia.ModiaMath: plot

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Modia.Test
end

@model LPfilter begin
    R = Resistor(R=100.0)
    C = Capacitor(C=2.5E-3, v=Float(start=0.0))
    V = ConstantVoltage(V=10.0)
    ground = Ground()
    @equations begin
        connect(V.n, ground.p)
        connect(V.p, R.p)
        connect(R.n, C.p)
        connect(C.n, V.n)
    end
end

# Array of components
@model TwoFilters begin
    F = [LPfilter(); LPfilter(C=Capacitor(C=10 * 2.5E-3))]
end
checkSimulation(TwoFilters, 1, "F[2].C.v", 3.2967995487381305)

nFilters = 10
# Array comprehensions of components
@model ManyFilters begin
    F = [LPfilter() for i in 1:nFilters]
end
checkSimulation(ManyFilters, 1, "F[1].C.v", 9.816405569531037)

nFilters = 10
# Array comprehensions of components
@model ManyDifferentFilters begin
    F = [LPfilter(C=Capacitor(C=i * 2.5E-3)) for i in 1:nFilters]
end
result = simulate(ManyDifferentFilters, 1)
plot(result, Tuple(["F[$i].C.v" for i in 1:nFilters]), heading="ManyDifferentFilters", figure=5)
@test result["F[1].C.v"][end] == 9.816405569531037
@test result["F[$nFilters].C.v"][end] == 3.2967995487381305

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

# Connecting arrays of components
@model AdvancedLPfilter begin
    M = [Resistor(R=1); Capacitor(C=1E-3)]
    V = ConstantVoltage(V=1)
    ground = Ground()
    @equations begin
        connect(V.n, ground.p)
        # connect(V.p, M[1].p)
        # connect(M[1].n, M[2].p)
        # connect(M[2].n, V.n)
    end
end

simulate(AdvancedLPfilter, 1)

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

# Experimental:
using Modia:@equation
#addEquation!(M, e) = begin @show e; push!(M.initializers, Modia.Instantiation.Equations([e])) end

@model LPfilterComponent begin
    R = Resistor(R=100.0)
    C = Capacitor(C=2.5E-3, v=Float(start=0.0))
    ground = Ground()
    in = Pin()
    out = Pin()
    @equations begin
        connect(in, R.p)
        connect(R.n, C.p)
        connect(C.n, ground.p)
        connect(C.p, out)
    end
end

@model ManyConnectedDifferentFilters begin
    V = ConstantVoltage(V=10.0)
    ground = Ground()
    F = [LPfilterComponent(C=Capacitor(C=i * 2.5E-3)) for i in 1:nFilters]
    @equations begin
        connect(V.n, ground.p)
#  connect(V.p, F[1].in)
    end
end

#=
# Not working yet
for i in 1:nFilters-1
  addEquation!(ManyConnectedDifferentFilters, @equation(connect(this.F[$i].out, this.F[$(i+1)].in)))
  # ERROR: LoadError: type Expr has no field out
end
simulate(ManyConnectedDifferentFilters, 1)
=#

#=
@model ManyConnectedDifferentFilters2 begin
    V = ConstantVoltage(V=10.0)
    ground = Ground()
    F = [LPfilterComponent(C=Capacitor(C=i * 2.5E-3)) for i in 1:nFilters]
    i = Integ()
    @equations begin
        connect(V.n, ground.p)
        #  connect(V.p, F[1].in)
        # Alternative connect constructs:

        for i in 1:length(F)                   # ERROR: LoadError: UndefVarError: i not defined
          connect(F[i].in, F[i+1].out)
        end

        #connect(F[1:end-1].in, F[2:end].out)     # ERROR: LoadError: UndefVarError: end not defined
        #connect(F[1:length(F)-1].in, F[2:length(F)].out)    # ERROR: LoadError: type Expr has no field in
        #[connect(F[i].in, F[i+1].out) for i in 1:length(F)] # ERROR: LoadError: UndefVarError: i not defined
        #connect(F[i].in, F[i+1].out) for i in length(F)     # ERROR: LoadError: syntax: extra token "for" after end of expression

    end
end

println("ManyConnectedDifferentFilters2")
simulate(ManyConnectedDifferentFilters2, 1)

=#
end


text:
using Documenter
include("../src/LightGraphs.jl")
using LightGraphs

# index is equal to the README for the time being
cp(normpath(@__FILE__, "../../README.md"), normpath(@__FILE__, "../src/index.md"); remove_destination=true)

# same for contributing and license
cp(normpath(@__FILE__, "../../CONTRIBUTING.md"), normpath(@__FILE__, "../src/contributing.md"); remove_destination=true)
cp(normpath(@__FILE__, "../../LICENSE.md"), normpath(@__FILE__, "../src/license.md"); remove_destination=true)

makedocs(modules=[LightGraphs], doctest = false)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo   = "github.com/JuliaGraphs/LightGraphs.jl.git",
    julia  = "0.6"
)

rm(normpath(@__FILE__, "../src/index.md"))
rm(normpath(@__FILE__, "../src/contributing.md"))
rm(normpath(@__FILE__, "../src/license.md"))


text:
function driver_version()
    ver = Ref{Cint}()
    rt.cudaDriverGetVersion(ver)
    return Int(ver[])
end

function runtime_version()
    ver = Ref{Cint}()
    rt.cudaRuntimeGetVersion(ver)
    return Int(ver[])
end


text:
# ---
# title: 16. 3Sum Closest
# id: problem16
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/3sum-closest/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers and an integer `target`, find three
# integers in `nums` such that the sum is closest to `target`. Return the sum of
# the three integers. You may assume that each input would have exactly one
# solution.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [-1,2,1,-4], target = 1
#     Output: 2
#     Explanation: The sum that is closest to the target is 2. (-1 + 2 + 1 = 2).
#
#
#
#
# **Constraints:**
#
#   * `3 <= nums.length <= 10^3`
#   * `-10^3 <= nums[i] <= 10^3`
#   * `-10^4 <= target <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Shapes
    abstract type AbstractShape end

    struct Separated <: AbstractShape
        Separated(sz...) = new()
    end
    struct Vectorized <: AbstractShape
        insize::Tuple
        Vectorized(sz::Integer...) = new(sz)
        Vectorized(sz::AbstractVector) = Vectorized(sz...)
        Vectorized(sz::Tuple) = Vectorized(sz...)
    end
    struct Combined <: AbstractShape
        Combined(sz...) = new()
    end

    Vec = Vectorized
end

isfixedsize(::A) where {A<:Shapes.AbstractShape} = isfixedsize(A)
isfixedsize(::Type{S}) where {S<:Shapes.AbstractShape}= false
isfixedsize(::Type{Shapes.Vec}) = true

function reshape_polyvec(::Shapes.Separated, ::AbstractNsolt, pvy::PolyphaseVector)
    [ reshape(@view(pvy.data[p,:]), pvy.nBlocks) for p in 1:size(pvy.data,1) ]
end

function reshape_polyvec(::Shapes.Combined, ::AbstractNsolt, pvy::PolyphaseVector)
    reshape(transpose(pvy.data), pvy.nBlocks..., size(pvy.data, 1))
end

function reshape_polyvec(::Shapes.Vec, ::AbstractNsolt, pvy::PolyphaseVector)
    vec(transpose(pvy.data))
end

function reshape_polyvec(::Shapes.Separated, ::ParallelFilters, y::AbstractArray)
    y
end

function reshape_polyvec(::Shapes.Combined, ::ParallelFilters{TF,D}, y::AbstractArray{TY,D}) where {TF,TY,D}
    cat(D+1, y...)
end

function reshape_polyvec(::Shapes.Vec, ::ParallelFilters, y::AbstractArray)
    vcat(vec.(y)...)
end

function reshape_coefs(::Shapes.Separated, ::AbstractNsolt, y::AbstractArray)
    PolyphaseVector(hcat(vec.(y)...) |> transpose |> Matrix, size(y[1]))
end

function reshape_coefs(::Shapes.Combined, ::AbstractNsolt, y::AbstractArray{T,D}) where {T,D}
    nBlocks = size(y)[1:D-1]
    outdata = transpose(reshape(y, prod(nBlocks), size(y,D)))
    PolyphaseVector(outdata, nBlocks)
end

function reshape_coefs(sv::Shapes.Vec, nsop::AbstractNsolt, y::AbstractArray)
    szout = fld.(sv.insize, decimations(nsop))
    ty = reshape(y, prod(szout), nchannels(nsop)) |> transpose
    PolyphaseVector(ty, szout)
end

function reshape_coefs(::Shapes.Separated, ::ParallelFilters, y::AbstractArray)
    y
end

function reshape_coefs(::Shapes.Combined, co::ParallelFilters{T,D}, y::AbstractArray) where {T,D}
    [ y[fill(:,D)..., p] for p in 1:nchannels(co)]
end

function reshape_coefs(sv::Shapes.Vec, co::ParallelFilters{T,D}, y::AbstractArray) where {T,D}
    ry = reshape(y, fld.(sv.insize, decimations(co))..., nchannels(co))
    [ ry[fill(:,D)..., p] for p in 1:nchannels(co) ]
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
export get_element,
       number_of_patterns,
       number_of_pixels,
       get_pixel,
       get_element!,
       spatial_position,
      pattern_dimension,
    get_element,
    assign!

using ColorTypes

function pattern_dimension(X)
  return length(get_element(X,1))
end

"""
```julia
function get_element!(o::Matrix,  img::Matrix{C}, i::Vector{Integer}) where C<:Colorant

```
"""
function get_element!(o::D, img::Matrix{C}, i::Vector{<:Integer}) where D<:AbstractArray  where C<:Colorant
    car_inds = CartesianIndices(img)[i]
    values = img[car_inds]
    for (j, pos) in enumerate(car_inds)
        @inbounds o[1, j] = pos[2]
        @inbounds o[2, j] = pos[1]
    end
    N = length(C)
    component = N >= 3 ? (comp1, comp2, comp3, alpha) : (comp1, alpha)
    for j=1:length(C)
        @inbounds o[2+j,:] = component[j].(values)
    end
end

function get_element!(o::D, data::Matrix{C}, i::Integer) where D<:AbstractArray  where C<:Number
    o[:] = data[:, i]
end

"""
```julia
function assign!(vec::T, val::C) where T<:AbstractArray where C<:Colorant
```
This function assigns the components of the color component val to a vector v
"""
function assign!(vec::T, val::C) where T<:AbstractArray where C<:Colorant
   N = length(C)
   component = N >= 3 ? (comp1, comp2, comp3, alpha) : (comp1, alpha)
    for j=1:length(C)
        @inbounds vec[j] = component[j](val)
    end
end

"""
```@julia
get_element!{T<:AbstractArray}(vec::T,  img::Matrix{Gray}, i::Integer)
```

Return through```vec``` the intensity image element  [x,y, i], where \$x,y\$ are the spatial
position of the pixel and the value i of the pixel \$(x,y)\$.
"""
function get_element!(vec::T,  img::Matrix{C}, i::Integer) where T<:AbstractArray where C<:Colorant
    ind    = spatial_position(img,i)
    @inbounds vec[1] = ind[2]
    @inbounds vec[2] = ind[1]
    assign!(view(vec,3:length(vec)), img[i])
end

"""
```
function get_element( img::Matrix{RGB}, i::Vector)
```
"""
function get_element(img::Matrix{T}, i::Vector)  where T<:Colorant
    m = zeros(length(T)+2,length(i))
    get_element!(m,img,i)
    return m
end

function get_element(img::Matrix{T}, i::Integer)  where T<:Colorant
   m = zeros(length(T)+2)
    get_element!(m,img,i)
    return m
end

"""
```@julia
number_of_patterns{T<:Any}(X::Matrix{T,3})
```

Return the number of pixels in the image
"""
number_of_patterns(X::Matrix{T}) where T<:Colorant = size(X,1)*size(X,2)

"""
```@julia
 spatial_position(X::Matrix, i::Int)
```
Returns the sub indexes from the linear index ```i```
"""
spatial_position(img::Matrix, i) = CartesianIndices(img)[i]

get_element(X::T,i) where T<:AbstractArray= view(X,:,i)
number_of_patterns(X::T) where T<:AbstractArray    = size(X,2)

get_element(X::Vector,i) = X[i]
number_of_patterns(X::Vector)     = length(X)

number_of_patterns(X::Array{T, 3}) where T = size(X, 2) * size(X, 3)
"""
```
function get_element(data::Array{T, 3}, i::Vector)  where T<:Number
```
"""
function get_element(data::Array{T, 3}, i::Vector)  where T<:Number
    m = zeros(size(data,1)+2,length(i))
    get_element!(m, data, i)
    return m
end

get_element(data::Array{T, 3}, i::Integer)  where T<:Number = get_element(data, [i])
get_element!(o::AbstractArray, data::Array{T,3}, i::Integer) where T<:Number = get_element!(o, data, [i])

function get_element!(o::AbstractArray, data::Array{T,3}, i::Vector) where T<:Number
    (_, nr, nc) = size(data)
    cart_indices = CartesianIndices((nr,nc))[i]
    @inbounds for (i, pos) in enumerate(cart_indices)
        o[3:end, i] = data[:, pos[1], pos[2]]
        o[1,i] = pos[2]
        o[2,i] = pos[1]
    end
end


text:
#######################
### Crop Single Way ###
#######################

function crop!(nodes::Dict, bounds::OpenStreetMapX.Bounds, way::OpenStreetMapX.Way)
    valid = falses(length(way.nodes)+2)
    n = 1
    while n <= length(way.nodes)
        if !haskey(nodes, way.nodes[n])
            splice!(way.nodes, n)
            splice!(valid, n+1)
        else
            valid[n+1] = OpenStreetMapX.inbounds(nodes[way.nodes[n]], bounds)
            n += 1
        end
    end
    if sum(valid) == 0
        return true
    elseif sum(valid) < (length(valid)-2)
		leave = trues(length(way.nodes))
        for i in 2:(length(valid)-1)
            if !valid[i]
                if valid[i-1] != valid[i]
                    if !OpenStreetMapX.onbounds(nodes[way.nodes[i-2]], bounds)
                        new_node = OpenStreetMapX.boundary_point(nodes[way.nodes[i-2]], nodes[way.nodes[i-1]],bounds)
                        new_id = OpenStreetMapX.add_new_node!(nodes, new_node)
                        way.nodes[i-1] = new_id
                    else
						leave[i-1] = false
                    end
                elseif valid[i] != valid[i+1]
                    if !OpenStreetMapX.onbounds(nodes[way.nodes[i]], bounds)
                        new_node = OpenStreetMapX.boundary_point(nodes[way.nodes[i-1]], nodes[way.nodes[i]],bounds)
                        new_id = OpenStreetMapX.add_new_node!(nodes, new_node)
                        way.nodes[i-1] = new_id
                    else
                        leave[i-1] = false
                    end
                else
                    leave[i-1] = false
                end
            end
        end
		way.nodes = way.nodes[leave]
        return false
    else
        return false
    end
end

#################
### Crop Ways ###
#################

function crop!(nodes::Dict, bounds::OpenStreetMapX.Bounds, ways::Vector{OpenStreetMapX.Way})
    leave = ways[[!OpenStreetMapX.crop!(nodes,bounds,way) for way in ways]]
	append!(empty!(ways),leave)
	return nothing
end

############################
### Crop Single Relation ###
############################

function crop!(nodes::Dict, bounds::OpenStreetMapX.Bounds, ways::Vector{OpenStreetMapX.Way},relations::Vector{OpenStreetMapX.Relation}, relation::OpenStreetMapX.Relation)
	valid = falses(length(relation.members))
	for i = 1:length(relation.members)
		ref = parse(Int,relation.members[i]["ref"])
		if relation.members[i]["type"] == "node" && haskey(nodes,ref)
			valid[i] = OpenStreetMapX.inbounds(nodes[ref],bounds)
		elseif relation.members[i]["type"] == "way"
			way_index = findfirst(way -> way.id == ref, ways)
			!isa(way_index,Nothing) && (valid[i] = !OpenStreetMapX.crop!(nodes, bounds, ways[way_index]))
		else
			relation_index = findfirst(relation -> relation.id == ref, relations)
			!isa(relation_index,Nothing) && (valid[i] = !OpenStreetMapX.crop!(nodes,bounds, ways, relations, relations[relation_index]))
		end
	end
	relation.members = relation.members[valid]
	if sum(valid) == 0
		return true
	else
		return false
	end
end

######################
### Crop Relations ###
######################

function crop!(nodes::Dict, bounds::OpenStreetMapX.Bounds, ways::Vector{OpenStreetMapX.Way}, relations::Vector{OpenStreetMapX.Relation})
    leave = relations[[!OpenStreetMapX.crop!(nodes,bounds,ways, relations, relation) for relation in relations]]
	append!(empty!(relations),leave)
	return nothing
end

####################################
### Crop Single Node and Feature ###
####################################

function crop!(nodes::Dict, bounds::OpenStreetMapX.Bounds, features::Dict, id::Int)
	if !OpenStreetMapX.inbounds(nodes[id], bounds)
		id in keys(features) && delete!(features, id)
		delete!(nodes,id)
	end
end

###############################
### Crop Nodes and Features ###
###############################

function crop!(nodes::Dict, bounds::OpenStreetMapX.Bounds, features::Dict)
    for (key, node) in nodes
        OpenStreetMapX.crop!(nodes,bounds,features,key)
    end
end

################
### Crop Map ###
################

function crop!(map::OpenStreetMapX.OSMData; crop_relations = true, crop_ways = true, crop_nodes = true)
	crop_relations && OpenStreetMapX.crop!(map.nodes, map.bounds, map.ways, map.relations)
	crop_ways && OpenStreetMapX.crop!(map.nodes, map.bounds, map.ways)
	crop_nodes && OpenStreetMapX.crop!(map.nodes, map.bounds, map.features)
end


text:
using Glob;
using DelimitedFiles;
using Quadmath;
using DecFP;

if length(ARGS) != 4
  println("usage: julia aggregate.jl <outfile> <indir> <pattern> <dielectric|polar>");
  exit(1);
end

outfile = open(ARGS[1], "w");
indir = ARGS[2];
pattern = Glob.GlobMatch(ARGS[3]);

if ARGS[4] == "dielectric"
  writedlm(outfile, ["E0" "K1" "K2" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "r3" "lambda1" "lambda2" "lambda3" "r1sq" "r2sq" "r3sq" "rsquared" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psquared" "U" "Usquared" "AR"], ',');
elseif ARGS[4] == "polar"
  writedlm(outfile, ["E0" "mu" "kT" "Fz" "Fx" "n" "b" "r1" "r2" "r3" "lambda1" "lambda2" "lambda3" "r1sq" "r2sq" "r3sq" "rsquared" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psquared" "U" "Usquared" "AR"], ',');
else
  println("I don't understand the last argument");
  exit(1);
end

for infile in readdir(pattern, indir)
  writedlm(outfile,
           hcat(
                transpose(map(pair -> eval(Meta.parse(split(pair, "-"; limit=2)[2]))*1e-3,
                              split(split(basename(infile), ".")[1], "_"))),
                map(line -> transpose(eval(Meta.parse(split(line, "=")[2]))),
                    readlines(infile))...
               ),
           ',');
end

close(outfile);

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
# The code to generate Fig 2 and SFig 1A
using PyPlot
using CSV
using SurNoR_2020
using Statistics
using MAT

PyPlot.svg(true)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["svg.fonttype"] = "none"
rcParams["pdf.fonttype"] = 42

# ------------------------------------------------------------------------------
# Behavioral data
# ------------------------------------------------------------------------------
BehavData = CSV.read("./data/BehavData_SwitchState.csv")

Sub_Num = 12
Sub_Set = [1,2,3,4,5,6,7,8,9,11,12,13]
Inputs = Array{SurNoR_2020.Str_Input,3}(undef,Sub_Num,2,5)
for n = 1:Sub_Num
  for Block = 1:2
    for Epi = 1:5
      Inputs[n,Block,Epi] = SurNoR_2020.Func_data_to_input(BehavData; Sub=Sub_Set[n],
                                                           Epi=Epi, Block=Block)
    end
  end
end

# ------------------------------------------------------------------------------
# Loading the action classes
# ------------------------------------------------------------------------------
ActionClass = matread("data/ActionClasses.mat")
AC1 = ActionClass["A1"]
AC2 = ActionClass["A2"]

# ------------------------------------------------------------------------------
# Fig 2A Trap-avoidance lenght
# ------------------------------------------------------------------------------
Epi = 1
K_movemean=1
N_thresh = 2

Trap_set = [1,2,6,7]

for Block = 1:2
  y, dy, N, Lenghts, y_med, y_Q25, y_Q75 = SurNoR_2020.Func_trap_av_lenght(Inputs;
                          Trap_set = Trap_set, Block=Block,K_movemean=K_movemean)

  if Block == 1
    Color = "#167CB7"
  else
    Color = "#B71616"
  end

  fig = figure(); ax = gca()
  ax.plot(N)
  ax.set_xlabel("Number of visit")
  ax.set_ylabel("Number of subjects")

  fig = figure(); ax = gca()

  x = 1:length(N)
  x = x[N .>= N_thresh]
  x_lim = length(x)

  dy_med = zeros(length(x),2)
  dy_med[:,1] .= y_med[x] .- y_Q25[x]
  dy_med[:,2] .= y_Q75[x] .- y_med[x]

  ax.errorbar(x,y_med[x],yerr = dy_med',color=Color,capsize=3)
  ax.scatter(x,y_med[x],s=10 .* N,c=Color)
  ax.set_xlim([0, x_lim + 1])

  ax.set_title("Number of actions to escape traps at Nth visit")
  ax.set_xlabel("Number of visit")
  ax.set_ylabel("Number of actions")

  ax.plot(-1:50, ones(52).*2, "--k",alpha=0.5)
  σ = 0.1
  for j = 1:Sub_Num
    y_new = Lenghts[j]
    if length(y_new) > x_lim
      y_new = y_new[1:x_lim]
    end
    x_1 = (1:length(y_new)) .+ 2*σ*(rand() - 0.5)
    ax.plot(x_1, y_new[:],".",color="k",alpha=0.5)
  end

  ax.set_ylim([1,18])
end

# ------------------------------------------------------------------------------
# Fig 2B, Fig 2C, and SFig 1A: Progress measures
# ------------------------------------------------------------------------------
Epi = 1
K_movemean = 1
scores = [-0.75,0.5,1.]
Map_S = [1,8,3,5,2,9,10,4,6,7]

for Block = 1:2
  fig = figure(figsize=(28,28))
  for State = 1:10
    if Block == 1
      AC = AC1
      Color = "#167CB7"
    else
      AC = AC2
      Color = "#B71616"
    end
    y, dy, N, Actions, y_med = SurNoR_2020.Func_action_progress_vec(Inputs, AC;
                                                        Block = Block, Epi = Epi,
                                                        State = State,
                                                        K_movemean=K_movemean,
                                                        scores=scores)

    N_thresh = 2
    ax = subplot(3,4,Map_S[State])

    x = 1:length(N)
    x = x[N .>= N_thresh]

    ax.errorbar(x,y[x],yerr = dy[x],color=Color,capsize=3)
    ax.scatter(x,y[x],s=10 .* N,c=Color)
    ax.set_ylim([-1.1,1.1])

    x_lim = length(x)
    ax.plot(-1:(x_lim+1), -ones(x_lim+3).*1.0, "--k",alpha=0.5)
    ax.plot(-1:(x_lim+1), ones(x_lim+3).*1.0, "--k",alpha=0.5)
    ax.plot(-1:(x_lim+1), ones(x_lim+3).*0.0, "--k",alpha=0.5)

    ax.set_title(string("Average progress at Nth visit - State:", string(Map_S[State])))
    ax.set_ylabel("Average progress")
    ax.set_xlim([0, x_lim + 1])
  end
end


text:
module FieldsBenchs

include("MockFieldsBenchs.jl")

include("ConstantFieldsBenchs.jl")

include("FieldApplyBenchs.jl")

include("FieldArraysBenchs.jl")

include("LincombBenchs.jl")

include("ComposeBenchs.jl")

include("VarinnerBenchs.jl")

include("AttachmapBenchs.jl")

include("IntegrateBenchs.jl")

end # module


text:
using Test

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    include("../src/Classic.jl")
    using .Classic
else
    # invoked during test
    using Classic
end

function test_Nucleotide()
    # Test Nucleotide(x::Char)
    @test Classic.A == Nucleotide('A')
    @test Classic.C == Nucleotide('C')
    @test Classic.G == Nucleotide('G')
    @test Classic.T == Nucleotide('T')
    @test_throws DomainError Nucleotide('Z')
    @test Nucleotide.(['A', 'C', 'G']) == [Classic.A, Classic.C, Classic.G]

    # Test Nucleotide(Int)
    @test Classic.A == Nucleotide(1)
    @test Classic.C == Nucleotide(2)
    @test Classic.G == Nucleotide(3)
    @test Classic.T == Nucleotide(4)
    @test_throws ArgumentError Nucleotide(5)
    @test Nucleotide.([1,2,3]) == [Classic.A, Classic.C, Classic.G]

    # Test Int(Nucleotide)
    @test 1 == Int(Nucleotide('A'))
    @test [1,2] == Int.([Nucleotide('A'), Nucleotide('C')])

    # convert()
    @test ['A', 'C'] == convert(Vector{Char}, [Nucleotide('A'), Nucleotide('C')])
    @test ['A', 'C'] == convert(Vector{Char}, [Classic.A, Classic.C])

    # Test <, ==, etc
    @test Nucleotide(1) == Nucleotide(1)
    @test Nucleotide(1) != Nucleotide(2)
    @test Nucleotide(1) < Nucleotide(2)
    @test Nucleotide(2) > Nucleotide(1)

    # isequal
    @test isequal(Nucleotide(1), Nucleotide(1))

    # hash
    @test hash(Nucleotide(1)) == hash(Nucleotide(1))
end

function test_Codon()
    # Test Codon(String)
    codon = Codon("AGT")
    @test [Classic.A, Classic.G, Classic.T] == codon.nucleotide
    codon = Codon("AGT1@#")
    @test [Classic.A, Classic.G, Classic.T] == codon.nucleotide
    @test_throws DomainError Codon("ZZZ")
    @test_throws BoundsError Codon("AG")

    # Test ==
    @test Codon("AGT") == Codon("AGT")
    @test Codon("ACT") != Codon("AGT")

    # isequal
    @test isequal(Codon("AGT"), Codon("AGT"))
    @test !isequal(Codon("ACT"), Codon("AGT"))

    # hash
    @test hash(Codon("AGT")) == hash(Codon("AGT"))
    @test hash(Codon("ACT")) != hash(Codon("AGT"))

    # <
    @test Codon("ACT") < Codon("AGT")
    @test Codon("AGT") > Codon("ACT")
    # isless
    @test isless(Codon("ACT"), Codon("AGT"))
    @test !isless(Codon("AGT"), Codon("ACT"))
end

function test_Gene()
    gene = Gene("AGTTGC")
    @test 2 == length(gene.codons)
    @test Codon("AGT") == gene.codons[1]
    @test Codon("TGC") == gene.codons[2]

    @test isempty(Gene("").codons)
    @test_throws BoundsError Gene("AGTT")
    @test_throws MethodError Gene()

    # linear search
    gene_str = "ACGTGGCTCTCTAACGTACGTACGTACGGGGTTTATATATACCCTAGGACTCCCTTT"
    @test true == linear_contains(Gene(gene_str), Codon("ACG"))
    @test false == linear_contains(Gene(gene_str), Codon("GAT"))

    # binary search
    gene_str = "ACGTGGCTCTCTAACGTACGTACGTACGGGGTTTATATATACCCTAGGACTCCCTTT"
    @test true == binary_contains(Gene(gene_str), Codon("ACT"))
    @test false == binary_contains(Gene(gene_str), Codon("GAT"))

end

@testset "Gene" begin
    test_Nucleotide()
    test_Codon()
    test_Gene()
end

nothing


text:
#!/usr/bin/env julia
using Plots

# sigmoid function
s(a,x)   = 1/(1+exp(-a*x))

# first derivative of sigmoid function
d1s(a,x) = a*exp(-a*x)/((1+exp(-a*x))^2)

# second derivative of sigmoid function
d2s(a,x) = a * (2*d1s(a,x)*exp(-a*x)/(1+exp(-a*x))  - d1s(a,x))

x  =  range(-5.0,5.0,step=0.01)
s1   =  s.(5,x)
d1s1 =  d1s.(5,x)
#d2s1 =  d2s.(5,x)

plot(x,[s1,d1s1],yscale = :log10)
# ,d2s1]
#label = ["sigmoid", "dsigmoid/dx", "d^2sigmoid/dx^2"])

max_stars_repo_name: joelandman/nlytiq-base

max_issues_repo_name: joelandman/nyltiq-base

max_forks_repo_name: joelandman/nyltiq-base


text:
using StatsBase
using Distributions
using CategoricalArrays
using DataFrames
using CSV
using RatPack

file = "../examples/test_ratings.csv"
(a,b) = read_ratings(file)

file = "../examples/test_competitions.csv"
(c,d) = read_results( file )

R = RatingsList( d )

file = "../data/test_competitions_1.csv"
(e,f) = read_results( file )
summarize(e, f)

input_ratings = RatingsList( f )
input_competitions = e

ratings = update_ratings(UpdateColley(),
                         input_ratings,
                         input_competitions)

ratings = update_ratings(UpdateKeenerScores(; skew = false),
                         input_ratings,
                         input_competitions)

ratings = update_ratings(UpdateElo(;r0 = 0.0),
                         input_ratings,
                         input_competitions)

# try out the NFL data
file = "../data/nfl_2009_regular.csv"
(e,f) = read_results( file )
summarize(e, f)

input_ratings = RatingsList( f )
input_competitions = e

ratings = update_ratings(UpdateKeenerScores(; skew=true, norm=false),
                         input_ratings,
                         input_competitions)

S = sort(collect(ratings.ratings), by = tuple -> last(tuple), rev=true)

# iterate Elo on the data, one by one
file = "../data/nfl_2009.csv"
(nfl_ext_competitions,  nfl_ext_player_list) = read_results( file )
nfl_ext_ratings = RatingsList( nfl_ext_player_list )

rule = UpdateElo(;r0 = 0.0, K = 32.0, θ=1000.0/log(10.0) )
r = nfl_ext_ratings
for i=1:size(nfl_ext_competitions,1)
    global r = update_ratings(rule, r, nfl_ext_competitions[ [i], :])
end
S = sort(collect(r.ratings), by = tuple -> last(tuple), rev=true)

mean(last.(S)) # should be close to r0 = 0.0


text:
@testset "chainrules.jl" begin
    test_rrule(NamedDimsArray, rand(2, 3), (:a, :b); check_inferred=VERSION >= v"1.6")
    test_rrule(NamedDimsArray{(:a, :b)}, (rand(2, 3)); check_inferred=VERSION >= v"1.6")
end


text:
##==============================================================================
## Abstract Types
##==============================================================================

abstract type InferenceType end
abstract type PackedInferenceType end

abstract type FunctorInferenceType <: Function end

# NOTE DF, Convolution is IIF idea, but DFG should know about "FactorOperationalMemory"
# DF, IIF.CommonConvWrapper <: FactorOperationalMemory #
abstract type FactorOperationalMemory <: Function end
# TODO to be removed from DFG,
# we can add to IIF or have IIF.CommonConvWrapper <: FactorOperationalMemory directly
# abstract type ConvolutionObject <: FactorOperationalMemory end

abstract type FunctorSingleton <: FunctorInferenceType end
abstract type FunctorPairwise <: FunctorInferenceType end
abstract type FunctorPairwiseMinimize <: FunctorInferenceType end

##==============================================================================
## GenericFunctionNodeData
##==============================================================================

"""
$(TYPEDEF)

Notes
- S::Symbol

Designing (WIP)
- T <: Union{FactorOperationalMemory, PackedInferenceType}
- in IIF.CCW{T <: DFG.InferenceType}
- in IIF.FunctorPairwiseMinimize <: InferenceType # DFG whatever, something, we'll figure it out
- in Main/User, SomeFactor <: FunctorPairwiseMinimize
"""
mutable struct GenericFunctionNodeData{T<:Union{PackedInferenceType, FunctorInferenceType, FactorOperationalMemory}}
    eliminated::Bool
    potentialused::Bool
    edgeIDs::Vector{Int}
    fnc::T
    multihypo::Vector{Float64} # likely to moved when GenericWrapParam is refactored
    certainhypo::Vector{Int}
    solveInProgress::Int
end

## Constructors

GenericFunctionNodeData{T}() where T =
    GenericFunctionNodeData{T}(false, false, Int[], T(), Float64[], Int[], 0)

function GenericFunctionNodeData(eliminated,
                                 potentialused,
                                 edgeIDs,
                                 fnc,
                                 multihypo=Float64[],
                                 certainhypo=Int[])
    return GenericFunctionNodeData(eliminated, potentialused, edgeIDs, fnc, multihypo, certainhypo, 0)
end

##------------------------------------------------------------------------------
## PackedFunctionNodeData and FunctionNodeData

const PackedFunctionNodeData{T} = GenericFunctionNodeData{T} where T <: PackedInferenceType
PackedFunctionNodeData(args...) = PackedFunctionNodeData{typeof(args[4])}(args...)

const FunctionNodeData{T} = GenericFunctionNodeData{T} where T <: Union{FunctorInferenceType, FactorOperationalMemory}
FunctionNodeData(args...) = FunctionNodeData{typeof(args[4])}(args...)

# PackedFunctionNodeData(x2, x3, x4, x6::T, multihypo::Vector{Float64}=[], certainhypo::Vector{Int}=Int[], x9::Int=0) where T <: PackedInferenceType =
#     GenericFunctionNodeData{T}(x2, x3, x4, x6, multihypo, certainhypo, x9)
# FunctionNodeData(x2, x3, x4, x6::T, multihypo::Vector{Float64}=[], certainhypo::Vector{Int}=Int[], x9::Int=0) where T <: Union{FunctorInferenceType, FactorOperationalMemory} =
#     GenericFunctionNodeData{T}(x2, x3, x4, x6, multihypo, certainhypo, x9)

##==============================================================================
## Factors
##==============================================================================
#
# |                   | label | tags | timestamp | solvable | solverData |
# |-------------------|:-----:|:----:|:---------:|:--------:|:----------:|
# | SkeletonDFGFactor |   X   |   x  |           |          |            |
# | DFGFactorSummary  |   X   |   X  |     X     |          |            |
# | DFGFactor         |   X   |   X  |     X     |     X    |      X     |

## DFGFactor lv2

"""
$(TYPEDEF)
Complete factor structure for a DistributedFactorGraph factor.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGFactor{T, N} <: AbstractDFGFactor
    """Factor label, e.g. :x1f1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref), [`setTimestamp`](@ref)"""
    timestamp::DateTime
    """Nano second time, for more resolution on timestamp (only subsecond information)"""
    nstime::Nanosecond
    """Factor tags, e.g [:FACTOR].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Solver data.
    Accessors: [`getSolverData`](@ref), [`setSolverData!`](@ref)"""
    solverData::Base.RefValue{GenericFunctionNodeData{T}}
    """Solvable flag for the factor.
    Accessors: [`getSolvable`](@ref), [`setSolvable!`](@ref)"""
    solvable::Base.RefValue{Int}
    """Internal cache of the ordering of the neighbor variables. Rather use getVariableOrder to get the list as this is an internal value.
    Accessors: [`getVariableOrder`](@ref)"""
    _variableOrderSymbols::NTuple{N,Symbol}
    # Inner constructor
    DFGFactor{T}(label::Symbol,
                 timestamp::DateTime,
                 nstime::Nanosecond,
                 tags::Set{Symbol},
                 solverData::GenericFunctionNodeData{T},
                 solvable::Int,
                 _variableOrderSymbols::NTuple{N,Symbol}) where {T,N} =
                 new{T,N}(label, timestamp, nstime, tags, Ref(solverData), Ref(solvable), _variableOrderSymbols)

end

##------------------------------------------------------------------------------
## Constructors

"""
$(SIGNATURES)

Construct a DFG factor given a label.
"""
DFGFactor(label::Symbol,
          timestamp::DateTime,
          nstime::Nanosecond,
          tags::Set{Symbol},
          solverData::GenericFunctionNodeData{T},
          solvable::Int,
          _variableOrderSymbols::Tuple) where {T} =
          DFGFactor{T}(label, timestamp, nstime, tags, solverData, solvable, _variableOrderSymbols)

DFGFactor{T}(label::Symbol, variableOrderSymbols::Vector{Symbol}, timestamp::DateTime=now(), data::GenericFunctionNodeData{T} = GenericFunctionNodeData{T}()) where {T} =
                DFGFactor(label, timestamp, Nanosecond(0), Set{Symbol}(), data, 1, Tuple(variableOrderSymbols))

DFGFactor(label::Symbol,
          variableOrderSymbols::Vector{Symbol},
          data::GenericFunctionNodeData{T};
          tags::Set{Symbol}=Set{Symbol}(),
          timestamp::DateTime=now(),
          solvable::Int=1,
          nstime::Nanosecond = Nanosecond(0)) where {T} =
                DFGFactor{T}(label, timestamp, nstime, tags, data, solvable, Tuple(variableOrderSymbols))

Base.getproperty(x::DFGFactor,f::Symbol) = begin
    if f == :solvable || f == :solverData
        getfield(x, f)[]
    elseif f == :_variableOrderSymbols
        [getfield(x,f)...]
    else
        getfield(x,f)
    end
end

Base.setproperty!(x::DFGFactor,f::Symbol, val) = begin
    if f == :solvable || f == :solverData
        getfield(x,f)[] = val
    else
        setfield!(x,f,val)
    end
end
##------------------------------------------------------------------------------
## DFGFactorSummary lv1
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Read-only summary factor structure for a DistributedFactorGraph factor.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGFactorSummary <: AbstractDFGFactor
    """Factor label, e.g. :x1f1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref)"""
    timestamp::DateTime
    """Factor tags, e.g [:FACTOR].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Internal cache of the ordering of the neighbor variables. Rather use getNeighbors to get the list as this is an internal value.
    Accessors: [`getVariableOrder`](@ref)"""
    _variableOrderSymbols::Vector{Symbol}
end

##------------------------------------------------------------------------------
## SkeletonDFGFactor lv0
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Skeleton factor structure for a DistributedFactorGraph factor.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct SkeletonDFGFactor <: AbstractDFGFactor
    """Factor label, e.g. :x1f1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Factor tags, e.g [:FACTOR].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Internal cache of the ordering of the neighbor variables. Rather use getNeighbors to get the list as this is an internal value.
    Accessors: [`getVariableOrder`](@ref)"""
    _variableOrderSymbols::Vector{Symbol}
end

##------------------------------------------------------------------------------
## Constructors

#NOTE I feel like a want to force a variableOrderSymbols
SkeletonDFGFactor(label::Symbol, variableOrderSymbols::Vector{Symbol} = Symbol[]) = SkeletonDFGFactor(label, Set{Symbol}(), variableOrderSymbols)

##==============================================================================
## Define factor levels
##==============================================================================
const FactorDataLevel0 = Union{DFGFactor, DFGFactorSummary, SkeletonDFGFactor}
const FactorDataLevel1 = Union{DFGFactor, DFGFactorSummary}
const FactorDataLevel2 = Union{DFGFactor}

##==============================================================================
## Convertion constructors
##==============================================================================

DFGFactorSummary(f::DFGFactor) =
    DFGFactorSummary(f.label, f.timestamp, deepcopy(f.tags), deepcopy(f._variableOrderSymbols))

SkeletonDFGFactor(f::FactorDataLevel1) =
    SkeletonDFGFactor(f.label, deepcopy(f.tags), deepcopy(f._variableOrderSymbols))


text:
#=
The function getweights() computes, for n in 0:N,
the integral I_n over [-1,1] of
dx * T_n(x) * (isodd(n) ? sin : cos)(w*x) where T_n
is the nth Chebyshev polynomial T_n(cos(x))=cos(n*x).

Supported data types are Float32, Float64, and
BigFloat with precision=256.

We use a three-term recurrence relation
found in, for example, (Evans and Webster, 1999).
The default, method=:thomas, is to use forward recursion
for n<=|w| and then to solve the recurrence for |w|<n<=N+P
for some padding length P using the Thomas tridiagonal solver algorithm,
the true boundary condition for I_{S-1} where S=ceil(|w|),
and the false boundary condition I_{N+P}=0.
The padding length depends only on the data type.

Keyword argument method=:forward uses forward recursion for all n,
even though this is unstable for n>|w|.

For testing purposes, I also implemented a slower algorithm
involving QR decomposition (method=:qr) with pivoting
to solve the recurrence relations with the true boundary condition
for I_0 the false boundary condition I_{N+P}=0.
(QR w/o pivoting and LU, even with pivoting,
are not inaccurate for some choices of w.)
=#
function getweights(N::Integer,w::Real,method::Symbol=:thomas,T::Type=Float64)::Vector
    out,work,paddedN=weights_alloc(N,method,T)
    getweights!(out,work,paddedN,w,method)
end
function weights_alloc(N::Integer,method::Symbol,T::Type=Float64)
    if method == :qr || method == :thomas
        N += padding(N,T)
    end
    out = Vector{T}(undef, N+1)
    if method == :qr
        work = Matrix{T}(undef,N+1,N+2)
    elseif method == :thomas
        work = Matrix{T}(undef,2,N+1)
    else
        work = Matrix{T}(undef,0,0)
    end
    out,work,N
end
function padding(N::Integer,T::Type=Float64)
    if T==Float64
        return 9+2round(Int64,0.5+sqrt(18.5N),RoundUp)
    elseif T==BigFloat #for unit roundoff 2.0^-256
        return 55+2round(Int64,sqrt(89N),RoundUp)
    elseif T==Float32
        return 3+2round(Int64,0.5+sqrt(8.5N),RoundUp)
    else
        throw("Unsupported type for Chebyshev weight computation")
    end
end
function getweights!(out::AbstractArray,work::AbstractArray,N::Integer,w::Real,method::Symbol)
    N>=0 || return out
    T = eltype(work)

    if iszero(w) #Clenshaw-Curtis weights
        out[1] = 2.0 #I_0
        for n in 2:2:N
            out[n] = zero(T) #I_{n-1}
            out[n+1] = T(2) / (1 - n^2) #I_n
        end
        if isodd(N)
            out[N+1] = zero(T) #I_N
        end
        return out
    end

    s,c,rw = sincos(w)...,inv(w)
    out[1] = 2s*rw #wI_0 = 2sin(w)
    N >= 1 || return out

    #slow QR decomposition with pivoting
    method == :qr && return getweightsQR!(out,work,N,w,s,c)

    #stable forward recursion for n < n0
    n0 = min(round(Int64,abs(w),RoundUp),N+1)
    n0 >=2 && getweightsForward!(out,1,n0-1,rw,s,c)
    n0 <= N || return out
    if n0 <= N #use method for n0 <= n <= N
        if method == :thomas
            getweightsThomas!(out,work,n0,N,w,s,c)
        elseif method == :forward
            #unstable forward recursion
            getweightsForward!(out,n0,N,rw,s,c)
        end
    end
    out
end

#=
Evans-Webster (1999) forward-recursive computation of I_n(w) for nmin <= n <= nmax.
Assumes w!=0 and 1 <= nmin <= nmax < length(out).
Assumes out[n+1]=I_n(w) initialized for max(0,nmin-2) <= n < nmin.
WARNING: unstable if nmax > |w|.
=#
function getweightsForward!(out::AbstractArray,nmin::Integer,nmax::Integer,
                            invw::Real,sinw::Real,cosw::Real)
    @assert 1 <= nmin <= nmax < length(out)

    if nmin <= 2
        if nmin <= 1
            out[2] = (out[1] - 2cosw)*invw #I_0 - wI_1 = 2cos(w)
            2 <= nmax || return out
        end
        out[3] = (2sinw - 4out[2])*invw #4I_1 + wI_2 = 2sin(w)
        3 <= nmax || return out
    end
    nstart=max(nmin,3)

    s4w,c4w = -4sinw*invw,4cosw*invw

    #n odd: wI_{n-2}/(2(n-2)) + I_{n-1} - wI_n/(2n) = -2cos(w)/(n(n-2))
    #out[n+1] = I_n = (nI_{n-2} + 4cos(w)/w)/(n-2) + 2nI_{n-1}/w
    getodd(n) = (n*out[n-1]+c4w)/(n-2) + 2n*out[n]*invw

    #n even: wI_{n-2}/(2(n-2)) - I_{n-1} - wI_n/(2n) = 2sin(w)/(n(n-2))
    #out[n+1] = I_n = (nI_{n-2} - 4sin(w)/w)/(n-2) - 2nI_{n-1}/w
    geteven(n) = (n*out[n-1]+s4w)/(n-2) - 2n*out[n]*invw

    isodd(nstart) && (out[nstart+1] = getodd(nstart))
    even_start = ((nstart+1)>>1)<<1 #round up
    for n in even_start:2:nmax-1
        out[n+1] = geteven(n)
        out[n+2] = getodd(n+1)
    end
    iseven(nmax) && (out[nmax+1] = geteven(nmax))
    out
end

#=
Evans-Webster (1999): Olver's method for I_n(w) for nmin <= n <= nmax,
specifically using Thomas' tridiagonal solver algorithm.
Assumes that 1 <= nmin <= nmax < length(out),
that work has at least 2 rows and nmax-nmin+1 columns, and
that out[nmin]=I_{nmin-1}(w) is initialized.
WARNING: unstable if nmin < |w|.
WARNING: uses Olver's method, making a tail of the output inaccurate.
=#
function getweightsThomas!(out::AbstractArray,work::AbstractArray,nmin::Integer,nmax::Integer,
                           w::Real,sinw::Real,cosw::Real)
    T = eltype(work)
    @assert 1 <= nmin <= nmax < length(out)
    sup = view(work,1,:) #superdiagonal after row operations
    rhs = view(work,2,:) #righthand side after row operations
    #1 == n odd:                      I_1 + wI_{1+1}/(2(1+1)) = 2sin(w)/4
    #2 <= n odd: -wI_{n-1}/(2(n-1)) + I_n + wI_{n+1}/(2(n+1)) = 2sin(w)/(1-n^2)
    #2 <= n even: wI_{n-1}/(2(n-1)) + I_n - wI_{n+1}/(2(n+1)) = 2cos(w)/(1-n^2)
    if nmin < 2
        sup[1] = sup_n = 0.25w
        rhs[1] = rhs_n = 0.5sinw
    else
        sup_n, rhs_n = zero(T), out[nmin]
    end
    p,q = 4cosw,-4sinw
    sign = 1
    nstart = max(nmin,2)
    if isodd(nstart)
        p,q = q,p
        sign = -sign
    end
    for n in nstart:nmax
        r = inv((n+1)*(2sign*(1-n) + w*sup_n))
        i = n-nmin+1
        sup[i] = sup_n = r*(n-1)*w;
        rhs[i] = rhs_n = r*(p + (n+1)*w*rhs_n);
        sign = -sign
        p,q = q,p
    end
    solution_n = zero(T)
    for n in nmax:-1:nmin
        i = n-nmin+1
        out[n+1] = solution_n = rhs[i] - sup[i]*solution_n
    end
    out
end

#=
Assumes N>=1, w!=0, that work has >=N+1 rows and >=N+2 columns, and that length(out)>=N+1.
Evans-Webster (1999) computation of I_n(w) for 0 <= n <= N,
but using QR with pivoting to solve the system. Stable but slow.
Assumes w!=0 and 2 <= N < length(out).
WARNING: uses Olver's method, making a tail of the output inaccurate.
=#
function getweightsQR!(out::AbstractArray,work::AbstractArray,N::Integer,
                       w::Real,sinw::Real,cosw::Real)
    T = eltype(work)
    #system of equations A*out = b; out=[I_0,...,I_N]
    A = view(work,1:N+1,1:N+1)
    A .= zero(T)
    b = view(work,1:N+1,N+2)

    #wI_0 = 2sin(w)
    #I_0 - wI_1 = 2cos(w)
    #4I_1 + wI_2 = 2sin(w)
    s2 = 2sinw
    c2 = 2cosw
    negw = -w
    A[1,1], b[1] = w, s2
    A[2,1], A[2,2], b[2] = 1, negw, c2
    A[3,2], A[3,3], b[3] = 4, w, s2

    for n in 2:2:N-2
        #n even: wI_{n-1}/(2(n-1)) + I_n - wI_{n+1}/(2(n+1)) = 2cos(w)/(1-n^2)
        #n+1 odd: -wI_n/(2n) + I_{n+1} + wI_{n+2}/(2(n+2)) = 2sin(w)/(1-(n+1)^2)
        A[n+2,n] = w/(2(n-1))
        A[n+2,n+1] = 1
        A[n+2,n+2] = negw/(2(n+1))
        b[n+2] = c2/(1-n^2)
        A[n+3,n+1] = negw/(2n)
        A[n+3,n+2] = 1
        A[n+3,n+3] = w/(2(n+2))
        b[n+3] = s2/(1-(n+1)^2)
    end
    if isodd(N)
        #N-1 even: wI_{N-2}/(2(N-2)) + I_{N-1} - wI_N/(2N) = 2cos(w)/(1-(N-1)^2)
        A[N+1,N-1] = w/(2(N-2))
        A[N+1,N] = 1
        A[N+1,N+1] = negw/(2N)
        b[N+1] = c2/(N*(2-N))
    end
    #println(cond(A))
    ldiv!(out,qr!(A,Val(true)),b) #QR with pivoting helps in ill-conditioned cases
    out
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
# Rationals Test
# Choose uniform collocation points in a cardinal basis (polynomial)
# and compute the derivative and the integral
#--------------------------------------------------------------------

function collocation{T<:Int}(i::T, N::T)::Rational
    return -(-1 + (2*(i-1)//N))
end

function derivative{T<:Int}(i::T, j::T, N::T)::Rational
    x = Rational[collocation(k, N) for k in 1:N+1]
    if i == j
        return sum((k==j ? 0 : (x[j] - x[k])^(-1)) for k in 1:N+1)
    else
        ai = prod((k==i ? 1 : (x[i] - x[k])) for k in 1:N+1)
        aj = prod((k==j ? 1 : (x[j] - x[k])) for k in 1:N+1)
        return ai/(aj*(x[i] - x[j]))
    end
end

N    = 5
grid = Rational[collocation(k, N) for k in 1:N+1]
D    = Rational[derivative(i, j, N) for i in 1:N+1, j in 1:N+1]
D2   = similar(D)

for i in 1:N+1, j in 1:N+1
    D2[i,j] = sum(D[i,k]*D[k,j] for k in 1:N+1)
end

#--------------------------------------------------------------------
# Test simple derivative
#--------------------------------------------------------------------
ϕ    = grid.^2
ψ    = similar(ϕ)

for i in 1:N+1
    ψ[i] = sum(D[i,k]*ϕ[k] for k in 1:N+1)
end

@test ψ == 2.*grid

#--------------------------------------------------------------------
# Solve a 1D elliptic PDE
#--------------------------------------------------------------------
x = grid
B = rationalize.(zeros(N+1, N+1))
B[1,1] = B[end, end] = 1//1
b = rationalize.(zeros(N+1))
b[1]   =  1//20
b[end] = -1//20
u = (D2 + B) \ (x.^3 + b)
@test u == x.^5/20

#--------------------------------------------------------------------
# Solve a 2D elliptic PDE
#--------------------------------------------------------------------

I   = rationalize.(eye(N+1))
Dx  = rationalize.(zeros((N+1), (N+1), (N+1), (N+1)))
Dy  = similar(Dx)
Φ   = Rational[xx^3 + yy^2 for xx in x, yy in x]
dΦx = similar(Φ)
dΦy = similar(Φ)

for index in CartesianRange(size(Dx))
    i,j,k,l   = index.I
    Dx[index] = D[j,l]*I[i,k]
    Dy[index] = I[j,l]*D[i,k]
end

for index in CartesianRange(size(Φ))
    i, ii      = index.I
    dΦx[index] = sum(Dy[i, ii, j, jj]*Φ[j, jj] for j in 1:N+1, jj in 1:N+1)
    dΦy[index] = sum(Dx[i, ii, j, jj]*Φ[j, jj] for j in 1:N+1, jj in 1:N+1)
end

@test dΦy == Rational[2*yy for xx in x, yy in x]
@test dΦx == Rational[3*xx^2 for xx in x, yy in x]


text:
using Documenter

makedocs(
    sitename= "CellMLToolkit",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md"
    ]
)


text:
#这种注释方式，不推荐大家使用，除非是需要用到数学公式之类的。
#因为此注释不能被julia插件解析，在vscode编辑器中无法悬停查看函数帮助。
@doc raw"""
    pt(a,b)

勾股定理，英文名为Pythagoras 也称为毕达哥拉斯定理。

在平面上的一个直角三角形中，两个直角边边长的平方加起来等于斜边长的平方。
如果设直角三角形的两条直角边长度分别是`a`和`b`，斜边长度是`c`，那么数学公式为：

``{\rm{c = }}\sqrt {{a^2} + {b^2}}``

返回 斜边长`c`
"""
function pt(a,b)
    c = sqrt(a^2+b^2)
end


text:
# this follows `markdown.jl` but spurred by bugs/issues

function inter(st::String)
    steps = explore_md_steps(st)
    return steps[:inter_md].inter_md, steps[:inter_html].inter_html
end

@testset "issue163" begin
    st = raw"""A _B `C` D_ E"""
    imd, ih = inter(st)
    @test imd == "A _B  ##FDINSERT## D_ E"
    @test isapproxstr(ih, "
A B  ##FDINSERT##  D E
")

    st = raw"""A _`B` C D_ E"""
    imd, ih = inter(st)
    @test isapproxstr(imd, "A _ ##FDINSERT##  C D_ E")
    @test isapproxstr(ih, "
A ##FDINSERT##  C D E
")

    st = raw"""A _B C `D`_ E"""
    imd, ih = inter(st)
    @test isapproxstr(imd, "A _B C  ##FDINSERT## _ E")
    @test isapproxstr(ih, "
A B C  ##FDINSERT## E
")

    st = raw"""A _`B` C `D`_ E"""
    imd, ih = inter(st)
    @test isapproxstr(imd, "A _ ##FDINSERT##  C  ##FDINSERT## _ E")
    @test isapproxstr(ih, "
A ##FDINSERT##  C  ##FDINSERT## E
")
end

@testset "TOC"  begin
    fs1()
    h = raw"""
        @def fd_rpath = "pages/ff/aa.md"
        \toc
        ## Hello `fd`
        #### weirdly nested
        ### Goodbye!
        ## Done
        done.
        """ |> seval
    @test isapproxstr(h, raw"""

              Hello fd

                      weirdly nested

                Goodbye!

            Done

Hello fd

weirdly nested

Goodbye!

Done

done.

        """)
end

@testset "TOC"  begin
    fs1()
    s = raw"""
        @def fd_rpath = "pages/ff/aa.md"
        @def mintoclevel = 2
        @def maxtoclevel = 3
        \toc
        # A
        ## B
        #### C
        ### D
        ## E
        ### F
        done.
        """ |> seval
    @test isapproxstr(s, raw"""

B

D

E

F

A

B

C

D

E

F

done.

        """)
end

@testset "TOC-fs2"  begin
    fs2()
    h = raw"""
        @def fd_rpath = "pages/ff/aa.md"
        \toc
        ## Hello `fd`
        #### weirdly nested
        ### Goodbye!
        ## Done
        done.
        """ |> seval
    @test isapproxstr(h, raw"""

              Hello fd

weirdly nested

Goodbye!

Done

Hello fd

weirdly nested

Goodbye!

Done

done.

        """)
end

@testset "TOC-fs2"  begin
    fs2()
    s = raw"""
        @def fd_rpath = "pages/ff/aa.md"
        @def mintoclevel = 2
        @def maxtoclevel = 3
        \toc
        # A
        ## B
        #### C
        ### D
        ## E
        ### F
        done.
        """ |> seval
    @test isapproxstr(s, raw"""

B

D

E

F

A

B

C

D

E

F

done.

        """)
end


text:
a_1 = 5
a_2 = -3

if (a_1 == 5) && (a_2 == -3)
a_1 = 10
end

println(a_1)


text:
using HDF5
using SparseArrays
using Statistics
include("../convergence.jl")
include("../hmc.jl")
include("../models.jl")

function main()
    R = 1001
    control = Dict(:iterations => 2000, :chains => 4)
    M = sparse([1.0 0.5; 0.5 4.0])
    d = 32

    S = Matrix(blockdiag([M for i in 1:d÷2]...))
    f = mvgaussian_(zeros(d), S)
    metric = Matrix(Diagonal(ones(d)))
    essbulk = zeros(R)
    esstail = zeros(R)
    esssq = zeros(R)
    leapfrog = zeros(R)
    time = zeros(R)
    Qs = zeros(R, 3)
    for r in 1:R
        time[r] = @elapsed samples, c = stan(f, d; M = metric, control = control)
        s = samples[:, :, 2]
        essbulk[r] = ess_bulk(s)
        esstail[r] = ess_tail(s)
        esssq[r] = ess_sq(s)
        leapfrog[r] = sum(c[:leapfrog])
        draws = reshape(s, :, 2)
        Qs[r, :] = quantile(draws[:, 2], [0.1 0.5 0.9])
    end
    g = "anisotropicgaussianC25/yastan/"
    h5open("data.h5", "w") do file
        write(file, g * "essbulk", essbulk)
        write(file, g * "esstail", esstail)
        write(file, g * "esssq", esssq)
        write(file, g * "leapfrog", leapfrog)
        write(file, g * "q10", Qs[:, 1])
        write(file, g * "q50", Qs[:, 2])
        write(file, g * "q90", Qs[:, 3])
        write(file, g * "seconds", time)
    end
end

main()

max_stars_repo_name: roualdes/yastan

max_issues_repo_name: roualdes/yastan

max_forks_repo_name: roualdes/yastan