text:
@testset "example2" begin

my_graph = "graph [
  node [id 1]
  node [id 2]
  node [id 3]
  edge [source 1 target 2]
  edge [source 2 target 3]
  edge [source 3 target 1]
]"

root = parse_dict(my_graph)

for graph in root[:graph]  # there could be multiple graphs
    for node in graph[:node]
        println("node $(node[:id])")
    end
    for edge in graph[:edge]
        println("edge $(edge[:source]) - $(edge[:target])")
    end
end

end


text:
using Mosek
using Polyopt
using JuMP
using LinearAlgebra

function elimination_equality(data_a)
 I,J,K,L,AI,AJ,AL,C_I,C_J,C_L,lCI,lCL,lCJ,UI,UJ,UL,Mu_max,Mu_min, Lambda,costI,costL,costJ,Demandcost=data_a;

 Ma_lamb=maximum(Lambda);
 Mi_lamb=minimum(Lambda);
 MO=maximum(C_J);
 UIindex=getindex.(findall(!iszero,1 ./ UI),2)
 ULindex=getindex.(findall(!iszero,1 ./ UL),2);
 UJindex=getindex.(findall(!iszero,1 ./ UJ),2);
 SizeUIindex=size(UIindex);
 SizeUJindex=size(UJindex);
 SizeULindex=size(ULindex);
 CIindex=getindex.(findall(!iszero,1 ./ C_I),2);
 CLindex=getindex.(findall(!iszero,1 ./ C_L),2);
 CJindex=getindex.(findall(!iszero,1 ./ C_J),2);
 SizeCIindex=size(CIindex);
 SizeCJindex=size(CJindex);
 SizeCLindex=size(CLindex);
 lCIindex=getindex.(findall(!iszero,lCI),2);
 lCLindex=getindex.(findall(!iszero,lCL),2);
 lCJindex=getindex.(findall(!iszero,lCJ),2);

 Poolinputsize=zeros(L,1); # number of inputs feed each pool
 Poolinput=Array{Any}(undef,L); # contains inputs' indeces that feed pool l, for each pool l.

 for l in 1:L
    a=findall(!iszero,AI[:,l]);
    Poolinput[l]=collect(a) ; #a cell whose elements are a vector including pools l's inputs

    Poolinputsize[l,:] .= Int(size(a,1));
 end

 Pooloutputsize=zeros(L,1); # number of outputs for each pool
 Pooloutput=Array{Any}(undef,L); # contains outputs' indeces of pool l, for each pool l.

 for l in 1:L
    a=findall(!iszero,AL[l,:]);
    Pooloutput[l]=a ; #a cell whose elements are a vector including pools' outputs
    Pooloutputsize[l,:] .= Int(size(a,1));
 end

 Inputoutputsize=zeros(I,1);
 Inputoutput=Array{Any}(undef,I);

 for i in 1:I
    a=findall(!iszero,AJ[i,:]);
    Inputoutput[i]=collect(a) ; #a cell whose elements are a vector including inputs' outputs
    Inputoutputsize[i,:] .= Int(size(a,1));
 end

 n=L*K+Int(sum(Pooloutputsize))+
      Int(sum(Inputoutputsize))+Int(sum(Poolinputsize))-L*(K+1);

 y = variables("y:\n",n);
 #y=collect(1:n);
 yout=y[1:Int(sum(Pooloutputsize))];# pool-output
 # print("yout:\n",yout);
 yin=y[Int(sum(Pooloutputsize))+1:Int(sum(Pooloutputsize))+Int(sum(Inputoutputsize))]; # input-output
 # print("yin:\n",yin);
 assignedNumber=Int(sum(Pooloutputsize))+Int(sum(Inputoutputsize));
 # print("assigned number sofar: \n", assignedNumber);
 # making the matrices in the polynomial{Float64,1} format
 z=variables("z" ,maximum([L*K,L*J,I*L,I*J]));
 #z=zeros(1,maximum([L*K,L*J,I*L,I*J]));

 P_lk=reshape(0.1*z[1:L*K],L,K);
 Y_lj=reshape(0.1*z[1:L*J],L,J);
 Y_il=reshape(0.1*z[1:I*L],I,L);
 Y_ij=reshape(0.1*z[1:I*J],I,J);

 #constructing the input-output variables as Y_ij
   for i=1:I
    for j=1:Int(Inputoutputsize[i])
      h=Inputoutput[i];
      Y_ij[i,h[j]]=yin[Int(sum(Inputoutputsize[1:i-1]))+j];
    end
    for j=1:J
      if !(j in Inputoutput[i])
        Y_ij[i,j]=0;
      end
    end
    if i==I
       # printstyled("=================Input-output variables are build!=================\n",color=:yellow)
       printstyled("=================Input-output variables are build!=================\n",color=:yellow)
    end
   end

 for l=1:L
  #constructing the pool-output variables as Y_lj
    for j=1:Int(Pooloutputsize[l])
      h=Pooloutput[l];
      Y_lj[l,h[j]]=yout[Int(sum(Pooloutputsize[1:l-1]))+j];
    end
    for j=1:J
      if !(j in Pooloutput[l])
        Y_lj[l,j]=0;
      end
    end
    if l==L
         # printstyled("=================Pool-output variables are build!=================\n",color=:yellow)
         printstyled("=================Pool-output variables are build!=================\n",color=:yellow)
    end
 end

 for l=1:L
  if Poolinputsize[l]>=K+1
    # the elimination will be on the input-pools variables when the incoming arcs are enough
    p=y[assignedNumber+1:assignedNumber+K];
    yinpool=y[assignedNumber+K+1: assignedNumber+K+Int(Poolinputsize[l])-(K+1)]; # input-pool
    assignedNumber=assignedNumber+K+Int(Poolinputsize[l])-(K+1);
    # print("assigned number sofar: \n", assignedNumber);
    #constructing the concentration variables as p_lk
      P_lk[l,:]=p[:];
      if l==L
        # printstyled("=================Concentration variables are build!=================\n",color=:yellow)
        printstyled("=================Concentration variables are build!=================\n",color=:yellow)
      end
    #constructing the input-pool variables as Y_il
    h=Lambda[Poolinput[l],:];
    A=[ones(1,Int(Poolinputsize[l]));
        h'] ; # A is the coeffiecient matrix of the input flow variables
       # print(A)
    (Q,R)=qr(A); #Q is full rank (K+1)*(K+1) and R is uppertriangular (K+1)*Inputsize(l)
    #Solving elimination (25) in Marandi, de Klerk, Dahl's paper
    Q =Q*Matrix(LinearAlgebra.I,K+1,K+1);

    h=sum(Y_lj[l,:]);
    Pvec=[];
    for k=1:K
      Pvec=[Pvec;h*((Ma_lamb-Mi_lamb)*P_lk[l,k]+Mi_lamb)];
    end
    Pvec=[h;Pvec];
      # R[:,1:K+1] yinpool[sum(Poolinputsize[1:l-1])+1:sum(Poolinputsize[1:l-1])+K+1]=h*Q'*[1;P_lk[l,:]]-R[:,K+2:end]*[Y_il[k+1:end]]
     #we use pseudoinverse for elimination
     A_plas=pinv(R[:,1:K+1]);

     y_help=A_plas*(Q'*Pvec - R[:,K+2:end]*yinpool);

      hh=Poolinput[l];
      for i=1:Int(Poolinputsize[l]) #elimination of the first K+1 feeding inputs
        if i<=K+1
          Y_il[Int(hh[i]),l]=y_help[i];
         #h itself is Any but h[1] is Poly. Also, right side is float
        else
          Y_il[Int(hh[i]),l]=yinpool[i-K-1];
        end
      end
      for i=1:I
        if !(i in hh)
          Y_il[i,l]=0;
        end
      end
    if l==L
      # printstyled("==========Input-pool variables are build, considering elimination of the equality constraints!=============\n \n",color=:yellow)
      printstyled("==========Input-pool variables are build, considering elimination of the equality constraints!=============\n \n",color=:yellow)
    end
  else
     # the elimination will be on the input-pools and concentration variables
    p=y[assignedNumber+1:assignedNumber+Int(Poolinputsize[l])-1];
    # yinpool=y[assignedNumber+K+1: assignedNumber+K+Int(Poolinputsize[l])-(K+1)]; # input-pool
    assignedNumber=assignedNumber+Int(Poolinputsize[l])-1;
    #Solving elimination (25) in Marandi, de Klerk, Dahl's paper
    #constructing the concentration variables as p_lk
    for k=1:Int(Poolinputsize[l])-1
      P_lk[l,k]=p[k];
    end
    h=Lambda[Poolinput[l],:];
    A=[ones(1,Int(Poolinputsize[l]));
        h'] ; # A is the coeffiecient matrix of the input flow variables
       # print(A)
    Q,R=qr(A); #Q is full rank (K+1)*(K+1) and R is uppertriangular (K+1)*Inputsize(l)

    # println(Matrix{Int}(LinearAlgebra.I,K+1,K+1))
    Q= Q * Matrix(LinearAlgebra.I,K+1,K+1);

    Q= Q';

    # 0=Q'*[1; P_lk]

    h=-pinv(Q[Int(Poolinputsize[l])+1:K+1,Int(Poolinputsize[l])+1:K+1])*Q[Int(Poolinputsize[l])+1:K+1,1:Int(Poolinputsize[l])]*[1/(Ma_lamb-Mi_lamb);p+(Mi_lamb/(Ma_lamb-Mi_lamb))*ones(Int(Poolinputsize[l])-1,1)]-(Mi_lamb/(Ma_lamb-Mi_lamb))*ones(-Int(Poolinputsize[l])+K+1,1);
    for k=Int(Poolinputsize[l]):K
      P_lk[l,k]=h[k-Int(Poolinputsize[l])+1];
    end
    if l==L
        # printstyled("=================Concentration variables are build!=================\n",color=:yellow)
        printstyled("=================Concentration variables are build!=================\n",color=:yellow)
    end
    h=sum(Y_lj[l,:]);
    Pvec=[];
    for k=1:K
      Pvec=[Pvec;h*((Ma_lamb-Mi_lamb)*P_lk[l,k]+Mi_lamb)];
    end
    Pvec=[h;Pvec];
      # R[:,1:K+1] yinpool[sum(Poolinputsize[1:l-1])+1:sum(Poolinputsize[1:l-1])+K+1]=h*Q'*[1;P_lk[l,:]]-R[:,K+2:end]*[Y_il[k+1:end]]
     #we use pseudoinverse for elimination
     A_plas=pinv(R);
      y_help=A_plas*Q[1:Int(Poolinputsize[l]),:]*Pvec;
      hh=Poolinput[l];
      for i=1:Int(Poolinputsize[l]) #elimination of the first K+1 feeding inputs
          Y_il[Int(hh[i]),l]=y_help[i];
      end
      for i=1:I
        if !(i in hh)
          Y_il[i,l]=0;
        end
      end
    if l==L
      printstyled("==========Input-pool variables are build, considering elimination of the equality constraints!=============\n \n",color=:yellow)
    end
  end
 end

 printstyled("==========The model is being constructed ... =============\n \n",color=:yellow)

 f= MO * ( sum(costI .* Y_il) +
   sum(costL .* Y_lj)+
   sum(costJ .* Y_ij)-
   sum((sum(Y_ij,dims=1) + sum(Y_lj, dims=1) ) * Demandcost) #the first sum is to convert one-lenght array to a polynomial
   );
 printstyled("==========Objective function is constructed! =============\n \n",color=:yellow);

 SizeCIindex=size(CIindex,1);
 SizeCJindex=size(CJindex,1);
 SizeCLindex=size(CLindex,1);
 SizelCIindex=size(lCIindex,1);
 SizelCJindex=size(lCJindex,1);
 SizelCLindex=size(lCLindex,1);

 g=Polyopt.Poly{Int64}[];
 #capacity ristrictions
 for i=1:SizeCIindex

  g=[g;
    1-(MO/ C_I[CIindex[i]])*(sum(Y_il[CIindex[i],:])+
                sum(Y_ij[CIindex[i],:]))];
 end
 for i=1:SizelCIindex
       g=[g;
           -(lCI[lCIindex[i]])/ (MO*(L+J))+(1/(L+J))*(sum(Y_il[lCIindex[i],:])+
                           sum(Y_ij[lCIindex[i],:])
                          )];
 end
 for l=1:SizeCLindex
       g=[g;
           1-MO/(C_L[CLindex[l]])*(sum(Y_lj[CLindex[l],:])) ];
 end

 for l=1:SizelCLindex
       g=[g;
           -(lCL[lCLindex[l]])/(MO*(J))+(1/J)*(sum(Y_lj[lCLindex[l],:])) ];
 end

 for j=1:SizeCJindex
      # print("CJindex:", CJindex[j]);
      # print("\n")
       g=[g;
           1-MO/(C_J[CJindex[j]])*(sum(Y_ij[:,CJindex[j]])+
                         sum(Y_lj[:,CJindex[j]]))];
 end
 for j=1:SizelCJindex
       g=[g;
           -(lCJ[lCJindex[j]])/(MO*(I+L))+(1/(I+L))*(sum(Y_ij[:,CJindex[j]])+
                                           sum(Y_lj[:,CJindex[j]]))];
 end

 #MU_max constraint
 for j=1:J
    for k=1:K
      # if Mi_lamb>=0
        g=[g;
            (1/(max(Mu_max[j,k],0)*(I+L)-sum(min.(Lambda[:,k],0))-min(0,Mi_lamb)*L))*(Mu_max[j,k]*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))-
                                        (sum(Lambda[:,k].*Y_ij[:,j])+
                                          (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
                                                    Mi_lamb*sum(Y_lj[:,j])
                                        ))];
      # elseif Mu_max[j,k]>=0
      #           g=[g;
      #       (1/((Mu_max[j,k])*(I+L)-sum(min(Lambda[:,k],0))-Mi_lamb*J))*(Mu_max[j,k]*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))-
      #                                   (sum(Lambda[:,k].*Y_ij[:,j])+
      #                                     (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
      #                                               Mi_lamb*sum(Y_lj[:,j])
      #                                   ))];
      # else
      #     g=[g;
      #       (1/(-sum(min(Lambda[:,k],0))-Mi_lamb*J))*(Mu_max[j,k]*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))-
      #                                   (sum(Lambda[:,k].*Y_ij[:,j])+
      #                                     (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
      #                                               Mi_lamb*sum(Y_lj[:,j])
      #                                   ))];
      # end
    end

 end
 #Mu_min constraint
  for j=1:J
    for k=1:K
        help=1/Mu_min[j,k];
        if help!=0
            # if min( Mu_min[j,k],Mi_lamb)>=0
             g=[g;
               (1/(sum(max.(Lambda[:,k],0))+L*(Ma_lamb-Mi_lamb)+L*max(0,Mi_lamb)-min(Mu_min[j,k],0)*(I+L)))*(sum(Lambda[:,k].*Y_ij[:,j])+
                            (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
                             Mi_lamb*sum(Y_lj[:,j])-
                             (Mu_min[j,k])*(sum(Y_ij[:,j])+
                                      sum(Y_lj[:,j])
                                      )
                             )];
            # elseif
            #     g=[g;
            #    (1/((Ma_lamb*I+L-Mu_min[j,k]*(I+L))))*(sum(Lambda[:,k].*Y_ij[:,j])+
            #                         (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
            #                         Mi_lamb*sum(Y_lj[:,j])-
            #                         (Mu_min[j,k])*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))
            #                         )];
            # end
        end
    end
  end

  #sign constraints
  for l=1:L
     for k=1:K
        g=[g;
            P_lk[l,k]];
      end
  end
  for i=1:I
     for l=1:L
        if(AI[i,l]==1)
          g=[g;
            (1/((J+(I-1))))*Y_il[i,l]];
          end
      end
  end

  for j=1:J
     for l=1:L
        if(AL[l,j]==1)
        g=[g;
            (1/(I+(J-1)))*Y_lj[l,j]];
        end
     end
  end
  for j=1:J
      for i=1:I
          if(AJ[i,j]==1)
          g=[g;
              Y_ij[i,j]];
         end
      end
  end
  #pipe capacity restriction
  g=[g;
       (ones(SizeUIindex[1],1)-((MO./UI[UIindex]).*Y_il[UIindex]))[:,1]];

   g=[g;
       (ones(SizeULindex[1],1)-((MO./UL[ULindex]).*Y_lj[ULindex]))[:,1]];

  g=[g;
       (ones(SizeUJindex[1],1)-((MO./UJ[UJindex]).*Y_ij[UJindex]))[:,1]];

  printstyled("====================Constraints are constructed!=============\n",color=:yellow)
	f,g ,n
end

function with_equality(data_a)
  I,J,K,L,AI,AJ,AL,C_I,C_J,C_L,lCI,lCL,lCJ,UI,UJ,UL,Mu_max,Mu_min, Lambda,costI,costL,costJ,Demandcost=data_a;
 Ma_lamb=maximum(Lambda);
 Mi_lamb=minimum(Lambda);
 MO=maximum(C_J);

 UIindex=getindex.(findall(!iszero,1 ./ UI),2);
 ULindex=getindex.(findall(!iszero,1 ./ UL),2);
 UJindex=getindex.(findall(!iszero,1 ./ UJ),2);
 SizeUIindex=size(UIindex);
 SizeUJindex=size(UJindex);
 SizeULindex=size(ULindex);
 CIindex=getindex.(findall(!iszero,1 ./ C_I),2);
 CLindex=getindex.(findall(!iszero,1 ./ C_L),2);
 CJindex=getindex.(findall(!iszero,1 ./ C_J),2);

 SizeCIindex=size(CIindex);
 SizeCJindex=size(CJindex);
 SizeCLindex=size(CLindex);
 lCIindex=getindex.(findall(!iszero,lCI),2);
 lCLindex=getindex.(findall(!iszero,lCL),2);
 lCJindex=getindex.(findall(!iszero,lCJ),2);

 Poolinputsize=zeros(L,1); # number of inputs feed each pool
 Poolinput=Array{Any}(undef,L); # contains inputs' indeces that feed pool l, for each pool l.

 for l in 1:L
    a=findall(!iszero,AI[:,l]);
    Poolinput[l]=collect(a) ; #a cell whose elements are a vector including pools l's inputs
    Poolinputsize[l,:] .=Int(size(a,1));
 end

 Pooloutputsize=zeros(L,1); # number of outputs for each pool
 Pooloutput=Array{Any}(undef,L); # contains outputs' indeces of pool l, for each pool l.

 for l in 1:L
    a=findall(!iszero,AL[l,:]);
    Pooloutput[l]=a ; #a cell whose elements are a vector including pools' outputs
    Pooloutputsize[l,:] .=Int(size(a,1));
 end

 Inputoutputsize=zeros(I,1);
 Inputoutput=Array{Any}(undef,I);
 for i in 1:I
    a=findall(!iszero,AJ[i,:]);
    Inputoutput[i]=collect(a) ; #a cell whose elements are a vector including inputs' outputs
    Inputoutputsize[i,:] .=Int(size(a,1));
 end

 n=L*K+Int(sum(Pooloutputsize))+
            Int(sum(Inputoutputsize))+Int(sum(Poolinputsize));

 y=variables("y" ,n);
 #y=collect(1:n);
 yout=y[1:Int(sum(Pooloutputsize))];# pool-output
 yin=y[Int(sum(Pooloutputsize))+1:Int(sum(Pooloutputsize))+Int(sum(Inputoutputsize))]; # input-output
 assignedNumber=Int(sum(Pooloutputsize))+Int(sum(Inputoutputsize));
 ypool=y[assignedNumber+1:assignedNumber+L*K]; # P_lk
 assignedNumber=assignedNumber+L*K;
 yinputpool=y[assignedNumber+1:assignedNumber+Int(sum(Poolinputsize))];
 # println("OK")
 #  making the matrices in the polynomial{Float64,1} format
 z=variables("z" ,maximum([L*K,L*J,I*L,I*J]));
 #z=zeros(1,maximum([L*K,L*J,I*L,I*J]));

 P_lk=reshape(0.1*z[1:L*K],L,K);
 Y_lj=reshape(0.1*z[1:L*J],L,J);
 Y_il=reshape(0.1*z[1:I*L],I,L);
 Y_ij=reshape(0.1*z[1:I*J],I,J);

   P_lk_index=zeros(L,K);
   Y_lj_index=zeros(L,J);
   Y_il_index=zeros(I,L);
   Y_ij_index=zeros(I,J);
 #constructing the input-output variables as Y_ij
     for i=1:I
        for j=1:Int(Inputoutputsize[i])
            h=Inputoutput[i];
            Y_ij[i,h[j]]=yin[Int(sum(Inputoutputsize[1:i-1]))+j];
            Y_ij_index[i,h[j]]=sort(indexin([yin[Int(sum(Inputoutputsize[1:i-1]))+j]],y))[1];
        end
        for j=1:J
            if !(j in Inputoutput[i])
                Y_ij[i,j]=0;
            end
        end
        if i==I
             printstyled("=================Input-output variables are build!=================\n",color=:yellow)
        end
     end

 for l=1:L
    #constructing the pool-output variables as Y_lj
        for j=1:Int(Pooloutputsize[l])
            h=Pooloutput[l];
            Y_lj[l,h[j]]=yout[Int(sum(Pooloutputsize[1:l-1]))+j];
            Y_lj_index[l,h[j]]=sort(indexin([yout[Int(sum(Pooloutputsize[1:l-1]))+j]],y))[1];
        end
        for j=1:J
            if !(j in Pooloutput[l])
                Y_lj[l,j]=0;
            end
        end
        if l==L
                 printstyled("=================Pool-output variables are build!=================\n",color=:yellow)
        end

 # constructing the input-pool
        for i=1:Int(Poolinputsize[l])
            h=Poolinput[l];
            Y_il[h[i],l]=yinputpool[Int(sum(Poolinputsize[1:l-1]))+i];
            Y_il_index[h[i],l]=sort(indexin([yinputpool[Int(sum(Poolinputsize[1:l-1]))+i]],y))[1];
        end
        for i=1:I
            if !(i in Poolinput[l])
                Y_il[i,l]=0;
            end
        end
        if l==L
                 printstyled("=================Input-Pool variables are build!=================\n",color=:yellow)
        end
        P_lk[l,:]=ypool[(l-1)*K+1:l*K];
        P_lk_index[l,:]=sort(indexin(ypool[(l-1)*K+1:l*K],y));
        if l==L
                 printstyled("=================Pool_specification variables are build!=================\n",color=:yellow)
        end
 end
 printstyled("==========The model is being constructed ... =============\n \n",color=:yellow)

 # display(P_lk)
 # display(Y_il)
 # display(Y_lj)
 # display(Y_ij)

 f=MO*(sum(costI .*Y_il)+
   sum(costL .*Y_lj)+
   sum(costJ .*Y_ij)-
   sum((sum(Y_ij,dims=1)+sum(Y_lj,dims=1))*Demandcost) #the first sum is to convert one-lenght array to a polynomial
   );
 printstyled("==========Objective function is constructed! =============\n \n",color=:yellow)

 SizeCIindex=size(CIindex,1);
 SizeCJindex=size(CJindex,1);
 SizeCLindex=size(CLindex,1);
 SizelCIindex=size(lCIindex,1);
 SizelCJindex=size(lCJindex,1);
 SizelCLindex=size(lCLindex,1);

 eq=Polyopt.Poly{Int64}[]; #equality constraints
 for l=1:L
    # balance constraint of input and output of a pool l
    eq=[eq;
     (1/I)*(sum(Y_il[:,l])-sum(Y_lj[l,:]));
    (1/J)*(-sum(Y_il[:,l])+sum(Y_lj[l,:]))]; #just for replacing with 2 inequalities
    # balance between input and outputs specifications of a pool l

    for k=1:K
        lhelp=Lambda[:,k];
        # println("lhelp'*Y_il[:,l]",lhelp'*Y_il[:,l])
        eq=[eq;
         (1/(Ma_lamb*I))*(lhelp'*Y_il[:,l]-((Ma_lamb-Mi_lamb)*P_lk[l,k]+Mi_lamb)*sum(Y_lj[l,:]));
        (1/(Ma_lamb*J))*(-lhelp'*Y_il[:,l]+((Ma_lamb-Mi_lamb)*P_lk[l,k]+Mi_lamb)*sum(Y_lj[l,:]))];
    end
 end

 g=Polyopt.Poly{Int64}[];#inequality constraints
 #capacity ristrictions
 for i=1:SizeCIindex

    g=[g;
        1-(MO/C_I[CIindex[i]])*(sum(Y_il[CIindex[i],:])+
                                sum(Y_ij[CIindex[i],:]))];
 end
 for i=1:SizelCIindex
       g=[g;
            - (lCI[lCIindex[i]])/(MO*(L+J))+(1/(L+J))*(sum(Y_il[lCIindex[i],:])+
                                                     sum(Y_ij[lCIindex[i],:])
                                                    )];
 end
 for l=1:SizeCLindex
       g=[g;
           1 - MO/(C_L[CLindex[l]])*(sum(Y_lj[CLindex[l],:])) ];
 end

 for l=1:SizelCLindex
       g=[g;
           - (lCL[lCLindex[l]])/(MO*(J))+(1/J)*(sum(Y_lj[lCLindex[l],:])) ];
 end

 for j=1:SizeCJindex
       g=[g;
           1 - MO/(C_J[CJindex[j]])*(sum(Y_ij[:,CJindex[j]])+
                                   sum(Y_lj[:,CJindex[j]]))];
 end
 for j=1:SizelCJindex
       g=[g;
           - (lCJ[lCJindex[j]])/(MO*(I+L))+(1/(I+L))*(sum(Y_ij[:,CJindex[j]])+
                                                     sum(Y_lj[:,CJindex[j]]))];
 end

#MU_max constraint
 for j=1:J
    for k=1:K
      # if Mi_lamb>=0
        g=[g;
            (1/(max(Mu_max[j,k],0)*(I+L)-sum(min.(Lambda[:,k],0))-min(0,Mi_lamb)*L))*(Mu_max[j,k]*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))-
                                        (sum(Lambda[:,k].*Y_ij[:,j])+
                                          (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
                                                    Mi_lamb*sum(Y_lj[:,j])
                                        ))];
      # elseif Mu_max[j,k]>=0
      #           g=[g;
      #       (1/((Mu_max[j,k])*(I+L)-sum(min(Lambda[:,k],0))-Mi_lamb*J))*(Mu_max[j,k]*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))-
      #                                   (sum(Lambda[:,k].*Y_ij[:,j])+
      #                                     (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
      #                                               Mi_lamb*sum(Y_lj[:,j])
      #                                   ))];
      # else
      #     g=[g;
      #       (1/(-sum(min(Lambda[:,k],0))-Mi_lamb*J))*(Mu_max[j,k]*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))-
      #                                   (sum(Lambda[:,k].*Y_ij[:,j])+
      #                                     (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
      #                                               Mi_lamb*sum(Y_lj[:,j])
      #                                   ))];
      # end
    end

 end
 #Mu_min constraint
  for j=1:J
    for k=1:K
        help=1/Mu_min[j,k];
        if help!=0
            # if min( Mu_min[j,k],Mi_lamb)>=0
             g=[g;
               (1/(sum(max(Lambda[:,k],0))+L*(Ma_lamb-Mi_lamb)+L*max(0,Mi_lamb)-min(Mu_min[j,k],0)*(I+L)))*(sum(Lambda[:,k].*Y_ij[:,j])+
                            (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
                             Mi_lamb*sum(Y_lj[:,j])-
                             (Mu_min[j,k])*(sum(Y_ij[:,j])+
                                      sum(Y_lj[:,j])
                                      )
                             )];
            # elseif
            #     g=[g;
            #    (1/((Ma_lamb*I+L-Mu_min[j,k]*(I+L))))*(sum(Lambda[:,k].*Y_ij[:,j])+
            #                         (Ma_lamb-Mi_lamb)*(sum(P_lk[:,k].*Y_lj[:,j]))+
            #                         Mi_lamb*sum(Y_lj[:,j])-
            #                         (Mu_min[j,k])*(sum(Y_ij[:,j])+sum(Y_lj[:,j]))
            #                         )];
            # end
        end
    end
  end

 #sign constraints
 for l=1:L
    for k=1:K
        g=[g;
            P_lk[l,k]];
    end
 end
 for i=1:I
    for l=1:L
        if(AI[i,l]==1)
        g=[g;
            Y_il[i,l]];
        end
    end
 end

 for j=1:J
    for l=1:L
        if(AL[l,j]==1)
        g=[g;
            Y_lj[l,j]];
        end
    end
 end
 for j=1:J
    for i=1:I
        if(AJ[i,j]==1)
        g=[g;
            Y_ij[i,j]];
        end
    end
 end
 #pipe capacity restriction
  g=[g;
       (ones(SizeUIindex[1],1) - ((MO./UI[UIindex]).*Y_il[UIindex]))[:,1]];

   g=[g;
       (ones(SizeULindex[1],1) - ((MO./UL[ULindex]).*Y_lj[ULindex]))[:,1]];

  g=[g;
       (ones(SizeUJindex[1],1)-((MO./UJ[UJindex]).*Y_ij[UJindex]))[:,1]];

  f,g,eq,n
end

function pooling_with_eq_BSOS(data_a , d::Int, k::Int)

	@time f_eq, g_eq, eq_eq, n_eq=with_equality(data_a);
	 f_eqscale = maximum(abs.(f_eq.c));
    # f_eqscale=1;
    f_eq=Polyopt.truncate(1/f_eqscale*f_eq);
    g_eq=Array{Polyopt.Poly{Float64},1}(g_eq);
    for i=1:length(g_eq)
        g_eq[i] = Polyopt.truncate(0.9*g_eq[i])
    end
    eq_eq=Array{Polyopt.Poly{Float64},1}(eq_eq);
    for i=1:length(eq_eq)
        eq_eq[i] = Polyopt.truncate(0.9*eq_eq[i])
    end
    	I = Array{Int,1}[ collect(1:n_eq) ];
      y=variables("y",n_eq);
   for j=1:size(I,1)
            z=variables("z",size(I[j],1));
            Ihelp=I[j];
            for i=1:size(Ihelp,1)
                z[i]=y[Ihelp[i]]^2;
            end
            g_eq=[g_eq;
                 1-sum(z)*(1/size(Ihelp,1))];
   end
   println("$(size([g_eq;eq_eq]))")
      prob = bsosprob_chordal(d, k, I, f_eq, [g_eq;eq_eq]);
      printstyled("====================MOSEK is solving the problem... =============\n",color=:yellow)
      time_solution=@elapsed X, t, l, y ,solsta = solve_mosek(prob, tolrelgap=1e-10);
      return( t*f_eqscale,solsta, time_solution)
end

function pooling_without_eq_BSOS(data_a , d::Int, k::Int)

	@time f, g, n=elimination_equality(data_a);
    println("f before scaling:",f)
    fscale = maximum(abs.(f.c));
     # fscale=1;
    f=Polyopt.truncate(1/fscale*f);

    g=Array{Polyopt.Poly{Float64},1}(g);
    for i=1:length(g)
        g[i] = Polyopt.truncate(0.9*g[i])
    end
    I = Array{Int,1}[ collect(1:n) ];
    y=variables("y",n);

      for j=1:size(I,1)
            z=variables("z",size(I[j],1));
            Ihelp=I[j];
            for i=1:size(Ihelp,1)
                z[i]=y[Ihelp[i]]^2;
            end
            g=[g;
                 1-sum(z)*(1/size(Ihelp,1))];
   end
   println("$(size(g))")
   println("$(fscale)")

      prob = bsosprob_chordal(d, k, I, f, g);
      #row uniqueness
      # q=size(Al,1);
      # Ahelp=prob.El[1];
      # A=prob.Al[1]*Ahelp';
      # for i=1:q
      #    Ahelp=prob.El[i];
      #     A=[A;
      #     prob.Al[i]*Ahelp'];
      # end
      # prob.Al=unique_ro(A,Al);
      # printstyled("====================MOSEK is solving the problem... =============\n",color=:yellow)
      # io = fileio +".txt";
      time_solution=@elapsed X, t, l,  y ,solsta= solve_mosek(prob, tolrelgap=1e-10);
      deg=Array{Int64}[];
      for j=1:size(g,1)
        deg=[deg;
          g[j].deg];
      end
      tau=maximum(deg);
      tau=maximum([f.deg, 2*k , d*tau]);
      return( t*fscale,solsta,time_solution)
end

function pooling_with_eq_Sparse_BSOS(data_a , d::Int, k::Int)

	@time f_eq, g_eq, eq_eq, n_eq=with_equality(data_a);
	 f_eqscale = maximum(abs.(f_eq.c));
    # f_eqscale=1;
    f_eq=Polyopt.truncate(1/f_eqscale*f_eq);
    g_eq=Array{Polyopt.Poly{Float64},1}(g_eq);
    for i=1:length(g_eq)
        g_eq[i] = Polyopt.truncate(0.9*g_eq[i])
    end
    eq_eq=Array{Polyopt.Poly{Float64},1}(eq_eq);
    for i=1:length(eq_eq)
        eq_eq[i] = Polyopt.truncate(0.9*eq_eq[i])
    end
    y=variables("y",n_eq);
    	 I_block = Polyopt.chordal_embedding(Polyopt.correlative_sparsity(f_eq,[g_eq;eq_eq]));

   for j=1:size(I_block,1)
            z=variables("z",size(I_block[j],1));
            Ihelp=I_block[j];
            for i=1:size(Ihelp,1)
                z[i]=y[Ihelp[i]]^2;
            end
            g_eq=[g_eq;
                 1-sum(z)*(1/size(Ihelp,1))];
   end
   println("$(size(g_eq))")
   println("$(size(I_block))")
      prob = bsosprob_chordal(d, k, I_block, f_eq, g_eq,eq_eq);
       printstyled("====================MOSEK is solving the problem... =============\n",color=:yellow)
     time_solution=@elapsed X, t, l,  y ,solsta = solve_mosek(prob, tolrelgap=1e-10);
      return( size(I_block),t*f_eqscale,solsta,time_solution)
   # return(I_block)
end

function pooling_without_eq_Sparse_BSOS(data_a , d::Int, k::Int)

	f, g, n=elimination_equality(data_a);
	 fscale = maximum(abs.(f.c));
    # fscale=1;
    f=Polyopt.truncate(1/fscale*f);
    g=Array{Polyopt.Poly{Float64},1}(g);
    for i=1:length(g)
        g[i] = Polyopt.truncate(0.9*g[i])
    end

    	I = Polyopt.chordal_embedding(Polyopt.correlative_sparsity(f,g));
      y=variables("y",n);
   for j=1:size(I,1)
            z=variables("z",size(I[j],1));
            Ihelp=I[j];
            for i=1:size(Ihelp,1)
                z[i]=y[Ihelp[i]]^2;
            end
            g=[g;
                 1-sum(z)*(1/size(Ihelp,1))];
   end
   println("$(size(g))")
      prob = bsosprob_chordal(d, k, I, f, g);
      printstyled("====================MOSEK is solving the problem... =============\n",color=:yellow)
      time_solution=@elapsed X, t, l,  y ,solsta = solve_mosek(prob, tolrelgap=1e-10);
      return( size(I),t*fscale,solsta, time_solution)
end

function pooling_with_eq_Merge_Sparse_BSOS(data_a , d::Int, k::Int)

  @time f_eq, g_eq, eq_eq, n_eq=with_equality(data_a);
   f_eqscale = maximum(abs.(f_eq.c));
    # f_eqscale=1;
    f_eq=Polyopt.truncate(1/f_eqscale*f_eq);
    g_eq=Array{Polyopt.Poly{Float64},1}(g_eq);
    for i=1:length(g_eq)
        g_eq[i] = Polyopt.truncate(0.9*g_eq[i])
    end
    eq_eq=Array{Polyopt.Poly{Float64},1}(eq_eq);
    for i=1:length(eq_eq)
        eq_eq[i] = Polyopt.truncate(0.9*eq_eq[i])
    end
      I_block = Polyopt.chordal_embedding(Polyopt.correlative_sparsity(f_eq,[g_eq;eq_eq]));
      printstyled("====================Checking the overlaps between the cliques... =============\n",color=:yellow)
      #        I=Array{Array{Int64,1},1}(size(I_block)[1]);
      # l=1;
      # while l<=size(I,1)
      #   I[l]=[];
      #   q=1;
      #   if size(I_block,1)>=1
      #     while q<=size(I_block,1)
      #      if size(intersect(I_block[q],Y_il[:,l]),1)>0
      #          I[l]=union(I_block[q],I[l]);
      #          deleteat!(I_block,q);
      #       else
      #          q=q+1;
      #      end
      #     end
      #     l=l+1;
      #   else
      #       deleteat!(I,l);
      #   end
      # end
      # l=1;
      # while l<=size(I,1)
      #   if size(I[l],1)==0
      #     deleteat!(I,l);
      #   else
      #     l=l+1;
      #   end
      # end
      # if size(I_block,1)>=1
      #   I_block=[I;I_block];
      # else
      #   I_block=I;
      # end
     I=I_block;
     Cmatr=collect(combinations(collect(1:size(I,1)),2));
     i=1;
     while i<=size(Cmatr,1)
          if (size(intersect(I[Cmatr[i][1]],I[Cmatr[i][2]]),1)>0.75*min(size(I[Cmatr[i][1]],1),size(I[Cmatr[i][2]],1)))
             I[Cmatr[i][1]]=sort(union(I[Cmatr[i][1]],I[Cmatr[i][2]]));
             deleteat!(I,Cmatr[i][2]);
             Cmatr=collect(combinations(collect(1:size(I,1)),2));
             i=0;
             printstyled("==================== 2 cliques are megred ...=============\n",color=:yellow)
          end
          i=i+1;
     end
     I_block=I;
     # I_block=Array{Int,1}[ collect(1:n_eq) ];
     y=variables("y",n_eq);
   for j=1:size(I_block,1)
            z=variables("z",size(I_block[j],1));
            Ihelp=I_block[j];
            for i=1:size(Ihelp,1)
                z[i]=y[Ihelp[i]]^2;
            end
            g_eq=[g_eq;
                 1-sum(z)*(1/size(Ihelp,1))];
   end
   println("$(size(g_eq))")
   println("$(size(I_block,1))")
   # return(I_block)
      prob = bsosprob_chordal(d, k, I_block, f_eq, g_eq,eq_eq);
      # printstyled("====================MOSEK is solving the problem... =============\n",color=:yellow)
     time_solution=@elapsed X, t, l,  y ,solsta = solve_mosek(prob, tolrelgap=1e-8);
      return( size(I_block),t*f_eqscale,solsta,time_solution)
end

function pooling_without_eq_Merge_Sparse_BSOS(data_a , d::Int, k::Int)

  f, g, n=elimination_equality(data_a);
   fscale = maximum(abs.(f.c));
    # fscale=1;
    f=Polyopt.truncate(1/fscale*f);
    g=Array{Polyopt.Poly{Float64},1}(g);
    for i=1:length(g)
        g[i] = Polyopt.truncate(0.9*g[i])
    end

      I = Polyopt.chordal_embedding(Polyopt.correlative_sparsity(f,g));
      Cmatr=collect(combinations(collect(1:size(I,1)),2));
     i=1;
     while i<=size(Cmatr,1)
          if (size(intersect(I[Cmatr[i][1]],I[Cmatr[i][2]]),1)>0.75*min(size(I[Cmatr[i][1]],1),size(I[Cmatr[i][2]],1)))
             I[Cmatr[i][1]]=sort(union(I[Cmatr[i][1]],I[Cmatr[i][2]]));
             deleteat!(I,Cmatr[i][2]);
             Cmatr=collect(combinations(collect(1:size(I,1)),2));
             i=0;
             printstyled("==================== 2 cliques are megred ...=============\n",color=:yellow)
          end
          i=i+1;
     end
     # I=Array{Int,1}[ collect(1:n) ];
     y=variables("y",n);
   for j=1:size(I,1)
            z=variables("z",size(I[j],1));
            Ihelp=I[j];
            for i=1:size(Ihelp,1)
                z[i]=y[Ihelp[i]]^2;
            end
            g=[g;
                 1-sum(z)*(1/size(Ihelp,1))];
   end
   println("$(size(g))")
      prob = bsosprob_chordal(d, k, I, f, g);
      # printstyled("====================MOSEK is solving the problem... =============\n",color=:yellow)
     time_solution=@elapsed X, t, l,  y ,solsta = solve_mosek(prob, tolrelgap=1e-8);
      return( size(I),t*fscale,solsta,time_solution)
end

ext: jl

max_stars_repo_name: shudianzhao/PoolOpt

max_issues_repo_name: shudianzhao/PoolOpt

max_forks_repo_name: shudianzhao/PoolOpt


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct ESST1A <: AVR
        UEL_flags::Int
        PSS_flags::Int
        Tr::Float64
        Vi_lim::Tuple{Float64, Float64}
        Tc::Float64
        Tb::Float64
        Tc1::Float64
        Tb1::Float64
        Ka::Float64
        Ta::Float64
        Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kc::Float64
        Kf::Float64
        Tf::Float64
        K_lr::Float64
        I_lr::Float64
        V_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

This excitation system supplies power through a transformer from the generator terminals and its regulated by a controlled rectifier (via thyristors).
Parameters of IEEE Std 421.5 Type ST1A Excitacion System. ESST1A in PSSE and PSLF

# Arguments
- `UEL_flags::Int`: Code input for Underexcitization limiter (UEL) entry. Not supported., validation range: `(1, 3)`, action if invalid: `warn`
- `PSS_flags::Int`: Code input for Power System Stabilizer (PSS) or (VOS) entry., validation range: `(1, 2)`
- `Tr::Float64`: Regulator input filter time constant in s, validation range: `(0, 0.1)`, action if invalid: `warn`
- `Vi_lim::Tuple{Float64, Float64}`: Voltage error limits (regulator input) (Vi_min, Vi_max)
- `Tc::Float64`: First regulator denominator (lead) time constant in s, validation range: `(0, 10)`, action if invalid: `warn`
- `Tb::Float64`: First regulator denominator (lag) time constant in s, validation range: `(0, 20)`
- `Tc1::Float64`: Second regulator denominator (lead) time constant in s, validation range: `(0, 10)`, action if invalid: `warn`
- `Tb1::Float64`: Second regulator denominator (lead) time constant in s, validation range: `(0, 20)`, action if invalid: `warn`
- `Ka::Float64`: Voltage regulator gain, validation range: `(50, 1000)`, action if invalid: `warn`
- `Ta::Float64`: Voltage regulator time constant in s, validation range: `(0, 0.5)`, action if invalid: `warn`
- `Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for regulator output `(Va_min, Va_max)`
- `Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for exciter output `(Vr_min, Vr_max)`
- `Kc::Float64`: Rectifier loading factor proportional to commutating reactance, validation range: `(0, 0.3)`, action if invalid: `warn`
- `Kf::Float64`: Rate feedback gain, validation range: `(0, 0.3)`, action if invalid: `warn`
- `Tf::Float64`: Rate feedback time constant in s, validation range: `(eps(), 1.5)`, action if invalid: `error`
- `K_lr::Float64`: Exciter output current limiter gain, validation range: `(0, 5)`, action if invalid: `warn`
- `I_lr::Float64`: Exciter output current limit reference, validation range: `(0, 5)`, action if invalid: `warn`
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vm: Sensed terminal voltage,
	Vr1: First Lead-lag state,
	Vr2: Second lead-lag state,
	Va: Regulator output state,
	Vr3: Feedback output state
- `n_states::Int`: ST1A has 5 states
- `states_types::Vector{StateTypes}`: ST1A has 5 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ESST1A <: AVR
    "Code input for Underexcitization limiter (UEL) entry. Not supported."
    UEL_flags::Int
    "Code input for Power System Stabilizer (PSS) or (VOS) entry."
    PSS_flags::Int
    "Regulator input filter time constant in s"
    Tr::Float64
    "Voltage error limits (regulator input) (Vi_min, Vi_max)"
    Vi_lim::Tuple{Float64, Float64}
    "First regulator denominator (lead) time constant in s"
    Tc::Float64
    "First regulator denominator (lag) time constant in s"
    Tb::Float64
    "Second regulator denominator (lead) time constant in s"
    Tc1::Float64
    "Second regulator denominator (lead) time constant in s"
    Tb1::Float64
    "Voltage regulator gain"
    Ka::Float64
    "Voltage regulator time constant in s"
    Ta::Float64
    "Limits for regulator output `(Va_min, Va_max)`"
    Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Limits for exciter output `(Vr_min, Vr_max)`"
    Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Rectifier loading factor proportional to commutating reactance"
    Kc::Float64
    "Rate feedback gain"
    Kf::Float64
    "Rate feedback time constant in s"
    Tf::Float64
    "Exciter output current limiter gain"
    K_lr::Float64
    "Exciter output current limit reference"
    I_lr::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    ext::Dict{String, Any}
    "The states are:
	Vm: Sensed terminal voltage,
	Vr1: First Lead-lag state,
	Vr2: Second lead-lag state,
	Va: Regulator output state,
	Vr3: Feedback output state"
    states::Vector{Symbol}
    "ST1A has 5 states"
    n_states::Int
    "ST1A has 5 states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ESST1A(UEL_flags, PSS_flags, Tr, Vi_lim, Tc, Tb, Tc1, Tb1, Ka, Ta, Va_lim, Vr_lim, Kc, Kf, Tf, K_lr, I_lr, V_ref=1.0, ext=Dict{String, Any}(), )
    ESST1A(UEL_flags, PSS_flags, Tr, Vi_lim, Tc, Tb, Tc1, Tb1, Ka, Ta, Va_lim, Vr_lim, Kc, Kf, Tf, K_lr, I_lr, V_ref, ext, [:Vm, :Vr1, :Vr2, :Va, :Vr3], 5, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function ESST1A(; UEL_flags, PSS_flags, Tr, Vi_lim, Tc, Tb, Tc1, Tb1, Ka, Ta, Va_lim, Vr_lim, Kc, Kf, Tf, K_lr, I_lr, V_ref=1.0, ext=Dict{String, Any}(), states=[:Vm, :Vr1, :Vr2, :Va, :Vr3], n_states=5, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    ESST1A(UEL_flags, PSS_flags, Tr, Vi_lim, Tc, Tb, Tc1, Tb1, Ka, Ta, Va_lim, Vr_lim, Kc, Kf, Tf, K_lr, I_lr, V_ref, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function ESST1A(::Nothing)
    ESST1A(;
        UEL_flags=1,
        PSS_flags=1,
        Tr=0,
        Vi_lim=(0.0, 0.0),
        Tc=0,
        Tb=0,
        Tc1=0,
        Tb1=0,
        Ka=0,
        Ta=0,
        Va_lim=(min=0.0, max=0.0),
        Vr_lim=(min=0.0, max=0.0),
        Kc=0,
        Kf=0,
        Tf=0,
        K_lr=0,
        I_lr=0,
        V_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`ESST1A`](@ref) `UEL_flags`."""
get_UEL_flags(value::ESST1A) = value.UEL_flags
"""Get [`ESST1A`](@ref) `PSS_flags`."""
get_PSS_flags(value::ESST1A) = value.PSS_flags
"""Get [`ESST1A`](@ref) `Tr`."""
get_Tr(value::ESST1A) = value.Tr
"""Get [`ESST1A`](@ref) `Vi_lim`."""
get_Vi_lim(value::ESST1A) = value.Vi_lim
"""Get [`ESST1A`](@ref) `Tc`."""
get_Tc(value::ESST1A) = value.Tc
"""Get [`ESST1A`](@ref) `Tb`."""
get_Tb(value::ESST1A) = value.Tb
"""Get [`ESST1A`](@ref) `Tc1`."""
get_Tc1(value::ESST1A) = value.Tc1
"""Get [`ESST1A`](@ref) `Tb1`."""
get_Tb1(value::ESST1A) = value.Tb1
"""Get [`ESST1A`](@ref) `Ka`."""
get_Ka(value::ESST1A) = value.Ka
"""Get [`ESST1A`](@ref) `Ta`."""
get_Ta(value::ESST1A) = value.Ta
"""Get [`ESST1A`](@ref) `Va_lim`."""
get_Va_lim(value::ESST1A) = value.Va_lim
"""Get [`ESST1A`](@ref) `Vr_lim`."""
get_Vr_lim(value::ESST1A) = value.Vr_lim
"""Get [`ESST1A`](@ref) `Kc`."""
get_Kc(value::ESST1A) = value.Kc
"""Get [`ESST1A`](@ref) `Kf`."""
get_Kf(value::ESST1A) = value.Kf
"""Get [`ESST1A`](@ref) `Tf`."""
get_Tf(value::ESST1A) = value.Tf
"""Get [`ESST1A`](@ref) `K_lr`."""
get_K_lr(value::ESST1A) = value.K_lr
"""Get [`ESST1A`](@ref) `I_lr`."""
get_I_lr(value::ESST1A) = value.I_lr
"""Get [`ESST1A`](@ref) `V_ref`."""
get_V_ref(value::ESST1A) = value.V_ref
"""Get [`ESST1A`](@ref) `ext`."""
get_ext(value::ESST1A) = value.ext
"""Get [`ESST1A`](@ref) `states`."""
get_states(value::ESST1A) = value.states
"""Get [`ESST1A`](@ref) `n_states`."""
get_n_states(value::ESST1A) = value.n_states
"""Get [`ESST1A`](@ref) `states_types`."""
get_states_types(value::ESST1A) = value.states_types
"""Get [`ESST1A`](@ref) `internal`."""
get_internal(value::ESST1A) = value.internal

"""Set [`ESST1A`](@ref) `UEL_flags`."""
set_UEL_flags!(value::ESST1A, val) = value.UEL_flags = val
"""Set [`ESST1A`](@ref) `PSS_flags`."""
set_PSS_flags!(value::ESST1A, val) = value.PSS_flags = val
"""Set [`ESST1A`](@ref) `Tr`."""
set_Tr!(value::ESST1A, val) = value.Tr = val
"""Set [`ESST1A`](@ref) `Vi_lim`."""
set_Vi_lim!(value::ESST1A, val) = value.Vi_lim = val
"""Set [`ESST1A`](@ref) `Tc`."""
set_Tc!(value::ESST1A, val) = value.Tc = val
"""Set [`ESST1A`](@ref) `Tb`."""
set_Tb!(value::ESST1A, val) = value.Tb = val
"""Set [`ESST1A`](@ref) `Tc1`."""
set_Tc1!(value::ESST1A, val) = value.Tc1 = val
"""Set [`ESST1A`](@ref) `Tb1`."""
set_Tb1!(value::ESST1A, val) = value.Tb1 = val
"""Set [`ESST1A`](@ref) `Ka`."""
set_Ka!(value::ESST1A, val) = value.Ka = val
"""Set [`ESST1A`](@ref) `Ta`."""
set_Ta!(value::ESST1A, val) = value.Ta = val
"""Set [`ESST1A`](@ref) `Va_lim`."""
set_Va_lim!(value::ESST1A, val) = value.Va_lim = val
"""Set [`ESST1A`](@ref) `Vr_lim`."""
set_Vr_lim!(value::ESST1A, val) = value.Vr_lim = val
"""Set [`ESST1A`](@ref) `Kc`."""
set_Kc!(value::ESST1A, val) = value.Kc = val
"""Set [`ESST1A`](@ref) `Kf`."""
set_Kf!(value::ESST1A, val) = value.Kf = val
"""Set [`ESST1A`](@ref) `Tf`."""
set_Tf!(value::ESST1A, val) = value.Tf = val
"""Set [`ESST1A`](@ref) `K_lr`."""
set_K_lr!(value::ESST1A, val) = value.K_lr = val
"""Set [`ESST1A`](@ref) `I_lr`."""
set_I_lr!(value::ESST1A, val) = value.I_lr = val
"""Set [`ESST1A`](@ref) `V_ref`."""
set_V_ref!(value::ESST1A, val) = value.V_ref = val
"""Set [`ESST1A`](@ref) `ext`."""
set_ext!(value::ESST1A, val) = value.ext = val
"""Set [`ESST1A`](@ref) `states_types`."""
set_states_types!(value::ESST1A, val) = value.states_types = val


text:
using HelloMain

@info("Starting HelloJuliaHub...")

HelloMain.serve()

@info("Finished HelloJuliaHub")


text:
facts("\n   ---testing KSP solvers---") do
# create vectors and matrices

  b = PetscVec(comm);
  PetscVecSetType(b, "mpi");
  PetscVecSetSizes(b,sys_size, PetscInt(comm_size*sys_size));

  low, high = PetscVecGetOwnershipRange(b)
  b_global_indices = Array(low:PetscInt(high - 1))

  x = PetscVec(comm);
  PetscVecSetType(x,"mpi");
  PetscVecSetSizes(x,sys_size, PetscInt(comm_size*sys_size));

  low, high = PetscVecGetOwnershipRange(x)
  x_global_indices = Array(low:PetscInt(high - 1))

  for i=1:sys_size
    idxm = [ b_global_indices[i] ]   # index
    val = [ rhs[i] ]  # value
    PetscVecSetValues(b, idxm, val, PETSC_INSERT_VALUES)
  end

  PetscVecAssemblyBegin(b)
  PetscVecAssemblyEnd(b)

  A = PetscMat(comm)
  PetscMatSetType(A, "mpiaij")

  PetscMatSetSizes(A,sys_size,sys_size,PetscInt(comm_size*sys_size),PetscInt(comm_size*sys_size));
  PetscSetUp(A);

  low, high = PetscMatGetOwnershipRange(A)
  mat_global_indices = Array(low:PetscInt(high - 1))

  for i=1:sys_size
    for j = 1:sys_size
      idxm = [ mat_global_indices[i] ]  # row index
      idxn = [ mat_global_indices[j] ]  # column index
      PetscMatSetValues(A,idxm, idxn, [A_julia[i,j]],PETSC_INSERT_VALUES);
    end
  end

  PetscMatAssemblyBegin(A,PETSC_MAT_FINAL_ASSEMBLY);
  PetscMatAssemblyEnd(A,PETSC_MAT_FINAL_ASSEMBLY);

# perform solve
ksp = KSP(comm)
KSPSetOperators(ksp, A, A)
KSPSetFromOptions(ksp)
KSPSetUp(ksp)
KSPSolve(ksp, b, x)
reason = KSPGetConvergedReason(ksp)

println("KSP convergence reason = ", reason)
@fact reason => greater_than(0)  # convergence

PetscView(ksp)

# copy solution back to Julia
x_copy = zeros(PetscScalar, sys_size)
#idx = Array(0:2)
#idx = zeros(PetscInt, sys_size)
#for i=1:sys_size
#  idx[i] = i-1
#end

PetscVecGetValues(x, sys_size, x_global_indices, x_copy)
println("x_copy = ", x_copy)
println("x_julia = ", x_julia)
for i=1:sys_size
    @fact x_copy[i] => roughly(x_julia[i], atol=1e-14)
end

#PetscView(x, 0)

PetscDestroy(ksp)

println("   \n--- Testing LGMRES --- ")
# test using non default KSP method

# perform solve
ksp = KSP(comm)
KSPSetOperators(ksp, A, A)

# test PC
println("   \n---Testing PC Options ---")
pc = KSPGetPC(ksp)

PCSetType(pc, PETSc.PCBJACOBI)
pctype = PCGetType(pc)
println("pctype = ", pctype)
@fact pctype => PETSc.PCBJACOBI
#PCFactorSetUseInPlace(pc, PetscBool(true))

### Do some KSP setup

KSPSetFromOptions(ksp)
KSPSetTolerances(ksp, rtol, abstol, dtol, maxits)
KSPSetInitialGuessNonzero(ksp, PetscBool(true))
KSPSetType(ksp, PETSc.KSPLGMRES)
KSPSetUp(ksp)

### More PC testing
println("getting sub KSP objects")
n_local, first_local, ksp_arr = PCBJacobiGetSubKSP(pc)
ksp2 = ksp_arr[1]

println("ksp2 = ", ksp2)
println("typeof(ksp2) = ", typeof(ksp2))
println("getting sub KSP PC")
pc2 = KSPGetPC(ksp2)
println("finsihed getting sub KSP PC")

println("setting ReusePreconditioner")
PCSetReusePreconditioner(pc2, PetscBool(true))
@fact PCGetReusePreconditioner(pc2) => true
println("finished setting ReusePreconditioner")

PCFactorSetAllowDiagonalFill(pc2, PetscBool(true))
@fact PCFactorGetAllowDiagonalFill(pc2) => true

PCFactorSetLevels(pc2, PetscInt(1))
@fact PCFactorGetLevels(pc2) => 1  # should be pc2

PCSetReusePreconditioner(pc2, PetscBool(true))
@fact PCGetReusePreconditioner(pc2) => true

PCFactorSetFill(pc, PetscReal(7.0))

#=
PCJacobiSetType(pc, PETSc.PC_JACOBI_ROWMAX)
@fact PCJacobiGetType(pc) => PETSc.PC_JACOBI_ROWMAX
=#

#=
tmp = PCFactorGetUseInPlace(pc)
println("tmp = ", tmp)
@fact PCFactorGetUseInPlace(pc) => true
=#

KSPSolve(ksp, b, x)
reason = KSPGetConvergedReason(ksp)
ksptype = KSPGetType(ksp)
println("finished calling KSPGetType")
println("typeof(ksptype) = ", typeof(ksptype))
println("ksptype = ", ksptype)

@fact ksptype => PETSc.KSPLGMRES

rtol_ret, abstol_ret, dtol_ret, maxits_ret = KSPGetTolerances(ksp)

@fact rtol_ret => roughly(rtol)
@fact abstol_ret => roughly(abstol)
@fact dtol_ret => roughly(dtol)
@fact maxits_ret => maxits
@fact KSPGetInitialGuessNonzero(ksp) => true
println("KSP convergence reason = ", reason)
@fact reason => greater_than(0)  # convergence

rnorm = KSPGetResidualNorm(ksp)
@fact rnorm => less_than(abstol)
PetscView(ksp)

# copy solution back to Julia
x_copy = zeros(PetscScalar, sys_size)
#idx = Array(0:2)
#idx = zeros(PetscInt, sys_size)
#for i=1:sys_size
#  idx[i] = i-1
#end

PetscVecGetValues(x, sys_size, x_global_indices, x_copy)
println("x_copy = ", x_copy)
println("x_julia = ", x_julia)
for i=1:sys_size
    @fact x_copy[i] => roughly(x_julia[i], atol=1e-14)
end

PetscDestroy(x)
PetscDestroy(b)
PetscDestroy(A)
PetscDestroy(ksp)
end


text:
@json struct Object
    type::String
    sha::String
    url::String
end

@json struct Tag
    ref::String
    node_id::String
    url::String
    object::Object
end

endpoint(::GitHubAPI, ::typeof(get_tags), owner::AStr, repo::AStr) =
    Endpoint(:GET, "/repos/$owner/$repo/git/refs/tags")
@not_implemented(::GitHubAPI, ::typeof(get_tags), ::Int64)
into(::GitHubAPI, ::typeof(get_tags)) = Vector{Tag}


text:
# doc.jl - infrastructure

"""
    process(::Array{AppliSales.Order, 1}; path="./test_invoices.sqlite")

    process(::Array{UnpaidInvoice, 1}, ::Array{AppliGeneralLedger.BankStatement, 1}; path="./test_invoices.sqlite")

- Creates UnpaidInvoice's from AppliSale.Order's, archive them, and creates AppliGeneralLedger.Entry's for the general ledger.
- Creates PaidInvoices's from UnpaidInvoices by using AppliGeneralLedger.BankStatement's, and creates AppliGeneralLedger.Entry's for the general ledger.

# Example

```
julia> using AppliSales

julia> using AppliGeneralLedger

julia> using AppliAR

julia> const PATH_CSV = "./bank.csv"

julia> orders = AppliSales.process()

julia> journal_entries_1 = AppliAR.process(orders)

julia> stms = AppliAR.read_bank_statements(PATH_CSV)

julia> unpaid_invoices = AppliAR.retrieve_unpaid_invoices()

julia> journal_entries_2 = AppliAR.process(unpaid_invoices, stms)

julia> cmd = `rm test_invoicing.txt test_invoicing_paid.txt invoicenbr.txt`

julia> run(cmd)
```
"""
function process end

"""
    read_bank_statements(path::String)

Retrieves bank statements from a CSV-file.

# Example

```
julia> const PATH_CSV = "./bank.csv"

julia> stms = AppliAR.read_bank_statements(PATH_CSV)
```
"""
function read_bank_statements end

"""
    retrieve_unpaid_invoices(;path="./test_invoicing.txt")::Array{UnpaidInvoice, 1}

Retrieves UnpaidInvoice's from a text file.

# Example

```
julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process()

julia> AppliAR.process(orders)

julia> unpaid_invoices = retrieve_unpaid_invoices()

julia> cmd = `rm test_invoicing.txt invoicenbr.txt`

julia> run(cmd)
```
"""
function retrieve_unpaid_invoices end

"""
    retrieve_paid_invoices(;path="./test_invoicing_paid.txt")::Array{PaidInvoice, 1}

Retrieves PaidInvoice's from a text file.

# Example

```
julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process()

julia> AppliAR.process(orders)

julia> unpaid_invoices = retrieve_unpaid_invoices()

julia> const PATH_CSV = "./bank.csv"

julia> stms = AppliAR.read_bank_statements(PATH_CSV)

julia> AppliAR.process(unpaid_invoices, stms)

julia> paid_invoices = AppliAR.retrieve_paid_invoices()

julia> cmd = `rm test_invoicing.txt test_invoicing_paid.txt invoicenbr.txt`

julia> run(cmd)
```
"""
function retrieve_paid_invoices end


text:
# Runs the forward model of all ITMIX glaciers, updates cached files and can plot all.
@time using BITEModel
const BM=BITEModel
verbose = false

itmixversion = 2
use_glogem = false
update_cache = true
plotting = false

# make them available in case the loop craps out
global glacier, gid, gl,gb,pp,pm,pn,pl,hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio

# these don't load yet
toskip = [] # [:Brewster, :Devon, :Synthetic2, :Synthetic3, :Unteraar, :Urumqi, :Washmawapta #=only v1=#]
for (ii,glacier) in enumerate(BM.ITMIXnames)
    #ii<9 && continue

    glacier in toskip && continue
    itmixversion==1 && glacier in BM.itmix2only && continue
    println(" ")
    println(glacier)
    gid = ITMIXGlacier(glacier, itmixversion)
    print("Loading ")
    # gid = ITMIXGlacier(glacier, 2); pl = BM.LoadPara(gid); dt, pl = BM.make_datatable(gid, pl);  ds = dt.outline;
    # outline, iscalving, isicecap, outline_proj =  BM.load!(ds, gid)
    @time gl,gb,pp,pm,pn,pl = BM.init_forward(gid, verbose;
                                              update_cache=update_cache, use_glogem=use_glogem);

    print("Running fwd-model ")
    @time hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb =
        BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp);
    @show prod(size(gl.dem))
    if plotting
        println("Plotting")
        display(BM.plot_run(gb, hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, pp, pm,
                            reuse=false, plot2d=[:h,:iv,:ivlog][1]))
    end
end


text:
"""
    struct FileManager <: Any
Hold a pointer to a `clang::FileManager` object.
"""
struct FileManager
    ptr::CXFileManager
end

Base.unsafe_convert(::Type{CXCodeGenOptions}, x::FileManager) = x.ptr
Base.cconvert(::Type{CXCodeGenOptions}, x::FileManager) = x


text:
module TestNullableVector
    using NullableArrays
    using Base.Test

    A = [1:10...]
    B = [1:5...]

    X = NullableArray(A)
    Y = NullableArray(B)

    #--- Base.push!

    # Base.push!{T, V}(X::NullableVector{T}, v::V)
    Z = NullableVector{Int}()
    push!(Z, 5)
    @test isequal(Z[1], Nullable(5))

    # Base.push!{T, V}(X::NullableVector{T}, v::Nullable{V})
    push!(Z, Nullable())
    @test isequal(Z, NullableArray([5, 1], [false, true]))
    push!(Z, Nullable(5))
    @test isequal(Z, NullableArray([5, 1, 5], [false, true, false]))

    ZZ = NullableVector{Int}()
    push!(ZZ, 5, Nullable(), Nullable(5))
    @test isequal(ZZ, Z)

    #--- Base.pop!

    # Base.pop!(X::NullableArray)
    @test isequal(pop!(Z), Nullable(5))
    @test isequal(pop!(Z), Nullable{Int}())
    @test isequal(pop!(Z), Nullable(5))

    #--- Base.unshift!

    # Base.unshift!(X::NullableVector, v)
    @test isequal(unshift!(X, 3), NullableArray(vcat(3, [1:10...])))

    # Base.unshift!(X::NullableVector, v::Nullable)
    @test isequal(unshift!(X, Nullable(2)),
                  NullableArray(vcat([2, 3], [1:10...])))
    @test isequal(unshift!(X, Nullable{Int}()),
                  NullableArray(vcat([1, 2, 3], [1:10...]),
                                vcat(true, fill(false, 12))
                  )
          )

    @test isequal(unshift!(Y, 1, Nullable(), Nullable(3)),
                  NullableArray([1, 2, 3, 1, 2, 3, 4, 5],
                                [false, true, false, false,
                                false, false, false, false]
                  )
          )

    #--- Base.shift!

    # Base.shift!{T}(X::NullableVector{T})
    Z = NullableArray([1:10...])

    @test isequal(shift!(Z), Nullable(1))
    @test isequal(Z, NullableArray([2:10...]))

    unshift!(Z, Nullable{Int}())

    @test isequal(shift!(Z), Nullable{Int}())

    #--- test Base.splice!

    # Base.splice!(X::NullableVector, i::Integer, ins=_default_splice)
    A = [1:10...]
    B = [1:10...]
    X = NullableArray(B)

    i = rand(1:10)
    @test isequal(splice!(X, i), Nullable(splice!(A, i)))
    @test isequal(X, NullableArray(A))

    i = rand(1:9)
    j = rand(1:9)
    @test isequal(splice!(X, i, j), Nullable(splice!(A, i, j)))
    @test isequal(X, NullableArray(A))

    A = [1:10...]
    B = [1:10...]
    X = NullableArray(B)
    i = rand(1:5)
    n = rand(3:5)
    @test isequal(splice!(X, i, [1:n...]), Nullable(splice!(A, i, [1:n...])))
    @test isequal(X, NullableArray(A))

    # Base.splice!{T<:Integer}(X::NullableVector,
    #                          rng::UnitRange{T},
    #                          ins=_default_splice)

    # test with length(rng) > length(ins)
    A = [1:20...]
    B = [1:20...]
    X = NullableArray(B)
    f = rand(1:7)
    d = rand(3:5)
    l = f + d
    ins = [1, 2]
    @test isequal(splice!(X, f:l, ins),
                  NullableArray(splice!(A, f:l, ins)))
    @test isequal(X, NullableArray(A))

    i = rand(1:length(X))
    @test isequal(splice!(X, 1:i), NullableArray(splice!(A, 1:i)))

    A = [1:20...]
    B = [1:20...]
    X = NullableArray(B)
    f = rand(10:15)
    d = rand(3:5)
    l = f + d
    ins = [1, 2]
    @test isequal(splice!(X, f:l, ins),
                  NullableArray(splice!(A, f:l, ins)))
    @test isequal(X, NullableArray(A))

    # test with length(rng) < length(ins)
    A = [1:20...]
    B = [1:20...]
    X = NullableArray(B)
    f = rand(1:7)
    d = rand(3:5)
    l = f + d
    ins = [1, 2, 3, 4, 5, 6, 7]
    @test isequal(splice!(X, f:l, ins),
                  NullableArray(splice!(A, f:l, ins)))
    @test isequal(X, NullableArray(A))

    A = [1:20...]
    B = [1:20...]
    X = NullableArray(B)
    f = rand(10:15)
    d = rand(3:5)
    l = f + d
    ins = [1, 2, 3, 4, 5, 6, 7]
    @test isequal(splice!(X, f:l, ins),
                  NullableArray(splice!(A, f:l, ins)))
    @test isequal(X, NullableArray(A))

    #--- test Base.deleteat!

    # Base.deleteat!(X::NullableArray, inds)
    X = NullableArray(1:10)
    @test isequal(deleteat!(X, 1), NullableArray(2:10))

    #--- test Base.append!

    # Base.append!(X::NullableVector, items::AbstractVector)
    @test isequal(append!(X, [11, 12]),
                  NullableArray(2:12))
    @test isequal(append!(X, [Nullable(13), Nullable(14)]),
                  NullableArray(2:14))
    @test isequal(append!(X, [Nullable(15), Nullable{Int}()]),
                  NullableArray(Nullable{Int}[2:15..., Nullable()]))

    #--- test Base.prepend!

    # Base.prepend!(X::NullableVector, items::AbstractVector)

    X = NullableArray(3:12)
    @test isequal(prepend!(X, [1, 2]),
                  NullableArray(1:12))
    @test isequal(prepend!(X, [Nullable(-1), Nullable(0)]),
                  NullableArray(-1:12))
    @test isequal(prepend!(X, [Nullable{Int}(), Nullable(-2)]),
                  NullableArray(Nullable{Int}[Nullable(), -2:12...]))

    #--- test Base.sizehint!

    # Base.sizehint!(X::NullableVector, newsz::Integer)
    sizehint!(X, 20)

    #--- test padnull!

    # padnull!{T}(X::NullableVector{T}, front::Integer, back::Integer)
    X = NullableArray(1:5)
    padnull!(X, 2, 3)
    @test length(X.values) == 10
    @test X.isnull == vcat(true, true, fill(false, 5), true, true, true)

    # padnull(X::NullableVector, front::Integer, back::Integer)
    X = NullableArray(1:5)
    Y = padnull(X, 2, 3)
    @test length(Y.values) == 10
    @test Y.isnull == vcat(true, true, fill(false, 5), true, true, true)

    #--- test Base.reverse!/Base.reverse

    y = NullableArray([nothing, 2, 3, 4, nothing, 6], Int, Void)
    @assert isequal(reverse(y),
                    NullableArray([6, nothing, 4, 3, 2, nothing], Int, Void))

    # check case where only nothing occurs in final position
    @assert isequal(unique(NullableArray([1, 2, 1, nothing], Int, Void)),
                    NullableArray([1, 2, nothing], Int, Void))

    # Base.reverse!(X::NullableVector, s=1, n=length(X))
    # check for case where isbits(eltype(X)) = false
    Z = NullableArray(Array{Int, 1}[[1, 2], [3, 4], [5, 6]])
    @test isequal(reverse!(Z),
                  NullableArray(Array{Int, 1}[[5, 6], [3, 4], [1, 2]]))

    # Base.reverse!(X::NullableVector, s=1, n=length(X))
    # check for case where isbits(eltype(X)) = false & anynull(X) = true
    A = fill([1,2], 20)
    Z = NullableArray(Array{Int, 1}, 20)
    i = rand(2:7)
    for i in [i-1, i, i+1, 20 - (i + 2), 20 - (i - 1)]
        Z[i] = [1, 2]
    end
    vals = Z.values
    nulls = Z.isnull
    @test isequal(reverse!(Z), NullableArray(A, reverse!(Z.isnull)))

    ## empty!(X::NullableVector)
    n = rand(1:1_000)
    A, M = rand(n), rand(Bool, n)
    X = NullableArray(A, M)
    empty!(X)
    @test isempty(X)
end


text:
function sv_ms_nested_assign()
    s1 = MatrixShape(Float64,5,10)
    s2 = VectorShape(Float64,5)
    s3 = ScalarShape(Float64)
    ms1 = MultiShape(s1 = s1, s2 = s2, s3=s3)

    s4 = VectorShape(Float64,3)
    s5 = ScalarShape(Float64)
    ms2 = MultiShape(s4=s4, s5=s5)

    s6 = ScalarShape(Float64)
    ms3 = MultiShape(s6=s6, ms2=ms2)

    ms = MultiShape(ms1=ms1, ms2=ms2, ms3=ms3)

    src = rand!(allocate(ms))
    shaped_src = ms(rand!(allocate(ms)))
    dst = zeros(ms)
    shaped_dst = ms(zeros(ms))

    bench = @benchmarkable begin
        sv_ms_nested_assign_bench($shaped_dst, $shaped_src)
    end teardown=(@assert $shaped_dst == $shaped_src)

    std = @benchmarkable begin
        sv_ms_nested_assign_std($dst, $src)
    end teardown = (@assert $dst == $src)

    bench, std
end

function sv_ms_nested_assign_bench(dst, src)
    @uviews dst src @inbounds begin
        dst.ms1.s1 .= src.ms1.s1
        dst.ms1.s2 .= src.ms1.s2
        dst.ms1.s3 = src.ms1.s3

        dst.ms2.s4 .= src.ms2.s4
        dst.ms2.s5 = src.ms2.s5

        dst.ms3.s6 = src.ms3.s6
        dst.ms3.ms2.s4 .= src.ms3.ms2.s4
        dst.ms3.ms2.s5 = src.ms3.ms2.s5
    end
end

function sv_ms_nested_assign_std(dst, src)
    @uviews dst src @inbounds begin
        dst[1:50] .= view(src, 1:50)
        dst[51:55] .= view(src, 51:55)
        dst[56] = src[56]
        dst[57:59] .= view(src, 57:59)
        dst[60] = src[60]
        dst[61] = src[61]
        dst[62:64] .= view(src, 62:64)
        dst[65] = src[65]
    end
end

function sv_ms_assign()
    s1 = MatrixShape(Float64, 5,10)
    s2 = VectorShape(Float64, 15)
    s3 = ScalarShape(Float64)
    ms = MultiShape(s1 = s1, s2 = s2, s3=s3)

    src = rand!(allocate(ms))
    shaped_src = ms(rand!(allocate(ms)))
    dst = zeros(ms)
    shaped_dst = ms(zeros(ms))

    bench = @benchmarkable begin
        sv_ms_assign_bench($shaped_dst, $shaped_src)
    end teardown=(@assert $shaped_src == $shaped_dst)

    std = @benchmarkable begin
        sv_ms_assign_std($dst, $src)
    end teardown = @assert $dst == $src

    bench, std
end

function sv_ms_assign_bench(src, dst)
    @uviews dst src @inbounds begin
        dst.s1 .= src.s1
        dst.s2 .= src.s2
        dst.s3 = src.s3
    end
end

function sv_ms_assign_std(dst, src)
    @uviews dst src @inbounds begin
        dst[1:50] .= view(src, 1:50)
        dst[51:65] .= view(src, 51:65)
        dst[66] = src[66]
    end
end

function sv_vs_assign()
    s = VectorShape(Float64, 15)

    src = rand!(allocate(s))
    shaped_src = ShapedView(rand!(allocate(s)), s)
    dst = zeros(s)
    shaped_dst = ShapedView(zeros(s), s)

    bench = @benchmarkable begin
        @inbounds $shaped_dst .= $shaped_src
    end teardown=(@assert $shaped_src == $shaped_dst)
    std = @benchmarkable begin
        @inbounds $dst .= $src
    end
    bench, std
end


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    public static final Condition c = new Condition();
    int{L /c H} x = 0;

    public static void main(String[] args) {
	H.c.x = 42;
	H.y = H.c.x;
    }
}

abstract class H {
    public static final C c = new C();
    public static int{L /c.c H} y;
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
function main()
	N = 4_000_000
	s = 0
	i1 = 1
	i2 = 1
	while true
		fib = i1 + i2
		i1 = i2
		i2 = fib
		if fib > N
			break
		end
		if fib % 2 == 0
			s += fib
		end
	end
	println(s)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
# Calculate sample means for each bin and dataset
function meansofbins(N, M, v)
    # M cannot be larger thn the size of v
    cuts = [N:N:M];         # Cut points for each bucket
    means = Float64[];      # Array holding mean values
    start = 1;
    for cut in cuts
        push!(means,mean(v[start:cut]));
        start = cut + 1;
    end
    return means;
end


text:
export build_xxz_mpo,
       build_ssh_mpo

"""
    build_xxz_mpo(::Type{T},
                  L::Int,
                  Δ::Float64,
                  has_pbc::Bool=false) where {T<:Number}

Build the XXZ Hamiltonian:
    H = ∑_i (s⁺_i s⁻_{i+1} + h.c.) + Δ ∑_i sᶻ_i sᶻ_{i+1}
"""
function build_xxz_mpo(::Type{T},
                       L::Int,
                       Δ::Float64,
                       has_pbc::Bool=false) where {T<:Number}
    # S+ S- terms.
    J = diagm(1 => ones(L-1))
    J[1, L] = has_pbc ? 1. : 0.
    J = Symmetric(J)
    # Sz Sz terms.
    V = diagm(1 => fill(Δ, L-1))
    V[1, L] = has_pbc ? Δ : 0.

    H = Mpo(T, L, 2)
    add_ops!(H, "s+", "s-", convert.(T, J))
    add_ops!(H, "sz", "sz", convert.(T, V))

    return H
end

build_xxz_mpo(L, Δ, has_pbc=false) = build_xxz_mpo(Float64, L, Δ, has_pbc)

"""
    build_ssh_mpo(::Type{T},
                  L::Int,
                  t::Float64,
                  ε::Float64,
                  w::Vector{Float64},
                  has_pbc::Bool=false) where {T<:Number}

Build the SSH Hamiltonian:
    H = t/2 ∑_i [(1 + ε*(-)^i)*c⁺_i c_{i+1} + h.c.] + ∑_i w_i n_i
"""
function build_ssh_mpo(::Type{T},
                       L::Int,
                       t::Float64,
                       ε::Float64,
                       w::Vector{Float64},
                       has_pbc::Bool=false) where {T<:Number}
    J = zeros(T, L, L)
    for i=1:L
        J[i, i] = w[i]
        i == L && continue
        J[i, i+1] = t/2*(1 + (iseven(i) ? ε : -ε))
    end
    J[L, 1] = has_pbc ? t/2*(1 + (iseven(L) ? ε : -ε)) : 0.
    J = Symmetric(J)

    H = Mpo(T, L, 2)
    add_ops!(H, "c+", "c", J, ferm_op="Z")

    return H
end

build_ssh_mpo(L, t, ε, w, has_pbc=false) = build_ssh_mpo(Float64, L, t, ε, w, has_pbc)


text:
function solve(lines, allow_diagonals)
    map = zeros(Int, (1001, 1001))
    for line in lines
        from, to = split(line, " -> ")
        toCoords(x) = (a -> parse(Int, a)).(split(x, ","))
        from = toCoords(from) .+ 1  # Shift all by one, so that Julia indices would work
        to = toCoords(to) .+ 1
        if any(from .== to) || allow_diagonals
            while from != to
                map[from[2], from[1]] += 1
                from += (to .> from) - (to .< from)
            end
            map[from[2], from[1]] += 1
        end
    end
    return sum(map .>= 2)
end

function main()
    lines = readlines()
    println(solve(lines, false))
    println(solve(lines, true))
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
mutable struct Boltzmann{T<:Integer} <: AbstractAgent
  id::T
  wealth::T
end

function boltzmann_model(; numagents)
  model = ABM(Boltzmann{Int64}, scheduler=random_activation)
  for i in 1:numagents
    model.agents[i] = Boltzmann(i, 1)
  end
  return model
end

function Boltzmann_step!(agent::AbstractAgent, model::ABM)
  if agent.wealth == 0
    return
  else
    random_agent = model.agents[rand(keys(model.agents))]
    agent.wealth -= 1
    random_agent.wealth += 1
  end
end

@testset "Model tests" begin

  numagents = 20
  n = 10
  model = boltzmann_model(numagents=numagents)
  agent_properties = [:wealth]
  data = step!(model, Boltzmann_step!, n, agent_properties, when=1:n)
  @test size(data) == (numagents*(n+1), length(agent_properties)+2)

  nreps = 3
  data = step!(model, Boltzmann_step!, n, agent_properties, when=1:n, replicates=nreps)
  @test size(data) == ((numagents*(n+1))*nreps, 4)
end


text:
refVals = []
refPrecs = []

#! format: off
# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/jmc/cliMa/cliMa_update/test/Ocean/SplitExplicit/simple_box_ivd.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [  "oce Q_3D",   "u[1]", -1.66327525125104291881e-01,  1.63223287218068308091e-01,  4.57054199141500999692e-03,  2.21420008866952226778e-02 ],
 [  "oce Q_3D",   "u[2]", -2.16296814964942768489e-01,  2.44105976363460708267e-01, -2.57418981108877616831e-03,  2.49256704872522875938e-02 ],
 [  "oce Q_3D",       :η, -8.06721661019456748321e-01,  3.47844179114868090608e-01, -3.03380420083551127861e-04,  2.92158779168222249023e-01 ],
 [  "oce Q_3D",       :θ,  3.58991077492186536763e-03,  9.91901429055192807027e+00,  2.49592172215209329167e+00,  2.17949780843674956188e+00 ],
 [   "oce aux",       :w, -1.66963788550371410252e-04,  1.64989216196191281925e-04,  5.43118005323253843601e-07,  1.44011197008218777164e-05 ],
 [   "oce aux",    :pkin, -8.84700789178211977060e+00,  0.00000000000000000000e+00, -3.26178716377690980366e+00,  2.50117471654691314598e+00 ],
 [   "oce aux",     :wz0, -2.12922215345145814233e-05,  3.08186566412866553883e-05, -1.53993475195290930982e-10,  7.82327443240865498639e-06 ],
 [   "oce aux", "u_d[1]", -1.57099271619415420398e-01,  1.22268940698781983234e-01, -7.84250724128123411372e-05,  1.15452721709468873745e-02 ],
 [   "oce aux", "u_d[2]", -2.12964763257990491452e-01,  2.30578860790304179806e-01,  2.76337873040526178152e-05,  1.75119371638010334902e-02 ],
 [   "oce aux", "ΔGu[1]", -2.33686060871737704252e-06,  3.34281015324516805469e-06, -4.07662146058456868485e-08,  3.37377005914367152159e-07 ],
 [   "oce aux", "ΔGu[2]", -2.46095864522509354306e-06,  2.23139964389116549296e-06,  1.29016375934756163675e-06,  7.46140899703298511983e-07 ],
 [   "oce aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573163901915703900e+06 ],
 [ "baro Q_2D",   "U[1]", -1.77642579751613922667e+01,  6.90024701711891168543e+01,  4.64758402961474637038e+00,  1.81125164737756669808e+01 ],
 [ "baro Q_2D",   "U[2]", -3.48668457221064898022e+01,  8.56043954389619301537e+01, -2.60107324488959568143e+00,  1.80557567005594421516e+01 ],
 [ "baro Q_2D",       :η, -8.06729084084832237522e-01,  3.47845108442749906263e-01, -3.03371430044830942430e-04,  2.92173862421732766226e-01 ],
 [  "baro aux",  "Gᵁ[1]", -3.34281015324516816989e-03,  2.33686060871737691716e-03,  4.07662146058456834074e-05,  3.37393707332951826861e-04 ],
 [  "baro aux",  "Gᵁ[2]", -2.23139964389116544213e-03,  2.46095864522509347530e-03, -1.29016375934756159609e-03,  7.46177836457347174598e-04 ],
 [  "baro aux", "U_c[1]", -1.77641349600227265171e+01,  6.90014815427443437557e+01,  4.64756464081321674087e+00,  1.81123414932683921563e+01 ],
 [  "baro aux", "U_c[2]", -3.48658691554727866446e+01,  8.56044062434442594167e+01, -2.60101037971774573521e+00,  1.80557220201615180599e+01 ],
 [  "baro aux",     :η_c, -8.06721661019456748321e-01,  3.47844179114868090608e-01, -3.03380420083561807253e-04,  2.92173242116136766544e-01 ],
 [  "baro aux", "U_s[1]", -1.77642579751613922667e+01,  6.90024701711891168543e+01,  4.64758402961474637038e+00,  1.81125164737756669808e+01 ],
 [  "baro aux", "U_s[2]", -3.48668457221064898022e+01,  8.56043954389619301537e+01, -2.60107324488959568143e+00,  1.80557567005594421516e+01 ],
 [  "baro aux",     :η_s, -8.06729084084832237522e-01,  3.47845108442749906263e-01, -3.03371430044830942430e-04,  2.92173862421732766226e-01 ],
 [  "baro aux",  "Δu[1]", -4.51957131810988945505e-04,  3.99130692959488437947e-04, -1.09363921539353876238e-05,  1.02327771113594037156e-04 ],
 [  "baro aux",  "Δu[2]", -2.42568998424025800446e-04,  3.82991529986384212410e-04,  5.18845385039475763891e-05,  7.47269801857268440920e-05 ],
 [  "baro aux",  :η_diag, -8.06666755113148337131e-01,  3.47953196151883525911e-01, -3.02575448959868274160e-04,  2.92168272060050748795e-01 ],
 [  "baro aux",      :Δη, -5.84450730272745300198e-04,  7.45283859099332701703e-04, -8.04971123704421912570e-07,  1.64194305918131369573e-04 ],
 [  "baro aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15578885204060329124e+06 ],
]
parr = [
 [  "oce Q_3D",   "u[1]",    12,    12,    12,    12 ],
 [  "oce Q_3D",   "u[2]",    12,    12,    12,    12 ],
 [  "oce Q_3D",       :η,    12,    12,     8,    12 ],
 [  "oce Q_3D",       :θ,    12,    12,    12,    12 ],
 [   "oce aux",       :w,    12,    12,     8,    12 ],
 [   "oce aux",    :pkin,    12,    12,    12,    12 ],
 [   "oce aux",     :wz0,    12,    12,     8,    12 ],
 [   "oce aux", "u_d[1]",    12,    12,    12,    12 ],
 [   "oce aux", "u_d[2]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[1]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[2]",    12,    12,    12,    12 ],
 [   "oce aux",       :y,    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[1]",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[2]",    12,    12,    12,    12 ],
 [ "baro Q_2D",       :η,    12,    12,     8,    12 ],
 [  "baro aux",  "Gᵁ[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Gᵁ[2]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_c,    12,    12,     8,    12 ],
 [  "baro aux", "U_s[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_s[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_s,    12,    12,     8,    12 ],
 [  "baro aux",  "Δu[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Δu[2]",    12,    12,    12,    12 ],
 [  "baro aux",  :η_diag,    12,    12,     8,    12 ],
 [  "baro aux",      :Δη,     9,     9,     6,    10 ],
 [  "baro aux",       :y,    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on

max_forks_repo_name: jm-c/CLIMA


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

export show, ideal

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::NfOrdGenIdlSet)
  print(io, "Set of ideals of $(order(a))\n")
end

function show(io::IO, a::NfOrdGenIdl)
  print(io, "Ideal of (")
  print(io, order(a), ")\n")
  print(io, "with basis matrix\n")
  print(io, basis_mat(a))
end

################################################################################
#
#  Construction
#
################################################################################

doc"""
***
    ideal(O::NfOrd, a::fmpz) -> NfOrdIdl

> Returns the ideal of $\mathcal O$ which is generated by $a$.
"""
ideal(O::NfOrdGen, a::fmpz)  = NfOrdGenIdl(O, deepcopy(a))

doc"""
***
    ideal(O::NfOrd, a::Int) -> NfOrdIdl

> Returns the ideal of $\mathcal O$ which is generated by $a$.
"""
ideal(O::NfOrdGen, a::Int) = NfOrdGenIdl(O, a)

doc"""
***
    ideal(O::NfOrd, A::fmpz_mat) -> NfOrdIdl

> Assuming that $A$ is a $d \times d$ nonsingular integer matrix,
> this function returns the ideal $I$ given by $I = \bigoplus \mathbf Z \alpha_i$,
> where $(\alpha_1,\dotsc,\alpha_d) = (\omega_1,\dotsc,\omega_d) A^t$.
"""
ideal(O::NfOrdGen, a::fmpz_mat) = NfOrdGenIdl(O, deepcopy(a))

doc"""
***
    ideal(O::NfOrd, a::NfOrdElem) -> NfOrdIdl

> Returns the ideal of $\mathcal O$ which is generated by $a$.
"""
ideal(O::NfOrdGen, a::NfOrdElem{NfOrdGen}) = NfOrdGenIdl(O, deepcopy(a))


text:
# ---
# title: 382. Linked List Random Node
# id: problem382
# author: zhwang
# date: 2022-02-15
# difficulty: Medium
# categories: Reservoir Sampling
# link: <https://leetcode.com/problems/linked-list-random-node/description/>
# hidden: true
# ---
#
# Given a singly linked list, return a random node's value from the linked list.
# Each node must have the **same probability** of being chosen.
#
# **Follow up:**
# What if the linked list is extremely large and its length is unknown to you?
# Could you solve this efficiently without using extra space?
#
# **Example:**
#
#
#
#     // Init a singly linked list [1,2,3].
#     ListNode head = new ListNode(1);
#     head.next = new ListNode(2);
#     head.next.next = new ListNode(3);
#     Solution solution = new Solution(head);
#
#     // getRandom() should return either 1, 2, or 3 randomly. Each element should have equal probability of returning.
#     solution.getRandom();
#
#
#
## @lc code=start
using LeetCode

function get_random(head::ListNode)::Int
    node, i, res = head, 1, 0
    while !isnothing(node)
        if rand(1:i) == 1  ## 1/i prob for value change
            res = node.val
        end
        node = node.next
        i += 1
    end
    res
end
## @lc code=end

## @lc test=start
@testset "382.linked-list-random-node.jl" begin
    node = ListNode{Int}([1, 2, 2, 3])
    @test begin
        freq = Dict{Int, Int}(1 => 0, 2 => 0, 3 => 0)
        for _ in 1:10000
            freq[get_random(node)] += 1
        end
        num = freq[1] + freq[2] + freq[3]
        !(freq[1] == 2 * freq[2] == freq[3]) &&
        abs(freq[1] / num - freq[3] / num) + abs(freq[1] / num - freq[2] / num / 2) <= 0.03
    end
end
## @lc test=end

## @lc info=start
# link: [solution to 382](https://leetcode-cn.com/problems/linked-list-random-node/solution/382-lian-biao-sui-ji-jie-dian-by-rex-rs-0g47/)
## @lc info=end


text:
module Bibi
    const Tera = BigInt(10)^12
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ d55bc09f-90f8-4a68-8f66-49a0538c682d
begin
	using PlutoUI
	using Test
	bigbreak = HTML("
"^4)
	html""""""
end

# ╔═╡ 67910dd6-000d-4d10-9473-be3a5944cf5d
md"
# Crash-course: _Understanding_ Julia

## The origin story of Julia. 🥚 🐔

The _'two language problem'_

> Scientific computing has traditionally required the highest performance, yet domain   experts have largely moved to slower dynamic languages for daily work.

Our answer
> The Julia programming language is a **flexible dynamic language**, appropriate for scientific and numerical computing, **with performance comparable to traditional statically-typed** languages.
"

# ╔═╡ 4ef295dd-11e0-4fe2-9c11-9a0e8a27f2ab
html""""""

# ╔═╡ b54a2e64-0303-4ff6-82f4-fb3f91f9b771
html"""
In this course we will learn why Julia is fast and the basics of programming in Julia.
 """

# ╔═╡ c5a7ba98-0f68-414d-8923-2c35caac34cc
bigbreak

# ╔═╡ 9138271c-ddfe-4181-9506-ba3922593163
md"""
### We will use material from:
0. [JuliaAcademy: Introduction to Julia](https://github.com/JuliaAcademy/Introduction-to-Julia)
1. [The manual from the official Julia documentation](https://docs.julialang.org/en/v1/manual/getting-started/),
2. [The Fast Track to Julia (cheat-sheet)](https://steffenpl.github.io/Julia-Cheat-Sheet/).
4. [MATLAB–Python–Julia (cheat-sheet)](https://cheatsheets.quantecon.org/)
5. [Makie.jl Plotting function examples](https://makie.juliaplots.org/stable/examples/plotting_functions/)

!!! correct "Hint: Where to find details?"
	If you need more details on a specific topic, read the **the official documentation**! It is _actually_ **well-written** and an **excellent introduction**!

"""

# ╔═╡ e1ffc859-b8f7-4f74-aa64-c83041e5c808
bigbreak

# ╔═╡ 9a7536a5-5c0f-4514-a64a-02504b1bad46
md"""
## Outline

**Today:**
- Intro
- The Julia RELP and Juno IDE
- Essential data types
- Control flow
- Arrays and functions

**Tuesday:**
- Scope of variables
- The type system and multiple dispatch
- Custom datatypes
- The package manager
- Plotting _(with Makie.jl)_

**Wednesday:**
- Linear Algebra
- More plotting
- Interfaces
- Parametric methods and types

**Thursday:**
- Modules
- Testing
- Workflow
- ... (let's see were we end up)

"""

# ╔═╡ cf47e80d-a42d-4df0-81b0-9360b9decae9
bigbreak

# ╔═╡ 6592ddea-7126-458d-8a38-f058ab3cb366
html"""
Disclaimer: Julia is not  simply "MATLAB but faster"!
"""

# ╔═╡ bfd077aa-8c54-44b0-bafb-67cb1f8e7176
md"Even if similar on the surface, those languages are very different."

# ╔═╡ 3545e4f9-df86-4261-943b-e79a48274716
md"""
Imagine this MATLAB code (taken from the official documentation):
```julia
# MATLAB code
x = -2:0.25:2;
y = x;
[X,Y] = meshgrid(x, y);

F = X.*exp(-X.^2-Y.^2);
```
"""

# ╔═╡ 273be1ff-58b0-4dab-9726-80b1fa8aa0c1
md"Typical naive translation to Julia *could* look like this:"

# ╔═╡ 1af43c9b-151f-4b22-aa88-110ad018e0ad
meshgrid(x,y) = x * one.(y'), one.(x) * y'   # Julia does not have meshgrid function, let's just make one for this example

# ╔═╡ ae2fc32c-a471-4e64-9d33-e9e6bf677166
let
	x = -2:0.25:2
	y = x
	X, Y = meshgrid(x, y)

	F = X .* exp.( -X.^2 - Y.^2 )

	F[1:3,1:3]
end  # hide output

# ╔═╡ 3608f7cf-86e3-4cdb-ada6-46c6beb3310a
md"...but there is a small problem here:"

# ╔═╡ b6406021-00e6-419f-8e15-9f80196433eb
md"""
!!! danger "MATLAB ≠ Julia"
	Good MATLAB code is often poor Julia code. *(And vice versa!)*

However, good Julia code is often still pretty (and usually needs much less memory)!
"""

# ╔═╡ ff3539e9-296f-403d-a031-86c4df8566f7
[i+j for i in 1:10, j in 1:2]

# ╔═╡ be19363a-dfd3-430e-b1c9-f8b69309411d
let
	xs = -2:0.25:2
	ys = xs

	F = [x * exp( -x^2-y^2 ) for x in xs, y in ys]

	F[1:3,1:3]
end

# ╔═╡ 0c50f149-4cea-4793-acd2-91b58595dc3b
md"""
### 👉 Do not blindly apply MATLAB or python tips to Julia!

Check out the [Noteworthy Differences from other Languages](https://docs.julialang.org/en/v1/manual/noteworthy-differences/) if you want to see more! *(We will see some of these later.)*

Most prominent difference: There is often no reason to avoid `for` loops in Julia!
"""

# ╔═╡ a725c129-1f57-47d3-810c-bfee4f600bcf
md"""#### ☝ Important: In Julia, the first index is "1" *and not zero*.
*(For all non-MATLAB users: I feel you, but that's how it is.)*"""

# ╔═╡ 80986d60-d993-4293-bd22-72cdef15988b
bigbreak

# ╔═╡ 9bccf94a-19d1-4d31-82d5-d882b1fd1e29
md"## Task #0: How to work with Julia in the terminal (aka REPL) and the Juno IDE."

# ╔═╡ add8900c-f7a6-450e-8e04-0e7237cf80f3
@bind n PlutoUI.Slider(1:100)

# ╔═╡ 250726a2-4c56-46fe-a2b1-8f5da0dc761e
n^2

# ╔═╡ 56a560dc-dfac-43a9-afa2-a3cd6733d14f
n^5

# ╔═╡ 711dc218-0d78-4a39-a380-bc461e249ba5
bigbreak

# ╔═╡ 7bcf336b-c2dc-40ad-bb43-82fbca1c8565
md"""## Task #1: Compute 1+2 with Julia

...let's call it the mathematician's Hello World 😜
"""

# ╔═╡ e42073de-0d42-42c3-9f1a-53f2ff54b91c
bigbreak

# ╔═╡ 52fb0f4d-c1c0-4ca2-9bd1-913414fcf9c8
md"""
## Task #2: Some basic math
!!! correct "Warm-up: Fibonacci numbers"
    Compute the first 100.000 Fibonacci numbers (using a `for` loop).
"""

# ╔═╡ 20b697f4-0bf2-426b-ae3b-a3864386cc19
bigbreak

# ╔═╡ 6714a575-dc6a-4926-bcfa-349b91aca16e
md"# This looks just like MATLAB! Why is that faster?"

# ╔═╡ 5c66ea8a-ca3e-4d07-b16c-8f4015f87671
html"
The answer is not as simple as it might seem...
"

# ╔═╡ 5efe07fc-b7e4-4069-83c6-b1534c39ac99
md"""
Let's consider a simple function, like this generic `add(x,y) = x + y`.

The first time we call `add` with new input types,
the Julia **just-in-time** (JIT) compiler generates a code with *static types* and compiles that!
"""

# ╔═╡ 172cf452-791b-4d7a-934a-0dfa05374e14
add(x,y) = x+y

# ╔═╡ 8d9ca351-c8f6-49bd-8a91-2c34e883d0c4
md"For `(Int64, Int64)`"

# ╔═╡ 3675a743-1ddc-4cd4-b954-8ca354bea723
with_terminal() do
	@code_llvm debuginfo=:none add(1,2)
end

# ╔═╡ d533573e-10fd-4480-9b17-93ef3dacb91e
md"For `(Float64, Float64)`"

# ╔═╡ 2e34f342-f0a8-4e62-a4bc-6bda81c5b0c6
with_terminal() do
	@code_llvm debuginfo=:none add(1.0,2.0)
end

# ╔═╡ 6c9df6b6-746a-4721-b77b-52718abab2e5
md"Doing the same in a for loop is still fast!"

# ╔═╡ 549a0344-a442-4dfd-b01f-f1b06723df7e
function add_4_numbers(x,y)
	z = zero(eltype(x))
	@inbounds for i in 1:4
		z += x[i] + y[i]
	end
	return z
end

# ╔═╡ ab6b4c8b-0b62-4c50-bf2f-f9b969d23aff
with_terminal() do
	@code_llvm add_4_numbers([1,2,3,4],[1,2,3,5])
end

# ╔═╡ dd6dff51-8765-48da-9ccc-6fb6ca228bff
bigbreak

# ╔═╡ 299f6a19-a805-4ff2-8108-95dd94a5bb6d

# ╔═╡ 92da736b-2da2-43bb-9c8c-c661a327acc0
md"""
### Meeting Julia's JIT compiler and type stability

"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[compat]
PlutoUI = "~0.7.30"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "92f91ba9e5941fc781fecf5494ac1da87bdac775"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "5c0eb9099596090bb3215260ceca687b888a1575"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.30"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─d55bc09f-90f8-4a68-8f66-49a0538c682d
# ╟─67910dd6-000d-4d10-9473-be3a5944cf5d
# ╟─4ef295dd-11e0-4fe2-9c11-9a0e8a27f2ab
# ╟─b54a2e64-0303-4ff6-82f4-fb3f91f9b771
# ╟─c5a7ba98-0f68-414d-8923-2c35caac34cc
# ╟─9138271c-ddfe-4181-9506-ba3922593163
# ╟─e1ffc859-b8f7-4f74-aa64-c83041e5c808
# ╟─9a7536a5-5c0f-4514-a64a-02504b1bad46
# ╟─cf47e80d-a42d-4df0-81b0-9360b9decae9
# ╟─6592ddea-7126-458d-8a38-f058ab3cb366
# ╟─bfd077aa-8c54-44b0-bafb-67cb1f8e7176
# ╟─3545e4f9-df86-4261-943b-e79a48274716
# ╟─273be1ff-58b0-4dab-9726-80b1fa8aa0c1
# ╠═1af43c9b-151f-4b22-aa88-110ad018e0ad
# ╠═ae2fc32c-a471-4e64-9d33-e9e6bf677166
# ╟─3608f7cf-86e3-4cdb-ada6-46c6beb3310a
# ╟─b6406021-00e6-419f-8e15-9f80196433eb
# ╠═ff3539e9-296f-403d-a031-86c4df8566f7
# ╠═be19363a-dfd3-430e-b1c9-f8b69309411d
# ╟─0c50f149-4cea-4793-acd2-91b58595dc3b
# ╟─a725c129-1f57-47d3-810c-bfee4f600bcf
# ╟─80986d60-d993-4293-bd22-72cdef15988b
# ╟─9bccf94a-19d1-4d31-82d5-d882b1fd1e29
# ╠═add8900c-f7a6-450e-8e04-0e7237cf80f3
# ╠═250726a2-4c56-46fe-a2b1-8f5da0dc761e
# ╠═56a560dc-dfac-43a9-afa2-a3cd6733d14f
# ╟─711dc218-0d78-4a39-a380-bc461e249ba5
# ╟─7bcf336b-c2dc-40ad-bb43-82fbca1c8565
# ╟─e42073de-0d42-42c3-9f1a-53f2ff54b91c
# ╟─52fb0f4d-c1c0-4ca2-9bd1-913414fcf9c8
# ╟─20b697f4-0bf2-426b-ae3b-a3864386cc19
# ╟─6714a575-dc6a-4926-bcfa-349b91aca16e
# ╟─5c66ea8a-ca3e-4d07-b16c-8f4015f87671
# ╟─5efe07fc-b7e4-4069-83c6-b1534c39ac99
# ╠═172cf452-791b-4d7a-934a-0dfa05374e14
# ╟─8d9ca351-c8f6-49bd-8a91-2c34e883d0c4
# ╠═3675a743-1ddc-4cd4-b954-8ca354bea723
# ╟─d533573e-10fd-4480-9b17-93ef3dacb91e
# ╠═2e34f342-f0a8-4e62-a4bc-6bda81c5b0c6
# ╠═6c9df6b6-746a-4721-b77b-52718abab2e5
# ╠═549a0344-a442-4dfd-b01f-f1b06723df7e
# ╠═ab6b4c8b-0b62-4c50-bf2f-f9b969d23aff
# ╟─dd6dff51-8765-48da-9ccc-6fb6ca228bff
# ╠═299f6a19-a805-4ff2-8108-95dd94a5bb6d
# ╟─92da736b-2da2-43bb-9c8c-c661a327acc0
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/00 - Intro.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/00 - Intro.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/00 - Intro.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
#=
Capture image from webcamera.
The original implementation is taken from VideoIO.jl(see __init__ function)
=#

using Makie
using VideoIO

function center_crop(img)
    imW, imH = size(img)
    sz = min(imH, imW)
    y_offset = round(Int,(imH - sz) / 2.)
    x_offset = round(Int,(imW - sz) / 2.)
    y_slice = range(1+y_offset, stop=y_offset+sz)
    x_slice = range(1+x_offset, stop=x_offset+sz)
    cropped = img[x_slice, y_slice]
    cropped
end

function play(f; flipx = false, flipy = false, pixelaspectratio = nothing)
    if pixelaspectratio ≡ nothing # if user did not specify the aspect ratio we'll try to use the one stored in the video file
        pixelaspectratio = VideoIO.aspect_ratio(f)
    end
    sz=min(f.height, f.width)
    scene = Makie.Scene(resolution = (sz,sz))
    buf = read(f)
    disp=center_crop(buf)
    makieimg = Makie.image!(
        scene,
        1:sz,
        1:sz,
        disp,
        show_axis = false,
        scale_plot = false,
    )[end]
    Makie.rotate!(scene, -0.5pi)
    if flipx && flipy
        Makie.scale!(scene, -1, -1, 1)
    else
        flipx && Makie.scale!(scene, -1, 1, 1)
        flipy && Makie.scale!(scene, 1, -1, 1)
    end
    display(scene)
    while !eof(f) && isopen(scene)
        read!(f, buf)
        disp=center_crop(buf)
        makieimg[3] = disp
        sleep(1 / f.framerate)
    end

end

function begindemo(pixelaspectratio = nothing)
    camera = opencamera()
    play(camera, flipx = true, pixelaspectratio = pixelaspectratio)
    close(camera)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
#=
TestInverseMatrices.jl:
- Julia version: 1.0
- Author: qz
- Date: Apr 22, 2019
=#
module TestInverseMatrices

using LinearAlgebra
using Test

using APMAE4001.InverseMatrices

@testset "test 1" begin
    x = Diagonal([1, 2.0, 0.001])
    @test inverse(x) == inv(x)
end

end